From d8f9f97f74944643330ac8cbc1223f9314f9391a Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 24 Apr 2013 21:47:48 +0000 Subject: [PATCH] Refactored class NasaThermo to remove implementations from header --- src/thermo/NasaThermo.cpp | 334 ++++++++++++++++++++++++++++ src/thermo/NasaThermo.h | 292 +++--------------------- src/thermo/SpeciesThermoFactory.cpp | 70 ------ 3 files changed, 365 insertions(+), 331 deletions(-) create mode 100644 src/thermo/NasaThermo.cpp diff --git a/src/thermo/NasaThermo.cpp b/src/thermo/NasaThermo.cpp new file mode 100644 index 000000000..981971754 --- /dev/null +++ b/src/thermo/NasaThermo.cpp @@ -0,0 +1,334 @@ +/*! + * @file NasaThermo.cpp Implementation of class Cantera::NasaThermo + */ +#include "NasaThermo.h" +#include "cantera/base/utilities.h" + +namespace Cantera +{ + +NasaThermo::NasaThermo() : + ID(NASA), + m_tlow_max(0.0), + m_thigh_min(1.e30), + m_p0(-1.0), + m_ngroups(0) { + m_t.resize(6); + } + +NasaThermo::NasaThermo(const NasaThermo& right) : + ID(NASA), + m_tlow_max(0.0), + m_thigh_min(1.e30), + m_p0(-1.0), + m_ngroups(0) { + *this = operator=(right); +} + +NasaThermo& NasaThermo::operator=(const NasaThermo& right) +{ + /* + * Check for self assignment. + */ + if (this == &right) { + return *this; + } + + m_high = right.m_high; + m_low = right.m_low; + m_index = right.m_index; + m_tmid = right.m_tmid; + m_tlow_max = right.m_tlow_max; + m_thigh_min = right.m_thigh_min; + m_tlow = right.m_tlow; + m_thigh = right.m_thigh; + m_p0 = right.m_p0; + m_ngroups = right.m_ngroups; + m_t = right.m_t; + m_group_map = right.m_group_map; + m_posInGroup_map = right.m_posInGroup_map; + m_name = right.m_name; + + return *this; +} + +void NasaThermo::install(const std::string& name, size_t index, int type, + const doublereal* c, + doublereal min_temp, doublereal max_temp, + doublereal ref_pressure) +{ + m_name[index] = name; + int imid = int(c[0]); // midpoint temp converted to integer + int igrp = m_index[imid]; // has this value been seen before? + if (igrp == 0) { // if not, prepare new group + std::vector v; + m_high.push_back(v); + m_low.push_back(v); + m_tmid.push_back(c[0]); + m_index[imid] = igrp = static_cast(m_high.size()); + m_ngroups++; + } + + m_group_map[index] = igrp; + m_posInGroup_map[index] = (int) m_low[igrp-1].size(); + + doublereal tlow = min_temp; + doublereal tmid = c[0]; + doublereal thigh = max_temp; + + vector_fp chigh(c+8, c+15); + vector_fp clow(c+1, c+8); + + ensureContinuity(name, tmid, &clow[0], &chigh[0]); + + m_high[igrp-1].push_back(NasaPoly1(index, tmid, thigh, + ref_pressure, &chigh[0])); + m_low[igrp-1].push_back(NasaPoly1(index, tlow, tmid, + ref_pressure, &clow[0])); + + if (tlow > m_tlow_max) { + m_tlow_max = tlow; + } + if (thigh < m_thigh_min) { + m_thigh_min = thigh; + } + if (m_tlow.size() < index + 1) { + m_tlow.resize(index + 1, tlow); + m_thigh.resize(index + 1, thigh); + } + m_tlow[index] = tlow; + m_thigh[index] = thigh; + if (m_p0 < 0.0) { + m_p0 = ref_pressure; + } else if (fabs(m_p0 - ref_pressure) > 0.1) { + std::string logmsg = " ERROR NasaThermo: New Species, " + name + ", has a different reference pressure, " + + fp2str(ref_pressure) + ", than existing reference pressure, " + fp2str(m_p0) + "\n"; + writelog(logmsg); + logmsg = " This is now a fatal error\n"; + writelog(logmsg); + throw CanteraError("install()", "species have different reference pressures"); + } + m_p0 = ref_pressure; +} + +void NasaThermo::update_one(size_t k, doublereal t, doublereal* cp_R, + doublereal* h_RT, doublereal* s_R) const +{ + m_t[0] = t; + m_t[1] = t*t; + m_t[2] = m_t[1]*t; + m_t[3] = m_t[2]*t; + m_t[4] = 1.0/t; + m_t[5] = log(t); + + size_t grp = m_group_map[k]; + size_t pos = m_posInGroup_map[k]; + const std::vector &mlg = m_low[grp-1]; + const NasaPoly1* nlow = &(mlg[pos]); + + doublereal tmid = nlow->maxTemp(); + if (t < tmid) { + nlow->updateProperties(&m_t[0], cp_R, h_RT, s_R); + } else { + const std::vector &mhg = m_high[grp-1]; + const NasaPoly1* nhigh = &(mhg[pos]); + nhigh->updateProperties(&m_t[0], cp_R, h_RT, s_R); + } +} + +void NasaThermo::update(doublereal t, doublereal* cp_R, + doublereal* h_RT, doublereal* s_R) const +{ + int i; + + // load functions of temperature into m_t vector + m_t[0] = t; + m_t[1] = t*t; + m_t[2] = m_t[1]*t; + m_t[3] = m_t[2]*t; + m_t[4] = 1.0/t; + m_t[5] = log(t); + + // iterate over the groups + std::vector::const_iterator _begin, _end; + for (i = 0; i != m_ngroups; i++) { + if (t > m_tmid[i]) { + _begin = m_high[i].begin(); + _end = m_high[i].end(); + } else { + _begin = m_low[i].begin(); + _end = m_low[i].end(); + } + for (; _begin != _end; ++_begin) { + _begin->updateProperties(&m_t[0], cp_R, h_RT, s_R); + } + } +} + +void NasaThermo::reportParams(size_t index, int& type, + doublereal* const c, + doublereal& minTemp, + doublereal& maxTemp, + doublereal& refPressure) const +{ + type = reportType(index); + if (type == NASA) { + size_t grp = m_group_map[index]; + size_t pos = m_posInGroup_map[index]; + const std::vector &mlg = m_low[grp-1]; + const std::vector &mhg = m_high[grp-1]; + const NasaPoly1* lowPoly = &(mlg[pos]); + const NasaPoly1* highPoly = &(mhg[pos]); + int itype = NASA; + doublereal tmid = lowPoly->maxTemp(); + c[0] = tmid; + size_t n; + double ttemp; + lowPoly->reportParameters(n, itype, minTemp, ttemp, refPressure, + c + 1); + if (n != index) { + throw CanteraError(" ", "confused"); + } + if (itype != NASA1) { + throw CanteraError(" ", "confused"); + } + highPoly->reportParameters(n, itype, ttemp, maxTemp, refPressure, + c + 8); + if (n != index) { + throw CanteraError(" ", "confused"); + } + if (itype != NASA1) { + throw CanteraError(" ", "confused"); + } + } else { + throw CanteraError(" ", "confused"); + } +} + +#ifdef H298MODIFY_CAPABILITY +doublereal NasaThermo::reportOneHf298(const int k) const +{ + int grp = m_group_map[k]; + int pos = m_posInGroup_map[k]; + const std::vector &mlg = m_low[grp-1]; + const NasaPoly1* nlow = &(mlg[pos]); + doublereal tmid = nlow->maxTemp(); + double h; + if (298.15 <= tmid) { + h = nlow->reportHf298(0); + } else { + const std::vector &mhg = m_high[grp-1]; + const NasaPoly1* nhigh = &(mhg[pos]); + h = nhigh->reportHf298(0); + } + return h; +} + +void NasaThermo::modifyOneHf298(const int k, const doublereal Hf298New) +{ + int grp = m_group_map[k]; + int pos = m_posInGroup_map[k]; + std::vector &mlg = m_low[grp-1]; + NasaPoly1* nlow = &(mlg[pos]); + std::vector &mhg = m_high[grp-1]; + NasaPoly1* nhigh = &(mhg[pos]); + doublereal tmid = nlow->maxTemp(); + + double hnow = reportOneHf298(k); + double delH = Hf298New - hnow; + if (298.15 <= tmid) { + nlow->modifyOneHf298(k, Hf298New); + double h = nhigh->reportHf298(0); + double hnew = h + delH; + nhigh->modifyOneHf298(k, hnew); + } else { + nhigh->modifyOneHf298(k, Hf298New); + double h = nlow->reportHf298(0); + double hnew = h + delH; + nlow->modifyOneHf298(k, hnew); + } +} +#endif + +doublereal NasaThermo::enthalpy_RT(double t, const doublereal* c) { + return c[2] + 0.5*c[3]*t + OneThird*c[4]*t*t + + 0.25*c[5]*t*t*t + 0.2*c[6]*t*t*t*t + + c[0]/t; +} + +doublereal NasaThermo::entropy_R(double t, const doublereal* c) { + return c[2]*log(t) + c[3]*t + 0.5*c[4]*t*t + + OneThird*c[5]*t*t*t + 0.25*c[6]*t*t*t*t + + c[1]; +} + +void NasaThermo::ensureContinuity(const std::string& name, double tmid, + doublereal* clow, doublereal* chigh) +{ + // heat capacity + doublereal cplow = poly4(tmid, clow + 2); + doublereal cphigh = poly4(tmid, chigh + 2); + doublereal delta = cplow - cphigh; + if (fabs(delta/(fabs(cplow)+1.0E-4)) > 0.001) { + writelog("\n\n**** WARNING ****\nFor species "+name+ + ", discontinuity in cp/R detected at Tmid = " + +fp2str(tmid)+"\n"); + writelog("\tValue computed using low-temperature polynomial: " + +fp2str(cplow)+".\n"); + writelog("\tValue computed using high-temperature polynomial: " + +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); + delta = hrtlow - hrthigh; + if (fabs(delta/(fabs(hrtlow)+cplow*tmid)) > 0.001) { + writelog("\n\n**** WARNING ****\nFor species "+name+ + ", discontinuity in h/RT detected at Tmid = " + +fp2str(tmid)+"\n"); + writelog("\tValue computed using low-temperature polynomial: " + +fp2str(hrtlow)+".\n"); + writelog("\tValue computed using high-temperature polynomial: " + +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); + delta = srlow - srhigh; + if (fabs(delta/(fabs(srlow)+cplow)) > 0.001) { + writelog("\n\n**** WARNING ****\nFor species "+name+ + ", discontinuity in s/R detected at Tmid = " + +fp2str(tmid)+"\n"); + writelog("\tValue computed using low-temperature polynomial: " + +fp2str(srlow)+".\n"); + 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"); +} + +} diff --git a/src/thermo/NasaThermo.h b/src/thermo/NasaThermo.h index 37250f15d..63c2c473c 100644 --- a/src/thermo/NasaThermo.h +++ b/src/thermo/NasaThermo.h @@ -48,66 +48,12 @@ namespace Cantera class NasaThermo : public SpeciesThermo { public: - //! Initialized to the type of parameterization - /*! - * Note, this value is used in some template functions - */ - const int ID; + NasaThermo(); - //! constructor - NasaThermo() : - ID(NASA), - m_tlow_max(0.0), - m_thigh_min(1.e30), - m_p0(-1.0), - m_ngroups(0) { - m_t.resize(6); - } + NasaThermo(const NasaThermo& right); - //! Copy constructor - /*! - * @param right NasaThermo object to be copied. - */ - NasaThermo(const NasaThermo& right) : - ID(NASA), - m_tlow_max(0.0), - m_thigh_min(1.e30), - m_p0(-1.0), - m_ngroups(0) { - *this = operator=(right); - } + NasaThermo& operator=(const NasaThermo& right); - //! Assignment operator - /*! - * @param right NasaThermo object to be copied. - */ - NasaThermo& operator=(const NasaThermo& right) { - /* - * Check for self assignment. - */ - if (this == &right) { - return *this; - } - - m_high = right.m_high; - m_low = right.m_low; - m_index = right.m_index; - m_tmid = right.m_tmid; - m_tlow_max = right.m_tlow_max; - m_thigh_min = right.m_thigh_min; - m_tlow = right.m_tlow; - m_thigh = right.m_thigh; - m_p0 = right.m_p0; - m_ngroups = right.m_ngroups; - m_t = right.m_t; - m_group_map = right.m_group_map; - m_posInGroup_map = right.m_posInGroup_map; - m_name = right.m_name; - - return *this; - } - - //! destructor virtual ~NasaThermo() {} virtual SpeciesThermo* duplMyselfAsSpeciesThermo() const { @@ -138,61 +84,7 @@ public: virtual void install(const std::string& name, size_t index, int type, const doublereal* c, doublereal min_temp, doublereal max_temp, - doublereal ref_pressure) { - - m_name[index] = name; - int imid = int(c[0]); // midpoint temp converted to integer - int igrp = m_index[imid]; // has this value been seen before? - if (igrp == 0) { // if not, prepare new group - std::vector v; - m_high.push_back(v); - m_low.push_back(v); - m_tmid.push_back(c[0]); - m_index[imid] = igrp = static_cast(m_high.size()); - m_ngroups++; - } - - m_group_map[index] = igrp; - m_posInGroup_map[index] = (int) m_low[igrp-1].size(); - - doublereal tlow = min_temp; - doublereal tmid = c[0]; - doublereal thigh = max_temp; - - vector_fp chigh(c+8, c+15); - vector_fp clow(c+1, c+8); - - ensureContinuity(name, tmid, &clow[0], &chigh[0]); - - m_high[igrp-1].push_back(NasaPoly1(index, tmid, thigh, - ref_pressure, &chigh[0])); - m_low[igrp-1].push_back(NasaPoly1(index, tlow, tmid, - ref_pressure, &clow[0])); - - if (tlow > m_tlow_max) { - m_tlow_max = tlow; - } - if (thigh < m_thigh_min) { - m_thigh_min = thigh; - } - if (m_tlow.size() < index + 1) { - m_tlow.resize(index + 1, tlow); - m_thigh.resize(index + 1, thigh); - } - m_tlow[index] = tlow; - m_thigh[index] = thigh; - if (m_p0 < 0.0) { - m_p0 = ref_pressure; - } else if (fabs(m_p0 - ref_pressure) > 0.1) { - std::string logmsg = " ERROR NasaThermo: New Species, " + name + ", has a different reference pressure, " - + fp2str(ref_pressure) + ", than existing reference pressure, " + fp2str(m_p0) + "\n"; - writelog(logmsg); - logmsg = " This is now a fatal error\n"; - writelog(logmsg); - throw CanteraError("install()", "species have different reference pressures"); - } - m_p0 = ref_pressure; - } + doublereal ref_pressure); virtual void install_STIT(SpeciesThermoInterpType* stit_ptr) { throw CanteraError("install_STIT", "not implemented"); @@ -207,56 +99,10 @@ public: * @param s_R Vector of Dimensionless entropies. (length m_kk). */ virtual void update_one(size_t k, doublereal t, doublereal* cp_R, - doublereal* h_RT, doublereal* s_R) const { - m_t[0] = t; - m_t[1] = t*t; - m_t[2] = m_t[1]*t; - m_t[3] = m_t[2]*t; - m_t[4] = 1.0/t; - m_t[5] = log(t); - - size_t grp = m_group_map[k]; - size_t pos = m_posInGroup_map[k]; - const std::vector &mlg = m_low[grp-1]; - const NasaPoly1* nlow = &(mlg[pos]); - - doublereal tmid = nlow->maxTemp(); - if (t < tmid) { - nlow->updateProperties(&m_t[0], cp_R, h_RT, s_R); - } else { - const std::vector &mhg = m_high[grp-1]; - const NasaPoly1* nhigh = &(mhg[pos]); - nhigh->updateProperties(&m_t[0], cp_R, h_RT, s_R); - } - } + doublereal* h_RT, doublereal* s_R) const; virtual void update(doublereal t, doublereal* cp_R, - doublereal* h_RT, doublereal* s_R) const { - int i; - - // load functions of temperature into m_t vector - m_t[0] = t; - m_t[1] = t*t; - m_t[2] = m_t[1]*t; - m_t[3] = m_t[2]*t; - m_t[4] = 1.0/t; - m_t[5] = log(t); - - // iterate over the groups - std::vector::const_iterator _begin, _end; - for (i = 0; i != m_ngroups; i++) { - if (t > m_tmid[i]) { - _begin = m_high[i].begin(); - _end = m_high[i].end(); - } else { - _begin = m_low[i].begin(); - _end = m_low[i].end(); - } - for (; _begin != _end; ++_begin) { - _begin->updateProperties(&m_t[0], cp_R, h_RT, s_R); - } - } - } + doublereal* h_RT, doublereal* s_R) const; virtual doublereal minTemp(size_t k=npos) const { if (k == npos) { @@ -300,86 +146,19 @@ public: doublereal* const c, doublereal& minTemp, doublereal& maxTemp, - doublereal& refPressure) const { - type = reportType(index); - if (type == NASA) { - size_t grp = m_group_map[index]; - size_t pos = m_posInGroup_map[index]; - const std::vector &mlg = m_low[grp-1]; - const std::vector &mhg = m_high[grp-1]; - const NasaPoly1* lowPoly = &(mlg[pos]); - const NasaPoly1* highPoly = &(mhg[pos]); - int itype = NASA; - doublereal tmid = lowPoly->maxTemp(); - c[0] = tmid; - size_t n; - double ttemp; - lowPoly->reportParameters(n, itype, minTemp, ttemp, refPressure, - c + 1); - if (n != index) { - throw CanteraError(" ", "confused"); - } - if (itype != NASA1) { - throw CanteraError(" ", "confused"); - } - highPoly->reportParameters(n, itype, ttemp, maxTemp, refPressure, - c + 8); - if (n != index) { - throw CanteraError(" ", "confused"); - } - if (itype != NASA1) { - throw CanteraError(" ", "confused"); - } - } else { - throw CanteraError(" ", "confused"); - } - } + doublereal& refPressure) const; #ifdef H298MODIFY_CAPABILITY - virtual doublereal reportOneHf298(const int k) const { - - int grp = m_group_map[k]; - int pos = m_posInGroup_map[k]; - const std::vector &mlg = m_low[grp-1]; - const NasaPoly1* nlow = &(mlg[pos]); - doublereal tmid = nlow->maxTemp(); - double h; - if (298.15 <= tmid) { - h = nlow->reportHf298(0); - } else { - const std::vector &mhg = m_high[grp-1]; - const NasaPoly1* nhigh = &(mhg[pos]); - h = nhigh->reportHf298(0); - } - return h; - } - - virtual void modifyOneHf298(const int k, const doublereal Hf298New) { - int grp = m_group_map[k]; - int pos = m_posInGroup_map[k]; - std::vector &mlg = m_low[grp-1]; - NasaPoly1* nlow = &(mlg[pos]); - std::vector &mhg = m_high[grp-1]; - NasaPoly1* nhigh = &(mhg[pos]); - doublereal tmid = nlow->maxTemp(); - - double hnow = reportOneHf298(k); - double delH = Hf298New - hnow; - if (298.15 <= tmid) { - nlow->modifyOneHf298(k, Hf298New); - double h = nhigh->reportHf298(0); - double hnew = h + delH; - nhigh->modifyOneHf298(k, hnew); - } else { - nhigh->modifyOneHf298(k, Hf298New); - double h = nlow->reportHf298(0); - double hnew = h + delH; - nlow->modifyOneHf298(k, hnew); - } - - } + virtual doublereal reportOneHf298(const int k) const; + virtual void modifyOneHf298(const int k, const doublereal Hf298New); #endif + //! Initialized to the type of parameterization + /*! + * Note, this value is used in some template functions + */ + const int ID; + protected: //! Vector of vector of NasaPoly1's for the high temp region. /*! @@ -459,13 +238,26 @@ protected: //! Species name as a function of the species index mutable std::map m_name; -private: +protected: + //! for internal use by ensureContinuity + /*! + * @param t temperature + * @param c coefficient array + */ + doublereal enthalpy_RT(double t, const doublereal* c); + + //! for internal use by ensureContinuity + /*! + * @param t temperature + * @param c coefficient array + */ + doublereal entropy_R(double t, const doublereal* c); + //! 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 @@ -473,29 +265,7 @@ private: * @param chigh coefficients for higher temperature region */ void ensureContinuity(const std::string& name, double tmid, - doublereal* clow, doublereal* chigh); - - //! for internal use by ensureContinuity - /*! - * @param t temperature - * @param c coefficient array - */ - doublereal enthalpy_RT(double t, const doublereal* c) { - return c[2] + 0.5*c[3]*t + OneThird*c[4]*t*t - + 0.25*c[5]*t*t*t + 0.2*c[6]*t*t*t*t - + c[0]/t; - } - - //! for internal use by ensureContinuity - /*! - * @param t temperature - * @param c coefficient array - */ - doublereal entropy_R(double t, const doublereal* c) { - return c[2]*log(t) + c[3]*t + 0.5*c[4]*t*t - + OneThird*c[5]*t*t*t + 0.25*c[6]*t*t*t*t - + c[1]; - } + doublereal* clow, doublereal* chigh); }; } diff --git a/src/thermo/SpeciesThermoFactory.cpp b/src/thermo/SpeciesThermoFactory.cpp index 4cb9ee6b0..e34ba504a 100644 --- a/src/thermo/SpeciesThermoFactory.cpp +++ b/src/thermo/SpeciesThermoFactory.cpp @@ -198,76 +198,6 @@ SpeciesThermo* SpeciesThermoFactory::newSpeciesThermoManager(std::string& stype) return (SpeciesThermo*) 0; } -void NasaThermo::ensureContinuity(const std::string& name, double tmid, - doublereal* clow, doublereal* chigh) -{ - // heat capacity - doublereal cplow = poly4(tmid, clow + 2); - doublereal cphigh = poly4(tmid, chigh + 2); - doublereal delta = cplow - cphigh; - if (fabs(delta/(fabs(cplow)+1.0E-4)) > 0.001) { - writelog("\n\n**** WARNING ****\nFor species "+name+ - ", discontinuity in cp/R detected at Tmid = " - +fp2str(tmid)+"\n"); - writelog("\tValue computed using low-temperature polynomial: " - +fp2str(cplow)+".\n"); - writelog("\tValue computed using high-temperature polynomial: " - +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); - delta = hrtlow - hrthigh; - if (fabs(delta/(fabs(hrtlow)+cplow*tmid)) > 0.001) { - writelog("\n\n**** WARNING ****\nFor species "+name+ - ", discontinuity in h/RT detected at Tmid = " - +fp2str(tmid)+"\n"); - writelog("\tValue computed using low-temperature polynomial: " - +fp2str(hrtlow)+".\n"); - writelog("\tValue computed using high-temperature polynomial: " - +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); - delta = srlow - srhigh; - if (fabs(delta/(fabs(srlow)+cplow)) > 0.001) { - writelog("\n\n**** WARNING ****\nFor species "+name+ - ", discontinuity in s/R detected at Tmid = " - +fp2str(tmid)+"\n"); - writelog("\tValue computed using low-temperature polynomial: " - +fp2str(srlow)+".\n"); - 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"); -} - - //! Install a NASA polynomial thermodynamic property parameterization for species k into a SpeciesThermo instance. /*! * This is called by method installThermoForSpecies if a NASA block is found in the XML input.