chemistryModel: rewrite main solver loop and update chemistrySolvers

accordingly to reuse the estimated sub-time-step more effectively

src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C

@@ -514,7 +514,7 @@ void Foam::chemistryModel<CompType, ThermoType>::jacobian
     }
 
     // calculate the dcdT elements numerically
-    const scalar delta = 1.0e-8;
+    const scalar delta = 1.0e-3;
     const scalarField dcdT0(omega(c2, T - delta, p));
     const scalarField dcdT1(omega(c2, T + delta, p));
 
@@ -772,9 +772,6 @@ Foam::scalar Foam::chemistryModel<CompType, ThermoType>::solve
         this->thermo().rho()
     );
 
-    tmp<volScalarField> thc = this->thermo().hc();
-    const scalarField& hc = thc();
-    const scalarField& he = this->thermo().he();
     const scalarField& T = this->thermo().T();
     const scalarField& p = this->thermo().p();
 
@@ -784,8 +781,7 @@ Foam::scalar Foam::chemistryModel<CompType, ThermoType>::solve
     forAll(rho, celli)
     {
         const scalar rhoi = rho[celli];
-        const scalar hi = he[celli] + hc[celli];
+        scalar pi = p[celli];
-        const scalar pi = p[celli];
         scalar Ti = T[celli];
 
         for (label i=0; i<nSpecie_; i++)
@@ -794,32 +790,18 @@ Foam::scalar Foam::chemistryModel<CompType, ThermoType>::solve
             c0[i] = c[i];
         }
 
-        // initialise timing parameters
+        // Initialise time progress
-        scalar t = 0;
         scalar timeLeft = deltaT[celli];
-        scalar tauC = this->deltaTChem_[celli];
-        scalar dt = min(timeLeft, tauC);
 
-        // calculate the chemical source terms
+        // Calculate the chemical source terms
         while (timeLeft > SMALL)
         {
-            tauC = this->solve(c, Ti, pi, t, dt);
+            scalar dt = timeLeft;
-            t += dt;
+            this->solve(c, Ti, pi, dt, this->deltaTChem_[celli]);
-
-            // update the temperature
-            const scalar cTot = sum(c);
-            ThermoType mixture((c[0]/cTot)*specieThermo_[0]);
-            for (label i=1; i<nSpecie_; i++)
-            {
-                mixture += (c[i]/cTot)*specieThermo_[i];
-            }
-            Ti = mixture.THa(hi, pi, Ti);
-
             timeLeft -= dt;
-            this->deltaTChem_[celli] = tauC;
-            dt = max(SMALL, min(timeLeft, tauC));
         }
-        deltaTMin = min(tauC, deltaTMin);
+
+        deltaTMin = min(this->deltaTChem_[celli], deltaTMin);
 
         for (label i=0; i<nSpecie_; i++)
         {
@@ -857,13 +839,13 @@ Foam::scalar Foam::chemistryModel<CompType, ThermoType>::solve
 
 
 template<class CompType, class ThermoType>
-Foam::scalar Foam::chemistryModel<CompType, ThermoType>::solve
+void Foam::chemistryModel<CompType, ThermoType>::solve
 (
     scalarField &c,
-    const scalar T,
+    scalar& T,
-    const scalar p,
+    scalar& p,
-    const scalar t0,
+    scalar& deltaT,
-    const scalar dt
+    scalar& subDeltaT
 ) const
 {
     notImplemented
@@ -871,14 +853,12 @@ Foam::scalar Foam::chemistryModel<CompType, ThermoType>::solve
         "chemistryModel::solve"
         "("
             "scalarField&, "
-            "const scalar, "
+            "scalar&, "
-            "const scalar, "
+            "scalar&, "
-            "const scalar, "
+            "scalar&, "
-            "const scalar"
+            "scalar&"
         ") const"
     );
-
-    return (0);
 }
 
 

src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H

@@ -78,6 +78,8 @@ class chemistryModel
 
 protected:
 
+    typedef ThermoType thermoType;
+
     // Private data
 
         //- Reference to the field of specie mass fractions
