`m_rxn_rates` type incorrect in template class MultiRate

[chengcli]

Problem description

In the template class MultiRate, its member m_rxn_rates has the following type:

vector<pair<size_t, RateType>> m_rxn_rates;

It causes an inconsistency between the ReactionRate object stored in Reaction and the ReactionRate object stored in MultiRate. Unless it is intentionally designed to be the case, the type RateType should be changed to RateType&.

Behavior

ReactionRate objects are stored in two places, one in the class Reaction and the other in the template class MultiRate. They should reflect the same object rather than different objects. In the current version, the ReactionRate object is first instantiated when a Reaction object is instantiated. If the Reaction object is instantiated by a derived class of Kinetics that has a MultiRate evaluator, such as the BulkKinetics class, the same ReactionRate object is passed to a MultiRate object. For example, the following code snippet in BulkKinetics.cpp shows the passing of the ReactionRate object to a MultiRate evaluator:

    // Add reaction rate to evaluator
    size_t index = m_bulk_types[rtype];
    m_bulk_rates[index]->add(nReactions() - 1, *rate);

The function signature and behavior of the add function in MultiRate is:

    void add(size_t rxn_index, ReactionRate& rate) override {
        m_indices[rxn_index] = m_rxn_rates.size();
        m_rxn_rates.emplace_back(rxn_index, dynamic_cast<RateType&>(rate));
        m_shared.invalidateCache();
    }

So, it accepts a reference to ReactionRate and casts it the correct derived class of ReactionRate. However, since the type of the second member of m_rxn_rates is RateType rather than RateType&, the copy constructor is called in the emplace_back call. As a result, the ReactionRate object stored in m_rxn_rates is a copy of the ReactionRate object stored in class Reaction. This causes an inconsistency between the ReactionRate object in MultiRate and the object in Reaction.

How to fix

Simply changing the type of RateType to RateType& in the declaration of m_rxn_rates will fix the problem. Changing the concrete type to a reference type will avoid copying the object. So the two ReactionRate objects in Reaction and in MultiRate are exactly the same. All tests shall pass without any issues.

[speth]

This isn't a typo -- we're aware that MultiRate makes copies of the ReactionRate objects. We've discussed this previously in #1211, but without coming to a resolution. This is, essentially, the reason that we require the use of the Kinetics.modifyReaction function to handle dynamic rate modifications.

I think the change suggested here is insufficient. For instance, it is possible to replace the Reaction.m_rate object outright, which would leave the reference in MultiRate.m_rxn_rates pointing to now-freed memory rather than the new ReactionRate object.

[luminoctum]

Thanks for pointing out #1211 and speth@b8f6132. It could be a potential memory problem if Reaction.m_rate is changed. However, upon inspecting the code, it seems to me that the only way a user can update Reaction.m_rate is through creating a new Reaction object. If a user can only change the Reaction object inside Kinetics via the modifyReaction method, the method itself can ensure that MutiRate.m_rxn_rates are updated correctly as has been done in BulkKinetics:

void BulkKinetics::modifyReaction(size_t i, shared_ptr<Reaction> rNew)
{
    Kinetics::modifyReaction(i, rNew);
    ...
    size_t index = m_bulk_types[rtype];
    rate->setRateIndex(i);
    rate->setContext(*rNew, *this);
    m_bulk_rates[index]->replace(i, *rate);
}

So I think the memory is fine if RateType is changed to RateType&.

Context

I came across this problem when I was augmenting Cantera for photochemical reactions. I'm trying to combine all photolysis yields of a species into a "net photolysis reaction" for performance. Since the quantum yields depend on the wavelength of the photon and the temperature, the stoichiometric coefficients of the "net reaction" are variable. MutiRate objects are capable of calculating the net stoichiometric coefficients after getting the reaction rate. However, when passing the net stoichiometric coefficients calculated by the MultiRate object back to the Kinetics object, I encountered the inconsistency between the rate in Reaction and the rate in MultiRate. I tried to access the ReactionRate in Reaction but found that it was not updated because it was the MultiRate object that carried out the calculation. The trick of changing RateType to RateType& works for my purpose, but I'm unsure about other side effects.

[ischoegl]

Within the context of #1211, the original motivation for using copies were performance considerations. Keeping things synchronized is a little bit of a pain (an understatement), but as @speth mentioned the current layout is indeed intentional and thus not a bug or issue. While I am closing this issue, starting a discussion in Cantera/enhancements about alternatives (with associated benchmarks) would certainly be appropriate.