Main changes:

- Injection is now performed in parallalel (an IF statement could be
  required to avoid overhead when number of injected particles is
  below a margin).
- Added the possibility for multiple injections.

src/DSMC_Neutrals.f95

@@ -37,11 +37,10 @@ PROGRAM DSMC_Neutrals
     !$OMP SINGLE
     tStep = omp_get_wtime()
     tPush = omp_get_wtime()
+    !$OMP END SINGLE
 
     CALL doInjects()
 
-    !$OMP END SINGLE NOWAIT
-
     !Push old particles
     CALL doPushes()
 

src/modules/moduleInject.f95

@@ -119,23 +119,19 @@ MODULE moduleInject
       CLASS(injectGeneric), INTENT(in):: self
       REAL(8):: randomX
       INTEGER:: j
-      INTEGER:: randomEdge
       REAL(8):: randomPos
-      !Edge nodes
-      INTEGER:: n1 = 0, n2 = 0
-      !Edge nodes coordinates
-      REAL(8):: p1(1:3) = 0.D0, p2(1:3) = 0.D0
       INTEGER:: nMin, nMax !Min and Max index in partInj array
       INTEGER:: n
+      CLASS(meshEdge), POINTER:: randomEdge
 
       !Insert particles
-      !TODO: Adjust for multiple injections
+      nMin = SUM(inject(1:(self%id-1))%nParticles) + 1
-      nMin = 1
+      nMax = nMin + self%nParticles
-      nMax = self%nParticles
       !Assign particle type
       partInj(nMin:nMax)%pt = self%pt
       !Assign weight to particle.
       partInj(nMin:nMax)%weight = species(self%pt)%obj%weight
+      !$OMP DO
       DO n = nMin, nMax
         !Select edge randomly from which inject particle
         randomX = RAND()*self%sumWeight
@@ -145,34 +141,16 @@ MODULE moduleInject
 
         END DO
 
-        randomEdge = self%edges(j)
+        randomEdge => mesh%edges(self%edges(j))%obj
-
+        !Random position in edge
-        !Get coordinates of edge nodes
+        randomPos = RAND()
-        SELECT TYPE(edge => mesh%edges(randomEdge)%obj)
+        partInj(n)%r = randomEdge%randPos(randomPos)
-        CLASS IS (meshEdgeCyl)
-          n1 = edge%n1%n
-          n2 = edge%n2%n
-
-        END SELECT
-        SELECT TYPE(node => mesh%nodes(n1)%obj)
-        TYPE IS (meshNodeCyl)
-          p1(1) = node%z
-          p1(2) = node%r
-
-        END SELECT
-        SELECT TYPE(node => mesh%nodes(n2)%obj)
-        TYPE IS (meshNodeCyl)
-          p2(1) = node%z
-          p2(2) = node%r
-
-        END SELECT
-
         !Volume associated to the edge:
-        IF (DOT_PRODUCT(self%n, mesh%edges(randomEdge)%obj%normal) >= 0.D0) THEN
+        IF (DOT_PRODUCT(self%n, randomEdge%normal) >= 0.D0) THEN
-          partInj(n)%e_p = mesh%edges(randomEdge)%obj%e1%n
+          partInj(n)%e_p = randomEdge%e1%n
         
         ELSE
-          partInj(n)%e_p = mesh%edges(randomEdge)%obj%e2%n
+          partInj(n)%e_p = randomEdge%e2%n
 
         END IF
 
@@ -180,20 +158,13 @@ MODULE moduleInject
                           vBC(self%v(2), self%vTh(2)), &
                           vBC(self%v(3), self%vTh(3)) /)
 
-        !Random position in edge
-        !TODO: Use edge coordinates and transformations for this process
-        randomPos = RAND()
-        partInj(n)%r(1) = p1(1) + randomPos*(p2(1) - p1(1))
-        partInj(n)%r(2) = p1(2) + randomPos*(p2(2) - p1(2))
-        partInj(n)%r(3) = p1(3) + randomPos*(p2(3) - p1(3))
-
-
         !Push new particle
         CALL push(partInj(n))
         !Assign cell to new particle
         CALL mesh%vols(partInj(n)%e_p)%obj%findCell(partInj(n))
 
       END DO
+      !$OMP END DO
 
     END SUBROUTINE addParticlesMaxwellian
 

src/modules/moduleList.f95

@@ -1,6 +1,7 @@
 !Linked list of particles
 MODULE moduleList
   USE moduleSpecies
+  USE OMP_LIB
   IMPLICIT NONE
 
   TYPE lNode 
@@ -20,6 +21,10 @@ MODULE moduleList
 
   END TYPE listNode
 
+  INTEGER(KIND=OMP_LOCK_KIND):: resetLock
+
+  TYPE(listNode):: partNew
+
   TYPE pointerArray
     TYPE(particle), POINTER:: part
 
@@ -49,6 +54,7 @@ MODULE moduleList
 
     END SUBROUTINE addToList
 
+    !converts list to array
     FUNCTION convert2Array(self) RESULT(partArray)
       IMPLICIT NONE
 

src/modules/moduleMesh.f95

@@ -59,6 +59,8 @@ MODULE moduleMesh
     INTEGER:: bt = 0
     CONTAINS
       PROCEDURE(initEdge_interface), DEFERRED, PASS:: init
