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>
#include <map>
 
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> gained_external, lost_external;
    std::vector<double>
        gained_bed; 
    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;
    std::map<int, std::vector<double>> runtime_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;
    std::vector<double> gained_exfiltration;
    std::vector<double>
        generated_age; 
    std::vector<double> gained_washoff;
    std::vector<double> cell_runoff_vol;
 
    // ---- 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);
        gained_external.assign(ns, 0.0);
        lost_external.assign(ns, 0.0);
        gained_bed.clear();
        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
        runtime_bc_conc.clear();
        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);
        gained_exfiltration.assign(ns, 0.0); // T7.1
        generated_age.clear();
        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] +
               lost_external[us] - gained_external[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) -
               (us < gained_washoff.size() ? gained_washoff[us] : 0.0) -
               (us < gained_exfiltration.size() ? gained_exfiltration[us] : 0.0) -
               (us < generated_age.size() ? generated_age[us] : 0.0);
    }
};
 
} // namespace openswmm::twoD
 
#endif // OPENSWMM_ENGINE_2D_SURFACE_TRANSPORT_STATE_HPP