DynamicWave.hpp

Go to the documentation of this file.

 
#ifndef OPENSWMM_DYNAMIC_WAVE_HPP
#define OPENSWMM_DYNAMIC_WAVE_HPP
 
#include "XSectBatch.hpp"
#include "../core/Constants.hpp"
#include "../core/SimulationOptions.hpp"
#include "../data/NodeData.hpp"
#include "../data/LinkData.hpp"
#include <cstdint>
#include <functional>
#include <vector>
 
namespace openswmm {
 
struct SimulationContext;
 
namespace dynwave {
 
// ============================================================================
// Constants — imported from global Constants.hpp
// ============================================================================
 
using constants::OMEGA;
using constants::DEFAULT_HEAD_TOL;
using constants::DEFAULT_MAX_TRIALS;
using constants::MAX_VELOCITY;
using constants::MIN_TIMESTEP;
using constants::EXTRAN_CROWN_CUTOFF;
using constants::SLOT_CROWN_CUTOFF;
using constants::SLOT_WIDTH_FACTOR;
using constants::FUDGE;
using constants::MIN_SURFAREA;
 
// ============================================================================
// Dynamic Preissmann Slot (DPS) configuration and per-link state
// Sharior, Hodges & Vasconcelos (2023), J. Hydraul. Eng. 149(11)
// ============================================================================
 
struct DPSConfig {
    double c_pT   = 25.0;    
    double alpha   = 3.0;    
    double r       = 0.5;    
    double c_pT_sq = 625.0;  
};
 
struct DPSLinkArrays {
    std::vector<double>  As;            
    std::vector<double>  hs;            
    std::vector<double>  hs_prev_iter;  
    std::vector<double>  P;             
    std::vector<double>  P_hat;         
    std::vector<double>  P_hat_0;       
    std::vector<double>  T_s_target;    
    std::vector<double>  t_s;           
    std::vector<uint8_t> surcharged;    
 
    void resize(std::size_t n) {
        As.assign(n, 0.0);
        hs.assign(n, 0.0);
        hs_prev_iter.assign(n, 0.0);
        P.assign(n, 1.0);
        P_hat.assign(n, 1.0);
        P_hat_0.assign(n, 1.0);
        T_s_target.assign(n, 0.0);
        t_s.assign(n, 0.0);
        surcharged.assign(n, 0);
    }
};
 
// ============================================================================
// Per-node extended state for DW iterations
// ============================================================================
 
struct DWNodeArrays {
    std::vector<double>  new_surf_area;
    std::vector<double>  old_surf_area;   
    std::vector<double>  sumdqdh;
    std::vector<double>  dYdT;
    std::vector<uint8_t> converged;
    std::vector<uint8_t> is_surcharged;   
 
    void resize(std::size_t n) {
        new_surf_area.assign(n, 0.0);
        old_surf_area.assign(n, 0.0);
        sumdqdh.assign(n, 0.0);
        dYdT.assign(n, 0.0);
        converged.assign(n, 0);
        is_surcharged.assign(n, 0);
    }
};
 
// ============================================================================
// DW solver — batch-oriented
// ============================================================================
 
enum class SurchargeMethod : int {
    EXTRAN = 0,  
    SLOT   = 1,  
    DYNAMIC_SLOT = 2   
};
 
enum class MomentumCategory : uint8_t {
    SKIP_DRY          = 0,  
    MANNING_OPEN      = 1,  
    MANNING_CLOSED_FS = 2,  
    MANNING_CLOSED_FULL = 3,
    FORCE_MAIN_HW     = 4,  
    FORCE_MAIN_DW     = 5,  
    N_CATEGORIES      = 6
};
 
class DWSolver {
public:
    void init(int n_nodes, int n_links, const XSectGroups& groups,
              const SimulationContext& ctx);
 
    void setNumThreads(int n);
 
    using NonConduitFlowFunc = std::function<void(SimulationContext&, double, int)>;
 
    int execute(SimulationContext& ctx, double dt,
                NonConduitFlowFunc non_conduit_fn = nullptr);
 
    double getRoutingStep(SimulationContext& ctx,
                          double fixed_step, double courant_factor);
 
    double head_tol   = DEFAULT_HEAD_TOL;
    int    max_trials = DEFAULT_MAX_TRIALS;
    double omega      = OMEGA;
    SurchargeMethod surcharge_method = SurchargeMethod::EXTRAN;
    NodeContinuity  node_continuity  = NodeContinuity::EXPLICIT;
    bool   anderson_accel = false;       
 
    double evap_rate  = 0.0;
 
private:
    int n_nodes_ = 0;
    int n_links_ = 0;
    int n_conduits_ = 0;
    int num_threads_ = 1;  
    const XSectGroups* groups_ = nullptr;
 
    // Pre-built conduit index list for skipping non-conduits in inner loops
    std::vector<int> conduit_idx_;
 
    // Pre-computed per-link invariants (populated once at first execute, reused)
    // Using uint8_t instead of bool to allow .data() pointer access for SIMD/restrict
    std::vector<uint8_t> is_open_;        
    std::vector<uint8_t> is_force_main_;  
    std::vector<uint8_t> has_losses_;     
    std::vector<double> barrels_d_;       
    std::vector<double> cached_length_;   
    std::vector<double> inv_length_;      
 
