Alphie grid case and issues

Output for the example ALPHIE_Grid.

Found an issue when multiple injections were used with species with
different time steps.

Modification to the way to compute the ionization boundary:
  The maximum number of ionizations is computed by eRel/eThreshold
  (relative energy / threshold of ionization)
  For each possible ionization, the probability of ionization is
  computed based on the density of neutrals, cross section and effective
  time divided by the number of maximum ionizations.
  If an ionization takes place, the ionization energy is substracted
  from the relative energy.

runs/ALPHIE_Grid/inputBaseCase.json

@@ -46,14 +46,15 @@
   "boundaryEM": [
     {"name": "Extraction Grid",    "type": "dirichlet", "potential": -150.0, "physicalSurface": 4},
     {"name": "Acceleration Grid",  "type": "dirichlet", "potential": -600.0, "physicalSurface": 5},
-    {"name": "Ionization Chamber", "type": "dirichlet", "potential":    0.0, "physicalSurface": 1}
+    {"name": "Ionization Chamber", "type": "dirichlet", "potential":    0.0, "physicalSurface": 1},
+    {"name": "Infinite",           "type": "dirichlet", "potential": -600.0, "physicalSurface": 2}
   ],
   "inject": [
-    {"name": "Ionization Argon+",   "species": "Argon+",   "flow": 27.0e-6, "units": "A", "v":    322.0, "T": [ 500.0,  500.0,  500.0],
+    {"name": "Ionization Argon+",   "species": "Argon+",   "flow":  1.0e-5, "units": "A", "v":   2500.0, "T": [ 500.0,  500.0,  500.0],
                                     "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [ 1, 0, 0], "physicalSurface": 1},
-    {"name": "Ionization Electron", "species": "Electron", "flow": 27.0e-6, "units": "A", "v":  87000.0, "T": [ 500.0,  500.0,  500.0],
+    {"name": "Ionization Electron", "species": "Electron", "flow":  1.0e-5, "units": "A", "v":  87000.0, "T": [30000.0, 30000.0, 30000.0],
                                     "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [ 1, 0, 0], "physicalSurface": 1},
-    {"name": "Cathode Electron",    "species": "Electron", "flow":  9.0e-5, "units": "A", "v":  87000.0, "T": [2500.0, 2500.0, 2500.0],
+    {"name": "Cathode Electron",    "species": "Electron", "flow":  1.0e-4, "units": "A", "v":  87000.0, "T": [30000.0, 30000.0, 30000.0],
                                     "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [-1, 0, 0], "physicalSurface": 2}
   ],
   "reference": {

runs/ALPHIE_Grid/inputIonization_0.10mA.json

@@ -20,7 +20,7 @@
   "boundary": [ 
     {"name": "Ionization Chanber", "physicalSurface": 1, "bTypes": [
                                                            {"type": "transparent"},
-                                                           {"type": "ionization", "neutral": {"ion": "Argon+", "mass": 6.633e-26, "density": 1.0e17, "velocity": [323, 0, 0], "temperature": 300},
+                                                           {"type": "ionization", "neutral": {"ion": "Argon+", "mass": 6.633e-26, "density": 5.0e16, "velocity": [2500, 0, 0], "temperature": 300},
                                                                                   "effectiveTime": 5.0e-6,"energyThreshold": 15.76, "crossSection": "./data/collisions/IO_e-Ar.dat"}
                                                           ]},
     {"name": "Vacuum Chamber",     "physicalSurface": 2, "bTypes": [
@@ -47,10 +47,11 @@
   "boundaryEM": [
     {"name": "Extraction Grid",    "type": "dirichlet", "potential": -150.0, "physicalSurface": 4},
     {"name": "Acceleration Grid",  "type": "dirichlet", "potential": -600.0, "physicalSurface": 5},
-    {"name": "Ionization Chamber", "type": "dirichlet", "potential":    0.0, "physicalSurface": 1}
+    {"name": "Ionization Chamber", "type": "dirichlet", "potential":    0.0, "physicalSurface": 1},
+    {"name": "Infinite",           "type": "dirichlet", "potential": -600.0, "physicalSurface": 2}
   ],
   "inject": [
-    {"name": "Cathode Electron",    "species": "Electron", "flow":  1.0e-4, "units": "A", "v":  87000.0, "T": [2500.0, 2500.0, 2500.0],
+    {"name": "Cathode Electron",    "species": "Electron", "flow":  1.0e-4, "units": "A", "v":  87000.0, "T": [30000.0, 30000.0, 30000.0],
                                     "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [-1, 0, 0], "physicalSurface": 2}
   ],
   "reference": {

runs/ALPHIE_Grid/mesh.geo

@@ -1,11 +1,11 @@
-zg1  = 0.0025;
+zg1  = 0.0020;
-tg1  = 0.0004;
+tg1  = 0.0003;
 rg1  = 0.0005;
-dg   = 0.0025;
+dg   = 0.0020;
 zg2  = zg1 + tg1 + dg;
 tg2  = tg1;
 rg2  = rg1;
-zEnd = 0.0042;
+zEnd = 0.0050;
 Lz   = zg2 + tg2 + zEnd;
 Lr   = rg1 + 0.0001;
 

src/modules/mesh/moduleMeshBoundary.f90

@@ -110,31 +110,34 @@ MODULE moduleMeshBoundary
       USE moduleMesh
       USE moduleRefParam
       USE moduleRandom
+      USE moduleMath
       IMPLICIT NONE
 
       CLASS(meshEdge), INTENT(inout):: edge
       CLASS(particle), INTENT(inout):: part
       REAL(8):: vRel, eRel, mRel !relative velocity, energy and mass
-      REAL(8):: ionizationRate
+      REAL(8):: nIonizations !Number of ionizations based on eRel
-      INTEGER:: ionizationPair, p
+      REAL(8):: pIonization !Probability of ionization of each event
+      INTEGER:: p
       REAL(8):: v0(1:3) !random velocity of neutral
       TYPE(particle), POINTER:: newElectron
       TYPE(particle), POINTER:: newIon
 
       SELECT TYPE(bound => edge%boundary%bTypes(part%species%n)%obj)
       TYPE IS(boundaryIonization)
-        mRel = (bound%m0*part%species%m)*(bound%m0+part%species%m)
+        mRel = reducedMass(bound%m0, part%species%m)
         vRel = SUM(DABS(part%v-bound%v0))
         eRel = mRel*vRel**2*5.D-1
 
-        IF (eRel > bound%eThreshold) THEN
+        !Maximum number of possible ionizations based on relative energy
-          ionizationRate = part%weight*bound%n0*bound%crossSection%get(eRel)*vRel
+        nIonizations = eRel/bound%eThreshold
 
-          !Rounds the number of particles up
+        DO p = 1, FLOOR(nIonizations)
-          ionizationPair = NINT(ionizationRate*bound%effectiveTime/bound%species%weight)
+          !Get probability of ionization
+          pIonization = 1.D0 - DEXP(-bound%n0*bound%crossSection%get(eRel)*vRel*bound%effectiveTime/nIonizations)
 
-          !Create the new pair of particles
+          !If a random number is below the probability of ionization, create new pair of ion-electron
-          DO p = 1, ionizationPair
+          IF (random() < pIonization) THEN
             !Assign random velocity to the neutral
             v0(1) = bound%v0(1) + bound%vTh*randomMaxwellian()
             v0(2) = bound%v0(2) + bound%vTh*randomMaxwellian()
@@ -159,7 +162,7 @@ MODULE moduleMeshBoundary
             newElectron%xi = mesh%vols(part%vol)%obj%phy2log(newElectron%r)
             newIon%xi      = newElectron%xi
 
-            newElectron%weight = bound%species%weight
+            newElectron%weight = part%weight
             newIon%weight      = newElectron%weight
 
             newElectron%n_in = .TRUE.
@@ -171,10 +174,14 @@ MODULE moduleMeshBoundary
             CALL partSurfaces%add(newIon)
             CALL OMP_UNSET_LOCK(lockSurfaces)
 
-          END DO
+            !Electron loses energy due to ionization
+            eRel = eRel - bound%eThreshold
+            vRel = 2.D0*DSQRT(eRel)/mRel
 
           END IF
 
+        END DO
+
       END SELECT
 
       !Removes ionizing electron regardless the number of pair created

src/modules/moduleInject.f90

@@ -252,12 +252,22 @@ MODULE moduleInject
       CLASS(injectGeneric), INTENT(in):: self
       INTEGER:: randomX
       INTEGER, SAVE:: nMin, nMax !Min and Max index in partInj array
+      INTEGER:: i
       INTEGER:: n, sp
       CLASS(meshEdge), POINTER:: randomEdge
 
       !Insert particles
       !$OMP SINGLE
-      nMin = SUM(inject(1:(self%id-1))%nParticles) + 1
+      nMin = 0
+      DO i = 1, self%id -1
+        IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN 
+          nMin = nMin + inject(i)%nParticles
+
+        END IF
+
+      END DO
+
+      nMin = nMin + 1
       nMax = nMin + self%nParticles - 1
       !Assign weight to particle.
       partInj(nMin:nMax)%weight = self%species%weight