Solving

solve(problem::PyrolysisProblem; kwargs...) -> PyrolysisSolution

Solve a pyrolysis problem using DifferentialEquations.jl.

Arguments

• problem: The PyrolysisProblem to solve

Keyword Arguments

• integrator=nothing: ODE integrator constructor (e.g. KenCarp4) or fully constructed algorithm (e.g. KenCarp4(linsolve=SparspakFactorization())). nothing chooses KenCarp4 with the default linear_solver.
• linear_solver=nothing: LinearSolve backend. nothing chooses SparspakFactorization() for sparse Jacobians. Pair with integrator only when integrator is a constructor.
• jacobian=nothing: Jacobian handle (JacobianSpec). nothing picks Structured(strategy = DirectSolve()). Pass an explicit JacobianSpec(Structured(strategy = ApproxSolve()), problem) for the Schur-style approximation that drops the structured ALE couplings.
• verify_jacobian=nothing: NamedTuple (rtol=, atol=, sample_times=) installs a DiscreteCallback that compares the active backend against a DenseAD_ForwardDiff reference at the listed sample times. Reports via @warn; does not throw.
• abstol=1e-8: Absolute tolerance
• reltol=1e-6: Relative tolerance
• sensealg=nothing: Optional SciMLSensitivity sense algorithm forwarded to OrdinaryDiffEq.solve. This is primarily used by PyrolysisSciMLSensitivityExt.
• maxiters=1000000: Maximum iterations
• dtmin=1e-12: Minimum timestep
• dtmax=Inf: Maximum timestep
• saveat=problem.output_times: Times to save output
• callback=nothing: Additional callbacks
• progress=true: Show progress bar during solve
• validate=true: Validate problem before solving
• use_ale=false: Enable ALE mesh motion for shrinking/swelling materials
• min_thickness=1e-6: Minimum thickness before termination (ALE only) [m]
• handle_depletion=false: Enable surface cell depletion handling (ALE only)
• depletion_threshold=0.05: Thickness fraction threshold for cell merging (0.05 = 5% of the cell's own initial thickness)
• min_cells=2: Minimum number of cells before terminating instead of merging
• use_transpiration_bc=false: Opt-in surface gas sensible-enthalpy transpiration. Off by default — the volumetric advection source Sconv already carries the gas enthalpy (ThermaKin Eq.17 convention); enabling this alongside Sconv double-counts and is rejected at problem construction. Retained only for closed-ledger experiments with S_conv disabled.
• use_blowing_correction=false: Opt-in Couette-film blowing correction to every convective BC in the surface energy balance: h → h·B/(e^B − 1) with B = ṁ″·c_p,g/h built from the outward gas mass flux at the face. The momentum-side companion of the transpiration term — unlike it, fully compatible with the Sconv-only convention (no enthalpy double-count). Suppresses convective cooling at peak burning (B = 4 → heff/h = 0.075); the back face is unaffected automatically (no outflow → B = 0). Calibrated examples should enable this only together with flame heat-flux feedback (see the physics roadmap).
• warm_start=nothing: Reuse a prior PyrolysisSolution's final state as the initial state (u0). Component count, mesh size, and ALE/χ modes must match. tspan[1] still defines the integrator's start time, so the warm-started solve treats the prior final state as fresh initial data.

Jacobian backend

Pass jacobian = JacobianSpec(backend, problem) to override the auto-default. The production backend is Structured; behavior is selected by its strategy:

Returns

A PyrolysisSolution containing:

• Temperature and concentration profiles at saved times
• Mass loss rate history
• Surface temperature history
• Thickness history (when use_ale=true)
• Solver return code and timing

ALE Mode

When use_ale=true, the solver uses the unified ProblemDef/Workspace residual. Node positions evolve simultaneously with temperature and concentration fields, maintaining proper conservation (GCL compliance). The thickness history is extracted from the state vector and stored in solution.extras.thickness_history.

Cell Depletion Handling

When handle_depletion=true (requires use_ale=true), the solver monitors every cell's thickness. When a cell shrinks below depletion_threshold fraction of its own initial thickness (a per-cell reference, so stretched meshes are judged cell-by-cell), it is merged with a neighbor. This prevents the ρcₚ→0 singularity that causes temperature blow-up in materials that fully decompose. The mesh is dynamically resized during integration.

Recommended Solvers

The default is KenCarp4(autodiff = AutoFiniteDiff(), linsolve = SparspakFactorization()) with Structured(strategy = DirectSolve()). This is the production baseline for stiff pyrolysis runs.

Other stiff ODE algorithms such as KenCarp3, TRBDF2, Rodas5P, and QNDF can be passed through integrator. Rosenbrock and BDF methods may reevaluate or refactor the Jacobian more often on these problems, so benchmark against the default before changing the solver family. Note: under OrdinaryDiffEq v7 only the default solver set is re-exported by the umbrella package, so passing e.g. QNDF/KenCarp3/TRBDF2 requires loading the relevant sublibrary (OrdinaryDiffEqBDF, OrdinaryDiffEqSDIRK, …) yourself.

ILUGMRESFactorization is available for large-ALE experiments where direct factorization cost dominates, but it is not the default because inexact linear solves can increase Newton rejections on strongly Arrhenius transients.

Default Jacobian

With jacobian = nothing (the default), every mode uses Structured(strategy = DirectSolve()). The Structured backend assembles per-operator contributions into a StructuredJacobian and densifies any structured couplings (ALE prefix-sum, χ rank-1, Beer- Lambert suffix-sum) at solve time so the linear-solve sees a true sparse Jacobian. For ALE problems where the dense cumulative C block dominates factorization cost, opt into the Schur-style approximation via Structured(strategy = ApproxSolve()) — this drops the structured couplings, at the cost of an approximate Newton search direction.

Linear Solver

Sparspak is the default sparse direct factorization for pyrolysis problems. Pass linear_solver = KLUFactorization(), another LinearSolve backend, or a fully constructed integrator when you need to compare factorization choices.

Example

problem = PyrolysisProblem(mesh, material, bcs, (0.0, 100.0))
solution = solve(problem)
using OrdinaryDiffEqBDF  # QNDF lives in a sublibrary under OrdinaryDiffEq v7
solution = solve(problem, integrator=QNDF, abstol=1e-10)
solution = solve(problem, progress=false)  # Disable progress bar

# With ALE mesh motion for shrinking materials:
solution = solve(problem, use_ale=true, min_thickness=0.1e-3)
thickness_history = solution.extras.thickness_history

# With ALE — auto-default Jacobian is Structured(strategy=DirectSolve()):
solution = solve(problem, use_ale=true,
    integrator=KenCarp4(linsolve=SparspakFactorization()))

# Opt into the ALE Schur-style approximation:
solution = solve(problem, use_ale=true,
    jacobian=JacobianSpec(Structured(strategy=ApproxSolve()), problem))

# With surface cell depletion handling (prevents blow-up):
solution = solve(problem, use_ale=true, handle_depletion=true, depletion_threshold=0.01)

# Differential-correctness verification at sample times:
solution = solve(problem,
    verify_jacobian=(rtol=1e-6, atol=1e-10, sample_times=[0.0, 1.0, 10.0]))

parameter_sweep(problem, grid;
                modify::Function = default_modify,
                extract::Function = identity,
                parallel::Bool = nworkers() > 1,
                solve_kwargs...) -> Vector

Sweep problem over grid and return a vector of results, one per grid point. Each grid point p ∈ grid is materialised into a fresh problem via modify(problem, p), solved with solve_kwargs, then reduced via extract(solution) before being returned.

Arguments

• problem::PyrolysisProblem: The base problem. Not mutated.
• grid: Any iterable of parameter points. Common shapes: Vector{<:NamedTuple}, Iterators.product(...), Vector{Vector{Float64}}.
• modify: (base_problem, point) -> new_problem. Must return a fresh PyrolysisProblem for each call; do not mutate base_problem. The default raises a helpful error pointing to this kwarg.
• extract: solution -> summary. Defaults to identity. For large grids prefer a small NamedTuple — full PyrolysisSolution instances are expensive to ship across worker boundaries.
• parallel: Use Distributed.pmap when true, otherwise fall back to a serial map. Defaults to parallel when more than one worker is available (nworkers() > 1).
• solve_kwargs: Forwarded to solve().

Notes

• Use Distributed.addprocs() (and @everywhere using Pyrolysis) before calling this function to make workers available.
• Failed solves are not retried. Wrap the body of extract in a try if you need fault tolerance.

Example

using Distributed
addprocs(4)
@everywhere using Pyrolysis

grid = [(A = a, E = e) for a in (1e12, 1e13, 1e14), e in (1.5e5, 2.0e5)]

function modify(prob, p)
	# Rebuild the material with the swept Arrhenius parameters.
	new_rxn = Reaction(:decomp, 1 => 2; A = p.A, E = p.E, h = 200e3)
	new_mat = Material(prob.material; reactions = (new_rxn,))
	return PyrolysisProblem(prob; material = new_mat)
end

extract(sol) = (peak_T = maximum(sol.surface_T),
                final_mass = sol.extras.total_mass[end])

results = parameter_sweep(problem, grid; modify, extract)