outram_foam_basic_lib/ode/
euler.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
22use super::{adaptive_step, integrate_interval, normalize_error, OdeError, OdeSolverConfig, OdeSystem};
23
24/// Explicit first-order Euler solver with adaptive step size.
25/// Maps to `Foam::Euler` (which inherits from `adaptiveSolver`).
26#[derive(Debug, Clone)]
27pub struct Euler {
28    /// Adaptive step-size controller settings (tolerances, scale limits).
29    pub config: OdeSolverConfig,
30    // scratch buffers (sized to n_eqns)
31    dydx0: Vec<f64>,
32    y_temp: Vec<f64>,
33    err: Vec<f64>,
34}
35
36impl Euler {
37    /// Create a solver for an `n`-equation system with the given absolute and
38    /// relative per-equation error tolerances; other controller parameters take
39    /// their [`OdeSolverConfig::default`] values.
40    pub fn new(n: usize, abs_tol: f64, rel_tol: f64) -> Self {
41        let mut cfg = OdeSolverConfig::default();
42        cfg.abs_tol = abs_tol;
43        cfg.rel_tol = rel_tol;
44        Self {
45            config: cfg,
46            dydx0: vec![0.0; n],
47            y_temp: vec![0.0; n],
48            err: vec![0.0; n],
49        }
50    }
51
52    /// Take one adaptive step. On return `x` and `y` are updated and
53    /// `dx_try` holds a suggested step size for the next call.
54    pub fn solve_step<Sys: OdeSystem + ?Sized>(
55        &mut self,
56        ode: &Sys,
57        x: &mut f64,
58        y: &mut Vec<f64>,
59        dx_try: &mut f64,
60    ) -> Result<(), OdeError> {
61        let cfg = self.config.clone();
62        let err_buf = &mut self.err;
63        let abs_tol = cfg.abs_tol;
64        let rel_tol = cfg.rel_tol;
65
66        adaptive_step(
67            &cfg,
68            |_x0, y0, dydx0, dx, y_out| {
69                for i in 0..y0.len() {
70                    err_buf[i] = dx * dydx0[i];
71                    y_out[i] = y0[i] + err_buf[i];
72                }
73                normalize_error(y0, y_out, err_buf, abs_tol, rel_tol)
74            },
75            ode,
76            x,
77            y,
78            &mut self.dydx0,
79            &mut self.y_temp,
80            dx_try,
81        )
82    }
83
84    /// Integrate from `x_start` to `x_end`.
85    pub fn integrate<Sys: OdeSystem + ?Sized>(
86        &mut self,
87        ode: &Sys,
88        x_start: f64,
89        x_end: f64,
90        y: &mut Vec<f64>,
91        dx_est: &mut f64,
92    ) -> Result<(), OdeError> {
93        let cfg = self.config.clone();
94        integrate_interval(
95            &cfg,
96            &mut |x, y, dx| self.solve_step(ode, x, y, dx),
97            x_start,
98            x_end,
99            y,
100            dx_est,
101        )
102    }
103}
104
105#[cfg(test)]
106mod tests {
107    use petir::mathf::RealMath;
108    use super::*;
109
110    struct DecayOde;
111    impl OdeSystem for DecayOde {
112        fn n_eqns(&self) -> usize {
113            1
114        }
115        fn derivatives(&self, _x: f64, y: &[f64], dydx: &mut Vec<f64>) {
116            dydx[0] = -y[0];
117        }
118    }
119
120    struct RotationOde;
121    impl OdeSystem for RotationOde {
122        fn n_eqns(&self) -> usize {
123            2
124        }
125        fn derivatives(&self, _x: f64, y: &[f64], dydx: &mut Vec<f64>) {
126            dydx[0] = -y[1];
127            dydx[1] = y[0];
128        }
129    }
130
131    #[test]
132    fn euler_exponential_decay() {
133        // y' = -y, y(0) = 1  →  y(1) = e^{-1} ≈ 0.36788
134        // Euler is 1st-order: use loose tolerance so adaptive step stays ~0.01 (≈100 steps)
135        let ode = DecayOde;
136        let mut solver = Euler::new(1, 1e-3, 1e-2);
137        let mut y = vec![1.0_f64];
138        let mut dx = 0.1;
139        solver.integrate(&ode, 0.0, 1.0, &mut y, &mut dx).unwrap();
140        let expected = (-1.0_f64).r_exp();
141        assert!(
142            (y[0] - expected).abs() < 1e-2,
143            "y={}, expected {}",
144            y[0],
145            expected
146        );
147    }
148
149    #[test]
150    fn euler_rotation() {
151        // y1' = -y2, y2' = y1,  y(0) = (1, 0)  →  y(π/2) ≈ (0, 1)
152        let ode = RotationOde;
153        let mut solver = Euler::new(2, 1e-3, 1e-2);
154        let mut y = vec![1.0_f64, 0.0];
155        let mut dx = 0.1;
156        solver
157            .integrate(&ode, 0.0, std::f64::consts::PI / 2.0, &mut y, &mut dx)
158            .unwrap();
159        assert!(y[0].abs() < 5e-2, "y0={}", y[0]);
160        assert!((y[1] - 1.0).abs() < 5e-2, "y1={}", y[1]);
161    }
162}
outram_foam_basic_lib/ode/euler.rs

outram_foam_basic_lib/ode/
euler.rs
