16. API and Parameter Reference

This chapter is the consolidated reference for the public API surface of Pyrolysis.jl and every user-facing parameter and its default. Earlier chapters teach workflow; this one is the lookup table you return to when you need the exact constructor keyword, the precise default value, or the right place to import a name from. Every type, function, keyword, and default shown here was checked against the package source — where a name is reachable only through a submodule (and not from using Pyrolysis alone), that is stated explicitly so you never get an UndefVarError from a copied snippet. For the physics behind any quantity, follow the "see Technical Reference §N" pointers.

Two coordinate/convention reminders that the whole API assumes (see Technical Reference §2 and §3):

• z = 0 is the bottom / substrate (boundary id 2, tag :bottom); z = L is the top / exposed surface (boundary id 1, tag :top). Heat enters from the top. Boundary heat and mass fluxes are positive into the domain.
• The primary species state is the mass concentration ξ_j in kg·m⁻³, not a mass fraction. The heat-of-reaction sign convention is h > 0 endothermic (cools), h < 0 exothermic (heats) — see Technical Reference §6.

16.1 What using Pyrolysis gives you, and what it does not

The top-level Pyrolysis module exports a curated surface. Many useful names live in submodules and are not re-exported at the top level. There are three ways to reach a name:

1. Top-level export — available after using Pyrolysis. This is the stable public API; prefer it in production code.

2. Submodule — e.g. Pyrolysis.Discretization.get_surface_temperature, or import it: using Pyrolysis.Discretization: get_surface_temperature.

3. Pyrolysis.Internal — an umbrella namespace that re-exports every submodule's public symbol at once:

   using Pyrolysis            # curated top-level surface
   using Pyrolysis.Internal   # everything every submodule exports

   Symbols reached via Internal are not part of the stable API and may move between minor releases. Pyrolysis.Experimental is not included in Internal; access it explicitly (Pyrolysis.Experimental.X).

Caveat that bites copy-pasters: the solution post-processors get_profile, get_cell_positions, interpolate_solution; the callbacks create_depletion_callback, create_thickness_termination_callback, create_step_update_callback; the iterative linear solvers ILUGMRESFactorization, ILUBiCGSTABFactorization; the surface-temperature accessors get_surface_temperature, get_back_face_temperature; and the geometry accessors (n_cells, cell_center, total_volume, …) are exported by their submodules but not at the top level. Use using Pyrolysis.Internal or a submodule-qualified import to call them.

16.1.1 Submodule map

The package is a layered DAG of submodules. The table lists each submodule and the kind of names it owns. Everything a submodule exports is reachable as Pyrolysis.<Submodule>.X and is bundled into Pyrolysis.Internal.

16.1.2 Extension stubs

Three families of functions are stubs in core Pyrolysis that gain concrete methods only when an optional dependency is loaded. Calling a stub without its extension raises an error telling you which package to add.

Plotting recipes (plot(sol; kind = …)) attach when using Plots (or RecipesBase) is loaded; see Chapter 9.

16.2 Materials, components, and property functions

See Technical Reference §4 (materials and property functions) and §5 (effective properties). All names in this section are top-level exported.

16.2.1 Components

Three concrete component types encode phase at the type level. Each ships a keyword convenience constructor.

SolidComponent(name::Symbol; ρ, c, k,
               ε = 0.9, α = 0.0, γ = 1.0, λ = 0.0, κ = 0.0,
               ρ_intrinsic = nothing, φ = nothing)

LiquidComponent(name::Symbol; ρ, c, k,
                λ = 2e-5, κ = 0.0, ε = 0.95, α = 0.0, γ = 0.0,
                ρ_intrinsic = nothing, φ = nothing)

GaseousComponent(name::Symbol; M, c, k, λ, ε = 0.0, α = 0.0, γ = 0.0)

ρ, c, k, λ, κ accept either a plain number or a property-function object (§16.2.2); they are converted with to_property automatically. Pass either ρ_intrinsic or φ, never both — supplying both throws an ArgumentError. The bulk density ρ is what mixing rules and mixture density use; ρ_intrinsic is used only in the porosity formula φ = clamp(1 − Σ ξ_j/ρ_{i,j}, 0, 1) (see §5).

virgin = SolidComponent(:virgin, ρ = 530, c = 2000, k = 0.12, ε = 0.9, κ = 1e-15)
char   = SolidComponent(:char,   ρ = 150, c = 1500, k = 0.20, ε = 0.95, α = 5.0, κ = 1e-12)
water  = LiquidComponent(:H2O_l, ρ = 1000, c = 4180, k = 0.6)
gas    = GaseousComponent(:gas,  M = 0.044, c = 1800, k = 0.03, λ = 1e-4)
vapor  = GaseousComponent(:H2O,  M = 0.018, c = 2000, k = 0.025, λ = 2e-5)

Phase predicates (top-level exported): is_solid(c), is_liquid(c), is_gas(c), phase(c), and molar_mass(c) (returns 0.0 for non-gas).

16.2.2 Property-function types

Every property field can be a constant number or one of these callable property objects (all top-level exported). They implement p(T) (and p(T, ξ) for the state-dependent one) and carry analytical or finite-difference derivatives used in Jacobian assembly (see §4.5).

TableProperty requires strictly increasing temperatures (checked at construction); it uses a linear search for ≤8 points and a binary search above that. StateDependentProperty is the type to use for moisture-dependent heat of sorption and similar composition-dependent quantities.

The convenience converter to_property (used implicitly by the constructors) maps Real → ConstantProperty, Tuple → LinearProperty, Function → CallableProperty, and passes property objects through.

# Temperature-dependent conductivity that rises with T:
char = SolidComponent(:char, ρ = 150, c = 1500, k = LinearProperty(0.10, 1.5e-4))

# Tabulated heat capacity:
cp_tab = TableProperty((300.0, 500.0, 800.0), (1500.0, 1800.0, 2100.0))
virgin = SolidComponent(:virgin, ρ = 530, c = cp_tab, k = 0.12)

16.2.3 Reactions

Reaction (top-level exported) is a single-reactant Arrhenius reaction. Four constructor forms cover numeric vs. symbolic component references and one vs. two products. Symbolic forms return a PendingReaction that the Material constructor resolves against the component names.

# one product, numeric indices
Reaction(name, reactant::Int => product::Int;
         A, E, h, n = 1.0, T_min = 0.0, T_max = Inf, validate_mass = true)

# two products, numeric indices
Reaction(name, reactant::Int => (p1::Int, p2::Int);
         A, E, h, yields::Tuple, n = 1.0, T_min = 0.0, T_max = Inf, validate_mass = true)

# one product, symbolic names
Reaction(name, reactant::Symbol => product::Symbol; A, E, h, n = 1.0, …)

# two products, symbolic names
Reaction(name, reactant::Symbol => (p1::Symbol, p2::Symbol); A, E, h, yields, n = 1.0, …)

Mass-balance convention: the single reactant is consumed with implicit mass 1.0; product yields must sum to 1.0 (the one-product form uses yield 1.0 automatically). Each yield must additionally be ≥ 0 and the order must satisfy n ≥ 0 — both enforced unconditionally at construction; validate_mass = false disables only the sum check. h accepts a number (ConstantProperty), a (a, b) tuple (LinearProperty, i.e. a + b·T), or any property object — use StateDependentProperty for moisture sorption. See §6 for the rate law, tanh gates (≈10 K leakage caveat), and the depletion limiter.

# endothermic charring: virgin -> 0.2 char + 0.8 gas
rxn = Reaction(:charring, :virgin => (:char, :gas);
               A = 1e15, E = 200e3, h = 500e3, yields = (0.2, 0.8),
               T_min = 473.0, T_max = 1273.0)

# exothermic char oxidation (negative h releases heat)
rxn2 = Reaction(:char_oxidation, 1 => 2; A = 1e10, E = 160e3, h = -25e6)

ReactionSet is a tuple-based collection type (top-level exported); the common path is to pass a plain tuple of reactions to Material.

16.2.4 Mixing rules: MixingRule, MixingSpec, MaterialMixing

Effective conductivity k, gas transfer λ, and permeability κ are each mixed with their own rule (see §5). The rule enum (top-level exported) is:

@enum MixingRule PARALLEL SERIES WEIGHTED BRUGGEMAN CARMAN_KOZENY

MixingSpec bundles a rule with the weighting β and (for Carman–Kozeny) the characteristic length:

MixingSpec(rule::MixingRule, β::Real = 0.5; length_scale::Real = NaN)

β must be in [0, 1] (default 0.5). For CARMAN_KOZENY you must pass a positive length_scale (the characteristic pore/particle diameter ℓ in metres) via the keyword — note it is length_scale, not a positional ℓ.

MaterialMixing builds the (k, λ, κ) NamedTuple that Material expects. Each channel accepts a MixingSpec or a bare MixingRule (in which case β = 0.5):

MaterialMixing(; k = MixingSpec(WEIGHTED, 0.5),
                 λ = MixingSpec(PARALLEL, 0.5),
                 κ = MixingSpec(WEIGHTED, 0.5))

mixing = MaterialMixing(
    k = MixingSpec(WEIGHTED, 0.6),
    λ = PARALLEL,                              # bare rule -> β = 0.5
    κ = MixingSpec(CARMAN_KOZENY; length_scale = 5e-5),
)

16.2.5 Material

Material (top-level exported) assembles components, reactions, mixing, and an optional lateral-shrinkage law.

Material(; name::Symbol,
           components::Tuple,
           reactions::Tuple = (),
           mixing = MaterialMixing(),
           lateral_shrinkage_law::Union{Function,Nothing} = nothing,
           depletion_limiter::DepletionLimiter = DepletionLimiter())

Providing a non-nothing lateral_shrinkage_law switches the simulation into the variable-cross-section (WithChi) mode under ALE (see §10 and §16.6). The law is any function of the column-average pyrolysis progress χ̄, e.g. χ̄ -> 1 - 0.3*χ̄ for 30 % lateral shrink at full conversion.

pmma = Material(
    name = :PMMA,
    components = (
        SolidComponent(:virgin, ρ = 1190, c = 1420, k = 0.21, ε = 0.86),
        GaseousComponent(:MMA, M = 0.100, c = 1500, k = 0.02, λ = 1e-5),
    ),
    reactions = (
        Reaction(:decomposition, :virgin => :MMA; A = 8.5e12, E = 188e3, h = 870e3),
    ),
    mixing = MaterialMixing(k = MixingSpec(WEIGHTED, 0.5)),
)

DepletionLimiter(; threshold = 1.0, enabled = true) — threshold [kg·m⁻³] is the concentration scale of the smooth tanh(ξ/threshold) rate roll-off applied per reactant; enabled = false removes the only guard that stops zero-order reactions at ξ = 0. Exported by the Physics and Materials submodules (not at the top level): using Pyrolysis.Physics: DepletionLimiter.

Pyrolysis.jl does not provide a "copy-with-overrides" Material(prob.material; …) or PyrolysisProblem(prob; …) constructor. To vary a parameter (e.g. in a sweep), rebuild the full object from scratch — see §16.9.

16.2.6 RadiationModel enum

The volumetric radiation absorption model (top-level exported, except the P1 value) is an enum:

@enum RadiationModel NO_RADIATION SURFACE_ABSORPTION BEER_LAMBERT P1_QUASI_STEADY

P1_QUASI_STEADY is not wired into the core solve path; passing it to solve raises an error (see §16.7 and §8). It lives in Pyrolysis.Experimental with its own residual.

16.3 Meshes and geometry

See Technical Reference §13 (discretization) and the geometry notes in §13. Only the two mesh creators and the topology/geometry view types are top-level exported; the geometry accessors require Pyrolysis.Internal or the Geometry submodule.

16.3.1 Mesh creators (top-level exported)

create_uniform_mesh(L::Real, n_cells::Int, ::Val{NC};
                    T_initial = 300.0, ξ_initial = nothing,
                    cross_section_area = 1.0)

create_stretched_mesh(L::Real, n_cells::Int, ::Val{NC}, stretch::Real;
                      T_initial = 300.0, ξ_initial = nothing,
                      cross_section_area = 1.0, stretch_from = :top)

The third argument is a Val: Val(NC) where NC is the component count — e.g. create_uniform_mesh(0.01, 100, Val(2)) for a 1 cm, 100-cell, 2-component mesh. The initial conditions seeded here are usually overwritten by the problem's T_initial/ξ_initial closures, so you typically leave T_initial/ξ_initial at their defaults and set them on PyrolysisProblem.

mesh   = create_uniform_mesh(0.01, 100, Val(2))                 # 1 cm, 100 cells
meshst = create_stretched_mesh(0.02, 80, Val(3), 1.03; stretch_from = :top)

16.3.2 Geometry accessors (not top-level; use Internal or Geometry)

After a solve you usually read geometry through the PyrolysisSolution (§16.5). When you need to inspect the mesh directly, import from the submodule:

