GSoC 2020 Blog Post 3

How Does Cantera's Reactor Network Time Integration Feature Work?

There's a great description of the science behind Cantera's reactor network simulation capabilities available on the Cantera website, here. This post will go into more developer-oriented detail about how the last step, ReactorNet's time integration methods, actually work. A ReactorNet object doesn't perform time integration on its own. It generates a system of ODE's based on the combined governing equations of all contained Reactors, which is then passed off to an Integrator object for solution. What is an Integrator? How does this work?

Reactor Network Time Integration, Explained.

First, let's take a look at a basic example to see how we might utilize Cantera's time integration functionality. We'll simulate an isolated reactor that is homogeneously filled by a gas mixture (the gas state used in this example is arbitrary, but interesting because it's explosive). Then we'll advance the simulation in time to an (arbitrary) absolute time of 1 second, noting the changes in the state of the gas. Follow along by typing this code into an interactive Python interpreter (like IPython):

>>> import cantera as ct                           #import Cantera's Python module
>>> gas = ct.Solution('gri30.yaml')                #create a default GRI-Mech 3.0 gas mixture
>>> gas.TPX = 1000.0, ct.one_atm, 'H2:2,O2:1,N2:4' #set gas to an interesting state
>>> reac = ct.IdealGasReactor(gas)                 #create a reactor containing the gas
>>> sim = ct.ReactorNet([reac])                    #add the reactor to a new ReactorNet simulator
>>> gas()                #view the initial state of the mixture (state summary will be printed to console)
>>> sim.advance(1)       #advance the simulation to the specified absolute time, t = 1 sec
>>> gas()                #view the updated state of the mixture, reflecting properties at t = 1 sec

Equivalently, the following can be compiled and run using Cantera's C++ interface:

