From cec66ce8cfc73301ec499510722ee27743074a75 Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Tue, 3 Aug 2021 14:15:38 +0200
Subject: [PATCH] z-dep kernel and mirror force terms
---
src/array_mod.F90 | 13 +++---
src/memory.F90 | 20 ++++++++--
src/moments_eq_rhs.F90 | 89 +++++++++++++++++++++++++++++-------------
src/numerics_mod.F90 | 83 ++++++++++++++++++++++++---------------
4 files changed, 138 insertions(+), 67 deletions(-)
diff --git a/src/array_mod.F90 b/src/array_mod.F90
index 16a1bdf3..ccdb267a 100644
--- a/src/array_mod.F90
+++ b/src/array_mod.F90
@@ -25,10 +25,13 @@ MODULE array
! lin rhs p,j coefficient storage (ip,ij)
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: xnepj,xnipj
- REAL(dp), DIMENSION(:), ALLOCATABLE :: xnepp1j, xnepm1j, xnepp2j, xnepm2j, xnepjp1, xnepjm1
+ REAL(dp), DIMENSION(:), ALLOCATABLE :: xnepp1j, xnepm1j, xnepp2j, xnepm2j, xnepjp1, xnepjm1
+ REAL(dp), DIMENSION(:,:), ALLOCATABLE :: ynepp1j, ynepm1j, ynepp1jm1, ynepm1jm1 ! mirror lin coeff for non adiab mom
+ REAL(dp), DIMENSION(:,:), ALLOCATABLE :: zNepm1j, zNepm1jp1, zNepm1jm1 ! mirror lin coeff for adiab mom
REAL(dp), DIMENSION(:), ALLOCATABLE :: xnipp1j, xnipm1j, xnipp2j, xnipm2j, xnipjp1, xnipjm1
+ REAL(dp), DIMENSION(:,:), ALLOCATABLE :: ynipp1j, ynipm1j, ynipp1jm1, ynipm1jm1 ! mirror lin coeff for non adiab mom
+ REAL(dp), DIMENSION(:,:), ALLOCATABLE :: zNipm1j, zNipm1jp1, zNipm1jm1 ! mirror lin coeff for adiab mom
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: xphij, xphijp1, xphijm1
-
! Geoemtrical operators
! Curvature
REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: Ckxky ! dimensions: kx, ky, z
@@ -49,9 +52,9 @@ MODULE array
REAL(dp), DIMENSION(:), allocatable :: Gamma3
! Some geometrical coefficients
REAL(dp), DIMENSION(:) , allocatable :: gradz_coeff ! 1 / [ J_{xyz} \hat{B} ]
- ! Kernel function evaluation (ij,ikx,iky)
- REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: kernel_e
- REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: kernel_i
+ ! Kernel function evaluation (ij,ikx,iky,iz)
+ REAL(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: kernel_e
+ REAL(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: kernel_i
!! Diagnostics
! Gyrocenter density for electron and ions (ikx,iky,iz)
diff --git a/src/memory.F90 b/src/memory.F90
index dab37da4..31d098e0 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -25,11 +25,11 @@ SUBROUTINE memory
CALL allocate_array(temp_e, ikxs,ikxe, ikys,ikye, izs,ize)
CALL allocate_array(temp_i, ikxs,ikxe, ikys,ikye, izs,ize)
- !___________________ 3D ARRAYS __________________________
+ !___________________ 4D ARRAYS __________________________
!! FLR kernels functions
! Kernel evaluation from j= -1 to jmax+1 for truncation
- CALL allocate_array(Kernel_e, ijsg_e,ijeg_e, ikxs,ikxe, ikys,ikye)
- CALL allocate_array(Kernel_i, ijsg_i,ijeg_i, ikxs,ikxe, ikys,ikye)
+ CALL allocate_array(Kernel_e, ijsg_e,ijeg_e, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array(Kernel_i, ijsg_i,ijeg_i, ikxs,ikxe, ikys,ikye, izs,ize)
!___________________ 5D ARRAYS __________________________
! Moments with ghost degrees for p+2 p-2, j+1, j-1 truncations
@@ -58,6 +58,13 @@ SUBROUTINE memory
CALL allocate_array( xnepm2j, ips_e,ipe_e)
CALL allocate_array( xnepjp1, ijs_e,ije_e)
CALL allocate_array( xnepjm1, ijs_e,ije_e)
+ CALL allocate_array( ynepp1j, ips_e,ipe_e, ijs_e,ije_e)
+ CALL allocate_array( ynepm1j, ips_e,ipe_e, ijs_e,ije_e)
+ CALL allocate_array( ynepp1jm1, ips_e,ipe_e, ijs_e,ije_e)
+ CALL allocate_array( ynepm1jm1, ips_e,ipe_e, ijs_e,ije_e)
+ CALL allocate_array( zNepm1j, ips_e,ipe_e, ijs_e,ije_e)
+ CALL allocate_array( zNepm1jp1, ips_e,ipe_e, ijs_e,ije_e)
+ CALL allocate_array( zNepm1jm1, ips_e,ipe_e, ijs_e,ije_e)
! ions
CALL allocate_array( xnipj, ips_i,ipe_i, ijs_i,ije_i)
CALL allocate_array( xnipp2j, ips_i,ipe_i)
@@ -66,6 +73,13 @@ SUBROUTINE memory
CALL allocate_array( xnipm2j, ips_i,ipe_i)
CALL allocate_array( xnipjp1, ijs_i,ije_i)
CALL allocate_array( xnipjm1, ijs_i,ije_i)
+ CALL allocate_array( ynipp1j, ips_i,ipe_i, ijs_i,ije_i)
+ CALL allocate_array( ynipm1j, ips_i,ipe_i, ijs_i,ije_i)
+ CALL allocate_array( ynipp1jm1, ips_i,ipe_i, ijs_i,ije_i)
+ CALL allocate_array( ynipm1jm1, ips_i,ipe_i, ijs_i,ije_i)
+ CALL allocate_array( zNipm1j, ips_i,ipe_i, ijs_i,ije_i)
+ CALL allocate_array( zNipm1jp1, ips_i,ipe_i, ijs_i,ije_i)
+ CALL allocate_array( zNipm1jm1, ips_i,ipe_i, ijs_i,ije_i)
! elect. pot.
CALL allocate_array( xphij, ips_i,ipe_i, ijs_i,ije_i)
CALL allocate_array( xphijp1, ips_i,ipe_i, ijs_i,ije_i)
diff --git a/src/moments_eq_rhs.F90 b/src/moments_eq_rhs.F90
index 9d41db73..58ccf7ff 100644
--- a/src/moments_eq_rhs.F90
+++ b/src/moments_eq_rhs.F90
@@ -14,11 +14,14 @@ SUBROUTINE moments_eq_rhs_e
USE model
USE prec_const
USE collision
+ use geometry
IMPLICIT NONE
INTEGER :: p_int, j_int ! loops indices and polynom. degrees
- REAL(dp) :: kx, ky, kperp2
- COMPLEX(dp) :: Tnepj, Tnepp2j, Tnepm2j, Tnepjp1, Tnepjm1, Tpare, Tphi
+ REAL(dp) :: kx, ky, kperp2, dzlnB_o_J
+ COMPLEX(dp) :: Tnepj, Tnepp2j, Tnepm2j, Tnepjp1, Tnepjm1, Tpare, Tphi ! Terms from b x gradB and drives
+ COMPLEX(dp) :: Tmir, Tnepp1j, Tnepm1j, Tnepp1jm1, Tnepm1jm1 ! Terms from mirror force with non adiab moments
+ COMPLEX(dp) :: UNepm1j, UNepm1jp1, UNepm1jm1 ! Terms from mirror force with adiab moments
COMPLEX(dp) :: TColl ! terms of the rhs
COMPLEX(dp) :: i_ky
REAL(dp) :: delta_p0, delta_p1, delta_p2
@@ -49,38 +52,52 @@ SUBROUTINE moments_eq_rhs_e
ky = kyarray(iky) ! toroidal wavevector
i_ky = imagu * ky ! toroidal derivative
IF (Nky .EQ. 1) i_ky = imagu * kxarray(ikx) ! If 1D simulation we put kx as ky
- kperp2 = kx**2 + ky**2 ! perpendicular wavevector
+ kperp2= kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
!! Compute moments mixing terms
! Perpendicular dynamic
! term propto n_e^{p,j}
Tnepj = xnepj(ip,ij) * (moments_e(ip,ij,ikx,iky,iz,updatetlevel) &
- +kernel_e(ij,ikx,iky)*qe_taue*phi(ikx,iky,iz)*delta_p0)
+ +kernel_e(ij,ikx,iky,iz)*qe_taue*phi(ikx,iky,iz)*delta_p0)
! term propto n_e^{p+2,j}
Tnepp2j = xnepp2j(ip) * moments_e(ip+2,ij,ikx,iky,iz,updatetlevel)
! term propto n_e^{p-2,j}
Tnepm2j = xnepm2j(ip) * (moments_e(ip-2,ij,ikx,iky,iz,updatetlevel) &
- +kernel_e(ij,ikx,iky)*qe_taue*phi(ikx,iky,iz)*delta_p2)
+ +kernel_e(ij,ikx,iky,iz)*qe_taue*phi(ikx,iky,iz)*delta_p2)
! term propto n_e^{p,j+1}
Tnepjp1 = xnepjp1(ij) * (moments_e(ip,ij+1,ikx,iky,iz,updatetlevel) &
- +kernel_e(ij+1,ikx,iky)*qe_taue*phi(ikx,iky,iz)*delta_p0)
+ +kernel_e(ij+1,ikx,iky,iz)*qe_taue*phi(ikx,iky,iz)*delta_p0)
! term propto n_e^{p,j-1}
Tnepjm1 = xnepjm1(ij) * (moments_e(ip,ij-1,ikx,iky,iz,updatetlevel) &
- +kernel_e(ij-1,ikx,iky)*qe_taue*phi(ikx,iky,iz)*delta_p0)
+ +kernel_e(ij-1,ikx,iky,iz)*qe_taue*phi(ikx,iky,iz)*delta_p0)
! Parallel dynamic
- ! term propto n_i^{p+1,j}, centered FDF
+ ! term propto centered FDF dz(n_a)
Tpare = xnepp1j(ip) * &
(moments_e(ip+1,ij,ikx,iky,iznext,updatetlevel)&
-moments_e(ip+1,ij,ikx,iky,izprev,updatetlevel)) &
+xnepm1j(ip) * &
- (moments_e(ip-1,ij,ikx,iky,iznext,updatetlevel)+kernel_e(ij,ikx,iky)*qe_taue*phi(ikx,iky,iznext)*delta_p1&
- -moments_e(ip-1,ij,ikx,iky,izprev,updatetlevel)-kernel_e(ij,ikx,iky)*qe_taue*phi(ikx,iky,izprev)*delta_p1)
+ (moments_e(ip-1,ij,ikx,iky,iznext,updatetlevel)+kernel_e(ij,ikx,iky,iz)*qe_taue*phi(ikx,iky,iznext)*delta_p1&
+ -moments_e(ip-1,ij,ikx,iky,izprev,updatetlevel)-kernel_e(ij,ikx,iky,iz)*qe_taue*phi(ikx,iky,izprev)*delta_p1)
+
+ ! Mirror terms
+ Tnepp1j = ynepp1j(ip,ij) * moments_e(ip+1, ij,ikx,iky,iz,updatetlevel)
+ Tnepp1jm1 = ynepp1jm1(ip,ij) * moments_e(ip+1,ij-1,ikx,iky,iz,updatetlevel)
+ Tnepm1j = ynepm1j(ip,ij) * (moments_e(ip-1, ij,ikx,iky,iz,updatetlevel) &
+ +kernel_e(ij,ikx,iky,iz)*qe_taue*phi(ikx,iky,iz)*delta_p1)
+ Tnepm1jm1 = ynepm1jm1(ip,ij) * (moments_e(ip-1,ij-1,ikx,iky,iz,updatetlevel) &
+ +kernel_e(ij-1,ikx,iky,iz)*qe_taue*phi(ikx,iky,iz)*delta_p1)
+
+ UNepm1j = zNepm1j(ip,ij) * moments_e(ip+1, ij,ikx,iky,iz,updatetlevel)
+ UNepm1jp1 = zNepm1jp1(ip,ij) * moments_e(ip-1,ij+1,ikx,iky,iz,updatetlevel)
+ UNepm1jm1 = zNepm1jm1(ip,ij) * moments_e(ip-1,ij-1,ikx,iky,iz,updatetlevel)
+
+ Tmir = Tnepp1j + Tnepp1jm1 + Tnepm1j + Tnepm1jm1 + UNepm1j + UNepm1jp1 + UNepm1jm1
!! Electrical potential term
IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2
- Tphi = phi(ikx,iky,iz) * (xphij(ip,ij)*kernel_e(ij, ikx, iky) &
- + xphijp1(ip,ij)*kernel_e(ij+1, ikx, iky) &
- + xphijm1(ip,ij)*kernel_e(ij-1, ikx, iky) )
+ Tphi = phi(ikx,iky,iz) * (xphij(ip,ij)*kernel_e(ij,ikx,iky,iz) &
+ + xphijp1(ip,ij)*kernel_e(ij+1,ikx,iky,iz) &
+ + xphijm1(ip,ij)*kernel_e(ij-1,ikx,iky,iz) )
ELSE
Tphi = 0._dp
ENDIF
@@ -96,8 +113,9 @@ SUBROUTINE moments_eq_rhs_e
!! Sum of all linear terms (the sign is inverted to match RHS)
moments_rhs_e(ip,ij,ikx,iky,iz,updatetlevel) = &
- - Ckxky(ikx,iky,iz) * (Tnepj + Tnepp2j + Tnepm2j + Tnepjp1 + Tnepjm1)&
- - Tpare/2._dp/deltaz/q0 &
+ - imagu*Ckxky(ikx,iky,iz) * (Tnepj + Tnepp2j + Tnepm2j + Tnepjp1 + Tnepjm1)&
+ - Tpare/2._dp/deltaz*gradz_coeff(iz) &
+ - gradzB(iz)* Tmir *gradz_coeff(iz) &
+ i_ky * Tphi &
- mu*kperp2**2 * moments_e(ip,ij,ikx,iky,iz,updatetlevel) &
+ TColl
@@ -141,6 +159,8 @@ SUBROUTINE moments_eq_rhs_i
INTEGER :: p_int, j_int ! loops indices and polynom. degrees
REAL(dp) :: kx, ky, kperp2
COMPLEX(dp) :: Tnipj, Tnipp2j, Tnipm2j, Tnipjp1, Tnipjm1, Tpari, Tphi
+ COMPLEX(dp) :: Tmir, Tnipp1j, Tnipm1j, Tnipp1jm1, Tnipm1jm1 ! Terms from mirror force with non adiab moments
+ COMPLEX(dp) :: UNipm1j, UNipm1jp1, UNipm1jm1 ! Terms from mirror force with adiab moments
COMPLEX(dp) :: TColl ! terms of the rhs
COMPLEX(dp) :: i_ky
REAL(dp) :: delta_p0, delta_p1, delta_p2
@@ -171,38 +191,52 @@ SUBROUTINE moments_eq_rhs_i
ky = kyarray(iky) ! toroidal wavevector
i_ky = imagu * ky ! toroidal derivative
IF (Nky .EQ. 1) i_ky = imagu * kxarray(ikx) ! If 1D simulation we put kx as ky
- kperp2 = kx**2 + ky**2 ! perpendicular wavevector
+ kperp2= kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
!! Compute moments mixing terms
! Perpendicular dynamic
! term propto n_i^{p,j}
Tnipj = xnipj(ip,ij) * (moments_i(ip,ij,ikx,iky,iz,updatetlevel) &
- +kernel_i(ij,ikx,iky)*qi_taui*phi(ikx,iky,iz)*delta_p0)
+ +kernel_i(ij,ikx,iky,iz)*qi_taui*phi(ikx,iky,iz)*delta_p0)
! term propto n_i^{p+2,j}
Tnipp2j = xnipp2j(ip) * moments_i(ip+2,ij,ikx,iky,iz,updatetlevel)
! term propto n_i^{p-2,j}
Tnipm2j = xnipm2j(ip) * (moments_i(ip-2,ij,ikx,iky,iz,updatetlevel) &
- +kernel_i(ij,ikx,iky)*qi_taui*phi(ikx,iky,iz)*delta_p2)
+ +kernel_i(ij,ikx,iky,iz)*qi_taui*phi(ikx,iky,iz)*delta_p2)
! term propto n_e^{p,j+1}
Tnipjp1 = xnipjp1(ij) * (moments_i(ip,ij+1,ikx,iky,iz,updatetlevel) &
- +kernel_i(ij+1,ikx,iky)*qi_taui*phi(ikx,iky,iz)*delta_p0)
+ +kernel_i(ij+1,ikx,iky,iz)*qi_taui*phi(ikx,iky,iz)*delta_p0)
! term propto n_e^{p,j-1}
Tnipjm1 = xnipjm1(ij) * (moments_i(ip,ij-1,ikx,iky,iz,updatetlevel) &
- +kernel_i(ij-1,ikx,iky)*qi_taui*phi(ikx,iky,iz)*delta_p0)
+ +kernel_i(ij-1,ikx,iky,iz)*qi_taui*phi(ikx,iky,iz)*delta_p0)
! Parallel dynamic
! term propto N_i^{p,j+1}, centered FDF
Tpari = xnipp1j(ip) * &
(moments_i(ip+1,ij,ikx,iky,iznext,updatetlevel)&
-moments_i(ip+1,ij,ikx,iky,izprev,updatetlevel)) &
+xnipm1j(ip) * &
- (moments_i(ip-1,ij,ikx,iky,iznext,updatetlevel)+qi_taui*kernel_i(ij,ikx,iky)*phi(ikx,iky,iznext)*delta_p1&
- -moments_i(ip-1,ij,ikx,iky,izprev,updatetlevel)-qi_taui*kernel_i(ij,ikx,iky)*phi(ikx,iky,izprev)*delta_p1)
+ (moments_i(ip-1,ij,ikx,iky,iznext,updatetlevel)+qi_taui*kernel_i(ij,ikx,iky,iz)*phi(ikx,iky,iznext)*delta_p1&
+ -moments_i(ip-1,ij,ikx,iky,izprev,updatetlevel)-qi_taui*kernel_i(ij,ikx,iky,iz)*phi(ikx,iky,izprev)*delta_p1)
+
+ ! Mirror terms
+ Tnipp1j = ynipp1j(ip,ij) * moments_i(ip+1, ij,ikx,iky,iz,updatetlevel)
+ Tnipp1jm1 = ynipp1jm1(ip,ij) * moments_i(ip+1,ij-1,ikx,iky,iz,updatetlevel)
+ Tnipm1j = ynipm1j(ip,ij) * (moments_i(ip-1, ij,ikx,iky,iz,updatetlevel) &
+ +kernel_i(ij,ikx,iky,iz)*qi_taui*phi(ikx,iky,iz)*delta_p1)
+ Tnipm1jm1 = ynipm1jm1(ip,ij) * (moments_i(ip-1,ij-1,ikx,iky,iz,updatetlevel) &
+ +kernel_i(ij-1,ikx,iky,iz)*qi_taui*phi(ikx,iky,iz)*delta_p1)
+
+ Unipm1j = znipm1j(ip,ij) * moments_i(ip+1,ij,ikx,iky,iz,updatetlevel)
+ Unipm1jp1 = znipm1jp1(ip,ij) * moments_i(ip-1,ij+1,ikx,iky,iz,updatetlevel)
+ Unipm1jm1 = znipm1jm1(ip,ij) * moments_i(ip-1,ij-1,ikx,iky,iz,updatetlevel)
+
+ Tmir = Tnipp1j + Tnipp1jm1 + Tnipm1j + Tnipm1jm1 + UNipm1j + UNipm1jp1 + UNipm1jm1
!! Electrical potential term
IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2
- Tphi = phi(ikx,iky,iz) * (xphij(ip,ij)*kernel_i(ij, ikx, iky) &
- + xphijp1(ip,ij)*kernel_i(ij+1, ikx, iky) &
- + xphijm1(ip,ij)*kernel_i(ij-1, ikx, iky) )
+ Tphi = phi(ikx,iky,iz) * (xphij(ip,ij)*kernel_i(ij,ikx,iky,iz) &
+ + xphijp1(ip,ij)*kernel_i(ij+1,ikx,iky,iz) &
+ + xphijm1(ip,ij)*kernel_i(ij-1,ikx,iky,iz) )
ELSE
Tphi = 0._dp
ENDIF
@@ -218,8 +252,9 @@ SUBROUTINE moments_eq_rhs_i
!! Sum of linear terms
moments_rhs_i(ip,ij,ikx,iky,iz,updatetlevel) = &
- - Ckxky(ikx,iky,iz) * (Tnipj + Tnipp2j + Tnipm2j + Tnipjp1 + Tnipjm1)&
- - Tpari/2._dp/deltaz/q0 &
+ - imagu*Ckxky(ikx,iky,iz) * (Tnipj + Tnipp2j + Tnipm2j + Tnipjp1 + Tnipjm1)&
+ - Tpari/2._dp/deltaz*gradz_coeff(iz) &
+ - gradzB(iz)* Tmir *gradz_coeff(iz) &
+ i_ky * Tphi &
- mu*kperp2**2 * moments_i(ip,ij,ikx,iky,iz,updatetlevel) &
+ TColl
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index f247cecb..b215be7b 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -53,9 +53,9 @@ END SUBROUTINE build_dnjs_table
!******************************************************************************!
SUBROUTINE evaluate_kernels
USE basic
- USE array, Only : kernel_e, kernel_i
+ USE array, Only : kernel_e, kernel_i, gxy, gyy
USE grid
- use model, ONLY : tau_e, tau_i, sigma_e, sigma_i, q_e, q_i, lambdaD, CLOS, sigmae2_taue_o2, sigmai2_taui_o2
+ USE model, ONLY : tau_e, tau_i, sigma_e, sigma_i, q_e, q_i, lambdaD, CLOS, sigmae2_taue_o2, sigmai2_taui_o2
IMPLICIT NONE
REAL(dp) :: factj, j_dp, j_int
@@ -80,11 +80,11 @@ SUBROUTINE evaluate_kernels
kx = kxarray(ikx)
DO iky = ikys,ikye
ky = kyarray(iky)
-
- be_2 = (kx**2 + ky**2) * sigmae2_taue_o2
-
- kernel_e(ij, ikx, iky) = be_2**j_int * exp(-be_2)/factj
-
+ DO iz = izs,ize
+ kperp2= kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
+ be_2 = kperp2 * sigmae2_taue_o2
+ kernel_e(ij, ikx, iky,iz) = be_2**j_int * exp(-be_2)/factj
+ ENDDO
ENDDO
ENDDO
ENDDO
@@ -93,15 +93,18 @@ SUBROUTINE evaluate_kernels
kx = kxarray(ikx)
DO iky = ikys,ikye
ky = kyarray(iky)
- be_2 = (kx**2 + ky**2) * sigmae2_taue_o2
- ! Kernel ghost + 1 with Kj+1 = y/(j+1) Kj (/!\ ij = j+1)
- kernel_e(ijeg_e,ikx,iky) = be_2/(real(ijeg_e-1,dp))*kernel_e(ije_e,ikx,iky)
- ! Kernel ghost - 1 with Kj-1 = j/y Kj(careful about the kperp=0)
- IF ( be_2 .NE. 0 ) THEN
- kernel_e(ijsg_e,ikx,iky) = (real(ijsg_e,dp))/be_2*kernel_e(ijs_e,ikx,iky)
- ELSE
- kernel_e(ijsg_e,ikx,iky) = 0._dp
- ENDIF
+ DO iz = izs,ize
+ kperp2= kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
+ be_2 = kperp2 * sigmae2_taue_o2
+ ! Kernel ghost + 1 with Kj+1 = y/(j+1) Kj (/!\ ij = j+1)
+ kernel_e(ijeg_e,ikx,iky,iz) = be_2/(real(ijeg_e-1,dp))*kernel_e(ije_e,ikx,iky,iz)
+ ! Kernel ghost - 1 with Kj-1 = j/y Kj(careful about the kperp=0)
+ IF ( be_2 .NE. 0 ) THEN
+ kernel_e(ijsg_e,ikx,iky,iz) = (real(ijsg_e,dp))/be_2*kernel_e(ijs_e,ikx,iky,iz)
+ ELSE
+ kernel_e(ijsg_e,ikx,iky,iz) = 0._dp
+ ENDIF
+ ENDDO
ENDDO
ENDDO
@@ -122,9 +125,11 @@ SUBROUTINE evaluate_kernels
kx = kxarray(ikx)
DO iky = ikys,ikye
ky = kyarray(iky)
- bi_2 = (kx**2 + ky**2) * sigmai2_taui_o2
- kernel_i(ij, ikx, iky) = bi_2**j_int * exp(-bi_2)/factj
-
+ DO iz = izs,ize
+ kperp2= kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
+ bi_2 = kperp2 * sigmai2_taui_o2
+ kernel_i(ij, ikx, iky,iz) = bi_2**j_int * exp(-bi_2)/factj
+ ENDDO
ENDDO
ENDDO
ENDDO
@@ -133,24 +138,25 @@ SUBROUTINE evaluate_kernels
kx = kxarray(ikx)
DO iky = ikys,ikye
ky = kyarray(iky)
- bi_2 = (kx**2 + ky**2) * sigmai2_taui_o2
- ! Kernel ghost + 1 with Kj+1 = y/(j+1) Kj
- kernel_i(ijeg_i,ikx,iky) = bi_2/(real(ijeg_i-1,dp))*kernel_i(ije_i,ikx,iky)
- ! Kernel ghost - 1 with Kj-1 = j/y Kj(careful about the kperp=0)
- IF ( bi_2 .NE. 0 ) THEN
- kernel_i(ijsg_i,ikx,iky) = (real(ijsg_i,dp))/bi_2*kernel_i(ijs_i,ikx,iky)
- ELSE
- kernel_i(ijsg_i,ikx,iky) = 0._dp
- ENDIF
+ DO iz = izs,ize
+ kperp2= kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
+ bi_2 = kperp2 * sigmai2_taui_o2
+ ! Kernel ghost + 1 with Kj+1 = y/(j+1) Kj
+ kernel_i(ijeg_i,ikx,iky,iz) = bi_2/(real(ijeg_i-1,dp))*kernel_i(ije_i,ikx,iky,iz)
+ ! Kernel ghost - 1 with Kj-1 = j/y Kj(careful about the kperp=0)
+ IF ( bi_2 .NE. 0 ) THEN
+ kernel_i(ijsg_i,ikx,iky,iz) = (real(ijsg_i,dp))/bi_2*kernel_i(ijs_i,ikx,iky,iz)
+ ELSE
+ kernel_i(ijsg_i,ikx,iky,iz) = 0._dp
+ ENDIF
+ ENDDO
ENDDO
ENDDO
END SUBROUTINE evaluate_kernels
!******************************************************************************!
SUBROUTINE compute_lin_coeff
- USE array, ONLY: xnepj, xnepp1j, xnepm1j, xnepp2j, xnepm2j, xnepjp1, xnepjm1, &
- xnipj, xnipp1j, xnipm1j, xnipp2j, xnipm2j, xnipjp1, xnipjm1, &
- xphij, xphijp1, xphijm1
+ USE array
USE model, ONLY: taue_qe, taui_qi, sqrtTaue_qe, sqrtTaui_qi, eta_T, eta_n
USE prec_const
USE grid, ONLY: parray_e, parray_i, jarray_e, jarray_i, &
@@ -159,7 +165,6 @@ SUBROUTINE compute_lin_coeff
INTEGER :: p_int, j_int ! polynom. degrees
REAL(dp) :: p_dp, j_dp
REAL(dp) :: kx, ky, z
-
!! Electrons linear coefficients for moment RHS !!!!!!!!!!
DO ip = ips_e, ipe_e
p_int= parray_e(ip) ! Hermite degree
@@ -168,6 +173,13 @@ SUBROUTINE compute_lin_coeff
j_int= jarray_e(ij) ! Laguerre degree
j_dp = REAL(j_int,dp) ! REAL of Laguerre degree
xnepj(ip,ij) = taue_qe * 2._dp * (p_dp + j_dp + 1._dp)
+ ynepp1j (ip,ij) = -SQRT(tau_e)/sigma_e * (j_dp+1)*SQRT(p_dp+1._dp)
+ ynepp1jm1(ip,ij) = +SQRT(tau_e)/sigma_e * j_dp*SQRT(p_dp+1._dp)
+ ynepm1j (ip,ij) = -SQRT(tau_e)/sigma_e * (j_dp+1)*SQRT(p_dp)
+ ynepm1jm1(ip,ij) = +SQRT(tau_e)/sigma_e * j_dp*SQRT(p_dp)
+ zNepm1j (ip,ij) = -SQRT(tau_e)/sigma_e * (2._dp*j_dp+1_dp)*SQRT(2._dp*p_dp)
+ zNepm1jp1(ip,ij) = +SQRT(tau_e)/sigma_e * (j_dp+1_dp)*SQRT(2._dp*p_dp)
+ zNepm1jm1(ip,ij) = +SQRT(tau_e)/sigma_e * j_dp*SQRT(2._dp*p_dp)
ENDDO
ENDDO
DO ip = ips_e, ipe_e
@@ -193,6 +205,13 @@ SUBROUTINE compute_lin_coeff
j_int= jarray_i(ij) ! Laguerre degree
j_dp = REAL(j_int,dp) ! REAL of Laguerre degree
xnipj(ip,ij) = taui_qi * 2._dp * (p_dp + j_dp + 1._dp)
+ ynipp1j (ip,ij) = -SQRT(tau_i)/sigma_i * (j_dp+1)*SQRT(p_dp+1._dp)
+ ynipp1jm1(ip,ij) = +SQRT(tau_i)/sigma_i * j_dp*SQRT(p_dp+1._dp)
+ ynipm1j (ip,ij) = -SQRT(tau_i)/sigma_i * (j_dp+1)*SQRT(p_dp)
+ ynipm1jm1(ip,ij) = +SQRT(tau_i)/sigma_i * j_dp*SQRT(p_dp)
+ zNipm1j (ip,ij) = -SQRT(tau_i)/sigma_i * (2._dp*j_dp+1_dp)*SQRT(2._dp*p_dp)
+ zNipm1jp1(ip,ij) = +SQRT(tau_i)/sigma_i * (j_dp+1_dp)*SQRT(2._dp*p_dp)
+ zNipm1jm1(ip,ij) = +SQRT(tau_i)/sigma_i * j_dp*SQRT(2._dp*p_dp)
ENDDO
ENDDO
DO ip = ips_i, ipe_i
--
GitLab