NetworkMeshData.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_NETWORK_MESH_DATA_HPP
#define OPENSWMM_ENGINE_FV_NETWORK_MESH_DATA_HPP
 
#include <cstdint>
#include <vector>
 
#include "../XSectBatch.hpp"
#include "../HydClosureKernels.hpp"
#include "../XSectKernels.hpp"
#include "FvClosureKernels.hpp"
 
namespace openswmm::fv {

inline constexpr int kNodeVolSamples = 129;

struct FvGeometry {
    XSectParams xs{}; 
 
    double y_full = 0.0; 
    double a_full = 0.0; 
    double w_max = 0.0;  
    double r_full = 0.0; 

    double barrel_scale = 1.0;

    double y_crown = 0.0;

    double t_slot = 0.0;

    double a_crown = 0.0;

    double i1_crown = 0.0;
 
    uint8_t is_open = 0; 
 
    int barrels = 1; 
 
    // Friction, losses, slope and the culvert curve are PER CONDUIT and live
    // in NetworkMeshData::conduit_* (plan Phase 1f): this block is one shared
    // entry per DISTINCT (section, barrels, open) so a 5000-pipe network with
    // a dozen sections keeps a dozen closure tables cache-resident instead of
    // streaming ~4 kB of geometry per cell every substep.

    FvClosure closure_tbl{};
};

struct NetworkMeshData {
 
    // -----------------------------------------------------------------------
    // Cells — indexed by cell index [0, n_cells)
    // -----------------------------------------------------------------------
 
    std::vector<int> cell_geom;    
    std::vector<int> cell_conduit; 
    std::vector<double> cell_dx;   
    std::vector<double> cell_zb;   

    std::vector<double> cell_dzdx;

    std::vector<int> cell_face0;
    std::vector<int> cell_face1;
    std::vector<int8_t> cell_side0;
    std::vector<int8_t> cell_side1;
 
    // -----------------------------------------------------------------------
    // Faces — indexed by face index [0, n_faces). Flux is positive LEFT → RIGHT.
    // -----------------------------------------------------------------------
    // A face with cell_l == -1 has a NODE on its left (a conduit's upstream
    // boundary); cell_r == -1 puts the node on its right (downstream boundary).
    // Interior faces (both cells ≥ 0) include virtual junctions, which are
    // nothing but a splice of two conduits' chains — mass and momentum flux
    // continuity are then properties of the scheme, not a special treatment.
 
    std::vector<int> face_cl;   
    std::vector<int> face_cr;   
    std::vector<int> face_node; 

    std::vector<int> face_geom; 

    void deriveFaceGeom() {
        if (face_geom.size() == static_cast<std::size_t>(n_faces()))
            return;
        face_geom.assign(static_cast<std::size_t>(n_faces()), 0);
        for (int f = 0; f < n_faces(); ++f) {
            const auto uf = static_cast<std::size_t>(f);
            const int c = (face_cl[uf] >= 0) ? face_cl[uf] : face_cr[uf];
            face_geom[uf] = (c >= 0 && static_cast<std::size_t>(c) < cell_geom.size())
                                ? cell_geom[static_cast<std::size_t>(c)]
                                : 0;
        }
    }

    std::vector<uint8_t> face_gate;

    std::vector<int> face_culvert;
    std::vector<double> face_zb; 
    std::vector<double> face_dx; 

    std::vector<int8_t> face_dir_l;
    std::vector<int8_t> face_dir_r;

    std::vector<uint8_t> face_virtual;

    std::vector<int> face_vj_node;

    std::vector<int> node_vj_face;
 
    // -----------------------------------------------------------------------
    // Conduit → cell map. Cells of a conduit are CONTIGUOUS by construction,
    // so a begin/count pair is a complete (and cheaper) CSR.
    // -----------------------------------------------------------------------
 
    std::vector<int> conduit_cell_begin; 
    std::vector<int> conduit_cell_count; 
    std::vector<int> conduit_link;       
    std::vector<int> conduit_section;    
 
