sll_m_poisson_1d_periodic.F90

Go to the documentation of this file.

!**************************************************************
!  Copyright INRIA
!  Authors :
!     CALVI project team
!
!  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".
!**************************************************************
 
module sll_m_poisson_1d_periodic
 
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#include "sll_assert.h"
#include "sll_memory.h"
#include "sll_working_precision.h"
#include "sll_errors.h"
 
   use sll_m_constants, only: &
      sll_p_pi
   use sll_m_poisson_1d_base, only: &
      sll_c_poisson_1d_base
 
   implicit none
 
   public :: &
      sll_o_initialize, &
      sll_o_new, &
      sll_t_poisson_1d_periodic, &
      sll_c_poisson_1d_periodic, &
      sll_o_solve, &
      sll_f_new_poisson_1d_periodic
 
   private
 
   type :: sll_t_poisson_1d_periodic
      sll_int32                         :: nc_eta1 
      sll_real64                        :: eta1_min 
      sll_real64                        :: eta1_max 
      sll_real64, dimension(:), pointer :: wsave 
      sll_real64, dimension(:), pointer :: work  
   end type sll_t_poisson_1d_periodic
 
   type, extends(sll_c_poisson_1d_base) :: sll_c_poisson_1d_periodic
 
      type(sll_t_poisson_1d_periodic) :: poiss
 
   contains
      procedure, pass(poisson) :: init => &
         initialize_poisson_1d_periodic_wrapper
      procedure, pass(poisson) :: compute_phi_from_rho => &
         compute_phi_from_rho_1d_periodic
      procedure, pass(poisson) :: compute_e_from_rho => &
         compute_e_from_rho_1d_periodic
 
   end type sll_c_poisson_1d_periodic
 
   interface sll_o_new
      module procedure new_poisson_1d_periodic
   end interface
 
   interface sll_o_initialize
      module procedure initialize_poisson_1d_periodic
   end interface
 
   interface sll_o_solve
      module procedure solve_poisson_1d_periodic
   end interface
 
