How Time Integration of Reactor Networks Works#

This section provides a developer-oriented description of how time integration of reactor networks works in Cantera.

A ReactorNet object doesn’t perform time integration on its own. Instead, the role of the ReactorNet object is to assemble a system of governing equations from multiple reactors. It then uses the CVODES or IDAS integrators from SUNDIALS to integrate the governing equations in time or space, respectively.

Creating a Reactor and Reactor Network#

First, let’s take a look at what happens when creating a reactor network by setting up an isolated reactor in Python.

import cantera as ct
gas = ct.Solution("gri30.yaml")
gas.TPX = 1000.0, ct.one_atm, "H2:2,O2:1,N2:4"
reac = ct.IdealGasReactor(gas)
sim = ct.ReactorNet([reac])

The __init__ method of the Python ReactorNet class calls ReactorNet::addReactor() for each Reactor object provided in the list supplied. When the first Reactor is added to the network, the ReactorNet creates a new Integrator used to integrate the governing equations.

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) provide implementations for.

The newIntegrator() factory function 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. A FlowReactor defines a DAE system and uses the IDAS integrator, while the other reactor types define ODE systems and use the CVODES integrator. In this guide, the implementation of CVodesIntegrator is described; IdasIntegrator works in a similar way, although the IDAS function names are different.

Communicating with SUNDIALS using Wrapper Functions#

Because SUNDIALS is written in C, the CVodesIntegrator C++ wrapper is used to access the CVODES solver and make the appropriate calls such as those to the integrator function CVode. For details on CVODES API, see the CVODES User Guide.

Calling the ReactorNet.advance Method#

After configuring a reactor network and its components in Cantera, a call to the ReactorNet.advance() method can be used to obtain 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.

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

Calling the ReactorNet::advance() method invokes the method CVodesIntegrator::integrate() to integrate the system to the specified time. It is most efficient to let CVODES determine the actual integration step sizes on its own. Therefore, we take individual time steps using CVODES using the CVode function until we have reached or passed the specified time. We then interpolate the system state back to the specified time using the CVodeGetDky function. With some additional handling of special cases and errors, the implementation of CVodesIntegrator::integrate() is:

void CVodesIntegrator::integrate(double tout)
{
    if (tout == m_time) {
        return;
    } else if (tout < m_time) {
        throw CanteraError("CVodesIntegrator::integrate",
                           "Cannot integrate backwards in time.\n"
                           "Requested time {} < current time {}",
                           tout, m_time);
    }
    int nsteps = 0;
    while (m_tInteg < tout) {
        if (nsteps >= m_maxsteps) {
            string f_errs = m_func->getErrors();
            if (!f_errs.empty()) {
                f_errs = "\nExceptions caught during RHS evaluation:\n" + f_errs;
            }
            throw CanteraError("CVodesIntegrator::integrate",
                "Maximum number of timesteps ({}) taken without reaching output "
                "time ({}).\nCurrent integrator time: {}{}",
                nsteps, tout, m_tInteg, f_errs);
        }
        int flag = CVode(m_cvode_mem, tout, m_y, &m_tInteg, CV_ONE_STEP);
        if (flag != CV_SUCCESS) {
            string f_errs = m_func->getErrors();
            if (!f_errs.empty()) {
                f_errs = "Exceptions caught during RHS evaluation:\n" + f_errs;
            }
            throw CanteraError("CVodesIntegrator::integrate",
                "CVodes error encountered. Error code: {}\n{}\n"
                "{}"
                "Components with largest weighted error estimates:\n{}",
                flag, m_error_message, f_errs, getErrorInfo(10));
        }
        nsteps++;
    }
    int flag = CVodeGetDky(m_cvode_mem, tout, 0, m_y);
    checkError(flag, "integrate", "CVodeGetDky");
    m_time = tout;
    m_sens_ok = false;
}

The arguments taken by the SUNDIALS CVode() function are:

  • CVodesIntegrator::m_cvode_mem, a pointer to the block of memory that was allocated and configured during initialization.

  • tout is the desired integrator output time. CVODES will not necessarily reach this time when operating in “one step” mode, but it is used in the selection of the initial step size.

  • After execution, CVodesIntegrator::m_y will contain the computed system state at the time reached by the integrator, and will later be used to update the ReactorNet to its time-integrated state.

  • After execution, CVodesIntegrator::m_tInteg will contain the time reached by the integrator.

  • The CV_ONE_STEP option tells the solver to take a single internal step.

The return value of the CVode() function is assigned to the flag variable. CVode() returns the constant CV_SUCCESS to indicate success an error code if integration was unsuccessful.

Specifying ODEs using Class FuncEval#

How does CVODES know what ODE system it should be solving?

In the example above, the ODE system was actually already specified using CVodeInit(), one of the functions automatically invoked by the ReactorNet::initialize() method. CVODES requires a C function with a specific signature that defines the ODE system by computing 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 the 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 and DAE right-hand-side function evaluators. Classes derived from FuncEval implement the evaluation of the provided ODE system.

An ODE right-hand-side evaluator is always needed in the ODE solution process to provide the definition of the system, and for that reason a FuncEval object is a required parameter when initializing any type of Integrator.

ReactorNet handles this requirement by inheriting the FuncEval class, meaning that it provides the implementation for the ODE function and actually specifies itself using the this pointer when calling Integrator::initialize() in ReactorNet::initialize().

To be a valid FuncEval object, a ReactorNet needs to provide implementations for several of FuncEval’s virtual methods, particularly the actual ODE right-hand-side computation method, FuncEval::eval():

virtual void eval(double t, double* y, double* ydot, double* p);

where the arguments are:

  • t, the current time in seconds.

  • y, a pointer to the start of the current state vector.

  • ydot, a pointer to the start of the array where the computed time derivatives should be stored.

  • p, a pointer to the start of an array containing the (potentially perturbed) sensitivity parameters, if any have been defined.

To take a single timestep, CVODES will call eval() one or more times and use the computed values of ydot to determine the value of y at the new time.

Within ReactorNet::eval(), the governing equations for each each reactor in the network are evaluated and assembled to form the full ydot vector for the system. To handle some complexities of the Reactor model, namely the fact that multiple Reactors can share the same ThermoPhase object while having different states, the governing equation evaluation takes place in two steps. First, the state of each Reactor is set according to the values in the global state vector y using the ReactorNet::updateState() method, which calls the Reactor::updateState() method for each reactor:

void ReactorNet::updateState(double* y)
{
    checkFinite("y", y, m_nv);
    for (size_t n = 0; n < m_reactors.size(); n++) {
        m_reactors[n]->updateState(y + m_start[n]);
    }
}

To simplify implementation of some modifications of the governing equations, for example, using the ExtensibleReactor class, the governing equations for each reactor are written in the form:

\[ \t{LHS}_i \frac{dy_i}{dt} = \t{RHS}_i \]

where the Reactor::eval() method or the eval() method of any class derived from Reactor calculates the values for the LHS (left-hand side) and RHS (right-hand side) vectors, whose values default to 1 and 0, respectively, by implementing a method with the signature:

void eval(double time, double* LHS, double* RHS);

These values are then assembled into the global ydot vector by ReactorNet::eval():

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);
}
On this page
Show Source