XSectBatch.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_XSECT_BATCH_HPP
#define OPENSWMM_XSECT_BATCH_HPP
 
#ifndef OPENSWMM_RESTRICT
#  if defined(_MSC_VER)
#    define OPENSWMM_RESTRICT __restrict
#  else
#    define OPENSWMM_RESTRICT __restrict__
#  endif
#endif
 
#include <vector>
#include <cstddef>
#include <cstdint>
#include <cmath>
 
#include "XSectLookup.hpp"
 
namespace openswmm {
 
// Forward declaration
struct SimulationContext;
 
// ============================================================================
// Cross-section shape codes (matches legacy enums.h XsectType)
// ============================================================================
 
enum class XSectShape : int {
    DUMMY              =  0,
    CIRCULAR           =  1,
    FILLED_CIRCULAR    =  2,
    RECT_CLOSED        =  3,
    RECT_OPEN          =  4,
    TRAPEZOIDAL        =  5,
    TRIANGULAR         =  6,
    PARABOLIC          =  7,
    POWERFUNC          =  8,
    RECT_TRIANG        =  9,
    RECT_ROUND         = 10,
    MOD_BASKET         = 11,
    HORIZ_ELLIPSE      = 12,
    VERT_ELLIPSE       = 13,
    ARCH               = 14,
    EGGSHAPED          = 15,
    HORSESHOE          = 16,
    GOTHIC             = 17,
    CATENARY           = 18,
    SEMIELLIPTICAL     = 19,
    BASKETHANDLE       = 20,
    SEMICIRCULAR       = 21,
    IRREGULAR          = 22,
    CUSTOM             = 23,
    FORCE_MAIN         = 24,
    STREET_XSECT       = 25
};
 
// ============================================================================
// Cross-section parameter struct (mirrors legacy TXsect)
// ============================================================================
 
struct XSectParams {
    int    type         = 0;
    int    culvert_code = 0;
    int    transect     = -1;
 
    double y_full  = 0.0;      
    double w_max   = 0.0;      
    double yw_max  = 0.0;      
    double a_full  = 0.0;      
    double r_full  = 0.0;      
    double s_full  = 0.0;      
    double s_max   = 0.0;      
 
    double y_bot   = 0.0;      
    double a_bot   = 0.0;      
    double s_bot   = 0.0;      
    double r_bot   = 0.0;      
 
    // Tabulated geometry for IRREGULAR / CUSTOM / STREET_XSECT shapes, whose
    // A/R/W vs depth come from a per-link transect table rather than a shared
    // static table.  These point into ctx.transect_tables (stable for the run)
    // and stay null for every self-contained shape.  Without them the
    // context-free per-element accessors return 0 for irregular sections (so
    // e.g. init-time getDepthFromFlow could not compute normal depth/storage).
    const double* area_tbl  = nullptr;   
    const double* hrad_tbl  = nullptr;   
    const double* width_tbl = nullptr;   
    int    transect_tbl_size = 0;        

    const xsect::LocateLut* area_lut = nullptr;
};
 
// ============================================================================
// Per-element functions (xsect:: namespace)
// ============================================================================
 
namespace xsect {
 
double getAofY(const XSectParams& xs, double y);
double getRofY(const XSectParams& xs, double y);
double getWofY(const XSectParams& xs, double y);
double getYofA(const XSectParams& xs, double a);
double getSofA(const XSectParams& xs, double a);
double getRofA(const XSectParams& xs, double a);
double getdSdA(const XSectParams& xs, double a);
double getAofS(const XSectParams& xs, double s_factor);
double getAmax(const XSectParams& xs);
double getYcrit(const XSectParams& xs, double q);
bool   isOpen(int type);
int    setParams(XSectParams& xs, int type, const double p[], double ucf);
 
// Lookup table helpers (exposed for batch kernels and testing)
double lookup(double x, const double* table, int n_items);
double invLookup(double y, const double* table, int n_items,
                 const LocateLut* lut = nullptr);
int    locate(double y, const double* table, int n);
double getYcircular(double alpha);
double getScircular(double alpha);
 
} // namespace xsect
 
// ============================================================================
// Shape group — contiguous SoA for all links sharing one shape type
// ============================================================================

struct ShapeGroup {
    XSectShape shape = XSectShape::DUMMY;
    int        count = 0;
 
    // Mapping back to global link arrays
    std::vector<int> link_idx;      
 
    // Geometry parameters (contiguous, aligned for SIMD)
    std::vector<double> y_full;     
    std::vector<double> a_full;     
    std::vector<double> r_full;     
    std::vector<double> s_full;     
    std::vector<double> w_max;      
 
    // Pre-computed reciprocal of y_full (avoids per-element division in kernels)
    std::vector<double> inv_y_full;    
 
    // Multi-purpose parameters (meaning depends on shape)
    std::vector<double> y_bot;
    std::vector<double> a_bot;
    std::vector<double> s_bot;
    std::vector<double> r_bot;
 
