Merge branch 'issue/injection' into 'feature/CoulombLinear'

Issue with injection of particles

See merge request JorgeGonz/fpakc!45

doc/user-manual/fpakc_UserManual.tex

@@ -619,6 +619,8 @@ make
                                   Type of distribution function used to obtain injected particle velocity:
                                   \begin{itemize}
                                     \item \textbf{Maxwellian}: Maxwellian distribution of temperature \textbf{T} and mean \textbf{v} times the value of \textbf{n} in the specified direction.
+                                    \item \textbf{Half-Maxwellian}: Half-Maxwellian distribution of temperature \textbf{T} and mean \textbf{v} times the value of \textbf{n} in the specified direction.
+                                                                    Only takes into account the positive part of the half-Maxwellian.
                                     \item \textbf{Delta}: Dirac's delta distribution function. All particles are injected with velocity \textbf{v} times the value of \textbf{n} in the specified direction.
                                   \end{itemize}
           \item \textbf{T}: Real.

src/modules/init/moduleInput.f90

@@ -1222,7 +1222,6 @@ MODULE moduleInput
       CHARACTER(:), ALLOCATABLE:: name
       REAL(8):: v
       REAL(8), ALLOCATABLE:: T(:), normal(:)
-      LOGICAL:: fixDirection
       REAL(8):: flow
       CHARACTER(:), ALLOCATABLE:: units
       INTEGER:: physicalSurface
@@ -1243,10 +1242,7 @@ MODULE moduleInput
         CALL config%get(object // '.v',               v,               found)
         CALL config%get(object // '.T',               T,               found)
         CALL config%get(object // '.n',               normal,          found)
-        IF (found) THEN
+        IF (.NOT. found) THEN
-          fixDirection = .TRUE.
-        ELSE
-          fixDirection = .FALSE.
           ALLOCATE(normal(1:3))
           normal = 0.D0
         END IF
@@ -1254,7 +1250,7 @@ MODULE moduleInput
         CALL config%get(object // '.units',           units,           found)
         CALL config%get(object // '.physicalSurface', physicalSurface, found)
 
-        CALL inject(i)%init(i, v, normal, fixDirection, T, flow, units, sp, physicalSurface)
+        CALL inject(i)%init(i, v, normal, T, flow, units, sp, physicalSurface)
 
         CALL readVelDistr(config, inject(i), object)
         
@@ -1333,6 +1329,10 @@ MODULE moduleInput
           T = inj%T(i)
           CALL initVelDistMaxwellian(inj%v(i)%obj, t, m)
 
+        CASE ("Half-Maxwellian")
+          T = inj%T(i)
+          CALL initVelDistHalfMaxwellian(inj%v(i)%obj, t, m)
+
         CASE ("Delta")
           v = inj%vMod*inj%n(i)
           CALL initVelDistDelta(inj%v(i)%obj)

src/modules/moduleInject.f90

@@ -35,6 +35,13 @@ MODULE moduleInject
   
   END TYPE velDistMaxwellian
 
+  TYPE, EXTENDS(velDistGeneric):: velDistHalfMaxwellian
+    REAL(8):: vTh !Thermal Velocity
+    CONTAINS
+      PROCEDURE, PASS:: randomVel => randomVelHalfMaxwellian
+  
+  END TYPE velDistHalfMaxwellian
+
   !Dirac's delta distribution function
   TYPE, EXTENDS(velDistGeneric):: velDistDelta
     CONTAINS
@@ -68,7 +75,7 @@ MODULE moduleInject
 
   CONTAINS
     !Initialize an injection of particles
-    SUBROUTINE initInject(self, i, v, n, fixDirection, T, flow, units, sp, physicalSurface)
+    SUBROUTINE initInject(self, i, v, n, T, flow, units, sp, physicalSurface)
       USE moduleMesh
       USE moduleRefParam
       USE moduleConstParam
@@ -80,7 +87,6 @@ MODULE moduleInject
       CLASS(injectGeneric), INTENT(inout):: self
       INTEGER, INTENT(in):: i
       REAL(8), INTENT(in):: v, n(1:3), T(1:3)
-      LOGICAL, INTENT(in):: fixDirection
       INTEGER, INTENT(in):: sp, physicalSurface
       REAL(8):: tauInject
       REAL(8), INTENT(in):: flow
@@ -91,12 +97,7 @@ MODULE moduleInject
 
       self%id         =  i
       self%vMod       =  v / v_ref
-      IF (.NOT. fixDirection) THEN
       self%n        =  n / NORM2(n)
-      ELSE
-        self%n        = 0.D0
-      END IF
-      self%fixDirection = fixDirection
       self%T          =  T / T_ref
       self%species    => species(sp)%obj
       tauInject = tau(self%species%n)
@@ -212,6 +213,16 @@ MODULE moduleInject
 
     END SUBROUTINE initVelDistMaxwellian
 
+    SUBROUTINE initVelDistHalfMaxwellian(velDist, T, m)
+      IMPLICIT NONE
+
+      CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
+      REAL(8), INTENT(in):: T, m
+
+      velDist = velDistHalfMaxwellian(vTh = DSQRT(T/m))
+
+    END SUBROUTINE initVelDistHalfMaxwellian
+
     SUBROUTINE initVelDistDelta(velDist)
       IMPLICIT NONE
 
@@ -234,6 +245,22 @@ MODULE moduleInject
 
     END FUNCTION randomVelMaxwellian
 
+    !Random velocity from Half Maxwellian distribution
+    FUNCTION randomVelHalfMaxwellian(self) RESULT (v)
+      USE moduleRandom
+      IMPLICIT NONE
+
+      CLASS(velDistHalfMaxwellian), INTENT(in):: self
+      REAL(8):: v
+      v = 0.D0
+
+      DO WHILE (v <= 0.D0)
+        v = self%vTh*randomMaxwellian()
+
+      END DO
+
+    END FUNCTION randomVelHalfMaxwellian
+
     !Random velocity from Dirac's delta distribution
     PURE FUNCTION randomVelDelta(self) RESULT(v)
       IMPLICIT NONE
@@ -305,18 +332,20 @@ MODULE moduleInject
         !Assign particle type
         partInj(n)%species => self%species
 
-        IF (self%fixDirection) THEN
         direction = self%n
 
-        ELSE
+        partInj(n)%v = 0.D0
-          direction = randomEdge%normal
-
-        END IF
 
         partInj(n)%v = self%vMod*direction +  (/ self%v(1)%obj%randomVel(), &
                                                  self%v(2)%obj%randomVel(), &
                                                  self%v(3)%obj%randomVel() /)
 
+        !If velocity is not in the right direction, invert it
+        IF (DOT_PRODUCT(direction, partInj(n)%v) < 0.D0) THEN
+          partInj(n)%v = - partInj(n)%v
+
+        END IF
+
         !Obtain natural coordinates of particle in cell
         partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r)
         !Push new particle with the minimum time step