SubsurfaceSolver.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_SUBSURFACE_SOLVER_HPP
#define OPENSWMM_ENGINE_2D_SUBSURFACE_SOLVER_HPP
 
#include "SigmaColumn.hpp"
#include "../solver/GwRechargeBackend.hpp"
#include "../../transport/components/ReactionModule/ReactionIntegrator.hpp"
#include "SurfaceExchange.hpp"
#include "../gw/GwSourceForcing.hpp"
#include "../gw/GwBoundaryForcing.hpp"
#include "SubsurfaceData.hpp"
#include "SubsurfaceTransportState.hpp"
#include "../coupling/RuntimeCoupling.hpp" // T7.1
 
#include <algorithm>
#include "../heat/GroundHeat.hpp"
#include <string>
#include <vector>
 
namespace openswmm::detail {
class CompletedSourceArchive;
}
 
namespace openswmm {
struct NodeData;
struct SimulationContext;
} // namespace openswmm
 
namespace openswmm::twoD {
struct RichardsBedStepper;
 
struct MeshData;
struct InertialEdges;
struct SurfaceStateData;
struct SolverOptions2D;
struct RichardsBedOwner;
struct RichardsBedAcceptance;
struct RichardsBedReceipt;
struct GwTransportData;
enum class GwScope : int8_t;

struct RowLayoutLite {
    int n_species = 0;
    int n_pollut = 0;
    int n_msx = 0;
    int age_row = -1;
    int temp_row = -1;
    std::vector<std::string> names;
};
 
using openswmm::NodeData;

class SubsurfaceSolver {
public:
    std::string initialize(const MeshData& mesh, const InertialEdges& edges,
                           const SolverOptions2D& opts, const GwUnitFactors& uf, int n_nodes,
                           SubsurfaceConfig& cfg, std::vector<std::string>& warnings);

    void bookSurfaceEt(int cell, double dt, double potential_rate, double actual_surface_volume);
    void bookAreaEt(int cell, double dt, double potential, double owner_evaporation,
                    double soil_demand);

    bool active() const noexcept { return state_.active; }
    void setRechargeBackend(GwRechargeBackend backend) { recharge_backend_ = std::move(backend); }
    const GwRechargeBackend& rechargeBackend() const noexcept { return recharge_backend_; }
    uint64_t rechargeBatches() const noexcept { return recharge_batches_; }
    uint64_t rechargeCells() const noexcept { return recharge_cells_; }
    int rechargeThreads() const noexcept {
        return recharge_backend_.threads ? recharge_backend_.threads() : 0;
    }

    const std::vector<uint8_t>& coverage() const noexcept { return covered_; }
    bool covers(int cell) const noexcept {
        return cell >= 0 && static_cast<std::size_t>(cell) < covered_.size() &&
               covered_[static_cast<std::size_t>(cell)] != 0;
    }

    SubsurfaceState& state() noexcept { return state_; }
    const SubsurfaceState& state() const noexcept { return state_; }
 
    // ---- T7.1: species / age / enthalpy -----------------------------------

    SubsurfaceTransportState& transport() noexcept { return tr_; }
    const SubsurfaceTransportState& transport() const noexcept { return tr_; }
    GroundHeatState& heat() noexcept { return heat_; }
    const GroundHeatState& heat() const noexcept { return heat_; }
    double pendingHeat(int cell) const noexcept;
    std::string initHeat(SimulationContext&, const GwTransportData&);
    std::string initPorewater();
    std::string refreshSigmaPorewater(int cell);
    void refreshSigmaPorewaterReservations(int cell) noexcept;
    std::string refreshPorewaterContactLedger();
    void initTransport(const RowLayoutLite& rows, const GwTransportData* gw,
                       const std::vector<double>& pollut_decay, std::vector<std::string>& warnings);
    double satVolume(int cell) const noexcept {
        const auto u = static_cast<std::size_t>(cell);
        return state_.hg[u] * state_.theta_s[u] * state_.area[u];
    }

