# 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 `Reactor`s, 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
std::cout << gas->thermo()->report(); //print the intitial state of the mixture to the console
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 `Integrator`s 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:

``````newIntegrator("CVODE")
``````

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)

``````m_integ(newIntegrator("CVODE"))
``````

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)

``````sim.advance(tnow);
``````

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) {
initialize();
} else if (!m_integrator_init) {
reinitialize();
``````

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)

``````m_integ->integrate(time);
``````

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 (d`y`/d`t`) 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 `Reactor`s 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.

@paulblum