diff --git a/src/thermo/NasaThermo.h b/src/thermo/NasaThermo.h index 618421d4b..b07a747fc 100644 --- a/src/thermo/NasaThermo.h +++ b/src/thermo/NasaThermo.h @@ -162,7 +162,7 @@ public: vector_fp chigh(c+8, c+15); vector_fp clow(c+1, c+8); - checkContinuity(name, tmid, &clow[0], &chigh[0]); + ensureContinuity(name, tmid, &clow[0], &chigh[0]); m_high[igrp-1].push_back(NasaPoly1(index, tmid, thigh, refPressure, &chigh[0])); @@ -460,17 +460,22 @@ protected: mutable std::map m_name; private: - //! see SpeciesThermoFactory.cpp for the definition + //! Adjust polynomials to be continuous at the midpoint temperature. /*! + * Check to see if the provided coefficients are nearly continuous. Adjust + * the values to get more precise contintinuity to avoid convergence + * issues with algorithms that expect these quantities to be continuous. + * See SpeciesThermoFactory.cpp for the definition. + * * @param name string name of species * @param tmid Mid temperature, between the two temperature regions * @param clow coefficients for lower temperature region * @param chigh coefficients for higher temperature region */ - void checkContinuity(const std::string& name, double tmid, - const doublereal* clow, doublereal* chigh); + void ensureContinuity(const std::string& name, double tmid, + doublereal* clow, doublereal* chigh); - //! for internal use by checkContinuity + //! for internal use by ensureContinuity /*! * @param t temperature * @param c coefficient array @@ -481,7 +486,7 @@ private: + c[0]/t; } - //! for internal use by checkContinuity + //! for internal use by ensureContinuity /*! * @param t temperature * @param c coefficient array diff --git a/src/thermo/SpeciesThermoFactory.cpp b/src/thermo/SpeciesThermoFactory.cpp index 5e6bf3d70..4cb9ee6b0 100644 --- a/src/thermo/SpeciesThermoFactory.cpp +++ b/src/thermo/SpeciesThermoFactory.cpp @@ -198,8 +198,8 @@ SpeciesThermo* SpeciesThermoFactory::newSpeciesThermoManager(std::string& stype) return (SpeciesThermo*) 0; } -void NasaThermo::checkContinuity(const std::string& name, double tmid, - const doublereal* clow, doublereal* chigh) +void NasaThermo::ensureContinuity(const std::string& name, double tmid, + doublereal* clow, doublereal* chigh) { // heat capacity doublereal cplow = poly4(tmid, clow + 2); @@ -215,6 +215,13 @@ void NasaThermo::checkContinuity(const std::string& name, double tmid, +fp2str(cphigh)+".\n"); } + // Adjust coefficients to eliminate any discontinuity + chigh[2] += 0.5 * delta; + clow[2] -= 0.5 * delta; + + AssertThrowMsg(std::abs(poly4(tmid, clow+2) - poly4(tmid, chigh+2)) < 1e-12, + "NasaThermo::ensureContinuity", "Cp/R does not match"); + // enthalpy doublereal hrtlow = enthalpy_RT(tmid, clow); doublereal hrthigh = enthalpy_RT(tmid, chigh); @@ -229,6 +236,14 @@ void NasaThermo::checkContinuity(const std::string& name, double tmid, +fp2str(hrthigh)+".\n"); } + // Adjust coefficients to eliminate any discontinuity + chigh[0] += 0.5 * delta * tmid; + clow[0] -= 0.5 * delta * tmid; + + AssertThrowMsg(std::abs(enthalpy_RT(tmid, clow) - + enthalpy_RT(tmid, chigh)) < 1e-12, + "NasaThermo::ensureContinuity", "H/RT does not match"); + // entropy doublereal srlow = entropy_R(tmid, clow); doublereal srhigh = entropy_R(tmid, chigh); @@ -242,6 +257,14 @@ void NasaThermo::checkContinuity(const std::string& name, double tmid, writelog("\tValue computed using high-temperature polynomial: " +fp2str(srhigh)+".\n"); } + + // Adjust coefficients to eliminate any discontinuity + chigh[1] += 0.5 * delta; + clow[1] -= 0.5 * delta; + + AssertThrowMsg(std::abs(entropy_R(tmid, clow) - + entropy_R(tmid, chigh)) < 1e-12, + "NasaThermo::ensureContinuity", "S/R does not match"); }