Merge branch 'feature/boundary' into 'development'

Ionization and symmetry axis fixin

See merge request JorgeGonz/fpakc!11

src/modules/mesh/moduleMeshBoundary.f90

@@ -111,6 +111,7 @@ MODULE moduleMeshBoundary
       USE moduleSpecies
       USE moduleMesh
       USE moduleRefParam
+      USE moduleRandom
       IMPLICIT NONE
 
       CLASS(meshEdge), INTENT(inout):: edge
@@ -118,6 +119,7 @@ MODULE moduleMeshBoundary
       REAL(8):: vRel, eRel, mRel !relative velocity, energy and mass
       REAL(8):: ionizationRate
       INTEGER:: ionizationPair, p
+      REAL(8):: v0(1:3) !random velocity of neutral
       TYPE(particle), POINTER:: newElectron
       TYPE(particle), POINTER:: newIon
 
@@ -128,21 +130,26 @@ MODULE moduleMeshBoundary
         eRel = mRel*vRel**2*5.D-1
 
         IF (eRel > bound%eThreshold) THEN
-          !TODO: Check units
+          ionizationRate = part%weight*bound%n0*bound%crossSection%get(eRel)*vRel
-          ionizationRate = bound%n0*bound%crossSection%get(eRel)
 
-          ionizationPair = INT(ionizationRate*tauMin*ti_ref/species(bound%sp)%obj%weight)
+          !Rounds the number of particles up
+          ionizationPair = NINT(ionizationRate*bound%effectiveTime/species(bound%sp)%obj%weight)
 
           !Create the new pair of particles
           DO p = 1, ionizationPair
             ALLOCATE(newElectron)
             ALLOCATE(newIon)
 
+            !Assign random velocity to the neutral
+            v0(1) = bound%v0(1) + bound%vTh*randomMaxwellian()
+            v0(2) = bound%v0(2) + bound%vTh*randomMaxwellian()
+            v0(3) = bound%v0(3) + bound%vTh*randomMaxwellian()
+
             newElectron%sp = part%sp
             newIon%sp      = bound%sp
 
-            newElectron%v = 10.D0*bound%v0 + (1.D0 + bound%deltaV/NORM2(bound%v0))
+            newElectron%v = v0 + (1.D0 + bound%deltaV*v0/NORM2(v0))
-            newIon%v      = bound%v0
+            newIon%v      = v0
 
             newElectron%r = edge%randPos()
             newIon%r      = newElectron%r
@@ -174,16 +181,18 @@ MODULE moduleMeshBoundary
 
           END DO
 
-          !Removes ionizing electron
-          part%n_in = .FALSE.
-
         END IF
 
       END SELECT
 
+      !Removes ionizing electron regardless the number of pair created
+      part%n_in = .FALSE.
+
     END SUBROUTINE ionization
 
-    !Symmetry axis. Dummy function
+    !Symmetry axis. Reflects particles.
+    !Although this function should never be called, it is set as a reflective boundary
+    !to properly deal with possible particles reaching a corner and selecting this boundary.
     SUBROUTINE symmetryAxis(edge, part)
       USE moduleSpecies
       IMPLICIT NONE
@@ -191,6 +200,8 @@ MODULE moduleMeshBoundary
       CLASS(meshEdge), INTENT(inout):: edge
       CLASS(particle), INTENT(inout):: part
 
+      CALL reflection(edge, part)
+
     END SUBROUTINE symmetryAxis
 
     !Points the boundary function to specific type

src/modules/moduleBoundary.f90

@@ -35,9 +35,10 @@ MODULE moduleBoundary
 
   !Ionization boundary
   TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryIonization
-    REAL(8):: m0, n0, v0(1:3), T0 !Properties of background neutrals.
+    REAL(8):: m0, n0, v0(1:3), vTh !Properties of background neutrals.
     INTEGER:: sp !Ion species
     TYPE(table1D):: crossSection
