Finally, some progress

I rewrote how particles are injected. Now the particles per edge and its
weight are calculated in the initialization. There is the possibility
for the user to select the particles per edge.

TODO: Write documentation for new feature.
TODO: Test in 2DCyl

src/modules/init/moduleInput.f90

@@ -1236,6 +1236,7 @@ MODULE moduleInput
       REAL(8):: flow
       CHARACTER(:), ALLOCATABLE:: units
       INTEGER:: physicalSurface
+      INTEGER:: particlesPerEdge
       INTEGER:: sp
 
       CALL config%info('inject', found, n_children = nInject)
@@ -1260,8 +1261,10 @@ MODULE moduleInput
         CALL config%get(object // '.flow',            flow,            found)
         CALL config%get(object // '.units',           units,           found)
         CALL config%get(object // '.physicalSurface', physicalSurface, found)
+        particlesPerEdge = 0
+        CALL config%get(object // '.particlesPerEdge', particlesPerEdge, found)
 
-        CALL inject(i)%init(i, v, normal, T, flow, units, sp, physicalSurface)
+        CALL inject(i)%init(i, v, normal, T, flow, units, sp, physicalSurface, particlesPerEdge)
 
         CALL readVelDistr(config, inject(i), object)
         

src/modules/moduleInject.f90

@@ -61,6 +61,8 @@ MODULE moduleInject
     CLASS(speciesGeneric), POINTER:: species !Species of injection
     INTEGER:: nEdges
     INTEGER, ALLOCATABLE:: edges(:) !Array with edges
+    INTEGER, ALLOCATABLE:: particlesPerEdge(:) ! Particles per edge
+    REAL(8), ALLOCATABLE:: weightPerEdge(:) ! Weight per edge
     REAL(8):: surface ! Total surface of injection
     TYPE(velDistCont):: v(1:3) !Velocity distribution function in each direction
     CONTAINS
@@ -74,7 +76,7 @@ MODULE moduleInject
 
   CONTAINS
     !Initialize an injection of particles
-    SUBROUTINE initInject(self, i, v, n, T, flow, units, sp, physicalSurface)
+    SUBROUTINE initInject(self, i, v, n, T, flow, units, sp, physicalSurface, particlesPerEdge)
       USE moduleMesh
       USE moduleRefParam
       USE moduleConstParam
@@ -86,52 +88,28 @@ MODULE moduleInject
       CLASS(injectGeneric), INTENT(inout):: self
       INTEGER, INTENT(in):: i
       REAL(8), INTENT(in):: v, n(1:3), T(1:3)
-      INTEGER, INTENT(in):: sp, physicalSurface
+      INTEGER, INTENT(in):: sp, physicalSurface, particlesPerEdge
       REAL(8):: tauInject
       REAL(8), INTENT(in):: flow
       CHARACTER(:), ALLOCATABLE, INTENT(in):: units
       INTEGER:: e, et
       INTEGER:: phSurface(1:mesh%numEdges)
       INTEGER:: nVolColl
+      REAL(8):: fluxPerStep = 0.D0
 
       self%id      =  i
       self%vMod    =  v / v_ref
       self%n       =  n / NORM2(n)
       self%T       =  T / T_ref
-      self%species    => species(sp)%obj
-      tauInject = tau(self%species%n)
-      SELECT CASE(units)
-      CASE ("sccm")
-        !Standard cubic centimeter per minute
-        self%nParticles = INT(flow*sccm2atomPerS*tauInject*ti_ref/species(sp)%obj%weight)
-
-      CASE ("A")
-        !Input current in Ampers
-        self%nParticles = INT(flow*tauInject*ti_ref/(qe*species(sp)%obj%weight))
-
-      CASE ("Am2")
-        !Input current in Ampers per square meter
-        self%nParticles = INT(flow*tauInject*ti_ref*L_ref**2/(qe*species(sp)%obj%weight))
-
-      CASE ("part/s")
-        !Input current in Ampers
-        self%nParticles = INT(flow*tauInject*ti_ref/species(sp)%obj%weight)
-
-      CASE DEFAULT
-        CALL criticalError("No support for units: " // units, 'initInject')
-        
-      END SELECT
-      !Scale particles for different species steps
-      IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject')
-
       !Gets the edge elements from which particles are injected
       DO e = 1, mesh%numEdges
         phSurface(e) = mesh%edges(e)%obj%physicalSurface
 
       END DO
-
       self%nEdges = COUNT(phSurface == physicalSurface)
-      ALLOCATE(inject(i)%edges(1:self%nEdges))
+      ALLOCATE(self%edges(1:self%nEdges))
+      ALLOCATE(self%particlesPerEdge(1:self%nEdges))
+      ALLOCATE(self%weightPerEdge(1:self%nEdges))
       et = 0
       DO e=1, mesh%numEdges
         IF (mesh%edges(e)%obj%physicalSurface == physicalSurface) THEN
@@ -170,6 +148,60 @@ MODULE moduleInject
 
       END DO
 
+      ! Information about species and flux
+      self%species => species(sp)%obj
+      tauInject    = tau(self%species%n)
+      SELECT CASE(units)
+      CASE ("sccm")
+        !Standard cubic centimeter per minute
+        fluxPerStep = flow*sccm2atomPerS*tauInject*ti_ref
+
+      CASE ("A")
+        !Current in Ampers
+        fluxPerStep = flow*tauInject*ti_ref/qe
+
+      CASE ("Am2")
+        !Input current in Ampers per square meter
+        fluxPerStep = flow*tauInject*ti_ref*self%surface*L_ref**2/qe
+
+      CASE ("part/s")
+        !Input current in Ampers
+        fluxPerStep = flow*tauInject*ti_ref
+
+      CASE DEFAULT
+        CALL criticalError("No support for units: " // units, 'initInject')
+
+      END SELECT
+      fluxPerStep = fluxPerStep / self%surface
+
+      !Assign particles per edge
+      IF (particlesPerEdge > 0) THEN
+        ! Particles per edge defined by the user
+        self%particlesPerEdge = particlesPerEdge
+        DO et = 1, self%nEdges
+          self%weightPerEdge(et) = fluxPerStep*mesh%edges(self%edges(et))%obj%surface / REAL(particlesPerEdge)
+
+        END DO
+
+      ELSE
+        ! No particles assigned per edge, use the species weight
+        self%weightPerEdge = self%species%weight
+        DO et = 1, self%nEdges
+          self%particlesPerEdge(et) = FLOOR(fluxPerStep*mesh%edges(self%edges(et))%obj%surface /self%species%weight)
+
+        END DO
+
+      END IF
+
+      print *, self%particlesPerEdge
+      print *, self%weightPerEdge
+
+      self%nParticles = SUM(self%particlesPerEdge)
+      print *, self%nParticles
+
+      !Scale particles for different species steps
+      IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject')
+
     END SUBROUTINE initInject
 
     !Injection of particles
@@ -285,9 +317,10 @@ MODULE moduleInject
       CLASS(injectGeneric), INTENT(in):: self
       INTEGER:: randomX
       INTEGER, SAVE:: nMin, nMax !Min and Max index in partInj array
-      INTEGER:: i
+      INTEGER:: i, e
       INTEGER:: n, sp
       CLASS(meshEdge), POINTER:: randomEdge
+      INTEGER:: particlesPerEdge
       REAL(8):: direction(1:3)
 
       !Insert particles
@@ -300,21 +333,23 @@ MODULE moduleInject
         END IF
 
       END DO
-
       nMin = nMin + 1
-      nMax = nMin + self%nParticles - 1
-      !Particle is considered to be outside the domain
-      partInj(nMin:nMax)%n_in = .FALSE.
       !$OMP END SINGLE
 
       !$OMP DO
-      DO n = nMin, nMax
+      DO e = 1, self%nEdges
-        randomX    =  random(1, self%nEdges)
+        ! Select edge for injection
-        randomEdge => mesh%edges(self%edges(randomX))%obj
+        randomEdge => mesh%edges(self%edges(e))%obj
+        ! Inject particles in edge
+        DO i = 1, self%particlesPerEdge(e)
+          ! Index in the global partInj array
+          n = nMin - 1 + SUM(self%particlesPerEdge(1:e-1)) + i
+          !Particle is considered to be outside the domain
+          partInj(n)%n_in = .FALSE.
           !Random position in edge
           partInj(n)%r = randomEdge%randPos()
           !Assign weight to particle.
-        partInj(n)%weight = self%species%weight * randomEdge%surface / self%surface
+          partInj(n)%weight = self%weightPerEdge(e)
           !Volume associated to the edge:
           IF (ASSOCIATED(randomEdge%e1)) THEN
             partInj(n)%cell = randomEdge%e1%n
@@ -354,6 +389,8 @@ MODULE moduleInject
           CALL solver%updateParticleCell(partInj(n))
 
         END DO
+
+      END DO
       !$OMP END DO
 
     END SUBROUTINE addParticles