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;
class Router;
 
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<double> conduit_loss; 
    std::vector<double> loss_factor; 
    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);
        conduit_loss.assign(n, 0.0);
        loss_factor.assign(n, 1.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);

    void accumulateLinkStepStats(SimulationContext& ctx, double dt, double fixed_step,
                                 double courant_factor) const;
 
    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:
    friend class ::openswmm::Router; // explicit internal checkpoint fields only
    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_crow_;
    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_; 
    std::vector<double> q1_;
 
    // 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;
        double rpt_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 limitConduitLosses(SimulationContext& ctx, double dt);
    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 traceLinkTerms(const SimulationContext& ctx, std::size_t uj, double qLast, double v,
                        double sig, double rho, double aWtd, double rWtd, double dq1, double dq2,
                        double dq3, double dq4, double dq5, double dq6, double qOld, double q,
                        double aMidConv);
    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