Function integrate_ensemble_mixed 

pub fn integrate_ensemble_mixed<S, G>(
    lanes: &[OdeLane<S>],
    solver_of: G,
    backend: ComputeBackend,
) -> OdeEnsemble
where
    S: OdeSystem + Sync,
    G: Fn(usize) -> OdeSolver + Sync,

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

The mixed-stiffness entry point. solver_of(i) is called once for lane i and its return value integrates that lane and no other, so a caller can put Rosenbrock23 on the handful of stiff cells and leave the rest on Rkf45 — paying the LU factorisation only where it is needed.

The natural way to use it is as a second pass: run integrate_ensemble with an explicit stepper, read OdeEnsemble::failures, and re-run just those lanes stiff.

Arguments

• lanes — one OdeLane per problem.
• solver_of — solver_of(i) returns the stepper for lane i. Called exactly once per lane, on whichever thread runs that lane, hence the Sync bound; the bound is present in both feature builds so that enabling parallel never changes a public signature. It must be a pure deterministic function of i — see the module-level “Determinism” section.
• backend — requested backend; see ensemble_backend_for.

Returns

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

Panics

As integrate_ensemble: a lane given OdeSolver::Rosenbrock23 whose system does not implement OdeSystem::jacobian panics.

Example — a stiff lane and a benign one in one ensemble

use outram_foam_basic_lib::compute::ComputeBackend;
use outram_foam_basic_lib::matrix::SquareMatrix;
use outram_foam_basic_lib::ode::{OdeSolver, OdeSystem};
use outram_foam_basic_lib::ode::parallel::{integrate_ensemble_mixed, OdeLane};

/// `dy/dx = -k y`, with an analytic Jacobian so the stiff stepper can run.
struct Decay { k: f64 }
impl OdeSystem for Decay {
    fn n_eqns(&self) -> usize { 1 }
    fn derivatives(&self, _x: f64, y: &[f64], d: &mut Vec<f64>) { d[0] = -self.k * y[0]; }
    fn jacobian(&self, _x: f64, _y: &[f64], dfdx: &mut Vec<f64>, dfdy: &mut SquareMatrix) {
        dfdx[0] = 0.0;
        dfdy.set(0, 0, -self.k);
    }
}

// Lane 0 is benign (k = 1), lane 1 is stiff (k = 5000).
let lanes = vec![
    OdeLane::new(Decay { k: 1.0 }, vec![1.0], 0.0, 1.0, 0.1),
    OdeLane::new(Decay { k: 5000.0 }, vec![1.0], 0.0, 1.0, 0.1),
];

let ensemble = integrate_ensemble_mixed(
    &lanes,
    |i| {
        if i == 1 {
            OdeSolver::rosenbrock23(1, 1e-10, 1e-8)
        } else {
            OdeSolver::rkf45(1, 1e-10, 1e-8)
        }
    },
    ComputeBackend::Serial,
);

assert!(ensemble.all_completed());
let states = ensemble.states().unwrap();
assert!((states[0][0] - (-1.0_f64).exp()).abs() < 1e-8);
assert!(states[1][0].abs() < 1e-8); // exp(-5000) underflows to ~0