SurfaceStateData.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_2D_SURFACE_STATE_DATA_HPP
#define OPENSWMM_ENGINE_2D_SURFACE_STATE_DATA_HPP
 
#include <vector>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <vector>
 
#include "SurfaceTransportState.hpp"   // S1 — species mass per cell
#include "MeshData.hpp"                 // kMaxCellVerts edge-slot stride
 
namespace openswmm { struct NodeData; }  // 1D node data (held during a 2D advance)
 
namespace openswmm::twoD {
 
struct CouplingPoint;  // fwd decl — 1D↔2D coupling descriptor (NodeCoupling.hpp)

struct BoundaryData;  // fwd decl — per-edge boundary conditions (BoundaryData.hpp)
 
struct SurfaceStateData {
 
    // -----------------------------------------------------------------------
    // State variables — per triangle [0, n_triangles)
    // -----------------------------------------------------------------------

    const BoundaryData* boundary = nullptr;

    const NodeData*                   nodes_1d        = nullptr;
    const std::vector<double>*        node_row_conc   = nullptr;
    const std::vector<CouplingPoint>* node_coupling   = nullptr;  

    SurfaceTransportState transport;
 
    std::vector<double> depth;          
    std::vector<double> head;           
    std::vector<double> volume;         
 
    // Gradient fields (per triangle)
    std::vector<double> grad_hx;        
    std::vector<double> grad_hy;        
    std::vector<double> grad_hx_lim;    
    std::vector<double> grad_hy_lim;    
 
    // Reconstructed head at vertices — [0, n_vertices)
    std::vector<double> vert_head;      

    std::vector<double> vert_depth_signed;
 
    // Cell-centred velocity (RT0 reconstruction from edge fluxes) — per triangle
    std::vector<double> face_vx;        
    std::vector<double> face_vy;        
 
    // Per-cell continuity residual — per triangle (m³/s, ≈0 when conservative)
    std::vector<double> cell_continuity_err;
 
    // Fluxes — flat 2D: [cell * kMaxCellVerts + edge]
    std::vector<double> edge_flux;      
 
    // Source/sink terms — per triangle
    std::vector<double> rainfall;       
    std::vector<double> evap_rate;      

    std::vector<double> infil_rate;

    std::vector<double> infil_applied;
 
    std::vector<double> coupling_flux;  

    std::vector<double> coupling_applied;
    std::vector<double> net_source;     
 
    // -----------------------------------------------------------------------
    // Forcing overrides (optional external control)
    // -----------------------------------------------------------------------
 
    std::vector<int8_t> rainfall_forced;      
    std::vector<int8_t> rainfall_persist;     
    std::vector<double> rainfall_force_val;   
    std::vector<int8_t> evap_forced;          
    std::vector<int8_t> evap_persist;         
    std::vector<double> evap_force_val;       
    std::vector<int8_t> coupling_forced;      
    std::vector<int8_t> coupling_persist;     
    std::vector<double> coupling_force_val;   

    bool forcing_dirty = false;

    bool forcing_ever_set = false;
 
    // -----------------------------------------------------------------------
    // Previous step state
    // -----------------------------------------------------------------------
 
    std::vector<double> old_depth;      
    std::vector<double> old_volume;     
 
    // -----------------------------------------------------------------------
    // Cumulative statistics
    // -----------------------------------------------------------------------
 
    std::vector<double> stat_max_depth;     
    std::vector<double> stat_max_velocity;  
    std::vector<double> stat_max_cont_err;  
    std::vector<double> stat_cum_volume;    

    double evap_loss_total = 0.0;
 
    // -----------------------------------------------------------------------
    // Lifecycle
    // -----------------------------------------------------------------------
 
