11. Solver and Jacobian Configuration; Extensions

This chapter is the reference for tuning how Pyrolysis.jl solves a problem, as opposed to what problem it solves (Chapters 3–7) or how to read the result (Chapter 9). It covers choosing a Jacobian backend and linear-solve strategy by problem size and simulation mode, configuring the integrator's linear solver (sparse-direct or ILU-preconditioned Krylov), loading the optional vendor and sensitivity extensions, debugging Jacobian correctness with verify_jacobian_at, and dropping down to the pure-residual layer (ProblemDef / Workspace / residual!) when you need to build a custom solver loop. The mathematics behind the structured Jacobian and the implicit-RK W matrix is derived in Technical Reference §15; this chapter is about the knobs.

Two distinct linear systems appear in a stiff implicit solve, and Pyrolysis.jl lets you configure each one separately. Keeping them straight is the key to this chapter:

1. The Jacobian backend (jacobian = ... keyword) decides how the matrix J = ∂(du/dt)/∂u is assembled — by structured operator assembly (the default: solve auto-builds Structured(strategy = DirectSolve()) when jacobian = nothing), by dense automatic differentiation (verification only), or by letting the integrator build it via colored finite differences (the fallback used when depletion handling disables the pluggable backend).
2. The integrator's linear solver (linear_solver = ... keyword) decides how the implicit-stage system W x = b, with W = I/(γ_RK Δt) − J, is solved once J exists — by sparse-direct factorization (Sparspak, KLU, Pardiso) or by a preconditioned Krylov iteration (ILU-GMRES, ILU-BiCGSTAB).

You can mix and match: a pluggable structured Jacobian with a sparse-direct linear solver, an integrator-built finite-difference Jacobian with KLU, and so on. The defaults (auto-built structured Jacobian; SparspakFactorization) work well for the typical 1D cone-calorimeter problem and you can ignore this chapter until a solve is too slow or you need sensitivities.

11.1 The two layers, and what the defaults do

When you call solve(problem) with no Jacobian or linear-solver keywords, the following happens (see src/Solver/solve.jl):

• The integrator defaults to KenCarp4(autodiff = AutoFiniteDiff(), linsolve = SparspakFactorization()). KenCarp4 is an SDIRK (singly diagonally implicit Runge–Kutta) method well suited to the moderate stiffness of Arrhenius pyrolysis. The AutoFiniteDiff() setting only matters on the rare paths where no pluggable Jacobian is active (depletion handling): there it makes the integrator finite-difference the residual itself.
• The Jacobian defaults to the structured analytical backend: with jacobian = nothing (the default), solve constructs Structured(strategy = DirectSolve()) internally — for both Eulerian and ALE problems — and hands the integrator an exact sparse Jacobian plus its sparsity prototype. The single exception is handle_depletion = true, which disables the pluggable backend (the spec cannot survive a mid-solve state resize) and falls back to the integrator's colored finite differences.
• The linear solver defaults to SparspakFactorization(), a sparse-direct LU that is typically faster than KLU for the block-tridiagonal sparsity pattern of a 1D problem.

In short: by default, Layer 1 is the structured analytical Jacobian (Structured(strategy = DirectSolve()), built for you) and Layer 2 is Sparspak. The jacobian keyword is for overriding the strategy (ApproxSolve, StructuredSolve, MatrixFreeSolve) or supplying a custom spec — not for switching the structured machinery on.

using Pyrolysis

# The default solve: auto-built structured Jacobian + Sparspak linear solve.
sol = solve(problem)

Defaults at a glance: integrator = nothing → KenCarp4, linear_solver = nothing → SparspakFactorization(), jacobian = nothing → auto-built Structured(strategy = DirectSolve()) (colored-FD fallback only under handle_depletion = true), abstol = 1e-8, reltol = 1e-6.

11.2 Choosing a Jacobian backend

A Jacobian backend is supplied as a JacobianSpec to the jacobian keyword of solve. Construct it from a backend tag plus your problem:

spec = JacobianSpec(Structured(strategy = DirectSolve()), problem)
sol  = solve(problem; jacobian = spec)

