Routing.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_ROUTING_HPP
#define OPENSWMM_ROUTING_HPP
 
#include "XSectBatch.hpp"
#include "../quality/FvTransportWater.hpp"
#include "KinematicWave.hpp"
#include "DynamicWave.hpp"
#include "Divider.hpp"
#include "fv/FvOptions.hpp"
#include "fv/INetworkSolver.hpp"
#include "fv/NetworkMeshData.hpp"
#include <array>
 
#include <memory>
#include <string>
#include <functional>
#include <vector>
 
namespace openswmm {
namespace twoD {
class SurfaceRouter2D;
}
 
struct SimulationContext;
struct OutfallData;
 
// ============================================================================
// Routing model enum
// ============================================================================
 
enum class RouteModel : int {
    STEADY = 0,  
    KINWAVE = 1, 
    DYNWAVE = 2, 
    FV = 3       
};

void applyConduitLengthening(SimulationContext& ctx, RouteModel model);
 
// ============================================================================
// Routing orchestrator
// ============================================================================

class Router {
public:
    void restoreFvClassicState(SimulationContext& ctx);
    void init(SimulationContext& ctx, RouteModel model);

    void acceptedQualityWater(const SimulationContext& ctx, double dt, quality::FvTransportWater& out) const;
 
    int step(SimulationContext& ctx, double dt, double evap_rate = 0.0,
             dynwave::DWSolver::NonConduitFlowFunc non_conduit_fn = nullptr);

    void enableCheckpointing(SimulationContext& ctx);
    void beginCheckpointStep(SimulationContext& ctx, double dt);
    void acceptCheckpointBoundary(SimulationContext& ctx);
    std::string saveCheckpoint(const SimulationContext& ctx, const std::string& path) const;
    std::string loadCheckpoint(SimulationContext& ctx, const std::string& path);

    double getAdaptiveStep(SimulationContext& ctx, double fixed_step, double courant);

    const XSectGroups& xsectGroups() const { return groups_; }

    bool hasCycle() const { return cycle_detected_; }

    void setDWNumThreads(int n, std::vector<std::string>* warnings = nullptr) {
        dw_solver_.setNumThreads(n, warnings);
    }

    dynwave::DWSolver& dwSolver() { return dw_solver_; }

    void setStructureSolver(hydstruct::StructureSolver* s) { structures_ = s; }
 
    // Internal water-only conduit window: copy persistent DW history and batch
    // geometry without retaining pointers into the live router. Checkpoint and
    // external structure sessions are not part of this transaction profile.
    std::unique_ptr<Router> privateDynamicWaveWindow() const;
    void publishDynamicWaveWindow(Router& r) noexcept;

    bool lastStepConverged() const {
        return (model_ == RouteModel::DYNWAVE) ? dw_solver_.lastConverged() : true;
    }

    bool enableLinkStepStats(SimulationContext& ctx);

    bool linkStepStatsEnabled() const { return lstep_enabled_; }

    void accumulateLinkStepStats(SimulationContext& ctx, double dt, bool in_window);

    const fv::INetworkSolver* fvSolver() const { return fv_solver_.get(); }

    const fv::NetworkMeshData& fvMesh() const { return fv_mesh_; }

    const std::string& fvBackend() const { return fv_backend_; }

