Real Parallel Execution of Domain-Decomposition Refinement

What is on this page. The projected-speedup analysis (companion page) quantified the parallelism the dynamic domain decomposition would buy — a projected per-region speedup built from eigensolve counts under a per-level-barrier model — but ran no parallel code. This page is the realisation: the per-level midpoint eigensolves now actually run on a bounded pool of worker processes, the run is proved numerically identical to the serial one, and the measured wall-clock is confronted with that projection. The parallel executor is opt-in — off by default, the refinement loop is the single-process serial path.

This is the answer to “the decomposition’s whole purpose is parallel refinement — does it actually deliver?” The headline is honest and now complete: the equivalence is exact, and parallelism genuinely wins once the problem is big enough — measured 1.3× at n_dof ≈ 16k, with the crossover near n_dof ≈ 5k (right where the §4.4 6D problem lives). At small sizes process overhead dominates (raw speedup < 1), and a naive worker pool collapses at large sizes from BLAS thread oversubscription — fixed by pinning each worker to one BLAS thread (see “Three regimes” below). Parallelism is a large-problem lever, correct at every size and fast at the sizes that matter.

Execution model — solve-grain task pool

The eigensolve (ARPACK shift-invert / LOBPCG over the full FE matrices) is the dominant cost and is exactly what the projection counts. Crucially, a solve needs only K(µ), M(µ) from the parametric problem — not the snapshot store. The store is consumed only by the cheap a-posteriori matching step, which stays in the main process. So the heavy unit is memory-light to ship to a worker, and every store + the graph stay single-process.

The per-level loop becomes:

Test edges (main process; reconstructs stored snapshots, runs matching).
Region decision → freeze / split (main process; the synchronisation barrier — the level-end rule that freezes certified regions of the parameter box and splits busy ones, unchanged from the domain-decomposition prototype).
Bisect flagged edges: the batch of midpoint µ’s is computed in the main process (cheap, pure 0.5·(a+b)), dispatched to a bounded ProcessPoolExecutor, and the returned eigenpairs are integrated back — store insertion + graph mutation — in flagged-edge order, not worker completion order.

The seed grid (level 0) always runs serially; only the per-level refinement batches — the work the companion projection scores — are parallelised.

                    main process                         workers
  level ℓ:
    test edges ───────────────┐
    region decide / freeze/split (barrier)
    midpoints = [0.5(a+b) …] ──┼──▶ eigensolve(µ₁)  ┐
                               │    eigensolve(µ₂)  ├ bounded pool
                               │    …               ┘ (worker cap)
    integrate in edge order ◀──┘   (results returned in input order)
      store insertion + graph mutation

The worker-pool configuration — the worker cap and the process start method — is the only knob. With the pool disabled (the default) the run is the single-process serial path, byte-for-byte unchanged.

Why it is OOM-safe (the laptop is the only host)

Workers carry one FE problem each and no store. Under the fork start method (Linux default) that problem is a copy-on-write share of the parent’s already-assembled affine cache — near-zero extra resident memory per worker.
The worker cap is a hard bound on concurrency; nothing in the batch path pickles a snapshot store or a cold tier.
The 6D validation point is explicitly bounded (capped n_dof, levels, and workers) and the caps are logged into the artifact.

The equivalence guarantee (HARD gate)

A parallel run is numerically identical to the serial run on the same config. Two ingredients make this exact rather than approximate:

Deterministic solve. The shift-invert eigensolver passes a fixed ARPACK start vector v0 = 𝟙/√n instead of ARPACK’s default random one. The eigenvector gauge in (near-)degenerate blocks is then a function of µ alone, so the same µ gives bit-identical eigenpairs no matter which process solves it. (This is the same fix the 6D reduced-order-model accuracy measurement needed; here it was promoted from a test fixture to a shipped solver option.)
Deterministic integration. Results are gathered for the whole batch and then applied in flagged-edge order — the same order the serial loop mutates the graph — so the graph that is built (node indices, edges, verdicts) is byte-identical. Worker completion order cannot leak in.

Because the parallel path feeds its precomputed eigenpairs back through the same bisection mutation the serial path uses, no graph-mutation logic is duplicated — only where the solve runs differs.

