diff --git a/src/closure_mod.F90 b/src/closure_mod.F90 index 82c42eb..4185a1a 100644 --- a/src/closure_mod.F90 +++ b/src/closure_mod.F90 @@ -1,107 +1,96 @@ module closure ! Contains the routines to define closures USE basic USE model, ONLY: CLOS, tau_e, tau_i, q_e, q_i, eta_B, nu USE grid USE array, ONLY: kernel_e, kernel_i USE fields, ONLY: moments_e, moments_i USE time_integration, ONLY: updatetlevel IMPLICIT NONE PUBLIC :: apply_closure_model CONTAINS ! Positive Oob indices are approximated with a model SUBROUTINE apply_closure_model IMPLICIT NONE - complex(dp) :: i_ky - real(dp) :: taue_qe_etaB_nu, taui_qi_etaB_nu - real(dp) :: sqpp2pp1_e, sqpp2pp1_i, sqpp1p_e, sqpp1p_i - real(dp) :: p_dp, j_dp - ! Spare some computations - taue_qe_etaB_nu = tau_e*eta_B/q_e/nu - taui_qi_etaB_nu = tau_i*eta_B/q_i/nu - sqpp2pp1_e = SQRT((pmaxe_dp+2)*(pmaxe_dp+1)) - sqpp2pp1_i = SQRT((pmaxi_dp+2)*(pmaxi_dp+1)) - sqpp1p_e = SQRT((pmaxe_dp+1)*(pmaxe_dp)) - sqpp1p_i = SQRT((pmaxi_dp+1)*(pmaxi_dp)) CALL cpu_time(t0_clos) ! zero truncation, An+1=0 for n+1>nmax IF (CLOS .EQ. 0) THEN DO ikx = ikxs,ikxe DO iky = ikys,ikye DO iz = izs,ize DO ip = ipsg_e,ipeg_e moments_e(ip,ijsg_e,ikx,iky,iz,updatetlevel) = 0._dp moments_e(ip,ijeg_e,ikx,iky,iz,updatetlevel) = 0._dp ENDDO DO ij = ijsg_e,ijeg_e moments_e(ipsg_e+1,ij,ikx,iky,iz,updatetlevel) = 0._dp moments_e(ipsg_e ,ij,ikx,iky,iz,updatetlevel) = 0._dp moments_e(ipeg_e-1,ij,ikx,iky,iz,updatetlevel) = 0._dp moments_e(ipeg_e ,ij,ikx,iky,iz,updatetlevel) = 0._dp ENDDO kernel_e(ijsg_e,ikx,iky) = 0._dp kernel_e(ijeg_e,ikx,iky) = 0._dp DO ip = ipsg_i,ipeg_i moments_i(ip,ijsg_i,ikx,iky,iz,updatetlevel) = 0._dp moments_i(ip,ijeg_i,ikx,iky,iz,updatetlevel) = 0._dp ENDDO DO ij = ijsg_i,ijeg_i moments_i(ipsg_i+1,ij,ikx,iky,iz,updatetlevel) = 0._dp moments_i(ipsg_i ,ij,ikx,iky,iz,updatetlevel) = 0._dp moments_i(ipeg_i-1,ij,ikx,iky,iz,updatetlevel) = 0._dp moments_i(ipeg_i ,ij,ikx,iky,iz,updatetlevel) = 0._dp ENDDO kernel_i(ijsg_i,ikx,iky) = 0._dp kernel_i(ijeg_i,ikx,iky) = 0._dp ENDDO ENDDO ENDDO ! zero truncation, An+1=0 for n+1>nmax ELSEIF (CLOS .EQ. 1) THEN DO ikx = ikxs,ikxe DO iky = ikys,ikye DO iz = izs,ize DO ip = ipsg_e,ipeg_e moments_e(ip,ijsg_e,ikx,iky,iz,updatetlevel) = 0._dp moments_e(ip,ijeg_e,ikx,iky,iz,updatetlevel) = 0._dp ENDDO DO ij = ijsg_e,ijeg_e moments_e(ipsg_e+1,ij,ikx,iky,iz,updatetlevel) = 0._dp moments_e(ipsg_e ,ij,ikx,iky,iz,updatetlevel) = 0._dp moments_e(ipeg_e-1,ij,ikx,iky,iz,updatetlevel) = 0._dp moments_e(ipeg_e ,ij,ikx,iky,iz,updatetlevel) = 0._dp ENDDO kernel_e(ijsg_e,ikx,iky) = 0._dp kernel_e(ijeg_e,ikx,iky) = 0._dp DO ip = ipsg_i,ipeg_i moments_i(ip,ijsg_i,ikx,iky,iz,updatetlevel) = 0._dp moments_i(ip,ijeg_i,ikx,iky,iz,updatetlevel) = 0._dp ENDDO DO ij = ijsg_i,ijeg_i moments_i(ipsg_i+1,ij,ikx,iky,iz,updatetlevel) = 0._dp moments_i(ipsg_i ,ij,ikx,iky,iz,updatetlevel) = 0._dp moments_i(ipeg_i-1,ij,ikx,iky,iz,updatetlevel) = 0._dp moments_i(ipeg_i ,ij,ikx,iky,iz,updatetlevel) = 0._dp ENDDO kernel_i(ijsg_i,ikx,iky) = 0._dp kernel_i(ijeg_i,ikx,iky) = 0._dp ENDDO ENDDO ENDDO ELSE if(my_id .EQ. 0) write(*,*) '! Closure scheme not found !' ENDIF CALL cpu_time(t1_clos) tc_clos = tc_clos + (t1_clos - t0_clos) END SUBROUTINE apply_closure_model END module closure diff --git a/src/model_mod.F90 b/src/model_mod.F90 index f06cf81..63a5e31 100644 --- a/src/model_mod.F90 +++ b/src/model_mod.F90 @@ -1,119 +1,119 @@ MODULE model ! Module for diagnostic parameters USE prec_const IMPLICIT NONE PRIVATE INTEGER, PUBLIC, PROTECTED :: CO = 0 ! Collision Operator 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) REAL(dp), PUBLIC, PROTECTED :: mu_p = 0._dp ! kinetic para hyperdiffusivity coefficient (for num. stability) REAL(dp), PUBLIC, PROTECTED :: mu_j = 0._dp ! kinetic perp hyperdiffusivity coefficient (for num. stability) REAL(dp), PUBLIC, PROTECTED :: nu = 1._dp ! Collision frequency REAL(dp), PUBLIC, PROTECTED :: tau_e = 1._dp ! Temperature REAL(dp), PUBLIC, PROTECTED :: tau_i = 1._dp ! REAL(dp), PUBLIC, PROTECTED :: sigma_e = 1._dp ! Mass REAL(dp), PUBLIC, PROTECTED :: sigma_i = 1._dp ! REAL(dp), PUBLIC, PROTECTED :: q_e = -1._dp ! Charge REAL(dp), PUBLIC, PROTECTED :: q_i = 1._dp ! REAL(dp), PUBLIC, PROTECTED :: eta_n = 1._dp ! Density gradient REAL(dp), PUBLIC, PROTECTED :: eta_T = 1._dp ! Temperature gradient REAL(dp), PUBLIC, PROTECTED :: eta_B = 1._dp ! Magnetic gradient REAL(dp), PUBLIC, PROTECTED :: lambdaD = 1._dp ! Debye length - REAL(dp), PUBLIC, PROTECTED :: taue_qe_etaB ! factor of the magnetic moment coupling - REAL(dp), PUBLIC, PROTECTED :: taui_qi_etaB ! + REAL(dp), PUBLIC, PROTECTED :: taue_qe ! factor of the magnetic moment coupling + REAL(dp), PUBLIC, PROTECTED :: taui_qi ! REAL(dp), PUBLIC, PROTECTED :: sqrtTaue_qe ! factor of parallel moment term REAL(dp), PUBLIC, PROTECTED :: sqrtTaui_qi ! REAL(dp), PUBLIC, PROTECTED :: qe_sigmae_sqrtTaue ! factor of parallel phi term REAL(dp), PUBLIC, PROTECTED :: qi_sigmai_sqrtTaui ! REAL(dp), PUBLIC, PROTECTED :: sigmae2_taue_o2 ! factor of the Kernel argument REAL(dp), PUBLIC, PROTECTED :: sigmai2_taui_o2 ! REAL(dp), PUBLIC, PROTECTED :: nu_e, nu_i ! electron-ion, ion-ion collision frequency REAL(dp), PUBLIC, PROTECTED :: nu_ee, nu_ie ! e-e, i-e coll. frequ. REAL(dp), PUBLIC, PROTECTED :: qe2_taue, qi2_taui ! factor of the gammaD sum REAL(dp), PUBLIC, PROTECTED :: INV_LIN_SYS = 1._dp ! Invert the sign of the linear RHS (for measuring damped eigenvalues) PUBLIC :: model_readinputs, model_outputinputs CONTAINS SUBROUTINE model_readinputs ! Read the input parameters USE basic, ONLY : lu_in USE prec_const IMPLICIT NONE 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) !Precompute species dependant factors IF( q_e .NE. 0._dp ) THEN - taue_qe_etaB = tau_e/q_e * eta_B ! factor of the magnetic moment coupling + taue_qe = tau_e/q_e ! factor of the magnetic moment coupling sqrtTaue_qe = sqrt(tau_e)/q_e ! factor of parallel moment term ELSE - taue_qe_etaB = 0._dp + taue_qe = 0._dp sqrtTaue_qe = 0._dp ENDIF - taui_qi_etaB = tau_i/q_i * eta_B ! factor of the magnetic moment coupling + taui_qi = tau_i/q_i ! factor of the magnetic moment coupling sqrtTaui_qi = sqrt(tau_i)/q_i ! factor of parallel moment term qe_sigmae_sqrtTaue = q_e/sigma_e/SQRT(tau_e) ! factor of parallel phi term qi_sigmai_sqrtTaui = q_i/sigma_i/SQRT(tau_i) qe2_taue = (q_e**2)/tau_e ! factor of the gammaD sum qi2_taui = (q_i**2)/tau_i sigmae2_taue_o2 = sigma_e**2 * tau_e/2._dp ! factor of the Kernel argument sigmai2_taui_o2 = sigma_i**2 * tau_i/2._dp !! We use the ion-ion collision as normalization with definition ! nu_ii = 4 sqrt(pi)/3 T_i^(-3/2) m_i^(-1/2) q^4 n_i0 ln(Lambda) ! nu_e = nu/sigma_e * (tau_e)**(3._dp/2._dp) ! electron-ion collision frequency (where already multiplied by 0.532) nu_i = nu ! ion-ion collision frequ. nu_ee = nu_e ! e-e coll. frequ. nu_ie = nu_i ! i-e coll. frequ. ! Old normalization (MOLI Jorge/Frei) ! nu_e = 0.532_dp*nu ! electron-ion collision frequency (where already multiplied by 0.532) ! nu_i = 0.532_dp*nu*sigma_e*tau_e**(-3._dp/2._dp)/SQRT2 ! ion-ion collision frequ. ! nu_ee = nu_e/SQRT2 ! e-e coll. frequ. ! nu_ie = 0.532_dp*nu*sigma_e**2 ! i-e coll. frequ. END SUBROUTINE model_readinputs SUBROUTINE model_outputinputs(fidres, str) ! Write the input parameters to the results_xx.h5 file USE futils, ONLY: attach USE prec_const IMPLICIT NONE INTEGER, INTENT(in) :: fidres CHARACTER(len=256), INTENT(in) :: str CALL attach(fidres, TRIM(str), "CO", CO) CALL attach(fidres, TRIM(str), "CLOS", CLOS) CALL attach(fidres, TRIM(str), "KERN", KERN) CALL attach(fidres, TRIM(str), "NON_LIN", NON_LIN) CALL attach(fidres, TRIM(str), "nu", nu) CALL attach(fidres, TRIM(str), "mu", mu) CALL attach(fidres, TRIM(str), "tau_e", tau_e) CALL attach(fidres, TRIM(str), "tau_i", tau_i) CALL attach(fidres, TRIM(str), "sigma_e", sigma_e) CALL attach(fidres, TRIM(str), "sigma_i", sigma_i) CALL attach(fidres, TRIM(str), "q_e", q_e) CALL attach(fidres, TRIM(str), "q_i", q_i) CALL attach(fidres, TRIM(str), "eta_n", eta_n) CALL attach(fidres, TRIM(str), "eta_T", eta_T) CALL attach(fidres, TRIM(str), "eta_B", eta_B) CALL attach(fidres, TRIM(str), "lambdaD", lambdaD) END SUBROUTINE model_outputinputs END MODULE model