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