A regression test pins this on a 1D and a 2D config: identical grid nodes, identical per-edge verdicts, identical stored eigenvalues, and eigenvector subspace angle < 1e-8. It also pins that the result is independent of the worker count (1 vs 4) and repeatable across runs.

Measured vs projected

Two distinct measurements, because they answer different questions.

1. Per-region measured model (comparable to the projection)

Bucket the real per-solve wall-times by region and apply the same per-level-barrier model the projection used (critical path = busiest region’s time). Speedup, efficiency, and load balance are dimensionless ratios, so the measured number (seconds) is directly comparable to the projected number (eigensolve counts). The gap between them is exactly the effect of solve-cost heterogeneity — the projection assumed every eigensolve costs the same; this shows whether that holds.

Case	regions (peak)	term L	projected speedup	measured-model speedup
1D (two-regime prototype)	2	2	1.20×	1.48×
2D (§3.3.3 diffusion)	3	2	1.09×	1.00×
6D (L-shape, bounded)	3	2	1.16×	1.77×

The projected column reproduces the companion page’s projection (1D 1.20×, 2D 1.09×) — confirming the two layers are scored on the same footing. The measured-model column diverges where solve cost is non-uniform: in 1D and 6D the busy region’s solves happen to be cheaper than average, so the measured model is more optimistic than the count-based projection; in 2D the run terminates before the regions carry refinement work, so the per-region model has a single active region per refinement level and reads 1.0 (the projection’s >1 comes from the seed-lattice distribution, which the measured model — refinement-levels only — does not see). All three are honest readings, not failures.

2. Raw end-to-end wall-clock (did multiprocessing pay off?)

The wall time of the whole refinement run, serial vs parallel:

Case	serial (s)	parallel (s)	raw speedup
1D	0.011	0.031	0.36×
2D	0.076	0.103	0.74×
6D (bounded)	0.228	0.356	0.64×

Raw speedup is below 1.0 in every bounded case — multiprocessing makes these small runs slower, because process spawn, argument/result pickling, and the per-level barrier together cost more than the eigensolves save when each solve is sub-millisecond-to-millisecond. These bounded sizes are the overhead-bound regime only; the next section measures the crossover where parallel wins. Eigenvalue max relative difference vs serial is 0.0 (bit-identical) in all three cases — the overhead buys nothing here because the compute is trivially small.

Overhead accounting

The wall-clock gap parallel − serial is dominated by:

Pool/worker spawn. Even with fork, creating the pool and its workers is a fixed cost paid once per run (the pool is created lazily on the first non-empty batch and reused across levels, so it is not paid per level).
Pickling. Each task ships a µ vector (tiny) and returns an (eigenvalues, eigenvectors) block (n_dof × k floats) — the return pickle grows with n_dof and is the per-solve tax that only pays off once the solve itself is expensive enough.
Barrier synchronisation. The region freeze/split decision is a hard per-level barrier; a level cannot start its solves until the previous level’s results are fully integrated.

The crossover where parallelism wins requires a larger n_dof (so each solve dwarfs its pickle), more flagged edges per level (so the bounded pool stays saturated), or deeper refinement (so the fixed spawn cost amortises). The bounded cases above sit entirely in the overhead-bound regime; the crossover itself is measured in the next section.

Three regimes, and the BLAS-oversubscription fix

Sweeping n_dof on the 2D problem (four workers) makes the full picture concrete:

n_dof	per-solve	serial	parallel	raw speedup
196	3 ms	48 ms	186 ms	0.26×
1 521	11 ms	200 ms	285 ms	0.70×
4 761	37 ms	553 ms	585 ms	0.95×
9 801	105 ms	1689 ms	1340 ms	1.26×
16 641	373 ms	5656 ms	4309 ms	1.31×

Three regimes: (1) overhead-bound (n_dof < ~2k) — pool spawn + result pickling exceed the sub-10 ms solves, raw speedup < 1; (2) crossover (n_dof ~ 5k) — solves grow heavy enough to amortise the overhead; (3) win (n_dof ≥ ~10k) — parallel is genuinely faster. The real §4.4 6D problem (n_dof ≈ 6000) sits right at the crossover.

