Jacobian backends and linear-solve strategies

JacobianBackend

Abstract supertype of every Jacobian-assembly strategy. The production subtype is Structured{S<:LinearSolveStrategy}; DenseAD_ForwardDiff is preserved as the verifier reference path.

Concrete subtypes should be lightweight bitstypes whenever possible — heavy state belongs in the cache field of JacobianSpec, not on the tag.

JacobianSpec{B<:JacobianBackend, S, C}

The single user-facing handle for Jacobian configuration. Carried as a keyword argument to solve(...).

Type parameters

• B: backend tag (concrete subtype of JacobianBackend)
• S: sparsity-pattern type. SparseMatrixCSC{Bool,Int} for the Structured backend, Nothing for the verifier-only DenseAD_ForwardDiff. The two sparsity fields share this type — they must both be present or both absent for a given backend.
• C: backend-specific cache type. Holds Dual workspaces, JacobianConfigs, reusable scratch arrays, etc. Each backend declares its own cache struct; sharing is forbidden (different backends would race over a shared build_ad_workspace Dual workspace).

Fields

• backend: instance of B
• sparsity_phys: physics-block (T, ξ, χ-row) sparsity. May be nothing.
• sparsity_geom: geometry-block (z-column) sparsity. ALE only; nothing otherwise.
• cache: backend-private workspace
• colorvec_phys: column coloring of sparsity_phys (for FD/sparse-AD). Int[] when not applicable.
• colorvec_geom: column coloring of sparsity_geom. Int[] when not applicable.

Structured{S<:LinearSolveStrategy} <: JacobianBackend

Backend that assembles a StructuredJacobian by walking an operator list. The strategy parameter S is propagated through is_approximate and consumed by the ldiv! chosen at solve time.

Constructor

Structured(; strategy = DirectSolve()) — the standard form. The strategy is a constructed instance, not a type, so future variants can carry per-strategy knobs in fields without touching this signature.

DenseAD_ForwardDiff <: JacobianBackend

Tag for the dense ForwardDiff backend. Carries no state.

compute_jacobian!(J, spec::JacobianSpec, u, def, t, ws)

Assemble the Jacobian into J using the backend carried by spec. This is the single dispatch entry point: backends specialize on the backend argument of _compute_jacobian!, never on compute_jacobian! itself.

def is expected to be a ProblemDef and ws a Workspace, but the dispatcher itself does not assert these types — it forwards to the backend, which is free to specialize on whatever problem-shape it understands.

LinearSolveStrategy

Abstract supertype of every strategy for inverting a StructuredJacobian. Concrete strategies override ldiv! on (out, strat, jac, b).

DirectSolve <: LinearSolveStrategy

Materialize A + Σⱼ densify(structuredⱼ) and factor with sparse-direct (KLU / UMFPACK). Predictable, robust, the verification anchor for correctness gates. Use this strategy first when proving a new backend contributes correctly; switch to StructuredSolve once correctness is established and perf becomes the goal.

StructuredSolve <: LinearSolveStrategy

Sparse-direct on A plus fast apply! on each StructuredCoupling, combined via preconditioned Richardson. Asymptotically O(N · n_structured) per iteration; the production target for ALE problems at large N when exact dense assembly is too expensive.

ApproxSolve <: LinearSolveStrategy

Drop every StructuredCoupling entirely; solve A⁻¹ b only. A Schur-style approximation — Newton sees a Jacobian that omits the dense ALE C block, takes a few extra iterations, but each linear solve is much cheaper.

This is the only strategy where missing-entry verification is permitted. The is_approximate trait extends to backends parameterized over this strategy, opting them into masked comparison in verify_jacobian_at.

MatrixFreeSolve <: LinearSolveStrategy

Use a GMRES iteration over the StructuredJacobian operator action. The structured couplings are applied by their fast apply! methods and are not densified into a concrete A + C matrix. Useful for strategy-level experiments and for very large N where materializing the coupled operator is itself expensive. Production callers should still prefer DirectSolve() unless a problem-specific benchmark shows the matrix-free path wins.

