DynamicWave.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_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 <string>
#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,  
    TPA    = 3   
};

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, std::vector<std::string>* warnings = nullptr);

    using NonConduitFlowFunc = std::function<void(SimulationContext&, double, int)>;

    int execute(SimulationContext& ctx, double dt,
                NonConduitFlowFunc non_conduit_fn = nullptr);

    bool lastConverged() const { return last_converged_; }

    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;       

    int    unsteady_friction = 0;
    double uf_k3 = 0.015;  

    double tpa_celerity = 100.0;

    std::vector<int>    uf_nb_up_, uf_nb_dn_;
    std::vector<int8_t> uf_sg_up_, uf_sg_dn_;   

    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;

    bool any_conduit_seep_ = false;

    bool losses_all_zero_ = true;
 
    // 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<double>  tile_loss_avg_;     
    std::vector<double>  tile_roughness_;    
    std::vector<double>  tile_fm_sbot_;      
    std::vector<double>  tile_fm_rbot_;      
    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;            
 
    bool         last_converged_   = false; 

    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_;
 
    // ------------------------------------------------------------------------
    // B2 threading — CSR node→incident-conduit adjacency for the parallel
    // node-centric flow gather (see gatherConduitNodeFlows). Built once in
    // init() from static topology.
    //
    // PROOF OF BIT-EXACTNESS vs the serial per-link scatter: legacy
    // findLinkFlows (dynwave.c:385-388) scatters conduits one link at a time
    // in ascending link index; each link updates node1's accumulators, then
    // node2's. For a given NODE, the projection of that global order is
    // simply its incident conduit links in ascending link index (node1-end
    // entry before node2-end entry when both ends touch the same node). CSR
    // entries are stored per node in exactly that order, so the per-node
    // gather performs the identical FP accumulation sequence on each
    // accumulator (inflow, outflow, sumdqdh, new_surf_area) — bit-exact at
    // any thread count, since each node is owned by exactly one thread.
    // ------------------------------------------------------------------------
    std::vector<int>     csr_row_;           
    std::vector<int32_t> csr_link_;          
    std::vector<uint8_t> csr_is_n2_;         
    std::vector<uint8_t> csr_other_outfall_; 
 
    // ------------------------------------------------------------------------
    // Virtual junctions — zero-storage, momentum-transmitting pair nodes
    // (see plans/VIRTUAL_JUNCTION_IMPLEMENTATION_PLAN.md). vjunc_ is empty for
    // models without [VIRTUAL_JUNCTIONS]; every hot-loop hook below is gated
    // on that emptiness so VJ-free models compile down to the original paths.
    // ------------------------------------------------------------------------
    struct VJuncPair {
        int node    = -1;             
        int link_a  = -1, link_b = -1;
        int up_link = -1, dn_link = -1;
        uint8_t through = 0;          
        double lambda = 0.0;          
 
        // Per-Picard-iteration cache (filled by vjPrepareIteration):
        double  sigma_j  = 1.0;       
        double  a_up_mid = 0.0;       
        double  r_up_mid = 0.0;       
        double  dq4j     = 0.0;       
        uint8_t active   = 0;         
        int8_t  dirn     = 0;         
 
        // Momentum-residual diagnostic accumulators (whole run):
        double resid_max = 0.0;
        double resid_sum = 0.0;
        long long resid_n = 0;
    };
    std::vector<VJuncPair> vjunc_;
    std::vector<int32_t>   vj_pair_n1_;
    std::vector<int32_t>   vj_pair_n2_;

    static constexpr double kVjWetSeedFrac = 0.02;
    std::vector<double>    vj_wet_floor_;

    void buildVirtualJunctionPairs(const SimulationContext& ctx);
    void vjPrepareIteration(const SimulationContext& ctx, double dt);
    void vjAccumulateResiduals(SimulationContext& ctx);
 
    // Node-dense tile of the per-node invariants setNodeDepth touches every
    // Picard iteration. setNodeDepth reads ~20 SoA arrays per node; the seven
    // step-invariant ones below otherwise cost seven separate cache streams
    // (legacy's AoS TNode record pays 1-2 lines per node for the same reads).
    // Rebuilt once per routing step in execute() — amortised over the Picard
    // iterations and automatically correct even if the editing API mutates
    // node geometry mid-run. y_crown pre-evaluates crownElev − invertElev
    // with the identical operands the per-call subtraction used, so the
    // value is bit-identical (legacy setNodeDepth recomputes it per call).
    struct NodeTile {
        double  full_depth;
        double  y_crown;       
        double  invert_elev;
        double  ponded_area;
        double  sur_depth;
        double  full_volume;
        int32_t degree;
        uint8_t is_storage;
        uint8_t is_outfall;
        uint8_t is_virtual;    
    };
    std::vector<NodeTile> node_tile_;
    // Unit system for node volume/surf-area table dispatch, hoisted from the
    // per-call ucf::getUnitSystem(options.flow_units) (options are fixed
    // during a run).
    int unit_sys_ = 0;
 
    // 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 recomputeConduitLossOne(SimulationContext& ctx, double dt, int ci);
    void momentumKernels(SimulationContext& ctx, double dt, int step);

    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);
    void buildConduitNodeCSR(const SimulationContext& ctx);
    void gatherConduitNodeFlows(SimulationContext& ctx);
    void computeAASkipFlags(const SimulationContext& ctx);
    void updateNodeDepthsTeam(SimulationContext& ctx, double dt, int step,
                              int& unconv_shared);
    void setNodeDepth(SimulationContext& ctx, int node_idx, double dt, int step);
    void commitNodeDepthState(SimulationContext& ctx, int node_idx,
                              double y_new, double dV, double dt);
    double getLinkStep(const SimulationContext& ctx, int link_idx) const;
 
public:
    bool isInitialized() const noexcept { return !xnode_.sumdqdh.empty(); }

    uint8_t& nodeSurchargedFlag(int idx) { return xnode_.is_surcharged[static_cast<std::size_t>(idx)]; }

    double* nodeNewSurfAreaDataMut() { return xnode_.new_surf_area.data(); }

    bool isBypassed(int j) const {
        const auto uj = static_cast<std::size_t>(j);
        return uj < bypassed_.size() && bypassed_[uj] != 0;
    }

    double& nodeSumDqdh(int n) { return xnode_.sumdqdh[static_cast<std::size_t>(n)]; }

    void setNodeDepthForTest(SimulationContext& ctx, int node_idx, double dt,
                             int step) {
        setNodeDepth(ctx, node_idx, dt, step);
    }

    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 getCrownCutoff() const;
 
    // Dynamic Preissmann Slot (DPS) state and methods
    DPSConfig dps_config_;
    DPSLinkArrays dps_;  
    double sim_time_ = 0.0;               

    void applyDPSGeometry(SimulationContext& ctx);
 
    // -- TPA (SurchargeMethod::TPA, issue #156 Phase 5) ----------------------
    std::vector<double> tpa_w_;
    std::vector<uint8_t> tpa_latch_;
    std::vector<uint8_t> tpa_latch_changed_;

    void updateTpaLatch(SimulationContext& ctx);

    void applyTpaGeometry(SimulationContext& ctx);

    void updateDPSState(SimulationContext& ctx, double dt);

    void spatialSmoothP(const SimulationContext& ctx);
};
 
} // namespace dynwave
} // namespace openswmm
 
#endif // OPENSWMM_DYNAMIC_WAVE_HPP