Intermediate push

I still have to figure out how to do this properly, but I am tired of
working on my laptop.

src/modules/init/moduleInput.f90

@@ -339,10 +339,9 @@ MODULE moduleInput
       !Mean velocity and temperature at particle position
       REAL(8):: velocityXi(1:3), temperatureXi
       INTEGER:: nNewPart = 0.D0
-      CLASS(meshCell), POINTER:: cell
       TYPE(particle), POINTER:: partNew
+      CLASS(meshCell), POINTER:: cell
       REAL(8):: vTh
-      TYPE(lNode), POINTER:: partCurr, partNext
 
       CALL config%info('solver.initial', found, n_children = nInitial)
 
@@ -371,9 +370,12 @@ MODULE moduleInput
             !Calculate number of particles
             nNewPart = INT(densityCen * (mesh%cells(e)%obj%volume*Vol_ref) / species(sp)%obj%weight)
 
-            !Allocate new particles
+            !Allocate array of particles for this species
+            ALLOCATE(partOld(sp)%p(1:nNewPart))
+
+            !Create new particles
             DO p = 1, nNewPart
-              ALLOCATE(partNew)
+              partNew = partOld(sp)%p(p)
               partNew%species => species(sp)%obj
               partNew%r       =  mesh%cells(e)%obj%randPos()
               partNew%Xi      =  mesh%cells(e)%obj%phy2log(partNew%r)
@@ -408,9 +410,6 @@ MODULE moduleInput
 
               partNew%weight = species(sp)%obj%weight
 
-              !Assign particle to temporal list of particles
-              CALL partInitial%add(partNew)
-
               !Assign particle to list in volume
               IF (listInCells) THEN
                 cell => meshforMCC%cells(partNew%cellColl)%obj
@@ -430,30 +429,10 @@ MODULE moduleInput
 
           END DO
 
+          nPartOld(sp) = SIZE(partOld(sp)%p)
+
         END DO
 
-        !Convert temporal list of particles into initial partOld array
-        !Deallocate the list of initial particles
-        nNewPart = partInitial%amount
-        IF (nNewPart > 0) THEN
-          ALLOCATE(partOld(1:nNewPart))
-          partCurr => partInitial%head
-          DO p = 1, nNewPart
-            partNext => partCurr%next
-            partOld(p) = partCurr%part
-            DEALLOCATE(partCurr)
-            partCurr => partNext
-
-          END DO
-
-          IF (ASSOCIATED(partInitial%head)) NULLIFY(partInitial%head)
-          IF (ASSOCIATED(partInitial%tail)) NULLIFY(partInitial%tail)
-          partInitial%amount = 0
-
-        END IF
-
-        nPartOld = SIZE(partOld)
-
       END IF
 
     END SUBROUTINE readInitial
@@ -600,7 +579,10 @@ MODULE moduleInput
 
       END DO
 
+      !Allocate the wrapper array that contains particles
+      ALLOCATE(partOld(1:nSpecies))
       !Set number of particles to 0 for init state
+      ALLOCATE(nPartOld(1:nSpecies))
       nPartOld = 0
 
       !Initialize the lock for the non-analogue (NA) list of particles

src/modules/mesh/moduleMesh.f90

@@ -884,7 +884,7 @@ MODULE moduleMesh
                   !Obtain the cross sections for the different processes
                   !TODO: From here it might be a procedure in interactionMatrix
                   vRel = NORM2(part_i%v-part_j%v)
-                  rMass = reducedMass(part_i%weight*part_i%species%m, part_j%weight*part_j%species%m)
+                  rMass = reducedMass(part_i%weight*part_i%species%mass, part_j%weight*part_j%species%mass)
                   eRel = rMass*vRel**2
                   CALL interactionMatrix(k)%getSigmaVrel(vRel, eRel, sigmaVrelTotal, sigmaVrel)
 
@@ -1079,7 +1079,7 @@ MODULE moduleMesh
 
             !Compute changes in velocity for each particle
             deltaV_ij(p,1:3) = MATMUL(rotation, W) + velocity - partTemp%part%v
