[Kinetics] Allow adding species after reactions for homogeneous kinetics

This enables incremental mechanism construction for gas phase kinetics. For
surface kinetics, adding new species changes the kinetics species index of
existing species in other phases, so this feature is disabled.

include/cantera/kinetics/BulkKinetics.h

@@ -35,7 +35,7 @@ public:
                                      bool doIrreversible = false);
 
     virtual bool addReaction(shared_ptr<Reaction> r);
-    virtual void init();
+    virtual void resizeSpecies();
 
     virtual void setMultiplier(size_t i, double f);
     virtual void invalidateCache();

include/cantera/kinetics/InterfaceKinetics.h

@@ -198,6 +198,7 @@ public:
     virtual void addPhase(thermo_t& thermo);
 
     virtual void init();
+    virtual void resizeSpecies();
     virtual bool addReaction(shared_ptr<Reaction> r);
     virtual void modifyReaction(size_t i, shared_ptr<Reaction> rNew);
     //! @}

include/cantera/kinetics/Kinetics.h

@@ -714,15 +714,21 @@ public:
     virtual void addPhase(thermo_t& thermo);
 
     /**
-     * Prepare the class for the addition of reactions. This method is called
-     * by importKinetics() after all phases have been added but before any
-     * reactions have been. The base class method does nothing, but derived
-     * classes may use this to perform any initialization (allocating arrays,
-     * etc.) that requires knowing the phases and species, but before any
-     * reactions are added.
+     * Prepare the class for the addition of reactions, after all phases have
+     * been added. This method is called automatically when the first reaction
+     * is added. It needs to be called directly only in the degenerate case
+     * where there are no reactions. The base class method does nothing, but
+     * derived classes may use this to perform any initialization (allocating
+     * arrays, etc.) that requires knowing the phases.
      */
     virtual void init() {}
 
