Merge branch 'fix.relativistic.chemistry.check' into 'master.dev'

[fix.relativistic.chemistry.check] Fixed bug in chemistry. Inconsistent consideration of relativistic effects in...

See merge request piclas/piclas!719

src/particles/dsmc/dsmc_chemical_reactions.f90

@@ -351,7 +351,7 @@ SUBROUTINE DSMC_Chemistry(iPair, iReac)
 ! Routine performs an exchange reaction of the type A + B + C -> D + E + F, where A, B, C, D, E, F can be anything
 !===================================================================================================================================
 ! MODULES
-USE MOD_Globals                ,ONLY: abort,DOTPRODUCT,StringBeginsWith
+USE MOD_Globals                ,ONLY: abort,DOTPRODUCT,StringBeginsWith,UNIT_StdOut,myrank
 USE MOD_Globals_Vars
 USE MOD_DSMC_Vars              ,ONLY: Coll_pData, DSMC, CollInf, SpecDSMC, DSMCSumOfFormedParticles, ElectronicDistriPart
 USE MOD_DSMC_Vars              ,ONLY: ChemReac, PartStateIntEn, PolyatomMolDSMC, VibQuantsPar, RadialWeighting, BGGas, ElecRelaxPart
@@ -414,7 +414,7 @@ SUBROUTINE DSMC_Chemistry(iPair, iReac)
         ChemReac%RecombParticle = 0
         ChemReac%nPairForRec = ChemReac%nPairForRec + 1
       ELSE
-        ! Utilize the second particle of the pair for the recombination        
+        ! Utilize the second particle of the pair for the recombination
         ChemReac%RecombParticle = Coll_pData(ChemReac%LastPairForRec)%iPart_p2
       END IF
     ELSE
@@ -522,7 +522,7 @@ SUBROUTINE DSMC_Chemistry(iPair, iReac)
       iPartIndx_NodeNewAmbi(newAmbiParts) = ReactInx(3)
     END IF
     IF (DSMC%ElectronicModel.EQ.4) THEN
-      newElecRelaxParts = newElecRelaxParts + 1 
+      newElecRelaxParts = newElecRelaxParts + 1
       iPartIndx_NodeNewElecRelax(newElecRelaxParts) = ReactInx(3)
     END IF
     Weight(3) = Weight(1)
@@ -562,7 +562,7 @@ SUBROUTINE DSMC_Chemistry(iPair, iReac)
     iPartIndx_NodeNewAmbi(newAmbiParts) = ReactInx(4)
   END IF
   IF (DSMC%ElectronicModel.EQ.4) THEN
-    newElecRelaxParts = newElecRelaxParts + 1 
+    newElecRelaxParts = newElecRelaxParts + 1
     iPartIndx_NodeNewElecRelax(newElecRelaxParts) = ReactInx(4)
   END IF
   Weight(4) = Weight(1)
@@ -595,6 +595,7 @@ SUBROUTINE DSMC_Chemistry(iPair, iReac)
                                                                 * PartState(4:6,ReactInx(3)) * Weight(3)
 #endif /* CODE_ANALYZE */
 
+! This might fall below zero
 Coll_pData(iPair)%Ec = 0.5 * MassRed * Coll_pData(iPair)%CRela2 + ChemReac%EForm(iReac)*SumWeightProd/REAL(NumProd)
 
 IF(RadialWeighting%DoRadialWeighting.OR.usevMPF) THEN
@@ -642,6 +643,24 @@ SUBROUTINE DSMC_Chemistry(iPair, iReac)
 IF (DSMC%ElectronicModel.GT.0) THEN
   Coll_pData(iPair)%Ec = Coll_pData(iPair)%Ec + PartStateIntEn(3,ReactInx(1))*Weight(1) + PartStateIntEn(3,ReactInx(2))*Weight(2)
 END IF
