Organized injection

src/modules/init/moduleInput.f90

@@ -1431,7 +1431,6 @@ MODULE moduleInput
 
       CALL config%info('inject', found, n_children = nInject)
       ALLOCATE(injects(1:nInject))
-      nPartInj = 0
       DO i = 1, nInject
         WRITE(iString, '(i2)') i
         object = 'inject(' // trim(iString) // ')'
@@ -1442,12 +1441,6 @@ MODULE moduleInput
         
       END DO
 
-      !Allocate array for injected particles
-      IF (nPartInj > 0) THEN
-        ALLOCATE(partInj(1:nPartInj))
-
-      END IF
-      
     END SUBROUTINE readInject
 
     SUBROUTINE readAverage(config)

src/modules/makefile

@@ -14,7 +14,7 @@ output.o: moduleSpecies.o common.o
 mesh.o: moduleList.o moduleSpecies.o moduleCollisions.o moduleCoulomb.o output.o common.o 
 	$(MAKE) -C mesh all
 
-solver.o: moduleSpecies.o moduleProbe.o common.o output.o mesh.o
+solver.o: moduleInject.o moduleSpecies.o moduleProbe.o common.o output.o mesh.o
 	$(MAKE) -C solver all
 
 init.o: common.o solver.o moduleInject.o

src/modules/moduleInject.f90

@@ -313,23 +313,21 @@ MODULE moduleInject
     end subroutine updateQuasiNeutral
 
     !Add particles for the injection
-    SUBROUTINE addParticles(self, partArray)
+    SUBROUTINE addParticles(self)
       USE moduleSpecies
       USE moduleMesh
       USE moduleRandom
       USE moduleErrors
+      use moduleList, only: partInj
       IMPLICIT NONE
 
       CLASS(injectGeneric), INTENT(in):: self
-      type(particle), allocatable, intent(inout):: partArray(:)
+      type(particle), pointer:: part
       INTEGER:: i, e
-      INTEGER:: n, sp
+      INTEGER:: sp
       CLASS(meshEdge), POINTER:: edge
-      REAL(8):: direction(1:3)
 
-      !$omp single
+      REAL(8):: direction(1:3)
-      n = 0
-      !$omp end single
 
       !Insert particles
       !$OMP DO
@@ -338,30 +336,30 @@ MODULE moduleInject
         edge => self%edges(e)%obj
         ! Inject particles in edge
         DO i = 1, self%particlesPerEdge(e)
+          allocate(part)
           ! Index in the partInj array
-          n = n + 1
+          !Particle is considered to be inside the domain
-          !Particle is considered to be outside the domain for now
+          part%n_in = .true.
-          partArray(n)%n_in = .FALSE.
           !Random position in edge
-          partArray(n)%r = edge%randPos()
+          part%r = edge%randPos()
           !Assign weight to particle.
-          partArray(n)%weight = self%weightPerEdge(e)
+          part%weight = self%weightPerEdge(e)
           !Volume associated to the edge:
           IF (ASSOCIATED(edge%e1)) THEN
-            partArray(n)%cell = edge%e1%n
+            part%cell = edge%e1%n
 
           ELSEIF (ASSOCIATED(edge%e2)) THEN
-            partArray(n)%cell = edge%e2%n
+            part%cell = edge%e2%n
 
           ELSE
             CALL criticalError("No Volume associated to edge", 'addParticles')
 
           END IF
-          partArray(n)%cellColl = edge%eColl%n
+          part%cellColl = edge%eColl%n
           sp = self%species%n
 
           !Assign particle type
-          partArray(n)%species => self%species
+          part%species => self%species
 
           if (all(self%n == 0.D0)) then
             direction =  edge%normal
@@ -371,23 +369,28 @@ MODULE moduleInject
 
           end if
 
-          partArray(n)%v = 0.D0
+          part%v = 0.D0
 
           ! do while(dot_product(partInj(n)%v, direction) <= 0.d0)
-            partArray(n)%v = self%vMod*direction +  (/ self%v(1)%obj%randomVel(), &
+            part%v = self%vMod*direction +  (/ self%v(1)%obj%randomVel(), &
                                                      self%v(2)%obj%randomVel(), &
                                                      self%v(3)%obj%randomVel() /)
             !If injecting a no-drift distribution and velocity is negative, reflect
             if ((self%vMod == 0.D0) .and. &
-                (dot_product(partArray(n)%v, direction) <= 0.D0)) then
+                (dot_product(part%v, direction) <= 0.D0)) then
-              partArray(n)%v = - partArray(n)%v
+              part%v = - part%v
 
             end if
 
           ! end do
 
           !Obtain natural coordinates of particle in cell
-          partArray(n)%Xi = mesh%cells(partArray(n)%cell)%obj%phy2log(partArray(n)%r)
+          part%Xi = mesh%cells(part%cell)%obj%phy2log(part%r)
+
+          ! Add particle to global list
+          call partInj%setLock()
+          call partInj%add(part)
+          call partInj%unsetLock()
 
         END DO
 

