No progress in fixing Coulomb collisions with mass ratio

I am starting to think that the only fix is to reduce the time step, but
that is too harsh.

src/modules/mesh/moduleMesh.f90

@@ -978,6 +978,7 @@ MODULE moduleMesh
       REAL(8):: rnd
 
 
+      !$OMP SINGLE
       DO e = 1, self%numCells
         cell => self%cells(e)%obj
         cellNodes = cell%getNodes(cell%nNodes)
@@ -1016,6 +1017,7 @@ MODULE moduleMesh
             lW = l * normW
             AW = coulombMatrix(k)%A_i/normW
 
+
             delta_par = -coulombMatrix(k)%A_i*coulombMatrix(k)%one_plus_massRatio_ij*density*l2*G(lW)
 
             delta_par_square = AW*density*G(lW)
@@ -1050,38 +1052,40 @@ MODULE moduleMesh
 
           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, 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
-          
+
+            !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)
-          
+
               l2 = coulombMatrix(k)%l2_i/temperature
               l  = SQRT(l2)
-          
+  
               W = partTemp%part%v - velocity
               normW = NORM2(W)
               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))
-          
+
               !System of reference for the velocity change
               !First one is parallel to the relative velocity
               e1 = normalize(W)
@@ -1093,20 +1097,24 @@ MODULE moduleMesh
               !Third one is perpendicular to the other two
               e3 = crossProduct(e2, e1)
               e3 = normalize(e3)
-          
+  
+              !Random number for direction
               rnd = PI2*random()
+  
+              !Change particle velocity
               partTemp%part%v = partTemp%part%v + dW(1)*e1 + dW(2)*(COS(rnd)*e2 + &
                                                                     SIN(rnd)*e3)
-          
+  
               partTemp => partTemp%next
-          
+  
             END DO
-          
+
           END IF
 
         END DO
 
       END DO
+      !$OMP END SINGLE
 
     END SUBROUTINE doCoulomb
 

src/modules/moduleCoulomb.f90

@@ -79,7 +79,7 @@ MODULE moduleCoulomb
 
       END SELECT
 
-      self%lnCoulomb = 12.0 !Make this function dependent
+      self%lnCoulomb = 10.0 !Make this function dependent
 
       self%A_i = 8.D0*PI*Z_i**2*Z_j**2*self%lnCoulomb / self%sp_i%m**2
       self%A_j = 8.D0*PI*Z_j**2*Z_i**2*self%lnCoulomb / self%sp_j%m**2