Merge branch 'feature/collisionPairs' into feature/electromagnetic

Merging branches and fixing a number of important issues:

- Initial particles were not being assigned to the list of particles.

- List of particles was being erased every iteration, even if species
  was not pushed.

These caused issues with the calculation of collisions when a species
was frozen.

Now, things should work properly. All particles are properly added to
the volume list and the list is erased ONLY if the species has been
updated.

I hope that collisions are now properly accounted for per species pair.

src/modules/mesh/0D/moduleMesh0D.f90

@@ -63,6 +63,7 @@ MODULE moduleMesh0D
     !Inits dummy 0D volume
     SUBROUTINE initVol0D(self, n, p, nodes)
       USE moduleRefParam
+      USE moduleSpecies
       IMPLICIT NONE
 
       CLASS(meshVol0D), INTENT(out):: self
@@ -80,6 +81,9 @@ MODULE moduleMesh0D
 
       CALL OMP_INIT_LOCK(self%lock)
 
+      ALLOCATE(self%listPart_in(1:nSpecies))
+      ALLOCATE(self%totalWeight(1:nSpecies))
+
     END SUBROUTINE initVol0D
 
     PURE FUNCTION getNodes0D(self) RESULT(n)

src/modules/mesh/1DCart/moduleMesh1DCart.f90

@@ -220,6 +220,9 @@ MODULE moduleMesh1DCart
 
       CALL OMP_INIT_LOCK(self%lock)
 
+      ALLOCATE(self%listPart_in(1:nSpecies))
+      ALLOCATE(self%totalWeight(1:nSpecies))
+
     END SUBROUTINE initVol1DCartSegm
 
     !Calculates a random position in 1D volume

src/modules/mesh/1DRad/moduleMesh1DRad.f90

@@ -222,6 +222,9 @@ MODULE moduleMesh1DRad
 
       CALL OMP_INIT_LOCK(self%lock)
 
+      ALLOCATE(self%listPart_in(1:nSpecies))
+      ALLOCATE(self%totalWeight(1:nSpecies))
+
     END SUBROUTINE initVol1DRadSegm
 
     !Calculates a random position in 1D volume

src/modules/mesh/2DCart/moduleMesh2DCart.f90

@@ -311,6 +311,9 @@ MODULE moduleMesh2DCart
 
       CALL OMP_INIT_LOCK(self%lock)
 
+      ALLOCATE(self%listPart_in(1:nSpecies))
+      ALLOCATE(self%totalWeight(1:nSpecies))
+
     END SUBROUTINE initVolQuad2DCart
 
     !Computes element area
@@ -638,6 +641,9 @@ MODULE moduleMesh2DCart
 
       CALL OMP_INIT_LOCK(self%lock)
 
+      ALLOCATE(self%listPart_in(1:nSpecies))
+      ALLOCATE(self%totalWeight(1:nSpecies))
+
     END SUBROUTINE initVolTria2DCart
 
     !Random position in quadrilateral volume

src/modules/mesh/2DCyl/moduleMesh2DCyl.f90

@@ -299,6 +299,9 @@ MODULE moduleMesh2DCyl
 
       CALL OMP_INIT_LOCK(self%lock)
 
+      ALLOCATE(self%listPart_in(1:nSpecies))
+      ALLOCATE(self%totalWeight(1:nSpecies))
+
     END SUBROUTINE initVolQuad2DCyl
 
     !Computes element area
@@ -659,6 +662,9 @@ MODULE moduleMesh2DCyl
 
       CALL OMP_INIT_LOCK(self%lock)
 
+      ALLOCATE(self%listPart_in(1:nSpecies))
+      ALLOCATE(self%totalWeight(1:nSpecies))
+
     END SUBROUTINE initVolTria2DCyl
 
     !Random position in quadrilateral volume

src/modules/mesh/3DCart/moduleMesh3DCart.f90

