Residual layer

SimulationMode{Mesh, Rad, Trans, Chi}

Bundled trait carrying all four orthogonal axes as type parameters. The residual takes SimulationMode as a single parameter so adding new axes does not multiply method-table size.

StateLayout — unified state-vector descriptor.

Describes the packing of the state vector u in terms of the bundled SimulationMode. Block layout (T, then ξ, then optional z, then optional χ) is the single source of truth for every pack / unpack / view operation in the package.

Block presence is strictly determined by SimulationMode:

• T and ξ blocks are always present.
• z block is present iff mesh_mode(mode) isa ALE.
• χ block is present iff chi_mode(mode) isa WithChi.

Type Parameters

• NC: number of components
• Mode <: SimulationMode: bundled trait

Block-major offsets (1-based, with n = n_cells)

• T: 1 .. n
• ξ: n + (j-1)*n + (1..n) for j ∈ 1:NC
• z: n*(1+NC) + (1..n_nodes) (ALE only)
• χ: trailing block of length n (only if WithChi)

ProblemDef — immutable problem definition.

A ProblemDef carries the unchanging contract of a pyrolysis problem: material, topology, boundary conditions, initial conditions, time span, options. Every mutable cache lives in Workspace instead.

This split is what makes the residual a pure function of (du, u, def, t, ws) and unlocks parameter sweeps, sensitivity analysis, and pluggable Jacobian backends in later phases.

Field access is direct: read def.material, def.topology, def.bcs, def.tspan, def.mode, def.options. The single helper is get_option for fallback semantics on the runtime-options NamedTuple.

Type Parameters

• M: Material type
• Top: MeshTopology type
• B: BoundaryConditionSet type
• IC: Initial-conditions container type
• Mode: bundled SimulationMode type
• Opts: NamedTuple of runtime options

Workspace — mutable container that owns every per-step cache and work array used by the unified residual.

The Workspace is the only place where mutable per-step state lives. Material, boundary conditions, and the radiation model are read from def::ProblemDef (the immutable side); the workspace does not duplicate them.

Every cache is parameterized so that Nothing can appear as a concrete type parameter rather than as a Union{Nothing, T}. Trait dispatch on the type — not runtime isnothing checks — decides whether to touch a cache.

Type Parameters

• NC: number of components
• Mode: bundled SimulationMode (matches the parent ProblemDef)
• Tu: scalar eltype carried by the live state (Float64 for the canonical workspace; ForwardDiff.Dual{Tag, Float64, N} for the AD workspace allocated by an AD Jacobian backend)
• MeshT: Reference to the mutable Mesh1D carrying connectivity and per-step geometry. The mesh is read-only during residual / Jacobian assembly; writes happen only between integrator-accepted steps via StepUpdateCallback. Kernels only read connectivity (cell IDs, face IDs, neighbor lookups, boundary tags) off the mesh; live state and state-derived geometry live in state_cache.
• GeomT: MeshGeometry view
• LayoutT: StateLayout descriptor
• GeomCacheT: Geometry.GeometryCache for ALE; Nothing for Eulerian
• RadCacheT: Properties.RadiationCache for P1; Nothing otherwise
• BCStateT: per-boundary-tag state container (NamedTuple)
• StateCacheT: StateCache{NC, Tu, Mode} — typed mirror of the live state populated at every residual entry

to_problem_def(problem::PyrolysisProblem; use_ale=false,
                radiation_model=NO_RADIATION,
                radiation_intensity=0.0,
                radiation_kirchhoff=true,
                energy_form=:standard,
                use_transpiration_bc=false,
                use_blowing_correction=false) -> ProblemDef

Lift a user-facing PyrolysisProblem into the immutable ProblemDef. The returned object shares the underlying material/BC/initial-condition data with the input but no longer carries any mutable mesh state.

radiation_kirchhoff records the surface-transmissivity convention for the Beer–Lambert injection: true (the default) means radiation_intensity is the bare incident flux and the live surface-cell ε_eff is applied at injection (Kirchhoff α = ε); false means the RadiativeFluxBC's explicit absorptivity is already folded into the intensity. solve() derives it from the top thermal BC.

use_blowing_correction (off by default) rescales every ConvectiveBC in the surface energy balance by the Couette-film blowing factor h → h·B/(e^B − 1), B = ṁ″·c_p,g/h, built from the outward gas mass flux at each boundary face. This is the momentum-side companion of the transpiration enthalpy term and is fully compatible with the default Sconv-only convention (it modifies the convective coefficient, not the gas enthalpy accounting). At peak burning (B ≳ 1) it suppresses convective cooling substantially — e.g. B = 4 gives heff/h = 0.075. See Discretization.blowing_factor.

