SurfaceRouter2D.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_SURFACE_ROUTER_HPP
#define OPENSWMM_ENGINE_2D_SURFACE_ROUTER_HPP
 
#include "data/MeshData.hpp"
#include "data/SurfaceStateData.hpp"
#include "data/SolverOptions2D.hpp"
#include "data/BoundaryData.hpp"
#include "data/PendingRows2D.hpp"
#include "coupling/NodeCoupling.hpp"
#include "coupling/OutfallExchange.hpp"
#include "gw/GwTransportData.hpp"
#include "quality/SurfaceQuality2D.hpp"    // S7
#include "subsurface/SubsurfaceSolver.hpp" // G1: the two-zone GW kernel
#include "heat/ChannelGround.hpp"
#include "infil/Infil2D.hpp"
#include "mesh/RainfallInterpolator.hpp"
 
#include <memory>
#include <functional>
#include <array>
#include <unordered_map>
#include <unordered_set>
#include <string>
#include <vector>
#ifdef OPENSWMM_HAS_2D
#include "solver/ISurfaceSolver.hpp"
#endif
 
namespace openswmm::runoff {
class SourceWaterDriver;
}
namespace openswmm {
struct SimulationContext;
struct HeatState;
namespace climate {
struct ClimateState;
}
class Router;
class SWMMEngine;
}
 
namespace openswmm::twoD {
class CommonSourceStage;
class SourceResultHistory;

class SurfaceRouter2D {
public:
    SurfaceRouter2D() = default;
    ~SurfaceRouter2D() = default;
 
    // Non-copyable, movable
    SurfaceRouter2D(const SurfaceRouter2D&) = delete;
    SurfaceRouter2D& operator=(const SurfaceRouter2D&) = delete;
    SurfaceRouter2D(SurfaceRouter2D&&) = default;
    SurfaceRouter2D& operator=(SurfaceRouter2D&&) = default;

    void initialize(SimulationContext& ctx);
    void initializeChannelGround(SimulationContext& ctx);
    const ChannelGroundState& channelGround() const noexcept { return channel_ground_; }
    ChannelGroundState& channelGround() noexcept { return channel_ground_; }
    void setCommonSourceStage(SimulationContext& ctx, CommonSourceStage& stage);
    void setJointSourceStage(SimulationContext& ctx, CommonSourceStage& stage, Router& hydraulic);
    void setSourceResultHistory(SourceResultHistory& history);
    void resumeSourceResultHistory(SourceResultHistory& history);
    std::string saveCheckpoint(const SimulationContext& ctx, const std::string& path) const;
    std::string loadCheckpoint(SimulationContext& ctx, const std::string& path);

    void step(SimulationContext& ctx, double dt, double t);

    void updateOutfallsPreRouting(SimulationContext& ctx, double dt = 0.0);

    void advancePostRouting(SimulationContext& ctx, double dt, double t);

    void prepareOneShotForcing(SimulationContext& ctx);

    void flushPendingBatch(SimulationContext& ctx);

    double reportWeight(double report_ms) const noexcept;

    void bookGageRates(const SimulationContext& ctx, double runoff_start_s);

    void reportRainfall(const SimulationContext& ctx, double report_date,
                        std::vector<double>& out) const;

    const std::vector<double>& reportRainCumulative(double report_ms) const noexcept {
        return report_rain_ms_ == report_ms ? report_rain_cum_ : rain_cum_;
    }

    const std::vector<double>& reportOldDepth() const noexcept { return report_old_depth_; }
    const std::vector<double>& reportOldHead() const noexcept { return report_old_head_; }

    void finalize(SimulationContext& ctx);

    bool isActive() const noexcept { return active_; }
    void prepareLidStep(SimulationContext&, double dt, double evaporation);

    const std::vector<std::string>& threadWarnings() const noexcept { return thread_warnings_; }

    void computeCouplingConductances(const SimulationContext& ctx,
                                     std::vector<std::pair<int, double>>& out) const;

    void prepareForEdit();

    const MeshData& mesh() const noexcept { return mesh_; }

    MeshData& mesh() noexcept { return mesh_; }

    const SurfaceStateData& state() const noexcept { return state_; }

    SurfaceStateData& state() noexcept { return state_; }

    double surfaceReactionLoss(int row) const noexcept;

    const SolverOptions2D& options() const noexcept { return options_; }

