Merge branch 'feature/injectionNormal' into 'development'

Injection of particles from surfaces direction

See merge request JorgeGonz/fpakc!36

doc/user-manual/fpakc_UserManual.tex

@@ -608,7 +608,8 @@ make
           \item \textbf{n}: Real.
                             Array dimension $3$.
                             Direction of injection.
-                            This vector is normalized to $1$.
+                            If no value is provided, the normal of the edge is used as the direction of injection.
+                            This vector is normalized to $1$ by the code.
           \item \textbf{velDist}: Character.
                                   Array dimension $3$.
                                   Type of distribution function used to obtain injected particle velocity:

src/modules/init/moduleInput.f90

@@ -1140,6 +1140,7 @@ MODULE moduleInput
       CHARACTER(:), ALLOCATABLE:: name
       REAL(8):: v
       REAL(8), ALLOCATABLE:: T(:), normal(:)
+      LOGICAL:: fixDirection
       REAL(8):: flow
       CHARACTER(:), ALLOCATABLE:: units
       INTEGER:: physicalSurface
@@ -1160,11 +1161,18 @@ 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
+          fixDirection = .TRUE.
+        ELSE
+          fixDirection = .FALSE.
+          ALLOCATE(normal(1:3))
+          normal = 0.D0
+        END IF
         CALL config%get(object // '.flow',            flow,            found)
         CALL config%get(object // '.units',           units,           found)
         CALL config%get(object // '.physicalSurface', physicalSurface, found)
 
-        CALL inject(i)%init(i, v, normal, T, flow, units, sp, physicalSurface)
+        CALL inject(i)%init(i, v, normal, fixDirection, T, flow, units, sp, physicalSurface)
 
         CALL readVelDistr(config, inject(i), object)
         
@@ -1240,13 +1248,12 @@ MODULE moduleInput
 
         SELECT CASE(fvType)
         CASE ("Maxwellian")
-          v = inj%vMod*inj%n(i)
           T = inj%T(i)
-          CALL initVelDistMaxwellian(inj%v(i)%obj, v, t, m)
+          CALL initVelDistMaxwellian(inj%v(i)%obj, t, m)
 
         CASE ("Delta")
           v = inj%vMod*inj%n(i)
-          CALL initVelDistDelta(inj%v(i)%obj, v)
+          CALL initVelDistDelta(inj%v(i)%obj)
 
         CASE DEFAULT
           CALL criticalError("No velocity distribution type " // fvType // " defined", 'readVelDistr')

src/modules/moduleInject.f90

@@ -29,7 +29,6 @@ MODULE moduleInject
 
   !Maxwellian distribution function
   TYPE, EXTENDS(velDistGeneric):: velDistMaxwellian
-    REAL(8):: v !Velocity
     REAL(8):: vTh !Thermal Velocity
     CONTAINS
       PROCEDURE, PASS:: randomVel => randomVelMaxwellian
@@ -38,7 +37,6 @@ MODULE moduleInject
 
   !Dirac's delta distribution function
   TYPE, EXTENDS(velDistGeneric):: velDistDelta
-    REAL(8):: v !Velocity
     CONTAINS
       PROCEDURE, PASS:: randomVel => randomVelDelta
 
@@ -51,6 +49,7 @@ MODULE moduleInject
     REAL(8):: vMod !Velocity (module)
     REAL(8):: T(1:3) !Temperature
     REAL(8):: n(1:3) !Direction of injection
+    LOGICAL:: fixDirection !The injection of particles has a fix direction defined by n
     INTEGER:: nParticles !Number of particles to introduce each time step
     CLASS(speciesGeneric), POINTER:: species !Species of injection
     INTEGER:: nEdges
@@ -69,7 +68,7 @@ MODULE moduleInject
 
   CONTAINS
     !Initialize an injection of particles
-    SUBROUTINE initInject(self, i, v, n, T, flow, units, sp, physicalSurface)
+    SUBROUTINE initInject(self, i, v, n, fixDirection, T, flow, units, sp, physicalSurface)
       USE moduleMesh
       USE moduleRefParam
       USE moduleConstParam
@@ -81,6 +80,7 @@ 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,7 +91,12 @@ 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)
@@ -197,23 +202,22 @@ MODULE moduleInject
 
     END SUBROUTINE doInjects
 
-    SUBROUTINE initVelDistMaxwellian(velDist, v, T, m)
+    SUBROUTINE initVelDistMaxwellian(velDist, T, m)
       IMPLICIT NONE
 
       CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
-      REAL(8), INTENT(in):: v, T, m
+      REAL(8), INTENT(in):: T, m
 
-      velDist = velDistMaxwellian(v = v, vTh = DSQRT(T/m))
+      velDist = velDistMaxwellian(vTh = DSQRT(T/m))
 
     END SUBROUTINE initVelDistMaxwellian
 
-    SUBROUTINE initVelDistDelta(velDist, v)
+    SUBROUTINE initVelDistDelta(velDist)
       IMPLICIT NONE
 
       CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
-      REAL(8), INTENT(in):: v
 
-      velDist = velDistDelta(v = v)
+      velDist = velDistDelta()
 
     END SUBROUTINE initVelDistDelta
 
@@ -226,7 +230,7 @@ MODULE moduleInject
       REAL(8):: v
       v = 0.D0
 
-      v = self%v + self%vTh*randomMaxwellian()
+      v = self%vTh*randomMaxwellian()
 
     END FUNCTION randomVelMaxwellian
 
@@ -237,7 +241,7 @@ MODULE moduleInject
       CLASS(velDistDelta), INTENT(in):: self
       REAL(8):: v
 
-      v = self%v
+      v = 0.D0
 
     END FUNCTION randomVelDelta
 
@@ -256,6 +260,7 @@ MODULE moduleInject
       INTEGER:: i
       INTEGER:: n, sp
       CLASS(meshEdge), POINTER:: randomEdge
+      REAL(8):: direction(1:3)
 
       !Insert particles
       !$OMP SINGLE
@@ -300,10 +305,21 @@ MODULE moduleInject
         !Assign particle type
         partInj(n)%species => self%species
 
-        partInj(n)%v = (/ self%v(1)%obj%randomVel(), &
+        IF (self%fixDirection) THEN
+          direction = self%n
+
+        ELSE
+          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() /)
 
+        print *, direction
+        print*, partInj(n)%v
+
         !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