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"
 
namespace openswmm::fv {

inline constexpr int kI1Samples = 129;

inline constexpr int kNodeVolSamples = 129;

struct FvGeometry {
    XSectParams xs{};              

    const xsect::XsectEval* eval = nullptr;
 
    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;           
 
    // -- Friction and losses (already lengthening-adjusted, see Router::init) --
 
    double roughness    = 0.01;    
    double rough_factor = 0.0;     
    double loss_inlet   = 0.0;     
    double loss_outlet  = 0.0;     
 
    int    barrels = 1;            

    int    culvert_code = 0;
    double slope        = 0.0;

    hydkernels::CulvertCurve culvert_curve{};
    uint8_t culvert_mitered = 0;

    double i1_tbl[2 * kI1Samples] = {};

    double h_tbl[kI1Samples] = {};
 
};

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;       
 
    // -----------------------------------------------------------------------
    // 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();
        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 (per BARREL — link reporting multiplies by
    // geom.barrels).
    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, cell-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