SurfaceTransportState.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_TRANSPORT_STATE_HPP
#define OPENSWMM_ENGINE_2D_SURFACE_TRANSPORT_STATE_HPP
 
#include <cstddef>
#include <string>
#include <string_view>
#include <vector>
 
namespace openswmm::twoD {
 
struct SurfaceTransportState {
    int n_species = 0;   
    int n_cells   = 0;
 
    // ---- S4 row layout — the 1D engines' convention, verbatim -------------
    // rows [0, n_pollut) pollutants; [n_pollut, n_pollut + n_msx) MSX species
    // in ReactionData order; then the reserved __WATER_AGE__ row; then the
    // reserved __TEMPERATURE__ row LAST (H1 fixed the reported column order
    // so adding heat to an age model never moves the age column). Set by the
    // router after resize(); -1 ⇒ the row is not carried.
    int n_pollut = 0;
    int n_msx    = 0;
    int age_row  = -1;
    int temp_row = -1;
    std::vector<std::string> row_names;
    bool signedRow(int s) const noexcept { return s == temp_row; }
    int rowIndex(std::string_view name) const noexcept {
        for (std::size_t i = 0; i < row_names.size(); ++i)
            if (row_names[i] == name) return static_cast<int>(i);
        return -1;
    }

    std::vector<double> cell_mass;
 
    // ---- Ledgers (per species), m³·conc, cumulative over the run. ----------
    // These are what S3's coupling tuple and the continuity table will read;
    // in S1 they are the record that mass left the surface and where.
    std::vector<double> lost_infiltration;   
    std::vector<double> lost_boundary;       
    std::vector<double> lost_coupling;       
    std::vector<double> exch_mass;
    std::vector<double> exch_spill;
 
    // ---- S2 source concentrations (species units) -------------------------
    std::vector<double> rain_conc;
    std::vector<double> bc_conc;
    struct BoundaryQualityRow { int slot = -1; int species = -1; double conc = 0.0; };
    std::vector<BoundaryQualityRow> bc_quality_rows;
 
    // ---- S2 dispersion limiter telemetry ----------------------------------
    long dispersion_limiter_binds = 0;

    std::vector<double> gained_rainfall;
    std::vector<double> gained_boundary;
    std::vector<double> gained_coupling;
 
    // ---- S3 outfall-discharge species source ------------------------------
    std::vector<double> coupling_src;
 
    bool active() const noexcept { return n_species > 0 && n_cells > 0; }
 
    void resize(int n_spec, int n_tri, int n_coupling_points) {
        n_species = (n_spec > 0) ? n_spec : 0;
        n_cells   = (n_tri  > 0) ? n_tri  : 0;
        n_pollut = n_species; n_msx = 0; age_row = -1; temp_row = -1;   // S4: router refines
        row_names.clear();
        const auto ns = static_cast<std::size_t>(n_species);
        cell_mass.assign(ns * static_cast<std::size_t>(n_cells), 0.0);
        lost_infiltration.assign(ns, 0.0);
        lost_boundary.assign(ns, 0.0);
        lost_coupling.assign(ns, 0.0);
        exch_mass.assign(ns * static_cast<std::size_t>(
                                  n_coupling_points > 0 ? n_coupling_points
                                                        : 0),
                         0.0);
        exch_spill.assign(static_cast<std::size_t>(
                              n_coupling_points > 0 ? n_coupling_points : 0),
                          0.0);
        rain_conc.clear();          // S2: router fills when pollutants carry one
        bc_conc.clear();            // S2: solver sizes at initialize
        bc_quality_rows.clear();
        dispersion_limiter_binds = 0;
        gained_rainfall.assign(ns, 0.0);
        gained_boundary.assign(ns, 0.0);
        gained_coupling.assign(ns, 0.0);
        coupling_src.clear();       // S3: router sizes when outfalls exist
    }
 
    void clear() { *this = SurfaceTransportState{}; }
 
    std::size_t idx(int s, int c) const noexcept {
        return static_cast<std::size_t>(s) * static_cast<std::size_t>(n_cells) +
               static_cast<std::size_t>(c);
    }

    double concentration(int s, int c, double volume_m3,
                         double dry_volume_m3) const noexcept {
        if (!(volume_m3 > dry_volume_m3)) return 0.0;
        return cell_mass[idx(s, c)] / volume_m3;
    }

    double totalIncludingLedgers(int s) const noexcept {
        double m = 0.0;
        const auto base = static_cast<std::size_t>(s) *
                          static_cast<std::size_t>(n_cells);
        for (int c = 0; c < n_cells; ++c)
            m += cell_mass[base + static_cast<std::size_t>(c)];
        const auto us = static_cast<std::size_t>(s);
        // Sources are SUBTRACTED so the quantity is "what was there at t=0":
        // surface + everything that left − everything that arrived.
        return m + lost_infiltration[us] + lost_boundary[us] +
               lost_coupling[us] -
               (us < gained_rainfall.size() ? gained_rainfall[us] : 0.0) -
               (us < gained_boundary.size() ? gained_boundary[us] : 0.0) -
               (us < gained_coupling.size() ? gained_coupling[us] : 0.0);
    }
};
 
}  // namespace openswmm::twoD
 
#endif  // OPENSWMM_ENGINE_2D_SURFACE_TRANSPORT_STATE_HPP