From 9ef07a7a803b1cfd5d05bac20f99510b477bb609 Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Tue, 15 Jun 2021 14:05:34 +0200
Subject: [PATCH] Added density and temperature diagnostics
---
src/diagnose.F90 | 93 ++++++++++++-------------------------
src/diagnostics_par_mod.F90 | 4 +-
src/memory.F90 | 49 +++++++++++--------
src/processing_mod.F90 | 84 ++++++++++++++++++++++++++++++---
4 files changed, 140 insertions(+), 90 deletions(-)
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index 33e8dbee..21f0a64f 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -126,6 +126,16 @@ SUBROUTINE diagnose(kstep)
CALL creatg(fidres, "/data/var2d/Ni00", "Ni00")
ENDIF
+ IF (write_dens) THEN
+ CALL creatg(fidres, "/data/var2d/dens_e", "dens_e")
+ CALL creatg(fidres, "/data/var2d/dens_i", "dens_i")
+ ENDIF
+
+ IF (write_temp) THEN
+ CALL creatg(fidres, "/data/var2d/temp_e", "temp_e")
+ CALL creatg(fidres, "/data/var2d/temp_i", "temp_i")
+ ENDIF
+
IF (cstep==0) THEN
iframe2d=0
ENDIF
@@ -184,6 +194,8 @@ SUBROUTINE diagnose(kstep)
CALL attach(fidres, TRIM(str), "write_Na00",write_Na00)
CALL attach(fidres, TRIM(str), "write_Napj",write_Napj)
CALL attach(fidres, TRIM(str), "write_Sapj",write_Sapj)
+ CALL attach(fidres, TRIM(str), "write_dens",write_dens)
+ CALL attach(fidres, TRIM(str), "write_temp",write_temp)
CALL grid_outputinputs(fidres, str)
@@ -227,7 +239,7 @@ SUBROUTINE diagnose(kstep)
! Terminal info
IF (MOD(cstep, INT(1.0/dt)) == 0 .AND. (my_id .EQ. 0)) THEN
- WRITE(*,"(F5.0,A,F5.0)") time,"/",tmax
+ WRITE(*,"(F6.0,A,F6.0)") time,"/",tmax
ENDIF
! 2.1 0d history arrays
@@ -321,11 +333,13 @@ SUBROUTINE diagnose_2d
USE basic
USE futils, ONLY: append, getatt, attach, putarrnd
USE fields
- USE array, ONLY: Ne00, Ni00
+ USE array, ONLY: Ne00, Ni00, dens_e, dens_i, temp_e, temp_i
USE grid, ONLY: ikrs,ikre, ikzs,ikze, nkr, nkz, local_nkr, ikr, ikz, ips_e, ips_i
USE time_integration
USE diagnostics_par
USE prec_const
+ USE processing
+
IMPLICIT NONE
COMPLEX(dp) :: buffer(ikrs:ikre,ikzs:ikze)
@@ -345,72 +359,25 @@ SUBROUTINE diagnose_2d
Ni00(ikrs:ikre,ikzs:ikze) = moments_i(ips_e,1,ikrs:ikre,ikzs:ikze,updatetlevel)
ENDIF
- root = 0
- !!!!! This is a manual way to do MPI_BCAST !!!!!!!!!!!
- CALL MPI_COMM_RANK(comm_p,world_rank,ierr)
- CALL MPI_COMM_SIZE(comm_p,world_size,ierr)
-
- IF (world_size .GT. 1) THEN
- !! Broadcast phi to the other processes on the same k range (communicator along p)
- IF (world_rank .EQ. root) THEN
- ! Fill the buffer
- DO ikr = ikrs,ikre
- DO ikz = ikzs,ikze
- buffer(ikr,ikz) = Ne00(ikr,ikz)
- ENDDO
- ENDDO
- ! Send it to all the other processes
- DO i_ = 0,num_procs_p-1
- IF (i_ .NE. world_rank) &
- CALL MPI_SEND(buffer, local_nkr * nkz , MPI_DOUBLE_COMPLEX, i_, 0, comm_p, ierr)
- ENDDO
- ELSE
- ! Recieve buffer from root
- CALL MPI_RECV(buffer, local_nkr * nkz , MPI_DOUBLE_COMPLEX, root, 0, comm_p, MPI_STATUS_IGNORE, ierr)
- ! Write it in phi
- DO ikr = ikrs,ikre
- DO ikz = ikzs,ikze
- Ne00(ikr,ikz) = buffer(ikr,ikz)
- ENDDO
- ENDDO
- ENDIF
- ENDIF
-
+ CALL manual_2D_bcast(Ne00(ikrs:ikre,ikzs:ikze))
CALL write_field2d(Ne00(ikrs:ikre,ikzs:ikze), 'Ne00')
- !!!!! This is a manual way to do MPI_BCAST !!!!!!!!!!!
- CALL MPI_COMM_RANK(comm_p,world_rank,ierr)
- CALL MPI_COMM_SIZE(comm_p,world_size,ierr)
-
- IF (world_size .GT. 1) THEN
- !! Broadcast phi to the other processes on the same k range (communicator along p)
- IF (world_rank .EQ. root) THEN
- ! Fill the buffer
- DO ikr = ikrs,ikre
- DO ikz = ikzs,ikze
- buffer(ikr,ikz) = Ni00(ikr,ikz)
- ENDDO
- ENDDO
- ! Send it to all the other processes
- DO i_ = 0,num_procs_p-1
- IF (i_ .NE. world_rank) &
- CALL MPI_SEND(buffer, local_nkr * nkz , MPI_DOUBLE_COMPLEX, i_, 0, comm_p, ierr)
- ENDDO
- ELSE
- ! Recieve buffer from root
- CALL MPI_RECV(buffer, local_nkr * nkz , MPI_DOUBLE_COMPLEX, root, 0, comm_p, MPI_STATUS_IGNORE, ierr)
- ! Write it in phi
- DO ikr = ikrs,ikre
- DO ikz = ikzs,ikze
- Ni00(ikr,ikz) = buffer(ikr,ikz)
- ENDDO
- ENDDO
- ENDIF
- ENDIF
-
+ CALL manual_2D_bcast(Ni00(ikrs:ikre,ikzs:ikze))
CALL write_field2d(Ni00(ikrs:ikre,ikzs:ikze), 'Ni00')
ENDIF
+ IF (write_dens) THEN
+ CALL compute_density
+ CALL write_field2d(dens_e(ikrs:ikre,ikzs:ikze), 'dens_e')
+ CALL write_field2d(dens_i(ikrs:ikre,ikzs:ikze), 'dens_i')
+ ENDIF
+
+ IF (write_temp) THEN
+ CALL compute_temperature
+ CALL write_field2d(temp_e(ikrs:ikre,ikzs:ikze), 'temp_e')
+ CALL write_field2d(temp_i(ikrs:ikre,ikzs:ikze), 'temp_i')
+ ENDIF
+
CONTAINS
SUBROUTINE write_field2d(field, text)
diff --git a/src/diagnostics_par_mod.F90 b/src/diagnostics_par_mod.F90
index 0f743896..011e09ed 100644
--- a/src/diagnostics_par_mod.F90
+++ b/src/diagnostics_par_mod.F90
@@ -7,8 +7,9 @@ MODULE diagnostics_par
LOGICAL, PUBLIC, PROTECTED :: write_doubleprecision = .FALSE.
LOGICAL, PUBLIC, PROTECTED :: write_gamma
- LOGICAL, PUBLIC, PROTECTED :: write_phi, write_Na00
+ LOGICAL, PUBLIC, PROTECTED :: write_phi, write_Na00
LOGICAL, PUBLIC, PROTECTED :: write_Napj, write_Sapj
+ LOGICAL, PUBLIC, PROTECTED :: write_dens, write_temp
INTEGER, PUBLIC, PROTECTED :: nsave_0d , nsave_1d , nsave_2d , nsave_5d, nsave_cp
@@ -36,6 +37,7 @@ CONTAINS
NAMELIST /OUTPUT_PAR/ nsave_0d , nsave_1d , nsave_2d , nsave_5d, nsave_cp
NAMELIST /OUTPUT_PAR/ write_doubleprecision, write_gamma, write_phi
NAMELIST /OUTPUT_PAR/ write_Na00, write_Napj, write_Sapj
+ NAMELIST /OUTPUT_PAR/ write_dens, write_temp
NAMELIST /OUTPUT_PAR/ resfile0, rstfile0, job2load
READ(lu_in,output_par)
diff --git a/src/memory.F90 b/src/memory.F90
index 8215f48a..4458499a 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -10,27 +10,47 @@ SUBROUTINE memory
USE prec_const
IMPLICIT NONE
-
- ! Moments with ghost degrees for p+2 p-2, j+1, j-1 truncations
- CALL allocate_array( moments_e, ipsg_e,ipeg_e, ijsg_e,ijeg_e, ikrs,ikre, ikzs,ikze, 1,ntimelevel )
- CALL allocate_array( moments_i, ipsg_i,ipeg_i, ijsg_i,ijeg_i, ikrs,ikre, ikzs,ikze, 1,ntimelevel )
-
- ! Moments right-hand-side (contains linear part of hierarchy)
- CALL allocate_array( moments_rhs_e, ips_e,ipe_e, ijs_e,ije_e, ikrs,ikre, ikzs,ikze, 1,ntimelevel )
- CALL allocate_array( moments_rhs_i, ips_i,ipe_i, ijs_i,ije_i, ikrs,ikre, ikzs,ikze, 1,ntimelevel )
-
+ !___________________ 2D ARRAYS __________________________
! Electrostatic potential
CALL allocate_array(phi, ikrs,ikre, ikzs,ikze)
+ !! Diagnostics arrays
! Gyrocenter density *for 2D output*
CALL allocate_array(Ne00, ikrs,ikre, ikzs,ikze)
CALL allocate_array(Ni00, ikrs,ikre, ikzs,ikze)
+ ! particle density *for 2D output*
+ CALL allocate_array(dens_e, ikrs,ikre, ikzs,ikze)
+ CALL allocate_array(dens_i, ikrs,ikre, ikzs,ikze)
+ ! particle temperature *for 2D output*
+ CALL allocate_array(temp_e, ikrs,ikre, ikzs,ikze)
+ CALL allocate_array(temp_i, ikrs,ikre, ikzs,ikze)
+ !___________________ 3D ARRAYS __________________________
+ !! FLR kernels functions
! Kernel evaluation from j= -1 to jmax+1 for truncation
CALL allocate_array(Kernel_e, ijsg_e,ijeg_e, ikrs,ikre, ikzs,ikze)
CALL allocate_array(Kernel_i, ijsg_i,ijeg_i, ikrs,ikre, ikzs,ikze)
- ! Collision matrix
+ !___________________ 5D ARRAYS __________________________
+ ! Moments with ghost degrees for p+2 p-2, j+1, j-1 truncations
+ CALL allocate_array( moments_e, ipsg_e,ipeg_e, ijsg_e,ijeg_e, ikrs,ikre, ikzs,ikze, 1,ntimelevel )
+ CALL allocate_array( moments_i, ipsg_i,ipeg_i, ijsg_i,ijeg_i, ikrs,ikre, ikzs,ikze, 1,ntimelevel )
+
+ ! Moments right-hand-side (contains linear part of hierarchy)
+ CALL allocate_array( moments_rhs_e, ips_e,ipe_e, ijs_e,ije_e, ikrs,ikre, ikzs,ikze, 1,ntimelevel )
+ CALL allocate_array( moments_rhs_i, ips_i,ipe_i, ijs_i,ije_i, ikrs,ikre, ikzs,ikze, 1,ntimelevel )
+
+ ! Collision term
+ CALL allocate_array( TColl_e, ips_e,ipe_e, ijs_e,ije_e , ikrs,ikre, ikzs,ikze)
+ CALL allocate_array( TColl_i, ips_i,ipe_i, ijs_i,ije_i , ikrs,ikre, ikzs,ikze)
+
+ ! Non linear terms and dnjs table
+ CALL allocate_array( Sepj, ips_e,ipe_e, ijs_e,ije_e, ikrs,ikre, ikzs,ikze )
+ CALL allocate_array( Sipj, ips_i,ipe_i, ijs_i,ije_i, ikrs,ikre, ikzs,ikze )
+ CALL allocate_array( dnjs, 1,maxj+1, 1,maxj+1, 1,maxj+1)
+
+ !___________________ 2x5D ARRAYS __________________________
+ !! Collision matrices
IF (CO .LT. -1) THEN !DK collision matrix (same for every k)
CALL allocate_array( Ceepj, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), 1,1, 1,1)
CALL allocate_array( CeipjT, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), 1,1, 1,1)
@@ -47,13 +67,4 @@ SUBROUTINE memory
CALL allocate_array( CiepjF, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxe+1)*(jmaxe+1), ikrs,ikre, ikzs,ikze)
ENDIF
- ! Collision term
- CALL allocate_array( TColl_e, ips_e,ipe_e, ijs_e,ije_e , ikrs,ikre, ikzs,ikze)
- CALL allocate_array( TColl_i, ips_i,ipe_i, ijs_i,ije_i , ikrs,ikre, ikzs,ikze)
-
- ! Non linear terms and dnjs table
- CALL allocate_array( Sepj, ips_e,ipe_e, ijs_e,ije_e, ikrs,ikre, ikzs,ikze )
- CALL allocate_array( Sipj, ips_i,ipe_i, ijs_i,ije_i, ikrs,ikre, ikzs,ikze )
- CALL allocate_array( dnjs, 1,maxj+1, 1,maxj+1, 1,maxj+1)
-
END SUBROUTINE memory
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index aa483103..0833d3f7 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -6,11 +6,12 @@ MODULE processing
implicit none
REAL(dp), PUBLIC, PROTECTED :: pflux_ri, gflux_ri
-
- PUBLIC :: compute_radial_ion_transport
-
+
+ PUBLIC :: compute_radial_ion_transport, compute_density, compute_temperature
+
CONTAINS
+! 1D diagnostic to compute the average radial particle transport <n_i v_ExB>_r
SUBROUTINE compute_radial_ion_transport
USE fields, ONLY : moments_i, phi
USE array, ONLY : kernel_i
@@ -20,8 +21,8 @@ SUBROUTINE compute_radial_ion_transport
REAL(dp) :: kz_, buffer(1:2)
INTEGER :: i_, world_rank, world_size, root
- pflux_local = 0._dp
- gflux_local = 0._dp
+ pflux_local = 0._dp ! particle flux
+ gflux_local = 0._dp ! gyrocenter flux
IF(ips_i .EQ. 1) THEN
! Loop to compute gamma_kr = sum_kz sum_j -i k_z Kernel_j Ni00 * phi
DO ikz = ikzs,ikze
@@ -58,8 +59,77 @@ SUBROUTINE compute_radial_ion_transport
ENDIF
ENDIF
ENDIF
+END SUBROUTINE compute_radial_ion_transport
+! Compute the 2D particle density for electron and ions (sum over Laguerre)
+SUBROUTINE compute_density
+ USE fields, ONLY : moments_i, moments_e
+ USE array, ONLY : dens_e, dens_i, kernel_e, kernel_i
+ USE time_integration, ONLY : updatetlevel
+ IMPLICIT NONE
-END SUBROUTINE compute_radial_ion_transport
+ IF( (ips_i .EQ. 1) .AND. (ips_e .EQ. 1) ) THEN
+ ! Loop to compute dens_i = sum_j kernel_j Ni0j at each k
+ DO ikz = ikzs,ikze
+ DO ikr = ikrs,ikre
+ ! electron density
+ dens_e(ikr,ikz) = 0._dp
+ DO ij = ijs_e, ije_e
+ dens_e(ikr,ikz) = dens_e(ikr,ikz) + &
+ kernel_e(ij,ikr,ikz) * moments_e(1,ij,ikr,ikz,updatetlevel)
+ ENDDO
+ ! ion density
+ dens_i(ikr,ikz) = 0._dp
+ DO ij = ijs_i, ije_i
+ dens_i(ikr,ikz) = dens_i(ikr,ikz) + &
+ kernel_i(ij,ikr,ikz) * moments_i(1,ij,ikr,ikz,updatetlevel)
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+ CALL manual_2D_bcast(dens_e(ikrs:ikre,ikzs:ikze))
+ CALL manual_2D_bcast(dens_i(ikrs:ikre,ikzs:ikze))
+END SUBROUTINE compute_density
+
+! Compute the 2D particle temperature for electron and ions (sum over Laguerre)
+SUBROUTINE compute_temperature
+ USE fields, ONLY : moments_i, moments_e
+ USE array, ONLY : temp_e, temp_i, kernel_e, kernel_i
+ USE time_integration, ONLY : updatetlevel
+ IMPLICIT NONE
+ REAL(dp) :: j_dp
+ COMPLEX(dp) :: Tperp, Tpar
+
+ IF( ((ips_i .EQ. 1) .AND. (ips_e .EQ. 1)) ) THEN
+ ! Loop to compute T = 1/3*(Tpar + 2Tperp)
+ DO ikz = ikzs,ikze
+ DO ikr = ikrs,ikre
+ ! electron temperature
+ Tpar = 0._dp; Tperp = 0._dp
+ DO ij = ijs_e, ije_e
+ j_dp = REAL(ij-1,dp)
+ Tpar = Tpar + kernel_e(ij,ikr,ikz)* &
+ (SQRT2*moments_e(3,ij,ikr,ikz,updatetlevel) + moments_e(1,ij,ikr,ikz,updatetlevel))
+ Tperp = Tperp + moments_e(1,ij,ikr,ikz,updatetlevel)*&
+ ((2._dp*j_dp+1)*kernel_e(ij,ikr,ikz) - (j_dp+1)*kernel_e(ij+1,ikr,ikz) - j_dp*kernel_e(ij-1,ikr,ikz))
+ ENDDO
+ temp_e(ikr,ikz) = (Tpar + 2._dp*Tperp)/3._dp
+
+ ! ion temperature
+ Tpar = 0._dp; Tperp = 0._dp
+ DO ij = ijs_i, ije_i
+ j_dp = REAL(ij-1,dp)
+ Tpar = Tpar + kernel_i(ij,ikr,ikz)* &
+ (SQRT2*moments_i(3,ij,ikr,ikz,updatetlevel) + moments_i(1,ij,ikr,ikz,updatetlevel))
+ Tperp = Tperp + moments_i(1,ij,ikr,ikz,updatetlevel)*&
+ ((2._dp*j_dp+1)*kernel_i(ij,ikr,ikz) - (j_dp+1)*kernel_i(ij+1,ikr,ikz) - j_dp*kernel_i(ij-1,ikr,ikz))
+ ENDDO
+ temp_i(ikr,ikz) = (Tpar + 2._dp*Tperp)/3._dp
+ ENDDO
+ ENDDO
+ ENDIF
+ CALL manual_2D_bcast(temp_e(ikrs:ikre,ikzs:ikze))
+ CALL manual_2D_bcast(temp_i(ikrs:ikre,ikzs:ikze))
+END SUBROUTINE compute_temperature
-END MODULE processing
\ No newline at end of file
+END MODULE processing
--
GitLab