First implementation of ionization boundary. Still some work to do.

src/modules/mesh/1DCart/moduleMesh1DCart.f90

@@ -7,6 +7,9 @@ MODULE moduleMesh1DCart
   USE moduleMeshBoundary
   IMPLICIT NONE
 
+  REAL(8), PARAMETER:: corSeg(1:3) = (/ -DSQRT(3.D0/5.D0), 0.D0,       DSQRT(3.D0/5.D0) /)
+  REAL(8), PARAMETER:: wSeg(1:3)  = (/        5.D0/9.D0 , 8.D0/9.D0,        5.D0/9.D0  /)
+
   TYPE, PUBLIC, EXTENDS(meshNode):: meshNode1DCart
     !Element coordinates
     REAL(8):: x = 0.D0
@@ -152,23 +155,7 @@ MODULE moduleMesh1DCart
       ALLOCATE(self%fboundary(1:nSpecies))
       !Assign functions to boundary
       DO s = 1, nSpecies
-        SELECT TYPE(obj => self%boundary%bTypes(s)%obj)
+        CALL pointBoundaryFunction(self, s)
-        TYPE IS(boundaryAbsorption)
-          self%fBoundary(s)%apply => absorption
-
-        TYPE IS(boundaryReflection)
-          self%fBoundary(s)%apply => reflection
-
-        TYPE IS(boundaryTransparent)
-          self%fBoundary(s)%apply => transparent
-
-        TYPE IS(boundaryWallTemperature)
-          self%fBoundary(s)%apply => wallTemperature
-
-        CLASS DEFAULT
-          CALL criticalError("Boundary type not defined in this geometry", 'initEdge1DCart')
-
-        END SELECT
 
       END DO
       
@@ -322,22 +309,29 @@ MODULE moduleMesh1DCart
     END SUBROUTINE partialDerSegm
 
     !Computes local stiffness matrix
-    PURE FUNCTION elemKSegm(self) RESULT(ke)
+    FUNCTION elemKSegm(self) RESULT(ke)
       IMPLICIT NONE
 
       CLASS(meshVol1DCartSegm), INTENT(in):: self
       REAL(8):: ke(1:2,1:2)
       REAL(8):: Xii(1:3)
       REAL(8):: dPsi(1:1, 1:2)
-      REAL(8):: invJ
+      REAL(8):: invJ(1), detJ
+      INTEGER:: l
 
       ke      = 0.D0
-      Xii     = (/ 0.D0, 0.D0, 0.D0 /)
+      Xii     = 0.D0
-      dPsi    = self%dPsi(Xii)
+
-      invJ    = self%invJac(Xii, dPsi)
+      DO l = 1, 3
-      ke(1,:) = (/ dPsi(1,1)*dPsi(1,1), dPsi(1,1)*dPsi(1,2) /) 
+        xii(1) = corSeg(l)
-      ke(2,:) = (/ dPsi(1,2)*dPsi(1,1), dPsi(1,2)*dPsi(1,2) /) 
+        dPsi    = self%dPsi(Xii)
-      ke      = 2.D0*ke*invJ
+        detJ    = self%detJac(Xii, dPsi)
+        invJ    = self%invJac(Xii, dPsi)
+        ke      = ke + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), &
+                              RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* &
+                       wSeg(l)/detJ
+
+      END DO
 
     END FUNCTION elemKSegm
 
@@ -348,14 +342,22 @@ MODULE moduleMesh1DCart
       REAL(8), INTENT(in):: source(1:)
       REAL(8), ALLOCATABLE:: localF(:)
       REAL(8):: fPsi(1:2)
-      REAL(8):: detJ
+      REAL(8):: detJ, f
       REAL(8):: Xii(1:3)
+      INTEGER:: l
 
-      Xii = 0.D0
-      fPsi = self%fPsi(Xii)
-      detJ = self%detJac(Xii)
       ALLOCATE(localF(1:2))
-      localF = 2.D0*DOT_PRODUCT(fPsi, source)*detJ
+      localF = 0.D0
+      Xii    = 0.D0
+
+      DO l = 1, 3
+        xii(1) = corSeg(l)
+        detJ   = self%detJac(Xii)
+        fPsi   = self%fPsi(Xii)
+        f      = DOT_PRODUCT(fPsi, source)
+        localF = localF + f*fPsi*wSeg(l)*detJ
+
+      END DO
 
     END FUNCTION elemFSegm
 

