Combining ij - ji collisions

In an attempt to make the operator fully conservarive I have combined ij
and ji collisions (when i/=j).

Now the matter is to find a way that makes this conserve momentum and
energy for intraspecies.

src/modules/mesh/moduleMesh.f90

@@ -967,6 +967,7 @@ MODULE moduleMesh
       INTEGER:: i, j
       INTEGER:: n
       INTEGER:: t
+      INTEGER:: p
       TYPE(lNode), POINTER:: partTemp
       INTEGER(8), ALLOCATABLE:: cellNodes(:)
       CLASS(meshNode), POINTER:: node
@@ -985,6 +986,8 @@ MODULE moduleMesh
       REAL(8):: deltaW_par, deltaW_par_square, deltaW_per_square !Increments of W
       REAL(8):: theta_per !Random angle for perpendicular direction
       REAL(8):: eps = 1.D-12
+      REAL(8):: preV(1:3), totalP_ij(1:3), totalP_ji(1:3)
+      REAL(8), ALLOCATABLE:: deltaV_ji(:,:)
 
 
       !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(partTemp)
@@ -1012,6 +1015,7 @@ MODULE moduleMesh
 
           END DO
 
+          totalP_ij = 0.D0
           !Loop over particles of species_i
           partTemp => cell%listPart_in(i)%head
           DO WHILE(ASSOCIATED(partTemp))
@@ -1071,7 +1075,9 @@ MODULE moduleMesh
               W(3) = normC + deltaW_par + deltaW_par_square
 
               !Update particle velocity and return to laboratory frame
+              preV = partTemp%part%v
               partTemp%part%v = MATMUL(rotation, W) + velocity
+              totalP_ij = totalP_ij + pair%sp_i%m*(partTemp%part%v - preV)
 
             END DO
                                                
@@ -1080,6 +1086,94 @@ MODULE moduleMesh
 
           END DO
 
+          !Do corresponding collisions
+          IF (i /= j) THEN
+            !Do scattering of particles from species_j due to species i
+            !Compute background properties of species_i
+            DO n = 1, cell%nNodes
+              node => self%nodes(cellNodes(n))%obj
+              CALL calculateOutput(node%output(i), output, node%v, pair%sp_i)
+              densityNodes(n)      = output%density/n_ref
+              velocityNodes(n,1:3) = output%velocity(1:3)/v_ref
+              temperatureNodes(n)  = output%temperature/T_ref
+
+            END DO
+
+            totalP_ji = 0.D0
+            ALLOCATE(deltaV_ji(1:cell%listPart_in(j)%amount,1:3))
+            !Loop over particles of species_j
+            partTemp => cell%listPart_in(j)%head
+            p = 1
+            DO WHILE(ASSOCIATED(partTemp))
+              density     = cell%gatherF(partTemp%part%Xi, cell%nNodes, densityNodes)
+              velocity    = cell%gatherF(partTemp%part%Xi, cell%nNodes, velocityNodes)
+              temperature = cell%gatherF(partTemp%part%Xi, cell%nNodes, temperatureNodes)
+
+              !If cell temperature is too low, skip particle to avoid division by zero
+              IF (temperature>eps) THEN
+                l2 = pair%l2_i/temperature
+                l  = SQRT(l2)
+
+              ELSE
+                partTemp => partTemp%next
+
+                CYCLE
+
+              END IF
+              A  = pair%A_j*density
+
+              !Do the required substeps
+              DO t = 1, pair%nSubSteps
+                C = partTemp%part%v - velocity
+                normC = NORM2(C)
+
+                !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
+
+                !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 /))
+
+                !W at start is = (0, 0, normC), so normW = normC
+                lW = l * normC
+                GlW = G(lW)
+                HlW = H(lW)
+                AW = A / normC
+
+                !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()
+
+                !Random angle to distribute perpendicular change in velocity
+                theta_per = PI2*random()
+
+                !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
+
+                preV = partTemp%part%v
+                partTemp%part%v = MATMUL(rotation, W) + velocity
+                totalP_ji = totalP_ji + pair%sp_j%m*(partTemp%part%v - preV)
+
+              END DO
+              
+              !Move to the next particle in the list
+              partTemp => partTemp%next
+
+            END DO
+
+          END IF
+
+          print *, k, NORM2(totalP_ij), NORM2(totalP_ji)
+
         END DO
 
         DEALLOCATE(densityNodes, velocityNodes, temperatureNodes, cellNodes)

src/modules/moduleCoulomb.f90

@@ -11,7 +11,9 @@ MODULE moduleCoulomb
     REAL(8):: one_plus_massRatio_ij
     REAL(8):: lnCoulomb !This can be done a function in the future
     REAL(8):: A_i
+    REAL(8):: A_j
     REAL(8):: l2_j
+    REAL(8):: l2_i
     REAL(8):: tauSub
     INTEGER:: nSubSteps
     CONTAINS
@@ -92,8 +94,10 @@ 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_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%l2_j = self%sp_j%m / 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
 
     END SUBROUTINE initInteractionCoulomb