@@ -285,6 +285,9 @@ MODULE moduleMesh3DCart
 
       CALL OMP_INIT_LOCK(self%lock)
 
+      ALLOCATE(self%listPart_in(1:nSpecies))
+      ALLOCATE(self%totalWeight(1:nSpecies))
+
     END SUBROUTINE initVolTetra3DCart
 
     !Random position in volume tetrahedron

src/modules/mesh/moduleMesh.f90

@@ -152,13 +152,13 @@ MODULE moduleMesh
     !Volume
     REAL(8):: volume = 0.D0
     !List of particles inside the volume
-    TYPE(listNode):: listPart_in
+    TYPE(listNode), ALLOCATABLE:: listPart_in(:)
     !Lock indicator for listPart_in
     INTEGER(KIND=OMP_LOCK_KIND):: lock
     !Number of collisions per volume
     INTEGER:: nColl = 0
     !Total weight of particles inside cell
-    REAL(8):: totalWeight = 0.D0
+    REAL(8), ALLOCATABLE:: totalWeight(:)
     CONTAINS
       PROCEDURE(initVol_interface),     DEFERRED, PASS::   init
       PROCEDURE(getNodesVol_interface), DEFERRED, PASS::   getNodes
@@ -516,6 +516,7 @@ MODULE moduleMesh
       CLASS(meshVol), OPTIONAL, INTENT(in):: oldCell
       REAL(8):: xi(1:3)
       CLASS(meshElement), POINTER:: nextElement
+      INTEGER:: sp
 
       xi = self%phy2log(part%r)
       !Checks if particle is inside 'self' cell
@@ -525,8 +526,9 @@ MODULE moduleMesh
         part%n_in = .TRUE.
         !Assign particle to listPart_in
         CALL OMP_SET_LOCK(self%lock)
-        CALL self%listPart_in%add(part)
+        sp = part%species%n
-        self%totalWeight = self%totalWeight + part%weight
+        CALL self%listPart_in(sp)%add(part)
+        self%totalWeight(sp) = self%totalWeight(sp) + part%weight
         CALL OMP_UNSET_LOCK(self%lock)
 
       ELSE
@@ -586,6 +588,7 @@ MODULE moduleMesh
       CLASS(meshVol), POINTER:: vol
       REAL(8), DIMENSION(1:3):: xii
       CLASS(meshElement), POINTER:: nextElement
+      INTEGER:: sp
 
       found = .FALSE.
       
@@ -595,8 +598,9 @@ MODULE moduleMesh
         IF (vol%inside(xii)) THEN
           part%volColl = vol%n
           CALL OMP_SET_LOCK(vol%lock)
-          CALL vol%listPart_in%add(part)
+          sp = part%species%n
-          vol%totalWeight = vol%totalWeight + part%weight
+          CALL vol%listPart_in(sp)%add(part)
+          vol%totalWeight(sp) = vol%totalWeight(sp) + part%weight
           CALL OMP_UNSET_LOCK(vol%lock)
           found = .TRUE.
 
@@ -660,13 +664,12 @@ MODULE moduleMesh
       INTEGER, INTENT(in):: t
       INTEGER:: e
       CLASS(meshVol), POINTER:: vol
-      INTEGER:: nPart !Number of particles inside the cell
+      INTEGER:: k, nPairs, i, j
+      INTEGER:: nPart_i, nPart_j, nPart!Number of particles inside the cell
       REAL(8):: pMax !Maximum probability of collision
-      TYPE(pointerArray), ALLOCATABLE:: partTemp(:)
+      TYPE(pointerArray), ALLOCATABLE:: partTemp_i(:), partTemp_j(:)
-      INTEGER:: i, j !random particle indexes
       TYPE(particle), POINTER:: part_i, part_j
-      INTEGER:: n !collision
+      INTEGER:: n, c
-      INTEGER:: ij, k
       REAL(8):: vRel, eRel
       REAL(8):: sigmaVrelTotal
       REAL(8), ALLOCATABLE:: sigmaVrel(:), probabilityColl(:)
