SpeciesTransportKernels.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_TRANSPORT_FVKERNELS_SPECIES_TRANSPORT_KERNELS_HPP
#define OPENSWMM_ENGINE_TRANSPORT_FVKERNELS_SPECIES_TRANSPORT_KERNELS_HPP
 
#include <algorithm>
#include <cmath>
#include <vector>
 
#include "../../hydraulics/fv/FvKernels.hpp"
#include "../../hydraulics/fv/FvOptions.hpp"
#include "../../hydraulics/fv/NetworkMeshData.hpp"
 
namespace openswmm::transport::fvkernels {
 
using openswmm::fv::Limiter;
using openswmm::fv::NetworkMeshData;
using openswmm::fv::NetworkStateData;
using openswmm::fv::ScalarScheme;

inline double limitSlope(double a, double b, Limiter lim) {
    if (a * b <= 0.0) return 0.0;                  // extremum ⇒ no slope
    switch (lim) {
        case Limiter::VANLEER:
            return 2.0 * a * b / (a + b);
        case Limiter::SUPERBEE: {
            const double s = (a > 0.0) ? 1.0 : -1.0;
            return s * std::max(std::min(2.0 * std::fabs(a), std::fabs(b)),
                                std::min(std::fabs(a), 2.0 * std::fabs(b)));
        }
        case Limiter::MINMOD:
        default:
            return (std::fabs(a) < std::fabs(b)) ? a : b;
    }
}

struct SpeciesKernelView {
    const NetworkMeshData* mesh = nullptr;
    NetworkStateData*      state = nullptr;
 
    // Scheme options (subset of FvOptions the kernels consume).
    ScalarScheme scalar_scheme = ScalarScheme::MUSCL;
    Limiter      limiter       = Limiter::MINMOD;
    double       dispersion    = 0.0;
    bool         hllc          = true;

    const std::vector<double>* cell_dispersion = nullptr;
 
    // Hydrodynamic face records of the current substep (read-only).
    const std::vector<double>* f_mass  = nullptr;  
    const std::vector<double>* f_sstar = nullptr;  
    const std::vector<openswmm::fv::kernels::FaceState>* f_state_l = nullptr;
    const std::vector<openswmm::fv::kernels::FaceState>* f_state_r = nullptr;
    const std::vector<openswmm::fv::kernels::FaceFlux>*  f_flux    = nullptr;
 
    // Derived cell fields (read-only).
    const std::vector<double>* cell_u      = nullptr;  
    const std::vector<char>*   cell_active = nullptr;
    const std::vector<int>*    active_faces = nullptr;
 
    // Caller-owned outputs / scratch.
    std::vector<double>* f_phi_l    = nullptr;  
    std::vector<double>* f_phi_r    = nullptr;  
    std::vector<double>* f_phi_flux = nullptr;  
    std::vector<double>* cell_slope = nullptr;
    std::vector<double>* lo_flux    = nullptr;
    std::vector<double>* anti_flux  = nullptr;
    std::vector<double>* td         = nullptr;
    std::vector<double>* anew       = nullptr;
    std::vector<double>* rplus      = nullptr;
    std::vector<double>* rminus     = nullptr;
};

void reconstructScalars(const SpeciesKernelView& v, double dt);

void limitSpeciesFluxes(const SpeciesKernelView& v, int species, double dt);

void dispersionSolve(const SpeciesKernelView& v, double dt);
 
}  // namespace openswmm::transport::fvkernels
 
#endif  // OPENSWMM_ENGINE_TRANSPORT_FVKERNELS_SPECIES_TRANSPORT_KERNELS_HPP