    const std::vector<std::string>& fvWarnings() const { return fv_warnings_; }
    const std::vector<std::string>& fvErrors() const { return fv_errors_; }
 
private:
    friend class twoD::SurfaceRouter2D;
    friend class SWMMEngine;
    struct RestartState;
    std::string checkpointBytes(const SimulationContext& c, std::string& result) const;
    std::shared_ptr<RestartState>
    prepareCheckpointBytes(const SimulationContext& c, std::string bytes,
                           const std::array<double, 6>* boundary = nullptr,
                           const double* delivery = nullptr,
                           const OutfallData* feedback = nullptr) const;
    static void validateJointDelivery(const NodeData& nodes, const std::vector<uint8_t>& coupled,
                                      double delivery);
    static void validateSessionStorage(const std::shared_ptr<RestartState>& s,
                                       const SimulationContext& c, double expected, bool routed);
    void installCheckpoint(SimulationContext& c, std::shared_ptr<RestartState> s) noexcept;
    SimulationContext* checkpoint_context_ = nullptr; // non-owning, stable lifetime
    bool checkpoint_owned_session_ = false;           // engine owns structure/control memory
    bool checkpoint_aquifer_session_ =
        false; // paired engine/surface recovery owns node-bed receipts
    std::string checkpoint_model_;
    double checkpoint_time_ = 0.0, checkpoint_step_end_ = 0.0;
    bool checkpoint_ready_ = false, checkpoint_prepared_ = false, checkpoint_step_finished_ = false;
    std::string checkpointModel(const SimulationContext& ctx) const;
    void checkCheckpointProfile(const SimulationContext& ctx) const;
    void checkCheckpointStep(SimulationContext& ctx, double dt, double evap_rate, bool callback,
                             bool owned_callback = false);
    int stepImpl(SimulationContext& ctx, double dt, double evap_rate,
                 dynwave::DWSolver::NonConduitFlowFunc non_conduit_fn, bool owned_callback);
 
    RouteModel model_ = RouteModel::DYNWAVE;
    XSectGroups groups_;

    hydstruct::StructureSolver* structures_ = nullptr;

    std::vector<char> steady_storage_updated_;
    std::vector<double> steady_y_;
 
    // --- Explicit finite-volume solver (FLOW_ROUTING FV) -------------------
    // The mesh and state are owned here, not by the solver, exactly as
    // SurfaceRouter2D owns MeshData/SurfaceStateData: the solver is swappable
    // (CPU or Kokkos plugin) and must not own what survives a backend change.
    fv::NetworkMeshData fv_mesh_;
    fv::NetworkStateData fv_state_;
    fv::FvOptions fv_opts_; 
    std::unique_ptr<fv::INetworkSolver> fv_solver_;
    std::string fv_backend_;
    std::vector<std::string> fv_warnings_, fv_errors_;
 
    // Per-conduit time step summary ([REPORT] LINK_STEPS). FV samples arrive
    // through the solver's sink during advance() and wait here, per link,
    // until accumulateLinkStepStats keeps or discards them.
    bool lstep_enabled_ = false;
    double lstep_rs_ = 0.0; 
    std::vector<double> lstep_pmin_, lstep_pmax_, lstep_pdt_time_, lstep_ptime_, lstep_pbin_;
    static void fvStepSink(void* user, int conduit, double dt, double dur);
 
    // Per-step forcing buffers, allocated once at init.
    std::vector<double> fv_lateral_, fv_fixed_head_, fv_struct_flow_, fv_cond_loss_;

    std::vector<double> fv_link_q_cap_;

    std::vector<double> fv_struct_int_;
    double fv_struct_t_prev_ = 0.0;

    void initFv(SimulationContext& ctx);
    int stepFv(SimulationContext& ctx, double dt,
               dynwave::DWSolver::NonConduitFlowFunc non_conduit_fn);
    void publishFv(SimulationContext& ctx, double dt);
    kinwave::KWSolver kw_solver_;
    dynwave::DWSolver dw_solver_;
    // Divider data moved to ctx.node_subtypes.dividers (relational side-table).
    std::vector<int> steady_sorted_links_; 
    bool cycle_detected_ = false;          

    void initNodeFlows(SimulationContext& ctx, double dt, double evap_rate);

    void computeConduitLosses(SimulationContext& ctx, double dt, double evap_rate);

    void refreshFvBoundaryFlows(SimulationContext& ctx, double t_elapsed,
                                const std::function<void()>& eval_structures);

    void accumulateStructureFlows(const SimulationContext& ctx, double t_now);

    void updateLinkStates(SimulationContext& ctx);

    int executeSteadyFlow(SimulationContext& ctx, double dt);
};
 
} // namespace openswmm
 
#endif // OPENSWMM_ROUTING_HPP