src/modules/mesh/1DRad/moduleMesh1DRad.f90

@@ -153,23 +153,7 @@ MODULE moduleMesh1DRad
       ALLOCATE(self%fboundary(1:nSpecies))
       !Assign functions to boundary
       DO s = 1, nSpecies
-        SELECT TYPE(obj => self%boundary%bTypes(s)%obj)
+        CALL pointBoundaryFunction(self, s)
-        TYPE IS(boundaryAbsorption)
-          self%fBoundary(s)%apply => absorption
-
-        TYPE IS(boundaryReflection)
-          self%fBoundary(s)%apply => reflection
-
-        TYPE IS(boundaryTransparent)
-          self%fBoundary(s)%apply => transparent
-
-        TYPE IS(boundaryWallTemperature)
-          self%fBoundary(s)%apply => wallTemperature
-
-        CLASS DEFAULT
-          CALL criticalError("Boundary type not defined in this geometry", 'initEdge1DRad')
-
-        END SELECT
 
       END DO
 

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

@@ -200,23 +200,7 @@ MODULE moduleMesh2DCart
       ALLOCATE(self%fboundary(1:nSpecies))
       !Assign functions to boundary
       DO s = 1, nSpecies
-        SELECT TYPE(obj => self%boundary%bTypes(s)%obj)
+        CALL pointBoundaryFunction(self, s)
-        TYPE IS(boundaryAbsorption)
-          self%fBoundary(s)%apply => absorption
-
-        TYPE IS(boundaryReflection)
-          self%fBoundary(s)%apply => reflection
-
-        TYPE IS(boundaryTransparent)
-          self%fBoundary(s)%apply => transparent
-
-        TYPE IS(boundaryWallTemperature)
-          self%fBoundary(s)%apply => wallTemperature
-
-        CLASS DEFAULT
-          CALL criticalError("Boundary type not defined in this geometry", 'initEdge2DCart')
-
-        END SELECT
 
       END DO
 

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

@@ -200,26 +200,7 @@ MODULE moduleMesh2DCyl
       ALLOCATE(self%fboundary(1:nSpecies))
       !Assign functions to boundary
       DO s = 1, nSpecies
-        SELECT TYPE(obj => self%boundary%bTypes(s)%obj)
+        CALL pointBoundaryFunction(self, s)
-        TYPE IS(boundaryAbsorption)
-          self%fBoundary(s)%apply => absorption
-
-        TYPE IS(boundaryReflection)
-          self%fBoundary(s)%apply => reflection
-
-        TYPE IS(boundaryTransparent)
-          self%fBoundary(s)%apply => transparent
-
-        TYPE IS(boundaryWallTemperature)
-          self%fBoundary(s)%apply => wallTemperature
-
-        TYPE IS(boundaryAxis)
-          self%fBoundary(s)%apply => symmetryAxis
-
-        CLASS DEFAULT
-          CALL criticalError("Boundary type not defined in this geometry", 'initEdge2DCyl')
-
-        END SELECT
 
       END DO
 
@@ -459,7 +440,9 @@ MODULE moduleMesh2DCyl
           detJ      = self%detJac(xi,dPsi)
           invJ      = self%invJac(xi,dPsi)
           r         = DOT_PRODUCT(fPsi,self%r)
-          ke        = ke + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*r*wQuad(l)*wQuad(m)/detJ
+          ke        = ke + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), &
+                                            MATMUL(invJ,dPsi))* &
+                           r*wQuad(l)*wQuad(m)/detJ
 
         END DO
       END DO

src/modules/mesh/3DCart/moduleMesh3DCart.f90

@@ -166,23 +166,7 @@ MODULE moduleMesh3DCart
       ALLOCATE(self%fBoundary(1:nSpecies))
       !Assign functions to boundary
       DO s = 1, nSpecies
