Module parallel 

Backend-dispatched kernels for the element-wise field algebra.

What this module is for

A finite-volume solver spends most of a timestep doing element-wise arithmetic on fields that carry one value per mesh cell or per mesh face: adding an explicit source to a residual, scaling by a relaxation factor, forming rho*U, advancing phi += dt*ddt(phi), and taking a norm to decide whether the outer loop has converged. Each of those touches every cell, several times per timestep, so on a large mesh they are the operations worth spreading across CPU cores.

The arithmetic is identical to the serial operators in crate::fields::field, crate::fields::vol_field and crate::fields::surface_field. Only the execution strategy differs — and, for the reductions only, the order in which floating-point values are summed (see “Reduction determinism” below).

One entry point per operation

Hybrid execution here means dispatch, not two APIs. There is exactly one public function per operation and the backend is a parameter:

use outram_foam_basic_lib::compute::ComputeBackend;
use outram_foam_basic_lib::fields::field::Field;
use outram_foam_basic_lib::fields::parallel;

let a = Field::uniform(1_000, 2.0_f64);
let b = Field::uniform(1_000, 3.0_f64);

// Same function, different backend — no `add_parallel` twin exists.
let c_ref = parallel::add(ComputeBackend::Serial, &a, &b);
let c_mt = parallel::add(ComputeBackend::CpuMulti, &a, &b);
assert_eq!(c_ref.as_slice(), c_mt.as_slice());

The size threshold that decides whether ComputeBackend::CpuMulti actually spreads the work is one named, documented function — should_parallelise, reading field_parallel_crossover — never an if n > … scattered through the operators. Retuning the policy is a one-place edit.

Cargo features

Multi-threading lives behind the crate’s parallel feature (rayon), which is off by default so the default build stays dependency-light and Android-clean. The public API in this module is the same either way: with the feature off, ComputeBackend::CpuMulti and ComputeBackend::Gpu transparently run the serial path and every result is unchanged. Nothing here is target-gated — rayon is pure Rust, so with --features parallel this module builds and runs on Android/Termux (target_os = "android") too; std::thread::available_parallelism works there and a phone simply reports fewer cores.

The Gpu backend

There is no GPU kernel for field algebra yet. ComputeBackend::Gpu is accepted and routed to the best available CPU path (multi-threaded when the parallel feature is on, serial otherwise). This is stated rather than silently pretended; when a GPU field kernel exists it is wired in here and nothing about the call sites changes.

Which thread pool

The kernels use rayon’s global pool. A caller that wants a dedicated, explicitly sized pool does not need a second API — bind one with rayon::ThreadPool::install and every call inside routes to it:

// Sketch — requires the `parallel` feature, so it is shown rather than run
// as a doctest (rayon is not a dependency of the default build at all).
let pool = rayon::ThreadPoolBuilder::new().num_threads(8).build().unwrap();
pool.install(|| parallel::axpy_assign(ComputeBackend::CpuMulti, &mut y, dt, &ddt));

The vv_reduction_is_independent_of_thread_count test in this module does exactly this and asserts the answer is unchanged.

The reductions are written so that the answer does not depend on which pool or how many threads it has — see below.

Units

This layer is deliberately unit-agnostic, exactly like the serial field layer it mirrors (see the crate::fields module docs): a Field carries bare f64 / Vector3 / Tensor / SymmTensor values in whatever SI units the caller assigned them. No uom quantity is stripped here because none is present at this layer — this crate’s uom discipline lives in the thermophysics layer and nothing in this module weakens it. Each function’s doc states the units it implies. The one place a physical unit is unavoidable is vol_integral, which multiplies by the mesh’s cell volumes in m^3 and therefore returns [phi]*m^3.

Field name strings never grow here

Every operation that returns a VolField or SurfaceField copies the left operand’s name verbatim and never composes a new one. This is not a style preference. A solver that reassigns a persistent field from an expression containing itself (rho = rho + div(phi)) would double the name string every timestep under compositional naming — 2^step growth, invisible in the field data. That exact bug once drove this crate’s compressible_lid_cavity test to 24 GB and a SIGTERM. See the crate CLAUDE.md “Critical translation gotcha”, the matching notes in crate::fields::vol_field, and the name_does_not_grow_* regression tests in this module, which reassign a field from an expression containing itself 64 times and assert the name length is unchanged.

Reduction determinism

The reductions (sum, l2_norm, dot, vol_integral, …) use a fixed-chunk tree reduction: the slice is cut into consecutive chunks of exactly REDUCTION_CHUNK elements, each chunk is summed sequentially in index order, and the per-chunk partial sums are then combined sequentially in index order on the calling thread.

That buys a strong and deliberate guarantee:

• the parallel reduction depends only on the data and on the compile-time constant REDUCTION_CHUNK;
• it is therefore bit-reproducible run to run and identical for any thread count — unlike a work-stealing par_iter().sum(), whose accumulation tree depends on how rayon happened to split the work that run;
• it is not bit-identical to the serial left-to-right sum, because floating-point addition is not associative. The two differ in the last bits; the measured worst-case deviation is recorded in the V&V test vv_parallel_sum_matches_serial_within_tolerance.

min and max are bit-identical to the serial fold, because min and max are associative.

ComputeBackend::Serial is the deterministic trusted reference — the oracle every parallel result in this module is checked against, matching the convention in outram-mc-libs (src/physics/compute.rs) and boon-lay (src/compute.rs).

Crossover: parallel is slower on small fields

Handing a 50-element addition to four threads loses — the dispatch costs more than the arithmetic. should_parallelise therefore falls back to the serial path below field_parallel_crossover, whose value was measured on real hardware, not guessed; see FIELD_PARALLEL_CROSSOVER for the table.

Constants

FIELD_PARALLEL_CROSSOVER

Element count at or above which ComputeBackend::CpuMulti actually spreads the work across threads. Below it, every operation runs on the calling thread.

REDUCTION_CHUNK

Number of consecutive elements summed sequentially inside one chunk of a parallel reduction.

Functions

add

Element-wise sum c[i] = a[i] + b[i], returning a new field.

add_assign

In-place accumulation y[i] += x[i].

add_surface

Sum of two surface fields — internal faces and every boundary patch.

add_vol

Sum of two volume fields — internal field and every boundary patch.

add_vol_assign

In-place accumulation on a volume field, y += x — internal field and every boundary patch.

axpy

Fused combination c[i] = y[i] + a * x[i] (“axpy”), returning a new field.

axpy_assign

In-place fused combination y[i] += a * x[i] (“axpy”). Allocation-free.

axpy_surface_assign

In-place fused update on a surface field, y += a * x — the flux-field counterpart of axpy_vol_assign (phi += dt * ddt(phi)).

axpy_vol_assign

In-place fused update on a volume field, y += a * x — internal field and every boundary patch.

dot

Inner product of two scalar fields, sum_i a[i]*b[i].

dot_field

Element-wise dot product of two vector fields → scalar field: c[i] = a[i] . b[i].

field_parallel_crossover

The element count at or above which multi-threading is worth its overhead for the field-algebra kernels.

l2_norm

Euclidean (L2) norm, sqrt(sum_i x[i]^2).

max

Largest element. Returns -inf for an empty field, matching Field::max.

mean

Arithmetic mean of all elements; returns 0.0 for an empty field, matching Field::mean.

min

Smallest element. Returns +inf for an empty field, matching Field::min.

pointwise_div

Element-wise quotient of two scalar fields, c[i] = a[i] / b[i].

pointwise_mul

Element-wise product of two scalar fields, c[i] = a[i] * b[i].

scale

Uniform scaling c[i] = a[i] * s, returning a new field.

scale_assign

In-place uniform scaling y[i] *= s. Allocation-free.

scale_by_field

Scale a vector/tensor field by a per-element scalar field: c[i] = v[i] * s[i].

scale_surface

Uniform scaling of a surface field, c = a * s — internal faces and every boundary patch.

scale_vol

Uniform scaling of a volume field, c = a * s — internal field and every boundary patch.

scale_vol_assign

In-place uniform scaling of a volume field, y *= s. Allocation-free, name-preserving. Units: the field becomes [y]*[s].

should_parallelise

The one place that decides Serial vs multi-threaded execution for an element-wise field operation over n elements.

sub

Element-wise difference c[i] = a[i] - b[i], returning a new field.

sub_assign

In-place subtraction y[i] -= x[i]. Allocation-free.

sub_surface

Difference of two surface fields — internal faces and every boundary patch.

sub_vol

Difference of two volume fields — internal field and every boundary patch.

sub_vol_assign

In-place subtraction on a volume field, y -= x. Allocation-free, name-preserving.

sum

Sum of all elements, sum_i x[i].

vol_average

Volume-weighted average over the internal cells, integral(phi dV) / integral(dV).

vol_integral

Volume integral over the internal cells, integral(phi dV) = sum_i phi[i] * V[i].

vol_l2_norm

L2 norm of a volume field’s internal values — the residual measure a solver’s outer loop tests for convergence.

vol_max

Largest internal value of a volume field (-inf on an empty mesh).

vol_min

Smallest internal value of a volume field (+inf on an empty mesh).