From a135252ae78b47b57d7e77d19001700b3bedac8b Mon Sep 17 00:00:00 2001 From: Codex Date: Sun, 18 Oct 2026 14:54:17 -0700 Subject: [PATCH] feat: add experimental projected explosion per refinement pair --- docs/development/projected_pair_explosion.md | 33 ++ src/CelerisLab/lbm/kernels/blocks.cu | 33 ++ src/CelerisLab/lbm/refinement_stepper.py | 29 +- src/CelerisLab/lbm/refinement_tables.py | 44 +++ tests/unit/test_projected_explosion.py | 394 +++++++++++++++++++ tests/unit/test_projected_pair_stepper.py | 198 ++++++++++ 6 files changed, 728 insertions(+), 3 deletions(-) create mode 100644 docs/development/projected_pair_explosion.md create mode 100644 tests/unit/test_projected_explosion.py create mode 100644 tests/unit/test_projected_pair_stepper.py diff --git a/docs/development/projected_pair_explosion.md b/docs/development/projected_pair_explosion.md new file mode 100644 index 0000000..b6a1105 --- /dev/null +++ b/docs/development/projected_pair_explosion.md @@ -0,0 +1,33 @@ +# Experimental projected G explosion + +`FixedPatchStepper(..., explosion="projected_linear")` changes only the existing +G-to-G explosion reconstruction. The omitted/default `"constant"` selector retains +the same kernel, launch arguments and device allocations. Projected mode adds two +compact tables per adjacent pair; recursion, collision, BFL, restriction, buffer +ownership and synchronized publication are unchanged. `remap_child` preserves it. + +Runtime eligibility is periodic-root, fully blocked D2Q9/SRT/FP32, body/Jet-free, +without LES; the fixed-patch backend has no forcing interface. Shifted and unshifted +storage use the same reconstruction. One generic builder handles scalar windows, +unique union rows and all adjacent hierarchy pairs without depth/rectangle caps. +Donors are original explosion donors intersected with owned, uncovered ACTIVE +cells. An axis uses centered differences only if both neighbors are legal; +otherwise both addresses repeat the center. Nested layout periodicity is explicit. +No covered cells, holes, inactive cells or stale halo copies enter the stencil. + +The 33-line kernel is unchanged from our `519473a` experiment and retains its +mathematical attribution. Host tests cover union decomposition, concavity, narrow +gaps, grazing packet paths, offsets, sparse/nested owner addresses, three/four-level +source timing, unchanged constant dispatch, remapping and shifted/unshifted actual +kernel arithmetic compiled as host C++. These are not CUDA or flow validation. + +The older `519473a` / `11901a7` N128 result passed its velocity threshold narrowly +(0.001995299 versus 0.002, 0.235% margin); the matched constant velocity result was +0.004418077. N64 projected failures and the false aggregate result remain relevant. +These observations justify only the next qualification: bounded GPU default-path +equivalence, union and multilevel checks on this port, then matched physics. +No convergence, long-time, body, boundary, MRT or general-accuracy claim is made. + +The older LES flow runner pins the entire `blocks.cu` file. Its three pin checks +correctly reject this added-kernel revision; their hashes are deliberately unchanged. +A matching, reviewed harness revision is required before using that evidence path. diff --git a/src/CelerisLab/lbm/kernels/blocks.cu b/src/CelerisLab/lbm/kernels/blocks.cu index 746a4b8..8e3764f 100644 --- a/src/CelerisLab/lbm/kernels/blocks.cu +++ b/src/CelerisLab/lbm/kernels/blocks.cu @@ -63,6 +63,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( diff --git a/src/CelerisLab/lbm/refinement_stepper.py b/src/CelerisLab/lbm/refinement_stepper.py index 883d59d..2cf63ab 100644 --- a/src/CelerisLab/lbm/refinement_stepper.py +++ b/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,8 @@ class _Pair: n_ghosts: int n_rows: int parent_cells: int # parent block window cells (explode donor stride) + neighbors: object = None # projected-only owner-core cells; constant allocates neither + offsets: object = None @dataclass @@ -206,6 +208,12 @@ class FixedPatchStepper: halo -> BFL -> step four times per base step. Runtime child Jet changes require ``set_jet_amplitude(index, amplitude, level=2)``. + ``explosion="projected_linear"`` selects experimental transverse linear + reconstruction on each adjacent pair. Default "constant" is unchanged. + Projected mode requires periodic root, all-blocked SRT levels without + bodies/Jets or LES; this backend supplies no forcing. Geometry is generic, + but union/multilevel GPU and matched-physics qualification remain pending. + Level lists (any depth): ``refinements[l]`` gives level l+1 in level-l integer cells (scalar window or rectangle list, see lbm.hierarchy); ``level_blocks``/``level_objects``/``level_jets`` hold one entry per level, @@ -219,7 +227,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_window_exchange=False): + observe_raw_exchange=False, observe_window_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"): @@ -248,6 +256,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 ( + coarse_boundary != "periodic" or any(b is None for b in level_blocks) + or any(level_objects) or any(level_jets) + or cfg.collision != "SRT" or cfg.les_enabled): + raise ValueError("projected explosion requires periodic all-blocked SRT levels " + "without bodies/Jets or LES") + 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") @@ -479,6 +496,11 @@ 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, self.hierarchy.kinds[index]) + result.explode = child.module.get_function("BlockExplodeProjected") + result.neighbors, result.offsets = 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: @@ -653,7 +675,7 @@ class FixedPatchStepper: child_rectangles_in_parent=new.child.bounds, child_objects=child_objects, child_jets=jets, observe_raw_exchange=self._observe_raw_exchange, - observe_window_exchange=self._observe_window) + observe_window_exchange=self._observe_window, explosion=self._explosion) try: if self._ddf_shifting: candidate._initialize_raw(_StoredPopulations(transferred.g, transferred.shifted)) @@ -746,6 +768,7 @@ class FixedPatchStepper: 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), + *((pair.neighbors, pair.offsets) if pair.neighbors is not None else ()), 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) diff --git a/src/CelerisLab/lbm/refinement_tables.py b/src/CelerisLab/lbm/refinement_tables.py index 96423c4..b2dd7eb 100644 --- a/src/CelerisLab/lbm/refinement_tables.py +++ b/src/CelerisLab/lbm/refinement_tables.py @@ -66,6 +66,50 @@ def blocked_base_tables(patch, layout, order, wall_aware=False, *, composed_kind blocked_kind(layout, dense_kind)) +def projected_explosion_tables(pair, parent_layout, composed_parent_kind): + """Experimental transverse linear explosion, in the pair's existing row order. + + Schornbaum & Rüde (2016), arXiv:1508.07982, Algorithm 2, following + Chen et al., doi:10.1016/j.physa.2005.09.036. Independently implemented + from the mathematics: centered gradients, perpendicular projection and + child offsets +/-1/4. The kernel applies the projection per population. + + Only original explosion donors that are owned, uncovered ACTIVE cells + may supply gradients. Both axis neighbors must be legal; otherwise repeat + the center twice for zero slope. Owner-core addresses avoid stale halos. + The same rule applies to scalar windows, unions and every adjacent pair; + nested layouts retain their own per-axis (usually nonperiodic) topology. + """ + ny, nx = pair.coarse_shape + kind = np.asarray(composed_parent_kind) + if parent_layout.shape != pair.coarse_shape or kind.shape != pair.coarse_shape: + raise ValueError("projected explosion parent layout/kind has the wrong shape") + donor = pair.ghost_coarse.astype(np.intp) + legal = np.zeros(ny * nx, bool) + legal[donor] = True + legal &= (parent_layout.owned.ravel() & ~pair.covered.ravel() + & (kind.ravel() == ACTIVE)) + if not legal[donor].all(): + raise ValueError("projected explosion donors need owned uncovered parent ACTIVE 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] + if not parent_layout.periodic[axis]: + coordinate, extent = (x, nx) if axis == 0 else (y, ny) + valid &= (coordinate > 0) & (coordinate < extent - 1) + neighbors[valid, 2 * axis] = minus[valid] + neighbors[valid, 2 * axis + 1] = plus[valid] + fy, fx = np.divmod(pair.ghost_fine, pair.fine_shape[1]) + ox, oy = pair.fine_origin + offsets = np.column_stack((ox + (fx + .5) / 2 - (x + .5), + oy + (fy + .5) / 2 - (y + .5))) + return (np.ascontiguousarray(parent_layout.owner.ravel()[neighbors], dtype=np.uint64), + np.ascontiguousarray(offsets, dtype=np.float32)) + + # 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") diff --git a/tests/unit/test_projected_explosion.py b/tests/unit/test_projected_explosion.py new file mode 100644 index 0000000..3c16634 --- /dev/null +++ b/tests/unit/test_projected_explosion.py @@ -0,0 +1,394 @@ +"""CPU stencil and actual closed-interface packet budgets; no CUDA claims.""" +import numpy as np +import unittest +from pathlib import Path +from dataclasses import replace +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 +from CelerisLab.lbm.hierarchy import build_hierarchy +from tests.unit.test_hierarchy import CFG, FOUR + +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, np.where(p.covered, OUTER_GHOST, ACTIVE)) + 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 + + + + +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) + + +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,np.where(p.covered,OUTER_GHOST,ACTIVE)) + 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 \n#include \n#include \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) + + +UNIONS = [ + [(8, 8, 24, 20)], + [(8, 8, 16, 20), (16, 8, 24, 20)], # abutting + [(8, 8, 18, 20), (14, 8, 24, 20)], # overlap + [(8, 8, 16, 20), (17, 8, 24, 20)], # gap 1 + [(8, 8, 16, 20), (18, 8, 24, 20)], # gap 2 + [(8, 8, 16, 20), (19, 8, 24, 20)], # gap 3 + [(8, 8, 20, 14), (8, 14, 14, 22)], # L concavity + [(4, 4, 10, 10), (18, 16, 26, 26)], # separated + [(4, 4, 5, 7), (4, 5, 7, 7)], # known concave exit/re-entry path +] + + +def decode(layout, neighbors): + inverse = np.full(layout.hole, -1, dtype=int) + inverse[layout.owner[layout.owned]] = np.flatnonzero(layout.owned) + return inverse[neighbors] + + +def check_pair_stencil(pair, layout, kind): + ids, offsets = projected_explosion_tables(pair, layout, kind) + neighbors = decode(layout, ids) + legal = set(map(int, pair.ghost_coarse)) + assert len(np.unique(pair.ghost_fine)) == len(pair.ghost_fine) + assert set(neighbors.ravel()) <= legal + assert layout.owned.ravel()[neighbors].all() + assert (kind.ravel()[neighbors] == ACTIVE).all() + assert not pair.covered.ravel()[neighbors].any() + ny, nx = pair.coarse_shape + for row, donor in enumerate(pair.ghost_coarse): + y, x = divmod(int(donor), nx) + for axis, (dx, dy) in enumerate(((1, 0), (0, 1))): + ends = [(x - dx, y - dy), (x + dx, y + dy)] + expected = [] + for xx, yy in ends: + outside = not (0 <= xx < nx and 0 <= yy < ny) + if outside and not layout.periodic[axis]: + break + cell = (yy % ny) * nx + xx % nx + if cell not in legal: + break + expected.append(cell) + if len(expected) != 2: + expected = [int(donor)] * 2 + np.testing.assert_array_equal(neighbors[row, 2*axis:2*axis+2], expected) + fy, fx = divmod(int(pair.ghost_fine[row]), pair.fine_shape[1]) + expected = [pair.fine_origin[0]+(fx+.5)/2-(x+.5), + pair.fine_origin[1]+(fy+.5)/2-(y+.5)] + np.testing.assert_array_equal(offsets[row], expected) + assert set(offsets.ravel()) == {-.25, .25} + assert ids.dtype == np.uint64 and ids.flags.c_contiguous + assert offsets.dtype == np.float32 and offsets.flags.c_contiguous + return neighbors, offsets + + +class GenericProjectedPairs(unittest.TestCase): + def test_single_window_and_singleton_union_exact_stencils(self): + for bounds in ((4,4,12,12), (2,3,9,10), (2,2,3,3)): + a, an, ao = setup(False, bounds) + b, bn, bo = setup(True, bounds) + np.testing.assert_array_equal(a.ghost_coarse, b.ghost_coarse) + np.testing.assert_array_equal(an, bn) + np.testing.assert_array_equal(ao, bo) + # Coordinates differ between scalar and union arrays; physical rows don't. + ay, ax = np.divmod(a.ghost_fine, a.fine_shape[1]) + by, bx = np.divmod(b.ghost_fine, b.fine_shape[1]) + np.testing.assert_array_equal(ax + 2*a.fine_origin[0], bx) + np.testing.assert_array_equal(ay + 2*a.fine_origin[1], by) + + def test_union_decomposition_permutation_and_topologies(self): + cfg = LBMConfig(nx=30, ny=30, viscosity=.1) + layout = BlockLayout2D.build(30,30,6,5,periodic=(True,True)) + stencils = [] + for rectangles in UNIONS: + pair = FixedPatchSet2D.from_config(cfg, rectangles) + kind = np.where(pair.covered, OUTER_GHOST, ACTIVE) + n, o = check_pair_stencil(pair, layout, kind) + other = FixedPatchSet2D.from_config(cfg, rectangles[::-1]) + on, oo = projected_explosion_tables(other, layout, kind) + np.testing.assert_array_equal(n, decode(layout, on)) + np.testing.assert_array_equal(o, oo) + stencils.append((pair.ghost_fine, pair.ghost_coarse, n, o)) + c = np.random.default_rng(92).normal(size=(9,30,30)) + c[:, pair.covered] = np.nan + values, gradients = reconstruct(c, pair, n, o) + self.assertTrue(np.isfinite(values).all()) + np.testing.assert_allclose(np.einsum('qra,qa->qr', gradients, E),0,atol=1e-15) + # q0 keeps its full gradient; no division by zero or forced zero. + expected_q0 = .5*(c[0].ravel()[n[:,1::2]]-c[0].ravel()[n[:,::2]]) + np.testing.assert_array_equal(gradients[0], expected_q0) + self.assertGreater(np.max(abs(gradients[0])), .1) + for donor in np.unique(pair.ghost_coarse): + rows = pair.ghost_coarse == donor + self.assertEqual(rows.sum(),4) + np.testing.assert_allclose(values[:,rows].mean(axis=1),c.reshape(9,-1)[:,donor],atol=1e-15) + for q in range(1,9): + delta = o[rows,None,:]-o[None,rows,:] + parallel = delta[...,0]*E[q,1] == delta[...,1]*E[q,0] + diff = values[q,rows][:,None]-values[q,rows][None,:] + np.testing.assert_allclose(diff[parallel],0,atol=1e-15) + for split in stencils[1:3]: + for expected, actual in zip(stencils[0],split): + np.testing.assert_array_equal(actual,expected) + + def test_nested_scalar_origins_sparse_holes_and_unequal_strides(self): + cases = [(CFG, FOUR), (LBMConfig(nx=16,ny=16,viscosity=.02), + [(4,4,12,12),(8,8,16,16),(8,8,16,16)])] + for cfg, refs in cases: + hierarchy = build_hierarchy(cfg,refs) + for index,pair in enumerate(hierarchy.pairs): + if index == 0: + layout = BlockLayout2D.build(cfg.ny,cfg.nx,4,4,periodic=(True,True)) + else: + prior = hierarchy.pairs[index-1] + layout = (prior.fine_layout((4,8)) if isinstance(prior,FixedPatchSet2D) + else BlockLayout2D.build(*prior.fine_shape,4,8,periodic=(False,False))) + kind = hierarchy.kinds[index] + n,o = check_pair_stencil(pair,layout,kind) + c = np.random.default_rng(index).normal(size=(9,*pair.coarse_shape)) + c[:,~layout.owned | (kind != ACTIVE)] = np.nan + self.assertTrue(np.isfinite(reconstruct(c,pair,n,o)[0]).all()) + if index and isinstance(pair,FixedPatchSet2D): + self.assertTrue((~layout.owned).any()) + self.assertEqual(layout.periodic, (index == 0, index == 0)) + + def test_periodicity_is_per_axis_and_missing_pairs_repeat_center(self): + cfg=LBMConfig(nx=16,ny=16,viscosity=.1) + p=FixedPatch2D.from_config(cfg,(2,2,14,14)) + kind=np.where(p.covered,OUTER_GHOST,ACTIVE) + results=[] + for periodic in ((False,False),(True,False),(False,True),(True,True)): + layout=BlockLayout2D.build(16,16,4,4,periodic=periodic) + n,_=check_pair_stencil(p,layout,kind) + results.append(n) + y,x=np.divmod(p.ghost_coarse,16) + for axis,coordinate in enumerate((x,y)): + edge=(coordinate==0)|(coordinate==15) + if not periodic[axis]: + np.testing.assert_array_equal(n[edge,2*axis:2*axis+2], + np.repeat(p.ghost_coarse[edge,None],2,axis=1)) + self.assertTrue((results[0]!=results[1]).any()) + self.assertTrue((results[0]!=results[2]).any()) + self.assertTrue((results[1]!=results[3]).any()) + + def test_invalid_parent_donors_reject_before_owner_lookup(self): + p,_,_=setup() + layout=BlockLayout2D.build(16,16,4,4) + kind=np.where(p.covered,OUTER_GHOST,ACTIVE) + for nonactive in (1,2,255): + bad=kind.copy(); bad.ravel()[p.ghost_coarse[0]]=nonactive + with self.assertRaisesRegex(ValueError,'ACTIVE'): + projected_explosion_tables(p,layout,bad) + with self.assertRaisesRegex(ValueError,'shape'): + projected_explosion_tables(p,layout,kind[:1]) + hole_layout=BlockLayout2D.build(16,16,4,4,origins=[(0,0)]) + with self.assertRaisesRegex(ValueError,'owned'): + projected_explosion_tables(p,hole_layout,kind) + + def test_generic_two_hop_mass_and_momentum_budget(self): + cfg=LBMConfig(nx=30,ny=30,viscosity=.1) + layout=BlockLayout2D.build(30,30,6,5) + grazing_counts=set() + for rectangles in UNIONS: + pair=FixedPatchSet2D.from_config(cfg,rectangles) + ids,o=projected_explosion_tables(pair,layout,np.where(pair.covered,OUTER_GHOST,ACTIVE)) + n=decode(layout,ids) + q,k=np.divmod(pair.coalesce_coarse.astype(np.int64),900) + y,x=np.divmod(k,30) + grazing=~pair.covered[(y-E[q,1])%30,(x-E[q,0])%30] + grazing_counts.add(int(grazing.sum())) + self.assertTrue(grazing.any()) + for collision in (False,True): + rng=np.random.default_rng(17) + c=W[:,None,None]*(1+rng.normal(0,.002,(9,30,30))) + f=W[:,None,None]*(1+rng.normal(0,.002,(9,*pair.fine_shape))) + start=totals(c,f,pair) + c[:,pair.covered]=np.nan + for _ in range(4): + f.reshape(9,-1)[:,pair.ghost_fine]=reconstruct(c,pair,n,o)[0] + for _ in range(2): + f=np.array([np.roll(f[q],(ey,ex),(0,1)) for q,(ex,ey) in enumerate(E)]) + f[:,pair.fine_kind==OUTER_GHOST]=np.nan + if collision: + f[:,pair.fine_kind==ACTIVE]=collide(f[:,pair.fine_kind==ACTIVE],pair.tau_fine) + out=np.array([np.roll(c[q],(ey,ex),(0,1)) for q,(ex,ey) in enumerate(E)]) + packets=f.reshape(-1)[pair.coalesce_fine].mean(axis=1) + self.assertTrue(np.isfinite(packets[grazing]).all()) + out.reshape(-1)[pair.coalesce_coarse]=packets + if collision: + out[:,~pair.covered]=collide(out[:,~pair.covered],pair.tau_coarse) + out[:,pair.covered]=np.nan + c=out + np.testing.assert_allclose(totals(c,f,pair),start,rtol=0,atol=4e-12) + self.assertGreater(len(grazing_counts),1) # not a hard-coded eight-corner oracle diff --git a/tests/unit/test_projected_pair_stepper.py b/tests/unit/test_projected_pair_stepper.py new file mode 100644 index 0000000..4dfcbee --- /dev/null +++ b/tests/unit/test_projected_pair_stepper.py @@ -0,0 +1,198 @@ +"""Production host dispatch with CUDA doubles; no CUDA compilation or flow runs.""" +from contextlib import nullcontext +from dataclasses import replace +from unittest.mock import Mock +import unittest + +import numpy as np + +from CelerisLab.config import LBMConfig +from tests.unit.test_hierarchy import CFG, FOUR +from tests.unit.test_refinement_three_level_multipatch import host_module + + +def backend(**options): + host = host_module() + calls, allocations, storage = [], [], {} + context = Mock(activate=nullcontext, sm_arch='sm_70') + host.cuda.mem_alloc.side_effect = lambda size: allocate(size) + + def allocate(size): + pointer = Mock(name=f'buffer{len(allocations)}', size=size) + allocations.append(pointer) + return pointer + + def load(_): + module = Mock() + module.get_global.return_value = (Mock(), 20) + module.get_function.side_effect = lambda name: Mock( + name=name, side_effect=lambda *a, **kw: calls.append((name, a, kw))) + return module + + host.compiler.load_module.side_effect = load + host.cuda.memcpy_htod.side_effect = lambda ptr, a: storage.__setitem__(id(ptr), np.array(a, copy=True)) + host.cuda.memcpy_dtoh.side_effect = lambda a, ptr: np.copyto(a, storage[id(ptr)]) + cfg = options.pop('cfg', CFG) + refinements = options.pop('refinements', FOUR) + blocks = options.pop('level_blocks', [(4, 4)] * (len(refinements) + 1)) + run = host.FixedPatchStepper(cfg, context=context, refinements=refinements, + level_blocks=blocks, **options) + run.initialize(*(np.full(level.shape, 1/9, np.float32) for level in run._levels)) + return run, host, calls, allocations, storage + + +def trace(run, calls, allocations, storage): + """Stable pointer identities expose allocations/arguments, not just kernel names.""" + pointer = {id(value): index for index, value in enumerate(allocations)} + launches = [(name, tuple(('ptr', pointer[id(a)]) if id(a) in pointer else int(a) + for a in args), kw['grid'], kw['block']) + for name, args, kw in calls] + uploads = [(index, value.dtype.str, value.shape, value.tobytes()) + for index, ptr in enumerate(allocations) if (value := storage.get(id(ptr))) is not None] + return [p.size for p in allocations], uploads, launches + + +class ProjectedPairDispatch(unittest.TestCase): + def test_omitted_and_constant_have_identical_allocations_uploads_and_launches(self): + results = [] + for options in ({}, {'explosion': 'constant'}): + run, host, calls, allocations, storage = backend(**options) + self.assertEqual(run._explosion, 'constant') + self.assertTrue(all(p.neighbors is None and p.offsets is None for p in run._pairs)) + run.step() + explosion = [args for name, args, _ in calls if name == 'BlockExplode'] + self.assertEqual(len(explosion), 7) + self.assertTrue(all(len(args) == 6 for args in explosion)) + self.assertFalse(any(name == 'BlockExplodeProjected' for name, _, _ in calls)) + results.append(trace(run, calls, allocations, storage)) + run.close() + self.assertEqual(*results) + + def test_projected_adds_exactly_two_tables_per_adjacent_pair(self): + for depth in (2, 3, 4): + constant, _, _, ca, _ = backend(refinements=FOUR[:depth-1]) + projected, _, calls, pa, storage = backend( + refinements=FOUR[:depth-1], explosion='projected_linear') + self.assertEqual(len(pa)-len(ca), 2*(depth-1)) + for pair, geometry in zip(projected._pairs, projected.hierarchy.pairs): + self.assertEqual(storage[id(pair.neighbors)].shape, (len(geometry.ghost_fine), 4)) + self.assertEqual(storage[id(pair.offsets)].shape, (len(geometry.ghost_fine), 2)) + self.assertEqual(pair.neighbors.size, len(geometry.ghost_fine)*4*8) + self.assertEqual(pair.offsets.size, len(geometry.ghost_fine)*2*4) + projected.step() + explosions = [args for name, args, _ in calls if name == 'BlockExplodeProjected'] + self.assertEqual(len(explosions), 2**(depth-1)-1) + self.assertTrue(all(len(args) == 8 for args in explosions)) + self.assertFalse(any(name == 'BlockExplode' for name, _, _ in calls)) + for args in explosions: + pair = next(p for p in projected._pairs if args[2] is p.ghosts) + self.assertIs(args[6], pair.neighbors) + self.assertIs(args[7], pair.offsets) + self.assertEqual(args[5], pair.parent_cells) + constant.close(); projected.close() + + def test_three_and_four_level_schedule_has_fresh_midpoint_sources(self): + for depth in (3, 4): + for explosion in ('constant', 'projected_linear'): + run, host, calls, allocations, storage = backend( + refinements=FOUR[:depth-1], explosion=explosion) + before = len(allocations) + events = [] + for index, pair in enumerate(run._pairs): + pair.explode.side_effect = lambda source, *a, index=index, **kw: events.append(('explode', index, source)) + pair.coalesce.side_effect = lambda child, source, *a, index=index, **kw: events.append(('restrict', index, source)) + for index, level in enumerate(run._levels): + level.kernel.side_effect = lambda source, *a, index=index, **kw: events.append(('step', index, source)) + level.halo_copy.side_effect = lambda source, *a, index=index, **kw: events.append(('halo', index, source)) + expected = [] + + def hop(index, source): + if index+1 < depth: + expected.append(('explode', index, source)) + child = run._levels[index+1] + hop(index+1, child.current) + hop(index+1, child.spare) + expected.append(('halo', index, source)) + if index+1 < depth: + expected.append(('restrict', index, source)) + expected.append(('step', index, source)) + + hop(0, run._levels[0].current) + run.step() + self.assertEqual(events, expected) + self.assertEqual([level.kernel.call_count for level in run._levels], [2**i for i in range(depth)]) + self.assertEqual(len(allocations), before) + self.assertEqual(run.step_count, 1) + for index in range(1, depth-1): + feeds = [source for name, pair, source in events if name == 'explode' and pair == index] + self.assertEqual(feeds, [run._levels[index].current, run._levels[index].spare] * (2**(index-1))) + run.close() + + def test_scalar_windows_follow_the_same_pair_path(self): + run, _, _, _, _ = backend(cfg=LBMConfig(nx=16,ny=16,viscosity=.02), + refinements=[(4,4,12,12),(8,8,16,16),(8,8,16,16)], + level_blocks=[(4,4),(4,8),(8,4),(4,4)], explosion='projected_linear') + self.assertTrue(all(pair.neighbors is not None for pair in run._pairs)) + self.assertEqual([pair.parent_cells for pair in run._pairs], [36,60,60]) + run.step() + self.assertEqual([level.kernel.call_count for level in run._levels], [1,2,4,8]) + run.close() + + def test_physics_gate_precedes_context_and_allocations(self): + host = host_module() + context = Mock() + options = dict(refinements=FOUR[:1], level_blocks=[(4,4)]*2, + explosion='projected_linear', context=context) + for cfg, changes in [(CFG, {'explosion':'unknown'}), + (CFG, {'coarse_boundary':'slot_reservoir'}), + (CFG, {'coarse_boundary':'channel'}), + (CFG, {'level_blocks':[None,None], 'refinements':[(6,6,26,26)]}), + (CFG, {'level_objects':[(),(object(),)]}), + (CFG, {'level_jets':[(),(object(),)]}), + (replace(CFG,collision='MRT'), {}), + (replace(CFG,les_enabled=True), {})]: + with self.subTest(changes=changes, cfg=cfg), self.assertRaisesRegex(ValueError,'explosion'): + host.FixedPatchStepper(cfg, **dict(options, **changes)) + context.activate.assert_not_called() + host.cuda.mem_alloc.assert_not_called() + # The legacy constant path still allows MRT and dense storage. + for options in ({'cfg':replace(CFG,collision='MRT')}, + {'level_blocks':[None,None], 'refinements':[(6,6,26,26)]}): + run, _, _, _, _ = backend(**options) + run.close() + + def test_remap_preserves_selector_tables_and_constant_default(self): + for explosion in ('constant','projected_linear'): + run, host, _, _, storage = backend(refinements=FOUR[:2], explosion=explosion) + run._step_count = 7 + old = run.download() + clone = run.remap_child(FOUR[1]) + self.assertEqual(clone._explosion, explosion) + self.assertEqual(clone.step_count, 7) + for pair in clone._pairs: + self.assertEqual(pair.neighbors is None, explosion == 'constant') + for expected, actual in zip(old, clone.download()): + np.testing.assert_array_equal(actual,expected) + clone.step() + self.assertEqual((run.step_count,clone.step_count), (7,8)) + clone.close(); run.close() + + def test_changed_remap_geometry_keeps_projected_in_both_representations(self): + changed = [(22,22,42,42), (72,68,76,76)] + for shifted in (False, True): + run, _, _, _, _ = backend(cfg=replace(CFG,ddf_shifting=shifted), + refinements=FOUR[:2], explosion='projected_linear') + run._step_count = 3 + before = run.download() + clone = run.remap_child(changed) + self.assertEqual(clone._explosion, 'projected_linear') + self.assertEqual(clone.hierarchy.pairs[1].bounds, tuple(changed)) + self.assertNotEqual(clone._pairs[1].n_ghosts, run._pairs[1].n_ghosts) + for old, new in zip(run._pairs, clone._pairs): + self.assertIsNot(old.neighbors, new.neighbors) + self.assertIsNotNone(new.neighbors) + for expected, actual in zip(before, run.download()): + np.testing.assert_array_equal(actual, expected) + clone.step() + self.assertEqual((run.step_count, clone.step_count), (3, 4)) + clone.close(); run.close()