LagrangianSolver.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: Apache-2.0
//
// Copyright 2026 Caleb Buahin
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

 
#ifndef OPENSWMM_QUALITY_LARD_LAGRANGIAN_SOLVER_HPP
#define OPENSWMM_QUALITY_LARD_LAGRANGIAN_SOLVER_HPP
 
#include <algorithm>
#include <cmath>
#include <limits>
#include <cstdint>
#include <vector>
 
#include "../../core/SimulationContext.hpp"
#include "../../transport/InitialQualitySeeds.hpp"
#include "../NegativeSources.hpp"
#include "../QualityRouting.hpp"
#include "../../transport/components/ReactionModule/ReactionLegacyBinding.hpp"
#include "../../transport/TransportPolicy.hpp"   // E2
#include "RwptDispersion.hpp"
#include "SegmentStore.hpp"
#include "../../core/UnitConversion.hpp"
#include "../../hydraulics/Node.hpp"
#include "../../transport/components/HeatFluxModules/BedExchange.hpp"
#include "../../transport/components/HeatFluxModules/HeatFluxes.hpp"
#include "../../transport/components/HeatFluxModules/HeatOverrides.hpp"
#include "../../transport/components/HeatFluxModules/SolarRadiation.hpp"
#include "../../transport/components/HeatFluxModules/SurfaceExchange.hpp"
#include "../QualityRouting.hpp"
 
namespace openswmm {
namespace lard {
 
constexpr double kTinyFlow = 1.0e-8;  

struct SpeciesRowLayout {
    int np        = 0;   
    int age_row   = -1;  
    int temp_row  = -1;  
    int msx_first = -1;  
    int ns        = 0;   
};

inline SpeciesRowLayout rowLayout(const SimulationContext& ctx) {
    // E2: the class ENABLES come from the one policy
    // (transport::network1DEnables — allocation-free, this runs per step).
    // Identical to the pre-E2 reads: LARD only runs with IGNORE_QUALITY NO,
    // where the policy's gate is a no-op. The row ORDER below is deliberately
    // LARD's own (pollutants, age, temperature, MSX) — see the struct note;
    // E2 unifies the decision, not the index arithmetic.
    const transport::ClassEnables en = transport::network1DEnables(ctx);
    SpeciesRowLayout L;
    L.np = en.n_pollut;
    L.ns = L.np;
    if (en.age)         L.age_row  = L.ns++;
    if (en.temperature) L.temp_row = L.ns++;   // H7b
    // L3: MSX species ride the segments after the reserved rows —
    // stable across the run (the compiled species table is fixed at
    // open). n_msx == 0 leaves every index above IDENTICAL to H7b's
    // layout, which is the bit-inertness claim the corpus checks.
    if (transport::legacyReactionsActive(ctx)) {
        L.msx_first = L.ns;
        L.ns += en.n_msx;
    }
    return L;
}
 
class LagrangianSolver {
public:
    void step(SimulationContext& ctx, double dt_routing) {
        // X4/H7b: the age and temperature rows ride the segments after
        // the pollutants — age published to water_age_state (seconds),
        // temperature to heat_state (degC) — rather than the np-strided
        // conc arrays.
        const SpeciesRowLayout L = rowLayout(ctx);
        const int np = L.np;
        const int ns = L.ns;
        if (ns <= 0) return;
        if (!initialized_) init(ctx);
 
        const int nl = ctx.n_links();
        auto& links = ctx.links;
 
        // ---- 0. Flow reversal (§4.4) — once per routing step: the flow
        //      solution is constant within it. -----------------------------
        bool topo_dirty = false;
        for (int l = 0; l < nl; ++l) {
            const auto ul = static_cast<std::size_t>(l);
            const double q = links.flow[ul];
            if (std::abs(q) <= kTinyFlow) continue;
            const std::int8_t sign = (q >= 0.0) ? 1 : -1;
            if (sign != flow_sign_[ul]) {
                if (links.type[ul] == LinkType::CONDUIT) store_.reverse(l);
                flow_sign_[ul] = sign;
                topo_dirty = true;
            }
        }
        if (topo_dirty || topo_.empty()) computeTopoOrder(ctx);
 
        // ---- X3a: QUALITY_STEP substepping (strategy §4.2). Equal
        //      substeps; mass/age external loads are RATES and scale
        //      through dt, the per-routing-step external VOLUME
        //      (qual_vol_in) scales through frac. dtq absent or >= the
        //      routing step degenerates to one substep — bit-identical to
        //      the pre-X3a engine by construction.
        const double dtq = ctx.options.quality_step;
        const int nsub = (dtq > 0.0 && dtq < dt_routing)
                             ? static_cast<int>(std::ceil(dt_routing / dtq))
                             : 1;
        const double dt = dt_routing / static_cast<double>(nsub);
        const double frac = 1.0 / static_cast<double>(nsub);
        for (int sub = 0; sub < nsub; ++sub) substep(ctx, dt, frac);
 
        // H2/H3/H6b under LARD: surface + radiative fluxes and the bed
        // pair, once per routing step (the ARD slot), applied to the
        // segment field before publication.
        applyFluxesAndBed(ctx, dt_routing);
 
        publish(ctx, dt_routing);
    }