    // Per-link transect table pointers (IRREGULAR shapes only)
    // Each pointer → a normalized table of N_TRANSECT_TBL entries.
    std::vector<const double*> area_tables;   
    std::vector<const double*> hrad_tables;   
    std::vector<const double*> width_tables;  
    int transect_tbl_size = 0;                
 
    // Pre-allocated working buffers (avoids per-call allocation in hot loop)
    mutable std::vector<double> buf_d;   
    mutable std::vector<double> buf_r;   
    mutable std::vector<double> buf_r2;  

    void resize(int n);
};
 
// ============================================================================
// XSectGroups — the shape-grouped index over all links
// ============================================================================

class XSectGroups {
public:
    void build(const SimulationContext& ctx);

    void attachTransectTables(const SimulationContext& ctx);

    void build(const XSectParams* params, int n_links);
 
    // ========================================================================
    // Batch compute — results scattered to global link arrays
    // ========================================================================

    void computeAreas(const double* depths, double* areas, int n_links) const;

    void computeHydRad(const double* depths, double* hydrad, int n_links) const;

    void computeWidths(const double* depths, double* widths, int n_links) const;

    void computeAreaAndHydRad(const double* depths, double* areas,
                              double* hydrad, int n_links) const;

    void computeAreaHydRadTriple(
        const double* d1, const double* d2, const double* dm,
        double* a1, double* a2, double* am,
        double* hrad1, double* hrad_mid, int n_links) const;

    void computeWidthsTriple(
        const double* d1, const double* d2, const double* dm,
        double* w1, double* w2, double* wm, int n_links) const;

    void computeAreaHydRadTripleTeam(
        const double* d1, const double* d2, const double* dm,
        double* a1, double* a2, double* am,
        double* hrad1, double* hrad_mid, int n_links,
        int tid, int nthreads) const;
 
    void computeWidthsTripleTeam(
        const double* d1, const double* d2, const double* dm,
        double* w1, double* w2, double* wm, int n_links,
        int tid, int nthreads) const;

    void setBypassMask(const std::uint8_t* bypassed_by_link) const;

    void computeSectionFactors(const double* areas, double* sfact, int n_links) const;

    void computeDepthsFromArea(const double* areas, double* depths, int n_links) const;
 
    // ========================================================================
    // Accessors
    // ========================================================================

    int numGroups() const { return static_cast<int>(groups_.size()); }

    const ShapeGroup& group(int i) const { return groups_[i]; }

    const ShapeGroup* findGroup(XSectShape shape) const;
 
private:
    std::vector<ShapeGroup> groups_;
 
    // Per-iteration bypass-mask state (see setBypassMask).  packed_groups_[gi]
    // is a packed view of groups_[gi] holding only the active links; it reuses
    // the ShapeGroup layout so the gather/kernel/scatter path runs unchanged.
    // packed_count_[gi]:  kMaskFullGroup → no link in the group is bypassed
    // (use the original group directly, no packing was done); 0..count →
    // number of active links in the packed view.  All mutable because the
    // mask is per-Picard-iteration scratch while the compute API is const.
    static constexpr int kMaskFullGroup = -1;
    mutable std::vector<ShapeGroup> packed_groups_;
    mutable std::vector<int>        packed_count_;
    mutable bool                    mask_active_ = false;

