Implementation of ionization process.

Now collisions can have a different time step.

Added species name to output names as it was starting to get confusing
in Gmsh for multiple species.

Output filenames adapted to match any number of iterations.

data/collisions/IO_e-Ar.dat

@@ -0,0 +1,34 @@
+# H. C. Straub et. al, Physical Review A, 55,2(1995)
+# Relative energy (eV) cross section (m^2)
+17 1.700E-22
+20 4.600E-21
+25 1.240E-20
+30 1.840E-20
+35 2.260E-20
+40 2.550E-20
+45 2.660E-20
+50 2.700E-20
+55 2.690E-20
+60 2.670E-20
+65 2.670E-20
+70 2.670E-20
+75 2.660E-20
+80 2.690E-20
+85 2.700E-20
+90 2.690E-20
+95 2.670E-20
+100 2.640E-20
+110 2.610E-20
+120 2.550E-20
+140 2.450E-20
+160 2.350E-20
+180 2.270E-20
+200 2.180E-20
+225 2.100E-20
+250 1.990E-20
+275 1.870E-20
+300 1.790E-20
+350 1.630E-20
+400 1.510E-20
+450 1.390E-20
+500 1.310E-20

runs/Argon_Expansion/CX_case.json

@@ -62,7 +62,7 @@
   },
   "parallel": {
     "OpenMP":{
-      "nThreads": 8
+      "nThreads": 24
     }
   }
 }

src/fpakc.f90

@@ -7,6 +7,7 @@ PROGRAM fpakc
   USE moduleSolver
   USE moduleOutput
   USE moduleCompTime
+  USE moduleCollisions
   USE moduleMesh
   IMPLICIT NONE
 
@@ -59,7 +60,7 @@ PROGRAM fpakc
     tColl = omp_get_wtime()
     !$OMP END SINGLE
 
-    CALL doCollisions()
+    CALL doCollisions(t)
 
     !$OMP SINGLE
     tColl = omp_get_wtime() - tColl

src/modules/mesh/moduleMesh.f90

@@ -125,7 +125,7 @@ MODULE moduleMesh
     !Volume index
     INTEGER:: n = 0
     !Maximum collision rate
-    REAL(8):: sigmaVrelMax = 1.D-17
+    REAL(8):: sigmaVrelMax = 0.D0
     !Volume
     REAL(8):: volume = 0.D0
     !List of particles inside the volume
@@ -392,7 +392,7 @@ MODULE moduleMesh
       self%nColl = 0
       nPart = self%listPart_in%amount
       IF (nPart > 1) THEN
-        pMax = self%totalWeight*self%sigmaVrelMax*tauMin/self%volume
+        pMax = self%totalWeight*self%sigmaVrelMax*tauCollisions/self%volume
         self%nColl = INT(REAL(nPart)*pMax*0.5D0)
 
         !Converts the list of particles to an array for easy access
@@ -428,9 +428,6 @@ MODULE moduleMesh
 
       self%totalWeight = 0.D0
 
-      !Erase the list of particles inside the cell
-      CALL self%listPart_in%erase()
-
     END SUBROUTINE collision
 
     SUBROUTINE printOutputGmsh(self, t)
@@ -445,12 +442,12 @@ MODULE moduleMesh
       TYPE(outputFormat):: output(1:self%numNodes)
       REAL(8):: time
       CHARACTER(:), ALLOCATABLE:: fileName
-      CHARACTER (LEN=6):: tstring !TODO: Review to allow any number of iterations
+      CHARACTER (LEN=iterationDigits):: tstring
 
       time = DBLE(t)*tauMin*ti_ref
 
       DO i = 1, nSpecies
