Added the possibility to have sub-steps

Now per each Coulomb collision process there is the possibility to do
sub-steps. This helps in improving accuracy without reducing the time
step of the problem.

src/modules/mesh/moduleMesh.f90

@@ -961,10 +961,12 @@ MODULE moduleMesh
 
       CLASS(meshParticles), INTENT(in), TARGET:: self
       CLASS(meshCell), POINTER:: cell
+      TYPE(interactionsCoulomb):: pair
       INTEGER:: e
       INTEGER:: k
       INTEGER:: i, j
       INTEGER:: n
+      INTEGER:: t
       TYPE(lNode), POINTER:: partTemp
       INTEGER(8), ALLOCATABLE:: cellNodes(:)
       CLASS(meshNode), POINTER:: node
@@ -995,14 +997,15 @@ MODULE moduleMesh
                  temperatureNodes(1:cell%nNodes))
 
         DO k=1, nCoulombPairs
-          i = coulombMatrix(k)%sp_i%n
-          j = coulombMatrix(k)%sp_j%n
+          pair = coulombMatrix(k)
+          i = pair%sp_i%n
+          j = pair%sp_j%n
 
           !Do scattering of particles from species_i due to species j
           !Compute background properties of species_j
           DO n = 1, cell%nNodes
             node => self%nodes(cellNodes(n))%obj
-            CALL calculateOutput(node%output(j), output, node%v, coulombMatrix(k)%sp_j)
+            CALL calculateOutput(node%output(j), output, node%v, pair%sp_j)
             densityNodes(n)      = output%density/n_ref
             velocityNodes(n,1:3) = output%velocity(1:3)/v_ref
             temperatureNodes(n)  = output%temperature/T_ref
@@ -1018,7 +1021,7 @@ MODULE moduleMesh
 
             !If cell temperature is too low, skip particle to avoid division by zero
             IF (temperature>eps) THEN
-              l2 = coulombMatrix(k)%l2_j/temperature
+              l2 = pair%l2_j/temperature
               l  = SQRT(l2)
 
             ELSE
@@ -1028,52 +1031,49 @@ MODULE moduleMesh
 
             END IF
 
-            C = partTemp%part%v - velocity
-            normC = NORM2(C)
+            A  = pair%A_i*density
 
-            !If relative velocity is too low, skip collision to avoid division by zero and move to next particle
-            IF (normC < eps) THEN
-              partTemp => partTemp%next
+            !Do the required substeps
+            DO t = 1, pair%nSubSteps
+              C = partTemp%part%v - velocity
+              normC = NORM2(C)
 
-              CYCLE
+              !C_3 = z; C_1, C2 = x, y (per)
+              C_per = NORM2(C(1:2))
+              cosPhi = C(1)  / C_per
+              sinPhi = C(2)  / C_per
+              cosThe = C(3)  / normC
+              sinThe = C_per / normC
 
-            END IF
+              !Rotation matrix to go from W to C
+              rotation = RESHAPE((/ cosThe*cosPhi, cosThe*sinPhi, -sinThe,    & !First column
+                                          -sinPhi,        cosPhi,    0.D0,    & !Second column
+                                    sinThe*cosPhi, sinThe*sinPhi,  cosThe /), & !Third column
+                                    (/ 3, 3 /))
 
-            !C_3 = z; C_1, C2 = x, y (per)
-            C_per = NORM2(C(1:2))
-            cosPhi = C(1)  / C_per
-            sinPhi = C(2)  / C_per
-            cosThe = C(3)  / normC
-            sinThe = C_per / normC
+              !W at start is = (0, 0, normC), so normW = normC
+              lW = l * normC
+              GlW = G(lW)
+              HlW = H(lW)
+              AW = A / normC
 
-            !Rotation matrix to go from W to C
-            rotation = RESHAPE((/ cosThe*cosPhi, cosThe*sinPhi, -sinThe,    & !First column
-                                        -sinPhi,        cosPhi,    0.D0,    & !Second column
-                                  sinThe*cosPhi, sinThe*sinPhi,  cosThe /), & !Third column
-                                  (/ 3, 3 /))
+              !Calculate changes in W due to collision process
+              deltaW_par = - A * pair%one_plus_massRatio_ij * l2 * GlW * pair%tauSub
+              deltaW_par_square = SQRT(AW * GlW * pair%tauSub)*randomMaxwellian()
+              deltaW_per_square = SQRT(AW * HlW * pair%tauSub)*randomMaxwellian()
 
-            !W at start is = (0, 0, normC), so normW = normC
-            lW = l * normC
-            GlW = G(lW)
-            HlW = H(lW)
-            A  = coulombMatrix(k)%A_i*density
-            AW = A / normC
+              !Random angle to distribute perpendicular change in velocity
+              theta_per = PI2*random()
 
-            !Calculate changes in W due to collision process
-            deltaW_par = - A * coulombMatrix(k)%one_plus_massRatio_ij * l2 * GlW * tauMin
-            deltaW_par_square = SQRT(AW * GlW * tauMin)*randomMaxwellian()
-            deltaW_per_square = SQRT(AW * HlW * tauMin)*randomMaxwellian()
+              !Change W
+              W(1) =         deltaW_per_square * COS(theta_per) 
+              W(2) =         deltaW_per_square * SIN(theta_per)
+              W(3) = normC + deltaW_par + deltaW_par_square
 
-            !Random angle to distribute perpendicular change in velocity
-            theta_per = PI2*random()
+              !Update particle velocity and return to laboratory frame
+              partTemp%part%v = MATMUL(rotation, W) + velocity
 
-            !Change W
-            W(1) =         deltaW_per_square * COS(theta_per) 
-            W(2) =         deltaW_per_square * SIN(theta_per)
-            W(3) = normC + deltaW_par + deltaW_par_square
-
-            !Update particle velocity and return to laboratory frame
-            partTemp%part%v = MATMUL(rotation, W) + velocity
+            END DO
                                                
             !Move to the next particle in the list
             partTemp => partTemp%next