+! Sanity check
+IF(Coll_pData(iPair)%Ec.LT.0)THEN
+  IPWRITE(UNIT_StdOut,*) "DSMC Chemisry: Added vibrational and electronic energy"
+  IPWRITE(UNIT_StdOut,*) "iReac       =", iReac
+  IPWRITE(UNIT_StdOut,*) "ReactInx(1) =", ReactInx(1), " ReactInx(2) =", ReactInx(2)
+  IPWRITE(UNIT_StdOut,*) "Weight(1)   =", Weight(1)  , " Weight(2)   =", Weight(2)
+  IPWRITE(UNIT_StdOut,*) "Coll_pData(iPair)%Ec     = ", Coll_pData(iPair)%Ec
+  IPWRITE(UNIT_StdOut,*) "MassRed                  = ", MassRed
+  IPWRITE(UNIT_StdOut,*) "ChemReac%EForm(iReac)    =", ChemReac%EForm(iReac)
+  IPWRITE(UNIT_StdOut,*) "Coll_pData(iPair)%CRela2 = ", Coll_pData(iPair)%CRela2
+  IPWRITE(UNIT_StdOut,*) "PartStateIntEn(1:2,ReactInx(1)) =", PartStateIntEn(1:2,ReactInx(1))
+  IPWRITE(UNIT_StdOut,*) "PartStateIntEn(1:2,ReactInx(2)) =", PartStateIntEn(1:2,ReactInx(2))
+  IF(DSMC%ElectronicModel.GT.0)THEN
+    IPWRITE(UNIT_StdOut,*) "PartStateIntEn(3,ReactInx(1)) =", PartStateIntEn(3,ReactInx(1))
+    IPWRITE(UNIT_StdOut,*) "PartStateIntEn(3,ReactInx(2)) =", PartStateIntEn(3,ReactInx(2))
+  END IF ! DSMC%ElectronicModel.GT.0
+  CALL abort(__STAMP__,'Coll_pData(iPair)%Ec < 0.)',RealInfoOpt = Coll_pData(iPair)%Ec)
+END IF ! Coll_pData(iPair)%Ec.LT.0
 
 IF(EductReac(3).NE.0) THEN
   ! If a third collision partner exists (recombination/exchange reactions with defined third educt, A + B+ C), calculation of
@@ -721,8 +740,14 @@ SUBROUTINE DSMC_Chemistry(iPair, iReac)
 END IF
 FakXi = 0.5*Xi_total - 1.0
 !-------------------------------------------------------------------------------------------------------------------------------
-! Substracting the zero-point energy of the products (is added back later)
+! Subtracting the zero-point energy of the products (is added back later)
 Coll_pData(iPair)%Ec = Coll_pData(iPair)%Ec - SUM(EZeroTempToExec(:))
+! Sanity check
+IF(Coll_pData(iPair)%Ec.LT.0.)THEN
+  IPWRITE(UNIT_StdOut,*) "DSMC Chemisry: Zero-point energy"
+  IPWRITE(UNIT_StdOut,*) "iReac =", iReac
+  CALL abort(__STAMP__,'Coll_pData(iPair)%Ec < 0.)',RealInfoOpt = Coll_pData(iPair)%Ec)
+END IF ! Coll_pData(iPair)%Ec.LT.0.
 !-------------------------------------------------------------------------------------------------------------------------------
 ! Set the new species of the products
 DO iProd = 1, NumProd
@@ -739,14 +764,14 @@ SUBROUTINE DSMC_Chemistry(iPair, iReac)
       END IF
       PartStateIntEn(3,ReactInx(iProd)) = 0.0
       IF (DSMC%ElectronicModel.EQ.4) THEN
-        nElecRelaxChemParts = nElecRelaxChemParts + 1 
+        nElecRelaxChemParts = nElecRelaxChemParts + 1
         iPartIndx_NodeElecRelaxChem(nElecRelaxChemParts) = ReactInx(iProd)
         ElecRelaxPart(ReactInx(iProd)) = .FALSE.
       END IF
     ELSE
       IF (DSMC%ElectronicModel.EQ.4) THEN
         PartStateIntEn(3,ReactInx(iProd)) = 0.0
