I hate Coulomb and his Scattering

I found no way to ensure conservation in the linear Coulomb operator.
Thus, now two collisions have to be declared if sp_i /= sp_j: collision
ij and collision ji.

This does not conserve energy so please use under your own risk, like
everything else.

Still, I think something is wrong with this implementation and I'm
really tired.

src/modules/mesh/moduleMesh.f90

@@ -975,7 +975,7 @@ MODULE moduleMesh
       REAL(8):: delta_par, delta_par_square, delta_per, delta_per_square
       REAL(8):: W(3), dW(2), normW !Relative velocity between particle and species and its increment
       REAL(8):: l2, l, lW, AW
-      REAL(8):: deltaV(1:3), totalDeltaV_ij, normDeltaV
+      REAL(8):: deltaV(1:3)
       REAL(8):: rnd
       REAL(8):: eps = 1.D-10
 
@@ -1005,7 +1005,6 @@ MODULE moduleMesh
           END DO
 
           !Loop over particles of species_i
-          totalDeltaV_ij = 0.D0
           partTemp => cell%listPart_in(i)%head
           DO WHILE(ASSOCIATED(partTemp))
             density     = cell%gatherF(partTemp%part%Xi, cell%nNodes, densityNodes)
@@ -1026,6 +1025,7 @@ MODULE moduleMesh
 
             W = partTemp%part%v - velocity
             normW = NORM2(W)
+
             !If relative velocity is too low, skip collision to avoid division by zero and move to next particle
             IF (normW < eps) THEN
               partTemp => partTemp%next
@@ -1063,8 +1063,6 @@ MODULE moduleMesh
 
             deltaV =  dW(1)*e1 + dW(2)*(COS(rnd)*e2 + SIN(rnd)*e3)
 
-            totalDeltaV_ij = totalDeltaV_ij + NORM2(deltaV)
-
             !Change particle velocity
             partTemp%part%v = partTemp%part%v + deltaV
                                                
@@ -1072,101 +1070,6 @@ MODULE moduleMesh
 
           END DO
 
-          !Do scattering of particles from species_j due to species i
-          IF (i /= j) THEN
-            !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, coulombMatrix(k)%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
-
-            !Divide total momentum exchanged among all the particles of species j
-            !TODO: This is a dirty trick to ensure conservation between species
-            normDeltaV = totalDeltaV_ij / REAL(cell%listPart_in(j)%amount) *       &
-                         (coulombMatrix(k)%sp_i%weight*coulombMatrix(k)%sp_i%m) /  &
-                         (coulombMatrix(k)%sp_j%weight*coulombMatrix(k)%sp_j%m)
-
-            !Loop over particles of species_j
-            partTemp => cell%listPart_in(j)%head
-            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 = coulombMatrix(k)%l2_i/temperature
-                l  = SQRT(l2)
-
-              ELSE
-                partTemp => partTemp%next
-
-                CYCLE
-
-              END IF
-
-              W = partTemp%part%v - velocity
-              normW = NORM2(W)
-              !If relative velocity is too low, skip collision and move to next particle
-              IF (normW < eps) THEN
-                partTemp => partTemp%next
-
-                CYCLE
-
-              END IF
-
-              lW = l * normW
-              AW = coulombMatrix(k)%A_j/normW
-
-              delta_par = -coulombMatrix(k)%A_j*coulombMatrix(k)%one_plus_massRatio_ji*density*l2*G(lW)
-
-              delta_par_square = AW*density*G(lW)
-
-              delta_per_square = AW*density*H(lW)
-
-              dW(1) = delta_par*tauMin + randomMaxwellian()*SQRT(delta_par_square*tauMin)
-              dW(2) =                ABS(randomMaxwellian()*SQRT(delta_per_square*tauMin))
-
-              !Normalize with average exchange per particle
-              !TODO: This is a dirty trick to ensure conservation between species
-              IF (NORM2(dW) > eps) THEN
-                dW = normDeltaV*dW/NORM2(dW)
-
-              ELSE
-                dW = 0.D0
-
-              END IF
-
-              !System of reference for the velocity change
-              !First one is parallel to the relative velocity
-              e1 = normalize(W)
-              !Second one is perpendicular to it
-              e2(1) = -e1(2)
-              e2(2) =  e1(1)
-              e2(3) =  0.D0
-              e2 = normalize(e2)
-              !Third one is perpendicular to the other two
-              e3 = crossProduct(e1, e2)
-              e3 = normalize(e3)
-
-              !Random number for direction
-              rnd = PI2*random()
-
-              deltaV =  dW(1)*e1 + dW(2)*(COS(rnd)*e2 + SIN(rnd)*e3)
-
-              !Change particle velocity
-              partTemp%part%v = partTemp%part%v + deltaV
-
-              partTemp => partTemp%next
-
-            END DO
-
-          END IF
-
         END DO
 
         DEALLOCATE(densityNodes, velocityNodes, temperatureNodes, cellNodes)

src/modules/moduleCoulomb.f90

@@ -8,10 +8,10 @@ MODULE moduleCoulomb
   TYPE:: interactionsCoulomb
     CLASS(speciesGeneric), POINTER:: sp_i
     CLASS(speciesGeneric), POINTER:: sp_j
-    REAL(8):: one_plus_massRatio_ij, one_plus_massRatio_ji
+    REAL(8):: one_plus_massRatio_ij
     REAL(8):: lnCoulomb !This can be done a function in the future
-    REAL(8):: A_i, A_j
+    REAL(8):: A_i
-    REAL(8):: l2_j, l2_i
+    REAL(8):: l2_j
     CONTAINS
       PROCEDURE, PASS:: init => initInteractionCoulomb
 
@@ -60,7 +60,6 @@ MODULE moduleCoulomb
       self%sp_j => species(j)%obj
 
       self%one_plus_massRatio_ij = 1.D0 + (self%sp_i%weight*self%sp_i%m)/(self%sp_j%weight*self%sp_j%m)
-      self%one_plus_massRatio_ji = 1.D0 + (self%sp_j%weight*self%sp_j%m)/(self%sp_i%weight*self%sp_i%m)
 
       SELECT TYPE(sp => self%sp_i)
       TYPE IS (speciesCharged)
@@ -85,10 +84,8 @@ MODULE moduleCoulomb
       scaleFactor = (n_ref * qe**4) / (eps_0**2 * m_ref**2 * v_ref**3) * ti_ref
 
       self%A_i = 2.D0*Z_i**2*Z_j**2*self%lnCoulomb / self%sp_i%m**2 * scaleFactor
-      self%A_j = 2.D0*Z_j**2*Z_i**2*self%lnCoulomb / self%sp_j%m**2 * scaleFactor
 
       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