[Python] Add access to SpeciesThermo(InterpType) objects

doc/sphinx/cython/thermo.rst

@@ -9,3 +9,14 @@ These classes are used to describe the thermodynamic state of a system.
 .. autoclass:: InterfacePhase(infile='', phaseid='')
 .. autoclass:: PureFluid(infile='', phaseid='')
 .. autoclass:: Mixture
+
+Species Thermodynamic Properties
+================================
+
+These classes are used to describe the reference-state thermodynamic properties
+of a pure species.
+
+.. autoclass:: SpeciesThermo
+.. autoclass:: ConstantCp(T_low, T_high, P_ref, coeffs)
+.. autoclass:: NasaPoly2(T_low, T_high, P_ref, coeffs)
+.. autoclass:: ShomatePoly2(T_low, T_high, P_ref, coeffs)

interfaces/cython/cantera/_cantera.pxd

@@ -43,17 +43,27 @@ cdef extern from "cantera/base/smart_ptr.h":
         T* get()
         void reset(T*)
 
+cdef extern from "cantera/thermo/SpeciesThermoInterpType.h":
+    cdef cppclass CxxSpeciesThermo "Cantera::SpeciesThermoInterpType":
+        CxxSpeciesThermo()
+        int reportType()
+        void updatePropertiesTemp(double, double*, double*, double*) except +
+
+cdef extern from "cantera/thermo/SpeciesThermoFactory.h":
+    cdef CxxSpeciesThermo* CxxNewSpeciesThermo "Cantera::newSpeciesThermoInterpType"\
+        (int, double, double, double, double*) except +
+
 cdef extern from "cantera/thermo/Species.h":
     cdef cppclass CxxSpecies "Cantera::Species":
         CxxSpecies()
         CxxSpecies(string, stdmap[string,double])
+        shared_ptr[CxxSpeciesThermo] thermo
 
         string name
         stdmap[string,double] composition
         double charge
         double size
 
-
 cdef extern from "cantera/thermo/ThermoPhase.h" namespace "Cantera":
     cdef cppclass CxxThermoPhase "Cantera::ThermoPhase":
         CxxThermoPhase()
@@ -611,6 +621,11 @@ cdef class Species:
 
     cdef _assign(self, shared_ptr[CxxSpecies] other)
 
+cdef class SpeciesThermo:
+    cdef shared_ptr[CxxSpeciesThermo] _spthermo
+    cdef CxxSpeciesThermo* spthermo
+    cdef _assign(self, shared_ptr[CxxSpeciesThermo] other)
+
 cdef class _SolutionBase:
     cdef CxxThermoPhase* thermo
     cdef CxxKinetics* kinetics
@@ -771,3 +786,4 @@ cdef np.ndarray get_reaction_array(Kinetics kin, kineticsMethod1d method)
 cdef np.ndarray get_transport_1d(Transport tran, transportMethod1d method)
 cdef np.ndarray get_transport_2d(Transport tran, transportMethod2d method)
 cdef CxxIdealGasPhase* getIdealGasPhase(ThermoPhase phase) except *
+cdef wrapSpeciesThermo(shared_ptr[CxxSpeciesThermo] spthermo)

interfaces/cython/cantera/_cantera.pyx

@@ -14,6 +14,7 @@ include "constants.pyx"
 include "func1.pyx"
 
 include "base.pyx"
+include "speciesthermo.pyx"
 include "thermo.pyx"
 include "kinetics.pyx"
 include "transport.pyx"

interfaces/cython/cantera/speciesthermo.pyx

