Merge branch 'fix.surfChemFixedProb' into 'master.dev'

[fix.surfChemFixedProb] Resolve "Surface Chemistry Model: fixed probability"

Closes #269

See merge request piclas/piclas!953

docs/documentation/userguide/features-and-models/BC-particle-solver.md

@@ -283,7 +283,7 @@ As additional output, the cell-local sticking coefficient will be added to the s
 
 ### Fixed probability surface chemistry
 
-This simple fixed-probability surface chemistry model allows the user to define arbitrary surface reactions, by defining the impacting species, the products and a fixed event probability. The reaction is then assigned to the boundaries by specifying their number and index as defined previously.
+This simple fixed-probability surface chemistry model allows the user to define arbitrary surface reactions, by defining the impacting species, the products and a fixed event probability. The reaction is then assigned to the boundaries by specifying their number and index as defined previously. This model corresponds to `Part-BoundaryX-SurfaceModel = 2`, which is set automatically when a reaction of this type is defined.
 
     Surface-Reaction1-Type               = P
     Surface-Reaction1-Reactants          = (/1,0/)

regressioncheck/NIG_DSMC/SURF_PROB_MultiReac/parameter.ini

@@ -8,7 +8,6 @@ TimeStampLength = 16
 ! VARIABLES
 ! =============================================================================== !
 CFLscale      = 0.2
-IniExactFunc  = 0
 N             = 1
 NAnalyze      = 1
 ! =============================================================================== !
@@ -66,6 +65,8 @@ Part-Boundary1-Condition    = reflective
 Part-Boundary1-WallTemp     = 300.
 Part-Boundary1-MomentumACC  = 0.
 Part-Boundary1-TransACC     = 1.
+! Definition of SurfaceModel = 2 is not required but added here as a test
+Part-Boundary1-SurfaceModel = 2
 Part-Boundary2-SourceName   = BC_Xminus
 Part-Boundary2-Condition    = reflective
 Part-Boundary2-WallTemp     = 300.
@@ -95,7 +96,7 @@ Part-FIBGMdeltas          = (/2E-4,2E-4,2E-4/)
 ! =============================================================================== !
 ! SURFACE CHEMISTRY
 ! =============================================================================== !
-Surface-NumOfReactions = 7
+Surface-NumOfReactions = 8
 
 ! H -> H2 (recombination at the wall)
 Surface-Reaction1-Type = P
@@ -149,10 +150,18 @@ Surface-Reaction6-ProductAccommodation = 1.
 Surface-Reaction7-Type = P
 Surface-Reaction7-Reactants = (/6,0/)
 Surface-Reaction7-Products  = (/2,2,2/)
-Surface-Reaction7-NumOfBoundaries = 6
-Surface-Reaction7-Boundaries = (/1,2,3,4,5,6/)
+Surface-Reaction7-NumOfBoundaries = 3
+Surface-Reaction7-Boundaries = (/1,2,3/)
 Surface-Reaction7-EventProbability = 0.5
 Surface-Reaction7-ProductAccommodation = 1.
+! H3Ion -> 3H (reaction is duplicated to test the correct handling of multiple reactions for the same reactant species at different BCs)
+Surface-Reaction8-Type = P
+Surface-Reaction8-Reactants = (/6,0/)
+Surface-Reaction8-Products  = (/2,2,2/)
+Surface-Reaction8-NumOfBoundaries = 3
+Surface-Reaction8-Boundaries = (/4,5,6/)
+Surface-Reaction8-EventProbability = 0.5
+Surface-Reaction8-ProductAccommodation = 1.
 ! =============================================================================== !
 ! SPECIES
 ! =============================================================================== !
@@ -171,14 +180,6 @@ Part-Species1-ChargeIC              = 0.0
 ! =============================================================================== !
 Part-Species2-MassIC                      = 1.67400E-27
 Part-Species2-ChargeIC                    = 0.0
