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/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index bbc72e2..a304c5e 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -104,7 +104,6 @@ MODULE moduleMesh1DCart USE moduleSpecies USE moduleBoundary USE moduleErrors - USE moduleRefParam, ONLY: L_ref IMPLICIT NONE CLASS(meshEdge1DCart), INTENT(out):: self diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index dcd2128..db1ef4f 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -144,7 +144,6 @@ MODULE moduleMesh2DCart USE moduleSpecies USE moduleBoundary USE moduleErrors - USE moduleRefParam, ONLY: L_ref IMPLICIT NONE CLASS(meshEdge2DCart), INTENT(out):: self @@ -164,7 +163,7 @@ MODULE moduleMesh2DCart r2 = self%n2%getCoordinates() self%x = (/r1(1), r2(1)/) self%y = (/r1(2), r2(2)/) - self%surface = SQRT((self%x(2) - self%x(1))**2 + (self%y(2) - self%y(1))**2) / L_ref + self%surface = SQRT((self%x(2) - self%x(1))**2 + (self%y(2) - self%y(1))**2) !Normal vector self%normal = (/ -(self%y(2)-self%y(1)), & self%x(2)-self%x(1) , & @@ -559,7 +558,6 @@ MODULE moduleMesh2DCart !Compute element volume PURE SUBROUTINE volumeQuad(self) - USE moduleRefParam, ONLY: L_ref IMPLICIT NONE CLASS(meshCell2DCartQuad), INTENT(inout):: self diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90 index 655e4df..deb4459 100644 --- a/src/modules/mesh/moduleMeshBoundary.f90 +++ b/src/modules/mesh/moduleMeshBoundary.f90 @@ -218,6 +218,28 @@ MODULE moduleMeshBoundary END SUBROUTINE ionization + subroutine outflowAdaptive(edge, part) + use moduleRandom + implicit none + + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part + + select type(bound => edge%boundary%bTypes(part%species%n)%obj) + type is(boundaryOutflowAdaptive) + + 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) USE moduleErrors @@ -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 278870f..67465c7 100644 --- a/src/modules/moduleBoundary.f90 +++ b/src/modules/moduleBoundary.f90 @@ -56,6 +56,7 @@ MODULE moduleBoundary !Boundary for quasi-neutral outflow adjusting reflection coefficient type, public, extends(boundaryGeneric):: boundaryOutflowAdaptive real(8):: outflowCurrent + real(8):: reflectionFraction contains end type boundaryOutflowAdaptive diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index 528cb91..cb1c7fb 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -387,9 +387,11 @@ 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. &