Expose C++ Units to Python

include/cantera/base/Units.h

@@ -35,7 +35,8 @@ class Units
                    double quantity=0);
 
     //! Create an object with the specified dimensions
-    explicit Units(const std::string& name);
+    //! @param force_unity  ensure that conversion factor is equal to one
+    explicit Units(const std::string& name, bool force_unity=false);
 
     //! Returns `true` if the specified Units are dimensionally consistent
     bool convertible(const Units& other) const;
@@ -49,7 +50,8 @@ class Units
     Units& operator*=(const Units& other);
 
     //! Provide a string representation of these Units
-    std::string str() const;
+    //! @param skip_unity  do not print '1' if conversion factor is equal to one
+    std::string str(bool skip_unity=true) const;
 
     //! Raise these Units to a power, changing both the conversion factor and
     //! the dimensions of these Units.
@@ -109,6 +111,9 @@ class UnitSystem
     //! recognize an optional argument with a default value)
     UnitSystem() : UnitSystem({}) {}
 
+    //! Return default units used by the unit system
+    std::map<std::string, std::string> defaults() const;
+
     //! Set the default units to convert from when explicit units are not
     //! provided. Defaults can be set for mass, length, time, quantity, energy,
     //! and pressure. Conversion using the pressure or energy units is done only

include/cantera/base/YamlWriter.h

@@ -59,6 +59,12 @@ class YamlWriter
     //!     corresponding units supported by the UnitSystem class.
     void setUnits(const std::map<std::string, std::string>& units={});
 
+    //! Set the units to be used in the output file. Dimensions not specified
+    //! will use Cantera's defaults.
+    //! @param units  A UnitSystem object specifying dimensions (mass, length, time,
+    //!     quantity, pressure, energy, activation-energy).
+    void setUnitSystem(const UnitSystem& units=UnitSystem());
+
 protected:
     std::vector<shared_ptr<Solution>> m_phases;
 

interfaces/cython/cantera/_cantera.pxd

@@ -49,6 +49,19 @@ cdef extern from "cantera/base/xml.h" namespace "Cantera":
         XML_Node* findID(string)
         int nChildren()
 
+cdef extern from "cantera/base/Units.h" namespace "Cantera":
+    cdef cppclass CxxUnits "Cantera::Units":
+        CxxUnits()
+        CxxUnits(CxxUnits)
+        CxxUnits(string, cbool) except +translate_exception
+        string str()
+        double factor()
+
+    cdef cppclass CxxUnitSystem "Cantera::UnitSystem":
+        CxxUnitSystem()
+        stdmap[string, string] defaults()
+        void setDefaults(stdmap[string, string]&) except +translate_exception
+
 cdef extern from "cantera/base/AnyMap.h" namespace "Cantera":
     cdef cppclass CxxAnyValue "Cantera::AnyValue"
 
@@ -208,6 +221,7 @@ cdef extern from "cantera/thermo/ThermoPhase.h" namespace "Cantera":
         size_t stateSize()
         void saveState(size_t, double*) except +translate_exception
         void restoreState(size_t, double*) except +translate_exception
+        CxxUnits standardConcentrationUnits() except +translate_exception
 
         # initialization
         void addUndefinedElements() except +translate_exception
@@ -420,6 +434,7 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera":
         cbool allow_nonreactant_orders
         cbool allow_negative_orders
         cbool usesLegacy()
+        CxxUnits rate_units
 
     cdef cppclass CxxElementaryReaction2 "Cantera::ElementaryReaction2" (CxxReaction):
         CxxElementaryReaction2()
@@ -632,7 +647,7 @@ cdef extern from "cantera/base/YamlWriter.h" namespace "Cantera":
         void toYamlFile(string&) except +translate_exception
         void setPrecision(int)
         void skipUserDefined(cbool)
-        void setUnits(stdmap[string, string]&) except +translate_exception
+        void setUnitSystem(CxxUnitSystem&) except +translate_exception
 
 cdef extern from "cantera/equil/MultiPhase.h" namespace "Cantera":
     cdef cppclass CxxMultiPhase "Cantera::MultiPhase":