-
-Part-Species2-nInits = 0
-Part-Species2-Init1-velocityDistribution  = maxwell_lpn
-Part-Species2-Init1-SpaceIC               = cell_local
-Part-Species2-Init1-VeloIC                = 0.
-Part-Species2-Init1-PartDensity           = 1E21
-Part-Species2-Init1-VeloVecIC             = (/0.,1.,0./)
-Part-Species2-Init1-MWTemperatureIC       = 50000.
 ! =============================================================================== !
 ! Species3 | e
 ! =============================================================================== !
@@ -197,8 +198,6 @@ Part-Species4-Init1-VeloIC                = 0.
 Part-Species4-Init1-PartDensity           = 1E21
 Part-Species4-Init1-VeloVecIC             = (/0.,1.,0./)
 Part-Species4-Init1-MWTemperatureIC       = 50000.
-Part-Species4-Init1-TempVib               = 50000.
-Part-Species4-Init1-TempRot               = 50000.
 ! =============================================================================== !
 ! Species5 | HIon
 ! =============================================================================== !
@@ -225,68 +224,33 @@ Part-Species6-Init1-VeloIC                = 0.
 Part-Species6-Init1-PartDensity           = 1E21
 Part-Species6-Init1-VeloVecIC             = (/0.,1.,0./)
 Part-Species6-Init1-MWTemperatureIC       = 50000.
-Part-Species6-Init1-TempVib               = 50000.
-Part-Species6-Init1-TempRot               = 50000.
-
 ! =============================================================================== !
 ! Species1, H2
 ! =============================================================================== !
 Part-Species1-SpeciesName       = H2
 Part-Species1-InteractionID     = 2
-Part-Species1-Tref              = 1000
-Part-Species1-dref              = 2.68E-10
-Part-Species1-omega             = 0.407
-Part-Species1-HeatOfFormation_K = 0.0
-Part-Species1-CharaTempVib      = 6332.37
-Part-Species1-Ediss_eV          = 4.47
 ! =============================================================================== !
 ! Species2, H
 ! =============================================================================== !
 Part-Species2-SpeciesName       = H
 Part-Species2-InteractionID     = 1
-Part-Species2-Tref              = 1000
-Part-Species2-dref              = 2.581E-10
-Part-Species2-omega             = 0.407
-Part-Species2-HeatOfFormation_K = 26159.76
 ! =============================================================================== !
 ! Species3, e
 ! =============================================================================== !
 Part-Species3-SpeciesName       = electron
 Part-Species3-InteractionID     = 4
-Part-Species3-Tref              = 1000
-Part-Species3-dref              = 2.817920E-15
-Part-Species3-omega             = 0.407
 ! =============================================================================== !
 ! Species4, H2Ion
 ! =============================================================================== !
 Part-Species4-SpeciesName      = H2Ion1
 Part-Species4-InteractionID    = 20
-Part-Species4-Tref             = 1000
-Part-Species4-dref             = 3.883E-10
-Part-Species4-omega            = 0.407
-Part-Species4-CharaTempVib     = 3341.01
-Part-Species4-Ediss_eV         = 2.65
 ! =============================================================================== !
 ! Species5, HIon
 ! =============================================================================== !
 Part-Species5-SpeciesName       = HIon1
 Part-Species5-InteractionID     = 10
-Part-Species5-Tref              = 1000
-Part-Species5-dref              = 3.912E-10
-Part-Species5-omega             = 0.407
 ! =============================================================================== !
 ! Species6, H3Ion
 ! =============================================================================== !
 Part-Species6-SpeciesName       = H3Ion1
 Part-Species6-InteractionID     = 20
-Part-Species6-PolyatomicMol     = T
-Part-Species6-NumOfAtoms        = 3
-Part-Species6-LinearMolec       = F
-Part-Species6-Tref              = 1000
-Part-Species6-dref              = 4.5E-10       ! Guess
-Part-Species6-omega             = 0.407
-Part-Species6-CharaTempVib1     = 4572.92
-Part-Species6-CharaTempVib2     = 3627.94
-Part-Species6-CharaTempVib3     = 3627.94
-Part-Species6-Ediss_eV          = 4.51
-Part-Species6-HeatOfFormation_K = 132803.52

