ImplicitReachSolver.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: Apache-2.0
// Copyright 2026 Caleb Buahin

 
#pragma once
 
#include <array>
#include <vector>
 
#include "NetworkMeshData.hpp"
#include "FvKernels.hpp"
#include "NativeGmres.hpp"
 
namespace openswmm::fv {

class ImplicitReachSolver {
public:
    using Pair = std::array<double, 2>;
    using Block = std::array<double, 4>; // row major
    struct Problem {
        const NetworkMeshData* mesh = nullptr;
        std::vector<double> old_area, old_discharge;
        // NaN denotes a wall; finite values are depths over the flat bed (ft).
        double left_depth = 0.0, right_depth = 0.0;
        double lower_depth = 0.0, upper_depth = 0.0;
    };
    struct Endpoint {
        std::vector<double> area, discharge, depth;
        std::vector<kernels::FaceFlux> faces;
    };
    enum class LinearMethod { StructuredThenGmres, StructuredOnly, GmresOnly };
    enum class Preconditioner { BlockLine, BlockJacobi };
    struct BlockDiagnostics {
        int failed_block = -1;
        // Elimination pivot/row-scale ratio, NOT a condition-number estimate.
        double minimum_pivot_ratio = 1.0;
    };
    struct LinearControls {
        LinearMethod method = LinearMethod::StructuredThenGmres;
        Preconditioner preconditioner = Preconditioner::BlockLine;
        NativeGmres::Controls gmres;
    };
    struct LinearResult {
        bool converged = false, used_gmres = false, fallback = false;
        int identity_preconditioner_blocks = 0;
        int shifted_preconditioner_blocks = 0;
        double maximum_preconditioner_shift = 0.0;
        std::size_t wrapper_bytes = 0;
        BlockDiagnostics direct;
        NativeGmres::Result gmres;
        double relative_residual = 0.0;
    };
    struct Controls {
        int max_newton = 30, max_backtracks = 24;
        double absolute_volume = 1.e-11, relative_volume = 1.e-12;
        double absolute_momentum = 1.e-11, relative_momentum = 1.e-12;
        // Local finite-difference Jacobian for the first reference prototype;
        // all trial residuals re-evaluate the real closure and HLL selections.
        double difference_step = 1.e-4;
        LinearControls linear;
    };
    enum class Status {
        Converged,
        InvalidInput,
        NonFinite,
        LinearFailure,
        LineSearchFailed,
        IterationLimit
    };
    struct Result {
        Status status = Status::InvalidInput;
        int newton_iterations = 0, backtracks = 0;
        double scaled_residual = 0.0;
        double max_volume_residual = 0.0, max_momentum_residual = 0.0;
        double linear_relative_residual = 0.0;
        int gmres_iterations = 0, gmres_solves = 0, linear_fallbacks = 0;
        std::size_t max_gmres_workspace_bytes = 0;
        LinearResult last_linear;
        bool converged() const { return status == Status::Converged; }
    };
    struct Linearization {
        std::vector<Block> lower, diagonal, upper;
        std::vector<Pair> residual, row_scale;
        Pair column_scale{};
    };

    static bool residual(const Problem& p, const std::vector<Pair>& trial, double dt,
                         std::vector<Pair>& output,
                         std::vector<kernels::FaceFlux>* faces = nullptr);
    static bool linearize(const Problem& p, const std::vector<Pair>& trial, double dt,
                          double difference_step, Linearization& out);
    static bool solveBlocks(const Linearization& l, const std::vector<Pair>& rhs,
                            std::vector<Pair>& output, double& relative_residual,
                            BlockDiagnostics* diagnostics = nullptr);
    static LinearResult solveLinear(const Linearization& l, const std::vector<Pair>& rhs,
                                    const LinearControls& controls, std::vector<Pair>& output);
    Result solve(const Problem& p, double dt, const Controls& controls, Endpoint& output);
};
} // namespace openswmm::fv