sll_m_maxwell_2d_periodic_cartesian_par.F90

Go to the documentation of this file.

!**************************************************************
!  Copyright INRIA, CNRS
!
!  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".
!**************************************************************
 
#define D_DX(field)                                                   \
plan%d_dx = field; \
call sll_s_fft_exec_r2c_1d(plan%fwx, plan%d_dx, plan%fft_x_array); \
plan%fft_x_array(2:ncx/2 + 1) = -cmplx(0.0_f64, plan%kx(2:ncx/2 + 1), kind=f64)*plan%fft_x_array(2:ncx/2 + 1); \
call sll_s_fft_exec_c2r_1d(plan%bwx, plan%fft_x_array, plan%d_dx); \
plan%d_dx = plan%d_dx/plan%ncx
 
#define D_DY(field)                                                   \
plan%d_dy = field; \
call sll_s_fft_exec_r2c_1d(plan%fwy, plan%d_dy, plan%fft_y_array); \
plan%fft_y_array(2:ncy/2 + 1) = -cmplx(0.0_f64, plan%ky(2:ncy/2 + 1), kind=f64)*plan%fft_y_array(2:ncy/2 + 1); \
call sll_s_fft_exec_c2r_1d(plan%bwy, plan%fft_y_array, plan%d_dy); \
plan%d_dy = plan%d_dy/plan%ncy
 
#define MPI_MASTER 0
 
module sll_m_maxwell_2d_periodic_cartesian_par
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#include "sll_memory.h"
#include "sll_working_precision.h"
#include "sll_maxwell_solvers_macros.h"
 
   use sll_m_collective, only: &
      sll_f_get_collective_rank, &
      sll_f_get_collective_size, &
      sll_v_world_collective
 
   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, &
      sll_s_fft_free
 
   use sll_m_remapper, only: &
      sll_o_apply_remap_2d, &
      sll_o_compute_local_sizes, &
      sll_t_layout_2d, &
      sll_o_new_remap_plan, &
      sll_t_remap_plan_2d_real64, &
      sll_o_delete
 
   use sll_m_utilities, only: &
      sll_f_is_power_of_two
 
   implicit none
 
   public :: &
      sll_s_ampere_te, &
      sll_s_delete_maxwell_2d_periodic_plan_cartesian_par, &
      sll_s_faraday_te, &
      sll_t_maxwell_2d_periodic_plan_cartesian_par, &
      sll_f_new_maxwell_2d_periodic_plan_cartesian_par
 
   private
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
   type sll_t_maxwell_2d_periodic_plan_cartesian_par
      sll_int32                                 :: ncx         
      sll_int32                                 :: ncy         
      type(sll_t_fft)                           :: fwx
      type(sll_t_fft)                           :: fwy
      type(sll_t_fft)                           :: bwx
      type(sll_t_fft)                           :: bwy
      sll_real64, dimension(:), pointer         :: d_dx        
      sll_real64, dimension(:), pointer         :: d_dy        
      sll_comp64, dimension(:), pointer         :: fft_x_array 
      sll_comp64, dimension(:), pointer         :: fft_y_array 
      sll_real64                                :: e_0         
      sll_real64                                :: mu_0        
      type(sll_t_layout_2d), pointer :: layout_x
      type(sll_t_layout_2d), pointer :: layout_y
      type(sll_t_remap_plan_2d_real64), pointer :: rmp_xy
      type(sll_t_remap_plan_2d_real64), pointer :: rmp_yx
      sll_real64, dimension(:, :), pointer :: fz_x        
      sll_real64, dimension(:, :), pointer :: fz_y        
      sll_real64, dimension(:), pointer :: kx          
      sll_real64, dimension(:), pointer :: ky          
   end type sll_t_maxwell_2d_periodic_plan_cartesian_par
 