+      PROCEDURE(boundary_interface), DEFERRED, PASS:: fBoundary
+      PROCEDURE(randPos_interface), DEFERRED, PASS::  randPos
 
   END TYPE meshEdge
 
@@ -73,6 +75,23 @@ MODULE moduleMesh
 
     END SUBROUTINE initEdge_interface
 
+    SUBROUTINE boundary_interface(self, part)
+      USE moduleSpecies
+
+      IMPORT:: meshEdge
+      CLASS (meshEdge), INTENT(inout):: self
+      CLASS (particle), INTENT(inout):: part
+
+    END SUBROUTINE
+
+    PURE FUNCTION randPos_interface(self, rnd) RESULT(r)
+      IMPORT:: meshEdge
+      CLASS(meshEdge), INTENT(in):: self
+      REAL(8), INTENT(in):: rnd
+      REAL(8):: r(1:3)
+
+    END FUNCTION randPos_interface
+
   END INTERFACE
 
   !Containers for edges in the mesh

src/modules/moduleMeshCyl.f95

@@ -25,23 +25,11 @@ MODULE moduleMeshCyl
     !Connectivity to nodes
     CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL()
     CONTAINS
-      PROCEDURE, PASS::                               init      => initEdgeCyl
+      PROCEDURE, PASS:: init    => initEdgeCyl
-      PROCEDURE(boundary_interface), PASS, DEFERRED:: fBoundary
+      PROCEDURE, PASS:: randPos => randPosCyl
 
   END TYPE meshEdgeCyl
 
-  ABSTRACT INTERFACE
-    SUBROUTINE boundary_interface(self, part)
-      USE moduleSpecies
-
-      IMPORT:: meshEdgeCyl
-      CLASS (meshEdgeCyl), INTENT(inout):: self
-      CLASS (particle), INTENT(inout):: part
-
-    END SUBROUTINE
-
-  END INTERFACE
-
   TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVolCyl
     CONTAINS
       PROCEDURE, PASS:: collision => collision2DCyl
@@ -140,6 +128,20 @@ MODULE moduleMeshCyl
 
     END SUBROUTINE initEdgeCyl
 
+
+    PURE FUNCTION randPosCyl(self, rnd) RESULT(r)
+      CLASS(meshEdgeCyl), INTENT(in):: self
+      REAL(8), INTENT(in):: rnd
+      REAL(8):: r(1:3)
+      REAL(8):: p1(1:2), p2(1:2)
+
+      p1 = (/self%z(1), self%r(1) /)
+      p2 = (/self%z(2), self%r(2) /)
+      r(1:2) = (1.D0 - rnd)*p1 + rnd*p2
+      r(3)   =  0.D0
+
+    END FUNCTION randPosCyl
+
     !Vol functions
     !Quad Volume
     SUBROUTINE initVolQuadCyl(self, n, p)

src/modules/moduleSolver.f95

@@ -74,31 +74,14 @@ MODULE moduleSolver
 
     END SUBROUTINE doCollisions
     
-    !Scatter particles in the grid
-    SUBROUTINE doScatter
-      USE moduleSpecies
-      USE moduleMesh
-
-      INTEGER:: n
-
-      !Loops over the particles to scatter them
-      !$OMP DO
-      DO n=1, nPartOld
-        CALL mesh%vols(partOld(n)%e_p)%obj%scatter(partOld(n))
-
-      END DO
-      !$OMP END DO
-
-    END SUBROUTINE doScatter
-
     SUBROUTINE doReset()
       USE moduleSpecies
       IMPLICIT NONE
 
       INTEGER:: nn, n
-      INTEGER:: nPartNew
+      INTEGER, SAVE:: nPartNew
       INTEGER, SAVE:: nInjIn, nOldIn
-      TYPE(particle), ALLOCATABLE:: partTemp(:)
+      TYPE(particle), ALLOCATABLE, SAVE:: partTemp(:)
 
       !$OMP SECTIONS
       !$OMP SECTION
@@ -121,7 +104,10 @@ MODULE moduleSolver
       CALL MOVE_ALLOC(partOld, partTemp)
       nPartNew = nInjIn + nOldIn
       ALLOCATE(partOld(1:nPartNew))
+      !$OMP END SINGLE
 
+      !$OMP SECTIONS
+      !$OMP SECTION
       nn = 0
       DO n = 1, nPartInj
         IF (partInj(n)%n_in) THEN
@@ -132,6 +118,8 @@ MODULE moduleSolver
 
       END DO
 
+      !$OMP SECTION
+      nn = nInjIn
       DO n = 1, nPartOld
         IF (partTemp(n)%n_in) THEN
           nn = nn + 1
@@ -140,12 +128,31 @@ MODULE moduleSolver
         END IF
 
       END DO
+      !$OMP END SECTIONS
 
+      !$OMP SINGLE
       nPartOld = nPartNew
       !$OMP END SINGLE
 
     END SUBROUTINE doReset
 
+    !Scatter particles in the grid
+    SUBROUTINE doScatter
+      USE moduleSpecies
+      USE moduleMesh
+
+      INTEGER:: n
+
+      !Loops over the particles to scatter them
+      !$OMP DO
+      DO n=1, nPartOld
+        CALL mesh%vols(partOld(n)%e_p)%obj%scatter(partOld(n))
+
+      END DO
+      !$OMP END DO
+
+    END SUBROUTINE doScatter
+
     SUBROUTINE doOutput(t)
       USE moduleMesh
       USE moduleOutput