PressurizedHeadSolver.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_PRESSURIZED_HEAD_SOLVER_HPP
#define OPENSWMM_ENGINE_FV_PRESSURIZED_HEAD_SOLVER_HPP
 
#include <cstdint>
#include <vector>
 
#include "FvKernels.hpp"
#include "INetworkSolver.hpp"
#include "NetworkMeshData.hpp"
#include "NonlinearPressureStorage.hpp"
 
namespace openswmm::fv {

struct PressurizedView {
    const NetworkMeshData* mesh = nullptr;
    NetworkStateData* state = nullptr;
    const FvOptions* opts = nullptr;
    const FvStepForcing* forcing = nullptr;
 
    // Face scratch (owned by ExplicitFvSolver).
    double* f_mass = nullptr;
    double* f_sstar = nullptr;
    kernels::FaceState* f_state_l = nullptr;
    kernels::FaceState* f_state_r = nullptr;
    kernels::FaceFlux* f_flux = nullptr;
 
    // Cell scratch.
    const double* cell_eta = nullptr;
    double* cell_u = nullptr;

    const std::vector<double>* uf_grad = nullptr;
 
    const std::vector<int>* active_faces = nullptr;
    const std::vector<char>* cell_active = nullptr;
 
    // Node machinery (owned by ExplicitFvSolver).
    const std::vector<double>* node_qstruct = nullptr;
    std::vector<double>* node_carry = nullptr;
    const std::vector<std::uint8_t>* node_alg = nullptr;
    const std::vector<double>* node_vfull = nullptr;
    const std::vector<std::uint8_t>* node_lat_div = nullptr;
    const std::vector<double>* cell_qlat = nullptr;
};
 
class PressurizedHeadSolver {
public:
    enum : std::uint8_t { kNone = 0, kFull = 1, kDelta = 2 };

    bool classify(const PressurizedView& v);

    void solve(const PressurizedView& v, double dt);

    NonlinearPressureStorage::Result
    solveCandidate(const PressurizedView& v, double dt,
                   const NonlinearPressureStorage::Controls& controls);
    bool verifyCandidateEndpoint(const PressurizedView& v, double dt,
                                 const NonlinearPressureStorage::Controls& controls);

    void finalizeCells(const PressurizedView& v);

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

    bool anyImplicit() const noexcept { return !imp_faces_.empty(); }
 
    // -- static predicates shared with the census edit ----------------------

    static bool cellPressurized(const NetworkMeshData& mesh, const NetworkStateData& state,
                                int cell) {
        const auto uc = static_cast<std::size_t>(cell);
        if (!state.cell_tpa.empty())
            return state.cell_tpa[uc] != 0;
        const FvGeometry& g = mesh.geom[static_cast<std::size_t>(mesh.cell_geom[uc])];
        return !g.is_open && state.cell_h[uc] >= g.y_crown;
    }

    static std::uint8_t faceModeOf(const NetworkMeshData& mesh, const NetworkStateData& state,
                                   int f) {
        const auto uf = static_cast<std::size_t>(f);
        if (mesh.face_gate[uf] != 0)
            return kNone;
        if (mesh.face_culvert[uf] >= 0)
            return kNone;
        const int cl = mesh.face_cl[uf];
        const int cr = mesh.face_cr[uf];
        const int nd = mesh.face_node[uf];
        if (nd < 0 && (cl < 0 || cr < 0))
            return kNone;
 
        bool pl = false, pr = false, cell_side_press = false;
        if (cl >= 0) {
            pl = cellPressurized(mesh, state, cl);
            cell_side_press = cell_side_press || pl;
        }
        if (cr >= 0) {
            pr = cellPressurized(mesh, state, cr);
            cell_side_press = cell_side_press || pr;
        }
        if (nd >= 0) {
            const bool gp = ghostPressurized(mesh, state, f);
            if (cl < 0)
                pl = gp;
            if (cr < 0)
                pr = gp;
        }
        if (!cell_side_press)
            return kNone;
        return (pl && pr) ? kFull : kDelta;
    }

