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