src/fpakc.f90

@@ -1,39 +1,15 @@
 !  FPAKC main program
 PROGRAM fpakc
-  use moduleCompTime,  only: tStep,    &
+  USE moduleCompTime
-                             tEMField, &
+  USE moduleCaseParam
-                             tCoul,    &
+  USE moduleInput
-                             tColl,    &
+  USE moduleInject
-                             tPush,    &
+  USE moduleSolver
-                             tReset,   &
+  USE moduleMesh
-                             tWeight
+  USE moduleProbe
-  use omp_lib,         only: omp_get_wtime
+  USE moduleErrors
-  use moduleErrors,    only: criticalError, &
+  USE OMP_LIB
-                             verboseError
+  IMPLICIT NONE
-  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,20 +67,19 @@ 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):: rnd
-      real(8):: v1
+      REAL(8):: x
 
       rnd = 0.D0
-      v1 = 0.D0
+      x = 0.D0
-      do while (v1 == 0.D0)
+      DO WHILE (x == 0.D0)
-        v1 = random()
+        CALL RANDOM_NUMBER(x)
+      END DO
 
-      end do
+      rnd = DSQRT(-DLOG(x))
-
-      rnd = sqrt(-2.d0*log(v1))
 
     END FUNCTION randomHalfMaxwellian
 

src/modules/common/velocityDistribution.f90

@@ -53,11 +53,12 @@ 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()
+      v = self%vTh*randomMaxwellian()/sqrt(2.d0)
 
     end function randomVelMaxwellian
 

src/modules/init/moduleInput.f90

@@ -1431,6 +1431,7 @@ 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) // ')'
@@ -1441,6 +1442,12 @@ 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: moduleInject.o moduleSpecies.o moduleProbe.o common.o output.o mesh.o
+solver.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 inject each time step
+    INTEGER:: nParticles !Number of particles to introduce each time step
     CLASS(speciesGeneric), POINTER:: species !Species of injection
     INTEGER:: nEdges
     type(meshEdgePointer), allocatable:: edges(:)
@@ -72,6 +72,7 @@ MODULE moduleInject
       USE moduleRefParam
       USE moduleConstParam
       USE moduleSpecies
+      USE moduleSolver
       USE moduleErrors
       use json_module
       IMPLICIT NONE
@@ -210,6 +211,8 @@ MODULE moduleInject
 
         END DO
 
+        self%nParticles = SUM(self%particlesPerEdge)
+
       ELSE
         ! No particles assigned per edge, use the species weight
         self%weightPerEdge = self%species%weight
@@ -217,14 +220,13 @@ 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')
 
@@ -312,85 +314,126 @@ 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
-      type(particle), pointer:: part
+      INTEGER, SAVE:: nMin
       INTEGER:: i, e
-      INTEGER:: sp
+      INTEGER:: n, sp
-      CLASS(meshEdge), POINTER:: edge
+      CLASS(meshEdge), POINTER:: randomEdge
-
       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
-        edge => self%edges(e)%obj
+        randomEdge => self%edges(e)%obj
         ! Inject particles in edge
         DO i = 1, self%particlesPerEdge(e)
-          allocate(part)
+          ! Index in the global partInj array
-          ! Index in the partInj array
+          n = nMin - 1 + SUM(self%particlesPerEdge(1:e-1)) + i
-          !Particle is considered to be inside the domain
+          !Particle is considered to be outside the domain
-          part%n_in = .true.
+          partInj(n)%n_in = .FALSE.
           !Random position in edge
-          part%r = edge%randPos()
+          partInj(n)%r = randomEdge%randPos()
           !Assign weight to particle.
-          part%weight = self%weightPerEdge(e)
+          partInj(n)%weight = self%weightPerEdge(e)
           !Volume associated to the edge:
-          IF (ASSOCIATED(edge%e1)) THEN
+          IF (ASSOCIATED(randomEdge%e1)) THEN
-            part%cell = edge%e1%n
+            partInj(n)%cell = randomEdge%e1%n
 
-          ELSEIF (ASSOCIATED(edge%e2)) THEN
+          ELSEIF (ASSOCIATED(randomEdge%e2)) THEN
-            part%cell = edge%e2%n
+            partInj(n)%cell = randomEdge%e2%n
 
           ELSE
             CALL criticalError("No Volume associated to edge", 'addParticles')
 
           END IF
