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