Conservation diagnostics

Diagnostics

First-class runtime conservation diagnostics for the Pyrolysis.jl energy formulation.

Two trackers + a callback:

• MassConservationTracker — per-component mass ledger.
• EnergyConservationTracker — both the discrete-form residual (closes to integrator tolerance by construction) and the physical-form residual (quantifies the cost of the quasi-steady-gas approximation that Formulation A makes by excluding gases from ρcp_eff).
• ReactionMassFlowTracker — per-reaction mass balance.
• DiagnosticsCallback — DiscreteCallback builder that fires every stride steps and at the final time, accumulating the substep contributions into the trackers and emitting a row to a Vector{NamedTuple} log.

The callback is the visibility mechanism that makes the FDS-convention approximations (gas-excluded ρcp_eff, dropped T·cp·∇·J term) auditable on every solve. The log is connected to PyrolysisSolution.extras.

DiagnosticsState

Mutable container threaded through the callback. Holds:

• the running trackers,
• a row log written every stride accepted steps,
• the last time at which the callback ran (for substep dt computation).

Built once at the start of solve() and lives on the Workspace (or as a free-standing object passed via cb_kwargs).

DiagnosticsState(mesh, material; stride=10) -> DiagnosticsState

Initialize the diagnostics state with fresh trackers seeded from the current mesh state (this becomes t=0 baseline).

EnergyConservationTracker

Tracks two complementary energy ledgers for the Formulation A pyrolysis residual:

1. Discrete-form ledger (closes to integrator tolerance)

Built from the terms the discretized RHS actually assembles:

ΔE_matrix(t) = ∫ q_cond_boundary · A · dt
             + ∫ Q_rxn dV · dt
             + ∫ S_conv dV · dt
             + ∫ S_gen dV · dt          (zero unless `:with_generation_sink`)

with E_matrix = Σ_i ρcp_eff,i · T_i · V_i using the gas-excluded ρcp_eff (quasi-steady-gas approximation; NOT the ThermaKin convention, which keeps gas storage — see Properties.effective_heat_capacity). This residual closes to roughly integrator tolerance × |E(0)| by construction; failures here indicate bugs in flux/source assembly or callbacks. NOTE: because it is built from the terms the RHS actually assembles, this ledger closes even if the RHS itself double-counts a mechanism — it audits discrete bookkeeping, not the physical model. cumulative_transpiration is zero under the default (transpiration-off) convention.

2. Physical-form ledger (auditable cost of the FDS-convention approximation)

Built from total sensible energy including gas storage. The physical-form residual quantifies the cost of the quasi-steady-gas approximation:

E_total = Σ_i ρcp_total,i · T_i · V_i     (matrix + gas)
physical_residual(t) = E_total(t) − E_total(0)
                      − cumulative_boundary_applied
                      − cumulative_reaction_heat
                      + cumulative_transpiration

For canonical pyrolysis cases (gas residence time ≪ thermal time scale) this should be < 1% × |E(0)|; if it is larger, the regime stresses the quasi-steady-gas approximation and the user should reconsider.

Fields

• initial_energy_matrix::Float64 — Σ ρcp_eff · T · V at t=0 (matrix only).
• initial_energy_total::Float64 — Σ ρcp_total · T · V at t=0 (matrix + gas).
• cumulative_boundary_conduction::Float64 — ∫ q_cond_boundary · A · dt.
• cumulative_boundary_applied::Float64 — ∫ q_BC_applied · A · dt.
• cumulative_reaction_heat::Float64 — ∫ Q_rxn dV · dt.
• cumulative_advection_sink::Float64 — ∫ ∫ S_conv dV · dt.
• cumulative_generation_sink::Float64 — ∫ ∫ S_gen dV · dt.
• cumulative_transpiration::Float64 — ∫ Σ_g ṁ_g · Δh_g(T_s, T_∞) · A · dt.
• last_time::Float64 — Time of last update [s].

EnergyConservationTracker(mesh::Mesh1D, material)

Initialize tracker with current sensible energies (matrix-only and total-including-gas).

MassConservationTracker{NC}

Tracks cumulative mass for each component throughout simulation.

Fields

• initial_mass::NTuple{NC, Float64}: Initial mass per component [kg]
• cumulative_efflux::NTuple{NC, Float64}: Cumulative mass leaving domain [kg]
• last_time::Float64: Time of last update [s]

MassConservationTracker(mesh::Mesh1D{NC}, material) where NC

Initialize tracker with current domain mass.

ReactionMassFlowTracker{NC,NR}

Tracks mass flow through each reaction over time.

For each reaction i, tracks:

• Cumulative reactant consumed: Σ∫ rᵢ × νᵢⱼ dt for reactants
• Cumulative product generated: Σ∫ rᵢ × νᵢⱼ dt for products

This allows validation that each reaction conserves mass independently.

Fields

• reactant_consumed::Matrix{Float64}: [NC × NR] cumulative consumption [kg]
• product_generated::Matrix{Float64}: [NC × NR] cumulative generation [kg]
• reaction_counts::Vector{Float64}: [NR] cumulative ∫ rᵢ dt (extent of reaction) [kg]
• last_time::Float64: Time of last update [s]

ReactionMassFlowTracker(material::Material{NC,NR}) where {NC,NR}

Initialize reaction mass flow tracker.

DiagnosticsCallback(record!::Function)

Construct a DiscreteCallback (DiffEqCallbacks-compatible) that fires every accepted step and runs record!(integrator). Call sites supply a record! closure that has access to the live workspace and pushes a row to the log when appropriate. This module ships snapshot_diagnostics_row as the simplest builder; a full ledger callback that calls record_substep! requires per-substep accumulators that the unified residual layer threads through.

check_energy_conservation(tracker, mesh, material; tol=1e-6) -> NamedTuple

Compute both the discrete-form residual (closes to integrator tolerance by construction) and the physical-form residual (quantifies the quasi-steady-gas approximation cost).

Returns

NamedTuple with:

• discrete_residual::Float64 — E_matrix(t) − E_matrix(0) − Σ contributions. Should be ≈ 0 to integrator tolerance.
• discrete_relative::Float64 — relative magnitude vs |E_matrix(0)|.
• physical_residual::Float64 — E_total(t) − E_total(0) − applied + transpiration − reaction.
• physical_relative::Float64 — relative magnitude vs |E_total(0)|.
• current_energy_matrix::Float64 — current matrix sensible heat.
• current_energy_total::Float64 — current total sensible heat (matrix + gas).
• gas_storage_gap::Float64 — (E_total − E_matrix), the energy held in gas storage that the discrete equation does not evolve.
• discrete_conserved::Bool — discrete_relative < tol.
• physical_conserved::Bool — physical_relative < tol.

check_mass_conservation(tracker::MassConservationTracker{NC},
                        mesh::Mesh1D{NC}, material;
                        tol::Float64=1e-10) -> NamedTuple

Check mass conservation and return detailed diagnostics.

Conservation equation: mcurrent(t) + mescaped(t) = m_initial

Or equivalently: mcurrent(t) = minitial - m_escaped(t)

Arguments

• tracker: Mass conservation tracker
• mesh: Current mesh state
• material: Material definition
• tol: Relative error tolerance (default: 1e-10)

ALE gas-species caveat

The default tol = 1e-10 is a bookkeeping-level tolerance and is attainable only where the discrete ledger telescopes exactly: condensed components (mesh-bound, flux-form), and all components on a fixed (Eulerian) mesh. Under ALE the gas-species mesh-motion correction is the pointwise chain-rule term +w·∇ξ_g (non-conservative form), which does not telescope cell-by-cell — gas-species closure is then truncation-order in (Δt, Δz), not roundoff. Loosen tol accordingly for gas components in ALE runs (guide §16.3.3/16.3.4).

Returns

NamedTuple with:

• conserved::Bool: True if mass is conserved within tolerance
• current_mass::NTuple{NC, Float64}: Current mass in domain [kg]
• expected_mass::NTuple{NC, Float64}: Expected mass (initial - efflux) [kg]
• relative_error::NTuple{NC, Float64}: Relative error per component
• max_error::Float64: Maximum relative error across components

check_reaction_mass_balance(tracker::ReactionMassFlowTracker{NC,NR};
                             tol::Float64=1e-10) where {NC,NR} -> NamedTuple

Check mass balance for each reaction.

Arguments

• tracker: ReactionMassFlowTracker
• tol: Relative tolerance for mass balance error

Returns

NamedTuple with:

• balanced::Vector{Bool}: Whether each reaction is balanced
• errors::Vector{Float64}: Relative mass error per reaction
• all_balanced::Bool: True if all reactions are balanced

compute_domain_mass(mesh::Mesh1D{NC}, material) -> NTuple{NC, Float64}

Compute total mass of each component in the domain.

Mass for component j: mⱼ = Σᵢ ξⱼ(i) × V(i)

Only includes active cells (excludes refined-away cells in AMR).

compute_sensible_energy(mesh::Mesh1D{NC}, material) -> Float64

Matrix-only sensible heat: E_matrix = Σ_i ρcp_eff,i · T_i · V_i. Gas components are excluded from ρcp_eff per Formulation A. Only active cells (AMR-aware).

