Skip to content
Snippets Groups Projects
Commit 93124e15 authored by Antoine Cyril David Hoffmann's avatar Antoine Cyril David Hoffmann
Browse files

Use precomputed Kernel and closed truncation at Jmax

parent 50486030
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/compute_Sapj.F90
+ 6
− 14
View file @ 93124e15
......@@ -3,7 +3,7 @@ SUBROUTINE compute_Sapj
!! In real space Sapj ~ b*(grad(phi) x grad(g)) which in moments in fourier becomes
!! Sapj = Sum_n (ikr Kn phi)#(ikz Sum_s d_njs Naps) - (ikz Kn phi)#(ikr Sum_s d_njs Naps)
!! where # denotes the convolution.
USE array, ONLY : dnjs, Sepj, Sipj
USE array, ONLY : dnjs, Sepj, Sipj, kernel_i, kernel_e
USE basic
USE fourier
USE fields!, ONLY : phi, moments_e, moments_i
......@@ -20,7 +20,7 @@ SUBROUTINE compute_Sapj
REAL(dp), DIMENSION(irs:ire,izs:ize) :: fz_real, gr_real, f_times_g
INTEGER :: in, is
REAL(dp):: kr, kz, kerneln, be_2, bi_2, factn
REAL(dp):: kr, kz, kerneln
REAL(dp):: sigmae2_taue_o2, sigmai2_taui_o2
! Execution time start
......@@ -32,7 +32,6 @@ SUBROUTINE compute_Sapj
ploope: DO ip = ips_e,ipe_e ! Loop over Hermite moments
jloope: DO ij = ijs_e, ije_e ! Loop over Laguerre moments
real_data_c = 0._dp ! initialize sum over real nonlinear term
factn = 1
nloope: DO in = 1,jmaxe+1 ! Loop over laguerre for the sum
......@@ -40,8 +39,8 @@ SUBROUTINE compute_Sapj
kzloope: DO ikz = ikzs,ikze ! Loop over kz
kr = krarray(ikr)
kz = kzarray(ikz)
be_2 = sigmae2_taue_o2 * (kr**2 + kz**2)
kerneln = be_2**(in-1)/factn * EXP(-be_2)
kerneln = kernel_e(in, ikr, ikz)
! First convolution terms
Fr_cmpx(ikr,ikz) = imagu*kr* phi(ikr,ikz) * kerneln
Fz_cmpx(ikr,ikz) = imagu*kz* phi(ikr,ikz) * kerneln
......@@ -85,9 +84,6 @@ SUBROUTINE compute_Sapj
real_data_c = real_data_c - real_data_f/Nz/Nr * real_data_g/Nz/Nr
IF ( in + 1 .LE. jmaxe+1 ) THEN
factn = real(in,dp) * factn ! compute (n+1)!
ENDIF
ENDDO nloope
! Put the real nonlinear product into k-space
......@@ -110,7 +106,6 @@ SUBROUTINE compute_Sapj
ploopi: DO ip = ips_e,ipe_e ! Loop over Hermite moments
jloopi: DO ij = ijs_e, ije_e ! Loop over Laguerre moments
real_data_c = 0._dp ! initialize sum over real nonlinear term
factn = 1
nloopi: DO in = 1,jmaxi+1 ! Loop over laguerre for the sum
......@@ -118,8 +113,8 @@ SUBROUTINE compute_Sapj
kzloopi: DO ikz = ikzs,ikze ! Loop over kz
kr = krarray(ikr)
kz = kzarray(ikz)
bi_2 = sigmai2_taui_o2 * (kr**2 + kz**2)
kerneln = bi_2**(in-1)/factn * EXP(-bi_2)
kerneln = kernel_i(in, ikr, ikz)
! First convolution terms
Fr_cmpx(ikr,ikz) = imagu*kr* phi(ikr,ikz) * kerneln
Fz_cmpx(ikr,ikz) = imagu*kz* phi(ikr,ikz) * kerneln
......@@ -163,9 +158,6 @@ SUBROUTINE compute_Sapj
real_data_c = real_data_c - real_data_f/Nz/Nr * real_data_g/Nz/Nr
IF ( in + 1 .LE. jmaxe+1 ) THEN
factn = real(in,dp) * factn ! compute (n+1)!
ENDIF
ENDDO nloopi
! Put the real nonlinear product into k-space
......
src/moments_eq_rhs.F90
+ 29
− 19
View file @ 93124e15
......@@ -15,14 +15,14 @@ SUBROUTINE moments_eq_rhs
REAL(dp) :: kr, kz, kperp2
REAL(dp) :: taue_qe_etaB, taui_qi_etaB
REAL(dp) :: sqrtTaue_qe, sqrtTaui_qi, qe_sigmae_sqrtTaue, qi_sigmai_sqrtTaui
REAL(dp) :: kernelj, kerneljp1, kerneljm1, be_2, bi_2 ! Kernel functions and variable
REAL(dp) :: kernelj, kerneljp1, kerneljm1 ! Kernel functions and variable
REAL(dp) :: factj, sigmae2_taue_o2, sigmai2_taui_o2 ! Auxiliary variables
REAL(dp) :: xNapj, xNapp1j, xNapm1j, xNapp2j, xNapm2j, xNapjp1, xNapjm1 ! Mom. factors depending on the pj loop
REAL(dp) :: xphij, xphijp1, xphijm1, xphijpar ! ESpot. factors depending on the pj loop
REAL(dp) :: xCapj, xCa20, xCa01, xCa10 ! Coll. factors depending on the pj loop
COMPLEX(dp) :: TNapj, TNapp1j, TNapm1j, TNapp2j, TNapm2j, TNapjp1, TNapjm1, Tphi
COMPLEX(dp) :: TColl, TColl20, TColl01, TColl10 ! terms of the rhs
COMPLEX(dp) :: test_nan
COMPLEX(dp) :: test_nan, i_kz
REAL(dp) :: nu_e, nu_i, nu_ee, nu_ie ! Species collisional frequency
! Measuring execution time
......@@ -127,10 +127,11 @@ SUBROUTINE moments_eq_rhs
kzloope : DO ikz = ikzs,ikze
kr = krarray(ikr) ! Poloidal wavevector
kz = kzarray(ikz) ! Toroidal wavevector
IF (Nkz .EQ. 1) kz = krarray(ikr) ! If 1D simulation we put kr as kz
i_kz = imagu * kz ! Ddz derivative
! If 1D simulation we put kr as kz since this dim is halved
IF (Nkz .EQ. 1) i_kz = imagu * krarray(ikr)
kperp2 = kr**2 + kz**2 ! perpendicular wavevector
be_2 = kperp2 * sigmae2_taue_o2 ! Kernel argument
!! Compute moments and mixing terms
! term propto N_e^{p,j}
......@@ -226,14 +227,20 @@ SUBROUTINE moments_eq_rhs
ENDIF
!! Electrical potential term
IF ( (p_int .LE. 2) ) THEN ! kronecker p0 p1 p2
kernelj = be_2**j_int * exp(-be_2)/factj
kerneljp1 = kernelj * be_2 /(j_dp + 1._dp)
IF ( be_2 .NE. 0 ) THEN
kerneljm1 = kernelj * j_dp / be_2
IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2
kernelj = kernel_e(ij, ikr, ikz)
IF ( j_int+1 .LE. jmaxe ) THEN
kerneljp1 = kernel_e(ij+1, ikr, ikz)
ELSE
kerneljp1 = 0._dp
ENDIF
IF ( j_int-1 .GE. 0 ) THEN
kerneljm1 = kernel_e(ij-1, ikr, ikz)
ELSE
kerneljm1 = 0._dp
ENDIF
Tphi = (xphij*kernelj + xphijp1*kerneljp1 + xphijm1*kerneljm1) * phi(ikr,ikz)
ELSE
Tphi = 0._dp
......@@ -241,7 +248,7 @@ SUBROUTINE moments_eq_rhs
! Sum of all linear terms
moments_rhs_e(ip,ij,ikr,ikz,updatetlevel) = &
-imagu * kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 - Tphi)&
-i_kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 - Tphi)&
-imagu * kpar*(TNapp1j + TNapm1j + xphijpar*kernelj*phi(ikr,ikz)) &
- mu*kperp2**2 * moments_e(ip,ij,ikr,ikz,updatetlevel) &
+ TColl
......@@ -340,10 +347,9 @@ SUBROUTINE moments_eq_rhs
kzloopi : DO ikz = ikzs,ikze
kr = krarray(ikr) ! Poloidal wavevector
kz = kzarray(ikz) ! Toroidal wavevector
IF (Nkz .EQ. 1) kz = krarray(ikr) ! If 1D simulation we put kr as kz
i_kz = imagu * kz ! Ddz derivative
IF (Nkz .EQ. 1) i_kz = imagu * krarray(ikr) ! If 1D simulation we put kr as kz
kperp2 = kr**2 + kz**2 ! perpendicular wavevector
bi_2 = kperp2 * sigmai2_taui_o2 ! Kernel argument
!! Compute moments and mixing terms
! term propto N_i^{p,j}
......@@ -439,11 +445,15 @@ SUBROUTINE moments_eq_rhs
ENDIF
!! Electrical potential term
IF ( (p_int .LE. 2) ) THEN ! kronecker p0 p1 p2
kernelj = bi_2**j_int * exp(-bi_2)/factj
kerneljp1 = kernelj * bi_2 /(j_dp + 1._dp)
IF ( bi_2 .NE. 0 ) THEN
kerneljm1 = kernelj * j_dp / bi_2
IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2
kernelj = kernel_i(ij, ikr, ikz)
IF (j_int+1 .LE. jmaxi) THEN
kerneljp1 = kernel_i(ij+1, ikr, ikz)
ELSE
kerneljp1 = 0.0_dp
ENDIF
IF ( j_int-1 .GE. 0 ) THEN
kerneljm1 = kernel_i(ij-1, ikr, ikz)
ELSE
kerneljm1 = 0._dp
ENDIF
......@@ -454,7 +464,7 @@ SUBROUTINE moments_eq_rhs
! Sum of linear terms
moments_rhs_i(ip,ij,ikr,ikz,updatetlevel) = &
-imagu * kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 - Tphi)&
-i_kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 - Tphi)&
-imagu * kpar*(TNapp1j + TNapm1j + xphijpar*kernelj*phi(ikr,ikz)) &
- mu*kperp2**2 * moments_i(ip,ij,ikr,ikz,updatetlevel) &
+ TColl
......
src/poisson.F90
+ 23
− 53
View file @ 93124e15
SUBROUTINE poisson
! Solve poisson equation to get phi
USE basic, ONLY: t0_poisson, t1_poisson, tc_poisson
USE basic
USE time_integration, ONLY: updatetlevel
USE array
USE fields
......@@ -11,12 +11,10 @@ SUBROUTINE poisson
USE prec_const
IMPLICIT NONE
INTEGER :: ini,ine, i0j
REAL(dp) :: ini_dp, ine_dp
REAL(dp) :: kr, kz, kperp2
INTEGER :: ini,ine
REAL(dp) :: Kne, Kni ! sub kernel factor for recursive build
REAL(dp) :: be_2, bi_2, alphaD
REAL(dp) :: sigmae2_taue_o2, sigmai2_taui_o2, qe2_taue, qi2_taui ! To avoid redondant computation
REAL(dp) :: alphaD
REAL(dp) :: qe2_taue, qi2_taui ! To avoid redondant computation
REAL(dp) :: sum_kernel2_e, sum_kernel2_i ! Store sum Kn^2
COMPLEX(dp) :: sum_kernel_mom_e, sum_kernel_mom_i ! Store sum Kn*Napn
REAL(dp) :: gammaD
......@@ -26,67 +24,39 @@ SUBROUTINE poisson
CALL cpu_time(t0_poisson)
!Precompute species dependant factors
sigmae2_taue_o2 = sigma_e**2 * tau_e/2._dp ! factor of the Kernel argument
sigmai2_taui_o2 = sigma_i**2 * tau_i/2._dp ! (b_a/2)^2 = (kperp sqrt(2 tau_a) sigma_a/2)^2
! = kperp^2 tau_a sigma_a^2/2
qe2_taue = (q_e**2)/tau_e ! factor of the gammaD sum
qi2_taui = (q_i**2)/tau_i
DO ikr=ikrs,ikre
DO ikz=ikzs,ikze
kr = krarray(ikr)
kz = kzarray(ikz)
kperp2 = kr**2 + kz**2
! non dim kernel argument (kperp2 sigma_a sqrt(2 tau_a))
be_2 = kperp2 * sigmae2_taue_o2
bi_2 = kperp2 * sigmai2_taui_o2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!! Electrons sum(Kernel * Ne0n) (skm) and sum(Kernel**2) (sk2)
!! Sum(kernel * Moments_0n)
! Initialization for n = 0 (ine = 1)
Kne = exp(-be_2)
sum_kernel_mom_e = Kne*moments_e(1, 1, ikr, ikz, updatetlevel)
sum_kernel2_e = Kne**2
! loop over n only if the max polynomial degree is 1 or more
if (jmaxe .GT. 0) then
DO ine=2,jmaxe+1 ! ine = n+1
ine_dp = REAL(ine-1,dp)
! We update iteratively the kernel functions (to spare factorial computations)
Kne = Kne * be_2/ine_dp
sum_kernel_mom_e = sum_kernel_mom_e + Kne * moments_e(1, ine, ikr, ikz, updatetlevel)
sum_kernel2_e = sum_kernel2_e + Kne**2 ! ... sum recursively ...
END DO
ENDIF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
sum_kernel_mom_e = 0._dp
sum_kernel2_e = 0._dp
! loop over n only if the max polynomial degree
DO ine=1,jmaxe+1 ! ine = n+1
Kne = kernel_e(ine, ikr, ikz)
sum_kernel_mom_e = sum_kernel_mom_e + Kne * moments_e(1, ine, ikr, ikz, updatetlevel)
sum_kernel2_e = sum_kernel2_e + Kne**2 ! ... sum recursively ...
END DO
!!!!!!!!!!!!!!!!! Ions sum(Kernel * Ni0n) (skm) and sum(Kernel**2) (sk2)
!! Sum(kernel * Moments_0n)
! Initialization for n = 0 (ini = 1)
Kni = exp(-bi_2)
sum_kernel_mom_i = Kni*moments_i(1, 1, ikr, ikz, updatetlevel)
sum_kernel2_i = Kni**2
! loop over n only if the max polynomial degree is 1 or more
if (jmaxi .GT. 0) then
DO ini=2,jmaxi+1
ini_dp = REAL(ini-1,dp) ! Real index (0 to jmax)
! We update iteratively to spare factorial computations
Kni = Kni * bi_2/ini_dp
sum_kernel_mom_i = sum_kernel_mom_i + Kni * moments_i(1, ini, ikr, ikz, updatetlevel)
sum_kernel2_i = sum_kernel2_i + Kni**2 ! ... sum recursively ...
END DO
ENDIF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
sum_kernel_mom_i = 0
sum_kernel2_i = 0
! loop over n only if the max polynomial degree
DO ini=1,jmaxi+1
Kni = kernel_i(ini, ikr, ikz)
sum_kernel_mom_i = sum_kernel_mom_i + Kni * moments_i(1, ini, ikr, ikz, updatetlevel)
sum_kernel2_i = sum_kernel2_i + Kni**2 ! ... sum recursively ...
END DO
!!!!!!!!!!!!!!! Assembling the poisson equation !!!!!!!!!!!!!!!!!!!!!!!!!!
alphaD = kperp2 * lambdaD**2
alphaD = (krarray(ikr)**2 + kzarray(ikz)**2) * lambdaD**2
gammaD = alphaD + qe2_taue * (1._dp - sum_kernel2_e) & ! Called Poisson_ in MOLI
+ qi2_taui * (1._dp - sum_kernel2_i)
gammaD_phi = q_e * sum_kernel_mom_e + q_i * sum_kernel_mom_i
IF ( (gammaD .EQ. 0) .AND. (abs(kr)+abs(kz) .NE. 0._dp) ) THEN
WRITE(*,*) 'Warning gammaD = 0', sum_kernel2_i
ENDIF
phi(ikr, ikz) = gammaD_phi/gammaD
END DO
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment