SegmentStore.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_QUALITY_LARD_SEGMENT_STORE_HPP
#define OPENSWMM_QUALITY_LARD_SEGMENT_STORE_HPP
#include "../../core/CompletedSourceArchive.hpp"
 
#include <algorithm>
#include <cstddef>
#include <vector>
 
namespace openswmm {
namespace lard {

constexpr int kMaxSegmentsPerLink = 100;

constexpr double kMergeAtol = 1.0e-6;
constexpr double kMergeRtol = 1.0e-4;
 
class SegmentStore {
public:
    void checkpoint(detail::CompletedSourceArchive& a) {
        a.fixed(head_, count_, vol_, conc_);
        for (int l = 0; l < nl_; ++l)
            if (head_[l] < 0 || head_[l] >= cap_ || count_[l] < 0 || count_[l] > cap_)
                throw std::runtime_error("Invalid LARD checkpoint ring.");
        for (double v : vol_)
            if (!std::isfinite(v) || v < 0)
                throw std::runtime_error("Invalid LARD checkpoint water volume.");
        for (double v : conc_)
            if (!std::isfinite(v))
                throw std::runtime_error("Nonfinite LARD checkpoint concentration.");
    }

    void resize(int n_links, int n_species, int cap = kMaxSegmentsPerLink) {
        nl_ = n_links;
        ns_ = n_species;
        cap_ = cap;
        const std::size_t slab = static_cast<std::size_t>(nl_) * static_cast<std::size_t>(cap_);
        head_.assign(static_cast<std::size_t>(nl_), 0);
        count_.assign(static_cast<std::size_t>(nl_), 0);
        vol_.assign(slab, 0.0);
        conc_.assign(slab * static_cast<std::size_t>(ns_ > 0 ? ns_ : 1), 0.0);
    }

    int count(int l) const { return count_[static_cast<std::size_t>(l)]; }
    int species() const { return ns_; }
    int capacity() const { return cap_; }

    double total_volume(int l) const {
        double v = 0.0;
        for (int i = 0; i < count(l); ++i)
            v += vol_[phys(l, i)];
        return v;
    }

    double seg_volume(int l, int i) const { return vol_[phys(l, i)]; }
    double seg_conc(int l, int i, int s) const { return conc_[cidx(l, i, s)]; }
    void set_seg_conc(int l, int i, int s, double c) { conc_[cidx(l, i, s)] = c; }

    void mean_conc(int l, double* out) const {
        for (int s = 0; s < ns_; ++s)
            out[s] = 0.0;
        double v = 0.0;
        for (int i = 0; i < count(l); ++i) {
            const double sv = vol_[phys(l, i)];
            v += sv;
            for (int s = 0; s < ns_; ++s)
                out[s] += sv * conc_[cidx(l, i, s)];
        }
        if (v > 0.0)
            for (int s = 0; s < ns_; ++s)
                out[s] /= v;
    }

    void clear_link(int l) {
        head_[static_cast<std::size_t>(l)] = 0;
        count_[static_cast<std::size_t>(l)] = 0;
    }

    double drain_back(int l, double vol_target, double* mass_out) {
        double drained = 0.0;
        while (vol_target - drained > 0.0 && count(l) > 0) {
            const int i = count(l) - 1; // back
            const std::size_t p = phys(l, i);
            const double take = std::min(vol_[p], vol_target - drained);
            for (int s = 0; s < ns_; ++s)
                mass_out[s] += take * conc_[cidx(l, i, s)];
            vol_[p] -= take;
            drained += take;
            if (vol_[p] <= 0.0)
                count_[static_cast<std::size_t>(l)] -= 1;
        }
        return drained;
    }

    double drain_front(int l, double vol_target, double* mass_out) {
        double drained = 0.0;
        while (vol_target - drained > 0.0 && count(l) > 0) {
            const std::size_t p = phys(l, 0);
            const double take = std::min(vol_[p], vol_target - drained);
            for (int s = 0; s < ns_; ++s)
                mass_out[s] += take * conc_[cidx(l, 0, s)];
            vol_[p] -= take;
            drained += take;
            if (vol_[p] <= 0.0) {
                head_[static_cast<std::size_t>(l)] =
                    (head_[static_cast<std::size_t>(l)] + 1) % cap_;
                count_[static_cast<std::size_t>(l)] -= 1;
            }
        }
        return drained;
    }

    void push_front(int l, double volume, const double* c) {
        if (volume <= 0.0)
            return;
        if (count(l) > 0 && mergeable(l, 0, c)) {
            merge_into(l, 0, volume, c);
            return;
        }
        if (count(l) >= cap_) {
            // D-L5: merge front pair, mass-conservative, then push.
            merge_front_pair(l);
        }
        head_[static_cast<std::size_t>(l)] = (head_[static_cast<std::size_t>(l)] - 1 + cap_) % cap_;
        count_[static_cast<std::size_t>(l)] += 1;
        const std::size_t p = phys(l, 0);
        vol_[p] = volume;
        for (int s = 0; s < ns_; ++s)
            conc_[cidx(l, 0, s)] = c[s];
    }

    void reverse(int l) {
        const int n = count(l);
        for (int i = 0, j = n - 1; i < j; ++i, --j) {
            std::swap(vol_[phys(l, i)], vol_[phys(l, j)]);
            for (int s = 0; s < ns_; ++s)
                std::swap(conc_[cidx(l, i, s)], conc_[cidx(l, j, s)]);
        }
    }

    void add_species(int s, double delta) {
        for (int l = 0; l < nl_; ++l)
            for (int i = 0; i < count(l); ++i)
                conc_[cidx(l, i, s)] += delta;
    }

    double decay_species(int s, double factor) {
        double removed = 0.0;
        for (int l = 0; l < nl_; ++l)
            for (int i = 0; i < count(l); ++i) {
                const std::size_t ci = cidx(l, i, s);
                removed += vol_[phys(l, i)] * conc_[ci] * (1.0 - factor);
                conc_[ci] *= factor;
            }
        return removed;
    }
 
private:
    std::size_t phys(int l, int i) const {
        const int slot = (head_[static_cast<std::size_t>(l)] + i) % cap_;
        return static_cast<std::size_t>(l) * static_cast<std::size_t>(cap_) +
               static_cast<std::size_t>(slot);
    }
    std::size_t cidx(int l, int i, int s) const {
        return static_cast<std::size_t>(s) * static_cast<std::size_t>(nl_) *
                   static_cast<std::size_t>(cap_) +
               phys(l, i);
    }

    bool mergeable(int l, int i, const double* c) const {
        for (int s = 0; s < ns_; ++s) {
            const double cf = conc_[cidx(l, i, s)];
            if (std::abs(c[s] - cf) > kMergeAtol + kMergeRtol * std::abs(cf))
                return false;
        }
        return true;
    }

    void merge_into(int l, int i, double volume, const double* c) {
        const std::size_t p = phys(l, i);
        const double v0 = vol_[p];
        const double vt = v0 + volume;
        if (vt <= 0.0)
            return;
        for (int s = 0; s < ns_; ++s) {
            const std::size_t ci = cidx(l, i, s);
            conc_[ci] = (conc_[ci] * v0 + c[s] * volume) / vt;
        }
        vol_[p] = vt;
    }

    void merge_front_pair(int l) {
        if (count(l) < 2)
            return;
        // Fold segment 1 into segment 0 (volume-weighted), drop slot 1 by
        // moving the head forward over it after copying 0 into 1.
        const std::size_t p0 = phys(l, 0);
        const std::size_t p1 = phys(l, 1);
        const double v0 = vol_[p0], v1 = vol_[p1], vt = v0 + v1;
        if (vt > 0.0)
            for (int s = 0; s < ns_; ++s) {
                const std::size_t c0 = cidx(l, 0, s);
                const std::size_t c1 = cidx(l, 1, s);
                conc_[c1] = (conc_[c0] * v0 + conc_[c1] * v1) / vt;
            }
        vol_[p1] = vt;
        head_[static_cast<std::size_t>(l)] = (head_[static_cast<std::size_t>(l)] + 1) % cap_;
        count_[static_cast<std::size_t>(l)] -= 1;
    }
 
    int nl_ = 0, ns_ = 0, cap_ = kMaxSegmentsPerLink;
    std::vector<int> head_, count_;
    std::vector<double> vol_;
    std::vector<double> conc_; 
};
 
} // namespace lard
} // namespace openswmm
 
#endif // OPENSWMM_QUALITY_LARD_SEGMENT_STORE_HPP