    static bool ghostPressurized(const NetworkMeshData& mesh, const NetworkStateData& state,
                                 int f) {
        const auto uf = static_cast<std::size_t>(f);
        const int nd = mesh.face_node[uf];
        if (nd < 0)
            return false;
        const int other = (mesh.face_cl[uf] >= 0) ? mesh.face_cl[uf] : mesh.face_cr[uf];
        if (other < 0)
            return false;
        const FvGeometry& g =
            mesh.geom[static_cast<std::size_t>(mesh.cell_geom[static_cast<std::size_t>(other)])];
        if (g.is_open)
            return false;
        // TPA (issue #156 Phase 5): a SEALED node against a flagged interior
        // cell carries the pressurized column at sub-atmospheric head — the
        // ghost is pressurized even though head < crown (the same rule
        // faceSide uses for its ghost regime). Without this the sealed
        // sub-crown junction's faces degrade to kDelta, the fold pass skips
        // the node, and the implicit path breaks the column (the P4 pinned
        // known-issue).
        if (!state.cell_tpa.empty() && state.cell_tpa[static_cast<std::size_t>(other)] != 0 &&
            mesh.node_sur_depth[static_cast<std::size_t>(nd)] > 0.0)
            return true;
        const double hg = state.node_head[static_cast<std::size_t>(nd)] - mesh.face_zb[uf];
        return hg >= g.y_crown;
    }
 
private:
    struct ImpFace {
        int face = -1;
        std::uint8_t mode = kNone;
        // Unknown index of each side (-1 = known/Dirichlet). Left/right in
        // the FACE frame, so the assembly's sign convention is the flux's.
        int unk_l = -1;
        int unk_r = -1;
        // Dirichlet head of a known side (valid where unk < 0).
        double dir_l = 0.0;
        double dir_r = 0.0;
        // Time-n heads of both sides (the DELTA form needs them).
        double eta_l = 0.0;
        double eta_r = 0.0;
        double alpha = 1.0; 
        double cond = 0.0;  
        double qstar = 0.0; 
    };
 
    // classify() products, persistent to avoid per-substep allocation.
    std::vector<std::uint8_t> cell_press_;
    std::vector<std::uint8_t> node_fold_;
    std::vector<int> cell_unk_; 
    std::vector<int> node_unk_; 
    std::vector<ImpFace> imp_faces_;
 
    // Unknown table: entity id (cell c → c; node n → n_cells + n), time-n
    // head, storage coefficient (T·Δx/Δt for cells, 0 for folded junctions —
    // filled with Δt at solve time), and component id.
    std::vector<int> unk_entity_;
    std::vector<double> unk_eta0_;
    std::vector<double> unk_store_;   
    std::vector<double> unk_rhs_ext_; 
    std::vector<int> unk_comp_;
    std::vector<double> unk_head_; 
 
    // Union-find over unknowns.
    std::vector<int> uf_parent_;
    int ufFind(int i);
    void ufUnion(int a, int b);
 
    // Per-component scratch.
    std::vector<int> comp_ptr_, comp_unks_;
    std::vector<double> thom_a_, thom_b_, thom_c_, thom_d_;
    std::vector<int> path_order_, degree_, adj_ptr_, adj_idx_;
    std::vector<double> adj_c_;
    std::vector<double> cg_r_, cg_p_, cg_ap_, cg_z_;

    void buildComponents();
    void prepare(const PressurizedView& v, double dt);
    void assembleRhs(const PressurizedView& v, double dt);
    void solveComponent(int comp, const PressurizedView& v, double dt);
    void backSubstitute(const PressurizedView& v);

    void linearSolve(const std::vector<int>& unks, const std::vector<double>& diag,
                     const std::vector<double>& rhs);
 
    // diag/rhs assembled per solve (indexed by unknown).
    std::vector<double> diag_, rhs_;
    std::vector<std::uint8_t> unk_dirichlet_; 
    std::vector<double> unk_dirvalue_;
    NonlinearPressureStorage candidate_;
    std::vector<NonlinearPressureStorage::Row> candidate_rows_;
    std::vector<NonlinearPressureStorage::Edge> candidate_edges_;
    std::vector<double> candidate_head_, candidate_residual_, candidate_budget_;
};
 
} // namespace openswmm::fv
 
#endif // OPENSWMM_ENGINE_FV_PRESSURIZED_HEAD_SOLVER_HPP