!------------------------------------------------------------------------------ ! EPFL/Swiss Plasma Center !------------------------------------------------------------------------------ ! ! MODULE: celldiag ! !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> Represent a diagnostic positioned at cell indices (rindex,zindex) that saves the individual particles !> position and velocity !------------------------------------------------------------------------------ MODULE celldiag ! USE constants use mpi USE mpihelper USE basic, ONLY: mpirank, mpisize, vnorm, rnorm, Zbounds, zgrid, & & nlclassical, nlmaxwellsource, phinorm, nbcelldiag USE beam USE futils IMPLICIT NONE PRIVATE INTEGER, SAVE, ALLOCATABLE :: specieid(:) !< position of the specie in partslist INTEGER, SAVE, ALLOCATABLE :: rindex(:) !< radial index for the diagnostic position INTEGER, SAVE, ALLOCATABLE :: zindex(:) !< axial index for the diagnostic position TYPE(particles), ALLOCATABLE, SAVE :: diagnosed_parts(:) !< Stores the particles properties at position (rindex,zindex) CHARACTER(len=20), SAVE, ALLOCATABLE :: groupname(:) !< Name of the group in the hdf5 file INTEGER, SAVE , ALLOCATABLE :: h5storelength(:) !< particles capacity of the hdf5 dataset NAMELIST /celldiagparams/ specieid, rindex, zindex PUBLIC:: celldiag_init, celldiag_save contains subroutine celldiag_init(lu_in, diagfile_id) implicit none INTEGER, INTENT(IN) :: lu_in INTEGER, INTENT(IN):: diagfile_id INTEGER:: i ALLOCATE(specieid(nbcelldiag), rindex(nbcelldiag), zindex(nbcelldiag)) ALLOCATE(diagnosed_parts(nbcelldiag), groupname(nbcelldiag), h5storelength(nbcelldiag)) Rewind(lu_in) if(nbcelldiag .gt. 0) then READ(lu_in, celldiagparams) if(mpirank .eq. 0) WRITE(*, celldiagparams) Do i=1,nbcelldiag CALL creat_parts(diagnosed_parts(i), 500) IF(mpirank .eq. 0) THEN WRITE(groupname(i),'(a,i2.2)') "/data/celldiag/",i If(.not. isgroup(diagfile_id, "/data/celldiag/")) THEN CALL creatg(diagfile_id, "/data/celldiag") CALL attach(diagfile_id, "/data/celldiag", "nbcelldiag", nbcelldiag) END IF CALL celldiag_createh5group(diagfile_id, groupname(i), rindex(i), zindex(i), specieid(i), diagnosed_parts(i), h5storelength(i)) END IF END DO END IF End subroutine celldiag_init subroutine celldiag_createh5group(diagfile_id, groupname, rindex, zindex, specid, diag_parts, h5strlength) INTEGER, INTENT(IN):: diagfile_id, rindex, zindex, specid CHARACTER(len=*), INTENT(IN):: groupname TYPE(particles), INTENT(IN):: diag_parts INTEGER, INTENT(INOUT):: h5strlength INTEGER:: partsrank, partsdim(2) If(.not. isgroup(diagfile_id, TRIM(groupname))) CALL creatg(diagfile_id, TRIM(groupname)) CALL attach(diagfile_id, TRIM(groupname), "rindex", rindex) CALL attach(diagfile_id, TRIM(groupname), "zindex", zindex) CALL attach(diagfile_id, trim(groupname), "q", partslist(specid)%q) CALL attach(diagfile_id, trim(groupname), "m", partslist(specid)%m) CALL attach(diagfile_id, trim(groupname), "weight", partslist(specid)%weight) If(.not. isdataset(diagfile_id, trim(groupname) // "/time")) CALL creatd(diagfile_id, 0, SHAPE(rindex), trim(groupname) // "/time", "time") If(.not. isdataset(diagfile_id, trim(groupname) // "/Nparts")) CALL creatd(diagfile_id, 0, SHAPE(rindex), trim(groupname) //"/Nparts", "number of remaining parts") If(.not. isdataset(diagfile_id, trim(groupname) // "/R")) CALL creatd(diagfile_id, 1, SHAPE(diag_parts%pos(1,:)), trim(groupname) // "/R", "radial pos") If(.not. isdataset(diagfile_id, trim(groupname) // "/Z")) CALL creatd(diagfile_id, 1, SHAPE(diag_parts%pos(1,:)), trim(groupname) // "/Z", "axial pos") If(.not. isdataset(diagfile_id, trim(groupname) // "/THET")) CALL creatd(diagfile_id, 1, SHAPE(diag_parts%pos(1,:)), trim(groupname) // "/THET", "azimuthal pos") If(.not. isdataset(diagfile_id, trim(groupname) // "/UZ")) CALL creatd(diagfile_id, 1, SHAPE(diag_parts%pos(1,:)), trim(groupname) // "/UZ", "axial beta*gamma") If(.not. isdataset(diagfile_id, trim(groupname) // "/UR")) CALL creatd(diagfile_id, 1, SHAPE(diag_parts%pos(1,:)), trim(groupname) // "/UR", "radial beta*gamma") If(.not. isdataset(diagfile_id, trim(groupname) // "/UTHET")) CALL creatd(diagfile_id, 1, SHAPE(diag_parts%pos(1,:)), trim(groupname) // "/UTHET", "azimuthal beta*gamma") If(.not. isdataset(diagfile_id, trim(groupname) // "/pot")) CALL creatd(diagfile_id, 1, SHAPE(diag_parts%pot), trim(groupname) // "/pot", "electric potential") CALL getdims(diagfile_id, trim(groupname) // '/R', partsrank, partsdim) h5strlength=partsdim(1) END subroutine celldiag_createh5group subroutine celldiag_save(time, diagfile_id) implicit none REAL(kind=db), INTENT(IN) :: time INTEGER, INTENT(IN) :: diagfile_id INTEGER :: Nbtosave, i ! check if source is on IF(.not. celldiag_on(time)) THEN RETURN END IF Do i=1,nbcelldiag CALL celldiag_save_specie(partslist(specieid(i)),rindex(i),zindex(i),diagnosed_parts(i)) END DO !$OMP BARRIER !$OMP MASTER Do i=1,nbcelldiag if(mpisize .gt. 1) then call collectparts(diagnosed_parts(i)) else diagnosed_parts(i)%Nptot=diagnosed_parts(i)%Nploc end if Nbtosave=min(diagnosed_parts(i)%Nptot,h5storelength(i)) CALL celldiag_write_specie(diagfile_id, diagnosed_parts(i), groupname(i), Nbtosave, time) END DO !$OMP END MASTER !$OMP BARRIER end subroutine celldiag_save SUBROUTINE celldiag_save_specie(p, rindex, zindex, savedp) Type(particles), INTENT(IN) :: p Type(particles), INTENT(INOUT) :: savedp INTEGER, INTENT(IN) :: rindex, zindex INTEGER:: i, destcopyindex savedp%Nploc=0 savedp%collected=.false. IF (p%Nploc .gt. 0 .and. zindex .ge. Zbounds(mpirank) .and. zindex .lt. Zbounds(mpirank+1)) THEN ! Boundary condition at z direction !$OMP DO DO i=1,p%Nploc ! If the particle is in the correct cell, it is saved IF (p%cellindex(3,i) .eq. zindex.and. p%cellindex(1,i) .eq. rindex ) THEN !$OMP CRITICAL (diagparts) savedp%Nploc=savedp%Nploc+1 destcopyindex= savedp%Nploc !$OMP END CRITICAL (diagparts) CALL copy_part(p,i,destcopyindex,savedp) END IF END DO !$OMP END DO NOWAIT END IF END subroutine celldiag_save_specie SUBROUTINE celldiag_write_specie(diagfile_id, savedp, groupname, Nbtosave, time) Type(particles), INTENT(IN) :: savedp INTEGER, INTENT(IN) :: diagfile_id CHARACTER(LEN=*), INTENT(IN) :: groupname INTEGER, INTENT(IN) :: Nbtosave REAL(kind=db), INTENT(IN) :: time IF(mpirank .eq. 0) THEN CALL append(diagfile_id, trim(groupname) // "/time", time) CALL append(diagfile_id, trim(groupname) // "/Nparts", REAL(savedp%Nptot,kind=db)) CALL append(diagfile_id, trim(groupname) // "/R", savedp%pos(1,1:Nbtosave)*rnorm) CALL append(diagfile_id, trim(groupname) // "/Z", savedp%pos(3,1:Nbtosave)*rnorm) CALL append(diagfile_id, trim(groupname) // "/THET", savedp%pos(2,1:Nbtosave)) CALL append(diagfile_id, trim(groupname) // "/UZ", savedp%U(3,1:Nbtosave)/savedp%gamma(1:Nbtosave)) CALL append(diagfile_id, trim(groupname) // "/UR", savedp%U(1,1:Nbtosave)/savedp%gamma(1:Nbtosave)) CALL append(diagfile_id, trim(groupname) // "/UTHET", savedp%U(2,1:Nbtosave)/savedp%gamma(1:Nbtosave)) CALL append(diagfile_id, trim(groupname) // "/pot", savedp%pot(1:Nbtosave)*phinorm) END IF END subroutine celldiag_write_specie logical function celldiag_on(time) REAL(kind=db), intent(in):: time celldiag_on=.true. end function End Module celldiag