@@ -681,80 +684,111 @@ MODULE moduleMesh
           realCollisions = 0
 
           vol => self%vols(e)%obj
-          nPart = vol%listPart_in%amount
 
-          !Resets the number of collisions
+          !TODO: Simplify this, to many sublevels
-          vol%nColl = 0
+          !Iterate over the number of pairs
+          nPairs = SIZE(interactionMatrix) !TODO: This does not change, make a variable in a module
+          DO k = 1, nPairs
+            IF (interactionMatrix(k)%amount > 0) THEN
+              !Select the species for the collision pair
+              i = interactionMatrix(k)%sp_i%n
+              j = interactionMatrix(k)%sp_j%n
 
-          !Calculates number of collisions if there is more than one particle in the cell
+              !Number of particles per species in the collision pair
-          IF (nPart > 1) THEN
+              nPart_i = vol%listPart_in(i)%amount
-            !Probability of collision
+              nPart_j = vol%listPart_in(j)%amount
-            pMax = vol%totalWeight*vol%sigmaVrelMax*tauColl/vol%volume
 
-            !Number of collisions in the cell
+              IF (nPart_i > 0 .AND. nPart_j > 0) THEN
-            vol%nColl = NINT(REAL(nPart)*pMax*0.5D0)
+                !Total number of particles for the collision pair
+                nPart = nPart_i + nPart_j
 
-            IF (vol%nColl > 0) THEN
+                !Resets the number of collisions in the cell
-              !Converts the list of particles to an array for easy access
+                vol%nColl = 0
-              partTemp = vol%listPart_in%convert2Array()
 
-            END IF
+                !Probability of collision for pair i-j
+                pMax = (vol%totalWeight(i) + vol%totalWeight(j))*vol%sigmaVrelMax*tauColl/vol%volume
 
-            DO n = 1, vol%nColl
+                !Number of collisions in the cell
-              !Select random different particles
+                vol%nColl = NINT(REAL(nPart)*pMax*0.5D0)
-              i = 0
-              j = 0
-              DO WHILE (i == j)
-                i = random(1, nPart)
-                j = random(1, nPart)
-              END DO
-              part_i => partTemp(i)%part
-              part_j => partTemp(j)%part
-              !TODO: I think that from here forward it can be passed to a procedure in interactionMatrix
-              ij = interactionIndex(part_i%species%n, part_j%species%n)
 
-              IF (interactionMatrix(ij)%amount > 0) THEN
+                !Converts the list of particles to an array for easy access
-                !Obtain the cross sections for the different processes
+                IF (vol%nColl > 0) THEN
-                vRel = NORM2(part_i%v-part_j%v)
+                  partTemp_i = vol%listPart_in(i)%convert2Array()
-                eRel = interactionMatrix(ij)%rMass*vRel**2
+                  partTemp_j = vol%listPart_in(j)%convert2Array()
-                CALL interactionMatrix(ij)%getSigmaVrel(vRel, eRel, sigmaVrelTotal, sigmaVrel)
-
-                !Update maximum sigma*v_rel
-                IF (sigmaVrelTotal > vol%sigmaVrelMax) THEN
-                  vol%sigmaVrelMax = sigmaVrelTotal
 
                 END IF
 
-                ALLOCATE(probabilityColl(0:interactionMatrix(ij)%amount))
+                DO n = 1, vol%nColl
-                probabilityColl(0) = 0.0
+                  !Select random particles
-                probabilityColl(1:interactionMatrix(ij)%amount) = sigmaVrel/vol%sigmaVrelMax
+                  part_i => NULL()
+                  part_j => NULL()
+                  rnd = random(1, nPart_i)
+                  part_i => partTemp_i(rnd)%part
+                  rnd = random(1, nPart_j)
+                  part_j => partTemp_j(rnd)%part
+                  !If they are the same particle, skip
+                  !TODO: Maybe try to improve this
+                  IF (ASSOCIATED(part_i, part_j)) THEN
+                    CYCLE
 
