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    /// Adaptive step-size controller settings (tolerances, scale limits).
30    pub config: OdeSolverConfig,
31    // scratch buffers (sized to n_eqns)
32    dydx0: Vec<f64>,
33    y_temp: Vec<f64>,
34    err: Vec<f64>,
35}
36
37impl Euler {
38    /// Create a solver for an `n`-equation system with the given absolute and
39    /// relative per-equation error tolerances; other controller parameters take
40    /// their [`OdeSolverConfig::default`] values.
41    pub fn new(n: usize, abs_tol: f64, rel_tol: f64) -> Self {
42        let mut cfg = OdeSolverConfig::default();
43        cfg.abs_tol = abs_tol;
44        cfg.rel_tol = rel_tol;
45        Self {
46            config: cfg,
47            dydx0: vec![0.0; n],
48            y_temp: vec![0.0; n],
49            err: vec![0.0; n],
50        }
51    }
52
53    /// Take one adaptive step. On return `x` and `y` are updated and
54    /// `dx_try` holds a suggested step size for the next call.
55    pub fn solve_step(
56        &mut self,
57        ode: &dyn OdeSystem,
58        x: &mut f64,
59        y: &mut Vec<f64>,
60        dx_try: &mut f64,
61    ) -> Result<(), OdeError> {
62        let cfg = self.config.clone();
63        let err_buf = &mut self.err;
64        let abs_tol = cfg.abs_tol;
65        let rel_tol = cfg.rel_tol;
66
67        adaptive_step(
68            &cfg,
69            |_x0, y0, dydx0, dx, y_out| {
70                for i in 0..y0.len() {
71                    err_buf[i] = dx * dydx0[i];
72                    y_out[i] = y0[i] + err_buf[i];
73                }
74                normalize_error(y0, y_out, err_buf, abs_tol, rel_tol)
75            },
76            ode,
77            x,
78            y,
79            &mut self.dydx0,
80            &mut self.y_temp,
81            dx_try,
82        )
83    }
84
85    /// Integrate from `x_start` to `x_end`.
86    pub fn integrate(
87        &mut self,
88        ode: &dyn OdeSystem,
89        x_start: f64,
90        x_end: f64,
91        y: &mut Vec<f64>,
92        dx_est: &mut f64,
93    ) -> Result<(), OdeError> {
94        let cfg = self.config.clone();
95        integrate_interval(
96            &cfg,
97            &mut |x, y, dx| self.solve_step(ode, x, y, dx),
98            x_start,
99            x_end,
100            y,
101            dx_est,
102        )
103    }
104}
105
106#[cfg(test)]
107mod tests {
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).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