+    REAL(8):: effectiveTime
     REAL(8):: eThreshold
     REAL(8):: deltaV
 
@@ -99,7 +100,7 @@ MODULE moduleBoundary
 
     END SUBROUTINE initWallTemperature
 
-    SUBROUTINE initIonization(boundary, m0, n0, v0, T0, speciesID, crossSection, eThreshold)
+    SUBROUTINE initIonization(boundary, me, m0, n0, v0, T0, speciesID, effTime, crossSection, eThreshold)
       USE moduleRefParam
       USE moduleSpecies
       USE moduleCaseParam
@@ -107,8 +108,10 @@ MODULE moduleBoundary
       IMPLICIT NONE
 
       CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary
-      REAL(8), INTENT(in):: m0, n0, v0(1:3), T0
+      REAL(8), INTENT(in):: me !Electron mass
+      REAL(8), INTENT(in):: m0, n0, v0(1:3), T0 !Neutral properties
       INTEGER:: speciesID
+      REAL(8):: effTime
       CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSection
       REAL(8), INTENT(in):: eThreshold
 
@@ -117,13 +120,15 @@ MODULE moduleBoundary
       SELECT TYPE(boundary)
       TYPE IS(boundaryIonization)
         boundary%m0 = m0 / m_ref
-        boundary%n0 = n0
+        boundary%n0 = n0 * Vol_ref
         boundary%v0 = v0 / v_ref
-        boundary%T0 = T0 / T_ref
+        boundary%vTh = DSQRT(kb*T0/m0)/v_ref
         boundary%sp = speciesID
+        boundary%effectiveTime = effTime / ti_ref
         CALL boundary%crossSection%init(crossSection)
         CALL boundary%crossSection%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2)
         boundary%eThreshold = eThreshold*eV2J/(m_ref*v_ref**2)
+        boundary%deltaV = DSQRT(boundary%eThreshold/me)
 
       END SELECT
 

src/modules/moduleInput.f90

@@ -600,9 +600,9 @@ MODULE moduleInput
       CHARACTER(:), ALLOCATABLE:: object, bType
       REAL(8):: Tw, cw !Wall temperature and specific heat
       !Neutral Properties
-      REAL(8):: m0, n0
+      REAL(8):: m0, n0, T0
       REAL(8), DIMENSION(:), ALLOCATABLE:: v0
-      REAL(8):: T0
+      REAL(8):: effTime
       REAL(8):: eThreshold !Energy threshold
       INTEGER:: speciesID
       CHARACTER(:), ALLOCATABLE:: speciesName, crossSection
@@ -647,12 +647,15 @@ MODULE moduleInput
             IF (.NOT. found) CALL criticalError("missing parameter 'velocity' for neutrals in ionization", 'readBoundary')
             CALL config%get(object // '.neutral.temperature', T0, found)
             IF (.NOT. found) CALL criticalError("missing parameter 'temperature' for neutrals in ionization", 'readBoundary')
+            CALL config%get(object // '.effectiveTime', effTime, found)
+            IF (.NOT. found) CALL criticalError("missing parameter 'effectiveTime' for neutrals in ionization", 'readBoundary')
             CALL config%get(object // '.energyThreshold', eThreshold, found)
             IF (.NOT. found) CALL criticalError("missing parameter 'eThreshold' in ionization", 'readBoundary')
             CALL config%get(object // '.crossSection', crossSection, found)
             IF (.NOT. found) CALL criticalError("missing parameter 'crossSection' for neutrals in ionization", 'readBoundary')
 
-            CALL initIonization(boundary(i)%bTypes(s)%obj, m0, n0, v0, T0, speciesID, crossSection, eThreshold)
+            CALL initIonization(boundary(i)%bTypes(s)%obj, species(s)%obj%m, m0, n0, v0, T0, &
+                                                           speciesID, effTime, crossSection, eThreshold)
 
           CASE('wallTemperature')
             CALL config%get(object // '.temperature', Tw, found)