-            mass_ij(p)       = pair%sp_i%m*partTemp%part%weight
+            mass_ij(p)       = pair%sp_i%mass*partTemp%part%weight
             p_ij(p,1:3)      = mass_ij(p)*partTemp%part%v
 
             !Move to the next particle in the list
@@ -1162,7 +1162,7 @@ MODULE moduleMesh
 
               !Compute changes in velocity for each particle
               deltaV_ji(p,1:3) = MATMUL(rotation, W) + velocity - partTemp%part%v
-              mass_ji(p)       = pair%sp_j%m*partTemp%part%weight
+              mass_ji(p)       = pair%sp_j%mass*partTemp%part%weight
               p_ji(p,1:3)      = mass_ji(p)*partTemp%part%v
 
               !Move to the next particle in the list

src/modules/mesh/moduleMeshBoundary.f90

@@ -125,7 +125,7 @@ MODULE moduleMeshBoundary
 
       SELECT TYPE(bound => edge%boundary%bTypes(part%species%n)%obj)
       TYPE IS(boundaryIonization)
-        mRel = reducedMass(bound%m0, part%species%m)
+        mRel = reducedMass(bound%m0, part%species%mass)
         vRel = SUM(DABS(part%v-bound%v0))
         eRel = mRel*vRel**2*5.D-1
 
@@ -238,7 +238,7 @@ MODULE moduleMeshBoundary
         !Get relative velocity
         vRel = NORM2(part%v)
         !Convert to relative energy
-        eRel = part%species%m*vRel**2*5.D-1
+        eRel = part%species%mass*vRel**2*5.D-1
 
         !If energy is abound threshold calculate secondary electrons
         IF (eRel >= bound%energyThreshold) THEN

src/modules/moduleCollisions.f90

@@ -164,8 +164,8 @@ MODULE moduleCollisions
 
       self%amount = amount
 
-      mass_i = species(i)%obj%m
+      mass_i = species(i)%obj%mass
-      mass_j = species(j)%obj%m
+      mass_j = species(j)%obj%mass
 
       ALLOCATE(self%collisions(1:self%amount))
 
@@ -227,8 +227,8 @@ MODULE moduleCollisions
       REAL(8):: m_i, m_j
       REAL(8), DIMENSION(1:3):: vCM, vp
 
-      m_i = part_i%species%m
+      m_i = part_i%species%mass
-      m_j = part_j%species%m
+      m_j = part_j%species%mass
       !Applies the collision
       vCM = velocityCM(part_i%weight*m_i, part_i%v, part_j%weight*m_j, part_j%v)
       vp  = vRel*randomDirectionVHS()
@@ -283,7 +283,7 @@ MODULE moduleCollisions
         END SELECT
 
         !momentum change per ionization process
-        collision%deltaV = sqrt(collision%eThreshold / collision%electron%m)
+        collision%deltaV = sqrt(collision%eThreshold / collision%electron%mass)
 
       END SELECT
 
@@ -307,7 +307,7 @@ MODULE moduleCollisions
       REAL(8), DIMENSION(1:3):: vChange
       TYPE(particle), POINTER:: newElectron => NULL(), remainingNeutral => NULL()
 
-      rMass = reducedMass(part_i%weight*part_i%species%m, part_j%weight*part_j%species%m)
+      rMass = reducedMass(part_i%weight*part_i%species%mass, part_j%weight*part_j%species%mass)
       eRel = rMass*vRel**2
       !Relative energy must be higher than threshold
       IF (eRel > self%eThreshold) THEN

src/modules/moduleCoulomb.f90

@@ -61,7 +61,7 @@ MODULE moduleCoulomb
       self%sp_i => species(i)%obj
       self%sp_j => species(j)%obj
 