src/particles/boundary/particle_boundary_init.f90

@@ -466,7 +466,7 @@ SUBROUTINE InitializeVariablesPartBoundary()
     IF (PartBound%SurfaceModel(iPartBound).GT.0)THEN
       IF (.NOT.useDSMC) CALL abort(__STAMP__,'Cannot use surfacemodel>0 with useDSMC=F for particle boundary: ',iPartBound)
       SELECT CASE (PartBound%SurfaceModel(iPartBound))
-      CASE (0)
+      CASE (0,2)
         PartBound%Reactive(iPartBound)        = .FALSE.
       CASE (1)
         PartBound%Reactive(iPartBound)        = .FALSE.
@@ -478,7 +478,7 @@ SUBROUTINE InitializeVariablesPartBoundary()
         ! SEE models require reactive BC
         PartBound%Reactive(iPartBound)        = .TRUE.
       CASE DEFAULT
-        CALL abort(__STAMP__,'Error in particle init: only allowed SurfaceModels: 0,1,20,SEE_MODELS_ID! SurfaceModel=',&
+        CALL abort(__STAMP__,'Error in particle init: only allowed SurfaceModels: 0,1,2,20,SEE_MODELS_ID! SurfaceModel=',&
                   IntInfoOpt=PartBound%SurfaceModel(iPartBound))
       END SELECT
     END IF
@@ -496,9 +496,7 @@ SUBROUTINE InitializeVariablesPartBoundary()
     PartBound%MaxTotalCoverage(iPartBound) = GETREAL('Part-Boundary'//TRIM(hilf)//'-MaxTotalCoverage', '1.')
     ! Check if the maximum of the coverage is reached
     IF (PartBound%TotalCoverage(iPartBound).GT.PartBound%MaxTotalCoverage(iPartBound)) THEN
-      CALL abort(&
-    __STAMP__&
-    ,'ERROR: Maximum surface coverage reached.', iPartBound)
+      CALL abort(__STAMP__,'ERROR: Maximum surface coverage reached.', iPartBound)
     END IF
     IF (PartBound%NbrOfSpeciesSwaps(iPartBound).GT.0) THEN
       !read Species to be changed at wall (in, out), out=0: delete

src/particles/surfacemodel/surfacemodel_chemistry.f90

@@ -137,6 +137,7 @@ SUBROUTINE InitializeVariablesSurfaceChemistry()
 INTEGER               :: iReac, iReac2, iBound, iVal, err
 INTEGER               :: ReadInNumOfReact
 INTEGER               :: iSpec, SpecID, ReactionPathPerSpecies(nSpecies), BCID
+REAL                  :: ReacProbTest(nPartBound)
 !===================================================================================================================================
 
 IF(SurfChem%NumOfReact.LE.0) RETURN
@@ -170,8 +171,11 @@ SUBROUTINE InitializeVariablesSurfaceChemistry()
   IF (SurfChemReac(iReac)%NumOfBounds.EQ.0) THEN
     CALL abort(__STAMP__,'ERROR: At least one boundary must be defined for each surface reaction!',iReac)
   END IF
-
   SurfChemReac(iReac)%Boundaries = GETINTARRAY('Surface-Reaction'//TRIM(hilf)//'-Boundaries',SurfChemReac(iReac)%NumOfBounds)
+  ! Sanity check if the boundary index is within the range of nPartBound
+  IF(ANY(SurfChemReac(iReac)%Boundaries.LT.1).OR.ANY(SurfChemReac(iReac)%Boundaries.GT.nPartBound)) THEN
+    CALL abort(__STAMP__,'ERROR: Boundary index out of range for surface reaction: ', iReac)
+  END IF
   ! Define the surface model
   SELECT CASE (TRIM(SurfChemReac(iReac)%ReactType))
   CASE('P')
@@ -213,9 +217,9 @@ SUBROUTINE InitializeVariablesSurfaceChemistry()
   END IF
 END DO
 
+! Probability based surface chemistry model: Allocate the species specific type with the number of the possible reaction paths
 DO iSpec = 1, nSpecies
   IF(SurfChem%EventProbInfo(iSpec)%NumOfReactionPaths.GT.0) THEN
-    ! Allocate the species specific type with the number of the possible reaction paths
     ALLOCATE(SurfChem%EventProbInfo(iSpec)%ReactionIndex(SurfChem%EventProbInfo(iSpec)%NumOfReactionPaths))
     SurfChem%EventProbInfo(iSpec)%ReactionIndex = 0
     ALLOCATE(SurfChem%EventProbInfo(iSpec)%ReactionProb(SurfChem%EventProbInfo(iSpec)%NumOfReactionPaths))
@@ -415,8 +419,7 @@ SUBROUTINE InitializeVariablesSurfaceChemistry()
   CALL H5CLOSE_F(err)
 END IF !SpeciesDatabase
 
-
-IF (SurfChem%OverwriteCatParameters) THEN
+IF(SurfChem%OverwriteCatParameters) THEN
   ! Loop over the surface reactions
   DO iReac = 1, ReadInNumOfReact
     WRITE(UNIT=hilf,FMT='(I0)') iReac
@@ -489,13 +492,24 @@ SUBROUTINE InitializeVariablesSurfaceChemistry()
   END DO
 END IF
 
-! Sanity check: Total reaction probability of a single species at the surface must not be above 1
+! Sanity check: Total reaction probability of a single species at a specific boundary must not be above 1
 DO iSpec = 1, nSpecies
-  IF(SurfChem%EventProbInfo(iSpec)%NumOfReactionPaths.GT.0) THEN
-    IF(SUM(SurfChem%EventProbInfo(iSpec)%ReactionProb(:)).GT.1.) THEN
-      CALL abort(__STAMP__,'ERROR: Total probability above unity for species: ', IntInfoOpt=iSpec)
-    END IF
-  END IF
+  IF(SurfChem%EventProbInfo(iSpec)%NumOfReactionPaths.EQ.0) CYCLE
+  ReacProbTest = 0.
+  ! Loop over the reaction paths
+  DO iVal = 1, SurfChem%EventProbInfo(iSpec)%NumOfReactionPaths
+    iReac = SurfChem%EventProbInfo(iSpec)%ReactionIndex(iVal)
+    ! Loop over the boundaries, where the reaction occurs
+    DO iBound = 1, SurfChemReac(iReac)%NumOfBounds
+      ! Get the boundary index
+      BCID = SurfChemReac(iReac)%Boundaries(iBound)
+      ! Sum up the reaction probabilities for each boundary
+      ReacProbTest(BCID) = ReacProbTest(BCID) + SurfChem%EventProbInfo(iSpec)%ReactionProb(iVal)
+      IF(ReacProbTest(BCID).GT.1.) THEN
+        CALL abort(__STAMP__,'ERROR: Total probability above unity for species: ', IntInfoOpt=iSpec)
+      END IF
+    END DO
+  END DO
 END DO
 
 END SUBROUTINE InitializeVariablesSurfaceChemistry
@@ -1009,7 +1023,11 @@ SUBROUTINE SurfaceModelEventProbability(PartID,SideID,GlobalElemID,n_loc,PartPos
 ! 1a.) Determine which reaction path to follow
   CALL RANDOM_NUMBER(RanNum)
   DO iPath = 1, SurfChem%EventProbInfo(SpecID)%NumOfReactionPaths
+    ! Check if the reaction is allowed at the boundary
+    IF(.NOT.ANY(SurfChemReac(SurfChem%EventProbInfo(SpecID)%ReactionIndex(iPath))%Boundaries(:).EQ.locBCID)) CYCLE
+    ! Sum up the probabilities of the reaction paths (sanity check during initialization to ensure that the sum is below 1.0)
     TotalProb = TotalProb + SurfChem%EventProbInfo(SpecID)%ReactionProb(iPath)
+    ! Decide which reaction path to follow
     IF(TotalProb.GT.RanNum) THEN
       PathTodo = iPath
       EXIT