Enable development of models written in high-level languages (e.g. Python)

[speth]

Abstract

The ultimate goal of this project is to allow users to write their own models for thermodynamics, kinetics, reactors, etc. in easy-to-use high level languages like Python, rather than requiring models to be implemented in C++.

For the sake of simplicity, this proposal is written with reference to models created in Python, although it could apply to other languages like Julia or anything else that exposes a C API to be able to call code written in that language.

Motivation

Currently, all of the models present in Cantera, representing everything from species thermodynamic properties to reaction rates to reactor governing equations, are implemented in C++. While this is computationally efficient, creating new models in C++ is a significant barrier to scientific users who are not necessarily expert software engineers, as well as being time-consuming even for those who are. This barrier discourages exploratory model implementations, which are valuable in many cases where those new models need to be evaluated and refined within the context of other Cantera capabilities, e.g., a new thermodynamic model that must be evaluated within the context of kinetic and transport calculations. Allowing new models to be implemented in high level languages like Python will reduce this barrier, providing scientific users with an environment that is well-suited for rapid prototyping and making it easy to share their developments as standalone modules. For a model that proves to be useful and needs to be used in performance-sensitive applications, a Python implementation could later be used as a “reference implementation” for the C++ version of the model.

Description

Cantera's interfaces to other languages mostly involve code in other languages calling the Cantera library through its C or C++ APIs. Implementing the capability described here requires the reverse - having Cantera call code written in another language which presents a C or C++ API. The Cantera Python interface already does this in one particular instance, for functions that provide reactor inputs as a function of time, where users can provide this information as native Python functions. While the implementation under the hood is not trivial (see func1.pyx and funcWrapper.h), in the case of Python, many of the difficulties such as translating C++ exceptions into Python has already been worked out, and the additions required are mostly just to handle functions with different forms besides double = f(double).

The next step is then to allow the user to define a Python class that overrides some virtual member functions of an existing Cantera C++ class, such as ThermoPhase. For this feature to be most useful, users should be able to start with any of the existing Cantera phase models. The implementation for this that I think would work is introduce a new templated class which derives from the base ThermoPhase type specified by the user, which holds a set of objects encapsulating the user-defined functions. A skeleton of this class might look something like the following:

template <class T>
class UserDefinedPhase : public T {
public:
    // An override of a ThermoPhase function
    virtual void getPartialMolarEnthalpies(double* hbar) const {
        if (m_getPartialMolarEnthalpiesFunc) {
            m_getPartialMolarEnthalpiesFunc(hbar);
        } else {
            T::getPartialMolarEnthalpies(hbar);
        }
    }

    // Setup function; this would need overrides for each distinct function signature handled
    void setOverride(const std::string& name, std::function<void(double*) func) {
        if (name == "getPartialMolarEnthalpies") {
            m_partialMolarEnthalpiesFunc = func;
        }
    }

private:
    std::function<void(double*)> m_getPartialMolarEnthalpiesFunc;
}

template class UserDefinedPhase<IdealGasPhase>;

This is of course a fairly "simple" example, and there are certainly some complexities to work out, like how to allow the user defined functions invoke the method of the base class.

On the Python side, I think the user would create a class that inherited from a Python class UserDefinedPhase. The UserDefinedPhase.__init__ method would then be responsible for encapsulating the methods in the user's class and calling the C++ setOverride method.

Some additional machinery that I haven't thought too deeply about yet would be required to allow Cantera to load the user's Python module based on some information defined in the input YAML file.

I think the approach outlined here would work well for the ThermoPhase, SpeciesThermoInterpType, Kinetics, Reactor, and Domain1D classes. One important exception to this at the moment is Reaction classes (or really, the reaction rate classes) which don't follow the same pattern of using (essentially) abstract base classes in a generic manner, which makes it difficult to add new reaction types even from C++. I think an overhaul there (per #63) is in order before this approach could be applied to reactions.

Alternatives

I don't know of an alternative that would provide the same flexibility. Perhaps something based on code generation would be an option, although I think that would require unique solutions for each supported language, and I don't think it would be able to provide the same coupling with Cantera's existing models. I suppose "rewrite Cantera from scratch in your own favorite language", but that's quite the time investment, and in the case of Python would come at a huge performance penalty as well.

References

• Enhancement proposal Make reactions extendable #63 - Make reactions extensible at the C++ level
• Posts to Users' Group that could be addressed by this capability:
  • Time dependent conditions/reservoir properties in Reactor Network
  • Simulating isentropic compression according to a known pressure trace
  • Simplified transport property relations for virtual species
• Other enhancements that could be implemented using this feature:
  • Tsang falloff format #113

[ischoegl]

@speth ... I am really happy to see this writeup also (I'm again tagging @jiweiqi here, as this is to some extent a spillover from a discussion on the user group).

[ischoegl]

@speth ... I decided to reopen Cantera/cantera#745 as Cantera/cantera#982 as I believe it would provide an initial 'hook' for custom reactions using the factory approach used elsewhere (it implements a ReactionFactory).

Regarding custom reactions (e.g. in Python), this may actually be the simplest case, as it mostly requires custom rates in RxnRates.h, similar to:

class PythonRate
{
public:
    PythonRate();

    // set custom rate
    void setPythonRate(Func1* f);

    // Update the value the natural logarithm of the rate constant.
    double updateLog(double logT, double recipT) const {
        return  std::log( m_ratefunc(T, recipT) );
    }

    // Update the value the rate constant.
    double updateRC(double logT, double recipT) const {
        return m_ratefunc(T, recipT);
    }

    //! Return NaN - pre-exponential factor does not apply
    double preExponentialFactor() const {
        return std::numeric_limits<double>::quiet_NaN();
    }

    // [...] other standard functions not shown (mostly ignoring A, b, E) 

protected:
    Func1* m_ratefunc;
};

Once Cantera/cantera#982 is adopted, adding a custom PythonReaction to Reaction.h plus a new entry in the ReactionFactory is trivial. Likewise, adding a user-defined rate function from Python should be straightforward. A similar approach may be possible for other high-level languages.

PS: Looked into this a little more and there are some gnarly bits in GasKinetics.h.

[ischoegl]

Adding links to embedding of Julia and how to call Julia from C via cxxWrap.jl here.

[ischoegl]

For custom Reactions in Python, I implemented a draft for a new CustomReaction class in Cantera/cantera#982, which can be instantiated as

rxn = CustomReaction(equation='H2 + O <=> H + OH',
                     rate=lambda T: 38.7 * T**2.7 * exp(-3150.15428/T),
                     kinetics=gas)

(which taps into Func1). As it is not formalized, there is no equivalent YAML syntax.

[ischoegl]

Cantera/cantera#982 (now merged) implements custom reaction rates specified by Func1 objects.