    void substep(SimulationContext& ctx, double dt, double frac) {
        const SpeciesRowLayout L = rowLayout(ctx);
        const int np = L.np;
        const int ns = L.ns;
        // Derived from the layout rather than re-read from options, so the
        // "is there a reserved row" question has exactly one answer per call.
        const bool age  = (L.age_row >= 0);
        const bool heat = (L.temp_row >= 0);
        const int nn = ctx.n_nodes();
        const int nl = ctx.n_links();
        auto& nodes = ctx.nodes;
        auto& links = ctx.links;
        auto& ws = ctx.water_age_state;
        auto& hs = ctx.heat_state;
 
        // ---- AGE (before transport): every parcel ages by exactly dt, and
        //      the aged value is this substep's "old" state — the plan §1
        //      convention routeLegacyAge follows (age then mix). ------------
        if (age) {
            store_.add_species(L.age_row, dt);
            for (int n = 0; n < nn; ++n)
                ws.node_age[static_cast<std::size_t>(n)] += dt;
        }
 
        // ---- 1. DRAIN (all conduits, before any node mixes) ---------------
        scratch_.assign(static_cast<std::size_t>(ns), 0.0);
        treat_cin_.assign(static_cast<std::size_t>(np > 0 ? np : 0), 0.0);
        for (int l = 0; l < nl; ++l) {
            const auto ul = static_cast<std::size_t>(l);
            if (links.type[ul] != LinkType::CONDUIT) continue;
            const double q = std::abs(links.flow[ul]);
            const int dn = downstreamNode(ctx, l);
            const int up = upstreamNode(ctx, l);
 
            if (q > kTinyFlow && dn >= 0) {
                std::fill(scratch_.begin(), scratch_.end(), 0.0);
                const double drained = store_.drain_back(l, q * dt,
                                                         scratch_.data());
                addToLedger(dn, drained, scratch_.data(), ns);
            }
            // Volume reconciliation: the slab must sum to links.volume.
            // A shortfall is filled at RELEASE with upstream water; an
            // excess left after the outflow drain leaves through the FRONT
            // to the upstream ledger — booked, not rescaled (see header).
            const double v_new = links.volume[ul];
            const double v_rem = store_.total_volume(l);
            if (v_rem > v_new && up >= 0) {
                std::fill(scratch_.begin(), scratch_.end(), 0.0);
                const double shed = store_.drain_front(l, v_rem - v_new,
                                                       scratch_.data());
                addToLedger(up, shed, scratch_.data(), ns);
            }
        }
 
        // ---- 2. MIX in topo order, passthrough zero-volume links ----------
        for (const int n : topo_) {
            const auto un = static_cast<std::size_t>(n);
            const double v_old = nodes.old_volume[un];
            const double v_in = node_vol_in_[un];
            // OUTFALL_BACKFLOW_QUALITY ZERO: an outfall taking no volume
            // inflow this substep is a fresh boundary — its held state
            // (pollutant rows AND the age row) reads zero, so the RELEASE
            // seeding and zero-volume passthrough below draw clean water.
            const bool zero_bf =
                ctx.options.outfall_backflow_zero &&
                nodes.type[un] == NodeType::OUTFALL &&
                v_in + nodes.qual_vol_in[un] * frac <= 0.0;
 
            for (int s = 0; s < ns; ++s) {
                // H7a: identity by INDEX, not by threshold. `s >= np` was
                // correct only while age was the sole reserved row; it would
                // silently capture the temperature row too (H7b).
                const bool is_age  = (s == L.age_row);
                const bool is_temp = (s == L.temp_row);
                // L3: species rows — state lives in msx_node_conc
                // ([node*nsp + sp]); U2 gave them the [INFLOWS] pathway
                // (msx_ext_mass_in, same index).
                const bool is_msx =
                    (L.msx_first >= 0 && s >= L.msx_first);
                const auto xi =
                    is_msx ? un * static_cast<std::size_t>(
                                      ns - L.msx_first) +
                                 static_cast<std::size_t>(s - L.msx_first)
                           : 0;
                const auto li = un * static_cast<std::size_t>(ns) +
                                static_cast<std::size_t>(s);  // ledger index
                // State and external load per row. Pollutants: nodes.conc +
                // qual_mass_in (rate × dt — the mixAtNodes convention). Age:
                // water_age_state.node_age (already aged +dt this step) +
                // node_age_vol_in (age·ft³/s rate, the D-UT10 parallel
                // accumulator filled by all seven loader pathways).
                const auto ci = un * static_cast<std::size_t>(np) +
                                static_cast<std::size_t>(s);
                // Temperature mirrors age one row over: state in
                // heat_state.node_temp (degC), external load from
                // node_temp_vol_in (degC.ft3/s, the D-UT10 twin filled by
                // the same seven loaders) -- H1's convention, consumed here
                // instead of by routeLegacyHeat.
                const double st_old =
                    is_age  ? ws.node_age[un]
                  : is_temp ? hs.node_temp[un]
                  : is_msx  ? ctx.reactions.msx_node_conc[xi]
                            : nodes.conc[ci];
                // U2: species rows take the [INFLOWS] species loads
                // (msx_ext_mass_in, a rate; empty when no such row).
                const double m_ext =
                    is_age  ? ws.node_age_vol_in[un] * dt
                  : is_temp ? hs.node_temp_vol_in[un] * dt
                  : is_msx  ? ((xi < ctx.reactions.msx_ext_mass_in.size())
                                   ? ctx.reactions.msx_ext_mass_in[xi] * dt : 0.0)
                            : nodes.qual_mass_in[ci] * dt;
                double m = st_old * v_old + node_mass_in_[li] + m_ext;
                // P2.3: stash the ARRIVING mass for the treatment cin below
                // (pollutant rows only; s == p for rows 0..np-1).
                if (!is_age && !is_temp && !is_msx &&
                    s < static_cast<int>(treat_cin_.size()))
                    treat_cin_[static_cast<std::size_t>(s)] =
                        node_mass_in_[li] + m_ext;
                // D-NS1 (X6, now observable): extraction beyond the
                // store's mass clamps to available — counted, warned
                // once, and (pollutant rows) un-booked so the ledger
                // carries what actually left.
                // The non-negativity clamp does NOT apply to temperature:
                // degC water below zero is an ordinary state (the freezing
                // gate in heat_watershed exists to keep it so), where a
                // negative mass or age is a defect. Same reasoning as the
                // report boundary's deliberate no-mask on temperature.
                if (m < 0.0 && !is_temp) {
                    if (is_age)
                        quality::bookNegativeAgeClamp(ctx, n);
                    else if (!is_msx)  // L3: species have no ledger row
                        quality::bookNegativeSourceClamp(ctx, n, s, -m);
                    m = 0.0;
                }
                const double denom =
                    v_old + v_in + nodes.qual_vol_in[un] * frac;
                // ALWAYS divide by the full denominator. m/denom is a convex
                // combination of st_old and the arriving values, so it can
                // never exceed its inputs; the fallback this replaces
                // divided a mass that included st_old*v_old by a divisor
                // that EXCLUDED v_old, and at a nearly-dry junction that
                // quotient amplified step over step -- measured on a
                // receding-flow deck (inflow stops at 1 h): node
                // concentrations reached 2.7e30, 2.0e281, then inf, and the
                // final-storage row went NaN. Below 1e-12 ft^3 there is no
                // meaningful water and the store keeps its value.
                const double st_new =
                    zero_bf ? 0.0
                            : ((denom > 1.0e-12) ? m / denom : st_old);
                if (is_age)       ws.node_age[un] = st_new;
                else if (is_temp) hs.node_temp[un] = st_new;
                else if (is_msx)  ctx.reactions.msx_node_conc[xi] = st_new;
                else              nodes.conc[ci] = st_new;
                node_mass_in_[li] = 0.0;  // consumed; cycle residue carries
            }
            node_vol_in_[un] = 0.0;
 
            // P2.3: [TREATMENT] applies HERE — to the freshly mixed node
            // state, BEFORE the passthrough and RELEASE below draw it: the
            // LEGACY ordering (mix → treat → links), transplanted. The
            // handoff's end-of-step application outside the solver was a
            // no-op twice over: the evaluator's cin read the LEGACY
            // accumulators (external loads only — 0 at an interior node,
            // and an R-typed expression keeps c_node when cin reads 0),
            // and a junction's ~zero stored volume gave the late write no
            // weight in the next mix. Fed with THIS substep's inflow
            // figures instead; nodes with no treatment pay one flag test.
            if (np > 0 && dt > 0.0 &&
                un < ctx.treatment.has_treatment.size() &&
                ctx.treatment.has_treatment[un]) {
                const double v_in_total =
                    v_in + nodes.qual_vol_in[un] * frac;
                for (int p = 0; p < np; ++p)
                    treat_cin_[static_cast<std::size_t>(p)] =
                        (v_in_total > 1.0e-12)
                            ? std::max(0.0,
                                       treat_cin_[
                                           static_cast<std::size_t>(p)]) /
                                  v_in_total
                            : 0.0;
                quality::applyNodeTreatment(ctx, n, dt, v_in_total / dt,
                                            treat_cin_.data());
            }
 
            // Zero-volume passthrough (§2.4): outgoing pump/orifice/weir/
            // outlet links deliver the node's NEW concentration downstream
            // within this step — the property the topo order exists for.
            for (const int l : node_out_links_[un]) {
                const auto ul = static_cast<std::size_t>(l);
                if (links.type[ul] == LinkType::CONDUIT) continue;
                const double q = std::abs(links.flow[ul]);
                if (q <= kTinyFlow) continue;
                const int dn = downstreamNode(ctx, l);
                if (dn < 0) continue;
                for (int p = 0; p < np; ++p) {
                    const auto ni = un * static_cast<std::size_t>(np) +
                                    static_cast<std::size_t>(p);
                    const auto lp = ul * static_cast<std::size_t>(np) +
                                    static_cast<std::size_t>(p);
                    scratch_[static_cast<std::size_t>(p)] = nodes.conc[ni];
                    links.conc[lp] = nodes.conc[ni];
                }
                if (age) {
                    scratch_[static_cast<std::size_t>(L.age_row)] =
                        ws.node_age[un];
                    ws.link_age[ul] = ws.node_age[un];
                }
                if (heat) {
                    scratch_[static_cast<std::size_t>(L.temp_row)] =
                        hs.node_temp[un];
                    hs.link_temp[ul] = hs.node_temp[un];
                }
                if (L.msx_first >= 0) {  // L3
                    const int nsp = ns - L.msx_first;
                    for (int sp = 0; sp < nsp; ++sp) {
                        const double c = ctx.reactions.msx_node_conc[
                            un * static_cast<std::size_t>(nsp) +
                            static_cast<std::size_t>(sp)];
                        scratch_[static_cast<std::size_t>(
                            L.msx_first + sp)] = c;
                        ctx.reactions.msx_link_conc[
                            ul * static_cast<std::size_t>(nsp) +
                            static_cast<std::size_t>(sp)] = c;
                    }
                }
                addToLedgerRate(dn, q * dt, scratch_.data(), ns);
            }
        }
 
        // ---- 3. RELEASE new front segments --------------------------------
        for (int l = 0; l < nl; ++l) {
            const auto ul = static_cast<std::size_t>(l);
            if (links.type[ul] != LinkType::CONDUIT) continue;
            const int up = upstreamNode(ctx, l);
            const double need = links.volume[ul] - store_.total_volume(l);
            release_vol_[ul] = (need > 0.0) ? need : 0.0;  // X3b: RWPT's V_in
            if (need <= 0.0 || up < 0) continue;
            const auto uu = static_cast<std::size_t>(up);
            for (int p = 0; p < np; ++p)
                scratch_[static_cast<std::size_t>(p)] =
                    nodes.conc[uu * static_cast<std::size_t>(np) +
                               static_cast<std::size_t>(p)];
            if (age)
                scratch_[static_cast<std::size_t>(L.age_row)] =
                    ws.node_age[uu];
            if (heat)
                scratch_[static_cast<std::size_t>(L.temp_row)] =
                    hs.node_temp[uu];
            if (L.msx_first >= 0) {  // L3
                const int nsp = ns - L.msx_first;
                for (int sp = 0; sp < nsp; ++sp)
                    scratch_[static_cast<std::size_t>(
                        L.msx_first + sp)] =
                        ctx.reactions.msx_node_conc[
                            uu * static_cast<std::size_t>(nsp) +
                            static_cast<std::size_t>(sp)];
            }
            store_.push_front(l, need, scratch_.data());
        }
 
        // ---- 3b. RWPT dispersion (X3b) — on the substep's FINAL segment
        //      field, per link, resolved vertical shear + walk. The age row
        //      disperses with the water like every other species — mixing
        //      moves age, physically. H7b: so does TEMPERATURE, at the same
        //      coefficient as a solute — the ARD engine's deliberate choice
        //      for its temperature row, adopted here for cross-engine
        //      consistency (decision 2026-08-30). ---------------------------
        if (ctx.options.lard_rwpt) {
            ++substep_counter_;
            const auto& cond = ctx.link_subtypes.conduits;
            for (int l = 0; l < nl; ++l) {
                const auto ul = static_cast<std::size_t>(l);
                if (links.type[ul] != LinkType::CONDUIT) continue;
                const double q = std::abs(links.flow[ul]);
                if (q <= kTinyFlow) continue;
                const int row = ctx.link_subtypes.conduit_row(l);
                if (row < 0) continue;
                const auto ur = static_cast<std::size_t>(row);
                const double len = cond.length[ur];
                const double vol = links.volume[ul];
                if (len <= 0.0 || vol <= 0.0) continue;
                const double a_flow = vol / len;
                const double ubar = q / a_flow;
                const double h = links.depth[ul];
                const bool circ =
                    links.xsect_shape[ul] == XsectShape::CIRCULAR;
                const double rh =
                    rwpt_hyd_radius(a_flow, h, links.xsect_geom1[ul], circ);
                rwpt_.disperse(ctx, store_, l, ubar, h, rh,
                               cond.roughness[ur], release_vol_[ul],
                               q * dt, dt, substep_counter_,
                               static_cast<std::uint64_t>(
                                   ctx.options.rwpt_seed),
                               scratch_);
            }
        }
 
        // ---- 4. DECAY (exact exponential; species-major stripes) ----------
        for (int p = 0; p < np; ++p) {
            const double k = ctx.pollutants.k_decay[static_cast<std::size_t>(p)];
            if (k == 0.0) continue;
            const double f = std::exp(-k * dt);
            double removed = store_.decay_species(p, f);
            for (int n = 0; n < nn; ++n) {
                const auto idx = static_cast<std::size_t>(n) *
                                     static_cast<std::size_t>(np) +
                                 static_cast<std::size_t>(p);
                const double v = nodes.volume[static_cast<std::size_t>(n)];
                removed += nodes.conc[idx] * (1.0 - f) * v;
                nodes.conc[idx] *= f;
            }
            if (static_cast<std::size_t>(p) <
                ctx.mass_balance.qual_routing_reacted.size())
                ctx.mass_balance.qual_routing_reacted[
                    static_cast<std::size_t>(p)] += removed;
        }
 
        // ---- 4b. REACT (L3): MSX species integrate per SEGMENT (pipe
        //      scope, the segment's own pollutant rows as context) and
        //      per NODE store (tank scope, the node's HRT) — D-L1's
        //      gather/scatter through the shared integrator. Pollutant
        //      kdecay is NOT re-applied here: stage 4 owns it (the
        //      integrator handles only MSX kinetics), so there is no
        //      double-decay the way the legacy binding had to avoid.
        if (L.msx_first >= 0) {
            auto& rx = ctx.reactions;
            const int nsp = ns - L.msx_first;
            if (static_cast<int>(msx_block_.size()) < nsp)
                msx_block_.assign(static_cast<std::size_t>(nsp), 0.0);
            if (static_cast<int>(msx_poll_.size()) < np && np > 0)
                msx_poll_.assign(static_cast<std::size_t>(np), 0.0);
            for (int l = 0; l < nl; ++l) {
                const auto ul = static_cast<std::size_t>(l);
                if (links.type[ul] != LinkType::CONDUIT) continue;
                const int cnt = store_.count(l);
                for (int i = 0; i < cnt; ++i) {
                    for (int sp = 0; sp < nsp; ++sp)
                        msx_block_[static_cast<std::size_t>(sp)] =
                            store_.seg_conc(l, i, L.msx_first + sp);
                    for (int p = 0; p < np; ++p)
                        msx_poll_[static_cast<std::size_t>(p)] =
                            store_.seg_conc(l, i, p);
                    // TEMP: the segment's own temperature row when heat
                    // rides the store (H7b); NaN defers to the
                    // [REACTION_OPTIONS] TEMPERATURE constant.
                    const double seg_temp_c =
                        (L.temp_row >= 0)
                            ? store_.seg_conc(l, i, L.temp_row)
                            : std::numeric_limits<double>::quiet_NaN();
                    transport::reactSpeciesBlock(
                        ctx, /*tank=*/false, dt, msx_block_.data(),
                        np > 0 ? msx_poll_.data() : nullptr, 0.0,
                        seg_temp_c);
                    for (int sp = 0; sp < nsp; ++sp)
                        store_.set_seg_conc(l, i, L.msx_first + sp,
                                            msx_block_[
                                                static_cast<std::size_t>(
                                                    sp)]);
                }
            }
            for (int n = 0; n < nn; ++n) {
                const auto un = static_cast<std::size_t>(n);
                for (int sp = 0; sp < nsp; ++sp)
                    msx_block_[static_cast<std::size_t>(sp)] =
                        rx.msx_node_conc[un * static_cast<std::size_t>(
                                                  nsp) +
                                         static_cast<std::size_t>(sp)];
                for (int p = 0; p < np; ++p)
                    msx_poll_[static_cast<std::size_t>(p)] =
                        nodes.conc[un * static_cast<std::size_t>(np) +
                                   static_cast<std::size_t>(p)];
                const double hrt =
                    (un < nodes.hrt.size()) ? nodes.hrt[un] : 0.0;
                const double node_temp_c =
                    (L.temp_row >= 0 &&
                     un < ctx.heat_state.node_temp.size())
                        ? ctx.heat_state.node_temp[un]
                        : std::numeric_limits<double>::quiet_NaN();
                transport::reactSpeciesBlock(
                    ctx, /*tank=*/true, dt, msx_block_.data(),
                    np > 0 ? msx_poll_.data() : nullptr, hrt,
                    node_temp_c);
                for (int sp = 0; sp < nsp; ++sp)
                    rx.msx_node_conc[un * static_cast<std::size_t>(nsp) +
                                     static_cast<std::size_t>(sp)] =
                        msx_block_[static_cast<std::size_t>(sp)];
            }
        }
    }

