JGonzalez/fpakc
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Restructuring the geometry and pushers

Browse source 
The geometry and push structure has been reworked to allow eassy adding
new pushers.

Documentation not updated yet.

Baseline for merging Cartesian pushers into one.
This commit is contained in:
Jorge Gonzalez 2022-04-08 19:06:12 +02:00
commit 5b5dadce39
18 changed files with 429 additions and 1052 deletions
Download patch file Download diff file Expand all files Collapse all files

296
src/modules/moduleInput.f90
View file

@ -40,16 +40,6 @@ MODULE moduleInput
CALL readSpecies(config)
CALL checkStatus(config, "readSpecies")
!Reads case parameters
CALL verboseError('Reading Case parameters...')
CALL readCase(config)
CALL checkStatus(config, "readCase")
!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)
@ -60,6 +50,16 @@ MODULE moduleInput
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)
@ -85,9 +85,12 @@ MODULE moduleInput
!Copies input file to output folder
CALL EXECUTE_COMMAND_LINE('cp ' // inputFile // ' ' // path // folder)
!Copies particle mesh
CALL EXECUTE_COMMAND_LINE('cp ' // pathMeshParticle // ' ' // path // folder)
IF (doubleMesh) THEN
CALL EXECUTE_COMMAND_LINE('cp ' // pathMeshColl // ' ' // path // folder)
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
@ -156,7 +159,7 @@ MODULE moduleInput
END SUBROUTINE readReference
!Reads the specific case parameters
SUBROUTINE readCase(config)
SUBROUTINE readSolver(config)
USE moduleRefParam
USE moduleErrors
USE moduleCaseParam
@ -164,6 +167,7 @@ MODULE moduleInput
USE moduleSpecies
USE moduleCollisions
USE moduleOutput
USE moduleMesh
USE json_module
IMPLICIT NONE
@ -172,47 +176,67 @@ MODULE moduleInput
CHARACTER(:), ALLOCATABLE:: object
!simulation final and initial times in [t]
REAL(8):: finalTime, initialTime
CHARACTER(:), ALLOCATABLE:: pusherType, WSType
CHARACTER(:), ALLOCATABLE:: pusherType, WSType, EMType
INTEGER:: nTau, nSolver
INTEGER:: i
CHARACTER(2):: iString
CHARACTER(1):: tString
object = 'case'
object = 'solver'
!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')
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) = MINVAL(tau(1:nTau))
tau(nTau+1:nSpecies) = tauMin
END IF
tauMin = MINVAL(tau)
!Gets the simulation final time
CALL config%get(object // '.finalTime', finalTime, found)
IF (.NOT. found) CALL criticalError('Required parameter finalTime not found','readCase')
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) tInitial = INT(initialTime / tauMin)
IF (found) THEN
tInitial = INT(initialTime / tauMin)
END IF
!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')
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, nSolver
DO i = 1, nSpecies
WRITE(iString, '(I2)') i
CALL config%get(object // '.pusher(' // TRIM(iString) // ')', pusherType, found)
@ -220,23 +244,39 @@ MODULE moduleInput
CALL solver%pusher(i)%init(pusherType, tauMin, tau(i))
END DO
!Gest the non-analogue scheme
!Get the non-analogue scheme
CALL config%get(object // '.WeightingScheme', WSType, found)
CALL solver%initWS(WSType)
!Makes tau(s) non-dimensional
!Make tau(s) non-dimensional
tau = tau / ti_ref
tauMin = tauMin / ti_ref
!Sets the format of output files accordint to iteration number
!Set the format of output files accordint to iteration number
iterationDigits = INT(LOG10(REAL(tFinal))) + 1
WRITE(tString, '(I1)') iterationDigits
iterationFormat = "(I" // tString // "." // tString // ")"
END SUBROUTINE readCase
!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)
!Reads the initial information for the species
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
@ -270,13 +310,13 @@ MODULE moduleInput
REAL(8):: vTh
TYPE(lNode), POINTER:: partCurr, partNext
CALL config%info('case.initial', found, n_children = nInitial)
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 = 'case.initial(' // iString // ')'
object = 'solver.initial(' // iString // ')'
CALL config%get(object // '.species', spName, found)
sp = speciesName2Index(spName)
CALL config%get(object // '.file', spFile, found)
@ -372,7 +412,7 @@ MODULE moduleInput
END DO
!Convert temporal list of particles into initial partOld array
!Deallocates the list of initial particles
!Deallocate the list of initial particles
nNewPart = partInitial%amount
IF (nNewPart > 0) THEN
ALLOCATE(partOld(1:nNewPart))
@ -504,13 +544,13 @@ MODULE moduleInput
END DO
!Reads relations between species
!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)
!Gets species linked ion
!Get species linked ion
CALL config%get(object // '.ion', linkName, found)
IF (found) THEN
linkID = speciesName2Index(linkName)
@ -519,7 +559,7 @@ MODULE moduleInput
END IF
TYPE IS (speciesCharged)
!Gets species linked neutral
!Get species linked neutral
CALL config%get(object // '.neutral', linkName, found)
IF (found) THEN
linkID = speciesName2Index(linkName)
@ -580,19 +620,6 @@ MODULE moduleInput
!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
!Reads collision time step
CALL config%info('interactions.timeStep', found)
IF (found) THEN
@ -792,11 +819,6 @@ MODULE moduleInput
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
@ -805,86 +827,122 @@ MODULE moduleInput
IMPLICIT NONE
TYPE(json_file), INTENT(inout):: config
CHARACTER(:), ALLOCATABLE:: object
LOGICAL:: found
CHARACTER(:), ALLOCATABLE:: meshFormat, meshFile, EMType
CHARACTER(:), ALLOCATABLE:: meshFormat, meshFile
REAL(8):: volume
!Firstly, indicates if a specific mesh for MC collisions is being use
doubleMesh = ASSOCIATED(meshForMCC, meshColl)
object = 'geometry'
!Selects the type of geometry.
CALL config%get('geometry.type', mesh%geometry, found)
IF (doubleMesh) meshColl%geometry = mesh%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.
!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)
ELSE
CALL config%info('interactions', found)
IF (found) THEN
!Points the meshForMCC pointer to the Particles Mesh
meshForMCC => mesh
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)
pathMeshParticle = path // meshFile
CALL mesh%readMesh(pathMeshParticle)
DEALLOCATE(meshFile)
IF (doubleMesh) THEN
!Reads the mesh file for collisions
CALL config%get('interactions.meshCollisions', meshFile, found)
pathMeshColl = path // meshFile
CALL meshColl%readMesh(pathMeshColl)
END IF
doubleMesh = .FALSE.
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)
!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 config%get(object // '.meshType', meshFormat, found)
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%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
CASE (1:3)
!Select the type of geometry
CALL config%get(object // '.type', mesh%geometry, found)
IF (doubleMesh) THEN
meshColl%geometry = mesh%geometry
!Creates the connectivity between elements
SELECT CASE(mesh%geometry)
CASE("3DCart")
mesh%connectMesh => connectMesh3DCart
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')
CASE("2DCyl")
mesh%connectMesh => connectMesh2DCyl
END IF
CASE("2DCart")
mesh%connectMesh => connectMesh2DCart
CASE(2)
IF (mesh%geometry /= 'Cart' .AND. &
mesh%geometry /= 'Cyl') THEN
CALL criticalError('No available geometry ' // mesh%geometry // ' in 2D', 'readGeometry')
CASE("1DRad")
mesh%connectMesh => connectMesh1DRad
END IF
CASE("1DCart")
mesh%connectMesh => connectMesh1DCart
CASE(1)
IF (mesh%geometry /= 'Cart' .AND. &
mesh%geometry /= 'Rad') THEN
CALL criticalError('No available geometry ' // mesh%geometry // ' in 1D', 'readGeometry')
CASE("0D")
mesh%connectMesh => NULL()
END IF
END SELECT
!Get the type 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
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
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()
@ -898,16 +956,6 @@ MODULE moduleInput
END IF
!Gest EM solver
CALL config%get('case.EMSolver', EMType, found)
CALL solver%initEM(EMType)
SELECT CASE(EMType)
CASE("Electrostatic")
!Reads BC
CALL readEMBoundary(config)
END SELECT
END SUBROUTINE readGeometry
SUBROUTINE readProbes(config)
@ -1010,7 +1058,6 @@ MODULE moduleInput
END DO
!TODO: Improve this
IF (ALLOCATED(boundEM)) THEN
DO e = 1, mesh%numEdges
IF (ANY(mesh%edges(e)%obj%physicalSurface == boundEM(:)%physicalSurface)) THEN
@ -1026,7 +1073,10 @@ MODULE moduleInput
!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')
IF (info /= 0) THEN
CALL criticalError('Factorization of K matrix failed', 'readEMBoundary')
END IF
END SUBROUTINE readEMBoundary