From cd641add21fc0f94dfd5dec8824a871266ae7238 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Wed, 8 Jun 2022 11:07:42 +0200
Subject: [PATCH] important bug correction in poisson equation for adiabatic
electron
---
src/array_mod.F90 | 1 +
src/memory.F90 | 1 +
src/numerics_mod.F90 | 5 +++--
src/poisson.F90 | 43 +++++++++++++++++++++----------------------
4 files changed, 26 insertions(+), 24 deletions(-)
diff --git a/src/array_mod.F90 b/src/array_mod.F90
index 3dca4e30..8e921396 100644
--- a/src/array_mod.F90
+++ b/src/array_mod.F90
@@ -62,6 +62,7 @@ MODULE array
REAL(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: kernel_i
! Poisson operator (ikx,iky,iz)
REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: inv_poisson_op
+ REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: inv_pol_ion
! Gyrocenter density for electron and ions (ikx,iky,iz)
COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: Ne00
diff --git a/src/memory.F90 b/src/memory.F90
index c9509278..afe8aebc 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -17,6 +17,7 @@ SUBROUTINE memory
CALL allocate_array( phi_ZF, ikxs,ikxe, izs,ize)
CALL allocate_array( phi_EM, ikys,ikye, izs,ize)
CALL allocate_array(inv_poisson_op, ikys,ikye, ikxs,ikxe, izs,ize)
+ CALL allocate_array( inv_pol_ion, ikys,ikye, ikxs,ikxe, izs,ize)
!Electrons arrays
IF(KIN_E) THEN
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index 77d0d357..7b2e7435 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -94,7 +94,7 @@ END SUBROUTINE evaluate_kernels
!******************************************************************************!
SUBROUTINE evaluate_poisson_op
USE basic
- USE array, Only : kernel_e, kernel_i, inv_poisson_op
+ USE array, Only : kernel_e, kernel_i, inv_poisson_op, inv_pol_ion
USE grid
USE model, ONLY : tau_e, tau_i, q_e, q_i, KIN_E
IMPLICIT NONE
@@ -125,9 +125,10 @@ SUBROUTINE evaluate_poisson_op
END DO
! Adiabatic model
ELSE
- pol_e = 1._dp - qe2_taue
+ pol_e = qe2_taue - 1._dp
ENDIF
inv_poisson_op(iky, ikx, iz) = 1._dp/(qe2_taue + qi2_taui - pol_i - pol_e)
+ inv_pol_ion (iky, ikx, iz) = 1._dp/(qi2_taui - pol_i)
ENDIF
END DO zloop
END DO kyloop
diff --git a/src/poisson.F90 b/src/poisson.F90
index 75e2b887..49f1db40 100644
--- a/src/poisson.F90
+++ b/src/poisson.F90
@@ -15,9 +15,9 @@ SUBROUTINE poisson
IMPLICIT NONE
INTEGER :: ini,ine, i_, root_bcast
- COMPLEX(dp) :: fa_phi, intf_ ! current flux averaged phi
+ COMPLEX(dp) :: fsa_phi, intf_ ! current flux averaged phi
INTEGER :: count !! mpi integer to broadcast the electric potential at the end
- COMPLEX(dp), DIMENSION(izs:ize) :: rho_i, rho_e ! charge density q_a n_a
+ COMPLEX(dp), DIMENSION(izs:ize) :: rho_i, rho_e, integrant ! charge density q_a n_a and aux var
! Execution time start
CALL cpu_time(t0_poisson)
@@ -27,36 +27,35 @@ SUBROUTINE poisson
kxloop: DO ikx = ikxs,ikxe
kyloop: DO iky = ikys,ikye
- !!!! Compute ion gyro charge density
+ !!!! Compute ion particle charge density q_i n_i
rho_i = 0._dp
DO ini=1,jmaxi+1
- DO iz = izs,ize
- rho_i(iz) = rho_i(iz) &
- +q_i*kernel_i(ini,iky,ikx,iz,0)*moments_i(ip0_i,ini,iky,ikx,iz,updatetlevel)
- ENDDO
+ rho_i(izs:ize) = rho_i(izs:ize) &
+ +q_i*kernel_i(ini,iky,ikx,izs:ize,0)*moments_i(ip0_i,ini,iky,ikx,izs:ize,updatetlevel)
END DO
- !!!! Compute electron gyro charge density
+ !!!! Compute electron particle charge density q_e n_e
rho_e = 0._dp
IF (KIN_E) THEN ! Kinetic electrons
DO ine=1,jmaxe+1
- DO iz = izs,ize
- rho_e(iz) = rho_e(iz) &
- +q_e*kernel_e(ine,iky,ikx,iz,0)*moments_e(ip0_e,ine,iky,ikx,iz,updatetlevel)
- ENDDO
+ rho_e(izs:ize) = rho_e(izs:ize) &
+ +q_e*kernel_e(ine,iky,ikx,izs:ize,0)*moments_e(ip0_e,ine,iky,ikx,izs:ize,updatetlevel)
END DO
ELSE ! Adiabatic electrons
- ! Adiabatic charge density (linked to flux averaged phi)
- fa_phi = 0._dp
- IF(kyarray(iky).EQ.0._dp) THEN ! take ky=0 mode (average)
- DO ini=1,jmaxi+1 ! some over FLR contributions
- rho_e(izs:ize) = Jacobian(izs:ize,0)*moments_i(ip0_i,ini,iky,ikx,izs:ize,updatetlevel)&
- *kernel_i(ini,iky,ikx,izs:ize,0)*(inv_poisson_op(iky,ikx,izs:ize)-1._dp)
- call simpson_rule_z(rho_e(izs:ize),intf_) ! integrate over z
- fa_phi = fa_phi + intf_
- ENDDO
- rho_e(izs:ize) = fa_phi*iInt_Jacobian !Normalize by 1/int(Jxyz)dz
+ ! Adiabatic charge density (linked to flux surface averaged phi)
+ fsa_phi = 0._dp
+ ! We compute the flux surface average solving a flux surface averaged
+ ! Poisson equation, i.e.
+ !
+ ! [qi^2(1-sum_j K_j^2)/tau_i] <phi>_psi = <q_i n_i >_psi
+ ! inv_pol_ion^-1 fsa_phi = simpson(Jacobian rho_i ) * iInt_Jacobian
+ IF(kyarray(iky).EQ.0._dp) THEN ! take ky=0 mode (y-average)
+ ! Prepare integrant for z-average
+ integrant(izs:ize) = Jacobian(izs:ize,0)*rho_i(izs:ize)*inv_pol_ion(iky,ikx,izs:ize)
+ call simpson_rule_z(integrant(izs:ize),intf_) ! get the flux averaged phi
+ fsa_phi = intf_
ENDIF
+ rho_e(izs:ize) = fsa_phi * iInt_Jacobian !Normalize by 1/int(Jxyz)dz
ENDIF
!!!!!!!!!!!!!!! Inverting the poisson equation !!!!!!!!!!!!!!!!!!!!!!!!!!
--
GitLab