Fixed an issue in which the logical coordinates of the particle (xi)

were not being initialized and was causing particles to get stuck in
  the iteration loop.

runs/ALPHIE_Grid/inputDiffTau.json

@@ -58,17 +58,18 @@
   "reference": {
     "density": 1.0e16,
     "mass": 9.109e-31,
-    "temperature": 2500.0
+    "temperature": 2500.0,
+    "radius": 1.88e-10
   },
   "case": {
     "tau":  [1.0e-9, 1.0e-11],
-    "time": 2.0e-6,
+    "time": 1.0e-6,
     "pusher": ["2DCylCharged", "2DCylCharged"],
     "EMSolver": "Electrostatic"
   },
   "parallel": {
     "OpenMP":{
-      "nThreads": 8
+      "nThreads": 24
     }
   }
 }

runs/ALPHIE_Grid/inputSameTau.json

@@ -14,8 +14,8 @@
     "meshFile": "mesh.msh"
   },
   "species": [
-    {"name": "Argon+",   "type": "charged", "mass": 6.633e-26, "charge": 1.0, "weight": 1.0e1},
+    {"name": "Argon+",   "type": "charged", "mass": 6.633e-26, "charge": 1.0, "weight": 1.0e2},
-    {"name": "Electron", "type": "charged", "mass": 9.109e-31, "charge":-1.0, "weight": 1.0e1}
+    {"name": "Electron", "type": "charged", "mass": 9.109e-31, "charge":-1.0, "weight": 1.0e2}
   ],
   "boundary": [ 
     {"name": "Ionization Chanber", "physicalSurface": 1, "bTypes": [
@@ -58,7 +58,8 @@
   "reference": {
     "density": 1.0e16,
     "mass": 9.109e-31,
-    "temperature": 2500.0
+    "temperature": 2500.0,
+    "radius": 1.88e-10
   },
   "case": {
     "tau":  [1.0e-11, 1.0e-11],
@@ -68,7 +69,7 @@
   },
   "parallel": {
     "OpenMP":{
-      "nThreads": 8
+      "nThreads": 4
     }
   }
 }

runs/Argon_Expansion/CX_case.json

@@ -3,7 +3,7 @@
     "path": "./runs/Argon_Expansion/",
 	"triggerOutput": 10,
 	"cpuTime": false,
-	"numColl": true,
+	"numColl": false,
 	"folder": "CX_case"
   },
   "geometry": {

runs/Argon_Expansion/elastic_case.json

@@ -4,7 +4,7 @@
 	"triggerOutput": 10,
 	"cpuTime": false,
 	"numColl": false,
-	"folder": "base_case"
+	"folder": "Elastic_case"
   },
   "geometry": {
     "type":     "2DCyl",

src/modules/mesh/2DCart/moduleMesh2DCart.f90

@@ -117,7 +117,6 @@ MODULE moduleMesh2DCart
     !Connectivity to adjacent elements
     CLASS(*), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL(), e4 => NULL()
     REAL(8):: arNodes(1:4) = 0.D0
-
     CONTAINS
       PROCEDURE, PASS::   init        => initVolQuad2DCart
       PROCEDURE, PASS::   randPos     => randPosVolQuad

src/modules/mesh/2DCyl/moduleMesh2DCyl.f90

@@ -274,28 +274,6 @@ MODULE moduleMesh2DCyl
 
     END SUBROUTINE initEdge2DCyl
 
-    !Random position in quadrilateral volume
-    FUNCTION randPosVolQuad(self) RESULT(r)
-      USE moduleRandom
-      IMPLICIT NONE
-
-      CLASS(meshVol2DCylQuad), INTENT(in):: self
-      REAL(8):: r(1:3)
-      REAL(8):: xii(1:3)
-      REAL(8), ALLOCATABLE:: fPsi(:)
-
-      xii(1) = random(-1.D0, 1.D0)
-      xii(2) = random(-1.D0, 1.D0)
-      xii(3) = 0.D0
-
-      fPsi = self%fPsi(xii)
-
-      r(1) = DOT_PRODUCT(fPsi, self%z)
-      r(2) = DOT_PRODUCT(fPsi, self%r)
-      r(3) = 0.D0
-
-    END FUNCTION randposVolQuad
-
     !Get nodes from edge
     PURE FUNCTION getNodes2DCyl(self) RESULT(n)
       IMPLICIT NONE
@@ -448,6 +426,7 @@ MODULE moduleMesh2DCyl
 
     END FUNCTION dPsiQuadXi2
 
+    !Partial derivative in global coordinates
     PURE SUBROUTINE partialDerQuad(self, dPsi, dz, dr)
       IMPLICIT NONE
 
@@ -462,6 +441,28 @@ MODULE moduleMesh2DCyl
 
     END SUBROUTINE partialDerQuad
 
+    !Random position in quadrilateral volume
+    FUNCTION randPosVolQuad(self) RESULT(r)
+      USE moduleRandom
+      IMPLICIT NONE
+
+      CLASS(meshVol2DCylQuad), INTENT(in):: self
+      REAL(8):: r(1:3)
+      REAL(8):: xii(1:3)
+      REAL(8), ALLOCATABLE:: fPsi(:)
+
+      xii(1) = random(-1.D0, 1.D0)
+      xii(2) = random(-1.D0, 1.D0)
+      xii(3) = 0.D0
+
+      fPsi = self%fPsi(xii)
+
+      r(1) = DOT_PRODUCT(fPsi, self%z)
+      r(2) = DOT_PRODUCT(fPsi, self%r)
+      r(3) = 0.D0
+
+    END FUNCTION randposVolQuad
+
     !Computes element local stiffness matrix
     PURE FUNCTION elemKQuad(self) RESULT(ke)
       USE moduleConstParam, ONLY: PI2

src/modules/mesh/moduleMesh.f90

@@ -359,8 +359,9 @@ MODULE moduleMesh
           CALL nextElement%findCell(part, self)
 
         CLASS IS (meshEdge)
-          !Particle encountered an edge, apply boundary
+          !Particle encountered a surface, apply boundary
           CALL nextElement%fBoundary(part%sp)%apply(nextElement,part)
+
           !If particle is still inside the domain, call findCell
           IF (part%n_in) THEN 
             IF(PRESENT(oldCell)) THEN
@@ -474,7 +475,7 @@ MODULE moduleMesh
         WRITE(60, "(A)") '$EndMeshFormat'
         WRITE(60, "(A)") '$NodeData'
         WRITE(60, "(A)") '1'
-        WRITE(60, "(A)") '"Density' // species(i)%obj%name // ' (m^-3)"'
+        WRITE(60, "(A)") '"Density ' // species(i)%obj%name // ' (m^-3)"'
         WRITE(60, *) 1
         WRITE(60, *) time
         WRITE(60, *) 3

src/modules/moduleInject.f90

@@ -242,6 +242,8 @@ MODULE moduleInject
                           self%v(2)%obj%randomVel(), &
                           self%v(3)%obj%randomVel() /)
 
+        !Obtain natural coordinates of particle in cell
+        partInj(n)%xi = mesh%vols(partInj(n)%vol)%obj%phy2log(partInj(n)%r)
         !Push new particle with the minimum time step
         CALL solver%pusher(self%sp)%pushParticle(partInj(n), tauMin)
         !Assign cell to new particle