-      self%one_plus_massRatio_ij = 1.D0 + self%sp_i%m/self%sp_j%m
+      self%one_plus_massRatio_ij = 1.D0 + self%sp_i%mass/self%sp_j%mass
 
       Z_i = 0.D0
       Z_j = 0.D0
@@ -87,11 +87,11 @@ MODULE moduleCoulomb
 
       scaleFactor = (n_ref * qe**4 * ti_ref) / (eps_0**2 * m_ref**2 * v_ref**3)
 
-      self%A_i = Z_i**2*Z_j**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_i%m**2) * scaleFactor !Missing density because it's cell dependent
+      self%A_i = Z_i**2*Z_j**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_i%mass**2) * scaleFactor !Missing density because it's cell dependent
-      self%A_j = Z_j**2*Z_i**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_j%m**2) * scaleFactor !Missing density because it's cell dependent
+      self%A_j = Z_j**2*Z_i**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_j%mass**2) * scaleFactor !Missing density because it's cell dependent
 
-      self%l2_j = self%sp_j%m / 2.D0 !Missing temperature because it's cell dependent
+      self%l2_j = self%sp_j%mass / 2.D0 !Missing temperature because it's cell dependent
-      self%l2_i = self%sp_i%m / 2.D0 !Missing temperature because it's cell dependent
+      self%l2_i = self%sp_i%mass / 2.D0 !Missing temperature because it's cell dependent
 
     END SUBROUTINE initInteractionCoulomb
 

src/modules/moduleInject.f90

@@ -179,16 +179,21 @@ MODULE moduleInject
       IMPLICIT NONE
 
       INTEGER:: i
+      INTEGER, DIMENSION(1:nInject):: nMin, nMax
       
       !$OMP SINGLE
       nPartInj = 0
       DO i = 1, nInject
         IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN 
+          nMin(i)  = nPartInj + 1 
           nPartInj = nPartInj + inject(i)%nParticles
+          nMax(i)  = nPartInj
 
         END IF
 
       END DO
+      PRINT *, nMin
+      PRINT *, nMax
 
       IF (ALLOCATED(partInj)) DEALLOCATE(partInj)
       ALLOCATE(partInj(1:nPartInj))
@@ -196,7 +201,7 @@ MODULE moduleInject
 
       DO i=1, nInject
         IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN 
-          CALL inject(i)%addParticles()
+          CALL inject(i)%addParticles(nMin(i), nMax(i))
 
         END IF
       END DO
@@ -279,8 +284,8 @@ MODULE moduleInject
       IMPLICIT NONE
 
       CLASS(injectGeneric), INTENT(in):: self
+      INTEGER, INTENT(in):: nMin, nMax !Min and Max index in partInj array
       INTEGER:: randomX
-      INTEGER, SAVE:: nMin, nMax !Min and Max index in partInj array
       INTEGER:: i
       INTEGER:: n, sp
       CLASS(meshEdge), POINTER:: randomEdge
@@ -288,17 +293,6 @@ MODULE moduleInject
 
       !Insert particles
       !$OMP SINGLE
-      nMin = 0
-      DO i = 1, self%id -1
-        IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN 
-          nMin = nMin + inject(i)%nParticles
-
-        END IF
-
-      END DO
-
-      nMin = nMin + 1
-      nMax = nMin + self%nParticles - 1
       !Assign weight to particle.
       partInj(nMin:nMax)%weight = self%species%weight
       !Particle is considered to be outside the domain

src/modules/moduleSpecies.f90

@@ -4,12 +4,56 @@ MODULE moduleSpecies
   USE OMP_LIB
   IMPLICIT NONE
   
