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

Time variable Dirichlet condition #52

Merged
JorgeGonz merged 8 commits from feature/temporalDirichlet into development 2024-07-13 13:01:11 +02:00
Conversation 2 Commits 8 Files changed 16 +340 -220
14 changed files with 340 additions and 220 deletions
Download patch file Download diff file Expand all files Collapse all files
Showing only changes of commit ac27725940 - Show all commits

Big one...

Prev Next 
I should've commited before, but I wanted to make things compile.

The big change is that I've added a global time step so the parameter
does not need to be passed in each function. This is useful as we are
moving towards using time profiles for boundary conditions and injection
of particles (not in this branch, but in the future and the procedure
will be quite similar)
Jorge Gonzalez 2024-07-12 23:08:19 +02:00
commit ac27725940

25
src/fpakc.f90
View file

@ -11,12 +11,12 @@ PROGRAM fpakc
USE OMP_LIB USE OMP_LIB
IMPLICIT NONE IMPLICIT NONE
! t = time step
INTEGER:: t
! arg1 = Input argument 1 (input file) ! arg1 = Input argument 1 (input file)
CHARACTER(200):: arg1 CHARACTER(200):: arg1
! inputFile = path+name of input file ! inputFile = path+name of input file
CHARACTER(:), ALLOCATABLE:: inputFile CHARACTER(:), ALLOCATABLE:: inputFile
! generic integer for time step
INTEGER:: t
tStep = omp_get_wtime() tStep = omp_get_wtime()
!Gets the input file !Gets the input file
@ -32,10 +32,13 @@ PROGRAM fpakc
CALL initOutput(inputFile) CALL initOutput(inputFile)
!Do '0' iteration !Do '0' iteration
t = tInitial timeStep = tInitial
!$OMP PARALLEL DEFAULT(SHARED) !$OMP PARALLEL DEFAULT(SHARED)
!$OMP SINGLE !$OMP SINGLE
! Initial reset of probes
CALL resetProbes()
CALL verboseError("Initial scatter of particles...") CALL verboseError("Initial scatter of particles...")
!$OMP END SINGLE !$OMP END SINGLE
CALL doScatter() CALL doScatter()
@ -49,19 +52,21 @@ PROGRAM fpakc
tStep = omp_get_wtime() - tStep tStep = omp_get_wtime() - tStep
!Output initial state !Output initial state
CALL doOutput(t) CALL doOutput()
CALL verboseError('Starting main loop...') CALL verboseError('Starting main loop...')
!$OMP PARALLEL DEFAULT(SHARED) !$OMP PARALLEL DEFAULT(SHARED)
DO t = tInitial + 1, tFinal DO t = tInitial + 1, tFinal
!Insert new particles and push them
!$OMP SINGLE !$OMP SINGLE
tStep = omp_get_wtime() tStep = omp_get_wtime()
! Update global time step index
timeStep = t
!Checks if a species needs to me moved in this iteration !Checks if a species needs to me moved in this iteration
CALL solver%updatePushSpecies(t) CALL solver%updatePushSpecies()
!Checks if probes need to be calculated this iteration !Checks if probes need to be calculated this iteration
CALL resetProbes(t) CALL resetProbes()
tPush = omp_get_wtime() tPush = omp_get_wtime()
!$OMP END SINGLE !$OMP END SINGLE
@ -79,7 +84,7 @@ PROGRAM fpakc
!$OMP END SINGLE !$OMP END SINGLE
IF (doMCCollisions) THEN IF (doMCCollisions) THEN
CALL meshForMCC%doCollisions(t) CALL meshForMCC%doCollisions()
END IF END IF
@ -124,12 +129,12 @@ PROGRAM fpakc
!$OMP SINGLE !$OMP SINGLE
tEMField = omp_get_wtime() - tEMField tEMField = omp_get_wtime() - tEMField
CALL doAverage(t) CALL doAverage()
tStep = omp_get_wtime() - tStep tStep = omp_get_wtime() - tStep
!Output data !Output data
CALL doOutput(t) CALL doOutput()
!$OMP END SINGLE !$OMP END SINGLE
END DO END DO

5
src/modules/common/moduleCaseParam.f90
View file

@ -2,8 +2,13 @@
MODULE moduleCaseParam MODULE moduleCaseParam
!Final and initial iterations !Final and initial iterations
INTEGER:: tFinal, tInitial = 0 INTEGER:: tFinal, tInitial = 0
! Global index of current iteration
INTEGER:: timeStep
! Time step for all species
REAL(8), ALLOCATABLE:: tau(:) REAL(8), ALLOCATABLE:: tau(:)
! Minimum time step
REAL(8):: tauMin REAL(8):: tauMin
! Time step for Monte-Carlo Collisions
REAL(8):: tauColl REAL(8):: tauColl
END MODULE moduleCaseParam END MODULE moduleCaseParam

134
src/modules/init/moduleInput.f90
View file

