diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index e1a5a2becf..469f7dcaf6 100644
--- a/.github/workflows/main.yml
+++ b/.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
diff --git a/doc/doxygen/cantera.bib b/doc/doxygen/cantera.bib
index f8988d599c..7954204b94 100644
--- a/doc/doxygen/cantera.bib
+++ b/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},
diff --git a/doc/sphinx/_static/custom.css b/doc/sphinx/_static/custom.css
index 6dda8f69a2..cc79008e90 100644
--- a/doc/sphinx/_static/custom.css
+++ b/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)
+}
diff --git a/doc/sphinx/conf.py b/doc/sphinx/conf.py
index 760679f8ba..2aa47f712a 100644
--- a/doc/sphinx/conf.py
+++ b/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'
diff --git a/doc/sphinx/develop/index.md b/doc/sphinx/develop/index.md
index 0216f6a18a..4915c28655 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 0000000000..d2bda4bdd3
--- /dev/null
+++ b/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++
+```
diff --git a/doc/sphinx/reference/index.md b/doc/sphinx/reference/index.md
index 2f6faff810..ef3f285722 100644
--- a/doc/sphinx/reference/index.md
+++ b/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
```
diff --git a/doc/sphinx/reference/kinetics/index.md b/doc/sphinx/reference/kinetics/index.md
new file mode 100644
index 0000000000..cc8e99bad4
--- /dev/null
+++ b/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
+```
diff --git a/doc/sphinx/reference/kinetics/rate-constants.md b/doc/sphinx/reference/kinetics/rate-constants.md
new file mode 100644
index 0000000000..5ada2ed91d
--- /dev/null
+++ b/doc/sphinx/reference/kinetics/rate-constants.md
@@ -0,0 +1,381 @@
+# Rate Constant Parameterizations
+
+This page describes the different parameterizations available in Cantera for
+calculating the forward rate constant $k_f$ for a reaction.
+
+(sec-arrhenius-rate)=
+## Arrhenius Rate Expressions
+
+An Arrhenius rate is described by the
+[modified Arrhenius function](https://en.wikipedia.org/wiki/Arrhenius_equation#Modified_Arrhenius_equation):
+
+$$ k_f = A T^b e^{-E_a / RT} $$
+
+where $A$ is the pre-exponential factor, $T$ is the temperature, $b$ is the temperature
+exponent, $E_a$ is the activation energy, and $R$ is the universal gas constant.
+
+:::{admonition} YAML Usage
+:class: tip
+An Arrhenius rate can be specified for a reaction in the YAML format by providing an
+[Arrhenius](sec-yaml-Arrhenius-rate) rate expression for the reaction's `rate-constant`
+field.
+:::
+
+(sec-falloff-rate)=
+## Falloff Reactions
+
+A falloff reaction is one that has a rate that is first-order in the total concentration
+of third-body colliders $[\t{M}]$ at low pressure, like a
+[three-body reaction](sec-three-body-reaction), but becomes zero-order in $[\t{M}]$ as
+$[\t{M}]$ increases. Dissociation/association reactions of polyatomic molecules often
+exhibit this behavior.
+
+The simplest expression for the rate coefficient for a falloff reaction is the Lindemann
+form {cite:p}`lindemann1922`:
+
+$$ k_f(T, [\t{M}]) = \frac{k_0 [\t{M}]}{1 + \frac{k_0 [\t{M}]}{k_\infty}} $$
+
+In the low-pressure limit, this approaches $k_0 [\t{M}]$, and in the high-pressure limit
+it approaches $k_\infty$.
+
+Defining the non-dimensional reduced pressure:
+
+$$ P_r = \frac{k_0 [\t{M}]}{k_\infty} $$
+
+The rate constant may be written as
+
+$$ k_f(T, P_r) = k_\infty \left(\frac{P_r}{1 + P_r}\right) $$
+
+More accurate models for unimolecular processes lead to other, more complex,
+forms for the dependence on reduced pressure. These can be accounted for by
+multiplying the Lindemann expression by a function $F(T, P_r)$:
+
+$$ k_f(T, P_r) = k_\infty \left(\frac{P_r}{1 + P_r}\right) F(T, P_r) $$
+
+This expression is used to compute the rate coefficient for falloff reactions. The
+function $F(T, P_r)$ is the falloff function.
+
+```{admonition} YAML Usage
+:class: tip
+A falloff reaction may be defined in the YAML format using the
+[`falloff`](sec-yaml-falloff) reaction `type`.
+```
+
+(sec-troe-falloff)=
+### The Troe Falloff Function
+
+A widely-used falloff function is the one proposed by {cite:t}`gilbert1983`:
+
+\begin{gather*}
+\log_{10} F(T, P_r) = \frac{\log_{10} F_\t{cent}(T)}{1 + f_1^2} \\
+
+F_\t{cent}(T) = (1-A) \exp(-T/T_3) + A \exp (-T/T_1) + \exp(-T_2/T) \\
+
+f_1 = (\log_{10} P_r + C) / (N - 0.14 (\log_{10} P_r + C)) \\
+
+C = -0.4 - 0.67\; \log_{10} F_\t{cent} \\
+
+N = 0.75 - 1.27\; \log_{10} F_\t{cent}
+\end{gather*}
+
+```{admonition} YAML Usage
+:class: tip
+A Troe falloff function may be specified in the YAML format using the
+[`Troe`](sec-yaml-falloff) field in the reaction entry. The first three parameters, $(A,
+T_3, T_1)$, are required. The fourth parameter, $T_2$, is optional; if omitted, the last
+term of the falloff function is not used.
+```
+
+(sec-tsang-falloff)=
+### Tsang's Approximation to $F_\t{cent}$
+
+Wing Tsang presented approximations for the value of $F_\t{cent}$ for Troe falloff in
+databases of reactions, for example, {cite:t}`tsang1991`. Tsang's approximations are
+linear in temperature:
+
+$$ F_\t{cent} = A + BT $$
+
+where $A$ and $B$ are constants. The remaining equations for $C$, $N$, $f_1$, and $F$
+from the [Troe](sec-troe-falloff) falloff function are not affected.
+
+```{admonition} YAML Usage
+:class: tip
+A Tsang falloff function may be specified in the YAML format using the
+[`Tsang`](sec-yaml-falloff) field in the reaction entry.
+```
+
+```{versionadded} 2.6
+```
+
+(sec-sri-falloff)=
+### The SRI Falloff Function
+
+This falloff function is based on the one originally due to {cite:t}`stewart1989`, which
+required three parameters $a$, $b$, and $c$. {cite:t}`kee1989` generalized this function
+slightly by adding two more parameters $d$ and $e$. The original form corresponds to
+$d = 1$ and $e = 0$. Cantera supports the extended 5-parameter form, given by:
+
+$$ F(T, P_r) = d \bigl[a \exp(-b/T) + \exp(-T/c)\bigr]^{1/(1+\log_{10}^2 P_r )} T^e $$
+
+In keeping with the nomenclature of {cite:t}`kee1989`, we will refer to this as the
+**SRI falloff function**.
+
+```{admonition} YAML Usage
+:class: tip
+An SRI falloff function may be specified in the YAML format using the
+[`SRI`](sec-yaml-falloff) field in the entry.
+```
+
+(sec-chemically-activated-rate)=
+## Chemically-Activated Reactions
+
+For these reactions, the rate falls off as the pressure increases, due to collisional
+stabilization of a reaction intermediate. Example:
+
+$$ \t{Si + SiH_4 (+M) \leftrightarrow Si_2H_2 + H_2 (+M)} $$
+
+which competes with:
+
+$$ \t{Si + SiH_4 (+M) \leftrightarrow Si_2H_4 (+M)} $$
+
+Like falloff reactions, chemically-activated reactions are described by blending between
+a low-pressure and a high-pressure rate expression. The difference is that the forward
+rate constant is written as proportional to the *low-pressure* rate constant:
+
+$$ k_f(T, P_r) = k_0 \left(\frac{1}{1 + P_r}\right) F(T, P_r) $$
+
+and the optional blending function $F$ may be described by any of the parameterizations
+allowed for falloff reactions.
+
+```{admonition} YAML Usage
+:class: tip
+Chemically-activated reactions can be defined in the YAML format using the
+[`chemically-activated`](sec-yaml-chemically-activated) reaction `type`.
+```
+
+(sec-plog-rate)=
+## Pressure-Dependent Arrhenius Rate Expressions (P-Log)
+
+This parameterization represents pressure-dependent reaction rates by logarithmically
+interpolating between Arrhenius rate expressions at various pressures. Given two rate
+expressions at two specific pressures:
+
+$$
+P_1: k_1(T) = A_1 T^{b_1} e^{-E_1 / RT}
+
+P_2: k_2(T) = A_2 T^{b_2} e^{-E_2 / RT}
+$$
+
+The rate at an intermediate pressure $P_1 < P < P_2$ is computed as
+
+$$
+\log k(T,P) = \log k_1(T) + \bigl(\log k_2(T) - \log k_1(T)\bigr)
+ \frac{\log P - \log P_1}{\log P_2 - \log P_1}
+$$
+
+Multiple rate expressions may be given at the same pressure, in which case the rate used
+in the interpolation formula is the sum of all the rates given at that pressure. For
+pressures outside the given range, the rate expression at the nearest pressure is used.
+
+```{caution}
+Negative A-factors can be used for any of the rate expressions at a given pressure.
+However, the sum of all of the rates at a given pressure **must** be positive, due to
+the logarithmic interpolation of the rate for intermediate pressures. When a P-log type
+reaction is initialized, Cantera does a validation check for a range of temperatures
+that the sum of the reaction rates at each pressure is positive. Unfortunately, if these
+checks fail, the only options are to remove the reaction or contact the author of the
+reaction/mechanism in question, because the reaction is mathematically unsound.
+```
+
+```{admonition} YAML Usage
+:class: tip
+P-log reactions can be defined in the YAML format using the
+[`pressure-dependent-Arrhenius`](sec-yaml-pressure-dependent-Arrhenius) reaction `type`.
+```
+
+(sec-chebyshev-rate)=
+## Chebyshev Reaction Rate Expressions
+
+Chebyshev rate expressions represent a phenomenological rate coefficient $k(T,P)$ in
+terms of a bivariate Chebyshev polynomial. The rate constant can be written as:
+
+$$
+\log k(T,P) = \sum_{t=1}^{N_T} \sum_{p=1}^{N_P} \alpha_{tp}
+ \phi_t(\tilde{T}) \phi_p(\tilde{P})
+$$
+
+where $N_T$ is the order of the polynomial in the temperature dimension, $N_P$ is the
+order of the polynomial in the pressure dimension, $\alpha_{tp}$ are the constants
+defining the rate, $\phi_n(x)$ is the Chebyshev polynomial of the first kind of degree
+$n$ evaluated at $x$, and
+
+$$
+\tilde{T} \equiv \frac{2T^{-1} - T_\t{min}^{-1} - T_\t{max}^{-1}}
+ {T_\t{max}^{-1} - T_\t{min}^{-1}}
+
+\tilde{P} \equiv \frac{2 \log P - \log P_\t{min} - \log P_\t{max}}
+ {\log P_\t{max} - \log P_\t{min}}
+$$
+
+are reduced temperatures and reduced pressures which map the ranges $(T_\t{min},
+T_\t{max})$ and $(P_\t{min}, P_\t{max})$ to $(-1, 1)$.
+
+A Chebyshev rate expression is specified in terms of the coefficient matrix $\alpha$ and
+the temperature and pressure ranges.
+
+```{caution}
+The Chebyshev polynomials are not defined outside the interval $(-1,1)$, and therefore
+extrapolation of rates outside the range of temperatures and pressure for which they are
+defined is strongly discouraged.
+```
+
+```{admonition} YAML Usage
+:class: tip
+Chebyshev reactions can be defined in the YAML format using the
+[`Chebyshev`](sec-yaml-Chebyshev) reaction `type`.
+```
+
+(sec-blowers-masel)=
+## Blowers-Masel Reactions
+
+In some circumstances like thermodynamic sensitivity analysis, or modeling heterogeneous
+reactions from one catalyst surface to another, the enthalpy change of a reaction
+($\Delta H$) can be modified. Due to the change in $\Delta H$, the activation energy of
+the reaction must be adjusted accordingly to provide accurate simulation results. To
+adjust the activation energy due to changes in the reaction enthalpy, the Blowers-Masel
+rate expression is available. This approximation was proposed by {cite:t}`blowers2000`
+to automatically scale activation energy as the reaction enthalpy is changed. The
+_intrinsic activation energy_ $E_a^0$ is defined as the activation energy when
+$\Delta H = 0$. The activation energy can then be written as a function of $\Delta H$:
+
+$$
+E_a = \begin{cases}
+ 0 & \text{if } \Delta H \leq -4 E_a^0 \\
+ \Delta H & \text{if } \Delta H \geq 4 E_a^0 \\
+ \frac{\left( w + \frac{\Delta H }{2} \right) (V_P - 2 w + \Delta H) ^2}
+ {V_P^2 - 4 w^2 + \Delta H^2} & \text{otherwise}
+ \end{cases}
+$$
+
+where
+
+$$ V_P = 2 w \frac{w + E_a^0}{w - E_a^0}, $$
+
+and $w$ is the average of the bond dissociation energy of the bond breaking and that
+being formed. Note that the expression is insensitive to $w$ as long as $w \ge 2 E_a^0$,
+so we can use an arbitrarily high value of $w = 1000\text{ kJ/mol}$.
+
+After $E_a$ is evaluated, the reaction rate can be calculated using the modified
+Arrhenius expression
+
+$$ k_f = A T^b e^{-E_a / RT}. $$
+
+```{versionadded} 2.6
+```
+
+```{admonition} YAML Usage
+:class: tip
+Blowers Masel reactions can be defined in the YAML format using the
+[Blowers-Masel](sec-yaml-Blowers-Masel-rate) reaction `type`.
+```
+
+(sec-surface-rate)=
+## Surface Reactions
+
+Heterogeneous reactions on surfaces are represented by an extended Arrhenius- like rate
+expression, which combines the modified Arrhenius rate expression with further
+corrections dependent on the fractional surface coverages $\theta_k$ of one or more
+surface species. The forward rate constant for a reaction of this type is:
+
+$$
+k_f = A T^b \exp \left( - \frac{E_a}{RT} \right)
+ \prod_k 10^{a_k \theta_k}
+ \theta_k^{m_k}
+ \exp \left( \frac{- E_k \theta_k}{RT} \right)
+$$
+
+where $A$, $b$, and $E_a$ are the modified Arrhenius parameters and $a_k$, $m_k$, and
+$E_k$ are the coverage dependencies from species $k$.
+
+::::{admonition} YAML Usage
+:class: tip
+In the YAML format, surface reactions are identified by the presence of surface species
+and support several [additional options](sec-yaml-interface-Arrhenius).
+
+The surface reaction `type` defaults to `interface-Arrhenius`, where the rate
+expression uses the [`Arrhenius`](sec-elementary) parameterization (see [YAML
+documentation](sec-yaml-interface-Arrhenius)).
+
+:::{versionadded} 2.6
+As an alternative, Cantera also supports the `interface-Blowers-Masel` surface reaction
+`type`, which uses the [`Blowers-Masel`](sec-Blowers-Masel) parameterization (see [YAML
+documentation](sec-yaml-interface-Blowers-Masel)).
+:::
+::::
+
+(sec-sticking-rate)=
+## Sticking Reactions
+
+Sticking reactions represent a special case of surface reactions, where collisions
+between gas-phase molecules and surfaces result in the gas-phase molecule sticking to
+the surface. This process can be described as a reaction which is parameterized by a
+sticking coefficient:
+
+$$ \gamma = a T^b e^{-c/RT} $$
+
+where $a$, $b$, and $c$ are constants specific to the reaction. The values of these
+constants must be specified so that the sticking coefficient $\gamma$ is between 0 and 1
+for all temperatures.
+
+The sticking coefficient is related to the forward rate constant by the formula:
+
+$$ k_f = \frac{\gamma}{\Gamma_\t{tot}^m} \sqrt{\frac{RT}{2 \pi W}} $$
+
+where $\Gamma_\t{tot}$ is the total molar site density, $m$ is the sum of all the
+surface reactant stoichiometric coefficients, and $W$ is the molecular weight of the gas
+phase species.
+
+::::{admonition} YAML Usage
+:class: tip
+Sticking reactions can be defined in the YAML format by specifying the rate constant
+in the reaction's
+[sticking-coefficient](sec-yaml-sticking-Arrhenius) field.
+
+The sticking reaction `type` defaults to `sticking-Arrhenius`, where the rate expression
+uses the [`Arrhenius`](sec-elementary) parameterization (see [YAML
+documentation](sec-yaml-sticking-Arrhenius)).
+
+:::{versionadded} 2.6
+As an alternative, Cantera also supports the `sticking-Blowers-Masel` surface reaction
+`type`, which uses the [`Blowers-Masel`](sec-Blowers-Masel) parameterization (see
+[YAML documentation](sec-yaml-sticking-Blowers-Masel)).
+:::
+::::
+
+(sec-two-temperature-plasma-rate)=
+## Two-Temperature-Plasma Reactions
+
+The two-temperature-plasma reaction is commonly used for non-equilibrium plasmas. The
+reaction rate of a two-temperature-plasma reaction depends on both gas and electron
+temperature {cite:p}`kossyi1992`, and can be expressed as:
+
+$$
+k_f = A T_e^b \exp \left( - \frac{E_{a,g}}{RT} \right)
+ \exp \left(\frac{E_{a,e}(T_e - T)}{R T T_e}\right),
+$$
+
+where $A$ is the pre-exponential factor, $T$ is the temperature, $T_e$ is the electron
+temperature, $b$ is the electron temperature exponent, $E_{a,g}$ is the activation
+energy for gas, $E_{a,e}$ is the activation energy for electron and $R$ is the gas
+constant.
+
+```{versionadded} 2.6
+```
+
+:::{admonition} YAML Usage
+:class: tip
+
+Two-temperature plasma reactions can be defined in the YAML format by specifying
+[`two-temperature-plasma`](sec-yaml-two-temperature-plasma) as the reaction `type` and
+providing the two activation energies as part of the `rate-constant`.
+:::
diff --git a/doc/sphinx/reference/kinetics/reaction-rates.md b/doc/sphinx/reference/kinetics/reaction-rates.md
new file mode 100644
index 0000000000..becff62574
--- /dev/null
+++ b/doc/sphinx/reference/kinetics/reaction-rates.md
@@ -0,0 +1,145 @@
+# Reaction Rates
+
+Here, we describe how Cantera calculates chemical reaction rates for various reaction
+types.
+
+(sec-elementary)=
+## Elementary Reactions
+
+The basic reaction type is a homogeneous reaction with a pressure-independent
+rate coefficient and mass action kinetics. For example:
+
+$$ a\t{A} + b\t{B} \rightleftharpoons c\t{C} + d\t{D} $$
+
+where A and B are reactant species, C and D are product species, and $a, b, c, $ and $d$
+are stoichiometric coefficients.
+
+The forward reaction rate is then calculated as:
+
+$$ R_f = [\t{A}]^a [\t{B}]^b k_f $$
+
+where $k_f$ is the forward rate constant, calculated using one of the available rate
+parameterizations such as the [modified Arrhenius](sec-arrhenius-rate) form.
+
+```{admonition} YAML Usage
+:class: tip
+An elementary reaction with an Arrhenius reaction rate can be defined in the YAML format
+using the [`elementary`](sec-yaml-elementary) reaction `type`, or by omitting the
+reaction `type` entry, as it represents the default. An exception to this default is
+when the same species occurs on both sides of the reaction equation, in which case the
+reaction is treated as a
+[three-body reaction for a specific collider](sec-three-body-specific-collider).
+```
+
+(sec-three-body-reaction)=
+## Three-Body Reactions
+
+A three-body reaction is a gas-phase reaction of the form:
+
+$$ \t{A + B + M \rightleftharpoons AB + M} $$
+
+Here $\t{M}$ is an unspecified collision partner that carries away excess energy to
+stabilize the $\t{AB}$ molecule (forward direction) or supplies energy to break the
+$\t{AB}$ bond (reverse direction). In addition to the generic collision partner
+$\t{M}$, it is also possible to explicitly specify a colliding species. In both
+cases, the reaction type can be automatically inferred by Cantera and does not need to
+be explicitly specified by the user.
+
+Different species may be more or less effective in acting as the collision partner. A
+species that is much lighter than $\t{A}$ and $\t{B}$ may not be able to
+transfer much of its kinetic energy, and so would be inefficient as a collision partner.
+On the other hand, a species with a transition from its ground state that is nearly
+resonant with one in the $\t{AB^*}$ activated complex may be much more effective at
+exchanging energy than would otherwise be expected.
+
+These effects can be accounted for by defining a collision efficiency $\epsilon$ for
+each species, defined such that the forward reaction rate is
+
+$$ R_f = [\t{A}][\t{B}][\t{M}] k_f(T) $$
+
+where
+
+$$ [\t{M}] = \sum_{k} \epsilon_k C_k $$
+
+where $C_k$ is the concentration of species $k$. Since any constant collision efficiency
+can be absorbed into the rate coefficient $k_f(T)$, the default collision efficiency is
+1.0.
+
+:::{versionadded} 3.0
+The rate coefficient $k_f(T)$ may be implemented using any rate parameterization
+supported by Cantera, not just the modified Arrhenius form.
+:::
+
+(sec-three-body-specific-collider)=
+### Collider-specific rate parameterizations
+
+Sometimes, accounting for a particular third body's collision efficiency may require an
+alternate set of rate parameters entirely. In this case, two reactions are written:
+
+$$
+\t{A + B + M \rightleftharpoons AB + M \quad (R1)}
+
+\t{A + B + C \rightleftharpoons AB + C \quad (R2)}
+$$
+
+where the third-body efficiency for C in the first reaction should be explicitly set to
+zero. For the second reaction, the efficiencies will automatically be set to one for C
+and zero for all other colliders.
+
+```{admonition} YAML Usage
+:class: tip
+A three-body reaction may be defined in the YAML format using the
+[`three-body`](sec-yaml-three-body) reaction `type` or, if no `type` is specified,
+identified automatically by the presence of the generic third body M or a specific
+non-reactive species (for example, C in R2 above).
+```
+
+## Pressure-dependent Reactions
+
+For pressure-dependent reactions where the behavior is more complex than described
+by the three-body form, the pressure dependency is folded into the calculation of the
+rate constant. Cantera supports several ways of representing pressure-dependent
+reactions:
+
+- [](sec-falloff-rate)
+- [](sec-chemically-activated-rate)
+- [](sec-plog-rate)
+- [](sec-chebyshev-rate)
+
+(sec-reaction-orders)=
+## Reaction Orders
+
+Explicit reaction orders different from the stoichiometric coefficients are sometimes
+used for non-elementary reactions. For example, consider the global reaction:
+
+$$
+\t{C_8H_{18} + 12.5 O_2 \rightarrow 8 CO_2 + 9 H_2O}
+$$
+
+the forward rate constant might be given as {cite:p}`westbrook1981`:
+
+$$
+k_f = 4.6 \times 10^{11} [\t{C_8H_{18}}]^{0.25} [\t{O_2}]^{1.5}
+ \exp\left(\frac{30.0\,\t{kcal/mol}}{RT}\right)
+$$
+
+Special care is required in this case since the units of the pre-exponential factor
+depend on the sum of the reaction orders, which may not be an integer.
+
+Note that you can change reaction orders only for irreversible reactions.
+
+Normally, reaction orders are required to be positive. However, in some cases negative
+reaction orders are found to be better fits for experimental data. In these cases, the
+default behavior may be overridden in the input file.
+
+````{admonition} YAML Usage
+:class: tip
+To include explicit orders for the reaction above, it can be written in the YAML format as:
+
+```yaml
+- equation: C8H18 + 12.5 O2 => 8 CO2 + 9 H2O
+ units: {length: cm, quantity: mol, activation-energy: kcal/mol}
+ rate-constant: {A: 4.5e+11, b: 0.0, Ea: 30.0}
+ orders: {C8H18: 0.25, O2: 1.5}
+```
+````
diff --git a/doc/sphinx/reference/onedim/discretization.md b/doc/sphinx/reference/onedim/discretization.md
new file mode 100644
index 0000000000..40abd8cea6
--- /dev/null
+++ b/doc/sphinx/reference/onedim/discretization.md
@@ -0,0 +1,5 @@
+# Discretization of 1D Equations
+
+```{caution}
+This page is a work in progress.
+```
diff --git a/doc/sphinx/reference/onedim/governing-equations.md b/doc/sphinx/reference/onedim/governing-equations.md
new file mode 100644
index 0000000000..142d64afb0
--- /dev/null
+++ b/doc/sphinx/reference/onedim/governing-equations.md
@@ -0,0 +1,223 @@
+# Governing Equations for One-dimensional Flow
+
+Cantera models flames that are stabilized in an axisymmetric stagnation flow, and
+computes the solution along the stagnation streamline ($r=0$), using a similarity
+solution to reduce the three-dimensional governing equations to a single dimension.
+
+## Axisymmetric Stagnation Flow
+
+The governing equations for a steady axisymmetric stagnation flow follow those derived
+in Section 7.2 of {cite:t}`kee2017` and are implemented by class {ct}`StFlow`.
+
+*Continuity*:
+
+$$ \pxpy{u}{z} + 2 \rho V = 0 $$
+
+*Radial momentum*:
+
+$$
+\rho u \pxpy{V}{z} + \rho V^2 = - \Lambda + \pxpy{}{z}\left(\mu \pxpy{V}{z}\right)
+$$
+
+*Energy*:
+
+$$
+\rho c_p u \pxpy{T}{z} = \pxpy{}{z}\left(\lambda \pxpy{T}{z}\right)
+ - \sum_k j_k \pxpy{h_k}{z} - \sum_k h_k W_k \dot{\omega}_k
+$$
+
+*Species*:
+
+$$
+\rho u \pxpy{Y_k}{z} = - \pxpy{j_k}{z} + W_k \dot{\omega}_k
+$$
+
+where the following variables are used:
+
+- $z$ is the axial coordinate
+- $r$ is the radial coordinate
+- $\rho$ is the density
+- $u$ is the axial velocity
+- $v$ is the radial velocity
+- $V = v/r$ is the scaled radial velocity
+- $\Lambda$ is the pressure eigenvalue (independent of $z$)
+- $\mu$ is the dynamic viscosity
+- $c_p$ is the heat capacity at constant pressure
+- $T$ is the temperature
+- $\lambda$ is the thermal conductivity
+- $Y_k$ is the mass fraction of species $k$
+- $j_k$ is the diffusive mass flux of species $k$
+- $c_{p,k}$ is the specific heat capacity of species $k$
+- $h_k$ is the enthalpy of species $k$
+- $W_k$ is the molecular weight of species $k$
+- $\dot{\omega}_k$ is the molar production rate of species $k$.
+
+The tangential velocity $w$ has been assumed to be zero. The model is applicable to both
+ideal and non-ideal fluids, which follow ideal-gas or real-gas (Redlich-Kwong and
+Peng-Robinson) equations of state.
+
+```{versionadded} 3.0
+Support for real gases in the flame models was introduced in Cantera 3.0.
+```
+
+To help in the solution of the discretized problem, it is useful to write a
+differential equation for the scalar $\Lambda$:
+
+$$ \frac{d\Lambda}{dz} = 0 $$
+
+When discretized, the Jacobian terms introduced by this equation match the block
+diagonal structure produced by the other governing equations, rather than creating a
+column of entries that would cause fill-in when factorizing as part of the Newton
+solver.
+
+## Diffusive Fluxes
+
+The species diffusive mass fluxes $j_k$ are computed according to either a
+mixture-averaged or multicomponent formulation. If the mixture-averaged formulation is
+used, the calculation performed is:
+
+$$
+j_k^* = - \rho \frac{W_k}{\overline{W}} D_{km}^\prime \pxpy{X_k}{z}
+
+j_k = j_k^* - Y_k \sum_i j_i^*
+$$
+
+where $\overline{W}$ is the mean molecular weight of the mixture, $D_{km}^\prime$ is the
+mixture-averaged diffusion coefficient for species $k$, and $X_k$ is the mole fraction
+for species $k$. The diffusion coefficients used here are those computed by the method
+{ct}`GasTransport::getMixDiffCoeffs`. The correction applied by the second equation
+ensures that the sum of the mass fluxes is zero, a condition which is not inherently
+guaranteed by the mixture-averaged formulation.
+
+When using the multicomponent formulation, the mass fluxes are computed according to:
+
+$$
+j_k = \frac{\rho W_k}{\overline{W}^2} \sum_i W_i D_{ki} \pxpy{X_i}{z}
+ - \frac{D_k^T}{T} \pxpy{T}{z}
+$$
+
+where $D_{ki}$ is the multicomponent diffusion coefficient and $D_k^T$ is the Soret
+diffusion coefficient. Inclusion of the Soret calculation must be explicitly enabled
+when setting up the simulation, on top of specifying a multicomponent transport model,
+for example by using the {ct}`StFlow::enableSoret` method (C++) or setting the
+{py:attr}`~cantera.FlameBase.soret_enabled` property (Python).
+
+## Boundary Conditions
+
+### Inlet boundary
+
+For a boundary located at a point $z_0$ where there is an inflow, values are supplied
+for the temperature $T_0$, the species mass fractions $Y_{k,0}$ the scaled radial
+velocity $V_0$, and the mass flow rate $\dot{m}_0$. In the case of the
+freely-propagating flame, the mass flow rate is not an input but is determined
+indirectly by holding the temperature fixed at an intermediate location within the
+domain; see [](discretization) for details.
+
+The following equations are solved at the point $z = z_\t{in}$:
+
+$$
+T(z_\t{in}) &= T_0
+
+V(z_\t{in}) &= V_0
+
+\dot{m}_0 Y_{k,\t{in}} - j_k(z_\t{in}) - \rho(z_\t{in}) u(z_\t{in}) Y_k(z_\t{in}) &= 0
+$$
+
+If the mass flow rate is specified, we also solve:
+
+$$
+\rho(z_\t{in}) u(z_\t{in}) = \dot{m}_0
+$$
+
+Otherwise, we solve:
+
+$$ \Lambda(z_\t{in}) = 0 $$
+
+These equations are implemented by class {ct}`Inlet1D`.
+
+### Outlet boundary
+
+For a boundary located at a point $z_\t{out}$ where there is an outflow, we solve:
+
+$$
+\Lambda(z_\t{out}) = 0
+
+\left.\pxpy{T}{z}\right|_{z_\t{out}} = 0
+
+\left.\pxpy{Y_k}{z}\right|_{z_\t{out}} = 0
+
+V(z_\t{out}) = 0
+$$
+
+These equations are implemented by class {ct}`Outlet1D`.
+
+### Symmetry boundary
+
+For a symmetry boundary located at a point $z_\t{symm}$, we solve:
+
+$$
+\rho(z_\t{symm}) u(z_\t{symm}) = 0
+
+\left.\pxpy{V}{z}\right|_{z_\t{symm}} = 0
+
+\left.\pxpy{T}{z}\right|_{z_\t{symm}} = 0
+
+j_k(z_\t{symm}) = 0
+$$
+
+These equations are implemented by class {ct}`Symm1D`.
+
+### Reacting surface
+
+For a surface boundary located at a point $z_\t{surf}$ on which reactions may
+occur, the temperature $T_\t{surf}$ is specified. We solve:
+
+$$
+\rho(z_\t{surf}) u(z_\t{surf}) &= 0
+
+V(z_\t{surf}) &= 0
+
+T(z_\t{surf}) &= T_\t{surf}
+
+j_k(z_\t{surf}) + \dot{s}_k W_k &= 0
+$$
+
+where $\dot{s}_k$ is the molar production rate of the gas-phase species $k$ on the
+surface. In addition, the surface coverages $\theta_i$ for each surface species $i$ are
+computed such that $\dot{s}_i = 0$.
+
+These equations are implemented by class {ct}`ReactingSurf1D`.
+
+## The Drift-Diffusion Model
+
+To account for the transport of charged species in a flame, class {ct}`IonFlow` adds the
+drift term to the diffusive fluxes of the mixture-average formulation according to
+{cite:t}`pedersen1993`,
+
+$$
+j_k^* = \rho \frac{W_k}{\overline{W}} D_{km}^\prime \pxpy{X_k}{z} + s_k \mu_k E Y_k,
+$$
+
+where $s_k$ is the sign of charge (1,-1, and 0 respectively for positive, negative, and
+neutral charge), $\mu_k$ is the mobility, and $E$ is the electric field. The diffusion
+coefficients and mobilities of charged species can be more accurately calculated by
+{ct}`IonGasTransport::getMixDiffCoeffs` and {ct}`IonGasTransport::getMobilities`. The
+following correction is applied instead to preserve the correct fluxes of charged
+species:
+
+$$
+j_k = j_k^* - \frac {1 - |s_k|} {1 - \sum_i |s_i| Y_i} Y_k \sum_i j_i^*.
+$$
+
+In addition, Gauss's law is solved simultaneously with the species and energy equations,
+
+$$
+\pxpy{E}{z} &= \frac{e}{\epsilon_0}\sum_k Z_k n_k ,
+
+n_k &= N_a \rho Y_k / W_k,
+
+E|_{z=0} &= 0,
+$$
+
+where $Z_k$ is the charge number, $n_k$ is the number density, and $N_a$ is the Avogadro
+number.
diff --git a/doc/sphinx/reference/onedim/index.md b/doc/sphinx/reference/onedim/index.md
new file mode 100644
index 0000000000..2e0f25ec03
--- /dev/null
+++ b/doc/sphinx/reference/onedim/index.md
@@ -0,0 +1,39 @@
+# One-dimensional Flames
+
+Cantera includes a set of models for representing steady-state, quasi-one-dimensional
+reacting flows.
+
+These models can be used to simulate a number of common flames, such as:
+
+- freely-propagating premixed laminar flames
+- burner-stabilized premixed flames
+- counterflow diffusion flames
+- counterflow (strained) premixed flames
+
+Additional capabilities include simulation of surface reactions, which can be used to
+represent processes such as combustion on a catalytic surface or chemical vapor
+deposition processes.
+
+All of these configurations are simulated using a common set of governing equations
+within a 1D flow domain, with the differences between the models being represented by
+differences in the boundary conditions applied.
+
+[](governing-equations)
+: This page describes the governing equations and the various boundary conditions that
+ can be applied.
+
+[](discretization)
+: This page describes the finite difference schemes used to discretize the 1D governing
+ equations and the criteria used for refining the grid.
+
+[](nonlinear-solver)
+: This page describes the hybrid time-stepping--steady-state damped Newton solver that
+ is used to solve the discretized governing equations.
+
+```{toctree}
+:hidden:
+
+governing-equations
+discretization
+nonlinear-solver
+```
diff --git a/doc/sphinx/reference/onedim/nonlinear-solver.md b/doc/sphinx/reference/onedim/nonlinear-solver.md
new file mode 100644
index 0000000000..77f6f3a949
--- /dev/null
+++ b/doc/sphinx/reference/onedim/nonlinear-solver.md
@@ -0,0 +1,5 @@
+# Nonlinear Solver for One-dimensional Flows
+
+```{caution}
+This page is a work in progress.
+```
diff --git a/doc/sphinx/reference/reactors/constant-pressure-mole-reactor.md b/doc/sphinx/reference/reactors/constant-pressure-mole-reactor.md
new file mode 100644
index 0000000000..422ea4f254
--- /dev/null
+++ b/doc/sphinx/reference/reactors/constant-pressure-mole-reactor.md
@@ -0,0 +1,55 @@
+```{py:currentmodule} cantera
+```
+
+# Constant Pressure Mole Reactor
+
+A constant pressure mole reactor is implemented by the C++ class
+{ct}`ConstPressureMoleReactor` and is available in Python as the
+{py:class}`ConstPressureMoleReactor` class. It is defined by the state variables:
+
+- $H$, the total enthalpy of the reactor's contents (in J)
+- $n_k$, the number of moles for each species (in kmol)
+
+Equations 1 and 2 below are the governing equations for a constant pressure mole
+reactor.
+
+## Species Equations
+
+The moles of each species in the reactor changes as a result of flow through the
+reactor's [inlets and outlets](sec-flow-device) and production of gas phase species
+through homogeneous reactions and reactions on the reactor
+[surfaces](sec-reactor-surface). The rate at which species $k$ is generated through
+homogeneous phase reactions is $V \dot{\omega}_k$, and the total rate at which moles of
+species $k$ changes is:
+
+$$
+\frac{dn_k}{dt} = V \dot{\omega}_k + \sum_\t{in} \dot{n}_{k, \t{in}}
+ - \sum_\t{out} \dot{n}_{k, \t{out}} + \dot{n}_{k, \t{wall}}
+$$ (const-pressure-mole-reactor-species)
+
+Where the subscripts *in* and *out* refer to the sum of the corresponding property over
+all inlets and outlets respectively. A dot above a variable signifies a time derivative.
+
+## Energy Equation
+
+Writing the first law for an open system gives:
+
+$$
+\frac{dU}{dt} = - p \frac{dV}{dt} + \dot{Q} + \sum_\t{in} \dot{n}_\t{in} \hat{h}_\t{in}
+ - \hat{h} \sum_\t{out} \dot{n}_\t{out}
+$$
+
+where positive $\dot{Q}$ represents heat addition to the system and $h$ is the specific
+enthalpy of the reactor's contents.
+
+Differentiating the definition of the total enthalpy, $H = U + pV$, with respect to time
+gives:
+
+$$ \frac{dH}{dt} = \frac{dU}{dt} + p \frac{dV}{dt} + V \frac{dp}{dt} $$
+
+Noting that $dp/dt = 0$ and substituting into the energy equation yields:
+
+$$
+\frac{dH}{dt} = \dot{Q} + \sum_\t{in} \dot{n}_\t{in} \hat{h}_\t{in}
+ - \hat{h} \sum_\t{out} \dot{n}_\t{out}
+$$ (const-pressure-mole-reactor-energy)
diff --git a/doc/sphinx/reference/reactors/constant-pressure-reactor.md b/doc/sphinx/reference/reactors/constant-pressure-reactor.md
new file mode 100644
index 0000000000..1067249de2
--- /dev/null
+++ b/doc/sphinx/reference/reactors/constant-pressure-reactor.md
@@ -0,0 +1,76 @@
+```{py:currentmodule} cantera
+```
+
+# Constant Pressure Reactor
+
+For this reactor model, the pressure is held constant and the energy equation is defined
+in terms of the total enthalpy. This model is implemented by the C++ class
+{ct}`ConstPressureReactor` and available in Python as the
+{py:class}`ConstPressureReactor` class. A constant pressure reactor is defined by the
+state variables:
+
+- $m$, the mass of the reactor's contents (in kg)
+- $H$, the total enthalpy of the reactor's contents (in J)
+- $Y_k$, the mass fractions for each species (dimensionless)
+
+Equations 1-3 below are the governing equations for a constant pressure reactor.
+
+## Mass Conservation
+
+The total mass of the reactor's contents changes as a result of flow through the
+reactor's [inlets and outlets](sec-flow-device), and production of homogeneous phase
+species on [surfaces](sec-reactor-surface):
+
+$$
+\frac{dm}{dt} = \sum_\t{in} \dot{m}_\t{in} - \sum_\t{out} \dot{m}_\t{out} +
+ \dot{m}_\t{wall}
+$$ (constpressurereactor-mass)
+
+Where the subscripts *in* and *out* refer to the sum of the superscripted property over
+all inlets and outlets respectively. A dot above a variable signifies a time derivative.
+
+## Species Equations
+
+The rate at which species $k$ is generated through homogeneous phase reactions is $V
+\dot{\omega}_k W_k$, and the total rate at which species $k$ is generated is:
+
+$$ \dot{m}_{k,\t{gen}} = V \dot{\omega}_k W_k + \dot{m}_{k,\t{wall}} $$
+
+The rate of change in the mass of each species is:
+
+$$
+\frac{d(mY_k)}{dt} = \sum_\t{in} \dot{m}_\t{in} Y_{k,\t{in}}
+ - \sum_\t{out} \dot{m}_\t{out} Y_k + \dot{m}_{k,\t{gen}}
+$$
+
+Expanding the derivative on the left hand side and substituting the equation for
+$dm/dt$, the equation for each homogeneous phase species is:
+
+$$
+m \frac{dY_k}{dt} = \sum_\t{in} \dot{m}_\t{in} (Y_{k,\t{in}} - Y_k)
+ + \dot{m}_{k,\t{gen}} - Y_k \dot{m}_\t{wall}
+$$ (constpressurereactor-species)
+
+## Energy Equation
+
+Writing the first law for an open system gives:
+
+$$
+\frac{dU}{dt} = - p \frac{dV}{dt} + \dot{Q}
+ + \sum_\t{in} \dot{m}_\t{in} h_\t{in} - h \sum_\t{out} \dot{m}_\t{out}
+$$
+
+where positive $\dot{Q}$ represents heat addition to the system and $h$ is the specific
+enthalpy of the reactor's contents.
+
+Differentiating the definition of the total enthalpy, $H = U + pV$, with respect to time
+gives:
+
+$$ \frac{dH}{dt} = \frac{dU}{dt} + p \frac{dV}{dt} + V \frac{dp}{dt} $$
+
+Noting that $dp/dt = 0$ and substituting into the energy equation yields:
+
+$$
+\frac{dH}{dt} = \dot{Q} + \sum_\t{in} \dot{m}_\t{in} h_\t{in}
+ - h \sum_\t{out} \dot{m}_\t{out}
+$$ (constpressurereactor-energy)
diff --git a/doc/sphinx/reference/science/reactors/controlreactor.md b/doc/sphinx/reference/reactors/controlreactor.md
similarity index 58%
rename from doc/sphinx/reference/science/reactors/controlreactor.md
rename to doc/sphinx/reference/reactors/controlreactor.md
index ea0535125f..d2cd0d2713 100644
--- a/doc/sphinx/reference/science/reactors/controlreactor.md
+++ b/doc/sphinx/reference/reactors/controlreactor.md
@@ -4,31 +4,34 @@
# Control Volume Reactor
This model represents a homogeneous zero-dimensional reactor, as implemented by the C++
-class {ct}`Reactor`. By default, these reactors are closed (no inlets or outlets), have
-fixed volume, and have adiabatic, chemically-inert walls. These properties may all be
-changed by adding appropriate components such as {py:class}`Wall`,
-{py:class}`ReactorSurface`, {py:class}`MassFlowController`, and
-{py:class}`Valve`.
-
-A control volume reactor is defined by the four state variables:
+class {ct}`Reactor` and available in Python as the {py:class}`Reactor` class. A control
+volume reactor is defined by the state variables:
- $m$, the mass of the reactor's contents (in kg)
-- $V$, the reactor volume (in m$^3$)
+- $V$, the reactor volume (in m{sup}`3`)
- $U$, the total internal energy of the reactors contents (in J)
- $Y_k$, the mass fractions for each species (dimensionless)
+Equations 1-4 below are the governing equations for a control volume reactor.
+
+## Mass Conservation
+
The total mass of the reactor's contents changes as a result of flow through the
-reactor's inlets and outlets, and production of homogeneous phase species on the reactor
-{py:class}`Wall`.
+reactor's [inlets and outlets](sec-flow-device), and production of homogeneous phase
+species on [surfaces](sec-reactor-surface):
$$
-\frac{dm}{dt} = \sum_{in} \dot{m}_{in} - \sum_{out} \dot{m}_{out} + \dot{m}_{wall}
+\frac{dm}{dt} = \sum_\t{in} \dot{m}_\t{in} - \sum_\t{out} \dot{m}_\t{out}
+ + \dot{m}_\t{wall}
$$ (mass)
Where the subscripts *in* and *out* refer to the sum of the corresponding property over
all inlets and outlets respectively. A dot above a variable signifies a time derivative.
-The reactor volume changes as a function of time due to the motion of one or more walls:
+## Volume Equation
+
+The reactor volume changes as a function of time due to the motion of one or more
+[walls](sec-wall):
$$
\frac{dV}{dt} = \sum_w f_w A_w v_w(t)
@@ -38,36 +41,38 @@ where $f_w = \pm 1$ indicates the facing of the wall (whether moving the wall in
or decreases the volume of the reactor), $A_w$ is the surface area of the wall, and
$v_w(t)$ is the velocity of the wall as a function of time.
-The equation for the total internal energy is found by writing the first law for an open
-system:
-
-$$
-\frac{dU}{dt} = - p \frac{dV}{dt} + \dot{Q} +
- \sum_{in} \dot{m}_{in} h_{in} - h \sum_{out} \dot{m}_{out}
-$$ (energy)
-
-Where $\dot{Q}$ is the net rate of heat addition to the system.
+## Species Equations
The rate at which species $k$ is generated through homogeneous phase reactions is
$V \dot{\omega}_k W_k$, and the total rate at which species $k$ is generated is:
$$
-\dot{m}_{k,gen} = V \dot{\omega}_k W_k + \dot{m}_{k,wall}
+\dot{m}_{k,\t{gen}} = V \dot{\omega}_k W_k + \dot{m}_{k,\t{wall}}
$$
The rate of change in the mass of each species is:
$$
-\frac{d(mY_k)}{dt} = \sum_{in} \dot{m}_{in} Y_{k,in} - \sum_{out} \dot{m}_{out} Y_k +
- \dot{m}_{k,gen}
+\frac{d(mY_k)}{dt} = \sum_\t{in} \dot{m}_\t{in} Y_{k,\t{in}}
+ - \sum_\t{out} \dot{m}_\t{out} Y_k + \dot{m}_{k,\t{gen}}
$$
Expanding the derivative on the left hand side and substituting the equation
for $dm/dt$, the equation for each homogeneous phase species is:
$$
-m \frac{dY_k}{dt} = \sum_{in} \dot{m}_{in} (Y_{k,in} - Y_k) +
- \dot{m}_{k,gen} - Y_k \dot{m}_{wall}
+m \frac{dY_k}{dt} = \sum_\t{in} \dot{m}_\t{in} (Y_{k,\t{in}} - Y_k) +
+ \dot{m}_{k,\t{gen}} - Y_k \dot{m}_\t{wall}
$$ (species)
-Equations 1-4 are the governing equations for a control volume reactor.
+## Energy Equation
+
+The equation for the total internal energy is found by writing the first law for an open
+system:
+
+$$
+\frac{dU}{dt} = - p \frac{dV}{dt} + \dot{Q} +
+ \sum_\t{in} \dot{m}_\t{in} h_\t{in} - h \sum_\t{out} \dot{m}_\t{out}
+$$ (cv-energy)
+
+Where $\dot{Q}$ is the net rate of heat addition to the system.
diff --git a/doc/sphinx/reference/reactors/ideal-gas-constant-pressure-mole-reactor.md b/doc/sphinx/reference/reactors/ideal-gas-constant-pressure-mole-reactor.md
new file mode 100644
index 0000000000..a0309a347a
--- /dev/null
+++ b/doc/sphinx/reference/reactors/ideal-gas-constant-pressure-mole-reactor.md
@@ -0,0 +1,71 @@
+```{py:currentmodule} cantera
+```
+
+# Ideal Gas Constant Pressure Mole Reactor
+
+An ideal gas constant pressure mole reactor, as implemented by the C++ class
+{ct}`IdealGasConstPressureMoleReactor` and available in Python as the
+{py:class}`IdealGasConstPressureMoleReactor` class. It is defined by the state
+variables:
+
+- $T$, the temperature (in K)
+- $n_k$, the number of moles for each species (in kmol)
+
+Equations 1 and 2 below are the governing equations for an ideal gas constant pressure
+mole reactor.
+
+## Species Equations
+
+The moles of each species in the reactor changes as a result of flow through the
+reactor's [inlets and outlets](sec-flow-device), and production of homogeneous gas phase
+species and reactions on the reactor [surfaces](sec-reactor-surface). The rate at which
+species $k$ is generated through homogeneous phase reactions is $V \dot{\omega}_k$, and
+the total rate at which moles of species $k$ changes is:
+
+$$
+\frac{dn_k}{dt} = V \dot{\omega}_k + \sum_\t{in} \dot{n}_{k, \t{in}}
+ - \sum_\t{out} \dot{n}_{k, \t{out}} + \dot{n}_{k, \t{wall}}
+$$ (ig-const-pressure-mole-reactor-species)
+
+Where the subscripts *in* and *out* refer to the sum of the corresponding property over
+all inlets and outlets respectively. A dot above a variable signifies a time derivative.
+
+## Energy Equation
+
+Writing the first law for an open system gives:
+
+$$
+\frac{dU}{dt} = - p \frac{dV}{dt} + \dot{Q} + \sum_\t{in} \dot{n}_\t{in} \hat{h}_\t{in}
+ - \hat{h} \sum_\t{out} \dot{n}_\t{out}
+$$
+
+where positive $\dot{Q}$ represents heat addition to the system and $h$ is the specific
+enthalpy of the reactor's contents.
+
+Differentiating the definition of the total enthalpy, $H = U + pV$, with respect to time
+gives:
+
+$$ \frac{dH}{dt} = \frac{dU}{dt} + p \frac{dV}{dt} + V \frac{dp}{dt} $$
+
+Noting that $dp/dt = 0$ and substituting into the energy equation yields:
+
+$$
+\frac{dH}{dt} = \dot{Q} + \sum_\t{in} \dot{n}_\t{in} \hat{h}_\t{in}
+ - \hat{h} \sum_\t{out} \dot{n}_\t{out}
+$$
+
+As for the [ideal gas mole reactor](ideal-gas-mole-reactor), we replace the total
+enthalpy as a state variable with the temperature by writing the total enthalpy in terms
+of the species moles and temperature:
+
+$$ H = \sum_k \hat{h}_k(T) n_k $$
+
+and differentiating with respect to time:
+
+$$ \frac{dH}{dt} = \frac{dT}{dt}\sum_k n_k \hat{c}_{p,k} + \sum \hat{h}_k \dot{n}_k $$
+
+Making this substitution and rearranging yields an equation for the temperature:
+
+$$
+\sum_k n_k \hat{c}_{p,k} \frac{dT}{dt} = \dot{Q} - \sum \hat{h}_k \dot{n}_k
+$$ (ig-const-pressure-mole-reactor-energy)
diff --git a/doc/sphinx/reference/reactors/ideal-gas-constant-pressure-reactor.md b/doc/sphinx/reference/reactors/ideal-gas-constant-pressure-reactor.md
new file mode 100644
index 0000000000..96ca482259
--- /dev/null
+++ b/doc/sphinx/reference/reactors/ideal-gas-constant-pressure-reactor.md
@@ -0,0 +1,72 @@
+```{py:currentmodule} cantera
+```
+
+# Ideal Gas Constant Pressure Reactor
+
+An ideal gas constant pressure reactor, as implemented by the C++ class
+{ct}`IdealGasConstPressureReactor` and available in Python as the
+{py:class}`IdealGasConstPressureReactor` class. It is defined by the state variables:
+
+- $m$, the mass of the reactor's contents (in kg)
+- $T$, the temperature (in K)
+- $Y_k$, the mass fractions for each species (dimensionless)
+
+Equations 1-3 below are the governing equations for an ideal gas constant pressure
+reactor.
+
+## Mass Conservation
+
+The total mass of the reactor's contents changes as a result of flow through the
+reactor's [inlets and outlets](sec-flow-device), and production of homogeneous phase
+species on [surfaces](sec-reactor-surface):
+
+$$
+\frac{dm}{dt} = \sum_\t{in} \dot{m}_\t{in} - \sum_\t{out} \dot{m}_\t{out}
+ + \dot{m}_\t{wall}
+$$ (igcpr-mass)
+
+Where the subscripts *in* and *out* refer to the sum of the corresponding property over
+all inlets and outlets respectively. A dot above a variable signifies a time derivative.
+
+## Species Equations
+
+The rate at which species $k$ is generated through homogeneous phase reactions is
+$V \dot{\omega}_k W_k$, and the total rate at which species $k$ is generated is:
+
+$$ \dot{m}_{k,\t{gen}} = V \dot{\omega}_k W_k + \dot{m}_{k,\t{wall}} $$
+
+The rate of change in the mass of each species is:
+
+$$
+\frac{d(mY_k)}{dt} = \sum_\t{in} \dot{m}_\t{in} Y_{k,\t{in}}
+ - \sum_\t{out} \dot{m}_\t{out} Y_k + \dot{m}_{k,gen}
+$$
+
+Expanding the derivative on the left hand side and substituting the equation
+for $dm/dt$, the equation for each homogeneous phase species is:
+
+$$
+m \frac{dY_k}{dt} = \sum_\t{in} \dot{m}_\t{in} (Y_{k,\t{in}} - Y_k)
+ + \dot{m}_{k,\t{gen}} - Y_k \dot{m}_\t{wall}
+$$ (igcpr-species)
+
+## Energy Equation
+
+As for the [ideal gas reactor](ideal-gas-reactor), we replace the total enthalpy as a
+state variable with the temperature by writing the total enthalpy in terms of the mass
+fractions and temperature and differentiating with respect to time:
+
+$$
+H &= m \sum_k Y_k h_k(T)
+
+\frac{dH}{dt} &= h \frac{dm}{dt} + m c_p \frac{dT}{dt}
+ + m \sum_k h_k \frac{dY_k}{dt}
+$$
+
+Substituting the corresponding derivatives into the constant pressure reactor energy
+equation {eq}`constpressurereactor-energy` yields an equation for the temperature:
+
+$$
+m c_p \frac{dT}{dt} = \dot{Q} - \sum_k h_k \dot{m}_{k,\t{gen}}
+ + \sum_\t{in} \dot{m}_\t{in} \left(h_\t{in} - \sum_k h_k Y_{k,\t{in}} \right)
+$$ (igcpr-energy)
diff --git a/doc/sphinx/reference/reactors/ideal-gas-mole-reactor.md b/doc/sphinx/reference/reactors/ideal-gas-mole-reactor.md
new file mode 100644
index 0000000000..fb26f00dca
--- /dev/null
+++ b/doc/sphinx/reference/reactors/ideal-gas-mole-reactor.md
@@ -0,0 +1,60 @@
+```{py:currentmodule} cantera
+```
+
+# Ideal Gas Control Volume Mole Reactor
+
+An ideal gas control volume mole reactor, as implemented by the C++ class
+{ct}`IdealGasMoleReactor` and available in Python as the {py:class}`IdealGasMoleReactor`
+class. It is defined by the state variables:
+
+- $T$, the temperature (in K)
+- $V$, the reactor volume (in m{sup}`3`)
+- $n_k$, the number of moles for each species (in kmol)
+
+Equations 1-3 are the governing equations for an ideal gas control volume mole reactor.
+
+## Volume Equation
+
+The reactor volume can change as a function of time due to the motion of one or more
+[walls](sec-wall):
+
+$$
+\frac{dV}{dt} = \sum_w f_w A_w v_w(t)
+$$ (ig-mole-reactor-volume)
+
+Where $f_w = \pm 1$ indicates the facing of the wall (whether moving the wall increases
+or decreases the volume of the reactor), $A_w$ is the surface area of the wall, and
+$v_w(t)$ is the velocity of the wall as a function of time.
+
+## Species Equations
+
+The moles of each species in the reactor changes as a result of flow through the
+reactor's [inlets and outlets](sec-flow-device), and production of homogeneous gas phase
+species and reactions on the reactor [surfaces](sec-reactor-surface). The rate at which
+species $k$ is generated through homogeneous phase reactions is $V \dot{\omega}_k$, and
+the total rate at which moles of species $k$ changes is:
+
+$$
+\frac{dn_k}{dt} = V \dot{\omega}_k + \sum_\t{in} \dot{n}_{k, \t{in}}
+ - \sum_\t{out} \dot{n}_{k, \t{out}} + \dot{n}_{k, \t{wall}}
+$$ (ig-mole-reactor-species)
+
+## Energy Equation
+
+In the case of the ideal gas control volume mole reactor model, the reactor temperature
+$T$ is used instead of the total internal energy $U$ as a state variable. For an ideal
+gas, we can rewrite the total internal energy in terms of the species moles and
+temperature:
+
+$$ U = \sum_k \hat{u}_k(T) n_k $$
+
+and differentiate it with respect to time to obtain:
+
+$$ \frac{dU}{dt} = \frac{dT}{dt}\sum_k n_k \hat{c}_{v,k} + \sum \hat{u}_k \dot{n}_k $$
+
+Substituting this into the energy equation for the control volume mole reactor
+{eq}`molereactor-energy` yields an equation for the temperature:
+
+$$
+\sum_k n_k \hat{c}_{v,k} \frac{dT}{dt} = \dot{Q} - \sum \hat{u}_k \dot{n}_k
+$$ (ig-mole-reactor-energy)
diff --git a/doc/sphinx/reference/reactors/ideal-gas-reactor.md b/doc/sphinx/reference/reactors/ideal-gas-reactor.md
new file mode 100644
index 0000000000..26cdf24d1c
--- /dev/null
+++ b/doc/sphinx/reference/reactors/ideal-gas-reactor.md
@@ -0,0 +1,90 @@
+```{py:currentmodule} cantera
+```
+
+# Ideal Gas Control Volume Reactor
+
+An ideal gas control volume reactor, as implemented by the C++ class
+{ct}`IdealGasReactor` and available in Python as the {py:class}`IdealGasReactor` class.
+It is defined by the state variables:
+
+- $m$, the mass of the reactor's contents (in kg)
+- $V$, the reactor volume (in m{sup}`3`)
+- $T$, the temperature (in K)
+- $Y_k$, the mass fractions for each species (dimensionless)
+
+Equations 1-4 below are the governing equations for an ideal gas reactor.
+
+## Mass Conservation
+
+The total mass of the reactor's contents changes as a result of flow through the
+reactor's [inlets and outlets](sec-flow-device), and production of gas phase species on
+[surfaces](sec-reactor-surface):
+
+$$
+\frac{dm}{dt} = \sum_\t{in} \dot{m}_\t{in} - \sum_\t{out} \dot{m}_\t{out}
+ + \dot{m}_\t{wall}
+$$ (igr-mass)
+
+where the subscripts *in* and *out* refer to the sum of the corresponding property over
+all inlets and outlets respectively. A dot above a variable signifies a time derivative.
+
+## Volume Equation
+
+The reactor volume can change as a function of time due to the motion of one or more
+[walls](sec-wall):
+
+$$
+\frac{dV}{dt} = \sum_w f_w A_w v_w(t)
+$$ (igr-volume)
+
+where $f_w = \pm 1$ indicates the facing of the wall (whether moving the wall increases
+or decreases the volume of the reactor), $A_w$ is the surface area of the wall, and
+$v_w(t)$ is the velocity of the wall as a function of time.
+
+## Species Equations
+
+The rate at which species $k$ is generated through homogeneous phase reactions is
+$V \dot{\omega}_k W_k$, and the total rate at which species $k$ is generated is:
+
+$$ \dot{m}_{k,\t{gen}} = V \dot{\omega}_k W_k + \dot{m}_{k,\t{wall}} $$
+
+The rate of change in the mass of each species is:
+
+$$
+\frac{d(mY_k)}{dt} = \sum_\t{in} \dot{m}_\t{in} Y_{k,\t{in}} - \sum_\t{out} \dot{m}_\t{out} Y_k +
+ \dot{m}_{k,\t{gen}}
+$$
+
+Expanding the derivative on the left hand side and substituting the equation
+for $dm/dt$, the equation for each homogeneous phase species is:
+
+$$
+m \frac{dY_k}{dt} = \sum_\t{in} \dot{m}_\t{in} (Y_{k,\t{in}} - Y_k)+ \dot{m}_{k,\t{gen}}
+ - Y_k \dot{m}_\t{wall}
+$$ (igr-species)
+
+## Energy Equation
+
+In the case of the ideal gas control volume reactor model, the reactor temperature $T$
+is used instead of the total internal energy $U$ as a state variable. For an ideal gas,
+we can rewrite the total internal energy in terms of the mass fractions and temperature:
+
+$$ U = m \sum_k Y_k u_k(T) $$
+
+and differentiate it with respect to time to obtain:
+
+$$
+\frac{dU}{dt} = u \frac{dm}{dt} + m c_v \frac{dT}{dt} + m \sum_k u_k \frac{dY_k}{dt}
+$$
+
+Substituting this into the energy equation for the control volume reactor
+{eq}`cv-energy` yields an equation for the temperature:
+
+$$
+m c_v \frac{dT}{dt} =& - p \frac{dV}{dt} + \dot{Q} + \sum_\t{in} \dot{m}_\t{in} \left( h_\t{in} - \sum_k u_k Y_{k,\t{in}} \right) \\
+ &- \frac{p V}{m} \sum_\t{out} \dot{m}_\t{out} - \sum_k \dot{m}_{k,\t{gen}} u_k
+$$ (igr-energy)
+
+While this form of the energy equation is somewhat more complicated, it significantly
+reduces the cost of evaluating the system Jacobian, since the derivatives of the species
+equations are taken at constant temperature instead of constant internal energy.
diff --git a/doc/sphinx/reference/reactors/index.md b/doc/sphinx/reference/reactors/index.md
new file mode 100644
index 0000000000..c2a67107fe
--- /dev/null
+++ b/doc/sphinx/reference/reactors/index.md
@@ -0,0 +1,166 @@
+```{py:currentmodule} cantera
+```
+
+# Reactors and Reactor Networks
+
+```{caution}
+This page is a work in progress. For more complete documentation of the current Cantera
+release (Cantera 3.0), please see this page.
+```
+
+In Cantera, a *reactor network* represents a set of one or more homogeneous reactors and
+reacting surfaces that may be connected to each other and to the environment through
+devices representing mass flow, heat transfer, and moving walls. The system is generally
+unsteady -- that is, all states are functions of time. In particular, transient state
+changes due to chemical reactions are possible.
+
+## Homogeneous Reactor Types and Governing Equations
+
+A Cantera *reactor* represents the simplest form of a chemically reacting system. It
+corresponds to an extensive thermodynamic control volume $V$, in which all state
+variables are homogeneously distributed. By default, these reactors are closed (no
+inlets or outlets) and have adiabatic, chemically-inert walls. These properties may all
+be changed by adding components such as walls, surfaces, mass flow controllers, and
+valves, as described [below](sec-reactor-interactions).
+
+The specific governing equations defining Cantera's reactor models are derived and
+described below. These models represent different combinations of whether pressure or
+volume are held constant, whether they support any equation of state or are specialized
+for ideal gas mixtures, and whether mass fractions or moles of each species are used as
+the state variables representing the composition.
+
+[Control Volume Reactor](controlreactor)
+: A reactor where the volume is prescribed by the motion of the reactor's walls.
+
+[Constant Pressure Reactor](constant-pressure-reactor)
+: A reactor where the pressure is held constant.
+
+[Ideal Gas Control Volume Reactor](ideal-gas-reactor)
+: A reactor where the volume is prescribed by the motion of the reactor's walls,
+ specialized for ideal gas mixtures.
+
+[Ideal Gas Constant Pressure Reactor](ideal-gas-constant-pressure-reactor)
+: A reactor where the pressure is held constant, specialized for ideal gas mixtures.
+
+[Control Volume Mole Reactor](mole-reactor)
+: A reactor where the volume is prescribed by the motion of the reactor's walls, with
+ the composition stored in moles.
+
+[Constant Pressure Mole Reactor](constant-pressure-mole-reactor)
+: A reactor where the pressure is held constant and the composition is stored in moles.
+
+[Ideal Gas Control Volume Mole Reactor](ideal-gas-mole-reactor)
+: A reactor where the volume is prescribed by the motion of the reactor's walls,
+ specialized for ideal gas mixtures and with the composition stored in moles.
+
+[Ideal Gas Constant Pressure Mole Reactor](ideal-gas-constant-pressure-mole-reactor)
+: A reactor where the pressure is held constant, specialized for ideal gas mixtures and
+ with the composition stored in moles.
+
+```{seealso}
+In some cases, Cantera's built-in reactor types are insufficient to model a problem.
+In this situation, the {py:class}`ExtensibleReactor` family of classes can be used to
+implement modified governing equations, starting from one of the built-in reactor types
+described above.
+
+The [Guide to Extending Reactor Models](/userguide/extensible-reactor) can help you get
+started with implementing your own customized reactor models.
+```
+
+## Plug Flow Reactors
+
+A *plug flow reactor* (PFR) represents a steady-state flow in a channel. The fluid is
+considered to be homogeneous perpendicular to the flow direction, while the state of the
+gas is allowed to change in the axial direction. However, all diffusion processes are
+neglected.
+
+These assumptions result in a system of equations that is similar to those used to model
+homogeneous reactors, with the time variable replaced by the axial coordinate. Because
+of this mathematical similarity, PFRs are also solved by Cantera's reactor network
+model. However, they can only be simulated alone, and not part of a network containing
+time-dependent reactors.
+
+[Plug Flow Reactor](pfr)
+: A reactor modeling one-dimensional steady-state flow in a channel that may contain
+ catalytically active surfaces where heterogeneous reactions occur.
+
+(sec-reactor-interactions)=
+## Reactor Interactions
+
+Reactors can interact with each other and the surrounding environment in multiple ways.
+Mass can flow from one reactor into another can be incorporated, heat can be
+transferred, and the walls between reactors can move. In addition, reactions can occur
+on surfaces within a reactor. The models used to establish these interconnections are
+described in the following sections:
+
+All of these interactions can vary as a function of time or system state. For example,
+heat transfer can be described as a function of the temperature difference between the
+reactor and the environment, or wall movement can be modeled as depending on the
+pressure difference. Interactions of the reactor with the environment are defined using
+the following models:
+
+[Reservoirs](sec-reservoir)
+: Reservoirs are used to represent constant conditions defining the inlets, outlets, and
+ surroundings of a reactor network.
+
+[Flow Devices](sec-flow-device)
+: Flow devices are used to define mass transfer between two reactors, or between
+ reactors and the surrounding environment as defined by a reservoir.
+
+[Walls](sec-wall)
+: Walls between reactors are used to allow heat transfer between reactors. By moving the
+ walls of the reactor, its volume can be changed and expansion or compression work can
+ be done by or on the reactor.
+
+[Reacting Surfaces](sec-reactor-surface)
+: Reactions may occur on the surfaces of a reactor. These reactions may include both
+ catalytic reactions and reactions resulting in net mass transfer between the surface
+ and the fluid.
+
+```{seealso}
+Cantera comes with a broad variety of well-commented example scripts for reactor
+networks. Please see the [Cantera Examples](/examples/python/reactors/index) for further
+information.
+```
+
+## 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
+:maxdepth: 1
+
+controlreactor
+constant-pressure-reactor
+ideal-gas-reactor
+ideal-gas-constant-pressure-reactor
+mole-reactor
+constant-pressure-mole-reactor
+ideal-gas-mole-reactor
+ideal-gas-constant-pressure-mole-reactor
+
+pfr
+interactions
+```
diff --git a/doc/sphinx/reference/reactors/interactions.md b/doc/sphinx/reference/reactors/interactions.md
new file mode 100644
index 0000000000..9075e6f649
--- /dev/null
+++ b/doc/sphinx/reference/reactors/interactions.md
@@ -0,0 +1,210 @@
+```{py:currentmodule} cantera
+```
+
+# Reactor Interactions
+
+(sec-reservoir)=
+## Reservoirs
+
+A reservoir can be thought of as an infinitely large volume, in which all states are
+predefined and never change from their initial values. Typically, it represents a vessel
+to define temperature and composition of a stream of mass flowing into a reactor, or the
+ambient fluid surrounding the reactor network. In addition, the fluid flow exiting a
+reactor network has to flow into a reservoir. In the latter case, the state of the
+reservoir, except for its pressure, is irrelevant.
+
+:::{admonition} Python Usage
+:class: tip
+
+In Python, a reservoir can be defined using the {py:class}`Reservoir` class.
+:::
+
+(sec-flow-device)=
+## Flow Devices
+
+A *flow device* connects two reactors and allows mass flow from an upstream reactor (1)
+to a downstream reactor (2). Different kinds of flow devices are defined based on how
+the mass flow rate is determined.
+
+A reactor can have multiple inlets and outlets. For the inlets, arbitrary states can be
+defined by setting a reservoir as the upstream reactor. Fluid with the current state of
+the reactor exits the reactor at the outlets.
+
+(sec-valve)=
+### Valves
+
+A *valve* is a flow device with mass flow rate that is a function of the pressure drop
+across it. The mass flow rate is computed as:
+
+$$ \dot m = K_v g(t) f(P_1 - P_2) $$
+
+with $K_v$ being a proportionality constant in kg/s/Pa and $f$ and $g$ being scalar
+functions that can be defined by the user. By default, $g(t) = 1$ and $f(P_1 - P_2) =
+P_1 - P_2$ such that the mass flow rate defaults to:
+
+$$ \dot m = K_v (P_1 - P_2) $$
+
+corresponding to a linear dependence of the mass flow rate on the pressure difference.
+The pressure difference between the upstream (1) and downstream (2) reactor is defined
+as $P_1 - P_2$. The flow is not allowed to reverse and go from the downstream to the
+upstream reactor/reservoir, which means that the flow rate is set to zero if the
+computed value of $\dot m$ is negative, for example if $P_1 < P_2$.
+
+Valves are often used between an upstream reactor and a downstream reactor or reservoir
+to maintain them both at nearly the same pressure. By setting the constant $K_v$ to a
+sufficiently large value, very small pressure differences will result in flow between
+the reactors that counteracts the pressure difference. However, excessively large values
+of $K_v$ may slow down the integrator or cause it to fail.
+
+:::{admonition} Python Usage
+:class: tip
+
+In Python, a valve is implemented by class {py:class}`Valve`; $K_v$ is set using the
+{py:attr}`Valve.valve_coeff` property; $f(P_1 - P_2)$ is set using the
+{py:attr}`Valve.time_function` property; and $g(t)$ is set using the
+{py:attr}`Valve.pressure_function` property.
+:::
+
+(sec-mass-flow-controller)=
+### Mass Flow Controllers
+
+A mass flow controller maintains a specified mass flow rate independent of upstream and
+downstream conditions. The equation used to compute the mass flow rate is
+
+$$ \dot m = m_0 g(t) $$
+
+where $m_0$ is a mass flow coefficient and $g(t)$ is user-specifiable function of time
+that defaults to $g(t) = 1$. Note that if $\dot m < 0$, the mass flow rate will be set
+to zero, since a reversal of the flow direction is not allowed.
+
+Unlike a real mass flow controller, a Cantera mass flow controller object will maintain
+the flow even if the downstream pressure is greater than the upstream pressure. This
+allows simple implementation of loops, in which exhaust gas from a reactor is fed back
+into it through an inlet. But note that this capability should be used with caution,
+since no account is taken of the work required to do this.
+
+:::{admonition} Python Usage
+:class: tip
+
+Mass flow controllers can be implemented in Python using the
+{py:class}`MassFlowController` class. The {py:attr}`MassFlowController.mass_flow_coeff`
+property can be used to set $m_0$ and the {py:attr}`MassFlowController.time_function`
+property can be used to set $g(t)$.
+:::
+
+(sec-pressure-controller)=
+### Pressure Controllers
+
+A pressure controller is designed to be used in conjunction with a primary flow
+controller, typically a mass flow controller. The primary flow controller is installed
+on the inlet of the reactor, and the corresponding pressure controller is installed on
+on outlet of the reactor. The mass flow rate of the pressure controller is equal to that
+of the primary mass flow rate, plus a small correction dependent on the pressure
+difference:
+
+$$ \dot m = \dot m_\t{primary} + K_v f(P_1 - P_2) $$
+
+where $K_v$ is a proportionality constant and $f$ is a function of the pressure drop
+that defaults to $f(P_1 - P_2) = P_1 - P_2$. If $\dot m < 0$, the mass flow rate will be
+set to zero, since a reversal of the flow direction is not allowed.
+
+:::{admonition} Python Usage
+:class: tip
+
+Pressure controllers can be defined in Python using the {py:class}`PressureController`
+class. The primary flow controller can be set using the
+{py:attr}`PressureController.primary` property; $K_v$ can be set using the
+{py:attr}`PressureController.pressure_coeff` property; and $f(P_1 - P_2)$ can be set
+using the {py:attr}`PressureController.pressure_function`.
+:::
+
+(sec-wall)=
+## Walls
+
+In Cantera, a wall separates two reactors or a reactor and a reservoir. A wall has a
+finite area, may conduct or radiate heat between the two reactors on either side, and
+may move like a piston.
+
+Walls are stateless objects in Cantera, meaning that no differential equation is
+integrated to determine any wall property. Since it is the wall, or piston, velocity
+that enters the energy equation, this means that it is the velocity, not the
+acceleration or displacement, that is specified. The wall velocity is computed from
+
+$$ v = K(P_\t{left} - P_\t{right}) + v_0(t) $$
+
+where $K$ is a non-negative constant, and $v_0(t)$ is a specified function of time. The
+velocity is positive if the wall is moving to the right.
+
+The total rate of heat transfer through all walls is:
+
+$$ \dot{Q} = \sum_w f_w \dot{Q}_w $$
+
+where $f_w = \pm 1$ indicates the facing of the wall (-1 for the reactor on the left, +1
+for the reactor on the right). The heat flux $\dot{Q}_w$ through a wall $w$ connecting
+reactors *left* and *right* is computed as:
+
+$$
+\dot{Q}_w = U A (T_\t{left} - T_\t{right})
+ + \epsilon\sigma A (T_\t{left}^4 - T_\t{right}^4) + A q_0(t)
+$$
+
+where $U$ is a user-specified heat transfer coefficient (W/m{sup}`2`-K), $A$ is the wall
+area (m{sup}`2`), $\epsilon$ is the user-specified emissivity, $\sigma$ is the
+Stefan-Boltzmann radiation constant, and $q_0(t)$ is a user-specified, time-dependent
+heat flux (W/m{sup}`2`). This definition is such that positive $q_0(t)$ implies heat
+transfer from the "left" reactor to the "right" reactor. Each of the user-specified
+terms defaults to 0.
+
+:::{admonition} Python Usage
+:class: tip
+
+In Python, walls are defined using the {py:class}`Wall` class.
+:::
+
+(sec-reactor-surface)=
+## Reacting Surfaces
+
+In case of surface reactions, there can be a net generation or destruction of
+homogeneous, gas phase species. The molar rate of production for each homogeneous phase
+species $k$ on surface $w$ is $\dot{s}_{k,w}$ (in kmol/s/m{sup}`2`).
+
+### Mass fraction-based reactors
+
+The total mass production rate for homogeneous phase species $k$ on all surfaces is:
+
+$$ \dot{m}_{k,\t{surf}} = W_k \sum_w A_w \dot{s}_{k,w} $$
+
+where $W_k$ is the molecular weight of species $k$ and $A_w$ is the area of each
+surface. The net mass flux from all reacting surfaces is then:
+
+$$ \dot{m}_\t{surf} = \sum_k \dot{m}_{k,\t{surf}} $$
+
+For each surface species $i$, the rate of change of the site fraction $\theta_{i,w}$ on
+each surface $w$ is integrated with time:
+
+$$ \frac{d\theta_{i,w}}{dt} = \frac{\dot{s}_{i,w} \sigma_i}{\Gamma_w} $$
+
+where $\Gamma_w$ is the total surface site density on surface $w$ and $\sigma_i$ is the
+number of surface sites occupied by a molecule of species $i$ and is sometimes referred
+to within Cantera as the molecule's *size*. The equations for $d\theta_{i,w}/dt$ are
+additional ODEs appended to the state vector for the corresponding reactor.
+
+### Mole-based reactors
+
+In the case of mole based reactors, $\dot{n}_\t{surf}$ is used instead, and is
+calculated as:
+
+$$ \dot{n}_{k,\t{surf}} = A_{w}\sum_{w}\dot{s}_{w, k} $$
+
+and the conservation equation for each surface species $i$ is
+
+$$ \frac{d n_{i,w}}{dt} = \dot{s}_{i,w} A_w $$
+
+These equations for $dn_{i,w}/dt$ are additional ODEs appended to the state
+vector for the corresponding reactor.
+
+:::{admonition} Python Usage
+:class: tip
+
+In Python, reacting surfaces are defined using the {py:class}`ReactorSurface` class.
+:::
diff --git a/doc/sphinx/reference/reactors/mole-reactor.md b/doc/sphinx/reference/reactors/mole-reactor.md
new file mode 100644
index 0000000000..f7ba0ef7cb
--- /dev/null
+++ b/doc/sphinx/reference/reactors/mole-reactor.md
@@ -0,0 +1,56 @@
+```{py:currentmodule} cantera
+```
+
+# Control Volume Mole Reactor
+
+A control volume mole reactor, as implemented by the C++ class {ct}`MoleReactor` and
+available in Python as the {py:class}`MoleReactor` class. It is defined by the state
+variables:
+
+- $U$, the total internal energy of the reactor's contents (in J)
+- $V$, the reactor volume (in m{sup}`3`)
+- $n_k$, the number of moles for each species (in kmol)
+
+Equations 1-3 are the governing equations for a control volume mole reactor.
+
+## Volume Equation
+
+The reactor volume changes as a function of time due to the motion of one or
+more [walls](sec-wall):
+
+$$
+\frac{dV}{dt} = \sum_w f_w A_w v_w(t)
+$$ (molereactor-volume)
+
+where $f_w = \pm 1$ indicates the facing of the wall (whether moving the wall increases
+or decreases the volume of the reactor), $A_w$ is the surface area of the wall, and
+$v_w(t)$ is the velocity of the wall as a function of time.
+
+## Species Equations
+
+The moles of each species in the reactor changes as a result of flow through the
+reactor's [inlets and outlets](sec-flow-device), and production of homogeneous gas phase
+species and reactions on the reactor [surfaces](sec-reactor-surface). The rate at which
+species $k$ is generated through homogeneous phase reactions is $V \dot{\omega}_k$, and
+the total rate at which moles of species $k$ changes is:
+
+$$
+\frac{dn_k}{dt} = V \dot{\omega}_k + \sum_\t{in} \dot{n}_{k, \t{in}}
+ - \sum_\t{out} \dot{n}_{k, \t{out}} + \dot{n}_{k, \t{wall}}
+$$ (molereactor-species)
+
+where the subscripts *in* and *out* refer to the sum of the corresponding property over
+all inlets and outlets respectively. A dot above a variable signifies a time derivative.
+
+## Energy Equation
+
+The equation for the total internal energy is found by writing the first law for an open
+system:
+
+$$
+\frac{dU}{dt} = - p \frac{dV}{dt} + \dot{Q} + \sum_\t{in} \dot{n}_\t{in} \hat{h}_\t{in}
+ - \hat{h} \sum_\t{out} \dot{n}_\t{out}
+$$ (molereactor-energy)
+
+where $\dot{Q}$ is the net rate of heat addition to the system and $\hat{h}$ is the
+molar enthalpy.
diff --git a/doc/sphinx/reference/reactors/pfr.md b/doc/sphinx/reference/reactors/pfr.md
new file mode 100644
index 0000000000..33ae427e2b
--- /dev/null
+++ b/doc/sphinx/reference/reactors/pfr.md
@@ -0,0 +1,127 @@
+```{py:currentmodule} cantera
+```
+
+# Plug Flow Reactor
+
+A plug flow reactor (PFR) represents a one-dimensional steady-state flow in a channel.
+Perpendicular to the flow direction, the gas is assumed to be homogenous.
+In the axial direction $z$, the state of the gas is allowed to change. However, all
+diffusion processes are neglected.
+
+In addition, the interior surface of the reactor may consist of one or more
+catalytically active surfaces where heterogeneous reactions occur.
+
+Plug-flow reactors are often used to simulate emission formation and catalytic
+processes.
+
+A plug flow reactor is defined by the state variables:
+
+- $\rho$, the density of the fluid phase (in kg/m{sup}`3`)
+- $u$, the velocity of the fluid phase (in m/s)
+- $p$, the pressure (in Pa)
+- $T$, the temperature (in K)
+- $Y_k$, the mass fractions for each fluid phase species (dimensionless)
+- $\theta_{i,j}$, the coverage of each surface species $i$ on each surface
+ $j$ (dimensionless)
+
+The reactor geometry is defined by the length $L$, total volume $V$, and the surface
+area $A_{s,j}$ for each surface, where the surface area per unit length is assumed to be
+a constant along the length of the reactor.
+
+The governing equations for a PFR are a system of differential-algebraic equations,
+which depend on the spatial derivatives of some but not all of the state variables. The
+plug flow reactor model in Cantera is implemented by class {ct}`FlowReactor` and
+available in Python as the {py:class}`FlowReactor` class.
+
+## Equation of State
+
+The fluid satisfies the ideal gas law:
+
+$$
+\rho = \frac{p \overline{W}}{R T}
+$$
+
+where $R$ is the universal gas constant and $\overline{W}$ is the mixture molecular
+weight.
+
+## Mass Conservation
+
+The net rate per unit cross sectional area at which fluid phase species are generated by
+reactions on the walls can be defined as
+
+$$ \dot{m}_s = \sum_j \frac{A_{s,j}}{V} \left(\sum_k \dot{s}_{k,j} W_k\right) $$
+
+where $\dot{s}_{k,j}$ is the net rate of production of species $k$ on surface $j$ and
+$W_k$ is the molecular weight of species $k$. The overall mass conservation equation for
+the reactor can then be written as:
+
+$$
+u \frac{d\rho}{dz} + \rho \frac{du}{dz} = \dot{m}_s
+$$ (pfr-mass)
+
+## Momentum Conservation in the Axial Direction
+
+$$
+\rho u \frac{du}{dz} = - u \dot{m}_s - \frac{dp}{dz}
+$$ (pfr-momentum)
+
+## Energy Equation
+
+$$
+\rho u c_p \frac{dT}{dz} =
+ - \sum_k \hat{h}_k \dot{\omega}_k
+ - \sum_j \frac{A_{s,j}}{V} \left(\sum_k \hat{h}_k \dot{s}_{k,j}\right)
+$$ (pfr-energy)
+
+where $c_p$ is the specific heat capacity at constant pressure of the mixture and
+$\hat{h}_k$ is the molar enthalpy of species $k$. Changes in kinetic energy and
+gravitational potential energy are neglected.
+
+## Gas Phase Species Equations
+
+$$
+\rho u \frac{d Y_k}{dz} = - Y_k \dot{m}_s
+ + \dot{\omega}_k W_k + \sum_j \frac{A_{s,j}}{V} \dot{s}_k W_k
+$$ (pfr-species)
+
+## Surface Phase Species Equations
+
+Because the PFR is modeled as steady state, net rate of production for each surface
+species must be zero.
+
+$$
+\dot{s}_{i,j} = 0,\quad i \ne 0
+$$ (pfr-surf-species-i)
+
+To satisfy the constraint that the total surface coverage is 1, the conservation
+equation for the first surface species on each surface is replaced by this constraint:
+
+$$
+\dot{s}_{0,j} = 1 - \sum_{i\ne 0} \theta_{i,j}
+$$ (pfr-surf-species-0)
+
+Without this constraint, the solver could find the trivial, non-physical solution
+$\theta_{i,j} = 0$ for all species, since this also satisfies the steady-state equations
+$\dot{s}_{i,j} = 0$.
+
+## Integrating the PFR Equations
+
+Because diffusion is neglected, downstream parts of the reactor have no influence on
+upstream parts. Therefore, PFRs can be integrated as initial value problems, starting
+from the composition at the inlet. Some care is required to determine initial values for
+the algebraic variables (the surface species coverages) and the time derivatives for the
+differential variables (the other state variables) that are consistent with the
+governing equations.
+
+To do this, we first solve the steady-state problem for each surface, holding the fluid
+phase composition, temperature, and pressure fixed at the inlet values to determine the
+values of $\theta_{i,j}(z=0)$. Then, we construct a linear system comprising the ideal
+gas law, differentiated with respect to $z$ and equations (1), (2), (3), and (4)
+evaluated at $z=0$. This system is then solved to obtain the values of $d\rho/dz$,
+$du/dz$, $dp/dz$, $dT/dz$, and $dY_k/dz$ at $z=0$.
+
+## Examples
+
+- [Partial oxidation of methane over a platinum catalyst](/examples/python/reactors/surf_pfr)
+- [silicon nitride (Si3N4) deposition from ammonia (NH3) and silicon tetrafluoride
+ (SiF4)](/examples/python/surface_chemistry/1D_pfr_surfchem)
diff --git a/doc/sphinx/reference/science/reactors/index.md b/doc/sphinx/reference/science/reactors/index.md
deleted file mode 100644
index fe80fddfd3..0000000000
--- a/doc/sphinx/reference/science/reactors/index.md
+++ /dev/null
@@ -1,67 +0,0 @@
-# Reactors and Reactor Networks
-
-```{caution}
-This page is a work in progress. For more complete documentation of the current Cantera
-release (Cantera 3.0), please see this page.
-```
-
-## Reactors
-
-A Cantera *reactor* represents the simplest form of a chemically reacting system. It
-corresponds to an extensive thermodynamic control volume $V$, in which all state
-variables are homogeneously distributed. The system is generally unsteady -- that is,
-all states are functions of time. In particular, transient state changes due to chemical
-reactions are possible. However, thermodynamic (but not chemical) equilibrium is assumed
-to be present throughout the reactor at all instants of time.
-
-Reactors can interact with the surrounding environment in multiple ways:
-
-- Expansion/compression work: By moving the walls of the reactor, its volume can be
- changed and expansion or compression work can be done by or on the reactor.
-
-- Heat transfer: An arbitrary heat transfer rate can be defined to cross the boundaries
- of the reactor.
-
-- Mass transfer: The reactor can have multiple inlets and outlets. For the inlets,
- arbitrary states can be defined. Fluid with the current state of the reactor exits the
- reactor at the outlets.
-
-- Surface interaction: One or multiple walls can influence the chemical reactions in the
- reactor. This is not just restricted to catalytic reactions, but mass transfer between
- the surface and the fluid can also be modeled.
-
-All of these interactions do not have to be constant, but can vary as a function of time
-or state. For example, heat transfer can be described as a function of the temperature
-difference between the reactor and the environment, or the wall movement can be modeled
-depending on the pressure difference. Interactions of the reactor with the environment
-are defined on one or more walls, inlets, and outlets.
-
-In addition to single reactors, Cantera is also able to interconnect reactors into a
-*reactor network*. Each reactor in a network may be connected so that the contents of one
-reactor flow into another. Reactors may also be in contact with one another or the
-environment via walls that conduct heat or move to do work.
-
-## Reactor Types and Governing Equations
-
-All reactor types are modelled using combinations of Cantera's governing equations of
-state. The specific governing equations defining Cantera's supported reactor models are
-derived and described below.
-
-
-````{grid} 2
-:gutter: 3
-
-```{grid-item-card} Control Volume Reactor
-:link: controlreactor
-:link-type: doc
-
-A reactor where the volume is prescribed by the motion of the reactor's walls. The state variables are the volume, mass, total internal energy, and species mass fractions.
-```
-
-````
-
-```{toctree}
-:hidden:
-
-controlreactor
-```
\ No newline at end of file
diff --git a/doc/sphinx/reference/thermo/index.md b/doc/sphinx/reference/thermo/index.md
new file mode 100644
index 0000000000..5861cb9e6a
--- /dev/null
+++ b/doc/sphinx/reference/thermo/index.md
@@ -0,0 +1,36 @@
+# Thermodynamic Properties
+
+In this section, we describe how Cantera uses species and phase information to calculate
+thermodynamic properties.
+
+Thermodynamic properties typically depend on information at both the species and phase
+levels. Generally, the species thermodynamic model and accompanying coefficient data
+specifies how the reference enthalpy and entropy values for each species are calculated
+as a function of temperature. The phase model then describes how the species interact
+with one another to determine phase properties and species specific properties for a
+given thermodynamic state. This includes both the mechanical equation of state
+($p$-$\hat{v}$-$T$ relationship) as well as how species-specific properties, such as
+internal energy, entropy, and others, depend on the state variables.
+
+The user must specify the thermodynamic models and provide input data to be used for
+both levels, and these selections must be compatible with one another. For instance: one
+cannot pair certain non-ideal species thermodynamic models with an ideal phase model.
+
+The following sections describe the species and phase thermodynamic models available
+in Cantera.
+
+[](species-thermo)
+: The models and equations that Cantera uses to calculate species thermodynamic
+ properties, such as the NASA 7-parameter polynomial form.
+
+[](phase-thermo)
+: The theory behind some of Cantera's phase models, such as the ideal gas law.
+
+```{toctree}
+:hidden:
+:maxdepth: 1
+:caption: Thermodynamic Models
+
+species-thermo
+phase-thermo
+```
diff --git a/doc/sphinx/reference/thermo/phase-thermo.md b/doc/sphinx/reference/thermo/phase-thermo.md
new file mode 100644
index 0000000000..85d50e3f61
--- /dev/null
+++ b/doc/sphinx/reference/thermo/phase-thermo.md
@@ -0,0 +1,164 @@
+# Phase Thermodynamic Models
+
+On this page, we list the phase thermodynamic models implemented in Cantera, with
+links to the documentation for their YAML input parameters and the documentation for
+the C++ classes which implement these models. This API documentation may also provide
+references or a mathematical description of the model.
+
+
+## Models for Gaseous Mixtures
+
+(sec-ideal-gas-phase)=
+Ideal Gas Mixture
+: A mixture which follows the ideal gas law. Defined in the YAML format by specifying
+ [`ideal-gas`](sec-yaml-ideal-gas) in the `thermo` field of the phase definition.
+ Implemented by class {ct}`IdealGasPhase`.
+
+(sec-Redlich-Kwong-phase)=
+Redlich-Kwong Real Gas Mixture
+: A multi-species real gas following the Redlich-Kwong equation of state. Defined in the
+ YAML format by specifying [`Redlich-Kwong`](sec-yaml-Redlich-Kwong) in the `thermo`
+ field of the phase definition. Implemented by class {ct}`RedlichKwongMFTP`.
+
+(sec-Peng-Robinson-phase)=
+Peng-Robinson Real Gas Mixture
+: A multi-species real gas following the Peng-Robinson equation of state. Defined in the
+ YAML format by specifying [`Peng-Robinson`](sec-yaml-Peng-Robinson) in the `thermo`
+ field of the phase definition. Implemented by class {ct}`PengRobinson`.
+
+(sec-plasma-phase)=
+Plasma
+: A phase that extends the ideal gas model to handle plasma properties such as the
+ electron energy distribution and electron temperature with different models. Defined
+ in the YAML format by specifying [`plasma`](sec-yaml-plasma) in the `thermo` field of
+ the phase definition. Implemented by class {ct}`PlasmaPhase`.
+
+
+## Models for Surfaces and Interfaces
+
+(sec-ideal-surface-phase)=
+Ideal Surface
+: An interface between two bulk phases where the species behave as an ideal solution and
+ the composition is described by the coverage of each species on the surface. Defined
+ in the YAML format by specifying [`ideal-surface`](sec-yaml-ideal-surface) in the
+ `thermo` field of the phase definition. Implemented by class {ct}`SurfPhase`.
+
+(sec-coverage-dependent-surface-phase)=
+Surface Phase with Coverage-Dependent Thermo
+: A coverage-dependent surface phase. That is, a surface phase where the enthalpy,
+ entropy, and heat capacity of each species may depend on its coverage and the coverage
+ of other species in the phase. Defined in the YAML format by specifying
+ [`coverage-dependent-surface`](sec-yaml-coverage-dependent-surface) in the `thermo`
+ field of the phase definition. Implemented by class {ct}`CoverageDependentSurfPhase`.
+
+(sec-edge-phase)=
+Edge
+: A one-dimensional edge between two surfaces, which defines a triple phase boundary.
+ Defined in the YAML format by specifying [`edge`](sec-yaml-edge) in the `thermo` field
+ of the phase definition. Implemented by class {ct}`EdgePhase`.
+
+
+## Single-species Phase Models
+
+(sec-fixed-stoichiometry-phase)=
+Stoichiometric Substance
+: A *stoichiometric substance* is one that is modeled as having a precise, fixed
+ composition, given by the composition of the one species present. Defined in the YAML
+ format by specifying [`fixed-stoichiometry`](sec-yaml-fixed-stoichiometry) in the
+ `thermo` field of the phase definition. Implemented by class {ct}`StoichSubstance`.
+
+(sec-electron-cloud-phase)=
+Electron Cloud
+: A phase representing an electron cloud, such as conduction electrons in a metal.
+ Defined in the YAML format by specifying [`electron-cloud`](sec-yaml-electron-cloud)
+ in the `thermo` field of the phase definition. Implemented by class {ct}`MetalPhase`.
+
+(sec-pure-fluid-phase)=
+Pure Fluid Phases
+: A phase representing a pure fluid equation of state for one of several pure substances
+ including liquid, vapor, two-phase, and supercritical regions. Defined in the YAML
+ format by specifying [`pure-fluid`](sec-yaml-pure-fluid) in the `thermo` field of the
+ phase definition. Implemented by class {ct}`PureFluidPhase`.
+
+(sec-liquid-water-IAPWS95-phase)=
+Liquid Water using the IAPWS95 Equation of State
+: An implementation of the IAPWS95 equation of state for water {cite:p}`wagner2002`, for
+ the liquid region only. Defined in the YAML format by specifying
+ [`liquid-water-IAPWS95`](sec-yaml-liquid-water-IAPWS95) in the `thermo` field of the
+ phase definition. Implemented by class {ct}`WaterSSTP`.
+
+
+## Ideal Solid and Liquid Solutions
+
+(sec-ideal-molal-solution-phase)=
+Ideal Molal Solution
+: An ideal solution based on the mixing-rule assumption that all molality-based activity
+ coefficients are equal to one. Defined in the YAML format by specifying
+ [`ideal-molal-solution`](sec-yaml-ideal-molal-solution) in the `thermo` field of the
+ phase definition. Implemented by class {ct}`IdealMolalSoln`.
+
+(sec-ideal-condensed-phase)=
+Ideal Condensed Phase
+: An ideal liquid or solid solution based on the mixing-rule assumption that all molar
+ concentration-based activity coefficients are equal to one. Defined in the YAML format
+ by specifying [`ideal-condensed`](sec-yaml-ideal-condensed) in the `thermo` field of
+ the phase definition. Implemented by class {ct}`IdealSolidSolnPhase`.
+
+(sec-ideal-solution-VPSS-phase)=
+Ideal Condensed Phase with VPSS Species
+: An ideal solution model using variable pressure standard state methods. This allows
+ the standard state molar volume of species to be specified as a function of
+ temperature. Defined in the YAML format by specifying
+ [`ideal-solution-VPSS`](sec-yaml-ideal-solution-VPSS) in the `thermo` field of the
+ phase definition. Implemented by class {ct}`IdealSolnGasVPSS`.
+
+(sec-lattice-phase)=
+Lattice Phase
+: A simple thermodynamic model for a bulk phase, assuming an incompressible lattice of
+ solid atoms. Defined in the YAML format by specifying [`lattice`](sec-yaml-lattice) in
+ the `thermo` field of the phase definition. Implemented by class {ct}`LatticePhase`.
+
+(sec-compound-lattice-phase)=
+Compound Lattice Phase
+: A phase that is comprised of a fixed additive combination of other lattice phases.
+ Defined in the YAML format by specifying [`compound-lattice`](sec-yaml-compound-lattice)
+ in the `thermo` field of the phase definition. Implemented by class
+ {ct}`LatticeSolidPhase`.
+
+
+## Non-ideal Solid and Liquid Solutions
+
+(sec-binary-solution-tabulated-phase)=
+Binary Solution with Tabulated Enthalpy and Entropy
+: A phase representing a non-ideal binary solution where the excess enthalpy and entropy
+ are interpolated between tabulated values as a function of mole fraction. Defined in
+ the YAML format by specifying
+ [`binary-solution-tabulated`](sec-yaml-binary-solution-tabulated) in the `thermo`
+ field of the phase definition. Implemented by class
+ {ct}`BinarySolutionTabulatedThermo`.
+
+(sec-Debye-Huckel-phase)=
+Debye-Huckel Solution
+: A dilute liquid electrolyte which obeys the Debye-Hückel formulation for nonideality.
+ Defined in the YAML format by specifying [`Debye-Huckel`](sec-yaml-Debye-Huckel) in
+ the `thermo` field of the phase definition. Implemented by class {ct}`DebyeHuckel`.
+
+(sec-HMW-electrolyte-phase)=
+Harvie--Møller--Weare electrolyte
+: A dilute or concentrated liquid electrolyte phase that obeys the Pitzer formulation
+ for nonideality. Defined in the YAML format by specifying
+ [`HMW-electrolyte`](sec-yaml-HMW-electrolyte) in the `thermo` field of the phase
+ definition. Implemented by class {ct}`HMWSoln`.
+
+(sec-Margules-phase)=
+Margules Solution
+: A condensed phase employing the Margules approximation for the excess Gibbs free
+ energy. Defined in the YAML format by specifying [`Margules`](sec-yaml-Margules) in
+ the `thermo` field of the phase definition. Implemented by class {ct}`MargulesVPSSTP`.
+
+(sec-Redlich-Kister-phase)=
+Redlich-Kister Solution
+: A phase employing the Redlich-Kister approximation for the excess Gibbs free energy.
+ Defined in the YAML format by specifying [`Redlich-Kister`](sec-yaml-Redlich-Kister)
+ in the `thermo` field of the phase definition. Implemented by class
+ {ct}`RedlichKisterVPSSTP`.
diff --git a/doc/sphinx/reference/thermo/species-thermo.md b/doc/sphinx/reference/thermo/species-thermo.md
new file mode 100644
index 0000000000..4752134d47
--- /dev/null
+++ b/doc/sphinx/reference/thermo/species-thermo.md
@@ -0,0 +1,341 @@
+# Species Thermodynamic Models
+
+The phase models discussed in the [](phase-thermo) section implement specific models for
+the thermodynamic properties appropriate for the type of phase or interface they
+represent. Although each one may use different expressions to compute the properties,
+they all require thermodynamic property information for the individual species.
+
+Generally, the phase models require a parameterization of the _standard state_ heat
+capacity, enthalpy, and entropy for each species at a fixed pressure $p^\circ$ as a
+function of $T$. In addition, phase models may require information describing how each
+species affects the equation of state, either in terms of the species standard molar
+volume, or through parameters that are specific to the equation of state.
+
+(sec-standard-state-species-thermo)=
+## Models for Species Standard State Enthalpy, Entropy, and Heat Capacity
+
+Many of Cantera's phase thermodynamic models are formulated to make use of the
+[standard state](https://goldbook.iupac.org/terms/view/S05925) thermodynamic properties
+for individual species, defined at a standard pressure $p^\circ$ and for the composition
+specified by the phase model. For example, this could include a pure gas in the case of
+the ideal gas model, or an ion at infinite dilution in water in the case of aqueous
+solutions. The value of $p^\circ$ is not fixed by Cantera, and may vary among different
+sources of thermodynamic data.
+
+```{caution}
+In some parts of the Cantera documentation, properties calculated at the user-specified
+standard pressure $p^\circ$ are referred to as *reference-state* thermodynamic
+properties, as they represent a well-known reference state and properties for all other
+states are calculated according to their departure from this known reference condition.
+In these same parts of the documentation, the term _standard-state properties_ refers to
+properties calculated using the composition defining the standard state but at the
+mixture's current pressure.
+
+This nomenclature is fairly unique to Cantera, based on the desire to distinguish these
+different steps in the calculation of the full thermodynamic properties, and is not
+often seen in other descriptions of solution thermodynamics.
+```
+
+The necessary properties are:
+
+1. $\hat{c}^\circ_p(T)$: the standard-state molar heat capacity at constant pressure as
+ a function of temperature;
+2. $\hat{h}^\circ(T)$, the standard-state molar enthalpy as a function of temperature;
+3. $\hat{s}^\circ(T)$: the standard-state molar entropy as a function of temperature.
+
+While each of these functions is implemented explicitly in Cantera for computational
+efficiency, $\hat{h}^\circ(T)$ and $\hat{s}^\circ(T)$ can be expressed in terms of
+$\hat{c}^\circ_p(T)$ using the relations
+
+$$ \hat{h}^\circ(T) = \hat{h}^\circ(T_\t{ref}) +
+ \int_{T_\t{ref}}^T \hat{c}^\circ_p(T) \; dT $$
+
+and
+
+$$ \hat{s}^\circ(T) = \hat{s}^\circ(T_\t{ref}) +
+ \int_{T_\t{ref}}^T \frac{\hat{c}^\circ_p(T)}{T} \; dT $$
+
+respectively. This means that a parameterization of $\hat{c}_p^\circ(T)$ plus the
+constants $\hat{h}^\circ(T_\t{ref})$ and $\hat{s}^\circ(T_\t{ref})$ at a reference
+temperature $T_\t{ref}$ is sufficient to define the standard state properties for a
+species.
+
+The models described in this section can be used to provide standard state thermodynamic
+data for each species in a phase. They are implemented by classes deriving from
+{ct}`SpeciesThermoInterpType`.
+
+```{note}
+There is no requirement that all species in a phase use the same parameterization; each
+species can use the one most appropriate to represent how the heat capacity depends on
+temperature.
+```
+
+(sec-thermo-nasa7)=
+### The NASA 7-Coefficient Polynomial Parameterization
+
+The NASA 7-coefficient polynomial parameterization is used to compute the species
+standard-state thermodynamic properties $\hat{c}^\circ_p(T)$, $\hat{h}^\circ(T)$, and
+$\hat{s}^\circ(T)$.
+
+The NASA parameterization represents $\hat{c}^\circ_p(T)$ with a fourth-order
+polynomial:
+
+$$
+\frac{\hat{c}_p^\circ(T)}{R} &= a_0 + a_1 T + a_2 T^2 + a_3 T^3 + a_4 T^4
+
+\frac{\hat{h}^\circ (T)}{R T} &= a_0 + \frac{a_1}{2} T + \frac{a_2}{3} T^2 +
+ \frac{a_3}{4} T^3 + \frac{a_4}{5} T^4 + \frac{a_5}{T}
+
+\frac{\hat{s}^\circ(T)}{R} &= a_0 \ln T + a_1 T + \frac{a_2}{2} T^2 +
+ \frac{a_3}{3} T^3 + \frac{a_4}{4} T^4 + a_6
+$$
+
+This model is implemented by the C++ classes {ct}`NasaPoly1` and {ct}`NasaPoly2`.
+
+:::{note}
+This is sometimes referred to as the *NASA7 polynomial* within Cantera. It was used in
+the original NASA equilibrium program and in Chemkin, which uses 7 coefficients in each
+of two temperature regions. It is not compatible with the form used in more recent
+versions of the NASA equilibrium program, which uses 9 coefficients for each temperature
+region.
+:::
+
+:::{admonition} YAML Usage
+:class: tip
+A NASA7 parameterization can be defined in the YAML format by specifying
+[`NASA7`](sec-yaml-nasa7) as the `model` in the species `thermo` field.
+:::
+
+(sec-thermo-nasa9)=
+### The NASA 9-Coefficient Polynomial Parameterization
+
+The NASA 9-coefficient polynomial parameterization {cite:p}`mcbride2002`, often called
+*NASA9* within Cantera for short, is an extension of the NASA 7-coefficient polynomial
+parameterization that includes two additional terms in each temperature region and
+supports an arbitrary number of temperature regions.
+
+The NASA9 parameterization represents the species thermodynamic properties with the
+following equations:
+
+$$
+\frac{\hat{c}_p^\circ(T)}{R} &= a_0 T^{-2} + a_1 T^{-1} + a_2 + a_3 T
+ + a_4 T^2 + a_5 T^3 + a_6 T^4
+
+\frac{\hat{h}^\circ(T)}{R T} &= - a_0 T^{-2} + a_1 \frac{\ln T}{T} + a_2
+ + \frac{a_3}{2} T + \frac{a_4}{3} T^2 + \frac{a_5}{4} T^3 +
+ \frac{a_6}{5} T^4 + \frac{a_7}{T}
+
+\frac{\hat{s}^\circ(T)}{R} &= - \frac{a_0}{2} T^{-2} - a_1 T^{-1} + a_2 \ln T
+ + a_3 T + \frac{a_4}{2} T^2 + \frac{a_5}{3} T^3 + \frac{a_6}{4} T^4 + a_8
+$$
+
+A common source for species data in the NASA9 format is the
+[NASA ThermoBuild](/userguide/thermobuild) tool. This model is implemented by the C++
+classes {ct}`Nasa9Poly1` and {ct}`Nasa9PolyMultiTempRegion`.
+
+:::{admonition} YAML Usage
+:class: tip
+A NASA-9 parameterization can be defined in the YAML format by specifying
+[`NASA9`](sec-yaml-nasa9) as the `model` in the species `thermo` field.
+:::
+
+(sec-thermo-shomate)=
+### The Shomate Parameterization
+
+The Shomate parameterization is:
+
+$$
+\hat{c}_p^\circ(T) &= A + Bt + Ct^2 + Dt^3 + \frac{E}{t^2}
+
+\hat{h}^\circ(T) &= At + \frac{Bt^2}{2} + \frac{Ct^3}{3} + \frac{Dt^4}{4} -
+ \frac{E}{t} + F
+
+\hat{s}^\circ(T) &= A \ln t + B t + \frac{Ct^2}{2} + \frac{Dt^3}{3} - \frac{E}{2t^2} + G
+$$
+
+where $t = T / 1000\textrm{ K}$. It requires 7 coefficients $A$, $B$, $C$, $D$, $E$,
+$F$, and $G$. This parameterization is used to represent standard-state properties in
+the [NIST Chemistry WebBook](http://webbook.nist.gov/chemistry). The values of the
+coefficients $A$ through $G$ should be entered precisely as shown there, with no units
+attached. Unit conversions to SI is handled internally. This model is implemented by the
+C++ classes {ct}`ShomatePoly` and {ct}`ShomatePoly2`.
+
+:::{admonition} YAML Usage
+:class: tip
+A Shomate parameterization can be defined in the YAML format by specifying
+[`Shomate`](sec-yaml-shomate) as the `model` in the species `thermo` field.
+:::
+
+(sec-thermo-const-cp)=
+### Constant Heat Capacity
+
+In some cases, species properties may only be required at a single temperature or over a
+narrow temperature range. In such cases, the heat capacity can be approximated as
+constant, and simple expressions can be used for the thermodynamic properties:
+
+$$
+\hat{c}_p^\circ(T) &= \hat{c}_p^\circ(T_\t{ref})
+
+\hat{h}^\circ(T) &= \hat{h}^\circ\left(T_\t{ref}\right) + \hat{c}_p^\circ \left(T-T_\t{ref}\right)
+
+\hat{s}^\circ(T) &= \hat{s}^\circ(T_\t{ref}) + \hat{c}_p^\circ \ln{\left(\frac{T}{T_\t{ref}}\right)}
+$$
+
+The parameterization uses four constants: $T_\t{ref}$,
+$\hat{c}_p^\circ(T_\t{ref})$, $\hat{h}^\circ(T_\t{ref})$, and
+$\hat{s}^\circ(T)$. The default value of $T_\t{ref}$ is 298.15 K; the default value
+for the other parameters is 0.0. This model is implemented by the C++ class
+{ct}`ConstCpPoly`.
+
+:::{admonition} YAML Usage
+:class: tip
+A constant heat capacity parameterization can be defined in the YAML format by
+specifying [`constant-cp`](sec-yaml-constcp) as the `model` in the species `thermo`
+field.
+:::
+
+(sec-thermo-piecewise-Gibbs)=
+### Piecewise Gibbs Free Energy
+
+This parameterization uses a piecewise interpolation of the Gibbs free energy between specified values, with a constant heat capacity between each pair of
+interpolation points. This model is implemented by the C++ class {ct}`Mu0Poly`.
+
+:::{admonition} YAML Usage
+:class: tip
+A piecewise Gibbs parameterization can be defined in the YAML format by specifying
+[`piecewise-Gibbs`](sec-yaml-constcp) as the `model` in the species `thermo` field.
+:::
+
+## Models for Species Contributions to the Equation of State
+
+Besides enthalpy, entropy, and heat capacity data, some phase models require additional
+parameters that describe how each species affects the mechanical equation of state.
+These models are used in combination with one of the above parameterizations for the
+standard state enthalpy, entropy, and heat capacity. This applies to most models that
+represent condensed phases.
+
+### Constant Volume Model
+
+This species equation of state model can be used for phase models that require standard
+state density information for each species, and simply specifies a constant specific
+volume, density or molar density for the species. Implemented by the C++ class
+{ct}`PDSS_ConstVol`.
+
+:::{admonition} YAML Usage
+:class: tip
+A constant volume parameterization can be defined in the YAML format by specifying
+[`constant-volume`](sec-yaml-eos-constant-volume) as the `model` in the species
+`equation-of-state` field.
+:::
+
+### Temperature-dependent Density
+
+This species equation of state model can be used for phase models that require the
+standard state molar volume $V_{m,k}^\circ$ for each species $k$. The standard state
+density $\rho^\circ_k$ is described using a third-order polynomial:
+
+$$ \rho^\circ_k(T) = \frac{W_k}{V_{m,k}^\circ(T)} = a_0 + a_1 T + a_2 T^2 + a_3 T^3 $$
+
+where $W_k$ is the molecular weight of the species. This model is implemented by class
+{ct}`PDSS_SSVol`.
+
+:::{admonition} YAML Usage
+:class: tip
+A temperature-dependent density parameterization can be defined in the YAML format by
+specifying [`density-temperature-polynomial`](sec-yaml-eos-density-temperature-polynomial)
+as the `model` in the species `equation-of-state` field.
+:::
+
+### Temperature-dependent Molar Volume
+
+This species equation of state model can be used for phase models that require the
+standard state molar volume $V_{m,k}^\circ$ for each species $k$. The standard state
+molar volume is described using a third-order polynomial:
+
+$$ V_{m,k}^\circ(T) = a_0 + a_1 T + a_2 T^2 + a_3 T^3 $$
+
+This model is implemented by class {ct}`PDSS_SSVol`.
+
+:::{admonition} YAML Usage
+:class: tip
+A temperature-dependent molar volume parameterization can be defined in the YAML format
+by specifying [`molar-volume-temperature-polynomial`](sec-yaml-eos-molar-volume-temperature-polynomial)
+as the `model` in the species `equation-of-state` field.
+:::
+
+### Peng-Robinson Equation of State
+
+This species equation of state model is parameterized using the coefficients $a$, $b$,
+and $\omega$ for a pure species that follows the Peng-Robinson equation of state:
+
+$$ P = \frac{RT}{V_m - b} - \frac{a\alpha}{V_m^2 + 2bV_m - b^2} $$
+
+where $V_m$ is the molar volume,
+
+$$ \alpha = \left[ 1 + \kappa \left(1 - \sqrt{T_r}\right) \right]^2 $$
+
+$$ \kappa =
+\begin{cases}
+0.37464 + 1.54226\omega - 0.26992\omega^2, & \omega \le 0.491 \\
+0.379642 + 1.487503\omega - 0.164423\omega^2 + 0.016667\omega^3 , & \omega > 0.491
+\end{cases}
+$$
+
+$T_r = T / T_c$ is the reduced temperature, and $T_c$ is the critical temperature.
+
+These pure-species properties are combined in the multi-species Peng-Robinson phase
+model, implemented by class {ct}`PengRobinson`.
+
+:::{admonition} YAML Usage
+:class: tip
+A Peng-Robinson parameterization for a species can be defined in the YAML format by
+specifying [`Peng-Robinson`](sec-yaml-eos-peng-robinson) as the `model` in the species
+`equation-of-state` field.
+:::
+
+### Redlich-Kwong Equation of State
+
+This species equation of state model is parameterized using the coefficients $a$ and $b$
+for a pure species that follows the Redlich-Kwong equation of state:
+
+$$ P = \frac{RT}{V_m-b} - \frac{a}{\sqrt{T} V_m (V_m + b) } $$
+
+where $V_m$ is the molar volume.
+
+These pure-species properties are combined in the multi-species Redlich-Kwong phase
+model, implemented by class {ct}`RedlichKwongMFTP`.
+
+:::{admonition} YAML Usage
+:class: tip
+A Redlich-Kwong parameterization for a species can be defined in the YAML format by
+specifying [`Redlich-Kwong`](sec-yaml-eos-redlich-kwong) as the `model` in the species
+`equation-of-state` field.
+:::
+
+## Other Species Equation of State Models
+
+The following models provide complete standard state parameterizations for a species,
+including molar volume, enthalpy, entropy, and heat capacity.
+
+### Helgeson-Kirkham-Flowers-Tanger Model
+
+The Helgeson-Kirkham-Flowers-Tanger model for aqueous species. Implemented by the C++
+class {ct}`PDSS_HKFT`.
+
+:::{admonition} YAML Usage
+:class: tip
+An HKFT parameterization can be defined in the YAML format by specifying
+[`HKFT`](sec-yaml-eos-hkft) as the `model` in the species `equation-of-state` field.
+:::
+
+### Liquid Water
+
+A model for liquid water that uses the IAPWS95 formulation {cite:p}`wagner2002`.
+Implemented by the C++ class {ct}`PDSS_Water`.
+
+:::{admonition} YAML Usage
+:class: tip
+A liquid water parameterization can be defined in the YAML format by specifying
+[`liquid-water-IAPWS95`](sec-yaml-eos-liquid-water-iapws95) as the `model` in the
+species `equation-of-state` field.
+:::
diff --git a/doc/sphinx/reference/transport/index.md b/doc/sphinx/reference/transport/index.md
new file mode 100644
index 0000000000..d8f2d1c06b
--- /dev/null
+++ b/doc/sphinx/reference/transport/index.md
@@ -0,0 +1,70 @@
+# Transport Properties
+
+Here, we describe how Cantera uses species and phase information to calculate transport
+properties and rates. Similar to Cantera's approach to
+[thermodynamic properties](/reference/thermo/index), transport property calculations in
+Cantera depend on information at both the species and phase levels. The user must
+specify transport models for both levels, and these selections must be compatible with
+one another.
+
+- The user must specify a transport model for each species and provide inputs that
+ inform how species properties are calculated. For example, the user provides inputs
+ that allow Cantera to calculate species collision integrals based on species-specific
+ Lennard-Jones parameters.
+- The user also selects a phase transport model. This model describes how the species
+ interact with one another to determine properties such as viscosity, thermal
+ conductivity, and diffusion coefficients for a given thermodynamic state.
+
+## Species Transport Parameters
+
+Transport property models in general require parameters that express the effect of each
+species on the transport properties of the phase. Currently, most transport models
+available in Cantera are applicable to gaseous phases.
+
+```{admonition} YAML Usage
+:class: tip
+Gas transport properties can be defined in the YAML format using the
+[`transport`](sec-yaml-species-transport) field of a `species` entry.
+```
+
+(sec-phase-transport-models)=
+## Phase Transport Models
+
+Multicomponent
+: A multicomponent transport model for ideal gases, based on the model described by
+ Dixon-Lewis {cite:t}`dixon-lewis1968`; See also Kee et al. {cite:t}`kee2017`. The
+ multicomponent transport model can be specified in the YAML format by setting the
+ [`transport`](sec-yaml-phase-transport) field of the phase entry to `multicomponent`.
+ Implemented by class {ct}`MultiTransport`.
+
+Mixture-averaged
+: A mixture-averaged transport model for ideal gases, as described in Kee et al.
+ {cite:t}`kee2017`. The mixture-averaged transport model can be specified in the YAML
+ format by setting the [`transport`](sec-yaml-phase-transport) field of the phase entry
+ to `mixture-averaged`. Implemented by class {ct}`MixTransport`.
+
+High-pressure Gas
+: A model for high-pressure gas transport properties based on a method of corresponding
+ states {cite:p}`takahashi1975,poling2001`. The high-pressure gas transport model can
+ be specified in the YAML format by setting the [`transport`](sec-yaml-phase-transport)
+ field of the phase entry to `high-pressure`. Implemented by
+ class {ct}`HighPressureGasTransport`.
+
+Ionized Gas
+: A model implementing the Stockmayer-(n,6,4) model for transport of ions in a gas. The
+ ionized gas transport model can be specified in the YAML format by setting the
+ [`transport`](sec-yaml-phase-transport) field of the phase entry to `ionized-gas`.
+ Implemented by class {ct}`IonGasTransport`.
+
+Unity Lewis Number
+: A transport model for ideal gases where diffusion coefficients for all species are set
+ so that the [Lewis number](https://en.wikipedia.org/wiki/Lewis_number) is 1. The unity
+ Lewis number transport model can be specified in the YAML format by setting the
+ [`transport`](sec-yaml-phase-transport) field of the phase entry to
+ `unity-Lewis-number`. Implemented by class {ct}`UnityLewisTransport`.
+
+Water
+: A transport model for pure water applicable in both liquid and vapor phases. The water
+ transport model can be specified in the YAML format by setting the
+ [`transport`](sec-yaml-phase-transport) field of the phase entry to `water`.
+ Implemented by class {ct}`WaterTransport`.
diff --git a/doc/sphinx/userguide/creating-mechanisms.md b/doc/sphinx/userguide/creating-mechanisms.md
index d66c199589..95b4c6932f 100644
--- a/doc/sphinx/userguide/creating-mechanisms.md
+++ b/doc/sphinx/userguide/creating-mechanisms.md
@@ -422,8 +422,12 @@ case the species is composed of nothing, and represents an empty surface site. T
also be done to represent vacancies in solids. A charged vacancy can be defined to be
composed solely of electrons.
-The number of atoms of an element must be non-negative, except for the special "element"
-`E` that represents an electron.
+The special pseudo-element `E` is used in representing charged species, where it
+specifies the net number of electrons compared to the number needed to form a neutral
+species. That is, negatively charged ions will have `E` > 0, while positively charged
+ions will have `E` \< 0.
+
+The number of atoms of an element must be non-negative, except for electrons.
Examples:
@@ -558,10 +562,10 @@ it is usually best not to specify units for $A$, in which case they will be comp
taking all of these factors into account.
```{note}
-If $b \ne 0$, then the term $T^b$ should have units of $\mathrm{K}^b$, which would
+If $b \ne 0$, then the term $T^b$ should have units of $\t{K}^b$, which would
change the units of $A$. This is not done, however, so the units associated with $A$
are really the units for $k_f$. One way to formally express this is to replace $T^b$
-by the non-dimensional quantity $[T/(1\;\mathrm{K})]^b$.
+by the non-dimensional quantity $[T/(1\;\t{K})]^b$.
```
The key `E` is used to specify $E_a$.
@@ -607,13 +611,13 @@ negative-A: true
Explicit reaction orders different from the stoichiometric coefficients are sometimes
used for non-elementary reactions. For example, consider the global reaction:
-$$ \mathrm{C_8H_{18} + 12.5 O_2 \rightarrow 8 CO_2 + 9 H_2O} $$
+$$ \t{C_8H_{18} + 12.5 O_2 \rightarrow 8 CO_2 + 9 H_2O} $$
the forward rate constant might be given as {cite:p}`westbrook1981`:
$$
-k_f = 4.6 \times 10^{11} [\mathrm{C_8H_{18}}]^{0.25} [\mathrm{O_2}]^{1.5}
- \exp\left(\frac{30.0\,\mathrm{kcal/mol}}{RT}\right)
+k_f = 4.6 \times 10^{11} [\t{C_8H_{18}}]^{0.25} [\t{O_2}]^{1.5}
+ \exp\left(\frac{30.0\,\t{kcal/mol}}{RT}\right)
$$
This reaction could be defined as:
@@ -659,9 +663,12 @@ case, the `nonreactant-orders` field must be added to the reaction entry:
(sec-yaml-guide-elements)=
## Elements
-Cantera provides built-in definitions for the chemical elements, including values for
-their atomic weights taken from IUPAC / CIAAW. These elements can be used by specifying
-the corresponding atomic symbols when specifying the composition of species.
+In Cantera, an *element* may refer to a chemical element or an isotope. Cantera provides
+built-in definitions for the
+chemical elements, including values for their atomic weights taken from
+[IUPAC / CIAAW](http://www.ciaaw.org/atomic-weights.htm). These elements can be used by
+specifying the corresponding atomic symbols when specifying the composition of species.
+Explicit element definitions are usually only needed for isotopes.
In order to give a name to a particular isotope or a virtual element representing a
surface site, a custom `element` entry can be used. The default location for `element`
diff --git a/doc/sphinx/userguide/extensible-reactor.md b/doc/sphinx/userguide/extensible-reactor.md
new file mode 100644
index 0000000000..e17ea82d22
--- /dev/null
+++ b/doc/sphinx/userguide/extensible-reactor.md
@@ -0,0 +1,163 @@
+---
+file_format: mystnb
+kernelspec:
+ name: python3
+---
+
+```{py:currentmodule} cantera
+```
+
+# Extending Reactor Models with Custom Equations
+
+In some cases, Cantera's built-in reactor types are insufficient to model a problem.
+In this situation, the {py:class}`ExtensibleReactor` family of classes can be used to
+modify and governing equations, starting from one of the built-in reactor types.
+
+These reactor types allow the methods that implement the governing equations and related
+reactor configuration to be augmented or replaced with user-defined Python functions.
+New state variables can be added to the reactor, and existing ones can be redefined.
+User-defined `ExtensibleReactor` implementations can be used alongside the built-in
+reactor types to build reactor networks that can be integrated using the
+{py:class}`ReactorNet` integrator provided by Cantera.
+
+This guide introduces the use of the `ExtensibleReactor` model to modify the governing
+equations for a reactor.
+
+## Modifying Existing Governing Equations
+
+The variables in the governing equations for a reactor that are differentiated with
+respect to time are known as the *state variables*. The state variables depend on the
+type of reactor base class chosen. For example, choosing an
+[ideal gas constant pressure reactor](/reference/reactors/ideal-gas-constant-pressure-reactor)
+allows the user to modify the governing equations for the following state variables:
+
+- $m$, the mass of the reactor's contents (in kg)
+- $T$, the temperature (in K)
+- $Y_k$, the mass fractions for each species (dimensionless)
+
+As shown in the derivations of the governing equations, the equations implemented by
+the {ct}`IdealGasConstPressureReactor` class are:
+
+$$
+\frac{dm}{dt} = \sum_\t{in} \dot{m}_\t{in} - \sum_\t{out} \dot{m}_\t{out}
+ + \dot{m}_\t{wall}
+$$
+
+$$
+m c_p \frac{dT}{dt} = \dot{Q} - \sum_k h_k \dot{m}_{k,\t{gen}}
+ + \sum_\t{in} \dot{m}_\t{in} \left(h_\t{in} - \sum_k h_k Y_{k,\t{in}} \right)
+$$
+
+$$
+\frac{d(mY_k)}{dt} = \sum_\t{in} \dot{m}_\t{in} Y_{k,\t{in}}
+ - \sum_\t{out} \dot{m}_\t{out} Y_k + \dot{m}_{k,gen}
+$$
+
+Each of these equations is written with an expression on the left-hand side (LHS)
+multiplied by the time derivative of a state variable which equals an expression on the
+right-hand side (RHS). The interface of the reactor network model is defined so that the
+reactor governing equations are evaluated by calculating these LHS and RHS expressions
+for each governing equation. The extensible reactor model then allows calculation of
+these LHS and RHS terms to be replaced or augmented by user-provided methods.
+
+For example, to add a term for a large mass, say a rock, inside the reactor that affects
+the thermal mass, the energy equation would become:
+
+$$
+\left(m c_p + m_\t{rock} c_{p,\t{rock}}\right) \frac{dT}{dt} = \dot{Q}
+ - \sum_k h_k \dot{m}_{k,\t{gen}}
+ + \sum_\t{in} \dot{m}_\t{in} \left(h_\t{in} - \sum_k h_k Y_{k,\t{in}} \right)
+$$
+
+Here, the LHS coefficient has changed from $m c_p$ to
+$m c_p + m_{\t{rock}} c_{p,\t{rock}}$. Since the rock does not change the
+composition of the species in the reactor and does not change the mass flow rate of any
+inlets or outlets, the other governing equations defining the ideal gas constant
+pressure reactor can be left unmodified. To implement this change, we define a new class
+derived from {py:class}`ExtensibleIdealGasConstPressureReactor`. In this new class, we
+then define an implementation for the `after_eval` method which has the signature
+
+```py
+def after_eval(self, t, LHS, RHS) -> None
+```
+
+where `t` is the current simulation time and `LHS` and `RHS` are arrays with one element
+for each state variable. In this case, the length of each array is two plus the number of species.
+
+Optional user-defined methods for the `ExtensibleReactor` class start with the prefix
+`before_`, `after_`, or `replace_`, indicating whether the this method will called
+before, after, or instead of the corresponding method from the base class. In this case,
+our method will be called after the `eval` method of {ct}`IdealGasConstPressureReactor`,
+with the values of `LHS` and `RHS` already set based on the original governing
+equations. Instead of `after_eval`, we could implement either `before_eval` or
+`replace_eval` methods. However, in the case of `before_eval`, any element of `LHS` that
+we set would be overwritten by the initial implementation, and if we use `replace_eval`,
+we would have needed to calculate values for all elements of `RHS` and `LHS` instead of
+just modifying a single value.
+
+To obtain the above modification to the governing equations, the implementation of
+`after_eval` inside a new reactor class is then:
+
+```{code-cell} python
+import cantera as ct
+
+class RockReactor(ct.ExtensibleIdealGasConstPressureReactor):
+ def __init__(self, *args, mass_rock, cp_rock, **kwargs):
+ super().__init__(*args, **kwargs)
+ self.mass_rock = mass_rock
+ self.cp_rock = cp_rock
+
+ def after_eval(self, t, LHS, RHS):
+ # The index 1 corresponds to the governing equation for the second
+ # state variable (temperature)
+ LHS[1] += self.mass_rock * self.cp_rock
+```
+
+We can then use this new reactor class in a reactor network with some different values
+of `mass_rock` to see the impact of the modified governing equations. To do this, we
+can define a reactor network with an imposed heat transfer rate.
+
+
+```{code-cell} python
+cp_rock = 790 # J/kg/K
+
+for mass_rock in [0.0, 1.0, 3.0]:
+ # Define objects, properties, and initial conditions.
+ gas = ct.Solution('h2o2.yaml')
+ gas.TPX = 300, ct.one_atm, {'O2': 1, 'N2': 3.76}
+ T0 = gas.T
+
+ r1 = RockReactor(gas, mass_rock=mass_rock, cp_rock=cp_rock)
+ r1.volume = 2.0 # volume (not including 'rock')
+
+ env = ct.Reservoir(gas)
+ wall = ct.Wall(env, r1, Q=1e4, A=1.0) # Heat transfer of 10 kW
+ net = ct.ReactorNet([r1])
+
+ # Integrate for 10 seconds
+ net.advance(10)
+
+ print(f'mass_rock = {mass_rock} kg; ΔT = {r1.thermo.T - T0:.2f} K')
+```
+
+As expected, we see that the temperature rise decreases as the thermal mass of the
+reactor's contents increases.
+
+Other modifications to the reactor's governing equations and behavior may be made by
+implementing `before_*`, `after_`, and `replace_*` versions of any of the reactor
+methods listed in the documentation for the {py:class}`ExtensibleReactor` class.
+
+## Additional Examples
+
+[Using `ExtensibleReactor` to implement wall inertia](/examples/python/reactors/custom2)
+: This example shows how to introduce additional governing equations to a reactor. Here,
+ in addition to the `eval` method, overrides are introduced that modify the behavior of
+ the `initialize`, `get_state`, `update_state`, `component_index`, and `component_name`
+ methods. The result is a reactor network with a wall between two reactors that has
+ inertia and takes time to respond to changes in the reactor's pressure.
+
+[Reactor cascade model for reactive flows in inert porous media](/examples/python/reactors/PorousMediaBurner)
+: This example showcases the use of {py:class}`ExtensibleReactor` to add a temperature
+ equation for a solid phase and custom heat transfer/radiation models. Each reactor in
+ the network represents a different spatial locations along the length of a porous
+ media burner.
diff --git a/doc/sphinx/userguide/index.md b/doc/sphinx/userguide/index.md
index 869c10d35d..4a9ae1c502 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
@@ -27,21 +28,32 @@ The tutorials in this section are designed to help you accomplish a specific tas
using Cantera, such as evaluating the ignition delay time for a fuel under different
conditions, or calculating the voltage of a Lithium-ion battery as it is discharged.
+### Working with Input Data
+
- [](ck2yaml-tutorial)
- [](creating-mechanisms)
- [](thermobuild)
- [](input-errors)
- [](legacy2yaml-tutorial)
+### Implementing Custom Models
+
+- [](extensible-reactor)
+
```{toctree}
:hidden:
python-tutorial
input-tutorial
+reactor-tutorial
+
faq
+
ck2yaml-tutorial
creating-mechanisms
thermobuild
input-errors
legacy2yaml-tutorial
+
+extensible-reactor
```
diff --git a/doc/sphinx/userguide/reactor-tutorial.md b/doc/sphinx/userguide/reactor-tutorial.md
new file mode 100644
index 0000000000..155a609167
--- /dev/null
+++ b/doc/sphinx/userguide/reactor-tutorial.md
@@ -0,0 +1,218 @@
+---
+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.
+# This mechanism is included with Cantera for use in examples, but is not
+# recommended for research use.
+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 absolute time of 1 second, and examine
+how the state of the gas has changed. There isn't anything special about choosing 1
+second as the end time for the integration; any end time can be chosen. We are choosing
+1 second here because after that amount of time, the mixture will have almost certainly
+ignited for a wide range of initial conditions.
+
+```{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_\t{max}$. The
+ new time $t_\t{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_\t{new}$. $t_\t{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_\t{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:
+
+1. the absolute tolerances
+2. the relative tolerances
+3. the maximum time step size
+
+Reducing the third value is particularly useful when dealing with abrupt changes in the
+boundary conditions. One example of this is opening or 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
+
+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. Using these solvers can greatly speed up integration for systems with
+hundreds or thousands of species {cite:p}`walker2023`.
+
+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 the time integration are the same as before:
+```{code-cell} python
+sim.advance(1)
+reac.thermo()
+```
+
+The approximate Jacobian matrices used to construct the preconditioners do not currently
+account for terms that link multiple reactors in the same network. For networks of this
+type, the iterative solver may exhibit convergence errors for systems where the default
+direct solver does not. If the solver converges, the solution will satisfy the error
+tolerances, even though the Jacobian matrix is not exact.
+
+## 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 mixture changes state in time when
+isolated from the surroundings, 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.
+
+The initial state of the solution can be specified by composition and a set of intensive
+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 {py:class}`ReactorNet` containing only the single
+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.
+
+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 control
+volume, $m$, divided by $\dot{m}$ defines the mean residence time of the fluid in the
+control volume.
+
+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. This method is demonstrated in the
+[`combustor.py`](/examples/python/reactors/combustor) example.
+:::
diff --git a/doc/sphinx/yaml/reactions.rst b/doc/sphinx/yaml/reactions.rst
index 932d550f2f..132e16560a 100644
--- a/doc/sphinx/yaml/reactions.rst
+++ b/doc/sphinx/yaml/reactions.rst
@@ -47,9 +47,8 @@ The fields common to all ``reaction`` entries are:
``orders``
An optional mapping of species to explicit reaction orders to use. Reaction
orders for reactant species not explicitly mentioned are taken to be their
- respective stoichiometric coefficients. See
- `Reaction orders `__
- for additional information.
+ respective stoichiometric coefficients. See :ref:`sec-reaction-orders` for
+ additional information.
``negative-orders``
Boolean indicating whether negative reaction orders are allowed. The
@@ -92,9 +91,8 @@ Blowers-Masel
-------------
Blowers-Masel rate expressions calculate the rate constant based on the Blowers Masel
-approximation as
-`described here `__.
-The rate parameters are specified as a mapping with fields:
+approximation as :ref:`described here `. The rate parameters are
+specified as a mapping with fields:
``A``
The pre-exponential factor :math:`A`
@@ -119,7 +117,7 @@ Two-Temperature Plasma
Two-temperature plasma reactions involve an electron as one of the reactants, where the
electron temperature may differ from the gas temperature as
-`described here `__.
+:ref:`described here `.
The rate parameters are specified as a mapping with fields:
``A``
@@ -163,9 +161,8 @@ Reaction types
``elementary``
--------------
-A homogeneous reaction with a pressure-independent rate coefficient and mass
-action kinetics, as
-`described here `__.
+A homogeneous reaction with a pressure-independent rate coefficient and mass action
+kinetics, as :ref:`described here `.
Additional fields are:
@@ -188,8 +185,7 @@ Example::
``three-body``
--------------
-A three body reaction as
-`described here `__.
+A three body reaction as :ref:`described here `.
The reaction equation should include the third body collision partner ``M``.
@@ -203,6 +199,29 @@ Example::
rate-constant: [1.20000E+17 cm^6/mol^2/s, -1, 0]
efficiencies: {AR: 0.83, H2O: 5}
+*Changed in Cantera 3.0*: The ``type`` field of the YAML entry may be omitted. Reactions
+containing the generic third body M are automatically identified as three-body
+reactions. Reactions are also identified as three-body reactions if all of the following
+conditions are met:
+
+- There is exactly one species appearing as both a reactant and product
+- All reactants and products have integral stoichiometric coefficients
+- The sum of the stoichiometric coefficients for either the reactants or products is 3.
+
+Examples::
+
+ - equation: H + O2 + M <=> HO2 + M # Reaction 33
+ rate-constant: {A: 2.8e+18, b: -0.86, Ea: 0.0}
+ efficiencies: {O2: 0.0, H2O: 0.0, CO: 0.75, CO2: 1.5, C2H6: 1.5, N2: 0.0,
+ AR: 0.0}
+ - equation: H + 2 O2 <=> HO2 + O2 # Reaction 34
+ rate-constant: {A: 2.08e+19, b: -1.24, Ea: 0.0}
+ - equation: H + O2 + N2 <=> HO2 + N2 # Reaction 36
+ rate-constant: {A: 2.6e+19, b: -1.24, Ea: 0.0}
+
+.. caution::
+ If the third body efficiency of O2 and N2 in Reaction 33 was not set to zero, these
+ would be considered duplicate reactions and would be required to be marked as such.
.. _sec-yaml-Blowers-Masel:
@@ -227,7 +246,7 @@ Example::
Includes the fields of an :ref:`elementary ` reaction, except that
the ``rate-constant`` field is a
-:ref:`Two-temperature-plasma-type ` list or
+:ref:`Two-temperature-plasma ` list or
mapping.
Example::
@@ -242,16 +261,14 @@ Example::
``falloff``
-----------
-A falloff reaction as
-`described here `__.
+A falloff reaction as :ref:`described here `.
The reaction equation should include the pressure-dependent third body collision
partner ``(+M)`` or ``(+name)`` where ``name`` is the name of a species. The
latter case is equivalent to setting the efficiency for ``name`` to 1 and the
efficiency for all other species to 0.
-Includes field for specifying :ref:`efficiencies ` as well
-as:
+Includes field for specifying :ref:`efficiencies ` as well as:
``high-P-rate-constant``
An :ref:`sec-yaml-Arrhenius-rate` expression for the high-pressure limit
@@ -260,23 +277,18 @@ as:
An :ref:`sec-yaml-Arrhenius-rate` expression for the low-pressure limit
``Troe``
- Parameters for the
- `Troe `__
- falloff function. A mapping containing the keys ``A``, ``T3``, ``T1`` and
- optionally ``T2``. The default value for ``T2`` is 0.
+ Parameters for the :ref:`Troe ` falloff function. A mapping
+ containing the keys ``A``, ``T3``, ``T1`` and optionally ``T2``. The default value
+ for ``T2`` is 0.
``SRI``
- Parameters for the
- `SRI `__
- falloff function. A mapping containing the keys ``A``, ``B``, ``C``, and
- optionally ``D`` and ``E``. The default values for ``D`` and ``E`` are 1.0
- and 0.0, respectively.
+ Parameters for the :ref:`SRI ` falloff function. A mapping
+ containing the keys ``A``, ``B``, ``C``, and optionally ``D`` and ``E``. The default
+ values for ``D`` and ``E`` are 1.0 and 0.0, respectively.
``Tsang``
- Parameters for the
- `Tsang `__
- falloff function. A mapping containing the keys ``A`` and ``B``. The default value
- for ``B`` is 0.0.
+ Parameters for the :ref:`Tsang ` falloff function. A mapping
+ containing the keys ``A`` and ``B``. The default value for ``B`` is 0.0.
Example::
@@ -293,7 +305,7 @@ Example::
------------------------
A chemically activated reaction as
-`described here `__.
+:ref:`described here `.
The parameters are the same as for :ref:`sec-yaml-falloff` reactions.
@@ -311,7 +323,7 @@ Example::
--------------------------------
A pressure-dependent reaction using multiple Arrhenius expressions as
-`described here `__.
+:ref:`described here `.
The only additional field in this reaction type is:
@@ -336,7 +348,7 @@ Example::
-------------
A reaction parameterized as a bivariate Chebyshev polynomial as
-`described here `__.
+:ref:`described here `.
Additional fields are:
@@ -370,9 +382,8 @@ Example::
``interface-Arrhenius``
-----------------------
-A reaction occurring on a surface between two bulk phases, or along an edge
-at the intersection of two surfaces, as
-`described here `__.
+A reaction occurring on a surface between two bulk phases, or along an edge at the
+intersection of two surfaces, as :ref:`described here `.
Includes the fields of an :ref:`sec-yaml-elementary` reaction plus:
@@ -443,7 +454,7 @@ Example::
----------------------
A sticking reaction occurring on a surface adjacent to a bulk phase, as
-`described here `__.
+:ref:`described here `.
Includes the fields of an :ref:`sec-yaml-interface-Arrhenius` reaction plus:
diff --git a/doc/sphinx/yaml/species.rst b/doc/sphinx/yaml/species.rst
index 71c25bc75a..1b12bd6287 100644
--- a/doc/sphinx/yaml/species.rst
+++ b/doc/sphinx/yaml/species.rst
@@ -76,7 +76,7 @@ Fields of a species ``thermo`` entry used by all models are:
NASA 7-coefficient polynomials
------------------------------
-The polynomial form `described here `__,
+The polynomial form :ref:`described here `,
given for one or two temperature regions. Additional fields of a ``NASA7``
thermo entry are:
@@ -109,7 +109,7 @@ Example::
NASA 9-coefficient polynomials
------------------------------
-The polynomial form `described here `__,
+The polynomial form :ref:`described here `,
given for any number of temperature regions. Additional fields of a ``NASA9``
thermo entry are:
@@ -146,7 +146,7 @@ Example::
Shomate polynomials
-------------------
-The polynomial form `described here `__,
+The polynomial form :ref:`described here `,
given for one or two temperature regions. Additional fields of a ``Shomate``
thermo entry are:
@@ -179,7 +179,7 @@ Example::
Constant heat capacity
----------------------
-The constant heat capacity model `described here `__.
+The constant heat capacity model :ref:`described here `.
Additional fields of a ``constant-cp`` thermo entry are:
``T0``
diff --git a/include/cantera/kinetics/BlowersMaselRate.h b/include/cantera/kinetics/BlowersMaselRate.h
index 80117b7a70..d10db0b723 100644
--- a/include/cantera/kinetics/BlowersMaselRate.h
+++ b/include/cantera/kinetics/BlowersMaselRate.h
@@ -40,7 +40,7 @@ struct BlowersMaselData : public ReactionData
//! Blowers Masel reaction rate type depends on the enthalpy of reaction
/**
- * The Blowers Masel approximation @cite blowers2004 adjusts the activation energy
+ * The Blowers Masel approximation @cite blowers2000 adjusts the activation energy
* based on enthalpy change of a reaction:
*
* @f{eqnarray*}{
diff --git a/include/cantera/numerics/CVodesIntegrator.h b/include/cantera/numerics/CVodesIntegrator.h
index afaf01846c..f954cd59da 100644
--- a/include/cantera/numerics/CVodesIntegrator.h
+++ b/include/cantera/numerics/CVodesIntegrator.h
@@ -94,7 +94,7 @@ class CVodesIntegrator : public Integrator
void checkError(long flag, const string& ctMethod, const string& cvodesMethod) const;
size_t m_neq = 0;
- void* m_cvode_mem = nullptr;
+ void* m_cvode_mem = nullptr; //!< CVODES internal memory object
SundialsContext m_sundials_ctx; //!< SUNContext object for Sundials>=6.0
void* m_linsol = nullptr; //!< Sundials linear solver object
void* m_linsol_matrix = nullptr; //!< matrix used by Sundials
diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx
index e11dba7e75..4f98256e79 100644
--- a/interfaces/cython/cantera/reactor.pyx
+++ b/interfaces/cython/cantera/reactor.pyx
@@ -1215,7 +1215,7 @@ cdef class Valve(FlowDevice):
.. math:: \dot m = K_v*(P_1 - P_2)
- where :math:`K_v` is a constant set by the `set_valve_coeff` method.
+ where :math:`K_v` is a constant set using the `valve_coeff` property.
Note that :math:`P_1` must be greater than :math:`P_2`; otherwise,
:math:`\dot m = 0`. However, an arbitrary function can also be specified,
such that
diff --git a/src/kinetics/PlogRate.cpp b/src/kinetics/PlogRate.cpp
index 832359506f..d8cd0d9b9f 100644
--- a/src/kinetics/PlogRate.cpp
+++ b/src/kinetics/PlogRate.cpp
@@ -149,7 +149,10 @@ void PlogRate::validate(const string& equation, const Kinetics& kin)
}
if (err_reactions.size()) {
throw InputFileError("PlogRate::validate", m_input,
- "\nInvalid rate coefficient for reaction '{}'\n{}",
+ "\nInvalid rate coefficient for reaction '{}'\n{}\n"
+ "To fix this error, remove this reaction or contact the author of the\n"
+ "reaction/mechanism in question, because the rate expression is\n"
+ "mathematically unsound at the temperatures and pressures noted above.\n",
equation, to_string(err_reactions));
}
}