From 906e0d9a0e12b4b29c0800c1c68f904ec1cb099d Mon Sep 17 00:00:00 2001 From: Frank14f <1515444314@qq.com> Date: Thu, 15 Oct 2026 16:25:26 +0800 Subject: [PATCH] Support blocked parent storage in fixed patch stepping --- src/CelerisLab/lbm/kernels/blocks.cu | 37 +- .../lbm/kernels/boundary/planar_velocity.cuh | 20 +- src/CelerisLab/lbm/refinement_stepper.py | 173 +++++++--- .../multigrid/test_refinement_blocked_base.py | 324 ++++++++++++++++++ .../multigrid/test_refinement_blocked_fine.py | 6 +- .../multigrid/test_refinement_followup.py | 2 +- 6 files changed, 496 insertions(+), 66 deletions(-) create mode 100644 tests/validation/multigrid/test_refinement_blocked_base.py diff --git a/src/CelerisLab/lbm/kernels/blocks.cu b/src/CelerisLab/lbm/kernels/blocks.cu index c504a34..55b1849 100644 --- a/src/CelerisLab/lbm/kernels/blocks.cu +++ b/src/CelerisLab/lbm/kernels/blocks.cu @@ -16,6 +16,7 @@ #include "operators/collision_srt.cuh" #include "streaming/pull_double_buffer.cuh" #include "boundary/curved_boundary.cuh" +#include "boundary/planar_velocity.cuh" #if DIM != 2 || NQ != 9 || COLLISION_MODEL != 0 || STORE_PRECISION != 0 || USE_DDF_SHIFTING || USE_LES || STREAMING_MODEL != 0 #error "Blocks require unshifted FP32 D2Q9 SRT double buffers without LES" @@ -24,6 +25,11 @@ // Fixed-patch cell kinds, matching refinement.py ACTIVE/INNER_GHOST/OUTER_GHOST. constexpr unsigned char KIND_ACTIVE = 0u, KIND_INNER = 1u, KIND_OUTER = 2u; +// Zero-initialized for ordinary blocks/fine levels. Set only for parent walls; +// NX/NY remain local storage strides, never physical wall coordinates. +struct BlockWall { unsigned int nx, ny, slot_width; float amplitude; }; +extern "C" { __constant__ BlockWall d_block_wall; } + // Halo refresh: block cell ids b*TOTAL_CELLS+c, sources are always core cells. extern "C" __global__ void BlockHaloCopy( fpxx* f, const uint64_t* dst, const uint64_t* src, unsigned int count) @@ -38,17 +44,32 @@ extern "C" __global__ void BlockHaloCopy( } // Fixed-patch explosion into blocked fine storage: target = owner block cell -// id of a patch-ghost cell, donor = dense coarse cell (coarse is Q-major dense). +// id of a patch-ghost cell; donors are parent block cell ids. Dense parent +// storage is the one-window case with coarse_cells = full domain size. extern "C" __global__ void BlockExplode( - const fpxx* coarse, fpxx* fine, const uint64_t* target, const uint32_t* donor, + const fpxx* coarse, fpxx* fine, const uint64_t* target, const uint64_t* donor, unsigned int count, unsigned long coarse_cells) { unsigned int row = blockIdx.x * blockDim.x + threadIdx.x; if (row >= count) return; const uint64_t t = target[row]; fpxx* to = fine + (t / TOTAL_CELLS) * NQ * TOTAL_CELLS; + const uint64_t d = donor[row]; + const fpxx* from = coarse + (d / coarse_cells) * NQ * coarse_cells; for (unsigned int q = 0; q < NQ; ++q) - to[index_f(t % TOTAL_CELLS, q)] = coarse[q * coarse_cells + donor[row]]; + to[index_f(t % TOTAL_CELLS, q)] = from[q * coarse_cells + d % coarse_cells]; +} + +// 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( + const fpxx* fine, fpxx* coarse, const uint64_t* target, const uint64_t* children, + unsigned int count) +{ + unsigned int row = blockIdx.x * blockDim.x + threadIdx.x; + if (row >= count) return; + const uint64_t* s = children + 4ul * row; + coarse[target[row]] = 0.25f * (fine[s[0]] + fine[s[1]] + fine[s[2]] + fine[s[3]]); } // Static cut links, receiver-local staging: cell is the receiver's local core @@ -92,6 +113,16 @@ extern "C" __global__ void BlockStep(const fpxx* in, fpxx* out, unsigned int hal compute_neighbors(x, y, j); // Core neighbors stay inside the window: no wrap. stream_pull_load(k, f, in + base, j); if (cell_kind == KIND_ACTIVE) { + if (d_block_wall.nx != 0) { + const unsigned int bx = NX - 2 * halo, by = NY - 2 * halo; + const unsigned int nbx = d_block_wall.nx / bx; + const unsigned int gx = (blockIdx.z % nbx) * bx + x - halo; + const unsigned int gy = (blockIdx.z / nbx) * by + y - halo; + if (gy == 1) + reconstruct_bottom_velocity(f, 0.0f, bottom_wall_normal_velocity( + gx, d_block_wall.nx, d_block_wall.slot_width, d_block_wall.amplitude)); + else if (gy == d_block_wall.ny - 2) reconstruct_top_reservoir(f, RHO); + } float rho, ux, uy, feq[NQ]; compute_rho_u(f, rho, ux, uy); compute_feq(rho, ux, uy, feq); diff --git a/src/CelerisLab/lbm/kernels/boundary/planar_velocity.cuh b/src/CelerisLab/lbm/kernels/boundary/planar_velocity.cuh index 87db746..820950d 100644 --- a/src/CelerisLab/lbm/kernels/boundary/planar_velocity.cuh +++ b/src/CelerisLab/lbm/kernels/boundary/planar_velocity.cuh @@ -3,16 +3,26 @@ // the caller owns collision/store. No force, shifting or off-lattice walls. #ifndef CELERIS_BOUNDARY_PLANAR_VELOCITY_CUH #define CELERIS_BOUNDARY_PLANAR_VELOCITY_CUH +#include + +__device__ __forceinline__ float bottom_wall_normal_velocity( + unsigned int x, unsigned int width, unsigned int slot_width, float amplitude) +{ + // Width is the physical domain width, even when NX is a block window. + if (slot_width > 0 && slot_width < width) { + const float s = (2.0f*(float)x-(float)(width-1))/(float)slot_width; + const float w = fmaxf(0.0f, 1.0f-s*s); + return amplitude*w*w; + } + return amplitude; // Disabled or explicitly uniform full width. +} __device__ __forceinline__ float bottom_wall_normal_velocity(unsigned int x) { #if WALL_SLOT_WIDTH > 0 && WALL_SLOT_WIDTH < NX - // C1 compact profile centered at (NX-1)/2; amplitude is its continuous peak. - const float s = (2.0f*(float)x-(float)(NX-1))/(float)WALL_SLOT_WIDTH; - const float w = fmaxf(0.0f, 1.0f-s*s); - return WALL_NORMAL_VELOCITY*w*w; + return bottom_wall_normal_velocity(x, NX, WALL_SLOT_WIDTH, WALL_NORMAL_VELOCITY); #else - return WALL_NORMAL_VELOCITY; // Disabled or explicitly uniform full width. + return WALL_NORMAL_VELOCITY; #endif } diff --git a/src/CelerisLab/lbm/refinement_stepper.py b/src/CelerisLab/lbm/refinement_stepper.py index cb649a6..a63be11 100644 --- a/src/CelerisLab/lbm/refinement_stepper.py +++ b/src/CelerisLab/lbm/refinement_stepper.py @@ -18,7 +18,7 @@ from .blocks import BlockLayout2D, blocked_link_tables from .curved_links import surface_link_velocity from .block_stepper import MAX_GRID_DIM from .refinement import ACTIVE, HALO, OUTER_GHOST, FixedPatch2D, prepare_fine_boundary -from .descriptors import SOLID +from .descriptors import SOLID, D2Q9_EX, D2Q9_EY def blocked_fine_tables(patch, layout, order): @@ -37,6 +37,27 @@ def blocked_fine_tables(patch, layout, order): layout.slot(owner[cell.astype(np.intp)], q), kind) +def blocked_base_tables(patch, layout, order, wall_aware=False): + """Parent donor cells and receiver-local arrival slots, in transfer row order. + + Stage arrivals only after child reads and parent halo gather, immediately + before parent pull. Source *owner* slots are wrong at block seams. + """ + owner = layout.owner.ravel() + q, receiver = np.divmod(patch.coalesce_coarse[order], np.uint64(owner.size)) + offset = (np.asarray(D2Q9_EY)[q.astype(np.intp)] * layout.local_shape[1] + + np.asarray(D2Q9_EX)[q.astype(np.intp)]) + source = owner[receiver.astype(np.intp)].astype(np.int64) - offset + kind = np.full(layout.nblocks * int(np.prod(layout.local_shape)), OUTER_GHOST, np.uint8) + leaves = np.ones(patch.coarse_shape, bool) + x0, y0, x1, y1 = patch.bounds + leaves[y0:y1, x0:x1] = False + if wall_aware: + leaves[[0, -1]] = False + kind[owner[leaves.ravel()].astype(np.intp)] = ACTIVE + return (owner[patch.ghost_coarse].astype(np.uint64), layout.slot(source, q), kind) + + @dataclass class _Level: shape: tuple @@ -45,6 +66,11 @@ class _Level: spare: object kernel: object grid: tuple + layout: object = None + kind: object = None + halo_copy: object = None + halo_dst: object = None + halo_src: object = None class FixedPatchStepper: @@ -75,11 +101,14 @@ class FixedPatchStepper: Blocked fine_objects are stationary circles, optionally surface-actuated: prepare_fine_boundary decides legality on the dense fine grid, then cut links are staged receiver-local per hop: halo -> BFL -> step, same buffer. + ``base_block=(by, bx)`` also blocks the parent; it requires ``fine_block``. + Parent halos precede receiver-local coalescence staging and pull/collision; + covered cells are explicitly poisoned on every step. Public shapes stay dense. Blocked mode applies prescribed surface velocity only; its per-link exchange diagnostic is unavailable (rows keep physical cut-link order). """ def __init__(self, cfg, bounds, *, context, coarse_boundary="periodic", fine_objects=(), - fine_block=None): + fine_block=None, base_block=None): if HALO != 4: raise RuntimeError("refinement.cu requires a four-cell fine halo") if coarse_boundary not in ("periodic", "slot_reservoir"): @@ -94,6 +123,17 @@ class FixedPatchStepper: if self.patch.bounds[1] < 4 or self.patch.bounds[3] > cfg.ny - 4: raise ValueError("slot_reservoir requires y0 >= 4 and y1 <= ny-4") fine_objects = tuple(fine_objects) + if base_block is not None and fine_block is None: + raise ValueError("base_block requires fine_block") + self.base_layout = None + if base_block is not None: + if len(base_block) != 2: + raise ValueError("base_block must be (by, bx)") + self.base_layout = BlockLayout2D.build( + cfg.ny, cfg.nx, *base_block, periodic=(True, not self._wall_aware)) + if (self.base_layout.nblocks > MAX_GRID_DIM + or self.base_layout.local_shape[0] > MAX_GRID_DIM): + raise ValueError("base blocks exceed launch limits") self.fine_layout = None if fine_block is not None: if len(fine_block) != 2: @@ -121,7 +161,7 @@ class FixedPatchStepper: self._levels = [] self._stream = None self._fine_flags = self._fine_curved = self._curved = None - self._kind = self._halo_copy = self._halo_dst = self._halo_src = None + self._coalesce = None self._link_tables = {} self._block = (cfg.threads_per_block, 1, 1) fine_cfg = replace(cfg, nx=self.patch.fine_shape[1], ny=self.patch.fine_shape[0], @@ -141,14 +181,15 @@ class FixedPatchStepper: with self.ctx.activate(): try: self._stream = cuda.Stream() - layout = self.fine_layout - fine = ((fine_cfg, "PatchFineStep", "refinement.cu") if layout is None else - (replace(fine_cfg, nx=layout.local_shape[1], ny=layout.local_shape[0], - y_wall_bc="bounce_back", top_wall_bc="velocity", wall_slot_width=0, - wall_normal_velocity=0.0, velocity=0.0), "BlockStep", "blocks.cu")) - for config, name, tau, source in ( - (replace(cfg), "PatchCoarseStep", self.patch.tau_coarse, "refinement.cu"), - (*fine[:2], self.patch.tau_fine, fine[2])): + for index, (config, layout, tau, name) in enumerate(( + (replace(cfg), self.base_layout, self.patch.tau_coarse, "PatchCoarseStep"), + (fine_cfg, self.fine_layout, self.patch.tau_fine, "PatchFineStep"))): + source = "refinement.cu" + if layout is not None: + config = replace(config, nx=layout.local_shape[1], ny=layout.local_shape[0], + y_wall_bc="bounce_back", top_wall_bc="velocity", + wall_slot_width=0, wall_normal_velocity=0.0, velocity=0.0) + name, source = "BlockStep", "blocks.cu" config.validate() if not np.isclose(config.omega, 1 / tau, rtol=1e-14, atol=0): raise RuntimeError("level omega/tau mismatch") @@ -165,33 +206,45 @@ class FixedPatchStepper: if size != struct.calcsize("=ffffI"): raise RuntimeError("LBMParams host/device layout mismatch") cuda.memcpy_htod(ptr, struct.pack("=ffffI", config.omega, 0, 0, 0, n_objects)) - shape = (9, config.ny, config.nx) - size = int(np.prod(shape)) * 4 - grid = ((config.nx + self._block[0] - 1) // self._block[0], config.ny, 1) - if name == "BlockStep": # public dense shape, private blocked allocation - shape = (9, *self.patch.fine_shape) - size = int(np.prod(layout.blocked_shape)) * 4 - grid = grid[:2] + (layout.nblocks,) - self._levels.append(_Level(shape, module, self._allocate(size), - self._allocate(size), module.get_function(name), grid)) + shape = (9, *(self.patch.coarse_shape if index == 0 else self.patch.fine_shape)) + size = int(np.prod(shape if layout is None else layout.blocked_shape)) * 4 + grid = ((config.nx + self._block[0] - 1) // self._block[0], config.ny, + 1 if layout is None else layout.nblocks) + level = _Level(shape, module, self._allocate(size), self._allocate(size), + module.get_function(name), grid, layout) + self._levels.append(level) + if layout is not None: + level.halo_copy = module.get_function("BlockHaloCopy") + if len(layout.halo_dst): + level.halo_dst = self._upload(layout.halo_dst) + level.halo_src = self._upload(layout.halo_src) + if index == 0 and self._wall_aware: + ptr, size = module.get_global("d_block_wall") + if size != struct.calcsize("=IIIf"): + raise RuntimeError("Block wall host/device layout mismatch") + cuda.memcpy_htod(ptr, struct.pack("=IIIf", cfg.nx, cfg.ny, + cfg.wall_slot_width, cfg.wall_normal_velocity)) p = self.patch + c, f = self._levels order = np.argsort(p.coalesce_coarse) - tables = None if layout is None else blocked_fine_tables(p, layout, order) - self._ghosts = self._upload(p.ghost_fine if layout is None else tables[0]) - self._donors = self._upload(p.ghost_coarse) - self._targets = self._upload(p.coalesce_coarse[order]) - if layout is None: - self._children = self._upload(p.coalesce_fine[order]) - self._explode = self._levels[1].module.get_function("PatchExplode") - else: - self._children = self._upload(tables[1]) - self._kind = self._upload(tables[2]) - if len(layout.halo_dst): # a 1x1 non-periodic layout has no pairs - self._halo_dst = self._upload(layout.halo_dst) - self._halo_src = self._upload(layout.halo_src) - module = self._levels[1].module - self._explode = module.get_function("BlockExplode") - self._halo_copy = module.get_function("BlockHaloCopy") + fine_tables = (None if f.layout is None else + blocked_fine_tables(p, f.layout, order)) + base_tables = (None if c.layout is None else + blocked_base_tables(p, c.layout, order, self._wall_aware)) + self._ghosts = self._upload(p.ghost_fine if fine_tables is None else fine_tables[0]) + donors = p.ghost_coarse if base_tables is None else base_tables[0] + self._donors = self._upload(donors if fine_tables is None else donors.astype(np.uint64)) + self._targets = self._upload(p.coalesce_coarse[order] if base_tables is None + else base_tables[1]) + self._children = self._upload(p.coalesce_fine[order] if fine_tables is None + else fine_tables[1]) + self._explode = f.module.get_function("PatchExplode" if fine_tables is None + else "BlockExplode") + if fine_tables is not None: + f.kind = self._upload(fine_tables[2]) + if base_tables is not None: + c.kind = self._upload(base_tables[2]) + self._coalesce = c.module.get_function("BlockCoalesce") if self._fine_tables is not None: self._fine_flags = self._upload(self._fine_tables["flags"]) self._link_tables = {k: self._upload(v) for k, v in self._fine_tables.items() @@ -257,8 +310,8 @@ class FixedPatchStepper: if not np.isfinite(array[:, mask]).all() or (array[:, mask].sum(axis=0) <= 0).any(): raise ValueError("Leaf populations must be finite with positive density") self._mask_nonleaves(*arrays) - if self.fine_layout is not None: - arrays[1] = self.fine_layout.to_blocked(arrays[1]) + arrays = [a if level.layout is None else level.layout.to_blocked(a) + for a, level in zip(arrays, self._levels)] try: for level, array in zip(self._levels, arrays): cuda.memcpy_htod(level.current, array) @@ -303,10 +356,10 @@ class FixedPatchStepper: np.uint32(count), velocity if velocity is not None else np.uintp(0), grid=((count + self._block[0] - 1) // self._block[0], 1), **kw) - def _fine_halo(self, buffer, **kw): - pairs = len(self.fine_layout.halo_dst) + def _gather_halo(self, level, buffer, **kw): + pairs = len(level.layout.halo_dst) if pairs: - self._halo_copy(buffer, self._halo_dst, self._halo_src, np.uint32(pairs), + level.halo_copy(buffer, level.halo_dst, level.halo_src, np.uint32(pairs), grid=((pairs + self._block[0] - 1) // self._block[0], 1), **kw) def step(self, n=1): @@ -322,7 +375,9 @@ class FixedPatchStepper: self._status = "advancing" for _ in range(n): self._explode(c.current, f.current, self._ghosts, self._donors, - np.uint32(len(p.ghost_fine)), np.uint64(np.prod(p.coarse_shape)), + np.uint32(len(p.ghost_fine)), + np.uint64(np.prod(p.coarse_shape if c.layout is None + else c.layout.local_shape)), grid=((len(p.ghost_fine) + self._block[0] - 1) // self._block[0], 1), **kw) if self.fine_layout is not None: @@ -330,10 +385,10 @@ class FixedPatchStepper: # Same buffer per hop: halo -> receiver-local BFL -> step. flags = self._fine_flags if self._fine_flags is not None else np.uintp(0) for src, dst in ((f.current, f.spare), (f.spare, f.current)): - self._fine_halo(src, **kw) + self._gather_halo(f, src, **kw) self._blocked_fine_boundary(src, **kw) f.kernel(src, dst, np.uint32(self.fine_layout.halo), flags, - self._kind, grid=f.grid, **kw) + f.kind, grid=f.grid, **kw) else: flags = self._fine_flags if self._fine_flags is not None else np.uintp(0) if self.fine_boundary is not None and self._fine_flags is None: @@ -342,10 +397,21 @@ class FixedPatchStepper: f.kernel(f.current, f.spare, flags, grid=f.grid, **kw) self._apply_fine_boundary(f.spare, **kw) f.kernel(f.spare, f.current, flags, grid=f.grid, **kw) - c.kernel(c.current, c.spare, f.current, self._targets, self._children, - np.uint32(len(p.coalesce_coarse)), - *(np.uint32(v) for v in p.bounds), np.uint32(self._wall_aware), - grid=c.grid, **kw) + if c.layout is None: + c.kernel(c.current, c.spare, f.current, self._targets, self._children, + np.uint32(len(p.coalesce_coarse)), + *(np.uint32(v) for v in p.bounds), np.uint32(self._wall_aware), + grid=c.grid, **kw) + else: + # Child donors must be read before staging mutates parent G. + # No halo gather may overwrite these receiver-local arrivals. + self._gather_halo(c, c.current, **kw) + self._coalesce(f.current, c.current, self._targets, self._children, + np.uint32(len(p.coalesce_coarse)), + grid=((len(p.coalesce_coarse) + self._block[0] - 1) + // self._block[0], 1), **kw) + c.kernel(c.current, c.spare, np.uint32(c.layout.halo), np.uintp(0), + c.kind, grid=c.grid, **kw) self._stream.synchronize() c.current, c.spare = c.spare, c.current self._step_count += 1 @@ -358,17 +424,16 @@ class FixedPatchStepper: """Return two synchronized G arrays, NaN on all non-leaves; no host cache.""" with self._lock, self.ctx.activate(): self._require("ready") - arrays = [np.empty(level.shape, np.float32) for level in self._levels] - if self.fine_layout is not None: - arrays[1] = np.empty(self.fine_layout.blocked_shape, np.float32) + arrays = [np.empty(level.shape if level.layout is None else level.layout.blocked_shape, + np.float32) for level in self._levels] try: for level, array in zip(self._levels, arrays): cuda.memcpy_dtoh(array, level.current) except BaseException: self._status = "failed" raise - if self.fine_layout is not None: - arrays[1] = self.fine_layout.to_dense(arrays[1]) + arrays = [a if level.layout is None else level.layout.to_dense(a) + for a, level in zip(arrays, self._levels)] self._mask_nonleaves(*arrays) return tuple(arrays) @@ -387,7 +452,7 @@ class FixedPatchStepper: except Exception as exc: error = exc self._levels.clear() - self._explode = self._halo_copy = None + self._explode = self._coalesce = None if getattr(self, "_fine_curved", None) is not None: try: self._fine_curved.free() diff --git a/tests/validation/multigrid/test_refinement_blocked_base.py b/tests/validation/multigrid/test_refinement_blocked_base.py new file mode 100644 index 0000000..45c98c9 --- /dev/null +++ b/tests/validation/multigrid/test_refinement_blocked_base.py @@ -0,0 +1,324 @@ +"""All-blocked fixed patch: transfer ordering, failure ownership and CUDA parity. + +Run CPU and GPU classes in separate processes. GPU callers must first check +resources/device availability; this file never selects or reserves a GPU. +""" +from contextlib import ExitStack, contextmanager, nullcontext +from dataclasses import replace +from pathlib import Path +import struct +import unittest +from unittest.mock import Mock, patch + +import numpy as np +import pycuda.driver as cuda + +from CelerisLab.cuda import compiler_v2 as compiler +from CelerisLab.lbm import refinement_stepper +from CelerisLab.lbm.blocks import BlockLayout2D +from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY +from CelerisLab.lbm.refinement import ACTIVE, OUTER_GHOST, FixedPatch2D +from CelerisLab.lbm.refinement_stepper import FixedPatchStepper, blocked_base_tables +from test_refinement_blocked_fine import NX, NY, BOUNDS, config, state, circle, actuation + +BASE_BLOCKS = ((24, 32), (12, 8), (3, 4), (1, 1)) +FINE_BLOCK = (16, 20) + + +def raw(backend): + arrays = [] + for level in backend._levels: + array = np.empty(level.shape if level.layout is None else level.layout.blocked_shape, + np.float32) + cuda.memcpy_dtoh(array, level.current) + arrays.append(array if level.layout is None else level.layout.to_dense(array)) + return tuple(arrays) + + +def exact(test, actual, expected, label): + """Compare the complete invalid mask and every finite FP32 bit (incl. -0).""" + np.testing.assert_array_equal(np.isnan(actual), np.isnan(expected), err_msg=label) + test.assertFalse(np.isinf(actual).any(), label) + test.assertFalse(np.isinf(expected).any(), label) + finite = np.isfinite(expected) + np.testing.assert_array_equal(actual[finite].view(np.uint32), + expected[finite].view(np.uint32), err_msg=label) + + +@contextmanager +def mocked_backend(mode="periodic", fail_alloc=None): + """Use the real constructor, uploads and launch logic, recording driver calls.""" + calls, allocations, uploads, functions = [], [], [], {} + stream = Mock() + context = Mock(activate=nullcontext, sm_arch="sm_70") + + def function(name): + if name not in functions: + functions[name] = Mock(side_effect=lambda *a, **kw: calls.append((name, a, kw))) + return functions[name] + + def allocate(size): + if len(allocations) == fail_alloc: + raise MemoryError("injected allocation failure") + if size <= 0: + raise AssertionError("empty allocation") + pointer = Mock() + allocations.append(pointer) + return pointer + + module = Mock(get_function=Mock(side_effect=function)) + module.get_global.side_effect = lambda name: (name, 16 if name == "d_block_wall" else 20) + with ExitStack() as stack: + for name in ("generate_config", "compile_kernel", "cleanup_build_dir"): + stack.enter_context(patch.object(compiler, name)) + stack.enter_context(patch.object(compiler, "load_module", return_value=module)) + stack.enter_context(patch.object(cuda, "Stream", return_value=stream)) + stack.enter_context(patch.object(cuda, "mem_alloc", side_effect=allocate)) + stack.enter_context(patch.object(cuda, "memcpy_htod", + side_effect=lambda p, a: uploads.append((p, a)))) + yield dict(context=context, calls=calls, allocations=allocations, uploads=uploads, + functions=functions, stream=stream) + + +class BlockedBaseCPU(unittest.TestCase): + def test_reject_before_cuda(self): + context = Mock(activate=Mock(side_effect=AssertionError("CUDA reached"))) + for kwargs in ({"base_block": (12, 8)}, + {"base_block": (7, 8), "fine_block": FINE_BLOCK}, + {"base_block": (12,), "fine_block": FINE_BLOCK}): + with self.subTest(kwargs=kwargs), self.assertRaises(ValueError): + FixedPatchStepper(config("periodic"), BOUNDS, context=context, **kwargs) + context.activate.assert_not_called() + + def test_all_arrival_rows_donors_and_kind_masks(self): + p = FixedPatch2D.from_config(config("periodic"), BOUNDS) + order = np.argsort(p.coalesce_coarse) + self.assertEqual(len(order), 172) + for mode in ("periodic", "slot_reservoir"): + leaves = np.ones((NY, NX), bool) + x0, y0, x1, y1 = BOUNDS + leaves[y0:y1, x0:x1] = False + if mode != "periodic": + leaves[[0, -1]] = False + for by, bx in BASE_BLOCKS: + layout = BlockLayout2D.build(NY, NX, by, bx, periodic=(True, mode == "periodic")) + donors, targets, kind = blocked_base_tables(p, layout, order, mode != "periodic") + ly, lx = layout.local_shape + L = ly * lx + expected, source_cells = [], [] + for slot in p.coalesce_coarse[order]: + q, receiver = divmod(int(slot), NX * NY) + y, x = divmod(receiver, NX) + block = (y // by) * (NX // bx) + x // bx + # Scalar coordinate construction independent of layout.owner / slot. + sy, sx = y % by + 1 - D2Q9_EY[q], x % bx + 1 - D2Q9_EX[q] + expected.append((block * 9 + q) * L + sy * lx + sx) + source_cells.append((y - D2Q9_EY[q], x - D2Q9_EX[q])) + np.testing.assert_array_equal(targets, expected) + self.assertEqual(len(set(targets.tolist())), 172) + self.assertEqual(sum(not (y0 <= y < y1 and x0 <= x < x1) + for y, x in source_cells), 8) + expected_donors = [] + for cell in p.ghost_coarse: + y, x = divmod(int(cell), NX) + block = (y // by) * (NX // bx) + x // bx + expected_donors.append(block * L + (y % by + 1) * lx + x % bx + 1) + np.testing.assert_array_equal(donors, expected_donors) + expected_kind = np.full(layout.nblocks * L, OUTER_GHOST, np.uint8) + expected_kind[layout.owner[leaves].astype(np.intp)] = ACTIVE + np.testing.assert_array_equal(kind, expected_kind) + + def test_staging_negative_controls_detect_seam_and_donor_corruption(self): + p = FixedPatch2D.from_config(config("periodic"), BOUNDS) + order = np.argsort(p.coalesce_coarse) + layout = BlockLayout2D.build(NY, NX, 3, 4) + donors, targets, _ = blocked_base_tables(p, layout, order) + dense = np.arange(9 * NY * NX, dtype=np.float32).reshape(9, NY, NX) + initial = layout.to_blocked(dense) + arrivals = np.arange(len(targets), dtype=np.float32) + 10000 + staged = initial.copy() + staged.ravel()[targets.astype(np.intp)] = arrivals + def pull(storage): + out = np.empty_like(dense) + by, bx = layout.block + ly, lx = layout.local_shape + for y in range(NY): + for x in range(NX): + b = (y // by) * (NX // bx) + x // bx + for q in range(9): + out[q, y, x] = storage[b, q, y % by + 1 - D2Q9_EY[q], + x % bx + 1 - D2Q9_EX[q]] + return out + + expected = np.array([np.roll(dense[q], (D2Q9_EY[q], D2Q9_EX[q]), axis=(0, 1)) + for q in range(9)]) + expected.ravel()[p.coalesce_coarse[order].astype(np.intp)] = arrivals + leaves = np.ones((NY, NX), bool) + leaves[6:18, 8:24] = False + # Compare ALL parent arrivals, not only the slots deliberately written. + np.testing.assert_array_equal(pull(staged)[:, leaves], expected[:, leaves]) + # A late halo gather erases receiver-local arrivals at seams. + late_gather = staged.copy() + layout.exchange(late_gather) + self.assertGreater(np.count_nonzero(pull(late_gather)[:, leaves] != expected[:, leaves]), 0) + # Canonical source owners are not the slots read by each receiver's pull. + q, receiver = np.divmod(p.coalesce_coarse[order], NX * NY) + y, x = np.divmod(receiver.astype(np.intp), NX) + owners = layout.owner[y - np.asarray(D2Q9_EY)[q.astype(np.intp)], + x - np.asarray(D2Q9_EX)[q.astype(np.intp)]] + wrong = initial.copy() + wrong.ravel()[layout.slot(owners, q).astype(np.intp)] = arrivals + self.assertGreater(np.count_nonzero(pull(wrong)[:, leaves] != expected[:, leaves]), 0) + # With one parent block, all eight grazing sources are owner cores. + # Each is read for four fine ghost children: late explosion corrupts32. + single = BlockLayout2D.build(NY, NX, NY, NX) + donors, targets, _ = blocked_base_tables(p, single, order) + initial = single.to_blocked(dense) + staged = initial.copy() + staged.ravel()[targets.astype(np.intp)] = arrivals + donor_slots = single.slot(donors[None, :], np.arange(9)[:, None]).astype(np.intp) + self.assertEqual(np.count_nonzero(staged.ravel()[donor_slots] + != initial.ravel()[donor_slots]), 32) + + def test_raw_mask_gate_detects_leak_hidden_by_download_mask(self): + valid = state(NY, NX, 1) + valid[:, 6:18, 8:24] = np.nan + leaked = valid.copy() + leaked[5, 6, 8] = 42 # skipped invalid-cell store leaves stale staging data + public = leaked.copy() + public[:, 6:18, 8:24] = np.nan + exact(self, public, valid, "masked download cannot detect stale staging") + with self.assertRaises(AssertionError): + exact(self, leaked, valid, "raw owner-core must detect stale staging") + + def test_launch_order_stream_wall_constants_and_buffer_swaps(self): + with mocked_backend("slot_reservoir") as m: + backend = FixedPatchStepper(config("slot_reservoir"), BOUNDS, context=m["context"], + coarse_boundary="slot_reservoir", fine_block=FINE_BLOCK, + base_block=(3, 4), fine_objects=[circle(actuation(.002))]) + try: + backend.initialize(state(NY, NX, 1), state(*backend.patch.fine_shape, 2)) + coarse, fine = backend._levels + current, spare, fine_current = coarse.current, coarse.spare, fine.current + backend.step(2) + names = [name for name, _, _ in m["calls"]] + interval = ["BlockExplode", "BlockHaloCopy", "BlockCurvedBoundary", "BlockStep", + "BlockHaloCopy", "BlockCurvedBoundary", "BlockStep", + "BlockHaloCopy", "BlockCoalesce", "BlockStep"] + self.assertEqual(names, interval * 2) + self.assertEqual(backend.step_count, 2) + self.assertIs(coarse.current, current) + self.assertIs(coarse.spare, spare) + self.assertIs(fine.current, fine_current) + for _, _, kw in m["calls"]: + self.assertIs(kw["stream"], m["stream"]) + for offset, src, dst in ((0, current, spare), (10, spare, current)): + explode, gather, stage, collide = [m["calls"][offset + i][1] for i in (0, 7, 8, 9)] + self.assertIs(explode[0], src) + self.assertIs(gather[0], src) + self.assertIs(stage[0], fine_current) + self.assertIs(stage[1], src) + self.assertIs(collide[0], src) + self.assertIs(collide[1], dst) + self.assertIs(collide[4], coarse.kind) + walls = [a for ptr, a in m["uploads"] if isinstance(ptr, str) and ptr == "d_block_wall"] + self.assertEqual(len(walls), 1) + width, height, slot, speed = struct.unpack("=IIIf", walls[0]) + self.assertEqual((width, height, slot), (NX, NY, 8)) + self.assertAlmostEqual(speed, .001) + finally: + backend.close() + for allocation in m["allocations"]: + allocation.free.assert_called_once_with() + + def test_allocation_failure_releases_every_owned_buffer(self): + # Fail during the first level, second level and transfer-table uploads. + for fail_at in (1, 6, 10): + with self.subTest(fail_at=fail_at), mocked_backend(fail_alloc=fail_at) as m: + with self.assertRaisesRegex(MemoryError, "injected allocation"): + FixedPatchStepper(config("periodic"), BOUNDS, context=m["context"], + fine_block=FINE_BLOCK, base_block=(3, 4)) + self.assertEqual(len(m["allocations"]), fail_at) + for allocation in m["allocations"]: + allocation.free.assert_called_once_with() + + def test_failed_parent_stage_cannot_publish_or_resume(self): + with mocked_backend() as m: + backend = FixedPatchStepper(config("periodic"), BOUNDS, context=m["context"], + fine_block=FINE_BLOCK, base_block=(3, 4)) + backend.initialize(state(NY, NX, 1), state(*backend.patch.fine_shape, 2)) + m["functions"]["BlockCoalesce"].side_effect = RuntimeError("injected parent stage failure") + current = backend._levels[0].current + with self.assertRaisesRegex(RuntimeError, "injected parent stage"): + backend.step() + self.assertEqual(backend.step_count, 0) + self.assertEqual(backend._status, "failed") + self.assertIs(backend._levels[0].current, current) + with self.assertRaisesRegex(RuntimeError, "failed"): + backend.download() + with self.assertRaisesRegex(RuntimeError, "failed"): + backend.step() + backend.close() + for allocation in m["allocations"]: + allocation.free.assert_called_once_with() + + +class BlockedBaseGPU(unittest.TestCase): + def test_bodyfree_exact_parity_and_raw_invalid_slots(self): + self.compare([(None, ((24, 32), (3, 4), (1, 1)), 8)]) + + def test_actuated_circle_exact_parity_and_raw_invalid_slots(self): + self.compare([(.002, ((12, 8),), 8), (-.002, ((12, 8),), 8)]) + + def test_non_power_two_slot_width_exact_parity(self): + self.compare([(None, ((3, 4),), 7), (None, ((3, 4),), 10)], ("slot_reservoir",)) + + def compare(self, cases, modes=("periodic", "slot_reservoir")): + from CelerisLab.cuda.context import CudaContext + expected_source = Path(__file__).resolve().parents[3] / "src/CelerisLab/lbm/refinement_stepper.py" + self.assertEqual(Path(refinement_stepper.__file__).resolve(), expected_source) + with CudaContext(0) as context: + for mode in modes: + for speed, layouts, width in cases: + cfg = config(mode) + if mode == "slot_reservoir": + cfg = replace(cfg, wall_slot_width=width) + p = FixedPatch2D.from_config(cfg, BOUNDS) + arrays = state(NY, NX, 1), state(*p.fine_shape, 2) + reference = None + for base_block in (None, *layouts): + label = f"{mode} width={width} speed={speed} base={base_block}" + snapshots = [] + with FixedPatchStepper(cfg, BOUNDS, context=context, coarse_boundary=mode, + fine_block=FINE_BLOCK, base_block=base_block, + fine_objects=() if speed is None else [circle(actuation(speed))]) as backend: + backend.initialize(*arrays) + for step in (0, 1, 2, 8): + if step == 1 and base_block is not None: + # Blocked staging dirties invalid input slots; its output + # must poison them. Dense parent does not stage there. + coarse = backend._levels[0] + spare = coarse.layout.to_blocked(np.full(coarse.shape, 42, np.float32)) + cuda.memcpy_htod(coarse.spare, spare) + backend.step(step - backend.step_count) + public, owner = backend.download(), raw(backend) + self.assertEqual(backend.step_count, step) + self.assertTrue(np.isnan(owner[0][:, 6:18, 8:24]).all(), label) + if mode != "periodic": + self.assertTrue(np.isnan(owner[0][:, (0, -1), :]).all(), label) + snapshots.append(public + owner) + backend.initialize(*arrays) + self.assertEqual(backend.step_count, 0) + for actual, expected in zip(backend.download() + raw(backend), snapshots[0]): + exact(self, actual, expected, label + " reset") + if reference is None: + reference = snapshots + else: + for step, actual, expected in zip((0, 1, 2, 8), snapshots, reference): + for index, (a, b) in enumerate(zip(actual, expected)): + exact(self, a, b, f"{label} step={step} array={index}") + print(f"PASS {label}: full finite bits/NaN masks at 0,1,2,8 + reset", flush=True) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/validation/multigrid/test_refinement_blocked_fine.py b/tests/validation/multigrid/test_refinement_blocked_fine.py index 06448ce..23e36e6 100644 --- a/tests/validation/multigrid/test_refinement_blocked_fine.py +++ b/tests/validation/multigrid/test_refinement_blocked_fine.py @@ -128,8 +128,8 @@ class BlockedFineCPU(unittest.TestCase): backend = FixedPatchStepper(config("periodic"), BOUNDS, context=context, fine_block=FINE_BLOCKS[0]) self.assertEqual(len(backend.fine_layout.halo_dst), 0) - self.assertIsNone(backend._halo_dst) - backend._fine_halo(Mock()) # no launch without pairs + self.assertIsNone(backend._levels[1].halo_dst) + backend._gather_halo(backend._levels[1], Mock()) # no launch without pairs module.get_function.return_value.assert_not_called() backend.close() self.assertTrue(sizes) @@ -187,7 +187,7 @@ class BlockedFineCPU(unittest.TestCase): backend, calls = self.mocked_backend([circle(spec)]) c, f = backend._levels coarse, cur, spare = c.current, f.current, f.spare - flags, kind = backend._fine_flags, backend._kind + flags, kind = backend._fine_flags, backend._levels[1].kind velocity = backend._link_tables.get("surface_velocity") self.assertEqual(velocity is None, spec is None) backend.step(1) diff --git a/tests/validation/multigrid/test_refinement_followup.py b/tests/validation/multigrid/test_refinement_followup.py index 9d65166..9af399f 100644 --- a/tests/validation/multigrid/test_refinement_followup.py +++ b/tests/validation/multigrid/test_refinement_followup.py @@ -111,7 +111,7 @@ class FaultGuardsCPU(unittest.TestCase): backend._lock = RLock() backend._status = "ready" backend._stream = Mock() - backend._levels = [Mock(shape=(9, 1, 1)), Mock(shape=(9, 1, 1))] + backend._levels = [Mock(shape=(9, 1, 1), layout=None), Mock(shape=(9, 1, 1), layout=None)] backend._allocations = [Mock(), Mock()] return backend