diff --git a/src/modules/common/moduleRandom.f90 b/src/modules/common/moduleRandom.f90 index 156f5e4..1b8ec67 100644 --- a/src/modules/common/moduleRandom.f90 +++ b/src/modules/common/moduleRandom.f90 @@ -48,24 +48,23 @@ MODULE moduleRandom END FUNCTION randomIntAB !Returns a random number in a Maxwellian distribution of mean 0 and width 1 with the Box-Muller Method - FUNCTION randomMaxwellian() RESULT(rnd) - USE moduleConstParam, ONLY: PI - IMPLICIT NONE + function randomMaxwellian() result(rnd) + USE moduleConstParam, only: pi + implicit none - REAL(8):: rnd - REAL(8):: v1, v2, Rsquare + real(8):: rnd + real(8):: v1, v2, Rsquare - Rsquare = 1.D0 - do while (Rsquare >= 1.D0 .and. Rsquare > 0.D0) - v1 = 2.D0 * random() - 1.D0 - v2 = 2.D0 * random() - 1.D0 - Rsquare = v1**2 + v2**2 + v1 = 0.d0 + do while (v1 <= 0.d0) + v1 = random() 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 FUNCTION randomHalfMaxwellian() RESULT(rnd) diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 74da279..9f50947 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -829,71 +829,77 @@ MODULE moduleInput IF (nTypes /= nSpecies) CALL criticalError('Not enough boundary types defined in ' // object, 'readBoundary') ALLOCATE(boundary(i)%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) - SELECT CASE(bType) - CASE('reflection') - ALLOCATE(boundaryReflection:: boundary(i)%bTypes(s)%obj) + associate(bound => boundary(i)%bTypes(s)%obj) + WRITE(sString,'(i2)') s + object = 'boundary(' // TRIM(iString) // ').bTypes(' // TRIM(sString) // ')' + CALL config%get(object // '.type', bType, found) + SELECT CASE(bType) + CASE('reflection') + ALLOCATE(boundaryReflection:: bound) - CASE('absorption') - ALLOCATE(boundaryAbsorption:: boundary(i)%bTypes(s)%obj) + CASE('absorption') + ALLOCATE(boundaryAbsorption:: bound) - CASE('transparent') - ALLOCATE(boundaryTransparent:: boundary(i)%bTypes(s)%obj) + CASE('transparent') + ALLOCATE(boundaryTransparent:: bound) - 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('axis') + ALLOCATE(boundaryAxis:: bound) - CALL config%get(object // '.effectiveTime', effTime, found) - IF (.NOT. found) CALL criticalError("missing parameter 'effectiveTime' for ionization", '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 config%get(object // '.energyThreshold', eThreshold, found) - IF (.NOT. found) CALL criticalError("missing parameter 'eThreshold' in ionization", 'readBoundary') + CALL initWallTemperature(bound, Tw, cw) - CALL config%get(object // '.crossSection', crossSection, found) - IF (.NOT. found) CALL criticalError("missing parameter 'crossSection' 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) 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) - 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 // '.effectiveTime', effTime, found) + IF (.NOT. found) CALL criticalError("missing parameter 'effectiveTime' for ionization", 'readBoundary') - ELSE - CALL initIonization(boundary(i)%bTypes(s)%obj, species(s)%obj%m, m0, n0, v0, T0, & - speciesID, effTime, crossSection, eThreshold) + CALL config%get(object // '.energyThreshold', eThreshold, found) + IF (.NOT. found) CALL criticalError("missing parameter 'eThreshold' in ionization", 'readBoundary') - 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 // '.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 // '.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 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') - ALLOCATE(boundaryAxis:: boundary(i)%bTypes(s)%obj) + END IF - CASE DEFAULT - CALL criticalError('Boundary type ' // bType // ' undefined', 'readBoundary') + case('outflowAdaptive') + allocate(boundaryOutflowAdaptive:: bound) - END SELECT + CASE DEFAULT + CALL criticalError('Boundary type ' // bType // ' undefined', 'readBoundary') + + END SELECT + + end associate END DO diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90 index 091e52e..deb4459 100644 --- a/src/modules/mesh/moduleMeshBoundary.f90 +++ b/src/modules/mesh/moduleMeshBoundary.f90 @@ -77,6 +77,20 @@ MODULE moduleMeshBoundary 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 SUBROUTINE wallTemperature(edge, part) USE moduleSpecies @@ -204,19 +218,27 @@ MODULE moduleMeshBoundary END SUBROUTINE ionization - !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 + subroutine outflowAdaptive(edge, part) + use moduleRandom + implicit none - CLASS(meshEdge), INTENT(inout):: edge - CLASS(particle), INTENT(inout):: part + class(meshEdge), intent(inout):: edge + 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 SUBROUTINE pointBoundaryFunction(edge, s) @@ -236,14 +258,17 @@ MODULE moduleMeshBoundary 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(boundaryAxis) - edge%fBoundary(s)%apply => symmetryAxis + type is(boundaryOutflowAdaptive) + edge%fBoundary(s)%apply => outflowAdaptive CLASS DEFAULT CALL criticalError("Boundary type not defined in this geometry", 'pointBoundaryFunction') diff --git a/src/modules/moduleBoundary.f90 b/src/modules/moduleBoundary.f90 index 0b76105..67465c7 100644 --- a/src/modules/moduleBoundary.f90 +++ b/src/modules/moduleBoundary.f90 @@ -26,6 +26,12 @@ MODULE moduleBoundary END TYPE boundaryTransparent + !Symmetry axis + TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryAxis + CONTAINS + + END TYPE boundaryAxis + !Wall Temperature boundary TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryWallTemperature !Thermal velocity of the wall: square root(Wall temperature X specific heat) @@ -47,11 +53,13 @@ MODULE moduleBoundary END TYPE boundaryIonization - !Symmetry axis - TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryAxis - CONTAINS + !Boundary for quasi-neutral outflow adjusting reflection coefficient + type, public, extends(boundaryGeneric):: boundaryOutflowAdaptive + real(8):: outflowCurrent + real(8):: reflectionFraction + contains - END TYPE boundaryAxis + end type boundaryOutflowAdaptive !Wrapper for boundary types (one per species) TYPE:: bTypesCont diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index 528cb91..ff8a694 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -257,7 +257,7 @@ MODULE moduleInject CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist REAL(8), INTENT(in):: temperature, m - velDist = velDistMaxwellian(vTh = DSQRT(temperature/m)) + velDist = velDistMaxwellian(vTh = DSQRT(2.d0*temperature/m)) END SUBROUTINE initVelDistMaxwellian @@ -267,7 +267,7 @@ MODULE moduleInject CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist REAL(8), INTENT(in):: temperature, m - velDist = velDistHalfMaxwellian(vTh = DSQRT(temperature/m)) + velDist = velDistHalfMaxwellian(vTh = DSQRT(2.d0*temperature/m)) END SUBROUTINE initVelDistHalfMaxwellian @@ -289,7 +289,7 @@ MODULE moduleInject REAL(8):: v v = 0.D0 - v = sqrt(2.0)*self%vTh*randomMaxwellian() + v = self%vTh*randomMaxwellian()/sqrt(2.d0) END FUNCTION randomVelMaxwellian @@ -302,7 +302,7 @@ MODULE moduleInject REAL(8):: v v = 0.D0 - v = sqrt(2.0)*self%vTh*randomHalfMaxwellian() + v = self%vTh*randomHalfMaxwellian()/sqrt(2.d0) END FUNCTION randomVelHalfMaxwellian @@ -387,16 +387,12 @@ MODULE moduleInject partInj(n)%v = 0.D0 - partInj(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), & - self%v(2)%obj%randomVel(), & - self%v(3)%obj%randomVel() /) + do while(dot_product(partInj(n)%v, direction) <= 0.d0) + partInj(n)%v = self%vMod*direction + (/ self%v(1)%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 partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r)