From 503fa0fa9f9600771c34c351f59f8e3157869f4e Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Wed, 3 Mar 2021 11:23:23 +0100
Subject: [PATCH] start of parallel along p

---
 src/basic_mod.F90 |  9 ++++--
 src/ghosts.F90    | 79 +++++++++++++++++++++++++++++++++++++++++++++++
 src/ppinit.F90    | 74 +++++++++++++++++++++++++++++++++++++++++++-
 3 files changed, 159 insertions(+), 3 deletions(-)
 create mode 100644 src/ghosts.F90

diff --git a/src/basic_mod.F90 b/src/basic_mod.F90
index 6625c8e3..5a7e4a80 100644
--- a/src/basic_mod.F90
+++ b/src/basic_mod.F90
@@ -10,7 +10,10 @@ MODULE basic
   real(dp) :: tmax   = 100000.0    ! Maximum simulation time
   real(dp) :: dt     = 1.0         ! Time step
   real(dp) :: time   = 0           ! Current simulation time (Init from restart file)
-
+  
+  INTEGER :: comm0                 ! Default communicator with a topology
+  INTEGER :: commp, commr          ! Communicators for 1-dim cartesian subgrids of comm0
+  
   INTEGER :: jobnum  = 0           ! Job number
   INTEGER :: step    = 0           ! Calculation step of this run
   INTEGER :: cstep   = 0           ! Current step number (Init from restart file)
@@ -19,8 +22,10 @@ MODULE basic
   LOGICAL :: crashed = .FALSE.     ! Signal end of crashed run
 
   INTEGER :: ierr                  ! flag for MPI error
-  INTEGER :: my_id                 ! identification number of current process
+  INTEGER :: my_id                 ! Rank in COMM_WORLD
   INTEGER :: num_procs             ! number of MPI processes
+  INTEGER :: ncp, ncr              ! Number of processes in p and r
+  INTEGER :: me_0, me_p, me_z      ! Ranks in comm0, commp, commz
 
   INTEGER :: iframe1d              ! counting the number of times 1d datasets are outputed (for diagnose)
   INTEGER :: iframe2d              ! counting the number of times 2d datasets are outputed (for diagnose)
diff --git a/src/ghosts.F90 b/src/ghosts.F90
new file mode 100644
index 00000000..f099c055
--- /dev/null
+++ b/src/ghosts.F90
@@ -0,0 +1,79 @@
+
+!Communicate two extra buffer zones needed for using left/right fd schemes in parallel 
+!gradients routines
+SUBROUTINE update_ghosts_p
+    USE basic
+    USE fields
+    USE grid
+    USE ppinit
+  
+    use prec_const
+    IMPLICIT NONE
+     
+    INTEGER :: status(MPI_STATUS_SIZE)
+    complex(dp):: buff_snd_L( 1: 2,ijs_e:ije_e,ikrs:ikre,ikzs:ikze)
+    complex(dp):: buff_snd_R(-2:-1,ijs_e:ije_e,ikrs:ikre,ikzs:ikze)
+    complex(dp):: buff_rcv_R( 1: 2,ijs_e:ije_e,ikrs:ikre,ikzs:ikze)
+    complex(dp):: buff_rcv_L(-2:-1,ijs_e:ije_e,ikrs:ikre,ikzs:ikze)
+  
+    !! Set up data to send to left neighbor
+    DO ij=ijs_e,ije_e
+        DO ikr=ikrs,ikre
+            DO ikz=ikzs,ikze
+                ! send to left ipe + 1 moments
+                buff_snd_L(1,ij,ikr,ikz) = moments_e(ips,ij,ikr,ikz,updatetlevel)
+                ! send to left ipe + 2 moments
+                buff_snd_L(2,ij,ikr,ikz) = moments_e(ips+1,ij,ikr,ikz,updatetlevel)
+            END DO
+        END DO
+    END DO
+
+    ! Exchange data with left neighbor
+    CALL mpi_sendrecv(buff_snd_L, count, MPI_DOUBLE_PRECISION, left_neighbor, 0, &
+    buff_rcv_L, count, MPI_DOUBLE_PRECISION, source, 0, &
+    commx, status, ierr)
+   
+    ! Write received data
+    DO ij=ijs_e,ije_e
+        DO ikr=ikrs,ikre
+            DO ikz=ikzs,ikze  
+                ! ip - 1 moments
+                moments_e(ips-1,ij,ikr,ikz,updatetlevel) = buff_rcv_L(-1,ij,ikr,ikz)
+                ! ip - 2 moments
+                moments_e(ips-2,ij,ikr,ikz,updatetlevel) = buff_rcv_L(-2,ij,ikr,ikz)
+            END DO
+        END DO
+    END DO 
+
+    !! Set up data to send to right neighbors
+    DO ij=ijs_e,ije_e
+        DO ikr=ikrs,ikre
+            DO ikz=ikzs,ikze
+                ! send to right ips - 1 moments
+                buff_snd_R(-1,ij,ikr,ikz) = moments_e(ipe,ij,ikr,ikz,updatetlevel)
+                ! send to right ips - 2 moments
+                buff_snd_R(-2,ij,ikr,ikz) = moments_e(ipe-1,ij,ikr,ikz,updatetlevel)
+            END DO
+        END DO
+    END DO
+   
+    ! Exchange data with right neighbor
+    CALL mpi_sendrecv(buff_snd_R, count, MPI_DOUBLE_PRECISION, dest, 0, &
+    buff_rcv_R, count, MPI_DOUBLE_PRECISION, source, 0, &
+    commx, status, ierr)
+
+    ! Write received data
+    DO ij=ijs_e,ije_e
+        DO ikr=ikrs,ikre
+            DO ikz=ikzs,ikze  
+                ! ipe + 1 moments
+                moments_e(ipe+1,ij,ikr,ikz,updatetlevel) = buff_rcv_R(1,ij,ikr,ikz)
+                ! ip + 2 moments
+                moments_e(ipe+2,ij,ikr,ikz,updatetlevel) = buff_rcv_R(2,ij,ikr,ikz)
+            END DO
+        END DO
+    END DO 
+
+    
+  END SUBROUTINE update_ghosts_p
+  
\ No newline at end of file
diff --git a/src/ppinit.F90 b/src/ppinit.F90
index 57d7dc39..9e8b8632 100644
--- a/src/ppinit.F90
+++ b/src/ppinit.F90
@@ -2,11 +2,18 @@ SUBROUTINE ppinit
   !   Parallel environment
 
   USE basic
