diff --git a/src/array_mod.F90 b/src/array_mod.F90
index 5268f0d9d5e4768dc8b985b1ff491755674cba87..7e2e8ddb14fbe5a62c3ea32d4e81346a5aaa1c53 100644
--- a/src/array_mod.F90
+++ b/src/array_mod.F90
@@ -8,10 +8,14 @@ MODULE array
COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: moments_rhs_i ! (ip,ij,ikr,ikz,updatetlevel)
! To load collision matrix
- REAL(dp), DIMENSION(:,:), ALLOCATABLE :: Ceepj, CeipjT
- REAL(dp), DIMENSION(:,:), ALLOCATABLE :: CeipjF
- REAL(dp), DIMENSION(:,:), ALLOCATABLE :: Ciipj, CiepjT
- REAL(dp), DIMENSION(:,:), ALLOCATABLE :: CiepjF
+ REAL(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: Ceepj, CeipjT
+ REAL(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: CeipjF
+ REAL(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: Ciipj, CiepjT
+ REAL(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: CiepjF
+
+ ! Collision term
+ COMPLEX(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: TColl_e, TColl_i
+ COMPLEX(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: TColl_e_local, TColl_i_local
! dnjs coefficient storage (in, ij, is)
COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: dnjs
diff --git a/src/auxval.F90 b/src/auxval.F90
index 56fa89db92885dd11ab3ac73d79a26f6cdc7a6f1..2f57e672c3860f1eb5a2bb0bfcb777163b9d11e9 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -23,6 +23,7 @@ subroutine auxval
CALL set_krgrid ! MPI Distributed dimension
CALL set_kzgrid
+ CALL set_kpgrid
CALL memory ! Allocate memory for global arrays
diff --git a/src/collision_mod.F90 b/src/collision_mod.F90
index 406ba4c846ab94b2d4b1874ad455441c0c3fccf0..26b2be7dd5b6b3d38e24f3ca9503bc57b05a911f 100644
--- a/src/collision_mod.F90
+++ b/src/collision_mod.F90
@@ -1,31 +1,28 @@
module collision
! contains the Hermite-Laguerre collision operators. Solved using COSOlver.
-
-implicit none
-
-PUBLIC :: load_FC_mat
+USE fields
+USE array
+USE basic
+USE grid
+USE prec_const
+USE time_integration
+USE model
+IMPLICIT NONE
+
+PUBLIC :: compute_TColl
PUBLIC :: DoughertyGK_e, DoughertyGK_i
-
+PUBLIC :: load_COSOlver_mat
+PUBLIC :: apply_COSOlver_mat_e, apply_COSOlver_mat_i
CONTAINS
- !******************************************************************************!
- SUBROUTINE LenardBernsteinDK
-
- END SUBROUTINE LenardBernsteinDK
!******************************************************************************!
!! Doughtery gyrokinetic collision operator for electrons
- SUBROUTINE DoughertyGK_e(ip,ij,ikr,ikz,TColl)
- USE fields, ONLY: moments_e, phi
- USE array, ONLY: Kernel_e
- USE basic
- USE grid, ONLY: jmaxe, parray_e, jarray_e, krarray, kzarray
- USE prec_const
- USE time_integration, ONLY: updatetlevel
- USE model, ONLY: sigmae2_taue_o2, tau_e, nu_e, q_e
+ !******************************************************************************!
+ SUBROUTINE DoughertyGK_e(ip_,ij_,ikr_,ikz_,TColl_)
IMPLICIT NONE
- INTEGER, INTENT(IN) :: ip,ij,ikr,ikz
- COMPLEX(dp), INTENT(INOUT) :: TColl
+ INTEGER, INTENT(IN) :: ip_,ij_,ikr_,ikz_
+ COMPLEX(dp), INTENT(INOUT) :: TColl_
COMPLEX(dp) :: n_,upar_,uperp_,Tpar_, Tperp_
COMPLEX(dp) :: Dpj, Ppj, T_, ibe_
@@ -34,20 +31,20 @@ CONTAINS
REAL(dp) :: n_dp, j_dp, p_dp, be_2, q_e_tau_e
!** Auxiliary variables **
- p_dp = REAL(parray_e(ip),dp)
- j_dp = REAL(jarray_e(ij),dp)
- be_2 = (krarray(ikr)**2 + kzarray(ikz)**2) * sigmae2_taue_o2
+ p_dp = REAL(parray_e(ip_),dp)
+ j_dp = REAL(jarray_e(ij_),dp)
+ be_2 = (krarray(ikr_)**2 + kzarray(ikz_)**2) * sigmae2_taue_o2
ibe_ = imagu*2._dp*SQRT(be_2)
q_e_tau_e = q_e/tau_e
!** Assembling collison operator **
! Velocity-space diffusion (similar to Lenhard Bernstein)
- TColl = -(p_dp + 2._dp*j_dp + be_2)*moments_e(ip,ij,ikr,ikz,updatetlevel)
+ TColl_ = -(p_dp + 2._dp*j_dp + be_2)*moments_e(ip_,ij_,ikr_,ikz_,updatetlevel)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
IF( p_dp .EQ. 0 ) THEN ! Kronecker p0
! Get adiabatic moment
- TColl = TColl - (p_dp + 2._dp*j_dp + be_2) * q_e_tau_e * Kernel_e(ij,ikr,ikz)*phi(ikr,ikz)
+ TColl_ = TColl_ - (p_dp + 2._dp*j_dp + be_2) * q_e_tau_e * Kernel_e(ij_,ikr_,ikz_)*phi(ikr_,ikz_)
!** build required fluid moments **
n_ = 0._dp
upar_ = 0._dp; uperp_ = 0._dp
@@ -55,27 +52,27 @@ CONTAINS
DO in_ = 1,jmaxe+1
n_dp = REAL(in_-1,dp)
! Nonadiabatic moments (only different from moments when p=0)
- nadiab_moment_0j = moments_e(1,in_ ,ikr,ikz,updatetlevel) + q_e_tau_e * kernel_e(in_ ,ikr,ikz)*phi(ikr,ikz)
- nadiab_moment_0jp1 = moments_e(1,in_+1,ikr,ikz,updatetlevel) + q_e_tau_e * kernel_e(in_+1,ikr,ikz)*phi(ikr,ikz)
- nadiab_moment_0jm1 = moments_e(1,in_-1,ikr,ikz,updatetlevel) + q_e_tau_e * kernel_e(in_-1,ikr,ikz)*phi(ikr,ikz)
+ nadiab_moment_0j = moments_e(1,in_ ,ikr_,ikz_,updatetlevel) + q_e_tau_e * kernel_e(in_ ,ikr_,ikz_)*phi(ikr_,ikz_)
+ nadiab_moment_0jp1 = moments_e(1,in_+1,ikr_,ikz_,updatetlevel) + q_e_tau_e * kernel_e(in_+1,ikr_,ikz_)*phi(ikr_,ikz_)
+ nadiab_moment_0jm1 = moments_e(1,in_-1,ikr_,ikz_,updatetlevel) + q_e_tau_e * kernel_e(in_-1,ikr_,ikz_)*phi(ikr_,ikz_)
! Density
- n_ = n_ + Kernel_e(in_,ikr,ikz) * nadiab_moment_0j
+ n_ = n_ + Kernel_e(in_,ikr_,ikz_) * nadiab_moment_0j
! Perpendicular velocity
- uperp_ = uperp_ + ibe_*0.5_dp*Kernel_e(in_,ikr,ikz) * (nadiab_moment_0j - nadiab_moment_0jp1)
+ uperp_ = uperp_ + ibe_*0.5_dp*Kernel_e(in_,ikr_,ikz_) * (nadiab_moment_0j - nadiab_moment_0jp1)
! Parallel temperature
- Tpar_ = Tpar_ + Kernel_e(in_,ikr,ikz) * (SQRT2*moments_e(3,in_,ikr,ikz,updatetlevel) + nadiab_moment_0j)
+ Tpar_ = Tpar_ + Kernel_e(in_,ikr_,ikz_) * (SQRT2*moments_e(3,in_,ikr_,ikz_,updatetlevel) + nadiab_moment_0j)
! Perpendicular temperature
- Tperp_ = Tperp_ + Kernel_e(in_,ikr,ikz) * ((2._dp*n_dp+1._dp)* nadiab_moment_0j &
+ Tperp_ = Tperp_ + Kernel_e(in_,ikr_,ikz_) * ((2._dp*n_dp+1._dp)* nadiab_moment_0j &
- n_dp * nadiab_moment_0jm1 &
- (n_dp+1)* nadiab_moment_0jp1)
ENDDO
T_ = (Tpar_ + 2._dp*Tperp_)/3._dp - n_
! Add energy restoring term
- TColl = TColl + T_* 4._dp * j_dp * Kernel_e(ij ,ikr,ikz)
- TColl = TColl - T_* 2._dp * (j_dp + 1._dp) * Kernel_e(ij+1,ikr,ikz)
- TColl = TColl - T_* 2._dp * j_dp * Kernel_e(ij-1,ikr,ikz)
- TColl = TColl + uperp_*ibe_*((j_dp + 1._dp)* Kernel_e(ij ,ikr,ikz) &
- -j_dp * Kernel_e(ij-1,ikr,ikz))
+ TColl_ = TColl_ + T_* 4._dp * j_dp * Kernel_e(ij_ ,ikr_,ikz_)
+ TColl_ = TColl_ - T_* 2._dp * (j_dp + 1._dp) * Kernel_e(ij_+1,ikr_,ikz_)
+ TColl_ = TColl_ - T_* 2._dp * j_dp * Kernel_e(ij_-1,ikr_,ikz_)
+ TColl_ = TColl_ + uperp_*ibe_*((j_dp + 1._dp)* Kernel_e(ij_ ,ikr_,ikz_) &
+ -j_dp * Kernel_e(ij_-1,ikr_,ikz_))
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -84,9 +81,9 @@ CONTAINS
upar_ = 0._dp
DO in_ = 1,jmaxe+1
! Parallel velocity
- upar_ = upar_ + Kernel_e(in_,ikr,ikz) * moments_e(2,in_,ikr,ikz,updatetlevel)
+ upar_ = upar_ + Kernel_e(in_,ikr_,ikz_) * moments_e(2,in_,ikr_,ikz_,updatetlevel)
ENDDO
- TColl = TColl + upar_*Kernel_e(ij,ikr,ikz)
+ TColl_ = TColl_ + upar_*Kernel_e(ij_,ikr_,ikz_)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -97,39 +94,32 @@ CONTAINS
DO in_ = 1,jmaxe+1
n_dp = REAL(in_-1,dp)
! Nonadiabatic moments (only different from moments when p=0)
- nadiab_moment_0j = moments_e(1,in_ ,ikr,ikz,updatetlevel) + q_e_tau_e * kernel_e(in_ ,ikr,ikz)*phi(ikr,ikz)/tau_e
- nadiab_moment_0jp1 = moments_e(1,in_+1,ikr,ikz,updatetlevel) + q_e_tau_e * kernel_e(in_+1,ikr,ikz)*phi(ikr,ikz)/tau_e
- nadiab_moment_0jm1 = moments_e(1,in_-1,ikr,ikz,updatetlevel) + q_e_tau_e * kernel_e(in_-1,ikr,ikz)*phi(ikr,ikz)/tau_e
+ nadiab_moment_0j = moments_e(1,in_ ,ikr_,ikz_,updatetlevel) + q_e_tau_e * kernel_e(in_ ,ikr_,ikz_)*phi(ikr_,ikz_)/tau_e
+ nadiab_moment_0jp1 = moments_e(1,in_+1,ikr_,ikz_,updatetlevel) + q_e_tau_e * kernel_e(in_+1,ikr_,ikz_)*phi(ikr_,ikz_)/tau_e
+ nadiab_moment_0jm1 = moments_e(1,in_-1,ikr_,ikz_,updatetlevel) + q_e_tau_e * kernel_e(in_-1,ikr_,ikz_)*phi(ikr_,ikz_)/tau_e
! Density
- n_ = n_ + Kernel_e(in_,ikr,ikz) * nadiab_moment_0j
+ n_ = n_ + Kernel_e(in_,ikr_,ikz_) * nadiab_moment_0j
! Parallel temperature
- Tpar_ = Tpar_ + Kernel_e(in_,ikr,ikz) * (SQRT2*moments_e(3,in_,ikr,ikz,updatetlevel) + nadiab_moment_0j)
+ Tpar_ = Tpar_ + Kernel_e(in_,ikr_,ikz_) * (SQRT2*moments_e(3,in_,ikr_,ikz_,updatetlevel) + nadiab_moment_0j)
! Perpendicular temperature
- Tperp_ = Tperp_ + Kernel_e(in_,ikr,ikz) * ((2._dp*n_dp+1._dp)*nadiab_moment_0j - n_dp*nadiab_moment_0jm1 - (n_dp+1)*nadiab_moment_0jp1)
+ Tperp_ = Tperp_ + Kernel_e(in_,ikr_,ikz_) * ((2._dp*n_dp+1._dp)*nadiab_moment_0j - n_dp*nadiab_moment_0jm1 - (n_dp+1)*nadiab_moment_0jp1)
ENDDO
T_ = (Tpar_ + 2._dp*Tperp_)/3._dp - n_
- TColl = TColl + T_*SQRT2*Kernel_e(ij,ikr,ikz)
+ TColl_ = TColl_ + T_*SQRT2*Kernel_e(ij_,ikr_,ikz_)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
ENDIF
! Multiply by electron-ion collision coefficient
- TColl = nu_e * TColl
+ TColl_ = nu_e * TColl_
END SUBROUTINE DoughertyGK_e
!******************************************************************************!
!! Doughtery gyrokinetic collision operator for electrons
- SUBROUTINE DoughertyGK_i(ip,ij,ikr,ikz,TColl)
- USE fields, ONLY: moments_i, phi
- USE array, ONLY: Kernel_i
- USE basic
- USE grid, ONLY: jmaxi, parray_i, jarray_i, krarray, kzarray
- USE prec_const
- USE time_integration, ONLY: updatetlevel
- USE model, ONLY: sigmai2_taui_o2, tau_i, nu_i, q_i
+ SUBROUTINE DoughertyGK_i(ip_,ij_,ikr_,ikz_,TColl_)
IMPLICIT NONE
- INTEGER, INTENT(IN) :: ip,ij,ikr,ikz
- COMPLEX(dp), INTENT(INOUT) :: TColl
+ INTEGER, INTENT(IN) :: ip_,ij_,ikr_,ikz_
+ COMPLEX(dp), INTENT(INOUT) :: TColl_
COMPLEX(dp) :: n_,upar_,uperp_,Tpar_, Tperp_
COMPLEX(dp) :: Dpj, Ppj, T_, ibi_
@@ -138,20 +128,20 @@ CONTAINS
REAL(dp) :: n_dp, j_dp, p_dp, bi_2, q_i_tau_i
!** Auxiliary variables **
- p_dp = REAL(parray_i(ip),dp)
- j_dp = REAL(jarray_i(ij),dp)
- bi_2 = (krarray(ikr)**2 + kzarray(ikz)**2) * sigmai2_taui_o2
+ p_dp = REAL(parray_i(ip_),dp)
+ j_dp = REAL(jarray_i(ij_),dp)
+ bi_2 = (krarray(ikr_)**2 + kzarray(ikz_)**2) * sigmai2_taui_o2
ibi_ = imagu*2._dp*SQRT(bi_2)
q_i_tau_i = q_i/tau_i
!** Assembling collison operator **
! Velocity-space diffusion (similar to Lenhard Bernstein)
- TColl = -(p_dp + 2._dp*j_dp + bi_2)*moments_i(ip,ij,ikr,ikz,updatetlevel)
+ TColl_ = -(p_dp + 2._dp*j_dp + bi_2)*moments_i(ip_,ij_,ikr_,ikz_,updatetlevel)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
IF( p_dp .EQ. 0 ) THEN ! Kronecker p0
! Get adiabatic moment
- TColl = TColl - (p_dp + 2._dp*j_dp + bi_2) * q_i_tau_i * Kernel_i(ij,ikr,ikz)*phi(ikr,ikz)
+ TColl_ = TColl_ - (p_dp + 2._dp*j_dp + bi_2) * q_i_tau_i * Kernel_i(ij_,ikr_,ikz_)*phi(ikr_,ikz_)
!** build required fluid moments **
n_ = 0._dp
upar_ = 0._dp; uperp_ = 0._dp
@@ -159,27 +149,27 @@ CONTAINS
DO in_ = 1,jmaxi+1
n_dp = REAL(in_-1,dp)
! Nonadiabatic moments (only different from moments when p=0)
- nadiab_moment_0j = moments_i(1,in_ ,ikr,ikz,updatetlevel) + q_i_tau_i * kernel_i(in_ ,ikr,ikz)*phi(ikr,ikz)
- nadiab_moment_0jp1 = moments_i(1,in_+1,ikr,ikz,updatetlevel) + q_i_tau_i * kernel_i(in_+1,ikr,ikz)*phi(ikr,ikz)
- nadiab_moment_0jm1 = moments_i(1,in_-1,ikr,ikz,updatetlevel) + q_i_tau_i * kernel_i(in_-1,ikr,ikz)*phi(ikr,ikz)
+ nadiab_moment_0j = moments_i(1,in_ ,ikr_,ikz_,updatetlevel) + q_i_tau_i * kernel_i(in_ ,ikr_,ikz_)*phi(ikr_,ikz_)
+ nadiab_moment_0jp1 = moments_i(1,in_+1,ikr_,ikz_,updatetlevel) + q_i_tau_i * kernel_i(in_+1,ikr_,ikz_)*phi(ikr_,ikz_)
+ nadiab_moment_0jm1 = moments_i(1,in_-1,ikr_,ikz_,updatetlevel) + q_i_tau_i * kernel_i(in_-1,ikr_,ikz_)*phi(ikr_,ikz_)
! Density
- n_ = n_ + Kernel_i(in_,ikr,ikz) * nadiab_moment_0j
+ n_ = n_ + Kernel_i(in_,ikr_,ikz_) * nadiab_moment_0j
! Perpendicular velocity
- uperp_ = uperp_ + ibi_*0.5_dp*Kernel_i(in_,ikr,ikz) * (nadiab_moment_0j -nadiab_moment_0jp1)
+ uperp_ = uperp_ + ibi_*0.5_dp*Kernel_i(in_,ikr_,ikz_) * (nadiab_moment_0j -nadiab_moment_0jp1)
! Parallel temperature
- Tpar_ = Tpar_ + Kernel_i(in_,ikr,ikz) * (SQRT2*moments_i(3,in_,ikr,ikz,updatetlevel) + nadiab_moment_0j)
+ Tpar_ = Tpar_ + Kernel_i(in_,ikr_,ikz_) * (SQRT2*moments_i(3,in_,ikr_,ikz_,updatetlevel) + nadiab_moment_0j)
! Perpendicular temperature
- Tperp_ = Tperp_ + Kernel_i(in_,ikr,ikz) * ((2._dp*n_dp+1._dp)* nadiab_moment_0j &
+ Tperp_ = Tperp_ + Kernel_i(in_,ikr_,ikz_) * ((2._dp*n_dp+1._dp)* nadiab_moment_0j &
- n_dp * nadiab_moment_0jm1 &
- (n_dp+1)* nadiab_moment_0jp1)
ENDDO
T_ = (Tpar_ + 2._dp*Tperp_)/3._dp - n_
! Add energy restoring term
- TColl = TColl + T_* 4._dp * j_dp * Kernel_i(ij ,ikr,ikz)
- TColl = TColl - T_* 2._dp * (j_dp + 1._dp) * Kernel_i(ij+1,ikr,ikz)
- TColl = TColl - T_* 2._dp * j_dp * Kernel_i(ij-1,ikr,ikz)
- TColl = TColl + uperp_*ibi_*((j_dp + 1._dp)* Kernel_i(ij ,ikr,ikz) &
- -j_dp * Kernel_i(ij-1,ikr,ikz))
+ TColl_ = TColl_ + T_* 4._dp * j_dp * Kernel_i(ij_ ,ikr_,ikz_)
+ TColl_ = TColl_ - T_* 2._dp * (j_dp + 1._dp) * Kernel_i(ij_+1,ikr_,ikz_)
+ TColl_ = TColl_ - T_* 2._dp * j_dp * Kernel_i(ij_-1,ikr_,ikz_)
+ TColl_ = TColl_ + uperp_*ibi_*((j_dp + 1._dp)* Kernel_i(ij_ ,ikr_,ikz_) &
+ -j_dp * Kernel_i(ij_-1,ikr_,ikz_))
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -188,9 +178,9 @@ CONTAINS
upar_ = 0._dp
DO in_ = 1,jmaxi+1
! Parallel velocity
- upar_ = upar_ + Kernel_i(in_,ikr,ikz) * moments_i(2,in_,ikr,ikz,updatetlevel)
+ upar_ = upar_ + Kernel_i(in_,ikr_,ikz_) * moments_i(2,in_,ikr_,ikz_,updatetlevel)
ENDDO
- TColl = TColl + upar_*Kernel_i(ij,ikr,ikz)
+ TColl_ = TColl_ + upar_*Kernel_i(ij_,ikr_,ikz_)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -201,128 +191,193 @@ CONTAINS
DO in_ = 1,jmaxi+1
n_dp = REAL(in_-1,dp)
! Nonadiabatic moments (only different from moments when p=0)
- nadiab_moment_0j = moments_i(1,in_ ,ikr,ikz,updatetlevel) + q_i_tau_i * kernel_i(in_ ,ikr,ikz)*phi(ikr,ikz)
- nadiab_moment_0jp1 = moments_i(1,in_+1,ikr,ikz,updatetlevel) + q_i_tau_i * kernel_i(in_+1,ikr,ikz)*phi(ikr,ikz)
- nadiab_moment_0jm1 = moments_i(1,in_-1,ikr,ikz,updatetlevel) + q_i_tau_i * kernel_i(in_-1,ikr,ikz)*phi(ikr,ikz)
+ nadiab_moment_0j = moments_i(1,in_ ,ikr_,ikz_,updatetlevel) + q_i_tau_i * kernel_i(in_ ,ikr_,ikz_)*phi(ikr_,ikz_)
+ nadiab_moment_0jp1 = moments_i(1,in_+1,ikr_,ikz_,updatetlevel) + q_i_tau_i * kernel_i(in_+1,ikr_,ikz_)*phi(ikr_,ikz_)
+ nadiab_moment_0jm1 = moments_i(1,in_-1,ikr_,ikz_,updatetlevel) + q_i_tau_i * kernel_i(in_-1,ikr_,ikz_)*phi(ikr_,ikz_)
! Density
- n_ = n_ + Kernel_i(in_,ikr,ikz) * nadiab_moment_0j
+ n_ = n_ + Kernel_i(in_,ikr_,ikz_) * nadiab_moment_0j
! Parallel temperature
- Tpar_ = Tpar_ + Kernel_i(in_,ikr,ikz) * (SQRT2*moments_i(3,in_,ikr,ikz,updatetlevel) + nadiab_moment_0j)
+ Tpar_ = Tpar_ + Kernel_i(in_,ikr_,ikz_) * (SQRT2*moments_i(3,in_,ikr_,ikz_,updatetlevel) + nadiab_moment_0j)
! Perpendicular temperature
- Tperp_ = Tperp_ + Kernel_i(in_,ikr,ikz) * ((2._dp*n_dp+1._dp)*nadiab_moment_0j - n_dp*nadiab_moment_0jm1 - (n_dp+1)*nadiab_moment_0jp1)
+ Tperp_ = Tperp_ + Kernel_i(in_,ikr_,ikz_) * ((2._dp*n_dp+1._dp)*nadiab_moment_0j - n_dp*nadiab_moment_0jm1 - (n_dp+1)*nadiab_moment_0jp1)
ENDDO
T_ = (Tpar_ + 2._dp*Tperp_)/3._dp - n_
- TColl = TColl + T_*SQRT2*Kernel_i(ij,ikr,ikz)
+ TColl_ = TColl_ + T_*SQRT2*Kernel_i(ij_,ikr_,ikz_)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
ENDIF
! Multiply by ion-ion collision coefficient
- TColl = nu_i * TColl
+ TColl_ = nu_i * TColl_
END SUBROUTINE DoughertyGK_i
!******************************************************************************!
- !!!!!!! Compute ion Full Coulomb collision operator
+ !! compute the collision terms in a (Np x Nj x Nkr x Nkz) matrix all at once
+ !******************************************************************************!
+ SUBROUTINE compute_TColl
+ IMPLICIT NONE
+ COMPLEX(dp), DIMENSION(1:pmaxe+1) :: local_sum_e, buffer_e, total_sum_e
+ COMPLEX(dp), DIMENSION(ips_e:ipe_e) :: TColl_distr_e
+ COMPLEX(dp), DIMENSION(1:pmaxi+1) :: local_sum_i, buffer_i, total_sum_i
+ COMPLEX(dp), DIMENSION(ips_i:ipe_i) :: TColl_distr_i
+ COMPLEX(dp) :: TColl
+ INTEGER :: ikrs_C, ikre_C, ikzs_C, ikze_C
+ IF (ABS(CO) .GE. 2) THEN !compute only if COSOlver matrices are used
+
+ DO ikr = ikrs,ikre
+ DO ikz = ikzs,ikze
+ ! Electrons
+ DO ij = 1,Jmaxe+1
+ ! Loop over all p to compute sub collision term
+ DO ip = 1,Pmaxe+1
+ CALL apply_COSOlver_mat_e(ip,ij,ikr,ikz,TColl)
+ local_sum_e(ip) = TColl
+ ENDDO
+ ! Sum up all the sub collision terms on root 0
+ CALL MPI_REDUCE(local_sum_e, buffer_e, pmaxe+1, MPI_DOUBLE_COMPLEX, MPI_SUM, 0, comm_p, ierr)
+ ! IF(rank_p .EQ. 0) THEN
+ ! DO ip = 1,Pmaxe+1
+ ! total_sum_e(ip) = buffer_e(ip)
+ ! ENDDO
+ ! ENDIF
+ ! distribute the sum over the process among p
+ CALL MPI_SCATTERV(buffer_e, counts_np_e, displs_np_e, MPI_DOUBLE_COMPLEX,&
+ TColl_distr_e, local_np_e, MPI_DOUBLE_COMPLEX,&
+ 0, comm_p, ierr)
+ DO ip = ips_e,ipe_e
+ TColl_e(ip,ij,ikr,ikz) = TColl_distr_e(ip)
+ ENDDO
+ ENDDO
+ ! Ions
+ DO ij = 1,Jmaxi+1
+ DO ip = 1,Pmaxi+1
+ CALL apply_COSOlver_mat_i(ip,ij,ikr,ikz,TColl)
+ local_sum_i(ip) = TColl
+ ENDDO
+ ! Reduce the local_sums to root = 0
+ CALL MPI_REDUCE(local_sum_i, buffer_i, pmaxi+1, MPI_DOUBLE_COMPLEX, MPI_SUM, 0, comm_p, ierr)
+ ! IF(rank_p .EQ. 0) THEN
+ ! DO ip = 1,Pmaxi+1
+ ! total_sum_i(ip) = buffer_i(ip)
+ ! ENDDO
+ ! ENDIF
+ ! buffer_e contains the entire collision term along p, scatter it among
+ ! the other processes (use of scatterv since Pmax/Np is not an integer)
+ CALL MPI_SCATTERV(buffer_i, counts_np_i, displs_np_i, MPI_DOUBLE_COMPLEX,&
+ TColl_distr_i, local_np_i, MPI_DOUBLE_COMPLEX, &
+ 0, comm_p, ierr)
+ DO ip = ips_i,ipe_i
+ TColl_i(ip,ij,ikr,ikz) = TColl_distr_i(ip)
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+ ! IF(cstep .EQ. 1) THEN
+ ! write(*,*) rank_p, ': local_sum = ', local_sum_e
+ ! write(*,*) rank_p, ': buffer = ', buffer_e
+ ! write(*,*) rank_p, ': dist_sum = ', TColl_distr_e
+ ! ENDIF
+
+ END SUBROUTINE compute_TColl
+
+ !******************************************************************************!
+ !!!!!!! Compute ion collision term
!******************************************************************************!
- SUBROUTINE FullCoulombDK_e(p_int,j_int,ikr,ikz,TColl)
- USE basic
- USE fields, ONLY: moments_i, moments_e
- USE array, ONLY: Ceepj, CeipjT, CeipjF
- USE grid, ONLY: pmaxe,jmaxe, pmaxi,jmaxi, parray_e, parray_i, &
- jarray_e, jarray_i, bari, bare
- USE prec_const
- USE time_integration, ONLY: updatetlevel
- USE model, ONLY: nu_e, nu_ee
+ SUBROUTINE apply_COSOlver_mat_e(ip_,ij_,ikr_,ikz_,TColl_)
IMPLICIT NONE
- INTEGER, INTENT(IN) :: p_int,j_int ,ikr,ikz
- COMPLEX(dp), INTENT(INOUT) :: TColl
+ INTEGER, INTENT(IN) :: ip_, ij_ ,ikr_, ikz_
+ COMPLEX(dp), INTENT(OUT) :: TColl_
- INTEGER :: ip2,ij2, p2_int,j2_int
+ INTEGER :: ip2,ij2, p_int,j_int, p2_int,j2_int, ikr_C, ikz_C
+ p_int = ip_-1; j_int = ij_-1
- TColl = 0._dp ! Initialization
+ IF (CO .GT. 0) THEN ! GK operator (k-dependant)
+ ikr_C = ikr_; ikz_C = ikz_
+ ELSEIF (CO .LT. 0) THEN ! DK operator (only one mat for every k)
+ ikr_C = 1; ikz_C = 1
+ ENDIF
+
+ TColl_ = 0._dp ! Initialization of the local sum
- ploopee: DO ip2 = 1,pmaxe+1 ! sum the electron-self and electron-ion test terms
+ ! sum the electron-self and electron-ion test terms
+ ploopee: DO ip2 = ips_e,ipe_e
p2_int = parray_e(ip2)
- jloopee: DO ij2 = 1,jmaxe+1
+ jloopee: DO ij2 = ijs_e,ije_e
j2_int = jarray_e(ij2)
- TColl = TColl + moments_e(ip2,ij2,ikr,ikz,updatetlevel) &
- *( nu_e * CeipjT(bare(p_int,j_int), bare(p2_int,j2_int)) &
- +nu_ee * Ceepj (bare(p_int,j_int), bare(p2_int,j2_int)))
+ TColl_ = TColl_ + moments_e(ip2,ij2,ikr_,ikz_,updatetlevel) &
+ *( nu_e * CeipjT(bare(p_int,j_int), bare(p2_int,j2_int),ikr_C, ikz_C) &
+ +nu_ee * Ceepj (bare(p_int,j_int), bare(p2_int,j2_int),ikr_C, ikz_C))
ENDDO jloopee
ENDDO ploopee
- ploopei: DO ip2 = 1,pmaxi+1 ! sum the electron-ion field terms
+
+ ! sum the electron-ion field terms
+ ploopei: DO ip2 = ips_i,ipe_i
p2_int = parray_i(ip2)
- jloopei: DO ij2 = 1,jmaxi+1
+ jloopei: DO ij2 = ijs_i,ije_i
j2_int = jarray_i(ij2)
- TColl = TColl + moments_i(ip2,ij2,ikr,ikz,updatetlevel) &
- *(nu_e * CeipjF(bare(p_int,j_int), bari(p2_int,j2_int)))
+ TColl_ = TColl_ + moments_i(ip2,ij2,ikr_,ikz_,updatetlevel) &
+ *(nu_e * CeipjF(bare(p_int,j_int), bari(p2_int,j2_int),ikr_C, ikz_C))
END DO jloopei
ENDDO ploopei
- END SUBROUTINE FullCoulombDK_e
+ END SUBROUTINE apply_COSOlver_mat_e
!******************************************************************************!
- !!!!!!! Compute ion Full Coulomb collision operator
+ !!!!!!! Compute ion collision term
!******************************************************************************!
- SUBROUTINE FullCoulombDK_i(p_int,j_int,ikr,ikz,TColl)
- USE basic
- USE fields, ONLY: moments_i, moments_e
- USE array, ONLY: Ciipj, CiepjT, CiepjF
- USE grid, ONLY: pmaxe,jmaxe, pmaxi,jmaxi, parray_e, parray_i, &
- jarray_e, jarray_i, bari, bare
- USE prec_const
- USE time_integration, ONLY: updatetlevel
- USE model, ONLY: nu_i, nu_ie
+ SUBROUTINE apply_COSOlver_mat_i(ip_,ij_,ikr_,ikz_,TColl_)
IMPLICIT NONE
- INTEGER, INTENT(IN) :: p_int,j_int,ikr,ikz
- COMPLEX(dp), INTENT(INOUT) :: TColl
+ INTEGER, INTENT(IN) :: ip_, ij_ ,ikr_, ikz_
+ COMPLEX(dp), INTENT(OUT) :: TColl_
- INTEGER :: ip2,ij2, p2_int,j2_int
+ INTEGER :: ip2,ij2, p_int,j_int, p2_int,j2_int, ikr_C, ikz_C
+ p_int = ip_-1; j_int = ij_-1
- TColl = 0._dp ! Initialization
+ IF (CO .GT. 0) THEN ! GK operator (k-dependant)
+ ikr_C = ikr_; ikz_C = ikz_
+ ELSEIF (CO .LT. 0) THEN ! DK operator (only one mat for every k)
+ ikr_C = 1; ikz_C = 1
+ ENDIF
- ploopii: DO ip2 = 1,pmaxi+1 ! sum the ion-self and ion-electron test terms
+ TColl_ = 0._dp ! Initialization
+ ! sum the ion-self and ion-electron test terms
+ ploopii: DO ip2 = ips_i,ipe_i
p2_int = parray_i(ip2)
- jloopii: DO ij2 = 1,jmaxi+1
+ jloopii: DO ij2 = ijs_i,ije_i
j2_int = jarray_i(ij2)
- TColl = TColl + moments_i(ip2,ij2,ikr,ikz,updatetlevel) &
- *( nu_ie * CiepjT(bari(p_int,j_int), bari(p2_int,j2_int)) &
- +nu_i * Ciipj (bari(p_int,j_int), bari(p2_int,j2_int)))
+ TColl_ = TColl_ + moments_i(ip2,ij2,ikr_,ikz_,updatetlevel) &
+ *( nu_ie * CiepjT(bari(p_int,j_int), bari(p2_int,j2_int), ikr_C, ikz_C) &
+ +nu_i * Ciipj (bari(p_int,j_int), bari(p2_int,j2_int), ikr_C, ikz_C))
ENDDO jloopii
ENDDO ploopii
- ploopie: DO ip2 = 1,pmaxe+1 ! sum the ion-electron field terms
+ ploopie: DO ip2 = ips_e,ipe_e ! sum the ion-electron field terms
p2_int = parray_e(ip2)
- jloopie: DO ij2 = 1,jmaxe+1
+ jloopie: DO ij2 = ijs_e,ije_e
j2_int = jarray_e(ij2)
- TColl = TColl + moments_e(ip2,ij2,ikr,ikz,updatetlevel) &
- *(nu_ie * CiepjF(bari(p_int,j_int), bare(p2_int,j2_int)))
+ TColl_ = TColl_ + moments_e(ip2,ij2,ikr_,ikz_,updatetlevel) &
+ *(nu_ie * CiepjF(bari(p_int,j_int), bare(p2_int,j2_int), ikr_C, ikz_C))
ENDDO jloopie
ENDDO ploopie
- END SUBROUTINE FullCoulombDK_i
+ END SUBROUTINE apply_COSOlver_mat_i
!******************************************************************************!
- !!!!!!! Load the Full coulomb coefficient table from COSOlver results
+ !!!!!!! Load the collision matrix coefficient table from COSOlver results
!******************************************************************************!
- SUBROUTINE load_FC_mat ! Load a sub matrix from iCa files (works for pmaxa,jmaxa<=P_full,J_full)
- use prec_const
- use fields
- use array
- use grid
- use basic
+ SUBROUTINE load_COSOlver_mat ! Load a sub matrix from iCa files (works for pmaxa,jmaxa<=P_full,J_full)
use futils
use initial_par
- use model
IMPLICIT NONE
INTEGER :: irow_sub, irow_full, icol_sub, icol_full
INTEGER :: fid1, fid2, fid3, fid4
- INTEGER :: ip_e, ij_e, il_e, ik_e
+ INTEGER :: ip_e, ij_e, il_e, ik_e, ikrs_C, ikre_C, ikzs_C, ikze_C
INTEGER :: pdime, jdime
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: Ceepj_full, CeipjT_full
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: CeipjF_full
@@ -331,141 +386,190 @@ CONTAINS
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: Ciipj_full, CiepjT_full
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: CiepjF_full
- !!!!!!!!!!!! Electron matrices !!!!!!!!!!!!
- ! get the self electron colision matrix
- CALL openf(selfmat_file,fid1, 'r', 'D',mpicomm=MPI_COMM_WORLD)
-
- CALL getatt(fid1,'/Caapj/Ceepj/','Pmaxe',pdime)
- CALL getatt(fid1,'/Caapj/Ceepj/','Jmaxe',jdime)
-
- IF ( ((pdime .LT. pmaxe) .OR. (jdime .LT. jmaxe)) .AND. (my_id .EQ. 0)) WRITE(*,*) '!! FC Matrix too small !!'
-
- CALL allocate_array( Ceepj_full, 1,(pdime+1)*(jdime+1), 1,(pdime+1)*(jdime+1))
- CALL getarr(fid1, '/Caapj/Ceepj', Ceepj_full) ! get array (moli format)
-
- ! Fill sub array with only usefull polynmials degree
- DO ip_e = 0,pmaxe ! Loop over rows
- DO ij_e = 0,jmaxe
- irow_sub = (jmaxe +1)*ip_e + ij_e +1
- irow_full = (jdime +1)*ip_e + ij_e +1
- DO il_e = 0,pmaxe ! Loop over columns
- DO ik_e = 0,jmaxe
- icol_sub = (jmaxe +1)*il_e + ik_e +1
- icol_full = (jdime +1)*il_e + ik_e +1
- Ceepj (irow_sub,icol_sub) = Ceepj_full (irow_full,icol_full)
- ENDDO
- ENDDO
- ENDDO
- ENDDO
-
- CALL closef(fid1)
- DEALLOCATE(Ceepj_full)
-
- ! get the Test and Back field electron ion collision matrix
- CALL openf(eimat_file,fid2, 'r', 'D');
-
- CALL getatt(fid2,'/Ceipj/CeipjT/','Pmaxi',pdimi)
- CALL getatt(fid2,'/Ceipj/CeipjT/','Jmaxi',jdimi)
- IF ( (pdimi .LT. pmaxi) .OR. (jdimi .LT. jmaxi) ) WRITE(*,*) '!! FC Matrix too small !!'
-
- CALL allocate_array( CeipjT_full, 1,(pdime+1)*(jdime+1), 1,(pdime+1)*(jdime+1))
- CALL allocate_array( CeipjF_full, 1,(pdime+1)*(jdime+1), 1,(pdimi+1)*(jdimi+1))
-
- CALL getarr(fid2, '/Ceipj/CeipjT', CeipjT_full)
- CALL getarr(fid2, '/Ceipj/CeipjF', CeipjF_full)
-
- ! Fill sub array with only usefull polynmials degree
- DO ip_e = 0,pmaxe ! Loop over rows
- DO ij_e = 0,jmaxe
- irow_sub = (jmaxe +1)*ip_e + ij_e +1
- irow_full = (jdime +1)*ip_e + ij_e +1
- DO il_e = 0,pmaxe ! Loop over columns
- DO ik_e = 0,jmaxe
- icol_sub = (jmaxe +1)*il_e + ik_e +1
- icol_full = (jdime +1)*il_e + ik_e +1
- CeipjT(irow_sub,icol_sub) = CeipjT_full(irow_full,icol_full)
- ENDDO
- ENDDO
- DO il_i = 0,pmaxi ! Loop over columns
- DO ik_i = 0,jmaxi
- icol_sub = (Jmaxi +1)*il_i + ik_i +1
- icol_full = (jdimi +1)*il_i + ik_i +1
- CeipjF(irow_sub,icol_sub) = CeipjF_full(irow_full,icol_full)
- ENDDO
- ENDDO
- ENDDO
- ENDDO
-
- CALL closef(fid2)
- DEALLOCATE(CeipjF_full)
- DEALLOCATE(CeipjT_full)
-
- !!!!!!!!!!!!!!! Ion matrices !!!!!!!!!!!!!!
- ! get the self electron colision matrix
- CALL openf(selfmat_file, fid3, 'r', 'D',mpicomm=MPI_COMM_WORLD);
-
- CALL allocate_array( Ciipj_full, 1,(pdimi+1)*(jdimi+1), 1,(pdimi+1)*(jdimi+1))
-
- IF ( (pmaxe .EQ. pmaxi) .AND. (jmaxe .EQ. jmaxi) ) THEN ! if same degrees ion and electron matrices are alike so load Ceepj
- CALL getarr(fid3, '/Caapj/Ceepj', Ciipj_full) ! get array (moli format)
- ELSE
- CALL getarr(fid3, '/Caapj/Ciipj', Ciipj_full) ! get array (moli format)
+ CHARACTER(len=256) :: mat_filename, kperp_string
+
+ !! Some terminal info
+ IF (CO .EQ. 2) THEN
+ IF (my_id .EQ. 0) WRITE(*,*) '=== Load GK Sugama matrix ==='
+ ELSEIF(CO .EQ. 3) THEN
+ IF (my_id .EQ. 0) WRITE(*,*) '=== Load GK Full Coulomb matrix ==='
+ ELSEIF(CO .EQ. -2) THEN
+ IF (my_id .EQ. 0) WRITE(*,*) '=== Load DK Sugama matrix ==='
+ ELSEIF(CO .EQ. -3) THEN
+ IF (my_id .EQ. 0) WRITE(*,*) '=== Load DK Full Coulomb matrix ==='
ENDIF
- ! Fill sub array with only usefull polynmials degree
- DO ip_i = 0,Pmaxi ! Loop over rows
- DO ij_i = 0,Jmaxi
- irow_sub = (Jmaxi +1)*ip_i + ij_i +1
- irow_full = (jdimi +1)*ip_i + ij_i +1
- DO il_i = 0,Pmaxi ! Loop over columns
- DO ik_i = 0,Jmaxi
- icol_sub = (Jmaxi +1)*il_i + ik_i +1
- icol_full = (jdimi +1)*il_i + ik_i +1
- Ciipj (irow_sub,icol_sub) = Ciipj_full (irow_full,icol_full)
- ENDDO
- ENDDO
- ENDDO
- ENDDO
-
- CALL closef(fid3)
- DEALLOCATE(Ciipj_full)
-
- ! get the Test and Back field electron ion collision matrix
- CALL openf(iemat_file,fid4, 'r', 'D');
-
- CALL allocate_array( CiepjT_full, 1,(pdimi+1)*(jdimi+1), 1,(pdimi+1)*(jdimi+1))
- CALL allocate_array( CiepjF_full, 1,(pdimi+1)*(jdimi+1), 1,(pdime+1)*(jdime+1))
-
- CALL getarr(fid4, '/Ciepj/CiepjT', CiepjT_full)
- CALL getarr(fid4, '/Ciepj/CiepjF', CiepjF_full)
+ IF (CO .GT. 0) THEN ! GK operator (k-dependant)
+ ikrs_C = ikrs; ikre_C = ikre
+ ikzs_C = ikzs; ikze_C = ikze
+ ELSEIF (CO .LT. 0) THEN ! DK operator (only one mat for every k)
+ ikrs_C = 1; ikre_C = 1
+ ikzs_C = 1; ikze_C = 1
+ ENDIF
- ! Fill sub array with only usefull polynmials degree
- DO ip_i = 0,Pmaxi ! Loop over rows
- DO ij_i = 0,Jmaxi
- irow_sub = (Jmaxi +1)*ip_i + ij_i +1
- irow_full = (jdimi +1)*ip_i + ij_i +1
- DO il_i = 0,Pmaxi ! Loop over columns
- DO ik_i = 0,Jmaxi
- icol_sub = (Jmaxi +1)*il_i + ik_i +1
- icol_full = (jdimi +1)*il_i + ik_i +1
- CiepjT(irow_sub,icol_sub) = CiepjT_full(irow_full,icol_full)
- ENDDO
- ENDDO
- DO il_e = 0,pmaxe ! Loop over columns
- DO ik_e = 0,jmaxe
- icol_sub = (jmaxe +1)*il_e + ik_e +1
- icol_full = (jdime +1)*il_e + ik_e +1
- CiepjF(irow_sub,icol_sub) = CiepjF_full(irow_full,icol_full)
- ENDDO
- ENDDO
+ DO ikr = ikrs_C,ikre_C ! Loop over everz kperp values
+ DO ikz = ikzs_C,ikze_C ! Loop over everz kperp values
+ write(kperp_string,'(f6.4)') kparray_2D(ikr,ikz)
+ !!!!!!!!!!!! Electron matrices !!!!!!!!!!!!
+ ! get the self electron colision matrix
+ IF (CO .GT. 0) THEN
+ WRITE(mat_filename,'(a,a6,f6.4,a3)') TRIM(selfmat_file),'kperp_',TRIM(ADJUSTL(kperp_string)),'.h5'
+ ELSE
+ mat_filename = selfmat_file
+ ENDIF
+
+ CALL openf(mat_filename,fid1, 'r', 'D', mpicomm=MPI_COMM_WORLD);
+ CALL getatt(fid1,'/Caapj/Ceepj/','Pmaxe',pdime)
+ CALL getatt(fid1,'/Caapj/Ceepj/','Jmaxe',jdime)
+ IF ( ((pdime .LT. pmaxe) .OR. (jdime .LT. jmaxe)) .AND. (my_id .EQ. 0)) WRITE(*,*) '!! Sugama Matrix too small !!'
+
+ CALL allocate_array( Ceepj_full, 1,(pdime+1)*(jdime+1), 1,(pdime+1)*(jdime+1))
+ CALL getarr(fid1, '/Caapj/Ceepj', Ceepj_full) ! get array (moli format)
+ ! Fill sub array with only usefull polynmials degree
+ DO ip_e = 0,pmaxe ! Loop over rows
+ DO ij_e = 0,jmaxe
+ irow_sub = (jmaxe +1)*ip_e + ij_e +1
+ irow_full = (jdime +1)*ip_e + ij_e +1
+ DO il_e = 0,pmaxe ! Loop over columns
+ DO ik_e = 0,jmaxe
+ icol_sub = (jmaxe +1)*il_e + ik_e +1
+ icol_full = (jdime +1)*il_e + ik_e +1
+ Ceepj (irow_sub,icol_sub,ikr,ikz) = Ceepj_full (irow_full,icol_full)
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+
+ CALL closef(fid1)
+ DEALLOCATE(Ceepj_full)
+
+ ! get the Test and Back field electron ion collision matrix
+ IF (CO .GT. 0) THEN
+ WRITE(mat_filename,'(a,a6,f6.4,a3)') TRIM(eimat_file),'kperp_',TRIM(ADJUSTL(kperp_string)),'.h5'
+ ELSE
+ mat_filename = eimat_file
+ ENDIF
+
+ CALL openf(mat_filename,fid2, 'r', 'D', mpicomm=MPI_COMM_WORLD);
+
+ CALL getatt(fid2,'/Ceipj/CeipjT/','Pmaxi',pdimi)
+ CALL getatt(fid2,'/Ceipj/CeipjT/','Jmaxi',jdimi)
+ IF ( (pdimi .LT. pmaxi) .OR. (jdimi .LT. jmaxi) ) WRITE(*,*) '!! Sugama Matrix too small !!'
+
+ CALL allocate_array( CeipjT_full, 1,(pdime+1)*(jdime+1), 1,(pdime+1)*(jdime+1))
+ CALL allocate_array( CeipjF_full, 1,(pdime+1)*(jdime+1), 1,(pdimi+1)*(jdimi+1))
+
+ CALL getarr(fid2, '/Ceipj/CeipjT', CeipjT_full)
+ CALL getarr(fid2, '/Ceipj/CeipjF', CeipjF_full)
+
+ ! Fill sub array with only usefull polynmials degree
+ DO ip_e = 0,pmaxe ! Loop over rows
+ DO ij_e = 0,jmaxe
+ irow_sub = (jmaxe +1)*ip_e + ij_e +1
+ irow_full = (jdime +1)*ip_e + ij_e +1
+ DO il_e = 0,pmaxe ! Loop over columns
+ DO ik_e = 0,jmaxe
+ icol_sub = (jmaxe +1)*il_e + ik_e +1
+ icol_full = (jdime +1)*il_e + ik_e +1
+ CeipjT(irow_sub,icol_sub,ikr,ikz) = CeipjT_full(irow_full,icol_full)
+ ENDDO
+ ENDDO
+ DO il_i = 0,pmaxi ! Loop over columns
+ DO ik_i = 0,jmaxi
+ icol_sub = (Jmaxi +1)*il_i + ik_i +1
+ icol_full = (jdimi +1)*il_i + ik_i +1
+ CeipjF(irow_sub,icol_sub,ikr,ikz) = CeipjF_full(irow_full,icol_full)
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+
+ CALL closef(fid2)
+ DEALLOCATE(CeipjF_full)
+ DEALLOCATE(CeipjT_full)
+
+ !!!!!!!!!!!!!!! Ion matrices !!!!!!!!!!!!!!
+ ! get the self electron colision matrix
+ IF (CO .GT. 0) THEN
+ WRITE(mat_filename,'(a,a6,f6.4,a3)') TRIM(selfmat_file),'kperp_',TRIM(ADJUSTL(kperp_string)),'.h5'
+ ELSE
+ mat_filename = selfmat_file
+ ENDIF
+
+ CALL openf(mat_filename, fid3, 'r', 'D', mpicomm=MPI_COMM_WORLD);
+
+ CALL allocate_array( Ciipj_full, 1,(pdimi+1)*(jdimi+1), 1,(pdimi+1)*(jdimi+1))
+
+ IF ( (pmaxe .EQ. pmaxi) .AND. (jmaxe .EQ. jmaxi) ) THEN ! if same degrees ion and electron matrices are alike so load Ceepj
+ CALL getarr(fid3, '/Caapj/Ceepj', Ciipj_full) ! get array (moli format)
+ ELSE
+ CALL getarr(fid3, '/Caapj/Ciipj', Ciipj_full) ! get array (moli format)
+ ENDIF
+
+ ! Fill sub array with only usefull polynmials degree
+ DO ip_i = 0,Pmaxi ! Loop over rows
+ DO ij_i = 0,Jmaxi
+ irow_sub = (Jmaxi +1)*ip_i + ij_i +1
+ irow_full = (jdimi +1)*ip_i + ij_i +1
+ DO il_i = 0,Pmaxi ! Loop over columns
+ DO ik_i = 0,Jmaxi
+ icol_sub = (Jmaxi +1)*il_i + ik_i +1
+ icol_full = (jdimi +1)*il_i + ik_i +1
+ Ciipj (irow_sub,icol_sub,ikr,ikz) = Ciipj_full (irow_full,icol_full)
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+
+ CALL closef(fid3)
+ DEALLOCATE(Ciipj_full)
+
+ ! get the Test and Back field electron ion collision matrix
+ IF (CO .GT. 0) THEN
+ WRITE(mat_filename,'(a,a6,f6.4,a3)') TRIM(iemat_file),'kperp_',TRIM(ADJUSTL(kperp_string)),'.h5'
+ ELSE
+ mat_filename = iemat_file
+ ENDIF
+ write(*,*) 'Loading GK COSOLVER mat : ', mat_filename
+
+ CALL openf(mat_filename,fid4, 'r', 'D', mpicomm=MPI_COMM_WORLD);
+
+ CALL allocate_array( CiepjT_full, 1,(pdimi+1)*(jdimi+1), 1,(pdimi+1)*(jdimi+1))
+ CALL allocate_array( CiepjF_full, 1,(pdimi+1)*(jdimi+1), 1,(pdime+1)*(jdime+1))
+
+ CALL getarr(fid4, '/Ciepj/CiepjT', CiepjT_full)
+ CALL getarr(fid4, '/Ciepj/CiepjF', CiepjF_full)
+
+ ! Fill sub array with only usefull polynmials degree
+ DO ip_i = 0,Pmaxi ! Loop over rows
+ DO ij_i = 0,Jmaxi
+ irow_sub = (Jmaxi +1)*ip_i + ij_i +1
+ irow_full = (jdimi +1)*ip_i + ij_i +1
+ DO il_i = 0,Pmaxi ! Loop over columns
+ DO ik_i = 0,Jmaxi
+ icol_sub = (Jmaxi +1)*il_i + ik_i +1
+ icol_full = (jdimi +1)*il_i + ik_i +1
+ CiepjT(irow_sub,icol_sub,ikr,ikz) = CiepjT_full(irow_full,icol_full)
+ ENDDO
+ ENDDO
+ DO il_e = 0,pmaxe ! Loop over columns
+ DO ik_e = 0,jmaxe
+ icol_sub = (jmaxe +1)*il_e + ik_e +1
+ icol_full = (jdime +1)*il_e + ik_e +1
+ CiepjF(irow_sub,icol_sub,ikr,ikz) = CiepjF_full(irow_full,icol_full)
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+
+ CALL closef(fid4)
+ DEALLOCATE(CiepjF_full)
+ DEALLOCATE(CiepjT_full)
ENDDO
- ENDDO
-
- CALL closef(fid4)
- DEALLOCATE(CiepjF_full)
- DEALLOCATE(CiepjT_full)
+ ENDDO
+ IF (my_id .EQ. 0) WRITE(*,*) '============DONE==========='
- END SUBROUTINE load_FC_mat
- !******************************************************************************!
+END SUBROUTINE load_COSOlver_mat
+ !******************************************************************************!
end module collision
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 03ad393aada84ae2536ef627192dd8a2361d1132..1ba800c84dbb3c15766e642150634382438cfd73 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -37,16 +37,23 @@ MODULE grid
INTEGER, PUBLIC :: ir,iz ! counters
integer(C_INTPTR_T), PUBLIC :: local_nkr, local_nz
integer(C_INTPTR_T), PUBLIC :: local_nkr_offset, local_nz_offset
+ INTEGER, PUBLIC :: local_nkp
INTEGER, PUBLIC :: local_np_e, local_np_i
integer(C_INTPTR_T), PUBLIC :: local_np_e_offset, local_np_i_offset
+ INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: counts_np_e, counts_np_i
+ INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: displs_np_e, displs_np_i
! Grids containing position in fourier space
- REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: krarray
- REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: kzarray
+ REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: krarray
+ REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: kzarray
+ REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: kparray ! kperp array
+ REAL(dp), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: kparray_2D ! kperp array in 2D
+ INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: ikparray ! kperp indices array
+ INTEGER, DIMENSION(:,:), ALLOCATABLE, PUBLIC :: ikparray_2D ! to convert (ikr,ikz) -> ikp
REAL(dp), PUBLIC, PROTECTED :: deltakr, deltakz
- INTEGER, PUBLIC, PROTECTED :: ikrs, ikre, ikzs, ikze, a
+ INTEGER, PUBLIC, PROTECTED :: ikrs, ikre, ikzs, ikze, ikps, ikpe
INTEGER, PUBLIC, PROTECTED :: ikr_0, ikz_0 ! Indices of k-grid origin
- INTEGER, PUBLIC :: ikr, ikz, ip, ij ! counters
+ INTEGER, PUBLIC :: ikr, ikz, ip, ij, ikp ! counters
! Grid containing the polynomials degrees
INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: parray_e
@@ -61,7 +68,7 @@ MODULE grid
! Public Functions
PUBLIC :: init_1Dgrid_distr
PUBLIC :: set_pgrid, set_jgrid
- PUBLIC :: set_krgrid, set_kzgrid
+ PUBLIC :: set_krgrid, set_kzgrid, set_kpgrid
PUBLIC :: grid_readinputs, grid_outputinputs
PUBLIC :: bare, bari
@@ -84,11 +91,30 @@ CONTAINS
SUBROUTINE set_pgrid
USE prec_const
IMPLICIT NONE
- INTEGER :: ip
+ INTEGER :: ip, istart, iend, in
+ ! Local data distribution
CALL decomp1D(pmaxe+1, num_procs_p, rank_p, ips_e, ipe_e)
CALL decomp1D(pmaxi+1, num_procs_p, rank_p, ips_i, ipe_i)
-
+ local_np_e = ipe_e - ips_e + 1
+ local_np_i = ipe_i - ips_i + 1
+ ! List of shift and local numbers between the different processes (used in scatterv and gatherv)
+ ALLOCATE(counts_np_e (1:num_procs_p))
+ ALLOCATE(counts_np_i (1:num_procs_p))
+ ALLOCATE(displs_np_e (1:num_procs_p))
+ ALLOCATE(displs_np_i (1:num_procs_p))
+ DO in = 0,num_procs_p-1
+ CALL decomp1D(pmaxe+1, num_procs_p, in, istart, iend)
+ counts_np_e(in+1) = iend-istart+1
+ displs_np_e(in+1) = istart-1
+ CALL decomp1D(pmaxi+1, num_procs_p, in, istart, iend)
+ counts_np_i(in+1) = iend-istart+1
+ displs_np_i(in+1) = istart-1
+ ENDDO
+ write(*,*) rank_p, ': counts = ', counts_np_e
+ write(*,*) rank_p, ': disps = ', displs_np_e
+
+ ! local grid computation
ALLOCATE(parray_e(ips_e:ipe_e))
ALLOCATE(parray_i(ips_i:ipe_i))
DO ip = ips_e,ipe_e; parray_e(ip) = (ip-1); END DO
@@ -100,7 +126,6 @@ CONTAINS
! Precomputations
pmaxe_dp = real(pmaxe,dp)
pmaxi_dp = real(pmaxi,dp)
-
END SUBROUTINE set_pgrid
SUBROUTINE set_jgrid
@@ -123,7 +148,6 @@ CONTAINS
! Precomputations
jmaxe_dp = real(jmaxe,dp)
jmaxi_dp = real(jmaxi,dp)
-
END SUBROUTINE set_jgrid
@@ -166,7 +190,7 @@ CONTAINS
SUBROUTINE set_kzgrid
USE prec_const
IMPLICIT NONE
- INTEGER :: i_
+ INTEGER :: i_, counter
Nkz = Nz;
! Start and END indices of grid
@@ -200,6 +224,41 @@ CONTAINS
END DO
END SUBROUTINE set_kzgrid
+ SUBROUTINE set_kpgrid !Precompute the grid of kperp
+ IMPLICIT NONE
+ INTEGER :: ikz_sym, tmp_, counter
+ ! 2D to 1D indices array convertor
+ ALLOCATE(ikparray_2D(ikrs:ikre,ikzs:ikze))
+ ALLOCATE( kparray_2D(ikrs:ikre,ikzs:ikze))
+ ! local number of different kperp
+ local_nkp = local_nkr * (local_nkr-1)/2 + 1
+ ! Allocate 1D array of kperp values and indices
+ ALLOCATE(ikparray(1:local_nkr))
+ ALLOCATE( kparray(1:local_nkr))
+
+ ! Fill the arrays
+ tmp_ = 0; counter = 1
+ DO ikz = ikzs,ikze
+ DO ikr = ikrs,ikre
+ ! Set a symmetry on kz
+ IF (ikz .LE. Nkz/2+1) THEN
+ ikz_sym = ikz
+ ELSE
+ ikz_sym = Nkz+2 - ikz
+ ENDIF
+ ! Formula to find the 2D to 1D kperp equivalences ordered as
+ ! 10
+ ! 6 9
+ ! 3 5 8
+ !1 2 4 7 etc...
+ ikp = MAX(ikr-1,ikz_sym-1)*MIN(ikr-1,ikz_sym-1)/2 + min(ikr-1,ikz_sym-1)
+ ikparray_2D(ikr,ikz) = ikp
+ kparray_2D(ikr,ikz) = SQRT(krarray(ikr)**2 + kzarray(ikz)**2)
+ ENDDO
+ ENDDO
+
+ END SUBROUTINE
+
SUBROUTINE grid_readinputs
! Read the input parameters
USE prec_const
diff --git a/src/inital.F90 b/src/inital.F90
index f2768f6241fac64ac663b9df17a5b46ee8cdd61e..7aa10a3497a839b44929a4eb2c66aeccdbe4b0ca 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -58,11 +58,9 @@ SUBROUTINE inital
CALL build_dnjs_table
ENDIF
- !!!!!! Load the full coulomb collision operator coefficients !!!!!!
- IF (CO .EQ. -1) THEN
- IF (my_id .EQ. 0) WRITE(*,*) '=== Load Full Coulomb matrix ==='
- CALL load_FC_mat
- IF (my_id .EQ. 0) WRITE(*,*) '..done'
+ !!!!!! Load the COSOlver collision operator coefficients !!!!!!
+ IF (ABS(CO) .GT. 1) THEN
+ CALL load_COSOlver_mat
ENDIF
END SUBROUTINE inital
diff --git a/src/memory.F90 b/src/memory.F90
index 60b52058384a8e6d99c8f2e889d9e1aaea5f860f..8215f48a59cff9fe2c9baf45efcfc3c716d7c4bd 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -21,7 +21,7 @@ SUBROUTINE memory
! Electrostatic potential
CALL allocate_array(phi, ikrs,ikre, ikzs,ikze)
-
+
! Gyrocenter density *for 2D output*
CALL allocate_array(Ne00, ikrs,ikre, ikzs,ikze)
CALL allocate_array(Ni00, ikrs,ikre, ikzs,ikze)
@@ -31,16 +31,26 @@ SUBROUTINE memory
CALL allocate_array(Kernel_i, ijsg_i,ijeg_i, ikrs,ikre, ikzs,ikze)
! Collision matrix
- IF (CO .EQ. -1) THEN
- CALL allocate_array( Ceepj, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1))
- CALL allocate_array( CeipjT, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1))
- CALL allocate_array( CeipjF, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxi+1)*(jmaxi+1))
-
- CALL allocate_array( Ciipj, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxi+1)*(jmaxi+1))
- CALL allocate_array( CiepjT, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxi+1)*(jmaxi+1))
- CALL allocate_array( CiepjF, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxe+1)*(jmaxe+1))
+ IF (CO .LT. -1) THEN !DK collision matrix (same for every k)
+ CALL allocate_array( Ceepj, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), 1,1, 1,1)
+ CALL allocate_array( CeipjT, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), 1,1, 1,1)
+ CALL allocate_array( CeipjF, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxi+1)*(jmaxi+1), 1,1, 1,1)
+ CALL allocate_array( Ciipj, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxi+1)*(jmaxi+1), 1,1, 1,1)
+ CALL allocate_array( CiepjT, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxi+1)*(jmaxi+1), 1,1, 1,1)
+ CALL allocate_array( CiepjF, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxe+1)*(jmaxe+1), 1,1, 1,1)
+ ELSEIF (CO .GT. 1) THEN !GK collision matrices (one for each kperp)
+ CALL allocate_array( Ceepj, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), ikrs,ikre, ikzs,ikze)
+ CALL allocate_array( CeipjT, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), ikrs,ikre, ikzs,ikze)
+ CALL allocate_array( CeipjF, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxi+1)*(jmaxi+1), ikrs,ikre, ikzs,ikze)
+ CALL allocate_array( Ciipj, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxi+1)*(jmaxi+1), ikrs,ikre, ikzs,ikze)
+ CALL allocate_array( CiepjT, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxi+1)*(jmaxi+1), ikrs,ikre, ikzs,ikze)
+ CALL allocate_array( CiepjF, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxe+1)*(jmaxe+1), ikrs,ikre, ikzs,ikze)
ENDIF
+ ! Collision term
+ CALL allocate_array( TColl_e, ips_e,ipe_e, ijs_e,ije_e , ikrs,ikre, ikzs,ikze)
+ CALL allocate_array( TColl_i, ips_i,ipe_i, ijs_i,ije_i , ikrs,ikre, ikzs,ikze)
+
! Non linear terms and dnjs table
CALL allocate_array( Sepj, ips_e,ipe_e, ijs_e,ije_e, ikrs,ikre, ikzs,ikze )
CALL allocate_array( Sipj, ips_i,ipe_i, ijs_i,ije_i, ikrs,ikre, ikzs,ikze )
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index 00fb21cbeca0ba33921e506c8264aa847ed72dc9..eda2577064393a05453882ce587e1996533bb392 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -6,7 +6,8 @@ MODULE model
PRIVATE
INTEGER, PUBLIC, PROTECTED :: CO = 0 ! Collision Operator
- INTEGER, PUBLIC, PROTECTED :: CLOS = 0 ! Truncation method
+ INTEGER, PUBLIC, PROTECTED :: CLOS = 0 ! linear truncation method
+ INTEGER, PUBLIC, PROTECTED :: NL_CLOS = 0 ! nonlinear truncation method
INTEGER, PUBLIC, PROTECTED :: KERN = 0 ! Kernel model
LOGICAL, PUBLIC, PROTECTED :: NON_LIN = .true. ! To turn on non linear bracket term
REAL(dp), PUBLIC, PROTECTED :: mu = 0._dp ! spatial Hyperdiffusivity coefficient (for num. stability)
@@ -42,12 +43,11 @@ CONTAINS
SUBROUTINE model_readinputs
! Read the input parameters
-
USE basic, ONLY : lu_in
USE prec_const
IMPLICIT NONE
- NAMELIST /MODEL_PAR/ CO, CLOS, KERN, NON_LIN, mu, mu_p, mu_j, nu, tau_e, tau_i, sigma_e, sigma_i, &
+ NAMELIST /MODEL_PAR/ CO, CLOS, NL_CLOS, KERN, NON_LIN, mu, mu_p, mu_j, nu, tau_e, tau_i, sigma_e, sigma_i, &
q_e, q_i, eta_n, eta_T, eta_B, lambdaD
READ(lu_in,model_par)
@@ -77,7 +77,6 @@ CONTAINS
nu_i = nu * sigma_e * (tau_i)**(-3._dp/2._dp)/SQRT2 ! ion-ion collision frequ.
nu_ee = nu_e/SQRT2 ! e-e coll. frequ.
nu_ie = nu*sigma_e**2 ! i-e coll. frequ.
-
END SUBROUTINE model_readinputs
diff --git a/src/moments_eq_rhs.F90 b/src/moments_eq_rhs.F90
index 5f3411dcfe1a6b08dd130b47ef4ebb0889e8129b..d667d1867a7c856fe3478f157202f0203a8aa6bb 100644
--- a/src/moments_eq_rhs.F90
+++ b/src/moments_eq_rhs.F90
@@ -115,10 +115,10 @@ SUBROUTINE moments_eq_rhs_e
TNapjm1 = xNapjm1 * moments_e(ip,ij-1,ikr,ikz,updatetlevel)
!! Collision
- IF (CO .EQ. -3) THEN ! GK Dougherty
- CALL DoughertyGK_e(ip,ij,ikr,ikz,TColl)
+ IF (CO .EQ. 0) THEN ! Lenard Bernstein
+ TColl = xCapj * moments_e(ip,ij,ikr,ikz,updatetlevel)
- ELSEIF (CO .EQ. -2) THEN ! DK Dougherty
+ ELSEIF (CO .EQ. -1) THEN ! DK Dougherty
TColl20 = 0._dp; TColl01 = 0._dp; TColl10 = 0._dp
IF ( (pmaxe .GE. 2) ) TColl20 = xCa20 * moments_e(3,1,ikr,ikz,updatetlevel)
IF ( (jmaxe .GE. 1) ) TColl01 = xCa01 * moments_e(1,2,ikr,ikz,updatetlevel)
@@ -127,13 +127,12 @@ SUBROUTINE moments_eq_rhs_e
TColl = xCapj* moments_e(ip,ij,ikr,ikz,updatetlevel)&
+ TColl20 + TColl01 + TColl10
- ELSEIF (CO .EQ. -1) THEN ! Full Coulomb for electrons (COSOlver matrix)
- CALL FullCoulombDK_e(p_int,j_int,ikr,ikz,TColl)
-
- ELSEIF (CO .EQ. 0) THEN ! Lenard Bernstein
- TColl = xCapj * moments_e(ip,ij,ikr,ikz,updatetlevel)
+ ELSEIF (CO .EQ. 1) THEN ! GK Dougherty
+ CALL DoughertyGK_e(ip,ij,ikr,ikz,TColl)
- ENDIF
+ ELSE ! COSOLver matrix
+ TColl = TColl_e(ip,ij,ikr,ikz)
+ ENDIF
!! Electrical potential term
IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2
@@ -303,24 +302,22 @@ SUBROUTINE moments_eq_rhs_i
TNapjm1 = xNapjm1 * moments_i(ip,ij-1,ikr,ikz,updatetlevel)
!! Collision
- IF (CO .EQ. -3) THEN ! Gyrokin. Dougherty Collision terms
- CALL DoughertyGK_i(ip,ij,ikr,ikz,TColl)
+ IF (CO .EQ. 0) THEN ! Lenard Bernstein
+ TColl = xCapj * moments_i(ip,ij,ikr,ikz,updatetlevel)
- ELSEIF (CO .EQ. -2) THEN ! Dougherty Collision terms
+ ELSEIF (CO .EQ. -1) THEN ! DK Dougherty
TColl20 = 0._dp; TColl01 = 0._dp; TColl10 = 0._dp
- IF ( (pmaxe .GE. 2) ) TColl20 = xCa20 * moments_i(3,1,ikr,ikz,updatetlevel)
- IF ( (jmaxe .GE. 1) ) TColl01 = xCa01 * moments_i(1,2,ikr,ikz,updatetlevel)
- IF ( (pmaxe .GE. 1) ) TColl10 = xCa10 * moments_i(2,1,ikr,ikz,updatetlevel)
+ IF ( (pmaxi .GE. 2) ) TColl20 = xCa20 * moments_i(3,1,ikr,ikz,updatetlevel)
+ IF ( (jmaxi .GE. 1) ) TColl01 = xCa01 * moments_i(1,2,ikr,ikz,updatetlevel)
+ IF ( (pmaxi .GE. 1) ) TColl10 = xCa10 * moments_i(2,1,ikr,ikz,updatetlevel)
! Total collisional term
TColl = xCapj* moments_i(ip,ij,ikr,ikz,updatetlevel)&
+ TColl20 + TColl01 + TColl10
- ELSEIF (CO .EQ. -1) THEN !!! Full Coulomb for ions (COSOlver matrix) !!!
- CALL FullCoulombDK_i(p_int,j_int,ikr,ikz,TColl)
-
- ELSEIF (CO .EQ. 0) THEN! Lenhard Bernstein
- TColl = xCapj * moments_i(ip,ij,ikr,ikz,updatetlevel)
-
+ ELSEIF (CO .EQ. 1) THEN ! GK Dougherty
+ CALL DoughertyGK_i(ip,ij,ikr,ikz,TColl)
+ ELSE! COSOLver matrix (Sugama, Coulomb)
+ TColl = TColl_i(ip,ij,ikr,ikz)
ENDIF
!! Electrical potential term
diff --git a/src/stepon.F90 b/src/stepon.F90
index f3933416a2ba905aaccbb78462a07537c8de731c..af5ee588c8ac22a503f678f200d9d548487bd66d 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -1,17 +1,17 @@
SUBROUTINE stepon
! Advance one time step, (num_step=4 for Runge Kutta 4 scheme)
-
+ USE advance_field_routine, ONLY: advance_time_level, advance_field
+ USE array , ONLY: moments_rhs_e, moments_rhs_i, Sepj, Sipj
USE basic
- USE time_integration
+ USE closure
+ USE collision, ONLY : compute_TColl
USE fields, ONLY: moments_e, moments_i, phi
- USE array , ONLY: moments_rhs_e, moments_rhs_i, Sepj, Sipj
+ USE ghosts
USE grid
- USE advance_field_routine, ONLY: advance_time_level, advance_field
USE model
- USE closure
- USE utility, ONLY: checkfield
use prec_const
- USE ghosts
+ USE time_integration
+ USE utility, ONLY: checkfield
IMPLICIT NONE
@@ -26,6 +26,9 @@ SUBROUTINE stepon
! Exchanges the ghosts values updated
CALL update_ghosts
+ ! Compute collision
+ CALL compute_TColl
+
! Compute right hand side of moments hierarchy equation
CALL moments_eq_rhs_e
CALL moments_eq_rhs_i