    const ShapeGroup* maskedGroup(std::size_t gi) const {
        const ShapeGroup* gp = &groups_[gi];
        if (mask_active_) {
            const int na = packed_count_[gi];
            if (na == 0) return nullptr;
            if (na != kMaskFullGroup) gp = &packed_groups_[gi];
        }
        return gp;
    }
};
 
// ============================================================================
// Shape-specific batch kernels (called by XSectGroups, also usable directly)
// ============================================================================
 
namespace xsect_batch {

void area_circular(
    const double* OPENSWMM_RESTRICT depth,
    const double* OPENSWMM_RESTRICT y_full,
    const double* OPENSWMM_RESTRICT a_full,
    double*       OPENSWMM_RESTRICT area,
    int count
);

void area_rect(
    const double* OPENSWMM_RESTRICT depth,
    const double* OPENSWMM_RESTRICT w_max,
    double*       OPENSWMM_RESTRICT area,
    int count
);

void area_trapezoidal(
    const double* OPENSWMM_RESTRICT depth,
    const double* OPENSWMM_RESTRICT y_bot,
    const double* OPENSWMM_RESTRICT s_bot,
    double*       OPENSWMM_RESTRICT area,
    int count
);

void area_triangular(
    const double* OPENSWMM_RESTRICT depth,
    const double* OPENSWMM_RESTRICT s_bot,
    double*       OPENSWMM_RESTRICT area,
    int count
);

void area_parabolic(
    const double* OPENSWMM_RESTRICT depth,
    const double* OPENSWMM_RESTRICT r_bot,
    double*       OPENSWMM_RESTRICT area,
    int count
);

void area_powerfunc(
    const double* OPENSWMM_RESTRICT depth,
    const double* OPENSWMM_RESTRICT s_bot,
    const double* OPENSWMM_RESTRICT r_bot,
    double*       OPENSWMM_RESTRICT area,
    int count
);

void area_tabulated(
    const double* OPENSWMM_RESTRICT depth,
    const double* OPENSWMM_RESTRICT y_full,
    const double* OPENSWMM_RESTRICT a_full,
    const double* table,
    int            table_size,
    double*       OPENSWMM_RESTRICT area,
    int count
);

void perlink_tabulated(
    const double* OPENSWMM_RESTRICT depth,
    const double* OPENSWMM_RESTRICT nrm,
    const double* OPENSWMM_RESTRICT scale,
    const double* const* tables,
    int            table_size,
    double*       OPENSWMM_RESTRICT result,
    int count
);

void perlink_tabulated_pair(
    const double* OPENSWMM_RESTRICT depth,
    const double* OPENSWMM_RESTRICT nrm,
    const double* OPENSWMM_RESTRICT scale_a,
    const double* OPENSWMM_RESTRICT scale_b,
    const double* const* tables_a,
    const double* const* tables_b,
    int            table_size,
    double*       OPENSWMM_RESTRICT out_a,
    double*       OPENSWMM_RESTRICT out_b,
    int count
);

void area_inv_tabulated(
    const double* OPENSWMM_RESTRICT depth,
    const double* OPENSWMM_RESTRICT y_full,
    const double* OPENSWMM_RESTRICT a_full,
    const double* table,
    int            table_size,
    const xsect::LocateLut* lut,
    double*       OPENSWMM_RESTRICT area,
    int count
);
 
// --- Hydraulic radius batch kernels ---
 
void hydrad_circular(
    const double* OPENSWMM_RESTRICT depth,
    const double* OPENSWMM_RESTRICT y_full,
    const double* OPENSWMM_RESTRICT r_full,
    double*       OPENSWMM_RESTRICT hydrad,
    int count
);

void area_hydrad_circular(
    const double* OPENSWMM_RESTRICT depth,
    const double* OPENSWMM_RESTRICT y_full,
    const double* OPENSWMM_RESTRICT a_full,
    const double* OPENSWMM_RESTRICT r_full,
    double*       OPENSWMM_RESTRICT area,
    double*       OPENSWMM_RESTRICT hydrad,
    int count
);
 
void hydrad_trapezoidal(
    const double* OPENSWMM_RESTRICT depth,
    const double* OPENSWMM_RESTRICT y_bot,
    const double* OPENSWMM_RESTRICT s_bot,
    const double* OPENSWMM_RESTRICT r_bot,
    double*       OPENSWMM_RESTRICT hydrad,
    int count
);
 
void hydrad_triangular(
    const double* OPENSWMM_RESTRICT depth,
    const double* OPENSWMM_RESTRICT s_bot,
    const double* OPENSWMM_RESTRICT r_bot,
    double*       OPENSWMM_RESTRICT hydrad,
    int count
);
 
void hydrad_rect(
    const double* OPENSWMM_RESTRICT depth,
    const double* OPENSWMM_RESTRICT w_max,
    double*       OPENSWMM_RESTRICT hydrad,
    int count
);
 
void hydrad_tabulated(
    const double* OPENSWMM_RESTRICT depth,
    const double* OPENSWMM_RESTRICT y_full,
    const double* OPENSWMM_RESTRICT r_full,
    const double* table,
    int            table_size,
    double*       OPENSWMM_RESTRICT hydrad,
    int count
);
 
// --- Top width batch kernels ---
 
void width_circular(
    const double* OPENSWMM_RESTRICT depth,
    const double* OPENSWMM_RESTRICT y_full,
    const double* OPENSWMM_RESTRICT w_max,
    double*       OPENSWMM_RESTRICT width,
    int count
);
 
void width_trapezoidal(
    const double* OPENSWMM_RESTRICT depth,
    const double* OPENSWMM_RESTRICT y_bot,
    const double* OPENSWMM_RESTRICT s_bot,
    double*       OPENSWMM_RESTRICT width,
    int count
);
 
void width_triangular(
    const double* OPENSWMM_RESTRICT depth,
    const double* OPENSWMM_RESTRICT s_bot,
    double*       OPENSWMM_RESTRICT width,
    int count
);
 
void width_rect(
    const double* OPENSWMM_RESTRICT w_max,
    double*       OPENSWMM_RESTRICT width,
    int count
);
 
void width_tabulated(
    const double* OPENSWMM_RESTRICT depth,
    const double* OPENSWMM_RESTRICT y_full,
    const double* OPENSWMM_RESTRICT w_max,
    const double* table,
    int            table_size,
    double*       OPENSWMM_RESTRICT width,
    int count
);
 
} // namespace xsect_batch
 
} // namespace openswmm
 
#endif // OPENSWMM_XSECT_BATCH_HPP