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

Merge pull request 'issue/edgesVolume' (#59) from issue/edgesVolume into development

Browse source 
Reviewed-on: https://plasmalab.aero.upm.es/git/git/JGonzalez/fpakc/pulls/59
This commit is contained in:
Jorge Gonzalez 2026-02-26 19:21:52 +01:00
commit 796176b4b4
5 changed files with 126 additions and 92 deletions
Download patch file Download diff file Expand all files Collapse all files

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

@ -48,24 +48,23 @@ MODULE moduleRandom
END FUNCTION randomIntAB END FUNCTION randomIntAB
!Returns a random number in a Maxwellian distribution of mean 0 and width 1 with the Box-Muller Method !Returns a random number in a Maxwellian distribution of mean 0 and width 1 with the Box-Muller Method
FUNCTION randomMaxwellian() RESULT(rnd) function randomMaxwellian() result(rnd)
USE moduleConstParam, ONLY: PI USE moduleConstParam, only: pi
IMPLICIT NONE implicit none
REAL(8):: rnd real(8):: rnd
REAL(8):: v1, v2, Rsquare real(8):: v1, v2, Rsquare
Rsquare = 1.D0 v1 = 0.d0
do while (Rsquare >= 1.D0 .and. Rsquare > 0.D0) do while (v1 <= 0.d0)
v1 = 2.D0 * random() - 1.D0 v1 = random()
v2 = 2.D0 * random() - 1.D0
Rsquare = v1**2 + v2**2
end do end do
v2 = random()
rnd = v2 * sqrt(-2.D0 * log(Rsquare) / Rsquare) rnd = sqrt(-2.d0*log(v1))*cos(2*pi*v2)
END FUNCTION randomMaxwellian end function randomMaxwellian
!Returns a random number in a Maxwellian distribution of mean 0 and width 1 !Returns a random number in a Maxwellian distribution of mean 0 and width 1
FUNCTION randomHalfMaxwellian() RESULT(rnd) FUNCTION randomHalfMaxwellian() RESULT(rnd)

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

@ -829,71 +829,77 @@ MODULE moduleInput
IF (nTypes /= nSpecies) CALL criticalError('Not enough boundary types defined in ' // object, 'readBoundary') IF (nTypes /= nSpecies) CALL criticalError('Not enough boundary types defined in ' // object, 'readBoundary')
ALLOCATE(boundary(i)%bTypes(1:nSpecies)) ALLOCATE(boundary(i)%bTypes(1:nSpecies))
DO s = 1, nSpecies DO s = 1, nSpecies
WRITE(sString,'(i2)') s associate(bound => boundary(i)%bTypes(s)%obj)
object = 'boundary(' // TRIM(iString) // ').bTypes(' // TRIM(sString) // ')' WRITE(sString,'(i2)') s
CALL config%get(object // '.type', bType, found) object = 'boundary(' // TRIM(iString) // ').bTypes(' // TRIM(sString) // ')'
SELECT CASE(bType) CALL config%get(object // '.type', bType, found)
CASE('reflection') SELECT CASE(bType)
ALLOCATE(boundaryReflection:: boundary(i)%bTypes(s)%obj) CASE('reflection')
ALLOCATE(boundaryReflection:: bound)
CASE('absorption') CASE('absorption')
ALLOCATE(boundaryAbsorption:: boundary(i)%bTypes(s)%obj) ALLOCATE(boundaryAbsorption:: bound)
CASE('transparent') CASE('transparent')
ALLOCATE(boundaryTransparent:: boundary(i)%bTypes(s)%obj) ALLOCATE(boundaryTransparent:: bound)
CASE('ionization') CASE('axis')
!Neutral parameters ALLOCATE(boundaryAxis:: bound)
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')
CALL config%get(object // '.effectiveTime', effTime, found) CASE('wallTemperature')
IF (.NOT. found) CALL criticalError("missing parameter 'effectiveTime' for ionization", 'readBoundary') 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 config%get(object // '.energyThreshold', eThreshold, found) CALL initWallTemperature(bound, Tw, cw)
IF (.NOT. found) CALL criticalError("missing parameter 'eThreshold' in ionization", 'readBoundary')
CALL config%get(object // '.crossSection', crossSection, found) CASE('ionization')
IF (.NOT. found) CALL criticalError("missing parameter 'crossSection' for neutrals in ionization", 'readBoundary') !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')
CALL config%get(object // '.electronSecondary', electronSecondary, found) CALL config%get(object // '.effectiveTime', effTime, found)
electronSecondaryID = speciesName2Index(electronSecondary) IF (.NOT. found) CALL criticalError("missing parameter 'effectiveTime' for ionization", 'readBoundary')
IF (found) THEN
CALL initIonization(boundary(i)%bTypes(s)%obj, species(s)%obj%m, m0, n0, v0, T0, &
speciesID, effTime, crossSection, eThreshold,electronSecondaryID)
ELSE CALL config%get(object // '.energyThreshold', eThreshold, found)
CALL initIonization(boundary(i)%bTypes(s)%obj, species(s)%obj%m, m0, n0, v0, T0, & IF (.NOT. found) CALL criticalError("missing parameter 'eThreshold' in ionization", 'readBoundary')
speciesID, effTime, crossSection, eThreshold)
END IF CALL config%get(object // '.crossSection', crossSection, found)
IF (.NOT. found) CALL criticalError("missing parameter 'crossSection' for neutrals in ionization", 'readBoundary')
CASE('wallTemperature') CALL config%get(object // '.electronSecondary', electronSecondary, found)
CALL config%get(object // '.temperature', Tw, found) electronSecondaryID = speciesName2Index(electronSecondary)
IF (.NOT. found) CALL criticalError("temperature not found for wallTemperature boundary type", 'readBoundary') IF (found) THEN
CALL config%get(object // '.specificHeat', cw, found) CALL initIonization(bound, species(s)%obj%m, m0, n0, v0, T0, &
IF (.NOT. found) CALL criticalError("specificHeat not found for wallTemperature boundary type", 'readBoundary') speciesID, effTime, crossSection, eThreshold,electronSecondaryID)
CALL initWallTemperature(boundary(i)%bTypes(s)%obj, Tw, cw) ELSE
CALL initIonization(bound, species(s)%obj%m, m0, n0, v0, T0, &
speciesID, effTime, crossSection, eThreshold)
CASE('axis') END IF
ALLOCATE(boundaryAxis:: boundary(i)%bTypes(s)%obj)
CASE DEFAULT case('outflowAdaptive')
CALL criticalError('Boundary type ' // bType // ' undefined', 'readBoundary') allocate(boundaryOutflowAdaptive:: bound)
END SELECT CASE DEFAULT
CALL criticalError('Boundary type ' // bType // ' undefined', 'readBoundary')
END SELECT
end associate
END DO END DO

49
src/modules/mesh/moduleMeshBoundary.f90
View file

@ -77,6 +77,20 @@ MODULE moduleMeshBoundary
END SUBROUTINE transparent END SUBROUTINE transparent
!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)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
CALL reflection(edge, part)
END SUBROUTINE symmetryAxis
!Wall with temperature !Wall with temperature
SUBROUTINE wallTemperature(edge, part) SUBROUTINE wallTemperature(edge, part)
USE moduleSpecies USE moduleSpecies
@ -204,19 +218,27 @@ MODULE moduleMeshBoundary
END SUBROUTINE ionization END SUBROUTINE ionization
!Symmetry axis. Reflects particles. subroutine outflowAdaptive(edge, part)
!Although this function should never be called, it is set as a reflective boundary use moduleRandom
!to properly deal with possible particles reaching a corner and selecting this boundary. implicit none
SUBROUTINE symmetryAxis(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge class(meshEdge), intent(inout):: edge
CLASS(particle), INTENT(inout):: part class(particle), intent(inout):: part
CALL reflection(edge, part) select type(bound => edge%boundary%bTypes(part%species%n)%obj)
type is(boundaryOutflowAdaptive)
END SUBROUTINE symmetryAxis if (random() < 0.844d0) then
call reflection(edge, part)
else
call transparent(edge, part)
end if
end select
end subroutine outflowAdaptive
!Points the boundary function to specific type !Points the boundary function to specific type
SUBROUTINE pointBoundaryFunction(edge, s) SUBROUTINE pointBoundaryFunction(edge, s)
@ -236,14 +258,17 @@ MODULE moduleMeshBoundary
TYPE IS(boundaryTransparent) TYPE IS(boundaryTransparent)
edge%fBoundary(s)%apply => transparent edge%fBoundary(s)%apply => transparent
TYPE IS(boundaryAxis)
edge%fBoundary(s)%apply => symmetryAxis
TYPE IS(boundaryWallTemperature) TYPE IS(boundaryWallTemperature)
edge%fBoundary(s)%apply => wallTemperature edge%fBoundary(s)%apply => wallTemperature
TYPE IS(boundaryIonization) TYPE IS(boundaryIonization)
edge%fBoundary(s)%apply => ionization edge%fBoundary(s)%apply => ionization
TYPE IS(boundaryAxis) type is(boundaryOutflowAdaptive)
edge%fBoundary(s)%apply => symmetryAxis edge%fBoundary(s)%apply => outflowAdaptive
CLASS DEFAULT CLASS DEFAULT
CALL criticalError("Boundary type not defined in this geometry", 'pointBoundaryFunction') CALL criticalError("Boundary type not defined in this geometry", 'pointBoundaryFunction')

16
src/modules/moduleBoundary.f90
View file

@ -26,6 +26,12 @@ MODULE moduleBoundary
END TYPE boundaryTransparent END TYPE boundaryTransparent
!Symmetry axis
TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryAxis
CONTAINS
END TYPE boundaryAxis
!Wall Temperature boundary !Wall Temperature boundary
TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryWallTemperature TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryWallTemperature
!Thermal velocity of the wall: square root(Wall temperature X specific heat) !Thermal velocity of the wall: square root(Wall temperature X specific heat)
@ -47,11 +53,13 @@ MODULE moduleBoundary
END TYPE boundaryIonization END TYPE boundaryIonization
!Symmetry axis !Boundary for quasi-neutral outflow adjusting reflection coefficient
TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryAxis type, public, extends(boundaryGeneric):: boundaryOutflowAdaptive
CONTAINS real(8):: outflowCurrent
real(8):: reflectionFraction
contains
END TYPE boundaryAxis end type boundaryOutflowAdaptive
!Wrapper for boundary types (one per species) !Wrapper for boundary types (one per species)
TYPE:: bTypesCont TYPE:: bTypesCont

22
src/modules/moduleInject.f90
View file

@ -257,7 +257,7 @@ MODULE moduleInject
CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
REAL(8), INTENT(in):: temperature, m REAL(8), INTENT(in):: temperature, m
velDist = velDistMaxwellian(vTh = DSQRT(temperature/m)) velDist = velDistMaxwellian(vTh = DSQRT(2.d0*temperature/m))
END SUBROUTINE initVelDistMaxwellian END SUBROUTINE initVelDistMaxwellian
@ -267,7 +267,7 @@ MODULE moduleInject
CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
REAL(8), INTENT(in):: temperature, m REAL(8), INTENT(in):: temperature, m
velDist = velDistHalfMaxwellian(vTh = DSQRT(temperature/m)) velDist = velDistHalfMaxwellian(vTh = DSQRT(2.d0*temperature/m))
END SUBROUTINE initVelDistHalfMaxwellian END SUBROUTINE initVelDistHalfMaxwellian
@ -289,7 +289,7 @@ MODULE moduleInject
REAL(8):: v REAL(8):: v
v = 0.D0 v = 0.D0
v = sqrt(2.0)*self%vTh*randomMaxwellian() v = self%vTh*randomMaxwellian()/sqrt(2.d0)
END FUNCTION randomVelMaxwellian END FUNCTION randomVelMaxwellian
@ -302,7 +302,7 @@ MODULE moduleInject
REAL(8):: v REAL(8):: v
v = 0.D0 v = 0.D0
v = sqrt(2.0)*self%vTh*randomHalfMaxwellian() v = self%vTh*randomHalfMaxwellian()/sqrt(2.d0)
END FUNCTION randomVelHalfMaxwellian END FUNCTION randomVelHalfMaxwellian
@ -387,16 +387,12 @@ MODULE moduleInject
partInj(n)%v = 0.D0 partInj(n)%v = 0.D0
partInj(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), & do while(dot_product(partInj(n)%v, direction) <= 0.d0)
self%v(2)%obj%randomVel(), & partInj(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), &
self%v(3)%obj%randomVel() /) self%v(2)%obj%randomVel(), &
self%v(3)%obj%randomVel() /)
end do
!If injecting a no-drift distribution and velocity is negative, reflect
if ((self%vMod == 0.D0) .and. &
(dot_product(direction, partInj(n)%v) < 0.D0)) then
partInj(n)%v = - partInj(n)%v
end if
!Obtain natural coordinates of particle in cell !Obtain natural coordinates of particle in cell
partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r) partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r)