    void resize(int n_triangles, int n_vertices) {
        auto nt = static_cast<std::size_t>(n_triangles);
        auto nv = static_cast<std::size_t>(n_vertices);
        auto n3 = nt * static_cast<std::size_t>(kMaxCellVerts);
 
        depth.assign(nt, 0.0);
        head.assign(nt, 0.0);
        volume.assign(nt, 0.0);
        grad_hx.assign(nt, 0.0);
        grad_hy.assign(nt, 0.0);
        grad_hx_lim.assign(nt, 0.0);
        grad_hy_lim.assign(nt, 0.0);
        vert_head.assign(nv, 0.0);
        vert_depth_signed.assign(nv, 0.0);
        face_vx.assign(nt, 0.0);
        face_vy.assign(nt, 0.0);
        cell_continuity_err.assign(nt, 0.0);
        edge_flux.assign(n3, 0.0);
        rainfall.assign(nt, 0.0);
        evap_rate.assign(nt, 0.0);
        infil_rate.assign(nt, 0.0);
        infil_applied.assign(nt, 0.0);
        coupling_applied.assign(nt, 0.0);
        coupling_flux.assign(nt, 0.0);
        net_source.assign(nt, 0.0);
 
        rainfall_forced.assign(nt, 0);
        rainfall_persist.assign(nt, 0);
        rainfall_force_val.assign(nt, 0.0);
        evap_forced.assign(nt, 0);
        evap_persist.assign(nt, 0);
        evap_force_val.assign(nt, 0.0);
        coupling_forced.assign(nt, 0);
        coupling_persist.assign(nt, 0);
        coupling_force_val.assign(nt, 0.0);
 
        old_depth.assign(nt, 0.0);
        old_volume.assign(nt, 0.0);
        stat_max_depth.assign(nt, 0.0);
        stat_max_velocity.assign(nt, 0.0);
        stat_max_cont_err.assign(nt, 0.0);
        stat_cum_volume.assign(nt, 0.0);
        evap_loss_total = 0.0;
    }
 
    void save_state() noexcept {
        std::memcpy(old_depth.data(), depth.data(),
                     depth.size() * sizeof(double));
        std::memcpy(old_volume.data(), volume.data(),
                     volume.size() * sizeof(double));
    }
 
    void reset_state() noexcept {
        std::memcpy(depth.data(), old_depth.data(),
                     old_depth.size() * sizeof(double));
        std::memcpy(volume.data(), old_volume.data(),
                     old_volume.size() * sizeof(double));
    }

    void clear_reset_forcings() noexcept {
        if (!forcing_ever_set) return;
        for (std::size_t i = 0; i < rainfall_forced.size(); ++i) {
            if (rainfall_persist[i] == 0 && rainfall_forced[i] != 0) {
                rainfall_forced[i] = 0;
                rainfall_force_val[i] = 0.0;
                forcing_dirty = true;
            }
            if (evap_persist[i] == 0 && evap_forced[i] != 0) {
                evap_forced[i] = 0;
                evap_force_val[i] = 0.0;
                forcing_dirty = true;
            }
            if (coupling_persist[i] == 0 && coupling_forced[i] != 0) {
                coupling_forced[i] = 0;
                coupling_force_val[i] = 0.0;
                forcing_dirty = true;
            }
        }
    }

    void update_statistics([[maybe_unused]] const std::vector<double>& tri_area,
                           double dt,
                           [[maybe_unused]] int nthreads = 1) noexcept {
        // Each cell updates only its own envelope slots (max/cum into [i]);
        // schedule(static) keeps this bit-identical to serial for any thread
        // count. int loop index for OpenMP canonical-loop form.
        const int n = static_cast<int>(depth.size());
#if defined(SWMM_USE_OPENMP)
#pragma omp parallel for schedule(static) num_threads(nthreads)
#endif
        for (int i = 0; i < n; ++i) {
            if (depth[i] > stat_max_depth[i])
                stat_max_depth[i] = depth[i];
 
            const double speed = std::sqrt(face_vx[i] * face_vx[i]
                                         + face_vy[i] * face_vy[i]);
            if (speed > stat_max_velocity[i])
                stat_max_velocity[i] = speed;
 
            const double aerr = std::abs(cell_continuity_err[i]);
            if (aerr > stat_max_cont_err[i])
                stat_max_cont_err[i] = aerr;
 
            // Cell water volume × dt (VFR: V is the integrated state, not h̄·A).
            stat_cum_volume[i] += volume[i] * dt;
        }
    }
};
 
} // namespace openswmm::twoD
 
#endif // OPENSWMM_ENGINE_2D_SURFACE_STATE_DATA_HPP