From a042d51948d61a0223a7889a40719f18533172b2 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Thu, 4 Apr 2024 09:17:38 +0200
Subject: [PATCH] Add new parameters tuners for mirror, trapping and Landau
damping terms
---
src/model_mod.F90 | 10 +++++++---
src/numerics_mod.F90 | 28 ++++++++++++++--------------
2 files changed, 21 insertions(+), 17 deletions(-)
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index 7edf6e8..da0357d 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -21,8 +21,11 @@ MODULE model
REAL(xp), PUBLIC, PROTECTED :: tau_i = 1.0 ! electron-ion temperature ratio for ion adiabatic model
REAL(xp), PUBLIC, PROTECTED :: q_i = 1.0 ! ion charge for ion adiabatic model
REAL(xp), PUBLIC, PROTECTED :: nu = 0._xp ! collision frequency parameter
- REAL(xp), PUBLIC, PROTECTED :: k_gB = 1._xp ! Magnetic gradient strength (L_ref/L_gB)
- REAL(xp), PUBLIC, PROTECTED :: k_cB = 1._xp ! Magnetic curvature strength (L_ref/L_cB)
+ REAL(xp), PUBLIC, PROTECTED :: k_gB = 1._xp ! artificial magnetic gradient tuner (L_ref/L_gB)
+ REAL(xp), PUBLIC, PROTECTED :: k_cB = 1._xp ! artificial magnetic curvature tuner (L_ref/L_cB)
+ REAL(xp), PUBLIC, PROTECTED :: k_mB = 1._xp ! artificial mirror force tuner (L_ref/L_cB)
+ REAL(xp), PUBLIC, PROTECTED :: k_tB = 1._xp ! artificial trapping term tuner (L_ref/L_cB)
+ REAL(xp), PUBLIC, PROTECTED :: k_ldB = 1._xp ! artificial Landau damping tuner (L_ref/L_cB)
REAL(xp), PUBLIC, PROTECTED :: lambdaD = 0._xp ! Debye length
REAL(xp), PUBLIC, PROTECTED :: beta = 0._xp ! electron plasma Beta (8piNT_e/B0^2)
REAL(xp), PUBLIC, PROTECTED :: ExBrate = 0._xp ! ExB background shearing rate (radially constant shear flow)
@@ -59,7 +62,8 @@ CONTAINS
NAMELIST /MODEL/ LINEARITY, RM_LD_T_EQ, FORCE_SYMMETRY, MHD_PD, &
Na, ADIAB_E, ADIAB_I, q_i, tau_i, &
mu_x, mu_y, N_HD, HDz_h, mu_z, mu_p, mu_j, HYP_V, &
- nu, k_gB, k_cB, lambdaD, beta, ExBrate, ExB_NL_CORRECTION,&
+ nu, k_gB, k_cB, k_mB, k_tB, k_ldB, &
+ lambdaD, beta, ExBrate, ExB_NL_CORRECTION,&
ikxZF, ZFrate, ZF_ONLY, KN_MODEL, ORDER, ORDER_NUM, ORDER_DEN
READ(lu_in,model)
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index ebab9ba..6101fcb 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -287,7 +287,7 @@ SUBROUTINE compute_lin_coeff
xphij, xphijp1, xphijm1,&
xpsij, xpsijp1, xpsijm1
USE species, ONLY: k_T, k_N, tau, q, sqrt_tau_o_sigma
- USE model, ONLY: k_cB, k_gB
+ USE model, ONLY: k_cB, k_gB, k_mB, k_tB, k_ldB
USE prec_const, ONLY: xp, SQRT2, SQRT3
USE grid, ONLY: parray, jarray, local_na, local_np, local_nj, ngj, ngp
INTEGER :: ia,ip,ij,p_int, j_int ! polynom. dagrees
@@ -305,32 +305,32 @@ SUBROUTINE compute_lin_coeff
xnapj(ia,ip,ij) = tau(ia)/q(ia)*(k_cB*(2._xp*p_xp + 1._xp) &
+k_gB*(2._xp*j_xp + 1._xp))
! Mirror force terms
- ynapp1j (ia,ip,ij) = -sqrt_tau_o_sigma(ia) * (j_xp+1._xp)*SQRT(p_xp+1._xp)
- ynapm1j (ia,ip,ij) = -sqrt_tau_o_sigma(ia) * (j_xp+1._xp)*SQRT(p_xp)
- ynapp1jm1(ia,ip,ij) = +sqrt_tau_o_sigma(ia) * j_xp*SQRT(p_xp+1._xp)
- ynapm1jm1(ia,ip,ij) = +sqrt_tau_o_sigma(ia) * j_xp*SQRT(p_xp)
+ ynapp1j (ia,ip,ij) = -sqrt_tau_o_sigma(ia) * (j_xp+1._xp)*SQRT(p_xp+1._xp) * k_mB
+ ynapm1j (ia,ip,ij) = -sqrt_tau_o_sigma(ia) * (j_xp+1._xp)*SQRT(p_xp) * k_mB
+ ynapp1jm1(ia,ip,ij) = +sqrt_tau_o_sigma(ia) * j_xp*SQRT(p_xp+1._xp) * k_mB
+ ynapm1jm1(ia,ip,ij) = +sqrt_tau_o_sigma(ia) * j_xp*SQRT(p_xp) * k_mB
! Trapping terms
- zNapm1j (ia,ip,ij) = +sqrt_tau_o_sigma(ia) *(2._xp*j_xp+1._xp)*SQRT(p_xp)
- zNapm1jp1(ia,ip,ij) = -sqrt_tau_o_sigma(ia) * (j_xp+1._xp)*SQRT(p_xp)
- zNapm1jm1(ia,ip,ij) = -sqrt_tau_o_sigma(ia) * j_xp*SQRT(p_xp)
+ zNapm1j (ia,ip,ij) = +sqrt_tau_o_sigma(ia) *(2._xp*j_xp+1._xp)*SQRT(p_xp) * k_tB
+ zNapm1jp1(ia,ip,ij) = -sqrt_tau_o_sigma(ia) * (j_xp+1._xp)*SQRT(p_xp) * k_tB
+ zNapm1jm1(ia,ip,ij) = -sqrt_tau_o_sigma(ia) * j_xp*SQRT(p_xp) * k_tB
ENDDO
ENDDO
DO ip = 1,local_np
p_int= parray(ip+ngp/2) ! Hermite degree
p_xp = REAL(p_int,xp) ! REAL of Hermite degree
! Landau damping coefficients (ddz napj term)
- xnapp1j(ia,ip) = sqrt_tau_o_sigma(ia) * SQRT(p_xp+1._xp)
- xnapm1j(ia,ip) = sqrt_tau_o_sigma(ia) * SQRT(p_xp)
+ xnapp1j(ia,ip) = sqrt_tau_o_sigma(ia) * SQRT(p_xp+1._xp) * k_ldB
+ xnapm1j(ia,ip) = sqrt_tau_o_sigma(ia) * SQRT(p_xp) * k_ldB
! Magnetic curvature term
- xnapp2j(ia,ip) = tau(ia)/q(ia) * k_cB * SQRT((p_xp+1._xp)*(p_xp + 2._xp))
- xnapm2j(ia,ip) = tau(ia)/q(ia) * k_cB * SQRT( p_xp *(p_xp - 1._xp))
+ xnapp2j(ia,ip) = tau(ia)/q(ia) * SQRT((p_xp+1._xp)*(p_xp + 2._xp)) * k_cB
+ xnapm2j(ia,ip) = tau(ia)/q(ia) * SQRT( p_xp *(p_xp - 1._xp)) * k_cB
ENDDO
DO ij = 1,local_nj
j_int= jarray(ij+ngj/2) ! Laguerre degree
j_xp = REAL(j_int,xp) ! REAL of Laguerre degree
! Magnetic gradient term
- xnapjp1(ia,ij) = -tau(ia)/q(ia) * k_gB * (j_xp + 1._xp)
- xnapjm1(ia,ij) = -tau(ia)/q(ia) * k_gB * j_xp
+ xnapjp1(ia,ij) = -tau(ia)/q(ia) * (j_xp + 1._xp) * k_gB
+ xnapjm1(ia,ij) = -tau(ia)/q(ia) * j_xp * k_gB
ENDDO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! ES linear coefficients for moment RHS !!!!!!!!!!
--
GitLab