    void publish(SimulationContext& ctx, double dt_routing) {
        // H7b: the fourth layout-aware site joins rowLayout() (H7a
        // converted step/substep/init and flagged this one).
        const SpeciesRowLayout L = rowLayout(ctx);
        const int np = L.np;
        const bool age  = (L.age_row >= 0);
        const bool heat = (L.temp_row >= 0);
        const int nl = ctx.n_links();
        auto& nodes = ctx.nodes;
        auto& links = ctx.links;
        auto& ws = ctx.water_age_state;
        auto& hs = ctx.heat_state;
 
        // ---- 5. PUBLISH ---------------------------------------------------
        for (int l = 0; l < nl; ++l) {
            const auto ul = static_cast<std::size_t>(l);
            if (links.type[ul] != LinkType::CONDUIT) continue;
            store_.mean_conc(l, scratch_.data());
            for (int p = 0; p < np; ++p)
                links.conc[ul * static_cast<std::size_t>(np) +
                           static_cast<std::size_t>(p)] =
                    scratch_[static_cast<std::size_t>(p)];
            // Age: volume-weighted mean over the same segments, seconds.
            // The dry-element report mask (A2b) is at the report boundary,
            // engine-side — state keeps aging here regardless. An EMPTY
            // slab holds no parcels to average, so the link's held age
            // ages in place instead of resetting to 0 — the state/report
            // separation the dry-mask round established (a dry element's
            // STATE keeps aging; only the report masks it).
            if (age) {
                if (store_.count(l) > 0)
                    ws.link_age[ul] =
                        scratch_[static_cast<std::size_t>(L.age_row)];
                else
                    ws.link_age[ul] += dt_routing;
            }
            // Temperature: volume-weighted mean over the same segments. An
            // EMPTY slab HOLDS its temperature — unlike age it does not
            // grow, and unlike the age report it is not masked when dry
            // (0 degC is an ordinary temperature; the no-mask call is
            // documented at the snapshot builder).
            if (heat) {
                if (store_.count(l) > 0)
                    hs.link_temp[ul] =
                        scratch_[static_cast<std::size_t>(L.temp_row)];
            }
            // L3: species rows publish like temperature — an EMPTY
            // slab HOLDS its concentration.
            if (L.msx_first >= 0 && store_.count(l) > 0) {
                const int nsp = L.ns - L.msx_first;
                for (int sp = 0; sp < nsp; ++sp)
                    ctx.reactions.msx_link_conc[
                        ul * static_cast<std::size_t>(nsp) +
                        static_cast<std::size_t>(sp)] =
                        scratch_[static_cast<std::size_t>(
                            L.msx_first + sp)];
            }
        }
        // conc_old bookkeeping matches the ARD/legacy convention.
        links.conc_old = links.conc;
        nodes.conc_old = nodes.conc;
    }
 
private:
    void applyFluxesAndBed(SimulationContext& ctx, double dt) {
        const SpeciesRowLayout L = rowLayout(ctx);
        const bool heat_on = (L.temp_row >= 0);
        const bool bed_on  = transport::heat::bedExchangeEnabled(ctx);
        const bool flux_on =
            heat_on && (ctx.heat_config.surface_exchange ||
                        ctx.heat_config.radiative_exchange);
        if ((!flux_on && !bed_on) || !(dt > 0.0)) return;
        namespace th = transport::heat;
 
        if (flux_on) th::updateSolarForcing(ctx);
        if (bed_on) th::seedBedTemperature(ctx);
        const double t_gr = bed_on ? th::groundTemperature(ctx) : 0.0;
        const double rho = ctx.options.water_density;
        const double cp  = ctx.options.water_specific_heat;
        constexpr double kSqFt = 0.09290304;
        constexpr double kCuFt = 0.028316846592;
 
        auto& hs = ctx.heat_state;
        const int nl = ctx.n_links();
        const int nn = ctx.n_nodes();
 
        // ---- Nodes: storage surfaces only, HeatFluxes.cpp's convention. --
        if (flux_on) {
            const int unit_sys = ucf::getUnitSystem(
                static_cast<int>(ctx.options.flow_units));
            for (int n = 0; n < nn; ++n) {
                const auto un = static_cast<std::size_t>(n);
                if (un >= hs.node_temp.size()) break;
                const double vol = ctx.nodes.volume[un];
                if (!(vol > 0.0)) continue;
                const double a = node::getSurfArea(
                    ctx.nodes, n, ctx.nodes.depth[un], &ctx.tables,
                    unit_sys, &ctx.node_subtypes);
                if (!(a > 0.0)) continue;
                const double t = hs.node_temp[un];
                const HeatElement ne = HeatElement::node(n);
                hs.node_temp[un] += th::relaxT(
                    th::netFluxOut(ctx, ne, t),
                    th::netFluxOut(ctx, ne, t + th::kProbeC), th::kProbeC,
                    a * kSqFt, vol * kCuFt, dt, rho, cp);
            }
        }
 
        // ---- Links: segment field, uniform increments. -------------------
        // Bed rows: pollutants then MSX — every row before age/temp.
        const int n_bed = L.np + ((L.msx_first >= 0) ? (L.ns - L.msx_first)
                                                     : 0);
        auto& bed = ctx.bed_state;
        if (bed_on && n_bed > 0 &&
            (bed.n_species != n_bed ||
             bed.link_conc.size() != static_cast<std::size_t>(nl) *
                                         static_cast<std::size_t>(n_bed))) {
            bed.link_conc.assign(static_cast<std::size_t>(nl) *
                                     static_cast<std::size_t>(n_bed),
                                 0.0);
            bed.n_species = n_bed;
        }
 
        for (int l = 0; l < nl; ++l) {
            const auto ul = static_cast<std::size_t>(l);
            if (ctx.links.type[ul] != LinkType::CONDUIT) continue;
            const int cnt = store_.count(l);
            if (cnt <= 0) continue;
            const double vol_ft3 = store_.total_volume(l);
            if (!(vol_ft3 > 0.0)) continue;
 
            store_.mean_conc(l, scratch_.data());
            // PE1: every parcel of a conduit shares the conduit's
            // attributes (D-PE1) — the bed and the shading belong to the
            // pipe, not to the water passing through it.
            const HeatElement le = HeatElement::link(l);
 
            // Heat: one pair/relaxation against the mean, uniform dT.
            if (heat_on) {
                const double t_mean =
                    scratch_[static_cast<std::size_t>(L.temp_row)];
                const double surf_m2 =
                    flux_on ? th::linkFreeSurfaceFt2(ctx, l) * kSqFt : 0.0;
                double d_tw = 0.0;
                bool bed_stepped = false;
                if (bed_on && ul < bed.link_temp.size()) {
                    // PE2: this link's own bed material and boundary.
                    const auto& sd = th::sedimentFor(ctx, le);
                    const th::BedCoupling g = th::bedCouplingFromContact(
                        ctx, sd, th::linkBedAreaM2(ctx, l), vol_ft3,
                        th::groundTempFor(ctx, sd, t_gr));
                    if (g.viable()) {
                        const th::PairStep ps = th::relaxPair(
                            g, t_mean, bed.link_temp[ul],
                            flux_on ? th::netFluxOut(ctx, le, t_mean) : 0.0,
                            flux_on ? th::netFluxOut(ctx, le,
                                                     t_mean + th::kProbeC)
                                    : 0.0,
                            flux_on ? th::kProbeC : 0.0, surf_m2, dt);
                        d_tw = ps.dt_w;
                        bed.link_temp[ul] += ps.dt_b;
                        bed_stepped = true;
                    }
                }
                if (!bed_stepped && flux_on && surf_m2 > 0.0) {
                    d_tw = th::relaxT(
                        th::netFluxOut(ctx, le, t_mean),
                        th::netFluxOut(ctx, le, t_mean + th::kProbeC),
                        th::kProbeC, surf_m2, vol_ft3 * kCuFt, dt, rho, cp);
                }
                if (d_tw != 0.0)
                    for (int i = 0; i < cnt; ++i)
                        store_.set_seg_conc(
                            l, i, L.temp_row,
                            store_.seg_conc(l, i, L.temp_row) + d_tw);
            }
 
            // Solutes: pair against the mean, uniform dc, per bed row.
            if (bed_on && n_bed > 0) {
                // PE2: same per-link config the heat pair above used, so the
                // two halves of one bed cannot disagree about its material.
                const auto& sd = th::sedimentFor(ctx, le);
                const double bed_m2 = th::linkBedAreaM2(ctx, l);
                const double vol_b = bed_m2 * sd.bed_thickness;
                const double q_exch = th::bedExchangeQ(sd, bed_m2);
                if (q_exch > 0.0 && vol_b > 0.0) {
                    for (int b = 0; b < n_bed; ++b) {
                        // Bed row b: pollutants 0..np-1 map directly, MSX
                        // rows at msx_first + (b - np).
                        const int row = (b < L.np)
                                            ? b
                                            : L.msx_first + (b - L.np);
                        const double c_mean =
                            scratch_[static_cast<std::size_t>(row)];
                        double& cb = bed.link_conc[
                            ul * static_cast<std::size_t>(n_bed) +
                            static_cast<std::size_t>(b)];
                        const th::SolutePairStep st = th::exchangePair(
                            c_mean, cb, vol_ft3 * kCuFt, vol_b, q_exch,
                            dt);
                        if (st.dc_w != 0.0)
                            for (int i = 0; i < cnt; ++i)
                                store_.set_seg_conc(
                                    l, i, row,
                                    store_.seg_conc(l, i, row) + st.dc_w);
                        cb += st.dc_b;
                    }
                }
            }
        }
    }
 
