[Thermo] Add ability to modify species data for existing Phase objects

include/cantera/thermo/GeneralSpeciesThermo.h

@@ -37,6 +37,9 @@ public:
     virtual void install_STIT(size_t index,
                               shared_ptr<SpeciesThermoInterpType> stit_ptr);
 
+    virtual void modifySpecies(size_t index,
+                               shared_ptr<SpeciesThermoInterpType> stit_ptr);
+
     //! Install a PDSS object to handle the reference state thermodynamics
     //! calculation
     /*!
@@ -101,6 +104,9 @@ protected:
     //! Temperature polynomials for each thermo parameterization
     mutable tpoly_map m_tpoly;
 
+    //! Map from species index to location within #m_sp, such that
+    //! `m_sp[m_speciesLoc[k].first][m_speciesLoc[k].second]` is the
+    //! SpeciesThermoInterpType object for species `k`.
     std::map<size_t, std::pair<int, size_t> > m_speciesLoc;
 
     //! Maximum value of the lowest temperature

include/cantera/thermo/Phase.h

@@ -689,6 +689,15 @@ public:
     //! @see ignoreUndefinedElements addUndefinedElements throwUndefinedElements
     virtual bool addSpecies(shared_ptr<Species> spec);
 
+    //! Modify the thermodynamic data associated with a species.
+    /*!
+     * The species name, elemental composition, and type of thermo
+     * parameterization must be unchanged. If there are Kinetics objects that
+     * depend on this phase, Kinetics::invalidateCache() should be called on
+     * those objects after calling this function.
+     */
+    virtual void modifySpecies(size_t k, shared_ptr<Species> spec);
+
     //! Return the Species object for the named species.
     shared_ptr<Species> species(const std::string& name) const;
 

include/cantera/thermo/SpeciesThermo.h

@@ -142,6 +142,14 @@ public:
     virtual void install_STIT(size_t index,
                               shared_ptr<SpeciesThermoInterpType> stit) = 0;
 
