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

Implementing output for injects. Cleaning code

Browse source
This commit is contained in:
Jorge Gonzalez 2026-04-01 16:19:16 +02:00
commit 04cd9dffc6
10 changed files with 283 additions and 169 deletions
Download patch file Download diff file Expand all files Collapse all files

2
src/fpakc.f90
View file

@ -72,7 +72,7 @@ PROGRAM fpakc
call boundariesEM_update()
! Update injects
call updateInjects
call updateInjects()
!Checks if a species needs to me moved in this iteration
CALL solver%updatePushSpecies()

2
src/modules/common/moduleRandom.f90
View file

@ -53,7 +53,7 @@ MODULE moduleRandom
implicit none
real(8):: rnd
real(8):: v1, v2, Rsquare
real(8):: v1, v2
v1 = 0.d0
do while (v1 <= 0.d0)

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

@ -628,6 +628,7 @@ MODULE moduleInput
USE moduleMesh
USE moduleCaseParam
USE moduleRefParam
use moduleOutput, only: collOutput
USE OMP_LIB
USE json_module
IMPLICIT NONE
@ -1426,50 +1427,17 @@ MODULE moduleInput
CHARACTER(2):: iString
CHARACTER(:), ALLOCATABLE:: object
LOGICAL:: found
CHARACTER(:), ALLOCATABLE:: speciesName
CHARACTER(:), ALLOCATABLE:: name
REAL(8):: v
REAL(8), ALLOCATABLE:: temperature(:), normal(:)
REAL(8):: flow
CHARACTER(:), ALLOCATABLE:: units
character(:), allocatable:: type
INTEGER:: physicalSurface
INTEGER:: particlesPerEdge
INTEGER:: sp
CALL config%info('inject', found, n_children = nInject)
ALLOCATE(inject(1:nInject))
ALLOCATE(injects(1:nInject))
nPartInj = 0
DO i = 1, nInject
WRITE(iString, '(i2)') i
object = 'inject(' // trim(iString) // ')'
!Find species
CALL config%get(object // '.species', speciesName, found)
sp = speciesName2Index(speciesName)
CALL initInject(injects(i)%obj, i, object, config)
CALL config%get(object // '.name', name, found)
CALL config%get(object // '.v', v, found)
CALL config%get(object // '.T', temperature, found)
CALL config%get(object // '.n', normal, found)
IF (.NOT. found) THEN
ALLOCATE(normal(1:3))
normal = 0.D0
END IF
CALL config%get(object // '.flow', flow, found)
CALL config%get(object // '.units', units, found)
CALL config%get(object // '.type', type, found)
if (.not. found) then
! If no type is found, assume constant injection of particles
type = 'constant'
end if
CALL config%get(object // '.physicalSurface', physicalSurface, found)
particlesPerEdge = 0
CALL config%get(object // '.particlesPerEdge', particlesPerEdge, found)
CALL inject(i)%init(i, v, normal, temperature, flow, units, type, sp, physicalSurface, particlesPerEdge)
CALL readVelDistr(config, inject(i), object)
CALL readVelDistr(config, injects(i)%obj, object)
END DO
@ -1528,7 +1496,7 @@ MODULE moduleInput
IMPLICIT NONE
TYPE(json_file), INTENT(inout):: config
TYPE(injectGeneric), INTENT(inout):: inj
class(injectGeneric), INTENT(inout):: inj
CHARACTER(:), ALLOCATABLE, INTENT(in):: object
INTEGER:: i
CHARACTER(2):: iString

11
src/modules/mesh/inout/text/moduleMeshOutputText.f90
View file

@ -1,9 +1,14 @@
module moduleMeshOutputText
use moduleOutput, only: prefix, &
informFileCreation, &
formatFileName, &
generateFilePath
contains
subroutine writeSpeciesOutput(self, fileID, speciesIndex)
use moduleMesh
use moduleOutput
use moduleOutput, only: calculateOutput, outputFormat
use moduleRefParam, only: L_ref
implicit none
@ -88,7 +93,6 @@ module moduleMeshOutputText
speciesIndex)
use moduleMesh
use moduleOutput
use moduleAverage
use moduleRefParam, only: L_ref
implicit none
@ -151,8 +155,8 @@ module moduleMeshOutputText
subroutine printCollText(self)
use moduleMesh
use moduleOutput
use moduleCaseParam, only: timeStep
use moduleOutput, only: collOutput
implicit none
class(meshGeneric), intent(in):: self
@ -190,6 +194,7 @@ module moduleMeshOutputText
subroutine printEMText(self)
use moduleMesh
use moduleCaseParam, only: timeStep
use moduleOutput, only: emOutput
implicit none
class(meshParticles), intent(in):: self

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

@ -1,9 +1,11 @@
MODULE moduleMeshOutputVTU
use moduleOutput, only: prefix, informFileCreation, formatFileName, generateFilePath
CONTAINS
SUBROUTINE writeHeader(nNodes, nCells, fileID)
USE moduleMesh
use moduleOutput, only: fmtInt
IMPLICIT NONE
INTEGER, INTENT(in):: nNodes, nCells
@ -77,6 +79,7 @@ MODULE moduleMeshOutputVTU
SUBROUTINE writeMesh(self, fileID)
USE moduleMesh
USE moduleRefParam
use moduleOutput, only: fmtReal, fmtInt
IMPLICIT NONE
CLASS(meshGeneric), INTENT(in):: self
@ -159,6 +162,7 @@ MODULE moduleMeshOutputVTU
SUBROUTINE writeCollOutput(self,fileID)
USE moduleMesh
USE moduleCollisions
use moduleOutput, only: fmtInt
IMPLICIT NONE
CLASS(meshGeneric), INTENT(in):: self
@ -185,6 +189,7 @@ MODULE moduleMeshOutputVTU
SUBROUTINE writeEM(self, fileID)
USE moduleMesh
USE moduleRefParam
use moduleOutput, only: fmtReal
IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self
@ -381,6 +386,7 @@ MODULE moduleMeshOutputVTU
SUBROUTINE printEMVTU(self)
USE moduleMesh
USE moduleCaseParam, ONLY: timeStep
use moduleOutput, only: emOutput
IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self
@ -414,7 +420,6 @@ MODULE moduleMeshOutputVTU
SUBROUTINE printAverageVTU(self)
USE moduleMesh
use moduleOutput
USE moduleSpecies
IMPLICIT NONE

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

@ -1,7 +1,7 @@
!moduleMesh: General module for Finite Element mesh
MODULE moduleMesh
USE moduleList
USE moduleOutput
USE moduleOutput, only: outputNode, emNode
USE moduleCollisions
use moduleSpecies, only: nSpecies

1
src/modules/mesh/moduleMesh@boundaryEM.f90
View file

@ -56,6 +56,7 @@ submodule(moduleMesh) boundaryEM
use json_module
use moduleRefParam, ONLY: Volt_ref, ti_ref
use moduleErrors
use moduleOutput, only: path
implicit none
class(boundaryEMGeneric), allocatable, intent(inout):: self

12
src/modules/mesh/moduleMesh@boundaryParticle.f90
View file

@ -342,6 +342,7 @@ submodule(moduleMesh) boundaryParticle
! Print
subroutine ionization_print(self, fileID)
use moduleOutput, only: fmtInt
implicit none
class(boundaryParticleGeneric), intent(inout):: self
@ -473,6 +474,7 @@ submodule(moduleMesh) boundaryParticle
! Print output
subroutine quasiNeutrality_print(self, fileID)
use moduleOutput, only: fmtColInt, fmtColReal
implicit none
class(boundaryParticleGeneric), intent(inout):: self
@ -613,6 +615,7 @@ submodule(moduleMesh) boundaryParticle
! Print
subroutine outflowAdaptive_print(self, fileID)
use moduleOutput, only: fmtColInt, fmtColReal
implicit none
class(boundaryParticleGeneric), intent(inout):: self
@ -691,7 +694,12 @@ submodule(moduleMesh) boundaryParticle
! Writes the output into the Step_XXXXX_boundaryParticles file when the %print procedure is associated
module subroutine boundariesParticle_write()
use moduleCaseparam, only: timeStep
use moduleOutput, only:fileID_boundaryParticle, formatFileName, informFileCreation
use moduleOutput, only: fileID_boundaryParticle, &
formatFileName, &
informFileCreation,&
boundaryParticleOutput,&
generateFilePath, &
prefix
implicit none
integer:: b
@ -700,7 +708,7 @@ submodule(moduleMesh) boundaryParticle
if (boundaryParticleOutput) then
fileName = formatFileName(prefix, 'boundariesParticle', 'csv', timeStep)
call informFileCreation(fileName)
open(fileID_boundaryParticle, file = path // folder // '/' // fileName)
open(fileID_boundaryParticle, file = generateFilePath(fileName))
do b = 1, nBoundariesParticle
if (associated(boundariesParticle(b)%obj%print)) then

371
src/modules/moduleInject.f90
View file

@ -5,7 +5,7 @@ MODULE moduleInject
use moduleMesh, only: meshEdgePointer
!Generic injection of particles
TYPE:: injectGeneric
TYPE, abstract:: injectGeneric
INTEGER:: id
CHARACTER(:), ALLOCATABLE:: name
REAL(8):: vMod !Velocity (module)
@ -14,46 +14,130 @@ MODULE moduleInject
LOGICAL:: fixDirection !The injection of particles has a fix direction defined by n
INTEGER:: nParticles !Number of particles to introduce each time step
CLASS(speciesGeneric), POINTER:: species !Species of injection
character(:), allocatable:: type
INTEGER:: nEdges
type(meshEdgePointer), allocatable:: edges(:)
INTEGER, ALLOCATABLE:: particlesPerEdge(:) ! Particles per edge
REAL(8), ALLOCATABLE:: weightPerEdge(:) ! Weight per edge
REAL(8):: surface ! Total surface of injection
TYPE(velDistCont):: v(1:3) !Velocity distribution function in each direction
procedure(updateInject_interface), pointer, pass:: update => null()
procedure(printInject_interface), pointer, pass:: print => null()
CONTAINS
PROCEDURE, PASS:: init => initInject
PROCEDURE, PASS:: addParticles
END TYPE injectGeneric
abstract interface
! Update the values of the particle boundary model
subroutine updateInject_interface(self)
import injectGeneric
class(injectGeneric), intent(inout):: self
end subroutine updateInject_interface
! Write the values of the particle boundary model
subroutine printInject_interface(self, fileID)
import injectGeneric
class(injectGeneric), intent(inout):: self
integer, intent(in):: fileID
end subroutine printInject_interface
end interface
! Default type, constant inject of particles
type, extends(injectGeneric):: injectConstant
end type injectconstant
type, extends(injectGeneric):: injectQuasiNeutral
end type injectQuasiNeutral
! Wrapper for injects
type:: injectWrapper
class(injectGeneric), allocatable:: obj
end type injectWrapper
INTEGER:: nInject
TYPE(injectGeneric), ALLOCATABLE, target:: inject(:)
TYPE(injectWrapper), ALLOCATABLE, target:: injects(:)
CONTAINS
!Initialize an injection of particles
SUBROUTINE initInject(self, i, v, n, temperature, flow, units, type, sp, ps, particlesPerEdge)
SUBROUTINE initInject(self, i, object, config)
USE moduleMesh
USE moduleRefParam
USE moduleConstParam
USE moduleSpecies
USE moduleSolver
USE moduleErrors
use json_module
IMPLICIT NONE
CLASS(injectGeneric), INTENT(inout):: self
CLASS(injectGeneric), allocatable, INTENT(inout):: self
INTEGER, INTENT(in):: i
REAL(8), INTENT(in):: v, n(1:3), temperature(1:3)
INTEGER, INTENT(in):: sp, ps, particlesPerEdge
character(:), allocatable, intent(in):: object
type(json_file), intent(inout):: config
character(:), allocatable:: type
logical:: found
REAL(8):: v
real(8), allocatable:: n(:), temperature(:)
INTEGER:: sp, ps, particlesPerEdge
character(:), allocatable:: speciesName
integer:: ps_index
REAL(8):: tauInject
REAL(8), INTENT(in):: flow
CHARACTER(:), ALLOCATABLE, INTENT(in):: units
character(:), allocatable, intent(in):: type
REAL(8):: flow
CHARACTER(:), ALLOCATABLE:: units
INTEGER:: e
REAL(8):: fluxPerStep = 0.D0
self%id = i
!Type of injection
call config%get(object // '.type', type, found)
if (.not. found) then
! If no type is found, assume constant injection of particles
type = 'constant'
end if
! Assign type to inject and allocates it
select case(type)
case ('constant')
allocate(injectConstant:: self)
case ('quasiNeutral')
allocate(injectQuasiNeutral:: self)
self%update => updateQuasiNeutral
case default
call criticalError('No injection type ' // type // ' defined', 'initInject')
end select
self%id = i
call config%get(object // '.name', self%name, found)
!Find species
CALL config%get(object // '.species', speciesName, found)
sp = speciesName2Index(speciesName)
! Velocity and direction of injection
call config%get(object // '.v', v, found)
call config%get(object // '.n', n, found)
if (.not. found) then
allocate(n(1:3))
n = 0.D0
end if
! Temperature of injection
call config%get(object // '.T', temperature, found)
! Flow
call config%get(object // '.flow', flow, found)
call config%get(object // '.units', units, found)
CALL config%get(object // '.physicalSurface', ps, found)
particlesPerEdge = 0
CALL config%get(object // '.particlesPerEdge', particlesPerEdge, found)
self%vMod = v / v_ref
self%n = n / NORM2(n)
self%temperature = temperature / T_ref
@ -145,15 +229,91 @@ MODULE moduleInject
!Scale particles for different species steps
IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject')
! Assign type to inject
select case(type)
case ('constant', 'quasiNeutral')
self%type = type
case default
call criticalError('No injection type ' // type // ' defined', 'initInject')
END SUBROUTINE initInject
! Update the value of injects
subroutine updateInjects()
implicit none
integer:: i
do i = 1, nInject
if (associated(injects(i)%obj%update)) then
call injects(i)%obj%update()
end if
end do
end subroutine updateInjects
! Updates the flow in the injection to maintain quasineutrality
subroutine updateQuasiNeutral(self)
use moduleMesh, only: meshEdge, meshNode, mesh, qSpecies
implicit none
class(injectGeneric), intent(inout):: self
integer:: e, s, n
integer, allocatable:: nodes(:)
class(meshEdge), pointer:: edge
real(8), allocatable:: density_nodes(:)
class(meshNode), pointer:: node
real(8):: den_center
real(8):: density_incident, density_rest
real(8):: alpha
select type(self)
type is (injectQuasiNeutral)
do e = 1, self%nEdges
edge => self%edges(e)%obj
density_incident = 0.d0
density_rest = 0.d0
nodes = edge%getNodes(edge%nNodes)
allocate(density_nodes(1:edge%nNodes))
do s = 1, nSpecies
do n = 1, edge%nNodes
node => mesh%nodes(nodes(n))%obj
density_nodes(n) = node%output(s)%den
end do
! Gather the density at the edge center and multiply by the species charge
den_center = qSpecies(s)*edge%gatherF(edge%centerXi(), edge%nNodes, density_nodes)
if (s == self%species%n) then
density_incident = den_center
else
density_rest = density_rest + den_center
end if
end do
! Correction for this time step
if (density_rest > 1.0d-10) then
! If there is a rest population, correct
alpha = 1.d0 + density_incident/density_rest
else
! If not, alpha is assumed unchaged
alpha = 0.d0
end if
print *, e, alpha
! Adjust the weight of particles to match the new current
self%weightPerEdge(e) = self%weightPerEdge(e) + 1.0d-1 * alpha
end do
end select
END SUBROUTINE initInject
end subroutine updateQuasiNeutral
!Injection of particles
SUBROUTINE doInjects()
@ -166,8 +326,8 @@ MODULE moduleInject
!$OMP SINGLE
nPartInj = 0
DO i = 1, nInject
IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN
nPartInj = nPartInj + inject(i)%nParticles
IF (solver%pusher(injects(i)%obj%species%n)%pushSpecies) THEN
nPartInj = nPartInj + injects(i)%obj%nParticles
END IF
@ -178,112 +338,14 @@ MODULE moduleInject
!$OMP END SINGLE
DO i=1, nInject
IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN
CALL inject(i)%addParticles()
IF (solver%pusher(injects(i)%obj%species%n)%pushSpecies) THEN
CALL injects(i)%obj%addParticles()
END IF
END DO
END SUBROUTINE doInjects
! Update the value of injects
subroutine updateInjects()
use moduleMesh, only: meshEdge, meshNode, mesh, qSpecies
implicit none
integer:: i, e, s, n
class(injectGeneric), pointer:: self
integer, allocatable:: nodes(:)
class(meshEdge), pointer:: edge
real(8), allocatable:: density_nodes(:)
class(meshNode), pointer:: node
real(8):: den_center
real(8):: density_incident, density_rest
real(8):: alpha
do i = 1, nInject
self => inject(i)
select case(inject(i)%type)
case ('quasiNeutral')
do e = 1, self%nEdges
edge => self%edges(e)%obj
density_incident = 0.d0
density_rest = 0.d0
nodes = edge%getNodes(edge%nNodes)
allocate(density_nodes(1:edge%nNodes))
do s = 1, nSpecies
do n = 1, edge%nNodes
node => mesh%nodes(nodes(n))%obj
density_nodes(n) = node%output(s)%den
end do
! Gather the density at the edge center and multiply by the species charge
den_center = qSpecies(s)*edge%gatherF(edge%centerXi(), edge%nNodes, density_nodes)
if (s == self%species%n) then
density_incident = den_center
else
density_rest = density_rest + den_center
end if
end do
! Correction for this time step
if (density_rest > 1.0d-10) then
! If there is a rest population, correct
alpha = 1.d0 + density_incident/density_rest
else
! If not, alpha is assumed unchaged
alpha = 0.d0
end if
! Adjust the weight of particles to match the new current
self%weightPerEdge(e) = self%weightPerEdge(e) + 1.0d-2 * alpha
end do
end select
end do
end subroutine updateInjects
SUBROUTINE initVelDistMaxwellian(velDist, temperature, m)
IMPLICIT NONE
CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
REAL(8), INTENT(in):: temperature, m
velDist = velDistMaxwellian(vTh = DSQRT(2.d0*temperature/m))
END SUBROUTINE initVelDistMaxwellian
SUBROUTINE initVelDistHalfMaxwellian(velDist, temperature, m)
IMPLICIT NONE
CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
REAL(8), INTENT(in):: temperature, m
velDist = velDistHalfMaxwellian(vTh = DSQRT(2.d0*temperature/m))
END SUBROUTINE initVelDistHalfMaxwellian
SUBROUTINE initVelDistDelta(velDist)
IMPLICIT NONE
CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
velDist = velDistDelta()
END SUBROUTINE initVelDistDelta
!Add particles for the injection
SUBROUTINE addParticles(self)
USE moduleSpecies
@ -304,8 +366,8 @@ MODULE moduleInject
!$OMP SINGLE
nMin = 0
DO i = 1, self%id -1
IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN
nMin = nMin + inject(i)%nParticles
IF (solver%pusher(injects(i)%obj%species%n)%pushSpecies) THEN
nMin = nMin + injects(i)%obj%nParticles
END IF
@ -381,4 +443,67 @@ MODULE moduleInject
END SUBROUTINE addParticles
! Writes the value of injects
subroutine writeInjects()
use moduleCaseparam, only: timeStep
use moduleOutput, only: fileID_inject, &
formatFileName, &
informFileCreation,&
injectOutput,&
generateFilePath, &
prefix
implicit none
integer:: i
character(:), allocatable:: fileName
if (injectOutput) then
fileName = formatFileName(prefix, 'injects', 'csv', timeStep)
call informFileCreation(fileName)
open(fileID_inject, file = generateFilePath(fileName))
do i = 1, nInject
if (associated(injects(i)%obj%print)) then
call injects(i)%obj%print(fileID_inject)
end if
end do
close(fileID_inject)
end if
end subroutine writeInjects
! Inits for different velocity distribution functions
SUBROUTINE initVelDistMaxwellian(velDist, temperature, m)
IMPLICIT NONE
CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
REAL(8), INTENT(in):: temperature, m
velDist = velDistMaxwellian(vTh = DSQRT(2.d0*temperature/m))
END SUBROUTINE initVelDistMaxwellian
SUBROUTINE initVelDistHalfMaxwellian(velDist, temperature, m)
IMPLICIT NONE
CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
REAL(8), INTENT(in):: temperature, m
velDist = velDistHalfMaxwellian(vTh = DSQRT(2.d0*temperature/m))
END SUBROUTINE initVelDistHalfMaxwellian
SUBROUTINE initVelDistDelta(velDist)
IMPLICIT NONE
CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
velDist = velDistDelta()
END SUBROUTINE initVelDistDelta
END MODULE moduleInject

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

@ -14,6 +14,7 @@ MODULE moduleOutput
LOGICAL:: collOutput = .FALSE.
LOGICAL:: emOutput = .FALSE.
logical:: boundaryParticleOutput = .false.
logical:: injectOutput = .false.
! Prefix for iteration files
character(len=*), parameter:: prefix = 'Step'
@ -33,6 +34,7 @@ MODULE moduleOutput
integer, parameter:: fileID_output = 20 ! Base id for species/collisions/EM output
integer, parameter:: fileID_boundaryParticle = 30 ! Particle boundaries
integer, parameter:: fileID_boundaryEM = 31 ! EM boundaries
integer, parameter:: fileID_inject = 32 ! Injects
integer, parameter:: fileID_reference = 40 ! Reference values
integer, parameter:: fileID_time =50 ! Computation time