ExplicitInertialSolver.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_ENGINE_2D_EXPLICIT_INERTIAL_SOLVER_HPP
#define OPENSWMM_ENGINE_2D_EXPLICIT_INERTIAL_SOLVER_HPP
 
#include <array>
#include <functional>
#include <cstdint>
#include <string>
#include <vector>
 
#include "ISurfaceSolver.hpp"
#include "InertialEdges.hpp"
#include "MomentumPhysics2D.hpp"
#include "../gw/GwSourceForcing.hpp"
#include "../../transport/components/ReactionModule/ReactionIntegrator.hpp"
#include "../data/SolverOptions2D.hpp" // Momentum2D
 
namespace openswmm::detail {
class CompletedSourceArchive;
}
namespace openswmm {
class SWMMEngine;
struct SimulationContext;
struct Table;
}
namespace openswmm::runoff {
class SourceWaterDriver;
}
namespace openswmm::twoD {
struct Rk2CheckpointBinding;
 
class SubsurfaceSolver; // G1: the two-zone groundwater kernel, if authored
class CommonSourceStage;
class ChannelGroundState;
class SourceResultHistory;
class EngineSourceResultHistory;
 
class ExplicitInertialSolver final : public ISurfaceSolver {
public:
    void publishMomentumVelocity();

    void initialize(MeshData& mesh, SurfaceStateData& state, SolverOptions2D& opts) override;
    double advance(double t_current, double t_target) override;
    double completedSourceTime() const noexcept override {
        return initialized_ ? t_last_sync_ : -1.0;
    }

    void reinitialize(double t0) override;
    void resyncFromVolumes(double t0) override;
    void pinExternalSourceCells(const std::vector<int>& cells);
    void setCommonSourceStage(CommonSourceStage& stage);
    void setChannelGround(ChannelGroundState&);
    void setSourceResultHistory(SourceResultHistory& history);
    std::string saveCheckpoint(const std::string& path) const;
    std::string loadCheckpoint(const std::string& path);
    void resumeSourceResultHistory(SourceResultHistory& history);
    std::string configureRk2Reactions(const SimulationContext&);
    std::string configureRk2BoundarySignals(SimulationContext&, double stage_scale,
                                            double flow_scale);
    void resolveRk2BoundarySignals(double time);
    double surfaceReactionLoss(int row) const noexcept {
        return row >= 0 && std::size_t(row) < rk_reacted_.size() ? rk_reacted_[row] : 0;
    }

    void finalize() override;
    void refreshBoundaries();

    long last_num_steps() const noexcept override { return last_steps_; }
    double last_step_size() const noexcept override { return last_dt_; }
    RunStats run_stats() const noexcept override;
    bool is_initialized() const noexcept override { return initialized_; }
    const std::vector<double>& last_coupling_exchange() const noexcept override { return exch_; }

    double spillDeficit() const noexcept { return spill_deficit_; }
 
private:
    friend struct ExplicitInertialSolverTestAccess;
    friend class SurfaceRouter2D;
    friend class openswmm::SWMMEngine;
    std::string checkpointBytes(std::string& bytes) const;
    std::string rk2CheckpointProfile() const;
    std::string rk2CheckpointModel() const;
    std::string rk2CheckpointModel(const Rk2CheckpointBinding*) const;
    std::string rk2CheckpointBytes(std::string&) const;
    std::string loadRk2CheckpointBytes(std::string);
    std::string loadRk2CheckpointBytes(
        std::string, const Rk2CheckpointBinding*,
        const std::function<void(const ExplicitInertialSolver&, const SurfaceStateData&,
                                 const SubsurfaceSolver&)>*,
        const std::function<void()>*, bool inspect);
    std::string loadCheckpointBytes(
        std::string bytes, const std::array<double, 10>* router_boundary = nullptr,
        const std::vector<double>* pending_outfalls = nullptr,
        const std::function<void(const runoff::SourceWaterDriver&)>* validate_source = nullptr,
        const std::function<void()>* before_install = nullptr,
        const std::vector<int>* boundary_slots = nullptr, double* inspect_source_end = nullptr);
    void refreshSourceResultBoundary();
    double advanceImpl(double current, double target);
    std::string rk2ThermalProfileError() const;
    void advanceRk2Thermal(double seconds, double sample_time);
    void advanceRk2BedContact(double seconds);
    void advanceRk2Snow(double seconds, double start);
    void advanceRk2SurfaceReactions(double seconds);
    void syncRk2ReactionReport(SimulationContext&);
    const ReactionData* rk_reactions_ = nullptr; // immutable model config
    transport::RxWorkspace rk_reaction_workspace_;
    std::vector<int> rk_reaction_rows_;
    std::vector<double> rk_decay_, rk_reacted_, rk_reacted_reported_;
    std::string rk_reaction_identity_;
 