    // Per-conduit friction and losses (already lengthening-adjusted, see
    // Router::init) and the culvert inlet-control data — everything that
    // differs between conduits sharing one section block (plan Phase 1f).
    std::vector<double> conduit_roughness;    
    std::vector<double> conduit_rough_factor; 
    std::vector<double> conduit_loss_inlet;   
    std::vector<double> conduit_loss_outlet;  
    std::vector<double> conduit_slope;        
    std::vector<int> conduit_culvert_code;    
    std::vector<hydkernels::CulvertCurve> conduit_culvert_curve;
    std::vector<uint8_t> conduit_culvert_mitered;
 
    // -----------------------------------------------------------------------
    // Cell chains (CSR). A chain is a maximal run of cells joined by INTERIOR
    // faces, so it spans virtual junctions — a spliced pair of conduits is one
    // chain, which is exactly the property that makes a virtual junction
    // indistinguishable from an interior cut.
    //
    // Chains give the solver an ordered 1D stencil: the second-order scalar
    // reconstruction reads neighbours along the chain, and the implicit
    // dispersion solve is one tridiagonal system per chain (D-FV1). `dir` is
    // the cell's own axis relative to the chain's walking direction.
    // -----------------------------------------------------------------------
 
    std::vector<int> chain_ptr;      
    std::vector<int> chain_cells;    
    std::vector<int8_t> chain_dir;   
    std::vector<int> cell_chain;     
    std::vector<int> cell_chain_pos; 

    int n_chains() const noexcept {
        return chain_ptr.empty() ? 0 : static_cast<int>(chain_ptr.size()) - 1;
    }
 
    // -----------------------------------------------------------------------
    // Non-conduit links (pumps, orifices, weirs, outlets). Evaluated by their
    // existing structure equations outside the solver and applied here as
    // source/sink pairs on the two node volumes.
    //
    // DUMMY-xsect conduits belong here too. Legacy isTrueConduit
    // (dynwave.c:411-414) is false for them: they carry no cross-section, so
    // there is nothing to march, and DW routes them through findNonConduitFlow
    // as a pure pass-through instead. They differ from a real structure in that
    // their discharge is not a head relation the engine can evaluate outside
    // the solver — it is whatever arrives at the upstream node — so the solver
    // computes it itself (`struct_is_dummy`, see ExplicitFvSolver::
    // refreshDummyFlows) rather than reading FvStepForcing::structure_flow.
    // -----------------------------------------------------------------------
 
    std::vector<int> struct_link; 
    std::vector<int> struct_n1;   
    std::vector<int> struct_n2;   

    std::vector<uint8_t> struct_is_dummy;

    std::vector<uint8_t> node_dummy_drain;
 
    // -----------------------------------------------------------------------
    // Node → face map (CSR). Two-pass gather with NO atomics: the node update
    // reads its own face list in a fixed order, so the result is deterministic
    // and identical across backends and thread counts.
    // -----------------------------------------------------------------------
 
    std::vector<int> node_face_ptr;     
    std::vector<int> node_face_idx;     
    std::vector<double> node_face_sign; 

    std::vector<double> node_face_zb;
 
    // -----------------------------------------------------------------------
    // Node properties
    // -----------------------------------------------------------------------
 
    std::vector<double> node_invert;      
    std::vector<double> node_full_depth;  
    std::vector<double> node_ponded_area; 
    std::vector<double> node_sur_depth;   

    std::vector<uint8_t> node_can_pond;

    std::vector<uint8_t> node_kind;

    std::vector<double> node_area;
 
    // -- Tabulated storage relation (STORAGE nodes only) ---------------------
    //
    // The solver must invert volume → depth every substep, and it must do so
    // without reaching back into SimulationContext (the GPU plugin has no
    // access to it). Storage geometry is therefore FLATTENED at mesh build into
    // a monotone depth→volume table sampled from node::getVolume, exactly as
    // the cross-sections are flattened for the device path.
    //
    // Junctions need none of this: their storage is linear in depth with the
    // surface area held over the routing step (plan §3.4).
 
    std::vector<int> node_vol_off;     
    std::vector<double> node_vol_dmax; 
    std::vector<double> node_vol_atop; 
    std::vector<double> node_vol_tbl;  
 
    // -----------------------------------------------------------------------
    // Geometry closures — one per conduit row
    // -----------------------------------------------------------------------
 
    std::vector<FvGeometry> geom;
 
