1489 lines
49 KiB
Fortran
1489 lines
49 KiB
Fortran
! moduleInput: Reads JSON configuration file
|
|
MODULE moduleInput
|
|
USE json_module
|
|
IMPLICIT NONE
|
|
|
|
CONTAINS
|
|
!Main routine to read the input JSON file
|
|
SUBROUTINE readConfig(inputFile)
|
|
USE json_module
|
|
USE moduleErrors
|
|
USE moduleBoundary
|
|
USE moduleOutput
|
|
USE moduleMesh
|
|
IMPLICIT NONE
|
|
|
|
CHARACTER(:), ALLOCATABLE, INTENT(in):: inputFile
|
|
TYPE(json_file):: config
|
|
|
|
!Initialize the json file variable
|
|
CALL config%initialize()
|
|
|
|
!Loads the config file
|
|
CALL verboseError('Loading input file...')
|
|
CALL config%load(filename = inputFile)
|
|
CALL checkStatus(config, "loading")
|
|
|
|
!Reads reference parameters
|
|
CALL verboseError('Reading Reference parameters...')
|
|
CALL readReference(config)
|
|
CALL checkStatus(config, "readReference")
|
|
|
|
!Reads output parameters
|
|
CALL verboseError('Reading Output parameters...')
|
|
CALL readOutput(config)
|
|
CALL checkStatus(config, "readOutput")
|
|
|
|
!Read species
|
|
CALL verboseError('Reading species information...')
|
|
CALL readSpecies(config)
|
|
CALL checkStatus(config, "readSpecies")
|
|
|
|
!Read boundaries
|
|
CALL verboseError('Reading boundary conditions...')
|
|
CALL readBoundary(config)
|
|
CALL checkStatus(config, "readBoundary")
|
|
|
|
!Read Geometry
|
|
CALL verboseError('Reading Geometry...')
|
|
CALL readGeometry(config)
|
|
CALL checkStatus(config, "readGeometry")
|
|
|
|
!Read solver parameters
|
|
CALL verboseError('Reading Solver parameters...')
|
|
CALL readSolver(config)
|
|
CALL checkStatus(config, "readSolver")
|
|
|
|
!Read interactions between species
|
|
CALL verboseError('Reading interaction between species...')
|
|
CALL readInteractions(config)
|
|
CALL checkStatus(config, "readInteractions")
|
|
|
|
!Read probes
|
|
CALL verboseError('Reading Probes...')
|
|
CALL readProbes(config)
|
|
CALL checkStatus(config, 'readProbes')
|
|
|
|
!Read initial state for species
|
|
CALL verboseError('Reading Initial state...')
|
|
CALL readInitial(config)
|
|
CALL checkStatus(config, "readInitial")
|
|
|
|
!Read injection of particles
|
|
CALL verboseError('Reading injection of particles from boundaries...')
|
|
CALL readInject(config)
|
|
CALL checkStatus(config, "readInject")
|
|
|
|
!Read average scheme
|
|
CALL verboseError('Reading average scheme...')
|
|
CALL readAverage(config)
|
|
CALL checkStatus(config, "readAverage")
|
|
|
|
!Read parallel parameters
|
|
CALL verboseError('Reading Parallel configuration...')
|
|
CALL readParallel(config)
|
|
CALL checkStatus(config, "readParallel")
|
|
|
|
END SUBROUTINE readConfig
|
|
|
|
!Checks the status of the JSON case file and, if failed, exits the execution.
|
|
SUBROUTINE checkStatus(config, step)
|
|
USE moduleErrors
|
|
USE json_module
|
|
IMPLICIT NONE
|
|
|
|
TYPE(json_file), INTENT(inout):: config
|
|
CHARACTER(LEN=*), INTENT(in):: step
|
|
|
|
IF (config%failed()) CALL criticalError("Error reading the JSON input file", TRIM(step))
|
|
|
|
END SUBROUTINE checkStatus
|
|
|
|
!Reads the reference parameters
|
|
SUBROUTINE readReference(config)
|
|
USE moduleRefParam
|
|
USE moduleConstParam
|
|
USE moduleErrors
|
|
USE json_module
|
|
IMPLICIT NONE
|
|
|
|
TYPE(json_file), INTENT(inout):: config
|
|
LOGICAL:: found
|
|
CHARACTER(:), ALLOCATABLE:: object
|
|
|
|
object = 'reference'
|
|
!Mandatory parameters that define the case and computes derived parameters
|
|
CALL config%get(object // '.density', n_ref, found)
|
|
IF (.NOT. found) CALL criticalError('Reference density not found','readReference')
|
|
|
|
CALL config%get(object // '.mass', m_ref, found)
|
|
IF (.NOT. found) CALL criticalError('Reference mass not found','readReference')
|
|
|
|
CALL config%get(object // '.temperature', T_ref, found)
|
|
IF (.NOT. found) CALL criticalError('Reference temperature not found','readReference')
|
|
|
|
!Derived parameters
|
|
L_ref = DSQRT(kb*T_ref*eps_0/n_ref)/qe !reference length
|
|
v_ref = DSQRT(kb*T_ref/m_ref) !reference velocity
|
|
ti_ref = L_ref/v_ref !reference time
|
|
Vol_ref = L_ref**3 !reference volume
|
|
EF_ref = qe*n_ref*L_ref/eps_0 !reference electric field
|
|
Volt_ref = EF_ref*L_ref !reference voltage
|
|
B_ref = m_ref / (ti_ref * qe) !reference magnetic field
|
|
|
|
!If a reference cross section is given, it is used
|
|
CALL config%get(object // '.sigmaVrel', sigmavRel_ref, found)
|
|
|
|
!If not, the reference radius is searched
|
|
IF (.NOT. found) THEN
|
|
CALL config%get(object // '.radius', r_ref, found)
|
|
IF (found) THEN
|
|
sigmaVrel_ref = PI*(r_ref+r_ref)**2*v_ref !reference cross section times velocity
|
|
|
|
ELSE
|
|
sigmaVrel_ref = L_ref**2 * v_ref
|
|
|
|
END IF
|
|
|
|
END IF
|
|
|
|
END SUBROUTINE readReference
|
|
|
|
!Reads the specific case parameters
|
|
SUBROUTINE readSolver(config)
|
|
USE moduleRefParam
|
|
USE moduleErrors
|
|
USE moduleCaseParam
|
|
USE moduleSolver
|
|
USE moduleSpecies
|
|
USE moduleCollisions
|
|
USE moduleOutput
|
|
USE moduleMesh
|
|
USE json_module
|
|
IMPLICIT NONE
|
|
|
|
TYPE(json_file), INTENT(inout):: config
|
|
LOGICAL:: found
|
|
CHARACTER(:), ALLOCATABLE:: object
|
|
!simulation final and initial times in [t]
|
|
REAL(8):: finalTime, initialTime
|
|
CHARACTER(:), ALLOCATABLE:: pusherType, WSType, EMType
|
|
REAL(8):: B(1:3)
|
|
INTEGER:: nTau, nSolver
|
|
INTEGER:: i, n
|
|
CHARACTER(2):: iString
|
|
CHARACTER(1):: tString
|
|
|
|
object = 'solver'
|
|
|
|
!Time parameters
|
|
CALL config%info(object // '.tau', found, n_children = nTau)
|
|
IF (.NOT. found .OR. nTau == 0) THEN
|
|
CALL criticalError('Required parameter tau not found','readSolver')
|
|
|
|
END IF
|
|
|
|
ALLOCATE(tau(1:nSpecies))
|
|
|
|
DO i = 1, nTau
|
|
WRITE(iString, '(I2)') i
|
|
CALL config%get(object // '.tau(' // TRIM(iString) // ')', tau(i), found)
|
|
|
|
END DO
|
|
|
|
tauMin = MINVAL(tau(1:nTau))
|
|
|
|
IF (nTau < nSpecies) THEN
|
|
CALL warningError('Using minimum time step for some species')
|
|
tau(nTau+1:nSpecies) = tauMin
|
|
|
|
END IF
|
|
|
|
!Gets the simulation final time
|
|
CALL config%get(object // '.finalTime', finalTime, found)
|
|
IF (.NOT. found) THEN
|
|
CALL criticalError('Required parameter finalTime not found','readSolver')
|
|
|
|
END IF
|
|
!Convert simulation time to number of iterations
|
|
tFinal = INT(finalTime / tauMin)
|
|
|
|
!Gets the simulation initial time
|
|
CALL config%get(object // '.initialTime', initialTime, found)
|
|
IF (found) THEN
|
|
tInitial = INT(initialTime / tauMin)
|
|
|
|
END IF
|
|
|
|
!Gest the pusher for each species
|
|
IF (mesh%dimen > 0) THEN
|
|
CALL config%info(object // '.pusher', found, n_children = nSolver)
|
|
IF (.NOT. found .OR. nSolver /= nSpecies) THEN
|
|
CALL criticalError('Required parameter pusher not found','readSolver')
|
|
|
|
END IF
|
|
|
|
END IF
|
|
|
|
!Allocates all the pushers for particles
|
|
ALLOCATE(solver%pusher(1:nSpecies))
|
|
!Initialize pushers
|
|
DO i = 1, nSpecies
|
|
WRITE(iString, '(I2)') i
|
|
CALL config%get(object // '.pusher(' // TRIM(iString) // ')', pusherType, found)
|
|
|
|
!Associate the type of solver
|
|
CALL solver%pusher(i)%init(pusherType, tauMin, tau(i))
|
|
|
|
END DO
|
|
|
|
!Get the non-analogue scheme
|
|
CALL config%get(object // '.WeightingScheme', WSType, found)
|
|
CALL solver%initWS(WSType)
|
|
|
|
!Make tau(s) non-dimensional
|
|
tau = tau / ti_ref
|
|
tauMin = tauMin / ti_ref
|
|
|
|
!Set the format of output files accordint to iteration number
|
|
iterationDigits = INT(LOG10(REAL(tFinal))) + 1
|
|
WRITE(tString, '(I1)') iterationDigits
|
|
iterationFormat = "(I" // tString // "." // tString // ")"
|
|
|
|
!Get EM solver
|
|
CALL config%get(object // '.EMSolver', EMType, found)
|
|
IF (found) THEN
|
|
CALL solver%initEM(EMType)
|
|
SELECT CASE(EMType)
|
|
CASE("Electrostatic")
|
|
!Read BC
|
|
CALL readEMBoundary(config)
|
|
|
|
CASE("ElectrostaticBoltzmann")
|
|
!Read BC
|
|
CALL readEMBoundary(config)
|
|
|
|
CASE("ConstantB")
|
|
!Read BC
|
|
CALL readEMBoundary(config)
|
|
!Read constant magnetic field
|
|
DO i = 1, 3
|
|
WRITE(iString, '(i2)') i
|
|
CALL config%get(object // '.B(' // iString // ')', B(i), found)
|
|
IF (.NOT. found) THEN
|
|
CALL criticalError('Constant magnetic field not provided in direction ' // iString, 'readSolver')
|
|
|
|
END IF
|
|
|
|
END DO
|
|
|
|
!Non-dimensional magnetic field
|
|
B = B / B_ref
|
|
|
|
!Assign it to the nodes
|
|
DO n =1, mesh%numNodes
|
|
mesh%nodes(n)%obj%emData%B(1) = B(1)
|
|
mesh%nodes(n)%obj%emData%B(2) = B(2)
|
|
mesh%nodes(n)%obj%emData%B(3) = B(3)
|
|
|
|
END DO
|
|
|
|
CASE DEFAULT
|
|
CALL criticalError('EM Solver ' // EMType // ' not found', 'readSolver')
|
|
|
|
END SELECT
|
|
|
|
END IF
|
|
|
|
END SUBROUTINE readSolver
|
|
|
|
!Read the initial information for the species
|
|
SUBROUTINE readInitial(config)
|
|
USE moduleSpecies
|
|
USE moduleMesh
|
|
USE moduleOutput
|
|
USE moduleRefParam
|
|
USE moduleRandom
|
|
USE moduleProbe
|
|
USE json_module
|
|
IMPLICIT NONE
|
|
|
|
TYPE(json_file), INTENT(inout):: config
|
|
LOGICAL:: found
|
|
CHARACTER(:), ALLOCATABLE:: object
|
|
INTEGER:: nInitial
|
|
INTEGER:: i, p, e
|
|
CHARACTER(LEN=2):: iString
|
|
CHARACTER(:), ALLOCATABLE:: spName
|
|
INTEGER:: sp
|
|
CHARACTER(:), ALLOCATABLE:: spFile
|
|
CHARACTER(:), ALLOCATABLE:: filename
|
|
REAL(8), ALLOCATABLE, DIMENSION(:):: density, temperature
|
|
REAL(8), ALLOCATABLE, DIMENSION(:,:):: velocity
|
|
INTEGER, ALLOCATABLE, DIMENSION(:):: nodes
|
|
INTEGER:: nNodes
|
|
REAL(8), ALLOCATABLE, DIMENSION(:):: sourceScalar
|
|
REAL(8), ALLOCATABLE, DIMENSION(:,:):: sourceArray
|
|
!Density at the volume centroid
|
|
REAL(8):: densityCen
|
|
!Mean velocity and temperature at particle position
|
|
REAL(8):: velocityXi(1:3), temperatureXi
|
|
INTEGER:: nNewPart = 0
|
|
REAL(8):: weight = 0.D0
|
|
CLASS(meshCell), POINTER:: cell
|
|
TYPE(particle), POINTER:: partNew
|
|
REAL(8):: vTh
|
|
TYPE(lNode), POINTER:: partCurr, partNext
|
|
|
|
CALL config%info('solver.initial', found, n_children = nInitial)
|
|
|
|
IF (found) THEN
|
|
!Reads the information from initial species
|
|
DO i = 1, nInitial
|
|
WRITE(iString, '(I2)') i
|
|
object = 'solver.initial(' // iString // ')'
|
|
CALL config%get(object // '.species', spName, found)
|
|
sp = speciesName2Index(spName)
|
|
CALL config%get(object // '.file', spFile, found)
|
|
!Reads node values at the nodes
|
|
filename = path // spFile
|
|
CALL mesh%readInitial(filename, density, velocity, temperature)
|
|
!Check if initial number of particles is given
|
|
CALL config%get(object // '.particlesPerCell', nNewPart, found)
|
|
|
|
!For each volume in the node, create corresponding particles
|
|
DO e = 1, mesh%numCells
|
|
!Scale variables
|
|
!Density at centroid of cell
|
|
nNodes = mesh%cells(e)%obj%nNodes
|
|
nodes = mesh%cells(e)%obj%getNodes(nNodes)
|
|
ALLOCATE(sourceScalar(1:nNodes))
|
|
ALLOCATE(sourceArray(1:nNodes, 1:3))
|
|
sourceScalar = density(nodes)
|
|
densityCen = mesh%cells(e)%obj%gatherF((/ 0.D0, 0.D0, 0.D0 /), nNodes, sourceScalar)
|
|
|
|
!Calculate number of particles
|
|
IF (.NOT. found) THEN
|
|
nNewPart = FLOOR(densityCen * (mesh%cells(e)%obj%volume*Vol_ref) / species(sp)%obj%weight)
|
|
|
|
END IF
|
|
weight = densityCen * (mesh%cells(e)%obj%volume*Vol_ref) / REAL(nNewPart)
|
|
|
|
!Allocate new particles
|
|
DO p = 1, nNewPart
|
|
ALLOCATE(partNew)
|
|
partNew%species => species(sp)%obj
|
|
partNew%r = mesh%cells(e)%obj%randPos()
|
|
partNew%Xi = mesh%cells(e)%obj%phy2log(partNew%r)
|
|
!Get mean velocity at particle position
|
|
sourceArray(:,1) = velocity((nodes), 1)
|
|
sourceArray(:,2) = velocity((nodes), 2)
|
|
sourceArray(:,3) = velocity((nodes), 3)
|
|
|
|
velocityXi = mesh%cells(e)%obj%gatherF(partNew%Xi, nNodes, sourceArray)
|
|
velocityXi = velocityXi / v_ref
|
|
|
|
!Get temperature at particle position
|
|
sourceScalar = temperature(nodes)
|
|
temperatureXi = mesh%cells(e)%obj%gatherF(partNew%Xi, nNodes, sourceScalar)
|
|
temperatureXi = temperatureXi / T_ref
|
|
|
|
vTh = DSQRT(temperatureXi / species(sp)%obj%m)
|
|
partNew%v(1) = velocityXi(1) + vTh*randomMaxwellian()
|
|
partNew%v(2) = velocityXi(2) + vTh*randomMaxwellian()
|
|
partNew%v(3) = velocityXi(3) + vTh*randomMaxwellian()
|
|
|
|
partNew%cell = e
|
|
IF (doubleMesh) THEN
|
|
partNew%cellColl = findCellBrute(meshColl, partNew%r)
|
|
|
|
ELSE
|
|
partNew%cellColl = partNew%cell
|
|
|
|
END IF
|
|
|
|
partNew%n_in = .TRUE.
|
|
|
|
partNew%weight = weight
|
|
|
|
!Assign particle to temporal list of particles
|
|
CALL partInitial%add(partNew)
|
|
|
|
!Assign particle to list in volume
|
|
IF (listInCells) THEN
|
|
cell => meshforMCC%cells(partNew%cellColl)%obj
|
|
CALL OMP_SET_LOCK(cell%lock)
|
|
CALL cell%listPart_in(sp)%add(partNew)
|
|
cell%totalWeight(sp) = cell%totalWeight(sp) + partNew%weight
|
|
CALL OMP_UNSET_LOCK(cell%lock)
|
|
|
|
END IF
|
|
|
|
!Assign particles to probes
|
|
CALL doProbes(partNew)
|
|
|
|
END DO
|
|
|
|
DEALLOCATE(sourceScalar, sourceArray)
|
|
|
|
END DO
|
|
|
|
END DO
|
|
|
|
!Convert temporal list of particles into initial partOld array
|
|
!Deallocate the list of initial particles
|
|
nNewPart = partInitial%amount
|
|
IF (nNewPart > 0) THEN
|
|
ALLOCATE(partOld(1:nNewPart))
|
|
partCurr => partInitial%head
|
|
DO p = 1, nNewPart
|
|
partNext => partCurr%next
|
|
partOld(p) = partCurr%part
|
|
DEALLOCATE(partCurr)
|
|
partCurr => partNext
|
|
|
|
END DO
|
|
|
|
IF (ASSOCIATED(partInitial%head)) NULLIFY(partInitial%head)
|
|
IF (ASSOCIATED(partInitial%tail)) NULLIFY(partInitial%tail)
|
|
partInitial%amount = 0
|
|
|
|
END IF
|
|
|
|
nPartOld = SIZE(partOld)
|
|
|
|
END IF
|
|
|
|
END SUBROUTINE readInitial
|
|
|
|
!Reads configuration for the output files
|
|
SUBROUTINE readOutput(config)
|
|
USE moduleErrors
|
|
USE moduleOutput
|
|
USE json_module
|
|
IMPLICIT NONE
|
|
|
|
TYPE(json_file), INTENT(inout):: config
|
|
LOGICAL:: found
|
|
CHARACTER(:), ALLOCATABLE:: object
|
|
CHARACTER(:), ALLOCATABLE:: baseName
|
|
CHARACTER(8) :: date_now=''
|
|
CHARACTER(10) :: time_now=''
|
|
|
|
object = 'output'
|
|
CALL config%get(object // '.path', path, found)
|
|
CALL config%get(object // '.triggerOutput', triggerOutput, found)
|
|
IF (.NOT. found) THEN
|
|
triggerOutput = 100
|
|
CALL warningError('Using default trigger for output file of 100 iterations')
|
|
|
|
END IF
|
|
|
|
!Gets actual date and time
|
|
CALL DATE_AND_TIME(date_now, time_now)
|
|
|
|
!Gets the basename of the folder
|
|
CALL config%get(object // '.folder', baseName, found)
|
|
IF (found) THEN
|
|
folder = baseName // '_'
|
|
|
|
END IF
|
|
|
|
!Compose the folder name
|
|
folder = folder // date_now(1:4) // '-' // date_now(5:6) // '-' // date_now(7:8) // '_' &
|
|
// time_now(1:2) // '.' // time_now(3:4) // '.' // time_now(5:6)
|
|
|
|
CALL config%get(object // '.cpuTime', timeOutput, found)
|
|
CALL config%get(object // '.numColl', collOutput, found)
|
|
CALL config%get(object // '.EMField', emOutput, found)
|
|
|
|
CALL config%get(object // '.triggerCPUTime', triggerCPUTime, found)
|
|
IF (.NOT. found) THEN
|
|
triggerCPUTime = triggerOutput
|
|
IF (timeOutput) CALL warningError('No triggerCPUTime found, using same vale as triggerOutput')
|
|
|
|
END IF
|
|
|
|
END SUBROUTINE readOutput
|
|
|
|
!Reads information about the case species
|
|
SUBROUTINE readSpecies(config)
|
|
USE moduleSpecies
|
|
USE moduleErrors
|
|
USE moduleRefParam
|
|
USE moduleList
|
|
USE json_module
|
|
IMPLICIT NONE
|
|
|
|
TYPE(json_file), INTENT(inout):: config
|
|
CHARACTER(2):: iString
|
|
CHARACTER(:), ALLOCATABLE:: object
|
|
CHARACTER(:), ALLOCATABLE:: speciesType
|
|
REAL(8):: mass, charge
|
|
LOGICAL:: found
|
|
INTEGER:: i
|
|
CHARACTER(:), ALLOCATABLE:: linkName
|
|
INTEGER:: linkID
|
|
|
|
!Gets the number of species
|
|
CALL config%info('species', found, n_children = nSpecies)
|
|
!Zero species means critical error
|
|
IF (nSpecies == 0) CALL criticalError("No species found", "configRead")
|
|
|
|
ALLOCATE(species(1:nSpecies))
|
|
|
|
!Reads information of individual species
|
|
DO i = 1, nSpecies
|
|
WRITE(iString, '(I2)') i
|
|
object = 'species(' // TRIM(iString) // ')'
|
|
CALL config%get(object // '.type', speciesType, found)
|
|
CALL config%get(object // '.mass', mass, found)
|
|
mass = mass/m_ref
|
|
|
|
!Allocate species depending on type and assign specific parameters
|
|
SELECT CASE(speciesType)
|
|
CASE ("neutral")
|
|
ALLOCATE(species(i)%obj, source=speciesNeutral())
|
|
|
|
CASE ("charged")
|
|
CALL config%get(object // '.charge', charge, found)
|
|
IF (.NOT. found) CALL criticalError("Required parameter charge not found for species " // object, 'readSpecies')
|
|
ALLOCATE(species(i)%obj, source=speciesCharged(q = charge, &
|
|
qm = charge/mass))
|
|
|
|
CASE DEFAULT
|
|
CALL warningError("Species " // speciesType // " not supported yet")
|
|
|
|
END SELECT
|
|
!Assign shared parameters for all species
|
|
CALL config%get(object // '.name', species(i)%obj%name, found)
|
|
CALL config%get(object // '.weight', species(i)%obj%weight, found)
|
|
species(i)%obj%n = i
|
|
species(i)%obj%m = mass
|
|
|
|
END DO
|
|
|
|
!Read relations between species
|
|
DO i = 1, nSpecies
|
|
WRITE(iString, '(I2)') i
|
|
object = 'species(' // TRIM(iString) // ')'
|
|
SELECT TYPE(sp => species(i)%obj)
|
|
TYPE IS (speciesNeutral)
|
|
!Get species linked ion
|
|
CALL config%get(object // '.ion', linkName, found)
|
|
IF (found) THEN
|
|
linkID = speciesName2Index(linkName)
|
|
sp%ion => species(linkID)%obj
|
|
|
|
END IF
|
|
|
|
TYPE IS (speciesCharged)
|
|
!Get species linked neutral
|
|
CALL config%get(object // '.neutral', linkName, found)
|
|
IF (found) THEN
|
|
linkID = speciesName2Index(linkName)
|
|
sp%neutral => species(linkID)%obj
|
|
|
|
END IF
|
|
|
|
!Gets species linked ion
|
|
CALL config%get(object // '.ion', linkName, found)
|
|
IF (found) THEN
|
|
linkID = speciesName2Index(linkName)
|
|
sp%ion => species(linkID)%obj
|
|
|
|
END IF
|
|
|
|
END SELECT
|
|
|
|
END DO
|
|
|
|
!Set number of particles to 0 for init state
|
|
nPartOld = 0
|
|
|
|
!Initialize the lock for the non-analogue (NA) list of particles
|
|
CALL OMP_INIT_LOCK(partWScheme%lock)
|
|
|
|
END SUBROUTINE readSpecies
|
|
|
|
!Reads information about interactions between species
|
|
SUBROUTINE readInteractions(config)
|
|
USE moduleSpecies
|
|
USE moduleList
|
|
USE moduleCollisions
|
|
USE moduleCoulomb
|
|
USE moduleErrors
|
|
USE moduleMesh
|
|
USE moduleCaseParam
|
|
USE moduleRefParam
|
|
USE OMP_LIB
|
|
USE json_module
|
|
IMPLICIT NONE
|
|
|
|
TYPE(json_file), INTENT(inout):: config
|
|
CHARACTER(2):: iString, kString
|
|
CHARACTER(:), ALLOCATABLE:: object
|
|
CHARACTER(:), ALLOCATABLE:: species_i, species_j
|
|
CHARACTER(:), ALLOCATABLE:: crossSecFile
|
|
CHARACTER(:), ALLOCATABLE:: crossSecFilePath
|
|
CHARACTER(:), ALLOCATABLE:: cType
|
|
LOGICAL:: found
|
|
INTEGER:: nPairs, nCollisions
|
|
INTEGER:: i, k, ij
|
|
INTEGER:: pt_i, pt_j
|
|
REAL(8):: energyThreshold, energyBinding
|
|
CHARACTER(:), ALLOCATABLE:: electron, electronSecondary
|
|
INTEGER:: e
|
|
CLASS(meshCell), POINTER:: cell
|
|
|
|
!Firstly, check if the object 'interactions' exists
|
|
CALL config%info('interactions', found)
|
|
IF (found) THEN
|
|
!Check if MC collisions have been defined
|
|
CALL config%info('interactions.collisions', found)
|
|
IF (found) THEN
|
|
doMCCollisions = .TRUE.
|
|
!Read collision time step
|
|
CALL config%info('interactions.timeStep', found)
|
|
IF (found) THEN
|
|
CALL config%get('interactions.timeStep', tauColl, found)
|
|
tauColl = tauColl / ti_ref
|
|
|
|
ELSE
|
|
tauColl = tauMin
|
|
|
|
END IF
|
|
|
|
IF (tauColl < tauMin) THEN
|
|
CALL criticalError('Collisional time step below minimum time step.', 'readInteractions')
|
|
|
|
END IF
|
|
|
|
everyColl = NINT(tauColl / tauMin)
|
|
|
|
!Inits the MCC matrix
|
|
CALL initInteractionMatrix(interactionMatrix)
|
|
|
|
!Path for collision cross-section data files
|
|
CALL config%get('interactions.folderCollisions', pathCollisions, found)
|
|
|
|
!Inits lock for list of particles
|
|
CALL OMP_INIT_LOCK(partCollisions%lock)
|
|
|
|
CALL config%info('interactions.collisions', found, n_children = nPairs)
|
|
DO i = 1, nPairs
|
|
WRITE(iString, '(I2)') i
|
|
object = 'interactions.collisions(' // TRIM(iString) // ')'
|
|
CALL config%get(object // '.species_i', species_i, found)
|
|
pt_i = speciesName2Index(species_i)
|
|
CALL config%get(object // '.species_j', species_j, found)
|
|
pt_j = speciesName2Index(species_j)
|
|
CALL config%info(object // '.cTypes', found, n_children = nCollisions)
|
|
ij = interactionIndex(pt_i,pt_j)
|
|
|
|
!Allocates the required number of collisions per each pair of species ij
|
|
CALL interactionMatrix(ij)%init(nCollisions, pt_i, pt_j)
|
|
|
|
DO k = 1, nCollisions
|
|
WRITE (kString, '(I2)') k
|
|
object = 'interactions.collisions(' // TRIM(iString) // ').cTypes(' // TRIM(kString) // ')'
|
|
!Reads the cross section file
|
|
CALL config%get(object // '.crossSection', crossSecFile, found)
|
|
crossSecFilePath = pathCollisions // crossSecFile
|
|
IF (.NOT. found) CALL criticalError('crossSection not found for ' // object, 'readInteractions')
|
|
!Reads the type of collision
|
|
CALL config%get(object // '.type', cType, found)
|
|
!Initialize collision type and reads required additional data
|
|
SELECT CASE(cType)
|
|
CASE ('elastic')
|
|
!Elastic collision
|
|
CALL initBinaryElastic(interactionMatrix(ij)%collisions(k)%obj, crossSecFilePath)
|
|
|
|
CASE ('chargeExchange')
|
|
!Resonant charge exchange
|
|
CALL initBinaryChargeExchange(interactionMatrix(ij)%collisions(k)%obj, crossSecFilePath)
|
|
|
|
CASE ('ionization')
|
|
!Electorn impact ionization
|
|
CALL config%get(object // '.energyThreshold', energyThreshold, found)
|
|
IF (.NOT. found) CALL criticalError('energyThreshold not found for collision' // object, 'readInteractions')
|
|
CALL config%get(object // '.electron', electron, found)
|
|
IF (.NOT. found) CALL criticalError('electron not found for collision' // object, 'readInteractions')
|
|
CALL config%get(object // '.electronSecondary', electronSecondary, found)
|
|
IF (found) THEN
|
|
CALL initBinaryIonization(interactionMatrix(ij)%collisions(k)%obj, &
|
|
crossSecFilePath, energyThreshold, electron, electronSecondary)
|
|
|
|
ELSE
|
|
CALL initBinaryIonization(interactionMatrix(ij)%collisions(k)%obj, &
|
|
crossSecFilePath, energyThreshold, electron)
|
|
|
|
END IF
|
|
|
|
CASE ('recombination')
|
|
!Electorn impact ionization
|
|
CALL config%get(object // '.energyBinding', energyBinding, found)
|
|
IF (.NOT. found) CALL criticalError('energyThreshold not found for collision' // object, 'readInteractions')
|
|
CALL config%get(object // '.electron', electron, found)
|
|
IF (.NOT. found) CALL criticalError('electron not found for collision' // object, 'readInteractions')
|
|
CALL initBinaryRecombination(interactionMatrix(ij)%collisions(k)%obj, &
|
|
crossSecFilePath, energyBinding, electron)
|
|
|
|
CASE DEFAULT
|
|
CALL criticalError('Collision type' // cType // 'not defined', 'readInteractions')
|
|
|
|
END SELECT
|
|
|
|
END DO
|
|
|
|
END DO
|
|
|
|
!Init the required arrays in each volume to account for MCC.
|
|
DO e = 1, meshForMCC%numCells
|
|
cell => meshForMCC%cells(e)%obj
|
|
|
|
!Allocate Maximum cross section per collision pair and assign the initial collision rate
|
|
ALLOCATE(cell%sigmaVrelMax(1:nCollPairs))
|
|
cell%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref)
|
|
|
|
IF (collOutput) THEN
|
|
ALLOCATE(cell%tallyColl(1:nCollPairs))
|
|
DO k = 1, nCollPairs
|
|
ALLOCATE(cell%tallyColl(k)%tally(1:interactionmatrix(k)%amount))
|
|
cell%tallyColl(k)%tally = 0
|
|
|
|
END DO
|
|
|
|
END IF
|
|
|
|
END DO
|
|
|
|
END IF
|
|
|
|
!Check if Coulomb scattering is defined
|
|
CALL config%info('interactions.Coulomb', found)
|
|
IF (found) THEN
|
|
|
|
CALL config%info('interactions.Coulomb', found, n_children = nPairs)
|
|
IF (nPairs > 0) THEN
|
|
nCoulombPairs = nPairs
|
|
doCoulombScattering = .TRUE.
|
|
ALLOCATE(coulombMatrix(1:nPairs))
|
|
DO i = 1, nPairs
|
|
WRITE(iString, '(I2)') i
|
|
object = 'interactions.Coulomb(' // TRIM(iString) // ')'
|
|
CALL config%get(object // '.species_i', species_i, found)
|
|
pt_i = speciesName2Index(species_i)
|
|
CALL config%get(object // '.species_j', species_j, found)
|
|
pt_j = speciesName2Index(species_j)
|
|
|
|
CALL coulombMatrix(i)%init(pt_i, pt_j)
|
|
|
|
END DO
|
|
|
|
END IF
|
|
|
|
END IF
|
|
|
|
END IF
|
|
|
|
listInCells = doMCCollisions .OR. doCoulombScattering
|
|
|
|
END SUBROUTINE readInteractions
|
|
|
|
!Reads boundary conditions for the mesh
|
|
SUBROUTINE readBoundary(config)
|
|
USE moduleBoundary
|
|
USE moduleErrors
|
|
USE moduleSpecies
|
|
USE moduleRefParam
|
|
USE moduleList, ONLY: partSurfaces
|
|
USE json_module
|
|
IMPLICIT NONE
|
|
|
|
TYPE(json_file), INTENT(inout):: config
|
|
INTEGER:: i, s
|
|
CHARACTER(2):: iString, sString
|
|
CHARACTER(:), ALLOCATABLE:: object, bType
|
|
REAL(8):: Tw, cw !Wall temperature and specific heat
|
|
!Neutral Properties
|
|
REAL(8):: m0, n0, T0
|
|
REAL(8), DIMENSION(:), ALLOCATABLE:: v0
|
|
REAL(8):: effTime
|
|
REAL(8):: eThreshold !Energy threshold
|
|
INTEGER:: speciesID, electronSecondaryID
|
|
CHARACTER(:), ALLOCATABLE:: speciesName, crossSection, electronSecondary
|
|
LOGICAL:: found
|
|
INTEGER:: nTypes
|
|
|
|
CALL config%info('boundary', found, n_children = nBoundary)
|
|
ALLOCATE(boundary(1:nBoundary))
|
|
DO i = 1, nBoundary
|
|
WRITE(iString, '(i2)') i
|
|
object = 'boundary(' // TRIM(iString) // ')'
|
|
|
|
boundary(i)%n = i
|
|
CALL config%get(object // '.name', boundary(i)%name, found)
|
|
CALL config%get(object // '.physicalSurface', boundary(i)%physicalSurface, found)
|
|
CALL config%info(object // '.bTypes', found, n_children = nTypes)
|
|
IF (nTypes /= nSpecies) CALL criticalError('Not enough boundary types defined in ' // object, 'readBoundary')
|
|
ALLOCATE(boundary(i)%bTypes(1:nSpecies))
|
|
DO s = 1, nSpecies
|
|
WRITE(sString,'(i2)') s
|
|
object = 'boundary(' // TRIM(iString) // ').bTypes(' // TRIM(sString) // ')'
|
|
CALL config%get(object // '.type', bType, found)
|
|
SELECT CASE(bType)
|
|
CASE('reflection')
|
|
ALLOCATE(boundaryReflection:: boundary(i)%bTypes(s)%obj)
|
|
|
|
CASE('absorption')
|
|
ALLOCATE(boundaryAbsorption:: boundary(i)%bTypes(s)%obj)
|
|
|
|
CASE('transparent')
|
|
ALLOCATE(boundaryTransparent:: boundary(i)%bTypes(s)%obj)
|
|
|
|
CASE('ionization')
|
|
!Neutral parameters
|
|
CALL config%get(object // '.neutral.ion', speciesName, found)
|
|
IF (.NOT. found) CALL criticalError("missing parameter 'ion' for neutrals in ionization", 'readBoundary')
|
|
speciesID = speciesName2Index(speciesName)
|
|
CALL config%get(object // '.neutral.mass', m0, found)
|
|
IF (.NOT. found) THEN
|
|
m0 = species(s)%obj%m*m_ref
|
|
END IF
|
|
CALL config%get(object // '.neutral.density', n0, found)
|
|
IF (.NOT. found) CALL criticalError("missing parameter 'density' for neutrals in ionization", 'readBoundary')
|
|
CALL config%get(object // '.neutral.velocity', v0, found)
|
|
IF (.NOT. found) CALL criticalError("missing parameter 'velocity' for neutrals in ionization", 'readBoundary')
|
|
CALL config%get(object // '.neutral.temperature', T0, found)
|
|
IF (.NOT. found) CALL criticalError("missing parameter 'temperature' for neutrals in ionization", 'readBoundary')
|
|
|
|
CALL config%get(object // '.effectiveTime', effTime, found)
|
|
IF (.NOT. found) CALL criticalError("missing parameter 'effectiveTime' for ionization", 'readBoundary')
|
|
|
|
CALL config%get(object // '.energyThreshold', eThreshold, found)
|
|
IF (.NOT. found) CALL criticalError("missing parameter 'eThreshold' in ionization", 'readBoundary')
|
|
|
|
CALL config%get(object // '.crossSection', crossSection, found)
|
|
IF (.NOT. found) CALL criticalError("missing parameter 'crossSection' for neutrals in ionization", 'readBoundary')
|
|
|
|
CALL config%get(object // '.electronSecondary', electronSecondary, found)
|
|
electronSecondaryID = speciesName2Index(electronSecondary)
|
|
IF (found) THEN
|
|
CALL initIonization(boundary(i)%bTypes(s)%obj, species(s)%obj%m, m0, n0, v0, T0, &
|
|
speciesID, effTime, crossSection, eThreshold,electronSecondaryID)
|
|
|
|
ELSE
|
|
CALL initIonization(boundary(i)%bTypes(s)%obj, species(s)%obj%m, m0, n0, v0, T0, &
|
|
speciesID, effTime, crossSection, eThreshold)
|
|
|
|
END IF
|
|
|
|
CASE('wallTemperature')
|
|
CALL config%get(object // '.temperature', Tw, found)
|
|
IF (.NOT. found) CALL criticalError("temperature not found for wallTemperature boundary type", 'readBoundary')
|
|
CALL config%get(object // '.specificHeat', cw, found)
|
|
IF (.NOT. found) CALL criticalError("specificHeat not found for wallTemperature boundary type", 'readBoundary')
|
|
|
|
CALL initWallTemperature(boundary(i)%bTypes(s)%obj, Tw, cw)
|
|
|
|
CASE('axis')
|
|
ALLOCATE(boundaryAxis:: boundary(i)%bTypes(s)%obj)
|
|
|
|
CASE DEFAULT
|
|
CALL criticalError('Boundary type ' // bType // ' undefined', 'readBoundary')
|
|
|
|
END SELECT
|
|
|
|
END DO
|
|
|
|
END DO
|
|
|
|
!Init the list of particles from surfaces
|
|
CALL OMP_INIT_LOCK(partSurfaces%lock)
|
|
|
|
END SUBROUTINE readBoundary
|
|
|
|
!Read the geometry (mesh) for the case
|
|
SUBROUTINE readGeometry(config)
|
|
USE moduleMesh
|
|
USE moduleMeshInputGmsh2, ONLY: initGmsh2
|
|
USE moduleMeshInputVTU, ONLY: initVTU
|
|
USE moduleMeshInput0D, ONLY: init0D
|
|
USE moduleMeshInputText, ONLY: initText
|
|
USE moduleMesh3DCart
|
|
USE moduleMesh2DCyl
|
|
USE moduleMesh2DCart
|
|
USE moduleMesh1DRad
|
|
USE moduleMesh1DCart
|
|
USE moduleErrors
|
|
USE moduleOutput
|
|
USE moduleRefParam
|
|
USE moduleSolver
|
|
USE json_module
|
|
IMPLICIT NONE
|
|
|
|
TYPE(json_file), INTENT(inout):: config
|
|
CHARACTER(:), ALLOCATABLE:: object
|
|
LOGICAL:: found
|
|
CHARACTER(:), ALLOCATABLE:: meshFormat, meshFile
|
|
REAL(8):: volume
|
|
|
|
object = 'geometry'
|
|
|
|
!Checks if a mesh for collisions has been defined
|
|
!The mesh will be initialized and reader in readGeometry
|
|
CALL config%info(object // '.meshCollisions', found)
|
|
IF (found) THEN
|
|
!Points meshForMCC to the specific mesh defined
|
|
meshForMCC => meshColl
|
|
doubleMesh = .TRUE.
|
|
|
|
ELSE
|
|
CALL config%info('interactions', found)
|
|
IF (found) THEN
|
|
!Points the meshForMCC pointer to the Particles Mesh
|
|
meshForMCC => mesh
|
|
|
|
END IF
|
|
doubleMesh = .FALSE.
|
|
|
|
END IF
|
|
|
|
!Get the dimension of the geometry
|
|
CALL config%get(object // '.dimension', mesh%dimen, found)
|
|
IF (.NOT. found) THEN
|
|
CALL criticalError('Required parameter geometry.dimension not found', 'readGeometry')
|
|
|
|
END IF
|
|
IF (doubleMesh) THEN
|
|
meshColl%dimen = mesh%dimen
|
|
|
|
END IF
|
|
|
|
|
|
SELECT CASE(mesh%dimen)
|
|
CASE (0)
|
|
CALL init0D(mesh)
|
|
|
|
!Read the 0D mesh
|
|
CALL mesh%readMesh(pathMeshParticle)
|
|
|
|
!Get the volume
|
|
CALL config%get(object // '.volume', volume, found)
|
|
!Rescale the volume
|
|
IF (found) THEN
|
|
mesh%cells(1)%obj%volume = mesh%cells(1)%obj%volume*volume / Vol_ref
|
|
mesh%nodes(1)%obj%v = mesh%cells(1)%obj%volume
|
|
|
|
END IF
|
|
|
|
CASE (1:3)
|
|
!Select the type of geometry
|
|
CALL config%get(object // '.type', mesh%geometry, found)
|
|
IF (doubleMesh) THEN
|
|
meshColl%geometry = mesh%geometry
|
|
|
|
END IF
|
|
|
|
!Check the consistency between geometry and dimension
|
|
SELECT CASE(mesh%dimen)
|
|
CASE(3)
|
|
IF (mesh%geometry /= 'Cart') THEN
|
|
CALL criticalError('No available geometry ' // mesh%geometry // ' in 3D', 'readGeometry')
|
|
|
|
END IF
|
|
|
|
CASE(2)
|
|
IF (mesh%geometry /= 'Cart' .AND. &
|
|
mesh%geometry /= 'Cyl') THEN
|
|
CALL criticalError('No available geometry ' // mesh%geometry // ' in 2D', 'readGeometry')
|
|
|
|
END IF
|
|
|
|
CASE(1)
|
|
IF (mesh%geometry /= 'Cart' .AND. &
|
|
mesh%geometry /= 'Rad') THEN
|
|
CALL criticalError('No available geometry ' // mesh%geometry // ' in 1D', 'readGeometry')
|
|
|
|
END IF
|
|
|
|
END SELECT
|
|
|
|
!Link the procedure to connect meshes
|
|
SELECT CASE(mesh%dimen)
|
|
CASE(3)
|
|
mesh%connectMesh => connectMesh3DCart
|
|
|
|
CASE(2)
|
|
SELECT CASE(mesh%geometry)
|
|
CASE("Cyl")
|
|
mesh%connectMesh => connectMesh2DCyl
|
|
|
|
CASE("Cart")
|
|
mesh%connectMesh => connectMesh2DCart
|
|
|
|
END SELECT
|
|
|
|
CASE(1)
|
|
SELECT CASE(mesh%geometry)
|
|
CASE("Rad")
|
|
mesh%connectMesh => connectMesh1DRad
|
|
|
|
CASE("Cart")
|
|
mesh%connectMesh => connectMesh1DCart
|
|
|
|
END SELECT
|
|
|
|
END SELECT
|
|
|
|
IF (doubleMesh) THEN
|
|
meshColl%connectMesh => mesh%connectMesh
|
|
|
|
END IF
|
|
|
|
!Get the format of mesh
|
|
CALL config%get(object // '.meshType', meshFormat, found)
|
|
SELECT CASE(meshFormat)
|
|
CASE ("gmsh2")
|
|
CALL initGmsh2(mesh)
|
|
IF (doubleMesh) THEN
|
|
CALL initGmsh2(meshColl)
|
|
|
|
END IF
|
|
|
|
CASE ("vtu")
|
|
CALL initVTU(mesh)
|
|
IF (doubleMesh) THEN
|
|
CALL initVTU(meshColl)
|
|
|
|
END IF
|
|
|
|
case ("text")
|
|
!Check if the geometry is right.
|
|
if (mesh%dimen /= 1) then
|
|
call criticalError("Text mesh is only allowed for 1D geometries", 'readGeometry')
|
|
|
|
end if
|
|
|
|
!Read the mesh
|
|
call initText(mesh)
|
|
if (doubleMesh) then
|
|
call initText(meshColl)
|
|
|
|
end if
|
|
|
|
CASE DEFAULT
|
|
CALL criticalError('Mesh format ' // meshFormat // ' not defined.', 'readGeometry')
|
|
|
|
END SELECT
|
|
|
|
!Reads the mesh file
|
|
CALL config%get(object // '.meshFile', meshFile, found)
|
|
pathMeshParticle = path // meshFile
|
|
CALL mesh%readMesh(pathMeshParticle)
|
|
DEALLOCATE(meshFile)
|
|
IF (doubleMesh) THEN
|
|
!Reads the mesh file for collisions
|
|
CALL config%get(object // '.meshCollisions', meshFile, found)
|
|
pathMeshColl = path // meshFile
|
|
CALL meshColl%readMesh(pathMeshColl)
|
|
|
|
END IF
|
|
|
|
CASE DEFAULT
|
|
CALL criticalError("Mesh dimension not recogniced", "readGeometry")
|
|
|
|
END SELECT
|
|
|
|
!Builds the K matrix for the Particles mesh
|
|
CALL mesh%constructGlobalK()
|
|
|
|
!Assign the procedure to find a volume for meshColl
|
|
IF (doubleMesh) THEN
|
|
findCellColl => findCellCollMesh
|
|
|
|
ELSE
|
|
findCellColl => findCellSameMesh
|
|
|
|
END IF
|
|
|
|
END SUBROUTINE readGeometry
|
|
|
|
SUBROUTINE readProbes(config)
|
|
USE moduleProbe
|
|
USE moduleErrors
|
|
USE json_module
|
|
IMPLICIT NONE
|
|
|
|
TYPE(json_file), INTENT(inout):: config
|
|
CHARACTER(:), ALLOCATABLE:: object
|
|
LOGICAL:: found
|
|
CHARACTER(2):: iString
|
|
INTEGER:: i
|
|
CHARACTER(:), ALLOCATABLE:: speciesName
|
|
REAL(8), ALLOCATABLE, DIMENSION(:):: r
|
|
REAL(8), ALLOCATABLE, DIMENSION(:):: v1, v2, v3
|
|
INTEGER, ALLOCATABLE, DIMENSION(:):: points
|
|
REAL(8):: everyTimeStep
|
|
|
|
CALL config%info('output.probes', found, n_children = nProbes)
|
|
|
|
IF (found) ALLOCATE(probe(1:nProbes))
|
|
|
|
DO i = 1, nProbes
|
|
WRITE(iString, '(I2)') i
|
|
object = 'output.probes(' // trim(iString) // ')'
|
|
|
|
CALL config%get(object // '.species', speciesName, found)
|
|
CALL config%get(object // '.position', r, found)
|
|
CALL config%get(object // '.velocity_1', v1, found)
|
|
CALL config%get(object // '.velocity_2', v2, found)
|
|
CALL config%get(object // '.velocity_3', v3, found)
|
|
CALL config%get(object // '.points', points, found)
|
|
CALL config%get(object // '.timeStep', everyTimeStep, found)
|
|
|
|
IF (ANY(points < 2)) CALL criticalError("Number of points in probe " // iString // " incorrect", 'readProbes')
|
|
|
|
CALL probe(i)%init(i, speciesName, r, v1, v2, v3, points, everyTimeStep)
|
|
|
|
END DO
|
|
|
|
END SUBROUTINE readProbes
|
|
|
|
SUBROUTINE readEMBoundary(config)
|
|
USE moduleMesh
|
|
USE moduleOutput
|
|
USE moduleErrors
|
|
USE moduleEM
|
|
USE moduleSpecies
|
|
USE json_module
|
|
IMPLICIT NONE
|
|
|
|
TYPE(json_file), INTENT(inout):: config
|
|
CHARACTER(:), ALLOCATABLE:: object
|
|
LOGICAL:: found
|
|
CHARACTER(:), ALLOCATABLE:: typeEM
|
|
REAL(8):: potential
|
|
INTEGER:: physicalSurface
|
|
CHARACTER(:), ALLOCATABLE:: temporalProfile, temporalProfilePath
|
|
INTEGER:: b, s, n, ni
|
|
CHARACTER(2):: bString
|
|
INTEGER:: info
|
|
EXTERNAL:: dgetrf
|
|
|
|
CALL config%info('boundaryEM', found, n_children = nBoundaryEM)
|
|
|
|
IF (found) THEN
|
|
ALLOCATE(boundaryEM(1:nBoundaryEM))
|
|
|
|
END IF
|
|
|
|
DO b = 1, nBoundaryEM
|
|
WRITE(bString, '(I2)') b
|
|
object = 'boundaryEM(' // TRIM(bString) // ')'
|
|
|
|
CALL config%get(object // '.type', typeEM, found)
|
|
|
|
SELECT CASE(typeEM)
|
|
CASE ("dirichlet")
|
|
CALL config%get(object // '.potential', potential, found)
|
|
IF (.NOT. found) THEN
|
|
CALL criticalError('Required parameter "potential" for Dirichlet boundary condition not found', 'readEMBoundary')
|
|
|
|
END IF
|
|
|
|
CALL config%get(object // '.physicalSurface', physicalSurface, found)
|
|
IF (.NOT. found) THEN
|
|
CALL criticalError('Required parameter "physicalSurface" for Dirichlet boundary condition not found', &
|
|
'readEMBoundary')
|
|
|
|
END IF
|
|
|
|
CALL initDirichlet(boundaryEM(b)%obj, physicalSurface, potential)
|
|
|
|
CASE ("dirichletTime")
|
|
CALL config%get(object // '.potential', potential, found)
|
|
IF (.NOT. found) THEN
|
|
CALL criticalError('Required parameter "potential" for Dirichlet Time boundary condition not found', &
|
|
'readEMBoundary')
|
|
|
|
END IF
|
|
|
|
CALL config%get(object // '.temporalProfile', temporalProfile, found)
|
|
IF (.NOT. found) THEN
|
|
CALL criticalError('Required parameter "temporalProfile" for Dirichlet Time boundary condition not found', &
|
|
'readEMBoundary')
|
|
|
|
END IF
|
|
temporalProfilePath = path // temporalProfile
|
|
|
|
CALL config%get(object // '.physicalSurface', physicalSurface, found)
|
|
IF (.NOT. found) THEN
|
|
CALL criticalError('Required parameter "physicalSurface" for Dirichlet Time boundary condition not found', &
|
|
'readEMBoundary')
|
|
|
|
END IF
|
|
|
|
CALL initDirichletTime(boundaryEM(b)%obj, physicalSurface, potential, temporalProfilePath)
|
|
|
|
CASE DEFAULT
|
|
CALL criticalError('Boundary type ' // typeEM // ' not yet supported', 'readEMBoundary')
|
|
|
|
END SELECT
|
|
|
|
END DO
|
|
|
|
ALLOCATE(qSpecies(1:nSpecies))
|
|
DO s = 1, nSpecies
|
|
SELECT TYPE(sp => species(s)%obj)
|
|
TYPE IS (speciesCharged)
|
|
qSpecies(s) = sp%q
|
|
|
|
CLASS DEFAULT
|
|
qSpecies(s) = 0.D0
|
|
|
|
END SELECT
|
|
|
|
END DO
|
|
|
|
! Modify K matrix due to boundary conditions
|
|
DO b = 1, nBoundaryEM
|
|
SELECT TYPE(boundary => boundaryEM(b)%obj)
|
|
TYPE IS(boundaryEMDirichlet)
|
|
DO n = 1, boundary%nNodes
|
|
ni = boundary%nodes(n)%obj%n
|
|
mesh%K(ni, :) = 0.D0
|
|
mesh%K(ni, ni) = 1.D0
|
|
|
|
END DO
|
|
|
|
TYPE IS(boundaryEMDirichletTime)
|
|
DO n = 1, boundary%nNodes
|
|
ni = boundary%nodes(n)%obj%n
|
|
mesh%K(ni, :) = 0.D0
|
|
mesh%K(ni, ni) = 1.D0
|
|
|
|
END DO
|
|
|
|
END SELECT
|
|
|
|
END DO
|
|
|
|
!Compute the PLU factorization of K once boundary conditions have been read
|
|
CALL dgetrf(mesh%numNodes, mesh%numNodes, mesh%K, mesh%numNodes, mesh%IPIV, info)
|
|
IF (info /= 0) THEN
|
|
CALL criticalError('Factorization of K matrix failed', 'readEMBoundary')
|
|
|
|
END IF
|
|
|
|
END SUBROUTINE readEMBoundary
|
|
|
|
!Reads the injection of particles from the boundaries
|
|
SUBROUTINE readInject(config)
|
|
USE moduleSpecies
|
|
USE moduleErrors
|
|
USE moduleInject
|
|
USE json_module
|
|
IMPLICIT NONE
|
|
|
|
TYPE(json_file), INTENT(inout):: config
|
|
INTEGER:: i
|
|
CHARACTER(2):: iString
|
|
CHARACTER(:), ALLOCATABLE:: object
|
|
LOGICAL:: found
|
|
CHARACTER(:), ALLOCATABLE:: speciesName
|
|
CHARACTER(:), ALLOCATABLE:: name
|
|
REAL(8):: v
|
|
REAL(8), ALLOCATABLE:: temperature(:), normal(:)
|
|
REAL(8):: flow
|
|
CHARACTER(:), ALLOCATABLE:: units
|
|
INTEGER:: physicalSurface
|
|
INTEGER:: particlesPerEdge
|
|
INTEGER:: sp
|
|
|
|
CALL config%info('inject', found, n_children = nInject)
|
|
ALLOCATE(inject(1:nInject))
|
|
nPartInj = 0
|
|
DO i = 1, nInject
|
|
WRITE(iString, '(i2)') i
|
|
object = 'inject(' // trim(iString) // ')'
|
|
|
|
!Find species
|
|
CALL config%get(object // '.species', speciesName, found)
|
|
sp = speciesName2Index(speciesName)
|
|
|
|
CALL config%get(object // '.name', name, found)
|
|
CALL config%get(object // '.v', v, found)
|
|
CALL config%get(object // '.T', temperature, found)
|
|
CALL config%get(object // '.n', normal, found)
|
|
IF (.NOT. found) THEN
|
|
ALLOCATE(normal(1:3))
|
|
normal = 0.D0
|
|
END IF
|
|
CALL config%get(object // '.flow', flow, found)
|
|
CALL config%get(object // '.units', units, found)
|
|
CALL config%get(object // '.physicalSurface', physicalSurface, found)
|
|
particlesPerEdge = 0
|
|
CALL config%get(object // '.particlesPerEdge', particlesPerEdge, found)
|
|
|
|
CALL inject(i)%init(i, v, normal, temperature, flow, units, sp, physicalSurface, particlesPerEdge)
|
|
|
|
CALL readVelDistr(config, inject(i), object)
|
|
|
|
END DO
|
|
|
|
!Allocate array for injected particles
|
|
IF (nPartInj > 0) THEN
|
|
ALLOCATE(partInj(1:nPartInj))
|
|
|
|
END IF
|
|
|
|
END SUBROUTINE readInject
|
|
|
|
SUBROUTINE readAverage(config)
|
|
USE moduleAverage
|
|
USE moduleCaseParam, ONLY: tauMin
|
|
USE moduleMesh, ONLY: mesh
|
|
USE moduleSpecies, ONLY: nSpecies
|
|
USE moduleRefParam, ONLY: ti_ref
|
|
IMPLICIT NONE
|
|
|
|
TYPE(json_file), INTENT(inout):: config
|
|
LOGICAL:: found
|
|
REAL(8):: tStart
|
|
INTEGER:: n
|
|
|
|
CALL config%info('average', found)
|
|
IF (found) THEN
|
|
useAverage = .TRUE.
|
|
CALL config%get('average.startTime', tStart, found)
|
|
|
|
IF (found) THEN
|
|
tAverageStart = INT(tStart / ti_ref / tauMin)
|
|
|
|
ELSE
|
|
tAverageStart = 0
|
|
|
|
END IF
|
|
|
|
ALLOCATE(averageScheme(1:mesh%numNodes))
|
|
|
|
DO n = 1, mesh%numNodes
|
|
ALLOCATE(averageScheme(n)%mean%output(1:nSpecies))
|
|
ALLOCATE(averageScheme(n)%deviation%output(1:nSpecies))
|
|
|
|
END DO
|
|
|
|
END IF
|
|
|
|
END SUBROUTINE readAverage
|
|
|
|
!Reads the velocity distribution functions for each inject
|
|
SUBROUTINE readVelDistr(config, inj, object)
|
|
USE moduleErrors
|
|
USE moduleInject
|
|
USE moduleSpecies
|
|
USE json_module
|
|
IMPLICIT NONE
|
|
|
|
TYPE(json_file), INTENT(inout):: config
|
|
TYPE(injectGeneric), INTENT(inout):: inj
|
|
CHARACTER(:), ALLOCATABLE, INTENT(in):: object
|
|
INTEGER:: i
|
|
CHARACTER(2):: iString
|
|
CHARACTER(:), ALLOCATABLE:: fvType
|
|
LOGICAL:: found
|
|
REAL(8):: v, temperature, m
|
|
|
|
!Reads species mass
|
|
m = inj%species%m
|
|
!Reads distribution functions for velocity
|
|
DO i = 1, 3
|
|
WRITE(iString, '(i2)') i
|
|
CALL config%get(object // '.velDist('// TRIM(iString) //')', fvType, found)
|
|
IF (.NOT. found) CALL criticalError("No velocity distribution in direction " // iString // &
|
|
" found for " // object, 'readVelDistr')
|
|
|
|
SELECT CASE(fvType)
|
|
CASE ("Maxwellian")
|
|
temperature = inj%temperature(i)
|
|
CALL initVelDistMaxwellian(inj%v(i)%obj, temperature, m)
|
|
|
|
CASE ("Half-Maxwellian")
|
|
temperature = inj%temperature(i)
|
|
CALL initVelDistHalfMaxwellian(inj%v(i)%obj, temperature, m)
|
|
|
|
CASE ("Delta")
|
|
v = inj%vMod*inj%n(i)
|
|
CALL initVelDistDelta(inj%v(i)%obj)
|
|
|
|
CASE DEFAULT
|
|
CALL criticalError("No velocity distribution type " // fvType // " defined", 'readVelDistr')
|
|
|
|
END SELECT
|
|
|
|
END DO
|
|
|
|
END SUBROUTINE readVelDistr
|
|
|
|
!Reads configuration for parallel run
|
|
SUBROUTINE readParallel(config)
|
|
USE moduleParallel
|
|
USE moduleErrors
|
|
USE json_module
|
|
IMPLICIT NONE
|
|
|
|
TYPE(json_file), INTENT(inout):: config
|
|
CHARACTER(:), ALLOCATABLE:: object
|
|
LOGICAL:: found
|
|
|
|
!Reads OpenMP parameters
|
|
object = 'parallel.OpenMP'
|
|
|
|
CALL config%get(object // '.nThreads', openMP%nThreads, found)
|
|
|
|
IF (.NOT. found) THEN
|
|
openMP%nthreads = 8
|
|
CALL warningError('No OpenMP configuration detected, using 8 threads as default.')
|
|
|
|
END IF
|
|
|
|
CALL initParallel()
|
|
|
|
END SUBROUTINE readParallel
|
|
|
|
SUBROUTINE initOutput(inputFile)
|
|
USE moduleRefParam
|
|
USE moduleMesh, ONLY: mesh, doubleMesh, pathMeshParticle, pathMeshColl
|
|
USE moduleOutput, ONLY: path, folder
|
|
IMPLICIT NONE
|
|
|
|
CHARACTER(:), ALLOCATABLE, INTENT(in):: inputFile
|
|
INTEGER:: fileReference = 30
|
|
!If everything is correct, creates the output folder
|
|
CALL EXECUTE_COMMAND_LINE('mkdir ' // path // folder )
|
|
!Copies input file to output folder
|
|
CALL EXECUTE_COMMAND_LINE('cp ' // inputFile // ' ' // path // folder)
|
|
!Copies particle mesh
|
|
IF (mesh%dimen > 0) THEN
|
|
CALL EXECUTE_COMMAND_LINE('cp ' // pathMeshParticle // ' ' // path // folder)
|
|
IF (doubleMesh) THEN
|
|
CALL EXECUTE_COMMAND_LINE('cp ' // pathMeshColl // ' ' // path // folder)
|
|
|
|
END IF
|
|
|
|
END IF
|
|
|
|
! Write commit of fpakc
|
|
CALL SYSTEM('git rev-parse HEAD > ' // path // folder // '/' // 'fpakc_commit.txt')
|
|
|
|
! Write file with reference values
|
|
OPEN (fileReference, file=path // folder // '/' // 'reference.txt')
|
|
WRITE(fileReference, "(7(1X,A20))") 'L_ref', 'v_ref', 'ti_ref', 'Vol_ref', 'EF_ref', 'Volt_ref', 'B_ref'
|
|
WRITE(fileReference, "(7(1X,ES20.6E3))") L_ref, v_ref, ti_ref, Vol_ref, EF_ref, Volt_ref, B_ref
|
|
CLOSE(fileReference)
|
|
|
|
END SUBROUTINE initOutput
|
|
|
|
END MODULE moduleInput
|