diff --git a/src/kinetics/AqueousKinetics.cpp b/src/kinetics/AqueousKinetics.cpp index 5301365cc..91ea9a39d 100644 --- a/src/kinetics/AqueousKinetics.cpp +++ b/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); + } } } diff --git a/src/kinetics/GasKinetics.cpp b/src/kinetics/GasKinetics.cpp index 186703cb7..fda70b1ea 100644 --- a/src/kinetics/GasKinetics.cpp +++ b/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); + } } } //==================================================================================================================== diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index ecdb411c0..e8cda8afa 100644 --- a/src/kinetics/InterfaceKinetics.cpp +++ b/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(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(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(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; } //================================================================================================ diff --git a/src/kinetics/importKinetics.cpp b/src/kinetics/importKinetics.cpp index 2515e0438..261d1c525 100644 --- a/src/kinetics/importKinetics.cpp +++ b/src/kinetics/importKinetics.cpp @@ -902,6 +902,7 @@ bool installReactionArrays(const XML_Node& p, Kinetics& kin, p.getChildren("reactionArray",rarrays); int na = static_cast(rarrays.size()); if (na == 0) { + kin.finalize(); return false; } for (int n = 0; n < na; n++) {