Jacobian backends and linear-solve strategies

JacobianBackend

Abstract supertype of every Jacobian-assembly strategy. The production subtype is Structured{S<:LinearSolveStrategy}; DenseAD_ForwardDiff is 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}

A backend bundled with the sparsity patterns and caches it assembles into, built for one (ProblemDef, Workspace) pair. solve builds its own from its jacobian backend; build one directly for compute_jacobian! or verify_jacobian_at.

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

Structured{S<:LinearSolveStrategy} <: JacobianBackend

Backend that assembles a StructuredJacobian by walking an operator list. The strategy parameter S fixes the assembly profile (sparsity_profile) and is reported by is_approximate.

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 the Structured backend's assembly strategies.

DirectSolve <: LinearSolveStrategy

Materialize A + Σⱼ densify(structuredⱼ) for a sparse direct factorization (KLU / UMFPACK). Predictable, robust, the verification anchor for correctness gates.

ApproxSolve <: LinearSolveStrategy

Drop every StructuredCoupling entirely and keep the local block A. 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.

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. A diagonal slot the profile stores but never writes (structural_zero_diagonal) counts as outside the support. Used by Structured(strategy = ApproxSolve()).

Without reference_spec the reference spec is built on every call. For repeated verification, build it 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.

Operator

Abstract supertype of every Jacobian operator. Operators are stateless tags: per-problem caches (AD workspaces, scratch buffers) live on the backend cache and are passed to the operator's methods.

StructuredJacobian

Structured representation of the Jacobian.

Fields

• A::SparseMatrixCSC{Float64,Int} — local sparse block. Block- tridiagonal physics for non-ALE, plus z-diagonal and χ block for ALE.
• structured::Vector{StructuredCoupling} — non-local blocks. Empty for Eulerian without Beer–Lambert; the node prefix sums for ALE, plus the rank-1 channels under WithChi.
• smap::FastSparsityMap — O(1) sparse-write helper for A.
• n::Int — full state dimension (rows == cols of A).

The operator is A + Σⱼ structuredⱼ. Assembly is performed by the backend (_compute_jacobian! in structured_backend.jl).

StructuredCoupling

Abstract supertype of every non-local block in a StructuredJacobian. Concrete subtypes implement densify!(M, c) — M .+= dense(c) — which the :full assembly profile uses to write the block into the Jacobian. Where the block acts in the global Jacobian is encoded in the coupling's offsets.

PrefixSumCoupling{T}

Per-cell factors v[i] = Δzᵢ · ∂θᵢ/∂xᵢ of the ALE cumulative mesh-velocity coupling (L · diag(v) · x)[k] = Σ_{i ≤ k-1} v[i] · x[i]. The mesh-velocity operator fills one per source block; compose_into! embeds it at the node rows as an EmbeddedNodePrefixSumCoupling.

EmbeddedScaledSuffixSumCoupling{T} <: StructuredCoupling

Embedded, per-row-scaled suffix sum 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 rank-1 block 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.