Had to remove charge, it is working now and it only accounts for charged species

src/modules/mesh/moduleMesh@boundaryParticle.f90

@@ -546,10 +546,11 @@ submodule(moduleMesh) boundaryParticle
       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):: mom_center
+      real(8):: den_center, mom_center, vel_center
-      real(8):: mom_incident, mom_rest
+      real(8):: vel_incident, vel_rest
       real(8):: alpha
       
       select type(self)
@@ -557,34 +558,52 @@ submodule(moduleMesh) boundaryParticle
         do e = 1, size(self%edges)
           edge => self%edges(e)%obj
 
-          mom_incident = 0.d0
+          vel_incident = 0.d0
-          mom_rest     = 0.d0
+          vel_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
+            den_center = 0.0d0
-              node => mesh%nodes(nodes(n))%obj
+            mom_center = 0.0d0
+            vel_center = 0.0d0
+            select type(sp => species(s)%obj)
+            type is(speciesCharged)
+              do n = 1, edge%nNodes
+                node => mesh%nodes(nodes(n))%obj
 
-              mom_nodes(n) = dot_product(node%output(s)%mom, edge%normal)
+                den_nodes(n) = node%output(s)%den
+                mom_nodes(n) = dot_product(node%output(s)%mom, edge%normal)
 
-            end do
+              end do
 
-            mom_center = qSpecies(s)*edge%gatherF(edge%centerXi(), edge%nNodes, mom_nodes)
+              den_center = edge%gatherF(edge%centerXi(), edge%nNodes, den_nodes)
+              if (den_center > 1.0d-10) then
+                mom_center = edge%gatherF(edge%centerXi(), edge%nNodes, mom_nodes)
+                vel_center = mom_center / den_center
 
-            if (s == self%s_incident) then
+              else
-              mom_incident = mom_center
+                vel_center = 0.0d0
-            else
+              
-              mom_rest = mom_rest + mom_center
+              end if
-            end if
+
+              if (s == self%s_incident) then
+                vel_incident = vel_center
+              else
+                vel_rest = vel_rest + vel_center
+              end if
+
+            end select
 
           end do
 
-          alpha = mom_rest + mom_incident
+          alpha = vel_rest - vel_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