From 5d9e3a380a14c07a3f53346f26d28410b5ad5da2 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Sun, 19 Mar 2023 18:49:32 +0100
Subject: [PATCH] Mpi parallelization looks fine, need to check adiab e
---
src/calculus_mod.F90 | 10 +++-------
src/geometry_mod.F90 | 8 ++++----
src/grid_mod.F90 | 19 ++++++++++++++-----
src/moments_eq_rhs_mod.F90 | 30 +++++++++++++++---------------
src/parallel_mod.F90 | 1 -
src/processing_mod.F90 | 10 +++++-----
src/solve_EM_fields.F90 | 4 ++--
7 files changed, 43 insertions(+), 39 deletions(-)
diff --git a/src/calculus_mod.F90 b/src/calculus_mod.F90
index f36ce8c2..613ec92e 100644
--- a/src/calculus_mod.F90
+++ b/src/calculus_mod.F90
@@ -170,7 +170,7 @@ CONTAINS
END SUBROUTINE interp_z
-SUBROUTINE simpson_rule_z(local_nz,dz,f,intf)
+SUBROUTINE simpson_rule_z(local_nz,zweights_SR,f,intf)
! integrate f(z) over z using the simpon's rule. Assume periodic boundary conditions (f(ize+1) = f(izs))
!from molix BJ Frei
USE prec_const, ONLY: dp, onethird
@@ -178,7 +178,7 @@ SUBROUTINE simpson_rule_z(local_nz,dz,f,intf)
USE mpi
implicit none
INTEGER, INTENT(IN) :: local_nz
- REAL(dp),INTENT(IN) :: dz
+ REAL(dp), dimension(local_nz), intent(in) :: zweights_SR
complex(dp),dimension(local_nz), intent(in) :: f
COMPLEX(dp), intent(out) :: intf
COMPLEX(dp) :: buffer, local_int
@@ -186,11 +186,7 @@ SUBROUTINE simpson_rule_z(local_nz,dz,f,intf)
! Buil local sum using the weights of composite Simpson's rule
local_int = 0._dp
DO iz = 1,local_nz
- IF(MODULO(iz,2) .EQ. 1) THEN ! odd iz
- local_int = local_int + 2._dp*onethird*dz*f(iz)
- ELSE ! even iz
- local_int = local_int + 4._dp*onethird*dz*f(iz)
- ENDIF
+ local_int = local_int + zweights_SR(iz)*f(iz)
ENDDO
buffer = local_int
root = 0
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index 6afb1fed..c354536f 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -93,7 +93,8 @@ CONTAINS
END SUBROUTINE geometry_readinputs
subroutine eval_magnetic_geometry
- USE grid, ONLY: total_nky, total_nz, local_nkx, local_nky, local_nz, Ngz, kxarray, kyarray, set_kparray, nzgrid, deltaz
+ USE grid, ONLY: total_nky, total_nz, local_nkx, local_nky, local_nz, Ngz,&
+ kxarray, kyarray, set_kparray, nzgrid, zweights_SR, ieven
USE basic, ONLY: speak
USE miller, ONLY: set_miller_parameters, get_miller
USE calculus, ONLY: simpson_rule_z
@@ -165,8 +166,8 @@ CONTAINS
CALL set_ikx_zBC_map
!
! Compute the inverse z integrated Jacobian (useful for flux averaging)
- integrant = Jacobian(1+ngz/2:local_nz+ngz/2,1) ! Convert into complex array
- CALL simpson_rule_z(local_nz,deltaz,integrant,iInt_Jacobian)
+ integrant = Jacobian((1+ngz/2):(local_nz+ngz/2),ieven) ! Convert into complex array
+ CALL simpson_rule_z(local_nz,zweights_SR,integrant,iInt_Jacobian)
iInt_Jacobian = 1._dp/iInt_Jacobian ! reverse it
END SUBROUTINE eval_magnetic_geometry
!
@@ -220,7 +221,6 @@ CONTAINS
ENDDO
ENDDO
-
END SUBROUTINE eval_salpha_geometry
!
!--------------------------------------------------------------------------------
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 513fd0a2..5a1604b6 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -504,10 +504,10 @@ CONTAINS
local_nz_offset = izs - 1
! Ghosts boundaries (depend on the order of z operators)
IF(Nz .EQ. 1) THEN
- Ngz = 0
+ ngz = 0
zarray_full(izs) = 0
ELSEIF(Nz .GE. 4) THEN
- Ngz =4
+ ngz =4
ELSE
ERROR STOP '>> ERROR << Nz is not appropriate!!'
ENDIF
@@ -520,11 +520,11 @@ CONTAINS
displs_nz(in+1) = istart-1
ENDDO
! Local z array
- ALLOCATE(zarray(local_nz+Ngz,nzgrid))
+ ALLOCATE(zarray(local_nz+ngz,nzgrid))
!! interior point loop
- DO iz = 1,total_nz
+ DO iz = 1,local_nz
DO eo = 1,nzgrid
- zarray(iz+ngz/2,eo) = zarray_full(iz) + REAL(eo-1,dp)*grid_shift
+ zarray(iz+ngz/2,eo) = zarray_full(iz+local_nz_offset) + REAL(eo-1,dp)*grid_shift
ENDDO
ENDDO
CALL allocate_array(local_zmax,1,nzgrid)
@@ -562,6 +562,15 @@ CONTAINS
IF(mod(Nz,2) .NE. 0 ) THEN
ERROR STOP '>> ERROR << Nz must be an even number for Simpson integration rule !!!!'
ENDIF
+ ! local weights for Simpson rule
+ ALLOCATE(zweights_SR(local_nz))
+ DO iz = 1,local_nz
+ IF(MODULO(iz+local_nz_offset,2) .EQ. 1) THEN ! odd iz
+ zweights_SR(iz) = onethird*deltaz*2._dp
+ ELSE ! even iz
+ zweights_SR(iz) = onethird*deltaz*4._dp
+ ENDIF
+ ENDDO
END SUBROUTINE set_zgrid
SUBROUTINE set_kparray(gxx, gxy, gyy,hatB)
diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90
index 29c658c5..eecef868 100644
--- a/src/moments_eq_rhs_mod.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -159,21 +159,21 @@ SUBROUTINE compute_moments_eq_rhs
CALL cpu_time(t1_rhs)
tc_rhs = tc_rhs + (t1_rhs-t0_rhs)
- ! print*,'sumabs moments i', SUM(ABS(moments(1,:,:,:,:,:,updatetlevel)))
- ! print*,'moment rhs i 221', moments_rhs(1,1,1,2,2,1,updatetlevel)
- ! print*,'sum real nadiabe', SUM(REAL(nadiab_moments(2,:,:,:,:,:)))
- ! print*,'sumreal momrhs i', SUM(REAL(moments_rhs(1,:,:,:,:,:,:)))
- ! print*,'sumimag momrhs i', SUM(IMAG(moments_rhs(1,:,:,:,:,:,:)))
- ! print*,'sumreal phi ', SUM(REAL(phi(:,:,(1+ngz/2):(local_nz+ngz/2))))
- ! print*,'sumimag phi ', SUM(IMAG(phi(:,:,(1+ngz/2):(local_nz+ngz/2))))
- ! print*,'phi(2,2,1) ', phi(2,2,1+ngz/2)
- ! print*,'sumreal ddznipj ', SUM(REAL(ddz_napj(1,:,:,:,:,:)))
- ! print*,'sumimag ddznipj ', SUM(IMAG(ddz_napj(1,:,:,:,:,:)))
- ! print*,' ddznipj ',(ddz_napj(1,2+ngp/2,2+ngj/2,2,2,1))
- ! print*,'sumreal Capj ', SUM(REAL(Capj(1,:,:,:,:,:)))
- ! print*,'---'
- ! IF(updatetlevel .EQ. 4) stop
- ! ! stop
+ ! print*,'sumabs moments i', SUM(ABS(moments(1,:,:,:,:,:,updatetlevel)))
+ ! print*,'moment rhs i 221', moments_rhs(1,1,1,2,2,1,updatetlevel)
+ ! print*,'sum real nadiabe', SUM(REAL(nadiab_moments(2,:,:,:,:,:)))
+ ! print*,'sumreal momrhs i', SUM(REAL(moments_rhs(1,:,:,:,:,:,:)))
+ ! print*,'sumimag momrhs i', SUM(IMAG(moments_rhs(1,:,:,:,:,:,:)))
+ ! print*,'sumreal phi ', SUM(REAL(phi(:,:,(1+ngz/2):(local_nz+ngz/2))))
+ ! print*,'sumimag phi ', SUM(IMAG(phi(:,:,(1+ngz/2):(local_nz+ngz/2))))
+ ! print*,'phi(2,2,1) ', phi(2,2,1+ngz/2)
+ ! print*,'sumreal ddznipj ', SUM(REAL(ddz_napj(1,:,:,:,:,:)))
+ ! print*,'sumimag ddznipj ', SUM(IMAG(ddz_napj(1,:,:,:,:,:)))
+ ! print*,' ddznipj ',(ddz_napj(1,2+ngp/2,2+ngj/2,2,2,1))
+ ! print*,'sumreal Capj ', SUM(REAL(Capj(1,:,:,:,:,:)))
+ ! print*,'---'
+ ! IF(updatetlevel .EQ. 4) stop
+ ! stop
END SUBROUTINE compute_moments_eq_rhs
diff --git a/src/parallel_mod.F90 b/src/parallel_mod.F90
index a9f3df34..e4469d0f 100644
--- a/src/parallel_mod.F90
+++ b/src/parallel_mod.F90
@@ -47,7 +47,6 @@ CONTAINS
INTEGER, PARAMETER :: ndims=3 ! p, kx and z
INTEGER, DIMENSION(ndims) :: dims=0, coords=0
LOGICAL :: periods(ndims) = .FALSE., reorder=.FALSE.
- LOGICAL :: com_log(ndims) = .FALSE.
CHARACTER(len=32) :: str
INTEGER :: nargs, i, l
CALL MPI_INIT(ierr)
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index 65bfea95..47ede3ef 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -7,7 +7,7 @@ MODULE processing
jarray,jmax,ij0, dmax,&
kyarray, AA_y,&
kxarray, AA_x,&
- zarray, deltaz, ieven, iodd, inv_deltaz
+ zarray, zweights_SR, ieven, iodd, inv_deltaz
USE fields, ONLY: moments, phi, psi
USE array, ONLY : kernel, nadiab_moments, &
ddz_napj, ddzND_Napj, interp_napj,&
@@ -81,7 +81,7 @@ CONTAINS
ENDDO
ENDIF
! Integrate over z
- call simpson_rule_z(local_nz,deltaz,integrant,integral)
+ call simpson_rule_z(local_nz,zweights_SR,integrant,integral)
! Get process local gyrocenter flux with a factor two to account for the negative ky modes
gflux_local = 2._dp*integral*iInt_Jacobian
!
@@ -120,7 +120,7 @@ CONTAINS
ENDDO
ENDIF
! Integrate over z
- call simpson_rule_z(local_nz,deltaz,integrant,integral)
+ call simpson_rule_z(local_nz,zweights_SR,integrant,integral)
! Get process local particle flux with a factor two to account for the negative ky modes
pflux_local = 2._dp*integral*iInt_Jacobian
!!!!---------- Sum over all processes ----------
@@ -194,7 +194,7 @@ CONTAINS
ini = in + ngj/2 !interior index for ghosted arrays
n_dp = jarray(ini)
integrant(iz) = integrant(iz) &
- +Jacobian(izi,iodd)*tau(ia)*sqrt_tau_o_sigma(ia)*imagu*kyarray(iky)*CONJG(psi(iky,ikx,izi))&
+ +Jacobian(izi,iodd)*tau(ia)*sqrt_tau_o_sigma(ia)*imagu*kyarray(iky)*CONJG(psi(iky,ikx,izi))&
*kernel(ia,ini,iky,ikx,izi,iodd)*(&
0.5_dp*SQRT2*SQRT3*moments(ia,ip3,ini ,iky,ikx,izi,updatetlevel)&
+1.5_dp*moments(ia,ip1,ini ,iky,ikx,izi,updatetlevel)&
@@ -210,7 +210,7 @@ CONTAINS
! integrant(iz) = integrant(iz) + tau_i*imagu*ky_&
! *CONJG(phi(iky,ikx,iz))*phi(iky,ikx,iz) * HF_phi_correction_operator(iky,ikx,iz)
! Integrate over z
- call simpson_rule_z(local_nz,deltaz,integrant,integral)
+ call simpson_rule_z(local_nz,zweights_SR,integrant,integral)
! Double it for taking into account the other half plane
hflux_local = 2._dp*integral*iInt_Jacobian
buffer(2) = REAL(hflux_local,dp)
diff --git a/src/solve_EM_fields.F90 b/src/solve_EM_fields.F90
index db3b3dd4..a365937d 100644
--- a/src/solve_EM_fields.F90
+++ b/src/solve_EM_fields.F90
@@ -18,7 +18,7 @@ CONTAINS
USE array, ONLY: kernel, inv_poisson_op, inv_pol_ion
USE fields, ONLY: phi, moments
USE grid, ONLY: local_na, local_nky, local_nkx, local_nz,ngz, SOLVE_POISSON,&
- kyarray, contains_kx0, contains_ky0,ikx0,iky0, deltaz, ieven,&
+ kyarray, contains_kx0, contains_ky0,ikx0,iky0, zweights_SR, ieven,&
ip0, total_nj, ngj
USE calculus, ONLY: simpson_rule_z
USE parallel, ONLY: manual_3D_bcast
@@ -61,7 +61,7 @@ CONTAINS
IF(kyarray(iky).EQ.0._dp) THEN ! take ky=0 mode (y-average)
! Prepare integrant for z-average
integrant(iz) = Jacobian(iz+ngz/2,ieven)*rho(iz)*inv_pol_ion(iky,ikx,iz)
- call simpson_rule_z(local_nz,deltaz,integrant,intf_) ! get the flux averaged phi
+ call simpson_rule_z(local_nz,zweights_SR,integrant,intf_) ! get the flux averaged phi
fsa_phi = intf_ * iInt_Jacobian !Normalize by 1/int(Jxyz)dz
ENDIF
rho(iz) = rho(iz) + fsa_phi
--
GitLab