ExplicitInertialSolver.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_2D_EXPLICIT_INERTIAL_SOLVER_HPP
#define OPENSWMM_ENGINE_2D_EXPLICIT_INERTIAL_SOLVER_HPP
 
#include <array>
#include <cstdint>
#include <string>
#include <vector>
 
#include "ISurfaceSolver.hpp"
#include "InertialEdges.hpp"
#include "../data/SolverOptions2D.hpp"   // Momentum2D
 
namespace openswmm::twoD {
 
class SubsurfaceSolver;   // G1: the two-zone groundwater kernel, if authored
 
class ExplicitInertialSolver final : public ISurfaceSolver {
public:
    void initialize(MeshData& mesh, SurfaceStateData& state,
                    SolverOptions2D& opts) override;
    double advance(double t_current, double t_target) override;
    void reinitialize(double t0) override;
    void resyncFromVolumes(double t0) override;
    void finalize() override;
 
    long   last_num_steps() const noexcept override { return last_steps_; }
    double last_step_size() const noexcept override { return last_dt_; }
    RunStats run_stats() const noexcept override;
    bool is_initialized() const noexcept override { return initialized_; }
    const std::vector<double>& last_coupling_exchange()
        const noexcept override {
        return exch_;
    }

    double spillDeficit() const noexcept { return spill_deficit_; }
 
private:
    // Recompute η/depth from state volumes for the whole mesh.
    void reconstructAll();
    // Flush every pending face accumulator into its cell (CSR gather over ALL
    // cells) and refresh the closure of touched cells. MUST run before any
    // tier/active-set reassignment: a cell deactivated or re-tiered with a
    // pending side-accumulator strands flux whose counterpart the other side
    // already applied — the backstop then realizes it as created volume.
    void settleAccumulators();
    // Apply lazily-accumulated sources (rain/held coupling) to INACTIVE cells
    // over [t_last_sync_, t], rebuild the active cell/edge lists, and assign
    // the LTS tiers (dt0_ = finest active CFL requirement).
    void syncAndRebuild(double t);
    // Tighten-only dt0_ refresh from CURRENT depths/speeds between rebuilds
    // (dt0_ may only grow at syncAndRebuild, which reassigns the tiers).
    void refreshDt0();
    // Fire one tier's faces over their Δt: inertial update + Froude cap +
    // face-cadence positivity share, booking ±ΔM into both side accumulators.
    // global_step = every active cell fires at this same Δt (the window tail),
    // so a face's exporter republishes its volume every time this face does
    // and the positivity share is NOT divided down by a refire ratio.
    void fireFaces(const std::vector<int>& faces, double dt_f,
                   bool global_step = false);
    // Fire one tier's cells over their Δt: gather + clear own-side face
    // accumulators, apply sources, refresh closure + Perot vector; tier-0
    // firings also evaluate the availability-clamped boundary edges.
    // tier0 = this firing carries the tier-0 cadence work (boundary edges and
    // the live junction exchange), which fires once per finest substep.
    void fireCells(const std::vector<int>& cells, double dt_c, bool tier0);
    template <bool kSwe>
    void fireCellsImpl(const std::vector<int>& cells, double dt_c, bool tier0);
    // One halving-order macro cycle of nsub base substeps: tier k fires every
    // 2^k substeps.
    void runMacroCycle(double dt0, int nsub);
 
    MeshData*         mesh_  = nullptr;
    SurfaceStateData* state_ = nullptr;
    SolverOptions2D*  opts_  = nullptr;
 
    InertialEdges edges_;
    std::vector<double>  q_;
    std::vector<double>  qcx_, qcy_;
 
    // -----------------------------------------------------------------------
    // Momentum closure (MOMENTUM_EQUATION, 2D_FULL_SWE plan §2.1). One
    // marcher, three face laws; everything below the face/cell kernels
    // (tiers, active sets, positivity, coupling, transport) is shared.
    // -----------------------------------------------------------------------
    Momentum2D mode_ = Momentum2D::LOCAL_INERTIAL;
    bool       has_quads_ = false;
    std::vector<int>     wall_ptr_, wall_slot_;
    std::vector<double>  macc_x_L_, macc_x_R_, macc_y_L_, macc_y_R_;
    double               swe_bc_dqx_ = 0.0, swe_bc_dqy_ = 0.0;
    bool                 second_order_ = false;
    std::vector<double>  gex_, gey_, gux_, guy_, gvx_, gvy_;
    std::vector<double>  rk_v0_, rk_qx0_, rk_qy0_;
    std::vector<double>  rk_inf0_, rk_bc0_, rk_ex0_, rk_drawn0_, rk_cpl0_;
    void computeLimitedGradientsSwe();
    void runRk2Step(double dt);
    bool                 front_rebuild_ = false;
    std::vector<uint8_t> frontier_;
    bool                 front_breach_  = false;
    double               h_on_front_    = 0.0;   
    std::vector<double>  dw_slope_;
    void fireFacesInertial(const std::vector<int>& faces, double dt_f, bool global_step);
    void fireFacesSwe(const std::vector<int>& faces, double dt_f, bool global_step);
    void fireFacesDiffusive(const std::vector<int>& faces, double dt_f, bool global_step);
    void bookFaceSpecies(int e, int a, int b, double dM, double hf, double dt_f,
                         bool global_step) noexcept;
    double boundaryFluxSwe(std::size_t k, int i, double dt_c);
    void refreshDiffusiveSlopes();
    std::vector<uint8_t> cell_active_;
    std::vector<int>     active_cells_;
    std::vector<int>     active_faces_;
    std::vector<uint8_t> pin_t0_;       
 