@@ -0,0 +1,96 @@
+cdef extern from "cantera/thermo/speciesThermoTypes.h" namespace "Cantera":
+    cdef int SPECIES_THERMO_CONSTANT_CP "CONSTANT_CP"
+    cdef int SPECIES_THERMO_NASA2 "NASA2"
+    cdef int SPECIES_THERMO_SHOMATE2 "SHOMATE2"
+
+
+cdef class SpeciesThermo:
+    """
+    Base class for representing the reference-state thermodynamic properties of
+    a pure species. These properties are a function of temperature. Derived
+    classes implement a parameterization of this temperature dependence. This is
+    a wrapper for the C++ class :ct:`SpeciesThermoInterpType`.
+    """
+    def __cinit__(self, T_low=None, T_high=None, P_ref=None, coeffs=None, *args,
+                  init=True, **kwargs):
+        if not init:
+            return
+
+        if len(coeffs) != self.n_coeffs:
+            raise ValueError("Coefficient array has incorrect length")
+        cdef np.ndarray[np.double_t, ndim=1] data = np.ascontiguousarray(
+            coeffs, dtype=np.double)
+        self._spthermo.reset(CxxNewSpeciesThermo(self.derived_type, T_low,
+                                                 T_high, P_ref, &data[0]))
+        self.spthermo = self._spthermo.get()
+
+    cdef _assign(self, shared_ptr[CxxSpeciesThermo] other):
+        self._spthermo = other
+        self.spthermo = self._spthermo.get()
+
+    def cp(self, T):
+        """ Molar heat capacity at constant pressure [J/kmol/K] """
+        cdef double cp_r, h_rt, s_r
+        self.spthermo.updatePropertiesTemp(T, &cp_r, &h_rt, &s_r)
+        return cp_r * gas_constant
+
+    def h(self, T):
+        """ Molar enthalpy [J/kmol] """
+        cdef double cp_r, h_rt, s_r
+        self.spthermo.updatePropertiesTemp(T, &cp_r, &h_rt, &s_r)
+        return h_rt * gas_constant * T
+
+    def s(self, T):
+        """ Molar entropy [J/kmol/K] """
+        cdef double cp_r, h_rt, s_r
+        self.spthermo.updatePropertiesTemp(T, &cp_r, &h_rt, &s_r)
+        return s_r * gas_constant
+
+
+cdef class ConstantCp(SpeciesThermo):
+    """
+    Thermodynamic properties for a species that has a constant specific heat
+    capacity. This is a wrapper for the C++ class :ct:`ConstCpPoly`.
+    """
+    derived_type = SPECIES_THERMO_CONSTANT_CP
+    n_coeffs = 4
+
+
+cdef class NasaPoly2(SpeciesThermo):
+    """
+    Thermodynamic properties for a species which is parameterized using the
+    7-coefficient NASA polynomial form in two temperature ranges. This is a
+    wrapper for the C++ class :ct:`NasaPoly2`.
+    """
+    derived_type = SPECIES_THERMO_NASA2
+    n_coeffs = 15
+
+
+cdef class ShomatePoly2(SpeciesThermo):
+    """
+    Thermodynamic properties for a species which is parameterized using the
+    Shomate equation in two temperature ranges. This is a wrapper for the C++
+    class :ct:`ShomatePoly2`.
+    """
+    derived_type = SPECIES_THERMO_SHOMATE2
+    n_coeffs = 15
+
+
+cdef wrapSpeciesThermo(shared_ptr[CxxSpeciesThermo] spthermo):
+    """
+    Wrap a C++ SpeciesThermoInterpType object with a Python object of the
+    correct derived type.
+    """
+    cdef int thermo_type = spthermo.get().reportType()
+
+    if thermo_type == SPECIES_THERMO_NASA2:
+        st = NasaPoly2(init=False)
+    elif thermo_type == SPECIES_THERMO_CONSTANT_CP:
+        st = ConstantCp(init=False)
+    elif thermo_type == SPECIES_THERMO_SHOMATE2:
+        st = ShomatePoly2(init=False)
+    else:
+        st = SpeciesThermo()
+
+    st._assign(spthermo)
+    return st

interfaces/cython/cantera/test/test_thermo.py

@@ -657,3 +657,39 @@ class TestSpecies(utilities.CanteraTest):
         for name in gas.species_names:
             s = gas.species(name)
             self.assertEqual(s.name, name)
+
+
+class TestSpeciesThermo(utilities.CanteraTest):
+    def setUp(self):
+        self.gas = ct.Solution('h2o2.xml')
+        self.gas.X = 'H2O:1.0'
+
+    def test_create(self):
+        st = ct.NasaPoly2(300, 3500, 101325,
+                [1000.0, 3.03399249E+00, 2.17691804E-03, -1.64072518E-07,
+                 -9.70419870E-11, 1.68200992E-14, -3.00042971E+04, 4.96677010E+00,
+                 4.19864056E+00, -2.03643410E-03, 6.52040211E-06, -5.48797062E-09,
+                 1.77197817E-12, -3.02937267E+04, -8.49032208E-01])
+
+        for T in [300, 500, 900, 1200, 2000]:
+            self.gas.TP = T, 101325
+            self.assertAlmostEqual(st.cp(T), self.gas.cp_mole)
+            self.assertAlmostEqual(st.h(T), self.gas.enthalpy_mole)
+            self.assertAlmostEqual(st.s(T), self.gas.entropy_mole)
+
+    def test_invalid(self):
+        with self.assertRaises(ValueError):
+            # not enough coefficients
+            st = ct.NasaPoly2(300, 3500, 101325,
+                              [1000.0, 3.03399249E+00, 2.17691804E-03])
+
+    def test_wrap(self):
+        st = self.gas.species('H2O').thermo
+
+        self.assertTrue(isinstance(st, ct.NasaPoly2))
+
+        for T in [300, 500, 900, 1200, 2000]:
+            self.gas.TP = T, 101325
+            self.assertAlmostEqual(st.cp(T), self.gas.cp_mole)
+            self.assertAlmostEqual(st.h(T), self.gas.enthalpy_mole)
+            self.assertAlmostEqual(st.s(T), self.gas.entropy_mole)

interfaces/cython/cantera/thermo.pyx

@@ -54,6 +54,13 @@ cdef class Species:
         def __get__(self):
             return self.species.size
 
+    property thermo:
+        def __get__(self):
+            return wrapSpeciesThermo(self.species.thermo)
+
+        def __set__(self, SpeciesThermo spthermo):
+            self.species.thermo = spthermo._spthermo
+
 
 cdef class ThermoPhase(_SolutionBase):
     """