That regime-3 win only holds because each worker is pinned to one BLAS thread (the worker-pool configuration’s default). Without it, each worker’s ARPACK shift-invert spawns an all-core BLAS, so several worker processes × all-core BLAS oversubscribe the CPU and collapse the speedup — measured 0.51× at n_dof = 16641, versus 1.31× with pinning. The pin is applied per-worker by runtime thread-pool limiting — the only mechanism that works under fork, where OMP_NUM_THREADS is read by BLAS at import, before the worker forks; the serial path’s multithreaded BLAS is untouched. Equivalence is unchanged by the thread count (pinned ≡ unpinned ≡ serial, bit-for-bit) — pinning trades only how each solve uses cores, never what it computes. The worker count defaults to two fewer than the CPU count, so the pinned single-thread workers fill the cores without starving the OS or the main process.

Provenance

The measurements were produced by an instrumented driver in the project repository. Per-case artifacts record the equivalence verdict, the raw wall-clock, the measured-vs-projected comparison, and the 6D caps. The wall-clock numbers are a point-in-time measurement on one machine, not a pinned invariant — the pinned guarantee is the equivalence, enforced by a regression test.

Out of scope. The region-aware reduced-order-model evaluation path remains unimplemented — the parallel executor never calls it (it parallelises refinement solves only). That is a separate follow-up to the domain-decomposition prototype.

Related. The projection this realises: region-parallel-speedup.md. The decomposition design and the freeze/split criterion are set out in the domain-decomposition design note, not yet published. The interface-scalability frontier the high-dimensional follow-up inherits — whether the region boundary, a (d−1)-dimensional surface of never-compressed full-fidelity snapshots, stays cheap as the parameter dimension grows — belongs to a companion catalogue of open scalability problems, not yet published.

