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

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

    int step(SimulationContext& ctx, double dt,
             double evap_rate = 0.0,
             dynwave::DWSolver::NonConduitFlowFunc non_conduit_fn = nullptr);

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

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

    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:
    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-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