First attempt at Coulomb collisions

First attemp for Coulomb collisions based on the moments distribtuions.
Still the method is not done and far from being complete but input
options and basic math are implemented.

src/fpakc.f90

@@ -74,7 +74,7 @@ PROGRAM fpakc
     tColl = omp_get_wtime()
     !$OMP END SINGLE
 
-    IF (doMCC) THEN
+    IF (doMCCollisions) THEN
       CALL meshForMCC%doCollisions(t)
 
     END IF
@@ -86,7 +86,7 @@ PROGRAM fpakc
     tCoul = omp_get_wTime()
     !$OMP END SINGLE
 
-    IF (ASSOCIATED(mesh%doCoulomb)) THEN
+    IF (doCoulombScattering) THEN
       CALL mesh%doCoulomb()
 
     END IF

src/makefile

@@ -4,7 +4,7 @@ OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshBoundary.o $(OBJDIR)/module
 			  	$(OBJDIR)/moduleBoundary.o $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \
 			  	$(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o $(OBJDIR)/moduleParallel.o \
 			  	$(OBJDIR)/moduleEM.o $(OBJDIR)/moduleRandom.o $(OBJDIR)/moduleMath.o \
-					$(OBJDIR)/moduleProbe.o $(OBJDIR)/moduleAverage.o \
+					$(OBJDIR)/moduleProbe.o $(OBJDIR)/moduleAverage.o $(OBJDIR)/moduleCoulomb.o \
 					$(OBJDIR)/moduleMeshInoutCommon.o \
 			  	$(OBJDIR)/moduleMeshInputVTU.o $(OBJDIR)/moduleMeshOutputVTU.o \
 			  	$(OBJDIR)/moduleMeshInputGmsh2.o $(OBJDIR)/moduleMeshOutputGmsh2.o \

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
+  REAL(8):: L_ref, v_ref, ti_ref, Vol_ref, EF_ref, Volt_ref, B_ref, e_ref
 
 END MODULE moduleRefParam
 

src/modules/init/moduleInput.f90

@@ -144,6 +144,7 @@ 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)
@@ -412,7 +413,7 @@ MODULE moduleInput
               CALL partInitial%add(partNew)
 
               !Assign particle to list in volume
-              IF (doMCC) THEN
+              IF (listInCells) THEN
                 cell => meshforMCC%cells(partNew%cellColl)%obj
                 CALL OMP_SET_LOCK(cell%lock)
                 CALL cell%listPart_in(sp)%add(partNew)
@@ -613,6 +614,7 @@ MODULE moduleInput
       USE moduleSpecies
       USE moduleList
       USE moduleCollisions
+      USE moduleCoulomb
       USE moduleErrors
       USE moduleMesh
       USE moduleCaseParam
@@ -640,10 +642,11 @@ MODULE moduleInput
       !Firstly, check if the object 'interactions' exists
       CALL config%info('interactions', found)
       IF (found) THEN
-        !Checks if MC collisions have been defined
+        !Check if MC collisions have been defined
         CALL config%info('interactions.collisions', found)
         IF (found) THEN
-          !Reads collision time step
+          doMCCollisions = .TRUE.
+          !Read collision time step
           CALL config%info('interactions.timeStep', found)
           IF (found) THEN
             CALL config%get('interactions.timeStep', tauColl, found)
@@ -752,8 +755,35 @@ MODULE moduleInput
 
         END IF
 
