OdeSolver.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_ODE_SOLVER_HPP
#define OPENSWMM_ODE_SOLVER_HPP
 
#include <functional>
#include <vector>
 
namespace openswmm {
namespace ode {
 
// ============================================================================
// Constants (matching legacy)
// ============================================================================
 
constexpr int MAXSTP = 10000;      
constexpr double TINY = 1.0e-30;   
constexpr double SAFETY = 0.9;     
constexpr double PGROW = -0.2;     
constexpr double PSHRNK = -0.25;   
constexpr double ERRCON = 1.89e-4; 
 
// ============================================================================
// Return codes
// ============================================================================
 
constexpr int ODE_OK = 0;        
constexpr int ODE_TOO_MANY = 1;  
constexpr int ODE_UNDERFLOW = 2; 
constexpr int ODE_MAX_STEPS = 3; 
 
// ============================================================================
// Derivative callback
// ============================================================================

using DerivFunc = std::function<void(double x, const double* y, double* dydx)>;
 
// ============================================================================
// Per-element integrator
// ============================================================================

int integrate(double* y, int n, double x1, double x2, double eps, double h1,
              const DerivFunc& derivs);

template <typename F>
int integrateCallable(double* y, int n, double x1, double x2, double eps, double h1, F&& derivs);
 
// ============================================================================
// Batch integrator (SoA, for vectorized subcatchment processing)
// ============================================================================

using BatchDerivFunc = std::function<void(double x, const double* y0, const double* y1, double* dy0,
                                          double* dy1, int n_sys)>;

int integrate_batch_2eq(double* y0, double* y1, int n_sys, double x1, double x2, double eps,
                        double h1, const BatchDerivFunc& derivs);
 
} // namespace ode
} // namespace openswmm
 
#include "OdeSolverDetail.hpp"
 
#endif // OPENSWMM_ODE_SOLVER_HPP