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 "gw/GwTransportData.hpp"   // U4
#include "subsurface/SubsurfaceSolver.hpp"   // G1: the two-zone GW kernel
#include "infil/Infil2D.hpp"
#include "mesh/RainfallInterpolator.hpp"
 
#include <memory>
#include <unordered_map>
#include <unordered_set>
#include <string>
#include <vector>
#ifdef OPENSWMM_HAS_2D
#include "solver/ISurfaceSolver.hpp"
#endif
 
namespace openswmm {
struct SimulationContext;
}
 
namespace openswmm::twoD {

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 step(SimulationContext& ctx, double dt, double t);

    void updateOutfallsPreRouting(SimulationContext& ctx);

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

    void prepareOneShotForcing(SimulationContext& ctx);

    void flushPendingBatch(SimulationContext& ctx);

    double reportWeight(double report_ms) const noexcept;

    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_; }

    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_; }

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

    SolverOptions2D& options() noexcept { return options_; }

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

    BoundaryData& boundary() noexcept { return boundary_; }

    void invalidateBoundaryIndex() noexcept { bc_nonwall_dirty_ = true; }

    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_; }

    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_;
    }

    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_;
    }
 
    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_;
    }

    using PendingBoundaryRow = twoD::PendingBoundaryRow;
    using PendingInitialQualityRow = twoD::PendingInitialQualityRow;
    using PendingBoundaryQualityRow = twoD::PendingBoundaryQualityRow;
    std::vector<PendingInitialQualityRow>& pendingInitialQualityRows() noexcept {
        return pending_iq_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:
    void drainPendingRows();
 
    MeshData         mesh_;
    SurfaceStateData state_;
    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_;

    SubsurfaceConfig         aquifer_cfg_;
    std::vector<std::string> aquifer_node_names_;
    SubsurfaceSolver         subsurface_;

    std::vector<double> rain_cum_;

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

    std::vector<PendingEdgeConveyanceRow> pending_edge_conveyance_rows_;
 
    std::vector<CouplingPoint> coupling_points_;
    std::vector<CouplingPoint> node_coupling_points_;
    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);
    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;

    std::vector<double> window_outfall_accum_;
    std::vector<double> window_avail_budget_;

    void resetWindowAccumulators();
    long outfall_clamp_windows_ = 0;

    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);

    RainfallInterpolator interp_;

    std::vector<double> rain_si_;

    void resolveBoundaryValues(SimulationContext& ctx, double t);

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