DynamicWave.hpp

Go to the documentation of this file.

 
#ifndef OPENSWMM_DYNAMIC_WAVE_HPP
#define OPENSWMM_DYNAMIC_WAVE_HPP
 
#include "XSectBatch.hpp"
#include "../data/NodeData.hpp"
#include "../data/LinkData.hpp"
#include <vector>
 
namespace openswmm {
 
struct SimulationContext;
 
namespace dynwave {
 
// ============================================================================
// Constants (matching legacy)
// ============================================================================
 
constexpr double OMEGA               = 0.5;       
constexpr double DEFAULT_HEADTOL     = 0.005;     
constexpr int    DEFAULT_MAXTRIALS   = 8;          
constexpr double MAXVELOCITY         = 50.0;      
constexpr double MINTIMESTEP         = 0.001;      
constexpr double EXTRAN_CROWN_CUTOFF = 0.96;       
constexpr double SLOT_CROWN_CUTOFF   = 0.985257;   
constexpr double SLOT_WIDTH_FACTOR   = 0.001;       
 
// ============================================================================
// Per-node extended state for DW iterations
// ============================================================================
 
struct DWNodeState {
    double new_surf_area = 0.0;
    double old_surf_area = 0.0;   
    double sumdqdh       = 0.0;
    double dYdT          = 0.0;
    bool   converged     = false;
    bool   is_surcharged = false;  
};
 
// ============================================================================
// DW solver — batch-oriented
// ============================================================================
 
enum class SurchargeMethod : int {
    EXTRAN = 0,  
    SLOT   = 1,  
    DYNAMIC_SLOT = 2   
};
 
class DWSolver {
public:
    void init(int n_nodes, int n_links, const XSectGroups& groups);
 
    void setNumThreads(int n);
 
    int execute(SimulationContext& ctx, double dt);
 
    double getRoutingStep(const SimulationContext& ctx,
                          double fixed_step, double courant_factor) const;
 
    double head_tol   = DEFAULT_HEADTOL;
    int    max_trials = DEFAULT_MAXTRIALS;
    double omega      = OMEGA;
    SurchargeMethod surcharge_method = SurchargeMethod::EXTRAN;
 
private:
    int n_nodes_ = 0;
    int n_links_ = 0;
    int num_threads_ = 1;  
    const XSectGroups* groups_ = nullptr;
 
    // Per-node working state
    std::vector<DWNodeState> xnode_;
 
    // Per-link pre-computed geometry (batch-filled by XSectGroups each iteration)
    std::vector<double> area1_;      
    std::vector<double> area2_;      
    std::vector<double> area_mid_;   
    std::vector<double> hrad_mid_;   
    std::vector<double> width_mid_;  
    std::vector<double> depth1_;     
    std::vector<double> depth2_;     
    std::vector<double> depth_mid_;  
 
    // Per-link momentum working arrays
    std::vector<double> velocity_;   
    std::vector<double> froude_;     
    std::vector<double> sigma_;      
    std::vector<double> dqdh_;       
    std::vector<double> new_flow_;   
 
    // Per-link area from previous iteration (for unsteady term)
    std::vector<double> area_old_;
 
    // Per-link surface area contributions to upstream/downstream nodes
    // (matching legacy Link[].surfArea1/surfArea2 from dwflow.c findSurfArea)
    std::vector<double> surf_area1_;   
    std::vector<double> surf_area2_;   
 
    // Per-link upstream geometry (for proper weighted hyd. radius)
    std::vector<double> hrad1_;        
    std::vector<double> width1_;       
    std::vector<double> width2_;       
 
    // Per-link head values (persisted from computeLinkGeometry for solveMomentumBatch,
    // may be modified by flow classification for UP_CRITICAL/DN_CRITICAL cases)
    std::vector<double> h1_;           
    std::vector<double> h2_;           
    std::vector<double> fasnh_;        
 
    // Internal methods
    void initNodeStates(SimulationContext& ctx);
    void computeLinkGeometry(SimulationContext& ctx);
    void solveMomentumBatch(SimulationContext& ctx, double dt, int step);
    void updateNodeFlows(SimulationContext& ctx);
    bool updateNodeDepths(SimulationContext& ctx, double dt, int step);
    void setNodeDepth(SimulationContext& ctx, int node_idx, double dt, int step);
    double getLinkStep(const SimulationContext& ctx, int link_idx) const;
 
    // Preissmann slot helpers (matching legacy dwflow.c)
    double getSlotWidth(double y, double y_full, double w_max, XsectShape shape) const;
    double getSlotArea(double y, double y_full, double a_full, double slot_width) const;
    double getSlotHydRad(double y, double y_full, double r_full) const;
    double getCrownCutoff() const;
};
 
} // namespace dynwave
} // namespace openswmm
 
#endif // OPENSWMM_DYNAMIC_WAVE_HPP