From 159f2e7620a355549ee64c8bbda579a0f94768c4 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Tue, 3 Feb 2026 10:15:12 +0100 Subject: [PATCH 1/4] Trying partial reflection --- src/modules/init/moduleInput.f90 | 108 ++++++++++--------- src/modules/mesh/1DCart/moduleMesh1DCart.f90 | 1 - src/modules/mesh/2DCart/moduleMesh2DCart.f90 | 4 +- src/modules/mesh/moduleMeshBoundary.f90 | 29 ++++- src/modules/moduleBoundary.f90 | 1 + src/modules/moduleInject.f90 | 8 +- 6 files changed, 91 insertions(+), 60 deletions(-) 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. & From 2b6ea0ff8494693d9f5f7e5bfb3e5c3ec80472fa Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Tue, 3 Feb 2026 11:48:43 +0100 Subject: [PATCH 2/4] First attempt at working with a distribution function with u --- src/modules/mesh/moduleMeshBoundary.f90 | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90 index deb4459..eaa2d52 100644 --- a/src/modules/mesh/moduleMeshBoundary.f90 +++ b/src/modules/mesh/moduleMeshBoundary.f90 @@ -220,6 +220,7 @@ MODULE moduleMeshBoundary subroutine outflowAdaptive(edge, part) use moduleRandom + use moduleRefParam, only: v_ref implicit none class(meshEdge), intent(inout):: edge @@ -228,14 +229,19 @@ MODULE moduleMeshBoundary select type(bound => edge%boundary%bTypes(part%species%n)%obj) type is(boundaryOutflowAdaptive) - if (random() < 0.844d0) then - call reflection(edge, part) - - else + ! if (random() < 0.844d0) then + call reflection(edge, part) + part%v = part%v + 40e3/v_ref + if (dot_product(part%v, edge%normal) <= 0.d0) then call transparent(edge, part) end if + ! else + ! call transparent(edge, part) + + ! end if + end select end subroutine outflowAdaptive From 8ae558f3c84dbe4b0044cc6cf85587dfe7f66d10 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Wed, 4 Feb 2026 11:23:03 +0100 Subject: [PATCH 3/4] Correction to ensure always correct direction in the injection --- src/modules/moduleInject.f90 | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index cb1c7fb..0d761ea 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -391,15 +391,15 @@ MODULE moduleInject 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 + 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) !Push new particle with the minimum time step From 4aa6e5aab2e60bffa1c6548869b6e2abf20feb63 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Wed, 4 Feb 2026 11:23:26 +0100 Subject: [PATCH 4/4] Try to impose a reflection with drag --- src/modules/mesh/moduleMeshBoundary.f90 | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90 index eaa2d52..7fc7798 100644 --- a/src/modules/mesh/moduleMeshBoundary.f90 +++ b/src/modules/mesh/moduleMeshBoundary.f90 @@ -225,23 +225,26 @@ MODULE moduleMeshBoundary 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 * - ! if (random() < 0.844d0) then - call reflection(edge, part) - part%v = part%v + 40e3/v_ref - if (dot_product(part%v, edge%normal) <= 0.d0) then + else call transparent(edge, part) end if - ! else - ! call transparent(edge, part) - - ! end if - end select end subroutine outflowAdaptive