sll_m_fdistribu4d_dk.F90

Go to the documentation of this file.

!-----------------------------------------------------------
! SELALIB
!-----------------------------------------------------------
!
! MODULE: sll_m_fdistribu4d_dk
!
!
! DESCRIPTION:
!
!-----------------------------------------------------------
module sll_m_fdistribu4d_dk
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#include "sll_assert.h"
#include "sll_memory.h"
#include "sll_working_precision.h"
 
   use sll_m_boundary_condition_descriptors, only: &
      sll_p_hermite, &
      sll_p_periodic
 
   use sll_m_common_coordinate_transformations, only: &
      sll_f_polar_eta1, &
      sll_f_polar_eta2
 
   use sll_m_constants, only: &
      sll_p_pi
 
   use sll_m_cubic_splines, only: &
      sll_s_cubic_spline_1d_compute_interpolant, &
      sll_s_cubic_spline_2d_compute_interpolant, &
      sll_f_cubic_spline_1d_eval, &
      sll_f_cubic_spline_2d_eval, &
      sll_s_cubic_spline_1d_init, &
      sll_s_cubic_spline_2d_init, &
      sll_t_cubic_spline_1d, &
      sll_t_cubic_spline_2d, &
      sll_s_cubic_spline_1d_free, &
      sll_s_cubic_spline_2d_free
 
   implicit none
 
   public :: &
      sll_s_function_xy_from_rtheta, &
      sll_s_init_brtheta, &
      sll_f_init_exact_profile_r, &
      sll_s_init_fequilibrium, &
      sll_s_init_fequilibrium_xy, &
      sll_f_profil_xy_exacte
 
   private
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
   sll_real64, dimension(1) :: whatever  ! dummy params array
 
