sll_m_maxwell_2d_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
 
module sll_m_maxwell_2d_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, only: sll_t_fft, &
                        sll_s_fft_free, &
                        sll_s_fft_init_c2r_1d, &
                        sll_s_fft_init_r2c_1d, &
                        sll_s_fft_exec_c2r_1d, &
                        sll_s_fft_exec_r2c_1d
 
   implicit none
 
   public :: sll_t_maxwell_2d_pstd, &
             sll_s_maxwell_2d_pstd_init, &
             sll_s_solve_maxwell_2d_pstd, &
             sll_s_solve_ampere_2d_pstd, &
             sll_s_solve_faraday_2d_pstd, &
             sll_s_free_maxwell_2d_pstd
 
   private
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
   type :: sll_t_maxwell_2d_pstd
 
      private
      sll_int32           :: nc_eta1      
      sll_int32           :: nc_eta2      
      sll_int32           :: polarization 
      sll_real64          :: e_0          
      sll_real64          :: mu_0         
      sll_real64          :: c            
      sll_real64          :: eta1_min     
      sll_real64          :: eta1_max     
      sll_real64          :: delta_eta1   
      sll_real64          :: eta2_min     
      sll_real64          :: eta2_max     
      sll_real64          :: delta_eta2   
      sll_real64, pointer :: d_dx(:)      
      sll_real64, pointer :: d_dy(:)      
      sll_real64, pointer :: kx(:)        
      sll_real64, pointer :: ky(:)        
      type(sll_t_fft)     :: fwx
      type(sll_t_fft)     :: fwy
      type(sll_t_fft)     :: bwx
      type(sll_t_fft)     :: bwy
      sll_comp64, pointer :: tmp_x(:)     
      sll_comp64, pointer :: tmp_y(:)     
 
   end type sll_t_maxwell_2d_pstd
 
 
 
 
 