-        SELECT TYPE(obj => self%boundary%bTypes(s)%obj)
+        CALL pointBoundaryFunction(self, s)
-        TYPE IS(boundaryAbsorption)
-          self%fBoundary(s)%apply => absorption
-
-        TYPE IS(boundaryReflection)
-          self%fBoundary(s)%apply => reflection
-
-        TYPE IS(boundaryTransparent)
-          self%fBoundary(s)%apply => transparent
-
-        TYPE IS(boundaryWallTemperature)
-          self%fBoundary(s)%apply => wallTemperature
-
-        CLASS DEFAULT
-          CALL criticalError("Boundary type not defined in this geometry", 'initEdge3DCart')
-
-        END SELECT
 
       END DO
 

src/modules/mesh/moduleMeshBoundary.f90

@@ -108,6 +108,85 @@ MODULE moduleMeshBoundary
 
     END SUBROUTINE wallTemperature
 
+    !Ionization surface: an electron will pass through the surface
+    ! and create an ion-electron pair based on a neutral background
+    SUBROUTINE ionization(edge, part)
+      USE moduleList
+      USE moduleSpecies
+      USE moduleMesh
+      USE moduleRefParam
+      IMPLICIT NONE
+
+      CLASS(meshEdge), INTENT(inout):: edge
+      CLASS(particle), INTENT(inout):: part
+      REAL(8):: vRel, eRel, mRel !relative velocity, energy and mass
+      REAL(8):: ionizationRate
+      INTEGER:: ionizationPair, p
+      TYPE(particle), POINTER:: newElectron
+      TYPE(particle), POINTER:: newIon
+
+      SELECT TYPE(bound => edge%boundary%bTypes(part%sp)%obj)
+      TYPE IS(boundaryIonization)
+        mRel = (bound%m0*species(part%sp)%obj%m)*(bound%m0+species(part%sp)%obj%m)
+        vRel = SUM(DABS(part%v-bound%v0))
+        eRel = mRel*vRel**2*5.D-1
+
+        IF (eRel > bound%eThreshold) THEN
+          !TODO: Check units
+          ionizationRate = bound%n0*bound%crossSection%get(eRel)
+
+          ionizationPair = INT(ionizationRate*tauMin*ti_ref/species(bound%sp)%obj%weight)
+
+          !Create the new pair of particles
+          DO p = 1, ionizationPair
+            ALLOCATE(newElectron)
+            ALLOCATE(newIon)
+
+            newElectron%sp = part%sp
+            newIon%sp      = bound%sp
+
+            newElectron%v = 10.D0*bound%v0 + (1.D0 + bound%deltaV/NORM2(bound%v0))
+            newIon%v      = bound%v0
+
+            newElectron%r = edge%randPos()
+            newIon%r      = newElectron%r
+
+            newElectron%vol = part%vol
+            newIon%vol      = part%vol
+
+            newElectron%xi = mesh%vols(part%vol)%obj%phy2log(newElectron%r)
+            newIon%xi      = newElectron%xi
+
+            newElectron%qm = part%qm
+            SELECT TYPE(spe => species(bound%sp)%obj)
+            TYPE IS(speciesCharged)
+              newIon%qm = spe%qm
+
+            END SELECT
+
+            newElectron%weight = species(bound%sp)%obj%weight
+            newIon%weight      = newElectron%weight
+
+            newElectron%n_in = .TRUE.
+            newIon%n_in      = .TRUE.
+
+            !Add particles to list
+            CALL OMP_SET_LOCK(lockSurfaces)
+            CALL partSurfaces%add(newElectron)
+            CALL partSurfaces%add(newIon)
+            CALL OMP_UNSET_LOCK(lockSurfaces)
+
+          END DO
+
+          !Removes ionizing electron
+          part%n_in = .FALSE.
+
+        END IF
+
+      END SELECT
+
+    END SUBROUTINE ionization
+
     !Symmetry axis. Dummy function
     SUBROUTINE symmetryAxis(edge, part)
       USE moduleSpecies
@@ -118,4 +197,38 @@ MODULE moduleMeshBoundary
 
     END SUBROUTINE symmetryAxis
 
