From 43d8c11f6d81fc3c8babf65cb612e31851e5a96a Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Thu, 15 Jun 2023 15:43:15 +0200
Subject: [PATCH] adiabatic ion model (not tested)
---
src/model_mod.F90 | 23 +++++++++++++++--------
src/numerics_mod.F90 | 15 +++++++++------
src/solve_EM_fields.F90 | 5 +----
src/species_mod.F90 | 4 +++-
4 files changed, 28 insertions(+), 19 deletions(-)
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index fffedaba..4364d580 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -23,7 +23,8 @@ MODULE model
REAL(xp), PUBLIC, PROTECTED :: lambdaD = 0._xp ! Debye length
REAL(xp), PUBLIC, PROTECTED :: beta = 0._xp ! electron plasma Beta (8piNT_e/B0^2)
LOGICAL, PUBLIC :: ADIAB_E = .false. ! adiabatic electron model
- REAL(xp), PUBLIC, PROTECTED :: tau_e = 1.0 ! electron temperature ratio for adiabatic electrons
+ LOGICAL, PUBLIC :: ADIAB_I = .false. ! adiabatic ion model
+ REAL(xp), PUBLIC, PROTECTED :: tau_i = 1.0 ! electron-ion temperature ratio for ion adiabatic model
! Auxiliary variable
LOGICAL, PUBLIC, PROTECTED :: EM = .false. ! Electromagnetic effects flag
LOGICAL, PUBLIC, PROTECTED :: MHD_PD = .false. ! MHD pressure drift
@@ -43,18 +44,23 @@ CONTAINS
NAMELIST /MODEL_PAR/ KERN, LINEARITY, RM_LD_T_EQ, &
mu_x, mu_y, N_HD, HDz_h, mu_z, mu_p, mu_j, HYP_V, Na,&
- nu, k_gB, k_cB, lambdaD, MHD_PD, beta, ADIAB_E, tau_e
+ nu, k_gB, k_cB, lambdaD, MHD_PD, beta, ADIAB_E, ADIAB_I, tau_i
READ(lu_in,model_par)
- IF((HYP_V .EQ. 'dvpar4') .AND. (num_procs_p .GT. 1)) THEN
+ IF (ADIAB_E .AND. ADIAB_I) &
+ ERROR STOP '>> ERROR << cannot have both adiab e and adiab i models'
+
+ IF((HYP_V .EQ. 'dvpar4') .AND. (num_procs_p .GT. 1)) &
ERROR STOP '>> ERROR << dvpar4 velocity dissipation is not compatible with current p parallelization'
- ENDIF
IF(Na .EQ. 1) THEN
- IF(.NOT. ADIAB_E) ERROR STOP "With one species, ADIAB_E must be set to .true. STOP"
- CALL speak('Adiabatic electron model -> beta = 0')
- beta = 0._xp
+ IF((.NOT. ADIAB_E) .AND. (.NOT. ADIAB_I)) ERROR STOP "With one species, ADIAB_E or ADIAB_I must be set to .true. STOP"
+ IF(ADIAB_E) THEN
+ CALL speak('Adiabatic electron model -> beta = 0')
+ beta = 0._xp
+ ENDIF
+ IF(ADIAB_I) CALL speak('Adiabatic ion model -> beta = 0')
ENDIF
IF(beta .GT. 0) THEN
@@ -91,7 +97,8 @@ CONTAINS
CALL attach(fid, TRIM(str), "MHD_PD", MHD_PD)
CALL attach(fid, TRIM(str), "beta", beta)
CALL attach(fid, TRIM(str), "ADIAB_E", ADIAB_E)
- CALL attach(fid, TRIM(str), "tau_e", tau_e)
+ CALL attach(fid, TRIM(str), "ADIAB_I", ADIAB_I)
+ CALL attach(fid, TRIM(str), "tau_i", tau_i)
END SUBROUTINE model_outputinputs
END MODULE model
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index 9111ea3f..e6d4ac94 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -134,10 +134,10 @@ SUBROUTINE evaluate_poisson_op
USE grid, ONLY : local_na, local_nkx, local_nky, local_nz,&
kxarray, kyarray, local_nj, ngj, ngz, ieven
USE species, ONLY : q2_tau
- USE model, ONLY : ADIAB_E
+ USE model, ONLY : ADIAB_E, ADIAB_I, tau_i
USE prec_const, ONLY: xp
IMPLICIT NONE
- REAL(xp) :: pol_tot, operator, operator_ion ! (Z^2/tau (1-sum_n kernel_na^2))
+ REAL(xp) :: pol_tot, operator_ion ! (Z^2/tau (1-sum_n kernel_na^2))
INTEGER :: in,ikx,iky,iz,ia
REAL(xp) :: sumker ! (Z_a^2/tau_a (1-sum_n kernel_na^2))
@@ -149,7 +149,7 @@ SUBROUTINE evaluate_poisson_op
IF( (kxarray(ikx).EQ.0._xp) .AND. (kyarray(iky).EQ.0._xp) ) THEN
inv_poisson_op(iky, ikx, iz) = 0._xp
inv_pol_ion (iky, ikx, iz) = 0._xp
-ELSE
+ ELSE
! loop over n only up to the max polynomial degree
pol_tot = 0._xp ! total polarisation term
a:DO ia = 1,local_na ! sum over species
@@ -161,12 +161,15 @@ ELSE
pol_tot = pol_tot + q2_tau(ia) - sumker
ENDDO a
operator_ion = pol_tot
- IF(ADIAB_E) THEN ! Adiabatic electron model
+
+ IF(ADIAB_E) & ! Adiabatic electron model
pol_tot = pol_tot + 1._xp
- ENDIF
- operator = pol_tot
+
inv_poisson_op(iky, ikx, iz) = 1._xp/pol_tot
inv_pol_ion (iky, ikx, iz) = 1._xp/operator_ion
+
+ IF(ADIAB_I) & ! adiabatic ions
+ inv_poisson_op(iky, ikx, iz) = tau_i
ENDIF
END DO zloop
END DO kyloop
diff --git a/src/solve_EM_fields.F90 b/src/solve_EM_fields.F90
index 871cc185..099bfad6 100644
--- a/src/solve_EM_fields.F90
+++ b/src/solve_EM_fields.F90
@@ -22,7 +22,7 @@ CONTAINS
ip0, total_nj, ngj
USE calculus, ONLY: simpson_rule_z
USE parallel, ONLY: manual_3D_bcast
- USE model, ONLY: lambdaD, ADIAB_E
+ USE model, ONLY: lambdaD, ADIAB_E, ADIAB_I, tau_i
use species, ONLY: q
USE processing, ONLY: compute_density
USE geometry, ONLY: iInt_Jacobian, Jacobian
@@ -67,9 +67,6 @@ CONTAINS
ENDIF
rho = rho + fsa_phi
ENDIF
- !!!!!!!!!!!!!!! adiabatic ions ?
- ! IF (ADIAB_I) THEN
- ! ENDIF
!!!!!!!!!!!!!!! Inverting the poisson equation
DO iz = 1,local_nz
phi(iky,ikx,iz+ngz/2) = inv_poisson_op(iky,ikx,iz)*rho(iz)
diff --git a/src/species_mod.F90 b/src/species_mod.F90
index 088f4019..f5635ecb 100644
--- a/src/species_mod.F90
+++ b/src/species_mod.F90
@@ -37,7 +37,7 @@ CONTAINS
SUBROUTINE species_readinputs
! Read the input parameters
USE basic, ONLY : lu_in
- USE model, ONLY : Na, nu, ADIAB_E, MHD_PD
+ USE model, ONLY : Na, nu, ADIAB_E, ADIAB_I, MHD_PD
USE prec_const
IMPLICIT NONE
INTEGER :: ia,ib
@@ -80,6 +80,8 @@ CONTAINS
SELECT CASE (name_)
CASE ('electrons','e','electron')
ADIAB_E = .FALSE.
+ CASE ('ions','i','ion')
+ ADIAB_I = .FALSE.
END SELECT
ENDDO
IF (.NOT. MHD_PD) Ptot = 0._dp
--
GitLab