Analytic Jacobian for 1D Flames

Overview

The nonlinear solver for 1D flames needs the Jacobian \(\mathbf{J}\) of the residual function to compute each Newton step. The Jacobian is reused across Newton steps and only re-evaluated when the steps stop making progress (see the nonlinear solver page). Even so, each re-evaluation is one of the more expensive operations in the solution process, so reducing its cost is worthwhile.

By default (domain.jacobian_mode = "auto"), the analytic Jacobian replaces the finite-difference perturbation for the species mass-fraction columns at interior grid points with formulas derived from the kinetics composition derivatives, falling back to finite differences for columns it does not cover. For a domain with \(K\) species and \(N\) grid points, the solution vector has length

\[ N_v = (K + c) N, \]

where \(c\) is the number of non-species components per grid point (the axial velocity \(u\), scaled radial velocity \(V\), temperature \(T\), pressure eigenvalue \(\Lambda\), and, for two-point-controlled flames, the axial mass flux). Of the \(N_v\) columns, the analytic mode handles the \(K\,(N-2)\) species columns at the \(N-2\) interior points directly. The remaining

\[ N_v - K (N-2) = c N + 2K \]

columns — every non-species column, plus all columns at the two boundary points — are still evaluated by finite differences. The relative savings therefore grow with the species fraction \(K/(K+c)\) of the system, i.e. with the size of the mechanism.

Those finite-difference columns are built by perturbation: each component is perturbed in turn and the resulting change in the residual gives one column of \(\mathbf{J}\). Because the discretization is block tridiagonal — the residual at grid point \(j\) depends only on the solution at points \(j-1\), \(j\), and \(j+1\) — perturbing a variable at point \(p\) changes the residual only at points \(p-1\), \(p\), and \(p+1\). Each such column therefore costs a local residual evaluation over three grid points, not a sweep of the whole domain.

The derivation below covers the interior column points \(p = 1, \ldots, N-2\) (the two boundary points always use finite differences). All transport coefficients are treated as frozen at the base state, matching the frozen-transport approximation of the finite-difference Jacobian, so the two modes make equivalent approximations.

Notation

The notation follows the discretization page: grid points are denoted by a subscript, species mass fractions by \(Y_{k,j}\) (species \(k\) at grid point \(j\)) and mole fractions by \(X_{k,j}\), with \(W_k\) the molecular weight of species \(k\) and \(\overline{W}_j\) the mean molecular weight at point \(j\). A Jacobian entry involves two grid points: the row (residual) point \(j\) and the column (perturbed-variable) point \(p\). We write the column we are computing as \(\partial / \partial Y_{m,p}\) — the derivative with respect to the mass fraction of species \(m\) at point \(p\). As in the finite-difference solver, this is an unnormalized perturbation: a single \(Y_{m,p}\) is varied while the other mass fractions at point \(p\) are held fixed (the mixture is not renormalized to sum to one). Summation indices over species are written as \(n\), and \(\delta_{km}\) is the Kronecker delta (\(1\) if \(k = m\), \(0\) otherwise).

Species residual

Using the mixture-averaged diffusive flux and an upwind convection scheme, the steady-state species residual that Cantera assembles at interior point \(j\) is

\[ F_{Y_k,j} = \frac{1}{\rho_j}\left[ W_k \dot\omega_{k,j} - \rho_j u_j \left.\frac{\partial Y_k}{\partial z}\right|_j - \zeta_j \bigl(j_{k,j+1/2} - j_{k,j-1/2}\bigr) \right], \]

where \(\zeta_j \equiv 2 / (z_{j+1} - z_{j-1})\) and \(j_{k,j+1/2}\) is the diffusive mass flux of species \(k\) at the midpoint between points \(j\) and \(j+1\), \(\dot\omega_{k,j}\) is the net molar production rate, and the convection derivative \(\partial Y_k/\partial z|_j\) is upwinded. (The transient term that the time-stepping solver adds to the diagonal is handled separately and is not part of the analytic species column.) Compared with the discretization page, this is the same residual divided through by \(\rho_j\), which is the form actually assembled in Flow1D::evalSpecies() and the reason the \(1/\rho_j\) factors appear in the Jacobian below.

