! 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("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, yield, 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('secondaryEmission') CALL config%get(object // '.yield', yield, found) IF (.NOT. found) CALL criticalError("missing parameter 'yield' for secondary emission", 'readBoundary') CALL config%get(object // '.electron', speciesName, found) IF (.NOT. found) CALL criticalError("missing parameter 'electron' for secondary emission", 'readBoundary') speciesID = speciesName2Index(speciesName) CALL initSEE(boundary(i)%bTypes(s)%obj, yield, speciesID) 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 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 volumne CALL config%get(object // '.volume', volume, found) !Rescale the volumne 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 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 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 / tauMin) 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