Clínica dental

# GSoC 2020 Blog Post 3

## [título_original] – Nueva actualización 2023

Somos un pequeño grupo de apasionados profesionales de la odontología que han estado escribiendo artículos para el público durante más de 10 años. Nuestra misión es proporcionar información precisa y actualizada sobre la salud bucal para que las personas puedan tomar decisiones informadas sobre su atención dental.
Sabemos que una excelente salud bucal es esencial para la salud y el bienestar general, y estamos comprometidos a ayudar a nuestros lectores a lograr y mantener dientes y encías saludables.

### 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

Next Post – GSoC 2020 Blog Post 4

Previous Post – GSoC 2020 Blog Post 2

Start from the Beginning – Introduction

## Lea más publicaciones relacionadas [título_original] en la misma categoría

Somos un pequeño grupo de personas apasionadas que han estado escribiendo artículos dentales durante los últimos años. Creemos que una buena salud oral es la clave para una vida feliz y saludable. Nuestro objetivo es proporcionar información precisa y actualizada sobre todos los aspectos de la odontología para que nuestros lectores puedan tomar decisiones informadas sobre su salud bucal.
Sabemos que elegir un dentista o un tratamiento dental puede ser abrumador, pero esperamos que nuestros artículos ayuden a que el proceso sea un poco más fácil.
Categorías principales