    // Recompute η/depth from state volumes for the whole mesh.
    void reconstructAll();
    // Flush every pending face accumulator into its cell (CSR gather over ALL
    // cells) and refresh the closure of touched cells. MUST run before any
    // tier/active-set reassignment: a cell deactivated or re-tiered with a
    // pending side-accumulator strands flux whose counterpart the other side
    // already applied — the backstop then realizes it as created volume.
    void settleAccumulators();
    // Apply lazily-accumulated sources (rain/held coupling) to INACTIVE cells
    // over [t_last_sync_, t], rebuild the active cell/edge lists, and assign
    // the LTS tiers (dt0_ = finest active CFL requirement).
    void syncAndRebuild(double t);
    // Tighten-only dt0_ refresh from CURRENT depths/speeds between rebuilds
    // (dt0_ may only grow at syncAndRebuild, which reassigns the tiers).
    void refreshDt0();
    // Fire one tier's faces over their Δt: inertial update + Froude cap +
    // face-cadence positivity share, booking ±ΔM into both side accumulators.
    // global_step = every active cell fires at this same Δt (the window tail),
    // so a face's exporter republishes its volume every time this face does
    // and the positivity share is NOT divided down by a refire ratio.
    void fireFaces(const std::vector<int>& faces, double dt_f, bool global_step = false);
    // Fire one tier's cells over their Δt: gather + clear own-side face
    // accumulators, apply sources, refresh closure + Perot vector; tier-0
    // firings also evaluate the availability-clamped boundary edges.
    // tier0 = this firing carries the tier-0 cadence work (boundary edges and
    // the live junction exchange), which fires once per finest substep.
    void fireCells(const std::vector<int>& cells, double dt_c, bool tier0);
    template <bool kSwe, bool kFuse>
    void fireCellsImpl(const std::vector<int>& cells, double dt_c, bool tier0);
    // One halving-order macro cycle of nsub base substeps: tier k fires every
    // 2^k substeps.
    void runMacroCycle(double dt0, int nsub, double time);
 
    MeshData* mesh_ = nullptr;
    SurfaceStateData* state_ = nullptr;
    SolverOptions2D* opts_ = nullptr;
 
    void refreshConveyance(); 
    InertialEdges edges_;
    std::vector<double> q_;
    std::vector<double> qcx_, qcy_;
 
    // -----------------------------------------------------------------------
    // Momentum closure (MOMENTUM_EQUATION, 2D_FULL_SWE plan §2.1). One
    // marcher, three face laws; everything below the face/cell kernels
    // (tiers, active sets, positivity, coupling, transport) is shared.
    // -----------------------------------------------------------------------
    Momentum2D mode_ = Momentum2D::LOCAL_INERTIAL;
    bool has_quads_ = false;
    std::vector<int> wall_ptr_, wall_slot_;
    // Momentum packets share face_accum_ slots 2..5 in FULL_SWE.
    double swe_bc_dqx_ = 0.0, swe_bc_dqy_ = 0.0;
    bool second_order_ = false;
    bool rk_transport_stage_ = false;
    std::vector<double> rk_stage_mass_, rk_stage_volume_; // scratch beginning-state donor views
    std::vector<double> gex_, gey_, gux_, guy_, gvx_, gvy_;
    std::vector<double> rk_v0_, rk_qx0_, rk_qy0_;
    std::vector<double> rk_inf0_, rk_bc0_, rk_ex0_, rk_drawn0_, rk_cpl0_;
    void computeLimitedGradientsSwe();
    MomentumPhysics2D momentum_physics_;
    std::vector<int> wall_drag_slots_;
    double wall_drag_energy_ = 0; 
    void applyMomentumPhysics(double seconds);
    double runMomentumStep(double limit, double time);

