[Kinetics] Fixed crashes when reaction mechanism contains no reactions

Make sure several arrays are never of length zero so that taking &v[0]
is always well-defined.

src/kinetics/AqueousKinetics.cpp

@@ -594,6 +594,16 @@ void AqueousKinetics::finalize()
 {
     if (!m_finalized) {
         m_finalized = true;
+
+        // Guarantee that these arrays can be converted to double* even in the
+        // special case where there are no reactions defined.
+        if (!m_ii) {
+            m_perturb.resize(1, 1.0);
+            m_ropf.resize(1, 0.0);
+            m_ropr.resize(1, 0.0);
+            m_ropnet.resize(1, 0.0);
+            m_rkcn.resize(1, 0.0);
+        }
     }
 }
 

src/kinetics/GasKinetics.cpp

@@ -838,6 +838,16 @@ void GasKinetics::finalize()
         concm_3b_values.resize(m_3b_concm.workSize());
         concm_falloff_values.resize(m_falloff_concm.workSize());
         m_finalized = true;
+
+        // Guarantee that these arrays can be converted to double* even in the
+        // special case where there are no reactions defined.
+        if (!m_ii) {
+            m_perturb.resize(1, 1.0);
+            m_ropf.resize(1, 0.0);
+            m_ropr.resize(1, 0.0);
+            m_ropnet.resize(1, 0.0);
+            m_rkcn.resize(1, 0.0);
+        }
     }
 }
 //====================================================================================================================

src/kinetics/InterfaceKinetics.cpp

@@ -356,10 +356,7 @@ void InterfaceKinetics::updateKc()
 void InterfaceKinetics::checkPartialEquil()
 {
     vector_fp dmu(nTotalSpecies(), 0.0);
-    vector_fp rmu(nReactions(), 0.0);
-    vector_fp frop(nReactions(), 0.0);
-    vector_fp rrop(nReactions(), 0.0);
-    vector_fp netrop(nReactions(), 0.0);
+    vector_fp rmu(std::max<size_t>(nReactions(), 1), 0.0);
     if (m_nrev > 0) {
         doublereal rt = GasConstant*thermo(0).temperature();
         cout << "T = " << thermo(0).temperature() << " " << rt << endl;
@@ -380,16 +377,14 @@ void InterfaceKinetics::checkPartialEquil()
 
         // compute Delta mu^ for all reversible reactions
         m_rxnstoich.getRevReactionDelta(m_ii, DATA_PTR(dmu), DATA_PTR(rmu));
-        getFwdRatesOfProgress(DATA_PTR(frop));
-        getRevRatesOfProgress(DATA_PTR(rrop));
-        getNetRatesOfProgress(DATA_PTR(netrop));
+        updateROP();
         for (size_t i = 0; i < m_nrev; i++) {
             size_t irxn = m_revindex[i];
             cout << "Reaction " << reactionString(irxn)
                  << "  " << rmu[irxn]/rt << endl;
             printf("%12.6e  %12.6e  %12.6e  %12.6e \n",
-                   frop[irxn], rrop[irxn], netrop[irxn],
-                   netrop[irxn]/(frop[irxn] + rrop[irxn]));
+                   m_ropf[irxn], m_ropr[irxn], m_ropnet[irxn],
+                   m_ropnet[irxn]/(m_ropf[irxn] + m_ropr[irxn]));
         }
     }
 }
@@ -1273,7 +1268,8 @@ void InterfaceKinetics::init()
 void InterfaceKinetics::finalize()
 {
     Kinetics::finalize();
-    m_rwork.resize(nReactions());
+    size_t safe_reaction_size = std::max<size_t>(nReactions(), 1);
+    m_rwork.resize(safe_reaction_size);
     size_t ks = reactionPhaseIndex();
     if (ks == npos) throw CanteraError("InterfaceKinetics::finalize",
                                            "no surface phase is present.");
@@ -1284,13 +1280,23 @@ void InterfaceKinetics::finalize()
                            +int2str(m_surf->nDim()));
 
     m_StandardConc.resize(m_kk, 0.0);
-    m_deltaG0.resize(m_ii, 0.0);
-    m_ProdStanConcReac.resize(m_ii, 0.0);
+    m_deltaG0.resize(safe_reaction_size, 0.0);
+    m_ProdStanConcReac.resize(safe_reaction_size, 0.0);
 
     if (m_thermo.size() != m_phaseExists.size()) {
         throw CanteraError("InterfaceKinetics::finalize", "internal error");
     }
 
+    // Guarantee that these arrays can be converted to double* even in the
+    // special case where there are no reactions defined.
+    if (!m_ii) {
+        m_perturb.resize(1, 1.0);
+        m_ropf.resize(1, 0.0);
+        m_ropr.resize(1, 0.0);
+        m_ropnet.resize(1, 0.0);
+        m_rkcn.resize(1, 0.0);
+    }
+
     m_finalized = true;
 }
 
@@ -1430,7 +1436,7 @@ void InterfaceKinetics::setPhaseStability(const size_t iphase, const int isStabl
 //================================================================================================
 void EdgeKinetics::finalize()
 {
-    m_rwork.resize(nReactions());
+    m_rwork.resize(std::max<size_t>(nReactions(), 1));
     size_t ks = reactionPhaseIndex();
     if (ks == npos) throw CanteraError("EdgeKinetics::finalize",
                                            "no edge phase is present.");
@@ -1439,6 +1445,17 @@ void EdgeKinetics::finalize()
         throw CanteraError("EdgeKinetics::finalize",
                            "expected interface dimension = 1, but got dimension = "
                            +int2str(m_surf->nDim()));
+
+    // Guarantee that these arrays can be converted to double* even in the
+    // special case where there are no reactions defined.
+    if (!m_ii) {
+        m_perturb.resize(1, 1.0);
+        m_ropf.resize(1, 0.0);
+        m_ropr.resize(1, 0.0);
+        m_ropnet.resize(1, 0.0);
+        m_rkcn.resize(1, 0.0);
+    }
+
     m_finalized = true;
 }
 //================================================================================================

src/kinetics/importKinetics.cpp

@@ -902,6 +902,7 @@ bool installReactionArrays(const XML_Node& p, Kinetics& kin,
     p.getChildren("reactionArray",rarrays);
     int na = static_cast<int>(rarrays.size());
     if (na == 0) {
+        kin.finalize();
         return false;
     }
     for (int n = 0; n < na; n++) {