Boundary ionization now outputs the number of ionization events during the iteration

runs/ALPHIE_Grid/inputIonization_0.10mA.json

@@ -5,6 +5,7 @@
     "cpuTime": false,
     "numColl": false,
     "EMField": true,
+    "boundariesParticle": true,
     "folder": "ionization_0.10mA"
   },
   "geometry": {

src/modules/mesh/moduleMesh.f90

@@ -694,7 +694,7 @@ MODULE moduleMesh
       use moduleSpecies
       import boundaryParticleGeneric, meshEdge
 
-      class(boundaryParticleGeneric), intent(in):: self
+      class(boundaryParticleGeneric), intent(inout):: self
       class(meshEdge), intent(inout):: edge
       class(particle), intent(inout):: part
 
@@ -712,7 +712,7 @@ MODULE moduleMesh
     subroutine printParticle_interface(self, fileID)
       import boundaryParticleGeneric
 
-      class(boundaryParticleGeneric), intent(in):: self
+      class(boundaryParticleGeneric), intent(inout):: self
       integer, intent(in):: fileID
 
     end subroutine printParticle_interface
@@ -765,6 +765,8 @@ MODULE moduleMesh
     REAL(8):: effectiveTime
     REAL(8):: eThreshold
     REAL(8):: deltaV
+    integer:: tally
+    integer(KIND=OMP_LOCK_KIND):: tally_lock
     CONTAINS
       procedure, pass:: apply => ionization
 
@@ -785,7 +787,7 @@ MODULE moduleMesh
     module subroutine reflection(self, edge, part)
       use moduleSpecies
 
-      class(boundaryReflection), intent(in):: self
+      class(boundaryReflection), intent(inout):: self
       class(meshEdge), intent(inout):: edge
       class(particle), intent(inout):: part
 
@@ -794,7 +796,7 @@ MODULE moduleMesh
     module subroutine absorption(self, edge, part)
       use moduleSpecies
 
-      class(boundaryAbsorption), intent(in):: self
+      class(boundaryAbsorption), intent(inout):: self
       class(meshEdge), intent(inout):: edge
       class(particle), intent(inout):: part
 
@@ -803,7 +805,7 @@ MODULE moduleMesh
     module subroutine transparent(self, edge, part)
       use moduleSpecies
 
-      class(boundaryTransparent), intent(in):: self
+      class(boundaryTransparent), intent(inout):: self
       class(meshEdge), intent(inout):: edge
       class(particle), intent(inout):: part
 
@@ -812,7 +814,7 @@ MODULE moduleMesh
     module subroutine symmetryAxis(self, edge, part)
       use moduleSpecies
 
-      class(boundaryAxis), intent(in):: self
+      class(boundaryAxis), intent(inout):: self
       class(meshEdge), intent(inout):: edge
       class(particle), intent(inout):: part
 
@@ -821,7 +823,7 @@ MODULE moduleMesh
     module subroutine wallTemperature(self, edge, part)
       use moduleSpecies
 
-      class(boundaryWallTemperature), intent(in):: self
+      class(boundaryWallTemperature), intent(inout):: self
       class(meshEdge), intent(inout):: edge
       class(particle), intent(inout):: part
 
@@ -830,7 +832,7 @@ MODULE moduleMesh
     module subroutine ionization(self, edge, part)
       use moduleSpecies
 
-      class(boundaryIonization), intent(in):: self
+      class(boundaryIonization), intent(inout):: self
       class(meshEdge), intent(inout):: edge
       class(particle), intent(inout):: part
 
@@ -839,7 +841,7 @@ MODULE moduleMesh
     module subroutine quasiNeutrality(self, edge, part)
       use moduleSpecies
 
-      class(boundaryQuasiNeutrality), intent(in):: self
+      class(boundaryQuasiNeutrality), intent(inout):: self
       class(meshEdge), intent(inout):: edge
       class(particle), intent(inout):: part
 

src/modules/mesh/moduleMesh@boundaryParticle.f90

@@ -1,6 +1,7 @@
 !moduleMeshBoundary: Boundary functions for the mesh edges
 submodule(moduleMesh) boundaryParticle
   contains
+    ! Returns the array index of the boundary based on the name
     module function boundaryParticleName_to_Index(boundaryName) result(bp)
       use moduleErrors
       implicit none
@@ -25,6 +26,101 @@ submodule(moduleMesh) boundaryParticle
     
     end function boundaryParticleName_to_Index
 
+    module SUBROUTINE reflection(self, edge, part)
+      USE moduleCaseParam
+      USE moduleSpecies
+      IMPLICIT NONE
+
+      class(boundaryReflection), intent(inout):: self
+      CLASS(meshEdge), INTENT(inout):: edge
+      CLASS(particle), INTENT(inout):: part
+      !rp  = intersection between particle and edge
+      !rpp = final position of particle
+      !vpp = final velocity of particle
+      REAL(8), DIMENSION(1:3):: rp, vpp
+
+      !Reflect particle velocity
+      vpp = part%v - 2.D0*DOT_PRODUCT(part%v, edge%normal)*edge%normal
+      part%v = vpp
+
+      rp = edge%intersection(part%r)
+
+      part%r = 2.D0*(rp - part%r) + part%r
+
+      !particle is assumed to be inside
+      part%n_in = .TRUE.
+
+    END SUBROUTINE reflection
+
+    !Absoption in a surface
+    module SUBROUTINE absorption(self, edge, part)
+      USE moduleCaseParam
+      USE moduleSpecies
+      IMPLICIT NONE
+
+      class(boundaryAbsorption), intent(inout):: self
+      CLASS(meshEdge), INTENT(inout):: edge
+      CLASS(particle), INTENT(inout):: part
+      REAL(8):: rpp(1:3) !Position of particle projected to the edge
+      REAL(8):: d !Distance from particle to edge
+
+      rpp = edge%intersection(part%r)
+
+      d = NORM2(rpp - part%r)
+
+      IF (d >= 0.D0) THEN
+        part%weight = part%weight/d
+
+      END IF
+
+      !Assign new position to particle
+      part%r = rpp
+      !Remove particle from the domain
+      part%n_in = .FALSE.
+
+      !Scatter particle in associated volume
+      IF (ASSOCIATED(edge%e1)) THEN
+        CALL edge%e1%scatter(edge%e1%nNodes, part)
+
+      ELSE
+        CALL edge%e2%scatter(edge%e2%nNodes, part)
+
+      END IF
+
+    END SUBROUTINE absorption
+
+    !Transparent boundary condition
+    module SUBROUTINE transparent(self, edge, part)
+      USE moduleSpecies
+      IMPLICIT NONE
+
+      class(boundaryTransparent), intent(inout):: self
+      CLASS(meshEdge), INTENT(inout):: edge
+      CLASS(particle), INTENT(inout):: part
+
+      !Removes particle from domain
+      part%n_in = .FALSE.
+
+    END SUBROUTINE transparent
+
+    !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.
+    module SUBROUTINE symmetryAxis(self, edge, part)
+      USE moduleSpecies
+      IMPLICIT NONE
+
+      class(boundaryAxis), intent(inout):: self
+      CLASS(meshEdge), INTENT(inout):: edge
+      CLASS(particle), INTENT(inout):: part
+
+      CALL genericReflection(edge, part)
+
+    END SUBROUTINE symmetryAxis
+
+    ! wallTemperature
+    ! During the collision, the boundary transfers part of the energy to the incident particle
+    ! Init
     module SUBROUTINE initWallTemperature(boundary, T, c)
       USE moduleRefParam
       IMPLICIT NONE
@@ -44,12 +140,37 @@ submodule(moduleMesh) boundaryParticle
 
     END SUBROUTINE initWallTemperature
 
+    ! Apply
+    module SUBROUTINE wallTemperature(self, edge, part)
+      USE moduleSpecies
+      USE moduleRandom
+      IMPLICIT NONE
+
+      class(boundaryWallTemperature), intent(inout):: self
+      CLASS(meshEdge), INTENT(inout):: edge
+      CLASS(particle), INTENT(inout):: part
+      INTEGER:: i
+
+      !Modifies particle velocity according to wall temperature
+      DO i = 1, 3
+        part%v(i) = part%v(i) + self%vTh*randomMaxwellian()
+
+      END DO
+          
+      CALL genericReflection(edge, part)
+
+    END SUBROUTINE wallTemperature
+
+    ! ionization
+    ! An electron will pass through the surface and create a series of ion-electron pairs based on a neutral background
+    !init
     module SUBROUTINE initIonization(boundary, mImpact, m0, n0, v0, T0, ion, effTime, crossSection, eThreshold, electronSecondary)
       use moduleRefParam
       use moduleSpecies
       use moduleCaseParam
       use moduleConstParam
       use moduleErrors, only: criticalError
+      use omp_lib, only: omp_init_lock
       implicit none
 
       class(boundaryParticleGeneric), allocatable, intent(inout):: boundary
@@ -91,147 +212,17 @@ submodule(moduleMesh) boundaryParticle
         boundary%eThreshold = eThreshold*eV2J/(m_ref*v_ref**2)
         boundary%deltaV = DSQRT(boundary%eThreshold/mImpact)
 
+        boundary%update => ionization_resetTally
+        boundary%print  => ionization_print
+
+        boundary%tally = 0
+        call omp_init_lock(boundary%tally_lock)
+
       END SELECT
 
     END SUBROUTINE initIonization
 
-    module subroutine initQuasiNeutrality(boundary, s_incident)
+    ! Apply
-      implicit none
-
-      class(boundaryParticleGeneric), allocatable, intent(inout):: boundary
-      integer, intent(in):: s_incident
-
-      allocate(boundaryQuasiNeutrality:: boundary)
-
-      select type(boundary)
-      type is(boundaryQuasiNeutrality)
-        allocate(boundary%edges(0))
-
-        boundary%s_incident = s_incident
-
-
-        boundary%update => quasiNeutrality_update
-        boundary%print  => quasiNeutrality_print
-
-      end select
-
-    end subroutine initQuasiNeutrality
-
-    module SUBROUTINE reflection(self, edge, part)
-      USE moduleCaseParam
-      USE moduleSpecies
-      IMPLICIT NONE
-
-      class(boundaryReflection), intent(in):: self
-      CLASS(meshEdge), INTENT(inout)::        edge
-      CLASS(particle), INTENT(inout)::        part
-      !rp  = intersection between particle and edge
-      !rpp = final position of particle
-      !vpp = final velocity of particle
-      REAL(8), DIMENSION(1:3):: rp, vpp
-
-      !Reflect particle velocity
-      vpp = part%v - 2.D0*DOT_PRODUCT(part%v, edge%normal)*edge%normal
-      part%v = vpp
-
-      rp = edge%intersection(part%r)
-
-      part%r = 2.D0*(rp - part%r) + part%r
-
-      !particle is assumed to be inside
-      part%n_in = .TRUE.
-
-    END SUBROUTINE reflection
-
-    !Absoption in a surface
-    module SUBROUTINE absorption(self, edge, part)
-      USE moduleCaseParam
-      USE moduleSpecies
-      IMPLICIT NONE
-
-      class(boundaryAbsorption), intent(in):: self
-      CLASS(meshEdge), INTENT(inout)::        edge
-      CLASS(particle), INTENT(inout)::        part
-      REAL(8):: rpp(1:3) !Position of particle projected to the edge
-      REAL(8):: d !Distance from particle to edge
-
-      rpp = edge%intersection(part%r)
-
-      d = NORM2(rpp - part%r)
-
-      IF (d >= 0.D0) THEN
-        part%weight = part%weight/d
-
-      END IF
-
-      !Assign new position to particle
-      part%r = rpp
-      !Remove particle from the domain
-      part%n_in = .FALSE.
-
-      !Scatter particle in associated volume
-      IF (ASSOCIATED(edge%e1)) THEN
-        CALL edge%e1%scatter(edge%e1%nNodes, part)
-
-      ELSE
-        CALL edge%e2%scatter(edge%e2%nNodes, part)
-
-      END IF
-
-    END SUBROUTINE absorption
-
-    !Transparent boundary condition
-    module SUBROUTINE transparent(self, edge, part)
-      USE moduleSpecies
-      IMPLICIT NONE
-
-      class(boundaryTransparent), intent(in):: self
-      CLASS(meshEdge), INTENT(inout)::         edge
-      CLASS(particle), INTENT(inout)::         part
-
-      !Removes particle from domain
-      part%n_in = .FALSE.
-
-    END SUBROUTINE transparent
-
-    !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.
-    module SUBROUTINE symmetryAxis(self, edge, part)
-      USE moduleSpecies
-      IMPLICIT NONE
-
-      class(boundaryAxis), intent(in):: self
-      CLASS(meshEdge), INTENT(inout)::  edge
-      CLASS(particle), INTENT(inout)::  part
-
-      CALL genericReflection(edge, part)
-
-    END SUBROUTINE symmetryAxis
-
-    !Wall with temperature
-    module SUBROUTINE wallTemperature(self, edge, part)
-      USE moduleSpecies
-      USE moduleRandom
-      IMPLICIT NONE
-
-      class(boundaryWallTemperature), intent(in):: self
-      CLASS(meshEdge), INTENT(inout)::               edge
-      CLASS(particle), INTENT(inout)::               part
-      INTEGER:: i
-
-      !Modifies particle velocity according to wall temperature
-      DO i = 1, 3
-        part%v(i) = part%v(i) + self%vTh*randomMaxwellian()
-
-      END DO
-          
-      CALL genericReflection(edge, part)
-
-    END SUBROUTINE wallTemperature
-
-    !Ionization surface: an electron will pass through the surface
-    ! and create an ion-electron pair based on a neutral background
     module SUBROUTINE ionization(self, edge, part)
       USE moduleList
       USE moduleSpecies
@@ -239,9 +230,10 @@ submodule(moduleMesh) boundaryParticle
       USE moduleRefParam
       USE moduleRandom
       USE moduleMath
+      use omp_lib, only: omp_set_lock, omp_unset_lock
       IMPLICIT NONE
 
-      class(boundaryIonization), intent(in):: self
+      class(boundaryIonization), intent(inout):: self
       CLASS(meshEdge), INTENT(inout)::        edge
       CLASS(particle), INTENT(inout)::        part
       REAL(8):: vRel, eRel, mRel !relative velocity, energy and mass
@@ -320,6 +312,11 @@ submodule(moduleMesh) boundaryParticle
           !Reduce number of possible ionizations
           nIonizations = nIonizations - 1
 
+          ! Increase the tally
+          call omp_set_lock(self%tally_lock)
+          self%tally = self%tally + 1
+          call omp_unset_lock(self%tally_lock)
+
         END IF
 
       END DO
@@ -329,11 +326,68 @@ submodule(moduleMesh) boundaryParticle
 
     END SUBROUTINE ionization
 
+    ! Update
+    subroutine ionization_resetTally(self)
+      implicit none
+
+      class(boundaryParticleGeneric), intent(inout):: self
+
+      select type(self)
+      type is(boundaryIonization)
+        self%tally = 0
+
+      end select
+
+    end subroutine ionization_resetTally
+
+    ! Print
+    subroutine ionization_print(self, fileID)
+      implicit none
+
+      class(boundaryParticleGeneric), intent(inout):: self
+      integer, intent(in):: fileID
+      integer:: e
+
+      write(fileID, '(A)') self%name
+      select type(self)
+      type is(boundaryIonization)
+        write(fileID, '(A)') 'Number of successful ionization events per iteration'
+        write(fileID, '('//fmtInt//')') self%tally
+
+      end select
+
+    end subroutine ionization_print
+
+    ! quasiNeutrality
+    ! Maintains n_incident = n_rest by modifying the reflection parameter of the edge.
+    ! Init
+    module subroutine initQuasiNeutrality(boundary, s_incident)
+      implicit none
+
+      class(boundaryParticleGeneric), allocatable, intent(inout):: boundary
+      integer, intent(in):: s_incident
+
+      allocate(boundaryQuasiNeutrality:: boundary)
+
+      select type(boundary)
+      type is(boundaryQuasiNeutrality)
+        allocate(boundary%edges(0))
+
+        boundary%s_incident = s_incident
+
+        boundary%update => quasiNeutrality_update
+        boundary%print  => quasiNeutrality_print
+
+      end select
+
+    end subroutine initQuasiNeutrality
+
+    ! Apply
     module subroutine quasiNeutrality(self, edge, part)
       use moduleRandom
       implicit none
 
-      class(boundaryQuasiNeutrality), intent(in):: self
+      class(boundaryQuasiNeutrality), intent(inout):: self
       class(meshEdge), intent(inout):: edge
       class(particle), intent(inout):: part
 
@@ -347,6 +401,7 @@ submodule(moduleMesh) boundaryParticle
 
     end subroutine quasiNeutrality
 
+    ! Update
     subroutine quasiNeutrality_update(self)
       implicit none
 
@@ -404,10 +459,11 @@ submodule(moduleMesh) boundaryParticle
 
     end subroutine quasiNeutrality_update
 
+    ! Print output
     subroutine quasiNeutrality_print(self, fileID)
       implicit none
 
-      class(boundaryParticleGeneric), intent(in):: self
+      class(boundaryParticleGeneric), intent(inout):: self
       integer, intent(in):: fileID
       integer:: e
 
@@ -416,7 +472,7 @@ submodule(moduleMesh) boundaryParticle
       type is(boundaryQuasiNeutrality)
         write(fileID, '(A,",",A)') '"Edge id"', '"alpha"'
         do e = 1, size(self%edges)
-          write(fileID, '(I0,",",ES0.6E3)') self%edges(e)%obj%n, self%alpha(self%edges(e)%obj%n)
+          write(fileID, '('//fmtInt//',",",'//fmtReal//')') self%edges(e)%obj%n, self%alpha(self%edges(e)%obj%n)
 
         end do
 
@@ -425,6 +481,7 @@ submodule(moduleMesh) boundaryParticle
     end subroutine quasiNeutrality_print
 
     ! Generic boundary conditions for internal use
+    ! reflection
     module subroutine genericReflection(edge, part)
       use moduleCaseParam
       use moduleSpecies
@@ -450,6 +507,7 @@ submodule(moduleMesh) boundaryParticle
 
     end subroutine genericReflection
 
+    ! transparent
     module subroutine genericTransparent(edge, part)
       use moduleSpecies
       implicit none
@@ -462,6 +520,7 @@ submodule(moduleMesh) boundaryParticle
 
     end subroutine genericTransparent
 
+    ! Updates the boundary condition values when %update procedure is associated
     module subroutine boundariesParticle_update()
       implicit none
 
@@ -477,6 +536,7 @@ submodule(moduleMesh) boundaryParticle
 
     end subroutine boundariesParticle_update
 
+    ! Writes the output into the Step_XXXXX_boundaryParticles file when the %print procedure is associated
     module subroutine boundariesParticle_write()
       use moduleCaseparam, only: timeStep
       use moduleOutput, only:fileID_boundaryParticle, formatFileName, informFileCreation