outram_foam_basic_lib/ode/

mod.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/>.

//! Ordinary differential equation solvers for systems `dy/dx = f(x, y)`.
//!
//! Ports the OpenFOAM `ODE` layer: user systems implement the [`OdeSystem`](crate::ode::OdeSystem)
//! trait, and one of the concrete steppers integrates them with adaptive step
//! control — [`Euler`](crate::ode::euler::Euler) (explicit 1st order),
//! [`Rkf45`](crate::ode::rkf45::Rkf45) (explicit Runge-Kutta-Fehlberg 4(5)), and
//! [`Rosenbrock23`](crate::ode::rosenbrock23::Rosenbrock23) (semi-implicit, for stiff systems,
//! requiring a Jacobian). The independent variable `x`, state `y`, and step
//! size are bare `f64` in the caller's own units; tolerances are set through
//! [`OdeSolverConfig`](crate::ode::OdeSolverConfig).
//!
//! # Storing an integrator: [`OdeIntegrator`](crate::ode::integrator::OdeIntegrator)
//!
//! The three steppers above take the system by reference on every call, which
//! is awkward for any caller that wants to *keep* "the integrator for this
//! material point" as a struct field — storing a borrow would force a lifetime
//! parameter, which the workspace design rules forbid.
//!
//! [`integrator`](crate::ode::integrator) solves that with two enums that own
//! what they integrate: [`OdeSolver`](crate::ode::integrator::OdeSolver) selects
//! the stepper, and
//! [`OdeIntegrator`](crate::ode::integrator::OdeIntegrator) selects how the
//! system is supplied —
//! [`OdeIntegrator::TypedState`](crate::ode::integrator::OdeIntegrator::TypedState)
//! (a concrete system owned
//! by value, statically dispatched, **preferred**) or
//! [`OdeIntegrator::DynSystem`](crate::ode::integrator::OdeIntegrator::DynSystem)
//! (an `Arc<dyn OdeSystem + Send + Sync>`, kept by
//! maintainer decision for flexibility). Neither borrows, so neither needs a
//! lifetime.

pub mod euler;
pub mod integrator;
pub mod parallel;
pub mod rkf45;
pub mod rosenbrock23;

use petir::mathf::RealMath;
pub use euler::Euler;
pub use integrator::{
    DynSystemIntegrator, NoTypedSystem, OdeIntegrator, OdeSolver, SharedOdeSystem,
    TypedStateIntegrator,
};
pub use parallel::{
    adaptive_quadrature_batch, ensemble_backend_for, integrate_ensemble, integrate_ensemble_mixed,
    quadrature_backend_for, quadrature_batch, AdaptiveSettings, GaussOrder, OdeEnsemble,
    OdeEnsembleFailure, OdeLane, OdeLaneSolution, OdeLaneStatus, QuadratureBatch,
    QuadratureBatchFailure, QuadratureInterval, QuadratureRule, QuadratureSolution,
    QuadratureStatus,
};
pub use rkf45::Rkf45;
pub use rosenbrock23::Rosenbrock23;

use crate::matrix::SquareMatrix;

// ── ODE system trait ─────────────────────────────────────────────────────────

/// Abstract ODE system `dy/dx = f(x, y)`. Maps to `Foam::ODESystem`.
pub trait OdeSystem {
    /// Number of coupled equations (the length of the state vector `y`).
    fn n_eqns(&self) -> usize;

    /// Fill `dydx` with the derivatives at `(x, y)`.
    fn derivatives(&self, x: f64, y: &[f64], dydx: &mut Vec<f64>);

    /// Fill `dfdx` and `dfdy` with the Jacobian at `(x, y)`.
    ///
    /// Required only by stiff solvers (Rosenbrock23). The default panics.
    fn jacobian(&self, _x: f64, _y: &[f64], _dfdx: &mut Vec<f64>, _dfdy: &mut SquareMatrix) {
        unimplemented!("jacobian not implemented for this ODE system");
    }
}