    void init(SimulationContext& ctx) {
 
        const SpeciesRowLayout L = rowLayout(ctx);
        const int np = L.np;
        const int ns = L.ns;
        // Derived from the layout rather than re-read from options, so the
        // "is there a reserved row" question has exactly one answer per call.
        const bool age  = (L.age_row >= 0);
        const bool heat = (L.temp_row >= 0);
        const int nn = ctx.n_nodes();
        const int nl = ctx.n_links();
        // X3a: slab capacity from [OPTIONS] MAX_SEGMENTS_PER_LINK, floored
        // at 2 (one segment to hold, one to receive).
        store_.resize(nl, ns,
                      std::max(2, ctx.options.max_segments_per_link));
        flow_sign_.assign(static_cast<std::size_t>(nl), 1);
        node_mass_in_.assign(
            static_cast<std::size_t>(nn) * static_cast<std::size_t>(ns), 0.0);
        node_vol_in_.assign(static_cast<std::size_t>(nn), 0.0);
        scratch_.assign(static_cast<std::size_t>(ns), 0.0);
        treat_cin_.assign(static_cast<std::size_t>(np > 0 ? np : 0), 0.0);
        node_out_links_.assign(static_cast<std::size_t>(nn), {});
        release_vol_.assign(static_cast<std::size_t>(nl), 0.0);
        if (ctx.options.lard_rwpt) rwpt_.resize(nl);
 
        // L3: size + seed the shared MSX element state (GLOBAL fill +
        // [REACTION_QUALITY] overrides) with the SAME spelling the
        // LEGACY dispatch uses, then let the seeded link values seed
        // the segments below like every other row.
        if (L.msx_first >= 0) transport::ensureMsxState(ctx);
 
        // X4 age seeding, the ARD precedent: a hotstart-loaded state wins
        // (node_age/link_age already carry the restored values, A2a);
        // otherwise a configured INITIAL_STATE age fills the network.
        // Restored/seeded link ages then seed the segments below — the
        // same within-link profile collapse A2a recorded for ARD
        // (lesson 37): continuous, not bit-continuous.
        if (age) {
            auto& ws = ctx.water_age_state;
            if (ws.node_age.size() != static_cast<std::size_t>(nn))
                ws.resize(nn, nl, ctx.n_subcatches());
            if (!ws.hotstart_loaded) {
                const double a0 = ctx.water_age_config.global_age[
                    static_cast<int>(WaterAgeSource::INITIAL_STATE)];
                if (a0 > 0.0) {
                    std::fill(ws.node_age.begin(), ws.node_age.end(), a0);
                    std::fill(ws.link_age.begin(), ws.link_age.end(), a0);
                }
            }
            // E-A3: [INITIAL_QUALITY] __WATER_AGE__ rows override the
            // global fill; the helper no-ops under hotstart (D-IQ7). The
            // per-link values then seed the segments below.
            transport::applyInitialAgeOverrides(ctx);
        }
 
        // H7b temperature seeding, the routeLegacyHeat convention: size the
        // state if the loaders have not already, fill with the configured
        // INITIAL_STATE once (legacy_seeded — shared with the LEGACY mirror
        // so a fallback path never re-seeds), then let per-element
        // [INITIAL_QUALITY] __TEMPERATURE__ rows override. The seeded link
        // temperatures seed the segments below, the same within-link
        // profile collapse the age row records (lesson 37): continuous,
        // not bit-continuous. Hotstart does NOT restore temperature — no
        // engine's does, the record has no field for it (owed; see the
        // H7b round record) — so a restarted run re-seeds from
        // INITIAL_STATE exactly as LEGACY and ARD do.
        if (heat) {
            auto& hstate = ctx.heat_state;
            const double t0 = ctx.heat_config.global_temp[
                static_cast<int>(HeatSource::INITIAL_STATE)];
            if (hstate.node_temp.size() != static_cast<std::size_t>(nn))
                hstate.resize(nn, nl, t0);
            if (!hstate.legacy_seeded) {
                std::fill(hstate.node_temp.begin(), hstate.node_temp.end(),
                          t0);
                std::fill(hstate.link_temp.begin(), hstate.link_temp.end(),
                          t0);
                transport::applyInitialTempOverrides(ctx);
                hstate.legacy_seeded = true;
            }
        }
 
        // Seed: one segment per conduit at the link's current volume and
        // (initQuality-seeded) concentration — a dry link seeds nothing.
        for (int l = 0; l < nl; ++l) {
            const auto ul = static_cast<std::size_t>(l);
            store_.clear_link(l);
            if (ctx.links.type[ul] != LinkType::CONDUIT) continue;
            const double v = ctx.links.volume[ul];
            if (v <= 0.0) continue;
            for (int p = 0; p < np; ++p)
                scratch_[static_cast<std::size_t>(p)] =
                    ctx.links.conc[ul * static_cast<std::size_t>(np) +
                                   static_cast<std::size_t>(p)];
            if (age)
                scratch_[static_cast<std::size_t>(L.age_row)] =
                    ctx.water_age_state.link_age[ul];
            if (heat)
                scratch_[static_cast<std::size_t>(L.temp_row)] =
                    ctx.heat_state.link_temp[ul];
            if (L.msx_first >= 0) {  // L3
                const int nsp = ns - L.msx_first;
                for (int sp = 0; sp < nsp; ++sp)
                    scratch_[static_cast<std::size_t>(
                        L.msx_first + sp)] =
                        ctx.reactions.msx_link_conc[
                            ul * static_cast<std::size_t>(nsp) +
                            static_cast<std::size_t>(sp)];
            }
            store_.push_front(l, v, scratch_.data());
            flow_sign_[ul] = (ctx.links.flow[ul] >= 0.0) ? 1 : -1;
        }
        computeTopoOrder(ctx);
        initialized_ = true;
    }
 
