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

outram_foam_basic_lib/ode/
euler.rs