-        nElecRelaxChemParts = nElecRelaxChemParts + 1 
+        nElecRelaxChemParts = nElecRelaxChemParts + 1
         iPartIndx_NodeElecRelaxChem(nElecRelaxChemParts) = ReactInx(iProd)
         ElecRelaxPart(ReactInx(iProd)) = .FALSE.
       ELSE
@@ -759,6 +784,12 @@ SUBROUTINE DSMC_Chemistry(iPair, iReac)
         FakXi = FakXi - 0.5*Xi_elec(iProd)
         CALL ElectronicEnergyExchange(iPair,ReactInx(iProd),FakXi, NewPart = .TRUE., XSec_Level = 0)
         Coll_pData(iPair)%Ec = Coll_pData(iPair)%Ec - PartStateIntEn(3,ReactInx(iProd))*Weight(iProd)
+        IF(Coll_pData(iPair)%Ec.LT.0.) THEn
+          IPWRITE(UNIT_StdOut,*) "DSMC Chemisry: Electronic energy exchange"
+          IPWRITE(UNIT_StdOut,*) "iProd =", iProd
+          IPWRITE(UNIT_StdOut,*) "iReac =", iReac
+          CALL abort(__STAMP__,'Coll_pData(iPair)%Ec < 0.)',RealInfoOpt = Coll_pData(iPair)%Ec)
+        END IF
       END IF
     END IF
   END DO
@@ -780,6 +811,13 @@ SUBROUTINE DSMC_Chemistry(iPair, iReac)
       Coll_pData(iPair)%Ec = Coll_pData(iPair)%Ec + EZeroTempToExec(iProd)
       CALL DSMC_VibRelaxDiatomic(iPair,ReactInx(iProd),FakXi)
       Coll_pData(iPair)%Ec = Coll_pData(iPair)%Ec - PartStateIntEn(1,ReactInx(iProd))*Weight(iProd)
+      ! Sanity check
+      IF(Coll_pData(iPair)%Ec.LT.0.) THEn
+        IPWRITE(UNIT_StdOut,*) "DSMC Chemisry: Vibrational energy exchange"
+        IPWRITE(UNIT_StdOut,*) "iProd =", iProd
+        IPWRITE(UNIT_StdOut,*) "iReac =", iReac
+        CALL abort(__STAMP__,'Coll_pData(iPair)%Ec < 0.)',RealInfoOpt = Coll_pData(iPair)%Ec)
+      END IF
     END IF
   END IF
 END DO
@@ -797,6 +835,13 @@ SUBROUTINE DSMC_Chemistry(iPair, iReac)
       FakXi = FakXi - 0.5*SpecDSMC(ProductReac(iProd))%Xi_Rot
     END IF
     Coll_pData(iPair)%Ec = Coll_pData(iPair)%Ec - PartStateIntEn(2,ReactInx(iProd))
+    ! Sanity check
+    IF(Coll_pData(iPair)%Ec.LT.0.) THEn
+      IPWRITE(UNIT_StdOut,*) "DSMC Chemisry: Rotational energy exchange"
+      IPWRITE(UNIT_StdOut,*) "iProd =", iProd
+      IPWRITE(UNIT_StdOut,*) "iReac =", iReac
+      CALL abort(__STAMP__,'Coll_pData(iPair)%Ec < 0.)',RealInfoOpt = Coll_pData(iPair)%Ec)
+    END IF
     PartStateIntEn(2,ReactInx(iProd)) = PartStateIntEn(2,ReactInx(iProd))/Weight(iProd)
   ELSE
     PartStateIntEn(1,ReactInx(iProd)) = 0.0
@@ -984,6 +1029,12 @@ SUBROUTINE DSMC_Chemistry(iPair, iReac)
   END IF
   ERel_React1_React3 = Coll_pData(iPair)%Ec
 