@@ -260,13 +262,13 @@ public:
                 scalarSquareMatrix& dfdc
             ) const;
 
-            virtual scalar solve
+            virtual void solve
             (
                 scalarField &c,
-                const scalar T,
+                scalar& T,
-                const scalar p,
+                scalar& p,
-                const scalar t0,
+                scalar& deltaT,
-                const scalar dt
+                scalar& subDeltaT
             ) const;
 };
 

src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C

@@ -38,7 +38,8 @@ Foam::EulerImplicit<ChemistryModel>::EulerImplicit
     chemistrySolver<ChemistryModel>(mesh),
     coeffsDict_(this->subDict("EulerImplicitCoeffs")),
     cTauChem_(readScalar(coeffsDict_.lookup("cTauChem"))),
-    eqRateLimiter_(coeffsDict_.lookup("equilibriumRateLimiter"))
+    eqRateLimiter_(coeffsDict_.lookup("equilibriumRateLimiter")),
+    cTp_(this->nEqns())
 {}
 
 
@@ -52,17 +53,16 @@ Foam::EulerImplicit<ChemistryModel>::~EulerImplicit()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class ChemistryModel>
-Foam::scalar Foam::EulerImplicit<ChemistryModel>::solve
+void Foam::EulerImplicit<ChemistryModel>::solve
 (
     scalarField& c,
-    const scalar T,
+    scalar& T,
-    const scalar p,
+    scalar& p,
-    const scalar t0,
+    scalar& deltaT,
-    const scalar dt
+    scalar& subDeltaT
 ) const
 {
-    scalar pf, cf, pr, cr;
+    deltaT = min(deltaT, subDeltaT);
-    label lRef, rRef;
 
     const label nSpecie = this->nSpecie();
     simpleMatrix<scalar> RR(nSpecie, 0, 0);
@@ -72,13 +72,28 @@ Foam::scalar Foam::EulerImplicit<ChemistryModel>::solve
         c[i] = max(0.0, c[i]);
     }
 
+    // Calculate the absolute enthalpy
+    scalar cTot = sum(c);
+    typename ChemistryModel::thermoType mixture
+    (
+        (c[0]/cTot)*this->specieThermo_[0]
+    );
+    for (label i=1; i<nSpecie; i++)
+    {
+        mixture += (c[i]/cTot)*this->specieThermo_[i];
+    }
+    scalar ha = mixture.Ha(p, T);
+
     for (label i=0; i<nSpecie; i++)
     {
-        RR.source()[i] = c[i]/dt;
+        RR.source()[i] = c[i]/deltaT;
     }
 
     forAll(this->reactions(), i)
     {
+        scalar pf, cf, pr, cr;
+        label lRef, rRef;
+
         scalar omegai = this->omegaI(i, c, T, p, pf, cf, lRef, pr, cr, rRef);
 
         scalar corr = 1.0;
@@ -86,11 +101,11 @@ Foam::scalar Foam::EulerImplicit<ChemistryModel>::solve
         {
             if (omegai < 0.0)
             {
-                corr = 1.0/(1.0 + pr*dt);
+                corr = 1.0/(1.0 + pr*deltaT);
             }
             else
             {
-                corr = 1.0/(1.0 + pf*dt);
+                corr = 1.0/(1.0 + pf*deltaT);
             }
         }
 
@@ -100,7 +115,7 @@ Foam::scalar Foam::EulerImplicit<ChemistryModel>::solve
 
     for (label i=0; i<nSpecie; i++)
     {
-        RR[i][i] += 1.0/dt;
+        RR[i][i] += 1.0/deltaT;
     }
 
     c = RR.LUsolve();
@@ -109,26 +124,35 @@ Foam::scalar Foam::EulerImplicit<ChemistryModel>::solve
         c[i] = max(0.0, c[i]);
     }
 
