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

First commit of code.

Browse source 
New functionality:
- DSMC module:
  - 2D cyl geometry
    - GMSH file format
    - Elastic cross-section for Argon-Argon collisions.
    - Basic boundary conditions: reflection, absorption and axis
      symmetry.

Bugs fixed:

Other comments:
- Still searching for name.
This commit is contained in:
Jorge Gonzalez 2020-10-09 08:45:07 +02:00
commit bd7e8b040b
29 changed files with 4069 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

49
src/modules/makefile Normal file
View file

@ -0,0 +1,49 @@
OBJS = moduleCaseParam.o moduleCompTime.o moduleList.o\
moduleMesh.o moduleMeshCyl.o moduleMeshCylBoundary.o\
moduleMeshCylRead.o moduleOutput.o moduleInput.o \
moduleSolver.o moduleCollisions.o moduleTable.o
all: $(OBJS)
moduleCollisions.o: moduleTable.o moduleSpecies.o moduleRefParam.o moduleConstParam.o moduleMeshCyl.f95
$(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@
moduleInput.o: moduleRefParam.o moduleCaseParam.o moduleSolver.o moduleInject.o moduleBoundary.o moduleMesh.o moduleMeshCylRead.o moduleErrors.o moduleSpecies.o moduleInput.f95
$(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@
moduleInject.o: moduleSpecies.o moduleSolver.o moduleMesh.o moduleMeshCyl.o moduleInject.f95
$(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@
moduleList.o: moduleErrors.o moduleList.f95
$(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@
moduleMesh.o: moduleOutput.o moduleList.o moduleSpecies.o moduleMesh.f95
$(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@
moduleMeshCyl.o: moduleRefParam.o moduleCollisions.o moduleOutput.o moduleMesh.o moduleMeshCyl.f95
$(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@
moduleMeshCylBoundary.o: moduleMeshCyl.o moduleMeshCylBoundary.f95
$(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@
moduleMeshCylRead.o: moduleBoundary.o moduleMeshCyl.o moduleMeshCylBoundary.o moduleMeshCylRead.f95
$(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@
moduleOutput.o: moduleSpecies.o moduleRefParam.o moduleOutput.f95
$(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@
moduleRefParam.o: moduleConstParam.o moduleRefParam.f95
$(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@
moduleSpecies.o: moduleCaseParam.o moduleList.o moduleSpecies.f95
$(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@
moduleSolver.o: moduleSpecies.o moduleRefParam.o moduleMesh.o moduleOutput.o moduleSolver.f95
$(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@
moduleTable.o: moduleErrors.o moduleTable.f95
$(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@
%.o: %.f95
$(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@

36
src/modules/moduleBoundary.f95 Normal file
View file

@ -0,0 +1,36 @@
MODULE moduleBoundary
TYPE, PUBLIC:: boundaryGeneric
INTEGER:: id = 0
CHARACTER(:), ALLOCATABLE:: name
INTEGER:: physicalSurface = 0
CHARACTER(:), ALLOCATABLE:: boundaryType !TODO: substitute for extended types
END TYPE boundaryGeneric
TYPE:: boundaryCont
CLASS(boundaryGeneric), ALLOCATABLE:: obj
END TYPE boundaryCont
INTEGER:: nBoundary = 0
TYPE(boundaryCont), ALLOCATABLE:: boundary(:)
CONTAINS
FUNCTION getBoundaryId(physicalSurface) RESULT(id)
IMPLICIT NONE
INTEGER:: physicalSurface
INTEGER:: id
INTEGER:: i
id = 0
DO i = 1, nBoundary
IF (physicalSurface == boundary(i)%obj%physicalSurface) id = boundary(i)%obj%id
END DO
END FUNCTION getBoundaryId
END MODULE moduleBoundary

8
src/modules/moduleCaseParam.f95 Normal file
View file

@ -0,0 +1,8 @@
!Problems of the case
MODULE moduleCaseParam
!Maximum number of iterations and number of species
INTEGER:: tmax
REAL(8):: tau
END MODULE moduleCaseParam

152
src/modules/moduleCollisions.f95 Normal file
View file

@ -0,0 +1,152 @@
MODULE moduleCollisions
USE moduleTable
TYPE, ABSTRACT:: collisionBinary
TYPE(table1D):: crossSec
CONTAINS
PROCEDURE(initBinary_interface), PASS, DEFERRED:: init
PROCEDURE(collideBinary_interface), PASS, DEFERRED:: collide
END TYPE collisionBinary
ABSTRACT INTERFACE
SUBROUTINE initBinary_interface(self, crossSectionFilename)
IMPORT:: collisionBinary
CLASS(collisionBinary), INTENT(inout):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename
END SUBROUTINE
SUBROUTINE collideBinary_interface(self, sigmaVRelMax, sigmaVrel, part_i, part_j)
USE moduleSpecies
IMPORT:: collisionBinary
CLASS(collisionBinary), INTENT(in):: self
REAL(8), INTENT(in):: sigmaVrelMax
REAL(8), INTENT(out):: sigmaVrel
TYPE(particle), INTENT(inout):: part_i, part_j
END SUBROUTINE
END INTERFACE
!Container for collisions
TYPE:: collisionCont
CLASS(collisionBinary), ALLOCATABLE:: obj
END TYPE collisionCont
TYPE, EXTENDS(collisionBinary):: collisionBinaryElastic
CONTAINS
PROCEDURE, PASS:: init => initBinaryElastic
PROCEDURE, PASS:: collide => collideBinaryElastic
END TYPE collisionBinaryElastic
TYPE:: interactionsBinary
INTEGER:: amount
TYPE(collisionCont), ALLOCATABLE:: collisions(:)
CONTAINS
PROCEDURE, PASS:: init => initInteractionBinary
END TYPE interactionsBinary
!Collision 'Matrix'. A symmetric 2D matrix put into a 1D array to save memory
TYPE(interactionsBinary), ALLOCATABLE:: interactionMatrix(:)
!Folder for collision cross section tables
CHARACTER(:), ALLOCATABLE:: pathCollisions
CONTAINS
SUBROUTINE initInteractionMatrix(interactionMatrix)
USE moduleSpecies
IMPLICIT NONE
TYPE(interactionsBinary), INTENT(inout), ALLOCATABLE:: interactionMatrix(:)
INTEGER:: nInteractions
nInteractions = (nSpecies*(nSpecies+1))/2
ALLOCATE(interactionMatrix(1:nInteractions))
END SUBROUTINE initInteractionMatrix
FUNCTION interactionIndex(i,j) RESULT(k)
INTEGER:: i, j
INTEGER:: p
INTEGER:: k
k = i + j
p = (k + ABS(i - j))/2
k = k + (p*(p-3))/2
END FUNCTION interactionIndex
SUBROUTINE initInteractionBinary(self, amount)
IMPLICIT NONE
CLASS(interactionsBinary), INTENT(inout):: self
INTEGER, INTENT(in):: amount
INTEGER:: k
self%amount = amount
ALLOCATE(self%collisions(1:self%amount))
DO k= 1, self%amount
ALLOCATE(collisionBinaryElastic:: self%collisions(k)%obj)
END DO
END SUBROUTINE initInteractionBinary
SUBROUTINE initBinaryElastic(self, crossSectionFilename)
USE moduleTable
USE moduleRefParam
USE moduleConstParam
IMPLICIT NONE
CLASS(collisionBinaryElastic), INTENT(inout):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename
!Reads data from file
CALL self%crossSec%init(crossSectionFilename)
!Convert to no-dimensional units
CALL self%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2)
END SUBROUTINE
SUBROUTINE collideBinaryElastic(self, sigmaVrelMax, sigmaVrel, part_i, part_j)
USE moduleConstParam
USE moduleSpecies
USE moduleTable
IMPLICIT NONE
CLASS(collisionBinaryElastic), INTENT(in):: self
REAL(8), INTENT(in):: sigmaVrelMax
REAL(8), INTENT(out):: sigmaVrel
TYPE(particle), INTENT(inout):: part_i, part_j
REAL(8):: vRel !relative velocity
REAL(8):: vp_i, vp_j, v_i, v_j !post and pre-collision velocities
REAL(8):: alpha
!v relative (in units of [m][L]^2[s]^-2
vRel = species(1)%obj%m*NORM2(part_i%v - part_j%v)**2
sigmaVrel = self%crossSec%get(vRel)*vRel
IF (sigmaVrel/sigmaVrelMax > RAND()) THEN
!Applies the collision
v_i = NORM2(part_i%v)
v_j = NORM2(part_j%v)
vp_j = (v_i + v_j)*(1.D0+DSQRT(3.D0))/2.D0
vp_i = (v_i + v_j)*(DSQRT(3.D0)-1.D0)/2.D0
alpha = PI*RAND()
part_i%v(1) = v_i*DCOS(alpha)
part_i%v(2) = v_i*DSIN(alpha)
alpha = PI*RAND()
part_j%v(1) = v_j*DCOS(alpha)
part_j%v(2) = v_j*DSIN(alpha)
END IF
END SUBROUTINE
END MODULE moduleCollisions

10
src/modules/moduleCompTime.f95 Normal file
View file

@ -0,0 +1,10 @@
!Information to calculate computation time
MODULE moduleCompTime
REAL(8):: tStep=0.D0
REAL(8):: tPush=0.D0
REAL(8):: tReset=0.D0
REAL(8):: tColl=0.D0
REAL(8):: tWeight=0.D0
END MODULE moduleCompTime

14
src/modules/moduleConstParam.f95 Normal file
View file

@ -0,0 +1,14 @@
!Physical and mathematical constants
MODULE moduleConstParam
IMPLICIT NONE
PUBLIC
REAL(8), PARAMETER:: PI = 4.D0*DATAN(1.D0) !number pi
REAL(8), PARAMETER:: sccm2atomPerS = 4.5D17 !sccm to atom s^-1
REAL(8), PARAMETER:: eV2J = 1.60217662D-19 !Electron volt to Joule (elementary charge)
REAL(8), PARAMETER:: kb = 1.38064852D-23 !Boltzmann constants SI
! REAL(8), PARAMETER:: eps_0 = 8.8542D-12 !Epsilon_0 (Not used in Neutrals)
END MODULE moduleConstParam

43
src/modules/moduleErrors.f95 Normal file
View file

@ -0,0 +1,43 @@
!moduleErrors: Manages the different type of errors the program can produce.
! By errors we understand critical errors (that stop the program),
! warnings (that only output a message with a WARNING tag),
! o verbose outputs that can be used for debugging porpouse.
MODULE moduleErrors
CONTAINS
SUBROUTINE criticalError(msg, pgr)
IMPLICIT NONE
CHARACTER(*), INTENT(in):: msg, pgr
CHARACTER(:), ALLOCATABLE:: errorMsg
errorMsg = "CRITICAL error in " // pgr // " with message:" // NEW_LINE('A') // msg
ERROR STOP errorMsg
END SUBROUTINE criticalError
SUBROUTINE warningError(msg)
IMPLICIT NONE
CHARACTER(*), INTENT(in):: msg
CHARACTER(:), ALLOCATABLE:: errorMsg
errorMsg = "WARNING: " // msg
WRITE (*, '(A)') errorMsg
END SUBROUTINE warningError
SUBROUTINE verboseError(msg)
IMPLICIT NONE
CHARACTER(*), INTENT(in):: msg
CHARACTER(:), ALLOCATABLE:: errorMsg
errorMsg = msg
WRITE (*, '(A)') errorMsg
END SUBROUTINE verboseError
END MODULE moduleErrors

155
src/modules/moduleInject.f95 Normal file
View file

@ -0,0 +1,155 @@
!injection of particles
MODULE moduleInject
TYPE:: injectGeneric
INTEGER:: id
CHARACTER(:), ALLOCATABLE:: name
REAL(8):: v !Velocity (module)
REAL(8):: T(1:3) !Temperature
REAL(8):: n(1:3) !Direction of injection
INTEGER:: nParticles !Number of particles to introduce each time step
INTEGER:: pt !Species of injection
INTEGER, ALLOCATABLE:: edges(:) !Array with edges
CONTAINS
PROCEDURE, PASS:: init => initInject
PROCEDURE, PASS:: addParticles => addParticlesMaxwellian
END TYPE injectGeneric
INTEGER:: nInject
TYPE(injectGeneric), ALLOCATABLE:: inject(:)
CONTAINS
SUBROUTINE initInject(self, i, v, n, T, flow, pt, physicalSurface)
USE moduleMesh
USE moduleRefParam
USE moduleConstParam
USE moduleSpecies
IMPLICIT NONE
CLASS(injectGeneric), INTENT(inout):: self
INTEGER, INTENT(in):: i
REAL(8), INTENT(in):: v, n(1:3), T(1:3)
INTEGER, INTENT(in):: pt, physicalSurface
REAL(8), INTENT(in):: flow
INTEGER:: nEdges, e, et
INTEGER:: phSurface(1:mesh%numEdges)
self%id = i
self%v = v/v_ref
self%n = n
self%T = T/T_ref
self%nParticles = INT(flow*sccm2atomPerS*tau*ti_ref/species(pt)%obj%weight)
self%pt = pt
!Gets the edge elements from which particles are injected
!TODO: Improve this A LOT
DO e = 1, mesh%numEdges
phSurface(e) = mesh%edges(e)%obj%physicalSurface
END DO
nEdges = COUNT(phSurface == physicalSurface)
ALLOCATE(inject(i)%edges(1:nEdges))
et = 0
DO e=1, mesh%numEdges
IF (mesh%edges(e)%obj%physicalSurface == physicalSurface) THEN
et = et + 1
self%edges(et) = mesh%edges(e)%obj%n
END IF
END DO
nPartInj = nPartInj + self%nParticles
END SUBROUTINE
!Random velocity from maxwellian distribution
REAL(8) FUNCTION vBC(u, vth)
USE moduleConstParam, ONLY: PI
REAL(8), INTENT (in):: u, vth
REAL(8):: x, y
vBC=0.D0
x=RAND()
y=RAND()
vBC = u + vth*DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y)
END FUNCTION vBC
SUBROUTINE addParticlesMaxwellian(self)
USE moduleSpecies
USE moduleSolver
USE moduleMesh
USE moduleMeshCyl
IMPLICIT NONE
CLASS(injectGeneric), INTENT(in):: self
INTEGER:: randomEdge
REAL(8):: randomPos
REAL(8):: vVec(1:3), vTh(1:3)
!Edge nodes
INTEGER:: n1 = 0, n2 = 0
!Edge nodes coordinates
REAL(8):: p1(1:3) = 0.D0, p2(1:3) = 0.D0
INTEGER:: n
vVec = self%v * self%n
vTh = DSQRT(self%T/species(self%pt)%obj%m)
!Insert particles
DO n = 1, self%nParticles
!Select edge randomly from which inject particle
randomEdge = self%edges(INT(RAND()*(SIZE(self%edges)-1)+1))
!Get coordinates of edge nodes
SELECT TYPE(edge => mesh%edges(randomEdge)%obj)
CLASS IS (meshEdgeCyl)
n1 = edge%n1%n
n2 = edge%n2%n
END SELECT
SELECT TYPE(node => mesh%nodes(n1)%obj)
TYPE IS (meshNodeCyl)
p1(1) = node%z
p1(2) = node%r
END SELECT
SELECT TYPE(node => mesh%nodes(n2)%obj)
TYPE IS (meshNodeCyl)
p2(1) = node%z
p2(2) = node%r
END SELECT
!Volume associated to the edge:
IF (DOT_PRODUCT(self%n, mesh%edges(randomEdge)%obj%normal) >= 0.D0) THEN
part_inj(n)%e_p = mesh%edges(randomEdge)%obj%e1%n
ELSE
part_inj(n)%e_p = mesh%edges(randomEdge)%obj%e2%n
END IF
part_inj(n)%pt = self%pt
part_inj(n)%n_in = .TRUE.
part_inj(n)%v = (/ vBC(vVec(1), vTh(1)), &
vBC(vVec(2), vTh(2)), &
vBC(vVec(3), vTh(3)) /)
!Random position in edge
!TODO: Use edge coordinates and transformations for this process
randomPos = RAND()
part_inj(n)%r(1) = p1(1) + randomPos*(p2(1) - p1(1))
part_inj(n)%r(2) = p1(2) + randomPos*(p2(2) - p1(2))
part_inj(n)%r(3) = p1(3) + randomPos*(p2(3) - p1(3))
!Push new particle
CALL push(part_inj(n))
END DO
END SUBROUTINE addParticlesMaxwellian
END MODULE moduleInject

366
src/modules/moduleInput.f95 Normal file
View file

@ -0,0 +1,366 @@
! 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
IMPLICIT NONE
CHARACTER(:), ALLOCATABLE, INTENT(in):: inputFile
TYPE(json_file):: config
!Initialize the json file variable
CALL config%initialize()
!Loads the config file
CALL config%load(filename = inputFile)
!Reads reference parameters
CALL readReference(config)
!Reads case parameters
CALL readCase(config)
!Reads output parameters
CALL readOutput(config)
!Read species
CALL readSpecies(config)
!Read interactions between species
CALL readInteractions(config)
!Read boundaries
CALL readBoundary(config)
!Read Geometry
CALL readGeometry(config)
!Read injection of particles
CALL readInject(config)
END SUBROUTINE readConfig
!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')
CALL config%get(object // '.radius', r_ref, found)
IF (.NOT. found) CALL criticalError('Reference radius not found','readReference')
!Derived parameters
sigma_ref = PI*(r_ref+r_ref)**2 !reference cross section
L_ref = 1.D0/(sigma_ref*n_ref) !mean free path
ti_ref = L_ref/v_ref !reference time
Vol_ref = L_ref**3 !reference volume
v_ref = DSQRT(kb*T_ref/m_ref) !reference velocity
END SUBROUTINE readReference
!Reads the specific case parameters
SUBROUTINE readCase(config)
USE moduleRefParam
USE moduleErrors
USE moduleCaseParam
USE json_module
IMPLICIT NONE
TYPE(json_file), INTENT(inout):: config
LOGICAL:: found
CHARACTER(:), ALLOCATABLE:: object
REAL(8):: time !simulation time in [t]
object = 'case'
!Time parameters
CALL config%get(object // '.tau', tau, found)
IF (.NOT. found) CALL criticalError('Required parameter tau not found','readCase')
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/(ti_ref*tau))
END SUBROUTINE readCase
!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(8) :: date_now=''
CHARACTER(10) :: time_now=''
object = 'output'
CALL config%get(object // '.path', path, found)
CALL config%get(object // '.trigger', triggerOutput, found)
IF (.NOT. found) THEN
triggerOutput = 100
CALL warningError('Using default trigger for output file of 100 iterations')
END IF
!Creates output folder
!TODO: Add option for custon name output_folder
CALL DATE_AND_TIME(date_now, time_now)
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 SYSTEM('mkdir ' // path // folder )
CALL config%get(object // '.cpuTime', timeOutput, found)
CALL config%get(object // '.numColl', collOutput, found)
END SUBROUTINE readOutput
!Reads information about the case species
SUBROUTINE readSpecies(config)
USE moduleSpecies
USE moduleErrors
USE moduleRefParam
USE json_module
IMPLICIT NONE
TYPE(json_file), INTENT(inout):: config
CHARACTER(2):: iString
CHARACTER(:), ALLOCATABLE:: object
CHARACTER(:), ALLOCATABLE:: speciesType
REAL(8):: mass
LOGICAL:: found
INTEGER:: i
!Gets the number of species
CALL config%info('species', 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)
SELECT CASE(speciesType)
CASE ("neutral")
ALLOCATE(species(i)%obj, source=speciesNeutral())
!TODO: move to subroutine
CALL config%get(object // '.name', species(i)%obj%name, found)
CALL config%get(object // '.mass', mass, found)
CALL config%get(object // '.weight', species(i)%obj%weight, found)
species(i)%obj%pt = i
species(i)%obj%m = mass/m_ref
CASE DEFAULT
CALL warningError("Species " // speciesType // " not supported yet")
END SELECT
END DO
!Set number of particles to 0 for init state
!TODO: In a future, this should include the particles from init states
n_part_old = 0
END SUBROUTINE readSpecies
!Reads information about interactions between species
SUBROUTINE readInteractions(config)
USE moduleSpecies
USE moduleCollisions
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
INTEGER:: nInteractions, nCollisions
INTEGER:: i, k, ij
INTEGER:: pt_i, pt_j
CALL initInteractionMatrix(interactionMatrix)
CALL config%get('interactions.folderCollisions', pathCollisions)
CALL config%info('interactions.collisions', n_children = nInteractions)
DO i = 1, nInteractions
WRITE(iString, '(I2)') i
object = 'interactions.collisions(' // TRIM(iString) // ')'
CALL config%get(object // '.species_i', species_i)
pt_i = speciesName2Index(species_i)
CALL config%get(object // '.species_j', species_j)
pt_j = speciesName2Index(species_j)
CALL config%info(object // '.crossSections', n_children = nCollisions)
ij = interactionIndex(pt_i,pt_j)
CALL interactionMatrix(ij)%init(nCollisions)
DO k = 1, nCollisions
WRITE (kString, '(I2)') k
CALL config%get(object // '.crossSections(' // TRIM(kString)// ')', crossSecFile)
crossSecFilePath = pathCollisions // crossSecFile
CALL interactionMatrix(ij)%collisions(k)%obj%init(crossSecFilePath)
END DO
END DO
END SUBROUTINE readInteractions
!Reads boundary conditions for the mesh
SUBROUTINE readBoundary(config)
USE moduleBoundary
USE moduleErrors
USE json_module
IMPLICIT NONE
TYPE(json_file), INTENT(inout):: config
CHARACTER(2):: istring
CHARACTER(:), ALLOCATABLE:: object
LOGICAL:: found
INTEGER:: i
CALL config%info('boundary', n_children = nBoundary)
ALLOCATE(boundary(1:nBoundary))
DO i = 1, nBoundary
WRITE(istring, '(i2)') i
object = 'boundary(' // trim(istring) // ')'
ALLOCATE(boundaryGeneric:: boundary(i)%obj)
CALL config%get(object // '.type', boundary(i)%obj%boundaryType, found)
CALL config%get(object // '.physicalSurface', boundary(i)%obj%physicalSurface, found)
boundary(i)%obj%id = i
END DO
END SUBROUTINE readBoundary
!Read the geometry (mesh) for the case
SUBROUTINE readGeometry(config)
USE moduleMesh
USE moduleMeshCylRead
USE moduleErrors
USE moduleOutput
USE json_module
IMPLICIT NONE
TYPE(json_file), INTENT(inout):: config
LOGICAL:: found
CHARACTER(:), ALLOCATABLE:: geometryType, meshType, meshFile
CHARACTER(:), ALLOCATABLE:: fullPath
!Selects the type of geometry.
CALL config%get('geometry.type', geometryType, found)
SELECT CASE(geometryType)
CASE ("2DCyl")
!Creates a 2D cylindrical mesh
ALLOCATE(meshCylGeneric:: mesh)
!Gets the type of mesh
CALL config%get('geometry.meshType', meshType, found)
SELECT CASE(meshType)
CASE ("gmsh")
!Gets the gmsh file
CALL config%get('geometry.meshFile', meshFile, found)
CASE DEFAULT
CALL criticalError("Mesh type " // meshType // " not supported.", "readGeometry")
END SELECT
!Reads the mesh
fullpath = path // meshFile
CALL mesh%readMesh(fullPath)
CASE DEFAULT
CALL criticalError("Geometry type " // geometryType // " not supported.", "readGeometry")
END SELECT
END SUBROUTINE readGeometry
!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
INTEGER:: physicalSurface
INTEGER:: pt
CALL config%info('inject', 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)
pt = 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 // '.physicalSurface', physicalSurface, found)
CALL inject(i)%init(i, v, normal, T, flow, pt, physicalSurface)
END DO
!Allocate array for injected particles
IF (nPartInj > 0) THEN
ALLOCATE(part_inj(1:nPartInj))
END IF
END SUBROUTINE readInject
END MODULE moduleInput

79
src/modules/moduleList.f95 Normal file
View file

@ -0,0 +1,79 @@
!Linked list of particles
MODULE moduleList
IMPLICIT NONE
TYPE lNode
INTEGER:: n = 0
TYPE(lNode), POINTER:: next => NULL()
END TYPE lNode
TYPE listNode
INTEGER:: amount = 0
TYPE(lNode),POINTER:: head => NULL()
TYPE(lNode),POINTER:: tail => NULL()
CONTAINS
PROCEDURE,PASS:: add => addToList
PROCEDURE,PASS:: get => getFromList
PROCEDURE,PASS:: erase => eraseList
END TYPE listNode
CONTAINS
!Adds element to list
SUBROUTINE addToList(this,n)
INTEGER,INTENT(in):: n
CLASS(listNode):: this
TYPE(lNode),POINTER:: temp
ALLOCATE(temp)
temp%n = n
NULLIFY(temp%next)
IF (.NOT. ASSOCIATED(this%head)) THEN
!First element
this%head => temp
this%tail => temp
this%amount = 1
ELSE
!Append element
this%tail%next => temp
this%tail => temp
this%amount = this%amount + 1
END IF
END SUBROUTINE addToList
FUNCTION getFromList(self, iObj) RESULT(n)
use moduleErrors
IMPLICIT NONE
CLASS(listNode):: self
INTEGER:: iObj
INTEGER:: n
INTEGER:: i
TYPE(lNode), POINTER:: tempNode
IF (iObj > self%amount) CALL criticalError('Accessing to element list outisde range','getFromList')
tempNode => self%head
DO i = 1, iObj - 1
tempNode => tempNode%next
END DO
n = tempNode%n
END FUNCTION getFromList
!Erase list
SUBROUTINE eraseList(this)
CLASS(listNode):: this
TYPE(lNode),POINTER:: current, next
current => this%head
DO WHILE (ASSOCIATED(current))
next => current%next
DEALLOCATE(current)
current => next
END DO
IF (ASSOCIATED(this%head)) NULLIFY(this%head)
IF (ASSOCIATED(this%tail)) NULLIFY(this%tail)
this%amount = 0
END SUBROUTINE eraseList
END MODULE moduleList

198
src/modules/moduleMesh.f95 Normal file
View file

@ -0,0 +1,198 @@
!moduleMesh: General module for Finite Element mesh
MODULE moduleMesh
USE moduleList
USE moduleOutput
IMPLICIT NONE
!Parent of Node element
TYPE, PUBLIC, ABSTRACT:: meshNode
!Node index
INTEGER:: n = 0
!Node volume
REAL(8):: v = 0.D0
!Output values
TYPE(outputNode), ALLOCATABLE:: output(:)
CONTAINS
PROCEDURE(initNode_interface), DEFERRED, PASS:: init
PROCEDURE(getCoord_interface), DEFERRED, PASS:: getCoordinates
END TYPE meshNode
ABSTRACT INTERFACE
!Interface of init a node (3D generic coordinates)
SUBROUTINE initNode_interface(self, n, r)
IMPORT:: meshNode
CLASS(meshNode), INTENT(out):: self
INTEGER, INTENT(in):: n
REAL(8), INTENT(in):: r(1:3)
END SUBROUTINE initNode_interface
!Interface to get coordinates from node
FUNCTION getCoord_interface(self) RESULT(r)
IMPORT:: meshNode
CLASS(meshNode):: self
REAL(8):: r(1:3)
END FUNCTION getCoord_interface
END INTERFACE
!Containers for nodes in the mesh
TYPE:: meshNodeCont
CLASS(meshNode), ALLOCATABLE:: obj
END TYPE meshNodeCont
!Parent of Edge element
TYPE, PUBLIC, ABSTRACT:: meshEdge
!Element index
INTEGER:: n = 0
!Connectivity to vols
CLASS(meshVol), POINTER:: e1 => NULL(), e2 => NULL()
!Normal vector
REAL(8):: normal(1:3)
!Physical surface in mesh
INTEGER:: physicalSurface
!id for boundary condition
INTEGER:: bt = 0
CONTAINS
PROCEDURE(initEdge_interface), DEFERRED, PASS:: init
END TYPE meshEdge
ABSTRACT INTERFACE
SUBROUTINE initEdge_interface(self, n, p, bt, physicalSurface)
IMPORT:: meshEdge
CLASS(meshEdge), INTENT(out):: self
INTEGER, INTENT(in):: n
INTEGER, INTENT(in):: p(:)
INTEGER, INTENT(in):: bt
INTEGER, INTENT(in):: physicalSurface
END SUBROUTINE initEdge_interface
END INTERFACE
!Containers for edges in the mesh
TYPE:: meshEdgeCont
CLASS(meshEdge), ALLOCATABLE:: obj
END TYPE meshEdgeCont
!Parent of Volume element
TYPE, PUBLIC, ABSTRACT:: meshVol
!Volume index
INTEGER:: n = 0
!Maximum collision rate
REAL(8):: sigmaVrelMax = 1.D-15
!Volume
REAL(8):: volume = 0.D0
!List of particles inside the volume
TYPE(listNode):: listPart_in
!Number of collisions per volume
INTEGER:: nColl = 0
CONTAINS
PROCEDURE(initVol_interface), DEFERRED, PASS:: init
PROCEDURE(scatter_interface), DEFERRED, PASS:: scatter
PROCEDURE(collision_interface), DEFERRED, PASS:: collision
PROCEDURE(findCell_interface), DEFERRED, PASS:: findCell
END TYPE meshVol
ABSTRACT INTERFACE
SUBROUTINE initVol_interface(self, n, p)
IMPORT:: meshVol
CLASS(meshVol), INTENT(out):: self
INTEGER, INTENT(in):: n
INTEGER, INTENT(in):: p(:)
END SUBROUTINE initVol_interface
SUBROUTINE scatter_interface(self, part)
USE moduleSpecies
IMPORT:: meshVol
CLASS(meshVol), INTENT(in):: self
CLASS(particle), INTENT(in):: part
END SUBROUTINE scatter_interface
SUBROUTINE collision_interface(self)
IMPORT:: meshVol
CLASS(meshVol), INTENT(inout):: self
END SUBROUTINE collision_interface
SUBROUTINE findCell_interface(self, part)
USE moduleSpecies
IMPORT:: meshVol
CLASS(meshVol), INTENT(in):: self
CLASS(particle), INTENT(inout):: part
END SUBROUTINE findCell_interface
END INTERFACE
!Containers for volumes in the mesh
TYPE:: meshVolCont
CLASS(meshVol), ALLOCATABLE:: obj
END TYPE meshVolCont
!Abstract type of mesh
TYPE, PUBLIC, ABSTRACT:: meshGeneric
INTEGER:: numEdges, numNodes, numVols
!Array of nodes
TYPE(meshNodeCont), ALLOCATABLE:: nodes(:)
!Array of boundary elements
TYPE(meshEdgeCont), ALLOCATABLE:: edges(:)
!Array of volume elements
TYPE(meshVolCont), ALLOCATABLE:: vols(:)
!Global stiffness matrix
REAL(8), ALLOCATABLE, DIMENSION(:,:):: K
!Global load vector
REAL(8), ALLOCATABLE, DIMENSION(:):: F
CONTAINS
PROCEDURE(readMesh_interface), PASS, DEFERRED:: readMesh
PROCEDURE(printOutput_interface), PASS, DEFERRED:: printOutput
PROCEDURE(printColl_interface), PASS, DEFERRED:: printColl
END TYPE meshGeneric
ABSTRACT INTERFACE
!Reads the mesh from a file
SUBROUTINE readMesh_interface(self, filename)
IMPORT meshGeneric
CLASS(meshGeneric), INTENT(out):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
END SUBROUTINE readMesh_interface
!Prints output variables
SUBROUTINE printOutput_interface(self, t)
IMPORT meshGeneric
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
END SUBROUTINE printOutput_interface
!Prints unmber of collisions
SUBROUTINE printColl_interface(self, t)
IMPORT meshGeneric
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
END SUBROUTINE printColl_interface
END INTERFACE
!Generic mesh
CLASS(meshGeneric), ALLOCATABLE, TARGET:: mesh
END MODULE moduleMesh

534
src/modules/moduleMeshCyl.f95 Normal file
View file

@ -0,0 +1,534 @@
!moduleMeshCyl: 2D axial symmetric extension of generic mesh from GMSH format.
! x == z
! y == r
! z == theta (unused)
MODULE moduleMeshCyl
USE moduleMesh
IMPLICIT NONE
!TODO: Move this, this is horrible
REAL(8), PARAMETER:: w(1:3) = (/ 5.D0/9.D0, 8.D0/9.D0, 5.D0/9.D0 /)
REAL(8), PARAMETER:: cor(1:3) = (/ -DSQRT(3.D0/5.D0), 0.D0, DSQRT(3.D0/5.D0) /)
TYPE, PUBLIC, EXTENDS(meshNode):: meshNodeCyl
!Element coordinates
REAL(8):: r = 0.D0, z = 0.D0
CONTAINS
PROCEDURE, PASS:: init => initNodeCyl
PROCEDURE, PASS:: getCoordinates => getCoordCyl
END TYPE meshNodeCyl
TYPE, PUBLIC, ABSTRACT, EXTENDS(meshEdge):: meshEdgeCyl
!Element coordinates
REAL(8):: r(1:2) = 0.D0, z(1:2) = 0.D0
!Connectivity to nodes
CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL()
CONTAINS
PROCEDURE, PASS:: init => initEdgeCyl
PROCEDURE(boundary_interface), PASS, DEFERRED:: fBoundary
END TYPE meshEdgeCyl
ABSTRACT INTERFACE
SUBROUTINE boundary_interface(self, part)
USE moduleSpecies
IMPORT:: meshEdgeCyl
CLASS (meshEdgeCyl), INTENT(inout):: self
CLASS (particle), INTENT(inout):: part
END SUBROUTINE
END INTERFACE
TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVolCyl
CONTAINS
PROCEDURE, PASS:: collision => collision2DCyl
END TYPE meshVolCyl
TYPE, PUBLIC, EXTENDS(meshVolCyl):: meshVolCylQuad
!Element coordinates
REAL(8):: r(1:4) = 0.D0, z(1:4) = 0.D0
!Connectivity to nodes
CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL(), n4 => NULL()
!Connectivity to adjacent elements
CLASS(*), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL(), e4 => NULL()
!Maximum collision rate
!TODO: Move to other more appropiate section
REAL(8):: phi(1:4) = 0.D0
REAL(8):: Ke(1:4,1:4) = 0.D0
REAL(8):: arNodes(1:4) = 0.D0
CONTAINS
PROCEDURE, PASS:: init => initVolQuadCyl
PROCEDURE, PASS:: locKe => localKeQuad
PROCEDURE, PASS:: detJac => detJQuad
PROCEDURE, PASS:: invJ => invJQuad
PROCEDURE, PASS:: area => areaQuad
PROCEDURE, NOPASS:: weight => weightQuad
PROCEDURE, NOPASS:: inside => insideQuad
PROCEDURE, PASS:: dVal => dValQuad
PROCEDURE, PASS:: scatter => scatterQuad
PROCEDURE, PASS:: phy2log => phy2logQuad
PROCEDURE, PASS:: findCell => findCellCylQuad
END TYPE meshVolCylQuad
CONTAINS
!Inits node element
SUBROUTINE initNodeCyl(self, n, r)
USE moduleSpecies
USE moduleRefParam
IMPLICIT NONE
CLASS(meshNodeCyl), INTENT(out):: self
INTEGER, INTENT(in):: n
REAL(8), INTENT(in):: r(1:3)
self%n = n
self%r = r(1)/L_ref
self%z = r(2)/L_ref
!Node volume, to be determined in mesh
self%v = 0.D0
!Allocates output:
ALLOCATE(self%output(1:nSpecies))
END SUBROUTINE initNodeCyl
FUNCTION getCoordCyl(self) RESULT(r)
IMPLICIT NONE
CLASS(meshNodeCyl):: self
REAL(8):: r(1:3)
r = (/self%z, self%r, 0.D0/)
END FUNCTION getCoordCyl
!Edge functions
SUBROUTINE initEdgeCyl(self, n, p, bt, physicalSurface)
IMPLICIT NONE
CLASS(meshEdgeCyl), INTENT(out):: self
INTEGER, INTENT(in):: n
INTEGER, INTENT(in):: p(:)
INTEGER, INTENT(in):: bt
INTEGER, INTENT(in):: physicalSurface
REAL(8):: r1(1:3), r2(1:3)
self%n = n
self%n1 => mesh%nodes(p(1))%obj
self%n2 => mesh%nodes(p(2))%obj
!Get element coordinates
r1 = self%n1%getCoordinates()
r2 = self%n2%getCoordinates()
self%z = (/r1(1), r2(1)/)
self%r = (/r1(2), r2(2)/)
!Normal vector
self%normal = (/ self%r(2)-self%r(1), &
self%z(2)-self%z(1), &
0.D0 /)
!Boundary index
self%bt = bt
!Phyiscal Surface
self%physicalSurface = physicalSurface
END SUBROUTINE initEdgeCyl
!Vol functions
!Quad Volume
SUBROUTINE initVolQuadCyl(self, n, p)
USE moduleRefParam
IMPLICIT NONE
CLASS(meshVolCylQuad), INTENT(out):: self
INTEGER, INTENT(in):: n
INTEGER, INTENT(in):: p(:)
REAL(8):: r1(1:3), r2(1:3), r3(1:3), r4(1:3)
self%n = n
self%n1 => mesh%nodes(p(1))%obj
self%n2 => mesh%nodes(p(2))%obj
self%n3 => mesh%nodes(p(3))%obj
self%n4 => mesh%nodes(p(4))%obj
!Get element coordinates
r1 = self%n1%getCoordinates()
r2 = self%n2%getCoordinates()
r3 = self%n3%getCoordinates()
r4 = self%n4%getCoordinates()
self%z = (/r1(1), r2(1), r3(1), r4(1)/)
self%r = (/r1(2), r2(2), r3(2), r4(2)/)
!Assign node volume
CALL self%area()
self%n1%v = self%n1%v + self%arNodes(1)
self%n2%v = self%n2%v + self%arNodes(2)
self%n3%v = self%n3%v + self%arNodes(3)
self%n4%v = self%n4%v + self%arNodes(4)
CALL self%locKe()
self%sigmaVrelMax = sigma_ref/L_ref**2
END SUBROUTINE initVolQuadCyl
FUNCTION fpsi(xi1,xi2) RESULT(psi)
IMPLICIT NONE
REAL(8),INTENT(in):: xi1,xi2
REAL(8):: psi(1:4)
psi(1) = (1.D0-xi1)*(1.D0-xi2)
psi(2) = (1.D0+xi1)*(1.D0-xi2)
psi(3) = (1.D0+xi1)*(1.D0+xi2)
psi(4) = (1.D0-xi1)*(1.D0+xi2)
psi = psi*0.25D0
END FUNCTION fpsi
!Derivative element function (xi1)
FUNCTION dpsiXi1(xi2) RESULT(psiXi1)
IMPLICIT NONE
REAL(8),INTENT(in):: xi2
REAL(8):: psiXi1(1:4)
psiXi1(1) = -(1.D0-xi2)
psiXi1(2) = (1.D0-xi2)
psiXi1(3) = (1.D0+xi2)
psiXi1(4) = -(1.D0+xi2)
psiXi1 = psiXi1*0.25D0
END FUNCTION dpsiXi1
!Derivative element function (xi2)
FUNCTION dpsiXi2(xi1) RESULT(psiXi2)
IMPLICIT NONE
REAL(8),INTENT(in):: xi1
REAL(8):: psiXi2(1:4)
psiXi2(1) = -(1.D0-xi1)
psiXi2(2) = -(1.D0+xi1)
psiXi2(3) = (1.D0+xi1)
psiXi2(4) = (1.D0-xi1)
psiXi2 = psiXi2*0.25D0
END FUNCTION dpsiXi2
!Transforms physical coordinates to element coordinates
FUNCTION phy2logQuad(this,r) RESULT(xN)
IMPLICIT NONE
CLASS(meshVolCylQuad):: this
REAL(8):: r(1:2)
REAL(8):: xN(1:2), xO(1:2), dPsi(1:2,1:4), detJ, j(1:2,1:2), psi(1:4), f(1:2)
REAL(8):: conv
!Iterative newton method to transform coordinates
conv=1.D0
xO=0.D0
DO WHILE(conv>1.D-4)
dPsi(1,:) = dpsiXi1(xO(2))
dPsi(2,:) = dpsiXi2(xO(1))
j = this%invJ(xO(1),xO(2),dPsi)
psi = fpsi(xO(1), xO(2))
f(1) = DOT_PRODUCT(psi,this%z)-r(1)
f(2) = DOT_PRODUCT(psi,this%r)-r(2)
detJ = this%detJac(xO(1),xO(2),dPsi)
xN=xO - MATMUL(j, f)/detJ
conv=MAXVAL(DABS(xN-xO),1)
xO=xN
END DO
END FUNCTION phy2logQuad
!Computes element local stiffness matrix
SUBROUTINE localKeQuad(this)
USE moduleConstParam
IMPLICIT NONE
CLASS(meshVolCylQuad):: this
REAL(8):: r, xi1, xi2, dPsi(1:2,1:4)
REAL(8):: j(1:2,1:2), detJ
INTEGER:: l, m
this%Ke=0.D0
!Start 2D Gauss Quad Integral
DO l=1, 3
DO m = 1, 3
xi1 = cor(l)
xi2 = cor(m)
dPsi(1,:) = dpsiXi1(xi2)
dPsi(2,:) = dpsiXi2(xi1)
detJ = this%detJac(xi1,xi2,dPsi)
j = this%invJ(xi1,xi2,dPsi)
r = DOT_PRODUCT(fpsi(xi1,xi2),this%r)
this%Ke = this%Ke + MATMUL(TRANSPOSE(MATMUL(j,dPsi)),MATMUL(j,dPsi))*r*w(l)*w(m)/detJ
END DO
END DO
this%Ke = this%Ke*2.D0*PI
END SUBROUTINE localKeQuad
!Computes element Jacobian determinant
FUNCTION detJQuad(this,xi1,xi2,dPsi_in) RESULT(dJ)
IMPLICIT NONE
REAL(8),OPTIONAL:: dPsi_in(1:2,1:4)
REAL(8):: dPsi(1:2,1:4)
REAL(8):: dz(1:2), dr(1:2)
REAL(8):: xi1,xi2, dJ
CLASS(meshVolCylQuad):: this
IF(PRESENT(dPsi_in)) THEN
dPsi=dPsi_in
ELSE
dPsi(1,:) = dpsiXi1(xi2)
dPsi(2,:) = dpsiXi2(xi1)
END IF
dz(1) = DOT_PRODUCT(dPsi(1,:),this%z)
dz(2) = DOT_PRODUCT(dPsi(2,:),this%z)
dr(1) = DOT_PRODUCT(dPsi(1,:),this%r)
dr(2) = DOT_PRODUCT(dPsi(2,:),this%r)
dJ = dz(1)*dr(2)-dz(2)*dr(1)
END FUNCTION detJQuad
FUNCTION invJQuad(this,xi1,xi2,dPsi_in) RESULT(j) !Computes element Jacobian inverse matrix (without determinant)
IMPLICIT NONE
REAL(8),OPTIONAL:: dpsi_in(1:2,1:4)
REAL(8):: dPsi(1:2,1:4)
REAL(8):: dz(1:2), dr(1:2)
REAL(8):: xi1,xi2, j(1:2,1:2)
CLASS(meshVolCylQuad):: this
IF(PRESENT(dPsi_in)) THEN
dPsi=dPsi_in
ELSE
dPsi(1,:) = dpsiXi1(xi2)
dPsi(2,:) = dpsiXi2(xi1)
END IF
dz(1) = DOT_PRODUCT(dPsi(1,:),this%z)
dz(2) = DOT_PRODUCT(dPsi(2,:),this%z)
dr(1) = DOT_PRODUCT(dPsi(1,:),this%r)
dr(2) = DOT_PRODUCT(dPsi(2,:),this%r)
j(1,:) = (/ dr(2), -dz(2) /)
j(2,:) = (/ -dr(1), dz(1) /)
END FUNCTION invJQuad
SUBROUTINE areaQuad(this)!Computes element area
USE moduleConstParam
IMPLICIT NONE
CLASS(meshVolCylQuad):: this
REAL(8):: r, xi1, xi2
REAL(8):: detJ, fpsi0(1:4)
this%volume = 0.D0
this%arNodes = 0.D0
!2D 1 point Gauss Quad Integral
xi1 = 0.D0
xi2 = 0.D0
detJ = this%detJac(xi1,xi2)*8.D0*PI !4*2*pi
fpsi0 = fpsi(xi1,xi2)
r = DOT_PRODUCT(fpsi0,this%r)
this%volume = r*detJ
this%arNodes = fpsi0*r*detJ
END SUBROUTINE areaQuad
FUNCTION weightQuad(xi1_p, xi2_p) RESULT(w) !Computes weights in the four vertices. '_p' references particle position in the logical space
IMPLICIT NONE
REAL(8):: xi1_p, xi2_p
REAL(8):: w(1:4)
w = fpsi(xi1_p, xi2_p)
END FUNCTION weightQuad
FUNCTION insideQuad(xi1_p, xi2_p) RESULT(ins) !Checks if a particle is inside a quad element
IMPLICIT NONE
REAL(8):: xi1_p, xi2_p
LOGICAL:: ins
ins = (xi1_p >= -1.D0 .AND. xi1_p <= 1.D0) .AND. &
(xi2_p >= -1.D0 .AND. xi2_p <= 1.D0)
END FUNCTION insideQuad
FUNCTION dValQuad(this,xi1,xi2) RESULT(EF)!Computes electric field in xi1,xi2
IMPLICIT NONE
CLASS(meshVolCylQuad):: this
REAL(8):: xi1, xi2, dPsi(1:2,1:4)
REAL(8):: j(1:2,1:2), detJ
REAL(8):: EF(1:3)
dPsi(1,:) = dpsiXi1(xi2)
dPsi(2,:) = dpsiXi2(xi1)
detJ = this%detJac(xi1,xi2,dPsi)
j = this%invJ(xi1,xi2,dPsi)
EF(1) = -DOT_PRODUCT(dPsi(1,:), this%phi)*j(2,2)/detJ
EF(2) = -DOT_PRODUCT(dPsi(2,:), this%phi)*j(1,1)/detJ
EF(3) = 0.D0
END FUNCTION dValQuad
SUBROUTINE scatterQuad(self, part)
USE moduleOutput
USE moduleSpecies
IMPLICIT NONE
CLASS(meshVolCylQuad), INTENT(in):: self
CLASS(particle), INTENT(in):: part
TYPE(outputNode), POINTER:: vertex
REAL(8):: w_p(1:4)
REAL(8):: tensorS(1:3,1:3)
w_p = self%weight(part%xLog(1), part%xLog(2))
tensorS = outerProduct(part%v, part%v)
vertex => self%n1%output(part%pt)
vertex%den = vertex%den + w_p(1)
vertex%mom(:) = vertex%mom(:) + w_p(1)*part%v(:)
vertex%tensorS(:,:) = vertex%tensorS(:,:) + w_p(1)*tensorS
vertex => self%n2%output(part%pt)
vertex%den = vertex%den + w_p(2)
vertex%mom(:) = vertex%mom(:) + w_p(2)*part%v(:)
vertex%tensorS(:,:) = vertex%tensorS(:,:) + w_p(2)*tensorS
vertex => self%n3%output(part%pt)
vertex%den = vertex%den + w_p(3)
vertex%mom(:) = vertex%mom(:) + w_p(3)*part%v(:)
vertex%tensorS(:,:) = vertex%tensorS(:,:) + w_p(3)*tensorS
vertex => self%n4%output(part%pt)
vertex%den = vertex%den + w_p(4)
vertex%mom(:) = vertex%mom(:) + w_p(4)*part%v(:)
vertex%tensorS(:,:) = vertex%tensorS(:,:) + w_p(4)*tensorS
END SUBROUTINE scatterQuad
RECURSIVE SUBROUTINE findCellCylQuad(self, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshVolCylQuad), INTENT(in):: self
CLASS(particle), INTENT(inout):: part
REAL(8):: xLog(1:2)
REAL(8):: xLogArray(1:4)
CLASS(*), POINTER:: nextElement
INTEGER:: nextInt
xLog = self%phy2log(part%r(1:2))
IF (self%inside(xLog(1), xLog(2))) THEN
!Checks if particle is inside of current cell
part%e_p = self%n
part%xLog = xLog
ELSE
!If not, searches for a neighbour and repeats the process.
xLogArray = (/ -xLog(2), xLog(1), xLog(2), -xLog(1) /)
nextInt = MAXLOC(xLogArray,1)
!Selects the higher value of directions and searches in that direction
NULLIFY(nextElement)
SELECT CASE (nextInt)
CASE (1)
nextElement => self%e1
CASE (2)
nextElement => self%e2
CASE (3)
nextElement => self%e3
CASE (4)
nextElement => self%e4
END SELECT
!Defines the next step
SELECT TYPE(nextElement)
CLASS IS(meshVolCyl)
!Particle moved to new cell, repeat find procedure
CALL nextElement%findCell(part)
CLASS IS (meshEdgeCyl)
!Particle encountered an edge, execute boundary
CALL nextElement%fBoundary(part)
CLASS DEFAULT
WRITE(*,*) "ERROR, CHECK findCellCylQuad"
END SELECT
END IF
END SUBROUTINE findCellCylQuad
SUBROUTINE collision2DCyl(self)
USE moduleRefParam
USE moduleConstParam
USE moduleCollisions
USE moduleSpecies
USE moduleList
IMPLICIT NONE
CLASS(meshVolCyl), INTENT(inout):: self
INTEGER:: Npart !Number of particles inside the cell
INTEGER:: Ncoll !Maximum number of collisions
REAL(8):: Fn !Specific weight
REAL(8):: Pmax !Maximum probability of collision
INTEGER:: i,j !random particles index
INTEGER:: rnd !random index
TYPE(particle), POINTER:: part_i, part_j
INTEGER:: n !collision
INTEGER:: ij, k
REAL(8):: sigmaVrel
REAL(8):: sigmaVrelMaxNew
Fn = species(1)%obj%weight!Check how to do this for multiple species
Npart = self%listPart_in%amount
Pmax = Fn*self%sigmaVrelMax*tau/self%volume
Ncoll = INT(REAL(Npart*(Npart-1))*Pmax*0.5D0)
self%nColl = Ncoll
DO n = 1, NColl
!Select random numbers
rnd = 1 + FLOOR(Npart*RAND())
i = self%listPart_in%get(rnd)
part_i => part_old(i)
rnd = 1 + FLOOR(Npart*RAND())
j = self%listPart_in%get(rnd)
part_j => part_old(j)
ij = interactionIndex(part_i%pt, part_j%pt)
sigmaVrelMaxNew = 0.D0
DO k = 1, interactionMatrix(ij)%amount
CALL interactionMatrix(ij)%collisions(k)%obj%collide(self%sigmaVrelMax, sigmaVrel, part_i, part_j)
sigmaVrelMaxNew = sigmaVrel + SigmaVrelMaxNew
END DO
!Update maximum cross section times vrelative per each collision
IF (sigmaVrelMaxNew > self%sigmaVrelMax) self%sigmaVrelMax = sigmaVrelMaxNew
END DO
CALL self%listPart_in%erase()
END SUBROUTINE collision2DCyl
END MODULE moduleMeshCyl

80
src/modules/moduleMeshCylBoundary.f95 Normal file
View file

@ -0,0 +1,80 @@
!moduleMeshCylBoundary: Edge elements for Cylindrical mesh.
MODULE moduleMeshCylBoundary
USE moduleMeshCyl
TYPE, PUBLIC, EXTENDS(meshEdgeCyl):: meshEdgeCylRef
CONTAINS
PROCEDURE, PASS:: fBoundary => reflection
END TYPE meshEdgeCylRef
TYPE, PUBLIC, EXTENDS(meshEdgeCyl):: meshEdgeCylAbs
CONTAINS
PROCEDURE, PASS:: fBoundary => absorption
END TYPE meshEdgeCylAbs
TYPE, PUBLIC, EXTENDS(meshEdgeCyl):: meshEdgeCylAxis
CONTAINS
PROCEDURE, PASS:: fBoundary => symmetryAxis
END TYPE meshEdgeCylAxis
CONTAINS
SUBROUTINE reflection(self, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdgeCylRef), INTENT(inout):: self
CLASS(particle), INTENT(inout):: part
REAL(8):: edgeNorm, cosT, sinT, rp(1:2), rpp(1:2), vpp(1:2)
edgeNorm = DSQRT((self%r(2)-self%r(1))**2 + (self%z(2)-self%z(1))**2)
cosT = (self%z(2)-self%z(1))/edgeNorm
sinT = DSQRT(1-cosT**2)
rp(1) = part%r(1) - self%z(1);
rp(2) = part%r(2) - self%r(1);
rpp(1) = cosT*rp(1) - sinT*rp(2)
rpp(2) = sinT*rp(1) + cosT*rp(2)
rpp(2) = -rpp(2)
vpp(1) = cosT*part%v(1) - sinT*part%v(2)
vpp(2) = sinT*part%v(1) + cosT*part%v(2)
vpp(2) = -vpp(2)
part%r(1) = cosT*rpp(1) + sinT*rpp(2) + self%z(1);
part%r(2) = -sinT*rpp(1) + cosT*rpp(2) + self%r(1);
part%v(1) = cosT*vpp(1) + sinT*vpp(2)
part%v(2) = -sinT*vpp(1) + cosT*vpp(2)
END SUBROUTINE reflection
!Absoption in a surface
SUBROUTINE absorption(self, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdgeCylAbs), INTENT(inout):: self
CLASS(particle), INTENT(inout):: part
!TODO: Add scatter to mesh nodes
part%n_in = .FALSE.
END SUBROUTINE absorption
SUBROUTINE symmetryAxis(self, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdgeCylAxis), INTENT(inout):: self
CLASS(particle), INTENT(inout):: part
!Do to the 2D Cyl pusher, this function should never be called.
END SUBROUTINE symmetryAxis
END MODULE moduleMeshCylBoundary

428
src/modules/moduleMeshCylRead.f95 Normal file
View file

@ -0,0 +1,428 @@
MODULE moduleMeshCylRead
USE moduleMesh
USE moduleMeshCyl
USE moduleMeshCylBoundary
TYPE, EXTENDS(meshGeneric):: meshCylGeneric
CONTAINS
PROCEDURE, PASS:: readMesh => readMeshCyl
PROCEDURE, PASS:: printOutput => printOutputCyl
PROCEDURE, PASS:: printColl => printCollisionsCyl
END TYPE
INTERFACE connected
MODULE PROCEDURE connectedVolVol, connectedVolEdge
END INTERFACE connected
CONTAINS
SUBROUTINE readMeshCyl(self, filename)
USE moduleRefParam
USE moduleBoundary
IMPLICIT NONE
CLASS(meshCylGeneric), INTENT(out):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
REAL(8):: r, z
INTEGER:: p(1:4)
INTEGER:: e=0, et=0, n=0, eTemp=0, elemType=0, bt = 0
INTEGER:: totalNumElem
INTEGER:: boundaryType
!Read msh
OPEN(10, FILE=TRIM(filename))
!Skip header
READ(10, *)
READ(10, *)
READ(10, *)
READ(10, *)
!Read number of nodes
READ(10, *) self%numNodes
!Allocate required matrices and vectors
ALLOCATE(self%nodes(1:self%numNodes))
! ALLOCATE(self%K(1:numNodes,1:numNodes))
! ALLOCATE(self%F(1:numNodes))
!Read nodes cartesian coordinates (x=z, y=r, z=null)
DO e=1, self%numNodes
READ(10, *) n, z, r
ALLOCATE(meshNodeCyl:: self%nodes(n)%obj)
CALL self%nodes(n)%obj%init(n, (/r, z, 0.D0 /))
END DO
!Skips comments
READ(10, *)
READ(10, *)
!Reads Totalnumber of elements
READ(10, *) TotalnumElem
!counts edges and volume elements
self%numEdges = 0
DO e=1, TotalnumElem
READ(10,*) eTemp, elemType
IF (elemType==1) THEN
self%numEdges=e
END IF
END DO
!Substract the number of edges to the total number of elements
!to obtain the number of volume elements
self%numVols = TotalnumElem - self%numEdges
!Allocates arrays
ALLOCATE(self%edges(1:self%numEdges))
ALLOCATE(self%vols(1:self%numVols))
!Go back to the beggining to read elements
DO e=1, totalNumElem
BACKSPACE(10)
END DO
!TODO: symplify to read all elements independently if they are edges or vols
!Reads edges
DO e=1, self%numEdges
READ(10,*) n, elemType, eTemp, boundaryType, eTemp, p(1:2)
!Associate boundary condition procedure.
!TODO: move to subroutine
bt = getBoundaryId(boundaryType)
SELECT CASE(boundary(bt)%obj%boundaryType)
CASE ('reflection')
ALLOCATE(meshEdgeCylRef:: self%edges(e)%obj)
CASE ('absorption')
ALLOCATE(meshEdgeCylAbs:: self%edges(e)%obj)
CASE ('axis')
ALLOCATE(meshEdgeCylAxis:: self%edges(e)%obj)
END SELECT
CALL self%edges(e)%obj%init(n, p, bt, boundaryType)
END DO
!Read and initialize volumes
!TODO: Extend to triangular elements
DO e=1, self%numVols
READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:4)
ALLOCATE(meshVolCylQuad:: self%vols(e)%obj)
CALL self%vols(e)%obj%init(n - self%numEdges, p)
END DO
CLOSE(10)
!Build connectivity between elements
DO e = 1, self%numVols
!Connectivity between volumes
DO et = 1, self%numVols
CALL connected(self%vols(e)%obj, self%vols(et)%obj)
END DO
!Connectivity between vols and edges
DO et = 1, self%numEdges
CALL connected(self%vols(e)%obj, self%edges(et)%obj)
END DO
END DO
! !Compute global stiffness matrix
! GlobalK=0.D0
! DO e=1, numElem
! DO i=1, 4
! DO j=1, 4
! GlobalK(elems(e)%p(i),elems(e)%p(j)) = GlobalK(elems(e)%p(i),elems(e)%p(j)) + elems(e)%Ke(i,j)
! END DO
! END DO
! END DO
! ! Apply Dirichlet boundary conditions to GlobalK
! DO n=1, numNodes
! IF (nodes(n)%bound == 1) THEN
! GlobalK(n,:)=0.D0
! GlobalK(n,n)=1.D0
! END IF
! END DO
END SUBROUTINE
SUBROUTINE connectedVolVol(elemA, elemB)
IMPLICIT NONE
CLASS(meshVol), INTENT(inout):: elemA
CLASS(meshVol), INTENT(inout):: elemB
SELECT TYPE(elemA)
TYPE IS (meshVolCylQuad)
SELECT TYPE(elemB)
TYPE IS(meshVolCylQuad)
CALL connectedQuadQuad(elemA, elemB)
END SELECT
END SELECT
END SUBROUTINE connectedVolVol
SUBROUTINE connectedVolEdge(elemA, elemB)
IMPLICIT NONE
CLASS(meshVol), INTENT(inout):: elemA
CLASS(meshEdge), INTENT(inout):: elemB
SELECT TYPE(elemA)
TYPE IS (meshVolCylQuad)
SELECT TYPE(elemB)
CLASS IS(meshEdgeCyl)
CALL connectedQuadEdge(elemA, elemB)
END SELECT
END SELECT
END SUBROUTINE connectedVolEdge
SUBROUTINE connectedQuadQuad(elemA, elemB)
IMPLICIT NONE
CLASS(meshVolCylQuad), INTENT(inout), TARGET:: elemA
CLASS(meshVolCylQuad), INTENT(inout), TARGET:: elemB
!Check direction 1
IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
elemA%n1%n == elemB%n4%n .AND. &
elemA%n2%n == elemB%n3%n) THEN
elemA%e1 => elemB
elemB%e3 => elemA
END IF
!Check direction 2
IF (.NOT. ASSOCIATED(elemA%e2) .AND. &
elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n4%n) THEN
elemA%e2 => elemB
elemB%e4 => elemA
END IF
!Check direction 3
IF (.NOT. ASSOCIATED(elemA%e3) .AND. &
elemA%n3%n == elemB%n2%n .AND. &
elemA%n4%n == elemB%n1%n) THEN
elemA%e3 => elemB
elemB%e1 => elemA
END IF
!Check direction 4
IF (.NOT. ASSOCIATED(elemA%e4) .AND. &
elemA%n4%n == elemB%n3%n .AND. &
elemA%n1%n == elemB%n2%n) THEN
elemA%e4 => elemB
elemB%e2 => elemA
END IF
END SUBROUTINE connectedQuadQuad
SUBROUTINE connectedQuadEdge(elemA, elemB)
IMPLICIT NONE
CLASS(meshVolCylQuad), INTENT(inout), TARGET:: elemA
CLASS(meshEdgeCyl), INTENT(inout), TARGET:: elemB
!Check direction 1
IF (.NOT. ASSOCIATED(elemA%e1)) THEN
IF (elemA%n1%n == elemB%n1%n .AND. &
elemA%n2%n == elemB%n2%n) THEN
elemA%e1 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
elemA%n2%n == elemB%n1%n) THEN
elemA%e1 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 2
IF (.NOT. ASSOCIATED(elemA%e2)) THEN
IF (elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n2%n) THEN
elemA%e2 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n1%n) THEN
elemA%e2 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 3
IF (.NOT. ASSOCIATED(elemA%e3)) THEN
IF (elemA%n3%n == elemB%n1%n .AND. &
elemA%n4%n == elemB%n2%n) THEN
elemA%e3 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
elemA%n4%n == elemB%n1%n) THEN
elemA%e3 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 4
IF (.NOT. ASSOCIATED(elemA%e4)) THEN
IF (elemA%n4%n == elemB%n1%n .AND. &
elemA%n1%n == elemB%n2%n) THEN
elemA%e4 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n4%n == elemB%n2%n .AND. &
elemA%n1%n == elemB%n1%n) THEN
elemA%e4 => elemB
elemB%e1 => elemA
END IF
END IF
END SUBROUTINE connectedQuadEdge
SUBROUTINE printOutputCyl(self, t)
USE moduleRefParam
USE moduleSpecies
USE moduleOutput
IMPLICIT NONE
CLASS(meshCylGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
INTEGER:: n, i
TYPE(outputFormat):: output(1:self%numNodes)
REAL(8):: time
CHARACTER(:), ALLOCATABLE:: fileName
CHARACTER (LEN=6):: tstring !TODO: Review to allow any number of iterations
time = DBLE(t)*tau*ti_ref
DO i = 1, nSpecies
WRITE(tstring, '(I6.6)') t
fileName='OUTPUT_' // tstring// '_' // species(i)%obj%name // '.msh'
PRINT *, "Creado archivo de datos:", fileName
OPEN (60, file = path // folder // '/' // fileName)
WRITE(60, "(A)") '$MeshFormat'
WRITE(60, "(A)") '2.2 0 8'
WRITE(60, "(A)") '$EndMeshFormat'
WRITE(60, "(A)") '$NodeData'
WRITE(60, "(A)") '1'
WRITE(60, "(A)") '"Density (m^-3)"'
WRITE(60, *) 1
WRITE(60, *) time
WRITE(60, *) 3
WRITE(60, *) t
WRITE(60, *) 1
WRITE(60, *) self%numNodes
DO n=1, self%numNodes
CALL calculateOutput(self%nodes(n)%obj%output(i), output(n), self%nodes(n)%obj%v, species(i)%obj)
WRITE(60, "(I6,ES20.6E3)") n, output(n)%density
END DO
WRITE(60, "(A)") '$EndNodeData'
WRITE(60, "(A)") '$NodeData'
WRITE(60, "(A)") '1'
WRITE(60, "(A)") '"Velocity (m/s)"'
WRITE(60, *) 1
WRITE(60, *) time
WRITE(60, *) 3
WRITE(60, *) t
WRITE(60, *) 3
WRITE(60, *) self%numNodes
DO n=1, self%numNodes
WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%velocity
END DO
WRITE(60, "(A)") '$EndNodeData'
WRITE(60, "(A)") '$NodeData'
WRITE(60, "(A)") '1'
WRITE(60, "(A)") '"Pressure (Pa)"'
WRITE(60, *) 1
WRITE(60, *) time
WRITE(60, *) 3
WRITE(60, *) t
WRITE(60, *) 1
WRITE(60, *) self%numNodes
DO n=1, self%numNodes
WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%pressure
END DO
WRITE(60, "(A)") '$EndNodeData'
WRITE(60, "(A)") '$NodeData'
WRITE(60, "(A)") '1'
WRITE(60, "(A)") '"Temperature (K)"'
WRITE(60, *) 1
WRITE(60, *) time
WRITE(60, *) 3
WRITE(60, *) t
WRITE(60, *) 1
WRITE(60, *) self%numNodes
DO n=1, self%numNodes
WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%temperature
END DO
WRITE(60, "(A)") '$EndNodeData'
CLOSE (60)
END DO
END SUBROUTINE printOutputCyl
SUBROUTINE printCollisionsCyl(self, t)
USE moduleRefParam
USE moduleCaseParam
USE moduleOutput
IMPLICIT NONE
CLASS(meshCylGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
INTEGER:: n
REAL(8):: time
CHARACTER(:), ALLOCATABLE:: fileName
CHARACTER (LEN=6):: tstring !TODO: Review to allow any number of iterations
IF (collOutput) THEN
time = DBLE(t)*tau*ti_ref
WRITE(tstring, '(I6.6)') t
fileName='OUTPUT_' // tstring// '_Collisions.msh'
PRINT *, "Creado archivo de datos:", fileName
OPEN (60, file = path // folder // '/' // fileName)
WRITE(60, "(A)") '$MeshFormat'
WRITE(60, "(A)") '2.2 0 8'
WRITE(60, "(A)") '$EndMeshFormat'
WRITE(60, "(A)") '$ElementData'
WRITE(60, "(A)") '1'
WRITE(60, "(A)") '"Collisions"'
WRITE(60, *) 1
WRITE(60, *) time
WRITE(60, *) 3
WRITE(60, *) t
WRITE(60, *) 1
WRITE(60, *) self%numVols
DO n=1, self%numVols
WRITE(60, "(I6,I10)") n + self%numEdges, self%vols(n)%obj%nColl
END DO
WRITE(60, "(A)") '$EndElementData'
CLOSE(60)
END IF
END SUBROUTINE printCollisionsCyl
END MODULE moduleMeshCylRead

118
src/modules/moduleOutput.f95 Normal file
View file

@ -0,0 +1,118 @@
!Contains information about output
MODULE moduleOutput
IMPLICIT NONE
!Output for each node
TYPE outputNode
REAL(8):: den, mom(1:3), tensorS(1:3,1:3)
END TYPE
!Output in dimensional units to print
TYPE outputFormat
REAL(8):: density, velocity(1:3), pressure, temperature
END TYPE
CHARACTER(:), ALLOCATABLE:: path
CHARACTER(:), ALLOCATABLE:: folder
INTEGER:: triggerOutput, counterOutput
LOGICAL:: timeOutput = .FALSE.
LOGICAL:: collOutput = .FALSE.
CONTAINS
FUNCTION outerProduct(a,b) RESULT(s)
IMPLICIT NONE
REAL(8), DIMENSION(1:3):: a, b
REAL(8):: s(1:3,1:3)
s = SPREAD(a, dim = 2, ncopies = 3)*SPREAD(b, dim = 1, ncopies = 3)
END FUNCTION outerProduct
FUNCTION tensorTrace(a) RESULT(t)
IMPLICIT NONE
REAL(8), DIMENSION(1:3,1:3):: a
REAL(8):: t
t = 0.D0
t = a(1,1)+a(2,2)+a(3,3)
END FUNCTION tensorTrace
SUBROUTINE calculateOutput(rawValues, formatValues, nodeVol, speciesIn)
USE moduleRefParam
USE moduleSpecies
IMPLICIT NONE
TYPE(outputNode), INTENT(in):: rawValues
TYPE(outputFormat), INTENT(out):: formatValues
REAL(8), INTENT(in):: nodeVol
CLASS(speciesGeneric), INTENT(in):: speciesIn
REAL(8), DIMENSION(1:3,1:3):: tensorTemp
REAL(8), DIMENSION(1:3):: tempVel
REAL(8):: tempVol
!Resets the node outputs
formatValues%density = 0.D0
formatValues%velocity = 0.D0
formatValues%pressure = 0.D0
formatValues%temperature = 0.D0
tempVol = speciesIn%weight/(nodeVol*Vol_ref)
IF (rawValues%den > 0.D0) THEN
tempVel = rawValues%mom(:)/rawValues%den
tensorTemp = (rawValues%tensorS(:,:) - rawValues%den*outerProduct(tempVel,tempVel))
formatValues%density = rawValues%den*tempVol
formatValues%velocity(:) = tempVel
IF (tensorTrace(tensorTemp) > 0.D0) THEN
formatValues%pressure = speciesIn%m*tensorTrace(tensorTemp)*tempVol/3.D0
formatValues%temperature = formatValues%pressure/(formatValues%density*kb)
END IF
END IF
formatValues%velocity = formatValues%velocity*v_ref
formatValues%pressure = formatValues%pressure*m_ref*v_ref**2
formatValues%temperature = formatValues%temperature*m_ref*v_ref**2
END SUBROUTINE calculateOutput
SUBROUTINE printTime(t, first)
USE moduleSpecies
USE moduleCompTime
IMPLICIT NONE
INTEGER, INTENT(in):: t
LOGICAL, INTENT(in), OPTIONAL:: first
CHARACTER(:), ALLOCATABLE:: filename
filename = 'cpuTime.dat'
IF (timeOutput) THEN
IF (PRESENT(first)) THEN
IF (first) THEN
OPEN(20, file = path // folder // '/' // filename, action = 'write')
WRITE(20, "(A1, 8X, A1, 9X, A1, 5(A20))") "#","t","n","total","push","reset","collision","weighting"
ELSE
END IF
OPEN(20, file = path // folder // '/' // filename, position = 'append', action = 'write')
ELSE
OPEN(20, file = path // folder // '/' // filename, position = 'append', action = 'write')
END IF
WRITE (20, "(I10, I10, 5(ES20.6E3))") t, n_part_old, tStep, tPush, tReset, tColl, tWeight
CLOSE(20)
END IF
END SUBROUTINE
END MODULE moduleOutput

10
src/modules/moduleRefParam.f95 Normal file
View file

@ -0,0 +1,10 @@
!Reference parameters
MODULE moduleRefParam
USE moduleConstParam
!Parameters that define the problem (inputs)
REAL(8):: n_ref, m_ref, T_ref, r_ref, sigma_ref!, q_ref=1.6022D-19
!Reference parameters for non-dimensional problem
REAL(8):: L_ref, v_ref, ti_ref, Vol_ref
END MODULE moduleRefParam

113
src/modules/moduleSolver.f95 Normal file
View file

@ -0,0 +1,113 @@
MODULE moduleSolver
CONTAINS
SUBROUTINE scatterGrid()
USE moduleSpecies
USE moduleRefParam
USE moduleMesh
USE moduleOutput
INTEGER:: n, e, k
!Cleans previous output
!$OMP DO PRIVATE(k)
DO e = 1, mesh%numNodes
DO k= 1, nSpecies
mesh%nodes(e)%obj%output(k)%den = 0.D0
mesh%nodes(e)%obj%output(k)%mom = 0.D0
mesh%nodes(e)%obj%output(k)%tensorS = 0.D0
END DO
END DO
!$OMP END DO
!Loops over the particles to scatter them
!$OMP DO
DO n=1, n_part_old
CALL mesh%vols(part_old(n)%e_p)%obj%scatter(part_old(n))
END DO
!$OMP END DO
END SUBROUTINE scatterGrid
SUBROUTINE push(part)
USE moduleSpecies
USE moduleMesh
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
REAL(8):: v_p_oh_star(2:3)
TYPE(particle):: part_temp
REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha, alpha
part_temp = part
!z
part_temp%v(1) = part%v(1)
part_temp%r(1) = part%r(1) + part_temp%v(1)*tau
!r,theta
v_p_oh_star(2) = part%v(2)
x_new = part%r(2) + v_p_oh_star(2)*tau
v_p_oh_star(3) = part%v(3)
y_new = v_p_oh_star(3)*tau
r = DSQRT(x_new**2+y_new**2)
part_temp%r(2) = r
alpha = 0.D0!ATAN2(y_new,x_new) !0 in 2D problem
part_temp%r(3) = part%r(3) + alpha
IF (r > 0.D0) THEN
sin_alpha = y_new/r
cos_alpha = x_new/r
ELSE
sin_alpha = 0.D0
cos_alpha = 1.D0
END IF
part_temp%v(2) = cos_alpha*v_p_oh_star(2)+sin_alpha*v_p_oh_star(3)
part_temp%v(3) = -sin_alpha*v_p_oh_star(2)+cos_alpha*v_p_oh_star(3)
part_temp%pt = part%pt
part_temp%e_p = part%e_p
!Assign cell to particle
part=part_temp
CALL mesh%vols(part%e_p)%obj%findCell(part)
END SUBROUTINE push
SUBROUTINE resetParticles()
USE moduleSpecies
USE moduleList
USE moduleMesh
IMPLICIT NONE
INTEGER:: nn, n
TYPE(particle), ALLOCATABLE:: partTemp(:)
IF (n_part_old > 0) THEN
n_part_new = COUNT(part_old%n_in) + COUNT(part_inj%n_in)
ELSE
n_part_new = COUNT(part_inj%n_in)
END IF
CALL MOVE_ALLOC(part_old, partTemp)
ALLOCATE(part_old(1:n_part_new))
nn = 0
DO n = 1, nPartInj
IF (part_inj(n)%n_in) THEN
nn = nn + 1
part_old(nn) = part_inj(n)
CALL mesh%vols(part_old(nn)%e_p)%obj%listPart_in%add(nn)
END IF
END DO
DO n = 1, n_part_old
IF (partTemp(n)%n_in) THEN
nn = nn + 1
part_old(nn) = partTemp(n)
CALL mesh%vols(part_old(nn)%e_p)%obj%listPart_in%add(nn)
END IF
END DO
n_part_old = n_part_new
END SUBROUTINE resetParticles
END MODULE moduleSolver

66
src/modules/moduleSpecies.f95 Normal file
View file

@ -0,0 +1,66 @@
!Contains the information about species (particles)
MODULE moduleSpecies
USE moduleCaseParam
USE moduleList
IMPLICIT NONE
TYPE, ABSTRACT:: speciesGeneric
CHARACTER(:), ALLOCATABLE:: name
REAL(8):: m=0.D0, weight=0.D0
INTEGER:: pt=0
END TYPE speciesGeneric
TYPE, EXTENDS(speciesGeneric):: speciesNeutral
END TYPE speciesNeutral
TYPE, EXTENDS(speciesGeneric):: speciesCharged
REAL(8):: q=0.D0, qm=0.D0
END TYPE speciesCharged
TYPE:: speciesCont
CLASS(speciesGeneric), ALLOCATABLE:: obj
END TYPE
INTEGER:: nSpecies
TYPE(speciesCont), ALLOCATABLE:: species(:)
TYPE particle
REAL(8):: r(1:3) !Position
REAL(8):: v(1:3) !Velocity
INTEGER:: pt !Particle species id
INTEGER:: e_p !Index of element in which the particle is located
REAL(8):: xLog(1:2) !Logical coordinates of particle in element e_p.
LOGICAL:: n_in !Flag that indicates if a particle is in the domain
END TYPE particle
INTEGER:: n_part_old, n_part_new
INTEGER:: nPartInj
!Arrays that define the particles
TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: part_old
TYPE(particle), ALLOCATABLE, DIMENSION(:):: part_new
TYPE(particle), ALLOCATABLE, DIMENSION(:):: part_inj !array of inject particles
CONTAINS
FUNCTION speciesName2Index(speciesName) RESULT(pt)
IMPLICIT NONE
CHARACTER(:), ALLOCATABLE:: speciesName
INTEGER:: pt
INTEGER:: n
pt = 0
DO n = 1, nSpecies
IF (speciesName == species(n)%obj%name) THEN
pt = species(n)%obj%pt
EXIT
END IF
END DO
!TODO: add error handling when species not found
END FUNCTION speciesName2Index
END MODULE moduleSpecies

121
src/modules/moduleTable.f95 Normal file
View file

@ -0,0 +1,121 @@
MODULE moduleTable
TYPE:: table1D
REAL(8):: xMin, xMax
REAL(8):: fMin, fMax
REAL(8), ALLOCATABLE, DIMENSION(:):: x, f, k
CONTAINS
PROCEDURE, PASS:: init => initTable1D
PROCEDURE, PASS:: get => getValueTable1D
PROCEDURE, PASS:: convert => convertUnits
END TYPE table1D
CONTAINS
SUBROUTINE initTable1D(self, tableFile)
USE moduleErrors
IMPLICIT NONE
CLASS(table1D), INTENT(inout):: self
CHARACTER(:), ALLOCATABLE, INTENT(IN):: tableFile
CHARACTER(100):: dummy
INTEGER:: amount
INTEGER:: i
INTEGER:: stat
INTEGER:: id = 20
OPEN(id, file = tableFile)
amount = 0
DO
READ(id, '(A)', iostat = stat) dummy
!Skip comment
IF (INDEX(dummy,'#') /= 0) CYCLE
!If EOF or error, exit file
IF (stat /= 0) EXIT
!Add row
amount = amount + 1
END DO
IF (amount == 0) CALL criticalError('Table ' // tableFile // ' is empty', 'initTable1D')
IF (amount == 1) CALL criticalError('Table ' // tableFile // ' has only one row', 'initTable1D')
!Go bback to initial point
REWIND(id)
!Allocate table arrays
ALLOCATE(self%x(1:amount), self%f(1:amount), self%k(1:amount))
self%x = 0.D0
self%f = 0.D0
self%k = 0.D0
i = 0
DO
READ(id, '(A)', iostat = stat) dummy
IF (INDEX(dummy,'#') /= 0) CYCLE
IF (stat /= 0) EXIT
!Add data
!TODO: substitute with extracting information from dummy
BACKSPACE(id)
i = i + 1
READ(id, *) self%x(i), self%f(i)
END DO
CLOSE(10)
self%xMin = self%x(1)
self%xMax = self%x(amount)
self%fMin = self%f(1)
self%fMax = self%f(amount)
DO i = 1, amount - 1
self%k(i) = (self%f(i+1) - self%f(i))/(self%x(i+1) - self%x(i))
END DO
END SUBROUTINE initTable1D
FUNCTION getValueTable1D(self, x) RESULT(f)
IMPLICIT NONE
CLASS(table1D), INTENT(in):: self
REAL(8):: x
REAL(8):: f
REAL(8):: deltaX
INTEGER:: i
IF (x <= self%xMin) THEN
f = self%fMin
ELSEIF (x >= self%xMax) THEN
f = self%fMax
ELSE
i = MINLOC(x - self%x, 1)
deltaX = x - self%x(i)
IF (deltaX < 0 ) THEN
i = i - 1
deltaX = x - self%x(i)
END IF
f = self%f(i) + self%k(i)*deltaX
END IF
END FUNCTION getValueTable1D
SUBROUTINE convertUnits(self, data_x, data_f)
IMPLICIT NONE
CLASS(table1D), INTENT(inout):: self
REAL(8):: data_x, data_f
self%x = self%x * data_x
self%f = self%f * data_f
self%k = self%k * data_f / data_x
END SUBROUTINE convertUnits
END MODULE moduleTable