From c400defa9af00c51f41696de159796abf01223d5 Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Wed, 11 May 2022 10:02:19 +0200
Subject: [PATCH] output napjz spectrum and corrected clos 1 bug
---
src/advance_field.F90 | 5 +-
src/array_mod.F90 | 4 ++
src/closure_mod.F90 | 4 +-
src/diagnose.F90 | 79 ++++++++++++++++++---------
src/memory.F90 | 2 +
src/nonlinear_mod.F90 | 118 ++++++++++++++++++++++-------------------
src/processing_mod.F90 | 91 +++++++++++++++++++++++++++----
7 files changed, 207 insertions(+), 96 deletions(-)
diff --git a/src/advance_field.F90 b/src/advance_field.F90
index ae48815e..dcad3c44 100644
--- a/src/advance_field.F90
+++ b/src/advance_field.F90
@@ -31,7 +31,8 @@ CONTAINS
DO ip=ips_e,ipe_e
p_int = parray_e(ip)
DO ij=ijs_e,ije_e
- IF((CLOS .NE. 1) .OR. (ip-1+2*(ij-1)+1 .LE. dmaxe))&
+ j_int = jarray_e(ij)
+ IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxe))&
CALL advance_field(moments_e(ip,ij,ikys:ikye,ikxs:ikxe,izs:ize,:), moments_rhs_e(ip,ij,ikys:ikye,ikxs:ikxe,izs:ize,:))
ENDDO
ENDDO
@@ -40,7 +41,7 @@ CONTAINS
p_int = parray_i(ip)
DO ij=ijs_i,ije_i
j_int = jarray_i(ij)
- IF((CLOS .NE. 1) .OR. (ip-1+2*(ij-1)+1 .LE. dmaxi))&
+ IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxi))&
CALL advance_field(moments_i(ip,ij,ikys:ikye,ikxs:ikxe,izs:ize,:), moments_rhs_i(ip,ij,ikys:ikye,ikxs:ikxe,izs:ize,:))
ENDDO
ENDDO
diff --git a/src/array_mod.F90 b/src/array_mod.F90
index 419d3c6c..96485a88 100644
--- a/src/array_mod.F90
+++ b/src/array_mod.F90
@@ -68,6 +68,10 @@ 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
+
! particle density for electron and ions (ikx,iky,iz)
COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: dens_e
COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: dens_i
diff --git a/src/closure_mod.F90 b/src/closure_mod.F90
index 5d48206c..ad7eac91 100644
--- a/src/closure_mod.F90
+++ b/src/closure_mod.F90
@@ -33,14 +33,14 @@ SUBROUTINE apply_closure_model
IF(KIN_E) THEN
DO ip = ipgs_e,ipge_e
DO ij = ijgs_e,ijge_e
- IF ( parray_e(ip)+2*jarray_e(ip) .GT. dmaxe) &
+ IF ( parray_e(ip)+2*jarray_e(ij) .GT. dmaxe) &
moments_e(ip,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) = 0._dp
ENDDO
ENDDO
ENDIF
DO ip = ipgs_i,ipge_i
DO ij = ijgs_i,ijge_i
- IF ( parray_i(ip)+2*jarray_i(ip) .GT. dmaxi) &
+ IF ( parray_i(ip)+2*jarray_i(ij) .GT. dmaxi) &
moments_i(ip,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) = 0._dp
ENDDO
ENDDO
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index 925f515b..ced695af 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -134,6 +134,9 @@ SUBROUTINE diagnose(kstep)
IF(KIN_E)&
CALL creatg(fidres, "/data/var3d/Ne00", "Ne00")
CALL creatg(fidres, "/data/var3d/Ni00", "Ni00")
+ IF(KIN_E)&
+ CALL creatg(fidres, "/data/var3d/Nepjz", "Nepjz")
+ CALL creatg(fidres, "/data/var3d/Nipjz", "Nipjz")
ENDIF
IF (write_dens) THEN
@@ -426,7 +429,7 @@ END SUBROUTINE diagnose_2d
SUBROUTINE diagnose_3d
USE basic
- USE futils, ONLY: append, getatt, attach, putarrnd
+ USE futils, ONLY: append, getatt, attach, putarrnd, putarr
USE fields
USE array
USE grid, ONLY: ikxs,ikxe, ikys,ikye, Nkx, Nky, local_nkx, ikx, iky, ips_e, ips_i
@@ -446,21 +449,24 @@ SUBROUTINE diagnose_3d
iframe3d=iframe3d+1
CALL attach(fidres,"/data/var3d/" , "frames", iframe3d)
- IF (write_phi) CALL write_field3d(phi (ikys:ikye,ikxs:ikxe,izs:ize), 'phi')
+ IF (write_phi) CALL write_field3d_kykxz(phi (ikys:ikye,ikxs:ikxe,izs:ize), 'phi')
IF (write_Na00) THEN
IF(KIN_E)THEN
IF (ips_e .EQ. 1) THEN
Ne00(ikys:ikye,ikxs:ikxe,izs:ize) = moments_e(ips_e,1,ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel)
ENDIF
- ! CALL manual_3D_bcast(Ne00(ikys:ikye,ikxs:ikxe,izs:ize))
- CALL write_field3d(Ne00(ikys:ikye,ikxs:ikxe,izs:ize), 'Ne00')
+ CALL write_field3d_kykxz(Ne00(ikys:ikye,ikxs:ikxe,izs:ize), 'Ne00')
ENDIF
IF (ips_i .EQ. 1) THEN
Ni00(ikys:ikye,ikxs:ikxe,izs:ize) = moments_i(ips_e,1,ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel)
ENDIF
- ! CALL manual_3D_bcast(Ni00(ikys:ikye,ikxs:ikxe,izs:ize))
- CALL write_field3d(Ni00(ikys:ikye,ikxs:ikxe,izs:ize), 'Ni00')
+ 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')
ENDIF
!! Fuid moments
@@ -469,44 +475,38 @@ SUBROUTINE diagnose_3d
IF (write_dens) THEN
IF(KIN_E)&
- CALL write_field3d(dens_e(ikys:ikye,ikxs:ikxe,izs:ize), 'dens_e')
- CALL write_field3d(dens_i(ikys:ikye,ikxs:ikxe,izs:ize), 'dens_i')
+ CALL write_field3d_kykxz(dens_e(ikys:ikye,ikxs:ikxe,izs:ize), 'dens_e')
+ CALL write_field3d_kykxz(dens_i(ikys:ikye,ikxs:ikxe,izs:ize), 'dens_i')
ENDIF
IF (write_fvel) THEN
IF(KIN_E)&
- CALL write_field3d(upar_e(ikys:ikye,ikxs:ikxe,izs:ize), 'upar_e')
- CALL write_field3d(upar_i(ikys:ikye,ikxs:ikxe,izs:ize), 'upar_i')
+ CALL write_field3d_kykxz(upar_e(ikys:ikye,ikxs:ikxe,izs:ize), 'upar_e')
+ CALL write_field3d_kykxz(upar_i(ikys:ikye,ikxs:ikxe,izs:ize), 'upar_i')
IF(KIN_E)&
- CALL write_field3d(uper_e(ikys:ikye,ikxs:ikxe,izs:ize), 'uper_e')
- CALL write_field3d(uper_i(ikys:ikye,ikxs:ikxe,izs:ize), 'uper_i')
+ CALL write_field3d_kykxz(uper_e(ikys:ikye,ikxs:ikxe,izs:ize), 'uper_e')
+ CALL write_field3d_kykxz(uper_i(ikys:ikye,ikxs:ikxe,izs:ize), 'uper_i')
ENDIF
IF (write_temp) THEN
IF(KIN_E)&
- CALL write_field3d(Tpar_e(ikys:ikye,ikxs:ikxe,izs:ize), 'Tpar_e')
- CALL write_field3d(Tpar_i(ikys:ikye,ikxs:ikxe,izs:ize), 'Tpar_i')
+ CALL write_field3d_kykxz(Tpar_e(ikys:ikye,ikxs:ikxe,izs:ize), 'Tpar_e')
+ CALL write_field3d_kykxz(Tpar_i(ikys:ikye,ikxs:ikxe,izs:ize), 'Tpar_i')
IF(KIN_E)&
- CALL write_field3d(Tper_e(ikys:ikye,ikxs:ikxe,izs:ize), 'Tper_e')
- CALL write_field3d(Tper_i(ikys:ikye,ikxs:ikxe,izs:ize), 'Tper_i')
+ CALL write_field3d_kykxz(Tper_e(ikys:ikye,ikxs:ikxe,izs:ize), 'Tper_e')
+ CALL write_field3d_kykxz(Tper_i(ikys:ikye,ikxs:ikxe,izs:ize), 'Tper_i')
IF(KIN_E)&
- CALL write_field3d(temp_e(ikys:ikye,ikxs:ikxe,izs:ize), 'temp_e')
- CALL write_field3d(temp_i(ikys:ikye,ikxs:ikxe,izs:ize), 'temp_i')
+ CALL write_field3d_kykxz(temp_e(ikys:ikye,ikxs:ikxe,izs:ize), 'temp_e')
+ CALL write_field3d_kykxz(temp_i(ikys:ikye,ikxs:ikxe,izs:ize), 'temp_i')
ENDIF
CONTAINS
- SUBROUTINE write_field3d(field, text)
- USE futils, ONLY: attach, putarr
- USE grid, ONLY: ikxs,ikxe, ikys,ikye, Nkx, Nky, local_nkx
- USE prec_const
- USE basic, ONLY : comm_ky, num_procs_p, rank_p
+ SUBROUTINE write_field3d_kykxz(field, text)
IMPLICIT NONE
-
COMPLEX(dp), DIMENSION(ikys:ikye,ikxs:ikxe, izs:ize), INTENT(IN) :: field
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(ikys:ikye,ikxs:ikxe, izs:ize), ionode=0)
@@ -514,8 +514,35 @@ CONTAINS
CALL putarrnd(fidres, dset_name, field(ikys:ikye,ikxs:ikxe, izs:ize), (/1, 1, 3/))
ENDIF
CALL attach(fidres, dset_name, "time", time)
+ END SUBROUTINE write_field3d_kykxz
+
+ SUBROUTINE write_field3d_pjz_i(field, text)
+ IMPLICIT NONE
+ REAL(dp), DIMENSION(ips_i:ipe_i,ijs_i:ije_i,izs:ize), INTENT(IN) :: field
+ 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/))
+ ENDIF
+ CALL attach(fidres, dset_name, "time", time)
+ END SUBROUTINE write_field3d_pjz_i
- END SUBROUTINE write_field3d
+ SUBROUTINE write_field3d_pjz_e(field, text)
+ IMPLICIT NONE
+ REAL(dp), DIMENSION(ips_e:ipe_e,ijs_e:ije_e,izs:ize), INTENT(IN) :: field
+ 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/))
+ ENDIF
+ CALL attach(fidres, dset_name, "time", time)
+ END SUBROUTINE write_field3d_pjz_e
END SUBROUTINE diagnose_3d
diff --git a/src/memory.F90 b/src/memory.F90
index a1ebae20..a3837b7c 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -29,6 +29,7 @@ SUBROUTINE memory
CALL allocate_array( temp_e, ikys,ikye, ikxs,ikxe, izs,ize)
CALL allocate_array( Kernel_e, ijgs_e,ijge_e, ikys,ikye, ikxs,ikxe, izgs,izge, 0,1)
CALL allocate_array( moments_e, ipgs_e,ipge_e, ijgs_e,ijge_e, ikys,ikye, ikxs,ikxe, izgs,izge, 1,ntimelevel )
+ CALL allocate_array( Nepjz, ips_e,ipe_e, ijs_e,ije_e, izs,ize)
CALL allocate_array( moments_rhs_e, ips_e,ipe_e, ijs_e,ije_e, ikys,ikye, ikxs,ikxe, izs,ize, 1,ntimelevel )
CALL allocate_array( nadiab_moments_e, ipgs_e,ipge_e, ijgs_e,ijge_e, ikys,ikye, ikxs,ikxe, izgs,izge)
CALL allocate_array( ddz_nepj, ipgs_e,ipge_e, ijgs_e,ijge_e, ikys,ikye, ikxs,ikxe, izgs,izge)
@@ -64,6 +65,7 @@ SUBROUTINE memory
CALL allocate_array( temp_i, ikys,ikye, ikxs,ikxe, izs,ize)
CALL allocate_array( Kernel_i, ijgs_i,ijge_i, ikys,ikye, ikxs,ikxe, izgs,izge, 0,1)
CALL allocate_array( moments_i, ipgs_i,ipge_i, ijgs_i,ijge_i, ikys,ikye, ikxs,ikxe, izgs,izge, 1,ntimelevel )
+ CALL allocate_array( Nipjz, ips_i,ipe_i, ijs_i,ije_i, izs,ize)
CALL allocate_array( moments_rhs_i, ips_i,ipe_i, ijs_i,ije_i, ikys,ikye, ikxs,ikxe, izs,ize, 1,ntimelevel )
CALL allocate_array( nadiab_moments_i, ipgs_i,ipge_i, ijgs_i,ijge_i, ikys,ikye, ikxs,ikxe, izgs,izge)
CALL allocate_array( ddz_nipj, ipgs_i,ipge_i, ijgs_i,ijge_i, ikys,ikye, ikxs,ikxe, izgs,izge)
diff --git a/src/nonlinear_mod.F90 b/src/nonlinear_mod.F90
index 78ffce96..4d83dfc1 100644
--- a/src/nonlinear_mod.F90
+++ b/src/nonlinear_mod.F90
@@ -71,39 +71,42 @@ SUBROUTINE compute_nonlinear
zloope: DO iz = izs,ize
ploope: DO ip = ips_e,ipe_e ! Loop over Hermite moments
eo = MODULO(parray_e(ip),2)
+ p_int = parray_e(ip)
jloope: DO ij = ijs_e, ije_e ! Loop over Laguerre moments
j_int=jarray_e(ij)
- ! Set non linear sum truncation
- IF (NL_CLOS .EQ. -2) THEN
- nmax = Jmaxe
- ELSEIF (NL_CLOS .EQ. -1) THEN
- nmax = Jmaxe-(ij-1)
- ELSE
- nmax = NL_CLOS
- ENDIF
- bracket_sum_r = 0._dp ! initialize sum over real nonlinear term
- nloope: DO in = 1,nmax+1 ! Loop over laguerre for the sum
- ! First convolution terms
- F_cmpx(:,:) = phi(:,:,iz) * kernel_e(in, :, :, iz, eo)
- ! Second convolution terms
- G_cmpx(:,:) = 0._dp ! initialization of the sum
- smax = MIN( (in-1)+(ij-1), jmaxe );
- DO is = 1, smax+1 ! sum truncation on number of moments
- G_cmpx(:,:) = G_cmpx(:,:) + &
- dnjs(in,ij,is) * moments_e(ip,is,:,:,iz,updatetlevel)
- ENDDO
- !/!\ this function add its result to bracket_sum_r (hard to read sorry) /!\
- CALL poisson_bracket_and_sum(F_cmpx,G_cmpx)
- ENDDO nloope
+ IF((CLOS .EQ. 1) .AND. (p_int+2*j_int .LE. dmaxe)) THEN !compute
+ ! Set non linear sum truncation
+ IF (NL_CLOS .EQ. -2) THEN
+ nmax = Jmaxe
+ ELSEIF (NL_CLOS .EQ. -1) THEN
+ nmax = Jmaxe-j_int
+ ELSE
+ nmax = NL_CLOS
+ ENDIF
+ bracket_sum_r = 0._dp ! initialize sum over real nonlinear term
+ nloope: DO in = 1,nmax+1 ! Loop over laguerre for the sum
+ ! First convolution terms
+ F_cmpx(ikxs:ikxe,ikxs:ikxe) = phi(ikxs:ikxe,ikxs:ikxe,iz) * kernel_e(in, ikxs:ikxe,ikxs:ikxe, iz, eo)
+ ! Second convolution terms
+ G_cmpx(ikxs:ikxe,ikxs:ikxe) = 0._dp ! initialization of the sum
+ smax = MIN( (in-1)+(ij-1), jmaxe );
+ DO is = 1, smax+1 ! sum truncation on number of moments
+ G_cmpx(ikxs:ikxe,ikxs:ikxe) = G_cmpx(ikxs:ikxe,ikxs:ikxe) + &
+ dnjs(in,ij,is) * moments_e(ip,is,ikxs:ikxe,ikxs:ikxe,iz,updatetlevel)
+ ENDDO
+ !/!\ this function add its result to bracket_sum_r (hard to read sorry) /!\
+ CALL poisson_bracket_and_sum(F_cmpx,G_cmpx)
+ ENDDO nloope
- ! Put the real nonlinear product into k-space
- call fftw_mpi_execute_dft_r2c(planf, bracket_sum_r, bracket_sum_c)
- ! Retrieve convolution in input format
- DO ikx = ikxs, ikxe
+ ! Put the real nonlinear product into k-space
+ call fftw_mpi_execute_dft_r2c(planf, bracket_sum_r, bracket_sum_c)
+ ! Retrieve convolution in input format
DO iky = ikys, ikye
- Sepj(ip,ij,iky,ikx,iz) = bracket_sum_c(ikx,iky-local_nky_offset)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
+ Sepj(ip,ij,iky,ikxs:ikxe,iz) = bracket_sum_c(ikxs:ikxe,iky-local_nky_offset)*AA_x(ikxs:ikxe)*AA_y(iky) !Anti aliasing filter
ENDDO
- ENDDO
+ ELSE
+ Sepj(ip,ij,:,:,iz) = 0._dp
+ ENDIF
ENDDO jloope
ENDDO ploope
ENDDO zloope
@@ -113,39 +116,42 @@ ENDIF
!!!!!!!!!!!!!!!!!!!! ION non linear term computation (Sipj)!!!!!!!!!!
zloopi: DO iz = izs,ize
ploopi: DO ip = ips_i,ipe_i ! Loop over Hermite moments
+ p_int = parray_i(ip)
eo = MODULO(parray_i(ip),2)
jloopi: DO ij = ijs_i, ije_i ! Loop over Laguerre moments
j_int=jarray_i(ij)
- ! Set non linear sum truncation
- IF (NL_CLOS .EQ. -2) THEN
- nmax = Jmaxi
- ELSEIF (NL_CLOS .EQ. -1) THEN
- nmax = Jmaxi-(ij-1)
- ELSE
- nmax = NL_CLOS
- ENDIF
- bracket_sum_r = 0._dp ! initialize sum over real nonlinear term
- nloopi: DO in = 1,nmax+1 ! Loop over laguerre for the sum
- ! First convolution terms
- F_cmpx(:,:) = phi(:,:,iz) * kernel_i(in, :, :, iz, eo)
- ! Second convolution terms
- G_cmpx(:,:) = 0._dp ! initialization of the sum
- smax = MIN( (in-1)+(ij-1), jmaxi );
- DO is = 1, smax+1 ! sum truncation on number of moments
- G_cmpx(:,:) = G_cmpx(:,:) + &
- dnjs(in,ij,is) * moments_i(ip,is,:,:,iz,updatetlevel)
- ENDDO
- !/!\ this function add its result to bracket_sum_r (hard to read sorry) /!\
- CALL poisson_bracket_and_sum(F_cmpx,G_cmpx)
- ENDDO nloopi
- ! Put the real nonlinear product into k-space
- call fftw_mpi_execute_dft_r2c(planf, bracket_sum_r, bracket_sum_c)
- ! Retrieve convolution in input format
- DO ikx = ikxs, ikxe
+ IF((CLOS .EQ. 1) .AND. (p_int+2*j_int .LE. dmaxi)) THEN !compute
+ ! Set non linear sum truncation
+ IF (NL_CLOS .EQ. -2) THEN
+ nmax = Jmaxi
+ ELSEIF (NL_CLOS .EQ. -1) THEN
+ nmax = Jmaxi-j_int
+ ELSE
+ nmax = NL_CLOS
+ ENDIF
+ bracket_sum_r = 0._dp ! initialize sum over real nonlinear term
+ nloopi: DO in = 1,nmax+1 ! Loop over laguerre for the sum
+ ! First convolution terms
+ F_cmpx(ikys:ikye,ikxs:ikxe) = phi(ikys:ikye,ikxs:ikxe,iz) * kernel_i(in, ikys:ikye,ikxs:ikxe, iz, eo)
+ ! Second convolution terms
+ G_cmpx(ikys:ikye,ikxs:ikxe) = 0._dp ! initialization of the sum
+ smax = MIN( (in-1)+(ij-1), jmaxi );
+ DO is = 1, smax+1 ! sum truncation on number of moments
+ G_cmpx(ikys:ikye,ikxs:ikxe) = G_cmpx(ikys:ikye,ikxs:ikxe) + &
+ dnjs(in,ij,is) * moments_i(ip,is,ikys:ikye,ikxs:ikxe,iz,updatetlevel)
+ ENDDO
+ !/!\ this function add its result to bracket_sum_r (hard to read sorry) /!\
+ CALL poisson_bracket_and_sum(F_cmpx,G_cmpx)
+ ENDDO nloopi
+ ! Put the real nonlinear product into k-space
+ call fftw_mpi_execute_dft_r2c(planf, bracket_sum_r, bracket_sum_c)
+ ! Retrieve convolution in input format
DO iky = ikys, ikye
- Sipj(ip,ij,iky,ikx,iz) = bracket_sum_c(ikx,iky-local_nky_offset)*AA_x(ikx)*AA_y(iky)
+ Sipj(ip,ij,iky,ikxs:ikxe,iz) = bracket_sum_c(ikxs:ikxe,iky-local_nky_offset)*AA_x(ikxs:ikxe)*AA_y(iky)
ENDDO
- ENDDO
+ ELSE
+ Sipj(ip,ij,:,:,iz) = 0._dp
+ ENDIF
ENDDO jloopi
ENDDO ploopi
ENDDO zloopi
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index d98334cb..5a47ed6f 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -1,11 +1,7 @@
MODULE processing
- ! contains the Hermite-Laguerre collision operators. Solved using COSOlver.
USE basic
USE prec_const
USE grid
- USE geometry
- USE utility
- USE calculus
implicit none
REAL(dp), PUBLIC, PROTECTED :: pflux_ri, gflux_ri
@@ -15,14 +11,16 @@ MODULE processing
PUBLIC :: compute_density, compute_upar, compute_uperp
PUBLIC :: compute_Tpar, compute_Tperp, compute_fluid_moments
PUBLIC :: compute_radial_ion_transport, compute_radial_heatflux
+ PUBLIC :: compute_Napjz_spectrum
CONTAINS
-! 1D diagnostic to compute the average radial particle transport <n_i v_ExB>_r
+! 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
+ 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)
@@ -91,13 +89,14 @@ SUBROUTINE compute_radial_ion_transport
! 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_i v_ExB>_r
+! 1D diagnostic to compute the average radial particle transport <T_i v_ExB_x>_xyz
SUBROUTINE compute_radial_heatflux
USE fields, ONLY : moments_i, moments_e, phi
USE array, ONLY : kernel_e, kernel_i
- USE geometry, ONLY : Jacobian
+ USE geometry, ONLY : Jacobian, iInt_Jacobian
USE time_integration, ONLY : updatetlevel
- USE model, ONLY : qe_taue, qi_taui, KIN_E
+ USE model, ONLY : qe_taue, qi_taui, KIN_E
+ USE calculus, ONLY : simpson_rule_z
IMPLICIT NONE
COMPLEX(dp) :: hflux_local, integral
REAL(dp) :: ky_, buffer(1:2), j_dp
@@ -228,7 +227,7 @@ SUBROUTINE compute_nadiab_moments_z_gradients_and_interp
ENDDO
ENDIF
-!------------- INTERP AND GRADIENTS ALONG Z ----------------------------------
+ !------------- INTERP AND GRADIENTS ALONG Z ----------------------------------
IF (KIN_E) THEN
DO ikx = ikxs,ikxe
@@ -271,6 +270,78 @@ SUBROUTINE compute_nadiab_moments_z_gradients_and_interp
tc_process = tc_process + (t1_process - t0_process)
END SUBROUTINE compute_nadiab_moments_z_gradients_and_interp
+SUBROUTINE compute_Napjz_spectrum
+ USE fields, ONLY : moments_e, moments_i
+ USE model, ONLY : KIN_E
+ 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
+ INTEGER :: i_, world_rank, world_size, root, count
+ root = 0
+ ! Electron moments spectrum
+ IF (KIN_E) THEN
+ ! build local sum
+ 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))
+ ENDDO
+ ENDDO
+ ! sum reduction
+ buffer_e = local_sum_e
+ global_sum_e = 0._dp
+ count = (ipe_e-ips_e+1)*(ije_e-ijs_e+1)*(ize-izs+1)
+ 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)
+ 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)
+ global_sum_e = global_sum_e + buffer_e
+ ENDDO
+ ENDIF
+ ELSE
+ global_sum_e = local_sum_e
+ ENDIF
+ Nepjz = global_sum_e
+ ENDIF
+ ! Ion moment spectrum
+ ! build local sum
+ 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))
+ ENDDO
+ ENDDO
+ ! sum reduction
+ buffer_i = local_sum_i
+ global_sum_i = 0._dp
+ count = (ipe_i-ips_i+1)*(ije_i-ijs_i+1)*(ize-izs+1)
+ 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)
+ 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)
+ global_sum_i = global_sum_i + buffer_i
+ ENDDO
+ ENDIF
+ ELSE
+ global_sum_i = local_sum_i
+ ENDIF
+ Nipjz = global_sum_i
+END SUBROUTINE compute_Napjz_spectrum
+
!_____________________________________________________________________________!
!!!!! FLUID MOMENTS COMPUTATIONS !!!!!
! Compute the 2D particle density for electron and ions (sum over Laguerre)
--
GitLab