outram_foam_basic_lib/ode/
mod.rs

1// SPDX-License-Identifier: GPL-3.0-only
2// Copyright (C) 2026 OUTRAM PARK contributors
3// Derived from OpenFOAM (www.openfoam.com)
4// Copyright (C) 2004-2023 OpenFOAM Foundation
5// Copyright (C) 2016-2023 OpenCFD Ltd.
6//
7// This file is part of OUTRAM PARK.
8//
9// OUTRAM PARK is free software: you can redistribute it and/or modify it
10// under the terms of the GNU General Public License as published by the
11// Free Software Foundation, either version 3 of the License, or (at your
12// option) any later version.
13//
14// OUTRAM PARK is distributed in the hope that it will be useful, but
15// WITHOUT ANY WARRANTY; without even the implied warranty of
16// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17// General Public License for more details.
18//
19// You should have received a copy of the GNU General Public License along
20// with OUTRAM PARK.  If not, see <https://www.gnu.org/licenses/>.
21
22//! Ordinary differential equation solvers for systems `dy/dx = f(x, y)`.
23//!
24//! Ports the OpenFOAM `ODE` layer: user systems implement the [`OdeSystem`]
25//! trait, and one of the concrete steppers integrates them with adaptive step
26//! control — [`Euler`] (explicit 1st order), [`Rkf45`] (explicit Runge-Kutta-
27//! Fehlberg 4(5)), and [`Rosenbrock23`] (semi-implicit, for stiff systems,
28//! requiring a Jacobian). The independent variable `x`, state `y`, and step
29//! size are bare `f64` in the caller's own units; tolerances are set through
30//! [`OdeSolverConfig`].
31
32pub mod euler;
33pub mod rkf45;
34pub mod rosenbrock23;
35
36pub use euler::Euler;
37pub use rkf45::Rkf45;
38pub use rosenbrock23::Rosenbrock23;
39
40use crate::matrix::SquareMatrix;
41
42// ── ODE system trait ─────────────────────────────────────────────────────────
43
44/// Abstract ODE system `dy/dx = f(x, y)`. Maps to `Foam::ODESystem`.
45pub trait OdeSystem {
46    /// Number of coupled equations (the length of the state vector `y`).
47    fn n_eqns(&self) -> usize;
48
49    /// Fill `dydx` with the derivatives at `(x, y)`.
50    fn derivatives(&self, x: f64, y: &[f64], dydx: &mut Vec<f64>);
51
52    /// Fill `dfdx` and `dfdy` with the Jacobian at `(x, y)`.
53    ///
54    /// Required only by stiff solvers (Rosenbrock23). The default panics.
55    fn jacobian(&self, _x: f64, _y: &[f64], _dfdx: &mut Vec<f64>, _dfdy: &mut SquareMatrix) {
56        unimplemented!("jacobian not implemented for this ODE system");
57    }
58}
59
60// ── Solver configuration ─────────────────────────────────────────────────────
61
62/// Common parameters for the adaptive step-size controller.
63/// Defaults match `Foam::adaptiveSolver` and `Foam::ODESolver`.
64#[derive(Debug, Clone)]
65pub struct OdeSolverConfig {
66    /// Absolute per-equation tolerance.
67    pub abs_tol: f64,
68    /// Relative per-equation tolerance.
69    pub rel_tol: f64,
70    /// Safety factor on the step-size scaling (0 < safeScale < 1).
71    pub safe_scale: f64,
72    /// Exponent for step *increase*.
73    pub alpha_inc: f64,
74    /// Exponent for step *decrease*.
75    pub alpha_dec: f64,
76    /// Minimum scale factor applied per step.
77    pub min_scale: f64,
78    /// Maximum scale factor applied per step.
79    pub max_scale: f64,
80    /// Maximum sub-steps for one `integrate()` call.
81    pub max_steps: usize,
82}
83
84impl Default for OdeSolverConfig {
85    fn default() -> Self {
86        Self {
87            abs_tol: 1e-6,
88            rel_tol: 1e-4,
89            safe_scale: 0.9,
90            alpha_inc: 0.2,
91            alpha_dec: 0.25,
92            min_scale: 0.2,
93            max_scale: 10.0,
94            max_steps: 10_000,
95        }
96    }
97}
98
99// ── Error type ───────────────────────────────────────────────────────────────
100
101/// Failure modes of an adaptive integration.
102#[derive(Debug, Clone, PartialEq)]
103pub enum OdeError {
104    /// The step size shrank below machine epsilon while trying to meet the
105    /// error tolerance — the system is too stiff for the chosen solver, or the
106    /// tolerances are unattainable.
107    StepSizeUnderflow,
108    /// The interval could not be spanned within `max_steps` sub-steps; carries
109    /// the number of steps taken.
110    MaxStepsExceeded(usize),
111}
112
113impl std::fmt::Display for OdeError {
114    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
115        match self {
116            Self::StepSizeUnderflow => write!(f, "ODE step size underflow"),
117            Self::MaxStepsExceeded(n) => write!(f, "ODE exceeded {n} steps"),
118        }
119    }
120}
121
122impl std::error::Error for OdeError {}
123
124// ── Shared utilities ─────────────────────────────────────────────────────────
125
126/// Normalised scalar error — max over all equations.
127/// Maps to `Foam::ODESolver::normalizeError`.
128pub(crate) fn normalize_error(
129    y0: &[f64],
130    y: &[f64],
131    err: &[f64],
132    abs_tol: f64,
133    rel_tol: f64,
134) -> f64 {
135    err.iter()
136        .zip(y0)
137        .zip(y)
138        .map(|((e, &a), &b)| {
139            let tol = abs_tol + rel_tol * a.abs().max(b.abs());
140            e.abs() / tol
141        })
142        .fold(0.0_f64, f64::max)
143}
144
145/// Adaptive step-size loop shared by all explicit solvers.
146///
147/// Calls `inner_step(x0, y0, dydx0, dx, y_out) -> err`, retrying with a
148/// smaller `dx` whenever `err > 1`. Updates `x`, `y`, and `dx_try`.
149/// Matches `Foam::adaptiveSolver::solve`.
150pub(crate) fn adaptive_step(
151    cfg: &OdeSolverConfig,
152    mut inner_step: impl FnMut(f64, &[f64], &[f64], f64, &mut Vec<f64>) -> f64,
153    ode: &dyn OdeSystem,
154    x: &mut f64,
155    y: &mut Vec<f64>,
156    dydx0: &mut Vec<f64>,
157    y_temp: &mut Vec<f64>,
158    dx_try: &mut f64,
159) -> Result<(), OdeError> {
160    let mut dx = *dx_try;
161    ode.derivatives(*x, y, dydx0);
162
163    let err = loop {
164        let err = inner_step(*x, y, dydx0, dx, y_temp);
165        if err <= 1.0 {
166            break err;
167        }
168        let scale = (cfg.safe_scale * err.powf(-cfg.alpha_dec)).max(cfg.min_scale);
169        dx *= scale;
170        if dx.abs() < f64::EPSILON {
171            return Err(OdeError::StepSizeUnderflow);
172        }
173    };
174
175    *x += dx;
176    std::mem::swap(y, y_temp);
177
178    let threshold = (cfg.max_scale / cfg.safe_scale).powf(-1.0 / cfg.alpha_inc);
179    *dx_try = if err > threshold {
180        let scale = (cfg.safe_scale * err.powf(-cfg.alpha_inc)).clamp(cfg.min_scale, cfg.max_scale);
181        dx * scale
182    } else {
183        dx * cfg.safe_scale * cfg.max_scale
184    };
185
186    Ok(())
187}
188
189/// Integrate from `x_start` to `x_end` using repeated adaptive steps.
190pub(crate) fn integrate_interval(
191    cfg: &OdeSolverConfig,
192    step_fn: &mut dyn FnMut(&mut f64, &mut Vec<f64>, &mut f64) -> Result<(), OdeError>,
193    x_start: f64,
194    x_end: f64,
195    y: &mut Vec<f64>,
196    dx_est: &mut f64,
197) -> Result<(), OdeError> {
198    let mut x = x_start;
199    let mut dx = *dx_est;
200    let mut n_steps = 0usize;
201
202    while x < x_end {
203        let dx_try = (x + dx).min(x_end) - x;
204        let mut dx_limited = dx_try;
205        step_fn(&mut x, y, &mut dx_limited)?;
206        dx = dx_limited;
207        n_steps += 1;
208        if n_steps > cfg.max_steps {
209            return Err(OdeError::MaxStepsExceeded(n_steps));
210        }
211    }
212
213    *dx_est = dx;
214    Ok(())
215}
outram_foam_basic_lib/ode/mod.rs

outram_foam_basic_lib/ode/
mod.rs