+    !Points the boundary function to specific type
+    SUBROUTINE pointBoundaryFunction(edge, s)
+      USE moduleErrors
+      IMPLICIT NONE
+
+      CLASS(meshEdge), INTENT(inout):: edge
+      INTEGER, INTENT(in):: s !Species index
+
+      SELECT TYPE(obj => edge%boundary%bTypes(s)%obj)
+      TYPE IS(boundaryAbsorption)
+        edge%fBoundary(s)%apply => absorption
+
+      TYPE IS(boundaryReflection)
+        edge%fBoundary(s)%apply => reflection
+
+      TYPE IS(boundaryTransparent)
+        edge%fBoundary(s)%apply => transparent
+
+      TYPE IS(boundaryWallTemperature)
+        edge%fBoundary(s)%apply => wallTemperature
+
+      TYPE IS(boundaryIonization)
+        edge%fBoundary(s)%apply => ionization
+
+      TYPE IS(boundaryAxis)
+        edge%fBoundary(s)%apply => symmetryAxis
+
+      CLASS DEFAULT
+        CALL criticalError("Boundary type not defined in this geometry", 'pointBoundaryFunction')
+
+      END SELECT
+
+    END SUBROUTINE pointBoundaryFunction
+
 END MODULE moduleMeshBoundary

src/modules/moduleBoundary.f90

@@ -1,4 +1,5 @@
 MODULE moduleBoundary
+  USE moduleTable
 
   !Generic type for boundaries
   TYPE, PUBLIC:: boundaryGeneric
@@ -24,7 +25,7 @@ MODULE moduleBoundary
 
   END TYPE boundaryTransparent
 
-  !Transparent boundary
+  !Wall Temperature boundary
   TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryWallTemperature
     !Thermal velocity of the wall: square root(Wall temperature X specific heat)
     REAL(8):: vTh
@@ -32,6 +33,18 @@ MODULE moduleBoundary
 
   END TYPE boundaryWallTemperature
 
+  !Ionization boundary
+  TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryIonization
+    REAL(8):: m0, n0, v0(1:3), T0 !Properties of background neutrals.
+    INTEGER:: sp !Ion species
+    TYPE(table1D):: crossSection
+    REAL(8):: eThreshold
+    REAL(8):: deltaV
+
+    CONTAINS
+
+  END TYPE boundaryIonization
+
   !Symmetry axis
   TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryAxis
     CONTAINS
@@ -86,4 +99,34 @@ MODULE moduleBoundary
 
     END SUBROUTINE initWallTemperature
 
+    SUBROUTINE initIonization(boundary, m0, n0, v0, T0, speciesID, crossSection, eThreshold)
+      USE moduleRefParam
+      USE moduleSpecies
+      USE moduleCaseParam
+      USE moduleConstParam
+      IMPLICIT NONE
+
+      CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary
+      REAL(8), INTENT(in):: m0, n0, v0(1:3), T0
+      INTEGER:: speciesID
+      CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSection
+      REAL(8), INTENT(in):: eThreshold
+
+      ALLOCATE(boundaryIonization:: boundary)
+
+      SELECT TYPE(boundary)
+      TYPE IS(boundaryIonization)
+        boundary%m0 = m0 / m_ref
+        boundary%n0 = n0
+        boundary%v0 = v0 / v_ref
+        boundary%T0 = T0 / T_ref
+        boundary%sp = speciesID
+        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)
+
+      END SELECT
+
+    END SUBROUTINE initIonization
+
 END MODULE moduleBoundary

src/modules/moduleInject.f90

@@ -98,6 +98,10 @@ MODULE moduleInject
         !Input current in Ampers
         self%nParticles = INT(flow*tauMin*ti_ref/(qe*species(sp)%obj%weight))
 
