From a213746492c1b4fffa3e0a934787b8e3e691cd8c Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Tue, 31 Jan 2023 17:22:00 +0100
Subject: [PATCH] Finished the diagnostic of N_a(p,j,z), cheaper way to obtain
 pj spectrum

---
 src/array_mod.F90      |  6 +++---
 src/diagnose.F90       | 22 ++++++++++++++--------
 src/processing_mod.F90 | 22 ++++++++++++----------
 3 files changed, 29 insertions(+), 21 deletions(-)

diff --git a/src/array_mod.F90 b/src/array_mod.F90
index 062c9deb..56fc14b9 100644
--- a/src/array_mod.F90
+++ b/src/array_mod.F90
@@ -72,9 +72,9 @@ MODULE array
   COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: Ne00
   COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: Ni00
 
-  ! Kinetic spectrum sum_kx,ky(|Napj(z)|^2), (ip,ij,iz)
-  REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: Nepjz
-  REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: Nipjz
+  ! Kinetic spectrum sum_kx,ky(|Napj(z)|^2), (ip,ij,iz) (should be real)
+  COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: Nepjz
+  COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: Nipjz
 
   ! particle density for electron and ions (iky,ikx,iz)
   COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: dens_e
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index cf1bfe27..ed872ad1 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -421,10 +421,10 @@ SUBROUTINE diagnose_3d
       Ni00(ikys:ikye,ikxs:ikxe,izs:ize) = moments_i(ip0_i,ij0_i,ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel)
     CALL write_field3d_kykxz(Ni00(ikys:ikye,ikxs:ikxe,izs:ize), 'Ni00')
 
-    ! CALL compute_Napjz_spectrum
-    ! IF(KIN_E) &
-    ! CALL write_field3d_pjz_e(Nepjz(ips_e:ipe_e,ijs_e:ije_e,izs:ize), 'Nepjz')
-    ! CALL write_field3d_pjz_i(Nipjz(ips_i:ipe_i,ijs_i:ije_i,izs:ize), 'Nipjz')
+    CALL compute_Napjz_spectrum
+    IF(KIN_E) &
+    CALL write_field3d_pjz_e(Nepjz(ips_e:ipe_e,ijs_e:ije_e,izs:ize), 'Nepjz')
+    CALL write_field3d_pjz_i(Nipjz(ips_i:ipe_i,ijs_i:ije_i,izs:ize), 'Nipjz')
   ENDIF
 
   !! Fuid moments
@@ -482,29 +482,35 @@ SUBROUTINE diagnose_3d
   END SUBROUTINE write_field3d_kykxz
 
   SUBROUTINE write_field3d_pjz_i(field, text)
+    USE parallel, ONLY : gather_pjz_i
     IMPLICIT NONE
-    REAL(dp), DIMENSION(ips_i:ipe_i,ijs_i:ije_i,izs:ize), INTENT(IN) :: field
+    COMPLEX(dp), DIMENSION(ips_i:ipe_i,ijs_i:ije_i,izs:ize), INTENT(IN) :: field
+    COMPLEX(dp), DIMENSION(1:pmaxi+1,1:jmaxi+1,1:Nz) :: field_full
     CHARACTER(*), INTENT(IN) :: text
     CHARACTER(LEN=50) :: dset_name
     WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var3d", TRIM(text), iframe3d
     IF (num_procs .EQ. 1) THEN ! no data distribution
       CALL putarr(fidres, dset_name, field(ips_i:ipe_i,ijs_i:ije_i,izs:ize), ionode=0)
     ELSE
-      CALL putarrnd(fidres, dset_name, field(ips_i:ipe_i,ijs_i:ije_i,izs:ize),  (/1, 0, 3/))
+      CALL gather_pjz_i(field(ips_i:ipe_i,ijs_i:ije_i,izs:ize),field_full(1:pmaxi+1,1:jmaxi+1,1:Nz))
+      CALL putarr(fidres, dset_name, field(1:pmaxi+1,1:jmaxi+1,1:Nz), ionode=0)
     ENDIF
     CALL attach(fidres, dset_name, "time", time)
   END SUBROUTINE write_field3d_pjz_i
 
   SUBROUTINE write_field3d_pjz_e(field, text)
+    USE parallel, ONLY : gather_pjz_e
     IMPLICIT NONE
-    REAL(dp), DIMENSION(ips_e:ipe_e,ijs_e:ije_e,izs:ize), INTENT(IN) :: field
+    COMPLEX(dp), DIMENSION(ips_e:ipe_e,ijs_e:ije_e,izs:ize), INTENT(IN) :: field
+    COMPLEX(dp), DIMENSION(1:pmaxe+1,1:jmaxe+1,1:Nz) :: field_full
     CHARACTER(*), INTENT(IN) :: text
     CHARACTER(LEN=50) :: dset_name
     WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var3d", TRIM(text), iframe3d
     IF (num_procs .EQ. 1) THEN ! no data distribution
       CALL putarr(fidres, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,izs:ize), ionode=0)
     ELSE
