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 <limits>
#include <algorithm>
#include <vector>
 
#include "../coupling/RuntimeCoupling.hpp"
#include "../../hydrology/surface/MeshSnow.hpp"
#include "SurfaceTransportState.hpp" // S1 — species mass per cell
#include "../heat/MeshHeat.hpp"
#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; 
    const std::vector<double>* node_pond_base = nullptr;
    const std::vector<int>* coupling_head_cells = nullptr;
    const std::vector<double>* node_pond_area_m2 = nullptr;
    double coupling_window_s = 0.0;

    SurfaceTransportState transport;
    MeshHeatState heat;
 
    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;
 
    // Output-only cell velocity: native FULL_SWE or wetted-edge reconstruction.
    // Paired NaNs mark an unresolved wet-cell vector; dry cells carry zero.
    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
    // Snow owns precipitation only on its resolved weather footprints. The
    // context is a non-owning forcing view, like nodes_1d; it is never archived.
    surface::MeshSnow snow;
    SimulationContext* snow_context = nullptr;
    std::vector<double> snow_catch;
    double liquidRainfall(int i) const {
        return snow.owns(i) ? snow.liquidRates()[i] : runtime_forcings.rate(1, i, rainfall[i]);
    }

    double exposedEvaporation(int i) const {
        return runtime_forcings.rate(2, i, evap_rate[i]) *
               (snow.owns(i) ? snow.exposedFractions()[i] : 1.);
    }
    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;
 
    // Plugin boundary transfers (m³), compacted in boundary_receipt_slots order. Allocated
    // only by plugin solvers; additive until the router drains them into its
    // engine-owned receipts after advance(). Separate directions preserve
    // reversals within a window. CPU solvers book their receipts directly.
    std::vector<int> boundary_receipt_slots; 
    std::vector<double> boundary_in_applied;
    std::vector<double> boundary_out_applied;
    std::vector<double> boundary_requested; 
    std::vector<double> boundary_last_flow; 
    std::vector<double> boundary_last_dt;   
 
    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; 

    // Unforced coupling is restored before each batch; ADD never compounds.
    RuntimeSources runtime_sources;
    RuntimeForcings runtime_forcings;
    std::vector<double> coupling_native;
    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;
    std::vector<double> rk2_rain_pending;
 
    // -----------------------------------------------------------------------
    // Lifecycle
    // -----------------------------------------------------------------------

    void resize(int n_triangles, int n_vertices) {
        heat.clear();
        snow = {};
        snow_context = nullptr;
        snow_catch.clear();
        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);
        boundary_receipt_slots.clear();
        boundary_in_applied.clear();
        boundary_out_applied.clear();
        boundary_requested.clear();
        boundary_last_flow.clear();
        boundary_last_dt.clear();
        coupling_flux.assign(nt, 0.0);
        coupling_native.clear();
        runtime_sources = {};
        runtime_forcings = {};
        rk2_rain_pending.clear();
        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, std::numeric_limits<double>::quiet_NaN());
        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 (std::isfinite(speed) &&
                (!std::isfinite(stat_max_velocity[i]) || 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