ExplicitFvSolver.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_ENGINE_FV_EXPLICIT_FV_SOLVER_HPP
#define OPENSWMM_ENGINE_FV_EXPLICIT_FV_SOLVER_HPP
 
#include <array>
#include <cstdint>
#include <vector>
 
#include "FvKernels.hpp"
#include "FvOptions.hpp"
#include "INetworkSolver.hpp"
#include "NetworkMeshData.hpp"
#include "PressurizedHeadSolver.hpp"
 
namespace openswmm::transport::fvkernels {
struct SpeciesKernelView;  // transport/fvkernels/SpeciesTransportKernels.hpp
}
 
namespace openswmm::fv {
 
class ExplicitFvSolver : public INetworkSolver {
public:
    ExplicitFvSolver() = default;
    ~ExplicitFvSolver() override = default;
 
    ExplicitFvSolver(const ExplicitFvSolver&) = delete;
    ExplicitFvSolver& operator=(const ExplicitFvSolver&) = delete;
 
    void   initialize(NetworkMeshData& mesh, NetworkStateData& state,
                      const FvOptions& opts) override;
    double advance(double t_current, double t_target,
                   const FvStepForcing& forcing) override;
    void   reinitialize(double t0) override;
    void   finalize() override;
 
    long   last_num_steps() const noexcept override { return last_nsteps_; }
    double last_step_size() const noexcept override { return last_h_; }
    double suggested_step() const noexcept override { return suggested_h_; }
    RunStats run_stats() const noexcept override;
    bool   is_initialized() const noexcept override { return mesh_ != nullptr; }
    bool   divergence(Divergence& d) const noexcept override {
        if (div_link_ < 0) return false;
        d.link  = div_link_;
        d.value = div_value_;
        d.what  = div_what_;
        return true;
    }

    const std::vector<double>& node_exchange() const noexcept { return node_exch_; }

    const std::vector<double>& node_inflow_volume() const noexcept { return node_in_; }
    const std::vector<double>& node_outflow_volume() const noexcept { return node_out_; }

    const std::vector<uint8_t>& inlet_control() const noexcept {
        return inlet_control_;
    }

    const std::vector<double>& node_flood_volume() const noexcept {
        return flood_vol_;
    }

    const std::vector<double>& dummy_volume() const noexcept {
        return dummy_vol_;
    }

    const std::vector<std::uint8_t>& node_passthrough() const noexcept {
        return node_pass_;
    }

    const std::vector<std::uint8_t>& node_publish_stage() const noexcept {
        return node_pub_;
    }

    const std::vector<double>& cell_flow_integral() const noexcept {
        return cell_q_int_;
    }

    double nodeVolumeFromDepth(int node, double depth) const;

    int lts_tiers() const noexcept { return lts_tiers_; }

    const std::array<long, kMaxLtsTiers>& tier_occupancy() const noexcept {
        return tier_occupancy_;
    }
 
private:
    // -- substep pipeline ---------------------------------------------------
    void   refreshDepths();
    void   refreshNodeAreas();
    void   rebuildActiveLists();
    double censusDt(bool press_edit = false) const;
    void   reconstructState();
    void   computeFaceFlux(int face);
    void   reconstructScalars(double dt);
    transport::fvkernels::SpeciesKernelView speciesKernelView();
    void   computeFluxes();
    void   limitPositivity(double dt);

    void   relaxNodeFluxes(double dt, const FvStepForcing& forcing);

    void   relaxOneNode(int node, double dt, const FvStepForcing& forcing);

    void   solveAlgebraicNode(int node, double dt, const FvStepForcing& forcing);

    void   settleAlgebraicNode(int node, double carry);

    std::vector<std::uint8_t> node_alg_;
    std::vector<std::uint8_t> node_pass_static_;
    std::vector<std::uint8_t> node_pass_;
    std::vector<std::uint8_t> node_pub_static_;
    std::vector<std::uint8_t> node_pub_;
    std::vector<double> node_carry_;
    std::vector<double> cell_qlat_;
    std::vector<std::uint8_t> node_lat_div_;
    double nodeLateral(const FvStepForcing& forcing,
                       std::size_t un) const noexcept {
        if (!node_lat_div_.empty() && node_lat_div_[un]) return 0.0;
        return forcing.node_lateral ? forcing.node_lateral[un] : 0.0;
    }
    std::vector<double> face_a_rim_;
    std::vector<double> face_lat_spill_;
    double injectDivertedLateral(std::size_t uc, double a_new, double dt,
                                 double inv_dx);
    void   bookLateralSpills(const std::vector<int>* nodes);
    std::vector<double> node_vfull_;

    bool algebraicActive(int n) const noexcept {
        const auto un = static_cast<std::size_t>(n);
        if (!node_alg_[un]) return false;
        if (mesh_->node_can_pond[un] &&
            state_->node_volume[un] > node_vfull_[un]) return false;
        return true;
    }

    bool incidentPressurized(std::size_t un) const noexcept {
        for (int p = mesh_->node_face_ptr[un];
             p < mesh_->node_face_ptr[un + 1]; ++p) {
            const auto uf = static_cast<std::size_t>(
                mesh_->node_face_idx[static_cast<std::size_t>(p)]);
            const int c = (mesh_->face_cl[uf] >= 0) ? mesh_->face_cl[uf]
                                                    : mesh_->face_cr[uf];
            if (c < 0) continue;
            const auto uc = static_cast<std::size_t>(c);
            if (tpaCell(uc)) return true;
            const FvGeometry& g =
                mesh_->geom[static_cast<std::size_t>(mesh_->cell_geom[uc])];
            if (!g.is_open && state_->cell_h[uc] >= g.y_crown) return true;
        }
        return false;
    }
    void   updateCells(double dt, const FvStepForcing& forcing);
    void   updateNodes(double dt, const FvStepForcing& forcing);
    void   dispersionSolve(double dt);

    void   refreshStructFlows(const FvStepForcing& forcing);
    std::vector<double> node_qstruct_;

    std::vector<double> node_qdummy_;

    std::vector<double> dummy_vol_;

    template <class ArrivedFn>
    void refreshDummyFlows(double dt, const FvStepForcing& forcing,
                           ArrivedFn&& arrived);

    double nodeQDummy(std::size_t un) const noexcept {
        return node_qdummy_.empty() ? 0.0 : node_qdummy_[un];
    }
 
    double nodeDepthFromVolume(int node, double volume) const;

    bool ltsEligible() const noexcept {
        return opts_.lts &&
               opts_.time_integration != TimeIntegration::RK2 &&
               state_ && state_->n_species == 0;
    }
 
    bool faceIsLive(int f) const noexcept {
        return all_faces_live_ ||
               (!live_stamp_.empty() &&
                live_stamp_[static_cast<std::size_t>(f)] == live_gen_);
    }

    std::vector<std::uint32_t> live_stamp_;
    std::uint32_t              live_gen_ = 0;
    bool all_faces_live_ = true;

    double algebraicNodeStableDt(int n) const noexcept;
 
    double frictionFor(const FvGeometry& g, double q, double u, double h,
                       double dt) const;

    void applyNodeCapacity(int node, double v_prev, double& vol, double& depth);
 
    // -- local time stepping (plan §3.3) ------------------------------------
    double cellStableDt(int cell) const;
    double nodeStableDt(int node) const;

    int    assignTiers(double& dt0);

    void   fireFaces(const std::vector<int>& faces, double dt0);
    void   fireCells(const std::vector<int>& cells, double dt0,
                     const FvStepForcing& forcing);
    void   fireNodes(const std::vector<int>& nodes, double dt0,
                     const FvStepForcing& forcing);
    void   runMacroCycle(double dt0, int nsub, const FvStepForcing& forcing);

    void   settleAccumulators();

    void   saveState();
    void   restoreState();

    void   takeSubstep(double dt, const FvStepForcing& forcing);
 
    // -- implicit pressurized head update (slot program R2a) ----------------
    bool   anyPressurizedCell() const;
    bool   nodePressFolded(int n) const;
    PressurizedView pressView(const FvStepForcing& forcing);
    PressurizedHeadSolver press_;
    bool   press_step_ = false;

    void   rkSave();
    void   rkAverage(const FvStepForcing& forcing);

    void   faceSide(int face, int cell, int node, double zstar, int dir,
                    double u_interior, kernels::FaceState& out,
                    double& i1_unreconstructed, double& z_side,
                    bool measure_only) const;
 
    NetworkMeshData*  mesh_  = nullptr;
    NetworkStateData* state_ = nullptr;

    const FvStepForcing* forcing_ = nullptr;
    FvOptions         opts_{};
 
    // Face scratch. corr_l/corr_r are the Audusse well-balanced corrections
    // g·(I₁(h_K) − I₁(h*_K)) — per-cell, not per-face, which is why they are
    // stored separately from the shared flux.
    std::vector<double> f_mass_, f_mom_, f_sstar_, f_corr_l_, f_corr_r_;

    std::vector<double> f_phi_l_, f_phi_r_;

    std::vector<double> f_phi_flux_;

    std::vector<double> lo_flux_, anti_flux_, td_, anew_, rplus_, rminus_;

    std::vector<kernels::FaceState> f_state_l_, f_state_r_;
    std::vector<kernels::FaceFlux>  f_flux_;
    bool hllc_ = true;

    std::vector<double> cell_slope_;

    std::vector<double> cell_eta_slope_, cell_u_slope_;

    std::vector<char> cell_ho2_;
 
 
    // Cell scratch.
    std::vector<double> cell_eta_;     
    std::vector<double> cell_u_;       

    std::vector<double> uf_grad_;

    void computeUfGradients(const std::vector<int>* cells);
 
    // -- TPA pressure closure (issue #156 Phase 4) ---------------------------
    bool tpa_ = false;                    
    std::vector<uint8_t> tpa_scratch_;    
    std::vector<uint8_t> save_tpa_;       

    void updateTpaFlags();

    bool tpaCell(std::size_t uc) const noexcept {
        return tpa_ && state_->cell_tpa[uc] != 0;
    }
    std::vector<double> cell_q_int_;   
 
