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

fennecs/src/celldiag_mod.f90
7a81ec0f73eegeom_testmumps

Last Change
Raw File

celldiag_mod.f90
View Options

!------------------------------------------------------------------------------
! 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
USE beam
IMPLICIT NONE
PRIVATE
INTEGER, SAVE :: specieid !< position of the specie in partslist
INTEGER, SAVE :: rindex !< radial index for the diagnostic position
INTEGER, SAVE :: zindex !< axial index for the diagnostic position
TYPE(particles), SAVE :: diagnosed_parts !< Stores the particles properties at position (rindex,zindex)
CHARACTER(len=20), SAVE :: groupname !< Name of the group in the hdf5 file
INTEGER, SAVE :: h5storelength !< particles capacity of the hdf5 dataset
NAMELIST /celldiagparams/ speciedid, rindex, zindex
PUBLIC:: celldiag_init, celldiag_save
contains
subroutine celldiag_init(lu_in, time, diagfile_id)
implicit none
INTEGER, INTENT(IN) :: lu_in
REAL(kind=db), INTENT(IN):: time
INTEGER, INTENT(IN):: diagfile_id
INTEGER:: ierr, cellnb=1
INTEGER:: partsrank, partsdim(2)
READ(lu_in, celldiagparams)
IF(mpirank .eq. 0) THEN
WRITE(*, celldiagparams)
END IF
CALL creat_parts(diagnosed_parts, 500)
WRITE(groupname,'(a,i2.2)') "/data/celldiag/",cellnb
If(.not. isgroup(diagfile_id, "/data/celldiag/")) CALL creatg(diagfile_id, "/data/celldiag")
If(.not. isgroup(diagfile_id, 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(specieid)%q)
CALL attach(diagfile_id,trim(groupname), "m", partslist(specieid)%m)
CALL attach(diagfile_id,trim(groupname), "weight", partslist(specieid)%weight)
If(.not. isdataset(diagfile_id, trim(groupname) // "/time")) CALL creatd(diagfile_id, 0, 1, trim(groupname) // "/time", "time")
If(.not. isdataset(diagfile_id, trim(groupname) // "/Nparts")) CALL creatd(diagfile_id, 0, 1, trim(groupname) //"/Nparts", "number of remaining parts")
If(.not. isdataset(diagfile_id, trim(groupname) // "/R")) CALL creatd(diagfile_id, 1, SHAPE(diagnosed_parts%R), trim(groupname) // "/R", "radial pos")
If(.not. isdataset(diagfile_id, trim(groupname) // "/Z")) CALL creatd(diagfile_id, 1, SHAPE(diagnosed_parts%R), trim(groupname) // "/Z", "axial pos")
If(.not. isdataset(diagfile_id, trim(groupname) // "/THET")) CALL creatd(diagfile_id, 1, SHAPE(diagnosed_parts%R), trim(groupname) // "/THET", "azimuthal pos")
If(.not. isdataset(diagfile_id, trim(groupname) // "/UZ")) CALL creatd(diagfile_id, 1, SHAPE(diagnosed_parts%R), trim(groupname) // "/UZ", "axial beta*gamma")
If(.not. isdataset(diagfile_id, trim(groupname) // "/UR")) CALL creatd(diagfile_id, 1, SHAPE(diagnosed_parts%R), trim(groupname) // "/UR", "radial beta*gamma")
If(.not. isdataset(diagfile_id, trim(groupname) // "/UTHET")) CALL creatd(diagfile_id, 1, SHAPE(diagnosed_parts%R), trim(groupname) // "/UTHET", "azimuthal beta*gamma")
If(.not. isdataset(diagfile_id, trim(groupname) // "/pot")) CALL creatd(diagfile_id, 1, SHAPE(diagnosed_parts%R), trim(groupname) // "/pot", "electric potential")
CALL getdims(fidres, trim(groupname) // 'R', partsrank, partsdim)
h5storelength=partsdim(1)
End subroutine celldiag_init
subroutine celldiag_save(time)
implicit none
REAL(kind=db), INTENT(IN) :: time
Type(particle), ALLOCATABLE, Dimension(:):: newparts
INTEGER:: npartsadd
INTEGER:: i
REAL(kind=db), ALLOCATABLE:: VR(:),VZ(:),VTHET(:),R(:),Z(:)
! check if source is on
IF(.not. celldiag_on(time)) THEN
RETURN
END IF
CALL celldiag_save_specie(partslist(specieid,rindex,zindex,diagnose_parts))
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:: partsdim
savedp%Nploc=0
savedp%collected=.false.
IF (p%Nploc .gt. 0) THEN
! Boundary condition at z direction
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i)
DO i=1,p%Nploc
! If the particle in the correct cell, it is saved
IF (p%Z(i) .eq. zindex .and. p%R(i) .eq. rindex) THEN
!$OMP CRITICAL (diagparts)
savedp%Nploc=savedp%Nploc+1
CALL copy_part(p,i,savedp%Nploc,savedp)
!$OMP END CRITICAL (diagparts)
END IF
END DO
!$OMP END PARALLEL DO
END IF
if(mpisize .gt. 1) call collectparts(savedp)
partsdim=min(savedp%Nptot,h5storelength)
IF(mpirank .eq. 0) THEN
CALL append(fidres, trim(groupname) // "/time", time)
CALL append(fidres, trim(groupname) // "/Nparts", REAL(savedp%Nptot,kind=db))
CALL append(fidres, trim(groupname) // "/R", savedp%R(1:partsdim)*rnorm)
CALL append(fidres, trim(groupname) // "/Z", savedp%Z(1:partsdim)*rnorm)
CALL append(fidres, trim(groupname) // "/THET", savedp%THET(1:partsdim))
CALL append(fidres, trim(groupname) // "/UZ", savedp%UZ(1:partsdim)/savedp%gamma(1:partsdim))
CALL append(fidres, trim(groupname) // "/UR", savedp%UR(1:partsdim)/savedp%gamma(1:partsdim))
CALL append(fidres, trim(groupname) // "/UTHET", savedp%UTHET(1:partsdim)/savedp%gamma(1:partsdim))
CALL append(fidres, trim(groupname) // "/pot", savedp%pot(1:partsdim)*phinorm)
END IF
END subroutine celldiag_save_specie
logical function celldiag_on(time)
REAL(kind=db), intent(in):: time
celldiag_on=.true.
end function
End Module celldiag
c4science · Help