    double bookLinkGain(int cell, int species, double vol) noexcept {
        if (!tr_.active() || !(vol > 0.0) || cell < 0 || cell >= tr_.n_cells || species < 0 ||
            species >= tr_.n_species)
            return 0.0;
        const double v = satVolume(cell);
        if (!(v > 0.0))
            return 0.0;
        const auto u = static_cast<std::size_t>(cell);
        const auto k = tr_.idx(species, cell);
        const double mob = tr_.dissolvedFraction(species, cell, state_.theta_s[u]);
        double m = vol * (tr_.sat_mass[k] / v) * mob;
        if (heat_.active() && species == tr_.temp_row) {
            heat_link_out_volume_[u] += vol;
            bookPresetHeatOut(cell, vol, m);
        }
        if (!tr_.signedRow(species)) {
            const double avail =
                std::max(0.0, tr_.sat_mass[k] * mob - pendingSaturatedMassExport(cell, species));
            m = std::min(m, avail);
        }
        if (m == 0.0)
            return 0.0;
        bookPresetPorewaterOut(cell, species, m);
        tr_.lout_mass[k] += m;
        if (tr_.n_layers > 0 && heat_.active())
            refreshSigmaPorewaterReservations(cell);
        return m;
    }

    double unsatVolume(int cell) const noexcept {
        const auto u = static_cast<std::size_t>(cell);
        return state_.hu[u] * state_.area[u];
    }

    void bookInfiltrationMass(int cell, int species, double mass) noexcept;
    bool sourceReceiptOwnershipBound() const {
        return std::any_of(source_owned_.begin(), source_owned_.end(), [](bool x) { return x; });
    }

    std::string configureSourceReceipts(const std::vector<int>& cells);
    void retainSourceReceipt(const SurfaceIntakeReceipt&);
    const std::vector<SurfaceIntakeReceipt>& pendingSourceReceipts(int cell) const {
        return source_receipts_.at(cell);
    }

    double sourceCompletedTime(int cell) const { return source_cell_end_.at(cell); }
    const std::vector<SurfaceReturnReceipt>& sourceReturns() const { return source_returns_; }
    void acknowledgeSourceReturns(double completed_through);
    double takeToSurfaceMass(int cell, int species) noexcept;
    const GwOptions& options() const noexcept { return options_; }
 
    // ---- LTS hooks -------------------------------------------------------

    void refreshDtCell(const MeshData& mesh, const InertialEdges& edges);
    double minDtCell() const noexcept {
        double m = 1.0e30;
        for (const double v : state_.dt_cell)
            if (v < m)
                m = v;
        return m;
    }

    void assignTiers(double dt0, int n_tiers);

    int requiredTiers(double dt0) const noexcept;

    void fireGwFaces(int tier, double dt);
    void fireGwCells(int tier, double dt, SurfaceStateData& surf, double time = 0.0,
                     double accepted_end = -1.0);
    RuntimeSources& runtimeSources() noexcept { return runtime_sources_; }
    const RuntimeSources& runtimeSources() const noexcept { return runtime_sources_; }
    void setBoundaries(std::vector<GwResolvedBoundary>);
    const std::vector<GwResolvedBoundary>& boundaries() const noexcept { return boundaries_; }
    // One forward-Euler evaluation within the owned shared RK2 transaction.
    // Prescribed forcing is integrated without additional water substeps.
    void fireGwRk2Cells(double dt, SurfaceStateData&, double time);
    void advanceRk2Heat(double seconds, double sample_time);
    void setSources(std::vector<GwResolvedSource> sources);
    void settle(SurfaceStateData& surf);
    void restorePendingSurface() {
        pending_surface_.clear();
        std::fill(pending_flag_.begin(), pending_flag_.end(), 0);
        for (int i = 0; i < state_.n_cells; ++i)
            if (state_.xacc_to_surface[i] > 0.0)
                markPendingSurface(i);
    }

    bool tierHasCells(int tier) const noexcept {
        return tier >= 0 && tier < static_cast<int>(cells_by_tier_.size()) &&
               !cells_by_tier_[static_cast<std::size_t>(tier)].empty();
    }