-      CALL putarrnd(fidres, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,izs:ize),  (/1, 0, 3/))
+      CALL gather_pjz_e(field(ips_e:ipe_e,ijs_e:ije_e,izs:ize),field_full(1:pmaxe+1,1:jmaxe+1,1:Nz))
+      CALL putarr(fidres, dset_name, field(1:pmaxi+1,1:jmaxi+1,1:Nz), ionode=0)
     ENDIF
     CALL attach(fidres, dset_name, "time", time)
   END SUBROUTINE write_field3d_pjz_e
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index b590e4fa..577e9f00 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -528,8 +528,8 @@ SUBROUTINE compute_Napjz_spectrum
   USE array,  ONLY : Nipjz, Nepjz
   USE time_integration, ONLY : updatetlevel
   IMPLICIT NONE
-  REAL(dp), DIMENSION(ips_e:ipe_e,ijs_e:ije_e,izs:ize) :: local_sum_e,global_sum_e, buffer_e
-  REAL(dp), DIMENSION(ips_i:ipe_i,ijs_i:ije_i,izs:ize) :: local_sum_i,global_sum_i, buffer_i
+  COMPLEX(dp), DIMENSION(ips_e:ipe_e,ijs_e:ije_e,izs:ize) :: local_sum_e,global_sum_e, buffer_e
+  COMPLEX(dp), DIMENSION(ips_i:ipe_i,ijs_i:ije_i,izs:ize) :: local_sum_i,global_sum_i, buffer_i
   INTEGER  :: i_, root, count
   root = 0
   ! Electron moments spectrum
@@ -538,8 +538,9 @@ SUBROUTINE compute_Napjz_spectrum
     local_sum_e = 0._dp
     DO ikx = ikxs,ikxe
       DO iky = ikys,ikye
-        local_sum_e(:,:,:)  = local_sum_e(:,:,:)  + &
-        moments_e(:,:,iky,ikx,:,updatetlevel) * CONJG(moments_e(:,:,iky,ikx,:,updatetlevel))
+        local_sum_e(ips_e:ipe_e,ijs_e:ije_e,izs:ize)  = local_sum_e(ips_e:ipe_e,ijs_e:ije_e,izs:ize)  + &
+        (moments_e(ips_e:ips_e,ijs_e:ije_e,iky,ikx,izs:ize,updatetlevel)&
+         * CONJG(moments_e(ips_e:ips_e,ijs_e:ije_e,iky,ikx,izs:ize,updatetlevel)))
       ENDDO
     ENDDO
     ! sum reduction
@@ -549,12 +550,12 @@ SUBROUTINE compute_Napjz_spectrum
     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_e, count , MPI_DOUBLE_PRECISION, root, 1234, comm_ky, ierr)
+            CALL MPI_SEND(buffer_e, count , MPI_DOUBLE_COMPLEX, 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_e, count , MPI_DOUBLE_PRECISION, i_, 1234, comm_ky, MPI_STATUS_IGNORE, ierr)
+                    CALL MPI_RECV(buffer_e, count , MPI_DOUBLE_COMPLEX, i_, 1234, comm_ky, MPI_STATUS_IGNORE, ierr)
                     global_sum_e = global_sum_e + buffer_e
             ENDDO
         ENDIF
@@ -568,8 +569,9 @@ SUBROUTINE compute_Napjz_spectrum
   local_sum_i = 0._dp
   DO ikx = ikxs,ikxe
     DO iky = ikys,ikye
-      local_sum_i(:,:,:)  = local_sum_i(:,:,:)  + &
-      moments_i(:,:,iky,ikx,:,updatetlevel) * CONJG(moments_i(:,:,iky,ikx,:,updatetlevel))
+      local_sum_i(ips_i:ipe_i,ijs_i:ije_i,izs:ize)  = local_sum_i(ips_i:ipe_i,ijs_i:ije_i,izs:ize)  + &
+      (moments_i(ips_i:ipe_i,ijs_i:ije_i,iky,ikx,izs:ize,updatetlevel) &
+      * CONJG(moments_i(ips_i:ipe_i,ijs_i:ije_i,iky,ikx,izs:ize,updatetlevel)))
     ENDDO
   ENDDO
   ! sum reduction
@@ -579,12 +581,12 @@ SUBROUTINE compute_Napjz_spectrum
   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_i, count , MPI_DOUBLE_PRECISION, root, 1234, comm_ky, ierr)
+          CALL MPI_SEND(buffer_i, count , MPI_DOUBLE_COMPLEX, 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_i, count , MPI_DOUBLE_PRECISION, i_, 1234, comm_ky, MPI_STATUS_IGNORE, ierr)
+                  CALL MPI_RECV(buffer_i, count , MPI_DOUBLE_COMPLEX, i_, 1234, comm_ky, MPI_STATUS_IGNORE, ierr)
                   global_sum_i = global_sum_i + buffer_i
           ENDDO
       ENDIF
-- 
GitLab