    int upstreamNode(const SimulationContext& ctx, int l) const {
        const auto ul = static_cast<std::size_t>(l);
        return (flow_sign_[ul] >= 0) ? ctx.links.node1[ul]
                                     : ctx.links.node2[ul];
    }
    int downstreamNode(const SimulationContext& ctx, int l) const {
        const auto ul = static_cast<std::size_t>(l);
        return (flow_sign_[ul] >= 0) ? ctx.links.node2[ul]
                                     : ctx.links.node1[ul];
    }
 
    void addToLedger(int n, double vol, const double* mass, int np) {
        const auto un = static_cast<std::size_t>(n);
        node_vol_in_[un] += vol;
        for (int p = 0; p < np; ++p)
            node_mass_in_[un * static_cast<std::size_t>(np) +
                          static_cast<std::size_t>(p)] +=
                mass[static_cast<std::size_t>(p)];
    }
    void addToLedgerRate(int n, double vol, const double* conc, int np) {
        const auto un = static_cast<std::size_t>(n);
        node_vol_in_[un] += vol;
        for (int p = 0; p < np; ++p)
            node_mass_in_[un * static_cast<std::size_t>(np) +
                          static_cast<std::size_t>(p)] +=
                vol * conc[static_cast<std::size_t>(p)];
    }

