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 "SubsurfaceData.hpp"
 
#include <string>
#include <vector>
 
namespace openswmm {
struct NodeData;
}  // namespace openswmm
 
namespace openswmm::twoD {
 
struct MeshData;
struct InertialEdges;
struct SurfaceStateData;
struct SolverOptions2D;
 
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);
 
    bool active() const noexcept { return state_.active; }
    SubsurfaceState&       state()       noexcept { return state_; }
    const SubsurfaceState& state() const noexcept { return state_; }
    const GwOptions&       options() const noexcept { return options_; }
 
    // ---- LTS hooks -------------------------------------------------------

    void refreshDtCell(const MeshData& mesh, const InertialEdges& edges);

    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);
    void settle(SurfaceStateData& surf);

    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;
    double takeToSurface(int cell) noexcept;
    const std::vector<int>& pendingSurfaceCells() const noexcept {
        return pending_surface_;
    }

    void compactPending() noexcept;

    void sampleNodeExchange(const NodeData* nodes, double dt);

    const std::vector<double>& nodeExchangeVolumes() const noexcept {
        return node_exchange_vol_;
    }

    void resetNodeExchangeVolumes() noexcept;

    double storage() const noexcept { return state_.storage(); }
 
private:
    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 gatherLateral(int i) noexcept;
    double gatherNode(int i) noexcept;
    void markPendingSurface(int i) noexcept;
 
    SubsurfaceState  state_;
    GwOptions        options_;
    std::vector<GwNodeBed> node_beds_;
 
    const MeshData*      mesh_  = nullptr;
    const InertialEdges* edges_ = nullptr;
    const SolverOptions2D* opts_ = nullptr;
 
    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_;  
    bool accumulators_pending_ = false;

    std::vector<int> subcatch_outlet_node_;
    std::vector<int> cell_subcatch_;
};
 
}  // namespace openswmm::twoD
 
#endif  // OPENSWMM_ENGINE_2D_SUBSURFACE_SOLVER_HPP