-        WRITE(tstring, '(I6.6)') t
+        WRITE(tstring, iterationFormat) t
         fileName='OUTPUT_' // tstring// '_' // species(i)%obj%name // '.msh'
         WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
         OPEN (60, file = path // folder // '/' //  fileName)
@@ -459,7 +456,7 @@ MODULE moduleMesh
         WRITE(60, "(A)") '$EndMeshFormat'
         WRITE(60, "(A)") '$NodeData'
         WRITE(60, "(A)") '1'
-        WRITE(60, "(A)") '"Density (m^-3)"'
+        WRITE(60, "(A)") '"Density' // species(i)%obj%name // ' (m^-3)"'
         WRITE(60, *) 1
         WRITE(60, *) time
         WRITE(60, *) 3
@@ -473,7 +470,7 @@ MODULE moduleMesh
         WRITE(60, "(A)") '$EndNodeData'
         WRITE(60, "(A)") '$NodeData'
         WRITE(60, "(A)") '1'
-        WRITE(60, "(A)") '"Velocity (m/s)"'
+        WRITE(60, "(A)") '"Velocity ' // species(i)%obj%name // ' (m/s)"'
         WRITE(60, *) 1
         WRITE(60, *) time
         WRITE(60, *) 3
@@ -486,7 +483,7 @@ MODULE moduleMesh
         WRITE(60, "(A)") '$EndNodeData'
         WRITE(60, "(A)") '$NodeData'
         WRITE(60, "(A)") '1'
-        WRITE(60, "(A)") '"Pressure (Pa)"'
+        WRITE(60, "(A)") '"Pressure ' // species(i)%obj%name // ' (Pa)"'
         WRITE(60, *) 1
         WRITE(60, *) time
         WRITE(60, *) 3
@@ -499,7 +496,7 @@ MODULE moduleMesh
         WRITE(60, "(A)") '$EndNodeData'
         WRITE(60, "(A)") '$NodeData'
         WRITE(60, "(A)") '1'
-        WRITE(60, "(A)") '"Temperature (K)"'
+        WRITE(60, "(A)") '"Temperature ' // species(i)%obj%name // ' (K)"'
         WRITE(60, *) 1
         WRITE(60, *) time
         WRITE(60, *) 3
@@ -519,6 +516,7 @@ MODULE moduleMesh
     SUBROUTINE printCollGmsh(self, t)
       USE moduleRefParam
       USE moduleCaseParam
+      USE moduleCollisions
       USE moduleOutput
       IMPLICIT NONE
 
@@ -527,12 +525,12 @@ MODULE moduleMesh
       INTEGER:: n
       REAL(8):: time
       CHARACTER(:), ALLOCATABLE:: fileName
-      CHARACTER (LEN=6):: tstring !TODO: Review to allow any number of iterations
+      CHARACTER (LEN=iterationDigits):: tstring
 
 
-      IF (collOutput) THEN
+      IF (collOutput .AND. MOD(t, everyCollisions) == 0) THEN
         time = DBLE(t)*tauMin*ti_ref
-        WRITE(tstring, '(I6.6)') t
+        WRITE(tstring, iterationFormat) t
 
         fileName='OUTPUT_' // tstring// '_Collisions.msh'
         WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
@@ -571,14 +569,14 @@ MODULE moduleMesh
       INTEGER:: n, e
       REAL(8):: time
       CHARACTER(:), ALLOCATABLE:: fileName
-      CHARACTER (LEN=6):: tstring !TODO: Review to allow any number of iterations
+      CHARACTER (LEN=iterationDigits):: tstring
       REAL(8):: xi(1:3)
 
       xi = (/ 0.D0, 0.D0, 0.D0 /)
 
       IF (emOutput) THEN
         time = DBLE(t)*tauMin*ti_ref
-        WRITE(tstring, '(I6.6)') t
+        WRITE(tstring, iterationFormat) t
 
         fileName='OUTPUT_' // tstring// '_EMField.msh'
         WRITE(*, "(6X,A15,A)") "Creating file: ", fileName

src/modules/moduleCollisions.f90

@@ -1,4 +1,5 @@
 MODULE moduleCollisions
+  USE moduleSpecies
   USE moduleTable
 
   !Abstract type for collision between two particles
@@ -11,14 +12,6 @@ MODULE moduleCollisions
   END TYPE collisionBinary
 
   ABSTRACT INTERFACE
-    SUBROUTINE initBinary_interface(self, crossSectionFilename, mass_i, mass_j)
-      IMPORT:: collisionBinary
-      CLASS(collisionBinary), INTENT(inout):: self
-      CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename
-      REAL(8), INTENT(in):: mass_i, mass_j
-
-    END SUBROUTINE
-
     SUBROUTINE collideBinary_interface(self, sigmaVRelMax, sigmaVrelMaxNew, part_i, part_j)
       USE moduleSpecies
       IMPORT:: collisionBinary
@@ -26,7 +19,7 @@ MODULE moduleCollisions
       CLASS(collisionBinary), INTENT(in):: self
       REAL(8), INTENT(in):: sigmaVrelMax
       REAL(8), INTENT(inout):: sigmaVrelMaxNew
-      TYPE(particle), INTENT(inout):: part_i, part_j
+      TYPE(particle), INTENT(inout), TARGET:: part_i, part_j
 
     END SUBROUTINE
 
@@ -51,6 +44,10 @@ MODULE moduleCollisions
   !Ionization binary interaction
   TYPE, EXTENDS(collisionBinary):: collisionBinaryIonization
     REAL(8):: eThreshold !Minimum energy (non-dimensional units) required for ionization
+    REAL(8):: deltaV !Change in velocity due to exchange of eThreshold
+    CLASS(speciesCharged), POINTER:: electron !Pointer to species considerer as electrons
+    REAL(8):: w_i = (1.D0+DSQRT(3.D0))/2.D0
+    REAL(8):: w_j = (DSQRT(3.D0)-1.D0)/2.D0
     CONTAINS
       PROCEDURE, PASS:: collide => collideBinaryIonization
 
@@ -76,6 +73,10 @@ MODULE moduleCollisions
   TYPE(interactionsBinary), ALLOCATABLE, TARGET:: interactionMatrix(:)
   !Folder for collision cross section tables
   CHARACTER(:), ALLOCATABLE:: pathCollisions
+  !Time step for collisional process
+  REAL(8):: tauCollisions
+  !Number of iterations between collisional updates
+  INTEGER:: everyCollisions
 
   CONTAINS
     !Inits the interaction matrix
@@ -152,7 +153,7 @@ MODULE moduleCollisions
       CLASS(collisionBinaryElastic), INTENT(in):: self
       REAL(8), INTENT(in):: sigmaVrelMax
       REAL(8), INTENT(inout):: sigmaVrelMaxNew
-      TYPE(particle), INTENT(inout):: part_i, part_j
+      TYPE(particle), INTENT(inout), TARGET:: part_i, part_j
       REAL(8):: sigmaVrel
       REAL(8):: vRel !relative velocity
       REAL(8):: eRel !relative energy
@@ -176,12 +177,12 @@ MODULE moduleCollisions
         vp_i = v_ij*self%w_j
         CALL RANDOM_NUMBER(rnd)
         alpha = PI*rnd
-        part_i%v(1) = v_i*DCOS(alpha)
+        part_i%v(1) = vp_i*DCOS(alpha)
-        part_i%v(2) = v_i*DSIN(alpha)
+        part_i%v(2) = vp_i*DSIN(alpha)
         CALL RANDOM_NUMBER(rnd)
         alpha = PI*rnd
-        part_j%v(1) = v_j*DCOS(alpha)
+        part_j%v(1) = vp_j*DCOS(alpha)
-        part_j%v(2) = v_j*DSIN(alpha)
+        part_j%v(2) = vp_j*DSIN(alpha)
 
       END IF
 
@@ -189,16 +190,20 @@ MODULE moduleCollisions
 
     !ELECTRON IMPACT IONIZATION
     !Inits electron impact ionization
-    SUBROUTINE initBinaryIonization(collision, crossSectionFilename, energyThreshold, mass_i, mass_j)
+    SUBROUTINE initBinaryIonization(collision, crossSectionFilename, energyThreshold, mass_i, mass_j, electron)
       USE moduleTable
       USE moduleRefParam
       USE moduleConstParam
+      USE moduleSpecies
+      USE moduleErrors
       IMPLICIT NONE
 
       CLASS(collisionBinary), INTENT(out), ALLOCATABLE:: collision
       CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename
       REAL(8), INTENT(in):: energyThreshold
       REAL(8), INTENT(in):: mass_i, mass_j
+      CHARACTER(:), ALLOCATABLE, INTENT(in):: electron
+      INTEGER:: electronIndex
 
       ALLOCATE(collisionBinaryIonization:: collision)
 
@@ -211,10 +216,22 @@ MODULE moduleCollisions
       !Calculates reduced mass
       collision%rMass = (mass_i*mass_j)/(mass_i+mass_j)
 
+      !Specific parameters for ionization collision
       SELECT TYPE(collision)
       TYPE IS(collisionBinaryIonization)
         !Assign the energy threshold
-        collision%eThreshold = energyThreshold/(m_ref*v_ref**2)
+        !Input energy is in eV. Convert to J with ev2J and then to
+        !non-dimensional units.
+        collision%eThreshold = energyThreshold*eV2J/(m_ref*v_ref**2)
+        electronIndex = speciesName2Index(electron)
+        SELECT TYPE(sp => species(electronIndex)%obj)
+        TYPE IS(speciesCharged)
+          collision%electron => sp
+
+        CLASS DEFAULT
+          CALL criticalError("Species " // sp%name // " chosen for ionization is not a charged species", 'initBinaryIonization')
+
+        END SELECT
 
       END SELECT
 
@@ -224,12 +241,84 @@ MODULE moduleCollisions
     SUBROUTINE collideBinaryIonization(self, sigmaVrelMax, sigmaVrelMaxNew, &
                                     part_i, part_j)
       USE moduleSpecies
+      USE moduleErrors
+      USE moduleList
+      USE moduleRefParam !TODO: DELETE
+      USE moduleConstParam !TODO: DELETE
+      USE OMP_LIB
       IMPLICIT NONE
 
       CLASS(collisionBinaryIonization), INTENT(in):: self
       REAL(8), INTENT(in):: sigmaVrelMax
       REAL(8), INTENT(inout):: sigmaVrelMaxNew
-      TYPE(particle), INTENT(inout):: part_i, part_j
+      TYPE(particle), INTENT(inout), TARGET:: part_i, part_j
+      TYPE(particle), POINTER:: electron, neutral
+      TYPE(particle), POINTER:: newElectron
+      REAL(8):: vRel, eRel
+      REAL(8):: sigmaVrel
+      REAL(8):: rnd
+      REAL(8), DIMENSION(1:3):: vp_e, vp_n
+
+      !eRel (in units of [m][L]^2[t]^-2
+      vRel = SUM(DABS(part_i%v-part_j%v)) !TODO make function of norm1
+      eRel = self%rMass*vRel**2
+      !Relative energy must be higher than threshold
+      IF (eRel > self%eThreshold) THEN
+        sigmaVrel = self%crossSec%get(eRel)*vRel
+        sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel
+        CALL RANDOM_NUMBER(rnd)
+        IF (sigmaVrelMaxNew/sigmaVrelMax > rnd) THEN
+          !Find which particle is the ionizing electron
+          IF    (part_i%sp == self%electron%sp) THEN
+            electron => part_i
+            neutral  => part_j
+
+          ELSEIF(part_j%sp == self%electron%sp) THEN
+            electron => part_j
+            neutral  => part_i
+
+          ELSE
+            CALL criticalError("No matching between input particles and ionizing species", 'collideBinaryIonization')
+
+          END IF
+
+          !Exchange energy between 
+          vp_e = electron%v*(1.D0 - self%deltaV/NORM2(electron%v))
+          vp_n = neutral%v* (1.D0 + self%deltaV/NORM2(neutral%v) )
+          
+          !Changes velocity of impacting electron
+          electron%v = vp_e
+
+          !Creates a new electron from ionization
+          ALLOCATE(newElectron)
+          newElectron%sp = electron%sp
+          newElectron%v  = vp_n
+          newElectron%r  = neutral%r
+          newElectron%xi = neutral%xi
+          newElectron%n_in = .TRUE.
+          newElectron%vol = neutral%vol
+          newElectron%weight = neutral%weight
+          newElectron%qm = electron%qm
+
+          !Ionize neutral particle
+          SELECT TYPE(sp => species(neutral%sp)%obj)
+          TYPE IS(speciesNeutral)
+            CALL sp%ionize(neutral)
+
+          CLASS DEFAULT
+            CALL criticalError(sp%name // " is not a neutral", 'collideBinaryIonization')
+
+          END SELECT
+
+          !Adds new electron to list of new particles from collisions
+          CALL OMP_SET_LOCK(lockCollisions)
+          CALL partCollisions%add(newElectron)
+          CALL OMP_UNSET_LOCK(lockCollisions)
+
+        END IF
+
+      END IF
+
 
     END SUBROUTINE collideBinaryIonization
 
@@ -266,7 +355,7 @@ MODULE moduleCollisions
       CLASS(collisionBinaryChargeExchange), INTENT(in):: self
       REAL(8), INTENT(in):: sigmaVrelMax
       REAL(8), INTENT(inout):: sigmaVrelMaxNew
-      TYPE(particle), INTENT(inout):: part_i, part_j
+      TYPE(particle), INTENT(inout), TARGET:: part_i, part_j
       REAL(8):: sigmaVrel
       REAL(8):: vRel !relative velocity
       REAL(8):: eRel !relative energy

src/modules/moduleInject.f90

@@ -88,24 +88,23 @@ MODULE moduleInject
       self%vMod       = v/v_ref
       self%n          = n
       self%T          = T/T_ref
+      self%sp  = sp
       SELECT CASE(units)
       CASE ("sccm")
         !Standard cubic centimeter per minute
-        self%nParticles = INT(flow*sccm2atomPerS*tauMin*ti_ref/species(sp)%obj%weight)
+        self%nParticles = INT(flow*sccm2atomPerS*tau(self%sp)*ti_ref/species(sp)%obj%weight)
 
       CASE ("A")
         !Input current in Ampers
-        self%nParticles = INT(flow*tauMin*ti_ref/(qe*species(sp)%obj%weight))
+        self%nParticles = INT(flow*tau(self%sp)*ti_ref/(qe*species(sp)%obj%weight))
 
       CASE DEFAULT
         CALL criticalError("No support for units: " // units, 'initInject')
         
       END SELECT
+      !Scale particles for different species steps
       IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject')
 
-      self%nParticles = self%nParticles * solver%pusher(sp)%every
-      self%sp  = sp
-
       !Gets the edge elements from which particles are injected
       !TODO: Improve this A LOT
       DO e = 1, mesh%numEdges

src/modules/moduleInput.f90

@@ -67,7 +67,7 @@ MODULE moduleInput
       IMPLICIT NONE
 
       TYPE(json_file), INTENT(inout):: config
-      LOGICAL:: found, found_r
+      LOGICAL:: found
       CHARACTER(:), ALLOCATABLE:: object
 
       object = 'reference'
@@ -81,20 +81,25 @@ MODULE moduleInput
       CALL config%get(object // '.temperature', T_ref, found)
       IF (.NOT. found) CALL criticalError('Reference temperature not found','readReference')
 
-      CALL config%get(object // '.radius', r_ref, found_r)
+      !If a reference cross section is given, it is used
+      CALL config%get(object // '.crossSection', sigma_ref, found)
 
-      !Derived parameters
+      !If not, the reference radius is searched
-      v_ref     = DSQRT(kb*T_ref/m_ref)  !reference velocity
+      IF (.NOT. found) THEN
-      IF (found_r) THEN
+        CALL config%get(object // '.radius', r_ref, found)
+        IF (found) THEN
           sigma_ref = PI*(r_ref+r_ref)**2    !reference cross section
-        L_ref = 1.D0/(sigma_ref*n_ref)     !mean free path
 
         ELSE
-        L_ref = DSQRT(kb*T_ref*eps_0/n_ref)/qe !Debye length
+          sigma_ref = 0.D0 !Assume no collisions
-        !TODO: Obtain right sigma_ref for PIC case
-        sigma_ref = PI*(4.D-10)**2        !reference cross section
 
         END IF
+
+      END IF
+      
+      !Derived parameters
+      L_ref     = DSQRT(kb*T_ref*eps_0/n_ref)/qe !reference length
+      v_ref     = DSQRT(kb*T_ref/m_ref)          !reference velocity
       ti_ref    = L_ref/v_ref                    !reference time
       Vol_ref   = L_ref**3                       !reference volume
       EF_ref    = qe*n_ref*L_ref/eps_0           !reference electric field
@@ -109,6 +114,8 @@ MODULE moduleInput
       USE moduleCaseParam
       USE moduleSolver
       USE moduleSpecies
+      USE moduleCollisions
+      USE moduleOutput
       USE json_module
       IMPLICIT NONE
 
@@ -120,6 +127,7 @@ MODULE moduleInput
       INTEGER:: nTau, nSolver
       INTEGER:: i
       CHARACTER(2):: iString
+      CHARACTER(1):: tString
 
       object = 'case'
 
@@ -139,6 +147,17 @@ MODULE moduleInput
       END IF
       tauMin = MINVAL(tau)
 
+      !Calculates iterations between collisions
+      IF (tauCollisions /= 0.D0) THEN
+        everyCollisions = INT(tauCollisions/tauMin)
+
+      ELSE
+        CALL warningError('Using minimum time step for collisions')
+        tauCollisions = tauMin
+        everyCollisions = 1
+
+      END IF
+
       !Gets the simulation time
       CALL config%get(object // '.time', time, found)
       IF (.NOT. found) CALL criticalError('Required parameter time not found','readCase')
@@ -174,9 +193,15 @@ MODULE moduleInput
       CALL solver%initWS(WSType)
 
 
-      !Makes tau non-dimensional
+      !Makes tau(s) non-dimensional
       tau = tau / ti_ref
       tauMin = tauMin / ti_ref
+      tauCollisions = tauCollisions / ti_ref
+
+      !Sets the format of output files accordint to iteration number
+      iterationDigits = INT(LOG10(REAL(tmax))) + 1
+      WRITE(tString, '(I1)') iterationDigits
+      iterationFormat = "(I" // tString // "." // tString // ")"
 
     END SUBROUTINE readCase
 
@@ -227,6 +252,7 @@ MODULE moduleInput
       USE moduleSpecies
       USE moduleErrors
       USE moduleRefParam
+      USE moduleList
       USE json_module
       IMPLICIT NONE
 
@@ -262,6 +288,7 @@ MODULE moduleInput
 
         CASE ("charged")
           CALL config%get(object // '.charge', charge,      found)
+          IF (.NOT. found) CALL criticalError("Required parameter charge not found for species " // object, 'readSpecies')
           ALLOCATE(species(i)%obj, source=speciesCharged(q  = charge, &
                                                          qm = charge/mass))
 
@@ -293,7 +320,7 @@ MODULE moduleInput
 
         TYPE IS (speciesCharged)
           !Gets species linked neutral
-          CALL config%get(object // '.neutral', linkName)
+          CALL config%get(object // '.neutral', linkName, found)
           IF (found) THEN
             linkID     =  speciesName2Index(linkName)
             sp%neutral => species(linkID)%obj
@@ -324,8 +351,10 @@ MODULE moduleInput
     !Reads information about interactions between species
     SUBROUTINE readInteractions(config)
       USE moduleSpecies
+      USE moduleList
       USE moduleCollisions
       USE moduleErrors
+      USE OMP_LIB
       USE json_module
       IMPLICIT NONE
 
@@ -341,10 +370,17 @@ MODULE moduleInput
       INTEGER:: i, k, ij
       INTEGER:: pt_i, pt_j
       REAL(8):: energyThreshold
+      CHARACTER(:), ALLOCATABLE:: electron
 
       CALL initInteractionMatrix(interactionMatrix)
 
+      !Path for collision cross-section data files
       CALL config%get('interactions.folderCollisions', pathCollisions, found)
+      !Collisional time step
+      CALL config%get('interactions.tauCollisions', tauCollisions, found)
+
+      !Inits lock for list of particles
+      CALL OMP_INIT_LOCK(lockCollisions)
 
       CALL config%info('interactions.collisions', found, n_children = nInteractions)
       DO i = 1, nInteractions
@@ -382,10 +418,12 @@ MODULE moduleInput
 
           CASE ('ionization')
             !Electorn impact ionization
-            CALL config%get(object // '.energyThreshold', found)
+            CALL config%get(object // '.energyThreshold', energyThreshold, found)
             IF (.NOT. found) CALL criticalError('energyThreshold not found for collision' // object, 'readInteractions')
+            CALL config%get(object // '.electron', electron, found)
+            IF (.NOT. found) CALL criticalError('electron not found for collision' // object, 'readInteractions')
             CALL initBinaryIonization(interactionMatrix(ij)%collisions(k)%obj, &
-                                      crossSecFilePath, energyThreshold, species(pt_i)%obj%m, species(pt_j)%obj%m)
+                                      crossSecFilePath, energyThreshold, species(pt_i)%obj%m, species(pt_j)%obj%m, electron)
 
           CASE DEFAULT
             CALL criticalError('Collision type' // cType // 'not defined yet', 'readInteractions')

src/modules/moduleList.f90

@@ -21,6 +21,9 @@ MODULE moduleList
   END TYPE listNode
 
   TYPE(listNode):: partWScheme !Particles comming from the nonAnalogue scheme
+  INTEGER(KIND=OMP_LOCK_KIND):: lockWScheme !Lock for the NA list of particles
+  TYPE(listNode):: partCollisions !Particles created in collisional process
+  INTEGER(KIND=OMP_LOCK_KIND):: lockCollisions !Lock for the NA list of particles
 
   TYPE pointerArray
     TYPE(particle), POINTER:: part

src/modules/moduleOutput.f90

@@ -22,6 +22,8 @@ MODULE moduleOutput
 
   CHARACTER(:), ALLOCATABLE:: path
   CHARACTER(:), ALLOCATABLE:: folder
+  INTEGER:: iterationDigits
+  CHARACTER(:), ALLOCATABLE:: iterationFormat
   INTEGER:: triggerOutput,  counterOutput  = 0
   INTEGER:: triggerCPUTime, counterCPUTime = 0
   LOGICAL:: timeOutput = .FALSE.
@@ -72,9 +74,6 @@ MODULE moduleOutput
       tempVol = 1.D0/(nodeVol*Vol_ref)
       IF (rawValues%den > 0.D0) THEN
         tempVel    =  rawValues%mom(:)/rawValues%den
-        IF ((tempVel(1) - 1.D0) .EQ. tempVel(1)) THEN
-          PRINT *, rawValues%mom
-        END IF
         tensorTemp = (rawValues%tensorS(:,:) - rawValues%den*outerProduct(tempVel,tempVel))
         formatValues%density     = rawValues%den*tempVol
         formatValues%velocity(:) = tempVel
@@ -106,7 +105,7 @@ MODULE moduleOutput
         IF (PRESENT(first)) THEN
           IF (first) THEN
             OPEN(20, file = path // folder // '/' // fileName, action = 'write')
-            WRITE(20, "(A1, 8X, A1, 9X, A1, 5(A20))") "#","t","n","total","push","reset","collision","weighting"
+            WRITE(20, "(A1, 8X, A1, 9X, A1, 6(A20))") "#","t","n","total","push","reset","collision","weighting","EMField"
             WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
 
           ELSE
@@ -119,7 +118,7 @@ MODULE moduleOutput
 
         END IF
 
-        WRITE (20, "(I10, I10, 5(ES20.6E3))") t, nPartOld, tStep, tPush, tReset, tColl, tWeight
+        WRITE (20, "(I10, I10, 6(ES20.6E3))") t, nPartOld, tStep, tPush, tReset, tColl, tWeight, tEMField
 
         CLOSE(20)
 

src/modules/moduleSolver.f90

@@ -295,18 +295,23 @@ MODULE moduleSolver
     END SUBROUTINE push1DRadCharged
 
     !Do the collisions in all the cells
-    SUBROUTINE doCollisions()
+    SUBROUTINE doCollisions(t)
       USE moduleMesh
+      USE moduleCollisions
       IMPLICIT NONE
 
+      INTEGER, INTENT(in):: t
       INTEGER:: e
 
+      IF (MOD(t, everyCollisions) == 0) THEN
         !$OMP DO SCHEDULE(DYNAMIC)
         DO e=1, mesh%numVols
           CALL mesh%vols(e)%obj%collision()
         END DO
         !$OMP END DO
 
+      END IF
+
     END SUBROUTINE doCollisions
     
     SUBROUTINE doReset()
@@ -317,7 +322,7 @@ MODULE moduleSolver
 
       INTEGER:: nn, n, e
       INTEGER, SAVE:: nPartNew
-      INTEGER, SAVE:: nInjIn, nOldIn, nWScheme
+      INTEGER, SAVE:: nInjIn, nOldIn, nWScheme, nCollisions
       TYPE(particle), ALLOCATABLE, SAVE:: partTemp(:)
       TYPE(lNode), POINTER:: partCurr, partNext
 
@@ -336,13 +341,15 @@ MODULE moduleSolver
       END IF
       !$OMP SECTION
       nWScheme = partWScheme%amount
+      !$OMP SECTION
+      nCollisions = partCollisions%amount
       !$OMP END SECTIONS
 
       !$OMP BARRIER
 
       !$OMP SINGLE
       CALL MOVE_ALLOC(partOld, partTemp)
-      nPartNew = nInjIn + nOldIn + nWScheme
+      nPartNew = nInjIn + nOldIn + nWScheme + nCollisions
       ALLOCATE(partOld(1:nPartNew))
       !$OMP END SINGLE
 
@@ -385,12 +392,35 @@ MODULE moduleSolver
       IF (ASSOCIATED(partWScheme%tail)) NULLIFY(partWScheme%tail)
       partWScheme%amount = 0
 
+      !$OMP SECTION
+      !Reset particles from collisional process
+      nn = nInjIn + nOldIn + nWScheme
+      partCurr => partCollisions%head
+      DO n = 1, nCollisions
+        partNext => partCurr%next
+        partOld(nn+n) = partCurr%part
+        DEALLOCATE(partCurr)
+        partCurr => partNext
+
+      END DO
+      IF (ASSOCIATED(partCollisions%head)) NULLIFY(partCollisions%head)
+      IF (ASSOCIATED(partCollisions%tail)) NULLIFY(partCollisions%tail)
+      partCollisions%amount = 0
+
       !$OMP SECTION
       !Reset output in nodes
       DO e = 1, mesh%numNodes
         CALL mesh%nodes(e)%obj%resetOutput()
 
       END DO
+
+      !$OMP SECTION
+      !Erase the list of particles inside the cell
+      DO e = 1, mesh%numVols
+        CALL mesh%vols(e)%obj%listPart_in%erase()
+
+      END DO
+
       !$OMP END SECTIONS
 
       !$OMP SINGLE

src/modules/moduleSpecies.f90

@@ -52,7 +52,6 @@ MODULE moduleSpecies
   !Arrays that contain the particles
   TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: partOld !array of particles from previous iteration
   TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: partInj !array of inject particles
-  INTEGER(KIND=OMP_LOCK_KIND):: lockWScheme !Lock for the NA list of particles
 
   CONTAINS
     FUNCTION speciesName2Index(speciesName) RESULT(sp)

src/modules/moduleTable.f90

@@ -113,7 +113,11 @@ MODULE moduleTable
       REAL(8):: data_x, data_f
 
       self%x = self%x * data_x
+      self%xMin = self%xMin * data_x
+      self%xMax = self%xMax * data_x
       self%f = self%f * data_f
+      self%fMin = self%fMin * data_f
+      self%fMax = self%fMax * data_f
       self%k = self%k * data_f / data_x
 
     END SUBROUTINE convertUnits