-          part%cellColl = edge%eColl%n
+          partInj(n)%cellColl = randomEdge%eColl%n
           sp = self%species%n
 
           !Assign particle type
-          part%species => self%species
+          partInj(n)%species => self%species
 
           if (all(self%n == 0.D0)) then
-            direction =  edge%normal
+            direction =  randomEdge%normal
 
           else
             direction = self%n
 
           end if
 
-          part%v = 0.D0
+          partInj(n)%v = 0.D0
 
-          ! do while(dot_product(partInj(n)%v, direction) <= 0.d0)
+          do while(dot_product(partInj(n)%v, direction) <= 0.d0)
-            part%v = self%vMod*direction +  (/ self%v(1)%obj%randomVel(), &
+            partInj(n)%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(part%v, direction) <= 0.D0)) then
+                (dot_product(partInj(n)%v, direction) <= 0.D0)) then
-              part%v = - part%v
+              partInj(n)%v = - partInj(n)%v
 
             end if
 
-          ! end do
+          end do
 
           !Obtain natural coordinates of particle in cell
-          part%Xi = mesh%cells(part%cell)%obj%phy2log(part%r)
+          partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r)
-
+          !Push new particle with the minimum time step
-          ! Add particle to global list
+          CALL solver%pusher(sp)%pushParticle(partInj(n), tau(sp))
-          call partInj%setLock()
+          !Assign cell to new particle
-          call partInj%add(part)
+          CALL solver%updateParticleCell(partInj(n))
-          call partInj%unsetLock()
 
         END DO
 
@@ -457,7 +500,7 @@ MODULE moduleInject
       CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
       REAL(8), INTENT(in):: temperature, m
 
-      velDist = velDistMaxwellian(vTh = DSQRT(temperature/m))
+      velDist = velDistMaxwellian(vTh = DSQRT(2.d0*temperature/m))
 
     END SUBROUTINE initVelDistMaxwellian
 
@@ -467,7 +510,7 @@ MODULE moduleInject
       CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
       REAL(8), INTENT(in):: temperature, m
 
-      velDist = velDistHalfMaxwellian(vTh = DSQRT(temperature/m))
+      velDist = velDistHalfMaxwellian(vTh = DSQRT(2.d0*temperature/m))
 
     END SUBROUTINE initVelDistHalfMaxwellian
 

src/modules/moduleList.f90

@@ -23,7 +23,6 @@ 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,8 +48,11 @@ 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,21 +235,23 @@ MODULE moduleSolver
 
       INTEGER:: nn, n, e
       INTEGER, SAVE:: nPartNew
-      INTEGER, SAVE:: nOldIn, nInj, nWScheme, nCollisions, nSurfaces
+      INTEGER, SAVE:: nInjIn, nOldIn, 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
 
@@ -265,14 +267,30 @@ MODULE moduleSolver
 
       !$OMP SINGLE
       CALL MOVE_ALLOC(partOld, partTemp)
-      nPartNew = nInj + nOldIn + nWScheme + nCollisions + nSurfaces
+      nPartNew = nInjIn + nOldIn + nWScheme + nCollisions + nSurfaces
       ALLOCATE(partOld(1:nPartNew))
       !$OMP END SINGLE
 
       !$OMP SECTIONS
       !$OMP SECTION
-      !Reset particles from previous iteration
+      !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
       DO n = 1, nPartOld
         IF (partTemp(n)%n_in) THEN
           nn = nn + 1
@@ -286,24 +304,19 @@ 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 = nOldIn + nInj
+      nn = nInjIn + nOldIn
       CALL resetList(partWScheme, partOld, nn)
 
       !$OMP SECTION
       !Reset particles from collisional process
-      nn = nOldIn + nInj + nWScheme
+      nn = nInjIn + nOldIn + nWScheme
       CALL resetList(partCollisions, partOld, nn)
 
       !$OMP SECTION
       !Reset particles from surface process
-      nn = nOldIn + nInj + nWScheme + nCollisions
+      nn = nInjIn + nOldIn + nWScheme + nCollisions
       CALL resetList(partSurfaces, partOld, nn)
 
       !$OMP SECTION
@@ -460,28 +473,6 @@ 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