Organized injection

src/fpakc.f90

@@ -1,15 +1,39 @@
 !  FPAKC main program
 PROGRAM fpakc
-  USE moduleCompTime
-  USE moduleCaseParam
-  USE moduleInput
-  USE moduleInject
-  USE moduleSolver
-  USE moduleMesh
-  USE moduleProbe
-  USE moduleErrors
-  USE OMP_LIB
-  IMPLICIT NONE
+  use moduleCompTime,  only: tStep,    &
+                             tEMField, &
+                             tCoul,    &
+                             tColl,    &
+                             tPush,    &
+                             tReset,   &
+                             tWeight
+  use omp_lib,         only: omp_get_wtime
+  use moduleErrors,    only: criticalError, &
+                             verboseError
+  use moduleInput,     only: readCOnfig, &
+                             initOutput
+  use moduleCaseParam, only: timeStep, &
+                             tInitial, &
+                             tFinal
+  use moduleProbe,     only: resetProbes
+  use moduleSolver,    only: doScatter, &
+                             doEMField, &
+                             doOutput,  &
+                             solver,    &
+                             doPushes,  &
+                             doReset,   &
+                             doAverage, &
+                             doInjects
+  use moduleMesh,      only: boundariesEM_update,       &
+                             boundariesEM_update,       &
+                             boundariesParticle_update, &
+                             mesh,                      &
+                             meshForMCC,                &
+                             doMCCollisions,            &
+                             doCollisions,              &
+                             doCoulombScattering
+  use moduleInject,    only: updateInjects
+  implicit none
 
   ! arg1 = Input argument 1 (input file)
   CHARACTER(200):: arg1

src/modules/common/moduleRandom.f90

@@ -67,19 +67,20 @@ MODULE moduleRandom
     end function randomMaxwellian
 
     !Returns a random number in a Maxwellian distribution of mean 0 and width 1
-    FUNCTION randomHalfMaxwellian() RESULT(rnd)
+    function randomHalfMaxwellian() result(rnd)
       IMPLICIT NONE
 
-      REAL(8):: rnd
-      REAL(8):: x
+      real(8):: rnd
+      real(8):: v1
 
       rnd = 0.D0
-      x = 0.D0
-      DO WHILE (x == 0.D0)
-        CALL RANDOM_NUMBER(x)
-      END DO
+      v1 = 0.D0
+      do while (v1 == 0.D0)
+        v1 = random()
       
-      rnd = DSQRT(-DLOG(x))
+      end do
+
+      rnd = sqrt(-2.d0*log(v1))
 
     END FUNCTION randomHalfMaxwellian
 

src/modules/common/velocityDistribution.f90

@@ -53,12 +53,11 @@ module velocityDistribution
     function randomVelMaxwellian(self) result (v)
       use moduleRandom
       implicit none
-
       class(velDistMaxwellian), intent(in):: self
       real(8):: v
       v = 0.D0
 
-      v = self%vTh*randomMaxwellian()/sqrt(2.d0)
+      v = self%vTh*randomMaxwellian()
 
     end function randomVelMaxwellian
 

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

@@ -12,7 +12,7 @@ MODULE moduleInject
     REAL(8):: temperature(1:3) !Temperature
     REAL(8):: n(1:3) !Direction of injection
     LOGICAL:: fixDirection !The injection of particles has a fix direction defined by n
-    INTEGER:: nParticles !Number of particles to introduce each time step
+    INTEGER:: nParticles !Number of particles to inject each time step
     CLASS(speciesGeneric), POINTER:: species !Species of injection
     INTEGER:: nEdges
     type(meshEdgePointer), allocatable:: edges(:)
@@ -72,7 +72,6 @@ MODULE moduleInject
       USE moduleRefParam
       USE moduleConstParam
       USE moduleSpecies
-      USE moduleSolver
       USE moduleErrors
       use json_module
       IMPLICIT NONE
@@ -211,8 +210,6 @@ MODULE moduleInject
 
         END DO
 
-        self%nParticles = SUM(self%particlesPerEdge)
-
       ELSE
         ! No particles assigned per edge, use the species weight
         self%weightPerEdge = self%species%weight
@@ -220,13 +217,14 @@ MODULE moduleInject
           self%particlesPerEdge(e) = max(1,FLOOR(fluxPerStep*self%edges(e)%obj%surface / self%species%weight))
         END DO
 
-        self%nParticles = SUM(self%particlesPerEdge)
-
         !Rescale weight to match flux
         self%weightPerEdge = fluxPerStep * self%surface / (real(self%nParticles))
 
       END IF
 
+      ! Total number of particles to inject
+      self%nParticles = SUM(self%particlesPerEdge)
+
       !Scale particles for different species steps
       IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject')
 
@@ -314,126 +312,85 @@ MODULE moduleInject
 
     end subroutine updateQuasiNeutral
 
-    !Injection of particles
-    SUBROUTINE doInjects()
-      USE moduleSpecies
-      USE moduleSolver
-      IMPLICIT NONE
-
-      INTEGER:: i
-      
-      !$OMP SINGLE
-      nPartInj = 0
-      DO i = 1, nInject
-        IF (solver%pusher(injects(i)%obj%species%n)%pushSpecies) THEN 
-          nPartInj = nPartInj + injects(i)%obj%nParticles
-
-        END IF
-
-      END DO
-
-      IF (ALLOCATED(partInj)) DEALLOCATE(partInj)
-      ALLOCATE(partInj(1:nPartInj))
-      !$OMP END SINGLE
-
-      DO i=1, nInject
-        IF (solver%pusher(injects(i)%obj%species%n)%pushSpecies) THEN 
-          CALL injects(i)%obj%addParticles()
-
-        END IF
-      END DO
-
-    END SUBROUTINE doInjects
-
     !Add particles for the injection
     SUBROUTINE addParticles(self)
       USE moduleSpecies
-      USE moduleSolver
       USE moduleMesh
       USE moduleRandom
       USE moduleErrors
+      use moduleList, only: partInj
       IMPLICIT NONE
 
       CLASS(injectGeneric), INTENT(in):: self
-      INTEGER, SAVE:: nMin
+      type(particle), pointer:: part
       INTEGER:: i, e
-      INTEGER:: n, sp
-      CLASS(meshEdge), POINTER:: randomEdge
+      INTEGER:: sp
+      CLASS(meshEdge), POINTER:: edge
+
       REAL(8):: direction(1:3)
 
       !Insert particles
-      !$OMP SINGLE
-      nMin = 0
-      DO i = 1, self%id -1
-        IF (solver%pusher(injects(i)%obj%species%n)%pushSpecies) THEN 
-          nMin = nMin + injects(i)%obj%nParticles
-
-        END IF
-
-      END DO
-      nMin = nMin + 1
-      !$OMP END SINGLE
-
       !$OMP DO
       DO e = 1, self%nEdges
         ! Select edge for injection
-        randomEdge => self%edges(e)%obj
+        edge => 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.
+          allocate(part)
+          ! Index in the partInj array
+          !Particle is considered to be inside the domain
+          part%n_in = .true.
           !Random position in edge
-          partInj(n)%r = randomEdge%randPos()
+          part%r = edge%randPos()
           !Assign weight to particle.
-          partInj(n)%weight = self%weightPerEdge(e)
+          part%weight = self%weightPerEdge(e)
           !Volume associated to the edge:
-          IF (ASSOCIATED(randomEdge%e1)) THEN
-            partInj(n)%cell = randomEdge%e1%n
+          IF (ASSOCIATED(edge%e1)) THEN
+            part%cell = edge%e1%n
 
-          ELSEIF (ASSOCIATED(randomEdge%e2)) THEN
-            partInj(n)%cell = randomEdge%e2%n
+          ELSEIF (ASSOCIATED(edge%e2)) THEN
+            part%cell = edge%e2%n
 
           ELSE
             CALL criticalError("No Volume associated to edge", 'addParticles')
 
           END IF
-          partInj(n)%cellColl = randomEdge%eColl%n
+          part%cellColl = edge%eColl%n
           sp = self%species%n
 
           !Assign particle type
-          partInj(n)%species => self%species
+          part%species => self%species
 
           if (all(self%n == 0.D0)) then
-            direction =  randomEdge%normal
+            direction =  edge%normal
 
           else
             direction = self%n
 
           end if
 
-          partInj(n)%v = 0.D0
+          part%v = 0.D0
 
-          do while(dot_product(partInj(n)%v, direction) <= 0.d0)
-            partInj(n)%v = self%vMod*direction +  (/ self%v(1)%obj%randomVel(), &
+          ! do while(dot_product(partInj(n)%v, direction) <= 0.d0)
+            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(partInj(n)%v, direction) <= 0.D0)) then
-              partInj(n)%v = - partInj(n)%v
+                (dot_product(part%v, direction) <= 0.D0)) then
+              part%v = - part%v
 
             end if
 
-          end do
+          ! end do
 
           !Obtain natural coordinates of particle in cell
-          partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r)
-          !Push new particle with the minimum time step
-          CALL solver%pusher(sp)%pushParticle(partInj(n), tau(sp))
-          !Assign cell to new particle
-          CALL solver%updateParticleCell(partInj(n))
+          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
 
@@ -500,7 +457,7 @@ MODULE moduleInject
       CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
       REAL(8), INTENT(in):: temperature, m
 
-      velDist = velDistMaxwellian(vTh = DSQRT(2.d0*temperature/m))
+      velDist = velDistMaxwellian(vTh = DSQRT(temperature/m))
 
     END SUBROUTINE initVelDistMaxwellian
 
@@ -510,7 +467,7 @@ MODULE moduleInject
       CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
       REAL(8), INTENT(in):: temperature, m
 
-      velDist = velDistHalfMaxwellian(vTh = DSQRT(2.d0*temperature/m))
+      velDist = velDistHalfMaxwellian(vTh = DSQRT(temperature/m))
 
     END SUBROUTINE initVelDistHalfMaxwellian
 

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