Module parallel 

Batched numerical integration on the hybrid execution backend — N independent initial-value problems, or N independent definite integrals, advanced at once, serially or across CPU cores.

The two things in here, and why they share a module

Operation	The batch is	Entry point
ODE ensemble, one stepper for every lane	N independent IVPs	integrate_ensemble
ODE ensemble, stepper chosen per lane	N independent IVPs	integrate_ensemble_mixed
Fixed-rule quadrature	N independent definite integrals	quadrature_batch
Adaptive quadrature	N independent definite integrals	adaptive_quadrature_batch

Both halves are “numerical integration” in the sense bead op-yvj.4.5 uses the phrase, and a definite integral genuinely is the initial-value problem dy/dx = f(x), y(a) = 0, read at x = b — so quadrature sitting beside the ODE steppers is not a filing accident. They also share every piece of machinery that is not the arithmetic itself: the lane vocabulary, the per-lane status reporting, the backend-degradation policy, and the determinism argument below. Keeping them together means the crate has one batched-integration dialect rather than two.

Reuse, not reimplementation

No new integrator is written here. The ODE half drives the crate’s existing steppers — Euler, Rkf45, Rosenbrock23, selected through the existing OdeSolver enum — over the existing adaptive interval loop. This module adds the ensemble: the outer loop over lanes, the per-lane outcome reporting, and the backend dispatch.

The quadrature half is new code, because the workspace had no general quadrature. The only prior art is outram-park-fork-dwsim-libs’s clean_energies::pem_fuel_cell::simpson_integrate, which is deliberately not reused: it is a verbatim port of OPEM’s rule including its documented flaw (composite Simpson weights applied without checking that the sample count is odd), and it integrates a slice of pre-sampled values rather than a callable integrand. QuadratureRule keeps the one convention that does carry over — composite rules over equal subintervals — and makes the sample-count error structurally impossible by counting Simpson panels of two subintervals each rather than raw samples. raffles’s distributions::special::integrate_open_unit is the other in-workspace quadrature: a composite 8-point Gauss-Legendre over geometrically graded panels, hard-wired to the open unit interval and to quantile-function moments. It is not a general routine and is private to that crate’s special module, so it is cited here as precedent for the Gauss-Legendre choice rather than reused.

Hybrid means dispatch, not two APIs

Every entry point takes a ComputeBackend parameter, and there is no *_parallel() sibling anywhere. With the parallel feature off, ComputeBackend::CpuMulti resolves down to ComputeBackend::Serial through ComputeBackend::resolve and the answer is unchanged, bit for bit. There is no GPU kernel in this module yet, so a ComputeBackend::Gpu request degrades to the best available CPU path — see “GPU” below for why that is not merely laziness for the adaptive paths.

Determinism — bitwise identical, and the summation-order question

Every kernel in this module returns bit-for-bit identical output on ComputeBackend::Serial and ComputeBackend::CpuMulti, at any thread count, on every run, provided the caller’s system or integrand is itself a deterministic pure function of its arguments.

For the ODE ensemble this is the same argument as crate::math::parallel’s: lane i’s trajectory is a pure function of lane i’s system, initial condition and stepper. No lane reads another lane’s state, so there is no cross-lane arithmetic whose association could change, and both backends call the very same per-lane kernel. Each lane also gets its own clone of the stepper prototype, so no scratch buffer is ever shared between lanes and the result cannot depend on which lanes a worker happened to run first.

Quadrature needs one extra sentence, because a quadrature rule is a sum and floating-point addition is not associative. The reason the answer is still bit-identical is the shape of the batch: one lane is one integral, and one integral is summed sequentially by a single thread. The parallelism is over lanes, never within a lane, so no partial sum is ever re-associated. This is a deliberate design choice with a real cost — it means a single very expensive integral gets no speed-up from this module at all — and it is taken because a reduction split across threads would give a different answer at every thread count, which for a verification oracle is disqualifying. Splitting one integral’s panels across threads is not offered here, and if it is ever added it must be a separate, separately-named entry point whose documentation says plainly that it is not bit-reproducible.

Verified by the bitwise_* tests in parallel/tests.rs, which compare serial against rayon pools of 1, 2, 4 and 8 workers on batches built to have wildly uneven per-lane cost.

For scale, a 4 096-lane deliberately-imbalanced ensemble (half the lanes decaying at k near 1, half at k near 60, giving 278 318 accepted steps in total — a mean of 67.9 per lane against a maximum of 126) under Rkf45, measured 2026-08-13 on 4 logical cores by the #[ignore]d ensemble_thread_scaling_benchmark, best of 7 samples, with a second independent run alongside:

Worker threads	Time	Speed-up	(repeat)	Bitwise vs serial
serial reference	35808.33 us	1.00x	1.00x	—
1	35839.25 us	1.00x	1.00x	identical
2	18283.82 us	1.96x	1.92x	identical
4	9276.65 us	3.86x	3.09x	identical
8	9128.96 us	3.92x	3.05x	identical

The “identical” column is the determinism claim above measured rather than argued, and it is asserted by the benchmark itself, not merely printed. Going through rayon with a single worker costs essentially nothing here (1.00x), unlike the batched root finder where it costs about 6% — the per-lane work is large enough that the iterator machinery disappears into it. Eight workers on four cores buy nothing further, the expected signature of a compute-bound kernel that already saturates its cores. This is one machine, one ensemble and two runs; it is not a scaling study, and nothing here has been measured on Android hardware or on a many-core server.

Load imbalance — why there is no hand-rolled partition

An adaptive stepper takes a different number of sub-steps in every lane, and the spread is not small. On the benchmark ensemble below it is a mean of 67.9 accepted steps against a maximum of 126; on the stiff pair in parallel/tests.rs a single lane under Rkf45 burns all 10 000 of its allotted steps and still spans only a third of the interval, beside decay lanes finishing in 81. Adaptive quadrature has the same property by construction — 469 to 1 225 integrand evaluations across the three verification integrands. A static equal split across P threads therefore ends up waiting on whichever chunk drew the hard lanes.

Every parallel path here uses rayon’s adaptive splitting with no min_len floor, so an idle worker can steal down to a single lane. That is the deliberate answer to the imbalance, not an oversight. No granularity floor is imposed even on the fixed-rule quadrature path, where the crate’s closed-form polynomial kernels do impose one (crate::math::parallel::POLY_BLOCK): there, every lane provably costs the same handful of flops, whereas here the per-lane cost is set by the caller through QuadratureRule’s panel count and by the cost of the caller’s integrand, so any fixed floor would be wrong for most callers. Work-stealing handles both ends without a number that cannot be justified.

Whatever the splitter does, it cannot change a value — every lane is computed independently of every other.

Stiffness — how a mixed ensemble is handled

The realistic per-cell chemistry ensemble is mixed: most cells are benign and a few are stiff. Three things follow, and all three are deliberate.

1. The stepper is not switched behind the caller’s back. An ensemble run with integrate_ensemble uses one stepper for every lane. A stiff lane handed an explicit stepper does not silently return a wrong answer: the adaptive controller shrinks the step until it either meets tolerance — correct but slow — or runs out, at which point the lane reports OdeLaneStatus::MaxStepsExceeded or OdeLaneStatus::StepSizeUnderflow. Stiffness therefore shows up as a named per-lane failure, not as silent garbage.
2. Per-lane stepper selection exists — integrate_ensemble_mixed takes a closure Fn(usize) -> OdeSolver, so a caller that knows which cells are stiff (or that has just been told by a failed integrate_ensemble pass) can give those lanes Rosenbrock23 and leave the rest on Rkf45. This is the intended recovery path and it costs nothing when unused.
3. Rosenbrock23 needs a Jacobian, and this module does not supply one. OdeSystem::jacobian’s default implementation panics, and a panic inside a rayon worker will propagate out of the batch. Check OdeSolver::requires_jacobian before selecting it for a system you did not write. Batched numerical Jacobians are bead op-yvj.4.6, not this one.

Non-convergence is reported, never swallowed

An ensemble of 10 000 cells in which 3 fail must say so, and must say which 3. Both halves report failure per lane, and make it hard to ignore by construction — the same shape crate::math::parallel uses:

• OdeLaneSolution::state and QuadratureSolution::value return Option, and hand back Some only for a lane that succeeded.
• The raw number is behind the deliberately-named OdeLaneSolution::last_state / QuadratureSolution::last_value, so using a failed lane’s partial answer is a visible decision in the calling code rather than an accident.
• OdeEnsemble::states and QuadratureBatch::values are all-or-nothing: they return Err naming the failure count and the first failing lane, rather than a plausible-looking Vec.

A lane that ran out of steps is never presented as if it had reached x_end. It reports OdeLaneStatus::MaxStepsExceeded, its genuine partial state, and the x it actually reached — which is exactly the information needed to decide whether to sub-cycle it or re-run it stiff. A lane whose inputs were unusable reports OdeLaneStatus::InvalidLane and a NaN state, because there is no honest number to return.

GPU

There is no wgpu kernel here yet; Gpu degrades. Two of these four kernels would map to one if written, and two would not, and the distinction is worth recording rather than rediscovering:

• Fixed-rule quadrature is GPU-shaped. Every lane evaluates the same fixed number of nodes with no data-dependent branching, so a warp never diverges.
• Adaptive quadrature is not, and is CPU-only by design. Its subdivision pattern is decided by the integrand, so neighbouring lanes take different paths through the code and different numbers of evaluations; on SIMT hardware that serialises the divergent branches and needs a per-lane work stack. It stays on the CPU deliberately.
• Adaptive ODE ensembles have the same problem, sharpened. A batch run lockstep at the smallest step any lane needs is simple and correct but wastes the whole ensemble’s budget on its worst member; per-lane step control is correct but divergent. Neither is attempted here.

Note also that WGSL has no f64, and the accuracy consequences of an f32 quadrature or an f32 error controller are unmeasured in this workspace.

Units

Everything here is dimensionless f64, for the same reason crate::math::parallel is: a general integrator has no single physical dimension. One lane’s abscissa may be a time in seconds and another’s a length in metres, and the state vector’s components routinely carry different dimensions from each other.

uom typing is not stripped to get here — it is applied at the boundary, by the caller, exactly as the hybrid-backend epic requires: convert into the batch, convert back out. The doctests on integrate_ensemble and quadrature_batch show that boundary explicitly, one recovering a ThermodynamicTemperature and the other an Energy.

Cargo features and portability

The rayon paths sit behind the crate’s parallel feature, which is off by default; with it off this module still compiles and every entry point still works. rayon is pure Rust with no system component, so everything here compiles and runs on aarch64-linux-android / Termux exactly as on desktop. Nothing in this module is target-gated.

Example — an ensemble of independent decays

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};

/// `dy/dx = -k y` — one lane per decay constant.
struct Decay { k: f64 }
impl OdeSystem for Decay {
    fn n_eqns(&self) -> usize { 1 }
    fn derivatives(&self, _x: f64, y: &[f64], dydx: &mut Vec<f64>) {
        dydx[0] = -self.k * y[0];
    }
}

let lanes: Vec<OdeLane<Decay>> = (1..=4)
    .map(|i| OdeLane::new(Decay { k: i as f64 }, vec![1.0], 0.0, 1.0, 0.1))
    .collect();

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

let states = ensemble.states().expect("all four lanes complete");
for (i, s) in states.iter().enumerate() {
    let exact = (-(i as f64 + 1.0)).exp();
    assert!((s[0] - exact).abs() < 1e-8, "lane {i}: {} vs {exact}", s[0]);
}

// Asking for multi-CPU gives a bit-for-bit identical answer, whether or not
// the `parallel` feature is compiled in.
let serial = integrate_ensemble(
    &lanes,
    &OdeSolver::rkf45(1, 1e-10, 1e-8),
    ComputeBackend::Serial,
);
assert_eq!(states, serial.states().unwrap());

Structs

AdaptiveSettings

Tolerances and work limit for adaptive_quadrature_batch.

OdeEnsemble

A batch of N lane outcomes, in the same order as the lanes handed in.

OdeEnsembleFailure

One or more lanes of an OdeEnsemble did not complete.

OdeLane

One member of an ODE ensemble: a system, its initial condition, and the interval to advance it over.

OdeLaneSolution

The outcome of a single ensemble lane: its status, its final state, how far it got and how much work it took.

QuadratureBatch

A batch of N quadrature outcomes, in the same order as the intervals handed in.

QuadratureBatchFailure

One or more lanes of a QuadratureBatch failed.

QuadratureInterval

One lane of a quadrature batch: the limits of one definite integral.

QuadratureSolution

The outcome of a single quadrature lane.

Enums

GaussOrder

Node count of a fixed-order Gauss-Legendre rule.

OdeLaneStatus

How one lane of an ODE ensemble ended.

QuadratureRule

Which fixed rule quadrature_batch applies to every lane.

QuadratureStatus

How one lane of a quadrature batch ended.

Constants

MAX_ADAPTIVE_DEPTH

Hard ceiling on adaptive bisection depth in adaptive_quadrature_batch, regardless of AdaptiveSettings::max_subdivisions.

ODE_ENSEMBLE_MIN_LANES

Lane count below which a ComputeBackend::CpuMulti request runs integrate_ensemble on the calling thread instead.

QUADRATURE_MIN_INTERVALS

Interval count below which a ComputeBackend::CpuMulti request runs quadrature_batch and adaptive_quadrature_batch on the calling thread instead.

Functions

adaptive_quadrature_batch

Evaluate N independent definite integrals adaptively, on the chosen backend.

ensemble_backend_for

The ComputeBackend the ODE ensemble would actually use for n lanes if asked for requested — without running anything.

integrate_ensemble

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

integrate_ensemble_mixed

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

quadrature_backend_for

The ComputeBackend the quadrature kernels would actually use for n intervals if asked for requested — without running anything.

quadrature_batch

Evaluate N independent definite integrals with a fixed rule, on the chosen backend.