compute_total_sensible_energy(mesh::Mesh1D{NC}, material) -> Float64

Total sensible heat (matrix + gas storage): E_total = Σ_i ρcp_total,i · T_i · V_i, where ρcp_total includes all components per Properties.total_heat_capacity. Used by the physical-form energy residual to quantify the cost of the quasi-steady-gas approximation that the matrix-only ρcp_eff makes.

energy_conservation_report(tracker::EnergyConservationTracker,
                           mesh::Mesh1D, material) -> String

Generate a human-readable energy conservation report.

Useful for post-processing and debugging.

mass_conservation_report(tracker::MassConservationTracker{NC},
                         mesh::Mesh1D{NC}, material) -> String

Generate a human-readable mass conservation report.

Useful for post-processing and debugging.

reaction_mass_flow_report(tracker::ReactionMassFlowTracker{NC,NR},
                          material::Material{NC,NR}) where {NC,NR} -> String

Generate detailed report of mass flow through reactions.

Example

tracker = ReactionMassFlowTracker(material)
# ... run simulation with update_reaction_tracker! calls ...
println(reaction_mass_flow_report(tracker, material))

record_substep!(state::DiagnosticsState, mesh, material;
                q_cond_boundary, q_bc_applied, Q_rxn_total,
                S_conv_total, S_gen_total, q_transpiration,
                J_boundary, A, t)

Update the trackers with one substep's contributions and append a typed row to the log if step_counter % stride == 0.

Substep size dt = t − state.last_t is computed automatically.

snapshot_diagnostics_row(t, mesh, material) -> DiagnosticsRow

Cumulative-free snapshot at a single point: matrix-only energy, total energy (includes gas), and the gas-storage gap. Other ledger fields (cumulative residuals, mass error) are filled with NaN to mark them as not-yet-computed. This is what the snapshot-mode callback writes; the cumulative form requires per-substep accumulation that the unified residual layer wires through record_substep!.

update_energy_tracker!(tracker; q_cond_boundary, q_bc_applied,
                       Q_rxn_total, S_conv_total, S_gen_total,
                       q_transpiration, A, dt)

Update tracker with one substep's contributions.

All *_total arguments are volume-integrals over the cell-source field (Σ_i field_i · V_i) computed from the active RHS state. q_*_boundary are face-flux values [W/m²] (multiplied by A and dt internally).

Keyword arguments

• q_cond_boundary::Float64 — net conducted boundary flux [W/m²], positive into slab.
• q_bc_applied::Float64 — applied external BC flux [W/m²], positive into slab.
• Q_rxn_total::Float64 — ∫ Q_rxn dV [W].
• S_conv_total::Float64 — ∫ S_conv dV [W] (typically negative).
• S_gen_total::Float64 — ∫ S_gen dV [W] (zero unless :with_generation_sink).
• q_transpiration::Float64 — Σ_g ṁ_g · Δh_g(T_s, T_∞) summed over boundary faces [W/m²].
• A::Float64 — Cross-sectional area at this substep [m²].
• dt::Float64 — Substep size [s].

update_mass_tracker!(tracker::MassConservationTracker{NC},
                     J_boundary::Matrix{Float64},
                     A::Float64, dt::Float64) where NC

Update tracker with boundary mass fluxes.

Arguments

• tracker: The mass conservation tracker
• J_boundary: Boundary mass fluxes [NC × 2] for top/bottom boundaries [kg/(m²·s)]
• A: Cross-sectional area [m²]. Under variable cross-section, pass the trapezoidal-in-A average over the substep, A = 0.5 * (A_t_pre + A_t_post). With identity law / non-ALE, A is constant and equals mesh.cross_section_area_initial.
• dt: Time step [s]

Notes

• Convention: Positive J = flux into domain (inflow)
• Mass leaving = -J × A × dt (negative J = outflow)
• Typical case: Jtop < 0 (gas escaping), Jbottom = 0 (impermeable)

update_reaction_tracker!(tracker::ReactionMassFlowTracker{NC,NR},
                         rates::Vector{Float64},
                         mesh::Mesh1D{NC},
                         material::Material{NC,NR},
                         dt::Float64) where {NC,NR}

Update reaction mass flow tracker with current rates.

Arguments

• tracker: ReactionMassFlowTracker
• rates: Current reaction rates [NR] in kg/(m³·s) (domain-averaged or from a specific cell)
• mesh: Mesh (to compute total domain volume)
• material: Material definition
• dt: Time step [s]

Notes

This function should be called once per time step with domain-integrated rates.