using Pyrolysis.Geometry: n_cells, n_nodes, cell_center, cell_size,
                          total_volume, domain_bounds, boundary_face_ids

n = n_cells(mesh)               # number of cells
zc = [cell_center(mesh, i) for i in 1:n]
(zmin, zmax) = domain_bounds(mesh)

MeshTopology and MeshGeometry (the immutable-connectivity and mutable-geometry views used by the low-level residual layer) are top-level exported. The component count of a material is n_components(material), exported from Pyrolysis.Materials (also in Internal).

16.3.3 State-vector layout

The flat ODE state is packed block-major as [T | ξ_1..ξ_NC | z | χ], where the z block exists only in ALE mode and the χ block only in WithChi mode (see §15 and §16.6). StateLayout (top-level exported) is the descriptor; state_length, pack_state!, create_state_vector, update_mesh_from_state!, block_view, and the accessors get_temperature, get_concentration, get_node_position, get_node_positions, get_thickness_from_state, column_chi_bar_from_state, cross_section_area_from_state are also top-level exported for users who index raw state. Most users never touch these — solve and the solution object handle packing for you. See §16.11 for the low-level residual surface.

16.4 Boundary and initial conditions

See Technical Reference §12. All concrete BC types listed here and BoundaryConditionSet are top-level exported. The P1 radiation BCs (Marshak, Reflective, Dirichlet) are not; they live in Pyrolysis.Experimental.

16.4.1 Thermal BCs

Notes and defaults:

• ConvectiveBC requires h and T_inf keywords; both accept a number or a function of t.
• RadiativeBC: ε = nothing (default) uses the composition-weighted effective emissivity ε_eff (Kirchhoff), tracking the surface as it chars; pass a number or T -> ε to override. F is the view factor, default 1.0.
• RadiativeFluxBC is incident external radiation (e.g. a cone heater), not the surface's own re-radiation. absorptivity = nothing (default) uses ε_eff (Kirchhoff). Its behaviour depends on the radiation_model passed to solve: for NO_RADIATION/SURFACE_ABSORPTION it is applied as a surface flux; for BEER_LAMBERT/P1_QUASI_STEADY it returns zero boundary flux because the incident energy is deposited volumetrically (avoids double-counting). The flux can also be a time-varying function; solve re-samples it at every RHS evaluation.
• Combine BCs with +: HeatFluxBC(50e3) + RadiativeBC(ε = 0.9, T_inf = 300.0) + ConvectiveBC(h = 10.0, T_inf = 300.0).

The surface temperature T_s is found by a Newton solve of the surface energy balance q_BC(T_s) = k_eff·(T_s − T_cell)/Δz_{1/2} (clamped to [200, 3000] K, tol 1e-10, ≤50 iterations) — see §12.

16.4.2 Mass BCs

DiffusiveMassBC() defaults to h_m = Inf (no external resistance → pure diffusive limit J = λξ/d); a finite h_m gives the combined-resistance form J = λ h_m ξ/(λ + h_m d). ConvectiveMassBC defaults ξ_inf = 0.0 (escape to clean atmosphere); ξ_inf also accepts t -> ξ_∞(t) or a per-component NTuple (length = NC, entries ≥ 0).

16.4.3 Pressure BCs (Darcy mode)

ConvectivePressureBC takes h_P (pressure-transfer coefficient, m·Pa⁻¹·s⁻¹) and P_ext positionally. Pressure BCs only act when the material enables Darcy flow (any solid/liquid κ > 0); see §9.

16.4.4 Mesh motion (ALE)

There are no mesh-motion BC types. ALE mesh motion is governed entirely by the residual: the bottom node is pinned (w_1 = 0) and node velocities integrate the volumetric strain upward, so the surface follows the material automatically (§7.3; Technical Reference §11.3). Enable it with use_ale = true — nothing to construct.

16.4.5 BoundaryConditionSet

The container that pairs :top and :bottom BCs (top-level exported):

BoundaryConditionSet(; thermal::NamedTuple = (top = AdiabaticBC(), bottom = AdiabaticBC()),
                       mass::NamedTuple    = (top = ImpermeableBC(), bottom = ImpermeableBC()))

User-supplied partial NamedTuples are merged with the defaults, then validated. Tags map to ids :top → 1, :bottom → 2.

