KinematicWave.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_KINEMATIC_WAVE_HPP
#define OPENSWMM_KINEMATIC_WAVE_HPP
 
#include "XSectBatch.hpp"
#include <vector>
 
namespace openswmm {
 
struct SimulationContext;
 
namespace hydstruct {
class StructureSolver;
}
 
namespace kinwave {
 
// ============================================================================
// Tree-layout routing helpers — shared by KINWAVE and STEADY
// ============================================================================
//
// Legacy keeps these in flowrout.c because Steady Flow and Kinematic Wave run
// the SAME outer loop (flowrout_execute) over a topologically sorted link
// array; only the per-conduit kernel differs. They live here so the KW loop
// and Router::executeSteadyFlow share one implementation.

double getLinkInflow(SimulationContext& ctx, hydstruct::StructureSolver* structures, int j,
                     double dt);

void updateStorageState(SimulationContext& ctx, hydstruct::StructureSolver* structures,
                        const std::vector<int>& order, int pos, int i, double dt);

void finishRouting(SimulationContext& ctx, hydstruct::StructureSolver* structures,
                   const std::vector<int>& order, const std::vector<char>& storage_updated,
                   const std::vector<double>& link_y1, const std::vector<double>& link_y2,
                   double dt);

double capConduitLoss(SimulationContext& ctx, std::size_t ucr, double q);
 
// ============================================================================
// Constants (matching legacy)
// ============================================================================
 
constexpr double WX = 0.6;      
constexpr double WT = 0.6;      
constexpr double EPSIL = 0.001; 
 
// ============================================================================
// KW solver — operates on entire conduit set
// ============================================================================

class KWSolver {
public:
    void init(int n_conduits, const XSectGroups& groups);

    void setLinkOrder(const std::vector<int>& sorted_links) { sorted_links_ = sorted_links; }

    int execute(SimulationContext& ctx, double dt,
                hydstruct::StructureSolver* structures = nullptr);

    std::vector<int> sorted_links_;

    std::vector<char> storage_updated_;

    int solveConduit(int idx, const XSectParams& xs, double q_full, double a_full, double s_full,
                     double beta, double length, double dt, double loss_rate);
 
    // Per-conduit SoA state (indexed by conduit-link index)
    // Public so that unit tests can set up / inspect working arrays directly.
    std::vector<double> q1_; 
    std::vector<double> a1_; 
    std::vector<double> q2_; 
    std::vector<double> a2_; 
 
    // Working buffers (reused each timestep)
    std::vector<double> q_in_;  
    std::vector<double> a_in_;  
    std::vector<double> q_out_; 
    std::vector<double> a_out_; 
    std::vector<double> sf_in_; 
    std::vector<double> y1_;    
    std::vector<double> y2_;    
 
private:
    int n_conduits_ = 0;
};
 
} // namespace kinwave
} // namespace openswmm
 
#endif // OPENSWMM_KINEMATIC_WAVE_HPP