Function quadrature_batch 

pub fn quadrature_batch<F>(
    intervals: &[QuadratureInterval],
    rule: QuadratureRule,
    backend: ComputeBackend,
    f: F,
) -> QuadratureBatch
where
    F: Fn(usize, f64) -> f64 + Sync,

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

The GPU-shaped half of this module: every lane performs exactly QuadratureRule::evaluations integrand calls with no data-dependent branching. Use it when the integrand is smooth and the panel count can be chosen once for the whole batch — a band-averaged cross section, a face-integrated flux, a cell-integrated source term.

Arguments

• intervals — one QuadratureInterval per lane.
• rule — the rule applied to every lane. QuadratureRule::GaussLegendre is the default worth reaching for on a smooth integrand.
• backend — requested backend. What actually runs is quadrature_backend_for applied to it; a batch below QUADRATURE_MIN_INTERVALS runs serially. None of the degradations changes the answer.
• f — the integrand. f(i, x) must return lane i’s integrand at abscissa x. It must be a pure deterministic function of its arguments — see the module-level “Determinism” section. It is called from multiple threads on the CpuMulti path, hence the Sync bound; the bound is present in both feature builds so that enabling parallel never changes a public signature.

Returns

A QuadratureBatch with one QuadratureSolution per interval, in input order.

Determinism

Bit-for-bit identical across backends and thread counts. The sum within a lane is sequential and is never split across threads — see the module-level “Determinism” section for why that restriction is deliberate.

Verification

Methodology. Three oracles, all exact rather than another implementation. (1) Polynomial exactness: an n-point Gauss-Legendre rule must integrate every monomial x^d for d <= 2n - 1 exactly, Simpson must be exact for d <= 3 and trapezoid for d <= 1; checked over [0, 1] and [-2, 3] against the closed form (b^(d+1) - a^(d+1)) / (d + 1). (2) Published nodes: the computed 8-point nodes and weights against the Abramowitz & Stegun 25.4.30 values already carried in this workspace by crates/raffles/src/distributions.rs. (3) Transcendental reference: integral of exp(-x) sin(x) from 0 to pi = (1 + exp(-pi)) / 2. Pass criteria: < 1e-13 relative for exactness, < 1e-15 absolute against A&S, < 1e-12 absolute for the transcendental.

Results, measured 2026-08-13 by gauss_legendre_is_exact_to_its_degree, gauss_nodes_match_the_in_workspace_abramowitz_stegun_values, simpson_and_trapezoid_are_exact_to_their_degree and fixed_rules_match_a_transcendental_reference in parallel/tests.rs, release build:

• Exactness sweep (G2..G8, degrees 0 to 15, both intervals): worst relative error 5.769990e-16. Simpson over degrees 0 to 3: 3.552714e-16. Trapezoid over degrees 0 to 1: 0.000000e0.
• Computed G8 against A&S 25.4.30: worst node difference 1.110223e-16, worst weight difference 1.249001e-16, weights summing to 2.00000000000000000.
• Transcendental reference, closed form 0.52160695913188615: G8 over 8 panels (64 evaluations) gave 0.52160695913188615, error 0.000000e0; Simpson over 64 panels (129 evaluations) error 4.206809e-9; trapezoid over 128 panels (129 evaluations) error 5.236767e-5.

Interpretation. The Gauss nodes and weights are correct to the last bit f64 holds — agreeing with a published table to within one unit in the last place while being computed independently of it — the composite mapping onto arbitrary intervals is correct, and the order hierarchy behaves as theory requires: G8 reaches the rounding floor on a smooth transcendental at half the evaluations where a 64-panel Simpson is still nine orders away and a trapezoid four orders beyond that.

Example — the uom boundary

The batch is dimensionless, and the caller converts at its edge. These lanes integrate a linearly-ramping electrical power over time to recover an Energy:

use outram_foam_basic_lib::compute::ComputeBackend;
use outram_foam_basic_lib::ode::parallel::{
    quadrature_batch, GaussOrder, QuadratureInterval, QuadratureRule,
};
use uom::si::f64::{Energy, Power, Time};
use uom::si::energy::joule;
use uom::si::power::watt;
use uom::si::time::second;

// P(t) = P0 (1 - t / T), shutting down linearly over the window.
let p0 = Power::new::<watt>(2.0e6);
let window = Time::new::<second>(30.0);

// Convert in: limits in seconds, as plain f64.
let intervals = [
    QuadratureInterval::new(0.0, window.get::<second>()),
    QuadratureInterval::new(0.0, 0.5 * window.get::<second>()),
];

let batch = quadrature_batch(
    &intervals,
    QuadratureRule::GaussLegendre { order: GaussOrder::G4, panels: 4 },
    ComputeBackend::CpuMulti,
    |_, t| p0.get::<watt>() * (1.0 - t / window.get::<second>()),
);

// Convert out: back to typed energies.
let energies: Vec<Energy> = batch
    .values()
    .expect("both lanes evaluate")
    .into_iter()
    .map(Energy::new::<joule>)
    .collect();

// Closed form over [0, T] is P0 T / 2; over [0, T/2] it is 3 P0 T / 8.
assert!((energies[0].get::<joule>() - 0.5 * 2.0e6 * 30.0).abs() < 1e-6);
assert!((energies[1].get::<joule>() - 0.375 * 2.0e6 * 30.0).abs() < 1e-6);