diff --git a/src/diagnose.f90 b/src/diagnose.f90 index 9599799..300d36b 100644 --- a/src/diagnose.f90 +++ b/src/diagnose.f90 @@ -1,405 +1,409 @@ SUBROUTINE diagnose(kstep) ! ! Diagnostics ! USE basic USE futils USE hashtable Use maxwsrce Use neutcol 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 Use splinebound Use weighttypes 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") CALL creatd(fidres, 1, (/3/), "/data/var0d/nudcol", "nudcol") ! Part and fields vectors ! Initialize time-dependent particle variables CALL creatd(fidres, 1, SHAPE(partslist(1)%R), "/data/part/R", "R",compress=.true.,chunking=SHAPE(partslist(1)%R)) CALL creatd(fidres, 1, SHAPE(partslist(1)%Z), "/data/part/Z", "Z",compress=.true.,chunking=SHAPE(partslist(1)%R)) CALL creatd(fidres, 1, SHAPE(partslist(1)%Z), "/data/part/THET", "THET",compress=.true.,chunking=SHAPE(partslist(1)%R)) CALL creatd(fidres, 1, SHAPE(partslist(1)%UZ), "/data/part/UR", "UZ",compress=.true.,chunking=SHAPE(partslist(1)%R)) CALL creatd(fidres, 1, SHAPE(partslist(1)%UR), "/data/part/UZ", "UR",compress=.true.,chunking=SHAPE(partslist(1)%R)) CALL creatd(fidres, 1, SHAPE(partslist(1)%UTHET), "/data/part/UTHET", "UTHET",compress=.true.,chunking=SHAPE(partslist(1)%R)) CALL creatd(fidres, 1, SHAPE(partslist(1)%Rindex), "/data/part/Rindex", "Rindex",compress=.true.,chunking=SHAPE(partslist(1)%R)) CALL creatd(fidres, 1, SHAPE(partslist(1)%Zindex), "/data/part/Zindex", "Zindex",compress=.true.,chunking=SHAPE(partslist(1)%R)) CALL creatd(fidres, 1, SHAPE(partslist(1)%partindex), "/data/part/partindex", "partindex",compress=.true.,chunking=SHAPE(partslist(1)%R)) CALL creatd(fidres, 1, SHAPE(partslist(1)%pot), "/data/part/pot", "pot",compress=.true.,chunking=SHAPE(partslist(1)%R)) 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 creatd(fidres, 1, (/7/), "/data/part/nbchange", "number of added parts, lost parts zm,zp,rm,rp, and collisions per type io, ela") 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) 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) CALL putarr(fidres, TRIM(str)//"/potxt", potxt) CALL putarr(fidres, TRIM(str)//"/Erxt", Erxt) CALL putarr(fidres, TRIM(str)//"/Ezxt", Ezxt) ! Save geometry parameters for non conforming boundary conditions Call geom_diag(fidres,str,rnorm) ! Save geometry parameters for non conforming boundary conditions using b-spline curves call splinebound_diag(fidres, str, the_domain) ! Save Maxwellsource parameters for the ad-hoc source Call maxwsrce_diag(fidres,str,vnorm) ! Save neutcol parameters for the electron collisions with neutrals Call neutcol_diag(fidres,str,vnorm) if(.not. isdataset(fidres,'/data/var0d/nudcol'))then CALL creatd(fidres, 1, (/3/), "/data/var0d/nudcol", "nudcol") end if ! 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)) ! Prepare hdf5 file for storing test particles DO i=2,nbspecies IF ( partslist(i)%is_test) THEN WRITE(grpname,'(a,i2)')'/data/part/',i call create_parts_group(partslist(i),trim(grpname),time) END IF END DO CALL attach(fidres, "/data/part", "nbspecies", nbspecies) ! ! 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, 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 IF(modulo(step,itrestart) .eq. 0 .or. modulo(step,itparts) .eq. 0 .or. nlend) THEN ! We gather the traced particles on the mpi host DO i=1,nbspecies CALL collectparts(partslist(i)) END DO END IF IF(modulo(step,ittext) .eq. 0 .or. nlend) THEN ! We gather the number of gained and lost particles on the mpi host do i=1,nbspecies IF(mpirank .eq.0 ) THEN CALL MPI_REDUCE(MPI_IN_PLACE, partslist(i)%nbadded, 1, MPI_INTEGER, MPI_SUM, & & 0, MPI_COMM_WORLD, ierr) CALL MPI_REDUCE(MPI_IN_PLACE, partslist(i)%nblost, 4, MPI_INTEGER, MPI_SUM, & & 0, MPI_COMM_WORLD, ierr) CALL MPI_REDUCE(MPI_IN_PLACE, partslist(i)%nbcolls, 2, MPI_INTEGER, MPI_SUM, & & 0, MPI_COMM_WORLD, ierr) + CALL MPI_REDUCE(partslist(i)%nploc, partslist(i)%nptot, 1, MPI_INTEGER, MPI_SUM, & + & 0, MPI_COMM_WORLD, ierr) ELSE CALL MPI_REDUCE(partslist(i)%nbadded, partslist(i)%nbadded, 1, MPI_INTEGER, MPI_SUM, & & 0, MPI_COMM_WORLD, ierr) CALL MPI_REDUCE(partslist(i)%nblost, partslist(i)%nblost, 4, MPI_INTEGER, MPI_SUM, & & 0, MPI_COMM_WORLD, ierr) CALL MPI_REDUCE(partslist(i)%nbcolls, partslist(i)%nbcolls, 2, MPI_INTEGER, MPI_SUM, & & 0, MPI_COMM_WORLD, ierr) + CALL MPI_REDUCE(partslist(i)%nploc, partslist(i)%nptot, 1, MPI_INTEGER, MPI_SUM, & + & 0, MPI_COMM_WORLD, ierr) partslist(i)%nbadded=0 partslist(i)%nblost=0 partslist(i)%nbcolls=0 END IF end do 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,6(1pe12.4),1x,i8.8,a1,i8.8)') 'Epot, Ekin, Etot, Etot0, Eerr, Eerr rel, Nbparts/Ntotal', epot, ekin, etot, etot0, 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 ) then WRITE(*,'(a,64i10.7)') 'Nbparts per proc', Nplocs_all end if Write(*,'(a)')"speci, added, zmloss, zploss, rmloss, rploss, iocoll, elacoll, tot var, tot" do i=1,nbspecies WRITE(*,'(i04,9i9.7)') i, partslist(i)%nbadded,-partslist(i)%nblost,partslist(i)%nbcolls,partslist(i)%nbadded-sum(partslist(i)%nblost), partslist(i)%Nptot partslist(i)%nbadded=0 partslist(i)%nblost=0 partslist(i)%nbcolls=0 end do 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)) CALL append(fidres, "/data/var0d/nudcol", partslist(1)%nudcol/dt) 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)) !CALL append(fidres, "/data/part/nbchange", REAL((/partslist(1)%nbadded,partslist(1)%nblost,partslist(1)%nbcolls/),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", 0.5*(partslist(1)%UZ(1:partsdim(1))/partslist(1)%gamma(1:partsdim(1))+partslist(1)%UZold(1:partsdim(1))/partslist(1)%gammaold(1:partsdim(1)))) CALL append(fidres, "/data/part/UR", 0.5*(partslist(1)%UR(1:partsdim(1))/partslist(1)%gamma(1:partsdim(1))+partslist(1)%URold(1:partsdim(1))/partslist(1)%gammaold(1:partsdim(1)))) CALL append(fidres, "/data/part/UTHET", 0.5*(partslist(1)%UTHET(1:partsdim(1))/partslist(1)%gamma(1:partsdim(1))+partslist(1)%UTHETold(1:partsdim(1))/partslist(1)%gammaold(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)) !CALL append(fidres, trim(grpname) //"nbchange", REAL((/partslist(i)%nbadded,partslist(i)%nblost,partslist(i)%nbcolls/),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", 0.5*(partslist(i)%UZ(1:partsdim(1))/partslist(i)%gamma(1:partsdim(1)) + partslist(i)%UZold(1:partsdim(1))/partslist(i)%gammaold(1:partsdim(1)))) CALL append(fidres, trim(grpname) // "UR", 0.5*(partslist(i)%UR(1:partsdim(1))/partslist(i)%gamma(1:partsdim(1)) + partslist(i)%URold(1:partsdim(1))/partslist(i)%gammaold(1:partsdim(1)))) CALL append(fidres, trim(grpname) // "UTHET", 0.5*(partslist(i)%UTHET(1:partsdim(1))/partslist(i)%gamma(1:partsdim(1)) + partslist(i)%UTHETold(1:partsdim(1))/partslist(i)%gammaold(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 ! CONTAINS SUBROUTINE create_parts_group(p,grpname, time) USE beam,ONLY: particles type(particles):: p real(kind=db):: time character(len=*):: grpname If(isgroup(fidres, trim(grpname))) return 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(p%R), trim(grpname) // "/R", "radial pos") CALL creatd(fidres, 1, SHAPE(p%R), trim(grpname) // "/Z", "axial pos") CALL creatd(fidres, 1, SHAPE(p%R), trim(grpname) // "/THET", "azimuthal pos") CALL creatd(fidres, 1, SHAPE(p%R), trim(grpname) // "/UZ", "axial beta*gamma") CALL creatd(fidres, 1, SHAPE(p%R), trim(grpname) // "/UR", "radial beta*gamma") CALL creatd(fidres, 1, SHAPE(p%R), trim(grpname) // "/UTHET", "azimuthal beta*gamma") CALL creatd(fidres, 1, SHAPE(p%R), trim(grpname) // "/pot", "electric potential") CALL creatd(fidres, 1, SHAPE(p%R), trim(grpname) // "/Rindex", "radial grid index") CALL creatd(fidres, 1, SHAPE(p%R), trim(grpname) // "/Zindex", "axial grid index") CALL creatd(fidres, 1, SHAPE(p%R), trim(grpname) // "/partindex", "particle index") CALL creatd(fidres, 1, (/7/), trim(grpname) // "nbchange", "number of added parts, lost parts zm,zp,rm,rp, and collisions per type io, ela") CALL attach(fidres,trim(grpname), "q", p%q) CALL attach(fidres,trim(grpname), "m", p%m) CALL attach(fidres,trim(grpname), "weight", p%weight) END SUBROUTINE create_parts_group END SUBROUTINE diagnose