+        !Check if Coulomb scattering is defined
+        CALL config%info('interactions.Coulomb', found)
+        IF (found) THEN
+
+          CALL config%info('interactions.Coulomb', found, n_children = nPairs)
+          IF (nPairs > 0) THEN
+            nCoulombPairs = nPairs
+            doCoulombScattering = .TRUE.
+            ALLOCATE(coulombMatrix(1:nPairs))
+            DO i = 1, nPairs
+              WRITE(iString, '(I2)') i
+              object = 'interactions.Coulomb(' // TRIM(iString) // ')'
+              CALL config%get(object // '.species_i', species_i, found)
+              pt_i = speciesName2Index(species_i)
+              CALL config%get(object // '.species_j', species_j, found)
+              pt_j = speciesName2Index(species_j)
+
+              CALL coulombMatrix(i)%init(pt_i, pt_j)
+
+            END DO
+
+          END IF
+
+        END IF
+
       END IF
 
+      listInCells = doMCCollisions .OR. doCoulombScattering
+
     END SUBROUTINE readInteractions
 
     !Reads boundary conditions for the mesh
@@ -906,8 +936,6 @@ MODULE moduleInput
 
       END IF
 
-      doMCC = ASSOCIATED(meshForMCC)
-
       !Get the dimension of the geometry
       CALL config%get(object // '.dimension', mesh%dimen, found)
       IF (.NOT. found) THEN

src/modules/makefile

@@ -1,6 +1,6 @@
 OBJS = common.o output.o mesh.o solver.o init.o \
 			 moduleBoundary.o moduleCollisions.o moduleInject.o \
-			 moduleList.o moduleProbe.o \
+			 moduleList.o moduleProbe.o moduleCoulomb.o \
 			 moduleSpecies.o
 
 all: $(OBJS)
@@ -11,7 +11,7 @@ common.o:
 output.o: moduleSpecies.o common.o
 	$(MAKE) -C output all
 
-mesh.o: moduleCollisions.o moduleBoundary.o output.o common.o 
+mesh.o: moduleCollisions.o moduleCoulomb.o moduleBoundary.o output.o common.o 
 	$(MAKE) -C mesh all
 
 solver.o: moduleSpecies.o moduleProbe.o common.o output.o mesh.o
@@ -26,6 +26,9 @@ moduleBoundary.o: common.o moduleBoundary.f90
 moduleCollisions.o: moduleList.o moduleSpecies.o common.o moduleCollisions.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 
+moduleCoulomb.o: moduleSpecies.o common.o moduleCoulomb.f90
+	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
+
 moduleList.o: common.o moduleSpecies.o moduleList.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 

src/modules/mesh/moduleMesh.f90

@@ -393,11 +393,11 @@ MODULE moduleMesh
     !PROCEDURES SPECIFIC OF FILE TYPE
     PROCEDURE(printOutput_interface),  POINTER, PASS:: printOutput  => NULL()
     PROCEDURE(printEM_interface),      POINTER, PASS:: printEM      => NULL()
-    PROCEDURE(doCoulomb_interface),    POINTER, PASS:: doCoulomb    => NULL()
     PROCEDURE(printAverage_interface), POINTER, PASS:: printAverage => NULL()
     CONTAINS
       !GENERIC PROCEDURES
       PROCEDURE, PASS:: constructGlobalK
+      PROCEDURE, PASS:: doCoulomb
 
   END TYPE meshParticles
 
@@ -481,7 +481,11 @@ MODULE moduleMesh
   !Logical to indicate if an specific mesh for MC Collisions is used
   LOGICAL:: doubleMesh
   !Logical to indicate if MCC collisions are performed
-  LOGICAL:: doMCC
+  LOGICAL:: doMCCollisions = .FALSE.
+  !Logical to indicate if Coulomb scattering is performed
+  LOGICAL:: doCoulombScattering = .FALSE.
+  !Logica to indicate if particles have to be listed in list inside the cells
+  LOGICAL:: listInCells = .FALSE.
   !Complete path for the two meshes
   CHARACTER(:), ALLOCATABLE:: pathMeshColl, pathMeshParticle
 
@@ -650,7 +654,7 @@ MODULE moduleMesh
         part%Xi   = Xi
         part%n_in = .TRUE.
         !Assign particle to listPart_in
-        IF (doMCC) THEN
+        IF (listInCells) THEN
           CALL OMP_SET_LOCK(self%lock)
           sp = part%species%n
           CALL self%listPart_in(sp)%add(part)
@@ -725,7 +729,7 @@ MODULE moduleMesh
         Xi = cell%phy2log(part%r)
         IF (cell%inside(Xi)) THEN
           part%cellColl = cell%n
-          IF (doMCC) THEN
+          IF (listInCells) THEN
             CALL OMP_SET_LOCK(cell%lock)
             sp = part%species%n
             CALL cell%listPart_in(sp)%add(part)
@@ -946,9 +950,152 @@ MODULE moduleMesh
     END SUBROUTINE doCollisions
 
     SUBROUTINE doCoulomb(self)
+      USE moduleCoulomb
+      USE moduleRandom
+      USE moduleOutput
+      USE moduleList
+      USE moduleMath
+      USE moduleRefParam
+      USE moduleConstParam
       IMPLICIT NONE
 
-      CLASS(meshParticles), INTENT(inout):: self
+      CLASS(meshParticles), INTENT(in), TARGET:: self
+      CLASS(meshCell), POINTER:: cell
+      INTEGER:: e
+      INTEGER:: k
+      INTEGER:: i, j
+      INTEGER:: n
+      TYPE(lNode), POINTER:: partTemp
+      REAL(8):: W(3), dW(2) !Relative velocity between particle and species and its increment
+      INTEGER(8), ALLOCATABLE:: cellNodes(:)
+      CLASS(meshNode), POINTER:: node
+      TYPE(outputFormat):: output
+      REAL(8), ALLOCATABLE:: densityNodes(:), velocityNodes(:,:), temperatureNodes(:) !values in node
+      REAL(8):: density, velocity(1:3), temperature!values at particle position
+      REAL(8), DIMENSION(1:3):: e1, e2, e3
+      REAL(8):: delta_par, delta_par_square, delta_per, delta_per_square
+      REAL(8):: l, lW, AW
+      REAL(8):: rnd
+
+
+      DO e = 1, self%numCells
+        cell => self%cells(e)%obj
+        cellNodes = cell%getNodes(cell%nNodes)
+
+        ALLOCATE(densityNodes(1:cell%nNodes),       &
+                 velocityNodes(1:cell%nNodes, 1:3), &
+                 temperatureNodes(1:cell%nNodes))
+
+        DO k=1, nCoulombPairs
+          i = coulombMatrix(k)%sp_i%n
+          j = coulombMatrix(k)%sp_j%n
+
+          !Do scattering of particles from species_i due to species j
+          !Compute background properties of species_j
+          DO n = 1, cell%nNodes
+            node => self%nodes(cellNodes(n))%obj
+            CALL calculateOutput(node%output(j), output, node%v, coulombMatrix(k)%sp_j)
+            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_i
+          partTemp => cell%listPart_in(i)%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)
+            l = coulombMatrix(k)%l_j/SQRT(temperature)
+
+            W = partTemp%part%v - velocity
+            lW = l * NORM2(W)
+            AW = coulombMatrix(k)%A_ij/NORM2(W)
+            !Axis of the relative velocity
+            !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(e2, e1)
+            e3 = normalize(e3)
+
+            delta_par = -coulombMatrix(k)%A_ij*coulombMatrix(k)%one_plus_massRatio_ij*density*l**2*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) =               DABS(randomMaxwellian()*SQRT(delta_per_square*tauMin))
+
+            rnd = random()
+            partTemp%part%v = partTemp%part%v + dW(1)*e1 + dW(2)*(COS(PI2*rnd)*e2 + &
+                                                                  SIN(PI2*rnd)*e3)
+
+            partTemp => partTemp%next
+
+          END DO
+
+          IF (i /= j) THEN
+            !Do scattering of particles from species_j due to 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
+
+            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)
+              l = coulombMatrix(k)%l_i/SQRT(temperature)
+
+              W = partTemp%part%v - velocity
+              lW = l * NORM2(W)
+              AW = coulombMatrix(k)%A_ji/NORM2(W)
+              !Axis of the relative velocity
+              !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(e2, e1)
+              e3 = normalize(e3)
+
+              delta_par = -coulombMatrix(k)%A_ji*coulombMatrix(k)%one_plus_massRatio_ji*density*l**2*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) =               DABS(randomMaxwellian()*SQRT(delta_per_square*tauMin))
+
+              rnd = random()
+              partTemp%part%v = partTemp%part%v + dW(1)*e1 + dW(2)*(COS(PI2*rnd)*e2 + &
+                                                                    SIN(PI2*rnd)*e3)
+
+              partTemp => partTemp%next
+
+            END DO
+
+          END IF
+
+        END DO
+
+      END DO
 
     END SUBROUTINE doCoulomb
 

