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
22pub mod euler;
23pub mod rkf45;
24pub mod rosenbrock23;
25
26pub use euler::Euler;
27pub use rkf45::Rkf45;
28pub use rosenbrock23::Rosenbrock23;
29
30use crate::matrix::SquareMatrix;
31
32// ── ODE system trait ─────────────────────────────────────────────────────────
33
34/// Abstract ODE system `dy/dx = f(x, y)`. Maps to `Foam::ODESystem`.
35pub trait OdeSystem {
36    fn n_eqns(&self) -> usize;
37
38    /// Fill `dydx` with the derivatives at `(x, y)`.
39    fn derivatives(&self, x: f64, y: &[f64], dydx: &mut Vec<f64>);
40
41    /// Fill `dfdx` and `dfdy` with the Jacobian at `(x, y)`.
42    ///
43    /// Required only by stiff solvers (Rosenbrock23). The default panics.
44    fn jacobian(&self, _x: f64, _y: &[f64], _dfdx: &mut Vec<f64>, _dfdy: &mut SquareMatrix) {
45        unimplemented!("jacobian not implemented for this ODE system");
46    }
47}
48
49// ── Solver configuration ─────────────────────────────────────────────────────
50
51/// Common parameters for the adaptive step-size controller.
52/// Defaults match `Foam::adaptiveSolver` and `Foam::ODESolver`.
53#[derive(Debug, Clone)]
54pub struct OdeSolverConfig {
55    /// Absolute per-equation tolerance.
56    pub abs_tol: f64,
57    /// Relative per-equation tolerance.
58    pub rel_tol: f64,
59    /// Safety factor on the step-size scaling (0 < safeScale < 1).
60    pub safe_scale: f64,
61    /// Exponent for step *increase*.
62    pub alpha_inc: f64,
63    /// Exponent for step *decrease*.
64    pub alpha_dec: f64,
65    /// Minimum scale factor applied per step.
66    pub min_scale: f64,
67    /// Maximum scale factor applied per step.
68    pub max_scale: f64,
69    /// Maximum sub-steps for one `integrate()` call.
70    pub max_steps: usize,
71}
72
73impl Default for OdeSolverConfig {
74    fn default() -> Self {
75        Self {
76            abs_tol: 1e-6,
77            rel_tol: 1e-4,
78            safe_scale: 0.9,
79            alpha_inc: 0.2,
80            alpha_dec: 0.25,
81            min_scale: 0.2,
82            max_scale: 10.0,
83            max_steps: 10_000,
84        }
85    }
86}
87
88// ── Error type ───────────────────────────────────────────────────────────────
89
90#[derive(Debug, Clone, PartialEq)]
91pub enum OdeError {
92    StepSizeUnderflow,
93    MaxStepsExceeded(usize),
94}
95
96impl std::fmt::Display for OdeError {
97    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
98        match self {
99            Self::StepSizeUnderflow => write!(f, "ODE step size underflow"),
100            Self::MaxStepsExceeded(n) => write!(f, "ODE exceeded {n} steps"),
101        }
102    }
103}
104
105impl std::error::Error for OdeError {}
106
107// ── Shared utilities ─────────────────────────────────────────────────────────
108
109/// Normalised scalar error — max over all equations.
110/// Maps to `Foam::ODESolver::normalizeError`.
111pub(crate) fn normalize_error(
112    y0: &[f64],
113    y: &[f64],
114    err: &[f64],
115    abs_tol: f64,
116    rel_tol: f64,
117) -> f64 {
118    err.iter()
119        .zip(y0)
120        .zip(y)
121        .map(|((e, &a), &b)| {
122            let tol = abs_tol + rel_tol * a.abs().max(b.abs());
123            e.abs() / tol
124        })
125        .fold(0.0_f64, f64::max)
126}
127
128/// Adaptive step-size loop shared by all explicit solvers.
129///
130/// Calls `inner_step(x0, y0, dydx0, dx, y_out) -> err`, retrying with a
131/// smaller `dx` whenever `err > 1`. Updates `x`, `y`, and `dx_try`.
132/// Matches `Foam::adaptiveSolver::solve`.
133pub(crate) fn adaptive_step(
134    cfg: &OdeSolverConfig,
135    mut inner_step: impl FnMut(f64, &[f64], &[f64], f64, &mut Vec<f64>) -> f64,
136    ode: &dyn OdeSystem,
137    x: &mut f64,
138    y: &mut Vec<f64>,
139    dydx0: &mut Vec<f64>,
140    y_temp: &mut Vec<f64>,
141    dx_try: &mut f64,
142) -> Result<(), OdeError> {
143    let mut dx = *dx_try;
144    ode.derivatives(*x, y, dydx0);
145
146    let err = loop {
147        let err = inner_step(*x, y, dydx0, dx, y_temp);
148        if err <= 1.0 {
149            break err;
150        }
151        let scale = (cfg.safe_scale * err.powf(-cfg.alpha_dec)).max(cfg.min_scale);
152        dx *= scale;
153        if dx.abs() < f64::EPSILON {
154            return Err(OdeError::StepSizeUnderflow);
155        }
156    };
157
158    *x += dx;
159    std::mem::swap(y, y_temp);
160
161    let threshold = (cfg.max_scale / cfg.safe_scale).powf(-1.0 / cfg.alpha_inc);
162    *dx_try = if err > threshold {
163        let scale = (cfg.safe_scale * err.powf(-cfg.alpha_inc)).clamp(cfg.min_scale, cfg.max_scale);
164        dx * scale
165    } else {
166        dx * cfg.safe_scale * cfg.max_scale
167    };
168
169    Ok(())
170}
171
172/// Integrate from `x_start` to `x_end` using repeated adaptive steps.
173pub(crate) fn integrate_interval(
174    cfg: &OdeSolverConfig,
175    step_fn: &mut dyn FnMut(&mut f64, &mut Vec<f64>, &mut f64) -> Result<(), OdeError>,
176    x_start: f64,
177    x_end: f64,
178    y: &mut Vec<f64>,
179    dx_est: &mut f64,
180) -> Result<(), OdeError> {
181    let mut x = x_start;
182    let mut dx = *dx_est;
183    let mut n_steps = 0usize;
184
185    while x < x_end {
186        let dx_try = (x + dx).min(x_end) - x;
187        let mut dx_limited = dx_try;
188        step_fn(&mut x, y, &mut dx_limited)?;
189        dx = dx_limited;
190        n_steps += 1;
191        if n_steps > cfg.max_steps {
192            return Err(OdeError::MaxStepsExceeded(n_steps));
193        }
194    }
195
196    *dx_est = dx;
197    Ok(())
198}
outram_foam_basic_lib/ode/mod.rs

outram_foam_basic_lib/ode/
mod.rs