-                !Selects random number between 0 and 1
+                  END IF
-                rnd = random()
 
-                !If the random number is below the total probability of collision, do collisions
+                  !If particles do not belong to the species, skip collision
-                IF (rnd < sigmaVrelTotal / vol%sigmaVrelMax) THEN
+                  !This can happen, for example, if particle has been previously ionized or removed
+                  !TODO: Try to find a way to no lose these collisions. Check new 'k' and use that for the collision, maybe?
+                  IF (part_i%species%n /= i .OR. &
+                      part_j%species%n /= j) THEN
+                      CYCLE
 
-                  !Loop over collisions
+                  END IF
-                  DO k = 1, interactionMatrix(ij)%amount
-                    IF (SUM(probabilityColl(0:k-1)) + rnd <= probabilityColl(k)) THEN
-                      CALL interactionMatrix(ij)%collisions(k)%obj%collide(part_i, part_j, vRel)
-                      realCollisions = realCollisions + 1
 
-                    END IF
+                  !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)
+                  eRel = interactionMatrix(k)%rMass*vRel**2
+                  CALL interactionMatrix(k)%getSigmaVrel(vRel, eRel, sigmaVrelTotal, sigmaVrel)
 
-                  END DO
+                  !Update maximum sigma*v_rel
+                  IF (sigmaVrelTotal > vol%sigmaVrelMax) THEN
+                    vol%sigmaVrelMax = sigmaVrelTotal
 
-                END IF
+                  END IF
+
+                  ALLOCATE(probabilityColl(0:interactionMatrix(k)%amount))
+                  probabilityColl(0) = 0.0
+                  probabilityColl(1:interactionMatrix(k)%amount) = sigmaVrel/vol%sigmaVrelMax
+
+                  !Selects random number between 0 and 1
+                  rnd = random()
+
+                  !If the random number is below the total probability of collision, collide particles
+                  IF (rnd < sigmaVrelTotal / vol%sigmaVrelMax) THEN
+
+                    !Loop over collisions
+                    DO c = 1, interactionMatrix(k)%amount
+                      IF (SUM(probabilityColl(0:c-1)) + rnd <= probabilityColl(c)) THEN
+                        CALL interactionMatrix(k)%collisions(c)%obj%collide(part_i, part_j, vRel)
+                        realCollisions = realCollisions + 1
+
+                        !A collision has ocurred, exit the loop
+                        EXIT
+
+                      END IF
+
+                    END DO
+
+                  END IF
                 
-                !Deallocate arrays for next collision
+                  !Deallocate arrays for next collision
-                DEALLOCATE(sigmaVrel, probabilityColl)
+                  DEALLOCATE(sigmaVrel, probabilityColl)
+
+                END DO
 
               END IF
 
-            END DO
+            END IF
 
-          END IF
+          END DO
 
           vol%nColl = realCollisions
 

src/modules/moduleCollisions.f90

@@ -67,6 +67,8 @@ MODULE moduleCollisions
 
   !Type for interaction matrix
   TYPE:: interactionsBinary
+    CLASS(speciesGeneric), POINTER:: sp_i
+    CLASS(speciesGeneric), POINTER:: sp_j
     INTEGER:: amount
     REAL(8):: rMass !Reduced mass
     TYPE(collisionCont), ALLOCATABLE:: collisions(:)
@@ -142,16 +144,24 @@ MODULE moduleCollisions
     END FUNCTION interactionIndex
 
     !Inits the binary interaction
-    SUBROUTINE initInteractionBinary(self, amount, mass_i, mass_j)
+    SUBROUTINE initInteractionBinary(self, amount,  i, j)
       USE moduleMath