src/modules/moduleList.f90

@@ -23,6 +23,7 @@ MODULE moduleList
 
   END TYPE listNode
 
+  TYPE(listNode):: partInj !Particles comming from injections
   TYPE(listNode):: partWScheme !Particles comming from the nonAnalogue scheme
   TYPE(listNode):: partCollisions !Particles created in collisional process
   TYPE(listNode):: partSurfaces !Particles created in surface interactions

src/modules/moduleSpecies.f90

@@ -48,11 +48,8 @@ MODULE moduleSpecies
 
   !Number of old particles
   INTEGER:: nPartOld
-  !Number of injected particles
-  INTEGER:: nPartInj
   !Arrays that contain the particles
   TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: partOld !array of particles from previous iteration
-  TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: partInj !array of inject particles
 
   CONTAINS
     FUNCTION speciesName2Index(speciesName) RESULT(sp)

src/modules/solver/moduleSolver.f90

@@ -235,23 +235,21 @@ MODULE moduleSolver
 
       INTEGER:: nn, n, e
       INTEGER, SAVE:: nPartNew
-      INTEGER, SAVE:: nInjIn, nOldIn, nWScheme, nCollisions, nSurfaces
+      INTEGER, SAVE:: nOldIn, nInj, nWScheme, nCollisions, nSurfaces
       TYPE(particle), ALLOCATABLE, SAVE:: partTemp(:)
       INTEGER:: s
 
       !$OMP SECTIONS
       !$OMP SECTION
-      nInjIn = 0
-      IF (ALLOCATED(partInj)) THEN
-        nInjIn = COUNT(partInj%n_in)
-
-      END IF
-      !$OMP SECTION
       nOldIn = 0
       IF (ALLOCATED(partOld)) THEN
         nOldIn = COUNT(partOld%n_in)
 
       END IF
+
+      !$OMP SECTION
+      nInj = partInj%amount
+
       !$OMP SECTION
       nWScheme = partWScheme%amount
 
@@ -267,30 +265,14 @@ MODULE moduleSolver
 
       !$OMP SINGLE
       CALL MOVE_ALLOC(partOld, partTemp)
-      nPartNew = nInjIn + nOldIn + nWScheme + nCollisions + nSurfaces
+      nPartNew = nInj + nOldIn + nWScheme + nCollisions + nSurfaces
       ALLOCATE(partOld(1:nPartNew))
       !$OMP END SINGLE
 
       !$OMP SECTIONS
-      !$OMP SECTION
-      !Reset particles from injection
-      nn = 0
-      DO n = 1, nPartInj
-        IF (partInj(n)%n_in) THEN
-          nn = nn + 1
-          partOld(nn) = partInj(n)
-          IF (nProbes > 0) THEN
-            CALL doProbes(partOld(nn))
-
-          END IF
-
-        END IF
-
-      END DO
-
       !$OMP SECTION
       !Reset particles from previous iteration
-      nn = nInjIn
+      nn = 0
       DO n = 1, nPartOld
         IF (partTemp(n)%n_in) THEN
           nn = nn + 1
@@ -304,19 +286,24 @@ MODULE moduleSolver
 
       END DO
 
+      !$OMP SECTION
+      !Reset particles from injection
+      nn = nOldIn
+      call resetList(partInj, partOld, nn)
+
       !$OMP SECTION
       !Reset particles from weighting scheme
-      nn = nInjIn + nOldIn
+      nn = nOldIn + nInj
       CALL resetList(partWScheme, partOld, nn)
 
       !$OMP SECTION
       !Reset particles from collisional process
-      nn = nInjIn + nOldIn + nWScheme
+      nn = nOldIn + nInj + nWScheme
       CALL resetList(partCollisions, partOld, nn)
 
       !$OMP SECTION
       !Reset particles from surface process
-      nn = nInjIn + nOldIn + nWScheme + nCollisions
+      nn = nOldIn + nInj + nWScheme + nCollisions
       CALL resetList(partSurfaces, partOld, nn)
 
       !$OMP SECTION
@@ -473,6 +460,28 @@ MODULE moduleSolver
 
     END SUBROUTINE splitParticle
 
+    !Injection of particles
+    SUBROUTINE doInjects()
+      USE moduleSpecies
+      use moduleInject
+      IMPLICIT NONE
+
+      INTEGER:: i, n, sp
+      
+      DO i=1, nInject
+        associate(inject => injects(i)%obj)
+        sp = inject%species%n
+        IF (solver%pusher(sp)%pushSpecies) THEN 
+
+          CALL inject%addParticles()
+
+        END IF
+        end associate
+
+      END DO
+
+    END SUBROUTINE doInjects
+
     SUBROUTINE updateParticleCell(self, part)
       USE moduleSpecies
       USE moduleMesh