diff --git a/src/modules/common/moduleRandom.f90 b/src/modules/common/moduleRandom.f90 index 1b8ec67..156f5e4 100644 --- a/src/modules/common/moduleRandom.f90 +++ b/src/modules/common/moduleRandom.f90 @@ -48,23 +48,24 @@ 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 - v1 = 0.d0 - do while (v1 <= 0.d0) - v1 = random() + 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 end do - v2 = random() - rnd = sqrt(-2.d0*log(v1))*cos(2*pi*v2) + rnd = v2 * sqrt(-2.D0 * log(Rsquare) / Rsquare) - 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 9f50947..74da279 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -829,77 +829,71 @@ 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 - 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) + 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) - CASE('absorption') - ALLOCATE(boundaryAbsorption:: bound) + CASE('absorption') + ALLOCATE(boundaryAbsorption:: boundary(i)%bTypes(s)%obj) - CASE('transparent') - ALLOCATE(boundaryTransparent:: bound) + CASE('transparent') + ALLOCATE(boundaryTransparent:: boundary(i)%bTypes(s)%obj) - CASE('axis') - ALLOCATE(boundaryAxis:: 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('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 // '.effectiveTime', effTime, found) + IF (.NOT. found) CALL criticalError("missing parameter 'effectiveTime' for ionization", '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') - !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 // '.crossSection', crossSection, found) + IF (.NOT. found) CALL criticalError("missing parameter 'crossSection' 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 // '.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 // '.energyThreshold', eThreshold, found) - IF (.NOT. found) CALL criticalError("missing parameter 'eThreshold' in 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 // '.crossSection', crossSection, found) - IF (.NOT. found) CALL criticalError("missing parameter 'crossSection' for neutrals in ionization", 'readBoundary') + END IF - 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) + 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') - ELSE - CALL initIonization(bound, species(s)%obj%m, m0, n0, v0, T0, & - speciesID, effTime, crossSection, eThreshold) + CALL initWallTemperature(boundary(i)%bTypes(s)%obj, Tw, cw) - END IF + CASE('axis') + ALLOCATE(boundaryAxis:: boundary(i)%bTypes(s)%obj) - case('outflowAdaptive') - allocate(boundaryOutflowAdaptive:: bound) + 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 diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90 index deb4459..091e52e 100644 --- a/src/modules/mesh/moduleMeshBoundary.f90 +++ b/src/modules/mesh/moduleMeshBoundary.f90 @@ -77,20 +77,6 @@ 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 @@ -218,27 +204,19 @@ MODULE moduleMeshBoundary END SUBROUTINE ionization - subroutine outflowAdaptive(edge, part) - use moduleRandom - implicit none + !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 + CLASS(meshEdge), INTENT(inout):: edge + CLASS(particle), INTENT(inout):: part - select type(bound => edge%boundary%bTypes(part%species%n)%obj) - type is(boundaryOutflowAdaptive) + CALL reflection(edge, part) - if (random() < 0.844d0) then - call reflection(edge, part) - - else - call transparent(edge, part) - - end if - - end select - - end subroutine outflowAdaptive + END SUBROUTINE symmetryAxis !Points the boundary function to specific type SUBROUTINE pointBoundaryFunction(edge, s) @@ -258,17 +236,14 @@ 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(boundaryOutflowAdaptive) - edge%fBoundary(s)%apply => outflowAdaptive + TYPE IS(boundaryAxis) + edge%fBoundary(s)%apply => symmetryAxis 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 67465c7..0b76105 100644 --- a/src/modules/moduleBoundary.f90 +++ b/src/modules/moduleBoundary.f90 @@ -26,12 +26,6 @@ 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) @@ -53,13 +47,11 @@ MODULE moduleBoundary END TYPE boundaryIonization - !Boundary for quasi-neutral outflow adjusting reflection coefficient - type, public, extends(boundaryGeneric):: boundaryOutflowAdaptive - real(8):: outflowCurrent - real(8):: reflectionFraction - contains + !Symmetry axis + TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryAxis + CONTAINS - end type boundaryOutflowAdaptive + END TYPE boundaryAxis !Wrapper for boundary types (one per species) TYPE:: bTypesCont diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index ff8a694..528cb91 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(2.d0*temperature/m)) + velDist = velDistMaxwellian(vTh = DSQRT(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(2.d0*temperature/m)) + velDist = velDistHalfMaxwellian(vTh = DSQRT(temperature/m)) END SUBROUTINE initVelDistHalfMaxwellian @@ -289,7 +289,7 @@ MODULE moduleInject REAL(8):: v v = 0.D0 - v = self%vTh*randomMaxwellian()/sqrt(2.d0) + v = sqrt(2.0)*self%vTh*randomMaxwellian() END FUNCTION randomVelMaxwellian @@ -302,7 +302,7 @@ MODULE moduleInject REAL(8):: v v = 0.D0 - v = self%vTh*randomHalfMaxwellian()/sqrt(2.d0) + v = sqrt(2.0)*self%vTh*randomHalfMaxwellian() END FUNCTION randomVelHalfMaxwellian @@ -387,12 +387,16 @@ MODULE moduleInject 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(), & - self%v(2)%obj%randomVel(), & - self%v(3)%obj%randomVel() /) - end do + partInj(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), & + self%v(2)%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 !Obtain natural coordinates of particle in cell partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r)