Homec4science
Diffusion FENNECS lite (master) Browse

fennecslite/src/basic_mod.f90
6ab2c92d5c8amaster

Last Change
Raw File

basic_mod.f90
View Options

MODULE basic
!
! Basic module for time dependent problems
!
USE hashtable
USE constants
USE bsplines
USE mumps_bsplines
USE futils
USE mpihelper
use random_mod
IMPLICIT NONE
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=0 !< Calculation step of this run
INTEGER :: cstep=0 !< Current step number (Init from restart file)
LOGICAL :: nlend=.false. !< 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 :: 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 where r,z are in m and Br, Bz are in T
CHARACTER(len=256) :: partfile(10)="" !< Particle loading file
CHARACTER(len=256) :: addedtestspecfile(10)="" !< Particle file list for added particles at restart
INTEGER :: fidres !< File ID for resfile
INTEGER :: fidrst !< File ID for restart file
TYPE(BUFFER_TYPE) :: hbuf0 !< Hashtable for 0d var
!
! Plasma parameters
LOGICAL :: nlclassical= .FALSE. !< If true, solves the equation of motion according to classical
!! dynamics
LOGICAL :: partperiodic= .TRUE. !< Sets if the particles boundary conditions are periodic or open
INTEGER :: nbaddtestspecies=0 !< On restart number of files to read to add test particles
INTEGER :: nplasma !< Number of macro-particles on initialisation
REAL(kind=db) :: qsim !< Charge of superparticles [C]
REAL(kind=db) :: msim !< Mass of superparticles [kg]
REAL(kind=db) :: partmass=me !< Mass of physical particle [kg]
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) :: n0 !< Physical plasma density parameter [m-3] used in distribtype=1 and for time scales normlisation
REAL(kind=db) :: plasmadim(4) !< Zmin Zmax Rmin Rmax values for plasma particle loading [m]
REAL(kind=db) :: H0=0 !< Initial value of Hamiltonian for distribution 2 [J]
REAL(kind=db) :: P0=0 !< Initial canonical angular momentum for distribution 2 [kg m^2/s]
REAL(kind=db) :: temp !< Initial temperature of plasma [K] for distribtype=1
REAL(kind=db) :: weights_scale=1.0 !< Scale factor for the particle weights on restart (only for newres=.true.)
REAL(kind=db) :: temprescale = -1.0 !< Factor used for temperature rescaling in case of a restart (<0 -> no rescaling) Currently not implemented
INTEGER :: samplefactor =-1 !< Factor used for the up-sampling of the particles number
!
! Fields parameters
REAL(kind=db) :: B0 !< Max magnitude of magnetic field [T] and normalisation factor for magnetic field
! If magnetfile ='' The magnetic field is one of a magnetic mirror with maximum amplitude on axis of B0
REAL(kind=db) :: Rcurv = 1.0 !< Magnetic field curvature coefficient
REAL(kind=db) :: Width = 1.0 !< Distance between two magnetic mirrors
REAL(kind=db), allocatable :: Bz(:), Br(:) !< Normalised magnetic field components
REAL(kind=db), allocatable :: Athet(:) !< Theta component of the magnetic vector potential Tm
INTEGER :: bscaling = -1 !< if >0 rescale the magnetic field read from h5 file before calculating value at grid points, if <0 rescale after interpolation, if = 0 doesn't rescale
INTEGER :: Dimensionality = 1 !< Defines the dimensionality of the solution 1: (r,z) 2: (r,theta) 3: (r,theta,z)
LOGICAL :: nlperiod(3)=(/.false.,.false.,.true./)!< Set periodic splines on or off (z,r,theta)
LOGICAL :: nlPhis= .TRUE. !< Calculate self consistent electric field flag
LOGICAL :: nlfreezephi= .FALSE. !< Freeze the Poisson solver to the field obtained at (re-)start
INTEGER :: femorder(3) = 2 !< FEM order (z,r,theta)
INTEGER :: ngauss(3) = 6 !< Number of gauss points for FEM integration (z,r,theta)
LOGICAL :: nlppform =.TRUE. !< Defines if spline evaluation is done using ppform (faster with true)
INTEGER, SAVE :: nrank(3) !< Number of splines in all directions (z,r,theta)
TYPE(spline2d), SAVE :: splrz_ext !< Spline at r and z for total electric field
TYPE(spline2d), ALLOCATABLE,SAVE :: splrz(:) !< Spline at r and z for external electric field
REAL(kind=db) :: potinn=0 !< Electric potential at the inner metallic wall
REAL(kind=db) :: potout=0 !< Electric potential at the outer metallic wall
REAL(kind=db), allocatable :: Ez(:,:), Er(:,:), Ethet(:,:) !< Normalised electric field components ( ext+self )
REAL(kind=db), allocatable :: pot(:,:) !< Normalised electrostatic potential ( ext+self )
REAL(kind=db), allocatable :: Ezxt(:), Erxt(:), Ethetxt(:) !< Normalised external Electric field components
REAL(kind=db), allocatable :: potxt(:) !< Normalised external Electro static potential
!REAL(kind=db), allocatable :: npartspercell(:)
REAL(kind=db) :: radii(11) !< Inner and outer radius of cylinder and radii of fine mesh region [m]
INTEGER :: distribtype=1 !< Type of distribution function used to load the particles
!!1: gaussian, 2: Stable as defined in 4.95 of Davidson, 7 use particle input file
REAL(kind=db) :: lz(11) !< Lower and upper cylinder limits in z direction [m]
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 !< Total Number of grid intervals in z
INTEGER :: nnz(10) !< Number of grid intervals in z
INTEGER :: nsubz=10 !< Number of sub-intervals in z
INTEGER :: nr !< Total number of grid intervals in r
INTEGER :: nnr(10) !< Number of grid intervals in r in each subdomain
INTEGER :: nsubr=10 !< Number of sub-intervals in r
INTEGER :: nthet=1 !< Total Number of grid intervals in theta
REAL(kind=db) :: dz(10) !< Cell size in z
REAL(kind=db) :: dr(10) !< Cell size in r for each region
Real(kind=db) :: dthet !< Cell size in theta
REAL(kind=db) :: invdz(10), invdr(10) !< inverse of the grid cell step
Real(kind=db) :: invdthet !< inverse of the theta grid cell step
REAL(kind=db), ALLOCATABLE :: zgrid(:) !< Nodes positions in longitudinal direction
REAL(kind=db), ALLOCATABLE :: rgrid(:) !< Nodes positions in radial direction
REAL(kind=db), ALLOCATABLE :: thetgrid(:) !< Nodes positions in azimuthal direction
REAL(kind=db), ALLOCATABLE :: gradthet(:,:) !< Operator for the derivative wrt theta
REAL(kind=db) :: bnorm,enorm,vnorm,tnorm,rnorm,phinorm,qnorm !< Normalization constants
REAL(kind=db) :: omegac !< cylotronic frequency at B0 [1/s]
REAL(kind=db) :: omegap !< Plasma frequency at n0 [1/s]
INTEGER, DIMENSION(:), ALLOCATABLE :: Zbounds !< Index of bounds for local processus in Z direction for MPI decomposition
CONTAINS
!
!================================================================================
SUBROUTINE basic_data
!
! Define basic data
!
use mpihelper
USE omp_lib
Use random_mod
IMPLICIT NONE
!
! Local vars and arrays
CHARACTER(len=256) :: inputfilename
INTEGER :: i, nbprocs
!
NAMELIST /BASIC/ job_time, extra_time, nrun, tmax, dt, nlres, nlsave, newres, nlxg, &
& nplasma, potinn, potout, B0, lz, n0, nz, nthet, nnz, 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, magnetfile, weights_scale, nlfreezephi, nbaddtestspecies, &
& addedtestspecfile, bscaling, Dimensionality
!________________________________________________________________________________
! 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
CALL mpitypes_init ! 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"
WRITE(*,*)"db given kind", db
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
! Total number of intervals
nr=sum(nnr)
if (any(nnz.gt.0)) then
nz=sum(nnz)
else
nnz(1)=nz
end if
! Normalisation constants
if(nplasma .gt. 0) then
qsim=pi*(plasmadim(2)-plasmadim(1))*(plasmadim(4)**2-plasmadim(3)**2)*n0*elchar/nplasma
else
qsim=sign(n0,elchar)
end if
msim=abs(qsim)/elchar*partmass
vnorm=vlight
omegac=sign(elchar,qsim)/partmass*B0
omegap=sqrt(elchar**2*abs(n0)/partmass/eps_0)
tnorm=min(abs(1/omegac),abs(1/omegap))
rnorm=vnorm*tnorm
bnorm=B0
enorm=vlight*bnorm
phinorm=enorm*rnorm
! Normalised boundary conditions
potinn=potinn/phinorm
potout=potout/phinorm
! Normalised dt
dt=dt/tnorm
! Characteristic frequencies and normalised volume
IF(mpirank .eq. 0) THEN
IF(abs(omegap).GT. abs(omegac)) THEN
WRITE(*,'(a,3(1pe12.3))') 'omegap, omegac, omegap/omegac', omegap, omegac, omegap/omegac
ELSE
WRITE(*,'(a,3(1pe12.3))') 'omegap, omegac, omegac/omegap', omegap, omegac, omegac/omegap
END IF
END IF
! Construction of the mesh rgrid in r and zgrid in z and its normalisation
CALL mesh
rgrid=rgrid/rnorm
zgrid=zgrid/rnorm
dz=dz/rnorm
dr=dr/rnorm
Where(dr.gt.0) invdr=1/dr
Where(dz.gt.0) invdz=1/dz
! Initialize random number generator
nbprocs = omp_get_max_threads()
allocate(seed(ran_s,nbprocs), ran_index(nbprocs), ran_array(ran_k,nbprocs))
IF(.false.) then
call date_and_time(time=random_seed_str)
CALL MPI_BCAST(random_seed_str,10,MPI_CHARACTER,0,MPI_COMM_WORLD,ierr)
write(*,*) "MPI seed:", mpirank, random_seed_str
end if
Do i=1,nbprocs
! Generate seed from the default seed-string in random module
CALL decimal_to_seed(random_seed_str, seed(:,i))
! Generate a different seed for each processor from the mother seed
CALL next_seed(mpirank*nbprocs+i,seed(:,i))
! Initialize the random array (first hundred numbers)
CALL random_init(seed(:,i), ran_index(i), ran_array(:,i))
end do
!
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
!================================================================================
!---------------------------------------------------------------------------
!> @author
!> Patryk Kaminski EPFL/SPC
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Creates the mesh in r, theta and z direction for calculating the electric and magnetic fields.
!---------------------------------------------------------------------------
SUBROUTINE mesh
INTEGER :: j,i,k
ALLOCATE(zgrid(0:nz),rgrid(0:nr),thetgrid(0:nthet-1))
! Mesh in z
k=0
nsubz=count(nnz.gt.0)
zgrid(0)=lz(1)
do i=1,nsubz
dz(i)=(lz(i+1)-lz(i))/nnz(i)
if (nnz(i).gt.0) then
DO j=1,nnz(i)
zgrid(j+k)=lz(i)+j*dz(i)
END DO
end if
k=k+nnz(i)
end do
! Mesh in r
nsubr=count(nnr.gt.0)
k=0
rgrid(0)=radii(1)
do i=1,nsubr
dr(i)=(radii(i+1)-radii(i))/nnr(i)
if (nnr(i).gt.0) then
DO j=1,nnr(i)
rgrid(j+k)=radii(i)+j*dr(i)
END DO
end if
k=k+nnr(i)
end do
! Mesh in theta
dthet=0.0_db
invdthet=0.0_db
if(nthet.gt.0) then
dthet=2*pi/nthet
invdthet=1/dthet
do i=0,nthet-1
thetgrid(i)=i*dthet
end do
end if
! Derivative operator wrt theta
ALLOCATE(gradthet(nthet,nthet))
gradthet=0
IF(nthet .gt. 1) THEN
DO k=2,nthet-1
gradthet(k,k-1)=-0.5*invdthet
gradthet(k,k+1)=0.5*invdthet
END DO
gradthet(1,2)=0.5*invdthet
gradthet(1,nthet)=-0.5*invdthet
gradthet(nthet,1)=0.5*invdthet
gradthet(nthet,nthet-1)=-0.5*invdthet
END IF
END SUBROUTINE mesh
END MODULE basic
c4science · Help