// ── Solver configuration ─────────────────────────────────────────────────────

/// Common parameters for the adaptive step-size controller.
/// Defaults match `Foam::adaptiveSolver` and `Foam::ODESolver`.
#[derive(Debug, Clone)]
pub struct OdeSolverConfig {
    /// Absolute per-equation tolerance.
    pub abs_tol: f64,
    /// Relative per-equation tolerance.
    pub rel_tol: f64,
    /// Safety factor on the step-size scaling (0 < safeScale < 1).
    pub safe_scale: f64,
    /// Exponent for step *increase*.
    pub alpha_inc: f64,
    /// Exponent for step *decrease*.
    pub alpha_dec: f64,
    /// Minimum scale factor applied per step.
    pub min_scale: f64,
    /// Maximum scale factor applied per step.
    pub max_scale: f64,
    /// Maximum sub-steps for one `integrate()` call.
    pub max_steps: usize,
}

impl Default for OdeSolverConfig {
    fn default() -> Self {
        Self {
            abs_tol: 1e-6,
            rel_tol: 1e-4,
            safe_scale: 0.9,
            alpha_inc: 0.2,
            alpha_dec: 0.25,
            min_scale: 0.2,
            max_scale: 10.0,
            max_steps: 10_000,
        }
    }
}

// ── Error type ───────────────────────────────────────────────────────────────

/// Failure modes of an adaptive integration.
#[derive(Debug, Clone, PartialEq)]
pub enum OdeError {
    /// The step size shrank below machine epsilon while trying to meet the
    /// error tolerance — the system is too stiff for the chosen solver, or the
    /// tolerances are unattainable.
    StepSizeUnderflow,
    /// The interval could not be spanned within `max_steps` sub-steps; carries
    /// the number of steps taken.
    MaxStepsExceeded(usize),
    /// The system produced a non-finite (NaN or infinite) error estimate, so
    /// the state cannot be trusted and integration stopped.
    ///
    /// The usual cause is `derivatives` or `jacobian` returning a non-finite
    /// value — for example a numerically-differenced Jacobian
    /// ([`crate::math::differentiate::NumericalJacobian`]) that could not form
    /// a column, or an evaluation that left the model's valid range.
    ///
    /// # Why this variant exists
    ///
    /// Before bead `op-ad6h` this case did not error at all: it returned
    /// `Ok(())` with a NaN state, because the per-equation error fold used
    /// `f64::max`, which discards NaN. A wrong answer reported as success is
    /// the worst failure mode available to a solver, so it is now a distinct
    /// error rather than being folded into
    /// [`StepSizeUnderflow`](Self::StepSizeUnderflow), which would have named
    /// the wrong cause.
    NonFiniteState,
}

impl std::fmt::Display for OdeError {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            Self::StepSizeUnderflow => write!(f, "ODE step size underflow"),
            Self::MaxStepsExceeded(n) => write!(f, "ODE exceeded {n} steps"),
            Self::NonFiniteState => write!(
                f,
                "ODE produced a non-finite error estimate (NaN or infinite state)"
            ),
        }
    }
}

impl std::error::Error for OdeError {}

// ── Shared utilities ─────────────────────────────────────────────────────────

