Implement cross-reactor Jacobian terms for preconditioners

doc/sphinx/reference/reactors/index.md

@@ -143,6 +143,102 @@ 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).
 
+For reactor networks containing multiple reactors, each reactor contributes the Jacobian
+rows corresponding to its own governing equations. Terms involving flow devices, walls,
+and reacting surfaces may therefore appear in columns for other reactors in the same
+network.
+
+The connector implementation is split between connector-specific scalar derivatives
+and reactor-specific state derivatives. Connectors such as valves, pressure
+controllers, and walls know simple derivatives such as
+$\partial \dot m / \partial (P_1 - P_2)$,
+$\partial \dot V / \partial (P_\t{left} - P_\t{right})$, or
+$\partial \dot Q / \partial T_\t{left}$. Reactor objects know how pressure,
+temperature, and flow-carried composition depend on their own state variables. Helper
+methods such as ``addPressureJacobian``, ``addTemperatureJacobian``,
+``addSpeciesMassFractionJacobian``, and ``addEnthalpyJacobian`` combine these pieces.
+
+For an {ct}`IdealGasMoleReactor`, the pressure derivatives used in connector terms are
+
+$$
+\frac{\partial P}{\partial T}\bigg|_V
+  = \frac{\pi_T + P}{T}, \qquad
+\frac{\partial P}{\partial V}\bigg|_T
+  = -\frac{1}{V \kappa_T}, \qquad
+\frac{\partial P}{\partial n_j}
+  \approx \frac{RT}{V},
+$$
+
+where $\pi_T = T(\partial P/\partial T)_V - P$ is the internal pressure and
+$\kappa_T = -(1/V)(\partial V/\partial P)_T$ is the isothermal compressibility. The
+first two derivatives are evaluated exactly using the equation of state via
+{ct}`ThermoPhase::internalPressure` and {ct}`ThermoPhase::isothermalCompressibility`,
+so they are correct for both ideal and non-ideal phases. The species-mole pressure
+derivative uses an ideal-gas approximation $\partial P/\partial n_j = RT/V$, which is
+exact for ideal gases. This approximation avoids EOS-specific partial-molar pressure
+evaluation and these terms are skipped by default (see ``skip-connector-pressure-composition-dependence``).
+
+For an {ct}`IdealGasConstPressureMoleReactor`, the reactor pressure is treated as fixed
+for these connector derivatives, so pressure-coupling terms from this reactor type are
+zero. Composition carried by a flow device is based on mass
+fractions because flow-device species fluxes are defined as $\dot m Y_k$. Mole reactors
+therefore convert those composition derivatives back to their native species-mole state
+variables using
+
+$$
+Y_k = \frac{W_k n_k}{m}, \qquad
+\frac{\partial Y_k}{\partial n_j}
+ = \frac{W_k \delta_{kj}}{m} - \frac{Y_k W_j}{m}
+$$
+
+where $\delta_{kj}$ is the Kronecker delta function.
+
+For inlets carrying enthalpy into a reactor, the connector terms also include
+derivatives of the upstream specific enthalpy with respect to the upstream reactor's
+state. The temperature contribution is
+
+$$
+\frac{\partial h_\t{in}}{\partial T_\t{in}} = c_{p,\t{in}},
+$$
+
+and the optional composition contribution is
+
+$$
+\frac{\partial h_\t{in}}{\partial n_{k,\t{in}}}
+ = \frac{\bar{h}_k - h_\t{in} W_k}{m_\t{in}},
+$$
+
+where $\bar{h}_k$ is the partial molar enthalpy of species $k$ in the upstream reactor.
+The composition contribution adds one entry per species to the preconditioner and is
+controlled by the ``skip-connector-composition-dependence`` setting.
+
+The wall heat-transfer terms illustrate the difference between the full Jacobian and
+the sparse approximation used for preconditioning. For a wall contribution to reactor
+$i$ written as
+
+$$
+\dot{T}_i = \frac{f_i \dot{Q}_w}{C_i},
+$$
+
+where $C_i = n_i \bar{c}_{v,i}$ for an {ct}`IdealGasMoleReactor` and
+$C_i = n_i \bar{c}_{p,i}$ for an {ct}`IdealGasConstPressureMoleReactor`, the full
+derivative is
+
+$$
+\frac{\partial \dot{T}_i}{\partial y_j}
+ = \frac{f_i}{C_i}\frac{\partial \dot{Q}_w}{\partial y_j}
+ - \frac{f_i \dot{Q}_w}{C_i^2}\frac{\partial C_i}{\partial y_j}.
+$$
+
+The second term includes temperature and composition derivatives of the reactor heat
+capacity. These terms are omitted from connector preconditioner entries because they
+add cost and, for composition derivatives, can add substantial fill-in. The implemented
+wall connector terms retain the numerator derivatives, for example
+$(f_i / C_i)\partial\dot{Q}_w/\partial T_\t{neighbor}$, which are the sparse
+cross-reactor couplings that most directly affect the iterative linear solve. Similar
+denominator-derivative terms are neglected for connector contributions involving
+enthalpy, internal energy, and pressure-expansion work.
+
 
 ```{toctree}
 :hidden:

doc/sphinx/reference/reactors/interactions.md

@@ -47,7 +47,7 @@ $$  \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
+upstream reactor/reservoir, which means that the flow rate is clamped 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
@@ -56,6 +56,15 @@ sufficiently large value, very small pressure differences will result in flow be
 the reactors that counteracts the pressure difference. However, excessively large values
 of $K_v$ may slow down the integrator or cause it to fail.
 
+For preconditioner construction, active (unclamped) valves contribute
+
+$$
+\pxpy{\dot m}{(P_1 - P_2)} = K_v g(t) f'(P_1 - P_2),
+$$
+
+and contribute zero when the flow is clamped to zero. The pressure derivatives with
+respect to reactor state variables are supplied by the adjacent reactor models.
+
 :::{admonition} Python Usage
 :class: tip
 
@@ -83,6 +92,11 @@ allows simple implementation of loops, in which exhaust gas from a reactor is fe
 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.
 
+For preconditioner construction, mass flow controllers have
+$\partial \dot m / \partial P = 0$. They may still contribute derivatives through the
+composition carried from the upstream reactor unless connector composition dependence is
+disabled by derivative settings.
+
 :::{admonition} Python Usage
 :class: tip
 
@@ -108,6 +122,15 @@ where $K_v$ is a proportionality constant and $f$ is a function of the pressure
 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.
 
+For preconditioner construction, pressure controllers include the derivative of their
+primary flow device and the pressure-correction derivative
+
+$$
+\pxpy{\dot m}{(P_1 - P_2)} = K_v f'(P_1 - P_2),
+$$
+
+when the resulting flow is active and unclamped.
+
 :::{admonition} Python Usage
 :class: tip
 
@@ -155,6 +178,27 @@ heat flux (W/m{sup}`2`). This definition is such that positive $q_0(t)$ implies
 transfer from the "left" reactor to the "right" reactor. Each of the user-specified
 terms defaults to 0.
 
+For sparse preconditioner construction, walls provide the derivatives:
+
+$$
+\pxpy{\dot V_w}{P_\t{left}} = K A,\qquad
+\pxpy{\dot V_w}{P_\t{right}} = -K A,
+$$
+
+where $\dot V_w = A v$ is the volumetric expansion rate [m³/s] associated with the
+wall's motion. Heat-transfer derivatives are:
+
+$$
+\pxpy{\dot Q_w}{T_\t{left}} = U A + 4 \epsilon \sigma A T_\t{left}^3,\qquad
+\pxpy{\dot Q_w}{T_\t{right}} = -U A - 4 \epsilon \sigma A T_\t{right}^3.
+$$
+
+The time-dependent velocity and heat-flux functions are treated as independent of the
+reactor state. When these wall derivatives are used in reactor energy equations for
+the sparse preconditioner, they are multiplied by the current inverse reactor heat
+capacity. Derivatives of that heat capacity with respect to temperature or composition
+are not included in the connector terms; see [](sec-reactor-preconditioning).
+
 :::{admonition} Python Usage
 :class: tip
 

doc/sphinx/userguide/reactor-tutorial.md

@@ -164,11 +164,13 @@ sim.advance(1)
 reac.phase()
 ```
 