+    /**
+     * Resize arrays with sizes that depend on the total number of species.
+     * Automatically called before adding each Reaction and Phase.
+     */
+    virtual void resizeSpecies();
+
     /**
      * Finish adding reactions and prepare for use. This method is called by
      * importKinetics() after all reactions have been entered into the

src/kinetics/BulkKinetics.cpp

@@ -137,9 +137,9 @@ void BulkKinetics::modifyElementaryReaction(size_t i, ElementaryReaction& rNew)
     m_rates.replace(i, rNew.rate);
 }
 
-void BulkKinetics::init()
+void BulkKinetics::resizeSpecies()
 {
-    m_kk = thermo().nSpecies();
+    Kinetics::resizeSpecies();
     m_conc.resize(m_kk);
     m_grt.resize(m_kk);
 }

src/kinetics/InterfaceKinetics.cpp

@@ -828,20 +828,6 @@ void InterfaceKinetics::addPhase(thermo_t& thermo)
 
 void InterfaceKinetics::init()
 {
-    m_kk = 0;
-    for (size_t n = 0; n < nPhases(); n++) {
-        m_kk += thermo(n).nSpecies();
-    }
-    m_actConc.resize(m_kk);
-    m_conc.resize(m_kk);
-    m_StandardConc.resize(m_kk, 0.0);
-    m_mu0.resize(m_kk);
-    m_mu.resize(m_kk);
-    m_mu0_Kc.resize(m_kk);
-    m_grt.resize(m_kk);
-    m_pot.resize(m_kk, 0.0);
-    m_phi.resize(nPhases(), 0.0);
-
     size_t ks = reactionPhaseIndex();
     if (ks == npos) throw CanteraError("InterfaceKinetics::finalize",
                                            "no surface phase is present.");
@@ -855,6 +841,25 @@ void InterfaceKinetics::init()
     }
 }
 
+void InterfaceKinetics::resizeSpecies()
+{
+    size_t kOld = m_kk;
+    Kinetics::resizeSpecies();
+    if (m_kk != kOld && nReactions()) {
+        throw CanteraError("InterfaceKinetics::resizeSpecies", "Cannot add"
+            " species to InterfaceKinetics after reactions have been added.");
+    }
+    m_actConc.resize(m_kk);
+    m_conc.resize(m_kk);
+    m_StandardConc.resize(m_kk, 0.0);
+    m_mu0.resize(m_kk);
+    m_mu.resize(m_kk);
+    m_mu0_Kc.resize(m_kk);
+    m_grt.resize(m_kk);
+    m_pot.resize(m_kk, 0.0);
+    m_phi.resize(nPhases(), 0.0);
+}
+
 doublereal InterfaceKinetics::electrochem_beta(size_t irxn) const
 {
     for (size_t i = 0; i < m_ctrxn.size(); i++) {

src/kinetics/Kinetics.cpp

@@ -483,16 +483,6 @@ void Kinetics::getNetProductionRates(doublereal* net)
 
 void Kinetics::addPhase(thermo_t& thermo)
 {
-    // if not the first thermo object, set the start position
-    // to that of the last object added + the number of its species
-    if (m_thermo.size() > 0) {
-        m_start.push_back(m_start.back()
-                          + m_thermo.back()->nSpecies());
-    } else {
-        // otherwise start at 0
-        m_start.push_back(0);
-    }
-
     // the phase with lowest dimensionality is assumed to be the
     // phase/interface at which reactions take place
     if (thermo.nDim() <= m_mindim) {
@@ -513,7 +503,19 @@ void Kinetics::addPhase(thermo_t& thermo)
     }
     m_thermo.push_back(&thermo);
     m_phaseindex[m_thermo.back()->id()] = nPhases();
-    m_kk += thermo.nSpecies();
+    resizeSpecies();
+}
+
+void Kinetics::resizeSpecies()
+{
+    m_kk = 0;
+    m_start.resize(nPhases());
+
+    for (size_t i = 0; i < m_thermo.size(); i++) {
+        m_start[i] = m_kk; // global index of first species of phase i
+        m_kk += m_thermo[i]->nSpecies();
+    }
+    invalidateCache();
 }
 
 void Kinetics::finalize()
@@ -525,6 +527,10 @@ void Kinetics::finalize()
 bool Kinetics::addReaction(shared_ptr<Reaction> r)
 {
     r->validate();
+    if (m_kk == 0) {
+        init();
+    }
+    resizeSpecies();
 
     // If reaction orders are specified, then this reaction does not follow
     // mass-action kinetics, and is not an elementary reaction. So check that it

test/kinetics/kineticsFromScratch.cpp

@@ -19,7 +19,6 @@ public:
         th.push_back(&p_ref);
         importKinetics(p_ref.xml(), th, &kin_ref);
         kin.addPhase(p);
-        kin.init();
     }
 
     IdealGasPhase p;
@@ -318,7 +317,6 @@ public:
         importKinetics(surf_ref.xml(), th, &kin_ref);
         kin.addPhase(surf);
         kin.addPhase(gas);
-        kin.init();
     }
 
     IdealGasPhase gas;
@@ -380,3 +378,109 @@ TEST_F(InterfaceKineticsFromScratch, add_sticking_reaction)
     kin.addReaction(R);
     check_rates(0);
 }
+
+class KineticsAddSpecies : public testing::Test
+{
+public:
+    KineticsAddSpecies()
+        : p_ref("../data/kineticsfromscratch.cti")
+    {
+        std::vector<ThermoPhase*> th;
+        th.push_back(&p_ref);
+        importKinetics(p_ref.xml(), th, &kin_ref);
+
+        p.addUndefinedElements();
+        kin.addPhase(p);
+
+        std::vector<shared_ptr<Species>> S = getSpecies(*get_XML_File("h2o2.cti"));
+        for (auto sp : S) {
+            species[sp->name] = sp;
+        }
+        reactions = getReactions(*get_XML_File("../data/kineticsfromscratch.cti"));
+    }
+
+    IdealGasPhase p;
+    IdealGasPhase p_ref;
+    GasKinetics kin;
+    GasKinetics kin_ref;
+    std::vector<shared_ptr<Reaction>> reactions;
+    std::map<std::string, shared_ptr<Species>> species;
+
+    void check_rates(size_t N, const std::string& X) {
+        for (size_t i = 0; i < kin_ref.nReactions(); i++) {
+            if (i >= N) {
+                kin_ref.setMultiplier(i, 0);
+            } else {
+                kin_ref.setMultiplier(i, 1);
+            }
+        }
+        p.setState_TPX(1200, 5*OneAtm, X);
+        p_ref.setState_TPX(1200, 5*OneAtm, X);
+
+        vector_fp k(kin.nReactions()), k_ref(kin_ref.nReactions());
+        vector_fp w(kin.nTotalSpecies()), w_ref(kin_ref.nTotalSpecies());
+
+        kin.getCreationRates(w.data());
+        kin_ref.getCreationRates(w_ref.data());
+        for (size_t i = 0; i < kin.nTotalSpecies(); i++) {
+            size_t iref = p_ref.speciesIndex(p.speciesName(i));
+            EXPECT_DOUBLE_EQ(w_ref[iref], w[i]) << "sp = " << p.speciesName(i) << "; N = " << N;
+        }
+
+        kin.getFwdRateConstants(k.data());
+        kin_ref.getFwdRateConstants(k_ref.data());
+        for (size_t i = 0; i < kin.nReactions(); i++) {
+            EXPECT_DOUBLE_EQ(k_ref[i], k[i]) << "i = " << i << "; N = " << N;
+        }
+
+        kin.getRevRateConstants(k.data());
+        kin_ref.getRevRateConstants(k_ref.data());
+        for (size_t i = 0; i < kin.nReactions(); i++) {
+            EXPECT_DOUBLE_EQ(k_ref[i], k[i]) << "i = " << i << "; N = " << N;
+        }
+    }
+};
+
+TEST_F(KineticsAddSpecies, add_species_sequential)
+{
+    ASSERT_EQ((size_t) 0, kin.nReactions());
+
+    for (auto s : {"AR", "O", "H2", "H", "OH"}) {
+        p.addSpecies(species[s]);
+    }
+    kin.addReaction(reactions[0]);
+    ASSERT_EQ(5, (int) kin.nTotalSpecies());
+    check_rates(1, "O:0.001, H2:0.1, H:0.005, OH:0.02, AR:0.88");
+
+    p.addSpecies(species["O2"]);
+    p.addSpecies(species["H2O"]);
+    kin.addReaction(reactions[1]);
+    ASSERT_EQ(7, (int) kin.nTotalSpecies());
+    ASSERT_EQ(2, (int) kin.nReactions());
+    check_rates(2, "O:0.001, H2:0.1, H:0.005, OH:0.02, O2:0.5, AR:0.38");
+
+    p.addSpecies(species["H2O2"]);
+    kin.addReaction(reactions[2]);
+    kin.addReaction(reactions[3]);
+    check_rates(4, "O:0.001, H2:0.1, H:0.005, OH:0.02, O2:0.5, AR:0.38"); // no change
+    check_rates(4, "O:0.001, H2:0.1, H:0.005, OH:0.02, O2:0.5, AR:0.35, H2O2:0.03");
+
+    p.addSpecies(species["HO2"]);
+    kin.addReaction(reactions[4]);
+    check_rates(5, "O:0.01, H2:0.1, H:0.02, OH:0.03, O2:0.4, AR:0.3, H2O2:0.03, HO2:0.01");
+}
+
+TEST_F(KineticsAddSpecies, add_species_err_first)
+{
+    for (auto s : {"AR", "O", "H2", "H"}) {
+        p.addSpecies(species[s]);
+    }
+    ASSERT_THROW(kin.addReaction(reactions[0]), CanteraError);
+    ASSERT_EQ((size_t) 0, kin.nReactions());
+
+    p.addSpecies(species["OH"]);
+    kin.addReaction(reactions[0]);
+    ASSERT_EQ(5, (int) kin.nTotalSpecies());
+    ASSERT_EQ((size_t) 1, kin.nReactions());
+    check_rates(1, "O:0.001, H2:0.1, H:0.005, OH:0.02, AR:0.88");
+}