First verstion that seems to work!

src/modules/init/moduleInput.f90

@@ -898,10 +898,10 @@ MODULE moduleInput
 
               s_incident = speciesName2Index(speciesName)
 
-              call initQuasiNeutrality(bound, s_incident)
+              call quasiNeutrality_init(bound, s_incident)
 
             CASE('outflowAdaptive')
-              ALLOCATE(boundaryOutflowAdaptive:: bound)
+              call outflowAdaptive_init(bound)
 
             CASE DEFAULT
               CALL criticalError('Boundary type ' // bType // ' undefined', 'readBoundary')
@@ -1295,6 +1295,7 @@ MODULE moduleInput
               ! If the boundary for the species is linked to the one analysing, add the edges
               if (associated(physicalSurfaces(ps)%particles(s)%obj, bound)) then
                 bound%edges = [bound%edges, physicalSurfaces(ps)%edges]
+
               end if
 
             end do
@@ -1304,6 +1305,24 @@ MODULE moduleInput
           allocate(bound%alpha(1:mesh%numEdges)) ! TODO: Change this so only the edges associated to the boundary are here
           bound%alpha = 0.d0
 
+        type is(boundaryOutflowAdaptive)
+          ! Loop over all physical surfaces
+          do ps = 1, nPhysicalSurfaces
+            ! Loop over all species
+            do s = 1, nSpecies
+              ! If the boundary for the species is linked to the one analysing, add the edges
+              if (associated(physicalSurfaces(ps)%particles(s)%obj, bound)) then
+                bound%edges = [bound%edges, physicalSurfaces(ps)%edges]
+
+              end if
+
+            end do
+
+          end do
+
+          allocate(bound%velocity_shift(1:mesh%numEdges)) ! TODO: Change this so only the edges associated to the boundary are here
+          bound%velocity_shift = 0.d0
+
         end select
 
       end do

src/modules/mesh/moduleMesh.f90

@@ -682,11 +682,17 @@ MODULE moduleMesh
 
     end subroutine initIonization
 
-    module subroutine initQuasiNeutrality(boundary, s_incident)
+    module subroutine quasiNeutrality_init(boundary, s_incident)
       class(boundaryParticleGeneric), allocatable, intent(inout):: boundary
       integer, intent(in):: s_incident
 
-    end subroutine initQuasiNeutrality
+    end subroutine quasiNeutrality_init
+
+    module subroutine outflowAdaptive_init(boundary)!, s_incident)
+      class(boundaryParticleGeneric), allocatable, intent(inout):: boundary
+      ! integer, intent(in):: s_incident
+
+    end subroutine outflowAdaptive_init
 
     module function boundaryParticleName_to_Index(boundaryName) result(bp)
       character(:), allocatable:: boundaryName
@@ -791,8 +797,9 @@ MODULE moduleMesh
 
   !Boundary for quasi-neutral outflow adjusting reflection coefficient
   type, public, extends(boundaryParticleGeneric):: boundaryOutflowAdaptive
-    real(8):: outflowCurrent
-    real(8):: reflectionFraction
+    real(8), allocatable:: velocity_shift(:)
+    integer:: s_incident ! species index of the incident species
+    type(meshEdgePointer), allocatable:: edges(:) !Array with edges
     contains
       procedure, pass:: apply => outflowAdaptive
 

src/modules/mesh/moduleMesh@boundaryParticle.f90

@@ -360,7 +360,7 @@ submodule(moduleMesh) boundaryParticle
     ! quasiNeutrality
     ! Maintains n_incident = n_rest by modifying the reflection parameter of the edge.
     ! Init
-    module subroutine initQuasiNeutrality(boundary, s_incident)
+    module subroutine quasiNeutrality_init(boundary, s_incident)
       implicit none
 
       class(boundaryParticleGeneric), allocatable, intent(inout):: boundary
@@ -379,7 +379,7 @@ submodule(moduleMesh) boundaryParticle
 
       end select
 
-    end subroutine initQuasiNeutrality
+    end subroutine quasiNeutrality_init
 
     ! Apply
     module subroutine quasiNeutrality(self, edge, part)
@@ -416,7 +416,7 @@ submodule(moduleMesh) boundaryParticle
       select type(self)
       type is(boundaryQuasiNeutrality)
         do e = 1, size(self%edges)
-          edge => mesh%edges(e)%obj
+          edge => self%edges(e)%obj
 
           density_incident = 0.d0
           density_rest     = 0.d0
@@ -425,7 +425,8 @@ submodule(moduleMesh) boundaryParticle
           allocate(density_nodes(1:edge%nNodes))
           do s = 1, nSpecies
             do n = 1, edge%nNodes
