From 9f046171f545d5a6833d5e618279289be8d9d153 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Fri, 5 Aug 2022 17:54:06 +0200
Subject: [PATCH] Added eta_T, eta_N parameters to change gradients according
to species
---
src/model_mod.F90 | 14 +++++++++-----
src/moments_eq_rhs_mod.F90 | 26 +++++++++++++-------------
src/numerics_mod.F90 | 35 ++++++++++++++++++-----------------
3 files changed, 40 insertions(+), 35 deletions(-)
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index 12d16a58..1f1406fe 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -24,8 +24,10 @@ MODULE model
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 :: K_n = 1._dp ! Density drive
- REAL(dp), PUBLIC, PROTECTED :: K_T = 0._dp ! Temperature drive
+ REAL(dp), PUBLIC, PROTECTED :: k_N = 1._dp ! Ion density drive
+ REAL(dp), PUBLIC, PROTECTED :: eta_N = 1._dp ! electron-ion density drive ratio (k_Ne/k_Ni)
+ REAL(dp), PUBLIC, PROTECTED :: k_T = 0._dp ! Temperature drive
+ REAL(dp), PUBLIC, PROTECTED :: eta_T = 0._dp ! electron-ion temperature drive ratio (k_Te/k_Ti)
REAL(dp), PUBLIC, PROTECTED :: K_E = 0._dp ! Backg. electric field drive
REAL(dp), PUBLIC, PROTECTED :: GradB = 1._dp ! Magnetic gradient
REAL(dp), PUBLIC, PROTECTED :: CurvB = 1._dp ! Magnetic curvature
@@ -62,7 +64,7 @@ CONTAINS
NAMELIST /MODEL_PAR/ CLOS, NL_CLOS, KERN, LINEARITY, KIN_E, &
mu_x, mu_y, N_HD, mu_z, mu_p, mu_j, nu,&
tau_e, tau_i, sigma_e, sigma_i, q_e, q_i,&
- K_n, K_T, K_E, GradB, CurvB, lambdaD, beta
+ k_N, eta_N, k_T, eta_T, K_E, GradB, CurvB, lambdaD, beta
READ(lu_in,model_par)
@@ -132,8 +134,10 @@ CONTAINS
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), "K_n", K_n)
- CALL attach(fidres, TRIM(str), "K_T", K_T)
+ CALL attach(fidres, TRIM(str), "k_N", k_N)
+ CALL attach(fidres, TRIM(str), "eta_N", eta_N)
+ CALL attach(fidres, TRIM(str), "k_T", k_T)
+ CALL attach(fidres, TRIM(str), "eta_T", eta_T)
CALL attach(fidres, TRIM(str), "K_E", K_E)
CALL attach(fidres, TRIM(str), "lambdaD", lambdaD)
END SUBROUTINE model_outputinputs
diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90
index 2fe44a54..092b390d 100644
--- a/src/moments_eq_rhs_mod.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -300,7 +300,7 @@ SUBROUTINE add_Maxwellian_background_terms
! 40, 01,02, 21 with background gradient dependences.
USE prec_const
USE time_integration, ONLY : updatetlevel
- USE model, ONLY: taue_qe, taui_qi, K_N, K_T, KIN_E
+ USE model, ONLY: taue_qe, taui_qi, k_N, eta_N, k_T, eta_T, KIN_E
USE array, ONLY: moments_rhs_e, moments_rhs_i
USE grid, ONLY: contains_kx0, contains_ky0, ikx_0, iky_0,&
ips_e,ipe_e,ijs_e,ije_e,ips_i,ipe_i,ijs_i,ije_i,&
@@ -320,28 +320,28 @@ SUBROUTINE add_Maxwellian_background_terms
SELECT CASE (ip-1)
CASE(0) ! Na00 term
moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- +taue_qe * sinz(izs:ize) * (1.5_dp*K_N - 1.125_dp*K_T)
+ +taue_qe * sinz(izs:ize) * (1.5_dp*eta_N*k_N - 1.125_dp*eta_N*k_T)
CASE(2) ! Na20 term
moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- +taue_qe * sinz(izs:ize) * (SQRT2*0.5_dp*K_N - 2.75_dp*K_T)
+ +taue_qe * sinz(izs:ize) * (SQRT2*0.5_dp*eta_N*k_N - 2.75_dp*eta_T*k_T)
CASE(4) ! Na40 term
moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- +taue_qe * sinz(izs:ize) * SQRT6*0.75_dp*K_T
+ +taue_qe * sinz(izs:ize) * SQRT6*0.75_dp*eta_T*k_T
END SELECT
CASE(1) ! j = 1
SELECT CASE (ip-1)
CASE(0) ! Na01 term
moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- -taue_qe * sinz(izs:ize) * (K_N + 3.5_dp*K_T)
+ -taue_qe * sinz(izs:ize) * (eta_N*k_N + 3.5_dp*eta_T*k_T)
CASE(2) ! Na21 term
moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- -taue_qe * sinz(izs:ize) * SQRT2*K_T
+ -taue_qe * sinz(izs:ize) * SQRT2*eta_T*k_T
END SELECT
CASE(2) ! j = 2
SELECT CASE (ip-1)
CASE(0) ! Na02 term
moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- +taue_qe * sinz(izs:ize) * 2._dp*K_T
+ +taue_qe * sinz(izs:ize) * 2._dp*eta_T*k_T
END SELECT
END SELECT
ENDDO
@@ -355,28 +355,28 @@ SUBROUTINE add_Maxwellian_background_terms
SELECT CASE (ip-1)
CASE(0) ! Na00 term
moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- +taui_qi * sinz(izs:ize) * (1.5_dp*K_N - 1.125_dp*K_T)
+ +taui_qi * sinz(izs:ize) * (1.5_dp*k_N - 1.125_dp*k_T)
CASE(2) ! Na20 term
moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- +taui_qi * sinz(izs:ize) * (SQRT2*0.5_dp*K_N - 2.75_dp*K_T)
+ +taui_qi * sinz(izs:ize) * (SQRT2*0.5_dp*k_N - 2.75_dp*k_T)
CASE(4) ! Na40 term
moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- +taui_qi * sinz(izs:ize) * SQRT6*0.75_dp*K_T
+ +taui_qi * sinz(izs:ize) * SQRT6*0.75_dp*k_T
END SELECT
CASE(1) ! j = 1
SELECT CASE (ip-1)
CASE(0) ! Na01 term
moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- -taui_qi * sinz(izs:ize) * (K_N + 3.5_dp*K_T)
+ -taui_qi * sinz(izs:ize) * (k_N + 3.5_dp*k_T)
CASE(2) ! Na21 term
moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- -taui_qi * sinz(izs:ize) * SQRT2*K_T
+ -taui_qi * sinz(izs:ize) * SQRT2*k_T
END SELECT
CASE(2) ! j = 2
SELECT CASE (ip-1)
CASE(0) ! Na02 term
moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- +taui_qi * sinz(izs:ize) * 2._dp*K_T
+ +taui_qi * sinz(izs:ize) * 2._dp*k_T
END SELECT
END SELECT
ENDDO
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index 1ede5ac6..ccfb5f4b 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -201,7 +201,8 @@ END SUBROUTINE evaluate_ampere_op
SUBROUTINE compute_lin_coeff
USE array
USE model, ONLY: taue_qe, taui_qi, sqrtTaue_qe, sqrtTaui_qi, &
- K_T, K_n, CurvB, GradB, KIN_E, tau_e, tau_i, sigma_e, sigma_i
+ k_T, eta_T, k_N, eta_N, CurvB, GradB, KIN_E,&
+ tau_e, tau_i, sigma_e, sigma_i
USE prec_const
USE grid, ONLY: parray_e, parray_i, jarray_e, jarray_i, &
ip,ij, ips_e,ipe_e, ips_i,ipe_i, ijs_e,ije_e, ijs_i,ije_i
@@ -297,11 +298,11 @@ SUBROUTINE compute_lin_coeff
j_dp = REAL(j_int,dp) ! REALof Laguerre degree
!! Electrostatic potential pj terms
IF (p_int .EQ. 0) THEN ! kronecker p0
- xphij_e(ip,ij) =+K_n + 2.*j_dp*K_T
- xphijp1_e(ip,ij) =-K_T*(j_dp+1._dp)
- xphijm1_e(ip,ij) =-K_T* j_dp
+ xphij_e(ip,ij) =+eta_N*k_N + 2.*j_dp*eta_T*k_T
+ xphijp1_e(ip,ij) =-eta_T*k_T*(j_dp+1._dp)
+ xphijm1_e(ip,ij) =-eta_T*k_T* j_dp
ELSE IF (p_int .EQ. 2) THEN ! kronecker p2
- xphij_e(ip,ij) =+K_T/SQRT2
+ xphij_e(ip,ij) =+eta_T*k_T/SQRT2
xphijp1_e(ip,ij) = 0._dp; xphijm1_e(ip,ij) = 0._dp;
ELSE
xphij_e(ip,ij) = 0._dp; xphijp1_e(ip,ij) = 0._dp
@@ -317,11 +318,11 @@ SUBROUTINE compute_lin_coeff
j_dp = REAL(j_int,dp) ! REALof Laguerre degree
!! Electrostatic potential pj terms
IF (p_int .EQ. 0) THEN ! kronecker p0
- xphij_i(ip,ij) =+K_n + 2._dp*j_dp*K_T
- xphijp1_i(ip,ij) =-K_T*(j_dp+1._dp)
- xphijm1_i(ip,ij) =-K_T* j_dp
+ xphij_i(ip,ij) =+k_N + 2._dp*j_dp*k_T
+ xphijp1_i(ip,ij) =-k_T*(j_dp+1._dp)
+ xphijm1_i(ip,ij) =-k_T* j_dp
ELSE IF (p_int .EQ. 2) THEN ! kronecker p2
- xphij_i(ip,ij) =+K_T/SQRT2
+ xphij_i(ip,ij) =+k_T/SQRT2
xphijp1_i(ip,ij) = 0._dp; xphijm1_i(ip,ij) = 0._dp;
ELSE
xphij_i(ip,ij) = 0._dp; xphijp1_i(ip,ij) = 0._dp
@@ -339,11 +340,11 @@ SUBROUTINE compute_lin_coeff
j_dp = REAL(j_int,dp) ! REALof Laguerre degree
!! Electrostatic potential pj terms
IF (p_int .EQ. 1) THEN ! kronecker p1
- xpsij_e(ip,ij) =+(K_n + (2._dp*j_dp+1._dp)*K_T) * SQRT(tau_e)/sigma_e
- xpsijp1_e(ip,ij) =- K_T*(j_dp+1._dp) * SQRT(tau_e)/sigma_e
- xpsijm1_e(ip,ij) =- K_T* j_dp * SQRT(tau_e)/sigma_e
+ xpsij_e(ip,ij) =+(eta_N*k_N + (2._dp*j_dp+1._dp)*eta_T*k_T) * SQRT(tau_e)/sigma_e
+ xpsijp1_e(ip,ij) =- eta_T*k_T*(j_dp+1._dp) * SQRT(tau_e)/sigma_e
+ xpsijm1_e(ip,ij) =- eta_T*k_T* j_dp * SQRT(tau_e)/sigma_e
ELSE IF (p_int .EQ. 3) THEN ! kronecker p3
- xpsij_e(ip,ij) =+ K_T*SQRT3/SQRT2 * SQRT(tau_e)/sigma_e
+ xpsij_e(ip,ij) =+ eta_T*k_T*SQRT3/SQRT2 * SQRT(tau_e)/sigma_e
xpsijp1_e(ip,ij) = 0._dp; xpsijm1_e(ip,ij) = 0._dp;
ELSE
xpsij_e(ip,ij) = 0._dp; xpsijp1_e(ip,ij) = 0._dp
@@ -359,11 +360,11 @@ SUBROUTINE compute_lin_coeff
j_dp = REAL(j_int,dp) ! REALof Laguerre degree
!! Electrostatic potential pj terms
IF (p_int .EQ. 1) THEN ! kronecker p1
- xpsij_i(ip,ij) =+(K_n + (2._dp*j_dp+1._dp)*K_T) * SQRT(tau_i)/sigma_i
- xpsijp1_i(ip,ij) =- K_T*(j_dp+1._dp) * SQRT(tau_i)/sigma_i
- xpsijm1_i(ip,ij) =- K_T* j_dp * SQRT(tau_i)/sigma_i
+ xpsij_i(ip,ij) =+(k_N + (2._dp*j_dp+1._dp)*k_T) * SQRT(tau_i)/sigma_i
+ xpsijp1_i(ip,ij) =- k_T*(j_dp+1._dp) * SQRT(tau_i)/sigma_i
+ xpsijm1_i(ip,ij) =- k_T* j_dp * SQRT(tau_i)/sigma_i
ELSE IF (p_int .EQ. 3) THEN ! kronecker p3
- xpsij_i(ip,ij) =+ K_T*SQRT3/SQRT2 * SQRT(tau_i)/sigma_i
+ xpsij_i(ip,ij) =+ k_T*SQRT3/SQRT2 * SQRT(tau_i)/sigma_i
xpsijp1_i(ip,ij) = 0._dp; xpsijm1_i(ip,ij) = 0._dp;
ELSE
xpsij_i(ip,ij) = 0._dp; xpsijp1_i(ip,ij) = 0._dp
--
GitLab