    SolverOptions2D& options() noexcept { return options_; }

    const BoundaryData& boundary() const noexcept { return boundary_; }

    BoundaryData& boundary() noexcept { return boundary_; }

    bool runtimeCpuAvailable() const noexcept;
    RuntimeSources& runtimeSources() noexcept { return state_.runtime_sources; }
    bool commonSourceAttached() const noexcept { return common_stage_ != nullptr; }
    const CommonSourceStage* commonSourceStage() const noexcept { return common_stage_; }
    bool runtimeWaterCheckpointing() const noexcept;
    bool runtimeBoundaryCheckpointing() const noexcept;
    bool outfallTailwaterCheckpointing() const noexcept {
        return joint_hydraulic_ && !outfall_exchange_.points_.empty();
    }

    const OutfallExchange& outfallExchange() const noexcept { return outfall_exchange_; }
    std::string validateRuntimeWater(const RuntimeSources&, const RuntimeSources&,
                                     const RuntimeForcings&, double time,
                                     const BoundaryData* boundary = nullptr) const;
    bool runtimeActive() const noexcept {
        return state_.runtime_sources.active() || subsurface_.runtimeSources().active() ||
               !boundary_.runtime_original.empty() || state_.runtime_forcings.active();
    }

    double nextRuntimeTime(double now) const noexcept {
        return std::min({state_.runtime_sources.nextTime(now),
                         subsurface_.runtimeSources().nextTime(now),
                         state_.runtime_forcings.nextTime(now)});
    }

    void expireRuntime(double now) {
        state_.runtime_sources.expire(now);
        subsurface_.runtimeSources().expire(now);
        state_.runtime_forcings.expire(now);
    }

    void invalidateBoundaryIndex();
    void invalidateBoundaryNames() noexcept { boundary_names_resolved_ = false; }

    void refreshRenderFieldsIfStale() {
        if (render_dirty_)
            refreshRenderFields();
    }

    Infil2D& infil() noexcept { return infil_; }
    const Infil2D& infil() const noexcept { return infil_; }

    GwTransportData& gwTransport() noexcept { return gw_; }
    const GwTransportData& gwTransport() const noexcept { return gw_; }
    SurfaceQuality2D& surfaceQuality() noexcept { return surface_quality_; }
    const SurfaceQuality2D& surfaceQuality() const noexcept { return surface_quality_; }

    SubsurfaceConfig& aquiferConfig() noexcept { return aquifer_cfg_; }
    const SubsurfaceConfig& aquiferConfig() const noexcept { return aquifer_cfg_; }
    std::vector<std::string>& aquiferNodeNames() noexcept { return aquifer_node_names_; }
    const std::vector<std::string>& aquiferNodeNames() const noexcept {
        return aquifer_node_names_;
    }

    std::vector<std::string>& aquiferLinkNames() noexcept { return aquifer_link_names_; }
    const std::vector<std::string>& aquiferLinkNames() const noexcept {
        return aquifer_link_names_;
    }

    SubsurfaceSolver& subsurface() noexcept { return subsurface_; }
    const SubsurfaceSolver& subsurface() const noexcept { return subsurface_; }

    const std::vector<double>& subcatchRecharge() const noexcept { return subcatch_recharge_; }
    void drainSubcatchRecharge(std::vector<double>& out) {
        out = subcatch_recharge_;
        std::fill(subcatch_recharge_.begin(), subcatch_recharge_.end(), 0.0);
    }

    const std::vector<int>& cellSubcatchment() const noexcept { return cell_subcatch_; }

    double infiltrationElapsed() const noexcept { return infil_elapsed_; }
    void restoreInfiltration(const std::vector<double>& cumulative, double elapsed) {
        infil_cum_applied_ = cumulative;
        infil_elapsed_ = elapsed;
    }

    const std::vector<double>& infilCumulative() const noexcept { return infil_cum_applied_; }

    const std::vector<double>& couplingCumulative() const noexcept { return coupling_cum_applied_; }

    const std::vector<double>& rainCumulative() const noexcept { return rain_cum_; }

    const RainfallInterpolator& rainfallInterpolator() const noexcept { return interp_; }