-    // estimate the next time step
+    // Update the temperature
+    cTot = sum(c);
+    mixture = (c[0]/cTot)*this->specieThermo_[0];
+    for (label i=1; i<nSpecie; i++)
+    {
+        mixture += (c[i]/cTot)*this->specieThermo_[i];
+    }
+    T = mixture.THa(ha, p, T);
+
+    // Estimate the next time step
     scalar tMin = GREAT;
     const label nEqns = this->nEqns();
-    scalarField c1(nEqns, 0.0);
 
     for (label i=0; i<nSpecie; i++)
     {
-        c1[i] = c[i];
+        cTp_[i] = c[i];
     }
-    c1[nSpecie] = T;
+    cTp_[nSpecie] = T;
-    c1[nSpecie+1] = p;
+    cTp_[nSpecie+1] = p;
 
     scalarField dcdt(nEqns, 0.0);
-    this->derivatives(0.0, c1, dcdt);
+    this->derivatives(0.0, cTp_, dcdt);
 
     const scalar sumC = sum(c);
 
     for (label i=0; i<nSpecie; i++)
     {
         scalar d = dcdt[i];
+
         if (d < -SMALL)
         {
             tMin = min(tMin, -(c[i] + SMALL)/d);
@@ -141,7 +165,7 @@ Foam::scalar Foam::EulerImplicit<ChemistryModel>::solve
         }
     }
 
-    return cTauChem_*tMin;
+    subDeltaT = cTauChem_*tMin;
 }
 
 

src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.H

@@ -65,6 +65,9 @@ class EulerImplicit
             //- Equilibrium rate limiter flag (on/off)
             Switch eqRateLimiter_;
 
+        // Solver data
+        mutable scalarField cTp_;
+
 
 public:
 
@@ -85,13 +88,13 @@ public:
     // Member Functions
 
         //- Update the concentrations and return the chemical time
-        virtual scalar solve
+        virtual void solve
         (
             scalarField& c,
-            const scalar T,
+            scalar& T,
-            const scalar p,
+            scalar& p,
-            const scalar t0,
+            scalar& deltaT,
-            const scalar dt
+            scalar& subDeltaT
         ) const;
 };
 

src/thermophysicalModels/chemistryModel/chemistrySolver/chemistrySolver/chemistrySolver.H

@@ -69,13 +69,13 @@ public:
     // Member Functions
 
         //- Update the concentrations and return the chemical time
-        virtual scalar solve
+        virtual void solve
         (
             scalarField &c,
-            const scalar T,
+            scalar& T,
-            const scalar p,
+            scalar& p,
-            const scalar t0,
+            scalar& deltaT,
-            const scalar dt
+            scalar& subDeltaT
         ) const = 0;
 };
 

src/thermophysicalModels/chemistryModel/chemistrySolver/chemistrySolver/makeChemistrySolverTypes.H

@@ -31,9 +31,9 @@ License
 #include "chemistryModel.H"
 
 #include "noChemistrySolver.H"
+#include "sequential.H"
 #include "EulerImplicit.H"
 #include "ode.H"
-#include "sequential.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.C

@@ -48,17 +48,15 @@ Foam::noChemistrySolver<ChemistryModel>::~noChemistrySolver()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class ChemistryModel>
-Foam::scalar Foam::noChemistrySolver<ChemistryModel>::solve
+void Foam::noChemistrySolver<ChemistryModel>::solve
 (
     scalarField&,
-    const scalar,
+    scalar&,
-    const scalar,
+    scalar&,
-    const scalar,
+    scalar&,
-    const scalar
+    scalar&
 ) const
-{
+{}
-    return GREAT;
-}
 
 
 // ************************************************************************* //

src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.H

@@ -72,13 +72,13 @@ public:
     // Member Functions
 
         //- Update the concentrations and return the chemical time
-        virtual scalar solve
+        virtual void solve
         (
             scalarField& c,
-            const scalar T,
+            scalar& T,
-            const scalar p,
+            scalar& p,
-            const scalar t0,
+            scalar& deltaT,
-            const scalar dt
+            scalar& subDeltaT
         ) const;
 };
 

src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C

