Solving

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

Solve a pyrolysis problem using DifferentialEquations.jl.

Arguments

• problem: The PyrolysisProblem to solve. solve does not modify it: the run integrates a copy with its own mesh, so every solve of a problem, and every parameter-sweep point that reuses its mesh, starts from the same initial geometry. solution.problem is that copy; its mesh holds the run's final geometry, live-cell count and boundary faces.

Keyword Arguments

• integrator=nothing: ODE integrator constructor (e.g. KenCarp4) or fully constructed algorithm (e.g. KenCarp4(linsolve=KLUFactorization())). nothing chooses KenCarp4 with the default linear_solver.
• linear_solver=nothing: LinearSolve backend. nothing chooses by the Jacobian pattern (see "Linear Solver"): UMFPACKFactorization() for the exact ALE Jacobian, KLUFactorization() for every other Jacobian. Pair with integrator only when integrator is a constructor.
• jacobian=nothing: Jacobian backend. nothing picks Structured(strategy = DirectSolve()) on a fixed mesh and Structured(strategy = ApproxSolve()) under ALE. A Structured backend such as Structured(strategy = DirectSolve()) overrides the default; the spec is built for this solve's own mesh copy, mode and options.
• 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
• 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 live thickness before termination (ALE only) [m]
• handle_depletion=false: Park depleted cells (see "Cell depletion" below)
• depletion_threshold=0.05: Depletion trigger. ALE: a cell is parked when its thickness falls below this fraction of its own initial thickness. Fixed mesh: the top live cell is parked when its condensed concentration falls below this fraction of its initial value.
• min_cells=2: Minimum number of live cells; the solve terminates instead of parking below it
• warm_start=nothing: Reuse a prior PyrolysisSolution's final state as the initial state (u0). Component count, mesh size, and ALE/χ modes must match, and the prior solve must not have parked cells. 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 = backend to override the auto-default; the spec is built for this solve's mode. The 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, in either mesh mode, the solver monitors every live cell. Under ALE a cell that has shrunk below depletion_threshold of its own initial thickness (a per-cell reference, so stretched meshes are judged cell-by-cell) is parked; on a fixed mesh the top live cell is parked once its condensed concentration has fallen below the same fraction of its initial value. A fixed-mesh solve stops with a warning (retcode :Terminated) when a cell below the top live cell becomes a void: the material has gasified from within, the model carries no heat across the void, and a regressing material belongs on the ALE mesh. A park merges the cell's mass and sensible enthalpy into its neighbour and moves the freed state slot above the live surface as an inert ghost: the state vector, mesh arrays and Jacobian keep their size, every kernel runs over the live range 1:mesh.n_active_cells, and the boundary conditions act at the live surface. Vanishing cells therefore never enter the nonlinear system (their 1/Δz² conditioning is what stalls an implicit solve), and the structured Jacobian stays valid through burnout. Ghost slots read as NaN in solution.T, solution.ξ and solution.z; solution.extras.live_cells and solution.extras.park_times record the events. Parameter sensitivities are valid only on horizons without park events.

Recommended Solvers

The default is KenCarp4(autodiff = AutoFiniteDiff()) with the default Jacobian (see "Default Jacobian") and the default linear solver (see "Linear Solver"). 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), a fixed mesh uses Structured(strategy = DirectSolve()) and ALE uses Structured(strategy = ApproxSolve()). The Structured backend assembles per-operator contributions into a StructuredJacobian. DirectSolve densifies the structured couplings (ALE prefix-sum, χ rank-1, Beer–Lambert suffix-sum) at solve time, so the linear solve sees the exact Jacobian. ApproxSolve drops them: under ALE the exact pattern couples every node to every cell below it, and the reduced one keeps the Newton matrix local while the modified Newton iteration still converges. Pass jacobian = Structured(strategy = DirectSolve()) for the exact ALE Jacobian.

Linear Solver

With linear_solver = nothing the sparse factorization follows the Jacobian pattern. Structured(strategy = DirectSolve()) under ALE materializes the geometry block, in which every cell couples to every node below it; that pattern uses UMFPACKFactorization(), whose supernodal factorization handles the dense block. Every other sparse Jacobian — every ApproxSolve Jacobian, and DirectSolve on an Eulerian mesh, Beer–Lambert suffix included — uses KLUFactorization(), which refactors the fixed pattern in place without allocating. Pass linear_solver = SparspakFactorization(), another LinearSolve backend, or a fully constructed integrator to choose the factorization.

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

# Choose the factorization instead of the pattern-based default:
solution = solve(problem, use_ale=true, linear_solver=SparspakFactorization())

# The exact ALE Jacobian instead of the reduced default (factored by UMFPACK):
solution = solve(problem, use_ale=true, jacobian=Structured(strategy=DirectSolve()))

# With depletion handling (cells park at 1% of their initial thickness):
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 new problem may reuse the base problem's mesh: solve works on its own copy of the mesh, so every point starts from the base geometry. 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.
	mat = prob.material
	new_rxn = Reaction(:decomp, 1 => 2; A = p.A, E = p.E, h = 200e3)
	new_mat = Material(name = mat.name, components = mat.components,
	                   reactions = (new_rxn,), mixing = mat.mixing,
	                   lateral_shrinkage_law = mat.lateral_shrinkage_law)
	return PyrolysisProblem(mesh = prob.mesh, material = new_mat,
	                        bc_set = prob.bc_set, T_initial = prob.T_initial,
	                        ξ_initial = prob.ξ_initial, tspan = prob.tspan,
	                        dt_initial = prob.dt_initial,
	                        output_times = prob.output_times)
end

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

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