First conservative implementation of Coulomb

I am doing a trick in which I ensure that energy is conserved for
Coulomb collisions. This was not happening and what an issue for
different mass ratios. Still, this can cause an issue on getting the
right relaxation rates, still necessary to check it.

src/modules/common/moduleRefParam.f90

@@ -3,7 +3,7 @@ MODULE moduleRefParam
   !Parameters that define the problem (inputs)
   REAL(8):: n_ref, m_ref, T_ref, r_ref, debye_ref, sigmaVrel_ref
   !Reference parameters for non-dimensional problem
-  REAL(8):: L_ref, v_ref, ti_ref, Vol_ref, EF_ref, Volt_ref, B_ref, e_ref
+  REAL(8):: L_ref, v_ref, ti_ref, Vol_ref, EF_ref, Volt_ref, B_ref
 
 END MODULE moduleRefParam
 

src/modules/init/moduleInput.f90

@@ -144,7 +144,6 @@ MODULE moduleInput
       EF_ref    = qe*n_ref*L_ref/eps_0           !reference electric field
       Volt_ref  = EF_ref*L_ref                   !reference voltage
       B_ref     = m_ref / (ti_ref * qe)          !reference magnetic field
-      e_ref     = qe/(m_ref*v_ref**2)            !reference charge
 
       !If a reference cross section is given, it is used
       CALL config%get(object // '.sigmaVrel', sigmavRel_ref, found)

src/modules/mesh/moduleMesh.f90

@@ -975,10 +975,11 @@ 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):: rnd
 
 
-      !$OMP SINGLE
+      !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(partTemp)
       DO e = 1, self%numCells
         cell => self%cells(e)%obj
         cellNodes = cell%getNodes(cell%nNodes)
@@ -1003,6 +1004,7 @@ 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)
@@ -1036,15 +1038,17 @@ MODULE moduleMesh
             e2(3) =  0.D0
             e2 = normalize(e2)
             !Third one is perpendicular to the other two
-            e3 = crossProduct(e2, e1)
+            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)
+
+            totalDeltaV_ij = totalDeltaV_ij + NORM2(deltaV)
             !Change particle velocity
-            partTemp%part%v = partTemp%part%v + dW(1)*e1 + dW(2)*(COS(rnd)*e2 + &
+            partTemp%part%v = partTemp%part%v + deltaV
-                                                                  SIN(rnd)*e3)
                                                
             partTemp => partTemp%next
 
@@ -1062,6 +1066,8 @@ MODULE moduleMesh
   
             END DO
 
+            !Divide total momentum exchanged among all the particles of species j
+            normDeltaV = totalDeltaV_ij / REAL(cell%listPart_in(j)%amount) * (self%sp_i%weight*self%sp_i%m)/(self%sp_j%weight*self%sp_j%m)
             !Loop over particles of species_j
             partTemp => cell%listPart_in(j)%head
             DO WHILE(ASSOCIATED(partTemp))
@@ -1086,6 +1092,10 @@ MODULE moduleMesh
               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, it should not be neccessary but I don't lnow why it is.
+              dW = normDeltaV*dW/NORM2(dW)
+
               !System of reference for the velocity change
               !First one is parallel to the relative velocity
               e1 = normalize(W)
@@ -1095,15 +1105,16 @@ MODULE moduleMesh
               e2(3) =  0.D0
               e2 = normalize(e2)
               !Third one is perpendicular to the other two
-              e3 = crossProduct(e2, e1)
+              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 + dW(1)*e1 + dW(2)*(COS(rnd)*e2 + &
+              partTemp%part%v = partTemp%part%v + deltaV
-                                                                    SIN(rnd)*e3)
   
               partTemp => partTemp%next
   
@@ -1114,7 +1125,7 @@ MODULE moduleMesh
         END DO
 
       END DO
-      !$OMP END SINGLE
+      !$OMP END DO
 
     END SUBROUTINE doCoulomb
 

src/modules/moduleCoulomb.f90

@@ -54,6 +54,7 @@ MODULE moduleCoulomb
       CLASS(interactionsCoulomb), INTENT(out):: self
       INTEGER, INTENT(in):: i, j
       REAL(8):: Z_i, Z_j
+      REAL(8):: scaleFactor
 
       self%sp_i => species(i)%obj
       self%sp_j => species(j)%obj
@@ -81,8 +82,10 @@ MODULE moduleCoulomb
 
       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
+      scaleFactor = (n_ref * qe**4) / (eps_0**2 * m_ref**2 * v_ref**3) * ti_ref
-      self%A_j = 8.D0*PI*Z_j**2*Z_i**2*self%lnCoulomb / self%sp_j%m**2
+
+      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