#include "cantera/zerodim.h" //include Cantera's 0D reactor simulation module
using namespace Cantera;     //activate Cantera namespace to identify scope of class and method references
int main() {
    auto gas = newSolution("gri30.yaml"); //create a default GRI-Mech 3.0 gas mixture
    gas->thermo()->setState_TPX(1000.0, OneAtm, "H2:2,O2:1,N2:4"); //set gas to an interesting state
    Reactor reac;                         //create an empty Reactor
    reac.insert(gas);                     //fill the reactor with the specified gas
    ReactorNet sim;                       //create an empty ReactorNet simulator
    sim.addReactor(reac);                 //add the reactor to the ReactorNet
    std::cout << gas->thermo()->report(); //print the intitial state of the mixture to the console
    sim.advance(1);                       //advance the simulation to absolute time t = 1 sec
    std::cout << gas->thermo()->report(); //print the updated state of the mixture to the console
    return 0;

For a more advanced example that adds inlets and outlets to the reactor, see Cantera's combustor example (Python | C++). Additional examples can be found in the Python Reactor Network Examples section of the Cantera website.

In any case, after properly configuring a reactor network and its components in Cantera, a call to the ReactorNet's advance() method can be used to predict the state of the network at a specified time. Transient physical and chemical interactions are simulated by integrating the network's system of ODE governing equations through time, a process that's actually performed by an external Integrator object. What is an Integrator?

The Integrator class is Cantera's interface for ODE system integrators. This general-purpose ODE system integration tool can be accessed in any Cantera project by including the Integrator.h header file in your code:

  • Class Integrator: A base class interface for ODE system integrators. (Documentation)
    • Declared and (virtually) implemented in Integrator.h, line 52 (see this on GitHub)

Integrator is a polymorphic base class; it defines a set of virtual functionalities that derived classes (the actual ODE system integrators) will provide implementations for. This is done so that different types of Integrators can be used interchangeably, without having to modify their references in code. Method implementations in different Integrator subclasses can be executed using the same call to the Integrator base class's virtual function—the base class will refer the call to the appropriate subclass implementation, based on the Integrator object's type. How do you set the type of an Integrator?

Conveniently, Integrator.h provides newIntegrator(), a factory method for creating Integrator instances of arbitrary type. This feature hides subclass implementation modules from users, who only know of the generic Integrator object that they received from the factory method:

  • Factory Method newIntegrator(): Creates and returns a pointer to an Integrator instance of type itype.
    • Integrator* newIntegrator(const std::string& itype);
    • Declared in Integrator.h, line 237 (see this on GitHub)
    • Implemented in ODE_integrators.cpp, line 13 (see this on GitHub)

The header files of different Integrator implementations are included near the top of ODE_integrators.cpp (see this on GitHub). Cantera only includes one type of Integrator by default, the CVODES ODE solver, which automatically installs alongside Cantera. CVODES is an externally developed and maintained solver, a part of @LLNL's SUNDIALS library. If you're interested in learning more specific details about how CVODES actually solves an ODE system, I recommend that you read through the CVODES User Guide for detailed documentation and explanation of the module. Note that ReactorNet configures CVODES to solve via Backward Differentiation Formulas (see CVODES User Guide, section 2.1), based on linear system solutions provided by the SUNDIALS SUNLinSol_LapackDense module (see CVODES User Guide, section 10.7).

Because CVODES is written in C, the CVodesIntegrator C++ wrapper is used to access the solver:

  • Class CVodesIntegrator: A C++ wrapper class for CVODES. (Documentation)
    • Declared in CVodesIntegrator.h, line 25 (see this on GitHub)
    • Implemented in CVodesIntegrator.cpp, line 79 (see this on GitHub)
    • Included in ODE_integrators.cpp, line 8 (see this on GitHub)

To create an instance of this CVODES-type Integrator, the factory method can be called with the "CVODE" keyword:


This call returns a generic Integrator pointer, whose virtual functions are overridden by those of a CVodesIntegrator. This is exactly how a ReactorNet creates its CVODES-type Integrator object, before storing it locally as m_integ for future reference:

ReactorNet.cpp, line 18 (see this on GitHub)


So, what actually happens when you call one of a ReactorNet's time integration functions? Let's follow the execution of a call to ReactorNet::advance() through the source code. Consider this call to 'sim', the ReactorNet object in Cantera's C++ combustor example:

combustor.cpp, line 122 (see this on GitHub)


The first thing that ReactorNet::advance() does is check for proper initialization, and initialize if needed:

ReactorNet.cpp, line 128 (see this on GitHub)

void ReactorNet::advance(doublereal time)
    if (!m_init) {
    } else if (!m_integrator_init) {

A ReactorNet always needs to be initialized before solving a new reactor network configuration, or after making any changes to the Integrator's settings. Initialization can be done with a call to ReactorNet::initialize(), which will allocate new memory, configure Integrator settings, initialize all substructures, and populate internal memory appropriately with required data and specifications about the current system.

After a certain ReactorNet configuration has been initialized, it can also be reinitialized to simulate the same network with modified initial conditions. This can be done by calling ReactorNet::reinitialize(), which updates the internal memory with the data and specifications of the modified system. ReactorNet::reinitialize() will be called automatically after changing the simulation's initial time with ReactorNet::setInitialTime(), or after modifying the contents of any contained reactor.

Once the ReactorNet has been properly initialized and its internal memory is up-to-date with the current network state, the problem is passed off to m_integ, the Integrator object, for solution:

ReactorNet.cpp, line 135 (see this on GitHub)


The CVodesIntegrator wrapper class will then make the appropriate call to the CVODES driver function, CVode():

  • Method CVode(): Main driver of the CVODES package. Integrates over a time interval defined by the user, by calling cvStep() to do internal time steps. (Documentation: see CVODES User Guide, section 4.5.6)
    • int CVode(void *cvode_mem, realtype tout, N_Vector yout, realtype *tret, int itask);
    • Implemented in cvodes.c, line 2771 (see this on GitHub)

CVodesIntegrator.cpp, line 458 (see this on GitHub)

int flag = CVode(m_cvode_mem, tout, m_y, &m_time, CV_NORMAL);

There are some interesting things to note about this call to CVode():

  • m_cvode_mem is a pointer to the block of memory that was allocated and configured during initialization.
  • After execution, m_y will contain the computed solution vector, and will later be used to update the ReactorNet to its time-integrated state .
  • The CV_NORMAL option tells the solver that it should continue taking internal timesteps until it has reached user-specified tout (or just passed it, in which case solutions are reached by interpolation). This provides the appropriate functionality for ReactorNet::advance(). The alternate option, CV_ONE_STEP, tells the solver to take a single internal step, and is used in ReactorNet::step().

How does CVODES know what ODE system it should be solving? The ODE system was actually already specified using CVodeInit(), one of the methods automatically invoked during the ReactorNet::initialize() routine. CVODES requires that its user provide a C function that defines their ODE, able to compute the right-hand side of the ODE system (dy/dt) for a given value of the independent variable, t, and the state vector, y. For more information about ODE right-hand side function requirements, see CVODES User Guide, section 4.6.1.

The CVodesIntegrator wrapper class provides a useful C++ interface for configuring this C function by pairing with FuncEval, an abstract base class for ODE right-hand-side function evaluators. Like the Integrator base class, FuncEval defines virtual functions that derived classes will provide the implementations for. In this case, classes derived from FuncEval will implement the actual evaluation of their particular ODE system.

  • Class FuncEval: An abstract base class for ODE right-hand-side function evaluators. (Documentation)
    • Declared in FuncEval.h, line 26 (see this on GitHub)
    • Implemented in FuncEval.cpp, line 7 (see this on GitHub)

An ODE right-hand-side evaluator is always needed in the ODE solution process (it's the only way to describe the system!), and for that reason a FuncEval object is a required parameter when initializing any type of Integrator:

Integrator.h, line 99 (see this on GitHub)

 * Initialize the integrator for a new problem. Call after all options have
 * been set.
 * @param t0   initial time
 * @param func RHS evaluator object for system of equations.
virtual void initialize(doublereal t0, FuncEval& func) {

Let's take a look at how ReactorNet implements this FuncEval object. Instead of creating an external FuncEval subclass object, it defines itself as a FuncEval derivative:

ReactorNet.h, line 23 (see this on GitHub)

class ReactorNet : public FuncEval

Then, it initializes the Integrator, using a reference to itself (as a FuncEval) from the 'this' pointer:

ReactorNet.cpp, line 112 (see this on GitHub)

m_integ->initialize(m_time, *this);

To be a valid FuncEval object, a ReactorNet needs to provide implementations for all of FuncEval's virtual functions, particularly the actual ODE right-hand-side computation function, FuncEval::eval(). Note that this is declared as a pure virtual function, which makes FuncEval an abstract class:

FuncEval.h, line 32 (see this on GitHub)

 * Evaluate the right-hand-side function. Called by the integrator.
 * @param[in] t time.
 * @param[in] y solution vector, length neq()
 * @param[out] ydot rate of change of solution vector, length neq()
 * @param[in] p sensitivity parameter vector, length nparams()
virtual void eval(double t, double* y, double* ydot, double* p)=0;

Along with the rest of FuncEval's virtual functions, an appropriate override is provided for FuncEval::eval() in ReactorNet:

ReactorNet.cpp, line 233 (see this on GitHub)

void ReactorNet::eval(doublereal t, doublereal* y, doublereal* ydot, doublereal* p)
    m_time = t; // This will be replaced at the end of the timestep
    for (size_t n = 0; n < m_reactors.size(); n++) {
        m_reactors[n]->evalEqs(t, y + m_start[n], ydot + m_start[n], p);
    checkFinite("ydot", ydot, m_nv);

ReactorNet's eval() method invokes calls to Reactor::evalEqs(), to evaluate the governing equations of all Reactors contained in the network. This brings us right back to where we started; for more information see Cantera's reactor network science page.

Hope you enjoyed the post.


Keep Reading:

Next Post - GSoC 2020 Blog Post 4

Previous Post - GSoC 2020 Blog Post 2

Start from the Beginning - Introduction