New and improved method to calculate collisions per iteration:

In each iteration, number of collisions are calculate as a REAL variable
(collFrac) and stored in each cell. The number of collisions is
calculated as FLOOR(collFrac) and, if it is >1 collisions are computed
as usual. Per each collision calculated, 1.0 is removed from collFrac

runs/Argon_Expansion/CX_case.json

@@ -41,14 +41,13 @@
     "radius": 1.88e-10
   },
   "case": {
-    "tau":  [1.0e-5, 1.0e-6],
+    "tau":  [1.0e-6, 1.0e-6],
     "time": 4.0e-3,
     "pusher": ["2DCylNeutral", "2DCylNeutral"],
     "WeightingScheme": "Volume"
   },
   "interactions": {
     "folderCollisions": "./data/collisions/",
-    "tauCollisions": 1e-4,
     "collisions": [
       {"species_i": "Argon", "species_j": "Argon", 
        "cTypes": [

runs/Argon_Expansion/base_case.json

@@ -11,8 +11,8 @@
     "meshFile": "mesh.msh"
   },
   "species": [
-    {"name": "Argon",  "type": "neutral", "mass": 6.633e-26, "weight": 1.0e8, "ion": "Argon+"},
-    {"name": "Argon+", "type": "charged", "mass": 6.633e-26, "weight": 1.0e8, "neutral": "Argon"}
+    {"name": "Argon",  "type": "neutral", "mass": 6.633e-26, "weight": 1.0e8,                "ion": "Argon+"},
+    {"name": "Argon+", "type": "charged", "mass": 6.633e-26, "weight": 1.0e8, "charge": 1.0, "neutral": "Argon"}
   ],
   "boundary": [ 
     {"name": "Injection", "physicalSurface": 1, "bTypes": [
@@ -29,9 +29,9 @@
                                                  ]}
   ],
   "inject": [
-    {"name": "Exhausts Ar", "species": "Argon", "flow": 1.7, "units": "sccm",     "v": 300.0, "T": [300.0, 300.0, 300.0],
+    {"name": "Exhausts Ar", "species": "Argon", "flow": 0.7, "units": "sccm",     "v": 300.0, "T": [300.0, 300.0, 300.0],
                             "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [1, 0, 0], "physicalSurface": 1},
-    {"name": "Exhausts Ar+", "species": "Argon+", "flow": 1.3, "units": "sccm",     "v": 40000.0, "T": [300.0, 300.0, 300.0],
+    {"name": "Exhausts Ar+", "species": "Argon+", "flow": 0.3, "units": "sccm",     "v": 40000.0, "T": [300.0, 300.0, 300.0],
                              "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [1, 0, 0], "physicalSurface": 1}
   ],
   "reference": {

src/fpakc.f90

@@ -67,7 +67,7 @@ PROGRAM fpakc
     tColl = omp_get_wtime()
     !$OMP END SINGLE
 
-    CALL doCollisions(t)
+    CALL doCollisions()
 
     !$OMP SINGLE
     tColl = omp_get_wtime() - tColl

src/modules/mesh/moduleMesh.f90

@@ -134,6 +134,8 @@ MODULE moduleMesh
     INTEGER(KIND=OMP_LOCK_KIND):: lock
     !Number of collisions per volume
     INTEGER:: nColl = 0
+    !Collisional fraction
+    REAL(8):: collFrac = 0.D0
     !Total weight of particles inside cell
     REAL(8):: totalWeight = 0.D0
     CONTAINS
@@ -379,7 +381,7 @@ MODULE moduleMesh
     END SUBROUTINE findCell
 
     !Computes collisions in element
-    SUBROUTINE collision(self, t)
+    SUBROUTINE collision(self)
       USE moduleCollisions
       USE moduleSpecies
       USE moduleList
@@ -388,7 +390,6 @@ MODULE moduleMesh
       IMPLICIT NONE
 
       CLASS(meshVol), INTENT(inout):: self
-      INTEGER, INTENT(in):: t
       INTEGER:: modCollisions !Remain of current iteration and everyCollisions
       INTEGER:: iterToCollisions !Number of iterations from current to next collision
       INTEGER:: nPart !Number of particles inside the cell
@@ -401,36 +402,25 @@ MODULE moduleMesh
       REAL(8):: sigmaVrelMaxNew
       TYPE(pointerArray), ALLOCATABLE:: partTemp(:)
 
-      modCollisions = MOD(t, everyCollisions)
-      iterToCollisions = everyCollisions - modCollisions
-
       nPart = self%listPart_in%amount
+      !Computes iterations if there is more than one particle in the cell
       IF (nPart > 1) THEN
-        IF (modCollisions == 0) THEN
-          !Collisional iteration, computes the number of iterations
-          pMax = self%totalWeight*self%sigmaVrelMax*tauCollisions/self%volume
-          self%nColl = INT(REAL(nPart)*pMax*0.5D0)
+        !Probability of collision
+        pMax = self%totalWeight*self%sigmaVrelMax*tauMin/self%volume
 
-        END IF
+        !Increases the collisional fraction of the cell
+        self%collFrac = self%collFrac + REAL(nPart)*pMax*0.5D0
         
-        IF (self%nColl > iterToCollisions) THEN
-          nCollIter = self%nColl / iterToCollisions
+        !Number of collisions in the cell
+        self%nColl = FLOOR(self%collFrac)
 
-        ELSE
-          nCollIter = self%nColl
-
-        END IF
-
-        IF (nCollIter > 0) THEN
+        IF (self%nColl > 0) THEN
           !Converts the list of particles to an array for easy access
           partTemp = self%listPart_in%convert2Array()
 
-          !Removes collisions from this iteration form the total in the cell
-          self%nColl = self%nColl - nCollIter
-
         END IF
 
-        DO n = 1, nCollIter
+        DO n = 1, self%nColl
           !Select random numbers
           rnd = random(1, nPart)
           part_i => partTemp(rnd)%part
@@ -449,7 +439,11 @@ MODULE moduleMesh
 
           END IF
 
+          !Removes one collision from the collisional fraction
+          self%collFrac = self%collFrac - 1.D0
+            
         END DO
+
       END IF
 
     END SUBROUTINE collision
@@ -552,7 +546,7 @@ MODULE moduleMesh
       CHARACTER (LEN=iterationDigits):: tstring
 
 
-      IF (collOutput .AND. MOD(t, everyCollisions) == 0) THEN
+      IF (collOutput) THEN
         time = DBLE(t)*tauMin*ti_ref
         WRITE(tstring, iterationFormat) t
 

src/modules/moduleCollisions.f90

@@ -78,10 +78,6 @@ MODULE moduleCollisions
   TYPE(interactionsBinary), ALLOCATABLE, TARGET:: interactionMatrix(:)
   !Folder for collision cross section tables
   CHARACTER(:), ALLOCATABLE:: pathCollisions
-  !Time step for collisional process
-  REAL(8):: tauCollisions
-  !Number of iterations between collisional updates
-  INTEGER:: everyCollisions
 
   CONTAINS
     !Velocity of center of mass of two particles

src/modules/moduleInput.f90

@@ -172,17 +172,6 @@ MODULE moduleInput
       END IF
       tauMin = MINVAL(tau)
 
-      !Calculates iterations between collisions
-      IF (tauCollisions /= 0.D0) THEN
-        everyCollisions = INT(tauCollisions/tauMin)
-
-      ELSE
-        CALL warningError('Using minimum time step for collisions')
-        tauCollisions = tauMin
-        everyCollisions = 1
-
-      END IF
-
       !Gets the simulation time
       CALL config%get(object // '.time', time, found)
       IF (.NOT. found) CALL criticalError('Required parameter time not found','readCase')
@@ -221,7 +210,6 @@ MODULE moduleInput
       !Makes tau(s) non-dimensional
       tau = tau / ti_ref
       tauMin = tauMin / ti_ref
-      tauCollisions = tauCollisions / ti_ref
 
       !Sets the format of output files accordint to iteration number
       iterationDigits = INT(LOG10(REAL(tmax))) + 1
@@ -520,8 +508,6 @@ MODULE moduleInput
 
       !Path for collision cross-section data files
       CALL config%get('interactions.folderCollisions', pathCollisions, found)
-      !Collisional time step
-      CALL config%get('interactions.tauCollisions', tauCollisions, found)
 
       !Inits lock for list of particles
       CALL OMP_INIT_LOCK(lockCollisions)

src/modules/moduleSolver.f90

@@ -295,16 +295,15 @@ MODULE moduleSolver
     END SUBROUTINE push1DRadCharged
 
     !Do the collisions in all the cells
-    SUBROUTINE doCollisions(t)
+    SUBROUTINE doCollisions()
       USE moduleMesh
       IMPLICIT NONE
 
-      INTEGER, INTENT(in):: t
       INTEGER:: e
 
       !$OMP DO SCHEDULE(DYNAMIC)
       DO e=1, mesh%numVols
-        CALL mesh%vols(e)%obj%collision(t)
+        CALL mesh%vols(e)%obj%collision()
       END DO
       !$OMP END DO