build_workspace(problem::PyrolysisProblem, def::ProblemDef;
                  min_thickness=1e-6,
                  ale_limiter=VANLEER, scalar_eltype::Type=Float64) -> Workspace

Build a Workspace that owns every per-step cache and work array used by the unified residual. Each cache is sized to the current mesh and component count.

scalar_eltype is the live state's eltype (Tu parameter) — Float64 for the canonical workspace solve() allocates; ForwardDiff.Dual{...} for the AD workspace allocated lazily by an AD Jacobian backend.

The mesh carried by the workspace is the same mutable Mesh1D referenced by problem. The mesh is read-only during residual / Jacobian assembly; live state lives in ws.state_cache and is written back to the mesh only between accepted timesteps via StepUpdateCallback.

build_ad_workspace(def::ProblemDef, ws::Workspace, ::Type{Tu}) -> Workspace

Build a parallel Workspace{NC, Mode, Tu} whose state cache, property cache, RHS cache, and physics work arrays carry scalar eltype Tu (typically ForwardDiff.Dual{Tag, Float64, N}), sharing the canonical workspace's mesh, layout, boundary IDs, and runtime settings.

Used by AD Jacobian backends to keep a Dual-typed workspace alongside the Float64 one — the residual then dispatches on eltype(ws.dT_work) == Tu end-to-end without per-call reallocation.

The mesh is shared (same identity), reflecting the read-only mesh contract during residual / Jacobian assembly.

Unified ODE right-hand side for the pyrolysis package.

function residual!(du, u, def::ProblemDef, t, ws::Workspace) -> nothing

This is the single RHS for all non-P1 simulations. The trait dispatch chain is:

residual!(du, u, def, t, ws)
	→ _residual!(du, u, def, t, ws, mesh_mode(def.mode))
		::Eulerian → Eulerian path
		::ALE      → ALE path

Inside each path, transport is dispatched via transport_mode(def.mode) (Fickian vs DarcyFick) and the χ block is included iff chi_mode(def.mode) isa WithChi.

P1 radiation has its own residual in Pyrolysis.Experimental.RadiationP1 and is not routed through this entry point.

Contract

• def is immutable — safe to share across threads / sweep workers.
• ws is mutable — one per worker.
• t is a plain Real; the kernel must accept whatever scalar type du/u carry (Float64, ForwardDiff.Dual, …).

state_length(l::StateLayout) -> Int

Return the total length of the state vector encoded by this layout.

create_state_vector(layout::StateLayout, mesh) -> Vector{Float64}

Allocate a fresh state vector sized to layout and pack the current mesh state into it. Used by solve() to build u0.

pack_state!(u, mesh, layout::StateLayout) -> u

Pack the mesh's cell_states (and node positions when ALE; χ when WithChi) into the flat ODE state vector u. Mode is read from the layout's type parameter.

update_mesh_from_state!(mesh, u::AbstractVector{Float64}, layout::StateLayout) -> mesh

Mirror an integrator-accepted Float64 state vector into the mutable mesh. Called by StepUpdateCallback between timesteps so AMR, depletion checks, and post-processing read a consistent mesh.

This function is not invoked during residual / Jacobian assembly — those are pure functions of u and the typed state_cache. The mesh is read-only during assembly.

block_view(u, l::StateLayout, ::Val{:T}) -> view of T block
block_view(u, l::StateLayout, ::Val{:ξ}, j::Int) -> view of component j
block_view(u, l::StateLayout, ::Val{:z}) -> view of node-position block (ALE)
block_view(u, l::StateLayout, ::Val{:χ}) -> view of pyrolysis-progress block (WithChi)

Compile-time-dispatched view into the named block of u. Returns an empty view when the block does not exist for the layout's mode.

get_temperature(u, layout, i) -> Real

get_concentration(u, layout, i, j) -> Real

get_node_position(u, layout, i) -> Real

Only valid for ALE layouts.

get_node_positions(u, layout) -> view

View of node positions for ALE layouts.

get_thickness_from_state(u, layout) -> Real

Compute current domain thickness z_top - z_bottom from the ALE state. Returns 0 for Eulerian layouts.

column_chi_bar_from_state(u, layout, mesh) -> Real

Mass-weighted column-averaged χ̄. AD-Dual safe (numerator carries duals, denominator from the immutable m_dry_solid_initial). Thin wrapper over the single implementation _chi_bar_from_state, which also returns the column dry-solid mass used by update_cross_section_area!.

cross_section_area_from_state(mesh, material, u, layout) -> Real

Current column cross-sectional area implied by the χ block. Returns mesh.cross_section_area_initial for identity-law materials.