Conservation diagnostics

Diagnostics

Conservation diagnostics.

• MassConservationTracker — per-component mass ledger.
• EnergyConservationTracker — sensible-energy ledger of the temperature equation (E = Σ ρcp·T·V over every component), balanced against the terms the discretized RHS assembles.
• DiagnosticsState and record_substep! — advance both ledgers by one step's boundary fluxes and volumetric sources and log a row every stride steps. The test suite drives them step by step to check conservation; solve does not.
• DiagnosticsCallback — DiscreteCallback builder that runs a record!(u, t) closure on every accepted step. solve uses it to log snapshot_diagnostics_rows (the sensible energy at the step) when diagnostics are requested.

DiagnosticsState

Ledger state advanced by record_substep!. Holds:

• the running mass and energy trackers,
• a row log written every stride recorded steps,
• the time of the last recorded step (for the substep dt).

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

Sensible-energy ledger of the temperature equation. The stored quantity is

E(t) = Σ_i ρcp_i · T_i · V_i,    ρcp = Σ_j ξ_j cp_j(T)  (every component)

and the ledger accumulates the terms the discretized RHS actually assembles:

ΔE(t) = ∫ q_cond_boundary · A · dt
      + ∫ Q_rxn dV · dt
      + ∫ S_conv dV · dt

q_cond_boundary is the conducted boundary flux — the energy the boundary cell receives. The surface energy balance has no gas-enthalpy term, so it equals the applied boundary flux. The residual closes to roughly integrator tolerance × |E(0)| for constant composition and volume; with active reactions or a moving mesh it also carries the composition-change term ∫ T ∂(ρcp)/∂t dV dt and the volume term ∫ ρcp T ∂V/∂t dt, which the per-substep accumulation below does not reconstruct. 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.

Fields

• initial_energy::Float64 — Σ ρcp · T · V at t=0.
• cumulative_boundary_conduction::Float64 — ∫ q_cond_boundary · A · dt.
• cumulative_reaction_heat::Float64 — ∫ Q_rxn dV · dt.
• cumulative_advection_sink::Float64 — ∫ ∫ S_conv dV · dt.
• last_time::Float64 — Time of last update [s].

EnergyConservationTracker(mesh::Mesh1D, material)

Initialize tracker with the current sensible energy.

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.

DiagnosticsCallback(record!)

Construct a DiscreteCallback that fires every accepted step and runs record!(u, t) on the integrator's state u::Vector{Float64} and time t::Float64. It only observes the state, so it clears the derivative-discontinuity flag. Call sites supply a record! closure that has access to the live workspace and pushes a row to the log when appropriate; solve pushes snapshot_diagnostics_rows.

record! is held behind a fixed-signature FunctionWrapper, so the callback's type does not depend on the closure's and the integrator it joins keeps a problem-independent type.

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

Compute the energy-ledger residual E(t) − E(0) − Σ contributions.

Returns

NamedTuple with:

• residual::Float64 — E(t) − E(0) − Σ contributions [J]. ≈ 0 to integrator tolerance for constant composition and volume.
• relative::Float64 — relative magnitude vs |E(0)|.
• conserved::Bool — relative < tol.
• current_energy::Float64 — current sensible heat [J].

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 note

Every component rides the mesh (its only mesh-motion term is the dilation −ξθ, which conserves the cell mass ξV), so under ALE the mass of every component changes only through the transport flux J at the boundary faces — exactly what expected_mass integrates. A fully closed (impermeable) column conserves each component to integrator tolerance (V-CPL-03 P4).

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

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

Sensible heat E = Σ_i ρcp_i · T_i · V_i with the mixture ρcp (every component, gas included — Properties.effective_heat_capacity). Only active cells (AMR-aware).

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.

record_substep!(state::DiagnosticsState, mesh, material;
                q_cond_boundary, Q_rxn_total, S_conv_total,
                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

Snapshot at a single point: the sensible energy. The ledger fields (cumulative residual, mass error) are NaN; only record_substep! fills them. This is the row solve logs.

update_energy_tracker!(tracker; q_cond_boundary, Q_rxn_total,
                       S_conv_total, 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_cond_boundary is a face-flux value [W/m²] (multiplied by A and dt internally).

Keyword arguments

• q_cond_boundary::Float64 — net conducted boundary flux [W/m²], positive into slab.
• Q_rxn_total::Float64 — ∫ Q_rxn dV [W].
• S_conv_total::Float64 — ∫ S_conv dV [W] (typically negative).
• 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)