contains
 
   function sll_f_new_maxwell_2d_periodic_plan_cartesian_par( &
      layout_x, layout_y, ncx, ncy, Lx, Ly) result(plan)
 
      type(sll_t_maxwell_2d_periodic_plan_cartesian_par), pointer :: plan
      type(sll_t_layout_2d), pointer :: layout_x
      type(sll_t_layout_2d), pointer :: layout_y
      sll_int32                :: ncx      
      sll_int32                :: ncy      
      sll_real64               :: lx       
      sll_real64               :: ly       
      sll_int32                :: prank    
      sll_int32                :: psize    
      sll_int32                :: error    
      sll_int32                :: nx_loc   
      sll_int32                :: ny_loc   
      sll_real64               :: kx0
      sll_real64               :: ky0
      sll_int32                :: i, j
 
      prank = sll_f_get_collective_rank(sll_v_world_collective)
      psize = sll_f_get_collective_size(sll_v_world_collective)
 
      if ((.not. sll_f_is_power_of_two(int(ncx, i64))) .and. &
          (.not. sll_f_is_power_of_two(int(ncy, i64)))) then
         print *, 'This test needs to run with numbers of cells which are', &
            'powers of 2.'
         print *, 'Exiting...'
         stop
      end if
 
      sll_allocate(plan, error)
      sll_clear_allocate(plan%d_dx(1:ncx), error)
      sll_clear_allocate(plan%d_dy(1:ncy), error)
 
      plan%ncx = ncx
      plan%ncy = ncy
 
      sll_allocate(plan%fft_x_array(ncx/2 + 1), error)
      sll_allocate(plan%fft_y_array(ncy/2 + 1), error)
 
      call sll_s_fft_init_r2c_1d(plan%fwx, ncx, plan%d_dx, plan%fft_x_array)
      call sll_s_fft_init_c2r_1d(plan%bwx, ncx, plan%fft_x_array, plan%d_dx)
      call sll_s_fft_init_r2c_1d(plan%fwy, ncy, plan%d_dy, plan%fft_y_array)
      call sll_s_fft_init_c2r_1d(plan%bwy, ncy, plan%fft_y_array, plan%d_dy)
 
      ! Layout and local sizes for FFTs in x-direction
      plan%layout_x => layout_x
      call sll_o_compute_local_sizes(plan%layout_x, nx_loc, ny_loc)
      sll_clear_allocate(plan%fz_x(1:nx_loc, 1:ny_loc), error)
 
      ! Layout and local sizes for FFTs in y-direction
      plan%layout_y => layout_y
      call sll_o_compute_local_sizes(plan%layout_y, nx_loc, ny_loc)
      sll_clear_allocate(plan%fz_y(1:nx_loc, 1:ny_loc), error)
 
      plan%rmp_xy => sll_o_new_remap_plan(plan%layout_x, plan%layout_y, plan%fz_x)
      plan%rmp_yx => sll_o_new_remap_plan(plan%layout_y, plan%layout_x, plan%fz_y)
 
      sll_allocate(plan%kx(ncx/2 + 1), error)
      sll_allocate(plan%ky(ncy/2 + 1), error)
 
      kx0 = 2._f64*sll_p_pi/lx
      ky0 = 2._f64*sll_p_pi/ly
 
      do i = 2, ncx/2 + 1
         plan%kx(i) = (i - 1)*kx0
      end do
      plan%kx(1) = 1.0_f64
      do j = 2, ncy/2 + 1
         plan%ky(j) = (j - 1)*ky0
      end do
      plan%ky(1) = 1.0_f64
 
   end function sll_f_new_maxwell_2d_periodic_plan_cartesian_par
 
!********************************************************************************
 
   subroutine sll_s_faraday_te(plan, dt, ex, ey)
 
      type(sll_t_maxwell_2d_periodic_plan_cartesian_par), pointer :: plan
 
      sll_real64, dimension(:, :) :: ex       
      sll_real64, dimension(:, :) :: ey       
      sll_int32                  :: ncx      
      sll_int32                  :: ncy      
      sll_int32                  :: nx_loc   
      sll_int32                  :: ny_loc   
      sll_int32                  :: prank
      sll_int32                  :: psize
      sll_real64, intent(in)     :: dt       
 
      sll_real64                 :: dt_mu
      sll_int32                  :: i, j
 
      prank = sll_f_get_collective_rank(sll_v_world_collective)
      psize = sll_f_get_collective_size(sll_v_world_collective)
 
      ncx = plan%ncx
      ncy = plan%ncy
 
#ifdef DEBUG
      call verify_argument_sizes_par(plan%layout_x, ey)
      call verify_argument_sizes_par(plan%layout_y, ex)