    void runRk2Step(double dt, double time);
    void runRk2Step(double dt, double time, unsigned thermal_depth);
    bool surfaceLatentCooling() const noexcept;
    struct Rk2BoundarySignal {
        int slot = -1, table = -1;
        bool stage = false;
        double scale = 1;
        std::string name;
        GwSourceSignal signal;
        std::shared_ptr<const Table> curve{};
        double x_scale = 1;
    };
    std::vector<Rk2BoundarySignal> rk_boundary_signals_;
    double rk_boundary_anchor_ = 0;

    bool front_rebuild_ = false;
    std::vector<uint8_t> frontier_;
    bool front_breach_ = false;
    double h_on_front_ = 0.0; 
    std::vector<double> dw_slope_;
    void fireFacesInertial(const std::vector<int>& faces, double dt_f, bool global_step);
    void fireFacesSwe(const std::vector<int>& faces, double dt_f, bool global_step);
    void fireFacesDiffusive(const std::vector<int>& faces, double dt_f, bool global_step);
    void bookFaceSpecies(int e, int a, int b, double dM, double hf, double dt_f,
                         bool global_step) noexcept;
    double boundaryFluxSwe(std::size_t k, int i, double dt_c,
                           const std::array<double, 4>* held = nullptr);
    void refreshDiffusiveSlopes();
    std::vector<uint8_t> cell_active_;
    std::vector<int> active_cells_;
    std::vector<int> active_faces_;
    std::vector<uint8_t> pin_t0_; 
 
    // Tiered LTS. tier_[i] = k means cell i updates every 2^k base substeps
    // with Δt = 2^k·dt0; face tier = min of its incident cells so a face
    // always integrates at the finer cadence. Every face firing books the
    // identical ±ΔM into the two logical water sides (single writer per face); each cell
    // applies + clears its own side at its own firing — conservation across
    // tier interfaces is exact by construction.
    std::vector<uint8_t> rebuild_seed_;
    std::vector<double> rebuild_dt_cell_;
    std::vector<double> rebuild_dt_partial_;
 
    std::vector<uint8_t> tier_;      
    std::vector<uint8_t> face_tier_; 
    std::vector<double> face_accum_;
    std::size_t face_accum_stride_ = 2; 
 
    // S1 — species mass rides the SAME face accumulators, with each face
    // holding adjacent rows at [e * n_species + s]. Booked in fireFaces
    // immediately after the volume ΔM, from the FINAL qn1 (after the Froude
    // cap and the positivity share),
    // at the exporting cell's concentration read at that same substep;
    // gathered and cleared in fireCells alongside the volume side. That is
    // D-2DT2: the species flux inherits the volume flux's tier cadence rather
    // than reproducing it, so conservation across tier interfaces is the
    // marcher's own property and not a second one to prove. Empty unless
    // `state_->transport.active()`.
    std::vector<double> sacc_L_;
    std::vector<double> sacc_R_;
    bool species_on_ = false;
    double donorConc(int s, int cell) const noexcept;
    void sinkMassAtCellConc(int i, double dv_m3, std::vector<double>& ledger,
                            double* per_point_ledger = nullptr) noexcept;
    // Sources share one water budget after face transfers have landed. Gross
    // constituent transfers are evaluated even when net water change is zero.
    // Returns false only when no source bookkeeping or volume change occurred.
    bool applyCellSources(int cell, double dt, double* boundary_donor = nullptr);
    bool applyPlainCellSources(int cell, double dt);
    // Selected only during an advance; external rebuild/settlement uses the full path.
    bool plain_sources_ = false;
    bool applyRuntimeSources(int cell, double dt);
    struct CellSourceStep {
        double area_dt = 0, rain = 0, coupling = 0, available = 0, infil = 0, evap = 0, out = 0,
               final_volume = 0, requested_infil = 0;
        double* ledger = nullptr;
    };
    bool prepareCellSources(int cell, double dt, CellSourceStep& source);
    void mixCellSourceRow(int cell, int row, double& mass, const CellSourceStep&);
    void prepareChannelGround(double seconds);
    void applyCellSourceRow(int cell, int row, double& mass, const CellSourceStep& source);
    enum class SourceLedger {
        Infiltration,
        CouplingOut,
        Rainfall,
        CouplingIn,
        Exfiltration,
        Boundary,
        ExternalIn,
        ExternalOut,
        BedExchange,
        Count
    };
    void bookSourceLedger(SourceLedger ledger, int species, double mass) noexcept;
    void flushSourceLedgers();
    int sourceThread() const noexcept;
    // Each worker owns a cache-line-padded slice: six species ledgers and
    // applied evaporation volume. Fold only after the parallel region ends.
    std::vector<double> source_ledgers_;
    std::size_t source_ledger_stride_ = 0;
    ChannelGroundState* channel_ground_ = nullptr;
    std::vector<CellSourceStep> channel_sources_;
    std::vector<uint8_t> channel_source_ready_, channel_temperature_mixed_;
    bool accumulators_pending_ = false;
    std::vector<std::vector<int>> cells_by_tier_;
    std::vector<std::vector<int>> edges_by_tier_;
    double dt0_ = 0.0;         
    std::vector<int> bc_cell_; 
    std::vector<int> bc_slot_; 
    std::vector<int> gw_row_map_;
    std::vector<double> bc_accum_; 
    std::vector<double> bc_q_;     
 
