From e54c86a952c3d3cd8a92200375c6ce696d5b9a4a Mon Sep 17 00:00:00 2001
From: JGonzalez <jorge.gonzalez@upm.es>
Date: Thu, 12 Mar 2026 20:09:34 +0100
Subject: [PATCH] Recovered the ouflowBoundary code

---
 src/modules/init/moduleInput.f90              |  3 ++
 src/modules/mesh/moduleMesh.f90               | 19 +++++++++++-
 .../mesh/moduleMesh@boundaryParticle.f90      | 30 +++++++++++++++++++
 src/modules/moduleInject.f90                  |  2 +-
 4 files changed, 52 insertions(+), 2 deletions(-)

diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90
index f3d2310..34f2466 100644
--- a/src/modules/init/moduleInput.f90
+++ b/src/modules/init/moduleInput.f90
@@ -900,6 +900,9 @@ MODULE moduleInput
 
               call initQuasiNeutrality(bound, s_incident)
 
+            CASE('outflowAdaptive')
+              ALLOCATE(boundaryOutflowAdaptive:: bound)
+
             CASE DEFAULT
               CALL criticalError('Boundary type ' // bType // ' undefined', 'readBoundary')
 
diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90
index 5c050e9..71c9bce 100644
--- a/src/modules/mesh/moduleMesh.f90
+++ b/src/modules/mesh/moduleMesh.f90
@@ -789,7 +789,15 @@ MODULE moduleMesh
 
   end type boundaryQuasiNeutrality
 
-  !Wrapper for boundary types (one per species)
+  !Boundary for quasi-neutral outflow adjusting reflection coefficient
+  type, public, extends(boundaryParticleGeneric):: boundaryOutflowAdaptive
+    real(8):: outflowCurrent
+    real(8):: reflectionFraction
+    contains
+      procedure, pass:: apply => outflowAdaptive
+
+  end type boundaryOutflowAdaptive
+
   interface
     module subroutine reflection(self, edge, part)
       use moduleSpecies
@@ -854,6 +862,15 @@ MODULE moduleMesh
 
     end subroutine quasiNeutrality
 
+    module subroutine outflowAdaptive(self, edge, part)
+      use moduleSpecies
+
+      class(boundaryOutflowAdaptive), intent(inout):: self
+      class(meshEdge), intent(inout):: edge
+      class(particle), intent(inout):: part
+
+    end subroutine outflowAdaptive
+
     ! Generic basic boundary conditions to use internally in the code
     module subroutine genericReflection(edge, part) 
       use moduleSpecies
diff --git a/src/modules/mesh/moduleMesh@boundaryParticle.f90 b/src/modules/mesh/moduleMesh@boundaryParticle.f90
index 2e2a743..963116b 100644
--- a/src/modules/mesh/moduleMesh@boundaryParticle.f90
+++ b/src/modules/mesh/moduleMesh@boundaryParticle.f90
@@ -487,6 +487,36 @@ submodule(moduleMesh) boundaryParticle
 
     end subroutine quasiNeutrality_print
 
+    ! outflowAdaptive
+    ! Adjust the reflection coefficient of the boundary to maintain a quasi-neutral outflow
+    module subroutine outflowAdaptive(self, edge, part)
+      use moduleSpecies
+      use moduleRefParam, only: v_ref
+      implicit none
+
+      class(boundaryOutflowAdaptive), intent(inout):: self
+      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
+
+      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 genericReflection(edge, part)
+        ! print *,part%v(1)
+        ! print *
+
+      else
+        call genericTransparent(edge, part)
+
+      end if
+
+    end subroutine outflowAdaptive
+
     ! Generic boundary conditions for internal use
     ! reflection
     module subroutine genericReflection(edge, part)
diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90
index 26ef921..d22258a 100644
--- a/src/modules/moduleInject.f90
+++ b/src/modules/moduleInject.f90
@@ -281,7 +281,7 @@ MODULE moduleInject
                                                      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
+                (dot_product(partInj(n)%v, direction) <= 0.D0)) then
               partInj(n)%v = - partInj(n)%v
 
             end if