Function integrate_ensemble 

pub fn integrate_ensemble<S>(
    lanes: &[OdeLane<S>],
    solver: &OdeSolver,
    backend: ComputeBackend,
) -> OdeEnsemble
where
    S: OdeSystem + Sync,

Integrate N independent initial-value problems with one stepper for every lane, on the chosen backend.

This is the entry point for the common case: one ODE per cell, per particle, or per material point, all of the same character. Each lane gets its own clone of solver, so no scratch buffer is shared and no lane can perturb another.

For an ensemble of mixed stiffness, use integrate_ensemble_mixed — see the module-level “Stiffness” section.

Arguments

• lanes — one OdeLane per problem. Lane i of the result corresponds to lanes[i].
• solver — the stepper prototype, cloned once per lane. Its OdeSolverConfig tolerances and max_steps apply to every lane. It must have been built for the same equation count the lanes’ systems report, because the steppers size their scratch buffers at construction.
• backend — requested execution backend. What actually runs is ensemble_backend_for applied to it: an unavailable backend degrades, Gpu degrades (no GPU kernel here), and an ensemble below ODE_ENSEMBLE_MIN_LANES runs serially. None of those changes the answer.

Returns

An OdeEnsemble with one OdeLaneSolution per lane, in input order.

Determinism

Bit-for-bit identical across backends and thread counts, for systems whose derivatives/jacobian are pure deterministic functions of their arguments. See the module-level “Determinism” section.

Cost note

solver is cloned once per lane, which clones its per-equation scratch buffers — a handful of small allocations per lane. This is deliberate: it is what makes each lane a pure function of its own inputs, and therefore what makes the bitwise-identity claim hold without depending on every stepper’s buffers being write-before-read. On the measured ensemble it is a small fraction of the per-lane cost — about 8.6 us of integration per lane over 68 adaptive steps, see ODE_ENSEMBLE_MIN_LANES — but a caller integrating a very short interval per lane would see it, and it has not been measured separately from the integration it accompanies.

Panics

Panics if solver is OdeSolver::Rosenbrock23 and a lane’s system does not override OdeSystem::jacobian, whose default implementation panics. On the CpuMulti path the panic propagates out of the rayon scope. Check OdeSolver::requires_jacobian first.

Verification

Methodology. Checked against the closed-form solution of dy/dx = -k y, y(0) = 1, namely y(x) = exp(-k x), over 64 lanes with k spread evenly across [0.5, 8], integrated to x = 1 by all three steppers; and against the harmonic oscillator y1' = -y2, y2' = y1 from y(0) = (1, 0), whose solution is (cos x, sin x), over 16 lanes ending at m * pi/2 for m = 1..=16. Tolerances abs_tol = 1e-10, rel_tol = 1e-8 for Rkf45 and Rosenbrock23 (1e-12 / 1e-10 for the oscillator) and 1e-3 / 1e-2 for Euler, which cannot reach the high-order tolerances inside the default 10 000-step budget. Pass criteria: < 1e-8 for the high-order steppers, < 5e-2 for first-order Euler.

Results, measured 2026-08-13 by ensemble_matches_analytic_decay and ensemble_matches_harmonic_oscillator in parallel/tests.rs, release build: worst absolute error over the 64 decay lanes 1.518896e-9 (Rkf45, 3 038 total accepted steps, 81 in the worst lane), 1.387431e-9 (Rosenbrock23, 46 535 steps, 1 186 in the worst lane) and 1.297370e-3 (Euler, 21 817 steps, 421 in the worst lane); worst absolute error on the oscillator 4.905929e-10 (Rkf45).

Interpretation. The two high-order steppers agree with the closed form to their requested tolerance and with each other to about 1e-9, while Euler is six orders coarser — the expected signature of three genuinely different steppers being reached through the ensemble, which is what rules out the wrapper silently routing every lane to one of them. The step counts are the other half of the story: Rosenbrock23 needs 15x the steps of Rkf45 on a non-stiff problem to reach the same accuracy, which is exactly why integrate_ensemble_mixed exists rather than “just use the stiff solver everywhere”.

Example — the uom boundary

The ensemble is dimensionless, and the caller converts at its edge. These lanes are lumped-capacitance bodies cooling towards ambient, dT/dt = -(T - T_inf) / tau, one lane per time constant, recovering a ThermodynamicTemperature.

On the tolerance asserted below. The stepper’s rel_tol is 1e-8 and the states are of order 400 K, so the achievable absolute accuracy is a few microkelvin — the error floor is set by the controller’s tolerance and the magnitude of the state, not by the ensemble. Measured 2026-08-13 (release) by lumped_body_accuracy_is_set_by_the_relative_tolerance in parallel/tests.rs: worst absolute error 1.150937e-6 K at abs_tol = 1e-10, rel_tol = 1e-8; 1.222958e-8 K at 1e-12/1e-10; 1.289209e-10 K at 1e-14/1e-12. The bound below is 5e-6 K, set from the first of those measurements — a tighter assertion would be asserting something this stepper at this tolerance does not deliver. A caller wanting sub-nanokelvin must ask for it through the tolerances, and the third row says what that buys.

use outram_foam_basic_lib::compute::ComputeBackend;
use outram_foam_basic_lib::ode::{OdeSolver, OdeSystem};
use outram_foam_basic_lib::ode::parallel::{integrate_ensemble, OdeLane};
use uom::si::f64::{ThermodynamicTemperature, Time};
use uom::si::thermodynamic_temperature::kelvin;
use uom::si::time::second;

/// Lumped body: `dT/dt = -(T - T_inf) / tau`, all in SI base units.
struct LumpedBody { tau_s: f64, t_inf_k: f64 }
impl OdeSystem for LumpedBody {
    fn n_eqns(&self) -> usize { 1 }
    fn derivatives(&self, _t: f64, y: &[f64], dydt: &mut Vec<f64>) {
        dydt[0] = -(y[0] - self.t_inf_k) / self.tau_s;
    }
}

// Convert in: typed quantities out to plain f64 in named units.
let t_inf = ThermodynamicTemperature::new::<kelvin>(300.0);
let t0 = ThermodynamicTemperature::new::<kelvin>(500.0);
let horizon = Time::new::<second>(10.0);

let lanes: Vec<OdeLane<LumpedBody>> = [5.0_f64, 20.0]
    .iter()
    .map(|&tau| {
        OdeLane::new(
            LumpedBody { tau_s: tau, t_inf_k: t_inf.get::<kelvin>() },
            vec![t0.get::<kelvin>()],
            0.0,
            horizon.get::<second>(),
            0.1,
        )
    })
    .collect();

let ensemble = integrate_ensemble(
    &lanes,
    &OdeSolver::rkf45(1, 1e-10, 1e-8),
    ComputeBackend::CpuMulti,
);

// Convert out: back to typed temperatures.
let temperatures: Vec<ThermodynamicTemperature> = ensemble
    .states()
    .expect("both lanes complete")
    .iter()
    .map(|s| ThermodynamicTemperature::new::<kelvin>(s[0]))
    .collect();

// Closed form: T(t) = T_inf + (T0 - T_inf) exp(-t / tau). The 5e-6 K bound
// is the measured floor at rel_tol = 1e-8 on a ~400 K state; see above.
for (tau, temperature) in [5.0_f64, 20.0].iter().zip(&temperatures) {
    let exact = 300.0 + 200.0 * (-10.0_f64 / tau).exp();
    assert!((temperature.get::<kelvin>() - exact).abs() < 5e-6);
}