-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.
+```{versionadded} 4.0
+The approximate Jacobian matrices used to construct the preconditioners now include
+sparse terms for supported flow devices, walls, and reacting surfaces connecting mole
+reactors. Some dense connector terms are omitted by default to preserve sparsity; see
+{ct}`Reactor::setDerivativeSettings` and [](#Kinetics.derivative_settings) for the
+derivative settings controlling these terms.
+```
 
 ## Common Reactor Types and their Implementation in Cantera
 

include/cantera/base/Delegator.h

@@ -10,6 +10,7 @@
 #include "cantera/base/Units.h"
 #include "cantera/base/ctexceptions.h"
 #include "cantera/base/ExtensionManager.h"
+#include "cantera/numerics/eigen_sparse.h"
 
 namespace Cantera
 {
@@ -267,6 +268,20 @@ class Delegator
         *m_funcs_sz_csr[name] = makeDelegate(name, func, when, m_base_sz_csr[name]);
     }
 
+    //! Set delegates for member functions with the signature `void(SparseTriplets&)`
+    //! @since New in %Cantera 4.0.
+    void setDelegate(const string& name,
+                     const function<void(SparseTriplets&)>& func,
+                     const string& when)
+    {
+        if (!m_funcs_v_vET.count(name)) {
+            throw NotImplementedError("Delegator::setDelegate",
+                "for function named '{}' with signature "
+                "'void(SparseTriplets&)'.", name);
+        }
+        *m_funcs_v_vET[name] = makeDelegate(func, when, *m_funcs_v_vET[name]);
+    }
+
     //! Store a handle to a wrapper for the delegate from an external language interface
     void holdExternalHandle(const string& name,
                             const shared_ptr<ExternalHandle>& handle) {
@@ -402,6 +417,17 @@ class Delegator
         m_base_sz_csr[name] = base;
     }
 
