sll_m_maxwell_3d_pstd.F90

Go to the documentation of this file.

!**************************************************************
!  Copyright INRIA
!
!  This code SeLaLib (for Semi-Lagrangian-Library)
!  is a parallel library for simulating the plasma turbulence
!  in a tokamak.
!
!  This software is governed by the CeCILL-B license
!  under French law and abiding by the rules of distribution
!  of free software.  You can  use, modify and redistribute
!  the software under the terms of the CeCILL-B license as
!  circulated by CEA, CNRS and INRIA at the following URL
!  "http://www.cecill.info".
!**************************************************************
 
#include "sll_fftw.h"
 
#define D_DX(field)                                           \
self%d_dx = field; \
call sll_s_fft_exec_r2c_1d(self%fwx, self%d_dx, self%tmp_x); \
self%tmp_x(2:nc_x/2 + 1) = -cmplx(0.0_f64, self%kx(2:nc_x/2 + 1), kind=f64)*self%tmp_x(2:nc_x/2 + 1); \
call sll_s_fft_exec_c2r_1d(self%bwx, self%tmp_x, self%d_dx); \
self%d_dx = self%d_dx/nc_x
 
#define D_DY(field)                                           \
self%d_dy = field; \
call sll_s_fft_exec_r2c_1d(self%fwy, self%d_dy, self%tmp_y); \
self%tmp_y(2:nc_y/2 + 1) = -cmplx(0.0_f64, self%ky(2:nc_y/2 + 1), kind=f64)*self%tmp_y(2:nc_y/2 + 1); \
call sll_s_fft_exec_c2r_1d(self%bwy, self%tmp_y, self%d_dy); \
self%d_dy = self%d_dy/nc_y
 
#define D_DZ(field)                                           \
self%d_dz = field; \
call sll_s_fft_exec_r2c_1d(self%fwz, self%d_dz, self%tmp_z); \
self%tmp_z(2:nc_z/2 + 1) = -cmplx(0.0_f64, self%kz(2:nc_z/2 + 1), kind=f64)*self%tmp_z(2:nc_z/2 + 1); \
call sll_s_fft_exec_c2r_1d(self%bwz, self%tmp_z, self%d_dz); \
self%d_dz = self%d_dz/nc_z
 
module sll_m_maxwell_3d_pstd
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#include "sll_memory.h"
#include "sll_working_precision.h"
#include "sll_maxwell_solvers_macros.h"
 
   use sll_m_constants, only: sll_p_pi
 
   use sll_m_fft
 
   implicit none
 
   type, public :: sll_t_maxwell_3d_pstd
      private
      sll_int32                          :: nc_x         
      sll_int32                          :: nc_y         
      sll_int32                          :: nc_z         
      sll_real64, dimension(:), pointer  :: d_dx         
      sll_real64, dimension(:), pointer  :: d_dy         
      sll_real64, dimension(:), pointer  :: d_dz         
      sll_real64, dimension(:), pointer  :: kx           
      sll_real64, dimension(:), pointer  :: ky           
      sll_real64, dimension(:), pointer  :: kz           
      type(sll_t_fft)                    :: fwx
      type(sll_t_fft)                    :: fwy
      type(sll_t_fft)                    :: fwz
      type(sll_t_fft)                    :: bwx
      type(sll_t_fft)                    :: bwy
      type(sll_t_fft)                    :: bwz
      sll_comp64, pointer :: tmp_x(:)     
      sll_comp64, pointer :: tmp_y(:)     
      sll_comp64, pointer :: tmp_z(:)     
      sll_int32                          :: polarization 
      sll_real64                         :: e_0          
      sll_real64                         :: mu_0         
   end type sll_t_maxwell_3d_pstd
 