    using PendingBoundaryRow = twoD::PendingBoundaryRow;
    using PendingInitialQualityRow = twoD::PendingInitialQualityRow;
    using PendingInitialBedQualityRow = twoD::PendingInitialBedQualityRow;
    using PendingBoundaryQualityRow = twoD::PendingBoundaryQualityRow;
    std::vector<PendingInitialQualityRow>& pendingInitialQualityRows() noexcept {
        return pending_iq_rows_;
    }

    std::vector<PendingInitialBedQualityRow>& pendingInitialBedQualityRows() noexcept {
        return pending_ibq_rows_;
    }

    std::vector<PendingBoundaryQualityRow>& pendingBoundaryQualityRows() noexcept {
        return pending_bq_rows_;
    }

    std::vector<PendingBoundaryRow>& pendingBCRows() noexcept { return pending_bc_rows_; }
    const std::vector<PendingBoundaryRow>& pendingBCRows() const noexcept {
        return pending_bc_rows_;
    }

    using PendingEdgeConveyanceRow = twoD::PendingEdgeConveyanceRow;
    std::vector<PendingEdgeConveyanceRow>& pendingEdgeConveyanceRows() noexcept {
        return pending_edge_conveyance_rows_;
    }

    const std::vector<PendingEdgeConveyanceRow>& pendingEdgeConveyanceRows() const noexcept {
        return pending_edge_conveyance_rows_;
    }

    double totalVolume() const;

    double totalExchangeFlow() const;
#ifdef OPENSWMM_HAS_2D
    long lastSolverSteps() const { return solver_ ? solver_->last_num_steps() : 0; }
    double lastSolverStepSize() const { return solver_ ? solver_->last_step_size() : 0.0; }
    ISurfaceSolver::RunStats runStats() const {
        return solver_ ? solver_->run_stats() : ISurfaceSolver::RunStats{};
    }
#else
    long lastSolverSteps() const { return 0; }
    double lastSolverStepSize() const { return 0.0; }
#endif
    const std::string& backendName() const noexcept { return backend_name_; }
 
private:
    friend class openswmm::SWMMEngine;
    struct RestartState;
    std::string checkpointBytes(const SimulationContext&, std::string&) const;
    std::string loadCheckpointBytes(
        SimulationContext&, std::string, const std::array<double, 6>* engine_boundary = nullptr,
        const std::vector<double>* engine_rainfall = nullptr,
        const double* engine_storage = nullptr, bool engine_routed = false,
        const std::function<void(const runoff::SourceWaterDriver&)>* validate_source = nullptr,
        const std::function<void()>* before_install = nullptr, double* inspect_source_end = nullptr,
        const climate::ClimateState* prepared_climate = nullptr,
        const HeatState* prepared_heat = nullptr);
    void attachCommonSourceStage(SimulationContext&, CommonSourceStage&, Router*);
    bool standaloneRk2Checkpointing() const noexcept;
    void attachStandaloneRk2Checkpoint(SimulationContext&);
    std::string standaloneRk2Boundary(const SimulationContext&) const;
    void configureSnowGroundContacts();
    void checkJointRoutingStep(const SimulationContext&, double dt, double start) const;
    Router* joint_hydraulic_ = nullptr; // externally owned, stable lifetime
    static double couplingSyncSpan(const SimulationContext&, const SolverOptions2D&);
    CommonSourceStage* common_stage_ = nullptr;
    bool checkpoint_owned_initialization_ = false; // scoped explicit native builder
    SimulationContext* source_context_ = nullptr;
    std::string source_router_model_;
    bool source_boundary_ready_ = false;
    std::string sourceBoundary(const SimulationContext& ctx) const;
    std::string sourceRouterModel(const SimulationContext& ctx) const;
    void drainPendingRows();
 
    MeshData mesh_;
    SurfaceStateData state_;
    ChannelGroundState channel_ground_;
    SolverOptions2D options_;
    BoundaryData boundary_;
    std::vector<std::string> thread_warnings_; 
    std::string backend_name_;                 

    Infil2D infil_;

    double infil_elapsed_ = 0.0;

    std::vector<double> infil_cum_applied_;

    std::vector<double> coupling_cum_applied_;

    std::vector<int> cell_subcatch_;
    std::vector<double> subcatch_recharge_;

    GwTransportData gw_;
    SurfaceQuality2D surface_quality_; 