There are two production-relevant backend tags, both exported at top level:

Structured is parameterized by a linear-solve strategy (§11.3) that determines how the assembled structured Jacobian is inverted. The default, selected when you pass jacobian = nothing, is Structured(strategy = DirectSolve()).

DenseAD_ForwardDiff() runs ForwardDiff over the entire residual vector and returns a dense matrix. It is correct by construction but slow; its only role is as the reference against which the structured backend is verified (see §11.7). Do not use it for production solves.

11.2.1 The JacobianSpec constructor

JacobianSpec has two forms (src/Jacobian/interface.jl):

# From a PyrolysisProblem (the common form). Pass use_ale = true for ALE.
JacobianSpec(backend, problem; use_ale = false, radiation_model = nothing,
             radiation_intensity = 0.0, energy_form = :standard,
             use_transpiration_bc = false, min_thickness = 1e-6)

# From an already-built (ProblemDef, Workspace) pair (low-level form).
JacobianSpec(backend, def, ws)

Important: a JacobianSpec builds its sparsity pattern and AD workspaces once, at construction time, and it bakes in the simulation mode. If your solve uses ALE, you must tell the spec by passing use_ale = true — otherwise the sparsity pattern omits the geometry (z-column) block and the Jacobian will be wrong:

# WRONG for an ALE solve — the spec is built for the Eulerian mode.
spec = JacobianSpec(Structured(strategy = ApproxSolve()), problem)
solve(problem; jacobian = spec, use_ale = true)   # mode mismatch

# RIGHT — the spec's mode matches the solve.
spec = JacobianSpec(Structured(strategy = ApproxSolve()), problem; use_ale = true)
solve(problem; jacobian = spec, use_ale = true)

If a problem also uses radiation, pass the matching radiation_model (e.g. radiation_model = BEER_LAMBERT) so the spec includes the Beer–Lambert non-local sparsity. Similarly, pass energy_form = :with_generation_sink if you solve with that energy form. For most users the simplest correct approach is to let solve build the spec for you — see §11.2.3.

11.2.2 When not to supply a Jacobian

A pluggable JacobianSpec is incompatible with surface-depletion handling. Depletion (handle_depletion = true) resizes the state vector mid-solve via resize!(integrator, ...), but the spec's cache is built once at setup time and cannot survive a resize. solve rejects the combination at validation:

# ERROR: handle_depletion + a pluggable jacobian is not allowed.
solve(problem; jacobian = spec, handle_depletion = true, use_ale = true)

When you need depletion, leave jacobian = nothing; the integrator falls back to its own colored finite-difference Jacobian, which tolerates the resize. This is the only supported configuration for depletion-tracking solves (see Chapter 10 and Technical Reference §14.7).

11.2.3 Letting solve pick the backend

For routine work, the cleanest pattern is to pass jacobian = nothing (the default) and let solve construct the default structured Jacobian for you, or to opt into a non-default backend explicitly only when you need it. The auto-default when jacobian is left at its default is Structured(strategy = DirectSolve()) constructed internally with the correct mode flags — you do not have to thread use_ale / radiation_model through JacobianSpec yourself. Supply your own spec only when you want a non-default strategy (ApproxSolve, StructuredSolve, MatrixFreeSolve) or a custom reference for verification.

11.3 Linear-solve strategies for the structured Jacobian

The strategy argument to Structured(...) is a LinearSolveStrategy. It governs how the structured Jacobian — a sparse local block A plus a list of fast-apply non-local couplings (prefix-sum, suffix-sum, rank-1) — is inverted inside the Newton step. All four strategies are top-level exports.

Construct any of these as a parameterless instance:

JacobianSpec(Structured(strategy = DirectSolve()),     problem)               # default
JacobianSpec(Structured(strategy = StructuredSolve()), problem; use_ale = true)
JacobianSpec(Structured(strategy = ApproxSolve()),     problem; use_ale = true)
JacobianSpec(Structured(strategy = MatrixFreeSolve()), problem; use_ale = true)

11.3.1 DirectSolve — the default and the anchor

