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

fennecs/src/basic_mod.f90
22b60c541a80geom_testmumps

Last Change
Raw File

basic_mod.f90
View Options

MODULE basic
!
USE hashtable
USE constants
USE bsplines
USE mumps_bsplines
USE futils
USE mpihelper
IMPLICIT NONE
!
! Basic module for time dependent problems
!
CHARACTER(len=128) :: label1, label2, label3, label4
!
! BASIC Namelist
!
LOGICAL :: nlres = .FALSE. !< Restart flag
LOGICAL :: nlsave = .TRUE. !< Checkpoint (save) flag
LOGICAL :: newres=.FALSE. !< New result HDF5 file
LOGICAL :: nlxg=.FALSE. !< Show graphical interface Xgrafix
LOGICAL :: nlmaxwellsource = .FALSE. !< Activate the maxwell source
INTEGER :: nrun=1 !< Number of time steps to run
REAL(kind=db) :: job_time=3600.0 !< Time allocated to this job in seconds
REAL(kind=db) :: tmax=100000.0 !< Maximum simulation time
REAL(kind=db) :: extra_time=60.0 !< Extra time allocated
REAL(kind=db) :: dt=1 !< Time step
REAL(kind=db) :: time=0 !< Current simulation time (Init from restart file)
!
! Other basic global vars and arrays
!
INTEGER :: jobnum !< Job number
INTEGER :: step !< Calculation step of this run
INTEGER :: cstep=0 !< Current step number (Init from restart file)
LOGICAL :: nlend !< Signal end of run
INTEGER :: ierr !< Integer used for MPI
INTEGER :: it0d=1 !< Number of iterations between 0d values writes to hdf5
INTEGER :: it2d=100 !< Number of iterations between 2d values writes to hdf5
INTEGER :: itparts=1000 !< Number of iterations between particles values writes to hdf5
INTEGER :: ittext=10 !< Number of iterations between text outputs in the console
INTEGER :: itrestart=10000 !< Number of iterations between save of restart.h5 file
INTEGER :: ittracer=100 !< Number of iterations between save of traced particles position and velocity
INTEGER :: itcelldiag=100000 !< Number of iterations between save of celldiag diagnostic
INTEGER :: nbcelldiag=0 !< Number of celldiagnostics
INTEGER :: itgraph !< Number of iterations between graphical interface updates
INTEGER :: mpirank !< MPIrank of the current processus
INTEGER :: mpisize !< Size of the MPI_COMM_WORLD communicator
INTEGER :: rightproc !< Rank of next processor in the z decomposition
INTEGER :: leftproc !< Rank of previous processor in the z decomposition
!
! List of logical file units
INTEGER :: lu_in = 90 !< File duplicated from STDIN
INTEGER :: lu_stop = 91 !< stop file, see subroutine TESEND
INTEGER :: lu_partfile = 120 !< particle loading file, see beam::loadpartfile
!
! HDF5 file
CHARACTER(len=256) :: resfile = "results.h5" !< Main result file
CHARACTER(len=256) :: rstfile = "restart.h5" !< Restart file
CHARACTER(len=256) :: magnetfile = "" !< H5 file containing the magnetic field definition
CHARACTER(len=256) :: partfile(10)="" !< Particle loading file
INTEGER :: fidres !< File ID for resfile
INTEGER :: fidrst !< File ID for restart file
TYPE(BUFFER_TYPE) :: hbuf0 !< Hashtable for 0d var
!
! Plasma parameters
LOGICAL :: nlPhis= .TRUE. !< Calculate self consistent electric field flag
LOGICAL :: nlclassical= .FALSE. !< If true, solves the equation of motion according to classical
!! dynamics
LOGICAL :: nlperiod(2)=(/.false.,.false./)!< Set periodic splines on or off
LOGICAL :: partperiodic= .TRUE. !< Sets if the particles boundary conditions are periodic or open
INTEGER :: nplasma !< Number of macro-particles on initialisation
INTEGER :: nbspecies = 1 !< Number of particles species also counting tracing particles
INTEGER :: npartsalloc = 0 !< Size of particle memory allocated at the begining of the simulation
INTEGER :: nblock !< Number of slices in Z for stable distribution initialisation
REAL(kind=db) :: potinn !< Electric potential at the inner metallic wall
REAL(kind=db) :: potout !< Electric potential at the outer metallic wall
REAL(kind=db) :: B0 !< Max magnitude of magnetic field
REAL(kind=db), allocatable :: Bz(:), Br(:) !< Magnetic field components
REAL(kind=db), allocatable :: Athet(:) !< Theta component of the magnetic vector potential
TYPE(spline2d), SAVE :: splrz !< Spline at r and z for total electric field
TYPE(spline2d), SAVE :: splrz_ext !< Spline at r and z for external electric field
REAL(kind=db), allocatable :: Ez(:), Er(:) !< Electric field components
REAL(kind=db), allocatable :: pot(:) !< Electro static potential
REAL(kind=db) :: radii(4) !< Inner and outer radius of cylinder and radii of fine mesh region
REAL(kind=db) :: plasmadim(4) !< Zmin Zmax Rmin Rmax values for plasma particle loading
INTEGER :: distribtype=1 !< Type of distribution function used to load the particles
!!1: gaussian, 2: Stable as defined in 4.95 of Davidson
REAL(kind=db) :: H0=0 !< Initial value of Hamiltonian for distribution 2
REAL(kind=db) :: P0=0 !< Initial canonical angular momentum for distribution 2
REAL(kind=db) :: temprescale = -1.0 !< Factor used for temperature rescaling in case of a restart (<0 -> no rescaling)
INTEGER :: samplefactor =-1 !< Factor used for the up-sampling of the particles number
REAL(kind=db) :: lz(2) !< Lower and upper cylinder limits in z direction
REAL(kind=db) :: n0 !< Physical plasma density parameter
REAL(kind=db), DIMENSION(:,:), ALLOCATABLE, SAVE:: moments !< Moments of the distribution function evaluated every it2d
REAL(kind=db), DIMENSION(:), ALLOCATABLE, SAVE:: rhs !< right hand side of the poisson equation solver
REAL(kind=db), DIMENSION(:), ALLOCATABLE, SAVE:: volume !< Volume covered by each spline for density calculation
INTEGER :: nz !< Number of grid intervals in z
INTEGER :: nr !< Total number of grid intervals in r
INTEGER :: nnr(3) !< Number of grid intervals in r in each subdomain
REAL(kind=db) :: dz !< Cell size in z
REAL(kind=db) :: dr(3) !< Cell size in r for each region
REAL(kind=db), ALLOCATABLE :: zgrid(:) !< Nodes positions in longitudinal direction
REAL(kind=db), ALLOCATABLE :: rgrid(:) !< Nodes positions in radial direction
REAL(kind=db) :: bnorm,enorm,vnorm,tnorm,rnorm,phinorm !< Normalization constants
REAL(kind=db) :: qsim !< Charge of superparticles
REAL(kind=db) :: msim !< Mass of superparticles
REAL(kind=db) :: partmass=me !< Mass of physical particle
INTEGER :: femorder(2) !< FEM order
INTEGER :: ngauss(2) !< Number of gauss points
LOGICAL :: nlppform =.TRUE. !< Argument of set_spline
INTEGER, SAVE :: nrank(2) !< Number of splines in both directions
REAL(kind=db) :: omegac !< Cyclotronic frequency
REAL(kind=db) :: omegap !< Plasma frequency
REAL(kind=db) :: temp !< Initial temperature of plasma
REAL(kind=db) :: Rcurv !< Magnetic field curvature coefficient
REAL(kind=db) :: Width !< Distance between two magnetic mirrors
INTEGER, DIMENSION(:), ALLOCATABLE :: Zbounds !< Index of bounds for local processus in Z direction
TYPE(BASICDATA) :: bdata !< Structure used for the mpi communication of the run parameters
CONTAINS
!
!================================================================================
SUBROUTINE basic_data
!
! Define basic data
!
use mpihelper
IMPLICIT NONE
!
! Local vars and arrays
CHARACTER(len=256) :: inputfilename
INTEGER, EXTERNAL :: OMP_GET_MAX_THREADS
!
NAMELIST /BASIC/ job_time, extra_time, nrun, tmax, dt, nlres, nlsave, newres, nlxg, &
& nplasma, potinn, potout, B0, lz, n0, nz, nnr, femorder, ngauss, &
& nlppform, plasmadim, radii, temp, Rcurv, width, it0d, it2d, itparts, ittext, &
& resfile, rstfile, itgraph, nlPhis, distribtype, nblock, nlclassical, H0, P0, partperiodic, &
& temprescale, samplefactor, nlmaxwellsource, npartsalloc, partfile, partmass, nbspecies, &
& ittracer, itcelldiag, nbcelldiag, magnetfile
!________________________________________________________________________________
! 1. Process Standard Input File
!
IF(COMMAND_ARGUMENT_COUNT().NE.1)THEN
WRITE(*,*)'ERROR, ONE COMMAND-LINE ARGUMENT REQUIRED, STOPPING'
STOP
ENDIF
CALL GET_COMMAND_ARGUMENT(1,inputfilename)
OPEN(UNIT=lu_in,FILE=trim(inputfilename),ACTION='READ')
IF(mpirank .eq. 0) THEN
!________________________________________________________________________________
! 1. Label the run
!
READ(lu_in,'(a)') label1
READ(lu_in,'(a)') label2
READ(lu_in,'(a)') label3
READ(lu_in,'(a)') label4
!
WRITE(*,'(12x,a/)') label1(1:len_trim(label1))
WRITE(*,'(12x,a/)') label2(1:len_trim(label2))
WRITE(*,'(12x,a/)') label3(1:len_trim(label3))
WRITE(*,'(12x,a/)') label4(1:len_trim(label4))
!________________________________________________________________________________
! 2. Read in basic data specific to run
!
READ(lu_in,basic)
WRITE(*,basic)
#if _DEBUG==1
WRITE(*,*) "Compiled in debug mode"
#endif
ELSE
READ(lu_in,basic)
END IF
! Distribute run parameters to all MPI workers
CALL init_mpitypes ! initialize all mpi types that will be needed in the simulation
WRITE(*,'(a,i4.2,a,i4.2,a)')"Running on ",mpisize," tasks with", omp_get_max_threads() ," openMP threads"
IF(samplefactor .gt. 1 .and. .not. newres) THEN
IF(mpirank.eq.0) WRITE(*,*)"To increase the number of particles, you need to create a new result file (set newres to 1)"
CALL MPI_abort(MPI_COMM_WORLD,-1,ierr)
END IF
IF (npartsalloc .lt. nplasma) THEN
npartsalloc=nplasma
END IF
!
END SUBROUTINE basic_data
!================================================================================
SUBROUTINE daytim(str)
!
! Print date and time
!
IMPLICIT NONE
!
CHARACTER(len=*), INTENT(in) :: str
!
! Local vars and arrays
CHARACTER(len=16) :: d, t, dat, functime
!________________________________________________________________________________
!
CALL DATE_AND_TIME(d,t)
dat=d(7:8) // '/' // d(5:6) // '/' // d(1:4)
functime=t(1:2) // ':' // t(3:4) // ':' // t(5:10)
WRITE(*,'(a,1x,a,1x,a)') str, dat(1:10), functime(1:12)
!
END SUBROUTINE daytim
!================================================================================
SUBROUTINE timera(cntrl, str, eltime)
!
! Timers (cntrl=0/1 to Init/Update)
!
IMPLICIT NONE
INTEGER, INTENT(in) :: cntrl
CHARACTER(len=*), INTENT(in) :: str
REAL(kind=db), OPTIONAL, INTENT(out) :: eltime
!
INTEGER, PARAMETER :: ncmax=128
INTEGER, SAVE :: icall=0, nc=0
REAL(kind=db), SAVE :: startt0=0.0
CHARACTER(len=16), SAVE :: which(ncmax)
INTEGER :: lstr, found, i
REAL(kind=db) :: seconds
REAL(kind=db), DIMENSION(ncmax), SAVE :: startt = 0.0, endt = 0.0
!________________________________________________________________________________
IF( icall .EQ. 0 ) THEN
icall = icall+1
startt0 = seconds()
END IF
lstr = LEN_TRIM(str)
IF( lstr .GT. 0 ) found = loc(str)
!________________________________________________________________________________
!
SELECT CASE (cntrl)
!
CASE(-1) ! Current wall time
IF( PRESENT(eltime) ) THEN
eltime = seconds() - startt0
ELSE
WRITE(*,'(/a,a,1pe10.3/)') "++ ", ' Wall time used so far = ', seconds() - startt0
END IF
!
CASE(0) ! Init Timer
IF( found .EQ. 0 ) THEN ! Called for the 1st time for 'str'
nc = nc+1
which(nc) = str(1:lstr)
found = nc
END IF
startt(found) = seconds()
!
CASE(1) ! Update timer
endt(found) = seconds() - startt(found)
IF( PRESENT(eltime) ) THEN
eltime = endt(found)
ELSE
WRITE(*,'(/a,a,1pe10.3/)') "++ "//str, ' wall clock time = ', endt(found)
END IF
!
CASE(2) ! Update and reset timer
endt(found) = endt(found) + seconds() - startt(found)
startt(found) = seconds()
IF( PRESENT(eltime) ) THEN
eltime = endt(found)
END IF
!
CASE(9) ! Display all timers
IF( nc .GT. 0 ) THEN
WRITE(*,'(a)') "Timer Summary"
WRITE(*,'(a)') "============="
DO i=1,nc
WRITE(*,'(a20,2x,2(1pe12.3))') TRIM(which(i))//":", endt(i)
END DO
END IF
!
END SELECT
!
CONTAINS
INTEGER FUNCTION loc(funcstr)
CHARACTER(len=*), INTENT(in) :: funcstr
INTEGER :: j, ind
loc = 0
DO j=1,nc
ind = INDEX(which(j), funcstr(1:lstr))
IF( ind .GT. 0 .AND. LEN_TRIM(which(j)) .EQ. lstr ) THEN
loc = j
EXIT
END IF
END DO
END FUNCTION loc
END SUBROUTINE timera
!================================================================================
SUBROUTINE readbdata
nlres = bdata%nlres
nlsave = bdata%nlsave
newres = bdata%newres
nlxg = bdata%nlxg
nlppform = bdata%nlppform
nlPhis = bdata%nlPhis
nlclassical = bdata%nlclassical
nplasma = bdata%nplasma
nz = bdata%nz
it0d = bdata%it0d
it2d = bdata%it2d
itparts = bdata%itparts
itgraph = bdata%itgraph
distribtype = bdata%distribtype
nblock = bdata%nblock
nrun = bdata%nrun
job_time = bdata%job_time
extra_time = bdata%extra_time
tmax = bdata%tmax
dt = bdata%dt
potinn = bdata%potinn
potout = bdata%potout
B0 = bdata%B0
n0 = bdata%n0
temp = bdata%temp
Rcurv = bdata%Rcurv
width = bdata%width
H0 = bdata%H0
P0 = bdata%P0
femorder = bdata%femorder
ngauss = bdata%ngauss
nnr = bdata%nnr
lz = bdata%lz
radii = bdata%radii
plasmadim = bdata%plasmadim
resfile = bdata%resfile
partperiodic = bdata%partperiodic
samplefactor = bdata%samplefactor
END SUBROUTINE readbdata
!================================================================================
END MODULE basic
c4science · Help