+    //! Modify the species thermodynamic property parameterization for a species
+    /*!
+     * @param index Index of the species being installed
+     * @param spec  Pointer to the SpeciesThermoInterpType object
+     */
+    virtual void modifySpecies(size_t index,
+                               shared_ptr<SpeciesThermoInterpType> spec) = 0;
+
     //! Compute the reference-state properties for all species.
     /*!
      * Given temperature T in K, this method updates the values of the non-

include/cantera/thermo/ThermoPhase.h

@@ -1384,6 +1384,8 @@ public:
     using Phase::addSpecies;
     virtual bool addSpecies(shared_ptr<Species> spec);
 
+    virtual void modifySpecies(size_t k, shared_ptr<Species> spec);
+
     //! Store a reference pointer to the XML tree containing the species data
     //! for this phase.
     /*!

interfaces/cython/cantera/_cantera.pxd

@@ -123,6 +123,7 @@ cdef extern from "cantera/thermo/ThermoPhase.h" namespace "Cantera":
         # initialization
         void addUndefinedElements() except +
         cbool addSpecies(shared_ptr[CxxSpecies]) except +
+        void modifySpecies(size_t, shared_ptr[CxxSpecies]) except +
         void initThermo() except +
         void invalidateCache() except +
 

interfaces/cython/cantera/test/test_kinetics.py

@@ -344,6 +344,32 @@ class KineticsRepeatability(utilities.CanteraTest):
     def test_pdep_pressure(self):
         self.check_rates_pressure('pdep-test.xml')
 
+    def test_modify_thermo(self):
+        # Make sure that thermo modifications propagate through to Kinetics
+
+        # Set a gas state that is near enough to equilibrium that changes in the
+        # reverse rate always show up in the net rate
+        gas = self.setup_gas('gri30.xml')
+        gas.TPX = self.T0, self.P0, self.X0
+        gas.equilibrate('TP')
+        gas.TP = gas.T + 20, None
+
+        S = {sp.name: sp for sp in ct.Species.listFromFile('gri30.xml')}
+        w1 = gas.net_rates_of_progress
+
+        OH = gas.species('OH')
+        OH.thermo = S['CO2'].thermo
+        gas.modify_species(gas.species_index('OH'), OH)
+        w2 = gas.net_rates_of_progress
+
+        for i,R in enumerate(gas.reactions()):
+            if ('OH' in R.reactants or 'OH' in R.products) and R.reversible:
+                # Rate should be different if reaction involves OH
+                self.assertNotAlmostEqual(w2[i] / w1[i], 1.0)
+            else:
+                # Rate should be the same if reaction does not involve OH
+                self.assertAlmostEqual(w2[i] / w1[i], 1.0)
+
 
 class TestEmptyKinetics(utilities.CanteraTest):
     def test_empty(self):

interfaces/cython/cantera/test/test_thermo.py

@@ -871,6 +871,44 @@ class TestSpecies(utilities.CanteraTest):
         self.assertEqual({sp.name for sp in S},
                          set(self.gas.species_names))
 
+    def test_modify_thermo(self):
+        S = {sp.name: sp for sp in ct.Species.listFromFile('h2o2.xml')}
+        self.gas.TPX = 400, 2*ct.one_atm, 'H2:1.0'
+        g0 = self.gas.gibbs_mole
+
+        self.gas.TPX = None, None, 'O2:1.0'
+        self.assertNotAlmostEqual(g0, self.gas.gibbs_mole)
+        # Replace O2 thermo with the data from H2
+        S['O2'].thermo = S['H2'].thermo
+        self.gas.modify_species(self.gas.species_index('O2'), S['O2'])
+        self.assertAlmostEqual(g0, self.gas.gibbs_mole)
+
+    def test_modify_thermo_invalid(self):
+        S = {sp.name: sp for sp in ct.Species.listFromFile('h2o2.xml')}
+
+        orig = S['H2']
+        thermo = orig.thermo
+        copy = ct.Species('foobar', orig.composition)
+        copy.thermo = thermo
+        with self.assertRaises(Exception):
+            self.gas.modify_species(self.gas.species_index('H2'), copy)
+
+        copy = ct.Species('H2', {'H': 3})
+        copy.thermo = thermo
+        with self.assertRaises(Exception):
+            self.gas.modify_species(self.gas.species_index('H2'), copy)
+
+        copy = ct.Species('H2', orig.composition)
+        copy.thermo = ct.ConstantCp(thermo.min_temp, thermo.max_temp,
+            thermo.reference_pressure, [300, 123, 456, 789])
+        with self.assertRaises(Exception):
+            self.gas.modify_species(self.gas.species_index('H2'), copy)
+
+        copy = ct.Species('H2', orig.composition)
+        copy.thermo = ct.NasaPoly2(thermo.min_temp+200, thermo.max_temp,
+            thermo.reference_pressure, thermo.coeffs)
+        with self.assertRaises(Exception):
+            self.gas.modify_species(self.gas.species_index('H2'), copy)
 
 
 class TestSpeciesThermo(utilities.CanteraTest):

interfaces/cython/cantera/thermo.pyx

@@ -432,6 +432,11 @@ cdef class ThermoPhase(_SolutionBase):
                             " Got {!r}.".format(k))
         return s
 
+    def modify_species(self, k, Species species):
+        self.thermo.modifySpecies(k, species._species)
+        if self.kinetics:
+            self.kinetics.invalidateCache()
+
     def n_atoms(self, species, element):
         """
         Number of atoms of element *element* in species *species*. The element

src/thermo/GeneralSpeciesThermo.cpp

@@ -96,6 +96,38 @@ void GeneralSpeciesThermo::install_STIT(size_t index,
     markInstalled(index);
 }
 
+void GeneralSpeciesThermo::modifySpecies(size_t index,
+                                         shared_ptr<SpeciesThermoInterpType> spthermo)
+{
+    if (!spthermo) {
+        throw CanteraError("GeneralSpeciesThermo::modifySpecies",
+                           "null pointer");
+    }
+    if (m_speciesLoc.find(index) == m_speciesLoc.end()) {
+        throw CanteraError("GeneralSpeciesThermo::modifySpecies",
+                           "Species with this index not previously added: {}",
+                           index);
+    }
+    int type = spthermo->reportType();
+    if (m_speciesLoc[index].first != type) {
+        throw CanteraError("GeneralSpeciesThermo::modifySpecies",
+                           "Type of parameterization changed: {} != {}", type,
+                           m_speciesLoc[index].first);
+    }
+    if (spthermo->minTemp() > m_tlow_max) {
+        throw CanteraError("GeneralSpeciesThermo::modifySpecies",
+            "Cannot increase minimum temperature for phase from {} to {}",
+            m_tlow_max, spthermo->minTemp());
+    }
+    if (spthermo->maxTemp() < m_thigh_min) {
+        throw CanteraError("GeneralSpeciesThermo::modifySpecies",
+            "Cannot increase minimum temperature for phase from {} to {}",
+            m_thigh_min, spthermo->maxTemp());
+    }
+
+    m_sp[type][m_speciesLoc[index].second] = {index, spthermo};
+}
+
 void GeneralSpeciesThermo::installPDSShandler(size_t k, PDSS* PDSS_ptr,
         VPSSMgr* vpssmgr_ptr)
 {

src/thermo/Phase.cpp

@@ -837,6 +837,23 @@ bool Phase::addSpecies(shared_ptr<Species> spec) {
     return true;
 }
 
+void Phase::modifySpecies(size_t k, shared_ptr<Species> spec)
+{
+    if (speciesName(k) != spec->name) {
+        throw CanteraError("Phase::modifySpecies",
+            "New species name '{}' does not match existing name '{}'",
+                           spec->name, speciesName(k));
+    }
+    const shared_ptr<Species>& old = m_species[spec->name];
+    if (spec->composition != old->composition) {
+        throw CanteraError("Phase::modifySpecies",
+            "New composition for '{}' does not match existing composition",
+            spec->name);
+    }
+    m_species[spec->name] = spec;
+    invalidateCache();
+}
+
 shared_ptr<Species> Phase::species(const std::string& name) const
 {
     return getValue(m_species, name);

src/thermo/ThermoPhase.cpp

@@ -715,6 +715,18 @@ bool ThermoPhase::addSpecies(shared_ptr<Species> spec)
     return added;
 }
 
+void ThermoPhase::modifySpecies(size_t k, shared_ptr<Species> spec)
+{
+    Phase::modifySpecies(k, spec);
+    if (speciesName(k) != spec->name) {
+        throw CanteraError("ThermoPhase::modifySpecies",
+            "New species '{}' does not match existing species '{}' at index {}",
+                           spec->name, speciesName(k), k);
+    }
+    spec->thermo->validate(spec->name);
+    m_spthermo->modifySpecies(k, spec->thermo);
+}
+
 void ThermoPhase::saveSpeciesData(const size_t k, const XML_Node* const data)
 {
     if (m_speciesData.size() < (k + 1)) {