txs auto checkin

doc/user-manual/fpakc_UserManual.tex

@@ -697,8 +697,7 @@ make
                                \begin{itemize}
                                \item \textbf{constant}: The flow does not change.
                                  This is the default value if none is provided.
-                               \item \textbf{quasiNeutral}: The flow of the species is changed over time to keep the charge density close to zero at the injection surface.
-                                 The weight of the injected particles is changed.
+                               \item \textbf{quasiNeutral}: The weight of the injected species' particles is changed over time to keep a zero electric field at the boundary.
                                \end{itemize}
           \item \textbf{physicalSurface}: Integer.
                                           Identification of the edge in the mesh file.

src/modules/moduleInject.f90

@@ -53,6 +53,7 @@ MODULE moduleInject
   end type injectconstant
 
   type, extends(injectGeneric):: injectQuasiNeutral
+    real(8):: qSign_incident
 
   end type injectQuasiNeutral
 
@@ -228,6 +229,14 @@ MODULE moduleInject
       !Scale particles for different species steps
       IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject')
 
+      ! Specific parameters per type
+      select type(self)
+      type is(injectQuasiNeutral)
+        ! Get the sign of the charged species
+        self%qSign_incident = sign(1.d0, self%species%qm)
+
+      end select
+
     END SUBROUTINE initInject
 
     ! Update the value of injects
@@ -248,17 +257,15 @@ MODULE moduleInject
 
     ! Updates the flow in the injection to maintain quasineutrality
     subroutine updateQuasiNeutral(self)
-      use moduleMesh, only: meshEdge, meshNode, mesh, qSpecies
+      use moduleMesh, only: meshEdge, meshCell
       implicit none
 
       class(injectGeneric), intent(inout):: self
-      integer:: e, s, n
-      integer, allocatable:: nodes(:)
+      integer:: e
       class(meshEdge), pointer:: edge
-      real(8), allocatable:: density_nodes(:)
-      class(meshNode), pointer:: node
-      real(8):: den_center
-      real(8):: density_incident, density_rest
+      class(meshCell), pointer:: cell
+      real(8):: Xi(1:3)
+      real(8):: EF_normal
       real(8):: alpha
 
       select type(self)
@@ -266,42 +273,21 @@ MODULE moduleInject
         do e = 1, self%nEdges
           edge => self%edges(e)%obj
 
-          density_incident = 0.d0
-          density_rest     = 0.d0
+          Xi = edge%centerXi()
 
-          nodes = edge%getNodes(edge%nNodes)
-          allocate(density_nodes(1:edge%nNodes))
-          do s = 1, nSpecies
-            do n = 1, edge%nNodes
-              node => mesh%nodes(nodes(n))%obj
-
-              density_nodes(n) = node%output(s)%den
-
-            end do
-
-            ! Gather the density at the edge center and multiply by the species charge
-            den_center = qSpecies(s)*edge%gatherF(edge%centerXi(), edge%nNodes, density_nodes)
-
-            if (s == self%species%n) then
-              density_incident = den_center
+          if (associated(edge%e1)) then
+            cell => edge%e1
 
           else
-              density_rest = density_rest + den_center
+            cell => edge%e2
 
           end if
 
-          end do
+          ! Projection of EF on the edge normal vector
+          EF_normal = dot_product(cell%gatherElectricField(Xi), edge%normal)
 
           ! Correction for this time step
-          if (density_rest > 1.0d-10) then
-            ! If there is a rest population, correct
-            alpha = 1.d0 + density_incident/density_rest
-
-          else
-            ! If not, alpha is assumed unchaged
-            alpha = 0.d0
-
-          end if
+          alpha = self%qSign_incident * EF_normal
 
           ! Adjust the weight of particles to match the new current
           self%weightPerEdge(e) = self%weightPerEdge(e) + 1.0d-1 * alpha

src/modules/solver/moduleSolver.f90

@@ -466,7 +466,7 @@ MODULE moduleSolver
       use moduleInject
       IMPLICIT NONE
 
-      INTEGER:: i, n, sp
+      INTEGER:: i, sp
       
       DO i=1, nInject
         associate(inject => injects(i)%obj)