CvodeSurfaceSolver.hpp

Go to the documentation of this file.

 
#ifndef OPENSWMM_ENGINE_2D_CVODE_SURFACE_SOLVER_HPP
#define OPENSWMM_ENGINE_2D_CVODE_SURFACE_SOLVER_HPP
 
#include "../data/MeshData.hpp"
#include "../data/SurfaceStateData.hpp"
#include "../data/SolverOptions2D.hpp"
#include "ISurfaceSolver.hpp"
 
#ifdef OPENSWMM_HAS_2D
 
#include <vector>
#include <memory>
 
// Forward declarations for SUNDIALS types (avoid pulling in full headers)
struct SUNContext_;
typedef struct SUNContext_* SUNContext;
typedef struct _generic_N_Vector* N_Vector;
typedef struct _generic_SUNLinearSolver* SUNLinearSolver;
 
namespace openswmm::twoD {
 
class CvodeSurfaceSolver;
#if defined(OPENSWMM_HAVE_HYPRE)
class HypreAmgPreconditioner;  // complete type only in the .cpp (guarded)
#endif

struct CvodeSolverContext {
    MeshData*           mesh   = nullptr;
    SurfaceStateData*   state  = nullptr;
    SolverOptions2D*    opts   = nullptr;
    CvodeSurfaceSolver* solver = nullptr;
    bool                amg_active = false;  
};

class CvodeSurfaceSolver : public ISurfaceSolver {
public:
    // Defined out-of-line in the .cpp (where HypreAmgPreconditioner is a
    // complete type) so the unique_ptr member's destructor can be instantiated
    // there rather than at every construction site.
    CvodeSurfaceSolver();
    ~CvodeSurfaceSolver() override;
 
    // Non-copyable
    CvodeSurfaceSolver(const CvodeSurfaceSolver&) = delete;
    CvodeSurfaceSolver& operator=(const CvodeSurfaceSolver&) = delete;
 
    // Movable
    CvodeSurfaceSolver(CvodeSurfaceSolver&& o) noexcept;
    CvodeSurfaceSolver& operator=(CvodeSurfaceSolver&& o) noexcept;

    void initialize(MeshData& mesh, SurfaceStateData& state,
                    SolverOptions2D& opts) override;

    double advance(double t_current, double t_target) override;

    void reinitialize(double t0) override;

    void finalize() override;

    long last_num_steps() const noexcept override { return last_nsteps_; }

    double last_step_size() const noexcept override { return last_h_; }

    const std::vector<double>& last_coupling_exchange() const noexcept override {
        return last_coupling_exchange_;
    }

    bool is_initialized() const noexcept override { return cvode_mem_ != nullptr; }
 
private:
    void*           cvode_mem_ = nullptr;  
    SUNLinearSolver ls_        = nullptr;  
    N_Vector        y_         = nullptr;  
    N_Vector        abstol_    = nullptr;  
    SUNContext      sun_ctx_   = nullptr;  
 
    CvodeSolverContext ctx_;               
 
    long   last_nsteps_ = 0;
    double last_h_      = 0.0;
 
    // ── Live node-coupling (macro-step) state augmentation ────────────────────
    // When state.node_coupling is set (COUPLING_INTERVAL > 1), the state vector
    // is augmented to nt + nc_: the extra nc_ entries A_k = ∫Q_k dt accumulate
    // the live orifice exchange per coupling point so the 1D↔2D booking is
    // conservative (CVODE integrates A_k with the same BDF method as V). nc_ == 0
    // (default) ⇒ no augmentation, byte-identical to the legacy held-flux path.
    int    nc_ = 0;                              
    std::vector<double> coupling_accum_start_;   
    std::vector<double> last_coupling_exchange_; 

    std::vector<double> precond_diag_;
 
#if defined(OPENSWMM_HAVE_HYPRE)
    std::unique_ptr<HypreAmgPreconditioner> amg_precond_;
    bool amg_used_last_setup_ = false;
#endif
 
    // ------------------------------------------------------------------
    // SUNDIALS callbacks. All three have C linkage requirements imposed
    // by SUNDIALS' function-pointer typedefs; we expose them as static
    // member functions to keep them inside the class scope.
    // ------------------------------------------------------------------

    static int rhs_fn(double t, N_Vector y, N_Vector ydot, void* user_data);

    static int psetup_fn(double t, N_Vector y, N_Vector fy,
                          int jok, int* jcurPtr, double gamma,
                          void* user_data);

    static int psolve_fn(double t, N_Vector y, N_Vector fy,
                          N_Vector r, N_Vector z,
                          double gamma, double delta, int lr,
                          void* user_data);
};
 
} // namespace openswmm::twoD
 
#endif // OPENSWMM_HAS_2D
 
#endif // OPENSWMM_ENGINE_2D_CVODE_SURFACE_SOLVER_HPP