experiment: add opt-in projected conservative explosion

docs/projected-explosion-experiment.md

@@ -0,0 +1,34 @@
+# Projected explosion: isolated spatial-interpolation experiment
+
+Historical baseline: `9805383f9c50d867e8555471baf5dbda8333a652`. This candidate does not replace or relabel its results. The TG runner records the actual candidate revision, dirty state, source hash and harness hashes, plus the historical baseline revision and selected explosion scheme.
+
+Default remains `constant`. Explicit `explosion="projected_linear"` supports only a periodic, body-free, two-level, fully blocked single rectangle (window or singleton union representation). Dense, multilevel, multi-rectangle and body/Jet cases are rejected. Constructor guards have been statically reviewed; GPU construction is not yet exercised.
+
+For each population, form a centered gradient only where both axis neighbors are original explosion donors, owned and uncovered. Otherwise set that axis derivative to zero. Project the gradient perpendicular to the population lattice direction (rest population retains its gradient), then evaluate at child offsets ±1/4 coarse cell. Read owner cores, never stale block halos. Repeating the donor for a missing axis avoids invalid reads. All neighbor reads remain within the original coarse transfer footprint.
+
+Timing, two ghost-only streaming hops, four-child coalescence, relaxation, shifted-storage convention, bulk collision and legacy constant kernel are unchanged. There is no NEQ rescaling, limiter, time interpolation or new framework. Projection commutes with subtracting the constant lattice weights.
+
+Method attribution: Schornbaum and Rüde, [arXiv:1508.07982](https://arxiv.org/pdf/1508.07982), Algorithm 2; Chen et al., Physica A 362 (2006) 158–167, [DOI 10.1016/j.physa.2005.09.036](https://doi.org/10.1016/j.physa.2005.09.036). Independently implemented mathematical stencil, with no GPL implementation copied. The Chen full text was not directly audited in this experiment.
+
+## CPU checks
+
+Run `PYTHONPATH=src python -m unittest discover -s tests/unit -p test_projected_explosion.py -v`.
+
+- Constant populations, affine available-axis gradients, characteristic constancy, quartet means and shift invariance.
+- Original donor footprint only; poison covered cells and check finite reconstruction.
+- Both window and singleton-union layouts, centered/off-axis/minimal near-edge rectangles.
+- Eight complete intervals, identity and SRT collision: actual fine pulls and coarse restriction tables, all eight grazing corner paths included; composite mass and both momentum components conserved within 2e-12 FP64 absolute tolerance.
+- Actual CUDA kernel body compiled with a CPU C++ shim: blocked address mapping, FP32 shifted/unshifted inputs, poisoned halo/covered storage, independent NumPy parity. This checks arithmetic and indexing, not CUDA compilation or execution.
+- Unprojected reconstruction is distinguished by failure of characteristic constancy. It need not fail global mass conservation in this geometry.
+
+## Next diagnostic, only on authorized Horsea GPU
+
+Use the same candidate checkout for paired Re100 N64 depth2 runs, differing only in `--explosion constant` versus `--explosion projected_linear`:
+
+`python tests/validation/taylor_green/run_case.py --repo "$PWD" --re 100 --n 64 --depth 2 --explosion projected_linear --output <new-result-path>`
+
+Select the physical CUDA device explicitly. Never overwrite historical reports. Compare frozen TG error gates and full trajectories, including mass, energy, pressure and velocity. Amplitude/nonmodal shape decomposition can inform mechanism but cannot substitute for the gates. First establish CUDA compilation and CPU/device parity, then paired end-to-end budgets. No accuracy improvement, production suitability or GPU pass is claimed by these CPU results.
+
+## Local verification limits
+
+Focused candidate/refinement/block-layout/build-contract checks pass (62 unittest tests). Full unittest discovery is not green: the attempted 200-test run reports 18 missing-dependency/import errors and three historical source-pin refusals (jet-circle and LES guards), plus one skipped test. Those pin refusals correctly reject the changed `blocks.cu`; their immutable fixtures were not weakened or rewritten. Python syntax checks and `git diff --check` pass. No CUDA build, device execution, GPU performance or TG improvement has been established.

src/CelerisLab/lbm/kernels/blocks.cu

@@ -60,6 +60,39 @@ extern "C" __global__ void BlockExplode(
         to[index_f(t % TOTAL_CELLS, q)] = from[q * coarse_cells + d % coarse_cells];
 }
 
+// Experimental Schornbaum--Ruede (arXiv:1508.07982, Algorithm 2)
+// and Chen et al. (doi:10.1016/j.physa.2005.09.036) population interpolation.
+// Independent implementation of the mathematical stencil; no upstream source copied.
+// Stencil cells are owner cores; a disabled axis repeats its center twice.
+extern "C" __global__ void BlockExplodeProjected(
+    const fpxx* coarse, fpxx* fine, const uint64_t* target, const uint64_t* donor,
+    unsigned int count, unsigned long coarse_cells,
+    const uint64_t* neighbors, const float* offsets)
+{
+    unsigned int row = blockIdx.x * blockDim.x + threadIdx.x;
+    if (row >= count) return;
+    const uint64_t t = target[row], d = donor[row];
+    fpxx* to = fine + (t / TOTAL_CELLS) * NQ * TOTAL_CELLS;
+    for (unsigned int q = 0; q < NQ; ++q) {
+        float v[4];
+        for (unsigned int a = 0; a < 4; ++a) {
+            const uint64_t n = neighbors[4ul * row + a];
+            v[a] = coarse[(n / coarse_cells) * NQ * coarse_cells
+                          + q * coarse_cells + n % coarse_cells];
+        }
+        float gx = .5f * (v[1] - v[0]), gy = .5f * (v[3] - v[2]);
+        const float ex = d_cx[q], ey = d_cy[q];
+        if (q != 0) {
+            const float along = (ex * gx + ey * gy) / (ex * ex + ey * ey);
+            gx -= ex * along;
+            gy -= ey * along;
+        }
+        to[index_f(t % TOTAL_CELLS, q)] =
+            coarse[(d / coarse_cells) * NQ * coarse_cells + q * coarse_cells + d % coarse_cells]
+            + offsets[2ul * row] * gx + offsets[2ul * row + 1] * gy;
+    }
+}
+
 // Targets are receiver-local SOURCE slots for the immediately following pull.
 // Keep child summation order identical to PatchCoarseStep; no collision here.
 extern "C" __global__ void BlockCoalesce(

src/CelerisLab/lbm/refinement_stepper.py

@@ -22,7 +22,7 @@ from .block_stepper import MAX_GRID_DIM
 from .hierarchy import build_hierarchy, is_rectangle
 from .refinement import ACTIVE, HALO, FixedPatch2D, FixedPatchSet2D, prepare_fine_boundary
 from .refinement_tables import (blocked_base_tables, blocked_fine_tables,  # noqa: F401 (re-export)
-                                blocked_kind, _prepare_blocked_boundary,
+                                blocked_kind, projected_explosion_tables, _prepare_blocked_boundary,
                                 BLOCKED_BOUNDARY_DEVICE_COLUMNS)
 from .descriptors import SOLID
 from .raw_exchange import summarize_exchange
@@ -41,6 +41,7 @@ class _Pair:
     n_ghosts: int
     n_rows: int
     parent_cells: int    # parent block window cells (explode donor stride)
+    interpolation: tuple = ()  # optional neighbor/offset device tables
 
 
 @dataclass
@@ -216,7 +217,7 @@ class FixedPatchStepper:
                  fine_block=None, base_block=None, jets=(), child=None, child_block=None,
                  child_objects=(), child_jets=(), child_rectangles_in_parent=None,
                  refinements=None, level_blocks=None, level_objects=None, level_jets=None,
-                 observe_raw_exchange=False):
+                 observe_raw_exchange=False, explosion="constant"):
         if HALO != 4:
             raise RuntimeError("refinement.cu requires a four-cell fine halo")
         if coarse_boundary not in ("periodic", "slot_reservoir", "channel"):
@@ -245,6 +246,15 @@ class FixedPatchStepper:
         level_blocks = [None if b is None else tuple(b) for b in level_blocks]
         _check_levels(refinements, level_blocks, level_objects, level_jets)
         depth = len(refinements) + 1
+        if explosion not in ("constant", "projected_linear"):
+            raise ValueError("explosion must be constant or projected_linear")
+        if explosion == "projected_linear" and (
+                depth != 2 or coarse_boundary != "periodic"
+                or any(b is None for b in level_blocks)
+                or any(level_objects) or any(level_jets)
+                or not (is_rectangle(refinements[0]) or len(refinements[0]) == 1)):
+            raise ValueError("projected explosion requires periodic two-level blocked single rectangle without bodies/Jets")
+        self._explosion = explosion
         self._ddf_shifting = bool(cfg.ddf_shifting)
         if self._ddf_shifting and any(b is None for b in level_blocks):
             raise ValueError("shifted populations require all-blocked levels; dense/mixed is unsupported")
@@ -467,6 +477,10 @@ class FixedPatchStepper:
                        self._upload(children), len(ghosts), len(targets),
                        int(np.prod(pair.coarse_shape if parent.layout is None
                                    else parent.layout.local_shape)))
+        if self._explosion == "projected_linear":
+            neighbors, offsets = projected_explosion_tables(pair, parent.layout)
+            result.explode = child.module.get_function("BlockExplodeProjected")
+            result.interpolation = (self._upload(neighbors), self._upload(offsets))
         if child.layout is not None:
             child.kind = self._upload(blocked_kind(child.layout, self.hierarchy.kinds[index + 1]))
         if index == 0 and base is not None:
@@ -731,7 +745,7 @@ class FixedPatchStepper:
         """Explode from the given parent buffer, then two hops of level ``index``."""
         pair, level = self._pairs[index - 1], self._levels[index]
         pair.explode(parent_buffer, level.current, pair.ghosts, pair.donors, np.uint32(pair.n_ghosts),
-                     np.uint64(pair.parent_cells),
+                     np.uint64(pair.parent_cells), *pair.interpolation,
                      grid=((pair.n_ghosts + self._block[0] - 1) // self._block[0], 1), **kw)
         self._hop(index, level.current, level.spare, **kw)
         self._hop(index, level.spare, level.current, **kw)

src/CelerisLab/lbm/refinement_tables.py

@@ -66,6 +66,40 @@ def blocked_base_tables(patch, layout, order, wall_aware=False, *, composed_kind
             blocked_kind(layout, dense_kind))
 
 
+def projected_explosion_tables(patch, layout):
+    """Experimental periodic single-rectangle stencil, owner-core cell IDs.
+
+    Schornbaum & Rüde (2016), arXiv:1508.07982, Algorithm 2 explosion
+    interpolation, following Chen et al., Physica A 362 (2006), 158-167,
+    doi:10.1016/j.physa.2005.09.036: centered gradients, perpendicular projection,
+    and child offsets +/-1/4. Independently implemented from the mathematics.
+    Only original explosion donors may supply gradients: this preserves the
+    existing coarse transfer footprint. Missing axis pairs use center twice.
+    Caller must restrict use to two periodic, body-free blocked levels.
+    """
+    if hasattr(patch, "patches") and len(patch.patches) != 1:
+        raise ValueError("projected explosion requires one rectangle")
+    ny, nx = patch.coarse_shape
+    donor = patch.ghost_coarse.astype(np.intp)
+    legal = np.zeros(ny * nx, bool)
+    legal[donor] = True
+    legal &= layout.owned.ravel() & ~patch.covered.ravel()
+    if not legal[donor].all():
+        raise ValueError("projected explosion donors must be owned uncovered cells")
+    y, x = np.divmod(donor, nx)
+    neighbors = np.repeat(donor[:, None], 4, axis=1)
+    for axis, (dx, dy) in enumerate(((1, 0), (0, 1))):
+        minus = ((y - dy) % ny) * nx + (x - dx) % nx
+        plus = ((y + dy) % ny) * nx + (x + dx) % nx
+        valid = legal[minus] & legal[plus]
+        neighbors[valid, 2 * axis] = minus[valid]
+        neighbors[valid, 2 * axis + 1] = plus[valid]
+    fy, fx = np.divmod(patch.ghost_fine, patch.fine_shape[1])
+    offsets = np.column_stack((fx % 2, fy % 2)).astype(np.float32) * .5 - .25
+    return (np.ascontiguousarray(layout.owner.ravel()[neighbors], dtype=np.uint64),
+            np.ascontiguousarray(offsets))
+
+
 # Required device columns plus optional actuation; observation arms are uploaded
 # separately only when enabled. body_id/aperture remain host reduction metadata.
 BLOCKED_BOUNDARY_DEVICE_COLUMNS = ("block", "cell", "dir", "q", "fallback", "surface_velocity")

tests/unit/test_projected_explosion.py

@@ -0,0 +1,204 @@
+"""CPU stencil and actual closed-interface packet budgets; no CUDA claims."""
+import numpy as np
+import unittest
+from pathlib import Path
+import shutil
+import subprocess
+import tempfile
+
+from CelerisLab.config import LBMConfig
+from CelerisLab.lbm.blocks import BlockLayout2D
+from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY
+from CelerisLab.lbm.refinement import FixedPatch2D, FixedPatchSet2D, ACTIVE, OUTER_GHOST
+from CelerisLab.lbm.refinement_tables import projected_explosion_tables
+
+E = np.array([D2Q9_EX, D2Q9_EY]).T
+W = np.array([4/9] + [1/9]*4 + [1/36]*4)
+
+
+def setup(union=False, bounds=(4, 4, 12, 12)):
+    cfg = LBMConfig(nx=16, ny=16, viscosity=.1)
+    p = (FixedPatchSet2D.from_config(cfg, [bounds]) if union else
+         FixedPatch2D.from_config(cfg, bounds))
+    layout = BlockLayout2D.build(16, 16, 4, 4, periodic=(True, True))
+    neighbors, offsets = projected_explosion_tables(p, layout)
+    inverse = np.full(layout.nblocks * np.prod(layout.local_shape), -1, dtype=int)
+    inverse[layout.owner.ravel()] = np.arange(256)
+    return p, inverse[neighbors], offsets
+
+
+def reconstruct(c, p, neighbors, offsets, project=True):
+    # Independent NumPy evaluation of the documented stencil; neighbors have
+    # been decoded from actual blocked owner IDs, not recomputed geometry.
+    values = c.reshape(9, -1)[:, neighbors]
+    g = .5 * (values[:, :, 1::2] - values[:, :, ::2])
+    if project:
+        g[1:] -= (np.einsum('qra,qa->qr', g[1:], E[1:]) /
+                  np.sum(E[1:]**2, axis=1)[:, None])[:, :, None] * E[1:, None, :]
+    out = c.reshape(9, -1)[:, p.ghost_coarse] + np.einsum('qra,ra->qr', g, offsets)
+    return out, g
+
+
+def check_legal_footprint_mean_and_characteristic(union, bounds):
+    p, neighbors, offsets = setup(union, bounds)
+    assert set(neighbors.ravel()) <= set(p.ghost_coarse)
+    assert not p.covered.ravel()[neighbors].any()
+    assert (neighbors >= 0).all()
+    c = np.random.default_rng(72).normal(size=(9,16,16))
+    c[:, p.covered] = np.nan  # poison every forbidden donor
+    out, g = reconstruct(c, p, neighbors, offsets)
+    assert np.isfinite(out).all()
+    np.testing.assert_allclose(np.einsum('qra,qa->qr', g, E), 0, atol=5e-16)
+    for d in np.unique(p.ghost_coarse):
+        rows = p.ghost_coarse == d
+        assert rows.sum() == 4
+        np.testing.assert_allclose(out[:, rows].mean(axis=1), c.reshape(9,-1)[:,d], atol=5e-16)
+        for q in range(1,9):
+            delta = offsets[rows, None, :] - offsets[None, rows, :]
+            # Parallel displacements have zero 2D cross product.
+            parallel = np.isclose(delta[...,0]*E[q,1]-delta[...,1]*E[q,0], 0)
+            diff = out[q, rows][:,None]-out[q, rows][None,:]
+            np.testing.assert_allclose(diff[parallel], 0, atol=5e-16)
+    const = np.broadcast_to(W[:,None,None], (9,16,16)).copy()
+    constant, _ = reconstruct(const,p,neighbors,offsets)
+    np.testing.assert_array_equal(constant, np.broadcast_to(W[:,None], constant.shape))
+
+
+def collide(f, tau):
+    rho = f.sum(axis=0)
+    u = np.einsum('qa,q...->a...', E, f) / rho
+    cu = np.einsum('qa,a...->q...', E, u)
+    eq = W[:,None] * rho * (1 + 3*cu + 4.5*cu**2 - 1.5*np.sum(u*u,axis=0))
+    return f + (eq-f)/tau
+
+
+def totals(c, f, p):
+    q = c[:,~p.covered].sum(axis=1)+.25*f[:,p.fine_kind==ACTIVE].sum(axis=1)
+    return np.r_[q.sum(), E.T@q]
+
+
+def check_full_two_hop_closed_corner_budget(union, collision, bounds):
+    p, neighbors, offsets = setup(union,bounds)
+    rng=np.random.default_rng(17)
+    c=W[:,None,None]*(1+rng.normal(0,.002,(9,16,16)))
+    f=W[:,None,None]*(1+rng.normal(0,.002,(9,*p.fine_shape)))
+    start=totals(c,f,p)
+    # Count exactly the eight coarse diagonal links whose source is also
+    # uncovered, but whose midpoint lies on the closed rectangular interface.
+    q, k=np.divmod(p.coalesce_coarse.astype(np.int64),256)
+    y,x=np.divmod(k,16)
+    grazing=~p.covered[y-E[q,1],x-E[q,0]]
+    assert grazing.sum()==8
+    c[:,p.covered]=np.nan
+    for _ in range(8):
+        f.reshape(9,-1)[:,p.ghost_fine]=reconstruct(c,p,neighbors,offsets)[0]
+        for _ in range(2):
+            f=np.array([np.roll(f[q],(ey,ex),(0,1)) for q,(ex,ey) in enumerate(E)])
+            f[:,p.fine_kind==OUTER_GHOST]=np.nan
+            if collision:
+                f[:,p.fine_kind==ACTIVE]=collide(f[:,p.fine_kind==ACTIVE],p.tau_fine)
+        out=np.array([np.roll(c[q],(ey,ex),(0,1)) for q,(ex,ey) in enumerate(E)])
+        packets=f.reshape(-1)[p.coalesce_fine].mean(axis=1)
+        assert np.isfinite(packets[grazing]).all()
+        out.reshape(-1)[p.coalesce_coarse]=packets
+        if collision:out[:,~p.covered]=collide(out[:,~p.covered],p.tau_coarse)
+        out[:,p.covered]=np.nan
+        c=out
+        np.testing.assert_allclose(totals(c,f,p),start,rtol=0,atol=2e-12)
+
+
+def check_projection_is_distinguished_by_characteristic_not_global_mass():
+    p,n,o=setup()
+    c=np.random.default_rng(8).normal(size=(9,16,16))
+    _,g=reconstruct(c,p,n,o,project=False)
+    assert np.max(abs(np.einsum('qra,qa->qr',g,E))) > .1
+
+
+def check_multiple_rectangles_rejected():
+    cfg=LBMConfig(nx=24,ny=16,viscosity=.1)
+    p=FixedPatchSet2D.from_config(cfg, [(4,4,8,12),(12,4,16,12)])
+    layout=BlockLayout2D.build(16,24,4,4,periodic=(True,True))
+    with unittest.TestCase().assertRaisesRegex(ValueError, 'one rectangle'):
+        projected_explosion_tables(p,layout)
+
+
+class ProjectedExplosionTests(unittest.TestCase):
+    def test_stencil(self):
+        for union in (False, True):
+            for bounds in ((4,4,12,12),(2,3,9,10),(2,2,3,3)):
+                with self.subTest(union=union,bounds=bounds):
+                    check_legal_footprint_mean_and_characteristic(union,bounds)
+
+    def test_closed_interface(self):
+        for union in (False,True):
+            for collision in (False,True):
+                for bounds in ((4,4,12,12),(2,3,14,12),(2,2,3,3)):
+                    with self.subTest(union=union,collision=collision,bounds=bounds):
+                        check_full_two_hop_closed_corner_budget(union,collision,bounds)
+
+    def test_projection(self):
+        check_projection_is_distinguished_by_characteristic_not_global_mass()
+
+    def test_affine_and_shift_invariance(self):
+        p,n,o=setup()
+        y,x=np.indices((16,16))
+        c=np.array([.7+q*.01+x*.03+y*.05 for q in range(9)])
+        out,g=reconstruct(c,p,n,o)
+        available=n[:,::2]!=n[:,1::2]
+        gradient=np.array([.03,.05])*available
+        expected=np.broadcast_to(gradient,(9,*gradient.shape)).copy()
+        expected[1:]-=(np.einsum('qra,qa->qr',expected[1:],E[1:])/
+                       np.sum(E[1:]**2,axis=1)[:,None])[:,:,None]*E[1:,None,:]
+        np.testing.assert_allclose(g,expected,atol=2e-16)
+        shifted,_=reconstruct(c-W[:,None,None],p,n,o)
+        np.testing.assert_allclose(out-W[:,None],shifted,atol=5e-16)
+
+    def test_reject(self):
+        check_multiple_rectangles_rejected()
+
+
+class ProjectedKernelTests(unittest.TestCase):
+    @unittest.skipUnless(shutil.which("g++"), "g++ is required for actual-kernel CPU arithmetic")
+    def test_actual_kernel_fp32_shifted_and_unshifted(self):
+        p, neighbors, offsets = setup(True)
+        parent = BlockLayout2D.build(16,16,4,4,periodic=(True,True))
+        child = p.fine_layout((8,4))  # deliberately unequal parent/child strides
+        ids, offsets = projected_explosion_tables(p,parent)
+        cells = int(np.prod(parent.local_shape))
+        targets = child.owner.ravel()[p.ghost_fine]
+        donors = parent.owner.ravel()[p.ghost_coarse]
+        source = (Path(__file__).resolve().parents[2] /
+                  "src/CelerisLab/lbm/kernels/blocks.cu").read_text()
+        kernel = source.split('extern "C" __global__ void BlockExplodeProjected(',1)[1]
+        kernel = 'void BlockExplodeProjected(' + kernel.split('// Targets are receiver-local SOURCE slots',1)[0]
+        def array(name, values, kind):
+            v = np.asarray(values).ravel()
+            text = ','.join(str(int(x)) if kind=='uint64_t' else float(x).hex()+'f' for x in v)
+            return f"{kind} {name}[] = {{{text}}};\n"
+        for shifted in (False,True):
+            c=np.random.default_rng(11).normal(0,.001,(9,16,16))
+            if not shifted:c+=W[:,None,None]
+            c=c.astype(np.float32)
+            packed=np.full(parent.nblocks*9*cells,np.nan,np.float32)
+            for q in range(9):packed[parent.slot(parent.owner.ravel(),q)]=c[q].ravel()
+            # All halo/covered inputs are poisoned; kernel must use legal cores.
+            for q in range(9):packed[parent.slot(parent.owner.ravel()[p.covered.ravel()],q)]=np.nan
+            # Hex initializer needs a portable expression for poisoned NaNs.
+            data=array('coarse',packed,'float').replace('nanf','NAN')
+            code = '#include <cstdint>\n#include <cmath>\n#include <cstdio>\nusing fpxx=float;\n'
+            code += f'constexpr unsigned NQ=9, TOTAL_CELLS={int(np.prod(child.local_shape))};\n'
+            code += 'struct Id {unsigned x;}; Id blockIdx{0},blockDim{1},threadIdx{0};\n'
+            code += 'unsigned index_f(unsigned k,unsigned q){return q*TOTAL_CELLS+k;}\n'
+            code += 'int d_cx[]={0,1,-1,0,0,1,-1,1,-1},d_cy[]={0,0,0,1,-1,1,-1,-1,1};\n'+kernel
+            code += data+array('target',targets,'uint64_t')+array('donor',donors,'uint64_t')
+            code += array('neighbor',ids,'uint64_t')+array('offset',offsets,'float')
+            code += f'float fine[{child.nblocks*9*int(np.prod(child.local_shape))}]={{}};\n'
+            code += f'int main(){{for(unsigned r=0;r<{len(targets)};++r){{threadIdx.x=r;BlockExplodeProjected(coarse,fine,target,donor,{len(targets)},{cells},neighbor,offset);}}'
+            code += f'for(unsigned q=0;q<9;++q)for(unsigned r=0;r<{len(targets)};++r){{auto t=target[r];printf("%.9g\\n",fine[(t/TOTAL_CELLS)*9*TOTAL_CELLS+q*TOTAL_CELLS+t%TOTAL_CELLS]);}}}}'
+            with tempfile.TemporaryDirectory() as folder:
+                src=Path(folder)/'kernel.cpp'; binary=Path(folder)/'kernel'
+                src.write_text(code)
+                subprocess.run(['g++','-std=c++17','-O2',str(src),'-o',str(binary)],check=True,capture_output=True)
+                result=np.fromstring(subprocess.check_output([str(binary)],text=True),sep=' ').reshape(9,-1)
+            expected,_=reconstruct(c.astype(float),p,neighbors,offsets)
+            np.testing.assert_allclose(result,expected,atol=6e-8,rtol=2e-7)

tests/validation/taylor_green/run_case.py

@@ -1,4 +1,4 @@
-"""One bounded GPU case. No production edits, plots, checkpoints or field dumps."""
+"""One bounded GPU case, with explicit opt-in experimental explosion."""
 
 import argparse
 import hashlib
@@ -75,7 +75,10 @@ def main():
     p.add_argument("--n", type=int, choices=[32, 64, 128, 256], required=True)
     p.add_argument("--depth", type=int, choices=[1, 2, 3], required=True)
     p.add_argument("--output", type=Path, required=True)
+    p.add_argument("--explosion", choices=["constant", "projected_linear"], default="constant")
     args = p.parse_args()
+    if args.explosion != "constant" and args.depth != 2:
+        p.error("projected_linear requires depth 2")
     if not os.environ.get("CUDA_VISIBLE_DEVICES"):
         p.error("Explicit physical CUDA_VISIBLE_DEVICES selection required")
     repo = args.repo.resolve()
@@ -99,6 +102,8 @@ def main():
     report = dict(
         schema=1,
         source=provenance(repo),
+        historical_baseline_revision="9805383f9c50d867e8555471baf5dbda8333a652",
+        explosion=args.explosion,
         n=args.n,
         re_L=args.re,
         depth=args.depth,
@@ -132,6 +137,7 @@ def main():
                 coarse_boundary="periodic",
                 refinements=refinements,
                 level_blocks=[(8, 8)] * args.depth,
+                explosion=args.explosion,
             )
         with run:
             run.initialize(