    // Node scratch.
    std::vector<double> node_exch_;    
    std::vector<double> node_in_, node_out_;
    std::vector<double> flood_vol_;
    std::vector<uint8_t> inlet_control_;
 
    // Step-rejection snapshot. A cell can cross the crown INSIDE a substep,
    // taking its celerity from the free-surface value to the slot value — a
    // factor of ~20 on a 3 ft pipe — so a step sized on the pre-step state can
    // violate CFL by more than an order of magnitude exactly when the model is
    // doing something interesting. Re-censusing the POST-step state and rolling
    // back when it disagrees is the principled fix; a heuristic fill-rate cap
    // is not, because the stiffness ratio depends on FV_SLOT_CELERITY.
    std::vector<double> save_cell_a_, save_cell_q_, save_cell_phi_;
    std::vector<double> save_node_vol_, save_node_head_;
    std::vector<double> save_exch_, save_in_, save_out_, save_flood_, save_qint_,
                        save_carry_;
 
    // SSP-RK2 stage-1 snapshot. Separate from the rejection snapshot above:
    // a rejected RK2 step has to roll back to Uⁿ, and Uⁿ must survive the
    // stage-1 save that the rejection machinery performs.
    std::vector<double> rk_cell_a_, rk_cell_q_, rk_cell_phi_;
    std::vector<double> rk_node_vol_, rk_node_head_;
    std::vector<double> rk_exch_, rk_in_, rk_out_, rk_flood_, rk_qint_,
                        rk_carry_;
    std::vector<uint8_t> rk_tpa_;   

    static constexpr double kStepAcceptRatio = 0.5;
    static constexpr int    kMaxStepRetries  = 8;

    static constexpr double kMinSubstep = 1.0e-6;   
 
    // Divergence guard (issue #156 R3). The retry loop above is sound for CFL
    // violations, which shrink away with dt — but a dt-INDEPENDENT
    // amplification (measured: RK2 x TPA on the study's e3 siphon grows the
    // crown head x20-40 per millisecond at FV_CFL 0.9, 0.45 and 0.25 alike)
    // fails all 8 retries identically and the loop then ACCEPTS the diverged
    // step silently. These bounds are physical absurdities, not tolerances: no
    // water transient carries particle velocity past 1000 ft/s (celerity is
    // not particle speed), and 10,000 ft of pressure head is 300 bar. A run
    // that crosses either has left the physics and must say so instead of
    // scoring OK (P6 finding F3). Checked once per routing step, O(cells).
    static constexpr double kDivergedVelocity = 1000.0;    
    static constexpr double kDivergedDepth    = 10000.0;   
    int         div_link_  = -1;       
    double      div_value_ = 0.0;
    const char* div_what_  = nullptr;
 
    // Work lists (plan §5.2.1). `halo_` is the compaction safety margin: a wet
    // front advances at most CFL cells per substep, so a halo of `rebuild_
    // interval_` cells keeps a stale list conservative between rebuilds.
    std::vector<int>  active_faces_;
    std::vector<char> cell_active_;
    std::vector<char> halo_prev_;   
    bool              lists_valid_  = false;
    int               since_rebuild_ = 0;
    static constexpr int kRebuildInterval = 8;
 
    // Local time stepping (plan §3.3). Tiers are assigned only at a
    // synchronisation point — never mid-cycle — because a cell re-tiered
    // between firings would either skip a flux it owes or drain one twice.
    std::vector<std::uint8_t>     cell_tier_, face_tier_, node_tier_;
    std::vector<std::vector<int>> cells_by_tier_, faces_by_tier_, nodes_by_tier_;

    std::vector<std::vector<int>> due_f_upto_, due_c_upto_, due_n_upto_;

    bool   lts_valid_    = false;
    int    lts_countdown_ = 0;
    double lts_dt0_      = 0.0;
    static constexpr int kRetierEveryCycles = 8;

    std::vector<double> acc_a_, acc_q_, acc_nvol_;

    int  lts_tiers_ = 1;
    std::array<long, kMaxLtsTiers> tier_occupancy_{};
 
    // Statistics.
    long   last_nsteps_  = 0;
    double last_h_       = 0.0;
    double suggested_h_  = 0.0;
    long   total_steps_  = 0;
    long   total_flux_   = 0;
    double min_h_        = 0.0;
    double sim_time_     = 0.0;
 
    // dt-argmin attribution (slot program R0): who owned the binding CFL
    // element, counted once per censusDt (global path) / assignTiers (LTS
    // path). `mutable` because censusDt is const; these are telemetry.
    mutable long dt_argmin_pressurized_ = 0;
    mutable long dt_argmin_band_        = 0;
    mutable long dt_argmin_free_        = 0;
    mutable long dt_argmin_node_        = 0;
    double active_sum_   = 0.0;
    double active_min_   = -1.0;
    double active_max_   = -1.0;
    long   active_n_     = 0;

    double dt_cache_     = 0.0;

    double dt_census_    = 0.0;

    int    census_count_ = 0;
};
 
} // namespace openswmm::fv
 
#endif // OPENSWMM_ENGINE_FV_EXPLICIT_FV_SOLVER_HPP