outram_foam_basic_lib/ode/

rkf45.rs

// SPDX-License-Identifier: GPL-3.0-only
// Copyright (C) 2026 OUTRAM PARK contributors
// Derived from OpenFOAM (www.openfoam.com)
// Copyright (C) 2004-2023 OpenFOAM Foundation
// Copyright (C) 2016-2023 OpenCFD Ltd.
//
// This file is part of OUTRAM PARK.
//
// OUTRAM PARK is free software: you can redistribute it and/or modify it
// under the terms of the GNU General Public License as published by the
// Free Software Foundation, either version 3 of the License, or (at your
// option) any later version.
//
// OUTRAM PARK is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with OUTRAM PARK.  If not, see <https://www.gnu.org/licenses/>.

use super::{adaptive_step, integrate_interval, normalize_error, OdeError, OdeSolverConfig, OdeSystem};

// Butcher tableau coefficients — Runge-Kutta-Fehlberg 4(5)
// Source: Foam::RKF45 constants
const C2: f64 = 1.0 / 4.0;
const C3: f64 = 3.0 / 8.0;
const C4: f64 = 12.0 / 13.0;
const C5: f64 = 1.0;
const C6: f64 = 1.0 / 2.0;

const A21: f64 = 1.0 / 4.0;
const A31: f64 = 3.0 / 32.0;
const A32: f64 = 9.0 / 32.0;
const A41: f64 = 1932.0 / 2197.0;
const A42: f64 = -7200.0 / 2197.0;
const A43: f64 = 7296.0 / 2197.0;
const A51: f64 = 439.0 / 216.0;
const A52: f64 = -8.0;
const A53: f64 = 3680.0 / 513.0;
const A54: f64 = -845.0 / 4104.0;
const A61: f64 = -8.0 / 27.0;
const A62: f64 = 2.0;
const A63: f64 = -3544.0 / 2565.0;
const A64: f64 = 1859.0 / 4104.0;
const A65: f64 = -11.0 / 40.0;

// 5th-order weights
const B1: f64 = 16.0 / 135.0;
const B3: f64 = 6656.0 / 12825.0;
const B4: f64 = 28561.0 / 56430.0;
const B5: f64 = -9.0 / 50.0;
const B6: f64 = 2.0 / 55.0;

// Error = 5th-order − 4th-order weights
const E1: f64 = 25.0 / 216.0 - B1;
const E3: f64 = 1408.0 / 2565.0 - B3;
const E4: f64 = 2197.0 / 4104.0 - B4;
const E5: f64 = -1.0 / 5.0 - B5;
const E6: f64 = -B6;

/// Runge-Kutta-Fehlberg 4(5) explicit solver with adaptive step size.
/// Maps to `Foam::RKF45`.
#[derive(Debug, Clone)]
pub struct Rkf45 {
    /// Adaptive step-size controller settings (tolerances, scale limits).
    pub config: OdeSolverConfig,
    dydx0: Vec<f64>,
    y_temp: Vec<f64>,
    k2: Vec<f64>,
    k3: Vec<f64>,
    k4: Vec<f64>,
    k5: Vec<f64>,
    k6: Vec<f64>,
    err: Vec<f64>,
    y_stage: Vec<f64>,
}

impl Rkf45 {
    /// Create a solver for an `n`-equation system with the given absolute and
    /// relative per-equation error tolerances; other controller parameters take
    /// their [`OdeSolverConfig::default`] values.
    pub fn new(n: usize, abs_tol: f64, rel_tol: f64) -> Self {
        let mut cfg = OdeSolverConfig::default();
        cfg.abs_tol = abs_tol;
        cfg.rel_tol = rel_tol;
        Self {
            config: cfg,
            dydx0: vec![0.0; n],
            y_temp: vec![0.0; n],
            k2: vec![0.0; n],
            k3: vec![0.0; n],
            k4: vec![0.0; n],
            k5: vec![0.0; n],
            k6: vec![0.0; n],
            err: vec![0.0; n],
            y_stage: vec![0.0; n],
        }
    }