# Real Parallel Execution of Domain-Decomposition Refinement > **What is on this page.** The [projected-speedup analysis (companion > page)](region-parallel-speedup.md) quantified the parallelism the dynamic > domain decomposition *would* buy — a **projected** per-region speedup built > from eigensolve counts under a per-level-barrier model — but ran no parallel > code. This page is the realisation: the per-level midpoint eigensolves now > actually run on a bounded pool of worker processes, the run is **proved > numerically identical** to the serial one, and the measured wall-clock is > confronted with that projection. The parallel executor is opt-in — off by > default, the refinement loop is the single-process serial path. This is the answer to *"the decomposition's whole purpose is parallel refinement — does it actually deliver?"* The headline is honest and now complete: **the equivalence is exact, and parallelism genuinely wins once the problem is big enough — measured 1.3× at `n_dof ≈ 16k`, with the crossover near `n_dof ≈ 5k` (right where the §4.4 6D problem lives).** At small sizes process overhead dominates (raw speedup < 1), and a naive worker pool *collapses* at large sizes from BLAS thread oversubscription — fixed by pinning each worker to one BLAS thread (see "Three regimes" below). Parallelism is a *large-problem* lever, correct at every size and fast at the sizes that matter. --- ## Execution model — solve-grain task pool The eigensolve (ARPACK shift-invert / LOBPCG over the full FE matrices) is the dominant cost and is exactly what the projection counts. Crucially, **a solve needs only `K(µ), M(µ)` from the parametric problem — not the snapshot store.** The store is consumed only by the cheap a-posteriori matching step, which stays in the main process. So the heavy unit is memory-light to ship to a worker, and every store + the graph stay single-process. The per-level loop becomes: 1. **Test edges** (main process; reconstructs stored snapshots, runs matching). 2. **Region decision → freeze / split** (main process; the synchronisation barrier — the level-end rule that freezes certified regions of the parameter box and splits busy ones, unchanged from the domain-decomposition prototype). 3. **Bisect flagged edges**: the batch of midpoint µ's is computed in the main process (cheap, pure `0.5·(a+b)`), dispatched to a **bounded** `ProcessPoolExecutor`, and the returned eigenpairs are integrated back — store insertion + graph mutation — **in flagged-edge order**, not worker completion order. The seed grid (level 0) always runs serially; only the per-level refinement batches — the work the companion projection scores — are parallelised. ``` main process workers level ℓ: test edges ───────────────┐ region decide / freeze/split (barrier) midpoints = [0.5(a+b) …] ──┼──▶ eigensolve(µ₁) ┐ │ eigensolve(µ₂) ├ bounded pool │ … ┘ (worker cap) integrate in edge order ◀──┘ (results returned in input order) store insertion + graph mutation ``` The worker-pool configuration — the worker cap and the process start method — is the only knob. With the pool disabled (the default) the run is the single-process serial path, byte-for-byte unchanged. ### Why it is OOM-safe (the laptop is the only host) * Workers carry **one FE problem each and no store**. Under the `fork` start method (Linux default) that problem is a copy-on-write share of the parent's already-assembled affine cache — near-zero extra resident memory per worker. * The worker cap is a hard bound on concurrency; nothing in the batch path pickles a snapshot store or a cold tier. * The 6D validation point is explicitly **bounded** (capped `n_dof`, levels, and workers) and the caps are logged into the artifact. --- ## The equivalence guarantee (HARD gate) A parallel run is **numerically identical** to the serial run on the same config. Two ingredients make this exact rather than approximate: 1. **Deterministic solve.** The shift-invert eigensolver passes a fixed ARPACK start vector `v0 = 𝟙/√n` instead of ARPACK's default *random* one. The eigenvector gauge in (near-)degenerate blocks is then a function of `µ` alone, so the same `µ` gives bit-identical eigenpairs no matter which process solves it. (This is the same fix the 6D reduced-order-model accuracy measurement needed; here it was promoted from a test fixture to a shipped solver option.) 2. **Deterministic integration.** Results are gathered for the whole batch and then applied in flagged-edge order — the *same* order the serial loop mutates the graph — so the graph that is built (node indices, edges, verdicts) is byte-identical. Worker completion order cannot leak in. Because the parallel path feeds its precomputed eigenpairs back through the *same* bisection mutation the serial path uses, no graph-mutation logic is duplicated — only *where the solve runs* differs. A regression test pins this on a 1D and a 2D config: identical grid nodes, identical per-edge verdicts, identical stored eigenvalues, and eigenvector subspace angle `< 1e-8`. It also pins that the result is independent of the worker count (1 vs 4) and repeatable across runs. --- ## Measured vs projected Two distinct measurements, because they answer different questions. ### 1. Per-region measured model (comparable to the projection) Bucket the **real per-solve wall-times** by region and apply the *same* per-level-barrier model the projection used (critical path = busiest region's time). Speedup, efficiency, and load balance are dimensionless **ratios**, so the measured number (seconds) is directly comparable to the projected number (eigensolve counts). The gap between them is exactly the effect of solve-cost **heterogeneity** — the projection assumed every eigensolve costs the same; this shows whether that holds. | Case | regions (peak) | term L | projected speedup | measured-model speedup | |------|---------------:|-------:|------------------:|-----------------------:| | 1D (two-regime prototype) | 2 | 2 | **1.20×** | 1.48× | | 2D (§3.3.3 diffusion) | 3 | 2 | **1.09×** | 1.00× | | 6D (L-shape, bounded) | 3 | 2 | **1.16×** | 1.77× | The projected column **reproduces the companion page's projection** (1D 1.20×, 2D 1.09×) — confirming the two layers are scored on the same footing. The measured-model column diverges where solve cost is non-uniform: in 1D and 6D the busy region's solves happen to be cheaper than average, so the measured model is *more* optimistic than the count-based projection; in 2D the run terminates before the regions carry refinement work, so the per-region model has a single active region per refinement level and reads 1.0 (the projection's >1 comes from the seed-lattice distribution, which the measured model — refinement-levels only — does not see). All three are honest readings, not failures. ### 2. Raw end-to-end wall-clock (did multiprocessing pay off?) The wall time of the whole refinement run, serial vs parallel: | Case | serial (s) | parallel (s) | raw speedup | |------|-----------:|-------------:|------------:| | 1D | 0.011 | 0.031 | **0.36×** | | 2D | 0.076 | 0.103 | **0.74×** | | 6D (bounded) | 0.228 | 0.356 | **0.64×** | **Raw speedup is below 1.0 in every bounded case** — multiprocessing makes these *small* runs slower, because process spawn, argument/result pickling, and the per-level barrier together cost more than the eigensolves save when each solve is sub-millisecond-to-millisecond. These bounded sizes are the **overhead-bound regime** only; the next section measures the crossover where parallel wins. Eigenvalue max relative difference vs serial is **0.0** (bit-identical) in all three cases — the overhead buys nothing here because the compute is trivially small. ### Overhead accounting The wall-clock gap `parallel − serial` is dominated by: * **Pool/worker spawn.** Even with `fork`, creating the pool and its workers is a fixed cost paid once per run (the pool is created lazily on the first non-empty batch and reused across levels, so it is *not* paid per level). * **Pickling.** Each task ships a `µ` vector (tiny) and returns an `(eigenvalues, eigenvectors)` block (`n_dof × k` floats) — the return pickle grows with `n_dof` and is the per-solve tax that only pays off once the solve itself is expensive enough. * **Barrier synchronisation.** The region freeze/split decision is a hard per-level barrier; a level cannot start its solves until the previous level's results are fully integrated. The crossover where parallelism wins requires a larger `n_dof` (so each solve dwarfs its pickle), more flagged edges per level (so the bounded pool stays saturated), or deeper refinement (so the fixed spawn cost amortises). The bounded cases above sit entirely in the overhead-bound regime; the crossover itself is **measured** in the next section. ### Three regimes, and the BLAS-oversubscription fix Sweeping `n_dof` on the 2D problem (four workers) makes the full picture concrete: | n_dof | per-solve | serial | parallel | raw speedup | |------:|----------:|-------:|---------:|------------:| | 196 | 3 ms | 48 ms | 186 ms | 0.26× | | 1 521 | 11 ms | 200 ms | 285 ms | 0.70× | | 4 761 | 37 ms | 553 ms | 585 ms | 0.95× | | 9 801 | 105 ms | 1689 ms | 1340 ms | **1.26×** | | 16 641 | 373 ms | 5656 ms | 4309 ms | **1.31×** | Three regimes: **(1) overhead-bound** (`n_dof < ~2k`) — pool spawn + result pickling exceed the sub-10 ms solves, raw speedup `< 1`; **(2) crossover** (`n_dof ~ 5k`) — solves grow heavy enough to amortise the overhead; **(3) win** (`n_dof ≥ ~10k`) — parallel is genuinely faster. The real §4.4 6D problem (`n_dof ≈ 6000`) sits right at the crossover. That regime-3 win only holds **because each worker is pinned to one BLAS thread** (the worker-pool configuration's default). Without it, each worker's ARPACK shift-invert spawns an all-core BLAS, so several worker processes × all-core BLAS oversubscribe the CPU and *collapse* the speedup — measured **0.51× at `n_dof = 16641`**, versus **1.31×** with pinning. The pin is applied per-worker by runtime thread-pool limiting — the only mechanism that works under `fork`, where `OMP_NUM_THREADS` is read by BLAS at import, before the worker forks; the serial path's multithreaded BLAS is untouched. Equivalence is unchanged by the thread count (pinned ≡ unpinned ≡ serial, bit-for-bit) — pinning trades only *how* each solve uses cores, never *what* it computes. The worker count defaults to two fewer than the CPU count, so the pinned single-thread workers fill the cores without starving the OS or the main process. --- ## Provenance The measurements were produced by an instrumented driver in the project repository. Per-case artifacts record the equivalence verdict, the raw wall-clock, the measured-vs-projected comparison, and the 6D caps. The wall-clock numbers are a point-in-time measurement on one machine, not a pinned invariant — the pinned guarantee is the *equivalence*, enforced by a regression test. > **Out of scope.** The region-aware reduced-order-model evaluation path > remains unimplemented — the parallel executor never calls it (it parallelises > *refinement* solves only). That is a separate follow-up to the > domain-decomposition prototype. --- **Related.** The projection this realises: [region-parallel-speedup.md](region-parallel-speedup.md). The decomposition design and the freeze/split criterion are set out in the domain-decomposition design note, not yet published. The interface-scalability frontier the high-dimensional follow-up inherits — whether the region boundary, a `(d−1)`-dimensional surface of never-compressed full-fidelity snapshots, stays cheap as the parameter dimension grows — belongs to a companion catalogue of open scalability problems, not yet published.