@ -264,8 +264,8 @@ MODULE moduleInput
CALL readEMBoundary(config) CALL readEMBoundary(config)
!Read constant magnetic field !Read constant magnetic field
DO i = 1, 3 DO i = 1, 3
WRITE(istring, '(i2)') i WRITE(iString, '(i2)') i
CALL config%get(object // '.B(' // istring // ')', B(i), found) CALL config%get(object // '.B(' // iString // ')', B(i), found)
IF (.NOT. found) THEN IF (.NOT. found) THEN
CALL criticalError('Constant magnetic field not provided in direction ' // iString, 'readSolver') CALL criticalError('Constant magnetic field not provided in direction ' // iString, 'readSolver')
@ -799,7 +799,7 @@ MODULE moduleInput
TYPE(json_file), INTENT(inout):: config TYPE(json_file), INTENT(inout):: config
INTEGER:: i, s INTEGER:: i, s
CHARACTER(2):: istring, sString CHARACTER(2):: iString, sString
CHARACTER(:), ALLOCATABLE:: object, bType CHARACTER(:), ALLOCATABLE:: object, bType
REAL(8):: Tw, cw !Wall temperature and specific heat REAL(8):: Tw, cw !Wall temperature and specific heat
!Neutral Properties !Neutral Properties
@ -815,8 +815,8 @@ MODULE moduleInput
CALL config%info('boundary', found, n_children = nBoundary) CALL config%info('boundary', found, n_children = nBoundary)
ALLOCATE(boundary(1:nBoundary)) ALLOCATE(boundary(1:nBoundary))
DO i = 1, nBoundary DO i = 1, nBoundary
WRITE(istring, '(i2)') i WRITE(iString, '(i2)') i
object = 'boundary(' // TRIM(istring) // ')' object = 'boundary(' // TRIM(iString) // ')'
boundary(i)%n = i boundary(i)%n = i
CALL config%get(object // '.name', boundary(i)%name, found) CALL config%get(object // '.name', boundary(i)%name, found)
@ -1100,13 +1100,13 @@ MODULE moduleInput
TYPE(json_file), INTENT(inout):: config TYPE(json_file), INTENT(inout):: config
CHARACTER(:), ALLOCATABLE:: object CHARACTER(:), ALLOCATABLE:: object
LOGICAL:: found LOGICAL:: found
CHARACTER(2):: istring CHARACTER(2):: iString
INTEGER:: i INTEGER:: i
CHARACTER(:), ALLOCATABLE:: speciesName CHARACTER(:), ALLOCATABLE:: speciesName
REAL(8), ALLOCATABLE, DIMENSION(:):: r REAL(8), ALLOCATABLE, DIMENSION(:):: r
REAL(8), ALLOCATABLE, DIMENSION(:):: v1, v2, v3 REAL(8), ALLOCATABLE, DIMENSION(:):: v1, v2, v3
INTEGER, ALLOCATABLE, DIMENSION(:):: points INTEGER, ALLOCATABLE, DIMENSION(:):: points
REAL(8):: timeStep REAL(8):: everyTimeStep
CALL config%info('output.probes', found, n_children = nProbes) CALL config%info('output.probes', found, n_children = nProbes)
@ -1114,7 +1114,7 @@ MODULE moduleInput
DO i = 1, nProbes DO i = 1, nProbes
WRITE(iString, '(I2)') i WRITE(iString, '(I2)') i
object = 'output.probes(' // trim(istring) // ')' object = 'output.probes(' // trim(iString) // ')'
CALL config%get(object // '.species', speciesName, found) CALL config%get(object // '.species', speciesName, found)
CALL config%get(object // '.position', r, found) CALL config%get(object // '.position', r, found)
@ -1122,16 +1122,14 @@ MODULE moduleInput
CALL config%get(object // '.velocity_2', v2, found) CALL config%get(object // '.velocity_2', v2, found)
CALL config%get(object // '.velocity_3', v3, found) CALL config%get(object // '.velocity_3', v3, found)
CALL config%get(object // '.points', points, found) CALL config%get(object // '.points', points, found)
CALL config%get(object // '.timeStep', timeStep, found) CALL config%get(object // '.timeStep', everyTimeStep, found)
IF (ANY(points < 2)) CALL criticalError("Number of points in probe " // iString // " incorrect", 'readProbes') 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, timeStep) CALL probe(i)%init(i, speciesName, r, v1, v2, v3, points, everyTimeStep)
END DO END DO
CALL resetProbes(tInitial)
END SUBROUTINE readProbes END SUBROUTINE readProbes
SUBROUTINE readEMBoundary(config) SUBROUTINE readEMBoundary(config)
@ -1147,40 +1145,56 @@ MODULE moduleInput
TYPE(json_file), INTENT(inout):: config TYPE(json_file), INTENT(inout):: config
CHARACTER(:), ALLOCATABLE:: object CHARACTER(:), ALLOCATABLE:: object
LOGICAL:: found LOGICAL:: found
CHARACTER(2):: istring CHARACTER(:), ALLOCATABLE:: typeEM
INTEGER:: i, e, s REAL(8):: potential
INTEGER:: physicalSurface
CHARACTER(:), ALLOCATABLE:: timeProfile
INTEGER:: b, e, s, n, ni
CHARACTER(2):: bString
INTEGER:: info INTEGER:: info
EXTERNAL:: dgetrf EXTERNAL:: dgetrf
CALL config%info('boundaryEM', found, n_children = nBoundaryEM) CALL config%info('boundaryEM', found, n_children = nBoundaryEM)
IF (found) ALLOCATE(boundEM(1:nBoundaryEM)) IF (found) ALLOCATE(boundaryEM(1:nBoundaryEM))
DO i = 1, nBoundaryEM DO b = 1, nBoundaryEM
WRITE(istring, '(I2)') i WRITE(bString, '(I2)') b
object = 'boundaryEM(' // trim(istring) // ')' object = 'boundaryEM(' // trim(bString) // ')'
CALL config%get(object // '.type', boundEM(i)%typeEM, found) CALL config%get(object // '.type', typeEM, found)
SELECT CASE(boundEM(i)%typeEM) SELECT CASE(typeEM)
CASE ("dirichlet") CASE ("dirichlet")
CALL config%get(object // '.potential', boundEM(i)%potential, found) CALL config%get(object // '.potential', potential, found)
IF (.NOT. found) & IF (.NOT. found) THEN
CALL criticalError('Required parameter "potential" for Dirichlet boundary condition not found', 'readEMBoundary') CALL criticalError('Required parameter "potential" for Dirichlet boundary condition not found', 'readEMBoundary')
boundEM(i)%potential = boundEM(i)%potential/Volt_ref
CALL config%get(object // '.physicalSurface', boundEM(i)%physicalSurface, found) END IF
IF (.NOT. found) &
CALL config%get(object // '.physicalSurface', physicalSurface, found)
IF (.NOT. found) THEN
CALL criticalError('Required parameter "physicalSurface" for Dirichlet boundary condition not found', 'readEMBoundary') CALL criticalError('Required parameter "physicalSurface" for Dirichlet boundary condition not found', 'readEMBoundary')
END IF
CALL initDirichlet(boundaryEM(b)%obj, physicalSurface, potential)
CASE ("dirichletTime") CASE ("dirichletTime")
CALL config%get(object // '.potential', boundEM(i)%potential, found) CALL config%get(object // '.potential', potential, found)
IF (.NOT. found) & IF (.NOT. found) THEN
CALL criticalError('Required parameter "potential" for Dirichlet boundary condition not found', 'readEMBoundary') CALL criticalError('Required parameter "potential" for Dirichlet boundary condition not found', 'readEMBoundary')
boundEM(i)%potential = boundEM(i)%potential/Volt_ref
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
CASE DEFAULT CASE DEFAULT
CALL criticalError('Boundary type ' // boundEM(i)%typeEM // ' not yet supported', 'readEMBoundary') CALL criticalError('Boundary type ' // typeEM // ' not yet supported', 'readEMBoundary')
END SELECT END SELECT
@ -1199,18 +1213,28 @@ MODULE moduleInput
END DO END DO
IF (ALLOCATED(boundEM)) THEN ! Modify K matrix due to boundary conditions
DO e = 1, mesh%numEdges DO b = 1, nBoundaryEM
IF (ANY(mesh%edges(e)%obj%physicalSurface == boundEM(:)%physicalSurface)) THEN SELECT TYPE(boundary => boundaryEM(b)%obj)
DO i = 1, nBoundaryEM TYPE IS(boundaryEMDirichlet)
IF (mesh%edges(e)%obj%physicalSurface == boundEM(i)%physicalSurface) THEN DO n = 1, boundary%nNodes
CALL boundEM(i)%apply(mesh%edges(e)%obj) ni = boundary%nodes(n)%obj%n
mesh%K(ni, :) = 0.D0
mesh%K(ni, ni) = 1.D0
END IF END DO
END DO
END IF TYPE IS(boundaryEMDirichletTime)
END DO DO n = 1, boundary%nNodes
END IF 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 !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) CALL dgetrf(mesh%numNodes, mesh%numNodes, mesh%K, mesh%numNodes, mesh%IPIV, info)
@ -1231,13 +1255,13 @@ MODULE moduleInput
TYPE(json_file), INTENT(inout):: config TYPE(json_file), INTENT(inout):: config
INTEGER:: i INTEGER:: i
CHARACTER(2):: istring CHARACTER(2):: iString
CHARACTER(:), ALLOCATABLE:: object CHARACTER(:), ALLOCATABLE:: object
LOGICAL:: found LOGICAL:: found
CHARACTER(:), ALLOCATABLE:: speciesName CHARACTER(:), ALLOCATABLE:: speciesName
CHARACTER(:), ALLOCATABLE:: name CHARACTER(:), ALLOCATABLE:: name
REAL(8):: v REAL(8):: v
REAL(8), ALLOCATABLE:: T(:), normal(:) REAL(8), ALLOCATABLE:: temperature(:), normal(:)
REAL(8):: flow REAL(8):: flow
CHARACTER(:), ALLOCATABLE:: units CHARACTER(:), ALLOCATABLE:: units
INTEGER:: physicalSurface INTEGER:: physicalSurface
@ -1248,8 +1272,8 @@ MODULE moduleInput
ALLOCATE(inject(1:nInject)) ALLOCATE(inject(1:nInject))
nPartInj = 0 nPartInj = 0
DO i = 1, nInject DO i = 1, nInject
WRITE(istring, '(i2)') i WRITE(iString, '(i2)') i
object = 'inject(' // trim(istring) // ')' object = 'inject(' // trim(iString) // ')'
!Find species !Find species
CALL config%get(object // '.species', speciesName, found) CALL config%get(object // '.species', speciesName, found)
@ -1257,7 +1281,7 @@ MODULE moduleInput
CALL config%get(object // '.name', name, found) CALL config%get(object // '.name', name, found)
CALL config%get(object // '.v', v, found) CALL config%get(object // '.v', v, found)
CALL config%get(object // '.T', T, found) CALL config%get(object // '.T', temperature, found)
CALL config%get(object // '.n', normal, found) CALL config%get(object // '.n', normal, found)
IF (.NOT. found) THEN IF (.NOT. found) THEN
ALLOCATE(normal(1:3)) ALLOCATE(normal(1:3))
@ -1269,7 +1293,7 @@ MODULE moduleInput
particlesPerEdge = 0 particlesPerEdge = 0
CALL config%get(object // '.particlesPerEdge', particlesPerEdge, found) CALL config%get(object // '.particlesPerEdge', particlesPerEdge, found)
CALL inject(i)%init(i, v, normal, T, flow, units, sp, physicalSurface, particlesPerEdge) CALL inject(i)%init(i, v, normal, temperature, flow, units, sp, physicalSurface, particlesPerEdge)
CALL readVelDistr(config, inject(i), object) CALL readVelDistr(config, inject(i), object)
@ -1329,28 +1353,28 @@ MODULE moduleInput
TYPE(injectGeneric), INTENT(inout):: inj TYPE(injectGeneric), INTENT(inout):: inj
CHARACTER(:), ALLOCATABLE, INTENT(in):: object CHARACTER(:), ALLOCATABLE, INTENT(in):: object
INTEGER:: i INTEGER:: i
CHARACTER(2):: istring CHARACTER(2):: iString
CHARACTER(:), ALLOCATABLE:: fvType CHARACTER(:), ALLOCATABLE:: fvType
LOGICAL:: found LOGICAL:: found
REAL(8):: v, T, m REAL(8):: v, temperature, m
!Reads species mass !Reads species mass
m = inj%species%m m = inj%species%m
!Reads distribution functions for velocity !Reads distribution functions for velocity
DO i = 1, 3 DO i = 1, 3
WRITE(istring, '(i2)') i WRITE(iString, '(i2)') i
CALL config%get(object // '.velDist('// TRIM(istring) //')', fvType, found) CALL config%get(object // '.velDist('// TRIM(iString) //')', fvType, found)
IF (.NOT. found) CALL criticalError("No velocity distribution in direction " // istring // & IF (.NOT. found) CALL criticalError("No velocity distribution in direction " // iString // &
" found for " // object, 'readVelDistr') " found for " // object, 'readVelDistr')
SELECT CASE(fvType) SELECT CASE(fvType)
CASE ("Maxwellian") CASE ("Maxwellian")
T = inj%T(i) temperature = inj%temperature(i)
CALL initVelDistMaxwellian(inj%v(i)%obj, t, m) CALL initVelDistMaxwellian(inj%v(i)%obj, temperature, m)
CASE ("Half-Maxwellian") CASE ("Half-Maxwellian")
T = inj%T(i) temperature = inj%temperature(i)
CALL initVelDistHalfMaxwellian(inj%v(i)%obj, t, m) CALL initVelDistHalfMaxwellian(inj%v(i)%obj, temperature, m)
CASE ("Delta") CASE ("Delta")
v = inj%vMod*inj%n(i) v = inj%vMod*inj%n(i)

21
src/modules/mesh/inout/0D/moduleMeshOutput0D.f90
View file

@ -1,22 +1,22 @@
MODULE moduleMeshOutput0D MODULE moduleMeshOutput0D
CONTAINS CONTAINS
SUBROUTINE printOutput0D(self, t) SUBROUTINE printOutput0D(self)
USE moduleMesh USE moduleMesh
USE moduleRefParam USE moduleRefParam
USE moduleSpecies USE moduleSpecies
USE moduleOutput USE moduleOutput
USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self CLASS(meshParticles), INTENT(in):: self
INTEGER, INTENT(in):: t
INTEGER:: i INTEGER:: i
TYPE(outputFormat):: output TYPE(outputFormat):: output
CHARACTER(:), ALLOCATABLE:: fileName CHARACTER(:), ALLOCATABLE:: fileName
DO i = 1, nSpecies DO i = 1, nSpecies
fileName='OUTPUT_' // species(i)%obj%name // '.dat' fileName='OUTPUT_' // species(i)%obj%name // '.dat'
IF (t == 0) THEN IF (timeStep == 0) THEN
OPEN(20, file = path // folder // '/' // fileName, action = 'write') OPEN(20, file = path // folder // '/' // fileName, action = 'write')
WRITE(20, "(A1, 14X, A5, A20, 40X, A20, 2(A20))") "#","t (s)","density (m^-3)", "velocity (m/s)", & WRITE(20, "(A1, 14X, A5, A20, 40X, A20, 2(A20))") "#","t (s)","density (m^-3)", "velocity (m/s)", &
"pressure (Pa)", "temperature (K)" "pressure (Pa)", "temperature (K)"
@ -27,14 +27,17 @@ MODULE moduleMeshOutput0D
OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write') OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write')
CALL calculateOutput(self%nodes(1)%obj%output(i), output, self%nodes(1)%obj%v, species(i)%obj) CALL calculateOutput(self%nodes(1)%obj%output(i), output, self%nodes(1)%obj%v, species(i)%obj)
WRITE(20, "(7(ES20.6E3))") REAL(t)*tauMin*ti_ref, output%density, output%velocity, output%pressure, output%temperature WRITE(20, "(7(ES20.6E3))") REAL(timeStep)*tauMin*ti_ref, output%density, &
output%velocity, &
output%pressure, &
output%temperature
CLOSE(20) CLOSE(20)
END DO END DO
END SUBROUTINE printOutput0D END SUBROUTINE printOutput0D
SUBROUTINE printColl0D(self, t) SUBROUTINE printColl0D(self)
USE moduleMesh USE moduleMesh
USE moduleRefParam USE moduleRefParam
USE moduleCaseParam USE moduleCaseParam
@ -43,12 +46,11 @@ MODULE moduleMeshOutput0D
IMPLICIT NONE IMPLICIT NONE
CLASS(meshGeneric), INTENT(in):: self CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
CHARACTER(:), ALLOCATABLE:: fileName CHARACTER(:), ALLOCATABLE:: fileName
INTEGER:: k INTEGER:: k
fileName='OUTPUT_Collisions.dat' fileName='OUTPUT_Collisions.dat'
IF (t == tInitial) THEN IF (timeStep == tInitial) THEN
OPEN(20, file = path // folder // '/' // fileName, action = 'write') OPEN(20, file = path // folder // '/' // fileName, action = 'write')
WRITE(20, "(A1, 14X, A5, A20)") "#","t (s)","collisions" WRITE(20, "(A1, 14X, A5, A20)") "#","t (s)","collisions"
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
@ -57,12 +59,12 @@ MODULE moduleMeshOutput0D
END IF END IF
OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write') OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write')
WRITE(20, "(ES20.6E3, 10I20)") REAL(t)*tauMin*ti_ref, (self%cells(1)%obj%tallyColl(k)%tally, k=1,nCollPairs) WRITE(20, "(ES20.6E3, 10I20)") REAL(timeStep)*tauMin*ti_ref, (self%cells(1)%obj%tallyColl(k)%tally, k=1,nCollPairs)
CLOSE(20) CLOSE(20)
END SUBROUTINE printColl0D END SUBROUTINE printColl0D
SUBROUTINE printEM0D(self, t) SUBROUTINE printEM0D(self)
USE moduleMesh USE moduleMesh
USE moduleRefParam USE moduleRefParam
USE moduleCaseParam USE moduleCaseParam
@ -70,7 +72,6 @@ MODULE moduleMeshOutput0D
IMPLICIT NONE IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self CLASS(meshParticles), INTENT(in):: self
INTEGER, INTENT(in):: t
END SUBROUTINE printEM0D END SUBROUTINE printEM0D

38
src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90
View file

@ -80,50 +80,50 @@ MODULE moduleMeshOutputGmsh2
END SUBROUTINE writeGmsh2FooterElementData END SUBROUTINE writeGmsh2FooterElementData
!Prints the scattered properties of particles into the nodes !Prints the scattered properties of particles into the nodes
SUBROUTINE printOutputGmsh2(self, t) SUBROUTINE printOutputGmsh2(self)
USE moduleMesh USE moduleMesh
USE moduleRefParam USE moduleRefParam
USE moduleSpecies USE moduleSpecies
USE moduleOutput USE moduleOutput
USE moduleMeshInoutCommon USE moduleMeshInoutCommon
USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self CLASS(meshParticles), INTENT(in):: self
INTEGER, INTENT(in):: t
INTEGER:: n, i INTEGER:: n, i
TYPE(outputFormat):: output(1:self%numNodes) TYPE(outputFormat):: output(1:self%numNodes)
REAL(8):: time REAL(8):: time
CHARACTER(:), ALLOCATABLE:: fileName CHARACTER(:), ALLOCATABLE:: fileName
time = DBLE(t)*tauMin*ti_ref time = DBLE(timeStep)*tauMin*ti_ref
DO i = 1, nSpecies DO i = 1, nSpecies
fileName = formatFileName(prefix, species(i)%obj%name, 'msh', t) fileName = formatFileName(prefix, species(i)%obj%name, 'msh', timeStep)
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (60, file = path // folder // '/' // fileName) OPEN (60, file = path // folder // '/' // fileName)
CALL writeGmsh2HeaderMesh(60) CALL writeGmsh2HeaderMesh(60)
CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' density (m^-3)', t, time, 1, self%numNodes) CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' density (m^-3)', timeStep, time, 1, self%numNodes)
DO n=1, 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) 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 WRITE(60, "(I6,ES20.6E3)") n, output(n)%density
END DO END DO
CALL writeGmsh2FooterNodeData(60) CALL writeGmsh2FooterNodeData(60)
CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' velocity (m s^-1)', t, time, 3, self%numNodes) CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' velocity (m s^-1)', timeStep, time, 3, self%numNodes)
DO n=1, self%numNodes DO n=1, self%numNodes
WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%velocity WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%velocity
END DO END DO
CALL writeGmsh2FooterNodeData(60) CALL writeGmsh2FooterNodeData(60)
CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Pressure (Pa)', t, time, 1, self%numNodes) CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Pressure (Pa)', timeStep, time, 1, self%numNodes)
DO n=1, self%numNodes DO n=1, self%numNodes
WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%pressure WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%pressure
END DO END DO
CALL writeGmsh2FooterNodeData(60) CALL writeGmsh2FooterNodeData(60)
CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Temperature (K)', t, time, 1, self%numNodes) CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Temperature (K)', timeStep, time, 1, self%numNodes)
DO n=1, self%numNodes DO n=1, self%numNodes
WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%temperature WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%temperature
END DO END DO
@ -135,7 +135,7 @@ MODULE moduleMeshOutputGmsh2
END SUBROUTINE printOutputGmsh2 END SUBROUTINE printOutputGmsh2
!Prints the number of collisions into the volumes !Prints the number of collisions into the volumes
SUBROUTINE printCollGmsh2(self, t) SUBROUTINE printCollGmsh2(self)
USE moduleMesh USE moduleMesh
USE moduleRefParam USE moduleRefParam
USE moduleCaseParam USE moduleCaseParam
@ -145,7 +145,6 @@ MODULE moduleMeshOutputGmsh2
IMPLICIT NONE IMPLICIT NONE
CLASS(meshGeneric), INTENT(in):: self CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
INTEGER:: numEdges INTEGER:: numEdges
INTEGER:: k, c INTEGER:: k, c
INTEGER:: n INTEGER:: n
@ -167,9 +166,9 @@ MODULE moduleMeshOutputGmsh2
END SELECT END SELECT
IF (collOutput) THEN IF (collOutput) THEN
time = DBLE(t)*tauMin*ti_ref time = DBLE(timeStep)*tauMin*ti_ref
fileName = formatFileName(prefix, 'Collisions', 'msh', t) fileName = formatFileName(prefix, 'Collisions', 'msh', timeStep)
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (60, file = path // folder // '/' // fileName) OPEN (60, file = path // folder // '/' // fileName)
@ -179,7 +178,7 @@ MODULE moduleMeshOutputGmsh2
DO c = 1, interactionMatrix(k)%amount DO c = 1, interactionMatrix(k)%amount
WRITE(cString, "(I2)") c WRITE(cString, "(I2)") c
title = '"Pair ' // interactionMatrix(k)%sp_i%name // '-' // interactionMatrix(k)%sp_j%name // ' collision ' // cString title = '"Pair ' // interactionMatrix(k)%sp_i%name // '-' // interactionMatrix(k)%sp_j%name // ' collision ' // cString
CALL writeGmsh2HeaderElementData(60, title, t, time, 1, self%numCells) CALL writeGmsh2HeaderElementData(60, title, timeStep, time, 1, self%numCells)
DO n=1, self%numCells DO n=1, self%numCells
WRITE(60, "(I6,I10)") n + numEdges, self%cells(n)%obj%tallyColl(k)%tally(c) WRITE(60, "(I6,I10)") n + numEdges, self%cells(n)%obj%tallyColl(k)%tally(c)
END DO END DO
@ -196,7 +195,7 @@ MODULE moduleMeshOutputGmsh2
END SUBROUTINE printCollGmsh2 END SUBROUTINE printCollGmsh2
!Prints the electrostatic EM properties into the nodes and volumes !Prints the electrostatic EM properties into the nodes and volumes
SUBROUTINE printEMGmsh2(self, t) SUBROUTINE printEMGmsh2(self)
USE moduleMesh USE moduleMesh
USE moduleRefParam USE moduleRefParam
USE moduleCaseParam USE moduleCaseParam
@ -205,7 +204,6 @@ MODULE moduleMeshOutputGmsh2
IMPLICIT NONE IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self CLASS(meshParticles), INTENT(in):: self
INTEGER, INTENT(in):: t
INTEGER:: n, e INTEGER:: n, e
REAL(8):: time REAL(8):: time
CHARACTER(:), ALLOCATABLE:: fileName CHARACTER(:), ALLOCATABLE:: fileName
@ -214,27 +212,27 @@ MODULE moduleMeshOutputGmsh2
Xi = (/ 0.D0, 0.D0, 0.D0 /) Xi = (/ 0.D0, 0.D0, 0.D0 /)
IF (emOutput) THEN IF (emOutput) THEN
time = DBLE(t)*tauMin*ti_ref time = DBLE(timeStep)*tauMin*ti_ref
fileName = formatFileName(prefix, 'EMField', 'msh', t) fileName = formatFileName(prefix, 'EMField', 'msh', timeStep)
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (20, file = path // folder // '/' // fileName) OPEN (20, file = path // folder // '/' // fileName)
CALL writeGmsh2HeaderMesh(20) CALL writeGmsh2HeaderMesh(20)
CALL writeGmsh2HeaderNodeData(20, 'Potential (V)', t, time, 1, self%numNodes) CALL writeGmsh2HeaderNodeData(20, 'Potential (V)', timeStep, time, 1, self%numNodes)
DO n=1, self%numNodes DO n=1, self%numNodes
WRITE(20, *) n, self%nodes(n)%obj%emData%phi*Volt_ref WRITE(20, *) n, self%nodes(n)%obj%emData%phi*Volt_ref
END DO END DO
CALL writeGmsh2FooterNodeData(20) CALL writeGmsh2FooterNodeData(20)
CALL writeGmsh2HeaderElementData(20, 'Electric Field (V m^-1)', t, time, 3, self%numCells) CALL writeGmsh2HeaderElementData(20, 'Electric Field (V m^-1)', timeStep, time, 3, self%numCells)
DO e=1, self%numCells DO e=1, self%numCells
WRITE(20, *) e+self%numEdges, self%cells(e)%obj%gatherElectricField(Xi)*EF_ref WRITE(20, *) e+self%numEdges, self%cells(e)%obj%gatherElectricField(Xi)*EF_ref
END DO END DO
CALL writeGmsh2FooterElementData(20) CALL writeGmsh2FooterElementData(20)
CALL writeGmsh2HeaderNodeData(20, 'Magnetic Field (T)', t, time, 3, self%numNodes) CALL writeGmsh2HeaderNodeData(20, 'Magnetic Field (T)', timeStep, time, 3, self%numNodes)
DO n=1, self%numNodes DO n=1, self%numNodes
WRITE(20, *) n, self%nodes(n)%obj%emData%B * B_ref WRITE(20, *) n, self%nodes(n)%obj%emData%B * B_ref
END DO END DO

8
src/modules/mesh/inout/moduleMeshInoutCommon.f90
View file

@ -3,17 +3,17 @@ MODULE moduleMeshInoutCommon
CHARACTER(LEN=4):: prefix = 'Step' CHARACTER(LEN=4):: prefix = 'Step'
CONTAINS CONTAINS
PURE FUNCTION formatFileName(prefix, suffix, extension, t) RESULT(fileName) PURE FUNCTION formatFileName(prefix, suffix, extension, timeStep) RESULT(fileName)
USE moduleOutput USE moduleOutput
IMPLICIT NONE IMPLICIT NONE
CHARACTER(*), INTENT(in):: prefix, suffix, extension CHARACTER(*), INTENT(in):: prefix, suffix, extension
INTEGER, INTENT(in), OPTIONAL:: t INTEGER, INTENT(in), OPTIONAL:: timeStep
CHARACTER (LEN=iterationDigits):: tString CHARACTER (LEN=iterationDigits):: tString
CHARACTER(:), ALLOCATABLE:: fileName CHARACTER(:), ALLOCATABLE:: fileName
IF (PRESENT(t)) THEN IF (PRESENT(timeStep)) THEN
WRITE(tString, iterationFormat) t WRITE(tString, iterationFormat) timeStep
fileName = prefix // '_' // tString // '_' // suffix // '.' // extension fileName = prefix // '_' // tString // '_' // suffix // '.' // extension
ELSE ELSE

2
src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90
View file

@ -167,7 +167,7 @@ MODULE moduleMeshInputVTU
CLASS(meshGeneric), INTENT(inout):: self CLASS(meshGeneric), INTENT(inout):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: filename CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
REAL(8):: r(1:3) !3 generic coordinates REAL(8):: r(1:3) !3 generic coordinates
INTEGER:: fileID, error, found INTEGER:: fileID
CHARACTER(LEN=256):: line CHARACTER(LEN=256):: line
INTEGER:: numNodes, numElements, numEdges INTEGER:: numNodes, numElements, numEdges
INTEGER, ALLOCATABLE, DIMENSION(:):: entitiesID, offsets, connectivity, types INTEGER, ALLOCATABLE, DIMENSION(:):: entitiesID, offsets, connectivity, types

35
src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90
View file

@ -215,17 +215,16 @@ MODULE moduleMeshOutputVTU
END SUBROUTINE writeEM END SUBROUTINE writeEM
SUBROUTINE writeCollection(fileID, t, fileNameStep, fileNameCollection) SUBROUTINE writeCollection(fileID, fileNameStep, fileNameCollection)
USE moduleCaseParam USE moduleCaseParam
USE moduleOutput USE moduleOutput
USE moduleRefParam USE moduleRefParam
IMPLICIT NONE IMPLICIT NONE
INTEGER:: fileID INTEGER:: fileID
INTEGER, INTENT(in):: t
CHARACTER(*):: fileNameStep, fileNameCollection CHARACTER(*):: fileNameStep, fileNameCollection
IF (t == tInitial) THEN IF (timeStep == tInitial) THEN
!Create collection file !Create collection file
WRITE(*, "(6X,A15,A)") "Creating file: ", fileNameCollection WRITE(*, "(6X,A15,A)") "Creating file: ", fileNameCollection
OPEN (fileID + 1, file = path // folder // '/' // fileNameCollection) OPEN (fileID + 1, file = path // folder // '/' // fileNameCollection)
@ -237,10 +236,11 @@ MODULE moduleMeshOutputVTU
!Write iteration file in collection !Write iteration file in collection
OPEN (fileID + 1, file = path // folder // '/' // fileNameCollection, ACCESS='APPEND') OPEN (fileID + 1, file = path // folder // '/' // fileNameCollection, ACCESS='APPEND')
WRITE(fileID + 1, "(4X, A, ES20.6E3, A, A, A)") '<DataSet timestep="', DBLE(t)*tauMin*ti_ref,'" file="', fileNameStep,'"/>' WRITE(fileID + 1, "(4X, A, ES20.6E3, A, A, A)") &
'<DataSet timestep="', DBLE(timeStep)*tauMin*ti_ref,'" file="', fileNameStep,'"/>'
!Close collection file !Close collection file
IF (t == tFinal) THEN IF (timeStep == tFinal) THEN
WRITE (fileID + 1, "(2X, A)") '</Collection>' WRITE (fileID + 1, "(2X, A)") '</Collection>'
WRITE (fileID + 1, "(A)") '</VTKFile>' WRITE (fileID + 1, "(A)") '</VTKFile>'
@ -307,21 +307,21 @@ MODULE moduleMeshOutputVTU
END SUBROUTINE writeAverage END SUBROUTINE writeAverage
SUBROUTINE printOutputVTU(self,t) SUBROUTINE printOutputVTU(self)
USE moduleMesh USE moduleMesh
USE moduleSpecies USE moduleSpecies
USE moduleMeshInoutCommon USE moduleMeshInoutCommon
USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self CLASS(meshParticles), INTENT(in):: self
INTEGER, INTENT(in):: t
INTEGER:: i, fileID INTEGER:: i, fileID
CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection
fileID = 60 fileID = 60
DO i = 1, nSpecies DO i = 1, nSpecies
fileName = formatFileName(prefix, species(i)%obj%name, 'vtu', t) fileName = formatFileName(prefix, species(i)%obj%name, 'vtu', timeStep)
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (fileID, file = path // folder // '/' // fileName) OPEN (fileID, file = path // folder // '/' // fileName)
@ -337,28 +337,27 @@ MODULE moduleMeshOutputVTU
!Write collection file for time plotting !Write collection file for time plotting
fileNameCollection = formatFileName('Collection', species(i)%obj%name, 'pvd') fileNameCollection = formatFileName('Collection', species(i)%obj%name, 'pvd')
CALL writeCollection(fileID, t, fileName, filenameCollection) CALL writeCollection(fileID, fileName, filenameCollection)
END DO END DO
END SUBROUTINE printOutputVTU END SUBROUTINE printOutputVTU
SUBROUTINE printCollVTU(self,t) SUBROUTINE printCollVTU(self)
USE moduleMesh USE moduleMesh
USE moduleOutput USE moduleOutput
USE moduleMeshInoutCommon USE moduleMeshInoutCommon
USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE IMPLICIT NONE
CLASS(meshGeneric), INTENT(in):: self CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
INTEGER:: fileID INTEGER:: fileID
CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection
CHARACTER (LEN=iterationDigits):: tstring
fileID = 62 fileID = 62
IF (collOutput) THEN IF (collOutput) THEN
fileName = formatFileName(prefix, 'Collisions', 'vtu', t) fileName = formatFileName(prefix, 'Collisions', 'vtu', timeStep)
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (fileID, file = path // folder // '/' // fileName) OPEN (fileID, file = path // folder // '/' // fileName)
@ -374,26 +373,26 @@ MODULE moduleMeshOutputVTU
!Write collection file for time plotting !Write collection file for time plotting
fileNameCollection = formatFileName('Collection', 'Collisions', 'pvd') fileNameCollection = formatFileName('Collection', 'Collisions', 'pvd')
CALL writeCollection(fileID, t, fileName, filenameCollection) CALL writeCollection(fileID, fileName, filenameCollection)
END IF END IF
END SUBROUTINE printCollVTU END SUBROUTINE printCollVTU
SUBROUTINE printEMVTU(self, t) SUBROUTINE printEMVTU(self)
USE moduleMesh USE moduleMesh
USE moduleMeshInoutCommon USE moduleMeshInoutCommon
USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self CLASS(meshParticles), INTENT(in):: self
INTEGER, INTENT(in):: t
INTEGER:: fileID INTEGER:: fileID
CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection
fileID = 64 fileID = 64
IF (emOutput) THEN IF (emOutput) THEN
fileName = formatFileName(prefix, 'EMField', 'vtu', t) fileName = formatFileName(prefix, 'EMField', 'vtu', timeStep)
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (fileID, file = path // folder // '/' // fileName) OPEN (fileID, file = path // folder // '/' // fileName)
@ -409,7 +408,7 @@ MODULE moduleMeshOutputVTU
!Write collection file for time plotting !Write collection file for time plotting
fileNameCollection = formatFileName('Collection', 'EMField', 'pvd') fileNameCollection = formatFileName('Collection', 'EMField', 'pvd')
CALL writeCollection(fileID, t, fileName, filenameCollection) CALL writeCollection(fileID, fileName, filenameCollection)
END IF END IF

15
src/modules/mesh/moduleMesh.f90
View file

@ -372,10 +372,9 @@ MODULE moduleMesh
END SUBROUTINE connectMesh_interface END SUBROUTINE connectMesh_interface
!Prints number of collisions in each cell !Prints number of collisions in each cell
SUBROUTINE printColl_interface(self, t) SUBROUTINE printColl_interface(self)
IMPORT meshGeneric IMPORT meshGeneric
CLASS(meshGeneric), INTENT(in):: self CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
END SUBROUTINE printColl_interface END SUBROUTINE printColl_interface
@ -403,18 +402,16 @@ MODULE moduleMesh
ABSTRACT INTERFACE ABSTRACT INTERFACE
!Prints Species data !Prints Species data
SUBROUTINE printOutput_interface(self, t) SUBROUTINE printOutput_interface(self)
IMPORT meshParticles IMPORT meshParticles
CLASS(meshParticles), INTENT(in):: self CLASS(meshParticles), INTENT(in):: self
INTEGER, INTENT(in):: t
END SUBROUTINE printOutput_interface END SUBROUTINE printOutput_interface
!Prints EM info !Prints EM info
SUBROUTINE printEM_interface(self, t) SUBROUTINE printEM_interface(self)
IMPORT meshParticles IMPORT meshParticles
CLASS(meshParticles), INTENT(in):: self CLASS(meshParticles), INTENT(in):: self
INTEGER, INTENT(in):: t
END SUBROUTINE printEM_interface END SUBROUTINE printEM_interface
@ -789,7 +786,7 @@ MODULE moduleMesh
END FUNCTION findCellBrute END FUNCTION findCellBrute
!Computes collisions in element !Computes collisions in element
SUBROUTINE doCollisions(self, t) SUBROUTINE doCollisions(self)
USE moduleCollisions USE moduleCollisions
USE moduleSpecies USE moduleSpecies
USE moduleList USE moduleList
@ -797,10 +794,10 @@ MODULE moduleMesh
USE moduleRandom USE moduleRandom
USE moduleOutput USE moduleOutput
USE moduleMath USE moduleMath
USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE IMPLICIT NONE
CLASS(meshGeneric), INTENT(inout), TARGET:: self CLASS(meshGeneric), INTENT(inout), TARGET:: self
INTEGER, INTENT(in):: t
INTEGER:: e INTEGER:: e
CLASS(meshCell), POINTER:: cell CLASS(meshCell), POINTER:: cell
INTEGER:: k, i, j INTEGER:: k, i, j
@ -816,7 +813,7 @@ MODULE moduleMesh
REAL(8):: rnd_real !Random number for collision REAL(8):: rnd_real !Random number for collision
INTEGER:: rnd_int !Random number for collision INTEGER:: rnd_int !Random number for collision
IF (MOD(t, everyColl) == 0) THEN IF (MOD(timeStep, everyColl) == 0) THEN
!Collisions need to be performed in this iteration !Collisions need to be performed in this iteration
!$OMP DO SCHEDULE(DYNAMIC) PRIVATE(part_i, part_j, partTemp_i, partTemp_j) !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(part_i, part_j, partTemp_i, partTemp_j)
DO e=1, self%numCells DO e=1, self%numCells

26
src/modules/moduleInject.f90
View file

@ -54,7 +54,7 @@ MODULE moduleInject
INTEGER:: id INTEGER:: id
CHARACTER(:), ALLOCATABLE:: name CHARACTER(:), ALLOCATABLE:: name
REAL(8):: vMod !Velocity (module) REAL(8):: vMod !Velocity (module)
REAL(8):: T(1:3) !Temperature REAL(8):: temperature(1:3) !Temperature
REAL(8):: n(1:3) !Direction of injection REAL(8):: n(1:3) !Direction of injection
LOGICAL:: fixDirection !The injection of particles has a fix direction defined by n LOGICAL:: fixDirection !The injection of particles has a fix direction defined by n
INTEGER:: nParticles !Number of particles to introduce each time step INTEGER:: nParticles !Number of particles to introduce each time step
@ -76,7 +76,7 @@ MODULE moduleInject
CONTAINS CONTAINS
!Initialize an injection of particles !Initialize an injection of particles
SUBROUTINE initInject(self, i, v, n, T, flow, units, sp, physicalSurface, particlesPerEdge) SUBROUTINE initInject(self, i, v, n, temperature, flow, units, sp, physicalSurface, particlesPerEdge)
USE moduleMesh USE moduleMesh
USE moduleRefParam USE moduleRefParam
USE moduleConstParam USE moduleConstParam
@ -87,7 +87,7 @@ MODULE moduleInject
CLASS(injectGeneric), INTENT(inout):: self CLASS(injectGeneric), INTENT(inout):: self
INTEGER, INTENT(in):: i INTEGER, INTENT(in):: i
REAL(8), INTENT(in):: v, n(1:3), T(1:3) REAL(8), INTENT(in):: v, n(1:3), temperature(1:3)
INTEGER, INTENT(in):: sp, physicalSurface, particlesPerEdge INTEGER, INTENT(in):: sp, physicalSurface, particlesPerEdge
REAL(8):: tauInject REAL(8):: tauInject
REAL(8), INTENT(in):: flow REAL(8), INTENT(in):: flow
@ -97,10 +97,10 @@ MODULE moduleInject
INTEGER:: nVolColl INTEGER:: nVolColl
REAL(8):: fluxPerStep = 0.D0 REAL(8):: fluxPerStep = 0.D0
self%id = i self%id = i
self%vMod = v / v_ref self%vMod = v / v_ref
self%n = n / NORM2(n) self%n = n / NORM2(n)
self%T = T / T_ref self%temperature = temperature / T_ref
!Gets the edge elements from which particles are injected !Gets the edge elements from which particles are injected
DO e = 1, mesh%numEdges DO e = 1, mesh%numEdges
phSurface(e) = mesh%edges(e)%obj%physicalSurface phSurface(e) = mesh%edges(e)%obj%physicalSurface
@ -232,23 +232,23 @@ MODULE moduleInject
END SUBROUTINE doInjects END SUBROUTINE doInjects
SUBROUTINE initVelDistMaxwellian(velDist, T, m) SUBROUTINE initVelDistMaxwellian(velDist, temperature, m)
IMPLICIT NONE IMPLICIT NONE
CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
REAL(8), INTENT(in):: T, m REAL(8), INTENT(in):: temperature, m
velDist = velDistMaxwellian(vTh = DSQRT(T/m)) velDist = velDistMaxwellian(vTh = DSQRT(temperature/m))
END SUBROUTINE initVelDistMaxwellian END SUBROUTINE initVelDistMaxwellian
SUBROUTINE initVelDistHalfMaxwellian(velDist, T, m) SUBROUTINE initVelDistHalfMaxwellian(velDist, temperature, m)
IMPLICIT NONE IMPLICIT NONE
CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
REAL(8), INTENT(in):: T, m REAL(8), INTENT(in):: temperature, m
velDist = velDistHalfMaxwellian(vTh = DSQRT(T/m)) velDist = velDistHalfMaxwellian(vTh = DSQRT(temperature/m))
END SUBROUTINE initVelDistHalfMaxwellian END SUBROUTINE initVelDistHalfMaxwellian

28
src/modules/moduleProbe.f90
View file

@ -27,7 +27,7 @@ MODULE moduleProbe
CONTAINS CONTAINS
!Functions for probeDistFunc type !Functions for probeDistFunc type
SUBROUTINE init(self, id, speciesName, r, v1, v2, v3, points, timeStep) SUBROUTINE init(self, id, speciesName, r, v1, v2, v3, points, everyTimeStep)
USE moduleCaseParam USE moduleCaseParam
USE moduleRefParam USE moduleRefParam
USE moduleSpecies USE moduleSpecies
@ -41,7 +41,7 @@ MODULE moduleProbe
REAL(8), INTENT(in):: r(1:3) REAL(8), INTENT(in):: r(1:3)
REAL(8), INTENT(in):: v1(1:2), v2(1:2), v3(1:2) REAL(8), INTENT(in):: v1(1:2), v2(1:2), v3(1:2)
INTEGER, INTENT(in):: points(1:3) INTEGER, INTENT(in):: points(1:3)
REAL(8), INTENT(in):: timeStep REAL(8), INTENT(in):: everyTimeStep
INTEGER:: sp, i INTEGER:: sp, i
REAL(8):: dv(1:3) REAL(8):: dv(1:3)
@ -91,11 +91,11 @@ MODULE moduleProbe
1:self%nv(3))) 1:self%nv(3)))
!Number of iterations between output !Number of iterations between output
IF (timeStep == 0.D0) THEN IF (everyTimeStep == 0.D0) THEN
self%every = 1 self%every = 1
ELSE ELSE
self%every = NINT(timeStep/ tauMin / ti_ref) self%every = NINT(everyTimeStep/ tauMin / ti_ref)
END IF END IF
@ -189,13 +189,13 @@ MODULE moduleProbe
END SUBROUTINE calculate END SUBROUTINE calculate
SUBROUTINE output(self, t) SUBROUTINE output(self)
USE moduleOutput USE moduleOutput
USE moduleRefParam USE moduleRefParam
USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE IMPLICIT NONE
CLASS(probeDistFunc), INTENT(inout):: self CLASS(probeDistFunc), INTENT(inout):: self
INTEGER, INTENT(in):: t
CHARACTER (LEN=iterationDigits):: tstring CHARACTER (LEN=iterationDigits):: tstring
CHARACTER (LEN=3):: pstring CHARACTER (LEN=3):: pstring
CHARACTER(:), ALLOCATABLE:: filename CHARACTER(:), ALLOCATABLE:: filename
@ -204,14 +204,14 @@ MODULE moduleProbe
!Divide by the velocity cube volume !Divide by the velocity cube volume
self%f = self%f * self%dvInv self%f = self%f * self%dvInv
WRITE(tstring, iterationFormat) t WRITE(tstring, iterationFormat) timeStep
WRITE(pstring, "(I3.3)") self%id WRITE(pstring, "(I3.3)") self%id
fileName='Probe_' // tstring// '_f_' // pstring // '.dat' fileName='Probe_' // tstring// '_f_' // pstring // '.dat'
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (10, file = path // folder // '/' // fileName) OPEN (10, file = path // folder // '/' // fileName)
WRITE(10, "(A1, 1X, A)") "# ", self%species%name WRITE(10, "(A1, 1X, A)") "# ", self%species%name
WRITE(10, "(A6, 3(ES15.6E3), A2)") "# r = ", self%r(:)*L_ref, " m" WRITE(10, "(A6, 3(ES15.6E3), A2)") "# r = ", self%r(:)*L_ref, " m"
WRITE(10, "(A6, ES15.6E3, A2)") "# t = ", REAL(t)*tauMin*ti_ref, " s" WRITE(10, "(A6, ES15.6E3, A2)") "# t = ", REAL(timeStep)*tauMin*ti_ref, " s"
WRITE(10, "(A1, A19, 3(A20))") "#", "v1 (m s^-1)", "v2 (m s^-1)", "v3 (m s^-1)", "f" WRITE(10, "(A1, A19, 3(A20))") "#", "v1 (m s^-1)", "v2 (m s^-1)", "v3 (m s^-1)", "f"
DO i = 1, self%nv(1) DO i = 1, self%nv(1)
DO j = 1, self%nv(2) DO j = 1, self%nv(2)
@ -252,15 +252,15 @@ MODULE moduleProbe
END SUBROUTINE doProbes END SUBROUTINE doProbes
SUBROUTINE outputProbes(t) SUBROUTINE outputProbes()
USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE IMPLICIT NONE
INTEGER, INTENT(in):: t
INTEGER:: i INTEGER:: i
DO i = 1, nProbes DO i = 1, nProbes
IF (probe(i)%update) THEN IF (probe(i)%update) THEN
CALL probe(i)%output(t) CALL probe(i)%output()
END IF END IF
@ -268,15 +268,15 @@ MODULE moduleProbe
END SUBROUTINE outputProbes END SUBROUTINE outputProbes
SUBROUTINE resetProbes(t) SUBROUTINE resetProbes()
USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE IMPLICIT NONE
INTEGER, INTENT(in):: t
INTEGER:: i INTEGER:: i
DO i = 1, nProbes DO i = 1, nProbes
probe(i)%f = 0.D0 probe(i)%f = 0.D0
probe(i)%update = t == tFinal .OR. t == tInitial .OR. MOD(t, probe(i)%every) == 0 probe(i)%update = timeStep == tFinal .OR. timeStep == tInitial .OR. MOD(timeStep, probe(i)%every) == 0
END DO END DO

6
src/modules/output/moduleOutput.f90
View file

@ -160,12 +160,12 @@ MODULE moduleOutput
END SUBROUTINE calculateOutput END SUBROUTINE calculateOutput
SUBROUTINE printTime(t, first) SUBROUTINE printTime(first)
USE moduleSpecies USE moduleSpecies
USE moduleCompTime USE moduleCompTime
USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE IMPLICIT NONE
INTEGER, INTENT(in):: t
LOGICAL, INTENT(in), OPTIONAL:: first LOGICAL, INTENT(in), OPTIONAL:: first
CHARACTER(:), ALLOCATABLE:: fileName CHARACTER(:), ALLOCATABLE:: fileName
@ -187,7 +187,7 @@ MODULE moduleOutput
OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write') OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write')
WRITE (20, "(I10, I10, 7(ES20.6E3))") t, nPartOld, tStep, tPush, tReset, tColl, tCoul, tWeight, tEMField WRITE (20, "(I10, I10, 7(ES20.6E3))") timeStep, nPartOld, tStep, tPush, tReset, tColl, tCoul, tWeight, tEMField
CLOSE(20) CLOSE(20)

180
src/modules/solver/electromagnetic/moduleEM.f90
View file

@ -1,6 +1,7 @@
!Module to solve the electromagnetic field !Module to solve the electromagnetic field
MODULE moduleEM MODULE moduleEM
USE moduleMesh USE moduleMesh
USE moduleTable
IMPLICIT NONE IMPLICIT NONE
! Array of pointers to nodes. ! Array of pointers to nodes.
@ -11,57 +12,155 @@ MODULE moduleEM
END TYPE meshNodePointer END TYPE meshNodePointer
! Generic type for electromagnetic boundary conditions
! TODO: Make this a derived type. TYPE, PUBLIC, ABSTRACT:: boundaryEMGeneric
TYPE:: boundaryEM INTEGER:: nNodes
CHARACTER(:), ALLOCATABLE:: typeEM
INTEGER:: physicalSurface
REAL(8):: potential
TYPE(meshNodePointer), ALLOCATABLE:: nodes(:) TYPE(meshNodePointer), ALLOCATABLE:: nodes(:)
CONTAINS CONTAINS
PROCEDURE, PASS:: apply PROCEDURE(applyEM_interface), DEFERRED, PASS:: apply
!PROCEDURE, PASS:: update !only for time dependent boundary conditions or maybe change apply????? That might be better. !PROCEDURE, PASS:: update !only for time dependent boundary conditions or maybe change apply????? That might be better.
END TYPE boundaryEM END TYPE boundaryEMGeneric
ABSTRACT INTERFACE
! Apply boundary condition to the load vector for the Poission equation
SUBROUTINE applyEM_interface(self, vectorF)
IMPORT boundaryEMGeneric
CLASS(boundaryEMGeneric), INTENT(in):: self
REAL(8), INTENT(inout):: vectorF(:)
END SUBROUTINE applyEM_interface
END INTERFACE
TYPE, EXTENDS(boundaryEMGeneric):: boundaryEMDirichlet
REAL(8):: potential
CONTAINS
! boundaryEMGeneric DEFERRED PROCEDURES
PROCEDURE, PASS:: apply => applyDirichlet
END TYPE boundaryEMDirichlet
TYPE, EXTENDS(boundaryEMGeneric):: boundaryEMDirichletTime
REAL(8):: potential
TYPE(table1D):: temporalProfile
CONTAINS
! boundaryEMGeneric DEFERRED PROCEDURES
PROCEDURE, PASS:: apply => applyDirichletTime
END TYPE boundaryEMDirichletTime
! Container for boundary conditions
TYPE:: boundaryEMCont
CLASS(boundaryEMGeneric), ALLOCATABLE:: obj
END TYPE boundaryEMCont
INTEGER:: nBoundaryEM INTEGER:: nBoundaryEM
TYPE(boundaryEM), ALLOCATABLE:: boundEM(:) TYPE(boundaryEMCont), ALLOCATABLE:: boundaryEM(:)
!Information of charge and reference parameters for rho vector !Information of charge and reference parameters for rho vector
REAL(8), ALLOCATABLE:: qSpecies(:) REAL(8), ALLOCATABLE:: qSpecies(:)
CONTAINS CONTAINS
!Apply boundary conditions to the K matrix for Poisson's equation ! Initialize Dirichlet boundary condition
SUBROUTINE apply(self, edge) SUBROUTINE initDirichlet(self, physicalSurface, potential)
USE moduleMesh
USE moduleRefParam, ONLY: Volt_ref
IMPLICIT NONE
CLASS(boundaryEMGeneric), ALLOCATABLE, INTENT(out):: self
INTEGER:: physicalSurface
REAL(8), INTENT(in):: potential
CLASS(meshEdge), POINTER:: edge
INTEGER, ALLOCATABLE:: nodes(:), nodesEdge(:)
INTEGER:: nNodes, nodesNew
INTEGER:: e, n
! Allocate boundary edge
ALLOCATE(boundaryEMDirichlet:: self)
SELECT TYPE(self)
TYPE IS(boundaryEMDirichlet)
self%potential = potential / Volt_ref
!TODO: This is going into a function
!Temporal array to hold nodes
ALLOCATE(nodes(0))
! Loop thorugh the edges and identify those that are part of the boundary
DO e = 1, mesh%numEdges
edge => mesh%edges(e)%obj
IF (edge%physicalSurface == physicalSurface) THEN
! Edge is of the right boundary index
! Get nodes in the edge
nNodes = edge%nNodes
nodesEdge = edge%getNodes(nNodes)
! Collect all nodes that are not already in the temporal array
DO n = 1, nNodes
IF (ANY(nodes == nodesEdge(n))) THEN
! Node already in array, skip
CYCLE
ELSE
! If not, add element to array of nodes
nodes = [nodes, nodesEdge(n)]
END IF
END DO
END IF
END DO
! Point boundary to nodes
nNodes = SIZE(nodes)
ALLOCATE(self%nodes(nNodes))
DO n = 1, nNodes
self%nodes(n)%obj => mesh%nodes(nodes(n))%obj
END DO
END SELECT
END SUBROUTINE initDirichlet
!Apply Dirichlet boundary condition to the poisson equation
SUBROUTINE applyDirichlet(self, vectorF)
USE moduleMesh USE moduleMesh
IMPLICIT NONE IMPLICIT NONE
CLASS(boundaryEM), INTENT(in):: self CLASS(boundaryEMDirichlet), INTENT(in):: self
CLASS(meshEdge):: edge REAL(8), INTENT(inout):: vectorF(:)
INTEGER:: nNodes INTEGER:: n, ni
INTEGER, ALLOCATABLE:: nodes(:)
INTEGER:: n
nNodes = 1 DO n = 1, self%nNodes
nNodes = edge%nNodes self%nodes(n)%obj%emData%phi = self%potential
nodes = edge%getNodes(nNodes)
DO n = 1, nNodes
SELECT CASE(self%typeEM)
CASE ("dirichlet")
mesh%K(nodes(n), :) = 0.D0
mesh%K(nodes(n), nodes(n)) = 1.D0
! TODO: Change this to pointer
mesh%nodes(nodes(n))%obj%emData%type = self%typeEM
mesh%nodes(nodes(n))%obj%emData%phi = self%potential
END SELECT
END DO END DO
END SUBROUTINE apply END SUBROUTINE applyDirichlet
!Apply Dirichlet boundary condition with time temporal profile
SUBROUTINE applyDirichletTime(self, vectorF)
USE moduleMesh
USE moduleCaseParam, ONLY: timeStep, tauMin
IMPLICIT NONE
CLASS(boundaryEMDirichletTime), INTENT(in):: self
REAL(8), INTENT(inout):: vectorF(:)
INTEGER:: n, ni
DO n = 1, self%nNodes
self%nodes(n)%obj%emData%phi = self%potential
END DO
END SUBROUTINE applyDirichletTime
!Assemble the source vector based on the charge density to solve Poisson's equation !Assemble the source vector based on the charge density to solve Poisson's equation
SUBROUTINE assembleSourceVector(vectorF) SUBROUTINE assembleSourceVector(vectorF)
@ -74,13 +173,14 @@ MODULE moduleEM
INTEGER, ALLOCATABLE:: nodes(:) INTEGER, ALLOCATABLE:: nodes(:)
REAL(8), ALLOCATABLE:: rho(:) REAL(8), ALLOCATABLE:: rho(:)
INTEGER:: nNodes INTEGER:: nNodes
INTEGER:: e, i, ni INTEGER:: e, i, ni, b
CLASS(meshNode), POINTER:: node CLASS(meshNode), POINTER:: node
!$OMP SINGLE !$OMP SINGLE
vectorF = 0.D0 vectorF = 0.D0
!$OMP END SINGLE !$OMP END SINGLE
! Calculate charge density in each node
!$OMP DO REDUCTION(+:vectorF) !$OMP DO REDUCTION(+:vectorF)
DO e = 1, mesh%numCells DO e = 1, mesh%numCells
nNodes = mesh%cells(e)%obj%nNodes nNodes = mesh%cells(e)%obj%nNodes
@ -112,18 +212,12 @@ MODULE moduleEM
!$OMP END DO !$OMP END DO
!Apply boundary conditions !Apply boundary conditions
!$OMP DO !$OMP SINGLE
DO i = 1, mesh%numNodes DO b = 1, nBoundaryEM
node => mesh%nodes(i)%obj CALL boundaryEM(b)%obj%apply(vectorF)
SELECT CASE(node%emData%type)
CASE ("dirichlet")
vectorF(i) = node%emData%phi
END SELECT
END DO END DO
!$OMP END DO !$OMP END SINGLE
END SUBROUTINE assembleSourceVector END SUBROUTINE assembleSourceVector

37
src/modules/solver/moduleSolver.f90
View file

@ -491,47 +491,46 @@ MODULE moduleSolver
END SUBROUTINE updateParticleCell END SUBROUTINE updateParticleCell
!Update the information about if a species needs to be moved this iteration !Update the information about if a species needs to be moved this iteration
SUBROUTINE updatePushSpecies(self, t) SUBROUTINE updatePushSpecies(self)
USE moduleSpecies USE moduleSpecies
USE moduleCaseparam, ONLY: timeStep
IMPLICIT NONE IMPLICIT NONE
CLASS(solverGeneric), INTENT(inout):: self CLASS(solverGeneric), INTENT(inout):: self
INTEGER, INTENT(in):: t
INTEGER:: s INTEGER:: s
DO s=1, nSpecies DO s=1, nSpecies
self%pusher(s)%pushSpecies = MOD(t, self%pusher(s)%every) == 0 self%pusher(s)%pushSpecies = MOD(timeStep, self%pusher(s)%every) == 0
END DO END DO
END SUBROUTINE updatePushSpecies END SUBROUTINE updatePushSpecies
!Output the different data and information !Output the different data and information
SUBROUTINE doOutput(t) SUBROUTINE doOutput()
USE moduleMesh USE moduleMesh
USE moduleOutput USE moduleOutput
USE moduleSpecies USE moduleSpecies
USE moduleCompTime USE moduleCompTime
USE moduleProbe USE moduleProbe
USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE IMPLICIT NONE
INTEGER, INTENT(in):: t CALL outputProbes()
CALL outputProbes(t)
counterOutput = counterOutput + 1 counterOutput = counterOutput + 1
IF (counterOutput >= triggerOutput .OR. & IF (counterOutput >= triggerOutput .OR. &
t == tFinal .OR. t == tInitial) THEN timeStep == tFinal .OR. timeStep == tInitial) THEN
!Resets output counter !Resets output counter
counterOutput=0 counterOutput=0
CALL mesh%printOutput(t) CALL mesh%printOutput()
IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl(t) IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl()
CALL mesh%printEM(t) CALL mesh%printEM()
WRITE(*, "(5X,A21,I10,A1,I10)") "t/tFinal: ", t, "/", tFinal WRITE(*, "(5X,A21,I10,A1,I10)") "t/tFinal: ", timeStep, "/", tFinal
WRITE(*, "(5X,A21,I10)") "Particles: ", nPartOld WRITE(*, "(5X,A21,I10)") "Particles: ", nPartOld
IF (t == 0) THEN IF (timeStep == 0) THEN
WRITE(*, "(5X,A21,F8.1,A2)") " init time: ", 1.D3*tStep, "ms" WRITE(*, "(5X,A21,F8.1,A2)") " init time: ", 1.D3*tStep, "ms"
ELSE ELSE
@ -549,34 +548,32 @@ MODULE moduleSolver
counterCPUTime = counterCPUTime + 1 counterCPUTime = counterCPUTime + 1
IF (counterCPUTime >= triggerCPUTime .OR. & IF (counterCPUTime >= triggerCPUTime .OR. &
t == tFinal .OR. t == tInitial) THEN timeStep == tFinal .OR. timeStep == tInitial) THEN
!Reset CPU Time counter !Reset CPU Time counter
counterCPUTime = 0 counterCPUTime = 0
CALL printTime(t, t == 0) CALL printTime(timeStep == 0)
END IF END IF
!Output average values !Output average values
IF (useAverage .AND. t == tFinal) THEN IF (useAverage .AND. timeStep == tFinal) THEN
CALL mesh%printAverage() CALL mesh%printAverage()
END IF END IF
END SUBROUTINE doOutput END SUBROUTINE doOutput
SUBROUTINE doAverage(t) SUBROUTINE doAverage()
USE moduleAverage USE moduleAverage
USE moduleMesh USE moduleMesh
IMPLICIT NONE IMPLICIT NONE
INTEGER, INTENT(in):: t
INTEGER:: tAverage, n INTEGER:: tAverage, n
IF (useAverage) THEN IF (useAverage) THEN
tAverage = t - tAverageStart tAverage = timeStep - tAverageStart
IF (tAverage == 1) THEN IF (tAverage == 1) THEN
!First iteration in which average scheme is used !First iteration in which average scheme is used