From c7b37e0ce355c784b6094d860214538d40b8c80d Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Mon, 22 Jun 2020 14:52:10 +0200
Subject: [PATCH] multiple typos correction
---
src/moments_eq_rhs.F90 | 177 +++++++++++++++++++++--------------------
src/poisson.F90 | 84 ++++++++++---------
2 files changed, 136 insertions(+), 125 deletions(-)
diff --git a/src/moments_eq_rhs.F90 b/src/moments_eq_rhs.F90
index 58171853..538f9124 100644
--- a/src/moments_eq_rhs.F90
+++ b/src/moments_eq_rhs.F90
@@ -12,74 +12,79 @@ SUBROUTINE moments_eq_rhs
INTEGER :: ip,ij, ikr,ikz ! loops indices
REAL(dp) :: ip_dp, ij_dp
- COMPLEX(dp) :: kr, kz, kperp2
- COMPLEX(dp) :: taue_qe_etaBi, taui_qi_etaBi
- COMPLEX(dp) :: kernelj, kerneljp1, kerneljm1, b_e2, b_i2 ! Kernel functions and variable
- COMPLEX(dp) :: factj, xb_e2, xb_i2 ! Auxiliary variables
- COMPLEX(dp) :: xNapj, xNapp2j, xNapm2j, xNapjp1, xNapjm1 ! factors depending on the pj loop
- COMPLEX(dp) :: xphij, xphijp1, xphijm1, xColl
+ REAL(dp) :: kr, kz, kperp2
+ REAL(dp) :: taue_qe_etaB, taui_qi_etaB
+ REAL(dp) :: kernelj, kerneljp1, kerneljm1, b_e2, b_i2 ! Kernel functions and variable
+ REAL(dp) :: factj, xb_e2, xb_i2 ! Auxiliary variables
+ REAL(dp) :: xNapj, xNapp2j, xNapm2j, xNapjp1, xNapjm1 ! factors depending on the pj loop
+ REAL(dp) :: xphij, xphijp1, xphijm1, xColl
COMPLEX(dp) :: TNapj, TNapp2j, TNapm2j, TNapjp1, TNapjm1, Tphi, TColl ! terms of the rhs
- !!!!!!!!! Electrons moments RHS !!!!!!!!!
- Tphi = 0 ! electrostatic potential term
- taue_qe_etaBi = tau_e/q_e * eta_B * imagu ! one-time computable factor
- xb_e2 = sigma_e**2 * tau_e/2. ! species dependant factor of the Kernel, squared
+ !Precompute species dependant factors
+ taue_qe_etaB = tau_e/q_e * eta_B
+ xb_e2 = sigma_e**2 * tau_e!/2.0 ! species dependant factor of the Kernel, squared
+ taui_qi_etaB = tau_i/q_i * eta_B
+ xb_i2 = sigma_i**2 * tau_i!/2.0 ! species dependant factor of the Kernel, squared
- ploope : DO ip = ips_e, ipe_e
- ip_dp = real(ip-1,dp) ! real index is one minus the loop index (0 to pmax)
+ !!!!!!!!! Electrons moments RHS !!!!!!!!!
+ Tphi = 0 ! electrostatic potential term
+ ploope : DO ip = ips_e, ipe_e ! This loop is from 1 to pmaxe+1
+ ip_dp = REAL(ip-1.,dp) ! REAL index is one minus the loop index (0 to pmax)
- ! x N_e^{p+2,j}
- IF (ip+2 .LE. pmaxe + 1) THEN
- xNapp2j = taue_qe_etaBi * SQRT(ip_dp + 1._dp) * SQRT(ip_dp + 2._dp)
+ ! N_e^{p+2,j} multiplicator
+ IF (ip+2 .LE. ipe_e) THEN
+ xNapp2j = -taue_qe_etaB * SQRT((ip_dp + 1.) * (ip_dp + 2.))
ELSE
xNapp2j = 0.
ENDIF
- ! x N_e^{p-2,j}
- IF (ip-2 .GE. 1) THEN
- xNapm2j = taue_qe_etaBi * SQRT(ip_dp) * SQRT(ip_dp - 1.)
+ ! N_e^{p-2,j} multiplicator
+ IF (ip-2 .GE. ips_e) THEN
+ xNapm2j = -taue_qe_etaB * SQRT(ip_dp*(ip_dp - 1.))
ELSE
xNapm2j = 0.
ENDIF
- jloope : DO ij = ijs_e, ije_e
- ij_dp = real(ij-1,dp) ! real index is one minus the loop index (0 to jmax)
+ jloope : DO ij = ijs_e, ije_e ! This loop is from 1 to jmaxe+1
+ ij_dp = REAL(ij-1.,dp) ! REAL index is one minus the loop index (0 to jmax)
- ! x N_e^{p,j+1}
- IF (ij+1 .LE. jmaxe+1) THEN
- xNapjp1 = -taue_qe_etaBi * (ij_dp + 1.)
+ ! N_e^{p,j+1} multiplicator
+ IF (ij+1 .LE. ije_e) THEN
+ xNapjp1 = taue_qe_etaB * (ij_dp + 1.)
ELSE
xNapjp1 = 0.
ENDIF
- ! x N_e^{p,j-1}
- IF (ij-1 .GE. 1) THEN
- xNapjm1 = -taue_qe_etaBi * ij_dp
+ ! N_e^{p,j-1} multiplicator
+ IF (ij-1 .GE. ijs_e) THEN
+ xNapjm1 = taue_qe_etaB * ij_dp
ELSE
xNapjm1 = 0.
ENDIF
- ! x N_e^{pj}
- xNapj = taue_qe_etaBi * 2.*(ip_dp + ij_dp + 1.)
+ ! N_e^{pj} multiplicator
+ xNapj = -taue_qe_etaB * 2.*(ip_dp + ij_dp + 1.)
! first part of collision operator (Lenard-Bernstein)
- xColl = -nu*(ip_dp + 2.*ij_dp)
-
- ! x phi
- IF (ip .eq. 1) THEN !(krokecker delta_p^0)
- xphij = imagu * (eta_n + 2.*ij_dp*eta_T + 2.*eta_B*(ij_dp+1.) )
- xphijp1 = -imagu * (eta_B + eta_T)*(ij_dp+1.)
- xphijm1 = -imagu * (eta_B + eta_T)* ij_dp
- factj = real(Factorial(ij-1),dp)
- ELSE IF (ip .eq. 3) THEN !(krokecker delta_p^2)
- xphij = imagu * (SQRT2*eta_B + eta_T/SQRT2)
- factj = real(Factorial(ij-1),dp)
+ xColl = -(ip_dp + 2.*ij_dp)
+
+ ! phi multiplicator for different kernel numbers
+ IF (ip .EQ. 1) THEN !(kronecker delta_p^0)
+ xphij = (eta_n + 2.*ij_dp*eta_T - 2.*eta_B*(ij_dp+1.) )
+ xphijp1 = -(eta_T - eta_B)*(ij_dp+1.)
+ xphijm1 = -(eta_T - eta_B)* ij_dp
+ factj = REAL(Factorial(ij-1),dp)
+ ELSE IF (ip .EQ. 3) THEN !(kronecker delta_p^2)
+ xphij = (eta_T/SQRT2 - SQRT2*eta_B)
+ factj = REAL(Factorial(ij-1),dp)
ELSE
xphij = 0.; xphijp1 = 0.; xphijm1 = 0.
factj = 1
ENDIF
- ! Loop on kspace
+ !write(*,*) '(ip,ij) = (', ip,',', ij,')'
+
+ ! Loop on kspace
krloope : DO ikr = ikrs,ikre
kzloope : DO ikz = ikzs,ikze
kr = krarray(ikr) ! Poloidal wavevector
@@ -91,33 +96,33 @@ SUBROUTINE moments_eq_rhs
! term propto N_e^{p,j}
TNapj = moments_e(ip,ij,ikr,ikz,updatetlevel) * xNapj
! term propto N_e^{p+2,j}
- IF (xNapp2j .NE. (0.,0.)) THEN
+ IF (ip+2 .LE. ipe_e) THEN
TNapp2j = moments_e(ip+2,ij,ikr,ikz,updatetlevel) * xNapp2j
ELSE
TNapp2j = 0.
ENDIF
! term propto N_e^{p-2,j}
- IF (xNapm2j .NE. (0.,0.)) THEN
+ IF (ip-2 .GE. ips_e) THEN
TNapm2j = moments_e(ip-2,ij,ikr,ikz,updatetlevel) * xNapm2j
ELSE
TNapm2j = 0.
ENDIF
! xterm propto N_e^{p,j+1}
- IF (xNapjp1 .NE. (0.,0.)) THEN
- TNapjp1 = moments_e(ip,ij+1,ikr,ikz,updatetlevel) * xNapp2j
+ IF (ij+1 .LE. ije_e) THEN
+ TNapjp1 = moments_e(ip,ij+1,ikr,ikz,updatetlevel) * xNapjp1
ELSE
TNapjp1 = 0.
ENDIF
! term propto N_e^{p,j-1}
- IF (xNapjm1 .NE. (0.,0.)) THEN
- TNapjm1 = moments_e(ip,ij-1,ikr,ikz,updatetlevel) * xNapp2j
+ IF (ij-1 .GE. ijs_e) THEN
+ TNapjm1 = moments_e(ip,ij-1,ikr,ikz,updatetlevel) * xNapjm1
ELSE
TNapjm1 = 0.
ENDIF
! Collision term completed (Lenard-Bernstein)
- IF (xNapp2j .NE. (0.,0.)) THEN
- TColl = (xColl - nu * b_e2/4.) * moments_e(ip+2,ij,ikr,ikz,updatetlevel)
+ IF (ip+2 .LE. ipe_e) THEN
+ TColl = nu*(xColl - b_e2) * moments_e(ip+2,ij,ikr,ikz,updatetlevel)
ELSE
TColl = 0.
ENDIF
@@ -128,13 +133,12 @@ SUBROUTINE moments_eq_rhs
kernelj = b_e2**(ij-1) * exp(-b_e2)/factj
kerneljp1 = kernelj * b_e2 /(ij_dp + 1.)
kerneljm1 = kernelj * ij_dp / b_e2
- Tphi = (xphij * Kernelj + xphijp1 * Kerneljp1 + xphijm1 * Kerneljm1)* kz * phi(ikr,ikz)
+ Tphi = (xphij*Kernelj + xphijp1*Kerneljp1 + xphijm1*Kerneljm1) * phi(ikr,ikz)
ENDIF
! Sum of all precomputed terms
moments_rhs_e(ip,ij,ikr,ikz,updatetlevel) = &
- kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1) &
- + Tphi + TColl
+ imagu * kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 + Tphi) + TColl
END DO kzloope
END DO krloope
@@ -142,65 +146,64 @@ SUBROUTINE moments_eq_rhs
END DO jloope
END DO ploope
-
!!!!!!!!! Ions moments RHS !!!!!!!!!
- Tphi = 0 ! electrostatic potential term
- taui_qi_etaBi = tau_i/real(q_i)*eta_B * imagu ! one-time computable factor
- xb_i2 = sigma_i**2 * tau_i/2.0 ! species dependant factor of the Kernel, squared
+ Tphi = 0 ! electrostatic potential term
ploopi : DO ip = ips_i, ipe_i
- ip_dp = real(ip-1,dp)
+ ip_dp = REAL(ip-1.,dp)
! x N_i^{p+2,j}
- IF (ip+2 .LE. pmaxi + 1) THEN
- xNapp2j = taui_qi_etaBi * SQRT(ip_dp + 1.) * SQRT(ip_dp + 2.)
+ IF (ip+2 .LE. ipe_i) THEN
+ xNapp2j = -taui_qi_etaB * SQRT((ip_dp + 1.) * (ip_dp + 2.))
ELSE
xNapp2j = 0.
ENDIF
! x N_i^{p-2,j}
- IF (ip-2 .GE. 1) THEN
- xNapm2j = taui_qi_etaBi * SQRT(ip_dp) * SQRT(ip_dp - 1.)
+ IF (ip-2 .GE. ips_i) THEN
+ xNapm2j = -taui_qi_etaB * SQRT(ip_dp * (ip_dp - 1.))
ELSE
xNapm2j = 0.
ENDIF
jloopi : DO ij = ijs_i, ije_i
- ij_dp = real(ij-1,dp)
+ ij_dp = REAL(ij-1.,dp)
! x N_i^{p,j+1}
- IF (ij+1 .LE. jmaxi+1) THEN
- xNapjp1 = -taui_qi_etaBi * (ij_dp + 1.)
+ IF (ij+1 .LE. ije_i) THEN
+ xNapjp1 = taui_qi_etaB * (ij_dp + 1.)
ELSE
xNapjp1 = 0.
ENDIF
! x N_i^{p,j-1}
- IF (ij-1 .GE. 1) THEN
- xNapjm1 = -taui_qi_etaBi * ij_dp
+ IF (ij-1 .GE. ijs_i) THEN
+ xNapjm1 = taui_qi_etaB * ij_dp
ELSE
xNapjm1 = 0.
ENDIF
! x N_i^{pj}
- xNapj = taui_qi_etaBi * 2.*(ip_dp + ij_dp + 1.)
+ xNapj = -taui_qi_etaB * 2.*(ip_dp + ij_dp + 1.)
! first part of collision operator (Lenard-Bernstein)
- xColl = -nu*(ip_dp + 2.0*ij_dp)
+ xColl = -(ip_dp + 2.0*ij_dp)
! x phi
IF (ip .EQ. 1) THEN !(krokecker delta_p^0)
- xphij = imagu * (eta_n + 2.*ij_dp*eta_T + 2.*eta_B*(ij_dp+1.) )
- xphijp1 = -imagu * (eta_B + eta_T)*(ij_dp+1.)
- xphijm1 = -imagu * (eta_B + eta_T)* ij_dp
- factj = REAL(Factorial(ij-1),dp)
+ xphij = (eta_n + 2.*ij_dp*eta_T - 2.*eta_B*(ij_dp+1.) )
+ xphijp1 = -(eta_T - eta_B)*(ij_dp+1.)
+ xphijm1 = -(eta_T - eta_B)* ij_dp
+ factj = REAL(Factorial(ij-1),dp)
ELSE IF (ip .EQ. 3) THEN !(krokecker delta_p^2)
- xphij = imagu * (SQRT2*eta_B + eta_T/SQRT2)
+ xphij = (eta_T/SQRT2 - SQRT2*eta_B)
factj = REAL(Factorial(ij-1),dp)
ELSE
xphij = 0.; xphijp1 = 0.; xphijm1 = 0.
+ factj = 1.
ENDIF
- ! Loop on kspace
+
+ ! Loop on kspace
krloopi : DO ikr = ikrs,ikre
kzloopi : DO ikz = ikzs,ikze
kr = krarray(ikr) ! Poloidal wavevector
@@ -211,31 +214,34 @@ SUBROUTINE moments_eq_rhs
!! Compute moments and mixing terms
! term propto N_i^{p,j}
TNapj = moments_i(ip,ij,ikr,ikz,updatetlevel) * xNapj
-
- IF (xNapp2j .NE. (0.,0.)) THEN ! term propto N_i^{p+2,j}
+ ! term propto N_i^{p+2,j}
+ IF (ip+2 .LE. ipe_i) THEN
TNapp2j = moments_i(ip+2,ij,ikr,ikz,updatetlevel) * xNapp2j
ELSE
TNapp2j = 0.
ENDIF
- IF (xNapm2j .NE. (0.,0.)) THEN ! term propto N_i^{p-2,j}
+ ! term propto N_i^{p-2,j}
+ IF (ip-2 .GE. ips_i) THEN
TNapm2j = moments_i(ip-2,ij,ikr,ikz,updatetlevel) * xNapm2j
ELSE
TNapm2j = 0.
ENDIF
- IF (xNapjp1 .NE. (0.,0.)) THEN ! xterm propto N_a^{p,j+1}
- TNapjp1 = moments_i(ip,ij+1,ikr,ikz,updatetlevel) * xNapp2j
+ ! xterm propto N_i^{p,j+1}
+ IF (ij+1 .LE. ije_i) THEN
+ TNapjp1 = moments_i(ip,ij+1,ikr,ikz,updatetlevel) * xNapjp1
ELSE
TNapjp1 = 0.
ENDIF
- IF (xNapjm1 .NE. (0.,0.)) THEN ! term propto N_a^{p,j-1}
- TNapjm1 = moments_i(ip,ij-1,ikr,ikz,updatetlevel) * xNapp2j
+ ! term propto N_i^{p,j-1}
+ IF (ij-1 .GE. ijs_i) THEN
+ TNapjm1 = moments_i(ip,ij-1,ikr,ikz,updatetlevel) * xNapjm1
ELSE
TNapjm1 = 0.
ENDIF
! Collision term completed (Lenard-Bernstein)
- IF (xNapp2j .NE. (0.,0.)) THEN
- TColl = (xColl - nu * b_i2/4.) * moments_i(ip+2,ij,ikr,ikz,updatetlevel)
+ IF (ip+2 .LE. ipe_i) THEN
+ TColl = nu*(xColl - b_i2) * moments_i(ip+2,ij,ikr,ikz,updatetlevel)
ELSE
TColl = 0.
ENDIF
@@ -244,15 +250,14 @@ SUBROUTINE moments_eq_rhs
Tphi = 0
IF ( (ip .eq. 1) .or. (ip .eq. 3) ) THEN ! 0 otherwise (krokecker delta_p^0, delta_p^2)
kernelj = b_i2**(ij-1) * exp(-b_i2)/factj
- kerneljp1 = kernelj * b_i2 /(ij_dp + 1._dp)
+ kerneljp1 = kernelj * b_i2 /(ij_dp + 1.)
kerneljm1 = kernelj * ij_dp / b_i2
- Tphi = (xphij * Kernelj + xphijp1 * Kerneljp1 + xphijm1 * Kerneljm1)* kz * phi(ikr,ikz)
+ Tphi = (xphij*Kernelj + xphijp1*Kerneljp1 + xphijm1*Kerneljm1) * phi(ikr,ikz)
ENDIF
! Sum of all precomputed terms
moments_rhs_i(ip,ij,ikr,ikz,updatetlevel) = &
- kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1) &
- + Tphi + TColl
+ imagu * kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 + Tphi) + TColl
END DO kzloopi
END DO krloopi
diff --git a/src/poisson.F90 b/src/poisson.F90
index f7ab0c58..1819af6f 100644
--- a/src/poisson.F90
+++ b/src/poisson.F90
@@ -12,16 +12,19 @@ SUBROUTINE poisson
INTEGER :: ikr,ikz, ini,ine, i0j
REAL(dp) :: ini_dp, ine_dp
- COMPLEX(dp) :: kr, kz, kperp2
- COMPLEX(dp) :: xkm_, xk2_ ! sub kernel factor for recursive build
- COMPLEX(dp) :: b_e2, b_i2, alphaD
- COMPLEX(dp) :: sigmaetaue2, sigmaitaui2 ! To avoid redondant computation
+ REAL(dp) :: kr, kz, kperp2
+ REAL(dp) :: Kn, Kn2 ! sub kernel factor for recursive build
+ REAL(dp) :: b_e2, b_i2, alphaD
+ REAL(dp) :: sigmae2_taue_o2, sigmai2_taui_o2, qe2_taue, qi2_taui ! To avoid redondant computation
COMPLEX(dp) :: gammaD, gammaphi
- COMPLEX(dp) :: sum_kernel2_e, sum_kernel2_i
- COMPLEX(dp) :: sum_kernel_mom_e, sum_kernel_mom_i
-
- sigmaetaue2 = sigma_e**2 * tau_e/2.
- sigmaitaui2 = sigma_i**2 * tau_i/2.
+ COMPLEX(dp) :: sum_kernel2_e, sum_kernel2_i ! Store to compute sum Kn^2
+ COMPLEX(dp) :: sum_kernel_mom_e, sum_kernel_mom_i ! Store to compute sum Kn*Napn
+
+ !Precompute species dependant factors
+ sigmae2_taue_o2 = sigma_e**2 * tau_e/2.
+ sigmai2_taui_o2 = sigma_i**2 * tau_i/2.
+ qe2_taue = (q_e**2)/tau_e
+ qi2_taui = (q_i**2)/tau_i
DO ikr=ikrs,ikre
DO ikz=ikzs,ikze
@@ -29,57 +32,60 @@ SUBROUTINE poisson
kz = kzarray(ikz)
kperp2 = kr**2 + kz**2
- ! Compute electrons sum(Kernel * Ne0n) (skm) and sum(Kernel**2) (sk2)
- b_e2 = kperp2 * sigmaetaue2 ! non dim kernel argument (kperp2 sigma_a sqrt(2 tau_a)/2)
+ !! Electrons sum(Kernel * Ne0n) (skm) and sum(Kernel**2) (sk2)
+ b_e2 = kperp2 * sigmae2_taue_o2 ! non dim kernel argument (kperp2 sigma_a sqrt(2 tau_a))
- xkm_ = 1. ! Initialization for n = 0
- xk2_ = 1.
-
+ ! Initialization for n = 0 (ine = 1)
+ Kn = 1.
+ Kn2 = 1.
sum_kernel_mom_e = moments_e(1,1,ikr,ikz,updatetlevel)
- sum_kernel2_e = xk2_
+ sum_kernel2_e = Kn2
+ ! loop over n only if the max polynomial degree is 1 or more
if (jmaxe .GT. 0) then
- DO ine=2,jmaxe+1
+ DO ine=1,jmaxe+1 ! ine = n+1
ine_dp = REAL(ine-1,dp)
- xkm_ = xkm_ * b_e2/ine_dp
- xk2_ = xk2_ *(b_e2/ine_dp)**2
+ Kn = Kn * b_e2/(ine_dp+1) ! We update iteratively the kernel functions (to spare factorial computations)
+ Kn2 = Kn2 *(b_e2/(ine_dp+1))**2 ! the exp term is still missing but does not depend on ine ...
- sum_kernel_mom_e = sum_kernel_mom_e + xkm_ * moments_e(1,ine,ikr,ikz,updatetlevel)
- sum_kernel2_e = sum_kernel2_e + xk2_
+ sum_kernel_mom_e = sum_kernel_mom_e + Kn * moments_e(1,ine,ikr,ikz,updatetlevel)
+ sum_kernel2_e = sum_kernel2_e + Kn2 ! ... sum recursively ...
END DO
endif
- sum_kernel2_e = sum_kernel2_e * exp(-b_e2)
- sum_kernel_mom_e = sum_kernel_mom_e * exp(-2.*b_e2)
+ sum_kernel_mom_e = sum_kernel_mom_e * exp(-b_e2) ! ... multiply the final sum with the missing exponential term
+ sum_kernel2_e = sum_kernel2_e * exp(-2.*b_e2) ! and its squared using exp(x)^2 = exp(2x).
- ! Compute ions sum(Kernel * Ni0n) (skm) and sum(Kernel**2) (sk2)
- b_i2 = kperp2 * sigmaitaui2
-
- xkm_ = 1.
- xk2_ = 1.
+ !! Ions sum(Kernel * Ni0n) (skm) and sum(Kernel**2) (sk2)
+ b_i2 = kperp2 * sigmai2_taui_o2
+ ! Initialization for n = 0 (ini = 1)
+ Kn = 1.
+ Kn2 = 1.
sum_kernel_mom_i = moments_i(1, 1, ikr, ikz, updatetlevel)
- sum_kernel2_i = xk2_
+ sum_kernel2_i = Kn2
+
+ ! loop over n only if the max polynomial degree is 1 or more
if (jmaxi .GT. 0) then
- DO ini=2,jmaxi + 1
+ DO ini=1,jmaxi+1
ini_dp = REAL(ini-1,dp) ! Real index (0 to jmax)
- xkm_ = xkm_ * b_i2/ini_dp
- xk2_ = xk2_ *(b_i2/ini_dp)**2
+ Kn = Kn * b_i2/(ini_dp+1) ! We update iteratively the kernel functions this time for ions ...
+ Kn2 = Kn2 *(b_i2/(ini_dp+1.))**2
- sum_kernel_mom_i = sum_kernel_mom_i + xkm_ * moments_i(1,ini,ikr,ikz,updatetlevel)
- sum_kernel2_i = sum_kernel2_i + xk2_
+ sum_kernel_mom_i = sum_kernel_mom_i + Kn * moments_i(1,ini,ikr,ikz,updatetlevel)
+ sum_kernel2_i = sum_kernel2_i + Kn2 ! ... sum recursively ...
END DO
endif
- sum_kernel2_i = sum_kernel2_i * exp(-b_i2)
- sum_kernel_mom_i = sum_kernel_mom_i * exp(-2.*b_i2)
+ sum_kernel_mom_i = sum_kernel_mom_i * exp(-b_i2) ! ... multiply the final sum with the missing exponential term.
+ sum_kernel2_i = sum_kernel2_i * exp(-2.*b_i2)
- ! Assembling the poisson equation
+ !!! Assembling the poisson equation
alphaD = kperp2 * lambdaD**2.
- gammaD = alphaD + (q_e**2)/tau_e * sum_kernel2_e &
- + (q_i**2)/tau_i * sum_kernel2_i
+ gammaD = alphaD + qe2_taue * (1. - sum_kernel2_e) &
+ + qi2_taui * (1. - sum_kernel2_i)
- gammaphi = q_e * sum_kernel_mom_e + q_i * sum_kernel_mom_i
+ gammaphi = q_e * sum_kernel_mom_e + q_i * sum_kernel_mom_i ! gamma_D * phi term
phi(ikr, ikz) = gammaphi/gammaD
--
GitLab