contains
 
   subroutine sll_s_maxwell_2d_pstd_init(self, xmin, xmax, nc_x, ymin, ymax, nc_y, polarization)
 
      type(sll_t_maxwell_2d_pstd) :: self
      sll_real64, intent(in)      :: xmin         
      sll_real64, intent(in)      :: xmax         
      sll_real64, intent(in)      :: ymin         
      sll_real64, intent(in)      :: ymax         
      sll_int32, intent(in)      :: nc_x         
      sll_int32, intent(in)      :: nc_y         
      sll_int32, intent(in)      :: polarization 
 
      sll_int32                   :: error        
      sll_real64                  :: dx           
      sll_real64                  :: dy           
      sll_real64                  :: kx0
      sll_real64                  :: ky0
      sll_int32                   :: i, j
 
      self%nc_eta1 = nc_x
      self%nc_eta2 = nc_y
      self%polarization = polarization
 
      self%e_0 = 1._f64
      self%mu_0 = 1._f64
 
      sll_allocate(self%d_dx(nc_x), error)
      sll_allocate(self%d_dy(nc_y), error)
 
      sll_allocate(self%tmp_x(nc_x/2 + 1), error)
      sll_allocate(self%tmp_y(nc_y/2 + 1), 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)
 
      sll_allocate(self%kx(nc_x/2 + 1), error)
      sll_allocate(self%ky(nc_y/2 + 1), error)
 
      dx = (xmax - xmin)/nc_x
      dy = (ymax - ymin)/nc_y
 
      kx0 = 2._f64*sll_p_pi/(nc_x*dx)
      ky0 = 2._f64*sll_p_pi/(nc_y*dy)
 
      do i = 2, nc_x/2 + 1
         self%kx(i) = (i - 1)*kx0
      end do
      self%kx(1) = 1.0_f64
      do j = 2, nc_y/2 + 1
         self%ky(j) = (j - 1)*ky0
      end do
      self%ky(1) = 1.0_f64
 
   end subroutine sll_s_maxwell_2d_pstd_init
 
   subroutine sll_s_solve_maxwell_2d_pstd(self, fx, fy, fz, dt)
 
      type(sll_t_maxwell_2d_pstd), intent(inout)          :: self
      sll_real64, intent(inout), dimension(:, :) :: fx   
      sll_real64, intent(inout), dimension(:, :) :: fy   
      sll_real64, intent(inout), dimension(:, :) :: fz   
      sll_real64, intent(in)                    :: dt   
 
      IF (self%polarization == tm_polarization) then
         call faraday_tm_2d_pstd(self, fx, fy, fz, 0.5*dt)
         call ampere_tm_2d_pstd(self, fx, fy, fz, dt)
         call faraday_tm_2d_pstd(self, fx, fy, fz, 0.5*dt)
      end if
 
      IF (self%polarization == te_polarization) then
         call faraday_te_2d_pstd(self, fx, fy, fz, 0.5*dt)
         call ampere_te_2d_pstd(self, fx, fy, fz, dt)
         call faraday_te_2d_pstd(self, fx, fy, fz, 0.5*dt)
      end if
 
   end subroutine sll_s_solve_maxwell_2d_pstd
 
   subroutine sll_s_solve_faraday_2d_pstd(self, fx, fy, fz, dt)
      type(sll_t_maxwell_2d_pstd), intent(inout) :: self
      sll_real64, dimension(:, :), intent(inout) :: fx      
      sll_real64, dimension(:, :), intent(inout) :: fy      
      sll_real64, dimension(:, :), intent(inout) :: fz      
      sll_real64, intent(in)                   :: dt      
 
      if (self%polarization == tm_polarization) then
         call faraday_tm_2d_pstd(self, fx, fy, fz, dt)
      end if
 
      if (self%polarization == te_polarization) then
         call faraday_te_2d_pstd(self, fx, fy, fz, dt)
      end if
 
   end subroutine sll_s_solve_faraday_2d_pstd
 
   subroutine sll_s_solve_ampere_2d_pstd(self, fx, fy, fz, dt, sx, sy)
      type(sll_t_maxwell_2d_pstd), intent(inout) :: self
      sll_real64, dimension(:, :), intent(inout) :: fx      
      sll_real64, dimension(:, :), intent(inout) :: fy      
      sll_real64, dimension(:, :), intent(inout) :: fz      
      sll_real64, intent(in)                   :: dt      
      sll_real64, dimension(:, :), optional      :: sx      
      sll_real64, dimension(:, :), optional      :: sy      
 
      if (self%polarization == tm_polarization .and. &
          present(sx) .and. .not. present(sy)) then
         call ampere_tm_2d_pstd(self, fx, fy, fz, dt, sx)
      else if (self%polarization == te_polarization .and. &
               present(sx) .and. present(sy)) then
         call ampere_te_2d_pstd(self, fx, fy, fz, dt, sx, sy)
      end if
 
   end subroutine sll_s_solve_ampere_2d_pstd
 
   subroutine faraday_tm_2d_pstd(self, hx, hy, ez, dt)
 
      type(sll_t_maxwell_2d_pstd), intent(inout) :: self
      sll_real64, dimension(:, :), intent(inout) :: hx      
      sll_real64, dimension(:, :), intent(inout) :: hy      
      sll_real64, dimension(:, :), intent(inout) :: ez      
      sll_int32                                 :: nc_x    
      sll_int32                                 :: nc_y    
      sll_real64, intent(in)                    :: dt      
 
      sll_real64                                :: dt_mu
      sll_int32                                 :: i, j
 
      nc_x = self%nc_eta1
      nc_y = self%nc_eta2
 
      dt_mu = dt/self%mu_0
 
      do i = 1, nc_x
         d_dy(ez(i, 1:nc_y))
         hx(i, 1:nc_y) = hx(i, 1:nc_y) - dt_mu*self%d_dy
      end do
 
      do j = 1, nc_y
         d_dx(ez(1:nc_x, j))
         hy(1:nc_x, j) = hy(1:nc_x, j) + dt_mu*self%d_dx
      end do
 
      if (size(hx, 2) == nc_y + 1) hx(:, nc_y + 1) = hx(:, 1)
      if (size(hy, 1) == nc_x + 1) hy(nc_x + 1, :) = hy(1, :)
 
   end subroutine faraday_tm_2d_pstd
 
   subroutine faraday_te_2d_pstd(self, ex, ey, hz, dt)
 
      type(sll_t_maxwell_2d_pstd), intent(inout) :: self
      sll_real64, dimension(:, :), intent(inout) :: ex     
      sll_real64, dimension(:, :), intent(inout) :: ey     
      sll_real64, dimension(:, :), intent(inout) :: hz     
      sll_int32                                 :: nc_x   
      sll_int32                                 :: nc_y   
      sll_real64, intent(in)                    :: dt     
 
      sll_real64                                :: dt_mu
      sll_int32                                 :: i, j
 
      nc_x = self%nc_eta1
      nc_y = self%nc_eta2
 
      dt_mu = dt/self%mu_0
 
      do i = 1, nc_x
         d_dy(ex(i, 1:nc_y))
         hz(i, 1:nc_y) = hz(i, 1:nc_y) + dt_mu*self%d_dy
      end do
 
      do j = 1, nc_y
         d_dx(ey(1:nc_x, j))
         hz(1:nc_x, j) = hz(1:nc_x, j) - dt_mu*self%d_dx
      end do
 
      if (size(hz, 1) == nc_x + 1) hz(nc_x + 1, :) = hz(1, :)
      if (size(hz, 2) == nc_y + 1) hz(:, nc_y + 1) = hz(:, 1)
 
   end subroutine faraday_te_2d_pstd
 
   subroutine ampere_tm_2d_pstd(self, hx, hy, ez, dt, jz)
 
      type(sll_t_maxwell_2d_pstd), intent(inout) :: self
      sll_int32                                 :: nc_x   
      sll_int32                                 :: nc_y   
      sll_real64, dimension(:, :)                :: hx     
      sll_real64, dimension(:, :)                :: hy     
      sll_real64, dimension(:, :)                :: ez     
      sll_real64                                :: dt     
      sll_real64, dimension(:, :), optional      :: jz     
 
      sll_real64                                :: dt_e
      sll_int32                                 :: i, j
 
      nc_x = self%nc_eta1
      nc_y = self%nc_eta2
 
      dt_e = dt/self%e_0
 
      do j = 1, nc_y
         d_dx(hy(1:nc_x, j))
         ez(1:nc_x, j) = ez(1:nc_x, j) + dt_e*self%d_dx
      end do
 
      do i = 1, nc_x
         d_dy(hx(i, 1:nc_y))
         ez(i, 1:nc_y) = ez(i, 1:nc_y) - dt_e*self%d_dy
      end do
 
      if (present(jz)) then
         ez = ez - dt_e*jz
      end if
 
      if (size(ez, 1) == nc_x + 1) ez(nc_x + 1, :) = ez(1, :)
      if (size(ez, 2) == nc_y + 1) ez(:, nc_y + 1) = ez(:, 1)
 
   end subroutine ampere_tm_2d_pstd
 
   subroutine ampere_te_2d_pstd(self, ex, ey, hz, dt, jx, jy)
 
      type(sll_t_maxwell_2d_pstd), intent(inout) :: self
      sll_real64, dimension(:, :)                :: ex     
      sll_real64, dimension(:, :)                :: ey     
      sll_real64, dimension(:, :)                :: hz     
      sll_real64                                :: dt     
      sll_real64, dimension(:, :), optional      :: jx     
      sll_real64, dimension(:, :), optional      :: jy     
 
      sll_real64                                :: dt_e
      sll_int32                                 :: nc_x
      sll_int32                                 :: nc_y
      sll_int32                                 :: i, j
 
      nc_x = self%nc_eta1
      nc_y = self%nc_eta2
 
      dt_e = dt/self%e_0
 
      do j = 1, nc_y
         d_dx(hz(1:nc_x, j))
         ey(1:nc_x, j) = ey(1:nc_x, j) - dt_e*self%d_dx
      end do
 
      do i = 1, nc_x
         d_dy(hz(i, 1:nc_y))
         ex(i, 1:nc_y) = ex(i, 1:nc_y) + dt_e*self%d_dy
      end do
 
      If (present(jx) .and. present(jy)) then
         ex = ex - dt_e*jx
         ey = ey - dt_e*jy
      end if
 
      if (size(ex, 1) == nc_x + 1) ex(nc_x + 1, :) = ex(1, :)
      if (size(ey, 2) == nc_y + 1) ey(:, nc_y + 1) = ey(:, 1)
 
   end subroutine ampere_te_2d_pstd
 
   subroutine sll_s_free_maxwell_2d_pstd(self)
      type(sll_t_maxwell_2d_pstd) :: self
 
      call sll_s_fft_free(self%fwx)
      call sll_s_fft_free(self%fwy)
      call sll_s_fft_free(self%bwx)
      call sll_s_fft_free(self%bwy)
 
   end subroutine sll_s_free_maxwell_2d_pstd
 
end module sll_m_maxwell_2d_pstd