diff --git a/src/compute_Sapj.F90 b/src/compute_Sapj.F90
index 7d1dfb6ab64666b7cc8f99fa066408790b4cb8bb..07a5a4dbe01b4b6bcc9a12de8deefba22c74acec 100644
--- a/src/compute_Sapj.F90
+++ b/src/compute_Sapj.F90
@@ -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
diff --git a/src/moments_eq_rhs.F90 b/src/moments_eq_rhs.F90
index bdbfc06e201db7fa1a7d7d5a2a030d1a368007bb..13b955ac321f6da385496cc0969b5b42c1b2cb7b 100644
--- a/src/moments_eq_rhs.F90
+++ b/src/moments_eq_rhs.F90
@@ -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
diff --git a/src/poisson.F90 b/src/poisson.F90
index fbc1a3eba03cfbb65f7e092fc23db20d7ec0e6b3..f126d60ba58f8e662d7ca241a8471ff008cb7e4b 100644
--- a/src/poisson.F90
+++ b/src/poisson.F90
@@ -1,7 +1,7 @@
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