--- a/src/CelerisLab/lbm/refinement_stepper.py +++ b/src/CelerisLab/lbm/refinement_stepper.py @@ -1,5 +1,5 @@ # CelerisLab/lbm/refinement_stepper.py -"""Opt-in fixed periodic patch backend; no Simulation/channel/objects adapter. +"""Opt-in fixed patch backend; no Simulation/objects adapter. The caller owns the CudaContext and must close this backend before that context. Only synchronized postcollision_g leaf populations are public. No checkpoint, @@ -27,8 +27,12 @@ grid: tuple -class PeriodicPatchStepper: - """D2Q9/SRT/FP32, one unforced 2:1 patch in a doubly periodic coarse box. +class FixedPatchStepper: + """D2Q9/SRT/FP32, one unforced internal 2:1 patch in an x-periodic box. + + ``coarse_boundary`` is ``periodic`` (default) or ``slot_reservoir``: + a bottom on-site slot and top fixed-density reservoir from ``cfg``. + Physical walls remain coarse; fine halo donors must be ordinary fluid. ``initialize(coarse_g, fine_g)`` accepts physical Q-major (9,ny,nx) arrays already in postcollision phase. It does not collide or build channel flags. @@ -36,10 +40,20 @@ On a device/launch failure the instance becomes unusable: close and recreate it, rather than pretending the partially advanced fine buffer rolled back. """ - def __init__(self, cfg, bounds, *, context): + def __init__(self, cfg, bounds, *, context, coarse_boundary="periodic"): if HALO != 4: raise RuntimeError("refinement.cu requires a four-cell fine halo") + if coarse_boundary not in ("periodic", "slot_reservoir"): + raise ValueError("coarse_boundary must be periodic or slot_reservoir") self.patch = FixedPatch2D.from_config(cfg, bounds) + self._wall_aware = coarse_boundary == "slot_reservoir" + if self._wall_aware: + if (cfg.y_wall_bc != "on_site_velocity" or cfg.top_wall_bc != "reservoir" + or cfg.velocity != 0 or not 2 <= cfg.wall_slot_width < cfg.nx): + raise ValueError("slot_reservoir requires on_site_velocity, reservoir, " + "velocity=0 and 2 <= wall_slot_width < nx") + 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") self.ctx = context self._lock = RLock() self._status = "new" @@ -50,6 +64,9 @@ self._block = (cfg.threads_per_block, 1, 1) fine_cfg = replace(cfg, nx=self.patch.fine_shape[1], ny=self.patch.fine_shape[0], viscosity=2 * cfg.viscosity) + if self._wall_aware: + fine_cfg = replace(fine_cfg, y_wall_bc="bounce_back", top_wall_bc="velocity", + wall_slot_width=0, wall_normal_velocity=0.0, velocity=0.0) with self.ctx.activate(): try: self._stream = cuda.Stream() @@ -102,7 +119,7 @@ def _require(self, *states): if self._status not in states: - raise RuntimeError(f"Periodic patch is {self._status}; expected {' or '.join(states)}") + raise RuntimeError(f"Fixed patch is {self._status}; expected {' or '.join(states)}") @property def step_count(self): @@ -111,10 +128,12 @@ def _mask_nonleaves(self, coarse, fine): x0, y0, x1, y1 = self.patch.bounds coarse[:, y0:y1, x0:x1] = np.nan + if self._wall_aware: + coarse[:, (0, -1), :] = np.nan fine[:, self.patch.fine_kind != ACTIVE] = np.nan def initialize(self, coarse_g, fine_g): - """Upload an explicit periodic G_0; invalid non-leaves are discarded.""" + """Upload an explicit G_0; invalid non-leaves are discarded.""" with self._lock, self.ctx.activate(): self._require("new", "ready") arrays = [np.array(a, dtype=np.float32, order="C", copy=True) @@ -124,6 +143,8 @@ x0, y0, x1, y1 = self.patch.bounds leaf = np.ones(self.patch.coarse_shape, dtype=bool) leaf[y0:y1, x0:x1] = False + if self._wall_aware: + leaf[[0, -1]] = False for array, mask in zip(arrays, (leaf, self.patch.fine_kind == ACTIVE)): 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") @@ -158,7 +179,8 @@ f.kernel(f.spare, f.current, 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), grid=c.grid, **kw) + *(np.uint32(v) for v in p.bounds), np.uint32(self._wall_aware), + grid=c.grid, **kw) self._stream.synchronize() c.current, c.spare = c.spare, c.current self._step_count += 1 @@ -216,3 +238,7 @@ except Exception: if exc[0] is None: raise + + +# Preserve the original entry point and its periodic default. +PeriodicPatchStepper = FixedPatchStepper --- a/src/CelerisLab/lbm/kernels/refinement.cu +++ b/src/CelerisLab/lbm/kernels/refinement.cu @@ -1,4 +1,4 @@ -// Fixed, unforced, doubly periodic 2:1 patch. Separate compilation entry; +// Fixed, unforced internal 2:1 patch. Separate compilation entry; // default OneStep/channel kernels are deliberately not included. #include #include @@ -10,6 +10,7 @@ #include "operators/macro.cuh" #include "operators/collision_srt.cuh" #include "streaming/pull_double_buffer.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 "Refinement requires unshifted FP32 D2Q9 SRT double buffers without LES" @@ -69,13 +70,15 @@ extern "C" __global__ void PatchCoarseStep( const float* in, float* out, const float* fine, const uint64_t* targets, const uint64_t* children, unsigned int count, - unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1) + unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1, + unsigned int wall_aware) { unsigned int x, y; unsigned long k; index_from_thread(x, y, k); if (x >= NX || y >= NY) return; // Both buffers were poisoned at initialization; covered cells stay invalid. if (x >= x0 && x < x1 && y >= y0 && y < y1) return; + if (wall_aware && (y == 0 || y == NY - 1)) return; unsigned long j[NQ]; compute_neighbors(x, y, j); float f[NQ]; @@ -100,6 +103,10 @@ } // Includes the eight grazing-corner links with UNCOVERED coarse donors. // All replacement arrivals must be in place before computing rho/u. + if (wall_aware) { + if (y == 1) reconstruct_bottom_velocity(f, 0.0f, bottom_wall_normal_velocity(x)); + else if (y == NY - 2) reconstruct_top_reservoir(f, RHO); + } patch_collide(f, k); stream_pull_store(k, f, out); } --- a/tests/validation/multigrid/test_refinement_followup.py +++ b/tests/validation/multigrid/test_refinement_followup.py @@ -1,10 +1,11 @@ -"""Only review gaps: CPU fault guards, normal-build seed1/8 and rectangle/1. +"""Bounded fixed-patch checks: fault guards, periodic geometry and coarse slot BC. Use the existing bounded runner/environment; no audit-build rerun. The scalar -rectangle oracle covers one interval, not a second general-purpose reference. +oracle covers the fixed tau pair and two boundary modes used by these cases. """ from contextlib import nullcontext from dataclasses import replace +import json from pathlib import Path import resource from threading import RLock @@ -21,8 +22,8 @@ PeriodicPatchStepper) -def scalar_interval(coarse, fine, bounds): - """Independent G→G single interval, geometric links, no production tables.""" +def scalar_interval(coarse, fine, bounds, slot=None): + """Independent G→G interval; optional (amplitude, width) bottom slot/top rho=1.""" _, ny, nx = coarse.shape _, fy, fx = fine.shape x0, y0, x1, y1 = bounds @@ -53,7 +54,7 @@ result = np.full(coarse.shape, np.nan) for y in range(ny): for x in range(nx): - if x0 <= x < x1 and y0 <= y < y1: + if (x0 <= x < x1 and y0 <= y < y1) or (slot is not None and y in (0, ny-1)): continue values = [] for q in range(9): @@ -64,11 +65,45 @@ else: value = coarse[q, (y-cy[q]) % ny, (x-cx[q]) % nx] values.append(value) + if slot is not None and y == 1: + amplitude, width = slot + s = 2*((x+.5)-nx/2)/width + v = amplitude * max(0.0, 1-s*s)**2 + rho = (values[0]+values[1]+values[2] + + 2*(values[4]+values[6]+values[7]))/(1-v) + values[3] = values[4] + 2*rho*v/3 + values[5] = values[6] + (values[2]-values[1])/2 + rho*v/6 + values[8] = values[7] + (values[1]-values[2])/2 + rho*v/6 + elif slot is not None and y == ny-2: + jy = values[0]+values[1]+values[2]+2*(values[3]+values[5]+values[8])-1 + values[4] = values[3] - 2*jy/3 + values[6] = values[5] + (values[1]-values[2])/2 - jy/6 + values[7] = values[8] + (values[2]-values[1])/2 - jy/6 result[:, y, x] = collide(values, .8) return result, current class FaultGuardsCPU(unittest.TestCase): + def test_boundary_mode_and_margin_contract(self): + cfg = LBMConfig(nx=32, ny=24, viscosity=.1, velocity=0, wall_normal_velocity=.001, + wall_slot_width=8, y_wall_bc="on_site_velocity", top_wall_bc="reservoir") + context = Mock() + context.activate.side_effect = RuntimeError("validated before CUDA") + for mode, bounds in (("unknown", (12, 4, 20, 10)), + ("slot_reservoir", (12, 3, 20, 10)), + ("slot_reservoir", (12, 4, 20, 21))): + with self.subTest(mode=mode, bounds=bounds): + with self.assertRaises(ValueError): + implementation.FixedPatchStepper(cfg, bounds, context=context, coarse_boundary=mode) + context.activate.assert_not_called() + # Both exact channel limits and the original periodic two-cell margin + # must reach context activation, with no allocation or device required. + for mode, bounds in (("slot_reservoir", (12, 4, 20, 20)), + ("periodic", (2, 2, 6, 5))): + with self.subTest(mode=mode, bounds=bounds): + with self.assertRaisesRegex(RuntimeError, "validated before CUDA"): + implementation.FixedPatchStepper(cfg, bounds, context=context, coarse_boundary=mode) + def fake_backend(self): backend = PeriodicPatchStepper.__new__(PeriodicPatchStepper) backend.ctx = Mock(activate=nullcontext) @@ -145,6 +180,126 @@ def tearDown(self): rss = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / 1024 print(f"host_peak_rss={rss:.3f} MiB", flush=True) + + def test_slot_reservoir_eight_steps(self): + nx, ny, bounds, amplitude, width = 32, 24, (12, 4, 20, 10), .001, 8 + x0, y0, x1, y1 = bounds + fx, fy = 2*(x1-x0)+8, 2*(y1-y0)+8 + cfg = LBMConfig(nx=nx, ny=ny, viscosity=.1, velocity=0, rho=1, + wall_normal_velocity=amplitude, wall_slot_width=width, + y_wall_bc="on_site_velocity", top_wall_bc="reservoir") + cx = np.array((0, 1, -1, 0, 0, 1, -1, 1, -1)) + cy = np.array((0, 0, 0, 1, -1, 1, -1, -1, 1)) + weights = np.array((4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36)) + + def slot_velocity(x): + return amplitude * np.maximum(0, 1-(2*(x-nx/2)/width)**2)**2 + + def equilibrium(shape, origin, dx): + ys = origin[1] + dx*(np.arange(shape[0])+.5) + xs = origin[0] + dx*(np.arange(shape[1])+.5) + v = slot_velocity(xs)[None, :] * (1-(ys[:, None]-1.5)/(ny-3)) + cu = cy[:, None, None] * v + return np.array(weights[:, None, None]*(1+3*cu+4.5*cu*cu-1.5*v*v), np.float32) + + coarse_leaf = np.ones((ny, nx), bool) + coarse_leaf[[0, -1]] = False + coarse_leaf[y0:y1, x0:x1] = False + fine_leaf = np.zeros((fy, fx), bool) + fine_leaf[4:-4, 4:-4] = True + initial = [equilibrium((ny, nx), (0, 0), 1), + equilibrium((fy, fx), (x0-2, y0-2), .5)] + initial[0][:, ~coarse_leaf] = np.nan + initial[1][:, ~fine_leaf] = np.nan + # A postcollision zero-mass/zero-momentum perturbation at a grazing + # source distinguishes closed-corner coupling from ordinary coarse pull. + initial[0][[5, 6], y0-1, x0] += 1e-5 + initial[0][[7, 8], y0-1, x0] -= 1e-5 + expected = tuple(a.astype(np.float64) for a in initial) + + def discrete_mass(arrays): + return (arrays[0][:, coarse_leaf].sum(dtype=np.float64) + + .25*arrays[1][:, fine_leaf].sum(dtype=np.float64)) + + # Observe the actual generated headers, without changing generation or + # compilation. In particular the fine module must have slot BC disabled. + generated = [] + generate = compiler.generate_config + + def record_config(*args, **kwargs): + build = generate(*args, **kwargs) + metadata = json.loads((Path(build)/".celeris-build").read_text()) + definitions = {} + for line in (Path(build)/"config/config_method.h").read_text().splitlines(): + parts = line.split() + if len(parts) == 3 and parts[0] == "#define": + definitions[parts[1]] = parts[2] + generated.append((metadata["config"], definitions)) + return build + + max_error = max_bc_error = 0.0 + with CudaContext(0) as context: + with patch.object(compiler, "generate_config", side_effect=record_config): + backend = implementation.FixedPatchStepper(cfg, bounds, context=context, + coarse_boundary="slot_reservoir") + with backend: + self.assertEqual(len(generated), 2) + fine_config, fine_defines = generated[1] + self.assertEqual((fine_config["nx"], fine_config["ny"]), (fx, fy)) + for name in ("Y_WALL_BC", "TOP_WALL_BC", "WALL_SLOT_WIDTH", + "WALL_NORMAL_VELOCITY", "TOP_WALL_NORMAL_VELOCITY"): + self.assertEqual(float(fine_defines[name].removesuffix("f")), 0, name) + self.assertEqual(int(generated[0][1]["Y_WALL_BC"]), 2) + self.assertEqual(int(generated[0][1]["TOP_WALL_BC"]), 1) + # Independently derive all target links, including the eight + # whose coarse sources are uncovered, from physical midpoints. + targets, corners = set(), set() + for y, x in np.argwhere(coarse_leaf): + for q in range(1, 9): + if x0 <= x+.5-cx[q]/2 <= x1 and y0 <= y+.5-cy[q]/2 <= y1: + target = q*nx*ny+y*nx+x + targets.add(target) + sx, sy = x-cx[q], y-cy[q] + if not (x0 <= sx < x1 and y0 <= sy < y1): + corners.add(target) + self.assertEqual(len(corners), 8) + self.assertEqual(set(backend.patch.coalesce_coarse.tolist()), targets) + backend.initialize(*initial) + before = discrete_mass(initial) + for step in range(1, 9): + expected = scalar_interval(*expected, bounds, slot=(amplitude, width)) + backend.step() + self.assertEqual(backend.step_count, step) + raw = [np.empty(a.shape, np.float32) for a in initial] + for level, array in zip(backend._levels, raw): + cuda.memcpy_dtoh(array, level.current) + # Read raw storage before masking: active NaNs cannot be + # hidden by the public download's non-leaf cleanup. + self.assertTrue(np.isfinite(raw[0][:, coarse_leaf]).all()) + self.assertTrue(np.isfinite(raw[1][:, fine_leaf]).all()) + self.assertTrue(np.isnan(raw[0][:, ~coarse_leaf]).all()) + for got, want in zip(raw, expected): + np.testing.assert_array_equal(np.isnan(got), np.isnan(want)) + max_error = max(max_error, float(np.nanmax(np.abs(got-want)))) + self.assertLessEqual(max_error, POP_ATOL) + public = backend.download() + np.testing.assert_array_equal(public[0], raw[0]) + np.testing.assert_array_equal(public[1][:, fine_leaf], raw[1][:, fine_leaf]) + self.assertTrue(np.isnan(public[1][:, ~fine_leaf]).all()) + walls = raw[0][:, [1, ny-2], :].astype(np.float64) + rho = walls.sum(axis=0) + ux = np.einsum("q,qyx->yx", cx, walls)/rho + vy = np.einsum("q,qyx->yx", cy, walls)/rho + max_bc_error = max(max_bc_error, float(np.abs(ux).max()), + float(np.abs(vy[0]-slot_velocity(np.arange(nx)+.5)).max()), + float(np.abs(rho[1]-1).max())) + self.assertLessEqual(max_bc_error, 2e-6) + change = discrete_mass(raw)-before + bottom_flux, top_flux = (rho*vy).sum(axis=1) + print(f"slot/reservoir steps=8 max_population={max_error:.9g} " + f"max_BC_error={max_bc_error:.9g} discrete_leaf_mass_change={change:+.9g} " + f"final_sum_rho_v_bottom={bottom_flux:.9g} top={top_flux:.9g}; " + "raw halos/corners/fine BC isolation PASS", flush=True) def test_nonequilibrium_seed1_eight_steps(self): spec = ref.make_spec(collision="srt")