#endif
 
      dt_mu = dt/plan%mu_0
 
      call sll_o_compute_local_sizes(plan%layout_x, nx_loc, ny_loc)
      do j = 1, ny_loc
         d_dx(ey(1:ncx, j))
         plan%fz_x(:, j) = plan%fz_x(:, j) - dt_mu*plan%d_dx
      end do
 
      call sll_o_apply_remap_2d(plan%rmp_xy, plan%fz_x, plan%fz_y)
 
      call sll_o_compute_local_sizes(plan%layout_y, nx_loc, ny_loc)
      do i = 1, nx_loc
         d_dy(ex(i, 1:ncy))
         plan%fz_y(i, :) = plan%fz_y(i, :) + dt_mu*plan%d_dy
      end do
 
   end subroutine sll_s_faraday_te
 
   subroutine sll_s_ampere_te(plan, dt, ex, ey, jx, jy)
 
      type(sll_t_maxwell_2d_periodic_plan_cartesian_par), pointer :: plan
 
      sll_real64, dimension(:, :), intent(inout) :: ex     
      sll_real64, dimension(:, :), intent(inout) :: ey     
      sll_real64, dimension(:, :), optional      :: jx     
      sll_real64, dimension(:, :), optional      :: jy     
      sll_real64, intent(in)                    :: dt     
      sll_int32                                 :: nx_loc 
      sll_int32                                 :: ny_loc 
      sll_int32                                 :: prank
      sll_int32                                 :: psize
 
      sll_real64                                :: dt_e
      sll_int32                                 :: i, j
      sll_int32                                 :: ncx, ncy
 
      prank = sll_f_get_collective_rank(sll_v_world_collective)
      psize = sll_f_get_collective_size(sll_v_world_collective)
 
#ifdef DEBUG
      call verify_argument_sizes_par(plan%layout_x, ey)
      call verify_argument_sizes_par(plan%layout_y, ex)
#endif
 
      dt_e = dt/plan%e_0
      ncx = plan%ncx
      ncy = plan%ncy
 
      call sll_o_compute_local_sizes(plan%layout_y, nx_loc, ny_loc)
      do i = 1, nx_loc
         d_dy(plan%fz_y(i, 1:ncy))
         ex(i, :) = ex(i, :) + dt_e*plan%d_dy
      end do
 
      if (present(jx)) then
#ifdef DEBUG
         call verify_argument_sizes_par(plan%layout_y, jx)
#endif
         ex(:, :) = ex(:, :) - dt_e*jx(:, :)
      end if
 
      call sll_o_apply_remap_2d(plan%rmp_xy, plan%fz_x, plan%fz_y)
 
      call sll_o_compute_local_sizes(plan%layout_x, nx_loc, ny_loc)
      do j = 1, ny_loc
         d_dx(plan%fz_x(1:ncx, j))
         ey(:, j) = ey(:, j) - dt_e*plan%d_dx
      end do
 
      if (present(jy)) then
#ifdef DEBUG
         call verify_argument_sizes_par(plan%layout_x, jy)
#endif
         ey(:, :) = ey(:, :) - dt*jy(:, :)/plan%e_0
      end if
 
   end subroutine sll_s_ampere_te
 
   subroutine ampere_tm(plan, dt, bx, by, jz)
 
      type(sll_t_maxwell_2d_periodic_plan_cartesian_par), pointer :: plan
 
      sll_real64, dimension(:, :) :: bx           
      sll_real64, dimension(:, :) :: by           
      sll_real64, dimension(:, :), optional :: jz 
      sll_int32                  :: ncx          
      sll_int32                  :: ncy          
      sll_int32                  :: nx_loc       
      sll_int32                  :: ny_loc       
      sll_int32                  :: prank
      sll_int32                  :: psize
      sll_real64, intent(in)     :: dt           
 
      sll_real64                 :: dt_e
      sll_int32                  :: i, j
 
      prank = sll_f_get_collective_rank(sll_v_world_collective)
      psize = sll_f_get_collective_size(sll_v_world_collective)
 
      ncx = plan%ncx
      ncy = plan%ncy
 
#ifdef DEBUG
      call verify_argument_sizes_par(plan%layout_x, by)
      call verify_argument_sizes_par(plan%layout_y, bx)