+  !Basic type that defines a macro-particle
+  TYPE:: particle
+    REAL(8):: r(1:3)                         !Position
+    REAL(8):: v(1:3)                         !Velocity
+    CLASS(speciesGeneric), POINTER:: species !Pointer to species associated with this particle
+    INTEGER:: cell                           !Index of element in which the particle is located. TODO: Make these pointers
+    INTEGER:: cellColl                       !Index of element in which the particle is located in the Collision Mesh
+    REAL(8):: Xi(1:3)                        !Logical coordinates of particle in element e_p.
+    LOGICAL:: n_in                           !Flag that indicates if a particle is in the domain
+    REAL(8):: weight=0.D0                    !weight of particle
+
+  END TYPE particle
+
+  !Wrapper to store the particles per species
+  TYPE:: particleArray
+    TYPE(particle), ALLOCATABLE, DIMENSION(:):: p
+
+  END TYPE particleArray
+
+  !Array of pointers for the species to be pushed
+  TYPE:: particleArray_pointer
+    TYPE(particle), POINTER, DIMENSION(:):: p
+
+  END TYPE particleArray_pointer
+
+  !Number of old particles
+  INTEGER, ALLOCATABLE, DIMENSION(:):: nPartOld
+  INTEGER:: nPartOldTotal
+  !Number of injected particles
+  INTEGER:: nPartInj
+  !Arrays that contain the particles
+  TYPE(particle),      ALLOCATABLE, TARGET, DIMENSION(:):: partInj  !array of inject particles
+  TYPE(particleArray), ALLOCATABLE, TARGET, DIMENSION(:):: partOld  !array of particles from previous iteration
+  TYPE(particleArray_pointer), ALLOCATABLE, DIMENSION(:):: particlesToPush !particles pushed in each iteration
+
+  INTEGER:: nSpeciesToPush
+  INTEGER, ALLOCATABLE, DIMENSION(:):: nPartOldToPush
+
+  !Generic species type
   TYPE, ABSTRACT:: speciesGeneric
-    CHARACTER(:), ALLOCATABLE:: name
+    INTEGER:: n=0 !Index of species
-    REAL(8):: m=0.D0, weight=0.D0, qm=0.D0
+    CHARACTER(:), ALLOCATABLE:: name !Name of species
-    INTEGER:: n=0
+    !Mass, default weight of species and charge over mass
+    REAL(8):: mass=0.D0, weight=0.D0, qm=0.D0
+    INTEGER:: every !How many interations between advancing the species
+    LOGICAL:: pushSpecies !Boolean to indicate if the species is moved in the iteration
+
   END TYPE speciesGeneric
 
+  !Neutral species
   TYPE, EXTENDS(speciesGeneric):: speciesNeutral
     CLASS(speciesGeneric), POINTER:: ion => NULL()
     CONTAINS
@@ -17,6 +61,7 @@ MODULE moduleSpecies
 
   END TYPE speciesNeutral
 
+  !Charged species
   TYPE, EXTENDS(speciesGeneric):: speciesCharged
     REAL(8):: q=0.D0
     CLASS(speciesGeneric), POINTER:: ion => NULL(), neutral => NULL()
@@ -26,34 +71,17 @@ MODULE moduleSpecies
 
   END TYPE speciesCharged
 
+  !Wrapper for species
   TYPE:: speciesCont
     CLASS(speciesGeneric), ALLOCATABLE:: obj
 
   END TYPE
 
+  !Number of species
   INTEGER:: nSpecies
+  !Array for species
   TYPE(speciesCont), ALLOCATABLE, TARGET:: species(:)
 
-  TYPE particle
-    REAL(8):: r(1:3) !Position
-    REAL(8):: v(1:3) !Velocity
-    CLASS(speciesGeneric), POINTER:: species !Pointer to species associated with this particle
-    INTEGER:: cell !Index of element in which the particle is located
-    INTEGER:: cellColl !Index of element in which the particle is located in the Collision Mesh
-    REAL(8):: Xi(1:3) !Logical coordinates of particle in element e_p.
-    LOGICAL:: n_in !Flag that indicates if a particle is in the domain
-    REAL(8):: weight=0.D0 !weight of particle
-
-  END TYPE particle
-
-  !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)
       USE moduleErrors

src/modules/output/moduleOutput.f90

