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

fennecs/src/basic_mod.f90
e3b71186ea57geom_testpardiso

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=80) :: 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
INTEGER :: nrun=1 !< Number of time steps to run
DOUBLE PRECISION :: job_time=3600.0 !< Time allocated to this job in seconds
DOUBLE PRECISION :: tmax=100000.0 !< Maximum simulation time
DOUBLE PRECISION :: extra_time=60.0 !< Extra time allocated
DOUBLE PRECISION :: dt=1 !< Time step
DOUBLE PRECISION :: 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=1 !< Number of iterations between 2d values writes to hdf5
INTEGER :: itparts=1 !< Number of iterations between particles values writes to hdf5
INTEGER :: ittext=1 !< Number of iterations between text outputs in the console
INTEGER :: itrestart=10000 !< Number of iterations between save of restart.h5 file
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
!
! HDF5 file
CHARACTER(len=64) :: resfile = "results.h5" !< Main result file
CHARACTER(len=64) :: rstfile = "restart.h5" !< Restart 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
INTEGER :: nblock !< Number of slices in Z for stable distribution initialisation
DOUBLE PRECISION :: potinn !< Electric potential at the inner metallic wall
DOUBLE PRECISION :: potout !< Electric potential at the outer metallic wall
DOUBLE PRECISION :: B0 !< Max magnitude of magnetic field
DOUBLE PRECISION, allocatable :: Bz(:), Br(:) !< Magnetic field components
DOUBLE PRECISION, allocatable :: Athet(:) !< Theta component of the magnetic vector potential
TYPE(spline2d) :: splrz !< Spline at r and z
DOUBLE PRECISION, allocatable :: Ez(:), Er(:) !< Electric field components
DOUBLE PRECISION, allocatable :: pot(:) !< Electro static potential
DOUBLE PRECISION :: radii(4) !< Inner and outer radius of cylinder and radii of fine mesh region
DOUBLE PRECISION :: 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
DOUBLE PRECISION :: H0=0 !< Initial value of Hamiltonian for distribution 2
DOUBLE PRECISION :: P0=0 !< Initial canonical angular momentum for distribution 2
DOUBLE PRECISION :: temprescale = -1.0 !< Factor used for temperature rescaling in case of a restart (<0 -> no rescaling)
DOUBLE PRECISION :: lz(2) !< Lower and upper cylinder limits in z direction
DOUBLE PRECISION :: n0 !< Physical plasma density parameter
DOUBLE PRECISION, ALLOCATABLE:: partdensity(:) !< Number of particles per gridcell
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE:: fluidur !< Radial fluid velocity at gridcell
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE:: fluiduz !< Longitudinal fluid velocity at gridcell
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE:: fluiduthet !< Azimuthal fluid velocity at gridcell
DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE:: Presstens !< Pressure tensor per gridcell computed every it2d
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
DOUBLE PRECISION :: dz !< Cell size in z
DOUBLE PRECISION :: dr(3) !< Cell size in r for each region
DOUBLE PRECISION, ALLOCATABLE :: zgrid(:) !< Nodes positions in longitudinal direction
DOUBLE PRECISION, ALLOCATABLE :: rgrid(:) !< Nodes positions in radial direction
DOUBLE PRECISION :: bnorm,enorm,vnorm,tnorm,rnorm,phinorm !< Normalization constants
DOUBLE PRECISION :: qsim !< Charge of superparticles
DOUBLE PRECISION :: msim !< Mass of superparticles
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
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE, SAVE:: rhs !< right hand side of the poisson equation solver
DOUBLE PRECISION :: omegac !< Cyclotronic frequency
DOUBLE PRECISION :: omegap !< Plasma frequency
DOUBLE PRECISION :: V !< Normalised volume
DOUBLE PRECISION :: temp !< Initial temperature of plasma
DOUBLE PRECISION :: Rcurv !< Magnetic field curvature coefficient
DOUBLE PRECISION :: Width !< Distance between two magnetic mirrors
INTEGER, DIMENSION(:), ALLOCATABLE :: Zbounds !< Index of bounds for local processus in Z direction
TYPE(BASICDATA) :: bdata
CONTAINS
!
!================================================================================
SUBROUTINE basic_data
!
! Define basic data
!
use mpihelper
IMPLICIT NONE
!
! Local vars and arrays
CHARACTER(len=256) :: line, inputfilename
!
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
!________________________________________________________________________________
! 1. Process Standard Input File
!
IF(mpirank .eq. 0) THEN
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))
!________________________________________________________________________________
! 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)
bdata=BASICDATA( nlres, nlsave, newres, nlxg, nlppform, nlPhis, nlclassical, &
& nplasma, nz, it0d, it2d, itparts, itgraph, distribtype, nblock, &
& nrun, job_time, extra_time, tmax, dt, potinn, potout, B0, n0, temp, &
& Rcurv, width, H0, P0, femorder, ngauss, nnr, lz, radii, plasmadim, resfile, partperiodic)
END IF
! Distribute run parameters to all MPI workers
CALL init_mpitypes ! initialize all mpi types that will be needed in the simulation
CALL MPI_Bcast(bdata, 1, basicdata_type, 0, MPI_COMM_WORLD, ierr)
CALL readbdata
!
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
DOUBLE PRECISION, OPTIONAL, INTENT(out) :: eltime
!
INTEGER, PARAMETER :: ncmax=128
INTEGER, SAVE :: icall=0, nc=0
DOUBLE PRECISION, SAVE :: startt0=0.0
CHARACTER(len=16), SAVE :: which(ncmax)
INTEGER :: lstr, found, i
DOUBLE PRECISION :: seconds
DOUBLE PRECISION, 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(str)
CHARACTER(len=*), INTENT(in) :: str
INTEGER :: i, ind
loc = 0
DO i=1,nc
ind = INDEX(which(i), str(1:lstr))
IF( ind .GT. 0 .AND. LEN_TRIM(which(i)) .EQ. lstr ) THEN
loc = i
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
END SUBROUTINE readbdata
!================================================================================
END MODULE basic
c4science · Help