I need to ensure quasi-neutrality at the inlet with a new kind of BC

src/modules/mesh/moduleMeshBoundary.f90

@@ -77,10 +77,23 @@ 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
-      USE moduleBoundary
       USE moduleRandom
       IMPLICIT NONE
 
@@ -204,19 +217,28 @@ MODULE moduleMeshBoundary
 
     END SUBROUTINE ionization
 
-    !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
+    subroutine quasiNeutrality(edge, part)
+      use moduleRandom
+      implicit none
 
-      CLASS(meshEdge), INTENT(inout):: edge
-      CLASS(particle), INTENT(inout):: part
+      class(meshEdge), intent(inout):: edge
+      class(particle), intent(inout):: part
 
-      CALL reflection(edge, part)
+      select type(bound => edge%boundary%bTypes(part%species%n)%obj)
+      type is(boundaryQuasiNeutrality)
+        bound%alpha = 1.0d-1
 
-    END SUBROUTINE symmetryAxis
+        if (random() < bound%alpha) then
+          call reflection(edge, part)
+
+        else
+          call transparent(edge, part)
+
+        end if
+
+      end select
+
+    end subroutine quasiNeutrality
 
     !Points the boundary function to specific type
     SUBROUTINE pointBoundaryFunction(edge, s)
@@ -236,17 +258,20 @@ 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(boundaryQuasiNeutrality)
+        edge%fBoundary(s)%apply => quasiNeutrality
 
       CLASS DEFAULT
-        CALL criticalError("Boundary type not defined in this geometry", 'pointBoundaryFunction')
+        CALL criticalError("Boundary type not defined", 'pointBoundaryFunction')
 
       END SELECT