+  ! Sanity check
+  IF(ERel_React1_React3.LT.0.) THEn
+    IPWRITE(UNIT_StdOut,*) "iReac =", iReac
+    CALL abort(__STAMP__,'DSMC Chemisry: Coll_pData(iPair)%Ec < 0.)',RealInfoOpt = Coll_pData(iPair)%Ec)
+  END IF
+
   IF (RadialWeighting%DoRadialWeighting.OR.usevMPF) THEN
     FracMassCent1 = Species(ProductReac(1))%MassIC *Weight(1) &
         /(Species(ProductReac(1))%MassIC* Weight(1) + Species(ProductReac(2))%MassIC * Weight(2))
@@ -1768,9 +1819,9 @@ SUBROUTINE PhotoIonization_InsertProducts(iPair, iReac, iInit, InitSpec, iLineOp
   cRelaNew(1:3) = PostCollVec(iPair)
 
   !deltaV particle 1
-  PartState(4:6,ReactInx(1)) = VeloCOM(1:3) + FracMassCent2*cRelaNew(1:3)  
+  PartState(4:6,ReactInx(1)) = VeloCOM(1:3) + FracMassCent2*cRelaNew(1:3)
   !deltaV particle 2
-  PartState(4:6,ReactInx(2)) = VeloCOM(1:3) - FracMassCent1*cRelaNew(1:3)  
+  PartState(4:6,ReactInx(2)) = VeloCOM(1:3) - FracMassCent1*cRelaNew(1:3)
 END IF
 
 IF(CalcPartBalance) THEN
@@ -1803,7 +1854,7 @@ PPURE FUNCTION GetRandomVectorInPlane(b1,b2,VeloVec,RandVal)
 REAL               :: GetRandomVectorInPlane(1:3) ! Output velocity vector
 !-----------------------------------------------------------------------------------------------------------------------------------
 ! LOCAL VARIABLES
-REAL               :: Vabs ! Absolute velocity 
+REAL               :: Vabs ! Absolute velocity
 REAL               :: phi ! random angle between 0 and 2*PI
 !===================================================================================================================================
 Vabs = VECNORM(VeloVec)

src/particles/dsmc/dsmc_init.f90

@@ -885,7 +885,7 @@ SUBROUTINE InitDSMC()
     IF (useRelaxProbCorrFactor.AND.(DSMC%ElectronicModel.EQ.1)) THEN
       DO iSpec = 1, nSpecies
         ALLOCATE(SpecDSMC(iSpec)%ElecRelaxCorrectFac(nSpecies))
-      END DO 
+      END DO
     END IF
   END IF
 
@@ -936,7 +936,7 @@ SUBROUTINE InitDSMC()
 !-----------------------------------------------------------------------------------------------------------------------------------
 ! Calculate vib collision numbers and characteristic velocity, according to Abe
 !-----------------------------------------------------------------------------------------------------------------------------------
-  ! (i) dref changed from dref = 0.5 * (dref_1+dref_2) 
+  ! (i) dref changed from dref = 0.5 * (dref_1+dref_2)
   !                  to   dref(iSpec,jSpec) which is identical to old definition (for averagedCollisionParameters=TRUE (DEFAULT))
   ! in case of averagedCollisionParameter=FALSE dref(iSpec,jSpec) contains collision specific dref see --help for details
   IF((DSMC%VibRelaxProb.EQ.2).AND.(CollisMode.GE.2)) THEN
@@ -1370,6 +1370,7 @@ SUBROUTINE FinalizeDSMC()
 SDEALLOCATE(DSMC%QualityFacSampVibSamp)
 SDEALLOCATE(DSMC%CalcVibProb)
 SDEALLOCATE(DSMC%CalcRotProb)
+SDEALLOCATE(DSMC%InstantTXiElec)
 SDEALLOCATE(SampDSMC)
 SDEALLOCATE(PartStateIntEn)
 SDEALLOCATE(ElecRelaxPart)
@@ -1517,7 +1518,7 @@ RECURSIVE SUBROUTINE DeleteNodeVolume(Node)
 !-----------------------------------------------------------------------------------------------------------------------------------
 ! LOCAL VARIABLES
 !===================================================================================================================================
-IF(ASSOCIATED(Node%SubNode1)) THEN 
+IF(ASSOCIATED(Node%SubNode1)) THEN
   CALL DeleteNodeVolume(Node%SubNode1)
   DEALLOCATE(Node%SubNode1)
 END IF

src/particles/mcc/mcc_xsec.f90

@@ -499,7 +499,7 @@ PPURE REAL FUNCTION InterpolateCrossSection_Elec(iCase,iLevel,CollisionEnergy)
 ! INPUT VARIABLES
 !-----------------------------------------------------------------------------------------------------------------------------------
 INTEGER,INTENT(IN)            :: iCase                            !< Case index
-INTEGER,INTENT(IN)            :: iLevel                           !< 
+INTEGER,INTENT(IN)            :: iLevel                           !<
 REAL,INTENT(IN)               :: CollisionEnergy                  !< Collision energy in [J]
 !-----------------------------------------------------------------------------------------------------------------------------------
 ! LOCAL VARIABLES
@@ -571,7 +571,7 @@ PPURE REAL FUNCTION InterpolateVibRelaxProb(iCase,CollisionEnergy)
 InterpolateVibRelaxProb = 0.
 MaxDOF = SIZE(SpecXSec(iCase)%VibXSecData,2)
 
-IF(CollisionEnergy.GT.SpecXSec(iCase)%VibXSecData(1,MaxDOF)) THEN 
+IF(CollisionEnergy.GT.SpecXSec(iCase)%VibXSecData(1,MaxDOF)) THEN
   ! If the collision energy is greater than the maximal value, get the cross-section of the last level and leave routine
   InterpolateVibRelaxProb = SpecXSec(iCase)%VibXSecData(2,MaxDOF)
   ! Leave routine
@@ -1377,6 +1377,7 @@ SUBROUTINE XSec_CalcReactionProb(iPair,iCase,iElem,SpecNum1,SpecNum2,MacroPartic
 !===================================================================================================================================
 ! MODULES
 USE MOD_Globals_Vars
+USE MOD_Globals               ,ONLY: abort
 USE MOD_DSMC_Vars             ,ONLY: SpecDSMC, Coll_pData, CollInf, BGGas, ChemReac, RadialWeighting, DSMC, PartStateIntEn
 USE MOD_MCC_Vars              ,ONLY: SpecXSec
 USE MOD_Particle_Vars         ,ONLY: PartSpecies, Species, VarTimeStep, usevMPF
@@ -1392,7 +1393,7 @@ SUBROUTINE XSec_CalcReactionProb(iPair,iCase,iElem,SpecNum1,SpecNum2,MacroPartic
 INTEGER                       :: iPath, ReacTest, EductReac(1:3), ProductReac(1:4), ReactInx(1:2), nPair, iProd
 INTEGER                       :: NumWeightProd, targetSpec, bgSpec
 REAL                          :: EZeroPoint_Prod, dtCell, Weight(1:4), ReducedMass, ReducedMassUnweighted, CollEnergy, GammaFac
-REAL                          :: EZeroPoint_Educt, SpecNumTarget, SpecNumSource, CrossSection
+REAL                          :: EZeroPoint_Educt, SpecNumTarget, SpecNumSource, CrossSection, EcRelativistic
 !===================================================================================================================================
 Weight = 0.; ReactInx = 0
 nPair = SIZE(Coll_pData)
@@ -1452,14 +1453,8 @@ SUBROUTINE XSec_CalcReactionProb(iPair,iCase,iElem,SpecNum1,SpecNum2,MacroPartic
       dtCell = dt
     END IF
 
-    IF(Coll_pData(iPair)%cRela2 .LT. RelativisticLimit) THEN
-      Coll_pData(iPair)%Ec = 0.5 * ReducedMass * Coll_pData(iPair)%cRela2
-    ELSE
-      ! Relativistic treatment under the assumption that the velocity of the background species is zero or negligible
-      GammaFac = Coll_pData(iPair)%cRela2*c2_inv
-      GammaFac = 1./SQRT(1.-GammaFac)
-      Coll_pData(iPair)%Ec = (GammaFac-1.) * ReducedMass * c2
-    END IF
+    ! No relativistic treatment here until the chemistry is done fully considering relativistic effects
+    Coll_pData(iPair)%Ec = 0.5 * ReducedMass * Coll_pData(iPair)%cRela2
 
     Coll_pData(iPair)%Ec = Coll_pData(iPair)%Ec + (PartStateIntEn(1,ReactInx(1)) + PartStateIntEn(2,ReactInx(1))) * Weight(1) &
                                                 + (PartStateIntEn(1,ReactInx(2)) + PartStateIntEn(2,ReactInx(2))) * Weight(2)
@@ -1469,8 +1464,30 @@ SUBROUTINE XSec_CalcReactionProb(iPair,iCase,iElem,SpecNum1,SpecNum2,MacroPartic
     END IF
     ! Check first if sufficient energy is available for the products after the reaction
     ASSOCIATE( ReactionProb => ChemReac%CollCaseInfo(iCase)%ReactionProb(iPath) )
-      IF(((Coll_pData(iPair)%Ec-EZeroPoint_Prod).GE.(-ChemReac%EForm(ReacTest)*SUM(Weight)/NumWeightProd))) THEN
-        CollEnergy = (Coll_pData(iPair)%Ec-EZeroPoint_Educt) * 2./(Weight(1)+Weight(2))
+      IF(((Coll_pData(iPair)%Ec-EZeroPoint_Prod).GE.(-ChemReac%EForm(ReacTest)*SUM(Weight)/REAL(NumWeightProd)))) THEN
+
+        ! Check relativistic limit for the interpolation of the XSec data but not for in inquiry above
+        IF(Coll_pData(iPair)%cRela2 .LT. RelativisticLimit) THEN
+          CollEnergy = (Coll_pData(iPair)%Ec-EZeroPoint_Educt) * 2./(Weight(1)+Weight(2))
+        ELSE
+          ! Relativistic treatment under the assumption that the velocity of the background species is zero or negligible
+          GammaFac = Coll_pData(iPair)%cRela2*c2_inv
+          ! Sanity check
+          IF(GammaFac.GE.1.0)THEN
+            CALL abort(__STAMP__,'X = Coll_pData(iPair)%cRela2*c2_inv >= 1.0, which results in SQRT(1-X) failing.')
+          ELSE
+            GammaFac = 1./SQRT(1.-GammaFac)
+          END IF
+          EcRelativistic = (GammaFac-1.) * ReducedMass * c2
+
+          EcRelativistic = EcRelativistic + (PartStateIntEn(1,ReactInx(1)) + PartStateIntEn(2,ReactInx(1))) * Weight(1) &
+                                                      + (PartStateIntEn(1,ReactInx(2)) + PartStateIntEn(2,ReactInx(2))) * Weight(2)
+          IF (DSMC%ElectronicModel.GT.0) EcRelativistic = EcRelativistic + PartStateIntEn(3,ReactInx(1))*Weight(1) &
+                                                                         + PartStateIntEn(3,ReactInx(2))*Weight(2)
+
+          CollEnergy = (EcRelativistic-EZeroPoint_Educt) * 2./(Weight(1)+Weight(2))
+        END IF
+
         CrossSection = InterpolateCrossSection_Chem(iCase,iPath,CollEnergy)
         IF(SpecXSec(iCase)%UseCollXSec) THEN
           ! Interpolate the reaction cross-section