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/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index a304c5e..bbc72e2 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -104,6 +104,7 @@ 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 db1ef4f..dcd2128 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -144,6 +144,7 @@ MODULE moduleMesh2DCart USE moduleSpecies USE moduleBoundary USE moduleErrors + USE moduleRefParam, ONLY: L_ref IMPLICIT NONE CLASS(meshEdge2DCart), INTENT(out):: self @@ -163,7 +164,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) + self%surface = SQRT((self%x(2) - self%x(1))**2 + (self%y(2) - self%y(1))**2) / L_ref !Normal vector self%normal = (/ -(self%y(2)-self%y(1)), & self%x(2)-self%x(1) , & @@ -558,6 +559,7 @@ 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 7fc7798..655e4df 100644 --- a/src/modules/mesh/moduleMeshBoundary.f90 +++ b/src/modules/mesh/moduleMeshBoundary.f90 @@ -218,37 +218,6 @@ MODULE moduleMeshBoundary 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 SUBROUTINE pointBoundaryFunction(edge, s) USE moduleErrors @@ -267,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..278870f 100644 --- a/src/modules/moduleBoundary.f90 +++ b/src/modules/moduleBoundary.f90 @@ -56,7 +56,6 @@ 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 0d761ea..528cb91 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -387,18 +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() /) - !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 + partInj(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), & + self%v(2)%obj%randomVel(), & + self%v(3)%obj%randomVel() /) - 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 partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r)