    SubsurfaceConfig aquifer_cfg_;
    std::vector<std::string> aquifer_node_names_;
    std::vector<std::string> aquifer_link_names_; 
    SubsurfaceSolver subsurface_;
    std::vector<GwLinkShare> aquifer_link_shares_;
    std::vector<int> aquifer_link_ids_;

    std::vector<double> rain_cum_;

    std::vector<PendingBoundaryRow> pending_bc_rows_;
    std::vector<PendingInitialQualityRow> pending_iq_rows_; 
    std::vector<PendingInitialBedQualityRow> pending_ibq_rows_;
    std::vector<PendingBoundaryQualityRow> pending_bq_rows_; 

    std::vector<PendingEdgeConveyanceRow> pending_edge_conveyance_rows_;
 
    std::vector<CouplingPoint> coupling_points_;
    OutfallExchange outfall_exchange_;
    bool outfall_state_changed_ = false;
    std::vector<CouplingPoint> node_coupling_points_;
    std::vector<int> coupling_head_cells_;
    std::vector<int> pond_track_ptr_;
    std::vector<int> pond_track_cell_;
    std::vector<double> pond_track_share_;
    std::vector<double> pond_chamber_ft2_;
    std::vector<double> node_pond_area_m2_;
    void updateTrackingPonds(SimulationContext& ctx, double dt);
    void publishEqualizeAreas(const SimulationContext& ctx, double window);
    std::vector<double> node_row_conc_;
    void publishNodeRows(const SimulationContext& ctx);
    void applyAgingAndReactions(SimulationContext& ctx, double dt);
    std::vector<double> react_phi_, react_a_, react_dx_; 

    double co_refresh_elapsed_ = 0.0;
    double co_render_elapsed_ = 0.0;
    bool render_dirty_ = false;
    std::vector<double> render_eta_;
    std::vector<int> bc_nonwall_slots_;
    bool bc_nonwall_dirty_ = true;
    std::unordered_set<int> mb_seen_nodes_;
    double co_forcing_elapsed_ = 0.0;
    bool co_forcing_first_ = true;
    void coAdvanceStep(SimulationContext& ctx, double dt, double t);
    void publishLinkCoupling(SimulationContext& ctx);
    double linkRowConc(const SimulationContext& ctx, int link, int s,
                       const SubsurfaceTransportState& tr) const;
    void flushAquiferNodeSpecies(SimulationContext& ctx);
    bool refreshDue(const SimulationContext& ctx, double dt) const;
    void refreshRenderFields();
 
    bool active_ = false;
    double sim_time_ = 0.0;
    double pending_dt_ = 0.0;

    std::vector<double> report_old_depth_;
    std::vector<double> report_old_head_;
    double report_window_old_ms_ = 0.0;
    double report_window_new_ms_ = 0.0;
    bool report_old_valid_ = false;

    std::vector<double> ramp_prev_head_;
    double ramp_prev_span_ = 0.0;

    double prev_boundary_cum_ = 0.0;
    double prev_boundary_in_ = 0.0, prev_boundary_out_ = 0.0;

    std::vector<double> window_outfall_accum_;
    void resetWindowAccumulators();

    double last_t_ = 0.0;

    bool boundary_names_resolved_ = false;

    double bc_stage_scale_ = 1.0;

    double bc_stage_ts_scale_ = 1.0;

    double bc_flow_scale_ = 1.0;
#ifdef OPENSWMM_HAS_2D
    std::unique_ptr<ISurfaceSolver> solver_;
#endif

    void updateRainfall(SimulationContext& ctx);

    void gageRatesToCells(const SimulationContext& ctx, const std::vector<double>& rates,
                          std::vector<double>& out) const;

    RainfallInterpolator interp_;

    std::vector<double> rain_gage_last_;

    std::vector<double> gage_rate_, gage_cum_, gage_window_depth_, gage_window_rate_;
    double gage_since_ = 0.0;
    // Preserve the integral at the next report boundary before the runoff
    // clock moves beyond it. Only one gage row and one cell row are retained.
    std::vector<double> gage_report_depth_, report_rain_cum_;
    double report_rain_ms_ = -1.0;

    void resolveBoundaryValues(SimulationContext& ctx, double t);

    void accumulateMassBalance(SimulationContext& ctx, double dt);
};
 
} // namespace openswmm::twoD
 
#endif // OPENSWMM_ENGINE_2D_SURFACE_ROUTER_HPP