bcs = BoundaryConditionSet(
    thermal = (
        top    = RadiativeFluxBC(50e3) + RadiativeBC(ε = 0.9, T_inf = 300.0) +
                 ConvectiveBC(h = 10.0, T_inf = 300.0),
        bottom = AdiabaticBC(),
    ),
    mass = (
        top    = DiffusiveMassBC(),     # gas escapes at the surface
        bottom = ImpermeableBC(),
    ),
)

16.4.6 Initial conditions

Initial conditions are set on PyrolysisProblem (§16.5.1) via the T_initial and ξ_initial closures (functions of z). When ξ_initial is omitted, the default seeds the first solid component at its reference density and all others to zero. T_initial defaults to z -> 300.0. Profile helpers (step_profile, linear_profile, gaussian_profile, tanh_profile) live in Pyrolysis.Problem.

Transpiration caveat. The surface gas-transpiration enthalpy term is incompatible with the default S_conv gas-convection convention (they double-count). solve rejects use_transpiration_bc = true. Leave it off.

16.5 The problem and solution objects

See Technical Reference §3 and §15. PyrolysisProblem and PyrolysisSolution are top-level exported.

16.5.1 PyrolysisProblem

PyrolysisProblem(; mesh, material, bc_set::BoundaryConditionSet,
                   T_initial = z -> 300.0,
                   ξ_initial = nothing,
                   tspan::Tuple{Real,Real},
                   dt_initial::Real = 1e-4,
                   output_times::Vector{<:Real} = Float64[])

# convenience positional form:
PyrolysisProblem(mesh, material, bc_set, tspan; kwargs...)

problem = PyrolysisProblem(
    mesh = mesh,
    material = pmma,
    bc_set = bcs,
    T_initial = z -> 300.0,
    ξ_initial = [z -> 1190.0, z -> 0.0],   # virgin density, gas zero
    tspan = (0.0, 100.0),
)

validate_problem(problem; strict = false) (in Pyrolysis.Problem) checks mesh, material stoichiometry, BC coverage, time span, and initial-condition counts; solve runs it automatically unless validate = false. validate_material(m) (in Pyrolysis.Materials) validates a material in isolation.

16.5.2 PyrolysisSolution

solve returns a PyrolysisSolution. Core fields (all read directly):

extras is a NamedTuple of vectors, written in place during extraction: gas_generation_rate, back_temperature, total_mass, thickness_history, surface_mass_flux, surface_heat_flux_net, surface_heat_flux_absorbed, surface_heat_flux_reradiation, surface_heat_flux_convection, surface_heat_flux_conduction, surface_heat_flux_transpiration, volumetric_absorption, surface_emissivity, surface_conductivity, back_heat_flux, MLR_gauge, cross_section_area_history, diagnostics_log.

sol = solve(problem)
times    = sol.t
Tsurf    = sol.surface_T                       # K, per saved time
mlr_top  = sol.extras.surface_mass_flux        # kg/(m^2 s), top face only
mlrgauge = sol.extras.MLR_gauge                # normalized to A_0 (variable cross-section)
Tprofile = sol.T[:, end]                       # final temperature profile
ξvirgin  = sol.ξ[1, :, end]                    # final virgin concentration profile

MLR_gauge equals mass_loss_rate scaled by A(t)/A_0; for constant cross-section the two coincide.

16.5.3 Solution post-processors (not top-level; use Internal/Solver)

using Pyrolysis.Internal   # or: using Pyrolysis.Solver: get_profile, …

get_profile(sol, field::Symbol, time_idx::Int)   # :T, :ξ, or :ξ1/:ξ2/…
get_cell_positions(sol)                          # cell-center positions from the mesh
interpolate_solution(sol, t::Real)               # full raw state vector at arbitrary t

get_profile(sol, :T, i) returns sol.T[:, i]; get_profile(sol, :ξ, i) returns the full NC × n_cells slice; get_profile(sol, :ξ1, i) returns component 1. The surface/back-face temperature accessors read the mesh's solved boundary state, not the solution:

using Pyrolysis.Discretization: get_surface_temperature, get_back_face_temperature,
                                check_boundary_energy_balance
Ts_final    = get_surface_temperature(problem.mesh)   # mesh holds last accepted step
Tback_final = get_back_face_temperature(problem.mesh)

For per-step surface temperature, prefer sol.surface_T.

16.6 Simulation modes (traits)

See Technical Reference §15. A simulation's behaviour is a bundled trait SimulationMode{Mesh,Rad,Trans,Chi} derived automatically by solve/to_problem_def from your options and material — you rarely construct it by hand. SimulationMode is top-level exported; the trait subtypes are deliberately not exported (to keep the namespace clean). Reach them via Pyrolysis.Problem:

