From 1115f416c1f3654ad3edbac61adec3684a05c502 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Fri, 8 Dec 2023 14:47:38 +0100
Subject: [PATCH] enable adiabatic ion model
---
src/numerics_mod.F90 | 17 +++++------------
src/solve_EM_fields.F90 | 11 ++++++++++-
2 files changed, 15 insertions(+), 13 deletions(-)
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index 3a1e4a98..53d48b6d 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -162,14 +162,12 @@ 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, ADIAB_I, tau_i
+ USE model, ONLY : ADIAB_E, ADIAB_I, tau_i, q_i
USE prec_const, ONLY: xp
IMPLICIT NONE
REAL(xp) :: pol_tot, operator_ion
INTEGER :: in,ikx,iky,iz,ia
- REAL(xp) :: sumker, q_i, sigma_i
- q_i = 1._xp ! we assume single charge ions for adiab. ion model
- sigma_i = 1._xp ! trivial singe sigma_a = sqrt(m_a/m_i)
+ REAL(xp) :: sumker
! This term has no staggered grid dependence. It is evalued for the
! even z grid since poisson uses p=0 moments and phi only.
kxloop: DO ikx = 1,local_nkx
@@ -194,13 +192,8 @@ SUBROUTINE evaluate_poisson_op
IF(ADIAB_E) & ! Adiabatic electron model
pol_tot = pol_tot + 1._xp
- IF(ADIAB_I) THEN ! adiabatic ions model (add ions polarisation)
- sumker = 0._xp ! sum of ion polarisation term
- DO in=1,local_nj
- sumker = sumker + q_i**2/tau_i * kernel_i(in+ngj/2,iky,ikx,iz+ngz/2,ieven)**2 ! ... sum recursively ...
- END DO
- pol_tot = pol_tot + q_i**2/tau_i - sumker
- ENDIF
+ IF(ADIAB_I) & ! adiabatic ions model, kernel_i = 0 and -q_i/tau_i*phi = rho_i
+ pol_tot = pol_tot + q_i**2/tau_i
inv_poisson_op(iky, ikx, iz) = 1._xp/pol_tot
inv_pol_ion (iky, ikx, iz) = 1._xp/operator_ion
@@ -246,7 +239,7 @@ SUBROUTINE evaluate_ampere_op
END DO a
IF(ADIAB_I) THEN ! Add ion contribution on the polarisation
DO in=1,total_nj
- sum_jpol = sum_jpol + q_i**2/(sigma_i**2)*kernel_i(in+ngj/2,iky,ikx,iz+ngz/2,iodd)**2 ! ... sum recursively ...
+ sum_jpol = sum_jpol + 0*q_i**2/(sigma_i**2)*kernel_i(in+ngj/2,iky,ikx,iz+ngz/2,iodd)**2 ! ... sum recursively ...
END DO
ENDIF
operator = 2._xp*kperp2*hatB(iz+ngz/2,iodd)**2 + beta*sum_jpol
diff --git a/src/solve_EM_fields.F90 b/src/solve_EM_fields.F90
index 2c1e20b7..d28a4aca 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 , tau_i
+ USE model, ONLY: lambdaD, ADIAB_E, ADIAB_I, q_i, tau_i
use species, ONLY: q
USE processing, ONLY: compute_density
USE geometry, ONLY: iInt_Jacobian, Jacobian
@@ -66,6 +66,15 @@ CONTAINS
ENDIF
rho = rho + fsa_phi
ENDIF
+ !!!!!!!!!!!!!!! adiabatic ion model
+ ! Candy et al. 2007, rho_i = -q_i/tau_i phi
+ IF (ADIAB_I) THEN
+ DO iz = 1,local_nz
+ izi = iz+ngz/2
+ rho(iz) = rho(iz) - 0*q_i/tau_i * phi(iky,ikx,izi)
+ ENDDO
+ ENDIF
+
!!!!!!!!!!!!!!! Inverting the poisson equation
DO iz = 1,local_nz
izi = iz+ngz/2
--
GitLab