    /// Take one adaptive step. On return `x` and `y` are updated and `dx_try`
    /// holds a suggested step size for the next call.
    pub fn solve_step<Sys: OdeSystem + ?Sized>(
        &mut self,
        ode: &Sys,
        x: &mut f64,
        y: &mut Vec<f64>,
        dx_try: &mut f64,
    ) -> Result<(), OdeError> {
        let cfg = self.config.clone();
        let abs_tol = cfg.abs_tol;
        let rel_tol = cfg.rel_tol;

        let k2 = &mut self.k2;
        let k3 = &mut self.k3;
        let k4 = &mut self.k4;
        let k5 = &mut self.k5;
        let k6 = &mut self.k6;
        let err = &mut self.err;
        let y_stage = &mut self.y_stage;

        adaptive_step(
            &cfg,
            |x0, y0, dydx0, dx, y_out| {
                let n = y0.len();

                // Stage 2
                for i in 0..n {
                    y_stage[i] = y0[i] + A21 * dx * dydx0[i];
                }
                ode.derivatives(x0 + C2 * dx, y_stage, k2);

                // Stage 3
                for i in 0..n {
                    y_stage[i] = y0[i] + dx * (A31 * dydx0[i] + A32 * k2[i]);
                }
                ode.derivatives(x0 + C3 * dx, y_stage, k3);

                // Stage 4
                for i in 0..n {
                    y_stage[i] = y0[i] + dx * (A41 * dydx0[i] + A42 * k2[i] + A43 * k3[i]);
                }
                ode.derivatives(x0 + C4 * dx, y_stage, k4);

                // Stage 5
                for i in 0..n {
                    y_stage[i] =
                        y0[i] + dx * (A51 * dydx0[i] + A52 * k2[i] + A53 * k3[i] + A54 * k4[i]);
                }
                ode.derivatives(x0 + C5 * dx, y_stage, k5);

                // Stage 6
                for i in 0..n {
                    y_stage[i] = y0[i]
                        + dx * (A61 * dydx0[i]
                            + A62 * k2[i]
                            + A63 * k3[i]
                            + A64 * k4[i]
                            + A65 * k5[i]);
                }
                ode.derivatives(x0 + C6 * dx, y_stage, k6);

                // 5th-order solution
                for i in 0..n {
                    y_out[i] = y0[i]
                        + dx * (B1 * dydx0[i] + B3 * k3[i] + B4 * k4[i] + B5 * k5[i] + B6 * k6[i]);
                }

                // Error = 5th − 4th
                for i in 0..n {
                    err[i] =
                        dx * (E1 * dydx0[i] + E3 * k3[i] + E4 * k4[i] + E5 * k5[i] + E6 * k6[i]);
                }

                normalize_error(y0, y_out, err, abs_tol, rel_tol)
            },
            ode,
            x,
            y,
            &mut self.dydx0,
            &mut self.y_temp,
            dx_try,
        )
    }

    /// Integrate from `x_start` to `x_end`, updating `y` in place and leaving
    /// the last accepted step size in `dx_est`.
    pub fn integrate<Sys: OdeSystem + ?Sized>(
        &mut self,
        ode: &Sys,
        x_start: f64,
        x_end: f64,
        y: &mut Vec<f64>,
        dx_est: &mut f64,
    ) -> Result<(), OdeError> {
        let cfg = self.config.clone();
        integrate_interval(
            &cfg,
            &mut |x, y, dx| self.solve_step(ode, x, y, dx),
            x_start,
            x_end,
            y,
            dx_est,
        )
    }
}

#[cfg(test)]
mod tests {
    use petir::mathf::RealMath;
    use super::*;

    struct DecayOde;
    impl OdeSystem for DecayOde {
        fn n_eqns(&self) -> usize {
            1
        }
        fn derivatives(&self, _x: f64, y: &[f64], dydx: &mut Vec<f64>) {
            dydx[0] = -y[0];
        }
    }

    struct VanDerPol {
        mu: f64,
    }
    impl OdeSystem for VanDerPol {
        fn n_eqns(&self) -> usize {
            2
        }
        fn derivatives(&self, _x: f64, y: &[f64], dydx: &mut Vec<f64>) {
            dydx[0] = y[1];
            dydx[1] = self.mu * (1.0 - y[0] * y[0]) * y[1] - y[0];
        }
    }

    /// `y' = -y`, but the derivative is NaN once `|y| > 2`: finite on the
    /// whole true trajectory, non-finite only at the overshot stage states an
    /// oversized trial step visits.
    struct DecayNanOutsideRange;
    impl OdeSystem for DecayNanOutsideRange {
        fn n_eqns(&self) -> usize {
            1
        }
        fn derivatives(&self, _x: f64, y: &[f64], dydx: &mut Vec<f64>) {
            dydx[0] = if y[0].abs() > 2.0 { f64::NAN } else { -y[0] };
        }
    }

    /// **Methodology.** A trial step that overflows while the accepted state is
    /// finite must be rejected and retried with a smaller step, as upstream
    /// `Foam::adaptiveSolver::solve` does (`err = inf` -> `scale = minScale`).
    /// Starting `y(0) = 1` with a first step `dx = 10` sends RKF45's later
    /// stages beyond `|y| = 2`, where the derivative is NaN. Pass criterion:
    /// `Ok`, and `y(1)` within `1e-6` of `e^{-1}`. Before 2026-10-21 (gh:#985)
    /// the op-ad6h guard returned `Err(NonFiniteState)` on the first trial,
    /// which is what broke offbeat's VISCOCHAB relaxation test.
    ///
    /// **Result (2026-10-21).** Passes; the solver recovers and the error is
    /// below `1e-6`.
    #[test]
    fn rkf45_recovers_from_a_trial_step_that_overflows() {
        let mut solver = Rkf45::new(1, 1e-8, 1e-6);
        let mut y = vec![1.0_f64];
        let mut dx = 10.0;
        solver
            .integrate(&DecayNanOutsideRange, 0.0, 1.0, &mut y, &mut dx)
            .unwrap();
        assert!((y[0] - (-1.0_f64).r_exp()).abs() < 1e-6, "y = {}", y[0]);
    }

