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/CompletedSourceArchive.hpp"
 
#include "../../core/SimulationContext.hpp"
#include "../../transport/InitialQualitySeeds.hpp"
#include "../NegativeSources.hpp"
#include "../QualityRouting.hpp"
#include "../../transport/components/ReactionModule/ReactionLegacyBinding.hpp"
#include "../../transport/components/ReactionModule/ReactionGeometry.hpp"
#include "../../transport/TransportPolicy.hpp" // E2
#include "RwptDispersion.hpp"
#include "SegmentStore.hpp"
#include "../../transport/components/ReactionModule/ReactionSources.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:
    std::string checkpointState() const {
        auto copy = *this;
        detail::CompletedSourceArchive a;
        copy.checkpoint(a);
        return a.bytes();
    }

    void restoreCheckpointState(const std::string& bytes, SimulationContext& ctx) {
        detail::CompletedSourceArchive a(bytes);
        bool ready = false;
        a(ready);
        if (ready)
            init(ctx);
        checkpointFields(a);
        a.finish();
    }

    bool initialized() const noexcept { return initialized_; }
    std::vector<ChannelThermalRegion> thermalRegions(const SimulationContext& ctx) const {
        std::vector<ChannelThermalRegion> out;
        const auto layout = rowLayout(ctx);
        if (!initialized_ || layout.temp_row < 0)
            return out;
        std::vector<double> mean(layout.ns);
        for (int link = 0; link < ctx.n_links(); ++link) {
            const double volume = store_.total_volume(link) * .028316846592;
            if (!(volume > 0))
                continue;
            store_.mean_conc(link, mean.data());
            out.push_back({link, 0, 0, 1, volume, mean[layout.temp_row]});
            transport::heat::fillChannelThermalGeometry(ctx, out.back());
        }
        return out;
    }

    std::vector<ChannelSpeciesRegion> speciesRegions(const SimulationContext& ctx) const {
        std::vector<ChannelSpeciesRegion> out;
        const auto layout = rowLayout(ctx);
        if (!initialized_ || layout.temp_row < 0)
            return out;
        std::vector<double> mean(layout.ns);
        for (int link = 0; link < ctx.n_links(); ++link) {
            const double volume = store_.total_volume(link) * .028316846592;
            if (!(volume > 0))
                continue;
            store_.mean_conc(link, mean.data());
            ChannelSpeciesRegion region;
            region.link = link;
            region.region = 0;
            region.volume = volume;
            for (int p = 0; p < layout.np; ++p)
                region.row_names.push_back(ctx.pollutant_names.name_of(p));
            if (layout.age_row >= 0)
                region.row_names.emplace_back("__WATER_AGE__");
            if (layout.temp_row >= 0)
                region.row_names.emplace_back("__TEMPERATURE__");
            if (layout.msx_first >= 0)
                for (int row = layout.msx_first; row < layout.ns; ++row)
                    region.row_names.push_back(
                        ctx.reactions.species_name[std::size_t(row - layout.msx_first)]);
            region.concentration = mean;
            out.push_back(std::move(region));
        }
        return out;
    }

    std::string applySpeciesIncrements(SimulationContext& ctx,
                                       const std::vector<ChannelSpeciesIncrement>& updates) {
        const auto layout = rowLayout(ctx);
        if (!initialized_ || layout.temp_row < 0)
            return "LARD channel species state is unavailable.";
        for (const auto& u : updates) {
            if (u.link < 0 || u.link >= ctx.n_links() || u.region != 0 ||
                u.concentration_delta.size() != std::size_t(layout.ns))
                return "LARD species increment has an invalid channel region or row count.";
            if (u.concentration_delta[layout.temp_row] != 0)
                return "LARD species update cannot modify the temperature row.";
            std::vector<double> mean(layout.ns);
            store_.mean_conc(u.link, mean.data());
            for (int k = 0; k < store_.count(u.link); ++k)
                for (int row = 0; row < layout.ns; ++row) {
                    const double delta = u.concentration_delta[row];
                    const double target = mean[row] + delta;
                    const double next = mean[row] > 0
                                            ? store_.seg_conc(u.link, k, row) * target / mean[row]
                                            : target;
                    if (!std::isfinite(delta) || !std::isfinite(next) ||
                        (row != layout.temp_row && next < 0))
                        return "LARD species increment is nonfinite or makes a carried mass negative.";
                }
        }
        // A channel region is a whole reach in LARD. Apply the finite-contact
        // change proportionally to its current segment concentrations, which
        // preserves their pattern and keeps every segment nonnegative. The
        // scaling also makes an inverse rollback restore the original state.
        for (const auto& u : updates) {
            std::vector<double> mean(layout.ns);
            store_.mean_conc(u.link, mean.data());
            for (int k = 0; k < store_.count(u.link); ++k)
                for (int row = 0; row < layout.ns; ++row)
                    if (row != layout.temp_row) {
                        const double target = mean[row] + u.concentration_delta[row];
                        const double next =
                            mean[row] > 0 ? store_.seg_conc(u.link, k, row) * target / mean[row]
                                          : target;
                        store_.set_seg_conc(u.link, k, row, next);
                    }
        }
        publish(ctx, 0);
        return {};
    }

    std::string applyThermalIncrements(SimulationContext& ctx,
                                       const std::vector<ChannelThermalIncrement>& updates) {
        const auto layout = rowLayout(ctx);
        if (!initialized_ || layout.temp_row < 0)
            return "LARD channel thermal state is unavailable.";
        for (const auto& u : updates) {
            if (u.link < 0 || u.link >= ctx.n_links() || u.region != 0 || !std::isfinite(u.celsius))
                return "LARD thermal increment has an invalid channel region.";
            for (int k = 0; k < store_.count(u.link); ++k)
                if (!std::isfinite(store_.seg_conc(u.link, k, layout.temp_row) + u.celsius))
                    return "LARD thermal increment is nonfinite.";
        }
        for (const auto& u : updates)
            for (int k = 0; k < store_.count(u.link); ++k)
                store_.set_seg_conc(u.link, k, layout.temp_row,
                                    store_.seg_conc(u.link, k, layout.temp_row) + u.celsius);
        publish(ctx, 0);
        return {};
    }

    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);
            // The node's mixed STORE (node::storeVolume): the routers book
            // nodes.volume in legacy's convention, 0 for a junction below
            // its rim, but the store this solver mixes into is the water
            // the continuity equation attributes to the node.
            const double v_old =
                node::storeVolume(nodes, n, nodes.old_depth[un], 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.
                if (is_msx && !zero_bf)
                    m += transport::reactionBoosterMass(ctx, n, s - L.msx_first, m, denom,
                                                        v_in + nodes.qual_vol_in[un] * frac);
                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 = node::storeVolume(nodes, 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) {
            quality::updateUntreatedStorageHrt(ctx, dt); // review A6
            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, l);
                    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);
                if (!quality::nodeIsReactor(ctx, n) || nodes.volume[un] <= quality::ZERO_VOLUME)
                    continue;
                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,
                                             n);
                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 checkpoint(detail::CompletedSourceArchive& a) {
        a(initialized_);
        checkpointFields(a);
    }
    void checkpointFields(detail::CompletedSourceArchive& a) {
        if (!initialized_)
            return;
        store_.checkpoint(a);
        rwpt_.checkpoint(a);
        a(substep_counter_, topo_, node_out_links_);
        a.fixed(flow_sign_, node_mass_in_, node_vol_in_);
        const int nn = int(node_vol_in_.size()), nl = int(flow_sign_.size());
        if (topo_.size() != std::size_t(nn) || node_out_links_.size() != std::size_t(nn))
            throw std::runtime_error("Invalid LARD topology shape.");
        std::vector<bool> seen(nn, false);
        for (int n : topo_) {
            if (n < 0 || n >= nn || seen[n])
                throw std::runtime_error("Invalid LARD checkpoint topology.");
            seen[n] = true;
        }
        for (const auto& links : node_out_links_)
            for (int l : links)
                if (l < 0 || l >= nl)
                    throw std::runtime_error("Invalid LARD checkpoint link.");
        for (auto v : flow_sign_)
            if (v < -1 || v > 1)
                throw std::runtime_error("Invalid LARD flow direction.");
        for (const auto* v : {&node_mass_in_, &node_vol_in_})
            for (double x : *v)
                if (!std::isfinite(x))
                    throw std::runtime_error("Nonfinite LARD inflow receipt.");
    }
    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, dt);
        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 = node::storeVolume(ctx.nodes, n);
                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 (ul < ctx.heat_state.mesh_contact_link.size() &&
                ctx.heat_state.mesh_contact_link[ul])
                continue;
            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