src/modules/moduleCoulomb.f90

@@ -0,0 +1,91 @@
+MODULE moduleCoulomb
+  USE moduleSpecies
+  IMPLICIT NONE
+
+  INTEGER:: nCoulombPairs = 0
+
+  !Type for Coulomb iteraction matrix
+  TYPE:: interactionsCoulomb
+    CLASS(speciesGeneric), POINTER:: sp_i
+    CLASS(speciesGeneric), POINTER:: sp_j
+    REAL(8):: one_plus_massRatio_ij, one_plus_massRatio_ji
+    REAL(8):: lnCoulomb !This can be done a function in the future
+    REAL(8):: A_ij, A_ji
+    REAL(8):: l_j, l_i
+    CONTAINS
+      PROCEDURE, PASS:: init => initInteractionCoulomb
+
+  END TYPE interactionsCoulomb
+
+  !Coulomb collision 'matrix'
+  TYPE(interactionsCoulomb), ALLOCATABLE:: coulombMatrix(:)
+
+  CONTAINS
+    PURE REAL(8) FUNCTION G(x)
+      USE moduleConstParam
+      IMPLICIT NONE
+
+      REAL(8), INTENT(in):: x
+
+      G = 0.D0
+      IF (x /= 0.D0) THEN
+        G = (ERF(x) - x*2.D0/SQRT(PI)*EXP(-x**2))/(2.D0*x**2)
+
+      END IF
+
+    END FUNCTION G
+
+    PURE REAL(8) FUNCTION H(x)
+      IMPLICIT NONE
+
+      REAL(8), INTENT(in):: x
+
+      H = ERF(x) - G(x)
+
+    END FUNCTION H
+
+    SUBROUTINE initInteractionCoulomb(self, i, j)
+      USE moduleSpecies
+      USE moduleErrors
+      USE moduleConstParam
+      USE moduleRefParam
+      IMPLICIT NONE
+
+      CLASS(interactionsCoulomb), INTENT(out):: self
+      INTEGER, INTENT(in):: i, j
+      REAL(8):: Z_i, Z_j
+
+      self%sp_i => species(i)%obj
+      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)
+        Z_i = sp%q
+
+      CLASS DEFAULT
+        CALL criticalError('Species ' // sp%name // ' for Coulomb scattering has no charge', 'initInteractionCoulomb')
+
+      END SELECT
+
+      SELECT TYPE(sp => self%sp_j)
+      TYPE IS (speciesCharged)
+        Z_j = sp%q
+
+      CLASS DEFAULT
+        CALL criticalError('Species ' // sp%name // ' for Coulomb scattering has no charge', 'initInteractionCoulomb')
+
+      END SELECT
+
+      self%lnCoulomb = 12.0
+      self%A_ij = 8.D0*PI*Z_i**2*Z_j**2*e_ref**4*self%lnCoulomb / self%sp_i%m
+      self%A_ji = 8.D0*PI*Z_j**2*Z_i**2*e_ref**4*self%lnCoulomb / self%sp_j%m
+
+      self%l_j = SQRT(self%sp_j%m / 2.D0) !Missing temperature because it's cell dependent
+      self%l_i = SQRT(self%sp_i%m / 2.D0) !Missing temperature because it's cell dependent
+
+    END SUBROUTINE initInteractionCoulomb
+
+END MODULE moduleCoulomb

src/modules/solver/moduleSolver.f90

@@ -322,7 +322,7 @@ MODULE moduleSolver
 
       !$OMP SECTION
       !Erase the list of particles inside the cell if particles have been pushed
-      IF (doMCC) THEN
+      IF (listInCells) THEN
         DO s = 1, nSpecies
           DO e = 1, mesh%numCells
             IF (solver%pusher(s)%pushSpecies) THEN 
@@ -456,7 +456,7 @@ MODULE moduleSolver
         CALL partWScheme%unsetLock()
         !Add particle to cell list
         sp = part%species%n
-        IF (doMCC) THEN
+        IF (listInCells) THEN
           CALL OMP_SET_lock(cell%lock)
           CALL cell%listPart_in(sp)%add(newPart)
           CALL OMP_UNSET_lock(cell%lock)