From 2b92585d4ee52b8ee4f5222f6f14ed2b0947fd24 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Fri, 17 Jun 2022 14:39:41 +0200
Subject: [PATCH] correction + 2 factor on the heat flux
---
src/array_mod.F90 | 2 ++
src/memory.F90 | 1 +
src/numerics_mod.F90 | 16 ++++++++++++++-
src/processing_mod.F90 | 46 +++++++++++++++++-------------------------
4 files changed, 37 insertions(+), 28 deletions(-)
diff --git a/src/array_mod.F90 b/src/array_mod.F90
index 8e921396..2db05d5e 100644
--- a/src/array_mod.F90
+++ b/src/array_mod.F90
@@ -60,9 +60,11 @@ MODULE array
! Kernel function evaluation (ij,ikx,iky,iz,odd/even p)
REAL(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: kernel_e
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
+ REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: HF_phi_correction_operator
! 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 afe8aebc..f39ca92d 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -18,6 +18,7 @@ SUBROUTINE memory
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)
+ CALL allocate_array(HF_phi_correction_operator, 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 7b2e7435..e019167a 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -49,7 +49,7 @@ END SUBROUTINE build_dnjs_table
!******************************************************************************!
SUBROUTINE evaluate_kernels
USE basic
- USE array, Only : kernel_e, kernel_i
+ USE array, Only : kernel_e, kernel_i, HF_phi_correction_operator
USE grid
USE model, ONLY : tau_e, tau_i, sigma_e, sigma_i, q_e, q_i, &
lambdaD, CLOS, sigmae2_taue_o2, sigmai2_taui_o2, KIN_E
@@ -85,6 +85,20 @@ ENDDO
ENDDO
ENDDO
ENDDO
+!! Correction term for the evaluation of the heat flux
+HF_phi_correction_operator(ikys:ikye,ikxs:ikxe,izs:ize) = &
+ 2._dp * Kernel_i(1,ikys:ikye,ikxs:ikxe,izs:ize,0) &
+ -1._dp * Kernel_i(2,ikys:ikye,ikxs:ikxe,izs:ize,0)
+
+DO ij = ijgs_i, ijge_i
+ j_int = jarray_i(ij)
+ j_dp = REAL(j_int,dp)
+ HF_phi_correction_operator(ikys:ikye,ikxs:ikxe,izs:ize) = HF_phi_correction_operator(ikys:ikye,ikxs:ikxe,izs:ize) &
+ - Kernel_i(ij,ikys:ikye,ikxs:ikxe,izs:ize,0) * (&
+ 2._dp*(j_dp+1.5_dp) * Kernel_i(ij ,ikys:ikye,ikxs:ikxe,izs:ize,0) &
+ - (j_dp+1.0_dp) * Kernel_i(ij+1,ikys:ikye,ikxs:ikxe,izs:ize,0) &
+ - j_dp * Kernel_i(ij-1,ikys:ikye,ikxs:ikxe,izs:ize,0))
+ENDDO
END SUBROUTINE evaluate_kernels
!******************************************************************************!
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index 048266cb..ef695649 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -92,52 +92,44 @@ END SUBROUTINE compute_radial_ion_transport
! 1D diagnostic to compute the average radial particle transport <T_i v_ExB_x>_xyz
SUBROUTINE compute_radial_heatflux
USE fields, ONLY : moments_i, moments_e, phi
- USE array, ONLY : kernel_e, kernel_i
+ USE array, ONLY : kernel_e, kernel_i, HF_phi_correction_operator
USE geometry, ONLY : Jacobian, iInt_Jacobian
USE time_integration, ONLY : updatetlevel
- USE model, ONLY : qe_taue, qi_taui, KIN_E
+ USE model, ONLY : qe_taue, qi_taui, KIN_E, tau_i, tau_e
USE calculus, ONLY : simpson_rule_z
IMPLICIT NONE
COMPLEX(dp) :: hflux_local, integral
- REAL(dp) :: ky_, buffer(1:2), j_dp
- INTEGER :: i_, world_rank, world_size, root
+ REAL(dp) :: ky_, buffer(1:2), n_dp
+ INTEGER :: i_, world_rank, world_size, root, in
COMPLEX(dp), DIMENSION(izgs:izge) :: integrant ! charge density q_a n_a
- hflux_local = 0._dp ! particle flux
+ hflux_local = 0._dp ! heat flux
IF(ips_i .EQ. 1 .AND. ips_e .EQ. 1) THEN
! Loop to compute gamma_kx = sum_ky sum_j -i k_z Kernel_j Ni00 * phi
DO iky = ikys,ikye
ky_ = kyarray(iky)
DO ikx = ikxs,ikxe
integrant = 0._dp
- DO ij = ijs_i, ije_i
- DO iz = izgs,izge
- integrant(iz) = integrant(iz) + Jacobian(iz,0)*imagu*ky_*CONJG(phi(iky,ikx,iz))&
- *(twothird * ( 2._dp*j_dp * kernel_i(ij ,iky,ikx,iz,0) &
- - (j_dp+1) * kernel_i(ij+1,iky,ikx,iz,0) &
- - j_dp * kernel_i(ij-1,iky,ikx,iz,0))&
- * (moments_i(ip0_i,ij,iky,ikx,iz,updatetlevel)+qi_taui*kernel_i(ij,iky,ikx,iz,0)*phi(iky,ikx,iz)) &
- + SQRT2*onethird * kernel_i(ij,iky,ikx,iz,0) * moments_i(ip2_i,ij,iky,ikx,iz,updatetlevel))
- ENDDO
- ENDDO
- IF(KIN_E) THEN
- DO ij = ijs_e, ije_e
- DO iz = izgs,izge
- integrant(iz) = integrant(iz) + Jacobian(iz,0)*imagu*ky_*CONJG(phi(iky,ikx,iz))&
- *(twothird * ( 2._dp*j_dp * kernel_e(ij ,iky,ikx,iz,0) &
- - (j_dp+1) * kernel_e(ij+1,iky,ikx,iz,0) &
- - j_dp * kernel_e(ij-1,iky,ikx,iz,0))&
- * (moments_e(ip0_e,ij,iky,ikx,iz,updatetlevel)+qe_taue*kernel_e(ij,iky,ikx,iz,0)*phi(iky,ikx,iz)) &
- + SQRT2*onethird * kernel_e(ij,iky,ikx,iz,0) * moments_e(ip2_e,ij,iky,ikx,iz,updatetlevel))
- ENDDO
+ DO in = ijs_i, ije_i
+ n_dp = jarray_i(in)
+
+ integrant(izs:ize) = integrant(izs:ize) + Jacobian(izs:ize,0)*tau_i*imagu*ky_*CONJG(phi(iky,ikx,izs:ize))&
+ *kernel_i(in,iky,ikx,iz,0)*(&
+ 0.5_dp*SQRT2*moments_i(ip2_i,in ,iky,ikx,izs:ize,updatetlevel)&
+ +(2._dp*n_dp + 1.5_dp)*moments_i(ip0_i,in ,iky,ikx,izs:ize,updatetlevel)&
+ -(n_dp+1._dp)*moments_i(ip0_i,in+1,iky,ikx,izs:ize,updatetlevel)&
+ -n_dp*moments_i(ip0_i,in-1,iky,ikx,izs:ize,updatetlevel))
ENDDO
- ENDIF
+ ! Add polarisation contribution
+ integrant(izs:ize) = integrant(izs:ize) + Jacobian(izs:ize,0)*tau_i*imagu*ky_&
+ *CONJG(phi(iky,ikx,izs:ize))*phi(iky,ikx,izs:ize) * HF_phi_correction_operator(iky,ikx,izs:ize)
! Integrate over z
call simpson_rule_z(integrant,integral)
hflux_local = hflux_local + integral*iInt_Jacobian
ENDDO
ENDDO
- buffer(2) = REAL(hflux_local)
+ ! Double it for taking into account the other half plane
+ buffer(2) = 2._dp*REAL(hflux_local)
root = 0
!Gather manually among the rank_p=0 processes and perform the sum
hflux_x = 0
--
GitLab