    // Live junction exchange (windowless coupling): state_->node_coupling
    // points, evaluated at tier-0 cadence against live 2D heads and the
    // routing-step 1D heads. exch_[k] = ∫Q_k dt (m³, + = 2D→1D), reset per
    // advance; node_drawn_ caps a step's total spill at the node's stored
    // volume so fill-and-spill thrash is structurally impossible.
    std::vector<double> exch_;
    std::vector<double> node_drawn_; 
    double spill_deficit_ = 0.0;
 
public:
    void setExchangeHeadSlopes(std::vector<double> slopes) { exch_head_slope_ = std::move(slopes); }
 
private:
    std::vector<double> exch_head_slope_;
    double exch_tau_ = 0.0; 
 
    double t_last_sync_ = 0.0; 
    CommonSourceStage* common_stage_ = nullptr;
    SourceResultHistory* source_history_ = nullptr;
    EngineSourceResultHistory* engine_source_history_ = nullptr;
    class EngineSourceResultStream* engine_source_stream_ = nullptr;
    bool checkpoint_ready_ = false;
    bool joint_checkpoint_ = false; // only SurfaceRouter2D may restore live 1D pairing
    std::string checkpointBoundary() const;
    void checkpointState(detail::CompletedSourceArchive& archive);
    void captureSourceResults();
    double source_phase_end_ = 0.0;
    int cycles_since_rebuild_ = 1000; 
    void lazySourcesOnly(double t);
    long substeps_run_ = 0; 
    long face_passes_ = 0;  
    long last_steps_ = 0;   
    double last_dt_ = 0.0;
    bool initialized_ = false;
 
    // Active-fraction telemetry: (sim time, active cells) sampled at every
    // rebuild; dumped as CSV at finalize() when OPENSWMM_2D_MARCHER_TELEMETRY
    // names a file. The Phase-1 gate reads this to verify the thin-film budget
    // assumption on the Bellinge storm slice.
    std::vector<std::pair<double, int>> telemetry_; 
    std::string telemetry_path_;
    double active_frac_min_ = 1.0e30;
    double active_frac_max_ = -1.0e30;
    double active_frac_sum_ = 0.0;
    long active_samples_ = 0;
    std::array<long, 8> tier_occupancy_{};

    SubsurfaceSolver* gw_ = nullptr;
 
public:
    void setSubsurface(SubsurfaceSolver* gw) noexcept { gw_ = gw; }
    void refreshSubsurfaceRows();

    const InertialEdges& inertialEdges() const noexcept { return edges_; }
};
 
} // namespace openswmm::twoD
 
#endif // OPENSWMM_ENGINE_2D_EXPLICIT_INERTIAL_SOLVER_HPP