From 130acad6b1b50be9f68b858b5ac16b6ac380fdd5 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Tue, 12 Sep 2023 16:33:40 +0200
Subject: [PATCH] need of an xgrid array and local_nx
---
src/grid_mod.F90 | 72 ++++++++++++++++++++++++++++++------------------
1 file changed, 45 insertions(+), 27 deletions(-)
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 70fd7e09..f5a30438 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -57,7 +57,8 @@ MODULE grid
INTEGER, PUBLIC, PROTECTED :: local_np
INTEGER, PUBLIC, PROTECTED :: local_nj
INTEGER, PUBLIC, PROTECTED :: local_nky
- INTEGER, PUBLIC, PROTECTED :: local_nkx
+ INTEGER, PUBLIC, PROTECTED :: local_nkx ! = total_Nkx = Nx, not parallel
+ INTEGER, PUBLIC, PROTECTED :: local_nx ! = local_nx_ptr from FFTW (used only for Fourier)
INTEGER, PUBLIC, PROTECTED :: local_nz
INTEGER, PUBLIC, PROTECTED :: local_nkp
INTEGER, PUBLIC, PROTECTED :: ngp, ngj, ngx, ngy, ngz ! number of ghosts points
@@ -68,10 +69,11 @@ MODULE grid
INTEGER, PUBLIC, PROTECTED :: local_nj_offset
INTEGER, PUBLIC, PROTECTED :: local_nky_offset
INTEGER, PUBLIC, PROTECTED :: local_nkx_offset
+ INTEGER, PUBLIC, PROTECTED :: local_nx_offset
INTEGER, PUBLIC, PROTECTED :: local_nz_offset
! C-pointer type for FFTW3
- integer(C_INTPTR_T), PUBLIC,PROTECTED :: local_nkx_ptr, local_nky_ptr
- integer(C_INTPTR_T), PUBLIC,PROTECTED :: local_nkx_ptr_offset, local_nky_ptr_offset
+ integer(C_INTPTR_T), PUBLIC,PROTECTED :: local_nx_ptr, local_nky_ptr
+ integer(C_INTPTR_T), PUBLIC,PROTECTED :: local_nx_ptr_offset, local_nky_ptr_offset
! Grid spacing and limits
REAL(xp), PUBLIC, PROTECTED :: deltap, deltaz, inv_deltaz, inv_dkx
REAL(xp), PUBLIC, PROTECTED :: deltakx, deltaky, deltax, kx_max, ky_max, kx_min, ky_min!, kp_max
@@ -198,7 +200,7 @@ CONTAINS
!! Parallel distribution of kx ky grid
IF (LINEARITY .NE. 'linear') THEN ! we let FFTW distribute if we use it
IF (my_id .EQ. 0) write(*,*) 'FFTW3 y-grid distribution'
- CALL init_grid_distr_and_plans(FFT2D,Nx,Ny,comm_ky,local_nkx_ptr,local_nkx_ptr_offset,local_nky_ptr,local_nky_ptr_offset)
+ CALL init_grid_distr_and_plans(FFT2D,Nx,Ny,comm_ky,local_nx_ptr,local_nx_ptr_offset,local_nky_ptr,local_nky_ptr_offset)
ELSE ! otherwise we distribute equally
IF (my_id .EQ. 0) write(*,*) 'Manual y-grid distribution'
! balanced distribution among the processes
@@ -207,30 +209,32 @@ CONTAINS
local_nky_ptr_offset = ikys - 1
ENDIF
!!----------------- BINORMAL KY INDICES (parallelized)
- Nky = Ny/2+1 ! Defined only on positive kx since fields are real
- total_nky = Nky
- Ngy = 0 ! no ghosts cells in ky
- ikys = local_nky_ptr_offset + 1
- ikye = ikys + local_nky_ptr - 1
- local_nky = ikye - ikys + 1
+ Nky = Ny/2+1 ! Defined only on positive kx since fields are real
+ total_nky = Nky
+ Ngy = 0 ! no ghosts cells in ky
+ ikys = local_nky_ptr_offset + 1
+ ikye = ikys + local_nky_ptr - 1
+ local_nky = ikye - ikys + 1
local_nky_offset = local_nky_ptr_offset
ALLOCATE(kyarray_full(Nky))
ALLOCATE(kyarray(local_nky))
- ALLOCATE(ikyarray(local_nky))
- ALLOCATE(inv_ikyarray(local_nky))
+ ALLOCATE(ikyarray(Nky))
+ ALLOCATE(inv_ikyarray(Nky))
ALLOCATE(AA_y(local_nky))
!!---------------- RADIAL KX INDICES (not parallelized)
- Nkx = Nx
- total_nkx = Nx
- ikxs = 1
- ikxe = total_nkx
- local_nkx_ptr = ikxe - ikxs + 1
- local_nkx = ikxe - ikxs + 1
+ Nkx = Nx
+ total_nkx = Nx
+ ikxs = 1
+ ikxe = total_nkx
+ local_nkx = ikxe - ikxs + 1
local_nkx_offset = ikxs - 1
ALLOCATE(kxarray_full(total_nkx))
ALLOCATE(kxarray(local_nky,local_Nkx))
- ALLOCATE(xarray(total_nkx))
ALLOCATE(AA_x(local_nkx))
+ !!---------------- RADIAL X GRID (only for Fourier routines)
+ local_nx = local_nx_ptr
+ local_nx_offset = local_nx_ptr_offset
+ ALLOCATE(xarray(Nx))
!!---------------- PARALLEL Z GRID (parallelized)
total_nz = Nz
IF (SG) THEN
@@ -299,6 +303,7 @@ CONTAINS
CALL set_jgrid
CALL set_kygrid(LINEARITY,N_HD)
CALL set_kxgrid(shear,Npol,1._xp,LINEARITY,N_HD) ! this will be redone after geometry if Cyq0_x0 .NE. 1
+ CALL set_xgrid
CALL set_zgrid (Npol)
END SUBROUTINE set_grids
@@ -428,6 +433,13 @@ CONTAINS
! Build the full grids on process 0 to diagnose it without comm
DO iky = 1,Nky
kyarray_full(iky) = REAL(iky-1,xp) * deltaky
+ ! full indices grid (for ExB NL factor)
+ ikyarray(iky) = REAL((iky)-1,xp)
+ IF(ikyarray(iky) .GT. 0) THEN
+ inv_ikyarray(iky) = 1._xp/ikyarray(iky)
+ ELSE
+ inv_ikyarray(iky) = 0._xp
+ ENDIF
END DO
local_kymax = 0._xp
! Creating a grid ordered as dk*(0 1 2 3)
@@ -442,12 +454,6 @@ CONTAINS
SINGLE_KY = .TRUE.
ELSE
kyarray(iky) = kyarray_full(iky+local_nky_offset)
- ikyarray(iky) = REAL((iky+local_nky_offset)-1,xp)
- IF(ikyarray(iky) .GT. 0) THEN
- inv_ikyarray(iky) = 1._xp/ikyarray(iky)
- ELSE
- inv_ikyarray(iky) = 0._xp
- ENDIF
ENDIF
! Finding kx=0
IF (kyarray(iky) .EQ. 0) THEN
@@ -503,13 +509,11 @@ CONTAINS
ENDIF
deltakx = 2._xp*PI/Lx
inv_dkx = 1._xp/deltakx
- deltax = Lx/REAL(Nx,xp) ! periodic donc pas Lx/(Nx-1)
IF(MODULO(total_nkx,2) .EQ. 0) THEN ! Even number of kx (-2 -1 0 1 2 3)
! Creating a grid ordered as dk*(0 1 2 3 -2 -1)
DO ikx = 1,total_nkx
kxarray_full(ikx) = deltakx*REAL(MODULO(ikx-1,total_nkx/2)-(total_nkx/2)*FLOOR(2.*real(ikx-1)/real(total_nkx)),xp)
IF (ikx .EQ. total_nkx/2+1) kxarray_full(ikx) = -kxarray_full(ikx)
- xarray(ikx) = REAL(ikx-1,xp)*deltax
END DO
kx_max = MAXVAL(kxarray_full)!(total_nkx/2)*deltakx
kx_min = MINVAL(kxarray_full)!-kx_max+deltakx
@@ -554,6 +558,20 @@ CONTAINS
ENDIF
END SUBROUTINE set_kxgrid
+ !----------- Radial x grid
+ ! used only for the computation of the ExB shear NL factor
+ SUBROUTINE set_xgrid
+ USE prec_const
+ IMPLICIT NONE
+ INTEGER :: ix
+ deltax = Lx/REAL(Nx,xp) ! periodic donc pas Lx/(Nx-1)
+ ! full xgrid
+ DO ix = 1,Nx
+ xarray(ix) = REAL(ix-1,xp)*deltax
+ ENDDO
+ END SUBROUTINE set_xgrid
+
+ !----------- Parallel z grid
SUBROUTINE set_zgrid(Npol)
USE prec_const
IMPLICIT NONE
--
GitLab