    bool tierHasFaces(int tier) const noexcept {
        return tier >= 0 && tier < static_cast<int>(faces_by_tier_.size()) &&
               !faces_by_tier_[static_cast<std::size_t>(tier)].empty();
    }

    int tierCount() const noexcept { return static_cast<int>(cells_by_tier_.size()); }
    const std::vector<long>& tierFirings() const noexcept { return tier_firings_; }
 
    // ---- surface ↔ subsurface --------------------------------------------

    void bookInfiltrationFromSurface(int cell, double vol_m3) noexcept;
    void publishInfiltration(SurfaceStateData& surf, double interval, double time);
    double acceptSurfaceInfiltration(int cell, double requested, bool prescribed = false) noexcept;
    double infiltrationHeadroom(int cell) const noexcept;
    double sourceInfiltrationCapacity(int cell, double pond_m, double interval_s) const noexcept;
    void bookLinkSeepage(int cell, double vol_m3) noexcept;
    void bookLinkContact(int cell, double volume, const std::vector<double>& concentration,
                         double elevation_begin, double elevation_end);
    int contactLayer(int cell, double elevation) const noexcept;
    double takeToSurface(int cell) noexcept;
    const std::vector<int>& pendingSurfaceCells() const noexcept { return pending_surface_; }
    void compactPending() noexcept;

    void publishNodeConcentrations(const SimulationContext& ctx);
    void setNodeRowConc(const double* rows, int n_species) noexcept {
        node_row_conc_ = rows;
        node_row_ns_ = n_species;
    }

    const std::vector<double>& nodeExchangeMass() const noexcept { return tr_.node_out_mass; }
    void bookLinkSeepageMass(int cell, int species, double mass) noexcept;

    void sampleNodeExchange(const NodeData* nodes, double dt, double end = -1);

    const std::vector<double>& nodeExchangeVolumes() const noexcept { return node_exchange_vol_; }
    void resetNodeExchangeVolumes() noexcept;
    const std::vector<GwNodeBed>& nodeBeds() const noexcept { return node_beds_; }
    const std::vector<double>& bedExchangeCumulative() const noexcept { return bed_exchange_cum_; }
 
    // ---- G-W1b: what a gate needs to prove its own deck is an instrument --
    //
    // A gate that measures a cap must first establish that the cap BOUND on
    // the deck it is measuring, or it can pass with the defect present. The
    // ratio between the two specific yields does not establish that — a deck
    // scoring 56.8× was measured bit-identical with the bug — because what
    // makes the cap bind is `h_g · Sy · A` being small in absolute terms.
    // Only the clamp knows. These three say so directly, with no threshold
    // to tune and no geometry that can make them lie.

    long drainCapBinds() const noexcept { return drain_cap_binds_; }
    long nodeRefunds() const noexcept { return node_refunds_; }
    double specificYieldOf(int cell) const noexcept;
    double unreservedDrainableWater(int cell) const noexcept;
    double unreservedDrainableWater(int cell, double specific_yield) const noexcept;
    double pendingSaturatedMassExport(int cell, int species) const noexcept;
    // Internal R5 trial budget in SI, before accepting a new node receipt.
    // Separate ceilings; incoming pending water never finances a withdrawal.
    // Call after atmospheric demand/other accepted channels have been booked.
    struct NodeWaterBudget {
        double head = 0, donor = 0, receiver = 0;
    };

    std::string nodeWaterBudget(int cell, double interval, NodeWaterBudget&) const;
    double cellCompletedTime(int cell) const { return cell_completed_.at(cell); }
    const std::vector<double>& richardsBedDownCumulative() const noexcept {
        return richards_down_cum_;
    }

    const std::vector<double>& richardsBedUpCumulative() const noexcept { return richards_up_cum_; }
    // Hydraulic columns book completed intervals; groundwater consumes them
    // only at its own firing. Gross reservations cannot cancel each other.
    struct PendingLidBed {
        double start = 0, end = 0, down = 0, up = 0, down_j = 0, up_j = 0;
        bool heat_applied = false;
    };

    std::string configureLidReceipts(const std::vector<int>& beds,
                                     const std::vector<int>& mobile = {});
    std::string retainLidReceipts(double start, double end, const std::vector<RichardsBedReceipt>&);
    std::string completeLidThermalInterval(SimulationContext&);
    bool lidReceiptsEnabled() const noexcept { return lid_receipts_enabled_; }
    double pendingLidDown(int cell) const noexcept;
    double pendingLidUp(int cell) const noexcept;
    double pendingLidWater() const noexcept;
    const std::vector<std::vector<PendingLidBed>>& pendingLidReceipts() const {
        return lid_receipts_;
    }
    struct LidBedDiagnostic {
        double start = 0, end = 0, groundwater_time = 0, lid_head = 0, groundwater_head = 0;
        double requested_down = 0, requested_up = 0, reserved_down = 0, reserved_up = 0,
               accepted_down = 0, accepted_up = 0;
    };

    const std::vector<LidBedDiagnostic>& lidBedDiagnostics() const { return lid_bed_diagnostics_; }
    double pendingBedDown(int bed) const noexcept;
    double pendingBedUp(int bed) const noexcept;
    void setLidBedClogging(int bed, double factor) { lid_bed_factor_.at(bed) = factor; }

    double storage() const noexcept { return state_.storage(); }
    void checkpoint(detail::CompletedSourceArchive&, bool model); // internal accepted-boundary IO
    void checkpoint(detail::CompletedSourceArchive&, bool model, const double* held_air);
    std::string initReactions(SimulationContext& ctx);
    std::string rk2ProfileError() const;
    void blendRk2Water(const SubsurfaceSolver& initial, double end);
    void publishRk2Infiltration(SurfaceStateData&, double dt, double time);
    bool hasRk2Reactions() const noexcept;
    void advanceRk2Reactions(double seconds);
    bool needsTransportRestart() const noexcept {
        return tr_.n_layers > 0 || !tr_.generated_age.empty() || !boundaries_.empty() ||
               !reactors_.empty() ||
               (tr_.active() && std::any_of(source_owned_.begin(), source_owned_.end(),
                                            [](bool owned) { return owned; }));
    }

    std::string transportRestartState() const;
    std::string restoreTransportRestart(const std::string& bytes, bool apply);
 
private:
    double infiltrationStorageDeficit(int cell) const noexcept;
    double infiltrationHeadroom(int cell, double deficit) const noexcept;
    friend class CommonSourceStage;
    friend class ExplicitInertialSolver;
    friend RichardsBedAcceptance acceptRichardsBedsAtFiring(SubsurfaceSolver&, SurfaceStateData&,
                                                            double, double,
                                                            const std::vector<RichardsBedOwner>&);
    friend RichardsBedAcceptance acceptRichardsBedsForTier(SubsurfaceSolver&, SurfaceStateData&,
                                                           int, double, double,
                                                           const std::vector<RichardsBedOwner>&);
    friend RichardsBedAcceptance acceptRichardsBedsAcrossSteps(SubsurfaceSolver&, SurfaceStateData&,
                                                               double, double,
                                                               const std::vector<double>&,
                                                               const std::vector<RichardsBedOwner>&,
                                                               const RichardsBedStepper*);
    friend RichardsBedAcceptance acceptRichardsBedsForTierAcrossSteps(
        SubsurfaceSolver&, SurfaceStateData&, int, double, double, const std::vector<double>&,
        const std::vector<RichardsBedOwner>&, const RichardsBedStepper*);
    void resolveRows(const MeshData& mesh, const GwUnitFactors& uf, SubsurfaceConfig& cfg,
                     std::vector<std::string>& warnings);
    void resolveClosures(std::vector<std::string>& warnings);
    soil::Params paramsOf(int i) const noexcept;
    void fireCell(int i, double dt, SurfaceStateData& surf, double time);
    void fireGwCellList(const std::vector<int>&, double dt, SurfaceStateData&, double time,
                        double accepted_end);
    std::vector<bool> source_owned_;
    std::vector<std::string> source_species_;
    std::vector<double> source_cell_end_;
    std::vector<std::vector<SurfaceIntakeReceipt>> source_receipts_;
    std::vector<SurfaceReturnReceipt> source_returns_;
    double source_firing_end_ = 0.0;
    double source_completed_end_ =
        0.0; // transient R5 readiness; old restores cannot activate Richards ownership
    double gather_fraction_ = 1.0;
    double consumePending(double& value) noexcept {
        const double take = value * gather_fraction_;
        value -= take;
        return take;
    }