@@ -148,8 +148,8 @@ MODULE moduleOutput
         formatValues%density     = rawValues%den*tempVol
         formatValues%velocity(:) = tempVel
         IF (tensorTrace(tensorTemp) > 0.D0) THEN
-         formatValues%pressure    = speciesIn%m*tensorTrace(tensorTemp)*tempVol/3.D0
+          formatValues%pressure    = speciesIn%mass*tensorTrace(tensorTemp)*tempVol/3.D0
-         formatValues%temperature = formatValues%pressure/(formatValues%density*kb)
+          formatValues%temperature = formatValues%pressure/(formatValues%density*kb)
 
         END IF
       END IF

src/modules/solver/moduleSolver.f90

@@ -3,7 +3,7 @@ MODULE moduleSolver
 
   !Generic type for pusher of particles
   TYPE, PUBLIC:: pusherGeneric
-    PROCEDURE(push_interafece), POINTER, NOPASS:: pushParticle => NULL()
+    PROCEDURE(push_interface), POINTER, NOPASS:: pushParticle => NULL()
     !Boolean to indicate if the species is moved in the iteration
     LOGICAL:: pushSpecies
     !How many interations between advancing the species
@@ -29,13 +29,13 @@ MODULE moduleSolver
 
   INTERFACE
     !Push a particle
-    PURE SUBROUTINE push_interafece(part, tauIn)
+    PURE SUBROUTINE push_interface(part, tauIn)
       USE moduleSpecies
 
       TYPE(particle), INTENT(inout):: part
       REAL(8), INTENT(in):: tauIn
 
-    END SUBROUTINE push_interafece
+    END SUBROUTINE push_interface
 
     !Solve the electromagnetic field
     SUBROUTINE solveEM_interface()
@@ -169,24 +169,75 @@ MODULE moduleSolver
       USE moduleMesh
       IMPLICIT NONE
 
-      INTEGER:: n
+      INTEGER:: p
-      INTEGER:: sp
+      INTEGER:: s, sp
+      INTEGER:: e
 
-      !$OMP DO
+      !$OMP SINGLE
-      DO n = 1, nPartOld
+      !Allocate the array of particles to push
-        !Select species type
+      nSpeciesToPush = COUNT(solver%pusher(:)%pushSpecies)
-        sp = partOld(n)%species%n
+      ALLOCATE(particlesToPush(1:nSpeciesToPush))
+      ALLOCATE(nPartOldToPush(1:nSpeciesToPush))
+      !Point the arrays to the location of the particles
+      sp = 0
+      DO s = 1, nSpecies
         !Checks if the species sp is update this iteration
         IF (solver%pusher(sp)%pushSpecies) THEN
-          !Push particle
+          sp = sp + 1
-          CALL solver%pusher(sp)%pushParticle(partOld(n), tau(sp))
+          particlesToPush(sp)%p = partOld(s)%p
-          !Find cell in wich particle reside
+          nPartOldToPush(sp)    = nPartOld(s)
-          CALL solver%updateParticleCell(partOld(n))
 
         END IF
 
       END DO
-      !$OMP END DO
+
+      !Delete list of particles in cells for collisions if the particle is pushed
+      IF (listInCells) THEN
+        DO e = 1, mesh%numCells
+          DO s = 1, nSpecies
+            IF (solver%pusher(s)%pushSpecies) THEN 
+              CALL mesh%cells(e)%obj%listPart_in(s)%erase()
+              mesh%cells(e)%obj%totalWeight(s) = 0.D0
+
+            END IF
+
+          END DO
+
+        END DO
+
+      END IF
+
+      !Erase the list of particles inside the cell in coll mesh if the particle is pushed
+      IF (doubleMesh) THEN
+        DO e = 1, meshColl%numCells
+          DO s = 1, nSpecies
+            IF (solver%pusher(s)%pushSpecies) THEN 
+              CALL meshColl%cells(e)%obj%listPart_in(s)%erase()
+              meshColl%cells(e)%obj%totalWeight(s) = 0.D0
+
+            END IF
+
+          END DO
+
+        END DO
+    
+      END IF
+      !$OMP END SINGLE
+
+      !Now, push particles
+      !$OMP DO
+      DO sp = 1, nSpeciesToPush
+        DO p = 1, nPartOldToPush(sp)
+          !Push particle
+          CALL solver%pusher(sp)%pushParticle(particlesToPush(sp)%p(p), tau(sp))
+          !Find cell in which particle reside
+          CALL solver%updateParticleCell(particlesToPush(sp)%p(p))
+
+        END DO
+
+      END DO
+      !$END OMP DO
+
 
     END SUBROUTINE doPushes
 