With frozen transport, the midpoint flux is

\[ j_{k,j+1/2} = \frac{\tilde{D}_{k,j+1/2}}{\Delta z_j}\bigl(X_{k,j} - X_{k,j+1}\bigr) + Y_{k,j} S_{j+1/2}, \]

where

\[ S_{j+1/2} = -\sum_n \frac{\tilde{D}_{n,\,j+1/2}}{\Delta z_j} \bigl(X_{n,j} - X_{n,j+1}\bigr), \]

and \(\Delta z_j = z_{j+1} - z_j\). The notation \(\tilde{D}_{k,\,j+1/2}\) is a reminder that this is not a bare diffusion coefficient but the composite flux prefactor that multiplies the mole-fraction gradient. For the default molar-gradient, mixture-averaged case it is

\[ \tilde{D}_{k,\,j+1/2} = \frac{\rho\, W_k}{\overline{W}}\, D_{k,\mathrm{mix}}, \]

evaluated at the midpoint state, where \(D_{k,\mathrm{mix}}\) is the mixture-averaged diffusion coefficient of species \(k\); it therefore carries units of a mass-flux density coefficient rather than \(\mathrm{m^2/s}\). This is the quantity stored in Flow1D::m_diff, and (being a transport property) it is held frozen here. The term \(Y_{k,j}\,S_{j+1/2}\) is the correction flux that enforces \(\sum_k j_{k,\,j+1/2} = 0\).

The mass-gradient form replaces each \(X_{n,j}\) by \(Y_{n,j}\), sets \(\tilde{D}_{k,\,j+1/2} = \rho\, D_{k,\mathrm{mix}}^{\,\mathrm{mass}}\), and drops the \(\overline{W}/W\) factors that appear in the derivatives below.

Jacobian column for \(Y_{m,p}\)

The residual \(F_{Y_k,j}\) depends on \(Y_{m,p}\) through the two midpoint fluxes adjacent to point \(p\) — \(j_{k,\,p-1/2}\) (which has point \(p\) as its right endpoint) and \(j_{k,\,p+1/2}\) (which has point \(p\) as its left endpoint) — and, when \(j = p\), additionally through the reaction term, the upwinded convection term, and the \(1/\rho_j\) prefactor. The column therefore has nonzero entries only in rows \(j = p-1,\, p,\, p+1\):

\[\begin{split} \frac{\partial F_{Y_k,j}}{\partial Y_{m,p}} = \begin{cases} \displaystyle -\frac{\zeta_j}{\rho_j} \frac{\partial j_{k,\,p-1/2}}{\partial Y_{m,p}} & j = p-1 \\[12pt] \displaystyle -\frac{\zeta_j}{\rho_j} \left(\frac{\partial j_{k,\,p+1/2}}{\partial Y_{m,p}} - \frac{\partial j_{k,\,p-1/2}}{\partial Y_{m,p}}\right) + \frac{W_k}{\rho_j}\frac{\partial \dot\omega_k}{\partial Y_m} + R_{km}^{(p)} + C_{km}^{(p)} & j = p \\[12pt] \displaystyle +\frac{\zeta_j}{\rho_j} \frac{\partial j_{k,\,p+1/2}}{\partial Y_{m,p}} & j = p+1 \end{cases} \end{split}\]

Here \(R_{km}^{(p)}\) is the density chain rule term: because \(\rho_j = P\overline{W}_j/(RT_j)\) depends on composition through \(\overline{W}_j\), the \(1/\rho_j\) prefactor contributes

\[ R_{km}^{(p)} = \bigl(W_k\dot\omega_{k,p} - d_{k,p}\bigr)\, \frac{\overline{W}_p}{W_m\,\rho_p}, \qquad d_{k,p} = \frac{2\bigl(j_{k,\,p+1/2} - j_{k,\,p-1/2}\bigr)}{z_{p+1} - z_{p-1}}, \]