@@ -38,7 +38,8 @@ Foam::ode<ChemistryModel>::ode
     coeffsDict_(this->subDict("odeCoeffs")),
     solverName_(coeffsDict_.lookup("solver")),
     odeSolver_(ODESolver::New(solverName_, *this)),
-    eps_(readScalar(coeffsDict_.lookup("eps")))
+    eps_(readScalar(coeffsDict_.lookup("eps"))),
+    cTp_(this->nEqns())
 {}
 
 
@@ -52,44 +53,41 @@ Foam::ode<ChemistryModel>::~ode()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class ChemistryModel>
-Foam::scalar Foam::ode<ChemistryModel>::solve
+void Foam::ode<ChemistryModel>::solve
 (
     scalarField& c,
-    const scalar T,
+    scalar& T,
-    const scalar p,
+    scalar& p,
-    const scalar t0,
+    scalar& deltaT,
-    const scalar dt
+    scalar& subDeltaT
 ) const
 {
     label nSpecie = this->nSpecie();
-    scalarField c1(this->nEqns(), 0.0);
 
-    // copy the concentration, T and P to the total solve-vector
+    // Copy the concentration, T and P to the total solve-vector
-    for (label i = 0; i < nSpecie; i++)
+    for (register int i=0; i<nSpecie; i++)
     {
-        c1[i] = c[i];
+        cTp_[i] = c[i];
     }
-    c1[nSpecie] = T;
+    cTp_[nSpecie] = T;
-    c1[nSpecie+1] = p;
+    cTp_[nSpecie+1] = p;
-
-    scalar dtEst = dt;
 
     odeSolver_->solve
     (
         *this,
-        t0,
+        0,
-        t0 + dt,
+        deltaT,
-        c1,
+        cTp_,
         eps_,
-        dtEst
+        subDeltaT
     );
 
-    forAll(c, i)
+    for (register int i=0; i<nSpecie; i++)
     {
-        c[i] = max(0.0, c1[i]);
+        c[i] = max(0.0, cTp_[i]);
     }
-
+    T = cTp_[nSpecie];
-    return dtEst;
+    p = cTp_[nSpecie+1];
 }
 
 

src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.H

@@ -62,6 +62,9 @@ class ode
 
             scalar eps_;
 
+        // Solver data
+        mutable scalarField cTp_;
+
 
 public:
 
@@ -81,13 +84,14 @@ public:
 
     // Member Functions
 
-        virtual scalar solve
+        //- Update the concentrations and return the chemical time
+        virtual void solve
         (
             scalarField& c,
-            const scalar T,
+            scalar& T,
-            const scalar p,
+            scalar& p,
-            const scalar t0,
+            scalar& deltaT,
-            const scalar dt
+            scalar& subDeltaT
         ) const;
 };
 

src/thermophysicalModels/chemistryModel/chemistrySolver/sequential/sequential.C

@@ -51,42 +51,67 @@ Foam::sequential<ChemistryModel>::~sequential()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class ChemistryModel>
-Foam::scalar Foam::sequential<ChemistryModel>::solve
+void Foam::sequential<ChemistryModel>::solve
 (
     scalarField& c,
-    const scalar T,
+    scalar& T,
-    const scalar p,
+    scalar& p,
-    const scalar t0,
+    scalar& deltaT,
-    const scalar dt
+    scalar& subDeltaT
 ) const
 {
-    scalar tChemInv = SMALL;
+    deltaT = min(deltaT, subDeltaT);
 
-    scalar pf, cf, pb, cb;
+    const label nSpecie = this->nSpecie();
-    label lRef, rRef;
+
+    // Calculate the absolute enthalpy
+    scalar cTot = sum(c);
+    typename ChemistryModel::thermoType mixture
+    (
+        (c[0]/cTot)*this->specieThermo_[0]
+    );
+    for (label i=1; i<nSpecie; i++)
+    {
+        mixture += (c[i]/cTot)*this->specieThermo_[i];
+    }
+    scalar ha = mixture.Ha(p, T);
+
+    scalar tChemInv = SMALL;
 
     forAll(this->reactions(), i)
     {
+        scalar pf, cf, pb, cb;
+        label lRef, rRef;
+
         scalar omega = this->omegaI(i, c, T, p, pf, cf, lRef, pb, cb, rRef);
 
         if (eqRateLimiter_)
         {
             if (omega < 0.0)
             {
-                omega /= 1.0 + pb*dt;
+                omega /= 1.0 + pb*deltaT;
             }
             else
             {
-                omega /= 1.0 + pf*dt;
+                omega /= 1.0 + pf*deltaT;
             }
         }
 
         tChemInv = max(tChemInv, mag(omega));
 
-        this->updateConcsInReactionI(i, dt, omega, p, T, c);
+        this->updateConcsInReactionI(i, deltaT, omega, p, T, c);
     }
 
-    return cTauChem_/tChemInv;
+    // Update the temperature
+    cTot = sum(c);
+    mixture = (c[0]/cTot)*this->specieThermo_[0];
+    for (label i=1; i<nSpecie; i++)
+    {
+        mixture += (c[i]/cTot)*this->specieThermo_[i];
+    }
+    T = mixture.THa(ha, p, T);
+
+    subDeltaT = cTauChem_/tChemInv;
 }
 
 

