From 0d94fc4bbeb5f584e50903f1bacb93278434a836 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Sun, 19 Nov 2023 22:10:01 -0500 Subject: [PATCH] [Doc] Transfer & Revise reactor network documentation Split among science reference, an introductory tutorial, and more detailed information in the developer guide. --- doc/doxygen/cantera.bib | 11 ++ doc/sphinx/develop/index.md | 10 ++ doc/sphinx/develop/reactor-integration.md | 188 ++++++++++++++++++++ doc/sphinx/reference/reactors/index.md | 24 +++ doc/sphinx/userguide/index.md | 2 + doc/sphinx/userguide/reactor-tutorial.md | 200 ++++++++++++++++++++++ 6 files changed, 435 insertions(+) create mode 100644 doc/sphinx/develop/reactor-integration.md create mode 100644 doc/sphinx/userguide/reactor-tutorial.md diff --git a/doc/doxygen/cantera.bib b/doc/doxygen/cantera.bib index 76662c924fc..7954204b943 100644 --- a/doc/doxygen/cantera.bib +++ b/doc/doxygen/cantera.bib @@ -475,6 +475,17 @@ @article{wagner2002 url = {https://dx.doi.org/10.1063/1.1461829}, volume = {31}, year = {2002}} +@article{walker2023, + author = {A.~S.~Walker and R.~L.~Speth and K.~E.~Niemeyer}, + journal = {Proceedings of the Combustion Institute}, + title = {Generalized preconditioning for accelerating simulations with large kinetic models}, + volume = {39}, + number = {4}, + pages = {5395--5403}, + url = {https://www.sciencedirect.com/science/article/pii/S1540748922003273}, + doi = {10.1016/j.proci.2022.07.256}, + month = jun, + year = {2023}} @article{westbrook1981, author = {C.~K.~Westbrook and F.~L.~Dryer}, journal = {Combustion Science and Technology}, diff --git a/doc/sphinx/develop/index.md b/doc/sphinx/develop/index.md index 0216f6a18af..4915c286559 100644 --- a/doc/sphinx/develop/index.md +++ b/doc/sphinx/develop/index.md @@ -35,4 +35,14 @@ compiling/special-cases ## How Cantera Works +- [](reactor-integration) + +```{toctree} +:caption: How Cantera Works +:hidden: +:maxdepth: 1 + +reactor-integration +``` + ## Adding New Models to Cantera diff --git a/doc/sphinx/develop/reactor-integration.md b/doc/sphinx/develop/reactor-integration.md new file mode 100644 index 00000000000..fb5a5881867 --- /dev/null +++ b/doc/sphinx/develop/reactor-integration.md @@ -0,0 +1,188 @@ +```{py:currentmodule} cantera +``` + +# 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 {ct}`ReactorNet` object doesn't perform time integration on its own. Instead, the role +of the {ct}`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. + +```python +# import the Cantera Python module +import cantera as ct + +# create a gas mixture using the GRI-Mech 3.0 mechanism +gas = ct.Solution("gri30.yaml") + +# set gas to an interesting state +gas.TPX = 1000.0, ct.one_atm, "H2:2,O2:1,N2:4" + +# create a reactor containing the gas +reac = ct.IdealGasReactor(gas) + +# Add the reactor to a new ReactorNet simulator +sim = ct.ReactorNet([reac]) +``` + +The `__init__` method of the Python {py:class}`ReactorNet` class calls +{ct}`ReactorNet::addReactor` for each {py:class}`Reactor` object provided in the list +supplied. When the first {ct}`Reactor` is added to the network, the {ct}`ReactorNet` +creates a new {ct}`Integrator` used to integrate the governing equations. + +The {ct}`Integrator` class is Cantera's interface for ODE/DAE system integrators. +{ct}`Integrator` is a [polymorphic base class](http://www.cplusplus.com/doc/tutorial/polymorphism/); +it defines a set of *virtual methods* that derived classes (the actual system +integrators) provide implementations for. + +The {ct}`newIntegrator` factory method creates and returns an {ct}`Integrator` object of +the specified type. Calling `newIntegrator("CVODE")` creates a new +{ct}`CVodesIntegrator` object for integrating an ODE system, while calling +`newIntegrator("IDA")` creates a new {ct}`IdasIntegrator` object for integrating a DAE +system. The appropriate integrator type is determined by the {ct}`ReactorNet` class +based on the types of the installed reactors -- a {ct}`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 {ct}`CVodesIntegrator` is +described; {ct}`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 {ct}`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](https://sundials.readthedocs.io/en/latest/cvodes/Introduction_link.html). + +## Calling the `ReactorNet.advance()` Method + +After configuring a reactor network and its components in Cantera, a call to the +{py:meth}`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 {ct}`Integrator` +when calling `advance()`. Transient physical and chemical interactions are simulated by +integrating the network's system of ODE governing equations through time. + +```python +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 {ct}`ReactorNet::advance` method invokes the method +{ct}`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 {ct}`CVodesIntegrator::integrate` is: + +```{literalinclude} ../../../../src/numerics/CVodesIntegrator.cpp +:language: c++ +:start-at: "void CVodesIntegrator::integrate(" +:end-before: " CVodesIntegrator::" +``` + +The arguments taken by the `CVode()` method are: +- `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, but it is used in the selection of the initial step size. +- After execution, `m_y` will contain the computed solution vector, and will later be + used to update the {ct}`ReactorNet` to its time-integrated state. +- After execution, `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()` method is assigned to the `flag` object. `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 methods automatically invoked during the {ct}`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](https://sundials.readthedocs.io/en/latest/cvodes/Usage/SIM.html#user-supplied-functions). + +The {ct}`CVodesIntegrator` wrapper class provides a useful C++ interface for configuring +this C function by pairing with {ct}`FuncEval`, an abstract base class for ODE and DAE +right-hand-side function evaluators. Classes derived from {ct}`FuncEval` 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 {ct}`FuncEval` object is a +required parameter when initializing any type of {ct}`Integrator`. + +{ct}`ReactorNet` handles this requirement by inheriting the {ct}`FuncEval` class, +meaning that it provides the implementation for the ODE function and actually specifies +itself (using the [`this`](https://en.cppreference.com/w/cpp/language/this) pointer) +when calling {ct}`Integrator::initialize` in {ct}`ReactorNet::initialize`. + +To be a valid {ct}`FuncEval` object, a {ct}`ReactorNet` needs to provide implementations +for several of {ct}`FuncEval`'s virtual functions, particularly the actual ODE +right-hand-side computation function, {ct}`FuncEval::eval`: + +```C++ +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 {ct}`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 +`Reactor`s 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 +{ct}`ReactorNet::updateState` method, which calls the {ct}`Reactor::updateState` method +for each reactor: + +```{literalinclude} ../../../../src/zeroD/ReactorNet.cpp +:start-at: "void ReactorNet::updateState" +:end-before: " ReactorNet::" +:language: c++ +``` + +To simplify implementation of some modifications of the governing equations, for +example, using the {py:class}`ExtensibleReactor` class, the governing equations for each +reactor are written in the form: + +$$ +\mathrm{LHS}_i \frac{dy_i}{dt} = \mathrm{RHS}_i +$$ + +where the {ct}`Reactor::eval` method (or the `eval` method of any class derived from +{ct}`Reactor`) calculates the values for the LHS and RHS vectors, whose values default +to 1 and 0, respectively, by implementing a function with the signature: + +```c++ +void eval(double time, double* LHS, double* RHS); +``` + +These values are then assembled into the global `ydot` vector by `ReactorNet::eval`: + +```{literalinclude} ../../../../src/zeroD/ReactorNet.cpp +:start-at: "void ReactorNet::eval(" +:end-before: " ReactorNet::" +:language: c++ +``` diff --git a/doc/sphinx/reference/reactors/index.md b/doc/sphinx/reference/reactors/index.md index d1f0aff912d..7ed416969fd 100644 --- a/doc/sphinx/reference/reactors/index.md +++ b/doc/sphinx/reference/reactors/index.md @@ -123,6 +123,30 @@ time-dependent reactors. : A reactor modeling one-dimensional steady-state flow in a channel that may contain catalytically active surfaces where heterogeneous reactions occur. +## Reactor Networks + +While reactors by themselves define the governing equations, the time integration is +performed by assembling reactors into a reactor network. A reactor network is therefore +necessary even if only a single reactor is considered. + +Cantera uses the CVODES and IDAS solvers from the +[SUNDIALS](https://computing.llnl.gov/projects/sundials) package to integrate the +governing equations for the reactor network, which are a system of stiff ODEs or DAEs. + +### Preconditioning + +Some of Cantera's reactor formulations (specifically, the +[Ideal Gas Control Volume Mole Reactor](ideal-gas-mole-reactor) and the +[Ideal Gas Constant Pressure Mole Reactor](ideal-gas-constant-pressure-mole-reactor)) +provide implementations of a sparse, approximate Jacobian matrix, which can be used by +{ct}`ReactorNet` to generate a preconditioner and use a sparse, iterative linear solver +within CVODES. This sparse, preconditioned method can significantly accelerate +integration for reactors containing many species. A derivation of the derivative terms +and benchmarks demonstrating the achievable performance gains can be found in +{cite:t}`walker2023`. An example demonstrating the use of this feature can be found in +[`preconditioned_integration.py`](/examples/python/reactors/preconditioned_integration). + + ```{toctree} :hidden: :caption: Reactor models diff --git a/doc/sphinx/userguide/index.md b/doc/sphinx/userguide/index.md index 21a7e1fe07f..4a9ae1c502c 100644 --- a/doc/sphinx/userguide/index.md +++ b/doc/sphinx/userguide/index.md @@ -14,6 +14,7 @@ troubleshooting. - [Getting Started with Python](python-tutorial) - [](input-tutorial) +- [](reactor-tutorial) ## Frequently asked questions @@ -44,6 +45,7 @@ conditions, or calculating the voltage of a Lithium-ion battery as it is dischar python-tutorial input-tutorial +reactor-tutorial faq diff --git a/doc/sphinx/userguide/reactor-tutorial.md b/doc/sphinx/userguide/reactor-tutorial.md new file mode 100644 index 00000000000..7e8c26c6ecd --- /dev/null +++ b/doc/sphinx/userguide/reactor-tutorial.md @@ -0,0 +1,200 @@ +--- +file_format: mystnb +kernelspec: + name: python3 +--- + +```{py:currentmodule} cantera +``` + +# Using Reactor Networks for Time Integration + +This guide explains the basics of setting up and integrating the governing equations of +a transient reactor network problem. + +## Setting up a 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 filled +with a homogeneous gas mixture (the gas state used in this example is arbitrary, but +interesting because it's explosive). + +```{code-cell} python +# import the Cantera Python module +import cantera as ct + +# create a gas mixture using the GRI-Mech 3.0 mechanism +gas = ct.Solution("gri30.yaml") + +# set gas to an interesting state +gas.TPX = 1000.0, ct.one_atm, "H2:2,O2:1,N2:4" + +# create a reactor containing the gas +reac = ct.IdealGasReactor(gas) + +# Add the reactor to a new ReactorNet simulator +sim = ct.ReactorNet([reac]) + +# View the initial state of the mixture +reac.thermo() +``` + +Now, let's advance the simulation in time to an (arbitrary) absolute time of 1 second, +and examine how the state of the gas has changed + +```{code-cell} python +# advance the simulation to the specified absolute time, t = 1 sec +sim.advance(1) + +# view the updated state of the mixture, reflecting properties at t = 1 sec +reac.thermo() +``` + +## Methods for Time Integration + +The state of a {py:class}`ReactorNet` can be advanced in time by one of the following +three methods: + +- `step()`: The `step()` method computes the state of the system after one time step. + The size of the step is determined by SUNDIALS when the method is called. SUNDIALS + determines the step size by estimating the local error, which must satisfy tolerance + conditions. The step is redone with reduced step size whenever that error test fails. + SUNDIALS also periodically checks if the maximum step size is being used. The time + step must not be larger than a predefined maximum time step $\Delta t_{\mathrm{max}}$. + The new time $t_{\mathrm{new}}$ at the end of the single step is returned by this + function. This method produces the highest time resolution in the output data of the + methods implemented in Cantera. + +- `advance(t_new)`: This method computes the state of the system at the user-provided + time $t_{\mathrm{new}}$. $t_{\mathrm{new}}$ is the absolute time from the initial time + of the system. Although the user specifies the time when integration should stop, + SUNDIALS chooses the time step size as the network is integrated. Many internal + SUNDIALS time steps are usually required to reach $t_{\mathrm{new}}$. As such, + `advance(t_new)` preserves the accuracy of using `step()` but allows consistent + spacing in the output data. + +- `advance_to_steady_state(max_steps, residual_threshold, atol, write_residuals)` + *Python interface only*: If the steady state solution of a reactor network is of + interest, this method can be used. Internally, the steady state is approached by time + stepping. The network is considered to be at steady state if the feature-scaled + residual of the state vector is below a given threshold value (which by default is 10 + times the time step `rtol`). + +### Limiting Timesteps + +The `advance(t_new)` approach is typically used when consistent, regular, time steps are +required in the output data. This usually simplifies comparisons among many simulation +results at a single time point. However, some detail, for example, a fast ignition +process, might not be resolved in the output data if the user-provided time step is not +small enough. + +To avoid losing this detail, the {ct}`Reactor::setAdvanceLimit` method (C++) or the +{py:func}`Reactor.set_advance_limit` method (Python) can be used to set the maximum +amount that a specified solution component can change between output times. For an +example of this feature's use, see the example +[reactor1.py](/examples/python/reactors/reactor1). However, as a tradeoff, the time step +sizes in the output data are no longer guaranteed to be uniform. + +### Setting Tolerances + +Even though Cantera comes pre-defined with reasonable default values for tolerances and +the maximum internal time step, the solver may not be correctly configured for the +specific system. In this case SUNDIALS may stop with an error. To solve this problem, +three parameters can be tuned: The absolute tolerances, the relative tolerances, and the +maximum time step. A reduction of the latter value is particularly useful when dealing +with abrupt changes in the boundary conditions (for example, opening/closing valves; see +also the [IC engine example](/examples/python/reactors/ic_engine)). + +## Adding Inlets and Outlets + +The following examples demonstrate the use of [flow devices](sec-flow-device) such +as valves and mass flow controllers to model the inlets and outlets of reactors: + +- [](/examples/python/reactors/mix1) +- [](/examples/python/reactors/fuel_injection) +- [](/examples/python/reactors/ic_engine) +- [](/examples/cxx/combustor) (C++) + +## Modeling Volume Changes and Heat Transfer + +The following examples demonstrate the use of [walls](sec-wall) to model volume changes +and heat transfer between reactors: + +- [](/examples/python/reactors/piston) +- [](/examples/python/reactors/reactor2) +- [](/examples/python/reactors/ic_engine) + +Additional examples can be found in the +[Python Reactor Network Examples](/examples/python/reactors/index) section. + +## Accelerating Integration with Preconditioned Sparse Iterative Methods + +Some reactor models (specifically, the [](/reference/reactors/ideal-gas-mole-reactor) +and the [](/reference/reactors/ideal-gas-constant-pressure-mole-reactor)) are able to be +solved using preconditioned sparse iterative methods, instead of the normal direct +linear methods. To enable the use of this solver, a preconditioner object needs to be +created and specified for use by the reactor network. The previous example can be +modified as follows: + +```{code-cell} 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.IdealGasMoleReactor(gas) # Use reactor type that supports preconditioning +sim = ct.ReactorNet([reac]) +precon = ct.AdaptivePreconditioner() # Create the preconditioner +sim.preconditioner = precon # Add it to the network +``` + +The results of time time integration are the same as before: +```{code-cell} python +sim.advance(1) +reac.thermo() +``` + +## Common Reactor Types and their Implementation in Cantera + +### Batch Reactor at Constant Volume or at Constant Pressure + +If you are interested in how a homogeneous chemical composition changes in time when it +is left to its own devices, a simple batch reactor can be used. Two versions are +commonly considered: A rigid vessel with fixed volume but variable pressure, or a system +idealized at constant pressure but varying volume. + +In Cantera, such a simulation can be performed very easily. The initial state of the +solution can be specified by composition and a set of thermodynamic parameters (like +temperature and pressure) as a standard Cantera {py:class}`Solution` object. From this +base, a {py:class}`Reactor` or a {py:class}`ConstPressureReactor` can be created, +depending on if a constant volume or constant pressure batch reactor should be +considered, respectively. The behavior of the solution in time can be simulated as a +very simple {py:class}`ReactorNet` containing only the single created reactor. + +An example for such a batch reactor is given in +[`reactor1.py`](/examples/python/reactors/reactor1). + +### Continuously Stirred Tank Reactor + +A Continuously Stirred Tank Reactor (CSTR), also often referred to as Well-Stirred +Reactor (WSR), Perfectly Stirred Reactor (PSR), or Longwell Reactor, is essentially a +single Cantera reactor with an inlet, an outlet, and constant volume. Therefore, the +governing equations for single reactors defined above apply accordingly. + +Steady state solutions to CSTRs are often of interest. In this case, the mass flow rate +$\dot{m}$ is constant and equal at inlet and outlet. The mass contained in the +confinement $m$ divided by $\dot{m}$ defines the mean residence time of the fluid in the +confinement. + +While Cantera always solves a transient problem, if you are interested in steady-state +conditions, you can run your simulation for a long time until the states are converged. +This integration process can be managed automatically using the +{py:meth}`ReactorNet.advance_to_steady_state` method. See +[`combustor.py`](/examples/python/reactors/combustor) for an example using this method. + +:::{tip} +A problem can be the ignition of a CSTR: If the reactants are not reactive enough, the +simulation can result in the trivial solution that inflow and outflow states are +identical. To solve this problem, the reactor can be initialized with a high temperature +and/or radical concentration. A good approach is to use the equilibrium composition of +the reactants (which can be computed using Cantera's {py:meth}`ThermoPhase.equilibrate` +function) as an initial guess. +:::