+      CASE ("part/s")
+        !Input current in Ampers
+        self%nParticles = INT(flow*tauMin*ti_ref/species(sp)%obj%weight)
+
       CASE DEFAULT
         CALL criticalError("No support for units: " // units, 'initInject')
         

src/modules/moduleInput.f90

@@ -599,6 +599,13 @@ MODULE moduleInput
       CHARACTER(2):: istring, sString
       CHARACTER(:), ALLOCATABLE:: object, bType
       REAL(8):: Tw, cw !Wall temperature and specific heat
+      !Neutral Properties
+      REAL(8):: m0, n0
+      REAL(8), DIMENSION(:), ALLOCATABLE:: v0
+      REAL(8):: T0
+      REAL(8):: eThreshold !Energy threshold
+      INTEGER:: speciesID
+      CHARACTER(:), ALLOCATABLE:: speciesName, crossSection
       LOGICAL:: found
       INTEGER:: nTypes
 
@@ -628,6 +635,25 @@ MODULE moduleInput
           CASE('transparent')
             ALLOCATE(boundaryTransparent:: boundary(i)%bTypes(s)%obj)
 
+          CASE('ionization')
+            CALL config%get(object // '.neutral.name', speciesName, found)
+            IF (.NOT. found) CALL criticalError("missing parameter 'name' for neutrals in ionization", 'readBoundary')
+            speciesID = speciesName2Index(speciesName)
+            CALL config%get(object // '.neutral.mass', m0, found)
+            IF (.NOT. found) CALL criticalError("missing parameter 'mass' for neutrals in ionization", 'readBoundary')
+            CALL config%get(object // '.neutral.density', n0, found)
+            IF (.NOT. found) CALL criticalError("missing parameter 'density' for neutrals in ionization", 'readBoundary')
+            CALL config%get(object // '.neutral.velocity', v0, found)
+            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 // '.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)
+
           CASE('wallTemperature')
             CALL config%get(object // '.temperature', Tw, found)
             IF (.NOT. found) CALL criticalError("temperature not found for wallTemperature boundary type", 'readBoundary')

src/modules/moduleList.f90

@@ -24,6 +24,8 @@ MODULE moduleList
   INTEGER(KIND=OMP_LOCK_KIND):: lockWScheme !Lock for the NA list of particles
   TYPE(listNode):: partCollisions !Particles created in collisional process
   INTEGER(KIND=OMP_LOCK_KIND):: lockCollisions !Lock for the NA list of particles
+  TYPE(listNode):: partSurfaces !Particles created in surface interactions
+  INTEGER(KIND=OMP_LOCK_KIND):: lockSurfaces !Lock for the NA list of particles
   TYPE(listNode):: partInitial !Initial distribution of particles
 
   TYPE pointerArray

src/modules/moduleSolver.f90

@@ -437,7 +437,7 @@ MODULE moduleSolver
 
       INTEGER:: nn, n, e
       INTEGER, SAVE:: nPartNew
-      INTEGER, SAVE:: nInjIn, nOldIn, nWScheme, nCollisions
+      INTEGER, SAVE:: nInjIn, nOldIn, nWScheme, nCollisions, nSurfaces
       TYPE(particle), ALLOCATABLE, SAVE:: partTemp(:)
       TYPE(lNode), POINTER:: partCurr, partNext
 
@@ -458,13 +458,15 @@ MODULE moduleSolver
       nWScheme = partWScheme%amount
       !$OMP SECTION
       nCollisions = partCollisions%amount
+      !$OMP SECTION
+      nSurfaces = partSurfaces%amount
       !$OMP END SECTIONS
 
       !$OMP BARRIER
 
       !$OMP SINGLE
       CALL MOVE_ALLOC(partOld, partTemp)
-      nPartNew = nInjIn + nOldIn + nWScheme + nCollisions
+      nPartNew = nInjIn + nOldIn + nWScheme + nCollisions + nSurfaces
       ALLOCATE(partOld(1:nPartNew))
       !$OMP END SINGLE
 
@@ -522,6 +524,21 @@ MODULE moduleSolver
       IF (ASSOCIATED(partCollisions%tail)) NULLIFY(partCollisions%tail)
       partCollisions%amount = 0
 
+      !$OMP SECTION
+      !Reset particles from surface process
+      nn = nInjIn + nOldIn + nWScheme + nCollisions
+      partCurr => partSurfaces%head
+      DO n = 1, nSurfaces
+        partNext => partCurr%next
+        partOld(nn+n) = partCurr%part
+        DEALLOCATE(partCurr)
+        partCurr => partNext
+
+      END DO
+      IF (ASSOCIATED(partSurfaces%head)) NULLIFY(partSurfaces%head)
+      IF (ASSOCIATED(partSurfaces%tail)) NULLIFY(partSurfaces%tail)
+      partSurfaces%amount = 0
+
       !$OMP SECTION
       !Reset output in nodes
       DO e = 1, mesh%numNodes