1080 lines
37 KiB
Fortran
1080 lines
37 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 moduleInject
|
|
USE moduleOutput
|
|
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, "load")
|
|
|
|
!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 interactions between species
|
|
CALL verboseError('Reading interaction between species...')
|
|
CALL readInteractions(config)
|
|
CALL checkStatus(config, "readInteractions")
|
|
|
|
!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")
|
|
|
|
!Reads case parameters
|
|
CALL verboseError('Reading Case parameters...')
|
|
CALL readCase(config)
|
|
CALL checkStatus(config, "readCase")
|
|
|
|
!Read injection of particles
|
|
CALL verboseError('Reading injection of particles from boundaries...')
|
|
CALL readInject(config)
|
|
CALL checkStatus(config, "readInject")
|
|
|
|
!Read parallel parameters
|
|
CALL verboseError('Reading Parallel configuration...')
|
|
CALL readParallel(config)
|
|
CALL checkStatus(config, "readParallel")
|
|
|
|
!If everything is correct, creates the output folder
|
|
CALL EXECUTE_COMMAND_LINE('mkdir ' // path // folder )
|
|
|
|
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')
|
|
|
|
!If a reference cross section is given, it is used
|
|
CALL config%get(object // '.crossSection', sigma_ref, found)
|
|
|
|
!If not, the reference radius is searched
|
|
IF (.NOT. found) THEN
|
|
CALL config%get(object // '.radius', r_ref, found)
|
|
IF (found) THEN
|
|
sigma_ref = PI*(r_ref+r_ref)**2 !reference cross section
|
|
|
|
ELSE
|
|
sigma_ref = 0.D0 !Assume no collisions
|
|
|
|
END IF
|
|
|
|
END IF
|
|
|
|
!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
|
|
|
|
END SUBROUTINE readReference
|
|
|
|
!Reads the specific case parameters
|
|
SUBROUTINE readCase(config)
|
|
USE moduleRefParam
|
|
USE moduleErrors
|
|
USE moduleCaseParam
|
|
USE moduleSolver
|
|
USE moduleSpecies
|
|
USE moduleCollisions
|
|
USE moduleOutput
|
|
USE json_module
|
|
IMPLICIT NONE
|
|
|
|
TYPE(json_file), INTENT(inout):: config
|
|
LOGICAL:: found
|
|
CHARACTER(:), ALLOCATABLE:: object
|
|
REAL(8):: time !simulation time in [t]
|
|
CHARACTER(:), ALLOCATABLE:: pusherType, EMType, WSType
|
|
INTEGER:: nTau, nSolver
|
|
INTEGER:: i
|
|
CHARACTER(2):: iString
|
|
CHARACTER(1):: tString
|
|
|
|
object = 'case'
|
|
|
|
!Time parameters
|
|
CALL config%info(object // '.tau', found, n_children = nTau)
|
|
IF (.NOT. found .OR. nTau == 0) CALL criticalError('Required parameter tau not found','readCase')
|
|
ALLOCATE(tau(1:nSpecies))
|
|
DO i = 1, nTau
|
|
WRITE(iString, '(I2)') i
|
|
CALL config%get(object // '.tau(' // TRIM(iString) // ')', tau(i), found)
|
|
|
|
END DO
|
|
IF (nTau < nSpecies) THEN
|
|
CALL warningError('Using minimum time step for some species')
|
|
tau(nTau+1:nSpecies) = MINVAL(tau(1:nTau))
|
|
|
|
END IF
|
|
tauMin = MINVAL(tau)
|
|
|
|
!Gets the simulation time
|
|
CALL config%get(object // '.time', time, found)
|
|
IF (.NOT. found) CALL criticalError('Required parameter time not found','readCase')
|
|
!Convert simulation time to number of iterations
|
|
tmax = INT(time/tauMin)
|
|
|
|
!Gest the pusher for each species
|
|
CALL config%info(object // '.pusher', found, n_children = nSolver)
|
|
IF (.NOT. found .OR. nSolver /= nSpecies) CALL criticalError('Required parameter pusher not found','readCase')
|
|
!Allocates all the pushers for particles
|
|
ALLOCATE(solver%pusher(1:nSpecies))
|
|
!Initialize pushers
|
|
DO i = 1, nSolver
|
|
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
|
|
|
|
!Gets the solver for the electromagnetic field
|
|
CALL config%get(object // '.EMSolver', EMType, found)
|
|
CALL solver%initEM(EMType)
|
|
SELECT CASE(EMType)
|
|
CASE("Electrostatic")
|
|
CALL readEMBoundary(config)
|
|
|
|
END SELECT
|
|
|
|
!Gest the non-analogue scheme
|
|
CALL config%get(object // '.WeightingScheme', WSType, found)
|
|
CALL solver%initWS(WSType)
|
|
|
|
!Makes tau(s) non-dimensional
|
|
tau = tau / ti_ref
|
|
tauMin = tauMin / ti_ref
|
|
|
|
!Sets the format of output files accordint to iteration number
|
|
iterationDigits = INT(LOG10(REAL(tmax))) + 1
|
|
WRITE(tString, '(I1)') iterationDigits
|
|
iterationFormat = "(I" // tString // "." // tString // ")"
|
|
|
|
!Read initial state for species
|
|
CALL verboseError('Reading Initial state...')
|
|
CALL readInitial(config)
|
|
CALL checkStatus(config, "readInitial")
|
|
|
|
END SUBROUTINE readCase
|
|
|
|
!Reads the initial information for the species
|
|
SUBROUTINE readInitial(config)
|
|
USE moduleSpecies
|
|
USE moduleMesh
|
|
USE moduleOutput
|
|
USE moduleRefParam
|
|
USE moduleRandom
|
|
USE json_module
|
|
IMPLICIT NONE
|
|
|
|
TYPE(json_file), INTENT(inout):: config
|
|
LOGICAL:: found
|
|
CHARACTER(:), ALLOCATABLE:: object
|
|
INTEGER:: nInitial
|
|
INTEGER:: i, j, 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(:):: source, fPsi
|
|
!Density at the volume centroid
|
|
REAL(8):: densityCen
|
|
!Mean velocity and temperature at particle position
|
|
REAL(8):: velocityXi(1:3), temperatureXi
|
|
INTEGER:: nNewPart = 0.D0
|
|
TYPE(particle), POINTER:: partNew
|
|
REAL(8):: vTh
|
|
TYPE(lNode), POINTER:: partCurr, partNext
|
|
|
|
CALL config%info('case.initial', found, n_children = nInitial)
|
|
|
|
IF (found) THEN
|
|
!Reads the information from initial species
|
|
DO i = 1, nInitial
|
|
WRITE(iString, '(I2)') i
|
|
object = 'case.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(sp, filename, density, velocity, temperature)
|
|
!For each volume in the node, create corresponding particles
|
|
DO e = 1, mesh%numVols
|
|
!Scale variables
|
|
!Density at centroid of cell
|
|
nodes = mesh%vols(e)%obj%getNodes()
|
|
nNodes = SIZE(nodes)
|
|
!TODO: Procedure to obtain centroid from element (also for printing Electric Field)
|
|
fPsi = mesh%vols(e)%obj%fPsi((/0.D0, 0.D0, 0.D0/))
|
|
ALLOCATE(source(1:nNodes))
|
|
DO j = 1, nNodes
|
|
source(j) = density(nodes(j))
|
|
|
|
END DO
|
|
densityCen = DOT_PRODUCT(fPsi, source)
|
|
DEALLOCATE(fPsi)
|
|
|
|
!Calculate number of particles
|
|
nNewPart = INT(densityCen * (mesh%vols(e)%obj%volume*Vol_ref) / species(sp)%obj%weight)
|
|
|
|
!Allocate new particles
|
|
DO p = 1, nNewPart
|
|
ALLOCATE(partNew)
|
|
partNew%species => species(sp)%obj
|
|
partNew%r = mesh%vols(e)%obj%randPos()
|
|
partNew%xi = mesh%vols(e)%obj%phy2log(partNew%r)
|
|
!Get mean velocity at particle position
|
|
fPsi = mesh%vols(e)%obj%fPsi(partNew%xi)
|
|
DO j = 1, nNodes
|
|
source(j) = velocity(nodes(j), 1)
|
|
|
|
END DO
|
|
velocityXi(1) = DOT_PRODUCT(fPsi, source)
|
|
DO j = 1, nNodes
|
|
source(j) = velocity(nodes(j), 2)
|
|
|
|
END DO
|
|
velocityXi(2) = DOT_PRODUCT(fPsi, source)
|
|
DO j = 1, nNodes
|
|
source(j) = velocity(nodes(j), 3)
|
|
|
|
END DO
|
|
velocityXi(3) = DOT_PRODUCT(fPsi, source)
|
|
|
|
velocityXi = velocityXi / v_ref
|
|
|
|
!Get temperature at particle position
|
|
DO j = 1, nNodes
|
|
source(j) = temperature(nodes(j))
|
|
|
|
END DO
|
|
temperatureXi = DOT_PRODUCT(fPsi, source)
|
|
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%vol = e
|
|
IF (ASSOCIATED(meshForMCC, mesh)) THEN
|
|
partNew%volColl = partNew%vol
|
|
|
|
ELSE
|
|
partNew%volColl = findCellBrute(meshColl, partNew%r)
|
|
|
|
END IF
|
|
partNew%n_in = .TRUE.
|
|
partNew%weight = species(sp)%obj%weight
|
|
!If charged species, add qm to particle
|
|
SELECT TYPE(sp => species(sp)%obj)
|
|
TYPE IS (speciesCharged)
|
|
partNew%qm = sp%qm
|
|
|
|
CLASS DEFAULT
|
|
partNew%qm = 0.D0
|
|
|
|
END SELECT
|
|
|
|
!Assign particle to temporal list of particles
|
|
CALL partInitial%add(partNew)
|
|
|
|
END DO
|
|
|
|
DEALLOCATE(source)
|
|
|
|
END DO
|
|
|
|
END DO
|
|
|
|
!Convert temporal list of particles into initial partOld array
|
|
!Deallocates 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
|
|
|
|
!Reads relations between species
|
|
DO i = 1, nSpecies
|
|
WRITE(iString, '(I2)') i
|
|
object = 'species(' // TRIM(iString) // ')'
|
|
SELECT TYPE(sp => species(i)%obj)
|
|
TYPE IS (speciesNeutral)
|
|
!Gets 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)
|
|
!Gets 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(lockWScheme)
|
|
|
|
END SUBROUTINE readSpecies
|
|
|
|
!Reads information about interactions between species
|
|
SUBROUTINE readInteractions(config)
|
|
USE moduleSpecies
|
|
USE moduleList
|
|
USE moduleCollisions
|
|
USE moduleErrors
|
|
USE moduleMesh
|
|
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:: nInteractions, nCollisions
|
|
INTEGER:: i, k, ij
|
|
INTEGER:: pt_i, pt_j
|
|
REAL(8):: energyThreshold, energyBinding
|
|
CHARACTER(:), ALLOCATABLE:: electron
|
|
|
|
!Firstly, checks if the object 'interactions' exists
|
|
CALL config%info('interactions', found)
|
|
IF (found) THEN
|
|
!Checks if MC collisions have been defined
|
|
CALL config%info('interactions.collisions', found)
|
|
IF (found) THEN
|
|
!Checks if a mesh for collisions has been defined
|
|
!The mesh will be initialized and reader in readGeometry
|
|
CALL config%info('interactions.meshCollisions', found)
|
|
IF (found) THEN
|
|
!Points meshForMCC to the specific mesh defined
|
|
meshForMCC => meshColl
|
|
|
|
ELSE
|
|
!Points the meshForMCC pointer to the Particles Mesh
|
|
meshForMCC => mesh
|
|
|
|
END IF
|
|
|
|
!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(lockCollisions)
|
|
|
|
CALL config%info('interactions.collisions', found, n_children = nInteractions)
|
|
DO i = 1, nInteractions
|
|
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)
|
|
|
|
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, species(pt_i)%obj%m, species(pt_j)%obj%m)
|
|
|
|
CASE ('chargeExchange')
|
|
!Resonant charge exchange
|
|
CALL initBinaryChargeExchange(interactionMatrix(ij)%collisions(k)%obj, &
|
|
crossSecFilePath, species(pt_i)%obj%m, species(pt_j)%obj%m)
|
|
|
|
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 initBinaryIonization(interactionMatrix(ij)%collisions(k)%obj, &
|
|
crossSecFilePath, energyThreshold, species(pt_i)%obj%m, species(pt_j)%obj%m, electron)
|
|
|
|
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, species(pt_i)%obj%m, species(pt_j)%obj%m, electron)
|
|
|
|
CASE DEFAULT
|
|
CALL criticalError('Collision type' // cType // 'not defined yet', 'readInteractions')
|
|
|
|
END SELECT
|
|
|
|
END DO
|
|
|
|
END DO
|
|
|
|
END IF
|
|
|
|
END IF
|
|
|
|
END SUBROUTINE readInteractions
|
|
|
|
!Reads boundary conditions for the mesh
|
|
SUBROUTINE readBoundary(config)
|
|
USE moduleBoundary
|
|
USE moduleErrors
|
|
USE moduleSpecies
|
|
USE moduleRefParam
|
|
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
|
|
CHARACTER(:), ALLOCATABLE:: speciesName, crossSection
|
|
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 initIonization(boundary(i)%bTypes(s)%obj, species(s)%obj%m, m0, n0, v0, T0, &
|
|
speciesID, effTime, crossSection, eThreshold)
|
|
|
|
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
|
|
|
|
END SUBROUTINE readBoundary
|
|
|
|
!Read the geometry (mesh) for the case
|
|
SUBROUTINE readGeometry(config)
|
|
USE moduleMesh
|
|
USE moduleMeshInputGmsh2, ONLY: initGmsh2
|
|
USE moduleMeshInput0D, ONLY: init0D
|
|
USE moduleMesh3DCart, ONLY: connectMesh3DCart
|
|
USE moduleMesh2DCyl, ONLY: connectMesh2DCyl
|
|
USE moduleMesh2DCart, ONLY: connectMesh2DCart
|
|
USE moduleMesh1DRad, ONLY: connectMesh1DRad
|
|
USE moduleMesh1DCart, ONLY: connectMesh1DCart
|
|
USE moduleErrors
|
|
USE moduleOutput
|
|
USE moduleRefParam
|
|
USE json_module
|
|
IMPLICIT NONE
|
|
|
|
TYPE(json_file), INTENT(inout):: config
|
|
LOGICAL:: found
|
|
LOGICAL:: doubleMesh
|
|
CHARACTER(:), ALLOCATABLE:: meshFormat, meshFile
|
|
CHARACTER(:), ALLOCATABLE:: fullPath
|
|
REAL(8):: volume
|
|
|
|
!Firstly, indicates if a specific mesh for MC collisions is being use
|
|
doubleMesh = ASSOCIATED(meshForMCC, meshColl)
|
|
|
|
!Selects the type of geometry.
|
|
CALL config%get('geometry.type', mesh%geometry, found)
|
|
IF (doubleMesh) meshColl%geometry = mesh%geometry
|
|
|
|
!Gets the type of mesh
|
|
CALL config%get('geometry.meshType', meshFormat, found)
|
|
SELECT CASE(meshFormat)
|
|
CASE ("gmsh2")
|
|
CALL initGmsh2(mesh)
|
|
IF (doubleMesh) CALL initGmsh2(meshColl)
|
|
|
|
CASE ("0D")
|
|
CALL config%get('geometry.meshType', meshFormat, found)
|
|
CALL init0D(mesh)
|
|
|
|
CASE DEFAULT
|
|
CALL criticalError("Mesh format " // meshFormat // " not recogniced", "readGeometry")
|
|
|
|
END SELECT
|
|
|
|
!Reads the mesh file
|
|
CALL config%get('geometry.meshFile', meshFile, found)
|
|
fullpath = path // meshFile
|
|
CALL mesh%readMesh(fullPath)
|
|
DEALLOCATE(fullPath, meshFile)
|
|
IF (doubleMesh) THEN
|
|
!Reads the mesh file for collisions
|
|
CALL config%get('interactions.meshCollisions', meshFile, found)
|
|
fullpath = path // meshFile
|
|
CALL meshColl%readMesh(fullPath)
|
|
|
|
END IF
|
|
|
|
!Gets the volume for a 0D mesh
|
|
!TODO: Try to constrain this to the inout for 0D
|
|
IF (meshFormat == "0D") THEN
|
|
CALL config%get('geometry.volume', volume, found)
|
|
IF (found) THEN
|
|
mesh%vols(1)%obj%volume = mesh%vols(1)%obj%volume*volume / Vol_ref
|
|
mesh%nodes(1)%obj%v = mesh%vols(1)%obj%volume
|
|
|
|
END IF
|
|
|
|
END IF
|
|
|
|
!Creates the connectivity between elements
|
|
SELECT CASE(mesh%geometry)
|
|
CASE("3DCart")
|
|
mesh%connectMesh => connectMesh3DCart
|
|
|
|
CASE("2DCyl")
|
|
mesh%connectMesh => connectMesh2DCyl
|
|
|
|
CASE("2DCart")
|
|
mesh%connectMesh => connectMesh2DCart
|
|
|
|
CASE("1DRad")
|
|
mesh%connectMesh => connectMesh1DRad
|
|
|
|
CASE("1DCart")
|
|
mesh%connectMesh => connectMesh1DCart
|
|
|
|
CASE("0D")
|
|
mesh%connectMesh => NULL()
|
|
|
|
END SELECT
|
|
IF (ASSOCIATED(mesh%connectMesh)) CALL mesh%connectMesh
|
|
|
|
IF (doubleMesh) THEN
|
|
meshColl%connectMesh => mesh%connectMesh
|
|
CALL meshColl%connectMesh
|
|
|
|
END IF
|
|
|
|
!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 readEMBoundary(config)
|
|
USE moduleMesh
|
|
USE moduleOutput
|
|
USE moduleErrors
|
|
USE moduleEM
|
|
USE moduleRefParam
|
|
USE moduleSpecies
|
|
USE json_module
|
|
IMPLICIT NONE
|
|
|
|
TYPE(json_file), INTENT(inout):: config
|
|
CHARACTER(:), ALLOCATABLE:: object
|
|
LOGICAL:: found
|
|
CHARACTER(2):: istring
|
|
INTEGER:: i, e, s
|
|
INTEGER:: info
|
|
EXTERNAL:: dgetrf
|
|
|
|
CALL config%info('boundaryEM', found, n_children = nBoundaryEM)
|
|
|
|
IF (found) ALLOCATE(boundEM(1:nBoundaryEM))
|
|
|
|
DO i = 1, nBoundaryEM
|
|
WRITE(istring, '(i2)') i
|
|
object = 'boundaryEM(' // trim(istring) // ')'
|
|
|
|
CALL config%get(object // '.type', boundEM(i)%typeEM, found)
|
|
|
|
SELECT CASE(boundEM(i)%typeEM)
|
|
CASE ("dirichlet")
|
|
CALL config%get(object // '.potential', boundEM(i)%potential, found)
|
|
IF (.NOT. found) &
|
|
CALL criticalError('Required parameter "potential" for Dirichlet boundary condition not found', 'readEMBoundary')
|
|
boundEM(i)%potential = boundEM(i)%potential/Volt_ref
|
|
|
|
CALL config%get(object // '.physicalSurface', boundEM(i)%physicalSurface, found)
|
|
IF (.NOT. found) &
|
|
CALL criticalError('Required parameter "physicalSurface" for Dirichlet boundary condition not found', 'readEMBoundary')
|
|
|
|
CASE DEFAULT
|
|
CALL criticalError('Boundary type ' // boundEM(i)%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
|
|
|
|
!TODO: Improve this
|
|
IF (ALLOCATED(boundEM)) THEN
|
|
DO e = 1, mesh%numEdges
|
|
IF (ANY(mesh%edges(e)%obj%physicalSurface == boundEM(:)%physicalSurface)) THEN
|
|
DO i = 1, nBoundaryEM
|
|
IF (mesh%edges(e)%obj%physicalSurface == boundEM(i)%physicalSurface) THEN
|
|
CALL boundEM(i)%apply(mesh%edges(e)%obj)
|
|
|
|
END IF
|
|
END DO
|
|
END IF
|
|
END DO
|
|
END IF
|
|
|
|
!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) CALL criticalError('Factorization of K matrix failed', 'readEMBoundary')
|
|
|
|
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:: T(:), normal(:)
|
|
REAL(8):: flow
|
|
CHARACTER(:), ALLOCATABLE:: units
|
|
INTEGER:: physicalSurface
|
|
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', T, found)
|
|
CALL config%get(object // '.n', normal, found)
|
|
CALL config%get(object // '.flow', flow, found)
|
|
CALL config%get(object // '.units', units, found)
|
|
CALL config%get(object // '.physicalSurface', physicalSurface, found)
|
|
|
|
CALL inject(i)%init(i, v, normal, T, flow, units, sp, physicalSurface)
|
|
|
|
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
|
|
|
|
!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, T, 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")
|
|
v = inj%vMod*inj%n(i)
|
|
T = inj%T(i)
|
|
CALL initVelDistMaxwellian(inj%v(i)%obj, v, t, m)
|
|
|
|
CASE ("Delta")
|
|
v = inj%vMod*inj%n(i)
|
|
CALL initVelDistDelta(inj%v(i)%obj, v)
|
|
|
|
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
|
|
|
|
|
|
END MODULE moduleInput
|