contains
 
   function new_poisson_1d_periodic(eta1_min, eta1_max, nc_eta1, error) &
      result(self)
      type(sll_t_poisson_1d_periodic), pointer :: self
      sll_int32, intent(in)              :: nc_eta1  
      sll_int32, intent(out)            :: error    
      sll_real64, intent(in)            :: eta1_min 
      sll_real64, intent(in)            :: eta1_max 
 
      sll_allocate(self, error)
      call initialize_poisson_1d_periodic(self, eta1_min, eta1_max, nc_eta1, error)
 
   end function new_poisson_1d_periodic
 
   subroutine initialize_poisson_1d_periodic(self, eta1_min, eta1_max, nc_eta1, error)
 
      type(sll_t_poisson_1d_periodic), intent(out) :: self
      sll_int32, intent(in)                  :: nc_eta1  
      sll_int32, intent(out)                :: error    
      sll_real64, intent(in)                :: eta1_min 
      sll_real64, intent(in)                :: eta1_max 
 
      error = 0
      ! geometry
      self%nc_eta1 = nc_eta1
      self%eta1_min = eta1_min
      self%eta1_max = eta1_max
 
      sll_allocate(self%wsave(2*self%nc_eta1 + 15), error)
 
      call dffti(self%nc_eta1, self%wsave)
 
      ! Allocate auxiliary arrays for fft in order to keep rhs unchanged
      sll_allocate(self%work(nc_eta1 + 1), error)
 
   end subroutine initialize_poisson_1d_periodic
 
   subroutine solve_poisson_1d_periodic(self, field, rhs)
 
      type(sll_t_poisson_1d_periodic), intent(inout) :: self
      sll_real64, dimension(:), intent(out)   :: field
      sll_real64, dimension(:), intent(in)    :: rhs
      sll_int32                               :: ik
      sll_real64                              :: kx0, kx, k2
 
      ! Check that field and rhs are both associated to the
      ! same mesh with the right number of cells
      ! that has been initialized in new_poisson_1d_periodic
      sll_assert(size(field) == self%nc_eta1 + 1)
      sll_assert(size(rhs) == self%nc_eta1 + 1)
 
      ! copy rhs into auxiliary array for fftpack
      ! in order to keep rhs unchanged
      self%work = rhs
 
      ! Forward FFT
      call dfftf(self%nc_eta1, self%work, self%wsave)
 
      self%work = self%work/real(self%nc_eta1, f64)      ! normalize FFT
 
      kx0 = 2_f64*sll_p_pi/(self%eta1_max - self%eta1_min)
 
      ! La moyenne de Ex est nulle donc les composantes de Fourier
      ! correspondant a k=0 sont nulles
      field(1) = 0.0_f64
 
      ! Calcul des autres composantes de Fourier
      do ik = 1, (self%nc_eta1 - 2)/2
         kx = real(ik, f64)*kx0
         k2 = kx*kx
         field(2*ik) = kx/k2*self%work(2*ik + 1)
         field(2*ik + 1) = -kx/k2*self%work(2*ik)
         self%work(2*ik) = 1.0_f64/k2*self%work(2*ik)
         self%work(2*ik + 1) = 1.0_f64/k2*self%work(2*ik + 1)
      end do
 
      field(self%nc_eta1) = 0.0_f64          ! because Im(rhs/2)=0
 
      ! Backward FFT
 
      call dfftb(self%nc_eta1, field, self%wsave)
 
      ! complete last term by periodicity
      field(self%nc_eta1 + 1) = field(1)
 
   end subroutine solve_poisson_1d_periodic
 
   function sll_f_new_poisson_1d_periodic( &
      eta1_min, &
      eta1_max, &
      nc_eta1) &
      result(poisson)
 
      type(sll_c_poisson_1d_periodic), pointer :: poisson
      sll_real64, intent(in) :: eta1_min
      sll_real64, intent(in) :: eta1_max
      sll_int32, intent(in) :: nc_eta1
      sll_int32 :: ierr
 
      sll_allocate(poisson, ierr)
      call initialize_poisson_1d_periodic_wrapper( &
         poisson, &
         eta1_min, &
         eta1_max, &
         nc_eta1)
   end function sll_f_new_poisson_1d_periodic
 
   subroutine initialize_poisson_1d_periodic_wrapper( &
      poisson, &
      eta1_min, &
      eta1_max, &
      nc_eta1)
      class(sll_c_poisson_1d_periodic) :: poisson
      sll_real64, intent(in) :: eta1_min
      sll_real64, intent(in) :: eta1_max
      sll_int32, intent(in) :: nc_eta1
      sll_int32 :: ierr
 
      call initialize_poisson_1d_periodic(poisson%poiss, eta1_min, eta1_max, nc_eta1, ierr)
 
   end subroutine initialize_poisson_1d_periodic_wrapper
 
   ! solves -\Delta phi = rho in 2d
   subroutine compute_phi_from_rho_1d_periodic(poisson, phi, rho)
      class(sll_c_poisson_1d_periodic)    :: poisson
      sll_real64, dimension(:), intent(in)  :: rho
      sll_real64, dimension(:), intent(out) :: phi
 
      sll_error('compute_phi_from_rho_1d_periodic', 'not implemented yet')
 
   end subroutine compute_phi_from_rho_1d_periodic
 
   ! solves E = -\nabla Phi with -\Delta phi = rho in 2d
   subroutine compute_e_from_rho_1d_periodic(poisson, E, rho)
      class(sll_c_poisson_1d_periodic) :: poisson
      sll_real64, dimension(:), intent(in) :: rho
      sll_real64, dimension(:), intent(out) :: e
 
      call sll_o_solve(poisson%poiss, e, rho)
 
   end subroutine compute_e_from_rho_1d_periodic
 
end module sll_m_poisson_1d_periodic