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 "AcceptedWaterVolumes.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 enablePressureStorageCandidate(const NonlinearPressureStorage::Controls& controls = {}) {
        pressure_candidate_ = true;
        pressure_controls_ = controls;
    }

    const NonlinearPressureStorage::Result& pressureCandidateResult() const noexcept {
        return pressure_result_;
    }

    long pressureCandidateRejections() const noexcept { return pressure_rejections_; }

    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;
    }

    using ConduitStepSink = void (*)(void* user, int conduit, double dt, double dur);
    void setConduitStepSink(ConduitStepSink fn, void* user) noexcept {
        step_sink_ = fn;
        step_sink_user_ = user;
    }

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

    void trackAcceptedWater(bool enabled) {
        water_tracking_ = enabled;
        if (enabled && mesh_)
            water_.resize(mesh_->n_faces(), mesh_->struct_link.size(), mesh_->n_cells(), mesh_->n_nodes());
        else water_ = {};
        save_water_ = water_; rk_water_ = water_;
    }
    const AcceptedWaterVolumes& acceptedWater() const noexcept { return water_; }
    void requireConservativeNodeBalance(bool enabled) noexcept {conservative_node_balance_=enabled;}
    std::vector<double>& withdrawalRates() noexcept { return withdrawal_rates_; }
    const std::vector<double>& withdrawalRates() const noexcept { return withdrawal_rates_; }
    const std::vector<double>& unmetWithdrawalVolumes() const noexcept { return unmet_withdrawal_; }
 
 

    const std::vector<double>& node_continuity_carry() const noexcept { return node_carry_; }

    const std::vector<std::uint8_t>& pressure_folded_nodes() const noexcept {
        return press_.foldedNodes();
    }

    const std::vector<double>& node_inflow_volume() const noexcept { return node_in_; }
    const std::vector<double>& node_outflow_volume() const noexcept { return node_out_; }
    void bookDivertedLateral(std::size_t un, double dt, const FvStepForcing& forcing) noexcept {
        if (node_lat_div_.empty() || !node_lat_div_[un] || !forcing.node_lateral)
            return;
        const double lat = forcing.node_lateral[un];
        if (lat > 0.0)
            node_out_[un] += lat * dt;
    }

    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:
    double advanceImpl(double t_current, double t_target, const FvStepForcing& forcing);
    void validatePressureCandidate(const FvStepForcing& forcing) const;
    // -- substep pipeline ---------------------------------------------------
    void refreshDepths();
    void refreshNodeAreas();
    void rebuildActiveLists();
    double censusDt(bool press_edit = false) const;
    void reconstructState();
    void computeFaceFlux(int face, bool mass_only = false);
    void reconstructScalars(double dt);
    transport::fvkernels::SpeciesKernelView speciesKernelView();
    void computeFluxes();
    void limitPositivity(double dt);

    void relaxNodeFluxes(double dt, const FvStepForcing& forcing);
    template <class DtOf>
    void solveNodes(const std::vector<int>& nodes, DtOf&& dt_of, const FvStepForcing& forcing);
    std::vector<int> node_solve_list_; 
    std::vector<int>
        node_heavy_; 
    std::vector<int> node_light_; 

    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, bool book_receipts = true);

    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 roughness, double rough_factor, double q,
                       double u, double h, double dt, double radius) 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, int k_cap_override = 0);

    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();

    bool 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;
    bool pressure_candidate_ = false;
    NonlinearPressureStorage::Controls pressure_controls_;
    NonlinearPressureStorage::Result pressure_result_;
    long pressure_rejections_ = 0;

    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> cell_ah_; 
    std::vector<double> cell_t_;  
    std::vector<double> cell_i1_; 
    void cacheClosure(std::size_t uc, const FvGeometry& g, double h,
                      double* radius = nullptr) noexcept;

    std::vector<double> uf_grad_;

    void computeUfGradients(const std::vector<int>* cells);
 
    // -- TPA pressure closure (issue #156 Phase 4) ---------------------------
    bool tpa_ = false;       
    bool mixed_wave_ = true; 
    bool mw_rh_ = true;      
    bool deg1_guard_slot_ =
        true; 
    bool census_pass_far_ =
        true; 
    bool lts_fit_ =
        false; 
    int passThroughFarCell(int face, int nd) const noexcept;
    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_; 
 
    bool water_tracking_ = false;
    bool conservative_node_balance_ = false;
    std::vector<double> withdrawal_rates_, unmet_withdrawal_, save_unmet_, rk_unmet_;
    void bookFaceWater(int face, double dt);
    void bookSpliceReceipt(std::size_t vj, double flux, double dt);
    double capWithdrawal(std::size_t node, double balance, double dt, const FvStepForcing& forcing);
 
    AcceptedWaterVolumes water_, save_water_, rk_water_;
    void bookStructureWater(double dt, const FvStepForcing& forcing);
    double algebraicAvailableWater(std::size_t node, double dt, double pending,
                                   const std::vector<double>& cell_scales) const;
 
    // 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_dummy_, rk_dummy_;
    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_;
 
    ConduitStepSink step_sink_ = nullptr; 
    void* step_sink_user_ = nullptr;
    void emitConduitSteps(double dt0, double dur, bool tiered) const;
    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;
    long n_macro_cycles_ = 0;
    long n_macro_rejected_ = 0;

    bool state_clean_ = false;
    bool acc_dirty_ = false;
    std::vector<double> relax_beta_;
 
    // 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