DirectSolve materializes A + Σⱼ densify(couplingⱼ) and factors it directly. It is exact and predictable. For Eulerian problems and modest ALE meshes it is the right choice and the one you get when you do nothing. It is also the backend you should verify against first when proving a configuration is correct.

spec = JacobianSpec(Structured(strategy = DirectSolve()), problem)
sol  = solve(problem; jacobian = spec)

11.3.2 StructuredSolve — exact, cheaper assembly at scale

StructuredSolve keeps the solve exact but avoids densifying the non-local ALE/radiation couplings. It factors only the sparse local block A, then runs a preconditioned Richardson iteration that applies each coupling through its O(N) apply! method:

F = lu(A); y = F \ b
repeat:  r = b − A·y − Σⱼ Cⱼ·y;  Δy = F \ r;  y ← y + Δy   (2–5 iterations)

Convergence is fast (2–5 iterations to machine precision) when the coupling magnitudes are small relative to A's diagonal — which is the usual ALE regime. Choose StructuredSolve when an ALE mesh is large enough that the dense assembly of DirectSolve dominates, but you still need an exact Newton step.

11.3.3 ApproxSolve — fastest, for large ALE

ApproxSolve drops the structured couplings entirely and solves A⁻¹ b. The Newton step is then a Schur-style approximation: it omits the dense ALE mesh-velocity coupling block, takes a few extra Newton iterations, but each linear solve and the assembly itself are far cheaper (the assembly profile drops to :ale_reduced, skipping the cumulative ALE widening, the geometry AD pass, the Beer–Lambert suffix, and the χ-column rank-1 term). This is the strategy the Douglas-fir ALE example uses for roughly an 8× speedup over the full Jacobian.

# Large ALE cone problem — fast approximate Jacobian.
spec = JacobianSpec(Structured(strategy = ApproxSolve()), problem; use_ale = true)
sol  = solve(problem; jacobian = spec, use_ale = true)

ApproxSolve is the only strategy that deliberately omits known-nonzero entries; it carries the is_approximate trait, which switches verify_jacobian_at (§11.7) into masked-comparison mode so the dropped entries are reported informationally rather than counted as failures.

11.3.4 MatrixFreeSolve — research / very large N

MatrixFreeSolve runs GMRES over the structured-Jacobian operator action without ever forming A + C. It is exact in the limit of GMRES convergence and useful when even materializing the coupled operator is expensive. It is marked experimental in the source; benchmark a specific problem before using it in production, and prefer DirectSolve/StructuredSolve unless the matrix-free path measurably wins.

11.3.5 A sizing rule of thumb

• Eulerian (fixed mesh), any size, no depletion → jacobian = nothing (default) is fine; it already gives you Structured(strategy = DirectSolve()).
• ALE, small/medium mesh → DirectSolve.
• ALE, large mesh, exactness required → StructuredSolve.
• ALE, large mesh, speed over exactness → ApproxSolve.
• Depletion tracking → no pluggable Jacobian (jacobian = nothing).

11.4 Configuring the integrator's linear solver

Independently of the Jacobian backend, you choose how the per-stage implicit system W x = b (W = I/(γ_RK Δt) − J, see Technical Reference §15.5) is solved. This is the linear_solver keyword. There are two families: sparse-direct factorizations and ILU-preconditioned Krylov iterations.

11.4.1 Sparse-direct factorizations

The default is SparspakFactorization() (about 28% faster than KLU for the block-tridiagonal 1D pattern). KLUFactorization() is a slightly more general-purpose alternative with the same interface. You can pass either by pairing it with a constructor integrator:

using LinearSolve: KLUFactorization, SparspakFactorization

# Default (no need to write this explicitly):
sol = solve(problem; linear_solver = SparspakFactorization())

# Switch to KLU:
sol = solve(problem; linear_solver = KLUFactorization())

Pairing rule: linear_solver is combined with integrator only when integrator is a constructor such as KenCarp4 (or left at its default). If you want a fully constructed algorithm with a non-default linear solver, build it yourself and pass it as the integrator: solve(problem; integrator = KenCarp4(linsolve = KLUFactorization())). Passing both a fully constructed integrator and a separate linear_solver raises an error.

