[Python] Prevent adding species to in-use Thermo objects

interfaces/cython/cantera/_cantera.pxd

@@ -874,6 +874,7 @@ cdef class ThermoPhase(_SolutionBase):
     cpdef int species_index(self, species) except *
     cdef np.ndarray _getArray1(self, thermoMethod1d method)
     cdef void _setArray1(self, thermoMethod1d method, values) except *
+    cdef public object _references
 
 cdef class InterfacePhase(ThermoPhase):
     cdef CxxSurfPhase* surf
@@ -906,6 +907,7 @@ cdef class DustyGasTransport(Transport):
 cdef class Mixture:
     cdef CxxMultiPhase* mix
     cdef list _phases
+    cdef object _weakref_proxy
     cpdef int element_index(self, element) except *
 
 cdef class Func1:
@@ -919,6 +921,7 @@ cdef class ReactorBase:
     cdef list _inlets
     cdef list _outlets
     cdef list _walls
+    cdef object _weakref_proxy
 
 cdef class Reactor(ReactorBase):
     cdef CxxReactor* reactor
@@ -982,6 +985,7 @@ cdef class ReactorNet:
 cdef class Domain1D:
     cdef CxxDomain1D* domain
     cdef _SolutionBase gas
+    cdef object _weakref_proxy
     cdef public pybool have_user_tolerances
 
 cdef class Boundary1D(Domain1D):

interfaces/cython/cantera/mixture.pyx

@@ -1,5 +1,9 @@
 import warnings
 
+# Need a pure-python class to store weakrefs to
+class _WeakrefProxy(object):
+    pass
+
 cdef class Mixture:
     """
 
@@ -32,6 +36,7 @@ cdef class Mixture:
     def __cinit__(self, phases):
         self.mix = new CxxMultiPhase()
         self._phases = []
+        self._weakref_proxy = _WeakrefProxy()
 
         cdef _SolutionBase phase
         if isinstance(phases[0], _SolutionBase):
@@ -39,6 +44,7 @@ cdef class Mixture:
             phases = [(p, 1 if i == 0 else 0) for i,p in enumerate(phases)]
 
         for phase,moles in phases:
+            phase._references[self._weakref_proxy] = True
             self.mix.addPhase(phase.thermo, moles)
             self._phases.append(phase)
 

interfaces/cython/cantera/onedim.pyx

@@ -1,11 +1,16 @@
 import interrupts
 
+# Need a pure-python class to store weakrefs to
+class _WeakrefProxy(object):
+    pass
+
 cdef class Domain1D:
     def __cinit__(self, *args, **kwargs):
         self.domain = NULL
 
     # The signature of this function causes warnings for Sphinx documentation
     def __init__(self, _SolutionBase phase, *args, name=None, **kwargs):
+        self._weakref_proxy = _WeakrefProxy()
         if self.domain is NULL:
             raise TypeError("Can't instantiate abstract class Domain1D.")
 
@@ -13,6 +18,7 @@ cdef class Domain1D:
             self.name = name
 
         self.gas = phase
+        self.gas._references[self._weakref_proxy] = True
         self.have_user_tolerances = False
 
     property index:
@@ -407,6 +413,8 @@ cdef class _FlowBase(Domain1D):
         """
         Set the `Solution` object used for calculating transport properties.
         """
+        self._weakref_proxy = _WeakrefProxy()
+        self.gas._references[self._weakref_proxy] = True
         self.gas = phase
         self.flow.setTransport(deref(self.gas.transport))
 

interfaces/cython/cantera/reactor.pyx

@@ -3,6 +3,10 @@ import numbers as _numbers
 
 _reactor_counts = _defaultdict(int)
 
+# Need a pure-python class to store weakrefs to
+class _WeakrefProxy(object):
+    pass
+
 cdef class ReactorBase:
     """
     Common base class for reactors and reservoirs.
@@ -13,6 +17,7 @@ cdef class ReactorBase:
 
     # The signature of this function causes warnings for Sphinx documentation
     def __init__(self, ThermoPhase contents=None, name=None, *, volume=None):
+        self._weakref_proxy = _WeakrefProxy()
         self._inlets = []
         self._outlets = []
         self._walls = []
@@ -38,6 +43,7 @@ cdef class ReactorBase:
         properties and kinetic rates for this reactor.
         """
         self._thermo = solution
+        self._thermo._references[self._weakref_proxy] = True
         self.rbase.setThermoMgr(deref(solution.thermo))
 
     property name:

interfaces/cython/cantera/test/test_thermo.py

@@ -1,5 +1,6 @@
 from .utilities import unittest
 import numpy as np
+import gc
 
 import cantera as ct
 from . import utilities
@@ -596,6 +597,30 @@ class TestThermoPhase(utilities.CanteraTest):
         self.assertArrayNear(ref[self.phase.species_names].partial_molar_cp,
                              self.phase.partial_molar_cp)
 
+    def test_add_species_disabled(self):
+        ref = ct.Solution('gri30.xml')
+
+        reactor = ct.IdealGasReactor(self.phase)
+        with self.assertRaises(RuntimeError):
+            self.phase.add_species(ref.species('CH4'))
+        del reactor
+        gc.collect()
+        self.phase.add_species(ref.species('CH4'))
+
+        flame = ct.FreeFlame(self.phase, width=0.1)
+        with self.assertRaises(RuntimeError):
+            self.phase.add_species(ref.species('CO'))
+        del flame
+        gc.collect()
+        self.phase.add_species(ref.species('CO'))
+
+        mix = ct.Mixture([(self.phase, 2.0)])
+        with self.assertRaises(RuntimeError):
+            self.phase.add_species(ref.species('CH2O'))
+        del mix
+        gc.collect()
+        self.phase.add_species(ref.species('CH2O'))
+
 
 class TestThermo(utilities.CanteraTest):
     def setUp(self):

interfaces/cython/cantera/thermo.pyx

@@ -1,4 +1,5 @@
 import warnings
+import weakref
 
 cdef enum ThermoBasis:
     mass_basis = 0
@@ -223,6 +224,7 @@ cdef class ThermoPhase(_SolutionBase):
         super().__init__(*args, **kwargs)
         if 'source' not in kwargs:
             self.thermo_basis = mass_basis
+        self._references = weakref.WeakKeyDictionary()
 
     def report(self, show_thermo=True, float threshold=1e-14):
         """
@@ -455,6 +457,9 @@ cdef class ThermoPhase(_SolutionBase):
         Add a new species to this phase. Missing elements will be added
         automatically.
         """
+        if self._references:
+            raise RuntimeError('Cannot add species to ThermoPhase object if it'
+                ' is linked to a Reactor, Domain1D (flame), or Mixture object.')
         self.thermo.addUndefinedElements()
         self.thermo.addSpecies(species._species)
         self.thermo.initThermo()