SubsurfaceTransportState.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_SUBSURFACE_TRANSPORT_STATE_HPP
#define OPENSWMM_ENGINE_2D_SUBSURFACE_TRANSPORT_STATE_HPP
 
#include <cmath>
#include <cstddef>
#include <string>
#include <string_view>
#include <vector>
 
namespace openswmm::twoD {

struct SubsurfaceTransportState {
    // ---- rows (TransportPolicy's layout, copied at initialize) ------------
    double c_diff = 0.4; 
    int n_species = 0;
    int n_cells = 0;
    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;
    }
 
    // ---- the two stores ---------------------------------------------------
    // T7.2: a store holds the cell's TOTAL mass of a species — dissolved
    // plus whatever is sorbed on the grains. Only the dissolved part
    // travels, and `dissolvedFraction` is what every flux multiplies by, so
    // retardation is a property of the CONCENTRATION the channels read
    // rather than a factor sprinkled through them.
    std::vector<double> sat_mass;   
    std::vector<double> unsat_mass; 
    // In LAYERED mode this is the canonical unsaturated inventory. The
    // public bulk field above is its derived sum, never a second owner.
    int n_layers = 0;
    std::vector<double> layer_mass; 
    // Above-table receiving slots are fixed at producer time. Water is a
    // routing split of nacc/lacc, not additional storage. The paired masses
    // live here instead of nacc_mass/lacc_mass until the owner fires.
    std::vector<double> contact_water;                      
    std::vector<double> contact_node_out, contact_link_out; 
    std::vector<double> contact_node_mass, contact_link_mass; 
    std::size_t layerIdx(int s, int j, int c) const noexcept {
        return (std::size_t(s) * n_cells + c) * n_layers + j;
    }

    void refreshUnsat(int s, int c) noexcept {
        double sum = 0;
        for (int j = 0; j < n_layers; ++j)
            sum += layer_mass[layerIdx(s, j, c)];
        unsat_mass[idx(s, c)] = sum;
    }
 
    // ---- T7.2 resolved parameters ----------------------------------------
    std::vector<double> alpha_L; 
    std::vector<double> alpha_T; 
    std::vector<double> D_m;     
    std::vector<double> rho_b;   
    std::vector<double> kd;
    std::vector<double> decay;
    bool dispersion_on = true;
    long dispersion_limiter_binds = 0;

    double dissolvedFraction(int s, int c, double theta) const noexcept {
        if (s == age_row || s == temp_row)
            return 1.0;
        if (kd.empty() || !(theta > 0.0))
            return 1.0;
        const double k = kd[idx(s, c)];
        if (!(k > 0.0))
            return 1.0;
        const double R = 1.0 + rho_b[static_cast<std::size_t>(c)] * k / theta;
        return (R > 1.0) ? 1.0 / R : 1.0;
    }
 
    // ---- cross-cadence accumulators (paired with the water's) -------------
    std::vector<double> sacc_L, sacc_R;
    std::vector<double> xacc_from_surface;
    std::vector<double> xacc_to_surface;
 
    // ---- T7.4: the 1D seams ----------------------------------------------
    std::vector<double> nacc_mass;
    std::vector<double> lacc_mass;
    std::vector<double> lout_mass;
    std::vector<double> link_gain_shortfall;
    std::vector<double> node_out_mass;
 
    // ---- ledgers (per species, cumulative) --------------------------------
    // Signs follow the water ledger they mirror: a "lost_" row is mass that
    // left the aquifer, a "gained_" row is mass that arrived.
    std::vector<double> lost_deep;    
    std::vector<double> lost_node;    
    std::vector<double> lost_link;    
    std::vector<double> lost_dunne;   
    std::vector<double> lost_et;      
    std::vector<double> gained_infil; 
    std::vector<double> gained_node;  
    std::vector<double> gained_link;  
    std::vector<double> net_lateral;  
    std::vector<double> gained_source, lost_source; 
    std::vector<double> boundary_in, boundary_out, boundary_flux_in,
        boundary_flux_out;             
    std::vector<double> lost_reaction; 
    std::vector<double> generated_age,
        gained_contact; 
    std::vector<double> init_mass;
    std::vector<double> internal_recharge;

