GSoC 2020: 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 initial 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 anIntegratorinstance of typeitype.
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 forCVODES. (Documentation)
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 theCVODESpackage. Integrates over a time interval defined by the user, by callingcvStep()to do internal time steps. (Documentation: seeCVODESUser 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_memis a pointer to the block of memory that was allocated and configured during initialization.After execution,
m_ywill contain the computed solution vector, and will later be used to update theReactorNetto its time-integrated state .The
CV_NORMALoption tells the solver that it should continue taking internal timesteps until it has reached user-specifiedtout(or just passed it, in which case solutions are reached by interpolation). This provides the appropriate functionality forReactorNet::advance(). The alternate option,CV_ONE_STEP, tells the solver to take a single internal step, and is used inReactorNet::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)
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
updateState(y);
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.