    struct CellFlux {
        double node_in = 0;  
        double link_out = 0; 
        const sigma::ColumnTransfers* column = nullptr;
        double boundary_in = 0, boundary_out = 0;
        const std::vector<GwBoundaryPulse>* boundaries = nullptr;
        double hg_final = 0;
        double runtime_dt = 0.0, requested_source_out = 0.0;
        double v_sat0 = 0.0, v_uns0 = 0.0; 
        double lateral = 0.0;              
        double node_out = 0.0;
        double infil_in = 0.0;
        double link = 0.0;
        double source_in = 0.0, source_out = 0.0;
        const std::vector<double>* source_mass = nullptr;
        double recharge = 0.0;
        double handover = 0.0;
        double grid_volume =
            0.0;               
        double lid_in_j = 0.0; 
        double et = 0.0;       
        double deep = 0.0;
        double dt = 0.0;    
        double dunne = 0.0; 
        double paired_refund = 0.0;
        const std::vector<double>* refund_mass = nullptr;
        std::vector<double>* dunne_mass = nullptr;
        std::vector<double>* excess_mass = nullptr;
        double reject = 0.0; 
    };
    std::vector<std::string> transport_units_;
    struct GroundReactor {
        std::vector<ReactionExprForm> forms;
        std::vector<RxExprSpan> expressions;
    };
    const ReactionData* reactions_ = nullptr;
    transport::RxWorkspace reaction_workspace_;
    std::vector<GroundReactor> reactors_;
    std::string reaction_identity_;
    void reactCell(int cell, double dt, double hg);

    void fireCellSpecies(int i, const CellFlux& f);
    void fireLayerSpecies(int i, const CellFlux& f);
    void reactPorewaterCell(int cell, double dt, double hg);
    void fireGroundHeatSpecies(int i, const CellFlux& f) noexcept;
    std::string groundWaterIncrement(int cell, const CellFlux&, std::vector<double>&) const;
    void fireGroundPorewaterSpecies(int i, const CellFlux& f) noexcept;
    double gatherLateral(int i) noexcept;
    double gatherLateralMass(int i, int s) noexcept;
    bool scopeCoversCell(int i, GwScope scope, const std::string& tag, int cell) const noexcept;
    double gatherNodeMass(int i, int s) noexcept;
    double pollutantDecay(int s, const std::vector<double>& pollut_decay) const noexcept;
    double gatherNode(int i) noexcept;

    double specificYield(std::size_t u, const soil::Params& p, double ts, double L0, GwClosure cl,
                         double* theta_bot_out) const noexcept;
    void markPendingSurface(int i) noexcept;
 