A note on the no-Jacobian path: with handle_depletion = true (or a dense backend such as DenseAD_ForwardDiff), there is no sparse jac_prototype, so the sparse Sparspak/KLU factorizations cannot run. solve then defaults the linear solver to a dense LUFactorization() — and if you explicitly request a sparse factorization it warns and rebuilds your algorithm with LUFactorization instead. On the ordinary default path (auto-built structured Jacobian) the sparse prototype is present and Sparspak applies.

11.4.2 ILU-preconditioned GMRES

For large ALE problems where direct factorization of W dominates the cost, Pyrolysis.jl ships an ILU-preconditioned GMRES solver, ILUGMRESFactorization (exported by the Solver submodule; qualify as Pyrolysis.ILUGMRESFactorization or using Pyrolysis.Solver). It builds an incomplete-LU preconditioner from the concrete sparse W and iterates GMRES against it, falling back to diagonal scaling if the ILU factorization fails on a stiff transient.

# ILUGMRESFactorization(; τ = 0.1, gmres_rtol = 1e-8, restart = 30)
alg = KenCarp4(linsolve = ILUGMRESFactorization(τ = 0.1))

sol = solve(problem;
    integrator = alg,
    jacobian   = JacobianSpec(Structured(strategy = DirectSolve()), problem; use_ale = true),
    use_ale    = true)

Tuning parameters and their defaults (the field name τ is the ILU drop tolerance, written τ_ILU in the nomenclature):

Caveat: inexact linear solves can increase Newton rejections on Arrhenius transients, which sometimes raises total wall time despite cheaper individual solves. Always benchmark an ILU-GMRES configuration end-to-end against the default SparspakFactorization() before adopting it for a production study.

11.4.3 ILU-preconditioned BiCGSTAB

ILUBiCGSTABFactorization (exported by the Solver submodule; qualify or using Pyrolysis.Solver) is the BiCGSTAB analogue. It uses the same ILU preconditioner but a fixed-memory BiCGSTAB iteration (no restart parameter), which can be faster for some matrix structures:

# ILUBiCGSTABFactorization(; τ = 0.1, rtol = 1e-8)
alg = KenCarp4(linsolve = ILUBiCGSTABFactorization(τ = 0.1, rtol = 1e-8))
sol = solve(problem; integrator = alg,
    jacobian = JacobianSpec(Structured(strategy = DirectSolve()), problem; use_ale = true),
    use_ale  = true)

11.4.4 The MKL ipiv caveat and depletion

When surface depletion resizes the state vector mid-solve, the integrator must also resize its linear-solver workspace. The MKL-backed pivot array (ipiv) used by some factorizations does not resize cleanly across a state-vector length change. Depletion forces jacobian = nothing (§11.2.2), so there is no sparse Jacobian prototype and the sparse factorizations cannot run on this path: solve defaults the linear solver to Julia's native, resize-safe LUFactorization() and rebuilds an explicitly requested SparspakFactorization/KLUFactorization to LUFactorization, with a warning. A robust depletion configuration is therefore simply:

# Depletion-safe: no pluggable Jacobian; solve picks the resize-safe dense LU.
sol = solve(problem;
    use_ale             = true,
    handle_depletion    = true,
    depletion_threshold = 0.05)

