diff --git a/src/advance_field.F90 b/src/advance_field.F90 index 9ebdba4..a26bda4 100644 --- a/src/advance_field.F90 +++ b/src/advance_field.F90 @@ -1,93 +1,82 @@ MODULE advance_field_routine USE prec_const implicit none CONTAINS SUBROUTINE advance_time_level USE basic USE time_integration use prec_const IMPLICIT NONE CALL set_updatetlevel(mod(updatetlevel,ntimelevel)+1) END SUBROUTINE advance_time_level SUBROUTINE advance_moments USE basic USE time_integration USE grid use prec_const USE model, ONLY: CLOS, KIN_E use fields, ONLY: moments_e, moments_i use array, ONLY: moments_rhs_e, moments_rhs_i IMPLICIT NONE INTEGER :: p_int, j_int CALL cpu_time(t0_adv_field) IF(KIN_E) THEN DO ip=ips_e,ipe_e p_int = parray_e(ip) DO ij=ijs_e,ije_e IF((CLOS .NE. 1) .OR. (ip-1+2*(ij-1)+1 .LE. dmaxe))& CALL advance_field(moments_e(ip,ij,ikxs:ikxe,ikys:ikye,izs:ize,:), moments_rhs_e(ip,ij,ikxs:ikxe,ikys:ikye,izs:ize,:)) ENDDO ENDDO ENDIF DO ip=ips_i,ipe_i p_int = parray_i(ip) DO ij=ijs_i,ije_i j_int = jarray_i(ij) IF((CLOS .NE. 1) .OR. (ip-1+2*(ij-1)+1 .LE. dmaxi))& CALL advance_field(moments_i(ip,ij,ikxs:ikxe,ikys:ikye,izs:ize,:), moments_rhs_i(ip,ij,ikxs:ikxe,ikys:ikye,izs:ize,:)) ENDDO ENDDO ! Execution time end CALL cpu_time(t1_adv_field) tc_adv_field = tc_adv_field + (t1_adv_field - t0_adv_field) END SUBROUTINE advance_moments SUBROUTINE advance_field( f, f_rhs ) USE basic USE time_integration USE array USE grid use prec_const use initial_par, ONLY: ACT_ON_MODES IMPLICIT NONE COMPLEX(dp), DIMENSION ( ikxs:ikxe, ikys:ikye, izs:ize, ntimelevel ) :: f COMPLEX(dp), DIMENSION ( ikxs:ikxe, ikys:ikye, izs:ize, ntimelevel ) :: f_rhs INTEGER :: istage SELECT CASE (updatetlevel) CASE(1) - DO iky=ikys,ikye - DO ikx=ikxs,ikxe - DO iz=izs,ize - DO istage=1,ntimelevel - f(ikx,iky,iz,1) = f(ikx,iky,iz,1) + dt*b_E(istage)*f_rhs(ikx,iky,iz,istage) - END DO - END DO - END DO - END DO + DO istage=1,ntimelevel + f(ikxs:ikxe,ikys:ikye,izs:ize,1) = f(ikxs:ikxe,ikys:ikye,izs:ize,1) & + + dt*b_E(istage)*f_rhs(ikxs:ikxe,ikys:ikye,izs:ize,istage) + END DO CASE DEFAULT - DO iky=ikys,ikye - DO ikx=ikxs,ikxe - DO iz=izs,ize - f(ikx,iky,iz,updatetlevel) = f(ikx,iky,iz,1); - DO istage=1,updatetlevel-1 - f(ikx,iky,iz,updatetlevel) = f(ikx,iky,iz,updatetlevel) + & - dt*A_E(updatetlevel,istage)*f_rhs(ikx,iky,iz,istage) - END DO - END DO - END DO - END DO + f(ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel) = f(ikxs:ikxe,ikys:ikye,izs:ize,1); + DO istage=1,updatetlevel-1 + f(ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel) = f(ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel) + & + dt*A_E(updatetlevel,istage)*f_rhs(ikxs:ikxe,ikys:ikye,izs:ize,istage) + END DO END SELECT END SUBROUTINE advance_field END MODULE advance_field_routine diff --git a/src/array_mod.F90 b/src/array_mod.F90 index 5ac91e6..419d3c6 100644 --- a/src/array_mod.F90 +++ b/src/array_mod.F90 @@ -1,81 +1,89 @@ MODULE array use prec_const implicit none ! Arrays to store the rhs, for time integration (ip,ij,ikx,iky,iz,updatetlevel) COMPLEX(dp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: moments_rhs_e COMPLEX(dp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: moments_rhs_i ! Arrays of non-adiabatique moments COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: nadiab_moments_e COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: nadiab_moments_i + ! Derivatives and interpolated moments + COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: ddz_nepj + COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: interp_nepj + COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: ddz2_Nepj + COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: ddz_nipj + COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: interp_nipj + COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: ddz2_Nipj + ! Arrays to store special initial modes (semi linear simulation) ! Zonal ones (ky=0) COMPLEX(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: moments_e_ZF COMPLEX(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: moments_i_ZF COMPLEX(dp), DIMENSION(:,:), ALLOCATABLE :: phi_ZF ! Entropy modes (kx=0) COMPLEX(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: moments_e_EM COMPLEX(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: moments_i_EM COMPLEX(dp), DIMENSION(:,:), ALLOCATABLE :: phi_EM ! Non linear term array (ip,ij,ikx,iky,iz) COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: Sepj ! electron COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: Sipj ! ion ! To load collision matrix (ip,ij,ikx,iky,iz) REAL(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: Ceepj, CeipjT REAL(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: CeipjF REAL(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: Ciipj, CiepjT REAL(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: CiepjF ! Collision term (ip,ij,ikx,iky,iz) 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 ! 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 :: 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_e, xphijp1_e, xphijm1_e REAL(dp), DIMENSION(:,:), ALLOCATABLE :: xphij_i, xphijp1_i, xphijm1_i - + ! Kernel function evaluation (ij,ikx,iky,iz,odd/even p) REAL(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: kernel_e REAL(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: kernel_i ! Poisson operator (ikx,iky,iz) REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: inv_poisson_op !! Diagnostics ! Gyrocenter density for electron and ions (ikx,iky,iz) COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: Ne00 COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: Ni00 ! particle density for electron and ions (ikx,iky,iz) COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: dens_e COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: dens_i ! particle fluid velocity for electron and ions (ikx,iky,iz) COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: upar_e COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: upar_i COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: uper_e COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: uper_i ! particle temperature for electron and ions (ikx,iky,iz) COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: Tpar_e COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: Tpar_i COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: Tper_e COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: Tper_i COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: temp_e COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: temp_i END MODULE array diff --git a/src/collision_mod.F90 b/src/collision_mod.F90 index eed589b..33954f7 100644 --- a/src/collision_mod.F90 +++ b/src/collision_mod.F90 @@ -1,898 +1,898 @@ module collision ! contains the Hermite-Laguerre collision operators. Solved using COSOlver. USE array USE basic USE fields USE futils USE grid USE model USE prec_const USE time_integration USE utility 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 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 LOGICAL, PUBLIC, PROTECTED :: cosolver_coll ! check if cosolver matrices are used 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 CONTAINS SUBROUTINE collision_readinputs ! Read the input parameters IMPLICIT NONE NAMELIST /COLLISION_PAR/ collision_model, gyrokin_CO, interspecies, mat_file READ(lu_in,collision_par) SELECT CASE(collision_model) CASE ('LB') ! Lenhard bernstein cosolver_coll = .false. interspecies = .false. CASE ('DG') ! Dougherty cosolver_coll = .false. interspecies = .false. CASE ('SG') ! Sugama cosolver_coll = .true. CASE ('LR') ! Lorentz (Pitch angle) cosolver_coll = .true. interspecies = .false. CASE ('LD') ! Landau cosolver_coll = .true. CASE ('none') cosolver_coll = .false. interspecies = .false. CASE DEFAULT ERROR STOP 'Error stop: collision model not recognized!!' END SELECT END SUBROUTINE collision_readinputs SUBROUTINE coll_outputinputs(fidres, str) ! Write the input parameters to the results_xx.h5 file IMPLICIT NONE INTEGER, INTENT(in) :: fidres CHARACTER(len=256), INTENT(in) :: str CHARACTER(len=2) :: gkco = 'dk' CHARACTER(len=2) :: abco = 'aa' CHARACTER(len=6) :: coname IF (gyrokin_CO) gkco = 'GK' IF (interspecies) abco = 'ab' WRITE(coname,'(a2,a2,a2)') collision_model,gkco,abco CALL attach(fidres, TRIM(str), "CO", coname) CALL attach(fidres, TRIM(str), "matfilename", mat_file) END SUBROUTINE coll_outputinputs SUBROUTINE compute_TColl USE basic USE model, ONLY : nu IMPLICIT NONE ! Execution time start CALL cpu_time(t0_coll) IF (nu .NE. 0) THEN SELECT CASE(collision_model) CASE ('LB') CALL compute_lenard_bernstein CASE ('DG') CALL compute_dougherty CASE ('SG','LR','LD') CALL compute_cosolver_coll CASE ('none') IF(KIN_E) & TColl_e = 0._dp TColl_i = 0._dp CASE DEFAULT ERROR STOP 'Error stop: collision operator not recognized!!' END SELECT ELSE IF(KIN_E) & TColl_e = 0._dp TColl_i = 0._dp ENDIF ! Execution time end CALL cpu_time(t1_coll) tc_coll = tc_coll + (t1_coll - t0_coll) END SUBROUTINE compute_TColl !******************************************************************************! !! Lenard Bernstein collision operator !******************************************************************************! SUBROUTINE compute_lenard_bernstein IMPLICIT NONE COMPLEX(dp) :: TColl_ IF (KIN_E) THEN DO ip = ips_e,ipe_e;DO ij = ijs_e,ije_e DO ikx = ikxs, ikxe;DO iky = ikys, ikye; DO iz = izs,ize CALL LenardBernstein_e(ip,ij,ikx,iky,iz,TColl_) TColl_e(ip,ij,ikx,iky,iz) = TColl_ ENDDO;ENDDO;ENDDO ENDDO;ENDDO ENDIF DO ip = ips_i,ipe_i;DO ij = ijs_i,ije_i DO ikx = ikxs, ikxe;DO iky = ikys, ikye; DO iz = izs,ize CALL LenardBernstein_i(ip,ij,ikx,iky,iz,TColl_) TColl_i(ip,ij,ikx,iky,iz) = TColl_ ENDDO;ENDDO;ENDDO ENDDO;ENDDO END SUBROUTINE compute_lenard_bernstein !******************************************************************************! !! for electrons !******************************************************************************! SUBROUTINE LenardBernstein_e(ip_,ij_,ikx_,iky_,iz_,TColl_) IMPLICIT NONE INTEGER, INTENT(IN) :: ip_,ij_,ikx_,iky_,iz_ COMPLEX(dp), INTENT(OUT) :: TColl_ REAL(dp) :: j_dp, p_dp, be_2, kp INTEGER :: eo_ !** Auxiliary variables ** eo_ = MODULO(parray_e(ip_),2) p_dp = REAL(parray_e(ip_),dp) j_dp = REAL(jarray_e(ij_),dp) kp = kparray(ikx_,iky_,iz_,eo_) be_2 = kp**2 * sigmae2_taue_o2 ! this is (be/2)^2 eo_ = MODULO(parray_e(ip_),2) !** Assembling collison operator ** ! -nuee (p + 2j) Nepj TColl_ = -nu_ee * (p_dp + 2._dp*j_dp)*nadiab_moments_e(ip_,ij_,ikx_,iky_,iz_) IF(gyrokin_CO) THEN TColl_ = TColl_ - nu_ee *2._dp*be_2*nadiab_moments_e(ip_,ij_,ikx_,iky_,iz_) ENDIF END SUBROUTINE LenardBernstein_e !******************************************************************************! !! for ions !******************************************************************************! SUBROUTINE LenardBernstein_i(ip_,ij_,ikx_,iky_,iz_,TColl_) USE fields, ONLY: moments_i USE grid, ONLY: parray_i, jarray_i, kxarray, kyarray USE basic USE model, ONLY: sigmai2_taui_o2, nu_i USE time_integration, ONLY : updatetlevel IMPLICIT NONE INTEGER, INTENT(IN) :: ip_,ij_,ikx_,iky_,iz_ COMPLEX(dp), INTENT(OUT) :: TColl_ REAL(dp) :: j_dp, p_dp, kp, bi_2 INTEGER :: eo_ !** Auxiliary variables ** eo_ = MODULO(parray_i(ip_),2) p_dp = REAL(parray_i(ip_),dp) j_dp = REAL(jarray_i(ij_),dp) kp = kparray(ikx_,iky_,iz_,eo_) bi_2 = kp**2 * sigmai2_taui_o2 ! this is (be/2)^2 !** Assembling collison operator ** ! -nuii (p + 2j) Nipj TColl_ = -nu_i * (p_dp + 2._dp*j_dp)*nadiab_moments_i(ip_,ij_,ikx_,iky_,iz_) IF(gyrokin_CO) THEN TColl_ = TColl_ - nu_i *2._dp*bi_2*nadiab_moments_i(ip_,ij_,ikx_,iky_,iz_) ENDIF END SUBROUTINE LenardBernstein_i !******************************************************************************! !! Doughtery collision operator !******************************************************************************! SUBROUTINE compute_dougherty IMPLICIT NONE COMPLEX(dp) :: TColl_ IF (KIN_E) THEN DO ip = ips_e,ipe_e;DO ij = ijs_e,ije_e DO ikx = ikxs, ikxe;DO iky = ikys, ikye; DO iz = izs,ize IF(gyrokin_CO) THEN CALL DoughertyGK_ee(ip,ij,ikx,iky,iz,TColl_) ELSE CALL DoughertyDK_ee(ip,ij,ikx,iky,iz,TColl_) ENDIF TColl_e(ip,ij,ikx,iky,iz) = TColl_ ENDDO;ENDDO;ENDDO ENDDO;ENDDO ENDIF DO ip = ips_i,ipe_i;DO ij = ijs_i,ije_i DO ikx = ikxs, ikxe;DO iky = ikys, ikye; DO iz = izs,ize IF(gyrokin_CO) THEN CALL DoughertyGK_ii(ip,ij,ikx,iky,iz,TColl_) ELSE CALL DoughertyDK_ii(ip,ij,ikx,iky,iz,TColl_) ENDIF TColl_i(ip,ij,ikx,iky,iz) = TColl_ ENDDO;ENDDO;ENDDO ENDDO;ENDDO END SUBROUTINE compute_dougherty !******************************************************************************! !! Doughtery driftkinetic collision operator for electrons !******************************************************************************! SUBROUTINE DoughertyDK_ee(ip_,ij_,ikx_,iky_,iz_,TColl_) IMPLICIT NONE INTEGER, INTENT(IN) :: ip_,ij_,ikx_,iky_,iz_ COMPLEX(dp), INTENT(OUT) :: TColl_ COMPLEX(dp) :: upar, Tpar, Tperp 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)*nadiab_moments_e(ip_,ij_,ikx_,iky_,iz_) IF( (p_dp .EQ. 1._dp) .AND. (j_dp .EQ. 0._dp)) THEN !Ce10 TColl_ = TColl_ + nadiab_moments_e(ip1_e,1,ikx_,iky_,iz_) ELSEIF( (p_dp .EQ. 2._dp) .AND. (j_dp .EQ. 0._dp)) THEN ! Ce20 TColl_ = TColl_ + twothird*nadiab_moments_e(ip2_e,1,ikx_,iky_,iz_) & - SQRT2*twothird*nadiab_moments_e(ip0_e,2,ikx_,iky_,iz_) ELSEIF( (p_dp .EQ. 0._dp) .AND. (j_dp .EQ. 1._dp)) THEN ! Ce01 TColl_ = TColl_ + 2._dp*twothird*nadiab_moments_e(ip0_e,2,ikx_,iky_,iz_) & - SQRT2*twothird*nadiab_moments_e(ip2_e,1,ikx_,iky_,iz_) ENDIF TColl_ = nu_ee * TColl_ END SUBROUTINE DoughertyDK_ee !******************************************************************************! !! Doughtery driftkinetic collision operator for ions !******************************************************************************! SUBROUTINE DoughertyDK_ii(ip_,ij_,ikx_,iky_,iz_,TColl_) IMPLICIT NONE INTEGER, INTENT(IN) :: ip_,ij_,ikx_,iky_,iz_ COMPLEX(dp), INTENT(OUT) :: TColl_ COMPLEX(dp) :: upar, Tpar, Tperp 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)*nadiab_moments_i(ip_,ij_,ikx_,iky_,iz_) IF( (p_dp .EQ. 1._dp) .AND. (j_dp .EQ. 0._dp)) THEN ! kronecker p1j0 TColl_ = TColl_ + nadiab_moments_i(ip1_i,1,ikx_,iky_,iz_) ELSEIF( (p_dp .EQ. 2._dp) .AND. (j_dp .EQ. 0._dp)) THEN ! kronecker p2j0 TColl_ = TColl_ + twothird*nadiab_moments_i(ip2_i,1,ikx_,iky_,iz_) & - SQRT2*twothird*nadiab_moments_i(ip0_i,2,ikx_,iky_,iz_) ELSEIF( (p_dp .EQ. 0._dp) .AND. (j_dp .EQ. 1._dp)) THEN ! kronecker p0j1 TColl_ = TColl_ + 2._dp*twothird*nadiab_moments_i(ip0_i,2,ikx_,iky_,iz_) & - SQRT2*twothird*nadiab_moments_i(ip2_i,1,ikx_,iky_,iz_) ENDIF TColl_ = nu_i * TColl_ END SUBROUTINE DoughertyDK_ii !******************************************************************************! !! Doughtery gyrokinetic collision operator for electrons !******************************************************************************! SUBROUTINE DoughertyGK_ee(ip_,ij_,ikx_,iky_,iz_,TColl_) IMPLICIT NONE INTEGER, INTENT(IN) :: ip_,ij_,ikx_,iky_,iz_ COMPLEX(dp), INTENT(OUT) :: TColl_ COMPLEX(dp) :: n_,upar_,uperp_,Tpar_, Tperp_, T_ COMPLEX(dp) :: nadiab_moment_0j REAL(dp) :: Knp0, Knp1, Knm1, kp INTEGER :: in_, eo_ REAL(dp) :: n_dp, j_dp, p_dp, be_, be_2 !** Auxiliary variables ** p_dp = REAL(parray_e(ip_),dp) eo_ = MODULO(parray_e(ip_),2) j_dp = REAL(jarray_e(ij_),dp) kp = kparray(ikx_,iky_,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 !** 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_,ikx_,iky_,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) ! Store the kernels for sparing readings Knp0 = Kernel_e(in_ ,ikx_,iky_,iz_,eo_) Knp1 = Kernel_e(in_+1,ikx_,iky_,iz_,eo_) Knm1 = Kernel_e(in_-1,ikx_,iky_,iz_,eo_) ! Nonadiabatic moments (only different from moments when p=0) nadiab_moment_0j = nadiab_moments_e(ip0_e,in_,ikx_,iky_,iz_) ! Density n_ = n_ + Knp0 * nadiab_moment_0j ! Perpendicular velocity uperp_ = uperp_ + be_*0.5_dp*(Knp0 - Knm1) * nadiab_moment_0j ! Parallel temperature Tpar_ = Tpar_ + Knp0 * (SQRT2*nadiab_moments_e(ip2_e,in_,ikx_,iky_,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_)/3._dp - n_ ! Add energy restoring term TColl_ = TColl_ + T_* 4._dp * j_dp * Kernel_e(ij_ ,ikx_,iky_,iz_,eo_) TColl_ = TColl_ - T_* 2._dp * (j_dp + 1._dp) * Kernel_e(ij_+1,ikx_,iky_,iz_,eo_) TColl_ = TColl_ - T_* 2._dp * j_dp * Kernel_e(ij_-1,ikx_,iky_,iz_,eo_) TColl_ = TColl_ + uperp_*be_* (Kernel_e(ij_,ikx_,iky_,iz_,eo_) - Kernel_e(ij_-1,ikx_,iky_,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 ! Parallel velocity upar_ = upar_ + Kernel_e(in_,ikx_,iky_,iz_,eo_) * nadiab_moments_e(ip1_e,in_,ikx_,iky_,iz_) ENDDO TColl_ = TColl_ + upar_*Kernel_e(ij_,ikx_,iky_,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) ! Store the kernels for sparing readings Knp0 = Kernel_e(in_ ,ikx_,iky_,iz_,eo_) Knp1 = Kernel_e(in_+1,ikx_,iky_,iz_,eo_) Knm1 = Kernel_e(in_-1,ikx_,iky_,iz_,eo_) ! Nonadiabatic moments (only different from moments when p=0) nadiab_moment_0j = nadiab_moments_e(ip0_e,in_,ikx_,iky_,iz_) ! Density n_ = n_ + Knp0 * nadiab_moment_0j ! Parallel temperature Tpar_ = Tpar_ + Knp0 * (SQRT2*nadiab_moments_e(ip2_e,in_,ikx_,iky_,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_)/3._dp - n_ TColl_ = TColl_ + T_*SQRT2*Kernel_e(ij_,ikx_,iky_,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_,ikx_,iky_,iz_,TColl_) IMPLICIT NONE INTEGER, INTENT(IN) :: ip_,ij_,ikx_,iky_,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(ikx_,iky_,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_,ikx_,iky_,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_ ,ikx_,iky_,iz_,eo_) Knp1 = Kernel_i(in_+1,ikx_,iky_,iz_,eo_) Knm1 = Kernel_i(in_-1,ikx_,iky_,iz_,eo_) ! Nonadiabatic moments (only different from moments when p=0) nadiab_moment_0j = nadiab_moments_i(ip0_i,in_,ikx_,iky_,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_,ikx_,iky_,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_ ,ikx_,iky_,iz_,eo_) TColl_ = TColl_ - T_* 2._dp * (j_dp + 1._dp) * Kernel_i(ij_+1,ikx_,iky_,iz_,eo_) TColl_ = TColl_ - T_* 2._dp * j_dp * Kernel_i(ij_-1,ikx_,iky_,iz_,eo_) TColl_ = TColl_ + uperp_*bi_* (Kernel_i(ij_,ikx_,iky_,iz_,eo_) - Kernel_i(ij_-1,ikx_,iky_,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_,ikx_,iky_,iz_,eo_) * nadiab_moments_i(ip1_i,in_,ikx_,iky_,iz_) ENDDO TColl_ = TColl_ + upar_*Kernel_i(ij_,ikx_,iky_,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_ ,ikx_,iky_,iz_,eo_) Knp1 = Kernel_i(in_+1,ikx_,iky_,iz_,eo_) Knm1 = Kernel_i(in_-1,ikx_,iky_,iz_,eo_) ! Nonadiabatic moments (only different from moments when p=0) nadiab_moment_0j = nadiab_moments_i(ip0_i,in_,ikx_,iky_,iz_) ! Density n_ = n_ + Knp0 * nadiab_moment_0j ! Parallel temperature Tpar_ = Tpar_ + Knp0 * (SQRT2*nadiab_moments_i(ip2_i,in_,ikx_,iky_,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_,ikx_,iky_,iz_,eo_) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ENDIF ! Multiply by ion-ion collision coefficient TColl_ = nu_i * TColl_ END SUBROUTINE DoughertyGK_ii !******************************************************************************! !! compute the collision terms in a (Np x Nj x Nkx x Nky) matrix all at once !******************************************************************************! SUBROUTINE compute_cosolver_coll IMPLICIT NONE COMPLEX(dp), DIMENSION(1:total_np_e) :: local_sum_e, buffer_e, total_sum_e COMPLEX(dp), DIMENSION(ips_e:ipe_e) :: TColl_distr_e COMPLEX(dp), DIMENSION(1:total_np_i) :: local_sum_i, buffer_i, total_sum_i COMPLEX(dp), DIMENSION(ips_i:ipe_i) :: TColl_distr_i COMPLEX(dp) :: TColl INTEGER :: ikxs_C, ikxe_C, ikys_C, ikye_C - DO ikx = ikxs,ikxe + DO iz = izs,ize DO iky = ikys,ikye - DO iz = izs,ize + DO ikx = ikxs,ikxe IF(KIN_E) THEN - ! Electrons - DO ij = 1,Jmaxe+1 - ! Loop over all p to compute sub collision term - DO ip = 1,total_np_e - CALL apply_COSOlver_mat_e(ip,ij,ikx,iky,iz,TColl) - local_sum_e(ip) = TColl - ENDDO - IF (num_procs_p .GT. 1) THEN - ! Sum up all the sub collision terms on root 0 - CALL MPI_REDUCE(local_sum_e, buffer_e, total_np_e, MPI_DOUBLE_COMPLEX, MPI_SUM, 0, comm_p, ierr) - ! 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) - ELSE - TColl_distr_e = local_sum_e - ENDIF - ! Write in output variable - DO ip = ips_e,ipe_e - TColl_e(ip,ij,ikx,iky,iz) = TColl_distr_e(ip) + DO ij = 1,Jmaxe+1 + ! Electrons + ! Loop over all p to compute sub collision term + DO ip = 1,total_np_e + CALL apply_COSOlver_mat_e(ip,ij,ikx,iky,iz,TColl) + local_sum_e(ip) = TColl + ENDDO + IF (num_procs_p .GT. 1) THEN + ! Sum up all the sub collision terms on root 0 + CALL MPI_REDUCE(local_sum_e, buffer_e, total_np_e, MPI_DOUBLE_COMPLEX, MPI_SUM, 0, comm_p, ierr) + ! 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) + ELSE + TColl_distr_e = local_sum_e + ENDIF + ! Write in output variable + DO ip = ips_e,ipe_e + TColl_e(ip,ij,ikx,iky,iz) = TColl_distr_e(ip) + ENDDO ENDDO - ENDDO ENDIF ! Ions DO ij = 1,Jmaxi+1 DO ip = 1,total_np_i CALL apply_COSOlver_mat_i(ip,ij,ikx,iky,iz,TColl) local_sum_i(ip) = TColl ENDDO IF (num_procs_p .GT. 1) THEN ! Reduce the local_sums to root = 0 CALL MPI_REDUCE(local_sum_i, buffer_i, total_np_i, MPI_DOUBLE_COMPLEX, MPI_SUM, 0, comm_p, ierr) ! buffer contains the entire collision term along p, we scatter it between ! 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) ELSE TColl_distr_i = local_sum_i ENDIF ! Write in output variable DO ip = ips_i,ipe_i TColl_i(ip,ij,ikx,iky,iz) = TColl_distr_i(ip) ENDDO ENDDO ENDDO ENDDO ENDDO END SUBROUTINE compute_cosolver_coll !******************************************************************************! - !!!!!!! Compute ion collision term + !!!!!!! Compute electron collision term !******************************************************************************! SUBROUTINE apply_COSOlver_mat_e(ip_,ij_,ikx_,iky_,iz_,TColl_) IMPLICIT NONE INTEGER, INTENT(IN) :: ip_, ij_ ,ikx_, iky_, iz_ COMPLEX(dp), INTENT(OUT) :: TColl_ INTEGER :: ip2,ij2, p_int,j_int, p2_int,j2_int, ikx_C, iky_C, iz_C p_int = parray_e_full(ip_); j_int = jarray_e_full(ij_); IF (gyrokin_CO) THEN ! GK operator (k-dependant) ikx_C = ikx_; iky_C = iky_; iz_C = iz_ ELSE ! DK operator (only one mat for every k) ikx_C = 1; iky_C = 1; iz_C = 1; ENDIF TColl_ = 0._dp ! Initialization of the local sum ! sum the electron-self and electron-ion test terms ploopee: DO ip2 = ips_e,ipe_e p2_int = parray_e(ip2) jloopee: DO ij2 = ijs_e,ije_e j2_int = jarray_e(ij2) IF((CLOS .NE. 1) .OR. (p2_int+2*j2_int .LE. dmaxe))& TColl_ = TColl_ + nadiab_moments_e(ip2,ij2,ikx_,iky_,iz_) & *( nu_e * CeipjT(bare(p_int,j_int), bare(p2_int,j2_int),ikx_C, iky_C, iz_C) & +nu_ee * Ceepj (bare(p_int,j_int), bare(p2_int,j2_int),ikx_C, iky_C, iz_C)) ENDDO jloopee ENDDO ploopee ! sum the electron-ion field terms ploopei: DO ip2 = ips_i,ipe_i p2_int = parray_i(ip2) jloopei: DO ij2 = ijs_i,ije_i j2_int = jarray_i(ij2) IF((CLOS .NE. 1) .OR. (p2_int+2*j2_int .LE. dmaxi))& TColl_ = TColl_ + nadiab_moments_i(ip2,ij2,ikx_,iky_,iz_) & *(nu_e * CeipjF(bare(p_int,j_int), bari(p2_int,j2_int),ikx_C, iky_C, iz_C)) END DO jloopei ENDDO ploopei END SUBROUTINE apply_COSOlver_mat_e !******************************************************************************! !!!!!!! Compute ion collision term !******************************************************************************! SUBROUTINE apply_COSOlver_mat_i(ip_,ij_,ikx_,iky_,iz_,TColl_) IMPLICIT NONE INTEGER, INTENT(IN) :: ip_, ij_ ,ikx_, iky_, iz_ COMPLEX(dp), INTENT(OUT) :: TColl_ INTEGER :: ip2,ij2, p_int,j_int, p2_int,j2_int, ikx_C, iky_C, iz_C p_int = parray_i_full(ip_); j_int = jarray_i_full(ij_); IF (gyrokin_CO) THEN ! GK operator (k-dependant) ikx_C = ikx_; iky_C = iky_; iz_C = iz_; ELSE ! DK operator (only one mat for every k) ikx_C = 1; iky_C = 1; iz_C = 1; ENDIF 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 = ijs_i,ije_i j2_int = jarray_i(ij2) IF((CLOS .NE. 1) .OR. (p2_int+2*j2_int .LE. dmaxi))& ! Ion-ion collision TColl_ = TColl_ + nadiab_moments_i(ip2,ij2,ikx_,iky_,iz_) & * nu_i * Ciipj (bari(p_int,j_int), bari(p2_int,j2_int), ikx_C, iky_C, iz_C) IF(KIN_E) & ! Ion-electron collision test TColl_ = TColl_ + nadiab_moments_i(ip2,ij2,ikx_,iky_,iz_) & * nu_ie * CiepjT(bari(p_int,j_int), bari(p2_int,j2_int), ikx_C, iky_C, iz_C) ENDDO jloopii ENDDO ploopii IF(KIN_E) THEN ! Ion-electron collision field ploopie: DO ip2 = ips_e,ipe_e ! sum the ion-electron field terms p2_int = parray_e(ip2) jloopie: DO ij2 = ijs_e,ije_e j2_int = jarray_e(ij2) IF((CLOS .NE. 1) .OR. (p2_int+2*j2_int .LE. dmaxe))& TColl_ = TColl_ + nadiab_moments_e(ip2,ij2,ikx_,iky_,iz_) & *(nu_ie * CiepjF(bari(p_int,j_int), bare(p2_int,j2_int), ikx_C, iky_C, iz_C)) ENDDO jloopie ENDDO ploopie ENDIF END SUBROUTINE apply_COSOlver_mat_i !******************************************************************************! !!!!!!! Load the collision matrix coefficient table from COSOlver results !******************************************************************************! SUBROUTINE load_COSOlver_mat ! Load a sub matrix from iCa files (works for pmaxa,jmaxa<=P_full,J_full) IMPLICIT NONE ! Indices for row and columns of the COSOlver matrix (4D compressed 2D matrices) INTEGER :: irow_sub, irow_full, icol_sub, icol_full INTEGER :: fid ! file indexation INTEGER :: ip_e, ij_e, il_e, ik_e, ikps_C, ikpe_C ! indices for electrons loops REAL(dp), DIMENSION(2) :: dims_e INTEGER :: pdime, jdime ! dimensions of the COSOlver matrices REAL(dp), DIMENSION(:,:), ALLOCATABLE :: Ceepj_full, CeipjT_full ! To load the entire matrix REAL(dp), DIMENSION(:,:), ALLOCATABLE :: CeipjF_full ! '' REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: Ceepj__kp, CeipjT_kp ! To store the coeff that will be used along kperp REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: CeipjF_kp ! '' INTEGER :: ip_i, ij_i, il_i, ik_i ! same for ions INTEGER, DIMENSION(2) :: dims_i INTEGER :: pdimi, jdimi ! dimensions of the COSOlver matrices REAL(dp), DIMENSION(:,:), ALLOCATABLE :: Ciipj_full, CiepjT_full ! . REAL(dp), DIMENSION(:,:), ALLOCATABLE :: CiepjF_full ! . REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: Ciipj__kp, CiepjT_kp ! . REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: CiepjF_kp ! . INTEGER :: NFLR REAL(dp), DIMENSION(:), ALLOCATABLE :: kp_grid_mat ! kperp grid of the matrices INTEGER :: ikp_next, ikp_prev, nkp_mat, ikp_mat REAL(dp) :: kp_max REAL(dp) :: kp_next, kp_prev, kperp_sim, kperp_mat, zerotoone, be_2, bi_2 CHARACTER(len=128) :: var_name, kperp_string, ikp_string !! Some terminal info SELECT CASE (collision_model) CASE ('SG') IF (my_id .EQ. 0) WRITE(*,*) '=== Load Sugama matrix ===' CASE ('LR') IF (my_id .EQ. 0) WRITE(*,*) '=== Load Lorentz matrix ===' CASE ('LD') IF (my_id .EQ. 0) WRITE(*,*) '=== Load Landau matrix ===' END SELECT SELECT CASE (gyrokin_CO) CASE (.true.) IF (my_id .EQ. 0) WRITE(*,*) '..gyrokinetic model..' CASE (.false.) IF (my_id .EQ. 0) WRITE(*,*) '..driftkinetic model..' END SELECT ! Opening the compiled cosolver matrices results if(my_id.EQ.0)write(*,*) mat_file CALL openf(mat_file,fid, 'r', 'D', mpicomm=comm_p); ! Get matrices dimensions (polynomials degrees and kperp grid) CALL getarr(fid, '/dims_e', dims_e) ! Get the electron polynomial degrees pdime = dims_e(1); jdime = dims_e(2); CALL getarr(fid, '/dims_i', dims_i) ! Get the ion polynomial degrees pdimi = dims_i(1); jdimi = dims_i(2); IF ( ((pdime .LT. pmaxe) .OR. (jdime .LT. jmaxe)) .AND. (my_id .EQ. 0)) WRITE(*,*) '!! Pe,Je Matrix too small !!' IF ( ((pdimi .LT. pmaxi) .OR. (jdimi .LT. jmaxi)) .AND. (my_id .EQ. 0)) WRITE(*,*) '!! Pi,Ji Matrix too small !!' CALL getsize(fid, '/coordkperp', nkp_mat) ! Get the dimension kperp grid of the matrices CALL allocate_array(kp_grid_mat, 1,nkp_mat) CALL getarr(fid, '/coordkperp', kp_grid_mat) ! check that we have enough kperps mat IF (LINEARITY .NE. 'linear') THEN IF ( (kp_grid_mat(nkp_mat) .LT. 2./3.*kp_max) .AND. (my_id .EQ. 0)) WRITE(*,*) '!! Matrix kperp grid too small !!' ELSE IF ( (kp_grid_mat(nkp_mat) .LT. kp_max) .AND. (my_id .EQ. 0)) WRITE(*,*) '!! Matrix kperp grid too small !!' ENDIF IF (gyrokin_CO) THEN ! GK operator (k-dependant) ikps_C = 1; ikpe_C = nkp_mat ELSE ! DK operator (only one mat for all k) ikps_C = 1; ikpe_C = 1 ENDIF CALL allocate_array( Ceepj__kp, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), ikps_C,ikpe_C) CALL allocate_array( CeipjT_kp, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), ikps_C,ikpe_C) CALL allocate_array( CeipjF_kp, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), ikps_C,ikpe_C) CALL allocate_array( Ciipj__kp, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), ikps_C,ikpe_C) CALL allocate_array( CiepjT_kp, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), ikps_C,ikpe_C) CALL allocate_array( CiepjF_kp, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), ikps_C,ikpe_C) DO ikp = ikps_C,ikpe_C ! Loop over everz kperp values ! Kperp value in string format to select in cosolver hdf5 file IF (gyrokin_CO) THEN write(ikp_string,'(i5.5)') ikp-1 ELSE write(ikp_string,'(i5.5)') 0 ENDIF !!!!!!!!!!!! E-E matrices !!!!!!!!!!!! ! get the self electron colision matrix ! Allocate space for storing full collision matrix CALL allocate_array( Ceepj_full, 1,(pdime+1)*(jdime+1), 1,(pdime+1)*(jdime+1)) ! Naming of the array to load (kperp dependant) WRITE(var_name,'(a,a)') TRIM(ADJUSTL(ikp_string)),'/Caapj/Ceepj' CALL getarr(fid, var_name, Ceepj_full) ! get array (moli format) ! Fill sub array with the usefull polynmial degrees only 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__kp (irow_sub,icol_sub,ikp) = Ceepj_full (irow_full,icol_full) ENDDO ENDDO ENDDO ENDDO DEALLOCATE(Ceepj_full) !!!!!!!!!!!!!!! I-I matrices !!!!!!!!!!!!!! ! get the self electron colision matrix CALL allocate_array( Ciipj_full, 1,(pdimi+1)*(jdimi+1), 1,(pdimi+1)*(jdimi+1)) WRITE(var_name,'(a,a,a)') TRIM(ADJUSTL(ikp_string)),'/Caapj/Ciipj' CALL getarr(fid, var_name, Ciipj_full) ! get array (moli format) ! 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__kp (irow_sub,icol_sub,ikp) = Ciipj_full (irow_full,icol_full) ENDDO ENDDO ENDDO ENDDO DEALLOCATE(Ciipj_full) IF(interspecies) THEN ! Pitch angle is only applied on like-species !!!!!!!!!!!!!!! E-I matrices !!!!!!!!!!!!!! ! Get test and field e-i collision matrices 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)) WRITE(var_name,'(a,a)') TRIM(ADJUSTL(ikp_string)),'/Ceipj/CeipjT' CALL getarr(fid, var_name, CeipjT_full) WRITE(var_name,'(a,a)') TRIM(ADJUSTL(ikp_string)),'/Ceipj/CeipjF' CALL getarr(fid, var_name, 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_kp(irow_sub,icol_sub,ikp) = 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_kp(irow_sub,icol_sub,ikp) = CeipjF_full(irow_full,icol_full) ENDDO ENDDO ENDDO ENDDO DEALLOCATE(CeipjF_full) DEALLOCATE(CeipjT_full) !!!!!!!!!!!!!!! I-E matrices !!!!!!!!!!!!!! ! get the Test and Back field electron ion collision matrix 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)) WRITE(var_name,'(a,a,a)') TRIM(ADJUSTL(ikp_string)),'/Ciepj/CiepjT' CALL getarr(fid, var_name, CiepjT_full) WRITE(var_name,'(a,a,a)') TRIM(ADJUSTL(ikp_string)),'/Ciepj/CiepjF' CALL getarr(fid, var_name, 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_kp(irow_sub,icol_sub,ikp) = 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_kp(irow_sub,icol_sub,ikp) = CiepjF_full(irow_full,icol_full) ENDDO ENDDO ENDDO ENDDO DEALLOCATE(CiepjF_full) DEALLOCATE(CiepjT_full) ELSE CeipjT_kp = 0._dp; CeipjF_kp = 0._dp; CiepjT_kp = 0._dp; CiepjF_kp = 0._dp; ENDIF ENDDO CALL closef(fid) IF (gyrokin_CO) THEN ! Interpolation of the kperp matrix values on kx ky grid IF (my_id .EQ. 0 ) WRITE(*,*) '...Interpolation from matrices kperp to simulation kx,ky...' DO ikx = ikxs,ikxe DO iky = ikys,ikye DO iz = izs,ize ! Check for nonlinear case if we are in the anti aliased domain or the filtered one kperp_sim = kparray(ikx,iky,iz,0) ! current simulation kperp ! Find the interval in kp grid mat where kperp_sim is contained ! Loop over the whole kp mat grid to find the smallest kperp that is ! larger than the current kperp_sim (brute force...) DO ikp=1,nkp_mat ikp_mat = ikp ! the first indice of the interval (k0) kperp_mat = kp_grid_mat(ikp) IF(kperp_mat .GT. kperp_sim) EXIT ! a matrix with kperp2 > current kx2 + ky2 has been found ENDDO ! Interpolation ! interval boundaries ikp_next = ikp_mat !index of k1 (larger than kperp_sim thanks to previous loop) ikp_prev = ikp_mat - 1 !index of k0 (smaller neighbour to interpolate inbetween) if ( (kp_grid_mat(ikp_prev) .GT. kperp_sim) .OR. (kp_grid_mat(ikp_next) .LT. kperp_sim) ) THEN ! write(*,*) 'Warning, linear interp of collision matrix failed!! ' ! write(*,*) kp_grid_mat(ikp_prev), '<', kperp_sim, '<', kp_grid_mat(ikp_next) ENDIF ! 0->1 variable for linear interp, i.e. zero2one = (k-k0)/(k1-k0) zerotoone = (kperp_sim - kp_grid_mat(ikp_prev))/(kp_grid_mat(ikp_next) - kp_grid_mat(ikp_prev)) ! Linear interpolation between previous and next kperp matrix values Ceepj (:,:,ikx,iky,iz) = (Ceepj__kp(:,:,ikp_next) - Ceepj__kp(:,:,ikp_prev))*zerotoone + Ceepj__kp(:,:,ikp_prev) Ciipj (:,:,ikx,iky,iz) = (Ciipj__kp(:,:,ikp_next) - Ciipj__kp(:,:,ikp_prev))*zerotoone + Ciipj__kp(:,:,ikp_prev) IF(interspecies) THEN CeipjT(:,:,ikx,iky,iz) = (CeipjT_kp(:,:,ikp_next) - CeipjT_kp(:,:,ikp_prev))*zerotoone + CeipjT_kp(:,:,ikp_prev) CeipjF(:,:,ikx,iky,iz) = (CeipjF_kp(:,:,ikp_next) - CeipjF_kp(:,:,ikp_prev))*zerotoone + CeipjF_kp(:,:,ikp_prev) CiepjT(:,:,ikx,iky,iz) = (CiepjT_kp(:,:,ikp_next) - CiepjT_kp(:,:,ikp_prev))*zerotoone + CiepjT_kp(:,:,ikp_prev) CiepjF(:,:,ikx,iky,iz) = (CiepjF_kp(:,:,ikp_next) - CiepjF_kp(:,:,ikp_prev))*zerotoone + CiepjF_kp(:,:,ikp_prev) ELSE CeipjT(:,:,ikx,iky,iz) = 0._dp CeipjF(:,:,ikx,iky,iz) = 0._dp CiepjT(:,:,ikx,iky,iz) = 0._dp CiepjF(:,:,ikx,iky,iz) = 0._dp ENDIF ENDDO ENDDO ENDDO ELSE ! DK -> No kperp dep, copy simply to final collision matrices Ceepj (:,:,1,1,1) = Ceepj__kp(:,:,1) CeipjT(:,:,1,1,1) = CeipjT_kp(:,:,1) CeipjF(:,:,1,1,1) = CeipjF_kp(:,:,1) Ciipj (:,:,1,1,1) = Ciipj__kp(:,:,1) CiepjT(:,:,1,1,1) = CiepjT_kp(:,:,1) CiepjF(:,:,1,1,1) = CiepjF_kp(:,:,1) ENDIF ! Deallocate auxiliary variables DEALLOCATE (Ceepj__kp); DEALLOCATE (CeipjT_kp); DEALLOCATE (CeipjF_kp) DEALLOCATE (Ciipj__kp); DEALLOCATE (CiepjT_kp); DEALLOCATE (CiepjF_kp) IF( .NOT. interspecies ) THEN IF(my_id.EQ.0) write(*,*) "--Like Species operator--" CeipjF = 0._dp; CeipjT = 0._dp; CiepjF = 0._dp; CiepjT = 0._dp; ENDIF IF (my_id .EQ. 0) WRITE(*,*) '============DONE===========' END SUBROUTINE load_COSOlver_mat !******************************************************************************! end module collision diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90 index ac081a8..2701295 100644 --- a/src/moments_eq_rhs_mod.F90 +++ b/src/moments_eq_rhs_mod.F90 @@ -1,269 +1,252 @@ MODULE moments_eq_rhs IMPLICIT NONE PUBLIC :: compute_moments_eq_rhs CONTAINS SUBROUTINE compute_moments_eq_rhs USE model, only: KIN_E IMPLICIT NONE IF(KIN_E) CALL moments_eq_rhs_e CALL moments_eq_rhs_i END SUBROUTINE compute_moments_eq_rhs !_____________________________________________________________________________! !_____________________________________________________________________________! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!! Electrons moments RHS !!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !_____________________________________________________________________________! SUBROUTINE moments_eq_rhs_e USE basic USE time_integration USE array USE fields USE grid USE model USE prec_const USE collision use geometry USE calculus, ONLY : interp_z, grad_z, grad_z2 IMPLICIT NONE INTEGER :: p_int, j_int ! loops indices and polynom. degrees REAL(dp) :: kx, ky, kperp2, dzlnB_o_J COMPLEX(dp) :: Tnepj, Tnepp2j, Tnepm2j, Tnepjp1, Tnepjm1 ! Terms from b x gradB and drives COMPLEX(dp) :: Tnepp1j, Tnepm1j, Tnepp1jm1, Tnepm1jm1 ! Terms from mirror force with non adiab moments COMPLEX(dp) :: Tperp, Tpar, Tmir, Tphi COMPLEX(dp) :: Unepm1j, Unepm1jp1, Unepm1jm1 ! Terms from mirror force with adiab moments - COMPLEX(dp) :: i_ky - ! To store derivatives and odd-even z grid interpolations - COMPLEX(dp), DIMENSION(izs:ize) :: ddznepp1j, ddznepm1j, & - nepp1j, nepp1jm1, nepm1j, nepm1jm1, nepm1jp1, ddz2Nepj + COMPLEX(dp) :: i_ky,phikxkyz ! Measuring execution time CALL cpu_time(t0_rhs) - ! Kinetic loops - ploope : DO ip = ips_e, ipe_e ! loop over Hermite degree indices - p_int = parray_e(ip) ! Hermite polynom. degree - eo = MODULO(p_int,2) ! Indicates if we are on even or odd z grid - jloope : DO ij = ijs_e, ije_e ! loop over Laguerre degree indices - j_int = jarray_e(ij) - ! Spatial loops - kxloope : DO ikx = ikxs,ikxe - kx = kxarray(ikx) ! radial wavevector - kyloope : DO iky = ikys,ikye + ! Spatial loops + zloope : DO iz = izs,ize + kyloope : DO iky = ikys,ikye 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 - ! Compute z derivatives and odd-even z interpolations - CALL grad_z(eo,nadiab_moments_e(ip+1,ij ,ikx,iky,izgs:izge), ddznepp1j(izs:ize)) - CALL grad_z(eo,nadiab_moments_e(ip-1,ij ,ikx,iky,izgs:izge), ddznepm1j(izs:ize)) - CALL interp_z(eo,nadiab_moments_e(ip+1,ij ,ikx,iky,izgs:izge), nepp1j (izs:ize)) - CALL interp_z(eo,nadiab_moments_e(ip+1,ij-1,ikx,iky,izgs:izge), nepp1jm1(izs:ize)) - CALL interp_z(eo,nadiab_moments_e(ip-1,ij ,ikx,iky,izgs:izge), nepm1j (izs:ize)) - CALL interp_z(eo,nadiab_moments_e(ip-1,ij+1,ikx,iky,izgs:izge), nepm1jp1(izs:ize)) - CALL interp_z(eo,nadiab_moments_e(ip-1,ij-1,ikx,iky,izgs:izge), nepm1jm1(izs:ize)) - ! Parallel hyperdiffusion - CALL grad_z2(moments_e(ip,ij,ikx,iky,izgs:izge,updatetlevel), ddz2Nepj(izs:ize)) + i_ky = imagu * ky ! toroidal derivative + + kxloope : DO ikx = ikxs,ikxe + kx = kxarray(ikx) ! radial wavevector + phikxkyz = phi(ikx,iky,iz)! tmp phi value + + ! Kinetic loops + jloope : DO ij = ijs_e, ije_e ! This loop is from 1 to jmaxi+1 + j_int = jarray_e(ij) + + ploope : DO ip = ips_e, ipe_e ! Hermite loop + p_int = parray_e(ip) ! Hermite degree + eo = MODULO(p_int,2) ! Indicates if we are on odd or even z grid + kperp2= kparray(ikx,iky,iz,eo)**2 - zloope : DO iz = izs,ize - ! kperp - kperp2= kparray(ikx,iky,iz,eo)**2 !! Compute moments mixing terms Tperp = 0._dp; Tpar = 0._dp; Tmir = 0._dp ! Perpendicular dynamic ! term propto n_e^{p,j} Tnepj = xnepj(ip,ij)* nadiab_moments_e(ip,ij,ikx,iky,iz) ! term propto n_e^{p+2,j} Tnepp2j = xnepp2j(ip) * nadiab_moments_e(ip+pp2,ij,ikx,iky,iz) ! term propto n_e^{p-2,j} Tnepm2j = xnepm2j(ip) * nadiab_moments_e(ip-pp2,ij,ikx,iky,iz) ! term propto n_e^{p,j+1} Tnepjp1 = xnepjp1(ij) * nadiab_moments_e(ip,ij+1,ikx,iky,iz) ! term propto n_e^{p,j-1} Tnepjm1 = xnepjm1(ij) * nadiab_moments_e(ip,ij-1,ikx,iky,iz) + ! Parallel dynamic ! ddz derivative for Landau damping term - Tpar = xnepp1j(ip) * ddznepp1j(iz) + xnepm1j(ip) * ddznepm1j(iz) - ! Mirror terms (trapping) - Tnepp1j = ynepp1j (ip,ij) * nepp1j (iz) - Tnepp1jm1 = ynepp1jm1(ip,ij) * nepp1jm1(iz) - Tnepm1j = ynepm1j (ip,ij) * nepm1j (iz) - Tnepm1jm1 = ynepm1jm1(ip,ij) * nepm1jm1(iz) + Tpar = xnepp1j(ip) * ddz_nepj(ip+1,ij,ikx,iky,iz) & + + xnepm1j(ip) * ddz_nepj(ip-1,ij,ikx,iky,iz) + ! Mirror terms + Tnepp1j = ynepp1j (ip,ij) * interp_nepj(ip+1,ij ,ikx,iky,iz) + Tnepp1jm1 = ynepp1jm1(ip,ij) * interp_nepj(ip+1,ij-1,ikx,iky,iz) + Tnepm1j = ynepm1j (ip,ij) * interp_nepj(ip-1,ij ,ikx,iky,iz) + Tnepm1jm1 = ynepm1jm1(ip,ij) * interp_nepj(ip-1,ij-1,ikx,iky,iz) ! Trapping terms - Unepm1j = znepm1j (ip,ij) * nepm1j (iz) - Unepm1jp1 = znepm1jp1(ip,ij) * nepm1jp1(iz) - Unepm1jm1 = znepm1jm1(ip,ij) * nepm1jm1(iz) + Unepm1j = znepm1j (ip,ij) * interp_nepj(ip-1,ij ,ikx,iky,iz) + Unepm1jp1 = znepm1jp1(ip,ij) * interp_nepj(ip-1,ij+1,ikx,iky,iz) + Unepm1jm1 = znepm1jm1(ip,ij) * interp_nepj(ip-1,ij-1,ikx,iky,iz) Tmir = Tnepp1j + Tnepp1jm1 + Tnepm1j + Tnepm1jm1 + Unepm1j + Unepm1jp1 + Unepm1jm1 !! Electrical potential term IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2 Tphi = (xphij_i (ip,ij)*kernel_e(ij ,ikx,iky,iz,eo) & + xphijp1_i(ip,ij)*kernel_e(ij+1,ikx,iky,iz,eo) & - + xphijm1_i(ip,ij)*kernel_e(ij-1,ikx,iky,iz,eo))*phi(ikx,iky,iz) + + xphijm1_i(ip,ij)*kernel_e(ij-1,ikx,iky,iz,eo))*phikxkyz ELSE Tphi = 0._dp ENDIF !! Sum of all RHS terms moments_rhs_e(ip,ij,ikx,iky,iz,updatetlevel) = & ! Perpendicular magnetic gradient/curvature effects - imagu*Ckxky(ikx,iky,iz,eo)*hatR(iz,eo)* (Tnepj + Tnepp2j + Tnepm2j + Tnepjp1 + Tnepjm1)& ! Parallel coupling (Landau Damping) - Tpar*gradz_coeff(iz,eo) & ! Mirror term (parallel magnetic gradient) - gradzB(iz,eo)* Tmir *gradz_coeff(iz,eo) & ! Drives (density + temperature gradients) - i_ky * Tphi & ! Numerical perpendicular hyperdiffusion (totally artificial, for stability purpose) - (mu_x*kx**4 + mu_y*ky**4)*moments_e(ip,ij,ikx,iky,iz,updatetlevel) & ! Numerical parallel hyperdiffusion "+ (mu_z*kz**4)" - + mu_z * diff_dz_coeff * ddz2Nepj(iz) & + + mu_z * diff_dz_coeff * ddz2_Nepj(ip,ij,ikx,iky,iz) & ! Collision term + TColl_e(ip,ij,ikx,iky,iz) & ! Nonlinear term - Sepj(ip,ij,ikx,iky,iz) - END DO zloope - END DO kyloope + END DO ploope + END DO jloope END DO kxloope - END DO jloope - END DO ploope + END DO kyloope + END DO zloope ! Execution time end CALL cpu_time(t1_rhs) tc_rhs = tc_rhs + (t1_rhs-t0_rhs) END SUBROUTINE moments_eq_rhs_e !_____________________________________________________________________________! !_____________________________________________________________________________! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!! Ions moments RHS !!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !_____________________________________________________________________________! SUBROUTINE moments_eq_rhs_i USE basic USE time_integration, ONLY: updatetlevel USE array USE fields USE grid USE model USE prec_const USE collision USE geometry USE calculus, ONLY : interp_z, grad_z, grad_z2 IMPLICIT NONE INTEGER :: p_int, j_int ! loops indices and polynom. degrees REAL(dp) :: kx, ky, kperp2 COMPLEX(dp) :: Tnipj, Tnipp2j, Tnipm2j, Tnipjp1, Tnipjm1 COMPLEX(dp) :: 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) :: Tperp, Tpar, Tmir, Tphi - COMPLEX(dp) :: i_ky - ! To store derivatives and odd-even z grid interpolations - COMPLEX(dp), DIMENSION(izs:ize) :: ddznipp1j, ddznipm1j, & - nipp1j, nipp1jm1, nipm1j, nipm1jm1, nipm1jp1, ddz2Nipj + COMPLEX(dp) :: i_ky, phikxkyz + ! Measuring execution time CALL cpu_time(t0_rhs) - ! Kinetic loops - ploopi : DO ip = ips_i, ipe_i ! Hermite loop - p_int = parray_i(ip) ! Hermite degree - eo = MODULO(p_int,2) ! Indicates if we are on odd or even z grid - jloopi : DO ij = ijs_i, ije_i ! This loop is from 1 to jmaxi+1 - j_int = jarray_i(ij) - ! Spatial loops - kxloopi : DO ikx = ikxs,ikxe - kx = kxarray(ikx) ! radial wavevector - kyloopi : DO iky = ikys,ikye - 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 - ! Compute z derivatives and odd-even z interpolations - CALL grad_z(eo,nadiab_moments_i(ip+1,ij ,ikx,iky,izgs:izge),ddznipp1j(izs:ize)) - CALL grad_z(eo,nadiab_moments_i(ip-1,ij ,ikx,iky,izgs:izge),ddznipm1j(izs:ize)) - CALL interp_z(eo,nadiab_moments_i(ip+1,ij ,ikx,iky,izgs:izge), nipp1j (izs:ize)) - CALL interp_z(eo,nadiab_moments_i(ip+1,ij-1,ikx,iky,izgs:izge), nipp1jm1(izs:ize)) - CALL interp_z(eo,nadiab_moments_i(ip-1,ij ,ikx,iky,izgs:izge), nipm1j (izs:ize)) - CALL interp_z(eo,nadiab_moments_i(ip-1,ij+1,ikx,iky,izgs:izge), nipm1jp1(izs:ize)) - CALL interp_z(eo,nadiab_moments_i(ip-1,ij-1,ikx,iky,izgs:izge), nipm1jm1(izs:ize)) - ! for hyperdiffusion - CALL grad_z2(moments_i(ip,ij,ikx,iky,izgs:izge,updatetlevel), ddz2Nipj(:)) - - zloopi : DO iz = izs,ize - kperp2= kparray(ikx,iky,iz,eo)**2 - - !! Compute moments mixing terms - Tperp = 0._dp; Tpar = 0._dp; Tmir = 0._dp - ! Perpendicular dynamic - ! term propto n_i^{p,j} - Tnipj = xnipj(ip,ij) * nadiab_moments_i(ip ,ij ,ikx,iky,iz) - ! term propto n_i^{p+2,j} - Tnipp2j = xnipp2j(ip) * nadiab_moments_i(ip+pp2,ij ,ikx,iky,iz) - ! term propto n_i^{p-2,j} - Tnipm2j = xnipm2j(ip) * nadiab_moments_i(ip-pp2,ij ,ikx,iky,iz) - ! term propto n_i^{p,j+1} - Tnipjp1 = xnipjp1(ij) * nadiab_moments_i(ip ,ij+1,ikx,iky,iz) - ! term propto n_i^{p,j-1} - Tnipjm1 = xnipjm1(ij) * nadiab_moments_i(ip ,ij-1,ikx,iky,iz) - ! Tperp - Tperp = Tnipj + Tnipp2j + Tnipm2j + Tnipjp1 + Tnipjm1 - ! Parallel dynamic - ! ddz derivative for Landau damping term - Tpar = xnipp1j(ip) * ddznipp1j(iz) + xnipm1j(ip) * ddznipm1j(iz) - ! Mirror terms - Tnipp1j = ynipp1j (ip,ij) * nipp1j (iz) - Tnipp1jm1 = ynipp1jm1(ip,ij) * nipp1jm1(iz) - Tnipm1j = ynipm1j (ip,ij) * nipm1j (iz) - Tnipm1jm1 = ynipm1jm1(ip,ij) * nipm1jm1(iz) - ! Trapping terms - Unipm1j = znipm1j (ip,ij) * nipm1j (iz) - Unipm1jp1 = znipm1jp1(ip,ij) * nipm1jp1(iz) - Unipm1jm1 = znipm1jm1(ip,ij) * nipm1jm1(iz) - - Tmir = Tnipp1j + Tnipp1jm1 + Tnipm1j + Tnipm1jm1 + Unipm1j + Unipm1jp1 + Unipm1jm1 - !! Electrical potential term - IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2 - Tphi = (xphij_i (ip,ij)*kernel_i(ij ,ikx,iky,iz,eo) & - + xphijp1_i(ip,ij)*kernel_i(ij+1,ikx,iky,iz,eo) & - + xphijm1_i(ip,ij)*kernel_i(ij-1,ikx,iky,iz,eo))*phi(ikx,iky,iz) - ELSE - Tphi = 0._dp - ENDIF - - !! Sum of all RHS terms - moments_rhs_i(ip,ij,ikx,iky,iz,updatetlevel) = & - ! Perpendicular magnetic gradient/curvature effects - - imagu*Ckxky(ikx,iky,iz,eo)*hatR(iz,eo) * Tperp & - ! Parallel coupling (Landau damping) - - gradz_coeff(iz,eo) * Tpar & - ! Mirror term (parallel magnetic gradient) - - gradzB(iz,eo) * gradz_coeff(iz,eo) * Tmir & - ! Drives (density + temperature gradients) - - i_ky * Tphi & - ! Numerical hyperdiffusion (totally artificial, for stability purpose) - - (mu_x*kx**4 + mu_y*ky**4)*moments_i(ip,ij,ikx,iky,iz,updatetlevel) & - ! Numerical parallel hyperdiffusion "+ (mu_z*kz**4)" - + mu_z * diff_dz_coeff * ddz2Nipj(iz) & - ! Collision term - + TColl_i(ip,ij,ikx,iky,iz)& - ! Nonlinear term - - Sipj(ip,ij,ikx,iky,iz) + ! Spatial loops + zloopi : DO iz = izs,ize + kyloopi : DO iky = ikys,ikye + ky = kyarray(iky) ! toroidal wavevector + i_ky = imagu * ky ! toroidal derivative - END DO zloopi - END DO kyloopi + kxloopi : DO ikx = ikxs,ikxe + kx = kxarray(ikx) ! radial wavevector + phikxkyz = phi(ikx,iky,iz)! tmp phi value + + ! Kinetic loops + jloopi : DO ij = ijs_i, ije_i ! This loop is from 1 to jmaxi+1 + j_int = jarray_i(ij) + + ploopi : DO ip = ips_i, ipe_i ! Hermite loop + p_int = parray_i(ip) ! Hermite degree + eo = MODULO(p_int,2) ! Indicates if we are on odd or even z grid + kperp2= kparray(ikx,iky,iz,eo)**2 + + !! Compute moments mixing terms + Tperp = 0._dp; Tpar = 0._dp; Tmir = 0._dp + ! Perpendicular dynamic + ! term propto n_i^{p,j} + Tnipj = xnipj(ip,ij) * nadiab_moments_i(ip ,ij ,ikx,iky,iz) + ! term propto n_i^{p+2,j} + Tnipp2j = xnipp2j(ip) * nadiab_moments_i(ip+pp2,ij ,ikx,iky,iz) + ! term propto n_i^{p-2,j} + Tnipm2j = xnipm2j(ip) * nadiab_moments_i(ip-pp2,ij ,ikx,iky,iz) + ! term propto n_i^{p,j+1} + Tnipjp1 = xnipjp1(ij) * nadiab_moments_i(ip ,ij+1,ikx,iky,iz) + ! term propto n_i^{p,j-1} + Tnipjm1 = xnipjm1(ij) * nadiab_moments_i(ip ,ij-1,ikx,iky,iz) + ! Tperp + Tperp = Tnipj + Tnipp2j + Tnipm2j + Tnipjp1 + Tnipjm1 + ! Parallel dynamic + ! ddz derivative for Landau damping term + Tpar = xnipp1j(ip) * ddz_nipj(ip+1,ij,ikx,iky,iz) & + + xnipm1j(ip) * ddz_nipj(ip-1,ij,ikx,iky,iz) + ! Mirror terms + Tnipp1j = ynipp1j (ip,ij) * interp_nipj(ip+1,ij ,ikx,iky,iz) + Tnipp1jm1 = ynipp1jm1(ip,ij) * interp_nipj(ip+1,ij-1,ikx,iky,iz) + Tnipm1j = ynipm1j (ip,ij) * interp_nipj(ip-1,ij ,ikx,iky,iz) + Tnipm1jm1 = ynipm1jm1(ip,ij) * interp_nipj(ip-1,ij-1,ikx,iky,iz) + ! Trapping terms + Unipm1j = znipm1j (ip,ij) * interp_nipj(ip-1,ij ,ikx,iky,iz) + Unipm1jp1 = znipm1jp1(ip,ij) * interp_nipj(ip-1,ij+1,ikx,iky,iz) + Unipm1jm1 = znipm1jm1(ip,ij) * interp_nipj(ip-1,ij-1,ikx,iky,iz) + + Tmir = Tnipp1j + Tnipp1jm1 + Tnipm1j + Tnipm1jm1 + Unipm1j + Unipm1jp1 + Unipm1jm1 + + !! Electrical potential term + IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2 + Tphi = (xphij_i (ip,ij)*kernel_i(ij ,ikx,iky,iz,eo) & + + xphijp1_i(ip,ij)*kernel_i(ij+1,ikx,iky,iz,eo) & + + xphijm1_i(ip,ij)*kernel_i(ij-1,ikx,iky,iz,eo))*phikxkyz + ELSE + Tphi = 0._dp + ENDIF + + !! Sum of all RHS terms + moments_rhs_i(ip,ij,ikx,iky,iz,updatetlevel) = & + ! Perpendicular magnetic gradient/curvature effects + - imagu*Ckxky(ikx,iky,iz,eo)*hatR(iz,eo) * Tperp & + ! Parallel coupling (Landau damping) + - gradz_coeff(iz,eo) * Tpar & + ! Mirror term (parallel magnetic gradient) + - gradzB(iz,eo) * gradz_coeff(iz,eo) * Tmir & + ! Drives (density + temperature gradients) + - i_ky * Tphi & + ! Numerical hyperdiffusion (totally artificial, for stability purpose) + - (mu_x*kx**4 + mu_y*ky**4)*moments_i(ip,ij,ikx,iky,iz,updatetlevel) & + ! Numerical parallel hyperdiffusion "+ (mu_z*kz**4)" + + mu_z * diff_dz_coeff * ddz2_Nipj(ip,ij,ikx,iky,iz) & + ! Collision term + + TColl_i(ip,ij,ikx,iky,iz)& + ! Nonlinear term + - Sipj(ip,ij,ikx,iky,iz) + + END DO ploopi + END DO jloopi END DO kxloopi - END DO jloopi - END DO ploopi + END DO kyloopi + END DO zloopi ! Execution time end CALL cpu_time(t1_rhs) tc_rhs = tc_rhs + (t1_rhs-t0_rhs) END SUBROUTINE moments_eq_rhs_i END MODULE moments_eq_rhs diff --git a/src/stepon.F90 b/src/stepon.F90 index c518a43..d928a4b 100644 --- a/src/stepon.F90 +++ b/src/stepon.F90 @@ -1,169 +1,166 @@ SUBROUTINE stepon ! Advance one time step, (num_step=4 for Runge Kutta 4 scheme) USE advance_field_routine, ONLY: advance_time_level, advance_field, advance_moments - USE array , ONLY: moments_rhs_e, moments_rhs_i, Sepj, Sipj USE basic USE closure - USE fields, ONLY: moments_e, moments_i, phi - USE initial_par, ONLY: ACT_ON_MODES USE ghosts USE grid USE model, ONLY : LINEARITY, KIN_E use prec_const USE time_integration USE numerics, ONLY: play_with_modes - USE processing, ONLY: compute_nadiab_moments, compute_fluid_moments USE utility, ONLY: checkfield IMPLICIT NONE INTEGER :: num_step LOGICAL :: mlend DO num_step=1,ntimelevel ! eg RK4 compute successively k1, k2, k3, k4 - !----- BEFORE: All fields are updated for step = n + !----- BEFORE: All fields+ghosts are updated for step = n ! Compute right hand side from current fields ! N_rhs(N_n, nadia_n, phi_n, S_n, Tcoll_n) CALL assemble_RHS ! ---- step n -> n+1 transition ! Advance from updatetlevel to updatetlevel+1 (according to num. scheme) CALL advance_time_level ! ---- ! Update moments with the hierarchy RHS (step by step) ! N_n+1 = N_n + N_rhs(n) CALL advance_moments ! Closure enforcement of moments CALL apply_closure_model ! Exchanges the ghosts values of N_n+1 CALL update_ghosts_p_moments CALL update_ghosts_z_moments ! Update electrostatic potential phi_n = phi(N_n+1) CALL poisson CALL update_ghosts_z_phi ! Numerical experiments ! Store or cancel/maintain zonal modes artificially CALL play_with_modes !- Check before next step CALL checkfield_all() IF( nlend ) EXIT ! exit do loop CALL MPI_BARRIER(MPI_COMM_WORLD,ierr) !----- AFTER: All fields are updated for step = n+1 END DO CONTAINS !!!! Basic structure to simplify stepon SUBROUTINE assemble_RHS USE moments_eq_rhs, ONLY: compute_moments_eq_rhs USE collision, ONLY: compute_TColl USE nonlinear, ONLY: compute_Sapj + USE processing, ONLY: compute_nadiab_moments_z_gradients_and_interp IMPLICIT NONE - ! compute auxiliary non adiabatic moments arrays - CALL compute_nadiab_moments + ! compute auxiliary non adiabatic moments array, gradients and interp + CALL compute_nadiab_moments_z_gradients_and_interp ! compute nonlinear term ("if linear" is included inside) CALL compute_Sapj ! compute collision term ("if coll, if nu >0" is included inside) CALL compute_TColl ! compute the moments equation rhs CALL compute_moments_eq_rhs END SUBROUTINE assemble_RHS SUBROUTINE checkfield_all ! Check all the fields for inf or nan ! Execution time start CALL cpu_time(t0_checkfield) IF(LINEARITY .NE. 'linear') CALL anti_aliasing ! ensure 0 mode for 2/3 rule IF(LINEARITY .NE. 'linear') CALL enforce_symmetry ! Enforcing symmetry on kx = 0 mlend=.FALSE. IF(.NOT.nlend) THEN mlend=mlend .or. checkfield(phi,' phi') IF(KIN_E) THEN - DO ip=ipgs_e,ipge_e - DO ij=ijgs_e,ijge_e + DO ij=ijgs_e,ijge_e + DO ip=ipgs_e,ipge_e mlend=mlend .or. checkfield(moments_e(ip,ij,:,:,:,updatetlevel),' moments_e') ENDDO ENDDO ENDIF - DO ip=ipgs_i,ipge_i - DO ij=ijgs_i,ijge_i + DO ij=ijgs_i,ijge_i + DO ip=ipgs_i,ipge_i mlend=mlend .or. checkfield(moments_i(ip,ij,:,:,:,updatetlevel),' moments_i') ENDDO ENDDO CALL MPI_ALLREDUCE(mlend, nlend, 1, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ierr) ENDIF ! Execution time end CALL cpu_time(t1_checkfield) tc_checkfield = tc_checkfield + (t1_checkfield - t0_checkfield) END SUBROUTINE checkfield_all SUBROUTINE anti_aliasing IF(KIN_E)THEN - DO ip=ipgs_e,ipge_e - DO ij=ijgs_e,ijge_e + DO iz=izgs,izge + DO iky=ikys,ikye DO ikx=ikxs,ikxe - DO iky=ikys,ikye - DO iz=izgs,izge + DO ij=ijgs_e,ijge_e + DO ip=ipgs_e,ipge_e moments_e( ip,ij,ikx,iky,iz,:) = AA_x(ikx)* AA_y(iky) * moments_e( ip,ij,ikx,iky,iz,:) END DO END DO END DO END DO END DO ENDIF - DO ip=ipgs_i,ipge_i - DO ij=ijs_i,ije_i + DO iz=izgs,izge + DO iky=ikys,ikye DO ikx=ikxs,ikxe - DO iky=ikys,ikye - DO iz=izgs,izge + DO ij=ijs_i,ije_i + DO ip=ipgs_i,ipge_i moments_i( ip,ij,ikx,iky,iz,:) = AA_x(ikx)* AA_y(iky) * moments_i( ip,ij,ikx,iky,iz,:) END DO END DO END DO END DO END DO END SUBROUTINE anti_aliasing SUBROUTINE enforce_symmetry ! Force X(k) = X(N-k)* complex conjugate symmetry IF ( contains_kx0 ) THEN ! Electron moments IF(KIN_E) THEN - DO ip=ipgs_e,ipge_e - DO ij=ijgs_e,ijge_e - DO iz=izs,ize - DO iky=2,Nky/2 !symmetry at kx = 0 - moments_e( ip,ij,ikx_0,iky,iz, :) = CONJG(moments_e( ip,ij,ikx_0,Nky+2-iky,iz, :)) - END DO + DO iz=izs,ize + DO ij=ijgs_e,ijge_e + DO ip=ipgs_e,ipge_e + DO iky=2,Nky/2 !symmetry at kx = 0 + moments_e( ip,ij,ikx_0,iky,iz, :) = CONJG(moments_e( ip,ij,ikx_0,Nky+2-iky,iz, :)) + END DO ! must be real at origin moments_e(ip,ij, ikx_0,iky_0,iz, :) = REAL(moments_e(ip,ij, ikx_0,iky_0,iz, :)) END DO END DO END DO ENDIF ! Ion moments - DO ip=ipgs_i,ipge_i + DO iz=izs,ize DO ij=ijgs_i,ijge_i - DO iz=izs,ize + DO ip=ipgs_i,ipge_i DO iky=2,Nky/2 !symmetry at kx = 0 moments_i( ip,ij,ikx_0,iky,iz, :) = CONJG(moments_i( ip,ij,ikx_0,Nky+2-iky,iz, :)) END DO ! must be real at origin and top right moments_i(ip,ij, ikx_0,iky_0,iz, :) = REAL(moments_i(ip,ij, ikx_0,iky_0,iz, :)) END DO END DO END DO ! Phi DO iky=2,Nky/2 !symmetry at kx = 0 phi(ikx_0,iky,izs:ize) = phi(ikx_0,Nky+2-iky,izs:ize) END DO ! must be real at origin phi(ikx_0,iky_0,izs:ize) = REAL(phi(ikx_0,iky_0,izs:ize)) ENDIF END SUBROUTINE enforce_symmetry END SUBROUTINE stepon