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

EM for transport

parent 4b618ae3
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 380 additions and 268 deletions
src/processing_mod.F90
+ 380
− 268
View file @ 96062c8b
......@@ -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
......
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