using Pyrolysis.Problem: Eulerian, ALE, Fickian, DarcyFick, NoChi, WithChi,
                         mesh_mode, radiation_mode, transport_mode, chi_mode

Rules of thumb for choosing modes via solve options:

• TGA-like / fixed slab: defaults (Eulerian, Fickian, NO_RADIATION).
• Cone calorimeter with surface recession: use_ale = true, radiation_model = SURFACE_ABSORPTION (opaque) or BEER_LAMBERT (semi-transparent).
• Pressure-driven char / pressure-treated wood: give components κ > 0 ⇒ Darcy–Fick automatically; add pressure BCs.
• Intumescent / lateral shrink: give the material a lateral_shrinkage_law and run with use_ale = true ⇒ WithChi.

16.7 solve: full keyword reference

See Technical Reference §15. solve is top-level exported. It rejects unknown keywords (typo protection), so only the keywords below are valid.

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

Key behaviours and constraints:

• Default integrator / linear solver: KenCarp4(autodiff = AutoFiniteDiff(), linsolve = SparspakFactorization()). KenCarp4 is an SDIRK method suited to the moderately stiff Arrhenius system. Pass a constructor (e.g. integrator = KenCarp4) to combine with linear_solver, or a fully-built instance (e.g. integrator = KenCarp4(linsolve = ILUGMRESFactorization())) on its own.
• Default Jacobian: Structured(strategy = DirectSolve()). Pass a JacobianSpec to opt into other strategies (§16.8).
• radiation_model = P1_QUASI_STEADY is rejected (P1 is experimental).
• energy_form must be :standard (Formulation A, the FDS/Gpyro/ThermaKin convention) or :with_generation_sink (adds the −Σ Δh_g ∇·J_g cell source).
• handle_depletion = true is incompatible with a pluggable jacobian: the state vector resizes mid-solve, which a pre-built JacobianSpec cannot survive. solve errors if both are set; with depletion alone it falls back to the integrator's colored finite-difference Jacobian and switches the linear solver to a dense LU.
• MLR gauge note: in variable-cross-section runs, sol.mass_loss_rate is per current area and sol.extras.MLR_gauge is normalized to A_0.

# Cone calorimeter, opaque, ALE surface recession, default Jacobian:
sol = solve(problem;
            use_ale = true,
            radiation_model = SURFACE_ABSORPTION,
            min_thickness = 0.5e-3,
            abstol = 1e-9, reltol = 1e-7,
            integrator = KenCarp4(linsolve = SparspakFactorization()))

16.8 Jacobian and linear-solver configuration

See Technical Reference §15. These names are top-level exported.

16.8.1 Backends

Build a spec from a problem (the common form) or from a (def, ws) pair:

JacobianSpec(backend::JacobianBackend, problem::PyrolysisProblem; kwargs...)
JacobianSpec(backend::JacobianBackend, def::ProblemDef, ws::Workspace; kwargs...)

spec = JacobianSpec(Structured(strategy = ApproxSolve()), problem)
sol  = solve(problem; jacobian = spec, use_ale = true)

16.8.2 Linear-solve strategies

Structured takes a LinearSolveStrategy (all top-level exported):

ApproxSolve is marked approximate via the is_approximate trait; the verifier compares it in masked (in-support) mode. For the Douglas-fir cone case, Structured(strategy = ApproxSolve()) gives roughly an 8× assembly speedup over the full Jacobian at solver-tolerance-equivalent results.

16.8.3 Iterative linear solvers (not top-level; use Internal/Solver)

using Pyrolysis.Solver: ILUGMRESFactorization, ILUBiCGSTABFactorization

ILUGMRESFactorization(; τ = 0.1, gmres_rtol = 1e-8, restart = 30)
ILUBiCGSTABFactorization(; τ = 0.1, rtol = 1e-8)

sol = solve(problem;
            integrator = KenCarp4(linsolve = ILUGMRESFactorization(τ = 0.1)))

Both fall back to diagonal preconditioning if the ILU factorization is singular during a transient. Use them for large ALE problems where direct factorization dominates cost (trade-off: inexact solves can cause Newton rejections on sharp Arrhenius transients).

16.8.4 Verification

