Homec4science
Diffusion FENNECS: Finite Element Non-Neutral Electron Cloud Simulator (geom_testmumps) Browse

fennecs/src/diagnose.f90
22b60c541a80geom_testmumps

Last Change
Raw File

diagnose.f90
View Options

SUBROUTINE diagnose(kstep)
!
! Diagnostics
!
USE basic
USE futils
USE hashtable
USE beam, ONLY : partslist, epot, ekin, etot, etot0, Nplocs_all, collectparts
#if USE_X == 1
USE xg, ONLY : initw, updt_xg_var
#endif
USE fields, ONLY: phi_spline
USE celldiag
USE mpi
Use geometry
IMPLICIT NONE
!
INTEGER, INTENT(in) :: kstep
!
! Local vars and arrays
INTEGER, PARAMETER :: BUFSIZE = 20
CHARACTER(len=128) :: str, fname, grpname
INTEGER:: Ntotal ! Total number of simulated particles remaining in the simulation
INTEGER:: partsrank, partsdim(2)
INTEGER:: Nplocs_all_save(64)
INTEGER:: i
!________________________________________________________________________________
! 1. Initial diagnostics
IF( kstep .EQ. 0 .and. mpirank .eq. 0) THEN
! Only process 0 should save on file
!
WRITE(*,'(a)') ' Initial diagnostics'
!
! 1.1 Initial run or when NEWRES set to .TRUE.
!
IF( .NOT. nlres .OR. newres) THEN
CALL creatf(resfile, fidres, 'Espic2d result file', real_prec='d')
WRITE(*,'(3x,a,a)') TRIM(resfile), ' created'
!
! Label the run
IF( LEN_TRIM(label1).GT.0 ) CALL attach(fidres, "/", "label1", TRIM(label1))
IF( LEN_TRIM(label2).GT.0 ) CALL attach(fidres, "/", "label2", TRIM(label2))
IF( LEN_TRIM(label3).GT.0 ) CALL attach(fidres, "/", "label3", TRIM(label3))
IF( LEN_TRIM(label4).GT.0 ) CALL attach(fidres, "/", "label4", TRIM(label4))
!
! Job number
jobnum = 0
!
! Data group
CALL creatg(fidres, "/data", "data")
CALL creatg(fidres, "/data/var0d", "0d history arrays")
CALL creatg(fidres, "/data/var1d","1d history arrays")
CALL creatg(fidres, "/data/part", "part phase space")
CALL creatg(fidres, "/data/fields", "Electro static potential and Er Ez fields")
!
! File group
CALL creatg(fidres, "/files", "files")
CALL attach(fidres, "/files", "jobnum", jobnum)
CALL putarr(fidres, "/data/var1d/rgrid", rgrid)
CALL putarr(fidres, "/data/var1d/zgrid", zgrid)
CALL creatd(fidres, 1, (/ 64 /), "/data/var0d/Nplocs_all", "Nplocs_all")
! Part and fields vectors
! Initialize time-dependent particle variables
CALL creatd(fidres, 1, SHAPE(partslist(1)%R), "/data/part/R", "R")
CALL creatd(fidres, 1, SHAPE(partslist(1)%Z), "/data/part/Z", "Z")
CALL creatd(fidres, 1, SHAPE(partslist(1)%Z), "/data/part/THET", "THET")
CALL creatd(fidres, 1, SHAPE(partslist(1)%UZ), "/data/part/UR", "UZ")
CALL creatd(fidres, 1, SHAPE(partslist(1)%UR), "/data/part/UZ", "UR")
CALL creatd(fidres, 1, SHAPE(partslist(1)%UTHET), "/data/part/UTHET", "UTHET")
CALL creatd(fidres, 1, SHAPE(partslist(1)%Rindex), "/data/part/Rindex", "Rindex")
CALL creatd(fidres, 1, SHAPE(partslist(1)%Zindex), "/data/part/Zindex", "Zindex")
CALL creatd(fidres, 1, SHAPE(partslist(1)%partindex), "/data/part/partindex", "partindex")
CALL creatd(fidres, 1, SHAPE(partslist(1)%pot), "/data/part/pot", "pot")
CALL creatd(fidres, 0, SHAPE(time), "/data/part/time", "time" )
CALL creatd(fidres, 0, SHAPE(Ntotal), "/data/part/Nparts", "number of remaining parts in the simulation space")
CALL attach(fidres,'/data/part', "q", partslist(1)%q)
CALL attach(fidres,'/data/part', "m", partslist(1)%m)
CALL attach(fidres,'/data/part', "weight", partslist(1)%weight)
DO i=2,nbspecies
IF ( .not. partslist(i)%is_test) CONTINUE
WRITE(grpname,'(a,i2)')'/data/part/',i
CALL creatg(fidres, grpname, "specific specie phase space")
CALL creatd(fidres, 0, SHAPE(time), trim(grpname) // "/time", "time")
CALL creatd(fidres, 0, SHAPE(time), trim(grpname) //"/Nparts", "number of remaining parts")
CALL creatd(fidres, 1, SHAPE(partslist(i)%R), trim(grpname) // "/R", "radial pos")
CALL creatd(fidres, 1, SHAPE(partslist(i)%R), trim(grpname) // "/Z", "axial pos")
CALL creatd(fidres, 1, SHAPE(partslist(i)%R), trim(grpname) // "/THET", "azimuthal pos")
CALL creatd(fidres, 1, SHAPE(partslist(i)%R), trim(grpname) // "/UZ", "axial beta*gamma")
CALL creatd(fidres, 1, SHAPE(partslist(i)%R), trim(grpname) // "/UR", "radial beta*gamma")
CALL creatd(fidres, 1, SHAPE(partslist(i)%R), trim(grpname) // "/UTHET", "azimuthal beta*gamma")
CALL creatd(fidres, 1, SHAPE(partslist(i)%R), trim(grpname) // "/pot", "electric potential")
CALL creatd(fidres, 1, SHAPE(partslist(i)%R), trim(grpname) // "/Rindex", "radial grid index")
CALL creatd(fidres, 1, SHAPE(partslist(i)%R), trim(grpname) // "/Zindex", "axial grid index")
CALL creatd(fidres, 1, SHAPE(partslist(i)%R), trim(grpname) // "/partindex", "particle index")
CALL attach(fidres,trim(grpname), "q", partslist(i)%q)
CALL attach(fidres,trim(grpname), "m", partslist(i)%m)
CALL attach(fidres,trim(grpname), "weight", partslist(i)%weight)
END DO
CALL creatd(fidres, 1, SHAPE(pot), "/data/fields/pot", "pot")
CALL creatd(fidres, 1, SHAPE(Er), "/data/fields/Er", "Er")
CALL creatd(fidres, 1, SHAPE(Ez), "/data/fields/Ez", "Ez")
CALL creatd(fidres, 1, SHAPE(phi_spline), "/data/fields/phi", "spline form of Phi")
CALL creatd(fidres, 2, SHAPE(moments), "/data/fields/moments", "moments",compress=.true.,chunking=(/1,nrank(1)*nrank(2)/))
!CALL creatd(fidres, 2, SHAPE(moments), "/data/fields/moments", "moments")
CALL creatd(fidres, 0, SHAPE(time), "/data/fields/time", "time" )
CALL putarr(fidres, "/data/fields/Br", Br)
CALL putarr(fidres, "/data/fields/Bz", Bz)
CALL putarr(fidres, "/data/fields/Athet", Athet)
CALL putarr(fidres, "/data/fields/volume", Volume)
!
! 1.2 Restart run
!
ELSE
CALL cp2bk(resfile) ! backup previous result file
CALL openf(resfile, fidres, real_prec='d')
WRITE(*,'(3x,a,a)') TRIM(resfile), ' open'
CALL getatt(fidres, "/files", "jobnum", jobnum)
jobnum = jobnum+1
WRITE(*,'(3x,a,i3)') "Current Job Number =", jobnum
CALL attach(fidres, "/files", "jobnum", jobnum)
END IF
!
! Add input namelist variables as attributes of /data/input
WRITE(str,'(a,i2.2)') "/data/input.",jobnum
CALL creatg(fidres, TRIM(str))
CALL attach(fidres, TRIM(str), "job_time", job_time)
CALL attach(fidres, TRIM(str), "extra_time", extra_time)
CALL attach(fidres, TRIM(str), "dt", dt)
CALL attach(fidres, TRIM(str), "tmax", tmax)
CALL attach(fidres, TRIM(str), "nrun", nrun)
CALL attach(fidres, TRIM(str), "nlres", nlres)
CALL attach(fidres, TRIM(str), "nlsave", nlsave)
CALL attach(fidres, TRIM(str), "newres", newres)
CALL attach(fidres, TRIM(str), "nz", nz)
CALL putarr(fidres, TRIM(str)//"/lz", lz)
CALL attach(fidres, TRIM(str), "nplasma", nplasma)
CALL attach(fidres, TRIM(str), "potinn", potinn)
CALL attach(fidres, TRIM(str), "potout", potout)
CALL attach(fidres, TRIM(str), "B0", B0)
CALL attach(fidres, TRIM(str), "Rcurv", Rcurv)
CALL attach(fidres, TRIM(str), "width", width)
CALL attach(fidres, TRIM(str), "n0", n0)
CALL attach(fidres, TRIM(str), "temp", partslist(1)%temperature)
CALL attach(fidres, TRIM(str), "it0d", it0d)
CALL attach(fidres, TRIM(str), "it2d", it2d)
CALL attach(fidres, TRIM(str), "itparts", itparts)
CALL attach(fidres, TRIM(str), "nlclassical", nlclassical)
CALL attach(fidres, TRIM(str), "nlPhis", nlPhis)
CALL attach(fidres, TRIM(str), "qsim", partslist(1)%q*partslist(1)%weight)
CALL attach(fidres, TRIM(str), "msim", partslist(1)%m*partslist(1)%weight)
CALL attach(fidres, TRIM(str), "startstep", cstep)
CALL attach(fidres, TRIM(str), "H0", partslist(1)%H0)
CALL attach(fidres, TRIM(str), "P0", partslist(1)%P0)
CALL putarr(fidres, TRIM(str)//"/femorder", femorder)
CALL putarr(fidres, TRIM(str)//"/ngauss", ngauss)
CALL putarr(fidres, TRIM(str)//"/nnr", nnr)
CALL putarr(fidres, TRIM(str)//"/radii", radii)
CALL putarr(fidres, TRIM(str)//"/plasmadim", plasmadim)
CALL attach(fidres, TRIM(str), "rawparts", .true.)
CALL attach(fidres, TRIM(str), "nbspecies", nbspecies)
! Save geometry parameters for non conforming boundary conditions
Call diag_geom(fidres,str,rnorm)
! Save STDIN of this run
WRITE(str,'(a,i2.2)') "/files/STDIN.",jobnum
INQUIRE(unit=lu_in, name=fname)
CALL putfile(fidres, TRIM(str), TRIM(fname))
!
! Initialize buffers for 0d history arrays
CALL htable_init(hbuf0, BUFSIZE)
CALL set_htable_fileid(hbuf0, fidres, "/data/var0d")
!
! Initialize Xgrafix
#if USE_X == 1
IF(nlxg) THEN
CALL initw
END IF
#endif
END IF
IF(kstep .EQ. 0) THEN
! Initialize particle cell diagnostic
CALL celldiag_init(lu_in,time, fidres)
CLOSE(lu_in)
END IF
!________________________________________________________________________________
! 2. Periodic diagnostics
IF( kstep .NE. -1) THEN
IF(modulo(step,ittracer) .eq. 0 .or. nlend) THEN
! We gather the traced particles on the mpi host
DO i=1,nbspecies
IF(partslist(i)%is_test) CALL collectparts(partslist(i))
END DO
END IF
! Cell diag quantities
IF(modulo(step,itcelldiag).eq. 0 .or. nlend) THEN
!
CALL celldiag_save(time, fidres)
END IF
!
! Only process 0 should save on file
IF(mpirank .ne. 0) RETURN
!
IF (mpisize .gt. 1) THEN
partslist(1)%Nptot=sum(Nplocs_all)
END IF
!
IF(modulo(step,ittext).eq. 0 .or. nlend) THEN
WRITE(*,'(a,1x,i8.8,a1,i8.8,20x,a,1pe10.3)') '*** Timestep (this run/total) =', &
& step, '/', cstep, 'Time =', time
if( abs(etot).gt. 0) then
WRITE(*,'(a,5(1pe12.4),1x,i8.8,a1,i8.8)') 'Epot, Ekin, Etot, Eerr, Eerr rel, Nbparts/Ntotal', epot, ekin, etot, etot-etot0,(etot-etot0)/etot, partslist(1)%Nptot,'/',nplasma
else
WRITE(*,'(a,4(1pe12.4),1x,i8.8,a1,i8.8)') 'Epot, Ekin, Etot, Eerr, Nbparts/Ntotal', epot, ekin, etot, etot-etot0, partslist(1)%Nptot,'/',nplasma
end if
IF(mpisize .gt. 1 ) WRITE(*,'(a,64i10.7)') 'Nbparts per proc', Nplocs_all
END IF
!________________________________________________________________________________
!
! 2.1 0d history arrays
!
IF(modulo(step,it0d).eq. 0 .or. nlend) THEN
CALL add_record(hbuf0, "time", "simulation time", time)
CALL add_record(hbuf0, "epot", "potential energy", epot)
CALL add_record(hbuf0, "ekin", "kinetic energy", ekin)
CALL add_record(hbuf0, "etot", "total energy", etot)
CALL add_record(hbuf0, "etot0", "theoretical total energy", etot0)
CALL add_record(hbuf0, "nbparts", "number of remaining parts in the simulation space", REAL(partslist(1)%Nptot,kind=db))
CALL htable_endstep(hbuf0)
Nplocs_all_save=0
Nplocs_all_save(1:mpisize)=Nplocs_all(0:mpisize-1)
CALL append(fidres, "/data/var0d/Nplocs_all", REAL(Nplocs_all_save,kind=db))
END IF
!
! 2.2 2d profiles
IF(modulo(step,it2d).eq. 0 .or. nlend) THEN
CALL append(fidres, "/data/fields/time", time)
CALL append(fidres, "/data/fields/pot", pot*phinorm)
CALL append(fidres, "/data/fields/Er", Er*enorm)
CALL append(fidres, "/data/fields/Ez", Ez*enorm)
CALL append(fidres, "/data/fields/phi", phi_spline*phinorm)
CALL append(fidres, "/data/fields/moments", moments)
END IF
!
! 2.3 main specie quantities
IF(modulo(step,itparts).eq. 0 .or. nlend) THEN
!PRINT*, 'write particles to file_____________________'
CALL append(fidres, "/data/part/time", time)
CALL append(fidres, "/data/part/Nparts", REAL(partslist(1)%Nptot,kind=db))
IF ( isdataset(fidres,'/data/part/R') ) THEN
CALL getdims(fidres, '/data/part/R', partsrank, partsdim)
partsdim(1)=min(size(partslist(1)%R,1), partsdim(1))
CALL append(fidres, "/data/part/R", partslist(1)%R(1:partsdim(1))*rnorm)
CALL append(fidres, "/data/part/Z", partslist(1)%Z(1:partsdim(1))*rnorm)
CALL append(fidres, "/data/part/THET", partslist(1)%THET(1:partsdim(1)))
CALL append(fidres, "/data/part/UZ", partslist(1)%UZ(1:partsdim(1))/partslist(1)%gamma(1:partsdim(1)))
CALL append(fidres, "/data/part/UR", partslist(1)%UR(1:partsdim(1))/partslist(1)%gamma(1:partsdim(1)))
CALL append(fidres, "/data/part/UTHET", partslist(1)%UTHET(1:partsdim(1))/partslist(1)%gamma(1:partsdim(1)))
CALL append(fidres, "/data/part/pot", partslist(1)%pot(1:partsdim(1))*phinorm)
CALL append(fidres, "/data/part/Rindex", REAL(partslist(1)%Rindex(1:partsdim(1)),kind=db))
CALL append(fidres, "/data/part/Zindex", REAL(partslist(1)%Zindex(1:partsdim(1)),kind=db))
CALL append(fidres, "/data/part/partindex", REAL(partslist(1)%partindex(1:partsdim(1)),kind=db))
END IF
END IF
!
! 2.4 Tracer quantities
IF(modulo(step,ittracer).eq. 0 .or. nlend) THEN
!PRINT*, 'write particles to file_____________________'
DO i=2,nbspecies
IF ( .not. partslist(i)%is_test) CYCLE
WRITE(grpname,'(a,i2,a)')'/data/part/',i,'/'
CALL append(fidres, trim(grpname) // "time", time)
CALL append(fidres, trim(grpname) //"Nparts", REAL(partslist(i)%Nptot,kind=db))
IF ( isdataset(fidres,trim(grpname)//'R') ) THEN
CALL getdims(fidres, trim(grpname) // 'R', partsrank, partsdim)
partsdim(1)=min(size(partslist(i)%R,1), partsdim(1))
CALL append(fidres, trim(grpname) // "R", partslist(i)%R(1:partsdim(1))*rnorm)
CALL append(fidres, trim(grpname) // "Z", partslist(i)%Z(1:partsdim(1))*rnorm)
CALL append(fidres, trim(grpname) // "THET", partslist(i)%THET(1:partsdim(1)))
CALL append(fidres, trim(grpname) // "UZ", partslist(i)%UZ(1:partsdim(1))/partslist(i)%gamma(1:partsdim(1)))
CALL append(fidres, trim(grpname) // "UR", partslist(i)%UR(1:partsdim(1))/partslist(i)%gamma(1:partsdim(1)))
CALL append(fidres, trim(grpname) // "UTHET", partslist(i)%UTHET(1:partsdim(1))/partslist(i)%gamma(1:partsdim(1)))
CALL append(fidres, trim(grpname) // "pot", partslist(i)%pot(1:partsdim(1))*phinorm)
CALL append(fidres, trim(grpname) // "Rindex", REAL(partslist(i)%Rindex(1:partsdim(1)),kind=db))
CALL append(fidres, trim(grpname) // "Zindex", REAL(partslist(i)%Zindex(1:partsdim(1)),kind=db))
CALL append(fidres, trim(grpname) // "partindex", REAL(partslist(i)%partindex(1:partsdim(1)),kind=db))
END IF
END DO
!
END IF
! 2.5 3d profiles
!
!
! 2.6 Xgrafix
!
#if USE_X == 1
IF(nlxg .AND. modulo(kstep,itgraph) .eq. 0) THEN
call xgevent
CALL updt_xg_var
CALL xgupdate
END IF
#endif
!________________________________________________________________________________
! 3. Final diagnostics
ELSE
! Only process 0 should save on file
IF(mpirank .ne. 0) RETURN
!
! Flush 0d history array buffers
CALL htable_hdf5_flush(hbuf0)
!
! Close all diagnostic files
CALL closef(fidres)
!________________________________________________________________________________
END IF
!
END SUBROUTINE diagnose
c4science · Help