Boundaries done. Now to change the input and edges

This commit is contained in:
Jorge Gonzalez 2026-02-09 18:43:59 +01:00
commit 0ccf455986
5 changed files with 239 additions and 221 deletions

View file

@ -1,29 +1,30 @@
!moduleMeshBoundary: Boundary functions for the mesh edges
submodule(moduleMesh) boundary
CONTAINS
FUNCTION getBoundaryId(physicalSurface) RESULT(id)
IMPLICIT NONE
! FUNCTION getBoundaryId(physicalSurface) RESULT(id)
! IMPLICIT NONE
INTEGER:: physicalSurface
INTEGER:: id
INTEGER:: i
! INTEGER:: physicalSurface
! INTEGER:: id
! INTEGER:: i
id = 0
DO i = 1, nBoundary
IF (physicalSurface == boundaries(i)%physicalSurface) id = boundaries(i)%n
! id = 0
! DO i = 1, nBoundary
! IF (physicalSurface == boundariesParticle(i)%obj%physicalSurface) id = boundariesParticle(i)%obj%n
END DO
! END DO
END FUNCTION getBoundaryId
! END FUNCTION getBoundaryId
module subroutine initBoundary(self, config, object, i)
module subroutine initBoundary(self, config, object)
use json_module
use moduleRefParam, only: m_ref
use moduleConstParam, only: me
implicit none
class(boundaryGeneric), intent(out):: self
class(boundaryGeneric), allocatable, intent(out):: self
type(json_file), intent(inout):: config
character(:), allocatable, intent(in):: object
integer, intent(in):: i
character(:), allocatable:: bType
logical:: found
real(8):: Tw, cw !Wall temperature and specific heat
@ -35,86 +36,73 @@ submodule(moduleMesh) boundary
integer:: speciesID, electronSecondaryID
character(:), allocatable:: speciesName, crossSection, electronSecondary
self%n = i
CALL config%get(object // '.name', self%name, found)
CALL config%get(object // '.physicalSurface', self%physicalSurface, found)
CALL config%info(object // '.bTypes', found, n_children = nTypes)
IF (nTypes /= nSpecies) CALL criticalError('Not enough boundary types defined in ' // object, 'readBoundary')
ALLOCATE(self%bTypes(1:nSpecies))
DO s = 1, nSpecies
WRITE(sString,'(i2)') s
object = 'boundary(' // TRIM(iString) // ').bTypes(' // TRIM(sString) // ')'
CALL config%get(object // '.type', bType, found)
associate(bound => self%bTypes(s)%obj)
SELECT CASE(bType)
CASE('reflection')
ALLOCATE(boundaryReflection:: bound)
CALL config%get(object // '.name', self%name, found)
CALL config%get(object // '.type', bType, found)
SELECT CASE(bType)
CASE('reflection')
ALLOCATE(boundaryReflection:: self)
CASE('absorption')
ALLOCATE(boundaryAbsorption:: bound)
CASE('absorption')
ALLOCATE(boundaryAbsorption:: self)
CASE('transparent')
ALLOCATE(boundaryTransparent:: bound)
CASE('transparent')
ALLOCATE(boundaryTransparent:: self)
CASE('axis')
ALLOCATE(boundaryAxis:: bound)
CASE('axis')
ALLOCATE(boundaryAxis:: self)
CASE('wallTemperature')
CALL config%get(object // '.temperature', Tw, found)
IF (.NOT. found) CALL criticalError("temperature not found for wallTemperature boundary type", 'readBoundary')
CALL config%get(object // '.specificHeat', cw, found)
IF (.NOT. found) CALL criticalError("specificHeat not found for wallTemperature boundary type", 'readBoundary')
CASE('wallTemperature')
CALL config%get(object // '.temperature', Tw, found)
IF (.NOT. found) CALL criticalError("temperature not found for wallTemperature boundary type", 'readBoundary')
CALL config%get(object // '.specificHeat', cw, found)
IF (.NOT. found) CALL criticalError("specificHeat not found for wallTemperature boundary type", 'readBoundary')
CALL initWallTemperature(bound, Tw, cw)
CALL initWallTemperature(self, Tw, cw)
CASE('ionization')
!Neutral parameters
CALL config%get(object // '.neutral.ion', speciesName, found)
IF (.NOT. found) CALL criticalError("missing parameter 'ion' for neutrals in ionization", 'readBoundary')
speciesID = speciesName2Index(speciesName)
CALL config%get(object // '.neutral.mass', m0, found)
IF (.NOT. found) THEN
m0 = species(s)%obj%m*m_ref
END IF
CALL config%get(object // '.neutral.density', n0, found)
IF (.NOT. found) CALL criticalError("missing parameter 'density' for neutrals in ionization", 'readBoundary')
CALL config%get(object // '.neutral.velocity', v0, found)
IF (.NOT. found) CALL criticalError("missing parameter 'velocity' for neutrals in ionization", 'readBoundary')
CALL config%get(object // '.neutral.temperature', T0, found)
IF (.NOT. found) CALL criticalError("missing parameter 'temperature' for neutrals in ionization", 'readBoundary')
CASE('ionization')
!Neutral parameters
CALL config%get(object // '.neutral.ion', speciesName, found)
IF (.NOT. found) CALL criticalError("missing parameter 'ion' for neutrals in ionization", 'readBoundary')
speciesID = speciesName2Index(speciesName)
CALL config%get(object // '.neutral.mass', m0, found)
IF (.NOT. found) CALL criticalError("missing parameter 'mass' for neutrals in ionization", 'readBoundary')
CALL config%get(object // '.neutral.density', n0, found)
IF (.NOT. found) CALL criticalError("missing parameter 'density' for neutrals in ionization", 'readBoundary')
CALL config%get(object // '.neutral.velocity', v0, found)
IF (.NOT. found) CALL criticalError("missing parameter 'velocity' for neutrals in ionization", 'readBoundary')
CALL config%get(object // '.neutral.temperature', T0, found)
IF (.NOT. found) CALL criticalError("missing parameter 'temperature' for neutrals in ionization", 'readBoundary')
CALL config%get(object // '.effectiveTime', effTime, found)
IF (.NOT. found) CALL criticalError("missing parameter 'effectiveTime' for ionization", 'readBoundary')
CALL config%get(object // '.effectiveTime', effTime, found)
IF (.NOT. found) CALL criticalError("missing parameter 'effectiveTime' for ionization", 'readBoundary')
CALL config%get(object // '.energyThreshold', eThreshold, found)
IF (.NOT. found) CALL criticalError("missing parameter 'eThreshold' in ionization", 'readBoundary')
CALL config%get(object // '.energyThreshold', eThreshold, found)
IF (.NOT. found) CALL criticalError("missing parameter 'eThreshold' in ionization", 'readBoundary')
CALL config%get(object // '.crossSection', crossSection, found)
IF (.NOT. found) CALL criticalError("missing parameter 'crossSection' for neutrals in ionization", 'readBoundary')
CALL config%get(object // '.crossSection', crossSection, found)
IF (.NOT. found) CALL criticalError("missing parameter 'crossSection' for neutrals in ionization", 'readBoundary')
CALL config%get(object // '.electronSecondary', electronSecondary, found)
electronSecondaryID = speciesName2Index(electronSecondary)
IF (found) THEN
CALL initIonization(bound, species(s)%obj%m, m0, n0, v0, T0, &
speciesID, effTime, crossSection, eThreshold,electronSecondaryID)
CALL config%get(object // '.electronSecondary', electronSecondary, found)
electronSecondaryID = speciesName2Index(electronSecondary)
IF (found) THEN
CALL initIonization(self, me/m_ref, m0, n0, v0, T0, &
speciesID, effTime, crossSection, eThreshold,electronSecondaryID)
ELSE
CALL initIonization(bound, species(s)%obj%m, m0, n0, v0, T0, &
speciesID, effTime, crossSection, eThreshold)
ELSE
CALL initIonization(self, me/m_ref, m0, n0, v0, T0, &
speciesID, effTime, crossSection, eThreshold)
END IF
END IF
case('quasiNeutrality')
call initQuasiNeutrality(bound)
case('quasiNeutrality')
call initQuasiNeutrality(self)
CASE DEFAULT
CALL criticalError('Boundary type ' // bType // ' undefined', 'readBoundary')
CASE DEFAULT
CALL criticalError('Boundary type ' // bType // ' undefined', 'readBoundary')
END SELECT
end associate
END SELECT
END DO
end subroutine initBoundary
SUBROUTINE initWallTemperature(boundary, T, c)
@ -130,7 +118,7 @@ submodule(moduleMesh) boundary
END SUBROUTINE initWallTemperature
SUBROUTINE initIonization(boundary, me, m0, n0, v0, T0, ion, effTime, crossSection, eThreshold, electronSecondary)
SUBROUTINE initIonization(boundary, mImpact, m0, n0, v0, T0, ion, effTime, crossSection, eThreshold, electronSecondary)
USE moduleRefParam
USE moduleSpecies
USE moduleCaseParam
@ -139,7 +127,7 @@ submodule(moduleMesh) boundary
IMPLICIT NONE
CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary
REAL(8), INTENT(in):: me !Electron mass
real(8), intent(in):: mImpact
REAL(8), INTENT(in):: m0, n0, v0(1:3), T0 !Neutral properties
INTEGER, INTENT(in):: ion
INTEGER, OPTIONAL, INTENT(in):: electronSecondary
@ -175,7 +163,7 @@ submodule(moduleMesh) boundary
CALL boundary%crossSection%init(crossSection)
CALL boundary%crossSection%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2)
boundary%eThreshold = eThreshold*eV2J/(m_ref*v_ref**2)
boundary%deltaV = DSQRT(boundary%eThreshold/me)
boundary%deltaV = DSQRT(boundary%eThreshold/mImpact)
END SELECT
@ -197,13 +185,14 @@ submodule(moduleMesh) boundary
end subroutine initQuasiNeutrality
module SUBROUTINE reflection(edge, part)
module SUBROUTINE reflection(edge, part, self)
USE moduleCaseParam
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
class(boundaryGeneric), optional, intent(in):: self
!rp = intersection between particle and edge
!rpp = final position of particle
!vpp = final velocity of particle
@ -223,13 +212,14 @@ submodule(moduleMesh) boundary
END SUBROUTINE reflection
!Absoption in a surface
SUBROUTINE absorption(edge, part)
SUBROUTINE absorption(edge, part, self)
USE moduleCaseParam
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
class(boundaryGeneric), optional, intent(in):: self
REAL(8):: rpp(1:3) !Position of particle projected to the edge
REAL(8):: d !Distance from particle to edge
@ -259,12 +249,13 @@ submodule(moduleMesh) boundary
END SUBROUTINE absorption
!Transparent boundary condition
SUBROUTINE transparent(edge, part)
SUBROUTINE transparent(edge, part, self)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
class(boundaryGeneric), optional, intent(in):: self
!Removes particle from domain
part%n_in = .FALSE.
@ -274,44 +265,46 @@ submodule(moduleMesh) boundary
!Symmetry axis. Reflects particles.
!Although this function should never be called, it is set as a reflective boundary
!to properly deal with possible particles reaching a corner and selecting this boundary.
SUBROUTINE symmetryAxis(edge, part)
SUBROUTINE symmetryAxis(edge, part, self)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
class(boundaryGeneric), optional, intent(in):: self
CALL reflection(edge, part)
END SUBROUTINE symmetryAxis
!Wall with temperature
SUBROUTINE wallTemperature(edge, part)
SUBROUTINE wallTemperature(edge, part, self)
USE moduleSpecies
USE moduleRandom
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
class(boundaryGeneric), optional, intent(in):: self
INTEGER:: i
!Modifies particle velocity according to wall temperature
SELECT TYPE(bound => edge%boundary%bTypes(part%species%n)%obj)
TYPE IS(boundaryWallTemperature)
select type(self)
type is(boundaryWallTemperature)
!Modifies particle velocity according to wall temperature
DO i = 1, 3
part%v(i) = part%v(i) + bound%vTh*randomMaxwellian()
part%v(i) = part%v(i) + self%vTh*randomMaxwellian()
END DO
END SELECT
CALL reflection(edge, part)
CALL reflection(edge, part)
end select
END SUBROUTINE wallTemperature
!Ionization surface: an electron will pass through the surface
! and create an ion-electron pair based on a neutral background
SUBROUTINE ionization(edge, part)
SUBROUTINE ionization(edge, part, self)
USE moduleList
USE moduleSpecies
USE moduleMesh
@ -320,8 +313,9 @@ submodule(moduleMesh) boundary
USE moduleMath
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
class(boundaryGeneric), optional, intent(in):: self
REAL(8):: vRel, eRel, mRel !relative velocity, energy and mass
INTEGER:: nIonizations !Number of ionizations based on eRel
REAL(8):: pIonization !Probability of ionization of each event
@ -330,40 +324,40 @@ submodule(moduleMesh) boundary
TYPE(particle), POINTER:: newElectron
TYPE(particle), POINTER:: newIon
SELECT TYPE(bound => edge%boundary%bTypes(part%species%n)%obj)
TYPE IS(boundaryIonization)
mRel = reducedMass(bound%m0, part%species%m)
vRel = SUM(DABS(part%v-bound%v0))
select type(self)
type is(boundaryIonization)
mRel = reducedMass(self%m0, part%species%m)
vRel = SUM(DABS(part%v-self%v0))
eRel = mRel*vRel**2*5.D-1
!Maximum number of possible ionizations based on relative energy
nIonizations = FLOOR(eRel/bound%eThreshold)
nIonizations = FLOOR(eRel/self%eThreshold)
DO p = 1, nIonizations
!Get probability of ionization
pIonization = 1.D0 - DEXP(-bound%n0*bound%crossSection%get(eRel)*vRel*bound%effectiveTime/REAL(nIonizations))
pIonization = 1.D0 - DEXP(-self%n0*self%crossSection%get(eRel)*vRel*self%effectiveTime/REAL(nIonizations))
!If a random number is below the probability of ionization, create new pair of ion-electron
IF (random() < pIonization) THEN
!Assign random velocity to the neutral
v0(1) = bound%v0(1) + bound%vTh*randomMaxwellian()
v0(2) = bound%v0(2) + bound%vTh*randomMaxwellian()
v0(3) = bound%v0(3) + bound%vTh*randomMaxwellian()
v0(1) = self%v0(1) + self%vTh*randomMaxwellian()
v0(2) = self%v0(2) + self%vTh*randomMaxwellian()
v0(3) = self%v0(3) + self%vTh*randomMaxwellian()
!Allocates the new particles
ALLOCATE(newElectron)
ALLOCATE(newIon)
IF (ASSOCIATED(bound%electronSecondary)) THEN
newElectron%species => bound%electronSecondary
IF (ASSOCIATED(self%electronSecondary)) THEN
newElectron%species => self%electronSecondary
ELSE
newElectron%species => part%species
END IF
newIon%species => bound%species
newIon%species => self%species
newElectron%v = v0 + (1.D0 + bound%deltaV*v0/NORM2(v0))
newElectron%v = v0 + (1.D0 + self%deltaV*v0/NORM2(v0))
newIon%v = v0
newElectron%r = edge%randPos()
@ -394,7 +388,7 @@ submodule(moduleMesh) boundary
CALL partSurfaces%unsetLock()
!Electron loses energy due to ionization
eRel = eRel - bound%eThreshold
eRel = eRel - self%eThreshold
vRel = 2.D0*DSQRT(eRel)/mRel
!Reduce number of possible ionizations
@ -404,83 +398,79 @@ submodule(moduleMesh) boundary
END DO
END SELECT
!Removes ionizing electron regardless the number of pair created
part%n_in = .FALSE.
!Removes ionizing electron regardless the number of pair created
part%n_in = .FALSE.
end select
END SUBROUTINE ionization
subroutine quasiNeutrality(edge, part)
subroutine quasiNeutrality(edge, part, self)
use moduleRandom
implicit none
class(meshEdge), intent(inout):: edge
class(particle), intent(inout):: part
class(meshEdge), intent(inout):: edge
class(particle), intent(inout):: part
class(boundaryGeneric), optional, intent(in):: self
real(8), allocatable:: density(:)
class(meshCell), pointer:: cell
real(8):: EF_dir
real(8):: alpha
select type(bound => edge%boundary%bTypes(part%species%n)%obj)
type is(boundaryQuasiNeutrality)
if (associated(edge%e1)) then
cell => edge%e1
if (associated(edge%e1)) then
cell => edge%e1
else
cell => edge%e2
else
cell => edge%e2
end if
if (random() <= alpha) then
call reflection(edge, part)
end if
if (random() <= alpha) then
call reflection(edge, part)
else
call transparent(edge, part)
else
call transparent(edge, part)
end if
end select
end if
end subroutine quasiNeutrality
!Points the boundary function to specific type
module SUBROUTINE pointBoundaryFunction(edge, s)
USE moduleErrors
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
INTEGER, INTENT(in):: s !Species index
SELECT TYPE(obj => edge%boundary%bTypes(s)%obj)
TYPE IS(boundaryAbsorption)
edge%fBoundary(s)%apply => absorption
TYPE IS(boundaryReflection)
edge%fBoundary(s)%apply => reflection
TYPE IS(boundaryTransparent)
edge%fBoundary(s)%apply => transparent
TYPE IS(boundaryAxis)
edge%fBoundary(s)%apply => symmetryAxis
TYPE IS(boundaryWallTemperature)
edge%fBoundary(s)%apply => wallTemperature
TYPE IS(boundaryIonization)
edge%fBoundary(s)%apply => ionization
type is(boundaryQuasiNeutrality)
edge%fBoundary(s)%apply => quasiNeutrality
CLASS DEFAULT
CALL criticalError("Boundary type not defined", 'pointBoundaryFunction')
END SELECT
END SUBROUTINE pointBoundaryFunction
! !Points the boundary function to specific type
! module SUBROUTINE pointBoundaryFunction(edge, s)
! USE moduleErrors
! IMPLICIT NONE
!
! CLASS(meshEdge), INTENT(inout):: edge
! INTEGER, INTENT(in):: s !Species index
!
! SELECT TYPE(obj => edge%boundary%bTypes(s)%obj)
! TYPE IS(boundaryAbsorption)
! edge%fBoundary(s)%apply => absorption
!
! TYPE IS(boundaryReflection)
! edge%fBoundary(s)%apply => reflection
!
! TYPE IS(boundaryTransparent)
! edge%fBoundary(s)%apply => transparent
!
! TYPE IS(boundaryAxis)
! edge%fBoundary(s)%apply => symmetryAxis
!
! TYPE IS(boundaryWallTemperature)
! edge%fBoundary(s)%apply => wallTemperature
!
! TYPE IS(boundaryIonization)
! edge%fBoundary(s)%apply => ionization
!
! type is(boundaryQuasiNeutrality)
! edge%fBoundary(s)%apply => quasiNeutrality
!
! CLASS DEFAULT
! CALL criticalError("Boundary type not defined", 'pointBoundaryFunction')
!
! END SELECT
!
! END SUBROUTINE pointBoundaryFunction
end submodule boundary