From cdac034fef81060cad4c3fa64db38917dbb0e595 Mon Sep 17 00:00:00 2001 From: Stephen Copplestone Date: Thu, 27 Oct 2022 23:54:14 +0200 Subject: [PATCH 1/2] Fixed bug in chemistry. Inconsistent consideration of relativistic effects in the chemistry module can lead to negative collision energies, which in turn create NaN in the particle velocity that are propagated to the field solver that finally crashes. Only consider relativistic effects during the interpolation of the XSec data for now. Also added more sanity checks during chemistry. Fixed deallocation of DSMC%InstantTXiElec. --- .../dsmc/dsmc_chemical_reactions.f90 | 93 +++++++++++++++++-- src/particles/dsmc/dsmc_init.f90 | 7 +- src/particles/mcc/mcc_xsec.f90 | 43 ++++++--- 3 files changed, 117 insertions(+), 26 deletions(-) diff --git a/src/particles/dsmc/dsmc_chemical_reactions.f90 b/src/particles/dsmc/dsmc_chemical_reactions.f90 index 982f53146..0a5819802 100644 --- a/src/particles/dsmc/dsmc_chemical_reactions.f90 +++ b/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)) @@ -1096,6 +1147,28 @@ SUBROUTINE DSMC_Chemistry(iPair, iReac) END DO #endif /* CODE_ANALYZE */ +associate( iPart1 => ReactInx(1), iPart2 => ReactInx(2)) + IF(ANY(ISNAN(PartState(:,iPart1))))THEN + IPWRITE(UNIT_StdOut,*) "after eight " + IPWRITE(UNIT_StdOut,*) "iPart1,PartState(:,iPart1) =", iPart1,PartState(:,iPart1) + CALL abort(__STAMP__,'after 8') + END IF ! ANY(ISNAN(PartState(:,iPart1))) + IF(ANY(ISNAN(PartState(:,iPart2))))THEN + IPWRITE(UNIT_StdOut,*) "after eight " + IPWRITE(UNIT_StdOut,*) "iPart2,PartState(:,iPart2) =", iPart2,PartState(:,iPart2) + CALL abort(__STAMP__,'after 8') + END IF ! ANY(ISNAN(PartState(:,iPart2))) +END ASSOCIATE + +IF(ProductReac(3).NE.0)THEN +IF(ANY(ISNAN(PartState(:,ReactInx(3)))))THEN + IPWRITE(UNIT_StdOut,*) "AFTER " + IPWRITE(UNIT_StdOut,*) "ReactInx(3),PartState(:,ReactInx(3)) =", ReactInx(3),PartState(:,ReactInx(3)) + CALL abort(__STAMP__,'after 3') +END IF ! ANY(ISNAN(PartState(:,ReactInx(3)))) + +END IF ! ProductReac(3).EQ.0 + END SUBROUTINE DSMC_Chemistry @@ -1768,9 +1841,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 +1876,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) diff --git a/src/particles/dsmc/dsmc_init.f90 b/src/particles/dsmc/dsmc_init.f90 index 0ae8f82ad..54f5d94ab 100644 --- a/src/particles/dsmc/dsmc_init.f90 +++ b/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 diff --git a/src/particles/mcc/mcc_xsec.f90 b/src/particles/mcc/mcc_xsec.f90 index 090446112..57f071a62 100644 --- a/src/particles/mcc/mcc_xsec.f90 +++ b/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 From dc8d3742012ca3dd504ddfe93efab7fd6b4fc437 Mon Sep 17 00:00:00 2001 From: Stephen Copplestone Date: Mon, 31 Oct 2022 14:19:56 +0100 Subject: [PATCH 2/2] Removed debugging output --- .../dsmc/dsmc_chemical_reactions.f90 | 22 ------------------- 1 file changed, 22 deletions(-) diff --git a/src/particles/dsmc/dsmc_chemical_reactions.f90 b/src/particles/dsmc/dsmc_chemical_reactions.f90 index 0a5819802..9024048a7 100644 --- a/src/particles/dsmc/dsmc_chemical_reactions.f90 +++ b/src/particles/dsmc/dsmc_chemical_reactions.f90 @@ -1147,28 +1147,6 @@ SUBROUTINE DSMC_Chemistry(iPair, iReac) END DO #endif /* CODE_ANALYZE */ -associate( iPart1 => ReactInx(1), iPart2 => ReactInx(2)) - IF(ANY(ISNAN(PartState(:,iPart1))))THEN - IPWRITE(UNIT_StdOut,*) "after eight " - IPWRITE(UNIT_StdOut,*) "iPart1,PartState(:,iPart1) =", iPart1,PartState(:,iPart1) - CALL abort(__STAMP__,'after 8') - END IF ! ANY(ISNAN(PartState(:,iPart1))) - IF(ANY(ISNAN(PartState(:,iPart2))))THEN - IPWRITE(UNIT_StdOut,*) "after eight " - IPWRITE(UNIT_StdOut,*) "iPart2,PartState(:,iPart2) =", iPart2,PartState(:,iPart2) - CALL abort(__STAMP__,'after 8') - END IF ! ANY(ISNAN(PartState(:,iPart2))) -END ASSOCIATE - -IF(ProductReac(3).NE.0)THEN -IF(ANY(ISNAN(PartState(:,ReactInx(3)))))THEN - IPWRITE(UNIT_StdOut,*) "AFTER " - IPWRITE(UNIT_StdOut,*) "ReactInx(3),PartState(:,ReactInx(3)) =", ReactInx(3),PartState(:,ReactInx(3)) - CALL abort(__STAMP__,'after 3') -END IF ! ANY(ISNAN(PartState(:,ReactInx(3)))) - -END IF ! ProductReac(3).EQ.0 - END SUBROUTINE DSMC_Chemistry