src/thermophysicalModels/chemistryModel/chemistrySolver/sequential/sequential.H

@@ -28,9 +28,7 @@ Description
     Foam::sequential
 
 SourceFiles
-    sequentialI.H
     sequential.C
-    sequentialIO.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -86,13 +84,13 @@ public:
     // Member Functions
 
         //- Update the concentrations and return the chemical time
-        virtual scalar solve
+        virtual void solve
         (
             scalarField& c,
-            const scalar T,
+            scalar& T,
-            const scalar p,
+            scalar& p,
-            const scalar t0,
+            scalar& deltaT,
-            const scalar dt
+            scalar& subDeltaT
         ) const;
 };
 

src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolverType.H

@@ -34,8 +34,8 @@ Description
 
 #include "noChemistrySolver.H"
 #include "EulerImplicit.H"
-#include "ode.H"
 #include "sequential.H"
+#include "ode.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

tutorials/combustion/LTSReactingFoam/ras/counterFlowFlame2D/constant/chemistryProperties

@@ -40,7 +40,6 @@ odeCoeffs
 {
     solver          KRR4;
     eps             0.05;
-    scale           1;
 }
 
 // ************************************************************************* //

tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/constant/chemistryProperties

@@ -17,7 +17,7 @@ FoamFile
 
 chemistryType
 {
-    chemistrySolver   ode;
+    chemistrySolver   EulerImplicit;
     chemistryThermo   psi;
 }
 
@@ -27,12 +27,12 @@ initialChemicalTimeStep 1e-07;
 
 sequentialCoeffs
 {
-    cTauChem        0.001;
+    cTauChem        0.05;
 }
 
 EulerImplicitCoeffs
 {
-    cTauChem        0.05;
+    cTauChem        1;
     equilibriumRateLimiter off;
 }
 
@@ -40,7 +40,6 @@ odeCoeffs
 {
     solver          KRR4;
     eps             0.05;
-    scale           1;
 }
 
 // ************************************************************************* //

tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/chemistryProperties

@@ -41,7 +41,6 @@ odeCoeffs
 {
     solver          SIBS;
     eps             0.05;
-    scale           1;
 }
 
 

tutorials/lagrangian/reactingParcelFoam/filter/constant/chemistryProperties

@@ -41,7 +41,6 @@ odeCoeffs
 {
     solver          RK;
     eps             0.05;
-    scale           1;
 }
 
 

tutorials/lagrangian/reactingParcelFoam/parcelInBox/constant/chemistryProperties

@@ -41,7 +41,6 @@ odeCoeffs
 {
     solver          RK;
     eps             0.05;
-    scale           1;
 }
 
 

tutorials/lagrangian/reactingParcelFoam/verticalChannel/constant/chemistryProperties

@@ -41,7 +41,6 @@ odeCoeffs
 {
     solver          RK;
     eps             0.05;
-    scale           1;
 }