From 5368ff2bf385aab033a58a5c967f8f55e8b841c6 Mon Sep 17 00:00:00 2001
From: Jorge Gonzalez <jorge.gonzalez@upm.es>
Date: Wed, 21 Oct 2020 10:13:43 +0200
Subject: [PATCH] Fixing a bug in injection of particles that was creating
 ghost particles assumed to be inside the domain but without cell assigned.

Now, particles are assumed to be outside the domain (n_in = .FALSE.)
  until findCell assign them a cell.
---
 src/modules/moduleInject.f95  | 20 ++++++++++----------
 src/modules/moduleMeshCyl.f95 |  1 +
 src/modules/moduleSolver.f95  |  1 +
 3 files changed, 12 insertions(+), 10 deletions(-)

diff --git a/src/modules/moduleInject.f95 b/src/modules/moduleInject.f95
index fa23b02..63f53fb 100644
--- a/src/modules/moduleInject.f95
+++ b/src/modules/moduleInject.f95
@@ -4,8 +4,10 @@ MODULE moduleInject
   TYPE:: injectGeneric
     INTEGER:: id
     CHARACTER(:), ALLOCATABLE:: name
-    REAL(8):: v !Velocity (module)
+    REAL(8):: vMod !Velocity (module)
     REAL(8):: T(1:3) !Temperature
+    REAL(8):: v(1:3) !Velocity(vector)
+    REAL(8):: vTh(1:3) !Thermal Velocity
     REAL(8):: n(1:3) !Direction of injection
     INTEGER:: nParticles !Number of particles to introduce each time step
     INTEGER:: pt !Species of injection
@@ -52,12 +54,15 @@ MODULE moduleInject
       INTEGER:: phSurface(1:mesh%numEdges)
 
       self%id         = i
-      self%v          = v/v_ref
+      self%vMod       = v/v_ref
       self%n          = n
       self%T          = T/T_ref
       self%nParticles = INT(flow*sccm2atomPerS*tau*ti_ref/species(pt)%obj%weight)
       self%pt         = pt
 
+      self%v   = self%v * self%n
+      self%vTh = DSQRT(self%T/species(self%pt)%obj%m)
+
       !Gets the edge elements from which particles are injected
       !TODO: Improve this A LOT
       DO e = 1, mesh%numEdges
@@ -116,7 +121,6 @@ MODULE moduleInject
       INTEGER:: j
       INTEGER:: randomEdge
       REAL(8):: randomPos
-      REAL(8):: vVec(1:3), vTh(1:3)
       !Edge nodes
       INTEGER:: n1 = 0, n2 = 0
       !Edge nodes coordinates
@@ -124,9 +128,6 @@ MODULE moduleInject
       INTEGER:: nMin, nMax !Min and Max index in partInj array
       INTEGER:: n
 
-      vVec = self%v * self%n
-      vTh  = DSQRT(self%T/species(self%pt)%obj%m)
-
       !Insert particles
       !TODO: Adjust for multiple injections
       nMin = 1
@@ -135,7 +136,6 @@ MODULE moduleInject
       partInj(nMin:nMax)%pt = self%pt
       !Assign weight to particle.
       partInj(nMin:nMax)%weight = species(self%pt)%obj%weight
-      partInj(nMin:nMax)%n_in = .TRUE.
       DO n = nMin, nMax
         !Select edge randomly from which inject particle
         randomX = RAND()*self%sumWeight
@@ -176,9 +176,9 @@ MODULE moduleInject
 
         END IF
 
-        partInj(n)%v = (/ vBC(vVec(1), vTh(1)), &
-                           vBC(vVec(2), vTh(2)), &
-                           vBC(vVec(3), vTh(3)) /)
+        partInj(n)%v = (/ vBC(self%v(1), self%vTh(1)), &
+                          vBC(self%v(2), self%vTh(2)), &
+                          vBC(self%v(3), self%vTh(3)) /)
 
         !Random position in edge
         !TODO: Use edge coordinates and transformations for this process
diff --git a/src/modules/moduleMeshCyl.f95 b/src/modules/moduleMeshCyl.f95
index 2696b05..17689c4 100644
--- a/src/modules/moduleMeshCyl.f95
+++ b/src/modules/moduleMeshCyl.f95
@@ -451,6 +451,7 @@ MODULE moduleMeshCyl
         END IF
         part%e_p  = self%n
         part%xLog = xLog
+        part%n_in = .TRUE.
         !Assign particle to listPart_in
         CALL OMP_SET_LOCK(self%lock)
         CALL self%listPart_in%add(part)
diff --git a/src/modules/moduleSolver.f95 b/src/modules/moduleSolver.f95
index 7566bb1..1676425 100644
--- a/src/modules/moduleSolver.f95
+++ b/src/modules/moduleSolver.f95
@@ -53,6 +53,7 @@ MODULE moduleSolver
       part_temp%v(3) = -sin_alpha*v_p_oh_star(2)+cos_alpha*v_p_oh_star(3)
       part_temp%pt  = part%pt
       part_temp%e_p = part%e_p
+      part_temp%n_in = .FALSE. !Assume particle is outside until cell is found
       !Copy temporal particle to particle
       part=part_temp