Compare commits

..

No commits in common. "4aa6e5aab2e60bffa1c6548869b6e2abf20feb63" and "4747673d0a81af1ee55dad7bb197f595028d0625" have entirely different histories.

6 changed files with 65 additions and 105 deletions

View file

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

View file

@ -104,6 +104,7 @@ MODULE moduleMesh1DCart
USE moduleSpecies USE moduleSpecies
USE moduleBoundary USE moduleBoundary
USE moduleErrors USE moduleErrors
USE moduleRefParam, ONLY: L_ref
IMPLICIT NONE IMPLICIT NONE
CLASS(meshEdge1DCart), INTENT(out):: self CLASS(meshEdge1DCart), INTENT(out):: self

View file

@ -144,6 +144,7 @@ MODULE moduleMesh2DCart
USE moduleSpecies USE moduleSpecies
USE moduleBoundary USE moduleBoundary
USE moduleErrors USE moduleErrors
USE moduleRefParam, ONLY: L_ref
IMPLICIT NONE IMPLICIT NONE
CLASS(meshEdge2DCart), INTENT(out):: self CLASS(meshEdge2DCart), INTENT(out):: self
@ -163,7 +164,7 @@ MODULE moduleMesh2DCart
r2 = self%n2%getCoordinates() r2 = self%n2%getCoordinates()
self%x = (/r1(1), r2(1)/) self%x = (/r1(1), r2(1)/)
self%y = (/r1(2), r2(2)/) self%y = (/r1(2), r2(2)/)
self%surface = SQRT((self%x(2) - self%x(1))**2 + (self%y(2) - self%y(1))**2) self%surface = SQRT((self%x(2) - self%x(1))**2 + (self%y(2) - self%y(1))**2) / L_ref
!Normal vector !Normal vector
self%normal = (/ -(self%y(2)-self%y(1)), & self%normal = (/ -(self%y(2)-self%y(1)), &
self%x(2)-self%x(1) , & self%x(2)-self%x(1) , &
@ -558,6 +559,7 @@ MODULE moduleMesh2DCart
!Compute element volume !Compute element volume
PURE SUBROUTINE volumeQuad(self) PURE SUBROUTINE volumeQuad(self)
USE moduleRefParam, ONLY: L_ref
IMPLICIT NONE IMPLICIT NONE
CLASS(meshCell2DCartQuad), INTENT(inout):: self CLASS(meshCell2DCartQuad), INTENT(inout):: self

View file

@ -218,37 +218,6 @@ MODULE moduleMeshBoundary
END SUBROUTINE ionization END SUBROUTINE ionization
subroutine outflowAdaptive(edge, part)
use moduleRandom
use moduleRefParam, only: v_ref
implicit none
class(meshEdge), intent(inout):: edge
class(particle), intent(inout):: part
real(8):: v_cut
v_cut = 2.d0 * 40.d3/v_ref !This will be the drag velocity of the ions in the future
select type(bound => edge%boundary%bTypes(part%species%n)%obj)
type is(boundaryOutflowAdaptive)
if (dot_product(part%v,-edge%normal) > v_cut) then
! print *,part%v(1), v_cut
! part%v = part%v + v_cut*edge%normal
part%v(1) = part%v(1) - v_cut
! print *,part%v(1)
call reflection(edge, part)
! print *,part%v(1)
! print *
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)
USE moduleErrors USE moduleErrors
@ -267,17 +236,14 @@ 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(boundaryOutflowAdaptive) TYPE IS(boundaryAxis)
edge%fBoundary(s)%apply => outflowAdaptive edge%fBoundary(s)%apply => symmetryAxis
CLASS DEFAULT CLASS DEFAULT
CALL criticalError("Boundary type not defined in this geometry", 'pointBoundaryFunction') CALL criticalError("Boundary type not defined in this geometry", 'pointBoundaryFunction')

View file

@ -56,7 +56,6 @@ MODULE moduleBoundary
!Boundary for quasi-neutral outflow adjusting reflection coefficient !Boundary for quasi-neutral outflow adjusting reflection coefficient
type, public, extends(boundaryGeneric):: boundaryOutflowAdaptive type, public, extends(boundaryGeneric):: boundaryOutflowAdaptive
real(8):: outflowCurrent real(8):: outflowCurrent
real(8):: reflectionFraction
contains contains
end type boundaryOutflowAdaptive end type boundaryOutflowAdaptive

View file

@ -387,18 +387,16 @@ MODULE moduleInject
partInj(n)%v = 0.D0 partInj(n)%v = 0.D0
do while(dot_product(partInj(n)%v, direction) <= 0.d0) partInj(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), &
partInj(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), & self%v(2)%obj%randomVel(), &
self%v(2)%obj%randomVel(), & self%v(3)%obj%randomVel() /)
self%v(3)%obj%randomVel() /)
!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 !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 do 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)