@@ -247,7 +298,10 @@ MODULE moduleSolver
       !$OMP SECTION
       nOldIn = 0
       IF (ALLOCATED(partOld)) THEN
-        nOldIn = COUNT(partOld%n_in)
+        DO s = 1, nSpecies
+          nOldIn = nOldin + COUNT(partOld(s)%p(:)%n_in)
+
+        END DO
 
       END IF
       !$OMP SECTION
@@ -324,40 +378,6 @@ MODULE moduleSolver
 
       END DO
 
-      !$OMP SECTION
-      !Erase the list of particles inside the cell if particles have been pushed
-      IF (listInCells) THEN
-        DO s = 1, nSpecies
-          DO e = 1, mesh%numCells
-            IF (solver%pusher(s)%pushSpecies) THEN 
-              CALL mesh%cells(e)%obj%listPart_in(s)%erase()
-              mesh%cells(e)%obj%totalWeight(s) = 0.D0
-
-            END IF
-
-          END DO
-
-        END DO
-
-      END IF
-
-      !$OMP SECTION
-      !Erase the list of particles inside the cell in coll mesh
-      IF (doubleMesh) THEN
-        DO s = 1, nSpecies
-          DO e = 1, meshColl%numCells
-            IF (solver%pusher(s)%pushSpecies) THEN 
-              CALL meshColl%cells(e)%obj%listPart_in(s)%erase()
-              meshColl%cells(e)%obj%totalWeight(s) = 0.D0
-
-            END IF
-
-          END DO
-
-        END DO
-      
-      END IF
-
       !$OMP END SECTIONS
 
       !$OMP SINGLE
@@ -372,14 +392,17 @@ MODULE moduleSolver
       USE moduleMesh
       IMPLICIT NONE
 
-      INTEGER:: n
+      INTEGER:: n, sp
       CLASS(meshCell), POINTER:: cell
 
       !Loops over the particles to scatter them
       !$OMP DO
-      DO n = 1, nPartOld
+      DO sp = 1, nSpeciesToPush
-        cell => mesh%cells(partOld(n)%cell)%obj
+        DO n = 1, nPartOldToPush(sp)
-        CALL cell%scatter(cell%nNodes, partOld(n))
+          cell => mesh%cells(particlesToPush(sp)%p(n)%cell)%obj
+          CALL cell%scatter(cell%nNodes, particlesToPush(sp)%p(n))
+
+        END DO
 
       END DO
       !$OMP END DO
@@ -549,8 +572,8 @@ MODULE moduleSolver
 
         END IF
 
-        IF (nPartOld > 0) THEN
+        IF (nPartOldTotal > 0) THEN
-          WRITE(*, "(5X,A21,F8.1,A2)") "avg t/particle: ", 1.D9*tStep/DBLE(nPartOld), "ns"
+          WRITE(*, "(5X,A21,F8.1,A2)") "avg t/particle: ", 1.D9*tStep/DBLE(nPartOldTotal), "ns"
 
         END IF
         WRITE(*,*)

src/modules/solver/pusher/modulePusher.f90

@@ -14,7 +14,7 @@ MODULE modulePusher
     END SUBROUTINE pushCartNeutral
 
     PURE SUBROUTINE pushCartElectrostatic(part, tauIn)
-      USE moduleSPecies
+      USE moduleSpecies
       USE moduleMesh
       IMPLICIT NONE