The MKL pitfall survives only when you pass a fully constructed integrator whose linsolve is unset (LinearSolve's auto-selection may then pick MKL LU); solve warns when it detects that combination on a large depleting state.

11.5 Tolerances and integrator selection

The ODE error tolerances are independent of both Jacobian and linear-solver choices:

To use a different stiff integrator, pass its constructor (the linear solver and AutoFiniteDiff autodiff are wired in for you) or a fully built instance:

using OrdinaryDiffEq: Rodas5P, FBDF

# Constructor form — linear_solver is paired in automatically:
sol = solve(problem; integrator = Rodas5P, linear_solver = SparspakFactorization())

# Fully constructed form (do not also pass linear_solver):
sol = solve(problem; integrator = FBDF(linsolve = KLUFactorization()))

Loosening tolerances (reltol = 1e-4) speeds up exploratory runs; tighten them for convergence and validation studies (Technical Reference §18.5). Tolerance choice does not affect mass/energy bookkeeping closure — that closes to the integrator tolerance by construction (Technical Reference §16) — but it does affect physical accuracy.

11.6 Loading and using extensions

Optional functionality is delivered through Julia's package-extension system (Julia 1.9+). The hooks (klu_factorization, pardiso_factorization, amg_preconditioner, forward_sensitivity, adjoint_sensitivity, verify_mixing_derivatives) are exported stubs in core Pyrolysis.jl. Each stub errors with a load hint until you using the corresponding optional package, at which point the extension attaches a concrete method. No extension depends on another; all are independent peers.

11.6.1 Linear-solver extensions: KLU, Pardiso, AMG

These three extensions each expose a vendor solver under a stable Pyrolysis-level name so the benchmark harness and out-of-tree code can request it without core Pyrolysis taking a hard dependency.

using Pyrolysis
using KLU                              # triggers PyrolysisKLUExt

ls  = Pyrolysis.klu_factorization()    # -> LinearSolve.KLUFactorization()
sol = solve(problem; linear_solver = ls)

using Pyrolysis
using Pardiso                          # triggers PyrolysisPardisoExt

# Returns LinearSolve.MKLPardisoFactorize(); for users with an Intel MKL license.
ls  = Pyrolysis.pardiso_factorization()
sol = solve(problem; integrator = KenCarp4(linsolve = ls))

using Pyrolysis
using AlgebraicMultigrid               # triggers PyrolysisAMGExt
using SparseArrays, LinearAlgebra

# Smoothed-aggregation AMG preconditioner for a sparse matrix A, for use as a
# preconditioner inside a Krylov linsolve.
A    = sprand(100, 100, 0.05) + 100I   # illustrative sparse matrix
prec = Pyrolysis.amg_preconditioner(A)

Note that KLUFactorization is already reachable as LinearSolve.KLUFactorization (LinearSolve is a hard dependency); the KLU extension simply confirms KLU itself is loaded and exposes it under the stable name Pyrolysis.klu_factorization. The AMG preconditioner is meant for the preconditioner slot of a Krylov solver, not as a linear_solver on its own.

Stub behavior without the trigger. If you call any of these before loading the optional package, you get a MethodError pointing at the unloaded extension:

using Pyrolysis
Pyrolysis.klu_factorization()   # MethodError — `using KLU` not yet loaded

The fix is always the same: add the corresponding using line.

11.6.2 Sensitivity extension

forward_sensitivity and adjoint_sensitivity are stubs until you load SciMLSensitivity and Zygote, which activate PyrolysisSciMLSensitivityExt. Both take a p_inject closure that builds a fresh problem from a parameter vector θ (it must not mutate the base problem). The closure rebuilds the problem from the base problem's mesh and BCs with a material that carries the injected parameters; the exact way you parameterize the material is up to you (Chapter 13 develops realistic p_inject closures in full):

using Pyrolysis
using SciMLSensitivity, Zygote          # triggers PyrolysisSciMLSensitivityExt

# p_inject(base, θ) -> fresh PyrolysisProblem with θ injected.
# `make_material(θ)` is your own function that builds a Material from the
# parameter vector (e.g. setting reaction A, E). The fresh problem reuses the
# base problem's mesh, BCs, time span, and initial conditions.
function p_inject(base, θ)
    mat = make_material(θ)
    return PyrolysisProblem(base.mesh, mat, base.bc_set, base.tspan;
                            T_initial = base.T_initial,
                            output_times = base.output_times)
end

# Forward mode (ForwardDiff) — small parameter counts:
θ = [1e13, 2.0e5]
output, ∂out_∂θ = forward_sensitivity(base_problem, θ, p_inject;
                                      output_fn = sol -> maximum(sol.surface_T))

# Adjoint mode (Zygote + SciMLSensitivity) — large parameter counts, scalar loss:
loss, grad = adjoint_sensitivity(base_problem, θ, p_inject, sol -> sum(sol.surface_T);
                                 sensealg = InterpolatingAdjoint(autojacvec = ReverseDiffVJP()))

This chapter only covers the loading mechanics; the full sensitivity workflow (forward vs. adjoint choice, inverse analysis, Sobol global sensitivity) is the subject of Chapter 13, and the theory is in Technical Reference §17.

11.6.3 Symbolics verification extension (test-only)

verify_mixing_derivatives is a test-only stub activated by using Symbolics. It compares the package's analytical mixing-rule derivatives against Symbolics-derived ground truth. Production code never calls it; it exists to catch drift bugs in the mixing-rule Jacobian terms.

using Pyrolysis
using Symbolics                         # triggers PyrolysisSymbolicsExt

# Returns a per-rule max relative error (NaN for BRUGGEMAN, which has no closed form).
report = verify_mixing_derivatives(material; rtol = 1e-8, T_test = 600.0)

11.7 Debugging Jacobian correctness with verify_jacobian_at

When a custom or approximate Jacobian configuration produces suspicious Newton behavior (excessive rejections, divergence, or wrong steady state), compare it against the dense ForwardDiff reference at a state point. There are two ways: a one-shot call, and an in-solve callback.

11.7.1 One-shot verification

verify_jacobian_at takes an active spec, a (def, ws) pair, a state u, and a time t, and reports the deviation from the DenseAD_ForwardDiff reference. Its signature is verify_jacobian_at(spec, def, ws, u, t; reference = DenseAD_ForwardDiff(), rtol = 1e-6, atol = 1e-10):

using Pyrolysis

# Build the (def, ws) pair for the mode you intend to solve in.
def = to_problem_def(problem; use_ale = true)
ws  = build_workspace(problem, def)
u   = create_state_vector(ws.layout, ws.mesh)

spec = JacobianSpec(Structured(strategy = ApproxSolve()), def, ws)

report = verify_jacobian_at(spec, def, ws, u, 0.0;
                            reference = DenseAD_ForwardDiff(),
                            rtol = 1e-6, atol = 1e-10)

report.passed                      # Bool — strict: per-entry |ΔJ_ij| ≤ atol + rtol·max(|J_ref[i,j]|, row max);
                                   #        masked: max_abs_dev ≤ atol || max_rel_dev ≤ rtol
report.max_abs_dev                 # largest element-wise deviation under the mode's rule
report.max_rel_dev
report.worst_index                 # (row, col) of the worst deviation
report.mode                        # :strict (exact backends) or :masked (approximate)
report.max_abs_dev_in_support      # in-support disagreement
report.max_abs_dev_outside_support # reference magnitude the backend's support omits

Comparison mode is chosen automatically by the backend's is_approximate trait (§11.3.3):

• Strict / support-complete (exact backends, e.g. DirectSolve, StructuredSolve): every entry of the dense reference is compared against the active backend; a non-zero reference entry that the backend's support omits surfaces as a deviation, and its magnitude appears in max_abs_dev_outside_support. The strict criterion is per-entry with a row-magnitude floor, so a missing coupling in a small-magnitude row cannot hide behind the global matrix maximum. This is the regime exact backends must satisfy.
• Masked (ApproxSolve): entries outside the backend's declared support are zeroed in the diff; the omitted-entry magnitudes are reported in max_abs_dev_outside_support but not asserted. Use is_approximate(spec.backend) to confirm which mode you are in.

verify_jacobian_at reports deviations; it does not throw. You inspect report.passed and the per-entry fields to localize a discrepancy.

11.7.2 In-solve verification callback

To sample Jacobian agreement at several times during an actual solve, pass the verify_jacobian keyword as a NamedTuple with keys drawn from {:rtol, :atol, :sample_times}. This installs a callback that compares the active backend against the dense reference each time the integrator crosses one of the listed times and warns (does not abort) on a mismatch:

spec = JacobianSpec(Structured(strategy = StructuredSolve()), problem; use_ale = true)
sol  = solve(problem;
    jacobian        = spec,
    use_ale         = true,
    verify_jacobian = (rtol = 1e-6, atol = 1e-10, sample_times = [0.0, 1.0, 10.0]))

The callback only runs when a pluggable jacobian is supplied (there is nothing to verify against the integrator's own finite-difference Jacobian). It warns rather than throws so a single transient outlier does not abort an otherwise healthy production run.

11.8 The pure-residual layer: building a custom solver loop

Beneath solve, Pyrolysis.jl exposes a decoupled residual API so you can drive the time integration yourself, build a non-standard solver, or feed the residual to external tooling. The three pieces are all top-level exports:

• ProblemDef — an immutable frozen problem definition (material, mesh topology, BCs, initial conditions, simulation-mode traits, options). Safe to share across threads and sweep workers.
• Workspace — a mutable bundle of per-step caches, work arrays, and the geometry-from-state machinery. Allocate one per worker.
• residual!(du, u, def, t, ws) — the pure ODE right-hand side. Note the argument order: du, u, def, t, ws.

The builders to_problem_def, build_workspace, and create_state_vector turn a PyrolysisProblem into this lower-level form:

using Pyrolysis
import OrdinaryDiffEq                                  # from your environment

# 1. Freeze the problem definition (pick the mode flags you want).
def = to_problem_def(problem;
    use_ale         = false,
    radiation_model = NO_RADIATION,
    energy_form     = :standard)

# 2. Build a mutable workspace (one per worker/thread).
ws = build_workspace(problem, def)

# 3. Initial state vector, packed [T, ξ, z (ALE), χ (WithChi)].
u0 = create_state_vector(ws.layout, ws.mesh)

# 4. Wrap residual! in the (du, u, p, t) shape OrdinaryDiffEq expects.
#    Note the residual! argument order is (du, u, def, t, ws).
rhs!(du, u, p, t) = residual!(du, u, def, t, ws)

# 5. Hand-roll your own integration.
prob = OrdinaryDiffEq.ODEProblem(rhs!, u0, problem.tspan)
raw  = OrdinaryDiffEq.solve(prob, OrdinaryDiffEq.KenCarp4(autodiff = OrdinaryDiffEq.AutoFiniteDiff()))

StateLayout (also exported) describes the packing and provides accessors so you can read fields out of a raw state vector without re-deriving the offsets:

layout = ws.layout

# Index helpers and block accessors:
state_length(layout)                       # total state-vector length
get_temperature(u, layout, i)              # T in cell i
get_concentration(u, layout, i, j)         # ξ of component j in cell i (cell index first)
get_node_positions(u, layout)              # z nodes (ALE only)
get_thickness_from_state(u, layout)        # current sample thickness
block_view(u, layout, ...)                 # a view of one block

# For variable-cross-section runs:
column_chi_bar_from_state(u, layout, material)        # column-average χ̄
cross_section_area_from_state(mesh, material, u, layout)

You can also attach a structured Jacobian to your hand-rolled loop. Build a spec from the same (def, ws) pair and call compute_jacobian!:

spec = JacobianSpec(Structured(strategy = DirectSolve()), def, ws)
J    = zeros(state_length(layout), state_length(layout))
compute_jacobian!(J, spec, u0, def, 0.0, ws)

When driving the solver yourself, remember the contract that solve normally enforces for you: the residual and Jacobian assembly are mesh-read-only. The mutable mesh state is only meant to be synchronized to the integrator's u between accepted steps (this is what the built-in StepUpdateCallback does inside solve). If your custom loop relies on the mesh geometry being current for post-processing or AMR, you must call update_mesh_from_state!(mesh, u, layout) between steps yourself.

11.9 The Internal umbrella and importing trait subtypes

The top-level Pyrolysis namespace deliberately exports only a curated public surface. Two escape hatches give you the rest.

11.9.1 Pyrolysis.Internal

Pyrolysis.Internal re-exports every symbol exported by the core submodules (Materials, BoundaryConditions, Geometry, Properties, Physics, Adaptivity, Discretization, Problem, Residual, Jacobian, Diagnostics, Solver) in one namespace:

using Pyrolysis.Internal   # everything the core submodules export, in one place

Use this when you need an internal building block — a mixing function, a flux kernel, a cache type — that is not on the curated top-level list. Caveat: symbols reached through Internal are not part of the stable public API and may move or change between minor releases; prefer top-level exports in production code. (The Experimental namespace — P1 radiation, h-AMR — is not included in Internal; reach it explicitly as Pyrolysis.Experimental.X.)

11.9.2 Importing the simulation-mode trait subtypes

The SimulationMode abstract trait is exported, but the concrete subtypes — Eulerian, ALE, Fickian, DarcyFick, NoChi, WithChi, NoRadiation, SurfaceAbsorption, BeerLambert — are not re-exported at top level, to keep the public surface clean and avoid shadowing generic names like ALE. If you need to dispatch on or construct one (e.g. for a custom residual specialization), reach it through the Problem submodule:

# Fully qualified:
mode_tag = Pyrolysis.Problem.Eulerian()

# Or import the ones you need explicitly:
using Pyrolysis.Problem: Eulerian, ALE, Fickian, DarcyFick, WithChi

In normal usage you never construct these directly — solve and to_problem_def select the right mode bundle from your keyword flags (use_ale, radiation_model, the presence of permeability for DarcyFick, a lateral-shrinkage law for WithChi). You only reach for the trait subtypes when building custom dispatch at the residual layer (see Chapter 7 for how the trait axes map to the four residual paths, and Technical Reference §15.1).

Note: P1 radiation (P1_QUASI_STEADY) is not part of the unified residual; it lives in Pyrolysis.Experimental with its own residual wrapper and is rejected by solve (radiation_model = P1_QUASI_STEADY raises an error). The supported radiation models for solve are NO_RADIATION, SURFACE_ABSORPTION, and BEER_LAMBERT (Technical Reference §8.12).

11.10 Quick configuration recipes

# (a) Default, simplest correct solve — Eulerian, no depletion.
sol = solve(problem)

# (b) Explicit exact structured Jacobian (Eulerian).
spec = JacobianSpec(Structured(strategy = DirectSolve()), problem)
sol  = solve(problem; jacobian = spec)

# (c) Large ALE, fast approximate Jacobian (~8× faster assembly).
spec = JacobianSpec(Structured(strategy = ApproxSolve()), problem; use_ale = true)
sol  = solve(problem; jacobian = spec, use_ale = true)

# (d) Large ALE, exact but cheap assembly via Richardson.
spec = JacobianSpec(Structured(strategy = StructuredSolve()), problem; use_ale = true)
sol  = solve(problem; jacobian = spec, use_ale = true)

# (e) ILU-GMRES linear solve for a large W matrix.
sol = solve(problem;
    integrator = KenCarp4(linsolve = ILUGMRESFactorization(τ = 0.1)),
    jacobian   = JacobianSpec(Structured(strategy = DirectSolve()), problem; use_ale = true),
    use_ale    = true)

# (f) Depletion tracking — no pluggable Jacobian; solve picks the resize-safe dense LU.
sol = solve(problem;
    use_ale          = true,
    handle_depletion = true)

# (g) Vendor solver via extension.
using KLU
sol = solve(problem; linear_solver = Pyrolysis.klu_factorization())

# (h) Verify a custom backend against the dense reference during a solve.
spec = JacobianSpec(Structured(strategy = StructuredSolve()), problem; use_ale = true)
sol  = solve(problem; jacobian = spec, use_ale = true,
             verify_jacobian = (rtol = 1e-6, atol = 1e-10, sample_times = [0.0, 5.0]))

11.11 Cross-references

• Structured-Jacobian architecture, the implicit-RK W = I/(γ_RK Δt) − J matrix, per-operator Jacobian formulas, and the prefix/suffix/rank-1 couplings: Technical Reference §15.
• Surface-depletion handling and the depletion ↔ pluggable-Jacobian incompatibility: Chapter 10 and Technical Reference §14.7.
• Building a problem and selecting simulation modes: Chapter 7.
• The solve keyword reference and post-processing: Chapter 8 and Chapter 9.
• Sensitivity workflows that build on the sensitivity extension: Chapter 13 and Technical Reference §17.
• The complete API and default-value tables: Chapter 16.