verify_jacobian_at(spec, def, ws, u, t; reference = DenseAD_ForwardDiff(), rtol, atol)
is_approximate(backend)            # trait: true for ApproxSolve
benchmark_jacobian_backends(...)   # timing harness

Or run verification during a solve at chosen times:

sol = solve(problem; verify_jacobian = (rtol = 1e-6, atol = 1e-10,
                                        sample_times = [0.0, 1.0, 10.0]))

It warns (does not abort) when the active backend deviates from the dense AD reference.

16.8.5 Solver extensions (vendor factorizations)

Load the dependency, then pass the factorization to a constructor integrator:

using Pyrolysis, KLU
sol = solve(problem; integrator = KenCarp4(linsolve = Pyrolysis.klu_factorization()))

Pyrolysis.pardiso_factorization() (Intel MKL Pardiso, using Pardiso) and Pyrolysis.amg_preconditioner(A) (using AlgebraicMultigrid) work the same way. Calling a stub without its extension loaded raises an error naming the package to add.

MKL ipiv caveat. When handle_depletion = true resizes the state vector, avoid factorizations with fixed pivot buffers (MKL Pardiso). With depletion the solver already uses a dense, resize-safe LUFactorization() and a finite-difference Jacobian (explicit sparse choices are rebuilt to dense LU with a warning), so this is only a concern if you pass a fully built integrator with a non-default linsolve.

16.9 Parameter sweeps

See Chapter 12. parameter_sweep is top-level exported.

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

Because there is no copy-with-overrides constructor, modify must rebuild the objects it changes from scratch:

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)
    comps = prob.material.components
    new_rxn = Reaction(:decomp, 1 => 2; A = p.A, E = p.E, h = 200e3)
    new_mat = Material(name = prob.material.name,
                       components = comps,
                       reactions = (new_rxn,),
                       mixing = prob.material.mixing)
    return PyrolysisProblem(mesh = prob.mesh, material = new_mat,
                            bc_set = prob.bc_set, tspan = prob.tspan,
                            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 = modify, extract = extract)

16.10 Sensitivity analysis (extension)

See Chapter 13 and Technical Reference §17. forward_sensitivity and adjoint_sensitivity are top-level exported stubs — they gain methods only under using SciMLSensitivity, Zygote (which loads PyrolysisSciMLSensitivityExt).

forward_sensitivity(base_problem, θ::AbstractVector, p_inject::Function;
                    output_fn::Function = identity, solve_kwargs...)
    -> (output, ∂output_∂θ)

adjoint_sensitivity(base_problem, θ::AbstractVector, p_inject::Function, loss_fn::Function;
                    sensealg = InterpolatingAdjoint(autojacvec = ReverseDiffVJP()),
                    solve_kwargs...)
    -> (loss, grad)

Forward mode (ForwardDiff) suits small θ; adjoint mode (Zygote + SciMLSensitivity) suits large θ / scalar loss (inverse problems). p_inject is the same pattern as modify in §16.9 — rebuild the perturbed problem from θ. Sensitivities are reliable thanks to the AD-clean design (tanh gates, midpoint enthalpy, smooth limiters); avoid warm starts and the depletion callback when differentiating (non-smooth).

using Pyrolysis, SciMLSensitivity, Zygote

θ0 = [8.5e12, 188e3]   # A, E
p_inject = (prob, θ) -> begin
    rxn = Reaction(:decomp, 1 => 2; A = θ[1], E = θ[2], h = 870e3)
    mat = Material(name = :PMMA, components = prob.material.components,
                   reactions = (rxn,), mixing = prob.material.mixing)
    PyrolysisProblem(mesh = prob.mesh, material = mat, bc_set = prob.bc_set,
                     tspan = prob.tspan, output_times = prob.output_times)
end

out, J = forward_sensitivity(problem, θ0, p_inject;
                             output_fn = sol -> [maximum(sol.surface_T)])

16.11 Adaptivity and callbacks

See Technical Reference §14. The depletion path is production-ready; h-refinement/coarsening and the AMRController are in Pyrolysis.Experimental.

16.11.1 Surface-depletion via solve

The simplest path is the solve keywords (§16.7): handle_depletion = true (needs use_ale = true), depletion_threshold (default 0.05), min_cells (default 2). Remember this falls back to a finite-difference Jacobian.

16.11.2 Callbacks (not top-level; use Internal/Solver)

using Pyrolysis.Solver: create_depletion_callback, create_thickness_termination_callback,
                        create_step_update_callback, DepletionCallbackState,
                        get_depletion_callback_stats

