sll_m_poisson_3d_periodic.F90

Go to the documentation of this file.

 
 
module sll_m_poisson_3d_periodic
 
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#include "sll_assert.h"
#include "sll_memory.h"
#include "sll_working_precision.h"
 
   use sll_m_constants, only: &
      sll_p_pi
 
   use sll_m_fft, only: &
      sll_s_fft_exec_c2c_3d, &
      sll_p_fft_backward, &
      sll_s_fft_free, &
      sll_p_fft_forward, &
      sll_s_fft_init_c2c_3d, &
      sll_t_fft
 
   implicit none
 
   public :: &
      sll_s_poisson_3d_periodic_init, &
      sll_s_poisson_3d_periodic_free, &
      sll_t_poisson_3d_periodic, &
      sll_s_poisson_3d_periodic_solve
 
   private
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
   type sll_t_poisson_3d_periodic
      sll_int32                             :: nx      
      sll_int32                             :: ny      
      sll_int32                             :: nz      
      sll_real64                            :: lx      
      sll_real64                            :: ly      
      sll_real64                            :: lz      
      type(sll_t_fft)                       :: fw
      type(sll_t_fft)                       :: bw
      sll_comp64, dimension(:, :, :), pointer :: hat_rho 
      sll_comp64, dimension(:, :, :), pointer :: hat_phi 
   end type sll_t_poisson_3d_periodic
 
contains
 
   subroutine sll_s_poisson_3d_periodic_init(self, nx, ny, nz, Lx, Ly, Lz)
 
      sll_int32                              :: nx   
      sll_int32                              :: ny   
      sll_int32                              :: nz   
      sll_int32                              :: ierr
      sll_real64                             :: lx   
      sll_real64                             :: ly   
      sll_real64                             :: lz   
      type(sll_t_poisson_3d_periodic)       :: self
 
      sll_allocate(self%hat_rho(1:nx, 1:ny, 1:nz), ierr)
      sll_allocate(self%hat_phi(1:nx, 1:ny, 1:nz), ierr)
 
      ! Geometry informations
      self%nx = nx
      self%ny = ny
      self%nz = nz
      self%Lx = lx
      self%Ly = ly
      self%Lz = lz
 
      call sll_s_fft_init_c2c_3d(self%fw, nx, ny, nz, self%hat_rho, self%hat_phi, sll_p_fft_forward)
      call sll_s_fft_init_c2c_3d(self%bw, nx, ny, nz, self%hat_phi, self%hat_rho, sll_p_fft_backward)
 
   end subroutine sll_s_poisson_3d_periodic_init
 
   subroutine sll_s_poisson_3d_periodic_solve(self, rho, phi)
 
      type(sll_t_poisson_3d_periodic) :: self
      sll_real64, dimension(:, :, :)     :: rho  
      sll_real64, dimension(:, :, :)     :: phi  
      sll_int32                        :: nx, ny, nz
      sll_int32                        :: i, j, k
      sll_real64                       :: lx, ly, lz
      sll_real64                       :: ind_x, ind_y, ind_z
 
      nx = self%nx
      ny = self%ny
      nz = self%nz
      lx = self%Lx
      ly = self%Ly
      lz = self%Lz
 
      sll_assert(nx == size(rho, 1))
      sll_assert(ny == size(rho, 2))
      sll_assert(nz == size(rho, 3))
      sll_assert(nx == size(phi, 1))
      sll_assert(ny == size(phi, 2))
      sll_assert(nz == size(phi, 3))
 
      self%hat_phi = cmplx(rho, 0_f64, kind=f64)
      call sll_s_fft_exec_c2c_3d(self%fw, self%hat_phi, self%hat_rho)
      self%hat_rho = self%hat_rho/cmplx(nx*ny*nz, 0.0, f64)
 
      ! Compute hat_phi, phi = inv_fft(hat_phi)
      do k = 1, nz
         do j = 1, ny
            do i = 1, nx
               if (i <= nx/2) then
                  ind_x = real(i - 1, f64)
               else
                  ind_x = real(nx - (i - 1), f64)
               end if
               if (j <= ny/2) then
                  ind_y = real(j - 1, f64)
               else
                  ind_y = real(ny - (j - 1), f64)
               end if
               if (k <= nz/2) then
                  ind_z = real(k - 1, f64)
               else
                  ind_z = real(nz - (k - 1), f64)
               end if
               if ((ind_x == 0) .and. (ind_y == 0) .and. (ind_z == 0)) then
                  self%hat_rho(i, j, k) = (0._f64, 0._f64)
               else
                  self%hat_rho(i, j, k) = self%hat_rho(i, j, k)/ &
                                          (4.0_f64*sll_p_pi**2* &
                                           ((real(ind_x, f64)/lx)**2 &
                                            + (real(ind_y, f64)/ly)**2 &
                                            + (real(ind_z, f64)/lz)**2))
               end if
            end do
         end do
      end do
 
      call sll_s_fft_exec_c2c_3d(self%bw, self%hat_rho, self%hat_phi)
 
      phi = real(self%hat_phi, f64)
 
   end subroutine sll_s_poisson_3d_periodic_solve
 
   subroutine sll_s_poisson_3d_periodic_free(self)
 
      type(sll_t_poisson_3d_periodic) :: self
      sll_int32                       :: ierr
 
      call sll_s_fft_free(self%fw)
      call sll_s_fft_free(self%bw)
 
      sll_deallocate_array(self%hat_rho, ierr)
      sll_deallocate_array(self%hat_phi, ierr)
 
   end subroutine sll_s_poisson_3d_periodic_free
 
end module sll_m_poisson_3d_periodic