    void publishCellInfiltration(int, SurfaceStateData&, double interval, double time);
    std::vector<std::vector<std::size_t>> cell_beds_;
    RuntimeSources runtime_sources_;
    std::vector<GwResolvedBoundary> boundaries_;
    std::vector<std::vector<std::size_t>> cell_boundaries_;
    std::vector<GwResolvedSource> sources_;
    std::vector<std::vector<std::pair<std::size_t, double>>> cell_sources_;
    SubsurfaceState state_;
    SubsurfaceTransportState tr_; 
    GroundHeatState heat_;
    std::vector<double> heat_link_out_volume_;
    std::vector<double>
        heat_preset_out_j_; // [cell*thermal_layers+layer], accepted producer-time exports
    std::vector<double>
        heat_preset_out_water_; // paired physical layer volumes, including source-interval splits
    std::vector<double> heat_node_out_volume_,
        heat_node_gathered_out_;               // gross accepted exports per hydraulic bed
    std::vector<double> porewater_preset_out_; // [species*n_cells*layers + cell*layers + layer]
    void bookPresetHeatOut(int cell, double volume, double temperature_volume) noexcept;
    std::string conductGroundWithLids(int cell, double seconds);
    void bookPresetPorewaterOut(int cell, int species, double mass) noexcept;
    double saturatedPorewaterMobileMass(int cell, int species) const noexcept;
    struct ThermalFaceSlice {
        int left = 0, right = 0;
        double low = 0, high = 0;
    };
    struct ThermalFace {
        double distance_left = 0, distance_right = 0;
        std::vector<ThermalFaceSlice> slices;
    };
    std::vector<ThermalFace> heat_faces_;
    std::vector<double> heat_pending_left_, heat_pending_right_, heat_dt_;
    std::string initLateralHeat();
    void fireLateralHeatFace(int edge, double seconds);
    void gatherLateralHeat(int cell) noexcept;
    GwOptions options_;
    std::vector<GwNodeBed> node_beds_;
    std::vector<double> bed_exchange_cum_; 
    // Internal R5 owner clocks and directional receipts; execution/checkpoint
    // activation remains guarded until the routing/restart adapter is qualified.
    std::vector<double> cell_completed_, richards_down_cum_, richards_up_cum_;
    std::vector<bool> richards_bed_owned_;
    bool lid_receipts_enabled_ = false;
    std::vector<std::vector<PendingLidBed>> lid_receipts_;
    std::vector<double> lid_booked_end_;
    std::vector<bool> mobile_lid_bed_;
    std::vector<double> lid_bed_factor_;
    std::vector<double> node_pending_down_, node_pending_up_;
    std::vector<LidBedDiagnostic> lid_bed_diagnostics_;
    bool richards_clock_valid_ = true; // old checkpoints do not restore the joint owner clocks
    std::vector<double> bed_gathered_;
    std::vector<double> node_drawn_gw_;
    std::vector<double> node_concentrations_;
    const double* node_row_conc_ = nullptr;
    int node_row_ns_ = 0;
 
    const MeshData* mesh_ = nullptr;
    const InertialEdges* edges_ = nullptr;
    const SolverOptions2D* opts_ = nullptr;
 
    GwRechargeBackend recharge_backend_;
    std::vector<OpenSwmmGwRechargeInputV1> recharge_inputs_;
    std::vector<double> recharge_values_;
    std::vector<int> recharge_slots_;
    uint64_t recharge_batches_ = 0, recharge_cells_ = 0;
    void prepareRechargeBatch(const std::vector<int>&);
    double rawRecharge(int cell, const soil::Params&, double length) const noexcept;
 
    std::vector<std::vector<int>> cells_by_tier_;
    std::vector<std::vector<int>> faces_by_tier_;
    std::vector<uint8_t> face_tier_;
    std::vector<long> tier_firings_;
    std::vector<int> pending_surface_;
    std::vector<uint8_t> pending_flag_;
    std::vector<int> cell_nface_;           
    std::vector<double> node_exchange_vol_; 
    long drain_cap_binds_ = 0;              
    long node_refunds_ = 0;                 
    bool accumulators_pending_ = false;

    std::vector<int> subcatch_outlet_node_;
    std::vector<int> cell_subcatch_;
    std::vector<uint8_t> covered_;
    int bed_gathered_cell_ = -1;
    bool et_boundary_ = false, et_caprise_ = false;
    std::vector<int> reaction_rows_;
 
    // Nonphysical, pointer-free work storage. Snapshots do not copy its capacity;
    // assignment keeps the destination's workspace. Every used slot is reset or
    // overwritten by its next firing; no checkpoint record includes these values.
    struct CellScratch {
        std::vector<double> receipt_mass;
        std::vector<double> refund_mass, dunne_mass, excess_mass;
        CellScratch() = default;
        CellScratch(const CellScratch&) noexcept {}
        CellScratch& operator=(const CellScratch&) noexcept { return *this; }
        CellScratch(CellScratch&&) noexcept = default;
        CellScratch& operator=(CellScratch&&) noexcept = default;
    } cell_scratch_;
};
 
} // namespace openswmm::twoD
 
#endif // OPENSWMM_ENGINE_2D_SUBSURFACE_SOLVER_HPP