+      USE moduleSpecies
       IMPLICIT NONE
 
       CLASS(interactionsBinary), INTENT(inout):: self
       INTEGER, INTENT(in):: amount
-      REAL(8), INTENT(in):: mass_i, mass_j
+      INTEGER, INTENT(in):: i, j
+      REAL(8):: mass_i, mass_j
+
+      self%sp_i => species(i)%obj
+      self%sp_j => species(j)%obj
 
       self%amount = amount
 
+      mass_i = species(i)%obj%m
+      mass_j = species(j)%obj%m
+
       self%rMass = reducedMass(mass_i, mass_j)
 
       ALLOCATE(self%collisions(1:self%amount))
@@ -165,14 +175,14 @@ MODULE moduleCollisions
       REAL(8), INTENT(in):: vRel, eRel
       REAL(8), INTENT(out):: sigmaVrelTotal
       REAL(8), INTENT(out), ALLOCATABLE:: sigmaVrel(:)
-      INTEGER:: k
+      INTEGER:: c
 
       sigmaVrelTotal = 0.D0
 
       ALLOCATE(sigmaVrel(1:self%amount))
 
-      DO k = 1, self%amount
+      DO c = 1, self%amount
-        sigmaVrel(k) = self%collisions(k)%obj%crossSec%get(eRel)*vRel
+        sigmaVrel(c) = self%collisions(c)%obj%crossSec%get(eRel)*vRel
 
       END DO
       sigmaVrelTotal = SUM(sigmaVrel)
@@ -222,8 +232,13 @@ MODULE moduleCollisions
       vp  = vRel*randomDirectionVHS()
 
       !Assign velocities to particles
+      PRINT *, part_i%v
       part_i%v = vCM + m_j*vp / (m_i + m_j)
+      PRINT *, part_i%v
+      PRINT *, part_j%v
       part_j%v = vCM - m_i*vp / (m_i + m_j)
+      PRINT *, part_j%v
+      PRINT *
 
     END SUBROUTINE collideBinaryElastic
 

src/modules/moduleInput.f90

@@ -333,6 +333,7 @@ MODULE moduleInput
       !Mean velocity and temperature at particle position
       REAL(8):: velocityXi(1:3), temperatureXi
       INTEGER:: nNewPart = 0.D0
+      CLASS(meshVol), POINTER:: vol
       TYPE(particle), POINTER:: partNew
       REAL(8):: vTh
       TYPE(lNode), POINTER:: partCurr, partNext
@@ -356,7 +357,6 @@ MODULE moduleInput
             !Density at centroid of cell
             nodes = mesh%vols(e)%obj%getNodes()
             nNodes = SIZE(nodes)
-            !TODO: Procedure to obtain centroid from element (also for printing Electric Field)
             fPsi = mesh%vols(e)%obj%fPsi((/0.D0, 0.D0, 0.D0/))
             ALLOCATE(source(1:nNodes))
             DO j = 1, nNodes
@@ -430,6 +430,13 @@ MODULE moduleInput
               !Assign particle to temporal list of particles
               CALL partInitial%add(partNew)
 
+              !Assign particle to list in volume
+              vol => meshforMCC%vols(partNew%volColl)%obj
+              CALL OMP_SET_LOCK(vol%lock)
+              CALL vol%listPart_in(sp)%add(partNew)
+              vol%totalWeight(sp) = vol%totalWeight(sp) + partNew%weight
+              CALL OMP_UNSET_LOCK(vol%lock)
+
             END DO
 
             DEALLOCATE(source)
