Step 1: Reactor(s) created in Reactor Network¶

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 in Python 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.

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

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.

Optional Preconditioning¶

Some mole reactors are capable of leveraging preconditioning to accelerate integration of the system. The former code can be modified as follows to use preconditioning. The preconditioner can also be assigned to a python object and tunable parameters can adjusted.

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.IdealGasMoleReactor(gas)  # create a reactor containing the gas
sim = ct.ReactorNet([reac])  # add the reactor to a new ReactorNet simulator
sim.preconditioner = ct.AdaptivePreconditioner() # add preconditioner to the network
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

Step 2: advance() method called¶

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. The initial condition information is passed off to the Integrator when calling advance(). 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.

Step 3: Information about current gas state provided to an Integrator ¶

The Integrator class is Cantera's interface for ODE/DAE system integrators.

Integrator is a polymorphic base class; it defines a set of virtual methods that derived classes (the actual system integrators) will provide implementations for.

The newIntegrator() factory method creates and returns an Integrator object of the specified type. Calling newIntegrator("CVODE") creates a new CVodesIntegrator object for integrating an ODE system, while calling newIntegrator("IDA") creates a new IdasIntegrator object for integrating a DAE system. The appropriate integrator type is determined by the ReactorNet class based on the types of the installed reactors. Below, the implementation of CvodesIntegrator is described; the IdasIntegrator works in a similar way, though the function names are different.

Step 4: Communicate with SUNDIALS using a wrapper function¶

Because SUNDIALS is written in C, the CVodesIntegrator C++ wrapper is used to access the solver. The CVodesIntegrator class is a C++ wrapper class for CVODES. (Documentation) The CVodesIntegrator class makes the appropriate call to the CVODES driver function, CVode().

Step 5: Cvode() driver function is called¶

Method CVode() is the main driver of the CVODES package. CVode() integrates over a time interval defined by the user, by calling cvStep() to do internal time steps (not specified by the user). (Documentation: see CVODES User Guide)

The arguments taken by the CVode() method is shown below:

int CVode(void *cvode_mem, realtype tout, N_Vector yout, realtype *tret, int itask);

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().

The result of the CVode() method is assigned to the flag object. CVode() returns 1 or 0, corresponding to a successful or unsuccessful integration, respectively.

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

Step 6: FuncEval class describes ODEs to solve¶

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.

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 (Documentation). Classes derived from FuncEval will implement the evaluation of the provided ODE system.

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.

Let's take a look at how ReactorNet implements this FuncEval object. ReactorNet actually points to itself when defining a FuncEval type, meaning it defines itself as a FuncEval derivative.

Then, ReactorNet initializes the Integrator, using a reference to itself (as a FuncEval) from the this pointer.

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.

To evaluate the reactor governing equations the following parameters must be known:

1. t: Current time in seconds.

2. LHS: pointer to start of vector of left-hand side coefficients for governing equations.

   Has length m_nv, default values 1.

3. RHS: pointer to start of vector of right-hand side coefficients for governing equations.

   Has length m_nv, default values 0.

virtual void eval(double t, double* LHS, double* RHS);

eval() is called by ReactorNet::eval.

The above code shows the necessary inputs for solving the ODEs using the eval() function. eval() takes in the value of each state variable derivative (ydot) at a time t, and will write the integrated values for each state variable to the solution vector (y).

Step 7: eval() is called to solve provided ODEs¶

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

void ReactorNet::eval(double t, double* y, double* ydot, double* p)
{
    m_time = t;
    updateState(y);
    m_LHS.assign(m_nv, 1);
    m_RHS.assign(m_nv, 0);
    for (size_t n = 0; n < m_reactors.size(); n++) {
        m_reactors[n]->applySensitivity(p);
        m_reactors[n]->eval(t, m_LHS.data() + m_start[n], m_RHS.data() + m_start[n]);
        size_t yEnd = 0;
        if (n == m_reactors.size() - 1) {
            yEnd = m_RHS.size();
        } else {
            yEnd = m_start[n + 1];
        }
        for (size_t i = m_start[n]; i < yEnd; i++) {
            ydot[i] = m_RHS[i] / m_LHS[i];
        }
        m_reactors[n]->resetSensitivity(p);
    }
    checkFinite("ydot", ydot, m_nv);
}

ReactorNet's eval() method evaluates 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.

This documentation is based off @paulblum's blog post.