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

    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 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_;
};
 
} // namespace openswmm::fv
 
#endif // OPENSWMM_ENGINE_FV_PRESSURIZED_HEAD_SOLVER_HPP