Skip to content
Snippets Groups Projects
Commit 2b92585d authored by Antoine Cyril David Hoffmann's avatar Antoine Cyril David Hoffmann :seedling:
Browse files

correction + 2 factor on the heat flux

parent 0f50fdc7
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/array_mod.F90
+ 2
− 0
View file @ 2b92585d
......@@ -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
......
src/memory.F90
+ 1
− 0
View file @ 2b92585d
......@@ -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
......
src/numerics_mod.F90
+ 15
− 1
View file @ 2b92585d
......@@ -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
!******************************************************************************!
......
src/processing_mod.F90
+ 19
− 27
View file @ 2b92585d
......@@ -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
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment