Skip to content
Snippets Groups Projects
Commit 9ef07a7a authored by Antoine Cyril David Hoffmann's avatar Antoine Cyril David Hoffmann
Browse files

Added density and temperature diagnostics

parent 3fe3d63c
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/diagnose.F90
+ 30
− 63
View file @ 9ef07a7a
......@@ -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)
......
src/diagnostics_par_mod.F90
+ 3
− 1
View file @ 9ef07a7a
......@@ -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)
......
src/memory.F90
+ 30
− 19
View file @ 9ef07a7a
......@@ -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
src/processing_mod.F90
+ 77
− 7
View file @ 9ef07a7a
......@@ -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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment