From 96062c8b894d26aa7bf9645a463302a0d7dc5e13 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Wed, 21 Sep 2022 17:44:17 +0200
Subject: [PATCH] EM for transport
---
src/processing_mod.F90 | 648 ++++++++++++++++++++++++-----------------
1 file changed, 380 insertions(+), 268 deletions(-)
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index 5b64b26e..291608c5 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -17,279 +17,389 @@ CONTAINS
! 1D diagnostic to compute the average radial particle transport <n_i v_ExB_x>_xyz
SUBROUTINE compute_radial_ion_transport
- USE fields, ONLY : moments_i, phi
- USE array, ONLY : kernel_i
- USE geometry, ONLY : Jacobian, iInt_Jacobian
- USE time_integration, ONLY : updatetlevel
- USE calculus, ONLY : simpson_rule_z
- IMPLICIT NONE
- COMPLEX(dp) :: pflux_local, gflux_local, integral
- REAL(dp) :: ky_, buffer(1:2)
- INTEGER :: i_, world_rank, world_size, root
- COMPLEX(dp), DIMENSION(izgs:izge) :: integrant
-
- pflux_local = 0._dp ! particle flux
- gflux_local = 0._dp ! gyrocenter flux
- integrant = 0._dp ! auxiliary variable for z integration
- IF(ips_i .EQ. 1) THEN
- !! Gyro center flux (drift kinetic)
- DO iky = ikys,ikye
- ky_ = kyarray(iky)
- DO ikx = ikxs,ikxe
- DO iz = izgs,izge
- integrant(iz) = Jacobian(iz,0)*imagu * ky_ * moments_i(ip0_i,1,iky,ikx,iz,updatetlevel) * CONJG(phi(iky,ikx,iz))
- ENDDO
- ! Integrate over z
- call simpson_rule_z(integrant,integral)
- ! add contribution
- gflux_local = gflux_local + integral*iInt_Jacobian
- ENDDO
+ USE fields, ONLY : moments_i, phi, psi
+ USE array, ONLY : kernel_i
+ USE geometry, ONLY : Jacobian, iInt_Jacobian
+ USE time_integration, ONLY : updatetlevel
+ USE calculus, ONLY : simpson_rule_z
+ USE model, ONLY : sqrt_tau_o_sigma_i, EM
+ IMPLICIT NONE
+ COMPLEX(dp) :: pflux_local, gflux_local, integral
+ REAL(dp) :: ky_, buffer(1:2)
+ INTEGER :: i_, world_rank, world_size, root
+ COMPLEX(dp), DIMENSION(izgs:izge) :: integrant
+
+ pflux_local = 0._dp ! particle flux
+ gflux_local = 0._dp ! gyrocenter flux
+ integrant = 0._dp ! auxiliary variable for z integration
+ !!---------- Gyro center flux (drift kinetic) ------------
+ ! Electrostatic part
+ IF(CONTAINS_ip0_i) THEN
+ DO iky = ikys,ikye
+ ky_ = kyarray(iky)
+ DO ikx = ikxs,ikxe
+ integrant(izgs:izge) = integrant(izgs:izge) &
+ +moments_i(ip0_i,ij0_i,iky,ikx,izgs:izge,updatetlevel)&
+ *imagu*ky_*CONJG(phi(iky,ikx,izgs:izge))
ENDDO
- !! Particle flux
- DO iky = ikys,ikye
- ky_ = kyarray(iky)
- DO ikx = ikxs,ikxe
- integrant = 0._dp ! auxiliary variable for z integration
- DO ij = ijs_i, ije_i
- DO iz = izgs,izge
- integrant(iz) = integrant(iz) + &
- Jacobian(iz,0)*imagu * ky_ * kernel_i(ij,iky,ikx,iz,0) &
- *moments_i(ip0_i,ij,iky,ikx,iz,updatetlevel) * CONJG(phi(iky,ikx,iz))
- ENDDO
- ENDDO
- ! Integrate over z
- call simpson_rule_z(integrant,integral)
- ! add contribution
- pflux_local = pflux_local + integral*iInt_Jacobian
+ ENDDO
+ ENDIF
+ ! Electromagnetic part
+ IF( EM .AND. CONTAINS_ip1_i ) THEN
+ DO iky = ikys,ikye
+ ky_ = kyarray(iky)
+ DO ikx = ikxs,ikxe
+ integrant(izgs:izge) = integrant(izgs:izge)&
+ -sqrt_tau_o_sigma_i*moments_i(ip1_i,ij0_i,iky,ikx,izgs:izge,updatetlevel)&
+ *imagu*ky_*CONJG(psi(iky,ikx,izgs:izge))
+ ENDDO
+ ENDDO
+ ENDIF
+ ! Integrate over z
+ integrant(izgs:izge) = Jacobian(izgs:izge,0)*integrant(izgs:izge)
+ call simpson_rule_z(integrant,integral)
+ ! Get process local gyrocenter flux
+ gflux_local = integral*iInt_Jacobian
+
+ !
+ integrant = 0._dp ! reset auxiliary variable
+ !!---------- Particle flux (gyrokinetic) ------------
+ ! Electrostatic part
+ IF(CONTAINS_ip0_i) THEN
+ DO iky = ikys,ikye
+ ky_ = kyarray(iky)
+ DO ikx = ikxs,ikxe
+ DO ij = ijs_i, ije_i
+ integrant(izgs:izge) = integrant(izgs:izge)&
+ +moments_i(ip0_i,ij,iky,ikx,izgs:izge,updatetlevel)&
+ *imagu*ky_*kernel_i(ij,iky,ikx,izgs:izge,0)*CONJG(phi(iky,ikx,izgs:izge))
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+ ! Electromagnetic part
+ IF( EM .AND. CONTAINS_ip1_i ) THEN
+ DO iky = ikys,ikye
+ ky_ = kyarray(iky)
+ DO ikx = ikxs,ikxe
+ integrant = 0._dp ! auxiliary variable for z integration
+ DO ij = ijs_i, ije_i
+ integrant(izgs:izge) = integrant(izgs:izge)&
+ -sqrt_tau_o_sigma_i*moments_i(ip1_i,ij,iky,ikx,izgs:izge,updatetlevel)&
+ *imagu*ky_*kernel_i(ij,iky,ikx,izgs:izge,0)*CONJG(psi(iky,ikx,izgs:izge))
ENDDO
ENDDO
- ENDIF
- buffer(1) = REAL(gflux_local)
- buffer(2) = REAL(pflux_local)
- root = 0
- !Gather manually among the rank_p=0 processes and perform the sum
- gflux_ri = 0
- pflux_ri = 0
- IF (num_procs_ky .GT. 1) THEN
- !! Everyone sends its local_sum to root = 0
- IF (rank_ky .NE. root) THEN
- CALL MPI_SEND(buffer, 2 , MPI_DOUBLE_PRECISION, root, 1234, comm_ky, ierr)
- ELSE
- ! Recieve from all the other processes
- DO i_ = 0,num_procs_ky-1
- IF (i_ .NE. rank_ky) &
- CALL MPI_RECV(buffer, 2 , MPI_DOUBLE_PRECISION, i_, 1234, comm_ky, MPI_STATUS_IGNORE, ierr)
- gflux_ri = gflux_ri + buffer(1)
- pflux_ri = pflux_ri + buffer(2)
- ENDDO
- ENDIF
- ELSE
- gflux_ri = gflux_local
- pflux_ri = pflux_local
- ENDIF
- ! if(my_id .eq. 0) write(*,*) 'pflux_ri = ',pflux_ri
+ ENDDO
+ ENDIF
+ ! Integrate over z
+ integrant(izgs:izge) = Jacobian(izgs:izge,0)*integrant(izgs:izge)
+ call simpson_rule_z(integrant,integral)
+ ! Get process local particle flux
+ pflux_local = integral*iInt_Jacobian
+
+ !!!!---------- Sum over all processes ----------
+ buffer(1) = 2._dp*REAL(gflux_local)
+ buffer(2) = 2._dp*REAL(pflux_local)
+ root = 0
+ !Gather manually among the rank_p=0 processes and perform the sum
+ gflux_ri = 0
+ pflux_ri = 0
+ IF (num_procs_ky .GT. 1) THEN
+ !! Everyone sends its local_sum to root = 0
+ IF (rank_ky .NE. root) THEN
+ CALL MPI_SEND(buffer, 2 , MPI_DOUBLE_PRECISION, root, 1234, comm_ky, ierr)
+ ELSE
+ ! Recieve from all the other processes
+ DO i_ = 0,num_procs_ky-1
+ IF (i_ .NE. rank_ky) &
+ CALL MPI_RECV(buffer, 2 , MPI_DOUBLE_PRECISION, i_, 1234, comm_ky, MPI_STATUS_IGNORE, ierr)
+ gflux_ri = gflux_ri + buffer(1)
+ pflux_ri = pflux_ri + buffer(2)
+ ENDDO
+ ENDIF
+ ELSE
+ gflux_ri = gflux_local
+ pflux_ri = pflux_local
+ ENDIF
+ ! if(my_id .eq. 0) write(*,*) 'pflux_ri = ',pflux_ri
END SUBROUTINE compute_radial_ion_transport
! 1D diagnostic to compute the average radial particle transport <n_e v_ExB_x>_xyz
SUBROUTINE compute_radial_electron_transport
- USE fields, ONLY : moments_e, phi
- USE array, ONLY : kernel_e
- USE geometry, ONLY : Jacobian, iInt_Jacobian
- USE time_integration, ONLY : updatetlevel
- USE calculus, ONLY : simpson_rule_z
- IMPLICIT NONE
- COMPLEX(dp) :: pflux_local, gflux_local, integral
- REAL(dp) :: ky_, buffer(1:2)
- INTEGER :: i_, world_rank, world_size, root
- COMPLEX(dp), DIMENSION(izgs:izge) :: integrant
-
- pflux_local = 0._dp ! particle flux
- gflux_local = 0._dp ! gyrocenter flux
- integrant = 0._dp ! auxiliary variable for z integration
- IF(ips_e .EQ. 1) THEN
- !! Gyro center flux (drift kinetic)
- DO iky = ikys,ikye
- ky_ = kyarray(iky)
- DO ikx = ikxs,ikxe
- DO iz = izgs,izge
- integrant(iz) = Jacobian(iz,0)*imagu * ky_ * moments_e(ip0_e,1,iky,ikx,iz,updatetlevel) * CONJG(phi(iky,ikx,iz))
- ENDDO
- ! Integrate over z
- call simpson_rule_z(integrant,integral)
- ! add contribution
- gflux_local = gflux_local + integral*iInt_Jacobian
- ENDDO
+ USE fields, ONLY : moments_e, phi, psi
+ USE array, ONLY : kernel_e
+ USE geometry, ONLY : Jacobian, iInt_Jacobian
+ USE time_integration, ONLY : updatetlevel
+ USE calculus, ONLY : simpson_rule_z
+ USE model, ONLY : sqrt_tau_o_sigma_e, EM
+ IMPLICIT NONE
+ COMPLEX(dp) :: pflux_local, gflux_local, integral
+ REAL(dp) :: ky_, buffer(1:2)
+ INTEGER :: i_, world_rank, world_size, root
+ COMPLEX(dp), DIMENSION(izgs:izge) :: integrant
+
+ pflux_local = 0._dp ! particle flux
+ gflux_local = 0._dp ! gyrocenter flux
+ integrant = 0._dp ! auxiliary variable for z integration
+ !!---------- Gyro center flux (drift kinetic) ------------
+ ! Electrostatic part
+ IF(CONTAINS_ip0_e) THEN
+ DO iky = ikys,ikye
+ ky_ = kyarray(iky)
+ DO ikx = ikxs,ikxe
+ integrant(izgs:izge) = integrant(izgs:izge) &
+ +moments_e(ip0_e,ij0_e,iky,ikx,izgs:izge,updatetlevel)&
+ *imagu*ky_*CONJG(phi(iky,ikx,izgs:izge))
ENDDO
- !! Particle flux
- DO iky = ikys,ikye
- ky_ = kyarray(iky)
- DO ikx = ikxs,ikxe
- integrant = 0._dp ! auxiliary variable for z integration
- DO ij = ijs_e, ije_e
- DO iz = izgs,izge
- integrant(iz) = integrant(iz) + &
- Jacobian(iz,0)*imagu * ky_ * kernel_e(ij,iky,ikx,iz,0) &
- *moments_e(ip0_e,ij,iky,ikx,iz,updatetlevel) * CONJG(phi(iky,ikx,iz))
- ENDDO
- ENDDO
- ! Integrate over z
- call simpson_rule_z(integrant,integral)
- ! add contribution
- pflux_local = pflux_local + integral*iInt_Jacobian
+ ENDDO
+ ENDIF
+ ! Electromagnetic part
+ IF( EM .AND. CONTAINS_ip1_e ) THEN
+ DO iky = ikys,ikye
+ ky_ = kyarray(iky)
+ DO ikx = ikxs,ikxe
+ integrant(izgs:izge) = integrant(izgs:izge)&
+ -sqrt_tau_o_sigma_e*moments_e(ip1_e,ij0_e,iky,ikx,izgs:izge,updatetlevel)&
+ *imagu*ky_*CONJG(psi(iky,ikx,izgs:izge))
+ ENDDO
+ ENDDO
+ ENDIF
+ ! Integrate over z
+ integrant(izgs:izge) = Jacobian(izgs:izge,0)*integrant(izgs:izge)
+ call simpson_rule_z(integrant,integral)
+ ! Get process local gyrocenter flux
+ gflux_local = integral*iInt_Jacobian
+ !
+ integrant = 0._dp ! reset auxiliary variable
+ !!---------- Particle flux (gyrokinetic) ------------
+ ! Electrostatic part
+ IF(CONTAINS_ip0_e) THEN
+ DO iky = ikys,ikye
+ ky_ = kyarray(iky)
+ DO ikx = ikxs,ikxe
+ DO ij = ijs_e, ije_e
+ integrant(izgs:izge) = integrant(izgs:izge)&
+ +moments_e(ip0_e,ij,iky,ikx,izgs:izge,updatetlevel)&
+ *imagu*ky_*kernel_e(ij,iky,ikx,izgs:izge,0)*CONJG(phi(iky,ikx,izgs:izge))
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+ ! Electromagnetic part
+ IF( EM .AND. CONTAINS_ip1_e ) THEN
+ DO iky = ikys,ikye
+ ky_ = kyarray(iky)
+ DO ikx = ikxs,ikxe
+ integrant = 0._dp ! auxiliary variable for z integration
+ DO ij = ijs_e, ije_e
+ integrant(izgs:izge) = integrant(izgs:izge)&
+ -sqrt_tau_o_sigma_e*moments_e(ip1_e,ij,iky,ikx,izgs:izge,updatetlevel)&
+ *imagu*ky_*kernel_e(ij,iky,ikx,izgs:izge,0)*CONJG(psi(iky,ikx,izgs:izge))
ENDDO
ENDDO
- ENDIF
- buffer(1) = REAL(gflux_local)
- buffer(2) = REAL(pflux_local)
- root = 0
- !Gather manually among the rank_p=0 processes and perform the sum
- gflux_re = 0
- pflux_re = 0
- IF (num_procs_ky .GT. 1) THEN
- !! Everyone sends its local_sum to root = 0
- IF (rank_ky .NE. root) THEN
- CALL MPI_SEND(buffer, 2 , MPI_DOUBLE_PRECISION, root, 1234, comm_ky, ierr)
- ELSE
- ! Recieve from all the other processes
- DO i_ = 0,num_procs_ky-1
- IF (i_ .NE. rank_ky) &
- CALL MPI_RECV(buffer, 2 , MPI_DOUBLE_PRECISION, i_, 1234, comm_ky, MPI_STATUS_IGNORE, ierr)
- gflux_re = gflux_re + buffer(1)
- pflux_re = pflux_re + buffer(2)
- ENDDO
- ENDIF
- ELSE
- gflux_re = gflux_local
- pflux_re = pflux_local
- ENDIF
+ ENDDO
+ ENDIF
+ ! Integrate over z
+ integrant(izgs:izge) = Jacobian(izgs:izge,0)*integrant(izgs:izge)
+ call simpson_rule_z(integrant,integral)
+ ! Get process local particle flux
+ pflux_local = integral*iInt_Jacobian
+
+ !!!!---------- Sum over all processes ----------
+ buffer(1) = 2._dp*REAL(gflux_local)
+ buffer(2) = 2._dp*REAL(pflux_local)
+ root = 0
+ !Gather manually among the rank_p=0 processes and perform the sum
+ gflux_re = 0
+ pflux_re = 0
+ IF (num_procs_ky .GT. 1) THEN
+ !! Everyone sends its local_sum to root = 0
+ IF (rank_ky .NE. root) THEN
+ CALL MPI_SEND(buffer, 2 , MPI_DOUBLE_PRECISION, root, 1234, comm_ky, ierr)
+ ELSE
+ ! Recieve from all the other processes
+ DO i_ = 0,num_procs_ky-1
+ IF (i_ .NE. rank_ky) &
+ CALL MPI_RECV(buffer, 2 , MPI_DOUBLE_PRECISION, i_, 1234, comm_ky, MPI_STATUS_IGNORE, ierr)
+ gflux_re = gflux_re + buffer(1)
+ pflux_re = pflux_re + buffer(2)
+ ENDDO
+ ENDIF
+ ELSE
+ gflux_re = gflux_local
+ pflux_re = pflux_local
+ ENDIF
END SUBROUTINE compute_radial_electron_transport
! 1D diagnostic to compute the average radial particle transport <T_i v_ExB_x>_xyz
SUBROUTINE compute_radial_ion_heatflux
- USE fields, ONLY : moments_i, phi
- USE array, ONLY : kernel_i, HF_phi_correction_operator
- USE geometry, ONLY : Jacobian, iInt_Jacobian
- USE time_integration, ONLY : updatetlevel
- USE model, ONLY : qi_taui, KIN_E, tau_i
- USE calculus, ONLY : simpson_rule_z
- IMPLICIT NONE
- COMPLEX(dp) :: hflux_local, integral
- 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 ! 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 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,izs:ize,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
- ! 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
+ USE fields, ONLY : moments_i, phi, psi
+ USE array, ONLY : kernel_i, HF_phi_correction_operator
+ USE geometry, ONLY : Jacobian, iInt_Jacobian
+ USE time_integration, ONLY : updatetlevel
+ USE model, ONLY : qi_taui, KIN_E, tau_i, EM
+ USE calculus, ONLY : simpson_rule_z
+ USE model, ONLY : sqrt_tau_o_sigma_i
+ IMPLICIT NONE
+ COMPLEX(dp) :: hflux_local, integral
+ 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 ! heat flux
+ integrant = 0._dp ! z integration auxiliary variable
+ !!----------------ELECTROSTATIC CONTRIBUTION---------------------------
+ IF(CONTAINS_ip0_i .AND. CONTAINS_ip2_i) 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
+ DO in = ijs_i, ije_i
+ n_dp = jarray_i(in)
+ integrant(izgs:izge) = integrant(izgs:izge) + tau_i*imagu*ky_*CONJG(phi(iky,ikx,izgs:izge))&
+ *kernel_i(in,iky,ikx,izgs:izge,0)*(&
+ 0.5_dp*SQRT2*moments_i(ip2_i,in ,iky,ikx,izgs:izge,updatetlevel)&
+ +(2._dp*n_dp + 1.5_dp)*moments_i(ip0_i,in ,iky,ikx,izgs:izge,updatetlevel)&
+ -(n_dp+1._dp)*moments_i(ip0_i,in+1,iky,ikx,izgs:izge,updatetlevel)&
+ -n_dp*moments_i(ip0_i,in-1,iky,ikx,izgs:izge,updatetlevel))
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+ IF(EM .AND. CONTAINS_ip1_i .AND. CONTAINS_ip3_i) THEN
+ !!----------------ELECTROMAGNETIC CONTRIBUTION--------------------
+ DO iky = ikys,ikye
+ ky_ = kyarray(iky)
+ DO ikx = ikxs,ikxe
+ DO in = ijs_i, ije_i
+ n_dp = jarray_i(in)
+ integrant(izgs:izge) = integrant(izgs:izge) &
+ +tau_i*sqrt_tau_o_sigma_i*imagu*ky_*CONJG(psi(iky,ikx,izgs:izge))&
+ *kernel_i(in,iky,ikx,izgs:izge,0)*(&
+ 0.5_dp*SQRT2*SQRT3*moments_i(ip3_i,in ,iky,ikx,izgs:izge,updatetlevel)&
+ +1.5_dp*CONJG(moments_i(ip1_i,in ,iky,ikx,izgs:izge,updatetlevel))& !?????
+ +(2._dp*n_dp+1._dp)*moments_i(ip1_i,in ,iky,ikx,izgs:izge,updatetlevel)&
+ -(n_dp+1._dp)*moments_i(ip1_i,in+1,iky,ikx,izgs:izge,updatetlevel)&
+ -n_dp*moments_i(ip1_i,in-1,iky,ikx,izgs:izge,updatetlevel))
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+ ! Add polarisation contribution
+ integrant(izs:ize) = integrant(izs:ize) + 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
+ integrant(izgs:izge) = Jacobian(izgs:izge,0)*integrant(izgs:izge)
+ call simpson_rule_z(integrant,integral)
+ hflux_local = hflux_local + integral*iInt_Jacobian
+ ! 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_xi = 0
+ IF (num_procs_ky .GT. 1) THEN
+ !! Everyone sends its local_sum to root = 0
+ IF (rank_ky .NE. root) THEN
+ CALL MPI_SEND(buffer, 2 , MPI_DOUBLE_PRECISION, root, 1234, comm_ky, ierr)
+ ELSE
+ ! Recieve from all the other processes
+ DO i_ = 0,num_procs_ky-1
+ IF (i_ .NE. rank_ky) &
+ CALL MPI_RECV(buffer, 2 , MPI_DOUBLE_PRECISION, i_, 1234, comm_ky, MPI_STATUS_IGNORE, ierr)
+ hflux_xi = hflux_xi + buffer(2)
ENDDO
- ENDDO
- ! 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_xi = 0
- IF (num_procs_ky .GT. 1) THEN
- !! Everyone sends its local_sum to root = 0
- IF (rank_ky .NE. root) THEN
- CALL MPI_SEND(buffer, 2 , MPI_DOUBLE_PRECISION, root, 1234, comm_ky, ierr)
- ELSE
- ! Recieve from all the other processes
- DO i_ = 0,num_procs_ky-1
- IF (i_ .NE. rank_ky) &
- CALL MPI_RECV(buffer, 2 , MPI_DOUBLE_PRECISION, i_, 1234, comm_ky, MPI_STATUS_IGNORE, ierr)
- hflux_xi = hflux_xi + buffer(2)
- ENDDO
- ENDIF
- ELSE
- hflux_xi = hflux_local
- ENDIF
- ENDIF
+ ENDIF
+ ELSE
+ hflux_xi = hflux_local
+ ENDIF
END SUBROUTINE compute_radial_ion_heatflux
! 1D diagnostic to compute the average radial particle transport <T_e v_ExB_x>_xyz
SUBROUTINE compute_radial_electron_heatflux
- USE fields, ONLY : moments_e, phi
- USE array, ONLY : kernel_e, HF_phi_correction_operator
- USE geometry, ONLY : Jacobian, iInt_Jacobian
- USE time_integration, ONLY : updatetlevel
- USE model, ONLY : qe_taue, KIN_E, tau_e
- USE calculus, ONLY : simpson_rule_z
- IMPLICIT NONE
- COMPLEX(dp) :: hflux_local, integral
- 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 ! heat flux
- IF(ips_e .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 in = ijs_e, ije_e
- n_dp = jarray_e(in)
-
- integrant(izs:ize) = integrant(izs:ize) + Jacobian(izs:ize,0)*tau_e*imagu*ky_*CONJG(phi(iky,ikx,izs:ize))&
- *kernel_e(in,iky,ikx,izs:ize,0)*(&
- 0.5_dp*SQRT2*moments_e(ip2_e,in ,iky,ikx,izs:ize,updatetlevel)&
- +(2._dp*n_dp + 1.5_dp)*moments_e(ip0_e,in ,iky,ikx,izs:ize,updatetlevel)&
- -(n_dp+1._dp)*moments_e(ip0_e,in+1,iky,ikx,izs:ize,updatetlevel)&
- -n_dp*moments_e(ip0_e,in-1,iky,ikx,izs:ize,updatetlevel))
- ENDDO
- ! ! Add polarisation contribution
- ! integrant(izs:ize) = integrant(izs:ize) + Jacobian(izs:ize,0)*tau_e*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
+ USE fields, ONLY : moments_e, phi, psi
+ USE array, ONLY : kernel_e, HF_phi_correction_operator
+ USE geometry, ONLY : Jacobian, iInt_Jacobian
+ USE time_integration, ONLY : updatetlevel
+ USE model, ONLY : qi_taui, KIN_E, tau_e, EM
+ USE calculus, ONLY : simpson_rule_z
+ USE model, ONLY : sqrt_tau_o_sigma_e
+ IMPLICIT NONE
+ COMPLEX(dp) :: hflux_local, integral
+ 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 ! heat flux
+ integrant = 0._dp ! z integration auxiliary variable
+ !!----------------ELECTROSTATIC CONTRIBUTION---------------------------
+ IF(CONTAINS_ip0_e .AND. CONTAINS_ip2_e) 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
+ DO in = ijs_e, ije_e
+ n_dp = jarray_e(in)
+ integrant(izgs:izge) = integrant(izgs:izge) + tau_e*imagu*ky_*CONJG(phi(iky,ikx,izgs:izge))&
+ *kernel_e(in,iky,ikx,izgs:izge,0)*(&
+ 0.5_dp*SQRT2*moments_e(ip2_e,in ,iky,ikx,izgs:izge,updatetlevel)&
+ +(2._dp*n_dp + 1.5_dp)*moments_e(ip0_e,in ,iky,ikx,izgs:izge,updatetlevel)&
+ -(n_dp+1._dp)*moments_e(ip0_e,in+1,iky,ikx,izgs:izge,updatetlevel)&
+ -n_dp*moments_e(ip0_e,in-1,iky,ikx,izgs:izge,updatetlevel))
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+ IF(EM .AND. CONTAINS_ip1_e .AND. CONTAINS_ip3_e) THEN
+ !!----------------ELECTROMAGNETIC CONTRIBUTION--------------------
+ DO iky = ikys,ikye
+ ky_ = kyarray(iky)
+ DO ikx = ikxs,ikxe
+ DO in = ijs_e, ije_e
+ n_dp = jarray_e(in)
+ integrant(izgs:izge) = integrant(izgs:izge) &
+ +tau_e*sqrt_tau_o_sigma_e*imagu*ky_*CONJG(psi(iky,ikx,izgs:izge))&
+ *kernel_e(in,iky,ikx,izgs:izge,0)*(&
+ 0.5_dp*SQRT2*SQRT3*moments_e(ip3_e,in ,iky,ikx,izgs:izge,updatetlevel)&
+ +1.5_dp*CONJG(moments_e(ip1_e,in ,iky,ikx,izgs:izge,updatetlevel))& !?????
+ +(2._dp*n_dp+1._dp)*moments_e(ip1_e,in ,iky,ikx,izgs:izge,updatetlevel)&
+ -(n_dp+1._dp)*moments_e(ip1_e,in+1,iky,ikx,izgs:izge,updatetlevel)&
+ -n_dp*moments_e(ip1_e,in-1,iky,ikx,izgs:izge,updatetlevel))
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+ ! Add polarisation contribution
+ integrant(izs:ize) = integrant(izs:ize) + tau_e*imagu*ky_&
+ *CONJG(phi(iky,ikx,izs:ize))*phi(iky,ikx,izs:ize) * HF_phi_correction_operator(iky,ikx,izs:ize)
+ ! Integrate over z
+ integrant(izgs:izge) = Jacobian(izgs:izge,0)*integrant(izgs:izge)
+ call simpson_rule_z(integrant,integral)
+ hflux_local = hflux_local + integral*iInt_Jacobian
+ ! 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_xi = 0
+ IF (num_procs_ky .GT. 1) THEN
+ !! Everyone sends its local_sum to root = 0
+ IF (rank_ky .NE. root) THEN
+ CALL MPI_SEND(buffer, 2 , MPI_DOUBLE_PRECISION, root, 1234, comm_ky, ierr)
+ ELSE
+ ! Recieve from all the other processes
+ DO i_ = 0,num_procs_ky-1
+ IF (i_ .NE. rank_ky) &
+ CALL MPI_RECV(buffer, 2 , MPI_DOUBLE_PRECISION, i_, 1234, comm_ky, MPI_STATUS_IGNORE, ierr)
+ hflux_xi = hflux_xi + buffer(2)
ENDDO
- ENDDO
- ! 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_xe = 0
- IF (num_procs_ky .GT. 1) THEN
- !! Everyone sends its local_sum to root = 0
- IF (rank_ky .NE. root) THEN
- CALL MPI_SEND(buffer, 2 , MPI_DOUBLE_PRECISION, root, 1234, comm_ky, ierr)
- ELSE
- ! Recieve from all the other processes
- DO i_ = 0,num_procs_ky-1
- IF (i_ .NE. rank_ky) &
- CALL MPI_RECV(buffer, 2 , MPI_DOUBLE_PRECISION, i_, 1234, comm_ky, MPI_STATUS_IGNORE, ierr)
- hflux_xe = hflux_xe + buffer(2)
- ENDDO
- ENDIF
- ELSE
- hflux_xe = hflux_local
- ENDIF
- ENDIF
+ ENDIF
+ ELSE
+ hflux_xi = hflux_local
+ ENDIF
END SUBROUTINE compute_radial_electron_heatflux
+
SUBROUTINE compute_nadiab_moments_z_gradients_and_interp
! evaluate the non-adiabatique ion moments
!
@@ -297,12 +407,12 @@ SUBROUTINE compute_nadiab_moments_z_gradients_and_interp
!
USE fields, ONLY : moments_i, moments_e, phi, psi
USE array, ONLY : kernel_e, kernel_i, nadiab_moments_e, nadiab_moments_i, &
- ddz_nepj, ddz4_Nepj, interp_nepj,&
- ddz_nipj, ddz4_Nipj, interp_nipj
+ ddz_nepj, ddzND_nepj, interp_nepj,&
+ ddz_nipj, ddzND_nipj, interp_nipj
USE time_integration, ONLY : updatetlevel
USE model, ONLY : qe_taue, qi_taui,q_o_sqrt_tau_sigma_e, q_o_sqrt_tau_sigma_i, &
- KIN_E, CLOS
- USE calculus, ONLY : grad_z, grad_z4, interp_z
+ KIN_E, CLOS, beta
+ USE calculus, ONLY : grad_z, grad_z2, grad_z4, interp_z
IMPLICIT NONE
INTEGER :: eo, p_int, j_int
CALL cpu_time(t0_process)
@@ -313,17 +423,17 @@ SUBROUTINE compute_nadiab_moments_z_gradients_and_interp
DO ip=ipgs_e,ipge_e
IF(parray_e(ip) .EQ. 0) THEN
DO ij=ijgs_e,ijge_e
- nadiab_moments_e(ip,ij,:,:,:) = moments_e(ip,ij,:,:,:,updatetlevel) &
- + qe_taue*kernel_e(ij,:,:,:,0)*phi(:,:,:)
+ nadiab_moments_e(ip,ij,ikys:ikye,ikxs:ikxe,izgs:izge) = moments_e(ip,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) &
+ + qe_taue*kernel_e(ij,ikys:ikye,ikxs:ikxe,izgs:izge,0)*phi(ikys:ikye,ikxs:ikxe,izgs:izge)
ENDDO
- ELSEIF(parray_e(ip) .EQ. 1) THEN
+ ELSEIF( (parray_e(ip) .EQ. 1) .AND. (beta .GT. 0) ) THEN
DO ij=ijgs_e,ijge_e
- nadiab_moments_e(ip,ij,:,:,:) = moments_e(ip,ij,:,:,:,updatetlevel) &
- - q_o_sqrt_tau_sigma_e*kernel_e(ij,:,:,:,0)*psi(:,:,:)
+ nadiab_moments_e(ip,ij,ikys:ikye,ikxs:ikxe,izgs:izge) = moments_e(ip,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) &
+ - q_o_sqrt_tau_sigma_e*kernel_e(ij,ikys:ikye,ikxs:ikxe,izgs:izge,0)*psi(ikys:ikye,ikxs:ikxe,izgs:izge)
ENDDO
ELSE
DO ij=ijgs_e,ijge_e
- nadiab_moments_e(ip,ij,:,:,:) = moments_e(ip,ij,:,:,:,updatetlevel)
+ nadiab_moments_e(ip,ij,ikys:ikye,ikxs:ikxe,izgs:izge) = moments_e(ip,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel)
ENDDO
ENDIF
ENDDO
@@ -332,17 +442,17 @@ SUBROUTINE compute_nadiab_moments_z_gradients_and_interp
DO ip=ipgs_i,ipge_i
IF(parray_i(ip) .EQ. 0) THEN
DO ij=ijgs_i,ijge_i
- nadiab_moments_i(ip,ij,:,:,:) = moments_i(ip,ij,:,:,:,updatetlevel) &
- + qi_taui*kernel_i(ij,:,:,:,0)*phi(:,:,:)
+ nadiab_moments_i(ip,ij,ikys:ikye,ikxs:ikxe,izgs:izge) = moments_i(ip,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) &
+ + qi_taui*kernel_i(ij,ikys:ikye,ikxs:ikxe,izgs:izge,0)*phi(ikys:ikye,ikxs:ikxe,izgs:izge)
ENDDO
- ELSEIF(parray_i(ip) .EQ. 1) THEN
+ ELSEIF( (parray_i(ip) .EQ. 1) .AND. (beta .GT. 0) ) THEN
DO ij=ijgs_i,ijge_i
- nadiab_moments_i(ip,ij,:,:,:) = moments_i(ip,ij,:,:,:,updatetlevel) &
- - q_o_sqrt_tau_sigma_i*kernel_i(ij,:,:,:,0)*psi(:,:,:)
+ nadiab_moments_i(ip,ij,ikys:ikye,ikxs:ikxe,izgs:izge) = moments_i(ip,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) &
+ - q_o_sqrt_tau_sigma_i*kernel_i(ij,ikys:ikye,ikxs:ikxe,izgs:izge,0)*psi(ikys:ikye,ikxs:ikxe,izgs:izge)
ENDDO
ELSE
DO ij=ijgs_i,ijge_i
- nadiab_moments_i(ip,ij,:,:,:) = moments_i(ip,ij,:,:,:,updatetlevel)
+ nadiab_moments_i(ip,ij,ikys:ikye,ikxs:ikxe,izgs:izge) = moments_i(ip,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel)
ENDDO
ENDIF
ENDDO
@@ -381,7 +491,8 @@ SUBROUTINE compute_nadiab_moments_z_gradients_and_interp
! Compute z derivatives
CALL grad_z(eo,nadiab_moments_e(ip,ij,iky,ikx,izgs:izge), ddz_nepj(ip,ij,iky,ikx,izs:ize))
! Parallel hyperdiffusion
- CALL grad_z4(moments_e(ip,ij,iky,ikx,izgs:izge,updatetlevel),ddz4_Nepj(ip,ij,iky,ikx,izs:ize))
+ CALL grad_z4(moments_e(ip,ij,iky,ikx,izgs:izge,updatetlevel),ddzND_nepj(ip,ij,iky,ikx,izs:ize))
+ ! CALL grad_z2(nadiab_moments_e(ip,ij,iky,ikx,izgs:izge),ddzND_nepj(ip,ij,iky,ikx,izs:ize))
! Compute even odd grids interpolation
CALL interp_z(eo,nadiab_moments_e(ip,ij,iky,ikx,izgs:izge), interp_nepj(ip,ij,iky,ikx,izs:ize))
ENDDO
@@ -398,8 +509,9 @@ SUBROUTINE compute_nadiab_moments_z_gradients_and_interp
eo = MODULO(p_int,2) ! Indicates if we are on even or odd z grid
! Compute z first derivative
CALL grad_z(eo,nadiab_moments_i(ip,ij,iky,ikx,izgs:izge), ddz_nipj(ip,ij,iky,ikx,izs:ize))
- ! Parallel hyperdiffusion
- CALL grad_z4(moments_i(ip,ij,iky,ikx,izgs:izge,updatetlevel),ddz4_Nipj(ip,ij,iky,ikx,izs:ize))
+ ! Parallel numerical diffusion
+ CALL grad_z4(moments_i(ip,ij,iky,ikx,izgs:izge,updatetlevel),ddzND_nipj(ip,ij,iky,ikx,izs:ize))
+ ! CALL grad_z2(nadiab_moments_i(ip,ij,iky,ikx,izgs:izge),ddzND_nipj(ip,ij,iky,ikx,izs:ize))
! Compute even odd grids interpolation
CALL interp_z(eo,nadiab_moments_i(ip,ij,iky,ikx,izgs:izge), interp_nipj(ip,ij,iky,ikx,izs:ize))
ENDDO
--
GitLab