#endif
 
      dt_e = dt/plan%e_0
 
      call sll_o_compute_local_sizes(plan%layout_x, nx_loc, ny_loc)
      do j = 1, ny_loc
         d_dx(by(1:ncx, j))
         plan%fz_x(:, j) = plan%fz_x(:, j) + dt_e*plan%d_dx
      end do
 
      call sll_o_apply_remap_2d(plan%rmp_xy, plan%fz_x, plan%fz_y)
 
      call sll_o_compute_local_sizes(plan%layout_y, nx_loc, ny_loc)
      do i = 1, nx_loc
         d_dy(bx(i, 1:ncy))
         plan%fz_y(i, :) = plan%fz_y(i, :) - dt_e*plan%d_dy
      end do
 
      if (present(jz)) then
#ifdef DEBUG
         call verify_argument_sizes_par(plan%layout_y, jz)
#endif
         plan%fz_y = plan%fz_y - dt_e*jz
      end if
 
   end subroutine ampere_tm
 
   subroutine faraday_tm(plan, dt, bx, by)
 
      type(sll_t_maxwell_2d_periodic_plan_cartesian_par), pointer :: plan
 
      sll_real64, dimension(:, :), intent(inout) :: bx     
      sll_real64, dimension(:, :), intent(inout) :: by     
      sll_real64, intent(in)                    :: dt     
      sll_int32                                 :: nx_loc 
      sll_int32                                 :: ny_loc 
      sll_int32                                 :: prank
      sll_int32                                 :: psize
 
      sll_real64                                :: dt_mu
      sll_int32                                 :: i, j
      sll_int32                                 :: ncx, ncy
 
      prank = sll_f_get_collective_rank(sll_v_world_collective)
      psize = sll_f_get_collective_size(sll_v_world_collective)
 
#ifdef DEBUG
      call verify_argument_sizes_par(plan%layout_x, by)
      call verify_argument_sizes_par(plan%layout_y, bx)
#endif
 
      dt_mu = dt/plan%mu_0
 
      ncx = plan%ncx
      ncy = plan%ncy
 
      call sll_o_compute_local_sizes(plan%layout_y, nx_loc, ny_loc)
      do i = 1, nx_loc
         d_dy(plan%fz_y(i, 1:ncy))
         bx(i, :) = bx(i, :) - dt_mu*plan%d_dy
      end do
 
      call sll_o_apply_remap_2d(plan%rmp_xy, plan%fz_x, plan%fz_y)
 
      call sll_o_compute_local_sizes(plan%layout_x, nx_loc, ny_loc)
      do j = 1, ny_loc
         d_dx(plan%fz_x(1:ncx, j))
         by(:, j) = by(:, j) + dt_mu*plan%d_dx
      end do
 
   end subroutine faraday_tm
 
 
   subroutine sll_s_delete_maxwell_2d_periodic_plan_cartesian_par(plan)
      type(sll_t_maxwell_2d_periodic_plan_cartesian_par), pointer :: plan
 
      if (.not. associated(plan)) then
         print *, 'ERROR, delete_maxwell_3d_periodic_plan_par(): ', &
            'passed plan is not associated.'
         stop
      end if
 
      call sll_s_fft_free(plan%fwx)
      call sll_s_fft_free(plan%fwy)
      call sll_s_fft_free(plan%bwx)
      call sll_s_fft_free(plan%bwy)
 
      call sll_o_delete(plan%rmp_xy)
      call sll_o_delete(plan%rmp_yx)
      call sll_o_delete(plan%layout_x)
      call sll_o_delete(plan%layout_y)
 
   end subroutine sll_s_delete_maxwell_2d_periodic_plan_cartesian_par
 
   subroutine verify_argument_sizes_par(layout, array)
      type(sll_t_layout_2d), pointer       :: layout
      sll_real64, dimension(:, :)     :: array  
      sll_int32, dimension(2)       :: n      
      sll_int32                      :: i
 
      call sll_o_compute_local_sizes(layout, n(1), n(2))
 
      do i = 1, 2
         if (n(i) /= size(array, i)) then
            print *, 'ERROR: solve_maxwell_2d_periodic_cartesian_par()', &
               'size of either ex,ey or bz does not match expected size. '
            if (i == 1) then
               print *, 'solve_maxwell_2d_periodic_cartesian_par(): ', &
                  'mismatch in direction x'
            elseif (i == 2) then
               print *, 'solve_maxwell_2d_periodic_cartesian_par(): ', &
                  'mismatch in direction y'
            end if
            print *, 'Exiting...'
            stop
         end if
      end do
 
   end subroutine verify_argument_sizes_par
 
end module sll_m_maxwell_2d_periodic_cartesian_par