Migrate "Science" docs

.github/workflows/main.yml

@@ -312,7 +312,7 @@ jobs:
       - name: Setup Python
         uses: actions/setup-python@v4
         with:
-          python-version: '3.8'
+          python-version: '3.11'
           architecture: x64
       - name: Install Apt dependencies
         run: |
@@ -323,7 +323,7 @@ jobs:
         run: python3 -m pip install -U pip setuptools wheel
       - name: Install Python dependencies
         run: |
-          python3 -m pip install ruamel.yaml scons numpy cython 'sphinx>=5.3.0' \
+          python3 -m pip install ruamel.yaml scons numpy cython 'sphinx>=7.2,<8' \
             sphinxcontrib-matlabdomain sphinxcontrib-doxylink sphinxcontrib-bibtex \
             pydata-sphinx-theme==0.14.1 sphinx-argparse sphinx_design myst-nb \
             sphinx-copybutton matplotlib pandas scipy pint

doc/doxygen/cantera.bib

@@ -20,17 +20,17 @@ @article{bisetti2012
     url = {https://doi.org/10.1016/j.combustflame.2012.08.002},
     doi = {10.1016/j.combustflame.2012.08.002},
     year = {2012}}
-@article{blowers2004,
+@article{blowers2000,
     author = {P.~Blowers and R.~Masel},
     journal = {AIChE Journal},
-    number = {5},
-    pages = {1151--1156},
+    number = {10},
+    pages = {2041--2052},
     title = {Engineering approximations for activation energies in hydrogen transfer
         reactions},
     doi = {10.1002/aic.690461015},
     url = {https://dx.doi.org/10.1002/aic.690461015},
     volume = {46},
-    year = {2004}}
+    year = {2000}}
 @article{chiflikian1995,
     author = {R.~V.~Chiflikian},
     journal = {Physics of Plasmas},
@@ -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},

doc/sphinx/_static/custom.css

@@ -32,3 +32,18 @@ p + div.math {
         object-position: left
     }
 }
+
+/* Make "tip" background color different from "versionadded" */
+.admonition.tip,
+div.admonition.tip {
+  border-color:var(--pst-color-info)
+}
+.admonition.tip>.admonition-title:before,
+div.admonition.tip>.admonition-title:before {
+  background-color:var(--pst-color-info-bg);
+}
+.admonition.tip>.admonition-title:after,
+div.admonition.tip>.admonition-title:after {
+  color:var(--pst-color-info);
+  content:var(--pst-icon-admonition-tip)
+}

doc/sphinx/conf.py

@@ -225,6 +225,15 @@ def escape_splats(app, what, name, obj, options, lines):
 
 myst_enable_extensions = ["dollarmath", "amsmath", "deflist", "colon_fence"]
 
+mathjax3_config = {
+    'tex': {
+        'macros': {
+            't': ['\\mathrm{#1}', 1],
+            'pxpy': ['\\frac{\\partial #1}{\\partial #2}', 2]
+        }
+    }
+}
+
 # Ensure that the primary domain is the Python domain, since we've added the
 # MATLAB domain with sphinxcontrib.matlab
 primary_domain = 'py'

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

doc/sphinx/develop/reactor-integration.md

@@ -0,0 +1,186 @@
+```{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 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 {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 function 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 SUNDIALS
+[`CVode()`](https://sundials.readthedocs.io/en/latest/cvode/Usage/index.html#cvode-solver-function)
+function are:
+
+- {ct}`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, {ct}`CVodesIntegrator::m_y` will contain the computed system state
+  at the time reached by the integrator, and will later be used to update the
+  {ct}`ReactorNet` to its time-integrated state.
+- After execution, {ct}`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 {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 to provide
+the definition of 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 methods, particularly the actual ODE
+right-hand-side computation method, {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:
+
+$$
+\t{LHS}_i \frac{dy_i}{dt} = \t{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 (left-hand side) and RHS (right-hand
+side) vectors, whose values default to 1 and 0, respectively, by implementing a method
+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++
+```

doc/sphinx/reference/index.md

@@ -64,13 +64,41 @@ releasenotes/index
 ## Science Reference
 
 These sections document the scientific theories, mathematical models, and numerical
-methods implemented by Cantera
+methods implemented by Cantera. This starts with some of the fundamental scientific
+theory underpinning the ways that Cantera models phases of matter, which involves
+calculations for thermodynamic and transport properties and chemical reaction rates.
+Cantera then builds on top of these models to provide means of representing and solving
+a number of zero- and one-dimensional systems.
 
 ````{grid} 2
 :gutter: 3
 
+```{grid-item-card} Thermodynamics
+:link: thermo/index
+:link-type: doc
+:text-align: center
+```
+
+```{grid-item-card} Chemical Reactions
+:link: kinetics/index
+:link-type: doc
+:text-align: center
+```
+
+```{grid-item-card} Transport Properties
+:link: transport/index
+:link-type: doc
+:text-align: center
+```
+
 ```{grid-item-card} Reactors
-:link: science/reactors/index
+:link: reactors/index
+:link-type: doc
+:text-align: center
+```
+
+```{grid-item-card} 1D Flames
+:link: onedim/index
 :link-type: doc
 :text-align: center
 ```
@@ -94,7 +122,11 @@ methods implemented by Cantera
 :maxdepth: 1
 :caption: Science Reference
 
-science/reactors/index
+thermo/index
+kinetics/index
+transport/index
+reactors/index
+onedim/index
 bibliography
 glossary
 ```

doc/sphinx/reference/kinetics/index.md

@@ -0,0 +1,24 @@
+# Chemical Reactions
+
+Calculation of reaction rates in Cantera is done in two steps. First, a [*forward rate
+constant*](https://en.wikipedia.org/wiki/Reaction_rate_constant) is calculated, which is
+typically a function of temperature, but can also depend on other properties of the
+mixture state such as pressure and composition. This rate constant is then used along
+with the reactant and product concentrations to determine the forward and reverse rates
+of reactions, as well as other quantities such as the net rate of production for each
+species.
+
+[](reaction-rates)
+: This page describes how forward and reverse reaction rates and other quantities are
+  calculated by Cantera.
+
+[](rate-constants)
+: This page describes the different parameterizations available in Cantera for
+  calculating forward rate constants.
+
+```{toctree}
+:hidden:
+
+reaction-rates
+rate-constants
+```