    /// The other half of op-ad6h is kept: a derivative that is NaN at the
    /// ACCEPTED state is still reported as `NonFiniteState`, never `Ok`.
    #[test]
    fn rkf45_reports_a_non_finite_derivative_at_the_accepted_state() {
        let mut solver = Rkf45::new(1, 1e-8, 1e-6);
        let mut y = vec![3.0_f64];
        let mut dx = 0.1;
        let outcome = solver.integrate(&DecayNanOutsideRange, 0.0, 1.0, &mut y, &mut dx);
        assert_eq!(outcome, Err(OdeError::NonFiniteState));
    }

    #[test]
    fn rkf45_exponential_decay() {
        // y' = -y, y(0) = 1  →  y(1) = e^{-1}
        let ode = DecayOde;
        let mut solver = Rkf45::new(1, 1e-8, 1e-6);
        let mut y = vec![1.0_f64];
        let mut dx = 0.1;
        solver.integrate(&ode, 0.0, 1.0, &mut y, &mut dx).unwrap();
        let expected = (-1.0_f64).r_exp();
        assert!(
            (y[0] - expected).abs() < 1e-6,
            "y={:.10}, expected={:.10}",
            y[0],
            expected
        );
    }

    #[test]
    fn rkf45_decay_to_5() {
        // y' = -y, y(0) = 1  →  y(5) = e^{-5}
        let ode = DecayOde;
        let mut solver = Rkf45::new(1, 1e-9, 1e-7);
        let mut y = vec![1.0_f64];
        let mut dx = 0.1;
        solver.integrate(&ode, 0.0, 5.0, &mut y, &mut dx).unwrap();
        let expected = (-5.0_f64).r_exp();
        assert!(
            (y[0] - expected).abs() < 1e-7,
            "y={:.12}, expected={:.12}",
            y[0],
            expected
        );
    }

    #[test]
    fn rkf45_nonstiff_vdp() {
        // Mildly non-stiff Van der Pol (mu=0.5): y(0) = (2, 0), integrate to t=2
        // Just check it completes without error and the solution is bounded
        let ode = VanDerPol { mu: 0.5 };
        let mut solver = Rkf45::new(2, 1e-7, 1e-5);
        let mut y = vec![2.0_f64, 0.0];
        let mut dx = 0.1;
        solver.integrate(&ode, 0.0, 2.0, &mut y, &mut dx).unwrap();
        assert!(y[0].abs() < 3.0, "y0={}", y[0]);
        assert!(y[1].abs() < 3.0, "y1={}", y[1]);
    }

    #[test]
    fn rkf45_convergence_order_is_fifth() {
        // Order verification: halving h should reduce the global error by ~2^5 = 32.
        // This RKF45 propagates with the 5th-order (B) solution (Dormand-Prince convention).
        // Force fixed-step by setting abs_tol=1 >> expected error so every step is accepted,
        // then manually supply dx_try = h at each iteration.
        let exact = (-1.0_f64).r_exp();

        let run = |n: usize| -> f64 {
            let ode = DecayOde;
            let mut solver = Rkf45::new(1, 1.0, 0.0);
            let mut y = vec![1.0_f64];
            let h = 1.0 / n as f64;
            let mut x = 0.0_f64;
            for _ in 0..n {
                let mut dx_try = h;
                solver
                    .solve_step(&ode, &mut x, &mut y, &mut dx_try)
                    .unwrap();
            }
            (y[0] - exact).abs()
        };

        let err_coarse = run(5); // h = 0.2
        let err_fine = run(10); // h = 0.1
        let ratio = err_coarse / err_fine;
        // 5th-order propagation: ratio ≈ 2^5 = 32; accept [16, 64] to allow higher-order terms
        assert!(
            ratio > 16.0 && ratio < 64.0,
            "order ratio = {ratio:.2} (expected ~32 for RKF45 5th-order propagation)"
        );
    }

    #[test]
    fn rkf45_exceeds_steps_on_stiff_vdp() {
        // Stiff Van der Pol (mu=1000): RKF45 needs h < ~2/mu² ≈ 2e-6 for stability,
        // so it quickly exhausts max_steps even on a short interval.
        let ode = VanDerPol { mu: 1000.0 };
        let mut solver = Rkf45::new(2, 1e-6, 1e-4);
        solver.config.max_steps = 200;
        let mut y = vec![2.0_f64, 0.0];
        let mut dx = 0.1;
        let result = solver.integrate(&ode, 0.0, 0.5, &mut y, &mut dx);
        assert!(
            matches!(result, Err(OdeError::MaxStepsExceeded(_))),
            "expected MaxStepsExceeded, got {:?}",
            result
        );
    }
}