diff --git a/include/cantera/thermo/ConstCpPoly.h b/include/cantera/thermo/ConstCpPoly.h index 2a81f7af9..7edcaa32b 100644 --- a/include/cantera/thermo/ConstCpPoly.h +++ b/include/cantera/thermo/ConstCpPoly.h @@ -43,7 +43,9 @@ namespace Cantera class ConstCpPoly: public SpeciesThermoInterpType { public: - //! Normal constructor + ConstCpPoly(); + + //! Constructor with all input data /*! * @param tlow Minimum temperature * @param thigh Maximum temperature @@ -58,6 +60,14 @@ public: */ ConstCpPoly(double tlow, double thigh, double pref, const double* coeffs); + /*! + * @param t0 \f$ T_0 \f$ [K] + * @param h0 \f$ h_k^o(T_0, p_{ref}) \f$ [J/kmol] + * @param s0 \f$ s_k^o(T_0, p_{ref}) \f$ [J/kmol/K] + * @param cp0 \f$ c_{p,k}^o(T_0, p_{ref}) \f$ [J/kmol/K] + */ + void setParameters(double t0, double h0, double s0, double cp0); + virtual int reportType() const { return CONSTANT_CP; } diff --git a/include/cantera/thermo/Mu0Poly.h b/include/cantera/thermo/Mu0Poly.h index 26293024c..509ce4c98 100644 --- a/include/cantera/thermo/Mu0Poly.h +++ b/include/cantera/thermo/Mu0Poly.h @@ -73,11 +73,10 @@ class XML_Node; class Mu0Poly: public SpeciesThermoInterpType { public: - //! Normal constructor + Mu0Poly(); + + //! Constructor with all input data /*! - * In the constructor, we calculate and store the piecewise linear - * approximation to the thermodynamic functions. - * * @param tlow Minimum temperature * @param thigh Maximum temperature * @param pref reference pressure (Pa). @@ -98,6 +97,18 @@ public: */ Mu0Poly(double tlow, double thigh, double pref, const double* coeffs); + //! Set parameters for \f$ \mu^o(T) \f$ + /*! + * Calculates and stores the piecewise linear approximation to the + * thermodynamic functions. + * + * @param h0 Enthalpy at the reference temperature of 298.15 K [J/kmol] + * @param T_mu Map with temperature [K] as the keys and the Gibbs free + * energy [J/kmol] as the values. Must contain one point at + * 298.15 K. + */ + void setParameters(double h0, const std::map& T_mu); + virtual int reportType() const { return MU0_INTERP; } @@ -145,25 +156,6 @@ protected: //! Heat capacity at the points vector_fp m_cp0_R_int; - -private: - //! process the coefficients - /*! - * In the constructor, we calculate and store the piecewise linear - * approximation to the thermodynamic functions. - * - * @param coeffs coefficients. These are defined as follows: - * - coeffs[0] = number of points (integer) - * - coeffs[1] = \f$ h^o(298.15 K) \f$ (J/kmol) - * - coeffs[2] = \f$ T_1 \f$ (Kelvin) - * - coeffs[3] = \f$ \mu^o(T_1) \f$ (J/kmol) - * - coeffs[4] = \f$ T_2 \f$ (Kelvin) - * - coeffs[5] = \f$ \mu^o(T_2) \f$ (J/kmol) - * - coeffs[6] = \f$ T_3 \f$ (Kelvin) - * - coeffs[7] = \f$ \mu^o(T_3) \f$ (J/kmol) - * - ........ - */ - void processCoeffs(const doublereal* coeffs); }; //! Install a Mu0 polynomial thermodynamic reference state diff --git a/include/cantera/thermo/Nasa9Poly1.h b/include/cantera/thermo/Nasa9Poly1.h index 8282bbd66..e653d2bf7 100644 --- a/include/cantera/thermo/Nasa9Poly1.h +++ b/include/cantera/thermo/Nasa9Poly1.h @@ -61,7 +61,9 @@ namespace Cantera class Nasa9Poly1 : public SpeciesThermoInterpType { public: - //! Normal constructor + Nasa9Poly1(); + + //! Constructor with all input data /*! * @param tlow Minimum temperature * @param thigh Maximum temperature @@ -71,6 +73,9 @@ public: */ Nasa9Poly1(double tlow, double thigh, double pref, const double* coeffs); + //! Set the array of 9 polynomial coefficients + void setParameters(const vector_fp& coeffs); + virtual int reportType() const; virtual size_t temperaturePolySize() const { return 7; } diff --git a/include/cantera/thermo/Nasa9PolyMultiTempRegion.h b/include/cantera/thermo/Nasa9PolyMultiTempRegion.h index 5e153ebad..74e1783e3 100644 --- a/include/cantera/thermo/Nasa9PolyMultiTempRegion.h +++ b/include/cantera/thermo/Nasa9PolyMultiTempRegion.h @@ -36,7 +36,9 @@ namespace Cantera class Nasa9PolyMultiTempRegion : public SpeciesThermoInterpType { public: - //! Constructor used in templated instantiations + Nasa9PolyMultiTempRegion(); + + //! Constructor with all input data /*! * @param regionPts Vector of pointers to Nasa9Poly1 objects. These objects * all refer to the temperature regions for the same species. The vector @@ -49,6 +51,14 @@ public: */ Nasa9PolyMultiTempRegion(std::vector ®ionPts); + //! Set the array of polynomial coefficients for each temperature region + /*! + * @param regions Map where each key is the minimum temperature for a + * region and each value is the array of 9 polynomial + * coefficients for that region. + */ + void setParameters(const std::map& regions); + virtual ~Nasa9PolyMultiTempRegion(); virtual int reportType() const; diff --git a/include/cantera/thermo/NasaPoly1.h b/include/cantera/thermo/NasaPoly1.h index 252e3132c..39e4b40af 100644 --- a/include/cantera/thermo/NasaPoly1.h +++ b/include/cantera/thermo/NasaPoly1.h @@ -45,7 +45,9 @@ namespace Cantera class NasaPoly1 : public SpeciesThermoInterpType { public: - //! Normal constructor + NasaPoly1() : m_coeff(7), m_coeff5_orig(0.0) {} + + //! Constructor with all input data /*! * @param tlow Minimum temperature * @param thigh Maximum temperature @@ -60,6 +62,16 @@ public: m_coeff5_orig = m_coeff[5]; } + //! Set array of 7 polynomial coefficients + void setParameters(const vector_fp& coeffs) { + if (coeffs.size() != 7) { + throw CanteraError("NasaPoly1::setParameters", "Array must contain " + "7 coefficients, but {} were given.", coeffs.size()); + } + m_coeff = coeffs; + m_coeff5_orig = m_coeff[5]; + } + virtual int reportType() const { return NASA1; } diff --git a/include/cantera/thermo/NasaPoly2.h b/include/cantera/thermo/NasaPoly2.h index 678c1adcc..d615acc30 100644 --- a/include/cantera/thermo/NasaPoly2.h +++ b/include/cantera/thermo/NasaPoly2.h @@ -48,7 +48,9 @@ namespace Cantera class NasaPoly2 : public SpeciesThermoInterpType { public: - //! Full Constructor + NasaPoly2(); + + //! Constructor with all input data /*! * @param tlow output - Minimum temperature * @param thigh output - Maximum temperature @@ -66,6 +68,30 @@ public: mnp_high(coeffs[0], thigh, pref, coeffs + 1) { } + virtual void setMinTemp(double Tmin) { + SpeciesThermoInterpType::setMinTemp(Tmin); + mnp_low.setMinTemp(Tmin); + } + + virtual void setMaxTemp(double Tmax) { + SpeciesThermoInterpType::setMaxTemp(Tmax); + mnp_high.setMaxTemp(Tmax); + } + + virtual void setRefPressure(double Pref) { + SpeciesThermoInterpType::setRefPressure(Pref); + mnp_low.setRefPressure(Pref); + mnp_high.setRefPressure(Pref); + } + + /*! + * @param Tmid Temperature [K] at the boundary between the low and high + * temperature polynomials + * @param low Vector of 7 coefficients for the low temperature polynomial + * @param high Vector of 7 coefficients for the high temperature polynomial + */ + void setParameters(double Tmid, const vector_fp& low, const vector_fp& high); + virtual int reportType() const { return NASA2; } diff --git a/include/cantera/thermo/ShomatePoly.h b/include/cantera/thermo/ShomatePoly.h index d6c88378b..e5d7eee8e 100644 --- a/include/cantera/thermo/ShomatePoly.h +++ b/include/cantera/thermo/ShomatePoly.h @@ -57,7 +57,9 @@ namespace Cantera class ShomatePoly : public SpeciesThermoInterpType { public: - //! Normal constructor + ShomatePoly() : m_coeff(7), m_coeff5_orig(0.0) {} + + //! Constructor with all input data /*! * @param tlow Minimum temperature * @param thigh Maximum temperature @@ -78,6 +80,19 @@ public: m_coeff5_orig = m_coeff[5]; } + //! Set array of 7 polynomial coefficients. Input values are assumed to be + //! on a kJ/mol basis. + void setParameters(const vector_fp& coeffs) { + if (coeffs.size() != 7) { + throw CanteraError("ShomatePoly::setParameters", "Array must " + "contain 7 coefficients, but {} were given.", coeffs.size()); + } + for (size_t i = 0; i < 7; i++) { + m_coeff[i] = coeffs[i] * 1000 / GasConstant; + } + m_coeff5_orig = m_coeff[5]; + } + virtual int reportType() const { return SHOMATE; } @@ -210,7 +225,9 @@ protected: class ShomatePoly2 : public SpeciesThermoInterpType { public: - //! Normal constructor + ShomatePoly2() : m_midT(0.0) {} + + //! Constructor with all input data /*! * @param tlow Minimum temperature * @param thigh Maximum temperature @@ -226,6 +243,36 @@ public: { } + virtual void setMinTemp(double Tmin) { + SpeciesThermoInterpType::setMinTemp(Tmin); + msp_low.setMinTemp(Tmin); + } + + virtual void setMaxTemp(double Tmax) { + SpeciesThermoInterpType::setMaxTemp(Tmax); + msp_high.setMaxTemp(Tmax); + } + + virtual void setRefPressure(double Pref) { + SpeciesThermoInterpType::setRefPressure(Pref); + msp_low.setRefPressure(Pref); + msp_high.setRefPressure(Pref); + } + + /*! + * @param Tmid Temperature [K] at the boundary between the low and high + * temperature polynomials + * @param low Vector of 7 coefficients for the low temperature polynomial + * @param high Vector of 7 coefficients for the high temperature polynomial + */ + void setParameters(double Tmid, const vector_fp& low, const vector_fp& high) { + m_midT = Tmid; + msp_low.setMaxTemp(Tmid); + msp_high.setMinTemp(Tmid); + msp_low.setParameters(low); + msp_high.setParameters(high); + } + virtual int reportType() const { return SHOMATE2; } diff --git a/include/cantera/thermo/SpeciesThermoInterpType.h b/include/cantera/thermo/SpeciesThermoInterpType.h index 329675c73..675155332 100644 --- a/include/cantera/thermo/SpeciesThermoInterpType.h +++ b/include/cantera/thermo/SpeciesThermoInterpType.h @@ -129,17 +129,32 @@ public: return m_lowT; } + //! Set the minimum temperature at which the thermo parameterization is valid + virtual void setMinTemp(double Tmin) { + m_lowT = Tmin; + } + //! Returns the maximum temperature that the thermo parameterization is //! valid virtual doublereal maxTemp() const { return m_highT; } + //! Set the maximum temperature at which the thermo parameterization is valid + virtual void setMaxTemp(double Tmax) { + m_highT = Tmax; + } + //! Returns the reference pressure (Pa) virtual doublereal refPressure() const { return m_Pref; } + //! Set the reference pressure [Pa] + virtual void setRefPressure(double Pref) { + m_Pref = Pref; + } + //! Check for problems with the parameterization, and generate warnings or //! throw and exception if any are found. virtual void validate(const std::string& name) {} diff --git a/src/thermo/ConstCpPoly.cpp b/src/thermo/ConstCpPoly.cpp index b2cb5b4f1..ae318c0a4 100644 --- a/src/thermo/ConstCpPoly.cpp +++ b/src/thermo/ConstCpPoly.cpp @@ -13,16 +13,30 @@ namespace Cantera { +ConstCpPoly::ConstCpPoly() + : m_t0(0.0) + , m_cp0_R(0.0) + , m_h0_R(0.0) + , m_s0_R(0.0) + , m_logt0(0.0) + , m_h0_R_orig(0.0) +{ +} + ConstCpPoly::ConstCpPoly(double tlow, double thigh, double pref, const double* coeffs) : SpeciesThermoInterpType(tlow, thigh, pref) { - m_t0 = coeffs[0]; - m_h0_R = coeffs[1] / GasConstant; - m_s0_R = coeffs[2] / GasConstant; - m_cp0_R = coeffs[3] / GasConstant; + setParameters(coeffs[0], coeffs[1], coeffs[2], coeffs[3]); +} + +void ConstCpPoly::setParameters(double t0, double h0, double s0, double cp0) +{ + m_t0 = t0; m_logt0 = log(m_t0); - m_h0_R_orig = m_h0_R; + m_cp0_R = cp0 / GasConstant; + m_h0_R = h0 / GasConstant; + m_s0_R = s0 / GasConstant; } void ConstCpPoly::updateProperties(const doublereal* tt, diff --git a/src/thermo/Mu0Poly.cpp b/src/thermo/Mu0Poly.cpp index 86b2dce41..3a10115f2 100644 --- a/src/thermo/Mu0Poly.cpp +++ b/src/thermo/Mu0Poly.cpp @@ -17,13 +17,89 @@ using namespace std; namespace Cantera { +Mu0Poly::Mu0Poly() + : m_numIntervals(0) + , m_H298(0.0) +{ +} Mu0Poly::Mu0Poly(double tlow, double thigh, double pref, const double* coeffs) : SpeciesThermoInterpType(tlow, thigh, pref), m_numIntervals(0), m_H298(0.0) { - processCoeffs(coeffs); + std::map T_mu; + size_t nPoints = (size_t) coeffs[0]; + for (size_t i = 0; i < nPoints; i++) { + T_mu[coeffs[2*i+2]] = coeffs[2*i+3]; + } + setParameters(coeffs[1], T_mu); +} + +void Mu0Poly::setParameters(double h0, const std::map& T_mu) +{ + size_t nPoints = T_mu.size(); + if (nPoints < 2) { + throw CanteraError("Mu0Poly::setParameters", "nPoints must be >= 2"); + } + m_numIntervals = nPoints - 1; + m_H298 = h0 / GasConstant; + + // Distribute the data into the internal arrays, and find the index of the + // point at 298.15 K. + size_t iT298 = npos; + for (const auto& row : T_mu) { + double T1 = row.first; + if (T1 == 298.15) { + iT298 = m_t0_int.size(); + } + m_t0_int.push_back(T1); + m_mu0_R_int.push_back(row.second / GasConstant); + } + if (iT298 == npos) { + throw CanteraError("Mu0Poly", "One temperature has to be 298.15"); + } + + // Resize according to the number of points + m_h0_R_int.resize(nPoints); + m_s0_R_int.resize(nPoints); + m_cp0_R_int.resize(nPoints); + + // Starting from the interval with T298, we go up + m_h0_R_int[iT298] = m_H298; + m_s0_R_int[iT298] = - (m_mu0_R_int[iT298] - m_h0_R_int[iT298]) / m_t0_int[iT298]; + for (size_t i = iT298; i < m_numIntervals; i++) { + double T1 = m_t0_int[i]; + double s1 = m_s0_R_int[i]; + double T2 = m_t0_int[i+1]; + double deltaMu = m_mu0_R_int[i+1] - m_mu0_R_int[i]; + double deltaT = T2 - T1; + double cpi = (deltaMu - T1 * s1 + T2 * s1) / (deltaT - T2 * log(T2/T1)); + m_cp0_R_int[i] = cpi; + m_h0_R_int[i+1] = m_h0_R_int[i] + cpi * deltaT; + m_s0_R_int[i+1] = s1 + cpi * log(T2/T1); + m_cp0_R_int[i+1] = cpi; + } + + // Starting from the interval with T298, we go down + if (iT298 != 0) { + m_h0_R_int[iT298] = m_H298; + m_s0_R_int[iT298] = - (m_mu0_R_int[iT298] - m_h0_R_int[iT298]) / m_t0_int[iT298]; + for (size_t i = iT298 - 1; i != npos; i--) { + double T1 = m_t0_int[i]; + double T2 = m_t0_int[i+1]; + double s2 = m_s0_R_int[i+1]; + double deltaMu = m_mu0_R_int[i+1] - m_mu0_R_int[i]; + double deltaT = T2 - T1; + double cpi = (deltaMu - T1 * s2 + T2 * s2) / (deltaT - T1 * log(T2/T1)); + m_cp0_R_int[i] = cpi; + m_h0_R_int[i] = m_h0_R_int[i+1] - cpi * deltaT; + m_s0_R_int[i] = s2 - cpi * log(T2/T1); + if (i == (m_numIntervals-1)) { + m_cp0_R_int[i+1] = cpi; + } + } + } } void Mu0Poly::updateProperties(const doublereal* tt, doublereal* cp_R, @@ -134,81 +210,4 @@ Mu0Poly* newMu0ThermoFromXML(const XML_Node& Mu0Node) fpValue(Mu0Node["Pref"]), &c[0]); } -void Mu0Poly::processCoeffs(const doublereal* coeffs) -{ - size_t nPoints = (size_t) coeffs[0]; - if (nPoints < 2) { - throw CanteraError("Mu0Poly", - "nPoints must be >= 2"); - } - m_numIntervals = nPoints - 1; - m_H298 = coeffs[1] / GasConstant; - size_t iT298 = 0; - - // Resize according to the number of points - m_t0_int.resize(nPoints); - m_h0_R_int.resize(nPoints); - m_s0_R_int.resize(nPoints); - m_cp0_R_int.resize(nPoints); - m_mu0_R_int.resize(nPoints); - - // Calculate the T298 interval and make sure that the temperatures are - // strictly monotonic. Also distribute the data into the internal arrays. - bool ifound = false; - for (size_t i = 0, iindex = 2; i < nPoints; i++) { - double T1 = coeffs[iindex]; - m_t0_int[i] = T1; - m_mu0_R_int[i] = coeffs[iindex+1] / GasConstant; - if (T1 == 298.15) { - iT298 = i; - ifound = true; - } - if (i < nPoints - 1 && coeffs[iindex+2] <= T1) { - throw CanteraError("Mu0Poly", - "Temperatures are not monotonic increasing"); - } - iindex += 2; - } - if (!ifound) { - throw CanteraError("Mu0Poly", - "One temperature has to be 298.15"); - } - - // Starting from the interval with T298, we go up - m_h0_R_int[iT298] = m_H298; - m_s0_R_int[iT298] = - (m_mu0_R_int[iT298] - m_h0_R_int[iT298]) / m_t0_int[iT298]; - for (size_t i = iT298; i < m_numIntervals; i++) { - double T1 = m_t0_int[i]; - double s1 = m_s0_R_int[i]; - double T2 = m_t0_int[i+1]; - double deltaMu = m_mu0_R_int[i+1] - m_mu0_R_int[i]; - double deltaT = T2 - T1; - double cpi = (deltaMu - T1 * s1 + T2 * s1) / (deltaT - T2 * log(T2/T1)); - m_cp0_R_int[i] = cpi; - m_h0_R_int[i+1] = m_h0_R_int[i] + cpi * deltaT; - m_s0_R_int[i+1] = s1 + cpi * log(T2/T1); - m_cp0_R_int[i+1] = cpi; - } - - // Starting from the interval with T298, we go down - if (iT298 != 0) { - m_h0_R_int[iT298] = m_H298; - m_s0_R_int[iT298] = - (m_mu0_R_int[iT298] - m_h0_R_int[iT298]) / m_t0_int[iT298]; - for (size_t i = iT298 - 1; i != npos; i--) { - double T1 = m_t0_int[i]; - double T2 = m_t0_int[i+1]; - double s2 = m_s0_R_int[i+1]; - double deltaMu = m_mu0_R_int[i+1] - m_mu0_R_int[i]; - double deltaT = T2 - T1; - double cpi = (deltaMu - T1 * s2 + T2 * s2) / (deltaT - T1 * log(T2/T1)); - m_cp0_R_int[i] = cpi; - m_h0_R_int[i] = m_h0_R_int[i+1] - cpi * deltaT; - m_s0_R_int[i] = s2 - cpi * log(T2/T1); - if (i == (m_numIntervals-1)) { - m_cp0_R_int[i+1] = cpi; - } - } - } -} - } diff --git a/src/thermo/Nasa9Poly1.cpp b/src/thermo/Nasa9Poly1.cpp index 9c2d309d4..028faf540 100644 --- a/src/thermo/Nasa9Poly1.cpp +++ b/src/thermo/Nasa9Poly1.cpp @@ -17,6 +17,11 @@ namespace Cantera { +Nasa9Poly1::Nasa9Poly1() + : m_coeff(9) +{ +} + Nasa9Poly1::Nasa9Poly1(double tlow, double thigh, double pref, const double* coeffs) : SpeciesThermoInterpType(tlow, thigh, pref), @@ -24,6 +29,15 @@ Nasa9Poly1::Nasa9Poly1(double tlow, double thigh, double pref, { } +void Nasa9Poly1::setParameters(const vector_fp &coeffs) +{ + if (coeffs.size() != 9) { + throw CanteraError("Nasa9Poly1::setParameters", "Array must contain " + "9 coefficients, but {} were given.", coeffs.size()); + } + m_coeff = coeffs; +} + int Nasa9Poly1::reportType() const { return NASA9; diff --git a/src/thermo/Nasa9PolyMultiTempRegion.cpp b/src/thermo/Nasa9PolyMultiTempRegion.cpp index 33844bacb..730984b2b 100644 --- a/src/thermo/Nasa9PolyMultiTempRegion.cpp +++ b/src/thermo/Nasa9PolyMultiTempRegion.cpp @@ -22,6 +22,11 @@ using namespace std; namespace Cantera { +Nasa9PolyMultiTempRegion::Nasa9PolyMultiTempRegion() + : m_currRegion(0) +{ +} + Nasa9PolyMultiTempRegion::Nasa9PolyMultiTempRegion(vector& regionPts) : m_currRegion(0) { @@ -52,6 +57,24 @@ Nasa9PolyMultiTempRegion::Nasa9PolyMultiTempRegion(vector& regionPt } } +void Nasa9PolyMultiTempRegion::setParameters(const std::map& regions) +{ + m_regionPts.clear(); + m_lowerTempBounds.clear(); + for (const auto& region : regions) { + m_lowerTempBounds.push_back(region.first); + Nasa9Poly1* poly = new Nasa9Poly1; + poly->setRefPressure(refPressure()); + poly->setMinTemp(region.first); + poly->setParameters(region.second); + if (!m_regionPts.empty()) { + m_regionPts.back()->setMaxTemp(region.first); + } + m_regionPts.emplace_back(poly); + } + m_regionPts.back()->setMaxTemp(maxTemp()); +} + Nasa9PolyMultiTempRegion::~Nasa9PolyMultiTempRegion() { } diff --git a/src/thermo/NasaPoly2.cpp b/src/thermo/NasaPoly2.cpp index f727430ca..1bea91cd3 100644 --- a/src/thermo/NasaPoly2.cpp +++ b/src/thermo/NasaPoly2.cpp @@ -7,6 +7,20 @@ namespace Cantera { +NasaPoly2::NasaPoly2() + : m_midT(0) +{ +} + +void NasaPoly2::setParameters(double Tmid, const vector_fp& low, + const vector_fp& high) { + m_midT = Tmid; + mnp_low.setMaxTemp(Tmid); + mnp_high.setMinTemp(Tmid); + mnp_low.setParameters(low); + mnp_high.setParameters(high); +} + void NasaPoly2::validate(const std::string& name) { if (thermo_warnings_suppressed()) {