    bool active() const noexcept { return n_species > 0 && n_cells > 0; }

    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);
    }

    std::size_t lateralIndex(int species, std::size_t edge) const noexcept {
        return edge * static_cast<std::size_t>(n_species) + static_cast<std::size_t>(species);
    }
 

    void resize(int n_spec, int n_cell, int n_edges) {
        n_species = (n_spec > 0) ? n_spec : 0;
        n_cells = (n_cell > 0) ? n_cell : 0;
        const auto ns = static_cast<std::size_t>(n_species);
        const auto nc = static_cast<std::size_t>(n_cells);
        const auto ne = static_cast<std::size_t>(n_edges > 0 ? n_edges : 0);
        sat_mass.assign(ns * nc, 0.0);
        unsat_mass.assign(ns * nc, 0.0);
        sacc_L.assign(ns * ne, 0.0);
        sacc_R.assign(ns * ne, 0.0);
        xacc_from_surface.assign(ns * nc, 0.0);
        xacc_to_surface.assign(ns * nc, 0.0);
        lacc_mass.assign(ns * nc, 0.0); // T7.4
        lout_mass.assign(ns * nc, 0.0);
        link_gain_shortfall.assign(ns, 0.0);
        // nacc_mass and node_out_mass are sized by the solver, which knows
        // the bed and node counts.
        lost_deep.assign(ns, 0.0);
        lost_node.assign(ns, 0.0);
        lost_link.assign(ns, 0.0);
        lost_dunne.assign(ns, 0.0);
        lost_et.assign(ns, 0.0);
        gained_infil.assign(ns, 0.0);
        gained_node.assign(ns, 0.0);
        gained_link.assign(ns, 0.0);
        net_lateral.assign(ns, 0.0);
        gained_source.assign(ns, 0.0);
        lost_source.assign(ns, 0.0);
        boundary_in.assign(ns, 0);
        boundary_out.assign(ns, 0);
        boundary_flux_in.assign(ns, 0);
        boundary_flux_out.assign(ns, 0);
        lost_reaction.assign(ns, 0.0); // T7.2
        generated_age.clear();
        gained_contact.clear();
        init_mass.assign(ns, 0.0);
        internal_recharge.assign(ns, 0.0);
        alpha_L.assign(nc, 1.0); // T7.2: the row defaults
        alpha_T.assign(nc, 0.1);
        D_m.assign(nc, 1.0e-9);
        rho_b.assign(nc, 0.0);
        kd.assign(ns * nc, 0.0);
        decay.assign(ns * nc, 0.0);
        dispersion_limiter_binds = 0;
    }

    void clear() { *this = SubsurfaceTransportState{}; }

    double inFlight1D(int s) const noexcept {
        double m = 0.0;
        const auto ns = static_cast<std::size_t>(n_species > 0 ? n_species : 1);
        if (!nacc_mass.empty()) {
            const auto nb = nacc_mass.size() / ns;
            for (std::size_t b = 0; b < nb; ++b)
                m += nacc_mass[static_cast<std::size_t>(s) * nb + b];
        }
        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 += lacc_mass[base + static_cast<std::size_t>(c)];
        for (int j = 0; j < n_layers; ++j)
            for (int c = 0; c < n_cells; ++c) {
                const auto k = layerIdx(s, j, c);
                m += contact_node_mass[k] + contact_link_mass[k];
            }
        return m;
    }

    double storage(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 += sat_mass[base + static_cast<std::size_t>(c)];
            m += unsat_mass[base + static_cast<std::size_t>(c)];
            m += xacc_from_surface[base + static_cast<std::size_t>(c)];
            m += xacc_to_surface[base + static_cast<std::size_t>(c)];
        }
        const auto ne = sacc_L.size() / (n_species > 0 ? static_cast<std::size_t>(n_species) : 1);
        for (std::size_t e = 0; e < ne; ++e) {
            const auto k = lateralIndex(s, e);
            m += sacc_L[k] + sacc_R[k];
        }
        return m + inFlight1D(s); // T7.4
    }

    double ledgeredStorage(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 += sat_mass[base + static_cast<std::size_t>(c)] +
                 unsat_mass[base + static_cast<std::size_t>(c)];
        const auto ne = sacc_L.size() / (n_species > 0 ? static_cast<std::size_t>(n_species) : 1);
        for (std::size_t e = 0; e < ne; ++e) {
            const auto k = lateralIndex(s, e);
            m += sacc_L[k] + sacc_R[k];
        }
        // T7.4: the two 1D arrival accumulators are booked into
        // `gained_node` / `gained_link` only when the CELL gathers them, so
        // until then they are storage the ledger has not yet named — the
        // same phase rule the surface accumulators follow, and the reason
        // they are counted here while `xacc_*` are not.
        return m + inFlight1D(s);
    }

    double continuityDenominator(int s) const noexcept {
        const auto us = static_cast<std::size_t>(s);
        double total = std::fabs(init_mass[us]) + std::fabs(gained_infil[us]) +
                       std::fabs(gained_node[us]) + std::fabs(gained_link[us]) +
                       std::fabs(net_lateral[us]) + std::fabs(gained_source[us]);
        if (!generated_age.empty())
            total += std::abs(generated_age[us]) + std::abs(gained_contact[us]);
        return total;
    }

    double residual(int s) const noexcept {
        const auto us = static_cast<std::size_t>(s);
        const double out = lost_deep[us] + lost_node[us] + lost_link[us] + lost_dunne[us] +
                           lost_et[us] + lost_reaction[us] + lost_source[us];
        double in = gained_infil[us] + gained_node[us] + gained_link[us] + net_lateral[us] +
                    gained_source[us];
        if (!generated_age.empty())
            in += generated_age[us] + gained_contact[us];
        return ledgeredStorage(s) + out - in - init_mass[us];
    }
};
 
} // namespace openswmm::twoD
 
#endif // OPENSWMM_ENGINE_2D_SUBSURFACE_TRANSPORT_STATE_HPP