    // -----------------------------------------------------------------------
    // Capacity queries
    // -----------------------------------------------------------------------

    int n_cells() const noexcept { return static_cast<int>(cell_dx.size()); }
    int n_faces() const noexcept { return static_cast<int>(face_cl.size()); }
    int n_nodes() const noexcept { return static_cast<int>(node_invert.size()); }
    int n_conduits() const noexcept { return static_cast<int>(conduit_cell_begin.size()); }

    void clear() {
        cell_geom.clear();
        cell_conduit.clear();
        cell_dx.clear();
        cell_zb.clear();
        cell_dzdx.clear();
        cell_face0.clear();
        cell_face1.clear();
        cell_side0.clear();
        cell_side1.clear();
        face_cl.clear();
        face_cr.clear();
        face_node.clear();
        face_gate.clear();
        face_geom.clear();
        face_culvert.clear();
        face_zb.clear();
        face_dx.clear();
        face_virtual.clear();
        face_vj_node.clear();
        node_vj_face.clear();
        face_dir_l.clear();
        face_dir_r.clear();
        conduit_cell_begin.clear();
        conduit_cell_count.clear();
        conduit_link.clear();
        conduit_section.clear();
        conduit_roughness.clear();
        conduit_rough_factor.clear();
        conduit_loss_inlet.clear();
        conduit_loss_outlet.clear();
        conduit_slope.clear();
        conduit_culvert_code.clear();
        conduit_culvert_curve.clear();
        conduit_culvert_mitered.clear();
        chain_ptr.clear();
        chain_cells.clear();
        chain_dir.clear();
        cell_chain.clear();
        cell_chain_pos.clear();
        struct_link.clear();
        struct_n1.clear();
        struct_n2.clear();
        struct_is_dummy.clear();
        node_dummy_drain.clear();
        node_face_ptr.clear();
        node_face_idx.clear();
        node_face_sign.clear();
        node_face_zb.clear();
        node_invert.clear();
        node_full_depth.clear();
        node_ponded_area.clear();
        node_sur_depth.clear();
        node_can_pond.clear();
        node_kind.clear();
        node_area.clear();
        node_vol_off.clear();
        node_vol_dmax.clear();
        node_vol_atop.clear();
        node_vol_tbl.clear();
        geom.clear();
    }
};

enum : uint8_t { kNodeJunction = 0, kNodeVirtual = 1, kNodeStorage = 2, kNodeOutfall = 3 };

struct NetworkStateData {
    // Conserved cell variables aggregate every parallel barrel. Link reporting
    // uses these directly; it must not multiply them by geom.barrels again.
    std::vector<double> cell_a; 
    std::vector<double> cell_q; 
 
    // Derived, refreshed after every cell update so faces read a consistent
    // depth without re-inverting A(h) per face.
    std::vector<double> cell_h; 
 
    // Node state. `volume` is the prognostic mass ledger — depth follows from
    // it through the node's storage relation — so node mass is conserved
    // exactly (plan §3.4).
    std::vector<double> node_volume;    
    std::vector<double> node_head;      
    std::vector<double> node_surf_area; 
    std::vector<double> node_overflow;  

    std::vector<double> node_overflow_vol;
 
    // Advected scalars, species-major: phi[s * n_cells + c]. Empty when transport
    // is off.
    std::vector<double> cell_phi;
    int n_species = 0;

    std::vector<uint8_t> cell_tpa;

    void resize(int n_cells, int n_nodes, int n_species_in) {
        const auto nc = static_cast<std::size_t>(n_cells);
        const auto nn = static_cast<std::size_t>(n_nodes);
        cell_a.assign(nc, 0.0);
        cell_q.assign(nc, 0.0);
        cell_h.assign(nc, 0.0);
        node_volume.assign(nn, 0.0);
        node_head.assign(nn, 0.0);
        node_surf_area.assign(nn, 0.0);
        node_overflow.assign(nn, 0.0);
        node_overflow_vol.assign(nn, 0.0);
        n_species = n_species_in;
        cell_phi.assign(nc * static_cast<std::size_t>(n_species_in), 0.0);
    }
};
 
} // namespace openswmm::fv
 
#endif // OPENSWMM_ENGINE_FV_NETWORK_MESH_DATA_HPP