From d43774e721b8af3a5f06356bd45fcccec497189a Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Fri, 20 Sep 2019 13:55:58 -0500 Subject: [PATCH] [SpeciesThermo] enable cython constructors for Nasa9PolyMultiTempRegion Further: * Add unit tests for ShomatePoly2 and Mu0Poly --- .../cantera/thermo/Nasa9PolyMultiTempRegion.h | 12 +++- .../cantera/thermo/SpeciesThermoInterpType.h | 2 - interfaces/cython/cantera/speciesthermo.pyx | 14 ++-- interfaces/cython/cantera/test/test_thermo.py | 65 +++++++++++++++++-- src/thermo/Mu0Poly.cpp | 2 +- src/thermo/Nasa9PolyMultiTempRegion.cpp | 37 +++++++++++ src/thermo/SpeciesThermoFactory.cpp | 2 + 7 files changed, 118 insertions(+), 16 deletions(-) diff --git a/include/cantera/thermo/Nasa9PolyMultiTempRegion.h b/include/cantera/thermo/Nasa9PolyMultiTempRegion.h index 55bc4a40e..f6dd79a76 100644 --- a/include/cantera/thermo/Nasa9PolyMultiTempRegion.h +++ b/include/cantera/thermo/Nasa9PolyMultiTempRegion.h @@ -51,7 +51,17 @@ public: */ Nasa9PolyMultiTempRegion(std::vector ®ionPts); - //! Set the array of polynomial coefficients for each temperature region + //! Constructor with all input data + /*! + * @param tlow Minimum temperature + * @param thigh Maximum temperature + * @param pref reference pressure (Pa). + * @param coeffs Vector of coefficients used to set the parameters for the + * standard state. + */ + Nasa9PolyMultiTempRegion(double tlow, double thigh, double pref, const double* coeffs); + + //! 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 diff --git a/include/cantera/thermo/SpeciesThermoInterpType.h b/include/cantera/thermo/SpeciesThermoInterpType.h index 0c56f3a3e..b49dc486d 100644 --- a/include/cantera/thermo/SpeciesThermoInterpType.h +++ b/include/cantera/thermo/SpeciesThermoInterpType.h @@ -206,12 +206,10 @@ public: doublereal* h_RT, doublereal* s_R) const; - //! This utility function reports back the number of coefficients //! for a given type of species parameterization virtual size_t nCoeffs() const; - //! This utility function reports back the type of parameterization and all //! of the parameters for the species. /*! diff --git a/interfaces/cython/cantera/speciesthermo.pyx b/interfaces/cython/cantera/speciesthermo.pyx index 0fa06df1e..f67577343 100644 --- a/interfaces/cython/cantera/speciesthermo.pyx +++ b/interfaces/cython/cantera/speciesthermo.pyx @@ -32,7 +32,7 @@ cdef class SpeciesThermo: if not init: return - if not self.check_n_coeffs(len(coeffs)): + if not self._check_n_coeffs(len(coeffs)): raise ValueError("Coefficient array has incorrect length") cdef np.ndarray[np.double_t, ndim=1] data = np.ascontiguousarray( coeffs, dtype=np.double) @@ -79,7 +79,7 @@ cdef class SpeciesThermo: T_high, P_ref, &data[0]) return data - def check_n_coeffs(self, n): + def _check_n_coeffs(self, n): raise NotImplementedError('Needs to be overloaded') def cp(self, T): @@ -118,7 +118,7 @@ cdef class ConstantCp(SpeciesThermo): """ derived_type = SPECIES_THERMO_CONSTANT_CP - def check_n_coeffs(self, n): + def _check_n_coeffs(self, n): return n == 4 @@ -141,7 +141,7 @@ cdef class Mu0Poly(SpeciesThermo): """ derived_type = SPECIES_THERMO_MU0_INTERP - def check_n_coeffs(self, n): + def _check_n_coeffs(self, n): return n > 3 and n % 2 == 0 @@ -166,7 +166,7 @@ cdef class NasaPoly2(SpeciesThermo): """ derived_type = SPECIES_THERMO_NASA2 - def check_n_coeffs(self, n): + def _check_n_coeffs(self, n): return n == 15 @@ -191,7 +191,7 @@ cdef class Nasa9PolyMultiTempRegion(SpeciesThermo): """ derived_type = SPECIES_THERMO_NASA9MULTITEMP - def check_n_coeffs(self, n): + def _check_n_coeffs(self, n): return n > 11 and ((n - 1) % 11) == 0 @@ -217,7 +217,7 @@ cdef class ShomatePoly2(SpeciesThermo): """ derived_type = SPECIES_THERMO_SHOMATE2 - def check_n_coeffs(self, n): + def _check_n_coeffs(self, n): return n == 15 diff --git a/interfaces/cython/cantera/test/test_thermo.py b/interfaces/cython/cantera/test/test_thermo.py index 5b0e49ad6..6aa74ab45 100644 --- a/interfaces/cython/cantera/test/test_thermo.py +++ b/interfaces/cython/cantera/test/test_thermo.py @@ -1185,14 +1185,69 @@ class TestSpeciesThermo(utilities.CanteraTest): self.assertEqual(st.reference_pressure, 101325) self.assertArrayNear(self.h2o_coeffs, st.coeffs) self.assertTrue(st.n_coeffs == len(st.coeffs)) - self.assertTrue(st.check_n_coeffs(st.n_coeffs)) + self.assertTrue(st._check_n_coeffs(st.n_coeffs)) - def test_nasa9(self): + def test_nasa9_load(self): gas = ct.Solution('airNASA9.cti') st = gas.species(3).thermo - self.assertTrue(isinstance(st, ct.Nasa9PolyMultiTempRegion)) - self.assertTrue(st.n_coeffs == len(st.coeffs)) - self.assertTrue(st.check_n_coeffs(st.n_coeffs)) + self.assertIsInstance(st, ct.Nasa9PolyMultiTempRegion) + self.assertEqual(st.n_coeffs, len(st.coeffs)) + self.assertTrue(st._check_n_coeffs(st.n_coeffs)) + + def test_nasa9_create(self): + gas = ct.Solution('airNASA9.cti') + st = gas.species(3).thermo + t_min = st.min_temp + t_max = st.max_temp + p_ref = st.reference_pressure + coeffs = st.coeffs + st2 = ct.Nasa9PolyMultiTempRegion(t_min, t_max, p_ref, coeffs) + self.assertIsInstance(st2, ct.Nasa9PolyMultiTempRegion) + self.assertEqual(st.min_temp, t_min) + self.assertEqual(st.max_temp, t_max) + self.assertEqual(st.reference_pressure, p_ref) + for T in range(300, 20000, 1000): + self.assertNear(st.cp(T), st2.cp(T)) + self.assertNear(st.h(T), st2.h(T)) + self.assertNear(st.s(T), st2.s(T)) + + def test_shomate_load(self): + sol = ct.Solution('thermo-models.yaml', 'molten-salt-Margules') + st = sol.species(0).thermo + self.assertIsInstance(st, ct.ShomatePoly2) + self.assertEqual(st.n_coeffs, len(st.coeffs)) + self.assertTrue(st._check_n_coeffs(st.n_coeffs)) + + def test_shomate_create(self): + sol = ct.Solution('thermo-models.yaml', 'molten-salt-Margules') + st = sol.species(0).thermo + t_min = st.min_temp + t_max = st.max_temp + p_ref = st.reference_pressure + coeffs = st.coeffs + st2 = ct.ShomatePoly2(t_min, t_max, p_ref, coeffs) + self.assertIsInstance(st2, ct.ShomatePoly2) + self.assertEqual(st.min_temp, t_min) + self.assertEqual(st.max_temp, t_max) + self.assertEqual(st.reference_pressure, p_ref) + for T in [300, 500, 700, 900]: + self.assertNear(st.cp(T), st2.cp(T)) + self.assertNear(st.h(T), st2.h(T)) + self.assertNear(st.s(T), st2.s(T)) + + def test_mu0poly_create(self): + # use OH- ion data from test/thermo/phaseConstructors.cpp + h298 = -230.015e6 + T1 = 298.15 + mu1 = -91.50963 + T2 = 333.15 + mu2 = -85 + pref = 101325 + coeffs = [2, h298, T1, mu1*ct.gas_constant*T1, T2, mu2*ct.gas_constant*T2] + st2 = ct.Mu0Poly(200, 3500, pref, coeffs) + self.assertIsInstance(st2, ct.Mu0Poly) + self.assertEqual(st2.n_coeffs, len(coeffs)) + self.assertEqual(st2.n_coeffs, len(st2.coeffs)) class TestQuantity(utilities.CanteraTest): diff --git a/src/thermo/Mu0Poly.cpp b/src/thermo/Mu0Poly.cpp index 2620ec5dc..ed4f5b7f7 100644 --- a/src/thermo/Mu0Poly.cpp +++ b/src/thermo/Mu0Poly.cpp @@ -131,7 +131,7 @@ void Mu0Poly::updatePropertiesTemp(const doublereal T, size_t Mu0Poly::nCoeffs() const { - return 2*m_numIntervals + 2; + return 2*m_numIntervals + 4; } void Mu0Poly::reportParameters(size_t& n, int& type, diff --git a/src/thermo/Nasa9PolyMultiTempRegion.cpp b/src/thermo/Nasa9PolyMultiTempRegion.cpp index 7ea0519a6..3c04c17d6 100644 --- a/src/thermo/Nasa9PolyMultiTempRegion.cpp +++ b/src/thermo/Nasa9PolyMultiTempRegion.cpp @@ -57,6 +57,43 @@ Nasa9PolyMultiTempRegion::Nasa9PolyMultiTempRegion(vector& regionPt } } +Nasa9PolyMultiTempRegion::Nasa9PolyMultiTempRegion(double tlow, double thigh, double pref, + const double* coeffs) +{ + size_t regions = (size_t)coeffs[0]; + + for (size_t i=0; isetRefPressure(pref); + poly->setMinTemp(coeffs[11*i+1]); + poly->setMaxTemp(coeffs[11*i+2]); + for (size_t j=11*i+3; j<=11*(i+1); j++) { + cs.push_back(coeffs[j]); + } + poly->setParameters(cs); + m_regionPts.emplace_back(poly); + } + + m_lowerTempBounds.resize(regions); + m_lowT = tlow; + m_highT = thigh; + m_Pref = pref; + for (size_t i = 0; i < m_regionPts.size(); i++) { + m_lowerTempBounds[i] = m_regionPts[i]->minTemp(); + if (i > 0) { + if (m_lowerTempBounds[i-1] >= m_lowerTempBounds[i]) { + throw CanteraError("Nasa9PolyMultiTempRegion::Nasa9PolyMultiTempRegion", + "minTemp bounds inconsistency"); + } + if (fabs(m_regionPts[i-1]->maxTemp() - m_lowerTempBounds[i]) > 0.0001) { + throw CanteraError("Nasa9PolyMultiTempRegion::Nasa9PolyMultiTempRegion", + "Temp bounds inconsistency"); + } + } + } +} + void Nasa9PolyMultiTempRegion::setParameters(const std::map& regions) { m_regionPts.clear(); diff --git a/src/thermo/SpeciesThermoFactory.cpp b/src/thermo/SpeciesThermoFactory.cpp index a9ceaad08..450eae854 100644 --- a/src/thermo/SpeciesThermoFactory.cpp +++ b/src/thermo/SpeciesThermoFactory.cpp @@ -44,6 +44,8 @@ SpeciesThermoInterpType* newSpeciesThermoInterpType(int type, double tlow, return new ShomatePoly2(tlow, thigh, pref, coeffs); case NASA2: return new NasaPoly2(tlow, thigh, pref, coeffs); + case NASA9MULTITEMP: + return new Nasa9PolyMultiTempRegion(tlow, thigh, pref, coeffs); default: throw CanteraError("newSpeciesThermoInterpType", "Unknown species thermo type: {}.", type);