Coulomb Scattering fully conservative

Coulomb scattering is now fully conservative thanks to the method in
lemos2009small.

The trick was to conserve the momentum and energy of ALL particles
involved in the scattering in each cell.

The substeps in Coulomb collisions have been removed as they are no
longer necessary.

Still some issues with e-i, but I don't know right now.

src/modules/init/moduleInput.f90

@@ -637,7 +637,6 @@ MODULE moduleInput
       CHARACTER(:), ALLOCATABLE:: electron
       INTEGER:: e
       CLASS(meshCell), POINTER:: cell
-      INTEGER:: subSteps
 
       !Firstly, check if the object 'interactions' exists
       CALL config%info('interactions', found)
@@ -771,13 +770,8 @@ MODULE moduleInput
               pt_i = speciesName2Index(species_i)
               CALL config%get(object // '.species_j', species_j, found)
               pt_j = speciesName2Index(species_j)
-              CALL config%get(object // '.subSteps', subSteps, found)
-              IF (.NOT. found) THEN
-                subSteps = 1
 
-              END IF
+              CALL coulombMatrix(i)%init(pt_i, pt_j)
-
-              CALL coulombMatrix(i)%init(pt_i, pt_j, subSteps)
 
             END DO
 

src/modules/mesh/moduleMesh.f90

@@ -966,7 +966,6 @@ MODULE moduleMesh
       INTEGER:: k
       INTEGER:: i, j
       INTEGER:: n
-      INTEGER:: t
       INTEGER:: p
       TYPE(lNode), POINTER:: partTemp
       INTEGER(8), ALLOCATABLE:: cellNodes(:)
@@ -978,7 +977,6 @@ MODULE moduleMesh
       REAL(8):: l, lW, l2
       REAL(8):: GlW, HlW
       REAL(8):: normC
-      REAL(8):: theta, phi !angles between w and c
       REAL(8):: cosThe, sinThe
       REAL(8):: cosPhi, sinPhi
       REAL(8):: rotation(1:3,1:3) !Rotation matrix to go back to laboratory frame
@@ -986,8 +984,13 @@ 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, DIMENSION(:,:):: deltaV_ij, p_ij
-      REAL(8), ALLOCATABLE:: deltaV_ji(:,:)
+      REAL(8), ALLOCATABLE, DIMENSION(:):: mass_ij
+      REAL(8):: massSum_ij
+      REAL(8), ALLOCATABLE, DIMENSION(:,:):: deltaV_ji, p_ji
+      REAL(8), ALLOCATABLE, DIMENSION(:):: mass_ji
+      REAL(8):: massSum_ji
+      REAL(8):: alpha_num, alpha_den, alpha, beta(1:3)
 
 
       !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(partTemp)
@@ -1015,9 +1018,12 @@ MODULE moduleMesh
 
           END DO
 
-          totalP_ij = 0.D0
+          ALLOCATE(deltaV_ij(1:cell%listPart_in(i)%amount, 1:3))
+          ALLOCATE(p_ij(1:cell%listPart_in(i)%amount, 1:3))
+          ALLOCATE(mass_ij(1:cell%listPart_in(i)%amount))
           !Loop over particles of species_i
           partTemp => cell%listPart_in(i)%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)
@@ -1037,8 +1043,89 @@ MODULE moduleMesh
 
             A  = pair%A_i*density
 
-            !Do the required substeps
+            C = partTemp%part%v - velocity
-            DO t = 1, pair%nSubSteps
+            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 * tauMin
+            deltaW_par_square = SQRT(AW * GlW * tauMin)*randomMaxwellian()
+            deltaW_per_square = SQRT(AW * HlW * tauMin)*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
+
+            !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
+            p_ij(p,1:3)      = mass_ij(p)*partTemp%part%v
+
+            !Move to the next particle in the list
+            partTemp => partTemp%next
+            p = p + 1
+
+          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
+
+            ALLOCATE(deltaV_ji(1:cell%listPart_in(j)%amount, 1:3))
+            ALLOCATE(p_ji(1:cell%listPart_in(j)%amount, 1:3))
+            ALLOCATE(mass_ji(1:cell%listPart_in(j)%amount))
+            !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
+
               C = partTemp%part%v - velocity
               normC = NORM2(C)
 
@@ -1062,9 +1149,9 @@ MODULE moduleMesh
               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 = - A * pair%one_plus_massRatio_ij * l2 * GlW * tauMin
-              deltaW_par_square = SQRT(AW * GlW * pair%tauSub)*randomMaxwellian()
+              deltaW_par_square = SQRT(AW * GlW * tauMin)*randomMaxwellian()
-              deltaW_per_square = SQRT(AW * HlW * pair%tauSub)*randomMaxwellian()
+              deltaW_per_square = SQRT(AW * HlW * tauMin)*randomMaxwellian()
 
               !Random angle to distribute perpendicular change in velocity
               theta_per = PI2*random()
@@ -1074,105 +1161,90 @@ MODULE moduleMesh
               W(2) =         deltaW_per_square * SIN(theta_per)
               W(3) = normC + deltaW_par + deltaW_par_square
 
-              !Update particle velocity and return to laboratory frame
+              !Compute changes in velocity for each particle
-              preV = partTemp%part%v
+              deltaV_ji(p,1:3) = MATMUL(rotation, W) + velocity - partTemp%part%v
-              partTemp%part%v = MATMUL(rotation, W) + velocity
+              mass_ji(p)       = pair%sp_j%m*partTemp%part%weight
-              totalP_ij = totalP_ij + pair%sp_i%m*(partTemp%part%v - preV)
+              p_ji(p,1:3)      = mass_ji(p)*partTemp%part%v
 
-            END DO
-                                               
-            !Move to the next particle in the list
-            partTemp => partTemp%next
-
-          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
+              p = p + 1
 
             END DO
 
           END IF
 
-          print *, k, NORM2(totalP_ij), NORM2(totalP_ji)
+          !Calculate correction
+          !Total mass
+          massSum_ij = SUM(mass_ij)
+          massSum_ji = 0.D0
+
+          !Beta
+          beta = 0.D0
+          DO p = 1, cell%listPart_in(i)%amount
+            beta = beta + mass_ij(p) * deltaV_ij(p,1:3)
+
+          END DO
+          
+          IF (i /= j) THEN
+            massSum_ji = SUM(mass_ji)
+            DO p = 1, cell%listPart_in(j)%amount
+              beta = beta + mass_ji(p) * deltaV_ji(p,1:3)
+
+            END DO
+
+          END IF
+
+          beta = beta / (massSum_ij + massSum_ji)
+
+          !Alpha
+          alpha_num = 0.D0
+          alpha_den = 0.D0
+          DO p =1, cell%listPart_in(i)%amount
+            alpha_num = alpha_num + DOT_PRODUCT(p_ij(p,1:3), deltav_ij(p,1:3) - beta(1:3))
+            alpha_den = alpha_den + mass_ij(p) * NORM2(deltav_ij(p,1:3) - beta(1:3))**2
+
+          END DO
+
+          IF (i /= j) THEN
+            DO p = 1, cell%listPart_in(j)%amount
+              alpha_num = alpha_num + DOT_PRODUCT(p_ji(p,1:3), deltav_ji(p,1:3) - beta(1:3))
+              alpha_den = alpha_den + mass_ji(p) * NORM2(deltav_ji(p,1:3) - beta(1:3))**2
+
+            END DO
+
+          END IF
+
+          alpha = -2.D0*alpha_num / alpha_den
+
+          !Apply correction to particles velocity
+          partTemp => cell%listPart_in(i)%head
+          p = 1
+          DO WHILE(ASSOCIATED(partTemp))
+            partTemp%part%v = partTemp%part%v + alpha * (deltaV_ij(p,1:3) - beta(1:3))
+            partTemp => partTemp%next
+            p = p + 1
+
+          END DO
+
+          IF (i /= j) THEN
+            partTemp => cell%listPart_in(j)%head
+            p = 1
+            DO WHILE(ASSOCIATED(partTemp))
+              partTemp%part%v = partTemp%part%v + alpha * (deltaV_ji(p,1:3) - beta(1:3))
+              partTemp => partTemp%next
+              p = p + 1
+
+            END DO
+
+          END IF
+
+          DEALLOCATE(deltaV_ij, p_ij, mass_ij)
+
+          IF (i /= j) THEN
+            DEALLOCATE(deltaV_ji, p_ji, mass_ji)
+
+          END IF
 
         END DO
 

src/modules/moduleCoulomb.f90

@@ -14,8 +14,6 @@ MODULE moduleCoulomb
     REAL(8):: A_j
     REAL(8):: l2_j
     REAL(8):: l2_i
-    REAL(8):: tauSub
-    INTEGER:: nSubSteps
     CONTAINS
       PROCEDURE, PASS:: init => initInteractionCoulomb
 
@@ -48,7 +46,7 @@ MODULE moduleCoulomb
 
     END FUNCTION H
 
-    SUBROUTINE initInteractionCoulomb(self, i, j, subSteps)
+    SUBROUTINE initInteractionCoulomb(self, i, j)
       USE moduleSpecies
       USE moduleErrors
       USE moduleConstParam
@@ -56,14 +54,10 @@ MODULE moduleCoulomb
       IMPLICIT NONE
 
       CLASS(interactionsCoulomb), INTENT(out):: self
-      INTEGER, INTENT(in):: subSteps
       INTEGER, INTENT(in):: i, j
       REAL(8):: Z_i, Z_j
       REAL(8):: scaleFactor
 
-      self%nSubSteps = subSteps
-      self%tauSub    = tauMin / REAL(self%nSubSteps)
-
       self%sp_i => species(i)%obj
       self%sp_j => species(j)%obj