@@ -1108,6 +1123,14 @@ ctypedef void (*transportMethod2d)(CxxTransport*, size_t, double*) except +trans
 ctypedef void (*kineticsMethod1d)(CxxKinetics*, double*) except +translate_exception
 
 # classes
+cdef class Units:
+    cdef CxxUnits units
+    @staticmethod
+    cdef copy(CxxUnits)
+
+cdef class UnitSystem:
+    cdef CxxUnitSystem unitsystem
+
 cdef class SpeciesThermo:
     cdef shared_ptr[CxxSpeciesThermo] _spthermo
     cdef CxxSpeciesThermo* spthermo
@@ -1199,6 +1222,8 @@ cdef class DustyGasTransport(Transport):
 cdef class YamlWriter:
     cdef shared_ptr[CxxYamlWriter] _writer
     cdef CxxYamlWriter* writer
+    @staticmethod
+    cdef CxxUnitSystem _get_unitsystem(UnitSystem units)
 
 cdef class Mixture:
     cdef CxxMultiPhase* mix

interfaces/cython/cantera/_cantera.pyx

@@ -12,6 +12,7 @@ from ._cantera cimport *
 include "utils.pyx"
 include "constants.pyx"
 include "func1.pyx"
+include "units.pyx"
 
 include "base.pyx"
 include "speciesthermo.pyx"

interfaces/cython/cantera/base.pyx

@@ -263,8 +263,8 @@ cdef class _SolutionBase:
             Additional ThermoPhase / Solution objects to be included in the
             output file
         :param units:
-            A dictionary of the units to be used for each dimension. See
-            `YamlWriter.output_units`.
+            A `UnitSystem` object or dictionary of the units to be used for
+            each dimension. See `YamlWriter.output_units`.
         :param precision:
             For output floating point values, the maximum number of digits to
             the right of the decimal point. The default is 15 digits.

interfaces/cython/cantera/reaction.pyx

@@ -823,6 +823,12 @@ cdef class Reaction:
         def __get__(self):
             return self.reaction.usesLegacy()
 