    // Tiered LTS. tier_[i] = k means cell i updates every 2^k base substeps
    // with Δt = 2^k·dt0; face tier = min of its incident cells so a face
    // always integrates at the finer cadence. Every face firing books the
    // identical ±ΔM into facc_L_/facc_R_ (single writer per face); each cell
    // applies + clears its own side at its own firing — conservation across
    // tier interfaces is exact by construction.
    std::vector<uint8_t> rebuild_seed_;
    std::vector<double>  rebuild_dt_cell_;
    std::vector<double>  rebuild_dt_partial_;
 
    std::vector<uint8_t> tier_;         
    std::vector<uint8_t> face_tier_;    
    std::vector<double>  facc_L_;       
    std::vector<double>  facc_R_;       
 
    // S1 — species mass rides the SAME face accumulators, one pair per
    // species, [s * ne + e]. Booked in fireFaces immediately after the volume
    // ΔM, from the FINAL qn1 (after the Froude cap and the positivity share),
    // at the exporting cell's concentration read at that same substep;
    // gathered and cleared in fireCells alongside the volume side. That is
    // D-2DT2: the species flux inherits the volume flux's tier cadence rather
    // than reproducing it, so conservation across tier interfaces is the
    // marcher's own property and not a second one to prove. Empty unless
    // `state_->transport.active()`.
    std::vector<double>  sacc_L_;
    std::vector<double>  sacc_R_;
    double donorConc(int s, int cell) const noexcept;
    void sinkMassAtCellConc(int i, double dv_m3, std::vector<double>& ledger,
                            double* per_point_ledger = nullptr) noexcept;
    void addRainMass(int i, double rain_m3) noexcept;
    void addCouplingSourceMass(int i, double area_dt) noexcept;
    void sinkTemperatureWithEvap(int i, double evap_m3) noexcept;
    bool accumulators_pending_ = false;
    std::vector<std::vector<int>> cells_by_tier_;
    std::vector<std::vector<int>> edges_by_tier_;
    double dt0_ = 0.0;                  
    std::vector<int>     bc_cell_;      
    std::vector<int>     bc_slot_;      
    std::vector<double>  bc_accum_;     
    std::vector<double>  bc_q_;         
 
    // Live junction exchange (windowless coupling): state_->node_coupling
    // points, evaluated at tier-0 cadence against live 2D heads and the
    // routing-step 1D heads. exch_[k] = ∫Q_k dt (m³, + = 2D→1D), reset per
    // advance; node_drawn_ caps a step's total spill at the node's stored
    // volume so fill-and-spill thrash is structurally impossible.
    std::vector<double>  exch_;
    std::vector<double>  node_drawn_;   
    double               spill_deficit_ = 0.0;
 
public:
    void setExchangeHeadSlopes(std::vector<double> slopes) {
        exch_head_slope_ = std::move(slopes);
    }
private:
    std::vector<double>  exch_head_slope_;
    double               exch_tau_ = 0.0;  
 
    double t_last_sync_ = 0.0;          
    int    cycles_since_rebuild_ = 1000; 
    void lazySourcesOnly(double t);
    long   substeps_run_ = 0;           
    long   face_passes_  = 0;           
    long   last_steps_   = 0;           
    double last_dt_      = 0.0;
    bool   initialized_  = false;
 
    // Active-fraction telemetry: (sim time, active cells) sampled at every
    // rebuild; dumped as CSV at finalize() when OPENSWMM_2D_MARCHER_TELEMETRY
    // names a file. The Phase-1 gate reads this to verify the thin-film budget
    // assumption on the Bellinge storm slice.
    std::vector<std::pair<double, int>> telemetry_;
    std::string telemetry_path_;
    std::array<long, 8> tier_occupancy_{};

    SubsurfaceSolver* gw_ = nullptr;
 
public:
    void setSubsurface(SubsurfaceSolver* gw) noexcept { gw_ = gw; }

    const InertialEdges& inertialEdges() const noexcept { return edges_; }
};
 
} // namespace openswmm::twoD
 
#endif // OPENSWMM_ENGINE_2D_EXPLICIT_INERTIAL_SOLVER_HPP