is_approximate(backend::JacobianBackend) -> Bool

Trait declaring whether backend deliberately omits known-non-zero entries from its declared support pattern. Default: false.

Set to true for backends that approximate the Jacobian by dropping couplings rather than computing them — e.g. Structured(strategy = ApproxSolve()), which drops the structured ALE mesh-velocity coupling. Used by verify_jacobian_at to decide whether to run the support-complete strict check (assert that the reference Jacobian is also zero outside the active backend's declared support) or fall back to the masked check (compare only inside support).

The default false makes the strict check the norm: any backend that "forgets" entries inside its declared support pattern AND has a sparsity pattern that omits the same entries will fail verification, instead of silently passing as both shortcuts cancel. Approximate backends opt in to the masked compare via this trait.

Adding the trait

struct MyApproximateBackend <: JacobianBackend end
Pyrolysis.Jacobian.is_approximate(::MyApproximateBackend) = true

verify_jacobian_at(spec_active::JacobianSpec, def, ws, u, t;
        reference::JacobianBackend = DenseAD_ForwardDiff(),
        rtol::Float64 = 1e-6,
        atol::Float64 = 1e-10) -> NamedTuple

Compute both spec_active's Jacobian and the reference-backend Jacobian at state (u, t) and report the maximum absolute and relative deviations.

Comparison modes (chosen by the is_approximate trait)

• Strict / support-complete (default — is_approximate(backend) == false). Every entry must satisfy the per-entry tolerance

  |ΔJ_ij| ≤ atol + rtol · max(|J_ref[i,j]|, maxⱼ|J_ref[i,:]|)

  The row-magnitude floor maxⱼ|J_ref[i,:]| keeps roundoff at true zeros of a stiff row from failing the check, while a wholly missing block in a small-magnitude row (e.g. the ALE z-rows, ~10 orders below the Arrhenius ξ-rows) still fails — the missing entry contributes to its own row's scale, so it cannot hide behind the global matrix maximum. Entries outside the backend's declared support are compared against an implicit zero, so a missing sparsity slot with a non-zero reference entry fails at that (row, col). This is the regime every "exact" backend must satisfy. The verifier additionally exposes max_abs_dev_outside_support so a missing-entry failure is distinguishable from an in-support disagreement.

• Masked (is_approximate(backend) == true). Entries outside the active backend's support are zeroed in the diff, and the pass criterion is the coarse global-scale rule max_abs_dev ≤ atol || max_abs_dev / max|J_ref| ≤ rtol. Approximate backends deliberately drop coupling contributions, so their in-support values legitimately deviate by the dropped magnitude; the masked compare is a sanity check, not an exactness gate. Used by Structured(strategy = ApproxSolve()) and any future backend tagged approximate.

The reference backend is built fresh each call (lazy AD workspace allocation on first hit, reused across calls within the same Julia session via JacobianSpec construction-time caching). For repeated verification, build the reference spec once and pass it via the reference_spec keyword (low-level form, see below).

Returned NamedTuple

• passed::Bool — strict: every entry within its per-entry tolerance; masked: global-scale criterion (see above).
• max_abs_dev::Float64 — largest absolute element-wise difference under the mode's comparison rule (strict: full diff, masked: in-support diff).
• max_rel_dev::Float64 — strict: |ΔJ| / max(|J_ref[i,j]|, row_max_i, 1e-12) at the worst (most tolerance-violating) entry; masked: max_abs_dev / max(maximum(abs, J_ref), 1e-12).
• worst_index::Tuple{Int,Int} — (row, col) of the worst entry (strict: largest violation ratio; masked: largest absolute deviation).
• max_abs_dev_in_support::Float64 — in-support disagreement.
• max_abs_dev_outside_support::Float64 — magnitude of reference entries the backend's support omits. For exact backends this should be within the per-entry tolerance; for approximate backends it is reported but not asserted.
• mode::Symbol — :strict or :masked.
• J_active::AbstractMatrix — the active backend's Jacobian (for inspection; do not retain across calls — backends may reuse buffers).
• J_ref::Matrix{Float64} — the dense reference Jacobian.

benchmark_jacobian_backends(def, ws, u;
        backends = [Structured(strategy = DirectSolve()),
                    Structured(strategy = ApproxSolve())],
        reference = DenseAD_ForwardDiff(),
        n_repeats::Int = 10,
        t::Float64 = 0.0) -> Vector{NamedTuple}

Run each backend on the given (def, ws, u) and report timing, allocations, and accuracy against the reference backend. Returns one NamedTuple per backend.

Arguments

• def, ws, u: as returned by canonical builders in perf_helpers.jl.
• backends: vector of JacobianBackend instances.
• reference: a JacobianBackend for the accuracy comparison. Defaults to DenseAD_ForwardDiff().
• n_repeats: number of timed repetitions; the minimum elapsed time is reported (filters scheduling jitter).
• t: simulation time at which to assemble.

Output schema

Per backend:

• name::String — backend type name.
• assembly_seconds::Float64 — minimum-over-n_repeats @elapsed.
• bytes::Int — @allocated after warm-up (zero is the goal for production backends).
• max_abs_err::Float64 — vs reference Jacobian on this backend's structural support.
• max_rel_err::Float64 — max_abs_err divided by max(maximum(abs, J_ref), 1e-12).

Operator

Abstract supertype of every Jacobian operator. Concrete subtypes implement the residual contribution (apply!) and either inherit the default AD-based contribute_to_jacobian! or override it with an analytical derivative.

Operators are stateless or near-stateless tags; per-problem caches (AD workspaces, scratch buffers) belong on the backend cache, not on the operator. This keeps the operator list cheap to construct and re-walk per Newton step.

ResidualContext

Holds the dependency graph for an operator list. Constructed via build_residual_context at problem-build time.

StructuredJacobian{Mode}

Structured representation of the Jacobian. The type parameter Mode records the simulation mode (Eulerian / ALE / WithChi) so dispatch on Mode selects between e.g. an empty structured list (Eulerian) and the prefix-sum + rank-1 list (ALE WithChi).

Fields

• A::SparseMatrixCSC{Float64,Int} — local sparse block. Block- tridiagonal physics for non-ALE, plus z-diagonal and χ block for ALE. May include the ALE limiter chain rule's enlarged stencil (within the 4-cell upwind window).
• structured::Vector{StructuredCoupling} — list of typed fast-apply linear operators. Empty for Eulerian; [PrefixSumCoupling] for ALE NoChi; [PrefixSumCoupling, Rank1Coupling] for ALE WithChi.
• implicit::Vector{ImplicitTerm} — IFT-resolved couplings. Empty in the AD-through-Newton path (the default).
• smap::FastSparsityMap — O(1) sparse-write helper for A.
• n::Int — full state dimension (rows == cols of A plus structured / implicit augmentation).

Apply:

(A + Σⱼ structuredⱼ + chain_rule_from_implicit) y = x

Assembly is performed by the backend (_compute_jacobian! in structured_backend.jl). Inversion is the responsibility of LinearSolveStrategy — see strategies.jl.

StructuredCoupling

Abstract supertype of every fast-apply linear-operator block in a StructuredJacobian. Concrete subtypes override:

• apply!(y, c, x) — y .+= c * x in-place (additive accumulation).
• apply_inverse!(y, c, x) — y = c⁻¹ x. Optional; only required for StructuredSolve. Concretes that have no analytical inverse may leave this unimplemented (errors at solve time, not assembly time).
• densify!(M, c) — M .+= dense(c). Used by DirectSolve to collapse the structured term into the sparse direct factorization.
• Base.size(c) — (rows, cols) of the operator's range / domain.

The coupling captures a derived field's chain-rule contribution; its sources and targets (which rows/columns of the global Jacobian it acts on) are encoded by how it is composed in the StructuredJacobian assembly driver, not as fields on the coupling itself.

PrefixSumCoupling{T} <: StructuredCoupling

Represents L · diag(v) where L is the lower-triangular all-ones matrix and v ∈ ℝ^N is a per-cell scaling. Equivalently: (L · diag(v) · x)[k] = Σ_{i ≤ k-1} v[i] · x[i] (exclusive prefix sum of v .* x).

For ALE this captures the cumulative mesh-velocity coupling: v[i] = Δzᵢ · ∂θ_i/∂x_i, and the apply produces ∂w_node/∂x for every node k from the cell-local strain-rate derivatives ∂θ_i/∂x_i.

Fields

• v::Vector{T} — per-cell scaling. Length N = number of cells.
• _buf::Vector{T} — internal scratch for in-place scans.

Apply / inverse cost

• apply!, apply_inverse!: O(N) (one pass).
• densify!: O(N²) (writes the full lower-triangular block).

The apply, inverse, and densify paths are all implemented and validated against dense references in test_structured_couplings.jl.

SuffixSumCoupling{T} <: StructuredCoupling

Represents U · diag(v) where U is the strict upper-triangular all- ones matrix (U[i, j] = 1 for j ≥ i+1, else 0) and v ∈ ℝ^N is a per-cell scaling. Equivalently: (U · diag(v) · x)[k] = Σ_{i ≥ k+1} v[i] · x[i] (exclusive suffix sum of v .* x).

This is the dual of PrefixSumCoupling. The right shape for non-local sources whose value at cell i depends on every cell with index greater than i. Beer-Lambert volumetric absorption is the motivating use: q_rad[i] depends on the optical depth integrated from the surface (cell N) down to cell i+1, so cell i's row has nonzero entries from cell i+1 through cell N.

Fields

• v::Vector{T} — per-cell scaling. Length N = number of cells.
• _buf::Vector{T} — internal scratch (reserved for future buffered scans; the present implementation is allocation-free without it).

Apply / inverse cost

• apply!, apply_inverse!: O(N) (one pass).
• densify!: O(N²) (writes the full upper-triangular block).

The matrix M = U · diag(v) is singular: its last row is zero and its first column is zero. So M y = x requires x[N] == 0 for consistency, and y[1] is in the right-nullspace. apply_inverse! sets y[1] = 0 by convention; downstream solvers compose this inverse with A, where the rank-deficient direction is pinned by the rest of the system.

Rank1Coupling{T} <: StructuredCoupling

Represents u · vᵀ — a dense rank-1 outer product. Captures genuinely low-rank couplings (e.g. column-global θ_A back-reaction in WithChi ALE mode).

Fields

• u::Vector{T} — left factor (columns of M = u * v').
• v::Vector{T} — right factor (rows of M).

Apply / inverse cost

• apply!: O(N) (one inner product + one scaled add).
• apply_inverse!: not analytically invertible on its own (the rank-1 matrix is singular); composes with the full system via Sherman- Morrison in StructuredJacobian's solve path. Leaving unimplemented at the coupling level — the SMW combination lives at the StructuredSolve strategy.
• densify!: O(N²).

EmbeddedScaledPrefixSumCoupling{T} <: StructuredCoupling

The ALE consumer-side wrapper that lifts a bare per-cell prefix-sum coupling (capturing ∂w_node[k]/∂x_j from ALEMeshVelocityOperator) into a global-sized linear operator. Captures one (target row block, source col block) channel of the indirect chain rule that flows from the cumulative-mesh-velocity coupling into a consumer's residual row (today: the ALEAdvectionOperator's T row and gas ξ rows).

What it represents

Given:

• v::Vector — per-cell Δzᵢ · ∂θ_i/∂x_i (the bare PrefixSumCoupling's scaling) for the source channel x (T or one ξ component).
• n_left[i], n_right[i] — the bounding nodes of cell i (so w_cell = 0.5·(w_node[n_left[i]] + w_node[n_right[i]])).
• row_scale[i] — the consumer's per-cell coefficient of w_cell (today: dT_dz_i for T-row contributions; -dξ_dz[j_species, i] for gas-row contributions).
• row_offset, col_offset — the global-state offsets that locate this block inside the full Jacobian. Cell i's contribution lands at [row_offset + i, col_offset + j] for source cell j.

The block represents

M[row_offset + i, col_offset + j] =
    row_scale[i] · 0.5 · (∂w_node[n_left[i]]/∂x_j + ∂w_node[n_right[i]]/∂x_j)
  = row_scale[i] · 0.5 · (𝟙[j ≤ n_left[i] - 1] + 𝟙[j ≤ n_right[i] - 1]) · v[j]

i.e. the densified prefix-sum + cell-centering averaging, scaled per target row.

Solver support

densify! writes the block into a dense matrix at the right offsets so DirectSolve factors the correct system. apply! performs the same action in O(N) through an exclusive prefix scan, which is what StructuredSolve and MatrixFreeSolve use to avoid densifying the coupled operator.

apply_inverse! is intentionally undefined for the embedded block by itself: the cell-centered embedding is rank-deficient, and inversion is well-defined only for the full (A + Σⱼ Cⱼ) operator.

Fields

• v::Vector{T} — per-cell scaling of the bare prefix sum.
• n_left::Vector{Int}, n_right::Vector{Int} — node ids bounding each cell. Stored as Vector{Int} (not NTuple) so the type is concrete.
• row_scale::Vector{T} — per-target-row consumer coefficient.
• row_offset::Int, col_offset::Int — global offsets.
• n_cells::Int — cached length(v). Helps densify! keep tight bounds.

EmbeddedScaledSuffixSumCoupling{T} <: StructuredCoupling

Mirror of EmbeddedScaledPrefixSumCoupling for non-local sources that read every cell above the current cell. Beer-Lambert volumetric absorption is the canonical case: q_rad[i] integrates the optical depth from the surface (cell n) down to cell i, so cell i's residual row depends on every cell k ≥ i.

What it represents

M[row_offset + i, col_offset + j] =
    row_scale[i] · (j ≥ i) · v[j]

where v[j] carries the per-cell scaling (typically ∂α_j/∂x_j × the optical-depth Jacobian factor) and row_scale[i] carries the consumer's per-cell coefficient (typically −q_rad[i] for the energy row).

Fields

• v::Vector{T} — per-cell scaling.
• row_scale::Vector{T} — per-target-row consumer coefficient.
• row_offset::Int, col_offset::Int — global offsets.
• n_cells::Int — cached length(v).

EmbeddedRank1Coupling{T} <: StructuredCoupling

Embeds a bare Rank1Coupling u · vᵀ at consumer-specified (row_offset, col_offset). Captures genuinely low-rank channels — most notably the column-global θ_A back-reaction in WithChi mode, where every target row's coefficient depends on a single scalar functional of the χ vector.

What it represents

M[row_offset + i, col_offset + j] += u[i] · v[j]

Fields

• u::Vector{T} — left factor (target rows). Length n_rows.
• v::Vector{T} — right factor (source columns). Length n_cols.
• row_offset::Int, col_offset::Int — global offsets.

ImplicitTerm

IFT-resolved coupling for a derived state that's the output of an internal Newton solve (today: only T_s from solve_boundary_temperature!). When the term's sources/targets are small and bounded (BC T_s touches O(NC) boundary entries), the contribution is folded into StructuredJacobian.A as auxiliary sparse rows during assembly. When sources or targets are dense, the implicit term is promoted to a StructuredCoupling instead.

Fields

• sources::Vector{Int} — global state indices whose values determine the implicit value (e.g. boundary cell index for T_s).
• rows::Vector{Int} — global Jacobian rows the chain rule flows into.
• dx_d_sources::Matrix{Float64} — ∂x_implicit / ∂(state[sources]) evaluated at the converged Newton fixed point.
• drhs_dx::Vector{Float64} — ∂(rhs[rows]) / ∂x_implicit.

No operator currently constructs an ImplicitTerm; the boundary T_s derivative is taken by the AD-through-Newton path instead.