Fix issues with Quantity class

interfaces/cython/cantera/composite.py

@@ -228,7 +228,7 @@ class Quantity:
     __slots__ = ("state", "_phase", "_id", "mass", "constant", "_weakref_proxy")
 
     def __init__(self, phase, mass=None, moles=None, constant='UV'):
-        self.state = phase.TDY
+        self.state = phase.state
         self._phase = phase
         self._weakref_proxy = _WeakrefProxy()
         phase._references[self._weakref_proxy] = True
@@ -253,7 +253,7 @@ def phase(self):
         Get the underlying `Solution` object, with the state set to match the
         wrapping `Quantity` object.
         """
-        self._phase.TDY = self.state
+        self._phase.state = self.state
         return self._phase
 
     @property
@@ -300,7 +300,21 @@ def equilibrate(self, XY=None, *args, **kwargs):
         if XY is None:
             XY = self.constant
         self.phase.equilibrate(XY, *args, **kwargs)
-        self.state = self._phase.TDY
+        self.state = self._phase.state
+
+    def set_equivalence_ratio(self, phi, fuel, oxidizer, basis="mole", *, diluent=None,
+                              fraction=None):
+        self._phase.state = self.state
+        self._phase.set_equivalence_ratio(phi, fuel, oxidizer, basis, diluent=diluent,
+                                          fraction=fraction)
+        self.state = self._phase.state
+    set_equivalence_ratio.__doc__ = Solution.set_equivalence_ratio.__doc__
+
+    def set_mixture_fraction(self, mixture_fraction, fuel, oxidizer, basis='mole'):
+        self._phase.state = self.state
+        self._phase.set_mixture_fraction(mixture_fraction, fuel, oxidizer, basis)
+        self.state = self._phase.state
+    set_mixture_fraction.__doc__ = Solution.set_mixture_fraction.__doc__
 
     def __imul__(self, other):
         self.mass *= other
@@ -317,14 +331,31 @@ def __iadd__(self, other):
             raise ValueError('Cannot add Quantities with different phase '
                 'definitions.')
         assert self.constant == other.constant
-        a1,b1 = getattr(self.phase, self.constant)
-        a2,b2 = getattr(other.phase, self.constant)
+
         m = self.mass + other.mass
-        a = (a1 * self.mass + a2 * other.mass) / m
-        b = (b1 * self.mass + b2 * other.mass) / m
-        self._phase.Y = (self.Y * self.mass + other.Y * other.mass) / m
-        setattr(self._phase, self.constant, (a,b))
-        self.state = self._phase.TDY
+        Y = (self.Y * self.mass + other.Y * other.mass)
+        if self.constant == 'UV':
+            U = self.int_energy + other.int_energy
+            V = self.volume + other.volume
+            if self.basis == 'mass':
+                self._phase.UVY = U / m, V / m, Y
+            else:
+                n = self.moles + other.moles
+                self._phase.UVY = U / n, V / n, Y
+        else:  # self.constant == 'HP'
+            dp_rel = 2 * abs(self.P - other.P) / (self.P + other.P)
+            if dp_rel > 1.0e-7:
+                raise ValueError('Cannot add Quantities at constant pressure when'
+                    f'pressure is not equal ({self.P} != {other.P})')
+
+            H = self.enthalpy + other.enthalpy
+            if self.basis == 'mass':
+                self._phase.HPY = H / m, None, Y
+            else:
+                n = self.moles + other.moles
+                self._phase.HPY = H / n, None, Y
+
+        self.state = self._phase.state
         self.mass = m
         return self
 
@@ -347,7 +378,7 @@ def getter(self):
 
     def setter(self, value):
         setattr(self.phase, attr, value)
-        self.state = self._phase.TDY
+        self.state = self._phase.state
 
     return property(getter, setter, doc=getattr(Solution, attr).__doc__)
 
@@ -360,6 +391,8 @@ def f(self, *args, **kwargs):
 for _attr in dir(Solution):
     if _attr.startswith('_') or _attr in Quantity.__dict__ or _attr == 'state':
         continue
+    if _attr.startswith('set_unnormalized'):
+        continue
     else:
         _orig = getattr(Solution, _attr)
         if hasattr(_orig, "__call__"):

test/python/test_thermo.py

@@ -5,7 +5,7 @@
 from . import utilities
 from .utilities import allow_deprecated
 import pytest
-
+from pytest import approx
 
 class TestThermoPhase(utilities.CanteraTest):
     def setUp(self):
@@ -1610,6 +1610,7 @@ def setUpClass(cls):
 
     def setUp(self):
         self.gas.TPX = 300, 101325, 'O2:1.0, N2:3.76'
+        self.gas.basis = 'mass'
 
     def test_mass_moles(self):
         q1 = ct.Quantity(self.gas, mass=5)
@@ -1638,6 +1639,33 @@ def test_extensive(self):
         self.assertNear(q1.entropy, q1.S)
         self.assertNear(q1.gibbs, q1.G)
 
+    def test_set_equivalence_ratio(self):
+        q1 = ct.Quantity(self.gas, mass=3)
+        T1, P1 = q1.TP
+        q1.set_equivalence_ratio(2.0, 'CH4:1.0', 'O2:1.0, N2:3.76')
+        assert q1.T == approx(T1)
+        assert q1.P == approx(P1)
+        assert q1.X[q1.species_index('CH4')] == approx(1.0 / (1 + 4.76))
+
+    def test_set_mixture_fraction(self):
+        q1 = ct.Quantity(self.gas, mass=3)
+        T1, P1 = q1.TP
+        q1.set_mixture_fraction(1.0, 'CH3OH:1.0', 'O2:1.0, N2:3.76')
+        assert q1.T == approx(T1)
+        assert q1.P == approx(P1)
+        assert q1.mass == 3
+        assert q1.X[q1.species_index('CH3OH')] == approx(1.0)
+
+    def test_basis(self):
+        q1 = ct.Quantity(self.gas, mass=5)
+        T1, P1 = q1.TP
+        h1 = q1.h  # J/kg
+
+        q1.basis = 'molar'
+        assert q1.T == approx(T1)
+        assert q1.P == approx(P1)
+        assert q1.h == approx(h1 * q1.mean_molecular_weight)
+
     def test_multiply(self):
         q1 = ct.Quantity(self.gas, mass=5)
         q2 = q1 * 2.5
@@ -1681,6 +1709,38 @@ def test_add(self):
         self.assertNear(q1.V + q2.V, q3.V)
         self.assertArrayNear(q1.X*q1.moles + q2.X*q2.moles, q3.X*q3.moles)
 
+    def test_add_molar(self):
+        q1 = ct.Quantity(self.gas, mass=5)
+        q2 = ct.Quantity(self.gas, mass=5)
+        q2.TPX = 900, 101325, 'CH4:1.0'
+
+        q1.basis = 'molar'
+        q2.basis = 'molar'
+
+        # addition at constant UV
+        q3 = q1 + q2
+        assert q1.mass + q2.mass == approx(q3.mass)
+
+        assert q1.U + q2.U == approx(q3.U)
+        assert q1.V + q2.V == approx(q3.V)
+
+        # addition at constant HP
+        q1.constant = q2.constant = 'HP'
+        q4 = q1 + q2
+        assert q1.mass + q2.mass == approx(q4.mass)
+
+        assert q1.H + q2.H == approx(q4.H)
+        assert q4.P == approx(q1.P)
+        self.assertArrayNear(q1.X*q1.moles + q2.X*q2.moles, q4.X*q4.moles)
+
+    def test_add_errors(self):
+        q1 = ct.Quantity(self.gas, mass=5)
+        q2 = ct.Quantity(self.gas, mass=5)
+        q1.constant = q2.constant = 'HP'
+        q2.TP = q1.T, 1.2 * q1.P
+        with pytest.raises(ValueError, match="pressure is not equal"):
+            q3 = q1 + q2
+
     def test_equilibrate(self):
         self.gas.TPX = 300, 101325, 'CH4:1.0, O2:0.2, N2:1.0'
         q1 = ct.Quantity(self.gas)
@@ -1696,6 +1756,12 @@ def test_invalid_setter(self):
         with self.assertRaises(AttributeError):
             q1.HPQ = self.gas.H, self.gas.P, 1
 
+        with pytest.raises(AttributeError):
+            q1.set_unnormalized_mass_fractions(np.ones(q1.n_species))
+
+        with pytest.raises(AttributeError):
+            q1.set_unnormalized_mole_fractions(np.ones(q1.n_species))
+
     def test_incompatible(self):
         gas2 = ct.Solution('h2o2.yaml', transport_model=None)
         q1 = ct.Quantity(self.gas)