+    //! Install a function with the signature `void(SparseTriplets&)`
+    //! as being delegatable
+    //! @since New in %Cantera 4.0.
+    void install(const string& name,
+                 function<void(SparseTriplets&)>& target,
+                 const function<void(SparseTriplets&)>& base)
+    {
+        target = base;
+        m_funcs_v_vET[name] = &target;
+    }
+
     //! Create a delegate for a function with no return value
     template <typename BaseFunc, class ... Args>
     function<void(Args ...)> makeDelegate(
@@ -508,6 +534,7 @@ class Delegator
     //! - `sz` for `size_t`
     //! - `AM` for `AnyMap`
     //! - `US` for `UnitStack`
+    //! - `ET` for `Eigen::Triplet<double>`
     //! - prefix `c` for `const` arguments
     //! - suffix `r` for reference arguments
     //! - suffix `p` for span arguments (for example, `dp` for `span<double>`)
@@ -526,6 +553,7 @@ class Delegator
         function<void(double, span<double>, span<double>)>*> m_funcs_v_d_dp_dp;
     map<string,
         function<void(span<double>, span<double>, span<double>)>*> m_funcs_v_dp_dp_dp;
+    map<string, function<void(SparseTriplets&)>*> m_funcs_v_vET;
 
     // Delegates with a return value
     map<string, function<double(void*)>> m_base_d_vp;

include/cantera/zeroD/FlowDevice.h

@@ -10,6 +10,7 @@
 #include "cantera/base/global.h"
 #include "cantera/base/ctexceptions.h"
 #include "ConnectorNode.h"
+#include "cantera/numerics/eigen_sparse.h"
 
 namespace Cantera
 {
@@ -51,6 +52,16 @@ class FlowDevice : public ConnectorNode
         throw NotImplementedError("FlowDevice::setMassFlowRate");
     }
 
+    //! Derivative of mass flow rate with respect to pressure difference.
+    //!
+    //! Returns `d(mdot)/d(P_in - P_out)` for the current state. A return value of zero
+    //! indicates that pressure coupling does not contribute Jacobian entries.
+    //!
+    //! @since New in %Cantera 4.0.
+    virtual double massFlowRate_ddP() const {
+        return 0.0;
+    }
+
     //! Get the device coefficient (defined by derived class).
     //! @since  New in %Cantera 3.2.
     double deviceCoefficient() const {
@@ -73,6 +84,61 @@ class FlowDevice : public ConnectorNode
     //! is not present in the upstream mixture.
     double outletSpeciesMassFlowRate(size_t k);
 
+    //! Add Jacobian terms proportional to derivatives of the mass flow rate.
+    //!
+    //! Adds entries for `coeff * d(mdot)/dy_j` to the specified global row of the
+    //! reactor network Jacobian. The flow device supplies the scalar derivative of
+    //! mass flow rate with respect to connector variables such as pressure drop while
+    //! the adjacent reactors supply derivatives of those variables with respect to
+    //! their state variables.
+    //!
+    //! @param[in,out] trips  Sparse Jacobian entries. Implementations append entries
+    //!     using global row and column indices in the reactor network.
+    //! @param row  Global row index receiving these chain-rule terms.
+    //! @param coeff  Multiplicative factor applied to `d(mdot)/dy_j`.
+    //! @param includePressureSpecies  Include pressure derivatives with respect to
+    //!     species state variables when pressure is a derived quantity. These terms
+    //!     may be dense and are controlled by preconditioner sparsity settings.
+    //! @since New in %Cantera 4.0.
+    virtual void addMassFlowRateJacobian(SparseTriplets& trips, size_t row,
+        double coeff, bool includePressureSpecies=true);
+
+    //! Add Jacobian terms proportional to derivatives of
+    //! `outletSpeciesMassFlowRate(k)`.
+    //!
+    //! Adds entries for `coeff * d(mdot * Y_k)/dy_j`, where `Y_k` is the upstream mass
+    //! fraction mapped to downstream species `k`.
+    //!
+    //! @param[in,out] trips  Sparse Jacobian entries. Implementations append entries
+    //!     using global row and column indices in the reactor network.
+    //! @param row  Global row index receiving these chain-rule terms.
+    //! @param k  Species index in the downstream reactor's phase.
+    //! @param coeff  Multiplicative factor applied to `d(mdot * Y_k)/dy_j`.
+    //! @param includeComposition  Include derivatives of upstream mass fraction `Y_k`
+    //!     with respect to upstream reactor state variables.
+    //! @param includePressureSpecies  Include pressure derivatives with respect to
+    //!     species state variables when pressure is a derived quantity.
+    //! @since New in %Cantera 4.0.
+    void addOutletSpeciesMassFlowRateJacobian(
+        SparseTriplets& trips, size_t row, size_t k, double coeff,
+        bool includeComposition=true, bool includePressureSpecies=true);
+
+    //! Add Jacobian terms for the inlet enthalpy dependence on upstream state.
+    //!
+    //! Adds entries for `coeff * mdot * d(h_in)/dy_j`, where `h_in` is the specific
+    //! enthalpy of the upstream reactor. This captures how the energy carried into a
+    //! downstream reactor changes when the upstream temperature or composition changes.
+    //!
+    //! @param[in,out] trips  Sparse Jacobian entries.
+    //! @param row  Global row index receiving these chain-rule terms.
+    //! @param coeff  Multiplicative factor applied to `mdot * d(h_in)/dy_j`.
+    //! @param includeComposition  Include derivatives of upstream enthalpy with respect
+    //!     to upstream species moles. These terms add composition-mediated fill-in and
+    //!     are controlled by preconditioner sparsity settings.
+    //! @since New in %Cantera 4.0.
+    void addInletEnthalpyJacobian(SparseTriplets& trips, size_t row, double coeff,
+                                   bool includeComposition=true);
+
     //! specific enthalpy
     double enthalpy_mass();
 
@@ -111,7 +177,7 @@ class FlowDevice : public ConnectorNode
     //! depends on the derived flow device class.
     //! @since  Changed in %Cantera 3.2. Previous version used a raw pointer.
     virtual void setPressureFunction(shared_ptr<Func1> f) {
-        m_pfunc = f.get();
+        m_pfunc = f;
     }
 
     //! Return current value of the time function.
@@ -129,7 +195,7 @@ class FlowDevice : public ConnectorNode
     //! depends on the derived flow device class.
     //! @since  Changed in %Cantera 3.2. Previous version used a raw pointer.
     virtual void setTimeFunction(shared_ptr<Func1> g) {
-        m_tfunc = g.get();
+        m_tfunc = g;
     }
 
     //! Set current reactor network time
@@ -141,13 +207,18 @@ class FlowDevice : public ConnectorNode
     }
 
 protected:
+    //! Return the derivative of the pressure function at `deltaP`.
+    //! @param deltaP  Pressure difference `P_in - P_out` [Pa].
+    //! @since New in %Cantera 4.0.
+    double pressureFunction_ddP(double deltaP) const;
+
     double m_mdot = Undef;
 
     //! Function set by setPressureFunction; used by updateMassFlowRate
-    Func1* m_pfunc = nullptr;
+    shared_ptr<Func1> m_pfunc;
 
     //! Function set by setTimeFunction; used by updateMassFlowRate
-    Func1* m_tfunc = nullptr;
+    shared_ptr<Func1> m_tfunc;
 
     //! Coefficient set by derived classes; used by updateMassFlowRate
     double m_coeff = 1.0;