/// Normalised scalar error — max over all equations.
/// Maps to `Foam::ODESolver::normalizeError`.
pub(crate) fn normalize_error(
    y0: &[f64],
    y: &[f64],
    err: &[f64],
    abs_tol: f64,
    rel_tol: f64,
) -> f64 {
    err.iter()
        .zip(y0)
        .zip(y)
        .map(|((e, &a), &b)| {
            let tol = abs_tol + rel_tol * a.abs().max(b.abs());
            e.abs() / tol
        })
        // NB: `fold(0.0, f64::max)` here was a silent-wrong-answer bug
        // (bead `op-ad6h`). Rust's `f64::max` follows IEEE-754 `maxNum` and
        // DISCARDS a NaN operand, so `f64::max(0.0, NaN) == 0.0`: a NaN error
        // read as *zero* error, the step was accepted as perfectly converged,
        // and `integrate` returned `Ok(())` with a NaN state. Returning NaN
        // instead would not fix it either, because the caller's accept test is
        // `err <= 1.0` and every comparison against NaN is false — the step
        // would be rejected, but `dx` would then shrink to underflow with no
        // indication of why.
        //
        // So a non-finite error is reported as `INFINITY`, which is both true
        // ("this step is unboundedly bad") and actionable: the caller rejects
        // it and, per the check in `adaptive_step`, returns
        // `OdeError::NonFiniteState` rather than a misleading
        // `StepSizeUnderflow`.
        .fold(0.0_f64, |acc, e| {
            if e.is_nan() || acc.is_nan() {
                f64::INFINITY
            } else {
                acc.max(e)
            }
        })
}

/// Adaptive step-size loop shared by all explicit solvers.
///
/// Calls `inner_step(x0, y0, dydx0, dx, y_out) -> err`, retrying with a
/// smaller `dx` whenever `err > 1`. Updates `x`, `y`, and `dx_try`.
/// Matches `Foam::adaptiveSolver::solve`.
pub(crate) fn adaptive_step<Sys: OdeSystem + ?Sized>(
    cfg: &OdeSolverConfig,
    mut inner_step: impl FnMut(f64, &[f64], &[f64], f64, &mut Vec<f64>) -> f64,
    ode: &Sys,
    x: &mut f64,
    y: &mut Vec<f64>,
    dydx0: &mut Vec<f64>,
    y_temp: &mut Vec<f64>,
    dx_try: &mut f64,
) -> Result<(), OdeError> {
    let mut dx = *dx_try;
    ode.derivatives(*x, y, dydx0);

    let err = loop {
        let err = inner_step(*x, y, dydx0, dx, y_temp);
        // A non-finite error means the system produced NaN or an infinite
        // value; shrinking `dx` cannot recover from that, so fail immediately
        // and name the real cause instead of grinding down to a misleading
        // `StepSizeUnderflow`. See `OdeError::NonFiniteState` and bead
        // `op-ad6h`.
        if !err.is_finite() {
            return Err(OdeError::NonFiniteState);
        }
        if err <= 1.0 {
            break err;
        }
        let scale = (cfg.safe_scale * err.r_powf(-cfg.alpha_dec)).max(cfg.min_scale);
        dx *= scale;
        if dx.abs() < f64::EPSILON {
            return Err(OdeError::StepSizeUnderflow);
        }
    };

    *x += dx;
    std::mem::swap(y, y_temp);

    let threshold = (cfg.max_scale / cfg.safe_scale).r_powf(-1.0 / cfg.alpha_inc);
    *dx_try = if err > threshold {
        let scale = (cfg.safe_scale * err.r_powf(-cfg.alpha_inc)).clamp(cfg.min_scale, cfg.max_scale);
        dx * scale
    } else {
        dx * cfg.safe_scale * cfg.max_scale
    };

    Ok(())
}

/// Integrate from `x_start` to `x_end` using repeated adaptive steps.
pub(crate) fn integrate_interval(
    cfg: &OdeSolverConfig,
    step_fn: &mut dyn FnMut(&mut f64, &mut Vec<f64>, &mut f64) -> Result<(), OdeError>,
    x_start: f64,
    x_end: f64,
    y: &mut Vec<f64>,
    dx_est: &mut f64,
) -> Result<(), OdeError> {
    let mut x = x_start;
    let mut dx = *dx_est;
    let mut n_steps = 0usize;

    while x < x_end {
        let dx_try = (x + dx).min(x_end) - x;
        let mut dx_limited = dx_try;
        step_fn(&mut x, y, &mut dx_limited)?;
        dx = dx_limited;
        n_steps += 1;
        if n_steps > cfg.max_steps {
            return Err(OdeError::MaxStepsExceeded(n_steps));
        }
    }

    *dx_est = dx;
    Ok(())
}