-  USE model, ONLY : NON_LIN
   use prec_const
   IMPLICIT NONE
 
   INTEGER :: version_prov=-1
+  ! Variables for cartesian domain decomposition
+  INTEGER, PARAMETER :: ndims=2 ! p and kr
+  INTEGER, DIMENSION(ndims) :: dims=0, coords=0, coords_L=0, coords_R=0
+  LOGICAL :: periods(ndims) = .FALSE., reorder=.FALSE.
+  CHARACTER(len=32) :: str
+  INTEGER :: nargs, i, l
+  INTEGER :: nghb_L, nghb_R ! Left and right neighbors along p
+  INTEGER :: source
 
   CALL MPI_INIT(ierr)
   ! CALL MPI_INIT_THREAD(MPI_THREAD_SINGLE,version_prov,ierr)
@@ -15,4 +22,69 @@ SUBROUTINE ppinit
   CALL MPI_COMM_RANK (MPI_COMM_WORLD,     my_id, ierr)
   CALL MPI_COMM_SIZE (MPI_COMM_WORLD, num_procs, ierr)
   
+  ! nargs = COMMAND_ARGUMENT_COUNT()
+  ! IF( nargs .NE. 0 .AND. nargs .NE. ndims ) THEN
+  !    IF(my_id .EQ. 0) WRITE(*, '(a,i4,a)') 'Number of arguments not equal to NDIMS =', ndims, '!'
+  !    CALL MPI_ABORT(MPI_COMM_WORLD, -1, ierr)
+  ! END IF
+  ! !
+  ! IF( nargs .NE. 0 ) THEN
+  !    DO i=1,nargs
+  !       CALL GET_COMMAND_ARGUMENT(i, str, l, ierr)
+  !       READ(str(1:l),'(i3)')  dims(i)
+  !       ncp = dims(1) ! Number of processes along p
+  !       ncr = dims(2)  ! Number of processes along kr
+  !    END DO
+  !    IF( PRODUCT(dims) .NE. num_procs ) THEN
+  !     IF(my_id .EQ. 0) WRITE(*, '(a,i4,a,i4)') 'Product of dims: ', PRODUCT(dims), " is not consistent WITH NPROCS=",num_procs
+  !       CALL MPI_ABORT(MPI_COMM_WORLD, -2, ierr)
+  !    END IF
+  ! ELSE
+  !    CALL MPI_DIMS_CREATE(num_procs, ndims, dims, ierr)
+  ! END IF
+  ! !
+  ! !periodicity in p
+  ! periods(1)=.FALSE.
+  ! !periodicity in kr
+  ! periods(2)=.FALSE.
+
+  ! CALL MPI_CART_CREATE(MPI_COMM_WORLD, ndims, dims, periods, reorder, comm0, ierr)
+
+  ! CALL MPI_COMM_RANK(comm0, me_0,  ierr)
+
+  ! CALL MPI_CART_COORDS(comm0,me_0,ndims,coords,ierr)
+
+  ! DO i=0,num_procs-1
+  !   CALL mpi_barrier(MPI_COMM_WORLD, ierr)
+  !   IF (my_id .EQ. i) THEN
+  !     WRITE(*,*) 'The coords of process ',me_0,' are: ', coords
+
+  !     IF( coords(1) .GT. 0 ) THEN
+  !       CALL MPI_CART_SHIFT(comm0, 0, -1, source , nghb_L , ierr)
+  !       CALL MPI_CART_COORDS(comm0,nghb_L,ndims,coords_L,ierr)
+  !       WRITE(*,*) coords_L,' are L from ', coords
+  !     ELSE
+  !       nghb_L = MPI_PROC_NULL
+  !     ENDIF
+
+
+  !     IF( coords(1) .LT. dims(1)-1 ) THEN
+  !       CALL MPI_CART_SHIFT(comm0, 0, +1, source , nghb_R , ierr)
+  !       CALL MPI_CART_COORDS(comm0,nghb_R,ndims,coords_R,ierr)
+  !       WRITE(*,*) coords_R,' are R from ', coords
+  !     ELSE
+  !       nghb_R = MPI_PROC_NULL
+  !     ENDIF
+  !   ENDIF
+  ! ENDDO
+
+  ! !
+  ! !  Partitions 2-dim cartesian topology of comm0 into 1-dim cartesian subgrids
+  ! !
+  ! CALL MPI_CART_SUB (comm0, (/.TRUE.,.FALSE./), commp, ierr)
+  ! CALL MPI_CART_SUB (comm0, (/.FALSE.,.TRUE./), commr, ierr)
+
+  ! CALL MPI_COMM_RANK(commp, me_p, ierr)
+  ! CALL MPI_COMM_RANK(commr, me_r, ierr)
+
 END SUBROUTINE ppinit
-- 
GitLab