@@ -685,7 +692,7 @@ MODULE moduleInput
             CALL config%info(object // '.cTypes', found, n_children = nCollisions)
             ij = interactionIndex(pt_i,pt_j)
             !Allocates the required number of collisions per each pair of species ij
-            CALL interactionMatrix(ij)%init(nCollisions, species(pt_i)%obj%m, species(pt_j)%obj%m)
+            CALL interactionMatrix(ij)%init(nCollisions, pt_i, pt_j)
 
             DO k = 1, nCollisions
               WRITE (kString, '(I2)') k

src/modules/moduleList.f90

@@ -71,6 +71,7 @@ MODULE moduleList
       DO n=1, self%amount
         !Point element in array to element in list
         partArray(n)%part => tempNode%part
+
         !Go to next element
         tempNode => tempNode%next
 

src/modules/moduleSolver.f90

@@ -457,6 +457,7 @@ MODULE moduleSolver
       INTEGER, SAVE:: nPartNew
       INTEGER, SAVE:: nInjIn, nOldIn, nWScheme, nCollisions, nSurfaces
       TYPE(particle), ALLOCATABLE, SAVE:: partTemp(:)
+      INTEGER:: s
 
       !$OMP SECTIONS
       !$OMP SECTION
@@ -540,18 +541,30 @@ MODULE moduleSolver
       END DO
 
       !$OMP SECTION
-      !Erase the list of particles inside the cell
+      !Erase the list of particles inside the cell if particles have been pushed
-      DO e = 1, mesh%numVols
+      DO s = 1, nSpecies
-        mesh%vols(e)%obj%totalWeight = 0.D0
+        DO e = 1, mesh%numVols
-        CALL mesh%vols(e)%obj%listPart_in%erase()
+          IF (solver%pusher(s)%pushSpecies) THEN 
+            CALL mesh%vols(e)%obj%listPart_in(s)%erase()
+            mesh%vols(e)%obj%totalWeight(s) = 0.D0
+
+          END IF
+
+        END DO
 
       END DO
 
       !$OMP SECTION
       !Erase the list of particles inside the cell in coll mesh
-      DO e = 1, meshColl%numVols
+      DO s = 1, nSpecies
-        meshColl%vols(e)%obj%totalWeight = 0.D0
+        DO e = 1, meshColl%numVols
-        CALL meshColl%vols(e)%obj%listPart_in%erase()
+          IF (solver%pusher(s)%pushSpecies) THEN 
+            CALL meshColl%vols(e)%obj%listPart_in(s)%erase()
+            meshColl%vols(e)%obj%totalWeight(s) = 0.D0
+
+          END IF
+
+        END DO
 
       END DO
 
@@ -685,6 +698,7 @@ MODULE moduleSolver
       REAL(8):: newWeight
       TYPE(particle), POINTER:: newPart
       INTEGER:: p
+      INTEGER:: sp
 
       newWeight = part%weight / nSplit
 
@@ -697,12 +711,14 @@ MODULE moduleSolver
         ALLOCATE(newPart)
         !Copy data from original particle
         newPart = part
+        !Add particle to list of new particles from weighting scheme
         CALL OMP_SET_LOCK(lockWScheme)
         CALL partWScheme%add(newPart)
         CALL OMP_UNSET_LOCK(lockWScheme)
         !Add particle to cell list
         CALL OMP_SET_lock(vol%lock)
-        CALL vol%listPart_in%add(newPart)
+        sp = part%species%n
+        CALL vol%listPart_in(sp)%add(newPart)
         CALL OMP_UNSET_lock(vol%lock)
 
       END DO

src/modules/moduleSpecies.f90

@@ -91,7 +91,7 @@ MODULE moduleSpecies
         part%species => self%ion
 
       ELSE
-        CALL criticalError('No ion defined for species' // self%name, 'ionizeNeutral')
+        CALL criticalError('No ion defined for species :' // self%name, 'ionizeNeutral')
 
       END IF
 
@@ -109,7 +109,7 @@ MODULE moduleSpecies
         part%species => self%ion
 
       ELSE
-        CALL criticalError('No ion defined for species' // self%name, 'ionizeCharged')
+        CALL criticalError('No ion defined for species :' // self%name, 'ionizeCharged')
 
       END IF