-              node => mesh%nodes(n)%obj
+              node => mesh%nodes(nodes(n))%obj
+
               density_nodes(n) = node%output(s)%den
 
             end do
@@ -489,6 +490,28 @@ submodule(moduleMesh) boundaryParticle
 
     ! outflowAdaptive
     ! Adjust the reflection coefficient of the boundary to maintain a quasi-neutral outflow
+    ! Init
+    subroutine outflowAdaptive_init(boundary)
+      implicit none
+
+      class(boundaryParticleGeneric), allocatable, intent(inout):: boundary
+
+      allocate(boundaryOutflowAdaptive:: boundary)
+      
+      select type(boundary)
+      type is(boundaryOutflowAdaptive)
+        allocate(boundary%edges(0))
+
+        boundary%s_incident = 1 !TODO Make this an input parameter
+
+        boundary%update => outflowAdaptive_update
+        boundary%print  => outflowAdaptive_print
+
+      end select
+
+    end subroutine outflowAdaptive_init
+
+    ! Apply
     module subroutine outflowAdaptive(self, edge, part)
       use moduleSpecies
       use moduleRefParam, only: v_ref
@@ -497,13 +520,10 @@ submodule(moduleMesh) boundaryParticle
       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
 
       call genericReflection(edge, part)
 
-      part%v(1) = part%v(1) + v_cut
+      part%v(1) = part%v(1) - self%velocity_shift(edge%n)
 
       if (dot_product(part%v,edge%normal) <= 0.d0) then
 
@@ -513,6 +533,94 @@ submodule(moduleMesh) boundaryParticle
 
     end subroutine outflowAdaptive
 
+    ! Update
+    subroutine outflowAdaptive_update(self)
+      implicit none
+
+      class(boundaryParticleGeneric), intent(inout):: self
+      integer:: e, n, s
+      integer, allocatable:: nodes(:)
+      class(meshEdge), pointer:: edge
+      real(8), allocatable:: den_nodes(:)
+      real(8), allocatable:: mom_nodes(:)
+      class(meshNode), pointer:: node
+      real(8):: den_center, mom_center, vel_center
+      real(8):: velocity_incident, velocity_rest
+      real(8):: alpha
+      
+      select type(self)
+      type is(boundaryOutflowAdaptive)
+        do e = 1, size(self%edges)
+          edge => self%edges(e)%obj
+
+          velocity_incident = 0.d0
+          velocity_rest     = 0.d0
+
+          nodes = edge%getNodes(edge%nNodes)
+          allocate(den_nodes(1:edge%nNodes))
+          allocate(mom_nodes(1:edge%nNodes))
+          do s = 1, nSpecies
+            do n = 1, edge%nNodes
+              node => mesh%nodes(nodes(n))%obj
+
+              den_nodes(n) = node%output(s)%den
+              mom_nodes(n) = dot_product(node%output(s)%mom, edge%normal)
+
+            end do
+
+            den_center = edge%gatherF(edge%centerXi(), edge%nNodes, den_nodes)
+            mom_center = edge%gatherF(edge%centerXi(), edge%nNodes, mom_nodes)
+
+            if (den_center > 1.0d-10) then
+              vel_center = mom_center / den_center
+
+            else
+              vel_center = 0.0d0
+            end if
+
+            if (s == self%s_incident) then
+              velocity_incident = vel_center
+            else
+              velocity_rest = velocity_rest + vel_center
+            end if
+
+          end do
+
+          alpha = velocity_rest - velocity_incident
+
+          ! Apply correction with a factor of 0.1 to avoid fast changes
+          self%velocity_shift(edge%n) = self%velocity_shift(edge%n) + 1.0d-2 * alpha
+
+          deallocate(den_nodes)
+          deallocate(mom_nodes)
+
+        end do
+
+      end select
+
+    end subroutine outflowAdaptive_update
+
+    ! Print
+    subroutine outflowAdaptive_print(self, fileID)
+      implicit none
+
+      class(boundaryParticleGeneric), intent(inout):: self
+      integer, intent(in):: fileID
+      integer:: e
+
+      write(fileID, '(A)') self%name
+      select type(self)
+      type is(boundaryOutflowAdaptive)
+        write(fileID, '(A,",",A)') '"Edge id"', '"vel shift non-dimensional"'
+        do e = 1, size(self%edges)
+          write(fileID, '('//fmtColInt//','//fmtColReal//')') self%edges(e)%obj%n, self%velocity_shift(self%edges(e)%obj%n)
+        end do
+
+      end select
+
+    end subroutine outflowAdaptive_print
+
+
     ! Generic boundary conditions for internal use
     ! reflection
     module subroutine genericReflection(edge, part)