obtained from \(\partial(1/\rho_p)/\partial Y_{m,p} = \overline{W}_p / (W_m\,\rho_p)\) acting on the reaction and diffusion parts of the residual numerator. (The convection part has no density chain rule contribution because the \(\rho_j\) in \(\rho_j u_j\) cancels the \(1/\rho_j\) prefactor.) The convection term \(C_{km}^{(p)}\) is the derivative of the upwinded \(-u_j\,\partial Y_k/\partial z|_j\): it is species-diagonal (\(k=m\)) and contributes to whichever rows have point \(p\) in their upwind difference stencil.

Boundary rows \(j = 0\) and \(j = N-1\) are skipped: those residuals are fixed-value constraints with no dependence on interior species values.

Diffusive-flux derivatives

The two flux derivatives above are evaluated at the endpoints of the relevant interval. With frozen transport, differentiating \(j_{k,\,j+1/2}\) with respect to the mass fraction at its left endpoint (point \(j\)) gives

\[ \frac{\partial j_{k,\,j+1/2}}{\partial Y_{m,j}} = \frac{\tilde{D}_{k,\,j+1/2}}{\Delta z_j}\, f_j\,\bigl(\delta_{km} - X_{k,j}\bigr) + Y_{k,j}\,\frac{\partial S_{j+1/2}}{\partial Y_{m,j}} + \delta_{km}\, S_{j+1/2}, \]

and with respect to the right endpoint (point \(j+1\)),

\[ \frac{\partial j_{k,\,j+1/2}}{\partial Y_{m,j+1}} = -\frac{\tilde{D}_{k,\,j+1/2}}{\Delta z_j}\, f_{j+1}\,\bigl(\delta_{km} - X_{k,j+1}\bigr) + Y_{k,j}\,\frac{\partial S_{j+1/2}}{\partial Y_{m,j+1}}, \]

where the \(f_j = \overline{W}_j / W_{m,j}\) is the \(\partial X/\partial Y\) scale factor and the correction-flux derivatives are

\[ \begin{align}\begin{aligned} \frac{\partial S_{j+1/2}}{\partial Y_{m,j}} & = -f_j\!\left(\frac{\tilde{D}_{m,\,j+1/2}}{\Delta z_j} - a_j\right),\\\frac{\partial S_{j+1/2}}{\partial Y_{m,j+1}} & = +f_{j+1}\!\left(\frac{\tilde{D}_{m,\,j+1/2}}{\Delta z_j} - a_{j+1}\right), \end{aligned}\end{align} \]

with \(a_j = \sum_n (\tilde{D}_{n,\,j+1/2}/\Delta z_j)\, X_{n,j}\). The right-endpoint form has no \(\delta_{km}\,S\) term because the \(Y_{k,j}\,S_{j+1/2}\) part of the flux is anchored at the left endpoint. The implementation computes only the endpoint actually needed for each row.

Reaction derivatives: \(\partial\dot\omega_k / \partial Y_m\)

At constant \(T\) and \(P\), the chain rule from molar concentrations \(C_n = \rho Y_n / W_n\) to the unnormalized mass fractions gives

\[ \frac{\partial \dot\omega_k}{\partial Y_m} = \frac{\rho}{W_m}\frac{\partial \dot\omega_k}{\partial C_m} - \frac{\overline{W}}{W_m} \sum_n \frac{\partial \dot\omega_k}{\partial C_n}\, C_n. \]

The first term is the direct response to perturbing \(C_m\); the second arises because \(\rho\) and \(\overline{W}\) themselves depend on the unnormalized \(Y_m\). The sparse matrix \(\partial \dot\omega / \partial C\) is supplied by Kinetics::netProductionRates_ddCi(), whose nonzero pattern is fixed for a given mechanism and is reused across all grid points.

Energy and flow rows

The energy-row block \(\partial F_{T,j} / \partial Y_{m,p}\) (for \(j \in \{p-1,\,p,\,p+1\}\)) reuses the same diffusive-flux derivatives through the enthalpy-flux term \(\sum_k h_k\, j_k/W_k\), and adds a density chain rule correction together with a reaction-heat term and (when enabled) a radiation term at \(j = p\). The continuity and radial-momentum rows depend on the species mass fractions only through the density, so they receive only a density chain rule correction (\(\partial \rho_p / \partial Y_{m,p} = -\rho_p\,\overline{W}_p / W_m\)) from the column at point \(p\).