create_thickness_termination_callback(sv_ref, min_thickness::Float64)
create_depletion_callback(ws, mesh, material, sv_or_state; threshold = 0.05, min_cells = 2)
create_step_update_callback(ws, material)
get_depletion_callback_stats(cb)                  # merge history

You normally do not build these by hand — solve assembles them from its keywords. They are exposed for custom solver loops.

16.11.3 Error indicators, monitors, quality metrics (Adaptivity submodule)

For manual r-refinement / mesh-quality work, import from Pyrolysis.Adaptivity:

using Pyrolysis.Adaptivity: GradientIndicator, CurvatureIndicator,
                            ReactionRateIndicator, InterfaceProximityIndicator,
                            CompositeIndicator, compute_indicator,
                            GradientMonitor, InterfaceMonitor, CompositeMonitor,
                            relocate_nodes!, AspectRatioMetric, SmoothnessMetric,
                            JacobianMetric, MinCellSizeMetric, CompositeQualityMetric,
                            check_quality, enforce_quality!, mesh_quality_summary

relocate_nodes!(mesh; monitor = GradientMonitor(), max_iter = 10, relaxation = 0.5, tol = 1e-3, h_min = 1e-8) moves interior nodes by equidistribution (boundary nodes are pinned). Quality metrics carry defaults (e.g. AspectRatioMetric(2.0), SmoothnessMetric(0.5), MinCellSizeMetric(1e-8)).

16.11.4 Experimental h-AMR

refine_cell!, perform_h_refinement!, coarsen_cells!, AMRController, adapt_mesh!, create_default_controller, etc. are in Pyrolysis.Experimental and are unstable. The active-cell bookkeeping requires iterating the mask, not eachindex; treat this path as research-grade.

16.12 Conservation diagnostics

See Technical Reference §16. Enable snapshots inline with solve(problem; diagnostics = true, diagnostics_stride = N) — rows land in sol.extras.diagnostics_log. The tracker types and check functions are in Pyrolysis.Diagnostics (and in Internal).

using Pyrolysis.Internal   # or: using Pyrolysis.Diagnostics

# trackers
MassConservationTracker, EnergyConservationTracker, ReactionMassFlowTracker
DiagnosticsState, DiagnosticsRow, DiagnosticsCallback

# instantaneous quantities
compute_domain_mass(mesh, material)             # NTuple of per-component mass [kg]
compute_sensible_energy(mesh, material)         # matrix-only E_matrix [J]
compute_total_sensible_energy(mesh, material)   # E_total (matrix + gas) [J]

# checks (return NamedTuples)
check_mass_conservation(tracker, mesh, material; tol = 1e-10)
check_energy_conservation(tracker, mesh, material; tol = 1e-6)
check_reaction_mass_balance(tracker; tol = 1e-10)

# human-readable reports (return String)
mass_conservation_report(tracker, mesh, material)
energy_conservation_report(tracker, mesh, material)
reaction_mass_flow_report(tracker, material)

The energy diagnostics keep two ledgers: a discrete form that closes to integrator tolerance by construction (audits bookkeeping) and a physical form that quantifies the gas-storage gap ΔE_gas = E_total − E_matrix — a runtime measure of the quasi-steady-gas approximation cost (canonically < 1 %). A DiagnosticsRow carries t, E_matrix, E_total, discrete_residual, discrete_relative, physical_residual, physical_relative, gas_storage_gap, mass_max_relative_error.

The low-level residual surface — ProblemDef, Workspace, residual!, to_problem_def, build_workspace, build_ad_workspace, StateLayout — is top-level exported for users building custom solver or sensitivity loops. The canonical call is residual!(du, u, def, t, ws); build (def, ws) with def = to_problem_def(problem; use_ale, radiation_model, energy_form) (plus, for direct Beer–Lambert use, radiation_intensity and radiation_kirchhoff — set radiation_kirchhoff = false if you have already folded an explicit absorptivity into the intensity, else the live surface ε_eff is applied at injection) and ws = build_workspace(problem, def). See Technical Reference §15.

16.13 Default-value reference

A single place to look up the defaults baked into the package. Physics is in the noted TR chapter.

solve keyword defaults

Component / material / reaction defaults

BC defaults

Problem / mesh defaults

Linear-solver / Jacobian defaults

Kernel constants and limits

16.14 Quick import cheat-sheet