contains
 
   !----------------------------------------
   ! sech = cosh^-1 definition
   !----------------------------------------
   function sech(x)
      sll_real64, intent(in) :: x
      sll_real64             :: sech
 
      sech = 1._f64/cosh(x)
   end function sech
 
   !---------------------------------------- -----------------------
   ! Initialization of the radial density profiles
   !---------------------------------------- -----------------------
   subroutine init_n0_r(r_peak, inv_Ln, deltarn, n0_rmin, &
                        r_grid, n0_1d)
      sll_real64, intent(in) :: r_peak
      sll_real64, intent(in) :: inv_ln
      sll_real64, intent(in) :: deltarn
      sll_real64, intent(in) :: n0_rmin
      sll_real64, dimension(:), intent(in)    :: r_grid
      sll_real64, dimension(:), intent(inout) :: n0_1d
 
      sll_int32  :: ir, nr
      sll_real64 :: dr, lr
      sll_real64 :: etai, etai_mid
      sll_real64 :: tmp
      sll_real64 :: n0norm_tmp
      sll_real64 :: deltarn_norm
 
      nr = size(r_grid, 1)
      lr = r_grid(nr) - r_grid(1)
      dr = r_grid(2) - r_grid(1)
 
      deltarn_norm = deltarn*lr
 
      !*** compute ns0 solution of :                           ***
      !***  2/(n0(r)+n0(r-1))*(n0(r)-n0(r-1))/dr               ***
      !***                  = -(1/Ln0)*cosh^-2(r-rpeak/deltar) ***
      !***  where (1/Ln) = kappa_n0/R                          ***
      n0_1d(1) = n0_rmin
      do ir = 2, nr
         etai = r_grid(ir - 1)
         etai_mid = etai + dr*0.5_f64
         tmp = -inv_ln* &
               sech((etai_mid - r_peak)/deltarn_norm)**2
         tmp = 0.5_f64*dr*tmp
         n0_1d(ir) = (1._f64 + tmp)/(1._f64 - tmp)*n0_1d(ir - 1)
      end do
 
      !*** normalisation of the density at int(n0(r)rdr)/int(rdr) ***
      ! -> computation of int(n0(r)rdr)
      n0norm_tmp = 0._f64
      do ir = 2, nr - 1
         n0norm_tmp = n0norm_tmp + n0_1d(ir)*r_grid(ir)
      end do
      n0norm_tmp = n0norm_tmp + 0.5_f64* &
                   (n0_1d(1)*r_grid(1) + n0_1d(nr)*r_grid(nr))
      ! -> division by int(rdr)
      n0norm_tmp = n0norm_tmp*2._f64*dr/ &
                   (r_grid(nr)**2 - r_grid(1)**2)
 
      n0_1d(1:nr) = n0_1d(1:nr)/n0norm_tmp
   end subroutine init_n0_r
 
   !---------------------------------------- -----------------------
   ! Initialization of the profiles
   ! solution of the following equation
   !
   !       1   d T(r)    = - kappa cosh^(-2)( (r- rx)/delta )
   !      T(r) dr
   !
   !  kappa = params(1)
   !  delta = params(2)
   !  rx    = params(3)
   !---------------------------------------- -----------------------
   function sll_f_profil_xy_exacte(x, y, params_profil)
 
      sll_real64, intent(in) :: x
      sll_real64, intent(in) :: y
      sll_real64, dimension(:), intent(in) ::params_profil
      sll_real64 :: sll_f_profil_xy_exacte
      sll_real64 :: kappa
      sll_real64 :: delta
      sll_real64 :: rx
      sll_real64 :: r
 
      kappa = params_profil(1)
      delta = params_profil(2)
      rx = params_profil(3)
      r = sqrt(x**2 + y**2)
      sll_f_profil_xy_exacte = exp(-kappa*delta*tanh((r - rx)/delta))
 
   end function sll_f_profil_xy_exacte
 
   !---------------------------------------- -----------------------
   ! Initialization of the profiles
   ! solution of the following equation
   !
   !       1   d T(r)    = - kappa cosh^(-2)( (r- rx)/delta )
   !      T(r) dr
   !
   !  kappa = params(1)
   !  delta = params(2)
   !  rx    = params(3)
   !---------------------------------------- -----------------------
   function sll_f_init_exact_profile_r(r, params_profil)
 
      sll_real64, intent(in) :: r
      sll_real64, dimension(:), intent(in) :: params_profil
 
      sll_real64 :: sll_f_init_exact_profile_r
      sll_real64 :: kappa
      sll_real64 :: delta
      sll_real64 :: rx
 
      kappa = params_profil(1)
      delta = params_profil(2)
      rx = params_profil(3)
      sll_f_init_exact_profile_r = exp(-kappa*delta*tanh((r - rx)/delta))
 
   end function sll_f_init_exact_profile_r
 
   !---------------------------------------------------------------
   ! Initialization of the radial temperature profiles
   !---------------------------------------------------------------
   subroutine init_t_r(r_peak, inv_LT, &
                       deltarT, T_rmin, T_scal, r_grid, T_1d)
      sll_real64, intent(in) :: r_peak
      sll_real64, intent(in) :: inv_lt
      sll_real64, intent(in) :: deltart
      sll_real64, intent(in) :: t_rmin
      sll_real64, intent(in) :: t_scal
      sll_real64, dimension(:), intent(in)    :: r_grid
      sll_real64, dimension(:), intent(inout) :: t_1d
 
      sll_int32  ::  ir, nr
      sll_real64 :: dr, lr
      sll_real64 :: etai, etai_mid
      sll_real64 :: tmp
      sll_real64 :: w0, w1, tnorm_tmp
      sll_real64 :: deltart_norm
 
      nr = size(r_grid, 1)
      lr = r_grid(nr) - r_grid(1)
      dr = r_grid(2) - r_grid(1)
      deltart_norm = deltart*lr
 
      !*** compute Ts solution of :                           ***
      !***  2/(Ts(r)+Ts(r-1))*(Ts(r)-Ts(r-1))/dr               ***
      !***                  = -(1/LTs)*cosh^-2(r-rpeak/deltar) ***
      !***  where (1/LT) = kappa_Ts/R                          ***
      t_1d(1) = t_rmin
      do ir = 2, nr
         etai = r_grid(ir - 1)
         etai_mid = etai + dr*0.5_f64
         tmp = -inv_lt* &
               sech((etai_mid - r_peak)/deltart_norm)**2
         tmp = 0.5_f64*dr*tmp
         t_1d(ir) = (1._f64 + tmp)/(1._f64 - tmp)*t_1d(ir - 1)
      end do
 
      !*** normalisation of the temperature to 1 at r=rpeak ***
      ir = int((r_peak - r_grid(1))/dr)
      w1 = (r_peak - r_grid(ir))/dr
      w0 = 1._f64 - w1
      tnorm_tmp = w0*t_1d(ir) + w1*t_1d(ir + 1)
      t_1d(1:nr) = (t_1d(1:nr)/tnorm_tmp)/t_scal
   end subroutine init_t_r
 
   !---------------------------------------------------------------
   ! Initialisation of the magnetic field B(r,theta)
   !---------------------------------------------------------------
   subroutine sll_s_init_brtheta(r_grid, theta_grid, B_rtheta)
      sll_real64, dimension(:), intent(in)    :: r_grid
      sll_real64, dimension(:), intent(in)    :: theta_grid
      sll_real64, dimension(:, :), intent(inout) :: b_rtheta
 
      sll_int32 :: ir, itheta
      sll_int32 :: nr, ntheta
 
      nr = size(r_grid, 1)
      ntheta = size(theta_grid, 1)
 
      do itheta = 1, ntheta
         do ir = 1, nr
            b_rtheta(ir, itheta) = 1._f64
         end do
      end do
   end subroutine sll_s_init_brtheta
 
   !---------------------------------------------------------------
   ! Computation of the equilibrium distribution function for
   !  drift-kinetic 4D simulation
   !   feq(r,vpar) = n0(r)/(2*pi*Ti(r))**(1/2) *
   !                    exp(-0.5*vpar**2/Ti(r))
   !  where n0 and Ti are respectively the initial density
   !  and temperature profiles
   !---------------------------------------------------------------
   function compute_feq_val(r, vpar, n0_r, Ti_r) &
      result(val)
 
      sll_real64             :: val      ! sll_DK_initializer_4d
      sll_real64, intent(in) :: r
      sll_real64, intent(in) :: vpar
      sll_real64, intent(in) :: n0_r
      sll_real64, intent(in) :: ti_r
 
      sll_assert(r >= 0.0)
      val = n0_r/sqrt(2._f64*sll_p_pi*ti_r)* &
            exp(-0.5_f64*vpar**2/ti_r)
   end function compute_feq_val
 
   !----------------------------------------------------
   ! Initialization of the 2D array for the equilibrium
   !  distribution function feq(r,vpar)
   !  feq(r,vpar) = n0(r)/(2*pi*Ti(r))**(1/2) *
   !                    exp(-0.5*vpar**2/Ti(r))
   !----------------------------------------------------
   subroutine sll_s_init_fequilibrium(Nr, Nvpar, &
                                      r_grid, vpar_grid, n0_1d, Ti_1d, feq_2d)
      sll_int32, intent(in) :: nr
      sll_int32, intent(in) :: nvpar
      sll_real64, dimension(:), intent(in)  :: r_grid
      sll_real64, dimension(:), intent(in)  :: vpar_grid
      sll_real64, dimension(:), intent(in)  :: n0_1d
      sll_real64, dimension(:), intent(in)  :: ti_1d
      sll_real64, dimension(:, :), intent(out) :: feq_2d
 
      sll_int32  :: ir, ivpar
      sll_real64 :: r, vpar, n0_r, ti_r
 
      do ivpar = 1, nvpar
         vpar = vpar_grid(ivpar)
         do ir = 1, nr
            r = r_grid(ir)
            n0_r = n0_1d(ir)
            ti_r = ti_1d(ir)
            feq_2d(ir, ivpar) = compute_feq_val(r, vpar, n0_r, ti_r)
         end do
      end do
   end subroutine sll_s_init_fequilibrium
 
   !----------------------------------------------------
   ! Compute func(x,y) from func(r)
   !----------------------------------------------------
   subroutine function_xy_from_r(r_grid, func_r, &
                                 xgrid_2d, ygrid_2d, func_xy)
      sll_real64, dimension(:), intent(in)  :: r_grid
      sll_real64, dimension(:), intent(in)  :: func_r
      sll_real64, dimension(:, :), intent(in)  :: xgrid_2d
      sll_real64, dimension(:, :), intent(in)  :: ygrid_2d
      sll_real64, dimension(:, :), intent(out) :: func_xy
 
      sll_int32  :: nr, npt1, npt2
      sll_int32  :: ix, iy
      sll_real64 :: r, x, y
 
      type(sll_t_cubic_spline_1d) :: sp1d_r
 
      nr = size(r_grid, 1)
      npt1 = size(xgrid_2d, 1)
      npt2 = size(xgrid_2d, 2)
 
      call sll_s_cubic_spline_1d_init(sp1d_r, npt1, &
                                      r_grid(1), r_grid(nr), sll_p_hermite)
      call sll_s_cubic_spline_1d_compute_interpolant(func_r, sp1d_r)
 
      do iy = 1, npt2
         do ix = 1, npt1
            x = xgrid_2d(ix, iy)
            y = ygrid_2d(ix, iy)
            r = sll_f_polar_eta1(x, y, (/0.0_f64/)) ! params doesn't matter for sll_f_polar_eta1
            r = min(max(r, r_grid(1)), r_grid(nr))
            func_xy(ix, iy) = sll_f_cubic_spline_1d_eval(r, sp1d_r)
         end do
      end do
      call sll_s_cubic_spline_1d_free(sp1d_r)
   end subroutine function_xy_from_r
 
   !----------------------------------------------------
   ! Compute func(x,y) from func(r,theta)
   !----------------------------------------------------
   subroutine sll_s_function_xy_from_rtheta(r_grid, theta_grid, &
                                            func_rtheta, xgrid_2d, ygrid_2d, func_xy)
      sll_real64, dimension(:), intent(in)  :: r_grid
      sll_real64, dimension(:), intent(in)  :: theta_grid
      sll_real64, dimension(:, :), intent(in)  :: func_rtheta
      sll_real64, dimension(:, :), intent(in)  :: xgrid_2d
      sll_real64, dimension(:, :), intent(in)  :: ygrid_2d
      sll_real64, dimension(:, :), intent(out) :: func_xy
 
      sll_int32  :: nr, ntheta
      sll_int32  :: npt1, npt2
      sll_int32  :: ix, iy
      sll_real64 :: r, theta, x, y
 
      type(sll_t_cubic_spline_2d) :: sp2d_rtheta
 
      nr = size(r_grid, 1)
      ntheta = size(theta_grid, 1)
      npt1 = size(xgrid_2d, 1)
      npt2 = size(xgrid_2d, 2)
 
      call sll_s_cubic_spline_2d_init(sp2d_rtheta, &
                                      npt1, npt2, &
                                      r_grid(1), r_grid(nr), &
                                      theta_grid(1), theta_grid(ntheta), &
                                      sll_p_hermite, sll_p_periodic)
      call sll_s_cubic_spline_2d_compute_interpolant(func_rtheta, sp2d_rtheta)
 
      do iy = 1, npt2
         do ix = 1, npt1
            x = xgrid_2d(ix, iy)
            y = ygrid_2d(ix, iy)
            r = sll_f_polar_eta1(x, y, whatever)
            r = min(max(r, r_grid(1)), r_grid(nr))
            theta = sll_f_polar_eta2(x, y, whatever)
            theta = modulo(theta, 2._f64*sll_p_pi)
            func_xy(ix, iy) = sll_f_cubic_spline_2d_eval(r, theta, sp2d_rtheta)
         end do
      end do
      call sll_s_cubic_spline_2d_free(sp2d_rtheta)
   end subroutine sll_s_function_xy_from_rtheta
 
   !----------------------------------------------------
   ! Initialization of the 3D array for the equilibrium
   !  distribution function feq(x,y,vpar)
   !  feq(x,y,vpar) = n0(x,y)/(2*pi*Ti(x,y))**(1/2) *
   !                    exp(-0.5*vpar**2/Ti(x,y))
   !----------------------------------------------------
   subroutine sll_s_init_fequilibrium_xy(xgrid_2d, ygrid_2d, &
                                         vpar_grid, n0_xy, Ti_xy, feq_xyvpar)
      sll_real64, dimension(:, :), intent(in)  :: xgrid_2d
      sll_real64, dimension(:, :), intent(in)  :: ygrid_2d
      sll_real64, dimension(:), intent(in)  :: vpar_grid
      sll_real64, dimension(:, :), intent(in)  :: n0_xy
      sll_real64, dimension(:, :), intent(in)  :: ti_xy
      sll_real64, dimension(:, :, :), intent(out) :: feq_xyvpar
 
      sll_int32 :: npt1, npt2, nvpar
      sll_int32 :: ix, iy, ivpar
 
      npt1 = size(xgrid_2d, 1)
      npt2 = size(ygrid_2d, 2)
      nvpar = size(vpar_grid, 1)
 
      do ivpar = 1, nvpar
         do iy = 1, npt2
            do ix = 1, npt1
               feq_xyvpar(ix, iy, ivpar) = n0_xy(ix, iy)* &
                                           exp(-0.5_f64*vpar_grid(ivpar)**2/ti_xy(ix, iy)) &
                                           /sqrt(2._f64*sll_p_pi*ti_xy(ix, iy))
            end do
         end do
      end do
   end subroutine sll_s_init_fequilibrium_xy
 
end module sll_m_fdistribu4d_dk