    void computeTopoOrder(SimulationContext& ctx) {
        const int nn = ctx.n_nodes();
        const int nl = ctx.n_links();
        std::vector<int> indeg(static_cast<std::size_t>(nn), 0);
        for (auto& v : node_out_links_) v.clear();
        for (int l = 0; l < nl; ++l) {
            const auto ul = static_cast<std::size_t>(l);
            const int up = upstreamNode(ctx, l);
            const int dn = downstreamNode(ctx, l);
            if (up < 0 || dn < 0) continue;
            node_out_links_[static_cast<std::size_t>(up)].push_back(l);
            if (std::abs(ctx.links.flow[ul]) > kTinyFlow)
                indeg[static_cast<std::size_t>(dn)] += 1;
        }
        topo_.clear();
        topo_.reserve(static_cast<std::size_t>(nn));
        std::vector<int> q;
        for (int n = 0; n < nn; ++n)
            if (indeg[static_cast<std::size_t>(n)] == 0) q.push_back(n);
        std::vector<char> seen(static_cast<std::size_t>(nn), 0);
        std::size_t qi = 0;
        while (qi < q.size()) {
            const int n = q[qi++];
            if (seen[static_cast<std::size_t>(n)]) continue;
            seen[static_cast<std::size_t>(n)] = 1;
            topo_.push_back(n);
            for (const int l : node_out_links_[static_cast<std::size_t>(n)]) {
                if (std::abs(ctx.links.flow[static_cast<std::size_t>(l)]) <=
                    kTinyFlow)
                    continue;
                const int dn = downstreamNode(ctx, l);
                if (dn >= 0 && --indeg[static_cast<std::size_t>(dn)] == 0)
                    q.push_back(dn);
            }
        }
        for (int n = 0; n < nn; ++n)  // cycle leftovers
            if (!seen[static_cast<std::size_t>(n)]) topo_.push_back(n);
    }
 
    SegmentStore store_;
    RwptDispersion rwpt_;                
    std::vector<double> release_vol_;    
    std::uint64_t substep_counter_ = 0;  
    std::vector<int> topo_;
    std::vector<std::vector<int>> node_out_links_;
    std::vector<std::int8_t> flow_sign_;
    std::vector<double> node_mass_in_;  
    std::vector<double> node_vol_in_;   
    std::vector<double> scratch_;       
    std::vector<double> treat_cin_;
    std::vector<double> msx_block_;     
    std::vector<double> msx_poll_;      
    bool initialized_ = false;
};
 
}  // namespace lard
}  // namespace openswmm
 
#endif  // OPENSWMM_QUALITY_LARD_LAGRANGIAN_SOLVER_HPP