contains
 
   subroutine sll_s_maxwell_3d_pstd_init(self, xmin, xmax, nc_x, &
                                         ymin, ymax, nc_y, &
                                         zmin, zmax, nc_z)
 
      type(sll_t_maxwell_3d_pstd) :: self
      sll_real64, intent(in)      :: xmin   
      sll_real64, intent(in)      :: xmax   
      sll_real64, intent(in)      :: ymin   
      sll_real64, intent(in)      :: ymax   
      sll_real64, intent(in)      :: zmin   
      sll_real64, intent(in)      :: zmax   
      sll_int32, intent(in)      :: nc_x   
      sll_int32, intent(in)      :: nc_y   
      sll_int32, intent(in)      :: nc_z   
 
      sll_int32                   :: error  
      sll_real64                  :: dx     
      sll_real64                  :: dy     
      sll_real64                  :: dz     
      sll_real64                  :: kx0
      sll_real64                  :: ky0
      sll_real64                  :: kz0
 
      sll_int32                   :: i, j, k
 
      self%nc_x = nc_x
      self%nc_y = nc_y
      self%nc_z = nc_z
 
      self%e_0 = 1.0_f64
      self%mu_0 = 1.0_f64
 
      sll_allocate(self%tmp_x(1:nc_x/2 + 1), error); self%tmp_x = (0.0_f64, 0.0_f64)
      sll_allocate(self%tmp_y(1:nc_y/2 + 1), error); self%tmp_y = (0.0_f64, 0.0_f64)
      sll_allocate(self%tmp_z(1:nc_z/2 + 1), error); self%tmp_z = (0.0_f64, 0.0_f64)
 
      sll_clear_allocate(self%d_dx(1:nc_x), error)
      sll_clear_allocate(self%d_dy(1:nc_y), error)
      sll_clear_allocate(self%d_dz(1:nc_z), error)
 
      call sll_s_fft_init_r2c_1d(self%fwx, nc_x, self%d_dx, self%tmp_x)
      call sll_s_fft_init_c2r_1d(self%bwx, nc_x, self%tmp_x, self%d_dx)
      call sll_s_fft_init_r2c_1d(self%fwy, nc_y, self%d_dy, self%tmp_y)
      call sll_s_fft_init_c2r_1d(self%bwy, nc_y, self%tmp_y, self%d_dy)
      call sll_s_fft_init_r2c_1d(self%fwz, nc_z, self%d_dz, self%tmp_z)
      call sll_s_fft_init_c2r_1d(self%bwz, nc_z, self%tmp_z, self%d_dz)
 
      sll_allocate(self%kx(nc_x/2 + 1), error)
      sll_allocate(self%ky(nc_y/2 + 1), error)
      sll_allocate(self%kz(nc_z/2 + 1), error)
 
      dx = (xmax - xmin)/real(nc_x,f64)
      dy = (ymax - ymin)/real(nc_y,f64)
      dz = (zmax - zmin)/real(nc_z,f64)
 
      kx0 = 2._f64*sll_p_pi/(real(nc_x,f64)*dx)
      ky0 = 2._f64*sll_p_pi/(real(nc_y,f64)*dy)
      kz0 = 2._f64*sll_p_pi/(real(nc_z,f64)*dz)
 
      self%kx(1) = 1.0_f64
      do i = 2, nc_x/2 + 1
         self%kx(i) = real(i - 1, f64)*kx0
      end do
      self%ky(1) = 1.0_f64
      do j = 2, nc_y/2 + 1
         self%ky(j) = real(j - 1, f64)*ky0
      end do
      self%kz(1) = 1.0_f64
      do k = 2, nc_z/2 + 1
         self%kz(k) = real(k - 1, f64)*kz0
      end do
 
   end subroutine sll_s_maxwell_3d_pstd_init
 
   subroutine sll_s_maxwell_3d_pstd_faraday(self, ex, ey, ez, hx, hy, hz, dt)
 
      type(sll_t_maxwell_3d_pstd), intent(inout) :: self
      sll_real64, dimension(:, :, :), intent(inout) :: ex    
      sll_real64, dimension(:, :, :), intent(inout) :: ey    
      sll_real64, dimension(:, :, :), intent(inout) :: ez    
      sll_real64, dimension(:, :, :), intent(inout) :: hx    
      sll_real64, dimension(:, :, :), intent(inout) :: hy    
      sll_real64, dimension(:, :, :), intent(inout) :: hz    
      sll_int32                                   :: nc_x  
      sll_int32                                   :: nc_y  
      sll_int32                                   :: nc_z  
      sll_real64, intent(in)                      :: dt    
 
      sll_real64 :: dt_mu
      sll_int32  :: i, j, k
 
      nc_x = self%nc_x
      nc_y = self%nc_y
      nc_z = self%nc_z
 
      dt_mu = dt/self%mu_0
 
      do k = 1, nc_z + 1
         do i = 1, nc_x + 1
            d_dy(ez(i, 1:nc_y, k))
            hx(i, 1:nc_y, k) = hx(i, 1:nc_y, k) - dt_mu*self%d_dy
            d_dy(ex(i, 1:nc_y, k))
            hz(i, 1:nc_y, k) = hz(i, 1:nc_y, k) + dt_mu*self%d_dy
         end do
      end do
 
      hx(:, nc_y + 1, :) = hx(:, 1, :)
      hz(:, nc_y + 1, :) = hz(:, 1, :)
 
      do j = 1, nc_y + 1
         do i = 1, nc_x + 1
            d_dz(ey(i, j, 1:nc_z))
            hx(i, j, 1:nc_z) = hx(i, j, 1:nc_z) + dt_mu*self%d_dz
            d_dz(ex(i, j, 1:nc_z))
            hy(i, j, 1:nc_z) = hy(i, j, 1:nc_z) - dt_mu*self%d_dz
         end do
      end do
 
      hx(:, :, nc_z + 1) = hx(:, :, 1)
      hy(:, :, nc_z + 1) = hy(:, :, 1)
 
      do k = 1, nc_z + 1
         do j = 1, nc_y + 1
            d_dx(ez(1:nc_x, j, k))
            hy(1:nc_x, j, k) = hy(1:nc_x, j, k) + dt_mu*self%d_dx
            d_dx(ey(1:nc_x, j, k))
            hz(1:nc_x, j, k) = hz(1:nc_x, j, k) - dt_mu*self%d_dx
         end do
      end do
 
      hy(nc_x + 1, :, :) = hy(1, :, :)
      hz(nc_x + 1, :, :) = hz(1, :, :)
 
   end subroutine sll_s_maxwell_3d_pstd_faraday
 
   subroutine sll_s_maxwell_3d_pstd_ampere(self, hx, hy, hz, ex, ey, ez, dt, jx, jy, jz)
 
      type(sll_t_maxwell_3d_pstd), intent(inout) :: self
      sll_real64, dimension(:, :, :)              :: hx   
      sll_real64, dimension(:, :, :)              :: hy   
      sll_real64, dimension(:, :, :)              :: hz   
      sll_real64, dimension(:, :, :)              :: ex   
      sll_real64, dimension(:, :, :)              :: ey   
      sll_real64, dimension(:, :, :)              :: ez   
      sll_real64                                :: dt   
      sll_real64, dimension(:, :, :), optional    :: jx   
      sll_real64, dimension(:, :, :), optional    :: jy   
      sll_real64, dimension(:, :, :), optional    :: jz   
 
      sll_int32                                 :: nc_x 
      sll_int32                                 :: nc_y 
      sll_int32                                 :: nc_z 
 
      sll_real64 :: dt_e
      sll_int32  :: i, j, k
 
      nc_x = self%nc_x
      nc_y = self%nc_y
      nc_z = self%nc_z
 
      dt_e = dt/self%e_0
 
      do k = 1, nc_z + 1
      do i = 1, nc_x + 1
         d_dy(hz(i, 1:nc_y, k))
         ex(i, 1:nc_y, k) = ex(i, 1:nc_y, k) + dt_e*self%d_dy
         d_dy(hx(i, 1:nc_y, k))
         ez(i, 1:nc_y, k) = ez(i, 1:nc_y, k) - dt_e*self%d_dy
      end do
      end do
 
      ex(:, nc_y + 1, :) = ex(:, 1, :)
      ez(:, nc_y + 1, :) = ez(:, 1, :)
 
      do j = 1, nc_y + 1
      do i = 1, nc_x + 1
         d_dz(hy(i, j, 1:nc_z))
         ex(i, j, 1:nc_z) = ex(i, j, 1:nc_z) - dt_e*self%d_dz
         d_dz(hx(i, j, 1:nc_z))
         ey(i, j, 1:nc_z) = ey(i, j, 1:nc_z) + dt_e*self%d_dz
      end do
      end do
 
      ex(:, :, nc_z + 1) = ex(:, :, 1)
      ey(:, :, nc_z + 1) = ey(:, :, 1)
 
      do k = 1, nc_z + 1
      do j = 1, nc_y + 1
         d_dx(hz(1:nc_x, j, k))
         ey(1:nc_x, j, k) = ey(1:nc_x, j, k) - dt_e*self%d_dx
         d_dx(hy(1:nc_x, j, k))
         ez(1:nc_x, j, k) = ez(1:nc_x, j, k) + dt_e*self%d_dx
      end do
      end do
 
      ey(nc_x + 1, :, :) = ey(1, :, :)
      ez(nc_x + 1, :, :) = ez(1, :, :)
 
      if (present(jx) .and. present(jy) .and. present(jz)) then
         ex = ex - dt_e*jx
         ey = ey - dt_e*jy
         ez = ez - dt_e*jz
      end if
 
   end subroutine sll_s_maxwell_3d_pstd_ampere
 
   subroutine sll_s_maxwell_3d_pstd_free(self)
 
      type(sll_t_maxwell_3d_pstd) :: self
 
      call sll_s_fft_free(self%fwx)
      call sll_s_fft_free(self%fwy)
      call sll_s_fft_free(self%fwz)
      call sll_s_fft_free(self%bwx)
      call sll_s_fft_free(self%bwy)
      call sll_s_fft_free(self%bwz)
 
   end subroutine sll_s_maxwell_3d_pstd_free
 
end module sll_m_maxwell_3d_pstd