    // ------------------------------------------------------------------------
    // Phase A — Conduit-dense "hot tile" of timestep-invariant data.
    //
    // Sized n_conduits_, accessed by ci (0..n_conduits_-1) for dense linear
    // memory access pattern. This replaces sparse `links.X[uj]` /
    // `nodes.X[un]` reads inside the Picard inner loops, where uj/un are
    // sparse-indexed via conduit_idx_ + links.node1/node2.  Each ci-indexed
    // read maps to one contiguous cache line that holds 8+ conduits' data,
    // versus the sparse pattern where each uj read can miss into a new line.
    //
    // All fields below are populated once in init() (or refreshConduitTile
    // when a hot-start changes invariants) and remain constant for the rest
    // of the simulation.  None of these change per Picard iter, per
    // timestep, or per outfall update.
    // ------------------------------------------------------------------------
    std::vector<int>    tile_uj_;            
    std::vector<int>    tile_n1_;            
    std::vector<int>    tile_n2_;            
    std::vector<double> tile_inv1_elev_;     
    std::vector<double> tile_inv2_elev_;     
    std::vector<double> tile_z1_off_;        
    std::vector<double> tile_z2_off_;        
    std::vector<double> tile_y_full_;        
    std::vector<double> tile_a_full_;        
    std::vector<double> tile_r_full_;        
    std::vector<double> tile_w_max_;         
    std::vector<double> tile_length_;        
    std::vector<double> tile_inv_length_;    
    std::vector<double> tile_links_length_;  
    std::vector<double> tile_beta_;          
    std::vector<double> tile_q_max_;         
    std::vector<double> tile_rough_factor_;  
    std::vector<double> tile_barrels_d_;     
    std::vector<uint8_t> tile_is_open_;      
    std::vector<uint8_t> tile_is_force_main_;
    std::vector<uint8_t> tile_is_closed_;    
    std::vector<uint8_t> tile_has_losses_;   
    std::vector<int>     tile_xsect_batch_shape_; 
    std::vector<XsectShape> tile_shape_;     
    std::vector<uint8_t> tile_has_offset_;   
    std::vector<int>     tile_uj_to_ci_;
    std::vector<int>     tile_culvert_code_;
    std::vector<double>  tile_slope_;
    std::vector<double>  tile_q_limit_;
    std::vector<double>  tile_loss_inlet_;
    std::vector<double>  tile_loss_outlet_;
    std::vector<uint8_t> tile_has_flap_gate_;
    std::vector<int8_t>  tile_direction_;
 
    void refreshConduitTile(const SimulationContext& ctx);
 
    // Per-timestep constants
    double dt_gravity_ = 0.0;            
 
    double min_surf_area_ = constants::MIN_SURFAREA;
 
    // Pre-allocated width-capping buffers (avoids thread_local per-call allocation)
    std::vector<double> wcap_d1_, wcap_d2_, wcap_dm_;
 
    // Variable timestep state (matching legacy VariableStep in dynwave.c)
    mutable double variable_step_ = 0.0;
 
    // Per-node working state (SoA)
    DWNodeArrays 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 bypass flag (true when both end nodes converged; skip momentum solve)
    // uint8_t instead of bool: avoids std::vector<bool> bit-packing overhead
    std::vector<uint8_t> bypassed_;
 
    // 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_;        
 
    // Anderson acceleration state (per-node, depth-2 mixing)
    std::vector<double> aa_y_prev_;     
    std::vector<double> aa_g_prev_;     
    std::vector<double> aa_r_prev_;     
    std::vector<uint8_t> aa_skip_;      
 
    // Per-conduit momentum category (rebuilt each Picard iteration).
    // solveMomentumBatch dispatches on category_[uj] inline — no auxiliary
    // per-category index list is needed.
    std::vector<MomentumCategory> category_;
 
    // Internal methods
    void initNodeStates(SimulationContext& ctx);
    void findBypassedLinks(const SimulationContext& ctx);
    void computeLinkGeometry(SimulationContext& ctx);
    void solveMomentumBatch(SimulationContext& ctx, double dt, int step);
    void classifyMomentumCategories(SimulationContext& ctx);
 
    void processDryLink(SimulationContext& ctx, double dt, std::size_t uj);
    void processManningLink(SimulationContext& ctx, double dt, int step,
                            std::size_t uj, MomentumCategory cat);
    void processForceMainLink(SimulationContext& ctx, double dt, int step,
                              std::size_t uj, MomentumCategory cat);
    void applyFlowLimits(SimulationContext& ctx, double dt, int step,
                         std::size_t uj, double& q, double qLast,
                         double barrels_d, bool isFull);
    void updateNodeFlows(SimulationContext& ctx, bool conduits_only = false);
    void computeAASkipFlags(const 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;
 
public:
    uint8_t& nodeSurchargedFlag(int idx) { return xnode_.is_surcharged[static_cast<std::size_t>(idx)]; }
 
    double* nodeNewSurfAreaDataMut() { return xnode_.new_surf_area.data(); }
 
    double& nodeSumDqdh(int n) { return xnode_.sumdqdh[static_cast<std::size_t>(n)]; }
 
    const std::vector<uint8_t>& aaSkipFlags() const { return aa_skip_; }
 
    const DPSLinkArrays& dpsState() const { return dps_; }
    const DPSConfig& dpsConfig() const { return dps_config_; }
    DPSLinkArrays& dpsStateMut() { return dps_; }
private:
 
    // 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;
 
    // Dynamic Preissmann Slot (DPS) state and methods
    DPSConfig dps_config_;
    DPSLinkArrays dps_;  
    double sim_time_ = 0.0;               
 
    void applyDPSGeometry(SimulationContext& ctx);
 
    void updateDPSState(SimulationContext& ctx, double dt);
 
    void spatialSmoothP(const SimulationContext& ctx);
};
 
} // namespace dynwave
} // namespace openswmm
 
#endif // OPENSWMM_DYNAMIC_WAVE_HPP