475 lines
14 KiB
Fortran
475 lines
14 KiB
Fortran
!moduleMeshBoundary: Boundary functions for the mesh edges
|
|
submodule(moduleMesh) boundaryParticle
|
|
contains
|
|
module function boundaryParticleName_to_Index(boundaryName) result(bp)
|
|
use moduleErrors
|
|
implicit none
|
|
|
|
character(:), allocatable:: boundaryName
|
|
integer:: bp
|
|
integer:: b
|
|
|
|
bp = 0
|
|
do b = 1, nBoundariesParticle
|
|
if (boundaryName == boundariesParticle(b)%obj%name) then
|
|
bp = boundariesParticle(b)%obj%n
|
|
|
|
end if
|
|
|
|
end do
|
|
|
|
if (bp == 0) then
|
|
call criticalError('Boundary ' // boundaryName // ' not found', 'boundaryParticleName_to_Index')
|
|
|
|
end if
|
|
|
|
end function boundaryParticleName_to_Index
|
|
|
|
module SUBROUTINE initWallTemperature(boundary, T, c)
|
|
USE moduleRefParam
|
|
IMPLICIT NONE
|
|
|
|
CLASS(boundaryParticleGeneric), ALLOCATABLE, INTENT(inout):: boundary
|
|
REAL(8), INTENT(in):: T, c !Wall temperature and specific heat
|
|
REAL(8):: vTh
|
|
|
|
vTh = DSQRT(c * T) / v_ref
|
|
allocate(boundaryWallTemperature:: boundary)
|
|
|
|
select type(boundary)
|
|
type is(boundaryWallTemperature)
|
|
boundary%vTh = vTh
|
|
|
|
end select
|
|
|
|
END SUBROUTINE initWallTemperature
|
|
|
|
module SUBROUTINE initIonization(boundary, mImpact, m0, n0, v0, T0, ion, effTime, crossSection, eThreshold, electronSecondary)
|
|
use moduleRefParam
|
|
use moduleSpecies
|
|
use moduleCaseParam
|
|
use moduleConstParam
|
|
use moduleErrors, only: criticalError
|
|
implicit none
|
|
|
|
class(boundaryParticleGeneric), allocatable, intent(inout):: boundary
|
|
real(8), intent(in):: mImpact
|
|
real(8), intent(in):: m0, n0, v0(1:3), T0 !Neutral properties
|
|
integer, intent(in):: ion
|
|
real(8), intent(in):: effTime
|
|
character(:), allocatable, intent(in):: crossSection
|
|
real(8), intent(in):: eThreshold
|
|
integer, optional, intent(in):: electronSecondary
|
|
|
|
ALLOCATE(boundaryIonization:: boundary)
|
|
|
|
SELECT TYPE(boundary)
|
|
TYPE IS(boundaryIonization)
|
|
boundary%m0 = m0 / m_ref
|
|
boundary%n0 = n0 * Vol_ref
|
|
boundary%v0 = v0 / v_ref
|
|
boundary%vTh = DSQRT(kb*T0/m0)/v_ref
|
|
boundary%species => species(ion)%obj
|
|
IF (PRESENT(electronSecondary)) THEN
|
|
SELECT TYPE(sp => species(electronSecondary)%obj)
|
|
TYPE IS(speciesCharged)
|
|
boundary%electronSecondary => sp
|
|
|
|
CLASS DEFAULT
|
|
CALL criticalError("Species " // sp%name // " chosen for " // &
|
|
"secondary electron is not a charged species", 'initIonization')
|
|
|
|
END SELECT
|
|
|
|
ELSE
|
|
boundary%electronSecondary => NULL()
|
|
|
|
END IF
|
|
boundary%effectiveTime = effTime / ti_ref
|
|
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/mImpact)
|
|
|
|
END SELECT
|
|
|
|
END SUBROUTINE initIonization
|
|
|
|
module subroutine initQuasiNeutrality(boundary, s_incident)
|
|
implicit none
|
|
|
|
class(boundaryParticleGeneric), allocatable, intent(inout):: boundary
|
|
integer, intent(in):: s_incident
|
|
|
|
allocate(boundaryQuasiNeutrality:: boundary)
|
|
|
|
select type(boundary)
|
|
type is(boundaryQuasiNeutrality)
|
|
allocate(boundary%edges(0))
|
|
|
|
boundary%s_incident = s_incident
|
|
|
|
|
|
boundary%update => quasiNeutrality_update
|
|
boundary%print => quasiNeutrality_print
|
|
|
|
end select
|
|
|
|
end subroutine initQuasiNeutrality
|
|
|
|
module SUBROUTINE reflection(self, edge, part)
|
|
USE moduleCaseParam
|
|
USE moduleSpecies
|
|
IMPLICIT NONE
|
|
|
|
class(boundaryReflection), intent(in):: self
|
|
CLASS(meshEdge), INTENT(inout):: edge
|
|
CLASS(particle), INTENT(inout):: part
|
|
!rp = intersection between particle and edge
|
|
!rpp = final position of particle
|
|
!vpp = final velocity of particle
|
|
REAL(8), DIMENSION(1:3):: rp, vpp
|
|
|
|
!Reflect particle velocity
|
|
vpp = part%v - 2.D0*DOT_PRODUCT(part%v, edge%normal)*edge%normal
|
|
part%v = vpp
|
|
|
|
rp = edge%intersection(part%r)
|
|
|
|
part%r = 2.D0*(rp - part%r) + part%r
|
|
|
|
!particle is assumed to be inside
|
|
part%n_in = .TRUE.
|
|
|
|
END SUBROUTINE reflection
|
|
|
|
!Absoption in a surface
|
|
module SUBROUTINE absorption(self, edge, part)
|
|
USE moduleCaseParam
|
|
USE moduleSpecies
|
|
IMPLICIT NONE
|
|
|
|
class(boundaryAbsorption), intent(in):: self
|
|
CLASS(meshEdge), INTENT(inout):: edge
|
|
CLASS(particle), INTENT(inout):: part
|
|
REAL(8):: rpp(1:3) !Position of particle projected to the edge
|
|
REAL(8):: d !Distance from particle to edge
|
|
|
|
rpp = edge%intersection(part%r)
|
|
|
|
d = NORM2(rpp - part%r)
|
|
|
|
IF (d >= 0.D0) THEN
|
|
part%weight = part%weight/d
|
|
|
|
END IF
|
|
|
|
!Assign new position to particle
|
|
part%r = rpp
|
|
!Remove particle from the domain
|
|
part%n_in = .FALSE.
|
|
|
|
!Scatter particle in associated volume
|
|
IF (ASSOCIATED(edge%e1)) THEN
|
|
CALL edge%e1%scatter(edge%e1%nNodes, part)
|
|
|
|
ELSE
|
|
CALL edge%e2%scatter(edge%e2%nNodes, part)
|
|
|
|
END IF
|
|
|
|
END SUBROUTINE absorption
|
|
|
|
!Transparent boundary condition
|
|
module SUBROUTINE transparent(self, edge, part)
|
|
USE moduleSpecies
|
|
IMPLICIT NONE
|
|
|
|
class(boundaryTransparent), intent(in):: self
|
|
CLASS(meshEdge), INTENT(inout):: edge
|
|
CLASS(particle), INTENT(inout):: part
|
|
|
|
!Removes particle from domain
|
|
part%n_in = .FALSE.
|
|
|
|
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.
|
|
module SUBROUTINE symmetryAxis(self, edge, part)
|
|
USE moduleSpecies
|
|
IMPLICIT NONE
|
|
|
|
class(boundaryAxis), intent(in):: self
|
|
CLASS(meshEdge), INTENT(inout):: edge
|
|
CLASS(particle), INTENT(inout):: part
|
|
|
|
CALL genericReflection(edge, part)
|
|
|
|
END SUBROUTINE symmetryAxis
|
|
|
|
!Wall with temperature
|
|
module SUBROUTINE wallTemperature(self, edge, part)
|
|
USE moduleSpecies
|
|
USE moduleRandom
|
|
IMPLICIT NONE
|
|
|
|
class(boundaryWallTemperature), intent(in):: self
|
|
CLASS(meshEdge), INTENT(inout):: edge
|
|
CLASS(particle), INTENT(inout):: part
|
|
INTEGER:: i
|
|
|
|
!Modifies particle velocity according to wall temperature
|
|
DO i = 1, 3
|
|
part%v(i) = part%v(i) + self%vTh*randomMaxwellian()
|
|
|
|
END DO
|
|
|
|
CALL genericReflection(edge, part)
|
|
|
|
END SUBROUTINE wallTemperature
|
|
|
|
!Ionization surface: an electron will pass through the surface
|
|
! and create an ion-electron pair based on a neutral background
|
|
module SUBROUTINE ionization(self, edge, part)
|
|
USE moduleList
|
|
USE moduleSpecies
|
|
USE moduleMesh
|
|
USE moduleRefParam
|
|
USE moduleRandom
|
|
USE moduleMath
|
|
IMPLICIT NONE
|
|
|
|
class(boundaryIonization), intent(in):: self
|
|
CLASS(meshEdge), INTENT(inout):: edge
|
|
CLASS(particle), INTENT(inout):: part
|
|
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
|
|
INTEGER:: p
|
|
REAL(8):: v0(1:3) !random velocity of neutral
|
|
TYPE(particle), POINTER:: newElectron
|
|
TYPE(particle), POINTER:: newIon
|
|
|
|
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/self%eThreshold)
|
|
|
|
DO p = 1, nIonizations
|
|
!Get probability of ionization
|
|
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) = 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(self%electronSecondary)) THEN
|
|
newElectron%species => self%electronSecondary
|
|
|
|
ELSE
|
|
newElectron%species => part%species
|
|
|
|
END IF
|
|
newIon%species => self%species
|
|
|
|
newElectron%v = v0 + (1.D0 + self%deltaV*v0/NORM2(v0))
|
|
newIon%v = v0
|
|
|
|
newElectron%r = edge%randPos()
|
|
newIon%r = newElectron%r
|
|
|
|
IF (ASSOCIATED(edge%e1)) THEN
|
|
newElectron%cell = edge%e1%n
|
|
|
|
ELSEIF (ASSOCIATED(edge%e2)) THEN
|
|
newElectron%cell = edge%e2%n
|
|
|
|
END IF
|
|
newIon%cell = newElectron%cell
|
|
|
|
newElectron%Xi = mesh%cells(part%cell)%obj%phy2log(newElectron%r)
|
|
newIon%Xi = newElectron%Xi
|
|
|
|
newElectron%weight = part%weight
|
|
newIon%weight = newElectron%weight
|
|
|
|
newElectron%n_in = .TRUE.
|
|
newIon%n_in = .TRUE.
|
|
|
|
!Add particles to list
|
|
CALL partSurfaces%setLock()
|
|
CALL partSurfaces%add(newElectron)
|
|
CALL partSurfaces%add(newIon)
|
|
CALL partSurfaces%unsetLock()
|
|
|
|
!Electron loses energy due to ionization
|
|
eRel = eRel - self%eThreshold
|
|
vRel = 2.D0*DSQRT(eRel)/mRel
|
|
|
|
!Reduce number of possible ionizations
|
|
nIonizations = nIonizations - 1
|
|
|
|
END IF
|
|
|
|
END DO
|
|
|
|
!Removes ionizing electron regardless the number of pair created
|
|
part%n_in = .FALSE.
|
|
|
|
END SUBROUTINE ionization
|
|
|
|
module subroutine quasiNeutrality(self, edge, part)
|
|
use moduleRandom
|
|
implicit none
|
|
|
|
class(boundaryQuasiNeutrality), intent(in):: self
|
|
class(meshEdge), intent(inout):: edge
|
|
class(particle), intent(inout):: part
|
|
|
|
if (random() <= self%alpha(edge%n)) then
|
|
call genericReflection(edge, part)
|
|
|
|
else
|
|
call genericTransparent(edge, part)
|
|
|
|
end if
|
|
|
|
end subroutine quasiNeutrality
|
|
|
|
subroutine quasiNeutrality_update(self)
|
|
implicit none
|
|
|
|
class(boundaryParticleGeneric), intent(inout):: self
|
|
integer:: e, n, s
|
|
integer, allocatable:: nodes(:)
|
|
class(meshEdge), pointer:: edge
|
|
real(8), allocatable:: density_nodes(:)
|
|
class(meshNode), pointer:: node
|
|
real(8):: density_incident, density_rest
|
|
real(8):: alpha
|
|
|
|
select type(self)
|
|
type is(boundaryQuasiNeutrality)
|
|
do e = 1, size(self%edges)
|
|
edge => mesh%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(n)%obj
|
|
density_nodes(n) = node%output(s)%den
|
|
|
|
end do
|
|
|
|
if (s == self%s_incident) then
|
|
density_incident = edge%gatherF((/0.d0,0.d0,0.d0/), edge%nNodes, density_nodes)
|
|
|
|
else
|
|
density_rest = density_rest + edge%gatherF((/0.d0,0.d0,0.d0/), edge%nNodes, density_nodes)
|
|
|
|
end if
|
|
|
|
end do
|
|
|
|
! Correction for this time step
|
|
alpha = 1.d0 - density_incident/density_rest
|
|
|
|
! Apply correction with a factor of 0.1 to avoid fast changes
|
|
self%alpha(edge%n) = self%alpha(edge%n) + 1.0d-2 * alpha
|
|
|
|
! Limit alpha between 0 and 1
|
|
self%alpha(edge%n) = min(self%alpha(edge%n), 1.d0)
|
|
self%alpha(edge%n) = max(self%alpha(edge%n), 0.d0)
|
|
|
|
deallocate(density_nodes)
|
|
|
|
end do
|
|
|
|
end select
|
|
|
|
end subroutine quasiNeutrality_update
|
|
|
|
subroutine quasiNeutrality_print(self)
|
|
implicit none
|
|
|
|
class(boundaryParticleGeneric), intent(inout):: self
|
|
|
|
print *, 'test'
|
|
|
|
select type(self)
|
|
type is(boundaryQuasiNeutrality)
|
|
print*, self%alpha
|
|
|
|
end select
|
|
|
|
end subroutine quasiNeutrality_print
|
|
|
|
! Generic boundary conditions for internal use
|
|
module subroutine genericReflection(edge, part)
|
|
use moduleCaseParam
|
|
use moduleSpecies
|
|
implicit none
|
|
|
|
class(meshEdge), intent(inout):: edge
|
|
class(particle), intent(inout):: part
|
|
!rp = intersection between particle and edge
|
|
!rpp = final position of particle
|
|
!vpp = final velocity of particle
|
|
real(8), dimension(1:3):: rp, vpp
|
|
|
|
!Reflect particle velocity
|
|
vpp = part%v - 2.D0*dot_product(part%v, edge%normal)*edge%normal
|
|
part%v = vpp
|
|
|
|
rp = edge%intersection(part%r)
|
|
|
|
part%r = 2.D0*(rp - part%r) + part%r
|
|
|
|
!particle is assumed to be inside
|
|
part%n_in = .TRUE.
|
|
|
|
end subroutine genericReflection
|
|
|
|
module subroutine genericTransparent(edge, part)
|
|
use moduleSpecies
|
|
implicit none
|
|
|
|
class(meshEdge), intent(inout):: edge
|
|
class(particle), intent(inout):: part
|
|
|
|
!Removes particle from domain
|
|
part%n_in = .FALSE.
|
|
|
|
end subroutine genericTransparent
|
|
|
|
module subroutine boundariesParticle_update()
|
|
implicit none
|
|
|
|
integer:: b
|
|
|
|
do b = 1, nBoundariesParticle
|
|
if (associated(boundariesParticle(b)%obj%update)) then
|
|
call boundariesParticle(b)%obj%update
|
|
|
|
end if
|
|
|
|
end do
|
|
|
|
end subroutine boundariesParticle_update
|
|
|
|
end submodule boundaryParticle
|