From c082a7a664a9a37acc12bf9aa8fe626177add223 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Mon, 27 Feb 2023 18:02:44 +0100
Subject: [PATCH] add dv4 framework
---
src/auxval.F90 | 2 +
src/collision_mod.F90 | 373 ++++++++++++++-----------------------
src/model_mod.F90 | 14 +-
src/moments_eq_rhs_mod.F90 | 42 +++--
src/processing_mod.F90 | 22 ++-
5 files changed, 188 insertions(+), 265 deletions(-)
diff --git a/src/auxval.F90 b/src/auxval.F90
index f9817aef..b236c48e 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -50,6 +50,8 @@ subroutine auxval
CALL build_dnjs_table ! precompute the Laguerre nonlin product coeffs
ENDIF
+ CALL build_dv4Hp_table ! precompute the hermite fourth derivative table
+
!! Display parallel settings
DO i_ = 0,num_procs-1
CALL mpi_barrier(MPI_COMM_WORLD, ierr)
diff --git a/src/collision_mod.F90 b/src/collision_mod.F90
index bdbf6f5c..c2286157 100644
--- a/src/collision_mod.F90
+++ b/src/collision_mod.F90
@@ -13,7 +13,7 @@ IMPLICIT NONE
PRIVATE
! Set the collision model to use
! (Lenard-Bernstein: 'LB', Dougherty: 'DG', Sugama: 'SG', Lorentz: 'LR', Landau: 'LD')
-CHARACTER(len=2),PUBLIC,PROTECTED :: collision_model
+CHARACTER(len=16),PUBLIC,PROTECTED :: collision_model
LOGICAL, PUBLIC, PROTECTED :: gyrokin_CO =.true. ! activates GK effects on CO
LOGICAL, PUBLIC, PROTECTED :: interspecies =.true. ! activates interpecies collision
CHARACTER(len=128), PUBLIC :: mat_file ! COSOlver matrix file names
@@ -23,8 +23,6 @@ LOGICAL, PUBLIC, PROTECTED :: cosolver_coll ! check if cosolver matrices are use
PUBLIC :: collision_readinputs, coll_outputinputs
PUBLIC :: compute_TColl
PUBLIC :: compute_lenard_bernstein, compute_dougherty
-PUBLIC :: LenardBernstein_e, LenardBernstein_i!, LenardBernstein GK
-PUBLIC :: DoughertyGK_ee, DoughertyGK_ii!, Dougherty GK
PUBLIC :: load_COSOlver_mat, compute_cosolver_coll
PUBLIC :: apply_COSOlver_mat_e, apply_COSOlver_mat_i
@@ -49,7 +47,7 @@ CONTAINS
interspecies = .false.
CASE ('LD') ! Landau
cosolver_coll = .true.
- CASE ('none')
+ CASE ('none','hypcoll','dvpar4')
cosolver_coll = .false.
interspecies = .false.
CASE DEFAULT
@@ -90,7 +88,7 @@ CONTAINS
CALL compute_dougherty
CASE ('SG','LR','LD')
CALL compute_cosolver_coll
- CASE ('none')
+ CASE ('none','hypcoll','dvpar4')
IF(KIN_E) &
TColl_e = 0._dp
TColl_i = 0._dp
@@ -193,284 +191,187 @@ CONTAINS
!******************************************************************************!
SUBROUTINE compute_dougherty
IMPLICIT NONE
- COMPLEX(dp) :: TColl_
- IF (KIN_E) THEN
- DO iz = izs,ize
- DO iky = ikys, ikye;
- DO ikx = ikxs, ikxe;
- DO ij = ijs_e,ije_e
- DO ip = ips_e,ipe_e;
- IF(gyrokin_CO) THEN
- CALL DoughertyGK_ee(ip,ij,iky,ikx,iz,TColl_)
- ELSE
- CALL DoughertyDK_ee(ip,ij,iky,ikx,iz,TColl_)
- ENDIF
- TColl_e(ip,ij,iky,ikx,iz) = TColl_
- ENDDO;ENDDO;ENDDO
- ENDDO;ENDDO
+ IF(gyrokin_CO) THEN
+ if(KIN_E) &
+ CALL DoughertyGK_aa(ips_e,ipe_e,ijs_e,ije_e,ijgs_e,ijge_e,ip0_e,ip1_e,ip2_e,Jmaxe,&
+ parray_e(ips_e:ipe_e),jarray_e(ijs_e:ije_e),&
+ kernel_e(ijgs_e:ijge_e,ikys:ikye,ikxs:ikxe,izs:ize,0:1),&
+ nadiab_moments_e(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikxs:ikxe,izs:ize),&
+ nu_ee,sigmae2_taue_o2,sqrt_sigmae2_taue_o2,TColl_e)
+ CALL DoughertyGK_aa(ips_i,ipe_i,ijs_i,ije_i,ijgs_i,ijge_i,ip0_i,ip1_i,ip2_i,Jmaxi,&
+ parray_i(ips_i:ipe_i),jarray_i(ijs_i:ije_i),&
+ kernel_i(ijgs_i:ijge_i,ikys:ikye,ikxs:ikxe,izs:ize,0:1),&
+ nadiab_moments_i(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize),&
+ nu_i,sigmai2_taui_o2,sqrt_sigmai2_taui_o2,TColl_i)
+ ELSE
+ if(KIN_E) &
+ CALL DoughertyDK_aa(ips_e,ipe_e,ijs_e,ije_e,ip0_e,ip1_e,ip2_e,&
+ parray_e(ips_e:ipe_e),jarray_e(ijs_e:ije_e),&
+ moments_e(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel),&
+ nu_ee,TColl_e)
+ CALL DoughertyDK_aa(ips_i,ipe_i,ijs_i,ije_i,ip0_i,ip1_i,ip2_i,&
+ parray_i(ips_i:ipe_i),jarray_i(ijs_i:ije_i),&
+ moments_i(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel),&
+ nu_i,TColl_i)
ENDIF
- DO iz = izs,ize
- DO iky = ikys, ikye;
- DO ikx = ikxs, ikxe;
- DO ij = ijs_i,ije_i
- DO ip = ips_i,ipe_i;
- IF(gyrokin_CO) THEN
- CALL DoughertyGK_ii(ip,ij,iky,ikx,iz,TColl_)
- ELSE
- CALL DoughertyDK_ii(ip,ij,iky,ikx,iz,TColl_)
- ENDIF
- TColl_i(ip,ij,iky,ikx,iz) = TColl_
- ENDDO;ENDDO;ENDDO
- ENDDO;ENDDO
END SUBROUTINE compute_dougherty
!******************************************************************************!
- !! Doughtery driftkinetic collision operator for electrons
- !******************************************************************************!
- SUBROUTINE DoughertyDK_ee(ip_,ij_,iky_,ikx_,iz_,TColl_)
- IMPLICIT NONE
- INTEGER, INTENT(IN) :: ip_,ij_,iky_,ikx_,iz_
- COMPLEX(dp), INTENT(OUT) :: TColl_
- REAL(dp) :: j_dp, p_dp
- !** Auxiliary variables **
- p_dp = REAL(parray_e(ip_),dp)
- j_dp = REAL(jarray_e(ij_),dp)
- !** Assembling collison operator **
- TColl_ = -(p_dp + 2._dp*j_dp)*moments_e(ip_,ij_,iky_,ikx_,iz_,updatetlevel)
- IF( (p_dp .EQ. 1._dp) .AND. (j_dp .EQ. 0._dp)) THEN !Ce10
- TColl_ = TColl_ + moments_e(ip1_e,1,iky_,ikx_,iz_,updatetlevel)
- ELSEIF( (p_dp .EQ. 2._dp) .AND. (j_dp .EQ. 0._dp)) THEN ! Ce20
- TColl_ = TColl_ + twothird*moments_e(ip2_e,1,iky_,ikx_,iz_,updatetlevel) &
- - SQRT2*twothird*moments_e(ip0_e,2,iky_,ikx_,iz_,updatetlevel)
- ELSEIF( (p_dp .EQ. 0._dp) .AND. (j_dp .EQ. 1._dp)) THEN ! Ce01
- TColl_ = TColl_ + 2._dp*twothird*moments_e(ip0_e,2,iky_,ikx_,iz_,updatetlevel) &
- - SQRT2*twothird*moments_e(ip2_e,1,iky_,ikx_,iz_,updatetlevel)
- ENDIF
- TColl_ = nu_ee * TColl_
- END SUBROUTINE DoughertyDK_ee
- !******************************************************************************!
- !! Doughtery driftkinetic collision operator for ions
+ !! Doughtery driftkinetic collision operator (like-species)
!******************************************************************************!
- SUBROUTINE DoughertyDK_ii(ip_,ij_,iky_,ikx_,iz_,TColl_)
+ SUBROUTINE DoughertyDK_aa(ips_,ipe_,ijs_,ije_,ip0_,ip1_,ip2_,&
+ parray_,jarray_,moments_,nu_,Tcoll_)
+
IMPLICIT NONE
- INTEGER, INTENT(IN) :: ip_,ij_,iky_,ikx_,iz_
- COMPLEX(dp), INTENT(OUT) :: TColl_
+ !! INPUTS
+ INTEGER, INTENT(IN) :: ips_, ipe_, ijs_, ije_
+ INTEGER, INTENT(IN) :: ip0_, ip1_, ip2_
+ INTEGER, DIMENSION(ips_:ipe_), INTENT(IN) :: parray_
+ INTEGER, DIMENSION(ijs_:ije_), INTENT(IN) :: jarray_
+ COMPLEX(dp), DIMENSION(ips_:ipe_,ijs_:ije_,ikys:ikye,ikxs:ikxe,izs:ize),INTENT(IN) :: moments_
+ REAL(dp), INTENT(IN) :: nu_
+ !! OUTPUT
+ COMPLEX(dp), DIMENSION(ips_:ipe_, ijs_:ije_,ikys:ikye,ikxs:ikxe, izs:ize), INTENT(OUT) :: TColl_
+ !! Local variables
REAL(dp) :: j_dp, p_dp
- !** Auxiliary variables **
- p_dp = REAL(parray_i(ip_),dp)
- j_dp = REAL(jarray_i(ij_),dp)
- !** Assembling collison operator **
- TColl_ = -(p_dp + 2._dp*j_dp)*moments_i(ip_,ij_,iky_,ikx_,iz_,updatetlevel)
- IF( (p_dp .EQ. 1._dp) .AND. (j_dp .EQ. 0._dp)) THEN ! kronecker p1j0
- TColl_ = TColl_ + moments_i(ip1_i,1,iky_,ikx_,iz_,updatetlevel)
- ELSEIF( (p_dp .EQ. 2._dp) .AND. (j_dp .EQ. 0._dp)) THEN ! kronecker p2j0
- TColl_ = TColl_ + twothird*moments_i(ip2_i,1,iky_,ikx_,iz_,updatetlevel) &
- - SQRT2*twothird*moments_i(ip0_i,2,iky_,ikx_,iz_,updatetlevel)
- ELSEIF( (p_dp .EQ. 0._dp) .AND. (j_dp .EQ. 1._dp)) THEN ! kronecker p0j1
- TColl_ = TColl_ + 2._dp*twothird*moments_i(ip0_i,2,iky_,ikx_,iz_,updatetlevel) &
- - SQRT2*twothird*moments_i(ip2_i,1,iky_,ikx_,iz_,updatetlevel)
- ENDIF
- TColl_ = nu_i * TColl_
- END SUBROUTINE DoughertyDK_ii
+ COMPLEX(dp) :: Tmp
+ DO iz = izs,ize
+ DO iky = ikys, ikye;
+ DO ikx = ikxs, ikxe;
+ DO ij = ijs_,ije_
+ DO ip = ips_,ipe_;
+ !** Auxiliary variables **
+ p_dp = REAL(parray_(ip),dp)
+ j_dp = REAL(jarray_(ij),dp)
+ !** Assembling collison operator **
+ Tmp = -(p_dp + 2._dp*j_dp)*moments_(ip,ij,iky,ikx,iz)
+ IF( (p_dp .EQ. 1._dp) .AND. (j_dp .EQ. 0._dp)) THEN !Ce10
+ Tmp = Tmp + moments_(ip1_,1,iky,ikx,iz)
+ ELSEIF( (p_dp .EQ. 2._dp) .AND. (j_dp .EQ. 0._dp)) THEN ! Ce20
+ Tmp = Tmp + twothird*moments_(ip2_,1,iky,ikx,iz) &
+ - SQRT2*twothird*moments_(ip0_,2,iky,ikx,iz)
+ ELSEIF( (p_dp .EQ. 0._dp) .AND. (j_dp .EQ. 1._dp)) THEN ! Ce01
+ Tmp = Tmp + 2._dp*twothird*moments_(ip0_,2,iky,ikx,iz) &
+ - SQRT2*twothird*moments_(ip2_,1,iky,ikx,iz)
+ ENDIF
+ TColl_(ip,ij,iky,ikx,iz) = nu_ * Tmp
+ ENDDO;ENDDO;ENDDO
+ ENDDO;ENDDO
+ END SUBROUTINE DoughertyDK_aa
+
!******************************************************************************!
- !! Doughtery gyrokinetic collision operator for electrons
+ !! Doughtery gyrokinetic collision operator (species like)
!******************************************************************************!
- SUBROUTINE DoughertyGK_ee(ip_,ij_,iky_,ikx_,iz_,TColl_)
+ SUBROUTINE DoughertyGK_aa(ips_,ipe_,ijs_,ije_,ijgs_,ijge_,ip0_,ip1_,ip2_,jmax_,&
+ parray_,jarray_,kernel_,nadiab_moments_,&
+ nu_,sigmaa2_taua_o2, sqrt_sigmaa2_taua_o2,Tcoll_)
IMPLICIT NONE
- INTEGER, INTENT(IN) :: ip_,ij_,iky_,ikx_,iz_
- COMPLEX(dp), INTENT(OUT) :: TColl_
-
- COMPLEX(dp) :: n_,upar_,uperp_,Tpar_, Tperp_, T_
- COMPLEX(dp) :: nadiab_moment_0j
+ !! INPUTS
+ INTEGER, INTENT(IN) :: ips_, ipe_, ijs_, ije_, ijgs_, ijge_
+ INTEGER, INTENT(IN) :: ip0_, ip1_, ip2_,jmax_
+ INTEGER, DIMENSION(ips_:ipe_), INTENT(IN) :: parray_
+ INTEGER, DIMENSION(ijs_:ije_), INTENT(IN) :: jarray_
+ REAL(dp), DIMENSION(ijgs_:ijge_,ikys:ikye,ikxs:ikxe,izgs:izge,0:1), INTENT(IN) :: kernel_
+ COMPLEX(dp), DIMENSION(ips_:ipe_,ijs_:ije_,ikys:ikye,ikxs:ikxe,izs:ize),INTENT(IN) :: nadiab_moments_
+ REAL(dp), INTENT(IN) :: nu_, sigmaa2_taua_o2, sqrt_sigmaa2_taua_o2
+ !! OUTPUT
+ COMPLEX(dp), DIMENSION(ips_:ipe_,ijs_:ije_,ikys:ikye,ikxs:ikxe,izs:ize), INTENT(OUT) :: TColl_
+ !! Local variables
+ COMPLEX(dp) :: dens,upar,uperp,Tpar, Tperp, Temp
+ COMPLEX(dp) :: nadiab_moment_0j, Tmp
REAL(dp) :: Knp0, Knp1, Knm1, kp
- INTEGER :: in_, eo_
- REAL(dp) :: n_dp, j_dp, p_dp, be_, be_2
-
+ INTEGER :: in, eo
+ REAL(dp) :: n_dp, j_dp, p_dp, ba, ba_2
+ DO iz = izs,ize
+ DO iky = ikys, ikye;
+ DO ikx = ikxs, ikxe;
+ DO ij = ijs_,ije_
+ DO ip = ips_,ipe_;
!** Auxiliary variables **
- p_dp = REAL(parray_e(ip_),dp)
- eo_ = MODULO(parray_e(ip_),2)
- j_dp = REAL(jarray_e(ij_),dp)
- kp = kparray(iky_,ikx_,iz_,eo_)
- be_2 = kp**2 * sigmae2_taue_o2 ! this is (be/2)^2
- be_ = 2_dp*kp * sqrt_sigmae2_taue_o2 ! this is be
+ p_dp = REAL(parray_(ip),dp)
+ eo = MODULO(parray_(ip),2)
+ j_dp = REAL(jarray_(ij),dp)
+ kp = kparray(iky,ikx,iz,eo)
+ ba_2 = kp**2 * sigmaa2_taua_o2 ! this is (l_a/2)^2
+ ba = 2_dp*kp * sqrt_sigmaa2_taua_o2 ! this is l_a
!** Assembling collison operator **
! Velocity-space diffusion (similar to Lenard Bernstein)
- ! -nuee (p + 2j + b^2/2) Nepj
- TColl_ = -(p_dp + 2._dp*j_dp + 2._dp*be_2)*nadiab_moments_e(ip_,ij_,iky_,ikx_,iz_)
+ ! -nu (p + 2j + b^2/2) Napj
+ Tmp = -(p_dp + 2._dp*j_dp + 2._dp*ba_2)*nadiab_moments_(ip,ij,iky,ikx,iz)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
IF( p_dp .EQ. 0 ) THEN ! Kronecker p0
!** build required fluid moments **
- n_ = 0._dp
- upar_ = 0._dp; uperp_ = 0._dp
- Tpar_ = 0._dp; Tperp_ = 0._dp
- DO in_ = 1,jmaxe+1
- n_dp = REAL(in_-1,dp)
+ dens = 0._dp
+ upar = 0._dp; uperp = 0._dp
+ Tpar = 0._dp; Tperp = 0._dp
+ DO in = 1,jmax_+1
+ n_dp = REAL(in-1,dp)
! Store the kernels for sparing readings
- Knp0 = Kernel_e(in_ ,iky_,ikx_,iz_,eo_)
- Knp1 = Kernel_e(in_+1,iky_,ikx_,iz_,eo_)
- Knm1 = Kernel_e(in_-1,iky_,ikx_,iz_,eo_)
+ Knp0 = Kernel_(in ,iky,ikx,iz,eo)
+ Knp1 = Kernel_(in+1,iky,ikx,iz,eo)
+ Knm1 = Kernel_(in-1,iky,ikx,iz,eo)
! Nonadiabatic moments (only different from moments when p=0)
- nadiab_moment_0j = nadiab_moments_e(ip0_e,in_,iky_,ikx_,iz_)
+ nadiab_moment_0j = nadiab_moments_(ip0_,in,iky,ikx,iz)
! Density
- n_ = n_ + Knp0 * nadiab_moment_0j
+ dens = dens + Knp0 * nadiab_moment_0j
! Perpendicular velocity
- uperp_ = uperp_ + be_*0.5_dp*(Knp0 - Knm1) * nadiab_moment_0j
+ uperp = uperp + ba*0.5_dp*(Knp0 - Knm1) * nadiab_moment_0j
! Parallel temperature
- Tpar_ = Tpar_ + Knp0 * (SQRT2*nadiab_moments_e(ip2_e,in_,iky_,ikx_,iz_) + nadiab_moment_0j)
+ Tpar = Tpar + Knp0 * (SQRT2*nadiab_moments_(ip2_,in,iky,ikx,iz) + nadiab_moment_0j)
! Perpendicular temperature
- Tperp_ = Tperp_ + ((2._dp*n_dp+1._dp)*Knp0 - (n_dp+1._dp)*Knp1 - n_dp*Knm1)*nadiab_moment_0j
+ Tperp = Tperp + ((2._dp*n_dp+1._dp)*Knp0 - (n_dp+1._dp)*Knp1 - n_dp*Knm1)*nadiab_moment_0j
ENDDO
- T_ = (Tpar_ + 2._dp*Tperp_)/3._dp - n_
+ Temp = (Tpar + 2._dp*Tperp)/3._dp - dens
! Add energy restoring term
- TColl_ = TColl_ + T_* 4._dp * j_dp * Kernel_e(ij_ ,iky_,ikx_,iz_,eo_)
- TColl_ = TColl_ - T_* 2._dp * (j_dp + 1._dp) * Kernel_e(ij_+1,iky_,ikx_,iz_,eo_)
- TColl_ = TColl_ - T_* 2._dp * j_dp * Kernel_e(ij_-1,iky_,ikx_,iz_,eo_)
- TColl_ = TColl_ + uperp_*be_* (Kernel_e(ij_,iky_,ikx_,iz_,eo_) - Kernel_e(ij_-1,iky_,ikx_,iz_,eo_))
+ Tmp = Tmp + Temp* 4._dp * j_dp * Kernel_(ij ,iky,ikx,iz,eo)
+ Tmp = Tmp - Temp* 2._dp * (j_dp + 1._dp) * Kernel_(ij+1,iky,ikx,iz,eo)
+ Tmp = Tmp - Temp* 2._dp * j_dp * Kernel_(ij-1,iky,ikx,iz,eo)
+ Tmp = Tmp + uperp*ba* (Kernel_(ij,iky,ikx,iz,eo) - Kernel_(ij-1,iky,ikx,iz,eo))
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
ELSEIF( p_dp .eq. 1 ) THEN ! kronecker p1
!** build required fluid moments **
- upar_ = 0._dp
- DO in_ = 1,jmaxe+1
+ upar = 0._dp
+ DO in = 1,jmax_+1
! Parallel velocity
- upar_ = upar_ + Kernel_e(in_,iky_,ikx_,iz_,eo_) * nadiab_moments_e(ip1_e,in_,iky_,ikx_,iz_)
+ upar = upar + Kernel_(in,iky,ikx,iz,eo) * nadiab_moments_(ip1_,in,iky,ikx,iz)
ENDDO
- TColl_ = TColl_ + upar_*Kernel_e(ij_,iky_,ikx_,iz_,eo_)
+ Tmp = Tmp + upar*Kernel_(ij,iky,ikx,iz,eo)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
ELSEIF( p_dp .eq. 2 ) THEN ! kronecker p2
!** build required fluid moments **
- n_ = 0._dp
- upar_ = 0._dp; uperp_ = 0._dp
- Tpar_ = 0._dp; Tperp_ = 0._dp
- DO in_ = 1,jmaxe+1
- n_dp = REAL(in_-1,dp)
+ dens = 0._dp
+ upar = 0._dp; uperp = 0._dp
+ Tpar = 0._dp; Tperp = 0._dp
+ DO in = 1,jmax_+1
+ n_dp = REAL(in-1,dp)
! Store the kernels for sparing readings
- Knp0 = Kernel_e(in_ ,iky_,ikx_,iz_,eo_)
- Knp1 = Kernel_e(in_+1,iky_,ikx_,iz_,eo_)
- Knm1 = Kernel_e(in_-1,iky_,ikx_,iz_,eo_)
+ Knp0 = Kernel_(in ,iky,ikx,iz,eo)
+ Knp1 = Kernel_(in+1,iky,ikx,iz,eo)
+ Knm1 = Kernel_(in-1,iky,ikx,iz,eo)
! Nonadiabatic moments (only different from moments when p=0)
- nadiab_moment_0j = nadiab_moments_e(ip0_e,in_,iky_,ikx_,iz_)
+ nadiab_moment_0j = nadiab_moments_(ip0_,in,iky,ikx,iz)
! Density
- n_ = n_ + Knp0 * nadiab_moment_0j
+ dens = dens + Knp0 * nadiab_moment_0j
! Parallel temperature
- Tpar_ = Tpar_ + Knp0 * (SQRT2*nadiab_moments_e(ip2_e,in_,iky_,ikx_,iz_) + nadiab_moment_0j)
+ Tpar = Tpar + Knp0 * (SQRT2*nadiab_moments_(ip2_,in,iky,ikx,iz) + nadiab_moment_0j)
! Perpendicular temperature
- Tperp_ = Tperp_ + ((2._dp*n_dp+1._dp)*Knp0 - (n_dp+1._dp)*Knp1 - n_dp*Knm1)*nadiab_moment_0j
+ Tperp = Tperp + ((2._dp*n_dp+1._dp)*Knp0 - (n_dp+1._dp)*Knp1 - n_dp*Knm1)*nadiab_moment_0j
ENDDO
- T_ = (Tpar_ + 2._dp*Tperp_)/3._dp - n_
- TColl_ = TColl_ + T_*SQRT2*Kernel_e(ij_,iky_,ikx_,iz_,eo_)
+ Temp = (Tpar + 2._dp*Tperp)/3._dp - dens
+ Tmp = Tmp + Temp*SQRT2*Kernel_(ij,iky,ikx,iz,eo)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
ENDIF
- ! Multiply by electron-electron collision coefficient
- TColl_ = nu_ee * TColl_
-
- END SUBROUTINE DoughertyGK_ee
- !******************************************************************************!
- !! Doughtery gyrokinetic collision operator for ions
- !******************************************************************************!
- SUBROUTINE DoughertyGK_ii(ip_,ij_,iky_,ikx_,iz_,TColl_)
- IMPLICIT NONE
- INTEGER, INTENT(IN) :: ip_,ij_,iky_,ikx_,iz_
- COMPLEX(dp), INTENT(OUT) :: TColl_
-
- COMPLEX(dp) :: n_,upar_,uperp_,Tpar_, Tperp_, T_
- COMPLEX(dp) :: bi_, bi_2
- COMPLEX(dp) :: nadiab_moment_0j
- REAL(dp) :: Knp0, Knp1, Knm1, kp
- INTEGER :: in_, eo_
- REAL(dp) :: n_dp, j_dp, p_dp
-
- !** Auxiliary variables **
- p_dp = REAL(parray_i(ip_),dp)
- eo_ = MODULO(parray_i(ip_),2)
- j_dp = REAL(jarray_i(ij_),dp)
- kp = kparray(iky_,ikx_,iz_,eo_)
- bi_2 = kp**2 *sigmai2_taui_o2 ! this is (bi/2)^2
- bi_ = 2_dp*kp*sqrt_sigmai2_taui_o2 ! this is bi
-
- !** Assembling collison operator **
- ! Velocity-space diffusion (similar to Lenard Bernstein)
- ! -nui (p + 2j + b^2/2) Nipj
- TColl_ = -(p_dp + 2._dp*j_dp + 2._dp*bi_2)*nadiab_moments_i(ip_,ij_,iky_,ikx_,iz_)
-
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- IF( p_dp .EQ. 0 ) THEN ! Kronecker p0
- !** build required fluid moments **
- n_ = 0._dp
- upar_ = 0._dp; uperp_ = 0._dp
- Tpar_ = 0._dp; Tperp_ = 0._dp
- DO in_ = 1,jmaxi+1
- n_dp = REAL(in_-1,dp)
- ! Store the kernels for sparing readings
- Knp0 = Kernel_i(in_ ,iky_,ikx_,iz_,eo_)
- Knp1 = Kernel_i(in_+1,iky_,ikx_,iz_,eo_)
- Knm1 = Kernel_i(in_-1,iky_,ikx_,iz_,eo_)
- ! Nonadiabatic moments (only different from moments when p=0)
- nadiab_moment_0j = nadiab_moments_i(ip0_i,in_,iky_,ikx_,iz_)
- ! Density
- n_ = n_ + Knp0 * nadiab_moment_0j
- ! Perpendicular velocity
- uperp_ = uperp_ + bi_*0.5_dp*(Knp0 - Knm1) * nadiab_moment_0j
- ! Parallel temperature
- Tpar_ = Tpar_ + Knp0 * (SQRT2*nadiab_moments_i(ip2_i,in_,iky_,ikx_,iz_) + nadiab_moment_0j)
- ! Perpendicular temperature
- Tperp_ = Tperp_ + ((2._dp*n_dp+1._dp)*Knp0 - (n_dp+1._dp)*Knp1 - n_dp*Knm1)*nadiab_moment_0j
- ENDDO
- T_ = (Tpar_ + 2._dp*Tperp_)*onethird - n_
- ! Add energy restoring term
- TColl_ = TColl_ + T_* 4._dp * j_dp * Kernel_i(ij_ ,iky_,ikx_,iz_,eo_)
- TColl_ = TColl_ - T_* 2._dp * (j_dp + 1._dp) * Kernel_i(ij_+1,iky_,ikx_,iz_,eo_)
- TColl_ = TColl_ - T_* 2._dp * j_dp * Kernel_i(ij_-1,iky_,ikx_,iz_,eo_)
- TColl_ = TColl_ + uperp_*bi_* (Kernel_i(ij_,iky_,ikx_,iz_,eo_) - Kernel_i(ij_-1,iky_,ikx_,iz_,eo_))
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ELSEIF( p_dp .eq. 1 ) THEN ! kxonecker p1
- !** build required fluid moments **
- upar_ = 0._dp
- DO in_ = 1,jmaxi+1
- ! Parallel velocity
- upar_ = upar_ + Kernel_i(in_,iky_,ikx_,iz_,eo_) * nadiab_moments_i(ip1_i,in_,iky_,ikx_,iz_)
- ENDDO
- TColl_ = TColl_ + upar_*Kernel_i(ij_,iky_,ikx_,iz_,eo_)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ELSEIF( p_dp .eq. 2 ) THEN ! kxonecker p2
- !** build required fluid moments **
- n_ = 0._dp
- upar_ = 0._dp; uperp_ = 0._dp
- Tpar_ = 0._dp; Tperp_ = 0._dp
- DO in_ = 1,jmaxi+1
- n_dp = REAL(in_-1,dp)
- ! Store the kernels for sparing readings
- Knp0 = Kernel_i(in_ ,iky_,ikx_,iz_,eo_)
- Knp1 = Kernel_i(in_+1,iky_,ikx_,iz_,eo_)
- Knm1 = Kernel_i(in_-1,iky_,ikx_,iz_,eo_)
- ! Nonadiabatic moments (only different from moments when p=0)
- nadiab_moment_0j = nadiab_moments_i(ip0_i,in_,iky_,ikx_,iz_)
- ! Density
- n_ = n_ + Knp0 * nadiab_moment_0j
- ! Parallel temperature
- Tpar_ = Tpar_ + Knp0 * (SQRT2*nadiab_moments_i(ip2_i,in_,iky_,ikx_,iz_) + nadiab_moment_0j)
- ! Perpendicular temperature
- Tperp_ = Tperp_ + ((2._dp*n_dp+1._dp)*Knp0 - (n_dp+1._dp)*Knp1 - n_dp*Knm1)*nadiab_moment_0j
- ENDDO
- T_ = (Tpar_ + 2._dp*Tperp_)*onethird - n_
- TColl_ = TColl_ + T_*SQRT2*Kernel_i(ij_,iky_,ikx_,iz_,eo_)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- ENDIF
- ! Multiply by ion-ion collision coefficient
- TColl_ = nu_i * TColl_
-
- END SUBROUTINE DoughertyGK_ii
-
+ ! Multiply by collision parameter
+ TColl_(ip,ij,iky,ikx,iz) = nu_ * Tmp
+ ENDDO;ENDDO;ENDDO
+ ENDDO;ENDDO
+ END SUBROUTINE DoughertyGK_aa
!******************************************************************************!
!! compute the collision terms in a (Np x Nj x Nkx x Nky) matrix all at once
!******************************************************************************!
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index 2bb73616..400ec152 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -1,6 +1,5 @@
MODULE model
! Module for diagnostic parameters
-
USE prec_const
IMPLICIT NONE
PRIVATE
@@ -11,12 +10,15 @@ MODULE model
CHARACTER(len=16), &
PUBLIC, PROTECTED ::LINEARITY= 'linear' ! To turn on non linear bracket term
LOGICAL, PUBLIC, PROTECTED :: KIN_E = .true. ! Simulate kinetic electron (adiabatic otherwise)
+ INTEGER, PUBLIC, PROTECTED :: N_HD = 4 ! order of numerical spatial diffusion
REAL(dp), PUBLIC, PROTECTED :: mu_x = 0._dp ! spatial x-Hyperdiffusivity coefficient (for num. stability)
REAL(dp), PUBLIC, PROTECTED :: mu_y = 0._dp ! spatial y-Hyperdiffusivity coefficient (for num. stability)
+ LOGICAL, PUBLIC, PROTECTED :: HDz_h = .false. ! to apply z-hyperdiffusion on non adiab part
REAL(dp), PUBLIC, PROTECTED :: mu_z = 0._dp ! spatial z-Hyperdiffusivity coefficient (for num. stability)
+ CHARACTER(len=16), &
+ PUBLIC, PROTECTED :: HYP_V = 'hypcoll' ! hyperdiffusion model for velocity space ('none','hypcoll','dvpar4')
REAL(dp), PUBLIC, PROTECTED :: mu_p = 0._dp ! kinetic para hyperdiffusivity coefficient (for num. stability)
REAL(dp), PUBLIC, PROTECTED :: mu_j = 0._dp ! kinetic perp hyperdiffusivity coefficient (for num. stability)
- INTEGER, PUBLIC, PROTECTED :: N_HD = 4 ! order of numerical spatial diffusion
REAL(dp), PUBLIC, PROTECTED :: nu = 0._dp ! Collision frequency
REAL(dp), PUBLIC, PROTECTED :: tau_e = 1._dp ! Temperature
REAL(dp), PUBLIC, PROTECTED :: tau_i = 1._dp !
@@ -59,17 +61,21 @@ CONTAINS
SUBROUTINE model_readinputs
! Read the input parameters
- USE basic, ONLY : lu_in, my_id
+ USE basic, ONLY : lu_in, my_id, num_procs_p
USE prec_const
IMPLICIT NONE
NAMELIST /MODEL_PAR/ CLOS, NL_CLOS, KERN, LINEARITY, KIN_E, &
- mu_x, mu_y, N_HD, mu_z, mu_p, mu_j, nu,&
+ mu_x, mu_y, N_HD, HDz_h, mu_z, mu_p, mu_j, HYP_V, nu,&
tau_e, tau_i, sigma_e, sigma_i, q_e, q_i,&
k_Ne, k_Ni, k_Te, k_Ti, k_gB, k_cB, lambdaD, beta
READ(lu_in,model_par)
+ IF((HYP_V .EQ. 'dvpar4') .AND. (num_procs_p .GT. 1)) THEN
+ ERROR STOP '>> ERROR << dvpar4 velocity dissipation is not compatible with current p parallelization'
+ ENDIF
+
IF(.NOT. KIN_E) THEN
IF(my_id.EQ.0) print*, 'Adiabatic electron model -> beta = 0'
beta = 0._dp
diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90
index 617d52e7..4e2b7de1 100644
--- a/src/moments_eq_rhs_mod.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -24,7 +24,7 @@ SUBROUTINE compute_moments_eq_rhs
xphij_i, xphijp1_i, xphijm1_i, xpsij_i, xpsijp1_i, xpsijm1_i,&
kernel_i, nadiab_moments_i, ddz_nipj, interp_nipj, Sipj,&
moments_i(ipgs_i:ipge_i,ijgs_i:ijge_i,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel),&
- TColl_i, ddzND_nipj, diff_pi_coeff, diff_ji_coeff,&
+ TColl_i, ddzND_Nipj, diff_pi_coeff, diff_ji_coeff,&
moments_rhs_i(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel))
!compute ion moments_eq_rhs
@@ -36,7 +36,7 @@ SUBROUTINE compute_moments_eq_rhs
xphij_e, xphijp1_e, xphijm1_e, xpsij_e, xpsijp1_e, xpsijm1_e,&
kernel_e, nadiab_moments_e, ddz_nepj, interp_nepj, Sepj,&
moments_e(ipgs_e:ipge_e,ijgs_e:ijge_e,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel),&
- TColl_e, ddzND_nepj, diff_pe_coeff, diff_je_coeff,&
+ TColl_e, ddzND_Nepj, diff_pe_coeff, diff_je_coeff,&
moments_rhs_e(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel))
CONTAINS
@@ -193,21 +193,30 @@ SUBROUTINE compute_moments_eq_rhs
-mu_y*diff_ky_coeff*ky**N_HD*moments_(ip,ij,iky,ikx,iz) &
! Numerical parallel hyperdiffusion "mu_z*ddz**4" see Pueschel 2010 (eq 25)
-mu_z*diff_dz_coeff*ddzND_napj_(ip,ij,iky,ikx,iz)
- ! GX like Hermite hypercollisions see Mandell et al. 2023 (eq 3.23), unadvised to use it
- IF (p_int .GT. 2) &
- moments_rhs_(ip,ij,iky,ikx,iz) = &
- moments_rhs_(ip,ij,iky,ikx,iz) - mu_p*diff_pe_coeff*p_int**6*moments_(ip,ij,iky,ikx,iz)
- IF (j_int .GT. 1) &
- moments_rhs_(ip,ij,iky,ikx,iz) = &
- moments_rhs_(ip,ij,iky,ikx,iz) - mu_j*diff_je_coeff*j_int**6*moments_(ip,ij,iky,ikx,iz)
- ! fourth order numerical diffusion in vpar
- ! IF( (ip-4 .GT. 0) .AND. (num_procs_p .EQ. 1) ) &
- ! ! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33)
- ! ! (not used often so not parallelized)
- ! moments_rhs_(ip,ij,iky,ikx,iz) = &
- ! moments_rhs_(ip,ij,iky,ikx,iz) &
- ! + mu_p * moments_(ip-4,ij,iky,ikx,iz)
+ !! Velocity space dissipation (should be implemented somewhere else)
+ SELECT CASE(HYP_V)
+ CASE('hypcoll') ! GX like Hermite hypercollisions see Mandell et al. 2023 (eq 3.23), unadvised to use it
+ IF (p_int .GT. 2) &
+ moments_rhs_(ip,ij,iky,ikx,iz) = &
+ moments_rhs_(ip,ij,iky,ikx,iz) - mu_p*diff_pe_coeff*p_int**6*moments_(ip,ij,iky,ikx,iz)
+ IF (j_int .GT. 1) &
+ moments_rhs_(ip,ij,iky,ikx,iz) = &
+ moments_rhs_(ip,ij,iky,ikx,iz) - mu_j*diff_je_coeff*j_int**6*moments_(ip,ij,iky,ikx,iz)
+ CASE('dvpar4')
+ ! fourth order numerical diffusion in vpar
+ IF(ip-4 .GT. 0) &
+ ! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33)
+ ! (not used often so not parallelized)
+ moments_rhs_(ip,ij,iky,ikx,iz) = &
+ moments_rhs_(ip,ij,iky,ikx,iz) &
+ + mu_p*dv4_Hp_coeff(p_int)*moments_(ip-4,ij,iky,ikx,iz)
+ ! + dummy Laguerre diff
+ IF (j_int .GT. 1) &
+ moments_rhs_(ip,ij,iky,ikx,iz) = &
+ moments_rhs_(ip,ij,iky,ikx,iz) - mu_j*diff_je_coeff*j_int**6*moments_(ip,ij,iky,ikx,iz)
+ CASE DEFAULT
+ END SELECT
ELSE
moments_rhs_(ip,ij,iky,ikx,iz) = 0._dp
ENDIF
@@ -216,7 +225,6 @@ SUBROUTINE compute_moments_eq_rhs
END DO kyloop
END DO kxloop
END DO zloop
-
! Execution time end
CALL cpu_time(t1_rhs)
tc_rhs = tc_rhs + (t1_rhs-t0_rhs)
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index e8589da9..e5ccd88c 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -405,11 +405,11 @@ SUBROUTINE compute_nadiab_moments_z_gradients_and_interp
!
USE fields, ONLY : moments_i, moments_e, phi, psi
USE array, ONLY : kernel_e, kernel_i, nadiab_moments_e, nadiab_moments_i, &
- ddz_nepj, ddzND_nepj, interp_nepj,&
- ddz_nipj, ddzND_nipj, interp_nipj!, ddz_phi
+ ddz_nepj, ddzND_Nepj, interp_nepj,&
+ ddz_nipj, ddzND_Nipj, interp_nipj!, ddz_phi
USE time_integration, ONLY : updatetlevel
USE model, ONLY : qe_taue, qi_taui,q_o_sqrt_tau_sigma_e, q_o_sqrt_tau_sigma_i, &
- KIN_E, CLOS, beta
+ KIN_E, CLOS, beta, HDz_h
USE calculus, ONLY : grad_z, grad_z2, grad_z4, interp_z
IMPLICIT NONE
INTEGER :: eo, p_int, j_int
@@ -489,9 +489,12 @@ SUBROUTINE compute_nadiab_moments_z_gradients_and_interp
! Compute z derivatives
CALL grad_z(eo,nadiab_moments_e(ip,ij,iky,ikx,izgs:izge), ddz_nepj(ip,ij,iky,ikx,izs:ize))
! Parallel hyperdiffusion
- CALL grad_z4(moments_e(ip,ij,iky,ikx,izgs:izge,updatetlevel),ddzND_nepj(ip,ij,iky,ikx,izs:ize))
- ! CALL grad_z2(nadiab_moments_e(ip,ij,iky,ikx,izgs:izge),ddzND_nepj(ip,ij,iky,ikx,izs:ize))
- ! Compute even odd grids interpolation
+ IF (HDz_h) THEN
+ CALL grad_z4(nadiab_moments_e(ip,ij,iky,ikx,izgs:izge),ddzND_Nepj(ip,ij,iky,ikx,izs:ize))
+ ELSE
+ CALL grad_z4(moments_e(ip,ij,iky,ikx,izgs:izge,updatetlevel),ddzND_Nepj(ip,ij,iky,ikx,izs:ize))
+ ENDIF
+ ! Compute even odd grids interpolation
CALL interp_z(eo,nadiab_moments_e(ip,ij,iky,ikx,izgs:izge), interp_nepj(ip,ij,iky,ikx,izs:ize))
ENDDO
ENDDO
@@ -508,8 +511,11 @@ SUBROUTINE compute_nadiab_moments_z_gradients_and_interp
! Compute z first derivative
CALL grad_z(eo,nadiab_moments_i(ip,ij,iky,ikx,izgs:izge), ddz_nipj(ip,ij,iky,ikx,izs:ize))
! Parallel numerical diffusion
- CALL grad_z4(moments_i(ip,ij,iky,ikx,izgs:izge,updatetlevel),ddzND_nipj(ip,ij,iky,ikx,izs:ize))
- ! CALL grad_z2(nadiab_moments_i(ip,ij,iky,ikx,izgs:izge),ddzND_nipj(ip,ij,iky,ikx,izs:ize))
+ IF (HDz_h) THEN
+ CALL grad_z4(nadiab_moments_i(ip,ij,iky,ikx,izgs:izge),ddzND_Nipj(ip,ij,iky,ikx,izs:ize))
+ ELSE
+ CALL grad_z4(moments_i(ip,ij,iky,ikx,izgs:izge,updatetlevel),ddzND_Nipj(ip,ij,iky,ikx,izs:ize))
+ ENDIF
! Compute even odd grids interpolation
CALL interp_z(eo,nadiab_moments_i(ip,ij,iky,ikx,izgs:izge), interp_nipj(ip,ij,iky,ikx,izs:ize))
ENDDO
--
GitLab