diff --git a/doc/user-manual/fpakc_UserManual.pdf b/doc/user-manual/fpakc_UserManual.pdf
index ef1494a..1666f42 100644
Binary files a/doc/user-manual/fpakc_UserManual.pdf and b/doc/user-manual/fpakc_UserManual.pdf differ
diff --git a/doc/user-manual/fpakc_UserManual.tex b/doc/user-manual/fpakc_UserManual.tex
index 13be4cd..e572331 100644
--- a/doc/user-manual/fpakc_UserManual.tex
+++ b/doc/user-manual/fpakc_UserManual.tex
@@ -585,12 +585,20 @@ make
Type of boundary.
Accepted values are:
\begin{itemize}
- \item \textbf{dirichlet}: Elastic reflection of particles.
+ \item \textbf{dirichlet}: Constant value of electric potential on the surface.
+ \item \textbf{dirichletTime}: Constant value of the electric potential with a time variable profile.
+ The value of \textbf{boundaryEM.potential} will be multiplied for the corresponding value in the file \textbf{boundaryEM.temporalProfile}.
\end{itemize}
- \item \textbf{potential}: Real.
- Fixed potential for Dirichlet boundary condition.
+ \item \textbf{potential}: Real.
+ Fixed potential for Dirichlet boundary condition.
\item \textbf{physicalSurface}: Integer.
Identification of the edge in the mesh file.
+ \item \textbf{temporalProfile}: Character.
+ Filename of the 2 column file containing the time variable profile.
+ File must be located in \textbf{output.path}.
+ The first column is the time in $\unit{s}$.
+ The second column is the factor that will multiply the value of the boundary.
+
\end{itemize}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -611,7 +619,7 @@ make
\begin{itemize}
\item \textbf{A}: Ampere.
\item \textbf{Am2}: Ampere per square meter.
- This value will be multiplied by the surface of injection.
+ This value will be multiplied by the area of injection.
\item \textbf{sccm}: Standard cubic centimetre.
\item \textbf{part/s}: Particles (real) per second.
\end{itemize}
@@ -717,7 +725,7 @@ make
Output file from previous run used as an initial state for the species.
The file format must be the same as in \textbf{geometry.meshType}
Initial particles are assumed to have a Maxwellian distribution.
- File must be located at \textbf{output.path}.
+ File must be located in \textbf{output.path}.
\item \textbf{particlesPerCell}: Integer.
Optional.
Initial number of particles per cell.
diff --git a/src/fpakc.f90 b/src/fpakc.f90
index e90a5e0..ae6e7bb 100644
--- a/src/fpakc.f90
+++ b/src/fpakc.f90
@@ -11,12 +11,12 @@ PROGRAM fpakc
USE OMP_LIB
IMPLICIT NONE
- ! t = time step
- INTEGER:: t
! arg1 = Input argument 1 (input file)
CHARACTER(200):: arg1
! inputFile = path+name of input file
CHARACTER(:), ALLOCATABLE:: inputFile
+ ! generic integer for time step
+ INTEGER:: t
tStep = omp_get_wtime()
!Gets the input file
@@ -32,10 +32,13 @@ PROGRAM fpakc
CALL initOutput(inputFile)
!Do '0' iteration
- t = tInitial
+ timeStep = tInitial
!$OMP PARALLEL DEFAULT(SHARED)
!$OMP SINGLE
+ ! Initial reset of probes
+ CALL resetProbes()
+
CALL verboseError("Initial scatter of particles...")
!$OMP END SINGLE
CALL doScatter()
@@ -49,19 +52,21 @@ PROGRAM fpakc
tStep = omp_get_wtime() - tStep
!Output initial state
- CALL doOutput(t)
+ CALL doOutput()
CALL verboseError('Starting main loop...')
!$OMP PARALLEL DEFAULT(SHARED)
DO t = tInitial + 1, tFinal
- !Insert new particles and push them
!$OMP SINGLE
tStep = omp_get_wtime()
+ ! Update global time step index
+ timeStep = t
+
!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
- CALL resetProbes(t)
+ CALL resetProbes()
tPush = omp_get_wtime()
!$OMP END SINGLE
@@ -79,7 +84,7 @@ PROGRAM fpakc
!$OMP END SINGLE
IF (doMCCollisions) THEN
- CALL meshForMCC%doCollisions(t)
+ CALL meshForMCC%doCollisions()
END IF
@@ -124,12 +129,12 @@ PROGRAM fpakc
!$OMP SINGLE
tEMField = omp_get_wtime() - tEMField
- CALL doAverage(t)
+ CALL doAverage()
tStep = omp_get_wtime() - tStep
!Output data
- CALL doOutput(t)
+ CALL doOutput()
!$OMP END SINGLE
END DO
diff --git a/src/modules/common/moduleCaseParam.f90 b/src/modules/common/moduleCaseParam.f90
index 8df3210..551d867 100644
--- a/src/modules/common/moduleCaseParam.f90
+++ b/src/modules/common/moduleCaseParam.f90
@@ -2,8 +2,13 @@
MODULE moduleCaseParam
!Final and initial iterations
INTEGER:: tFinal, tInitial = 0
+ ! Global index of current iteration
+ INTEGER:: timeStep
+ ! Time step for all species
REAL(8), ALLOCATABLE:: tau(:)
+ ! Minimum time step
REAL(8):: tauMin
+ ! Time step for Monte-Carlo Collisions
REAL(8):: tauColl
END MODULE moduleCaseParam
diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90
index 6f1d5bb..b5f6881 100644
--- a/src/modules/init/moduleInput.f90
+++ b/src/modules/init/moduleInput.f90
@@ -264,8 +264,8 @@ MODULE moduleInput
CALL readEMBoundary(config)
!Read constant magnetic field
DO i = 1, 3
- WRITE(istring, '(i2)') i
- CALL config%get(object // '.B(' // istring // ')', B(i), found)
+ WRITE(iString, '(i2)') i
+ CALL config%get(object // '.B(' // iString // ')', B(i), found)
IF (.NOT. found) THEN
CALL criticalError('Constant magnetic field not provided in direction ' // iString, 'readSolver')
@@ -799,7 +799,7 @@ MODULE moduleInput
TYPE(json_file), INTENT(inout):: config
INTEGER:: i, s
- CHARACTER(2):: istring, sString
+ CHARACTER(2):: iString, sString
CHARACTER(:), ALLOCATABLE:: object, bType
REAL(8):: Tw, cw !Wall temperature and specific heat
!Neutral Properties
@@ -815,8 +815,8 @@ MODULE moduleInput
CALL config%info('boundary', found, n_children = nBoundary)
ALLOCATE(boundary(1:nBoundary))
DO i = 1, nBoundary
- WRITE(istring, '(i2)') i
- object = 'boundary(' // TRIM(istring) // ')'
+ WRITE(iString, '(i2)') i
+ object = 'boundary(' // TRIM(iString) // ')'
boundary(i)%n = i
CALL config%get(object // '.name', boundary(i)%name, found)
@@ -1100,13 +1100,13 @@ MODULE moduleInput
TYPE(json_file), INTENT(inout):: config
CHARACTER(:), ALLOCATABLE:: object
LOGICAL:: found
- CHARACTER(2):: istring
+ CHARACTER(2):: iString
INTEGER:: i
CHARACTER(:), ALLOCATABLE:: speciesName
REAL(8), ALLOCATABLE, DIMENSION(:):: r
REAL(8), ALLOCATABLE, DIMENSION(:):: v1, v2, v3
INTEGER, ALLOCATABLE, DIMENSION(:):: points
- REAL(8):: timeStep
+ REAL(8):: everyTimeStep
CALL config%info('output.probes', found, n_children = nProbes)
@@ -1114,7 +1114,7 @@ MODULE moduleInput
DO i = 1, nProbes
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 // '.position', r, found)
@@ -1122,16 +1122,14 @@ MODULE moduleInput
CALL config%get(object // '.velocity_2', v2, found)
CALL config%get(object // '.velocity_3', v3, found)
CALL config%get(object // '.points', points, found)
- CALL config%get(object // '.timeStep', timeStep, found)
+ CALL config%get(object // '.timeStep', everyTimeStep, found)
IF (ANY(points < 2)) CALL criticalError("Number of points in probe " // iString // " incorrect", 'readProbes')
- CALL probe(i)%init(i, speciesName, r, v1, v2, v3, points, timeStep)
+ CALL probe(i)%init(i, speciesName, r, v1, v2, v3, points, everyTimeStep)
END DO
- CALL resetProbes(tInitial)
-
END SUBROUTINE readProbes
SUBROUTINE readEMBoundary(config)
@@ -1139,7 +1137,6 @@ MODULE moduleInput
USE moduleOutput
USE moduleErrors
USE moduleEM
- USE moduleRefParam
USE moduleSpecies
USE json_module
IMPLICIT NONE
@@ -1147,34 +1144,72 @@ MODULE moduleInput
TYPE(json_file), INTENT(inout):: config
CHARACTER(:), ALLOCATABLE:: object
LOGICAL:: found
- CHARACTER(2):: istring
- INTEGER:: i, e, s
+ CHARACTER(:), ALLOCATABLE:: typeEM
+ REAL(8):: potential
+ INTEGER:: physicalSurface
+ CHARACTER(:), ALLOCATABLE:: temporalProfile, temporalProfilePath
+ INTEGER:: b, s, n, ni
+ CHARACTER(2):: bString
INTEGER:: info
EXTERNAL:: dgetrf
CALL config%info('boundaryEM', found, n_children = nBoundaryEM)
- IF (found) ALLOCATE(boundEM(1:nBoundaryEM))
+ IF (found) THEN
+ ALLOCATE(boundaryEM(1:nBoundaryEM))
+
+ END IF
- DO i = 1, nBoundaryEM
- WRITE(istring, '(I2)') i
- object = 'boundaryEM(' // trim(istring) // ')'
+ DO b = 1, nBoundaryEM
+ WRITE(bString, '(I2)') b
+ 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")
- CALL config%get(object // '.potential', boundEM(i)%potential, found)
- IF (.NOT. found) &
+ CALL config%get(object // '.potential', potential, found)
+ IF (.NOT. found) THEN
CALL criticalError('Required parameter "potential" for Dirichlet boundary condition not found', 'readEMBoundary')
- boundEM(i)%potential = boundEM(i)%potential/Volt_ref
- CALL config%get(object // '.physicalSurface', boundEM(i)%physicalSurface, found)
- IF (.NOT. found) &
- CALL criticalError('Required parameter "physicalSurface" for Dirichlet boundary condition not found', 'readEMBoundary')
+ END IF
+
+ CALL config%get(object // '.physicalSurface', physicalSurface, found)
+ IF (.NOT. found) THEN
+ CALL criticalError('Required parameter "physicalSurface" for Dirichlet boundary condition not found', &
+ 'readEMBoundary')
+
+ END IF
+
+ CALL initDirichlet(boundaryEM(b)%obj, physicalSurface, potential)
+
+ CASE ("dirichletTime")
+ CALL config%get(object // '.potential', potential, found)
+ IF (.NOT. found) THEN
+ CALL criticalError('Required parameter "potential" for Dirichlet Time boundary condition not found', &
+ 'readEMBoundary')
+
+ END IF
+
+ CALL config%get(object // '.temporalProfile', temporalProfile, found)
+ IF (.NOT. found) THEN
+ CALL criticalError('Required parameter "temporalProfile" for Dirichlet Time boundary condition not found', &
+ 'readEMBoundary')
+
+ END IF
+ temporalProfilePath = path // temporalProfile
+
+ CALL config%get(object // '.physicalSurface', physicalSurface, found)
+ IF (.NOT. found) THEN
+ CALL criticalError('Required parameter "physicalSurface" for Dirichlet Time boundary condition not found', &
+ 'readEMBoundary')
+
+ END IF
+
+ CALL initDirichletTime(boundaryEM(b)%obj, physicalSurface, potential, temporalProfilePath)
CASE DEFAULT
- CALL criticalError('Boundary type ' // boundEM(i)%typeEM // ' not yet supported', 'readEMBoundary')
+ CALL criticalError('Boundary type ' // typeEM // ' not yet supported', 'readEMBoundary')
END SELECT
@@ -1193,18 +1228,28 @@ MODULE moduleInput
END DO
- IF (ALLOCATED(boundEM)) THEN
- DO e = 1, mesh%numEdges
- IF (ANY(mesh%edges(e)%obj%physicalSurface == boundEM(:)%physicalSurface)) THEN
- DO i = 1, nBoundaryEM
- IF (mesh%edges(e)%obj%physicalSurface == boundEM(i)%physicalSurface) THEN
- CALL boundEM(i)%apply(mesh%edges(e)%obj)
+ ! Modify K matrix due to boundary conditions
+ DO b = 1, nBoundaryEM
+ SELECT TYPE(boundary => boundaryEM(b)%obj)
+ TYPE IS(boundaryEMDirichlet)
+ DO n = 1, boundary%nNodes
+ ni = boundary%nodes(n)%obj%n
+ mesh%K(ni, :) = 0.D0
+ mesh%K(ni, ni) = 1.D0
- END IF
- END DO
- END IF
- END DO
- END IF
+ END DO
+
+ TYPE IS(boundaryEMDirichletTime)
+ DO n = 1, boundary%nNodes
+ ni = boundary%nodes(n)%obj%n
+ mesh%K(ni, :) = 0.D0
+ mesh%K(ni, ni) = 1.D0
+
+ END DO
+
+ END SELECT
+
+ END DO
!Compute the PLU factorization of K once boundary conditions have been read
CALL dgetrf(mesh%numNodes, mesh%numNodes, mesh%K, mesh%numNodes, mesh%IPIV, info)
@@ -1225,13 +1270,13 @@ MODULE moduleInput
TYPE(json_file), INTENT(inout):: config
INTEGER:: i
- CHARACTER(2):: istring
+ CHARACTER(2):: iString
CHARACTER(:), ALLOCATABLE:: object
LOGICAL:: found
CHARACTER(:), ALLOCATABLE:: speciesName
CHARACTER(:), ALLOCATABLE:: name
REAL(8):: v
- REAL(8), ALLOCATABLE:: T(:), normal(:)
+ REAL(8), ALLOCATABLE:: temperature(:), normal(:)
REAL(8):: flow
CHARACTER(:), ALLOCATABLE:: units
INTEGER:: physicalSurface
@@ -1242,8 +1287,8 @@ MODULE moduleInput
ALLOCATE(inject(1:nInject))
nPartInj = 0
DO i = 1, nInject
- WRITE(istring, '(i2)') i
- object = 'inject(' // trim(istring) // ')'
+ WRITE(iString, '(i2)') i
+ object = 'inject(' // trim(iString) // ')'
!Find species
CALL config%get(object // '.species', speciesName, found)
@@ -1251,7 +1296,7 @@ MODULE moduleInput
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 // '.T', temperature, found)
CALL config%get(object // '.n', normal, found)
IF (.NOT. found) THEN
ALLOCATE(normal(1:3))
@@ -1263,7 +1308,7 @@ MODULE moduleInput
particlesPerEdge = 0
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)
@@ -1323,28 +1368,28 @@ MODULE moduleInput
TYPE(injectGeneric), INTENT(inout):: inj
CHARACTER(:), ALLOCATABLE, INTENT(in):: object
INTEGER:: i
- CHARACTER(2):: istring
+ CHARACTER(2):: iString
CHARACTER(:), ALLOCATABLE:: fvType
LOGICAL:: found
- REAL(8):: v, T, m
+ REAL(8):: v, temperature, m
!Reads species mass
m = inj%species%m
!Reads distribution functions for velocity
DO i = 1, 3
- WRITE(istring, '(i2)') i
- CALL config%get(object // '.velDist('// TRIM(istring) //')', fvType, found)
- IF (.NOT. found) CALL criticalError("No velocity distribution in direction " // istring // &
+ WRITE(iString, '(i2)') i
+ CALL config%get(object // '.velDist('// TRIM(iString) //')', fvType, found)
+ IF (.NOT. found) CALL criticalError("No velocity distribution in direction " // iString // &
" found for " // object, 'readVelDistr')
SELECT CASE(fvType)
CASE ("Maxwellian")
- T = inj%T(i)
- CALL initVelDistMaxwellian(inj%v(i)%obj, t, m)
+ temperature = inj%temperature(i)
+ CALL initVelDistMaxwellian(inj%v(i)%obj, temperature, m)
CASE ("Half-Maxwellian")
- T = inj%T(i)
- CALL initVelDistHalfMaxwellian(inj%v(i)%obj, t, m)
+ temperature = inj%temperature(i)
+ CALL initVelDistHalfMaxwellian(inj%v(i)%obj, temperature, m)
CASE ("Delta")
v = inj%vMod*inj%n(i)
diff --git a/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 b/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90
index 97ec729..c0dcfbb 100644
--- a/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90
+++ b/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90
@@ -1,22 +1,22 @@
MODULE moduleMeshOutput0D
CONTAINS
- SUBROUTINE printOutput0D(self, t)
+ SUBROUTINE printOutput0D(self)
USE moduleMesh
USE moduleRefParam
USE moduleSpecies
USE moduleOutput
+ USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self
- INTEGER, INTENT(in):: t
INTEGER:: i
TYPE(outputFormat):: output
CHARACTER(:), ALLOCATABLE:: fileName
DO i = 1, nSpecies
fileName='OUTPUT_' // species(i)%obj%name // '.dat'
- IF (t == 0) THEN
+ IF (timeStep == 0) THEN
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)", &
"pressure (Pa)", "temperature (K)"
@@ -27,14 +27,17 @@ MODULE moduleMeshOutput0D
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)
- 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)
END DO
END SUBROUTINE printOutput0D
- SUBROUTINE printColl0D(self, t)
+ SUBROUTINE printColl0D(self)
USE moduleMesh
USE moduleRefParam
USE moduleCaseParam
@@ -43,12 +46,11 @@ MODULE moduleMeshOutput0D
IMPLICIT NONE
CLASS(meshGeneric), INTENT(in):: self
- INTEGER, INTENT(in):: t
CHARACTER(:), ALLOCATABLE:: fileName
INTEGER:: k
fileName='OUTPUT_Collisions.dat'
- IF (t == tInitial) THEN
+ IF (timeStep == tInitial) THEN
OPEN(20, file = path // folder // '/' // fileName, action = 'write')
WRITE(20, "(A1, 14X, A5, A20)") "#","t (s)","collisions"
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
@@ -57,12 +59,12 @@ MODULE moduleMeshOutput0D
END IF
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)
END SUBROUTINE printColl0D
- SUBROUTINE printEM0D(self, t)
+ SUBROUTINE printEM0D(self)
USE moduleMesh
USE moduleRefParam
USE moduleCaseParam
@@ -70,7 +72,6 @@ MODULE moduleMeshOutput0D
IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self
- INTEGER, INTENT(in):: t
END SUBROUTINE printEM0D
diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90
index ccdcf3d..8176bb5 100644
--- a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90
+++ b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90
@@ -80,50 +80,50 @@ MODULE moduleMeshOutputGmsh2
END SUBROUTINE writeGmsh2FooterElementData
!Prints the scattered properties of particles into the nodes
- SUBROUTINE printOutputGmsh2(self, t)
+ SUBROUTINE printOutputGmsh2(self)
USE moduleMesh
USE moduleRefParam
USE moduleSpecies
USE moduleOutput
USE moduleMeshInoutCommon
+ USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self
- INTEGER, INTENT(in):: t
INTEGER:: n, i
TYPE(outputFormat):: output(1:self%numNodes)
REAL(8):: time
CHARACTER(:), ALLOCATABLE:: fileName
- time = DBLE(t)*tauMin*ti_ref
+ time = DBLE(timeStep)*tauMin*ti_ref
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
OPEN (60, file = path // folder // '/' // fileName)
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
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
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
WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%velocity
END DO
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
WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%pressure
END DO
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
WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%temperature
END DO
@@ -135,7 +135,7 @@ MODULE moduleMeshOutputGmsh2
END SUBROUTINE printOutputGmsh2
!Prints the number of collisions into the volumes
- SUBROUTINE printCollGmsh2(self, t)
+ SUBROUTINE printCollGmsh2(self)
USE moduleMesh
USE moduleRefParam
USE moduleCaseParam
@@ -145,7 +145,6 @@ MODULE moduleMeshOutputGmsh2
IMPLICIT NONE
CLASS(meshGeneric), INTENT(in):: self
- INTEGER, INTENT(in):: t
INTEGER:: numEdges
INTEGER:: k, c
INTEGER:: n
@@ -167,9 +166,9 @@ MODULE moduleMeshOutputGmsh2
END SELECT
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
OPEN (60, file = path // folder // '/' // fileName)
@@ -179,7 +178,7 @@ MODULE moduleMeshOutputGmsh2
DO c = 1, interactionMatrix(k)%amount
WRITE(cString, "(I2)") c
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
WRITE(60, "(I6,I10)") n + numEdges, self%cells(n)%obj%tallyColl(k)%tally(c)
END DO
@@ -196,7 +195,7 @@ MODULE moduleMeshOutputGmsh2
END SUBROUTINE printCollGmsh2
!Prints the electrostatic EM properties into the nodes and volumes
- SUBROUTINE printEMGmsh2(self, t)
+ SUBROUTINE printEMGmsh2(self)
USE moduleMesh
USE moduleRefParam
USE moduleCaseParam
@@ -205,7 +204,6 @@ MODULE moduleMeshOutputGmsh2
IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self
- INTEGER, INTENT(in):: t
INTEGER:: n, e
REAL(8):: time
CHARACTER(:), ALLOCATABLE:: fileName
@@ -214,27 +212,27 @@ MODULE moduleMeshOutputGmsh2
Xi = (/ 0.D0, 0.D0, 0.D0 /)
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
OPEN (20, file = path // folder // '/' // fileName)
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
WRITE(20, *) n, self%nodes(n)%obj%emData%phi*Volt_ref
END DO
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
WRITE(20, *) e+self%numEdges, self%cells(e)%obj%gatherElectricField(Xi)*EF_ref
END DO
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
WRITE(20, *) n, self%nodes(n)%obj%emData%B * B_ref
END DO
diff --git a/src/modules/mesh/inout/moduleMeshInoutCommon.f90 b/src/modules/mesh/inout/moduleMeshInoutCommon.f90
index e4a6c72..7dc2e84 100644
--- a/src/modules/mesh/inout/moduleMeshInoutCommon.f90
+++ b/src/modules/mesh/inout/moduleMeshInoutCommon.f90
@@ -3,17 +3,17 @@ MODULE moduleMeshInoutCommon
CHARACTER(LEN=4):: prefix = 'Step'
CONTAINS
- PURE FUNCTION formatFileName(prefix, suffix, extension, t) RESULT(fileName)
+ PURE FUNCTION formatFileName(prefix, suffix, extension, timeStep) RESULT(fileName)
USE moduleOutput
IMPLICIT NONE
CHARACTER(*), INTENT(in):: prefix, suffix, extension
- INTEGER, INTENT(in), OPTIONAL:: t
+ INTEGER, INTENT(in), OPTIONAL:: timeStep
CHARACTER (LEN=iterationDigits):: tString
CHARACTER(:), ALLOCATABLE:: fileName
- IF (PRESENT(t)) THEN
- WRITE(tString, iterationFormat) t
+ IF (PRESENT(timeStep)) THEN
+ WRITE(tString, iterationFormat) timeStep
fileName = prefix // '_' // tString // '_' // suffix // '.' // extension
ELSE
diff --git a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90
index 184d645..e07db01 100644
--- a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90
+++ b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90
@@ -167,7 +167,7 @@ MODULE moduleMeshInputVTU
CLASS(meshGeneric), INTENT(inout):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
REAL(8):: r(1:3) !3 generic coordinates
- INTEGER:: fileID, error, found
+ INTEGER:: fileID
CHARACTER(LEN=256):: line
INTEGER:: numNodes, numElements, numEdges
INTEGER, ALLOCATABLE, DIMENSION(:):: entitiesID, offsets, connectivity, types
diff --git a/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90
index da90b6b..81e4bbf 100644
--- a/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90
+++ b/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90
@@ -215,17 +215,16 @@ MODULE moduleMeshOutputVTU
END SUBROUTINE writeEM
- SUBROUTINE writeCollection(fileID, t, fileNameStep, fileNameCollection)
+ SUBROUTINE writeCollection(fileID, fileNameStep, fileNameCollection)
USE moduleCaseParam
USE moduleOutput
USE moduleRefParam
IMPLICIT NONE
INTEGER:: fileID
- INTEGER, INTENT(in):: t
CHARACTER(*):: fileNameStep, fileNameCollection
- IF (t == tInitial) THEN
+ IF (timeStep == tInitial) THEN
!Create collection file
WRITE(*, "(6X,A15,A)") "Creating file: ", fileNameCollection
OPEN (fileID + 1, file = path // folder // '/' // fileNameCollection)
@@ -237,10 +236,11 @@ MODULE moduleMeshOutputVTU
!Write iteration file in collection
OPEN (fileID + 1, file = path // folder // '/' // fileNameCollection, ACCESS='APPEND')
- WRITE(fileID + 1, "(4X, A, ES20.6E3, A, A, A)") ''
+ WRITE(fileID + 1, "(4X, A, ES20.6E3, A, A, A)") &
+ ''
!Close collection file
- IF (t == tFinal) THEN
+ IF (timeStep == tFinal) THEN
WRITE (fileID + 1, "(2X, A)") ''
WRITE (fileID + 1, "(A)") ''
@@ -307,21 +307,21 @@ MODULE moduleMeshOutputVTU
END SUBROUTINE writeAverage
- SUBROUTINE printOutputVTU(self,t)
+ SUBROUTINE printOutputVTU(self)
USE moduleMesh
USE moduleSpecies
USE moduleMeshInoutCommon
+ USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self
- INTEGER, INTENT(in):: t
INTEGER:: i, fileID
CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection
fileID = 60
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
OPEN (fileID, file = path // folder // '/' // fileName)
@@ -337,28 +337,27 @@ MODULE moduleMeshOutputVTU
!Write collection file for time plotting
fileNameCollection = formatFileName('Collection', species(i)%obj%name, 'pvd')
- CALL writeCollection(fileID, t, fileName, filenameCollection)
+ CALL writeCollection(fileID, fileName, filenameCollection)
END DO
END SUBROUTINE printOutputVTU
- SUBROUTINE printCollVTU(self,t)
+ SUBROUTINE printCollVTU(self)
USE moduleMesh
USE moduleOutput
USE moduleMeshInoutCommon
+ USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE
CLASS(meshGeneric), INTENT(in):: self
- INTEGER, INTENT(in):: t
INTEGER:: fileID
CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection
- CHARACTER (LEN=iterationDigits):: tstring
fileID = 62
IF (collOutput) THEN
- fileName = formatFileName(prefix, 'Collisions', 'vtu', t)
+ fileName = formatFileName(prefix, 'Collisions', 'vtu', timeStep)
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (fileID, file = path // folder // '/' // fileName)
@@ -374,26 +373,26 @@ MODULE moduleMeshOutputVTU
!Write collection file for time plotting
fileNameCollection = formatFileName('Collection', 'Collisions', 'pvd')
- CALL writeCollection(fileID, t, fileName, filenameCollection)
+ CALL writeCollection(fileID, fileName, filenameCollection)
END IF
END SUBROUTINE printCollVTU
- SUBROUTINE printEMVTU(self, t)
+ SUBROUTINE printEMVTU(self)
USE moduleMesh
USE moduleMeshInoutCommon
+ USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self
- INTEGER, INTENT(in):: t
INTEGER:: fileID
CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection
fileID = 64
IF (emOutput) THEN
- fileName = formatFileName(prefix, 'EMField', 'vtu', t)
+ fileName = formatFileName(prefix, 'EMField', 'vtu', timeStep)
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (fileID, file = path // folder // '/' // fileName)
@@ -409,7 +408,7 @@ MODULE moduleMeshOutputVTU
!Write collection file for time plotting
fileNameCollection = formatFileName('Collection', 'EMField', 'pvd')
- CALL writeCollection(fileID, t, fileName, filenameCollection)
+ CALL writeCollection(fileID, fileName, filenameCollection)
END IF
diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90
index ae53aa0..7ab3914 100644
--- a/src/modules/mesh/moduleMesh.f90
+++ b/src/modules/mesh/moduleMesh.f90
@@ -59,6 +59,13 @@ MODULE moduleMesh
END TYPE meshNodeCont
+ ! Array of pointers to nodes.
+ TYPE:: meshNodePointer
+ CLASS(meshNode), POINTER:: obj
+ CONTAINS
+
+ END TYPE meshNodePointer
+
!Type for array of boundary functions (one per species)
TYPE, PUBLIC:: fBoundaryGeneric
PROCEDURE(boundary_interface), POINTER, NOPASS:: apply => NULL()
@@ -372,10 +379,9 @@ MODULE moduleMesh
END SUBROUTINE connectMesh_interface
!Prints number of collisions in each cell
- SUBROUTINE printColl_interface(self, t)
+ SUBROUTINE printColl_interface(self)
IMPORT meshGeneric
CLASS(meshGeneric), INTENT(in):: self
- INTEGER, INTENT(in):: t
END SUBROUTINE printColl_interface
@@ -403,18 +409,16 @@ MODULE moduleMesh
ABSTRACT INTERFACE
!Prints Species data
- SUBROUTINE printOutput_interface(self, t)
+ SUBROUTINE printOutput_interface(self)
IMPORT meshParticles
CLASS(meshParticles), INTENT(in):: self
- INTEGER, INTENT(in):: t
END SUBROUTINE printOutput_interface
!Prints EM info
- SUBROUTINE printEM_interface(self, t)
+ SUBROUTINE printEM_interface(self)
IMPORT meshParticles
CLASS(meshParticles), INTENT(in):: self
- INTEGER, INTENT(in):: t
END SUBROUTINE printEM_interface
@@ -789,7 +793,7 @@ MODULE moduleMesh
END FUNCTION findCellBrute
!Computes collisions in element
- SUBROUTINE doCollisions(self, t)
+ SUBROUTINE doCollisions(self)
USE moduleCollisions
USE moduleSpecies
USE moduleList
@@ -797,10 +801,10 @@ MODULE moduleMesh
USE moduleRandom
USE moduleOutput
USE moduleMath
+ USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE
CLASS(meshGeneric), INTENT(inout), TARGET:: self
- INTEGER, INTENT(in):: t
INTEGER:: e
CLASS(meshCell), POINTER:: cell
INTEGER:: k, i, j
@@ -816,7 +820,7 @@ MODULE moduleMesh
REAL(8):: rnd_real !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
!$OMP DO SCHEDULE(DYNAMIC) PRIVATE(part_i, part_j, partTemp_i, partTemp_j)
DO e=1, self%numCells
diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90
index 18c7fbb..5b96be2 100644
--- a/src/modules/moduleInject.f90
+++ b/src/modules/moduleInject.f90
@@ -54,7 +54,7 @@ MODULE moduleInject
INTEGER:: id
CHARACTER(:), ALLOCATABLE:: name
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
LOGICAL:: fixDirection !The injection of particles has a fix direction defined by n
INTEGER:: nParticles !Number of particles to introduce each time step
@@ -76,7 +76,7 @@ MODULE moduleInject
CONTAINS
!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 moduleRefParam
USE moduleConstParam
@@ -87,7 +87,7 @@ MODULE moduleInject
CLASS(injectGeneric), INTENT(inout):: self
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
REAL(8):: tauInject
REAL(8), INTENT(in):: flow
@@ -97,10 +97,10 @@ MODULE moduleInject
INTEGER:: nVolColl
REAL(8):: fluxPerStep = 0.D0
- self%id = i
- self%vMod = v / v_ref
- self%n = n / NORM2(n)
- self%T = T / T_ref
+ self%id = i
+ self%vMod = v / v_ref
+ self%n = n / NORM2(n)
+ self%temperature = temperature / T_ref
!Gets the edge elements from which particles are injected
DO e = 1, mesh%numEdges
phSurface(e) = mesh%edges(e)%obj%physicalSurface
@@ -232,23 +232,23 @@ MODULE moduleInject
END SUBROUTINE doInjects
- SUBROUTINE initVelDistMaxwellian(velDist, T, m)
+ SUBROUTINE initVelDistMaxwellian(velDist, temperature, m)
IMPLICIT NONE
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
- SUBROUTINE initVelDistHalfMaxwellian(velDist, T, m)
+ SUBROUTINE initVelDistHalfMaxwellian(velDist, temperature, m)
IMPLICIT NONE
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
diff --git a/src/modules/moduleProbe.f90 b/src/modules/moduleProbe.f90
index c29eeda..754d56a 100644
--- a/src/modules/moduleProbe.f90
+++ b/src/modules/moduleProbe.f90
@@ -27,7 +27,7 @@ MODULE moduleProbe
CONTAINS
!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 moduleRefParam
USE moduleSpecies
@@ -41,7 +41,7 @@ MODULE moduleProbe
REAL(8), INTENT(in):: r(1:3)
REAL(8), INTENT(in):: v1(1:2), v2(1:2), v3(1:2)
INTEGER, INTENT(in):: points(1:3)
- REAL(8), INTENT(in):: timeStep
+ REAL(8), INTENT(in):: everyTimeStep
INTEGER:: sp, i
REAL(8):: dv(1:3)
@@ -91,11 +91,11 @@ MODULE moduleProbe
1:self%nv(3)))
!Number of iterations between output
- IF (timeStep == 0.D0) THEN
+ IF (everyTimeStep == 0.D0) THEN
self%every = 1
ELSE
- self%every = NINT(timeStep/ tauMin / ti_ref)
+ self%every = NINT(everyTimeStep/ tauMin / ti_ref)
END IF
@@ -189,13 +189,13 @@ MODULE moduleProbe
END SUBROUTINE calculate
- SUBROUTINE output(self, t)
+ SUBROUTINE output(self)
USE moduleOutput
USE moduleRefParam
+ USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE
CLASS(probeDistFunc), INTENT(inout):: self
- INTEGER, INTENT(in):: t
CHARACTER (LEN=iterationDigits):: tstring
CHARACTER (LEN=3):: pstring
CHARACTER(:), ALLOCATABLE:: filename
@@ -204,14 +204,14 @@ MODULE moduleProbe
!Divide by the velocity cube volume
self%f = self%f * self%dvInv
- WRITE(tstring, iterationFormat) t
+ WRITE(tstring, iterationFormat) timeStep
WRITE(pstring, "(I3.3)") self%id
fileName='Probe_' // tstring// '_f_' // pstring // '.dat'
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (10, file = path // folder // '/' // fileName)
WRITE(10, "(A1, 1X, A)") "# ", self%species%name
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"
DO i = 1, self%nv(1)
DO j = 1, self%nv(2)
@@ -252,15 +252,15 @@ MODULE moduleProbe
END SUBROUTINE doProbes
- SUBROUTINE outputProbes(t)
+ SUBROUTINE outputProbes()
+ USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE
- INTEGER, INTENT(in):: t
INTEGER:: i
DO i = 1, nProbes
IF (probe(i)%update) THEN
- CALL probe(i)%output(t)
+ CALL probe(i)%output()
END IF
@@ -268,15 +268,15 @@ MODULE moduleProbe
END SUBROUTINE outputProbes
- SUBROUTINE resetProbes(t)
+ SUBROUTINE resetProbes()
+ USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE
- INTEGER, INTENT(in):: t
INTEGER:: i
DO i = 1, nProbes
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
diff --git a/src/modules/output/moduleOutput.f90 b/src/modules/output/moduleOutput.f90
index 18fbb7f..e6dc91f 100644
--- a/src/modules/output/moduleOutput.f90
+++ b/src/modules/output/moduleOutput.f90
@@ -160,12 +160,12 @@ MODULE moduleOutput
END SUBROUTINE calculateOutput
- SUBROUTINE printTime(t, first)
+ SUBROUTINE printTime(first)
USE moduleSpecies
USE moduleCompTime
+ USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE
- INTEGER, INTENT(in):: t
LOGICAL, INTENT(in), OPTIONAL:: first
CHARACTER(:), ALLOCATABLE:: fileName
@@ -187,7 +187,7 @@ MODULE moduleOutput
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)
diff --git a/src/modules/solver/electromagnetic/moduleEM.f90 b/src/modules/solver/electromagnetic/moduleEM.f90
index d5d0793..126f9a6 100644
--- a/src/modules/solver/electromagnetic/moduleEM.f90
+++ b/src/modules/solver/electromagnetic/moduleEM.f90
@@ -1,53 +1,200 @@
!Module to solve the electromagnetic field
MODULE moduleEM
+ USE moduleMesh
+ USE moduleTable
IMPLICIT NONE
- TYPE:: boundaryEM
- CHARACTER(:), ALLOCATABLE:: typeEM
- INTEGER:: physicalSurface
+ ! Generic type for electromagnetic boundary conditions
+ TYPE, PUBLIC, ABSTRACT:: boundaryEMGeneric
+ INTEGER:: nNodes
+ TYPE(meshNodePointer), ALLOCATABLE:: nodes(:)
+
+ CONTAINS
+ PROCEDURE(applyEM_interface), DEFERRED, PASS:: apply
+ !PROCEDURE, PASS:: update !only for time dependent boundary conditions or maybe change apply????? That might be better.
+
+ 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
- PROCEDURE, PASS:: apply
+ ! boundaryEMGeneric DEFERRED PROCEDURES
+ PROCEDURE, PASS:: apply => applyDirichlet
- END TYPE boundaryEM
+ 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
- TYPE(boundaryEM), ALLOCATABLE:: boundEM(:)
+ TYPE(boundaryEMCont), ALLOCATABLE:: boundaryEM(:)
!Information of charge and reference parameters for rho vector
REAL(8), ALLOCATABLE:: qSpecies(:)
CONTAINS
- !Apply boundary conditions to the K matrix for Poisson's equation
- SUBROUTINE apply(self, edge)
+ SUBROUTINE findNodes(self, physicalSurface)
USE moduleMesh
IMPLICIT NONE
- CLASS(boundaryEM), INTENT(in):: self
- CLASS(meshEdge):: edge
- INTEGER:: nNodes
- INTEGER, ALLOCATABLE:: nodes(:)
- INTEGER:: n
+ CLASS(boundaryEMGeneric), INTENT(inout):: self
+ INTEGER, INTENT(in):: physicalSurface
+ CLASS(meshEdge), POINTER:: edge
+ INTEGER, ALLOCATABLE:: nodes(:), nodesEdge(:)
+ INTEGER:: nNodes, nodesNew
+ INTEGER:: e, n
- nNodes = 1
- nNodes = edge%nNodes
- nodes = edge%getNodes(nNodes)
+ !Temporal array to hold nodes
+ ALLOCATE(nodes(0))
- DO n = 1, nNodes
- SELECT CASE(self%typeEM)
- CASE ("dirichlet")
- mesh%K(nodes(n), :) = 0.D0
- mesh%K(nodes(n), nodes(n)) = 1.D0
-
- mesh%nodes(nodes(n))%obj%emData%type = self%typeEM
- mesh%nodes(nodes(n))%obj%emData%phi = self%potential
+ ! 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
- END SELECT
+ ELSE
+ ! If not, add element to array of nodes
+ nodes = [nodes, nodesEdge(n)]
+
+ END IF
+
+ END DO
+
+ END IF
END DO
- END SUBROUTINE
+ ! Point boundary to nodes
+ nNodes = SIZE(nodes)
+ ALLOCATE(self%nodes(nNodes))
+ self%nNodes = nNodes
+ DO n = 1, nNodes
+ self%nodes(n)%obj => mesh%nodes(nodes(n))%obj
+
+ END DO
+
+ END SUBROUTINE findNodes
+
+ ! Initialize Dirichlet boundary condition
+ SUBROUTINE initDirichlet(self, physicalSurface, potential)
+ USE moduleRefParam, ONLY: Volt_ref
+ IMPLICIT NONE
+
+ CLASS(boundaryEMGeneric), ALLOCATABLE, INTENT(out):: self
+ INTEGER, INTENT(in):: physicalSurface
+ REAL(8), INTENT(in):: potential
+
+ ! Allocate boundary edge
+ ALLOCATE(boundaryEMDirichlet:: self)
+
+ SELECT TYPE(self)
+ TYPE IS(boundaryEMDirichlet)
+ self%potential = potential / Volt_ref
+
+ CALL findNodes(self, physicalSurface)
+
+ END SELECT
+
+ END SUBROUTINE initDirichlet
+
+ ! Initialize Dirichlet boundary condition
+ SUBROUTINE initDirichletTime(self, physicalSurface, potential, temporalProfile)
+ USE moduleRefParam, ONLY: Volt_ref, ti_ref
+ IMPLICIT NONE
+
+ CLASS(boundaryEMGeneric), ALLOCATABLE, INTENT(out):: self
+ INTEGER, INTENT(in):: physicalSurface
+ REAL(8), INTENT(in):: potential
+ CHARACTER(:), ALLOCATABLE, INTENT(in):: temporalProfile
+
+ ! Allocate boundary edge
+ ALLOCATE(boundaryEMDirichletTime:: self)
+
+ SELECT TYPE(self)
+ TYPE IS(boundaryEMDirichletTime)
+ self%potential = potential / Volt_ref
+
+ CALL findNodes(self, physicalSurface)
+
+ CALL self%temporalProfile%init(temporalProfile)
+
+ CALL self%temporalProfile%convert(1.D0/ti_ref, 1.D0)
+
+ END SELECT
+
+ END SUBROUTINE initDirichletTime
+
+ !Apply Dirichlet boundary condition to the poisson equation
+ SUBROUTINE applyDirichlet(self, vectorF)
+ USE moduleMesh
+ IMPLICIT NONE
+
+ CLASS(boundaryEMDirichlet), INTENT(in):: self
+ REAL(8), INTENT(inout):: vectorF(:)
+ INTEGER:: n, ni
+
+ DO n = 1, self%nNodes
+ self%nodes(n)%obj%emData%phi = self%potential
+ vectorF(self%nodes(n)%obj%n) = self%nodes(n)%obj%emData%phi
+
+ END DO
+
+ 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(:)
+ REAL(8):: timeFactor
+ INTEGER:: n, ni
+
+ timeFactor = self%temporalProfile%get(DBLE(timeStep)*tauMin)
+
+ DO n = 1, self%nNodes
+ self%nodes(n)%obj%emData%phi = self%potential * timeFactor
+ vectorF(self%nodes(n)%obj%n) = self%nodes(n)%obj%emData%phi
+
+ END DO
+
+ END SUBROUTINE applyDirichletTime
!Assemble the source vector based on the charge density to solve Poisson's equation
SUBROUTINE assembleSourceVector(vectorF)
@@ -60,13 +207,14 @@ MODULE moduleEM
INTEGER, ALLOCATABLE:: nodes(:)
REAL(8), ALLOCATABLE:: rho(:)
INTEGER:: nNodes
- INTEGER:: e, i, ni
+ INTEGER:: e, i, ni, b
CLASS(meshNode), POINTER:: node
!$OMP SINGLE
vectorF = 0.D0
!$OMP END SINGLE
+ ! Calculate charge density in each node
!$OMP DO REDUCTION(+:vectorF)
DO e = 1, mesh%numCells
nNodes = mesh%cells(e)%obj%nNodes
@@ -98,18 +246,12 @@ MODULE moduleEM
!$OMP END DO
!Apply boundary conditions
- !$OMP DO
- DO i = 1, mesh%numNodes
- node => mesh%nodes(i)%obj
-
- SELECT CASE(node%emData%type)
- CASE ("dirichlet")
- vectorF(i) = node%emData%phi
-
- END SELECT
+ !$OMP SINGLE
+ DO b = 1, nBoundaryEM
+ CALL boundaryEM(b)%obj%apply(vectorF)
END DO
- !$OMP END DO
+ !$OMP END SINGLE
END SUBROUTINE assembleSourceVector
diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90
index 0b2837f..5928508 100644
--- a/src/modules/solver/moduleSolver.f90
+++ b/src/modules/solver/moduleSolver.f90
@@ -491,47 +491,46 @@ MODULE moduleSolver
END SUBROUTINE updateParticleCell
!Update the information about if a species needs to be moved this iteration
- SUBROUTINE updatePushSpecies(self, t)
+ SUBROUTINE updatePushSpecies(self)
USE moduleSpecies
+ USE moduleCaseparam, ONLY: timeStep
IMPLICIT NONE
CLASS(solverGeneric), INTENT(inout):: self
- INTEGER, INTENT(in):: t
INTEGER:: s
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 SUBROUTINE updatePushSpecies
!Output the different data and information
- SUBROUTINE doOutput(t)
+ SUBROUTINE doOutput()
USE moduleMesh
USE moduleOutput
USE moduleSpecies
USE moduleCompTime
USE moduleProbe
+ USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE
- INTEGER, INTENT(in):: t
-
- CALL outputProbes(t)
+ CALL outputProbes()
counterOutput = counterOutput + 1
IF (counterOutput >= triggerOutput .OR. &
- t == tFinal .OR. t == tInitial) THEN
+ timeStep == tFinal .OR. timeStep == tInitial) THEN
!Resets output counter
counterOutput=0
- CALL mesh%printOutput(t)
- IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl(t)
- CALL mesh%printEM(t)
- WRITE(*, "(5X,A21,I10,A1,I10)") "t/tFinal: ", t, "/", tFinal
+ CALL mesh%printOutput()
+ IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl()
+ CALL mesh%printEM()
+ WRITE(*, "(5X,A21,I10,A1,I10)") "t/tFinal: ", timeStep, "/", tFinal
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"
ELSE
@@ -549,34 +548,32 @@ MODULE moduleSolver
counterCPUTime = counterCPUTime + 1
IF (counterCPUTime >= triggerCPUTime .OR. &
- t == tFinal .OR. t == tInitial) THEN
+ timeStep == tFinal .OR. timeStep == tInitial) THEN
!Reset CPU Time counter
counterCPUTime = 0
- CALL printTime(t, t == 0)
+ CALL printTime(timeStep == 0)
END IF
!Output average values
- IF (useAverage .AND. t == tFinal) THEN
+ IF (useAverage .AND. timeStep == tFinal) THEN
CALL mesh%printAverage()
END IF
END SUBROUTINE doOutput
- SUBROUTINE doAverage(t)
+ SUBROUTINE doAverage()
USE moduleAverage
USE moduleMesh
IMPLICIT NONE
- INTEGER, INTENT(in):: t
INTEGER:: tAverage, n
-
IF (useAverage) THEN
- tAverage = t - tAverageStart
+ tAverage = timeStep - tAverageStart
IF (tAverage == 1) THEN
!First iteration in which average scheme is used