+    property rate_coeff_units:
+        """Get reaction rate coefficient units"""
+        def __get__(self):
+            cdef CxxUnits rate_units = self.reaction.rate_units
+            return Units.copy(rate_units)
+
 
 cdef class Arrhenius:
     r"""

interfaces/cython/cantera/test/test_composite.py

@@ -580,7 +580,7 @@ def test_yaml_simple(self):
                              gas2.forward_rate_constants)
         self.assertArrayNear(gas.mix_diff_coeffs, gas2.mix_diff_coeffs)
 
-    def test_yaml_outunits(self):
+    def test_yaml_outunits1(self):
         gas = ct.Solution('h2o2.yaml')
         gas.TPX = 500, ct.one_atm, 'H2: 1.0, O2: 1.0'
         gas.equilibrate('HP')
@@ -604,6 +604,30 @@ def test_yaml_outunits(self):
                              gas2.forward_rate_constants)
         self.assertArrayNear(gas.mix_diff_coeffs, gas2.mix_diff_coeffs)
 
+    def test_yaml_outunits2(self):
+        gas = ct.Solution('h2o2.yaml')
+        gas.TPX = 500, ct.one_atm, 'H2: 1.0, O2: 1.0'
+        gas.equilibrate('HP')
+        gas.TP = 1500, ct.one_atm
+        units = {'length': 'cm', 'quantity': 'mol', 'energy': 'cal'}
+        system = ct.UnitSystem(units)
+        gas.write_yaml('h2o2-generated.yaml', units=system)
+        generated = utilities.load_yaml("h2o2-generated.yaml")
+        original = utilities.load_yaml(self.cantera_data_path / "h2o2.yaml")
+
+        for r1, r2 in zip(original['reactions'], generated['reactions']):
+            if 'rate-constant' in r1:
+                self.assertNear(r1['rate-constant']['A'], r2['rate-constant']['A'])
+                self.assertNear(r1['rate-constant']['Ea'], r2['rate-constant']['Ea'])
+
+        gas2 = ct.Solution("h2o2-generated.yaml")
+        self.assertArrayNear(gas.concentrations, gas2.concentrations)
+        self.assertArrayNear(gas.partial_molar_enthalpies,
+                             gas2.partial_molar_enthalpies)
+        self.assertArrayNear(gas.forward_rate_constants,
+                             gas2.forward_rate_constants)
+        self.assertArrayNear(gas.mix_diff_coeffs, gas2.mix_diff_coeffs)
+
     def test_yaml_surface(self):
         gas = ct.Solution('ptcombust.yaml', 'gas')
         surf = ct.Interface('ptcombust.yaml', 'Pt_surf', [gas])

interfaces/cython/cantera/test/test_utils.py

@@ -5,6 +5,63 @@
 
 from cantera._cantera import _py_to_any_to_py
 
+
+class TestUnitSystem(utilities.CanteraTest):
+
+    def test_default(self):
+        units = ct.UnitSystem().units
+        checks = {
+            "activation-energy": "J / kmol",
+            "current": "A",
+            "energy": "J",
+            "length": "m",
+            "mass": "kg",
+            "pressure": "Pa",
+            "quantity": "kmol",
+            "temperature": "K",
+            "time": "s",
+        }
+        for dim, unit in units.items():
+            self.assertIn(dim, checks)
+            self.assertEqual(checks[dim], unit)
+
+    def test_cgs(self):
+        system = ct.UnitSystem({
+            "length": "cm", "mass": "g", "time": "s",
+            "quantity": "mol", "pressure": "dyn / cm^2", "energy": "erg",
+            "activation-energy": "cal / mol"})
+        units = system.units
+        checks = {
+            "activation-energy": "cal / mol",
+            "current": "A",
+            "energy": "erg",
+            "length": "cm",
+            "mass": "g",
+            "pressure": "dyn / cm^2",
+            "quantity": "mol",
+            "temperature": "K",
+            "time": "s",
+        }
+        for dim, unit in units.items():
+            self.assertIn(dim, checks)
+            self.assertEqual(checks[dim], unit)
+
+    def test_activation_energy(self):
+        system = ct.UnitSystem({"activation-energy": "eV"})
+        units = system.units
+        self.assertEqual(units["activation-energy"], "eV")
+
+        system = ct.UnitSystem({"activation-energy": "K"})
+        units = system.units
+        self.assertEqual(units["activation-energy"], "K")
+
+    def test_raises(self):
+        with self.assertRaisesRegex(ct.CanteraError, "non-unity conversion factor"):
+            ct.UnitSystem({"temperature": "2 K"})
+        with self.assertRaisesRegex(ct.CanteraError, "non-unity conversion factor"):
+            ct.UnitSystem({"current": "2 A"})
+
+
 class TestPyToAnyValue(utilities.CanteraTest):
 
     def check_conversion(self, value, check_type=None):

interfaces/cython/cantera/thermo.pyx

@@ -1570,6 +1570,12 @@ cdef class ThermoPhase(_SolutionBase):
         def __set__(self, double value):
             self.thermo.setElectricPotential(value)
 
+    property standard_concentration_units:
+        """Get standard concentration units for this phase."""
+        def __get__(self):
+            cdef CxxUnits units = self.thermo.standardConcentrationUnits()
+            return Units.copy(units)
+
 
 cdef class InterfacePhase(ThermoPhase):
     """ A class representing a surface or edge phase"""

interfaces/cython/cantera/units.pyx

@@ -0,0 +1,80 @@
+# This file is part of Cantera. See License.txt in the top-level directory or
+# at https://cantera.org/license.txt for license and copyright information.
+
+cdef class Units:
+    """
+    A representation of the units associated with a dimensional quantity.
+
+    This class is a light-weight interface to internal Cantera capabilities that are
+    used for converting quantities between unit systems and checking for dimensional
+    consistency. Internally, `Units` objects are mainly used within the `UnitSystem`
+    class to convert values from a user-specified Unit system to Cantera's base units
+    (SI + kmol).
+
+    The Python API handles display of `Units` that do not require a conversion factor,
+    with other functions not enabled.
+    """
+    def __cinit__(self, name=None):
+        if name:
+            self.units = CxxUnits(stringify(name), True)
+
+    def __repr__(self):
+        return f"<Units({pystr(self.units.str())}) at {id(self):0x}>"
+
+    @staticmethod
+    cdef copy(CxxUnits other):
+        """Copy a C++ Units object to a Python object."""
+        cdef Units units = Units()
+        units.units = CxxUnits(other)
+        return units
+
+
+cdef class UnitSystem:
+    """
+    Unit system used for YAML input and output.
+
+    The `UnitSystem` class is used to specify dimensional values for a given unit
+    system. The main use is for converting values specified in input files to Cantera's
+    native unit system, which is SI units except for the use of kmol as the base
+    unit of quantity, i.e. kilogram, meter, second, kelvin, ampere, and kmol.
+
+    The default unit system used by Cantera is SI+kmol::
+
+        ct.UnitSystem({
+            "length": "m", "mass": "kg", "time": "s",
+            "quantity": "kmol", "pressure": "Pa", "energy": "J",
+            "temperature": "K", "current": "A", "activation-energy": "J / kmol"})
+
+    A CGS+mol unit system with activation energy units of cal/mol is created as::
+
+        ct.UnitSystem({
+            "length": "cm", "mass": "g", "time": "s",
+            "quantity": "mol", "pressure": "dyn / cm^2", "energy": "erg",
+            "temperature": "K", "current": "A", "activation-energy": "cal / mol"})
+
+    Defaults for dimensions not specified will be left unchanged. Accordingly,
+    the default unit system is retrieved as::
+
+        ct.UnitSystem()
+    """
+    def __cinit__(self, units=None):
+        self.unitsystem = CxxUnitSystem()
+        if units:
+            self.units = units
+
+    def __repr__(self):
+        return f"<UnitSystem at {id(self):0x}>"
+
+    property units:
+        """
+        Units used by the unit system
+        """
+        def __get__(self):
+            cdef stdmap[string, string] cxxunits = self.unitsystem.defaults()
+            cdef pair[string, string] item
+            return {pystr(item.first): pystr(item.second) for item in cxxunits}
+        def __set__(self, units):
+            cdef stdmap[string, string] cxxunits
+            for dimension, unit in units.items():
+                cxxunits[stringify(dimension)] = stringify(unit)
+            self.unitsystem.setDefaults(cxxunits)

interfaces/cython/cantera/yamlwriter.pyx

@@ -50,15 +50,19 @@ cdef class YamlWriter:
         will use Cantera's defaults.
 
         :param units:
-            A map where keys are dimensions (mass, length, time, quantity,
-            pressure, energy, activation-energy), and the values are
+            A `UnitSystem` object or map where keys are dimensions (mass, length, time,
+            quantity, pressure, energy, activation-energy), and the values are
             corresponding units such as kg, mm, s, kmol, Pa, cal, and eV.
         """
         def __set__(self, units):
-            cdef stdmap[string, string] cxxunits
-            for dimension, unit in units.items():
-                cxxunits[stringify(dimension)] = stringify(unit)
-            self.writer.setUnits(cxxunits)
+            if not isinstance(units, UnitSystem):
+                units = UnitSystem(units)
+            cdef CxxUnitSystem cxxunits = YamlWriter._get_unitsystem(units)
+            self.writer.setUnitSystem(cxxunits)
+
+    @staticmethod
+    cdef CxxUnitSystem _get_unitsystem(UnitSystem units):
+        return units.unitsystem
 
     def __reduce__(self):
         raise NotImplementedError('YamlWriter object is not picklable')