diff --git a/src/CelerisLab/lbm/block_stepper.py b/src/CelerisLab/lbm/block_stepper.py index 1bbed04..c4c0cfe 100644 --- a/src/CelerisLab/lbm/block_stepper.py +++ b/src/CelerisLab/lbm/block_stepper.py @@ -1,12 +1,14 @@ # CelerisLab/lbm/block_stepper.py -"""Opt-in body-free periodic backend for one BlockLayout2D level. +"""Opt-in periodic backend for one BlockLayout2D level, optional static circles. D2Q9/SRT/FP32 unforced double buffers. Each step refreshes every halo from its owner core, then pulls/collides core cells into the spare buffer, then swaps. Only halo pair tables are uploaded; dense ownership stays on the host. The caller owns the CudaContext and must close this backend before it. -No bodies, walls, refinement, checkpoint or runtime parameter mutation; -wall/slot/inlet fields in cfg are ignored (fully periodic level). +Optional ``objects`` are stationary passive circles (static linear BFL via +receiver-local staging: halo refresh -> cut links -> pull, no copy between). +No force/actuation readout, motion, walls, refinement, checkpoint or runtime +parameter mutation; wall/slot/inlet fields in cfg are ignored. """ from dataclasses import replace from numbers import Integral @@ -18,17 +20,64 @@ import pycuda.driver as cuda from ..cuda import compiler_v2 as compiler from .blocks import NQ, BlockLayout2D +from .descriptors import D2Q9_EX, D2Q9_EY, FLUID, SOLID MAX_GRID_DIM = 65535 # CUDA gridDim.y/z limit (rows, blocks) +def _static_circles(layout, objects): + """Host tables for stationary circles; rejects before any CUDA work. + + Every periodic neighbor of every fluid core is inspected, including domain + edges that CircleGeometry does not enumerate, so clipped or incomplete + geometry cannot pass. Each physical link has one receiver-owning block. + """ + from ..body.geometry.circle import CircleGeometry + + ny, nx = layout.shape + flags = np.full(ny * nx, FLUID, np.uint16) + for obj in objects: + if (type(obj.geometry) is not CircleGeometry or obj.is_sensor or obj.is_force_region + or obj.surface_actuation is not None + or any(getattr(obj.state, name) != 0 for name in ("vx", "vy", "omega")) + or any(vars(obj.control).values())): + raise ValueError("Block bodies must be stationary passive circles") + mask = obj.geometry.build_flag_mask(nx, ny) + solid = (mask & SOLID) != 0 + if not solid.any() or np.any(flags[solid] & SOLID): + raise ValueError("Block circles must be resolved and disjoint") + flags[solid] = mask[solid] + solid = ((flags & SOLID) != 0).reshape(ny, nx) + if solid[[0, -1]].any() or solid[:, [0, -1]].any(): + raise ValueError("Block circles must not touch the domain edge") + links = [link for obj in objects + for link in obj.geometry.build_cut_links(nx, ny, domain_flags=flags)] + def source(r, q): + y, x = divmod(r, nx) + return (y + D2Q9_EY[q]) % ny, (x + D2Q9_EX[q]) % nx + receivers = {(int(link.fluid_idx), int(link.dir)) for link in links} + required = {(int(r), q) for r in np.flatnonzero(~solid.ravel()) + for q in range(1, NQ) if solid[source(r, q)]} + if len(receivers) != len(links) or receivers != required: + raise ValueError("Cut links must cover every fluid pull from a solid exactly once") + mapped = layout.map_cut_links(links) # unique receiver owner and local writer + cells = int(np.prod(layout.local_shape)) + blocked_flags = np.zeros(layout.nblocks * cells, np.uint16) # halo entries unused + blocked_flags[layout.owner.ravel().astype(np.intp)] = flags + return dict(solid=solid, flags=blocked_flags, block=mapped.block, + cell=(mapped.receiver % cells).astype(np.uint32), + dir=np.array([link.dir for link in links], np.uint8), + q=np.array([link.q for link in links], np.float32), + fallback=np.array([link.fallback for link in links], np.uint8)) + + class PeriodicBlockStepper: """G->G stepping of Q-major (9, ny, nx) postcollision populations. On a device/launch failure the instance becomes unusable: close and recreate it rather than assuming a partially advanced buffer rolled back. """ - def __init__(self, cfg, layout, *, context): + def __init__(self, cfg, layout, *, context, objects=()): if not isinstance(layout, BlockLayout2D) or layout.periodic != (True, True): raise ValueError("PeriodicBlockStepper requires a doubly periodic BlockLayout2D") checked = replace(cfg) @@ -44,8 +93,14 @@ class PeriodicBlockStepper: tpb = checked.threads_per_block if layout.nblocks > MAX_GRID_DIM or ly > MAX_GRID_DIM or not 0 < pairs < 2**31: raise ValueError("block count, block rows or halo pairs exceed launch limits") + self.boundary = _static_circles(layout, tuple(objects)) if objects else None + links = 0 if self.boundary is None else len(self.boundary["cell"]) + if links >= 2**31: + raise ValueError("cut links exceed launch limits") self.layout = layout self.ctx = context + self._links = links + self._link_grid = ((links + tpb - 1) // tpb, 1) # Wall/slot fields are not applied; neutralize them like the fine patch level. self._local_cfg = replace(checked, nx=lx, ny=ly, y_wall_bc="bounce_back", top_wall_bc="velocity", wall_slot_width=0, wall_normal_velocity=0.0, velocity=0.0) @@ -57,11 +112,13 @@ class PeriodicBlockStepper: self._status = "new" self._step_count = 0 self._allocations = [] - self._stream = self._module = self._copy = self._step = None + self._stream = self._module = self._copy = self._step = self._curved = None + self._tables = {} with self.ctx.activate(): try: free, _ = cuda.mem_get_info() - if 2 * self._bytes + layout.halo_dst.nbytes + layout.halo_src.nbytes > free: + tables = sum(a.nbytes for k, a in (self.boundary or {}).items() if k != "solid") + if 2 * self._bytes + layout.halo_dst.nbytes + layout.halo_src.nbytes + tables > free: raise MemoryError("blocked buffers exceed free device memory") self._stream = cuda.Stream() config = self._local_cfg @@ -85,6 +142,9 @@ class PeriodicBlockStepper: self._src = self._upload(layout.halo_src) self._copy = module.get_function("BlockHaloCopy") self._step = module.get_function("BlockStep") + if self.boundary is not None: + self._tables = {k: self._upload(a) for k, a in self.boundary.items() if k != "solid"} + self._curved = module.get_function("BlockCurvedBoundary") except BaseException: try: self.close() @@ -118,8 +178,10 @@ class PeriodicBlockStepper: dense = np.array(g, dtype=np.float32, order="C", copy=True) if dense.shape != (NQ, *self.layout.shape): raise ValueError("Expected a Q-major (9, ny, nx) array") - if not np.isfinite(dense).all() or (dense.sum(axis=0) <= 0).any(): - raise ValueError("Populations must be finite with positive density") + fluid = np.ones(self.layout.shape, bool) if self.boundary is None else ~self.boundary["solid"] + if not np.isfinite(dense[:, fluid]).all() or (dense[:, fluid].sum(axis=0) <= 0).any(): + raise ValueError("Fluid populations must be finite with positive density") + dense[:, ~fluid] = np.nan # Solid cores are deliberately invalid. blocked = self.layout.to_blocked(dense) spare = np.full_like(blocked, np.nan) try: @@ -143,8 +205,12 @@ class PeriodicBlockStepper: for _ in range(n): self._copy(self._current, self._dst, self._src, np.uint32(len(self.layout.halo_dst)), grid=self._copy_grid, **kw) + if self._curved is not None: + t = self._tables + self._curved(self._current, t["block"], t["cell"], t["dir"], t["q"], + t["fallback"], np.uint32(self._links), grid=self._link_grid, **kw) self._step(self._current, self._spare, np.uint32(self.layout.halo), - grid=self._grid, **kw) + self._tables.get("flags", np.uintp(0)), grid=self._grid, **kw) self._stream.synchronize() self._current, self._spare = self._spare, self._current self._step_count += 1 @@ -179,7 +245,8 @@ class PeriodicBlockStepper: self._stream.synchronize() except Exception as exc: error = exc - self._copy = self._step = self._module = None + self._copy = self._step = self._curved = self._module = None + self._tables = {} for pointer in self._allocations: try: pointer.free() diff --git a/src/CelerisLab/lbm/kernels/blocks.cu b/src/CelerisLab/lbm/kernels/blocks.cu index b1d0813..da3f6ad 100644 --- a/src/CelerisLab/lbm/kernels/blocks.cu +++ b/src/CelerisLab/lbm/kernels/blocks.cu @@ -4,14 +4,18 @@ // and linear_index address one block; blockIdx.z selects the block, whose // populations start at b*NQ*TOTAL_CELLS (block-contiguous, Q-major inside). #include +#include #include "config.h" #include "core/precision.cuh" #include "core/descriptors.cuh" #include "core/layout.cuh" #include "core/params.cuh" +#include "core/flags.cuh" +#include "core/obs.cuh" #include "operators/macro.cuh" #include "operators/collision_srt.cuh" #include "streaming/pull_double_buffer.cuh" +#include "boundary/curved_boundary.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" @@ -30,20 +34,38 @@ extern "C" __global__ void BlockHaloCopy( to[index_f(d % TOTAL_CELLS, q)] = from[index_f(s % TOTAL_CELLS, q)]; } +// Static cut links, receiver-local staging: cell is the receiver's local core +// cell, so the existing BFL writes the solid-source window image that this +// block's next pull reads (no halo copy between this kernel and BlockStep). +extern "C" __global__ void BlockCurvedBoundary( + fpxx* f, const uint32_t* block, const uint32_t* cell, const unsigned char* dir, + const float* q, const unsigned char* fallback, unsigned int count) +{ + unsigned int row = blockIdx.x * blockDim.x + threadIdx.x; + if (row >= count) return; + apply_bouzidi_link(dir[row], q[row], cell[row], f + (unsigned long)block[row] * NQ * TOTAL_CELLS, + 0.0f, 0.0f, 0.0f, 0.0f, fallback[row], nullptr, -1); +} + // Pull/collide/store core cells only; halo cells of out are left untouched. -extern "C" __global__ void BlockStep(const fpxx* in, fpxx* out, unsigned int halo) { +// flags (block cell ids, nullable): solid cores store NaN and never collide. +extern "C" __global__ void BlockStep(const fpxx* in, fpxx* out, unsigned int halo, + const uint16_t* flags) { unsigned int x, y; unsigned long k; index_from_thread(x, y, k); if (x < halo || y < halo || x >= NX - halo || y >= NY - halo) return; const unsigned long base = (unsigned long)blockIdx.z * NQ * TOTAL_CELLS; - in += base; - out += base; - unsigned long j[NQ]; - compute_neighbors(x, y, j); // Core neighbors stay inside the window: no wrap. - float f[NQ], rho, ux, uy, feq[NQ]; - stream_pull_load(k, f, in, j); - compute_rho_u(f, rho, ux, uy); - compute_feq(rho, ux, uy, feq); - collide_srt_no_force(f, feq, d_params.omega); - stream_pull_store(k, f, out); + float f[NQ]; + if (flags != nullptr && is_solid(flags[(unsigned long)blockIdx.z * TOTAL_CELLS + k])) { + for (int q = 0; q < NQ; ++q) f[q] = NAN; + } else { + unsigned long j[NQ]; + compute_neighbors(x, y, j); // Core neighbors stay inside the window: no wrap. + float rho, ux, uy, feq[NQ]; + stream_pull_load(k, f, in + base, j); + compute_rho_u(f, rho, ux, uy); + compute_feq(rho, ux, uy, feq); + collide_srt_no_force(f, feq, d_params.omega); + } + stream_pull_store(k, f, out + base); } diff --git a/tests/validation/multigrid/test_block_stepper.py b/tests/validation/multigrid/test_block_stepper.py index c2fb6c1..24e5634 100644 --- a/tests/validation/multigrid/test_block_stepper.py +++ b/tests/validation/multigrid/test_block_stepper.py @@ -2,7 +2,11 @@ Frozen GPU gates (32x24, 8 steps, positive low-Mach non-equilibrium state): finite cores; bit-identical 1x1, 2x2 and 4x2 splits; max |G - FP64 oracle| -<= 3e-6. Mass/momentum drift is reported, not gated. +<= 3e-6. Static circle (centre (15.5, 11.5), R=4.3, on the 2x2 corner seam), same +case: halo writes and q<.5 halo donors asserted present; finite fluid / NaN +solid cores; bit-identical splits; max |G - FP64 BFL+pull/SRT| <= 3e-6 on +fluid cells, from the same FP32-packed link data and initial populations. +Mass/momentum drift is reported, not gated. """ from contextlib import nullcontext from pathlib import Path @@ -14,10 +18,12 @@ from unittest.mock import MagicMock, Mock, patch import numpy as np import pycuda.driver as cuda +from CelerisLab.body.geometry.circle import CircleGeometry +from CelerisLab.body.objects import SimObject from CelerisLab.config import LBMConfig from CelerisLab.cuda import compiler_v2 as compiler from CelerisLab.lbm import block_stepper as implementation -from CelerisLab.lbm.blocks import BlockLayout2D +from CelerisLab.lbm.blocks import NO_DONOR, BlockLayout2D from CelerisLab.lbm.block_stepper import PeriodicBlockStepper from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY @@ -59,6 +65,48 @@ def oracle(g, steps): return g +CENTER, RADIUS = (15.5, 11.5), 4.3 + + +def circle(center=CENTER, radius=RADIUS): + return SimObject(0, geometry=CircleGeometry(*center, radius), center=center, radius=radius) + + +def dense_flags(): + """FLUID everywhere, circle flags on its solids (as the stepper builds them).""" + mask = circle().get_flag_mask(NX, NY) + flags = np.full(NX * NY, 1, np.uint16) + flags[(mask & 2) != 0] = mask[(mask & 2) != 0] + return flags + + +def circle_oracle(g, steps, solid, links): + """FP64 periodic static linear BFL -> pull -> SRT on fluid; solids NaN.""" + g = g.astype(np.float64) + g[:, solid] = np.nan + omega = 1 / (3 * NU + .5) + rows = [(*divmod(int(link.fluid_idx), NX), int(link.dir), float(np.float32(link.q)), + int(link.fallback)) for link in links] + for _ in range(steps): + staged = g.copy() + for y, x, d, q, fallback in rows: + ey, ex, o = D2Q9_EY[d], D2Q9_EX[d], d + 1 if d % 2 else d - 1 + toward = g[d, y, x] + if fallback: + value = toward + elif q < .5: + value = 2 * q * toward + (1 - 2 * q) * g[d, (y - ey) % NY, (x - ex) % NX] + else: + value = toward / (2 * q) + (1 - 1 / (2 * q)) * g[o, y, x] + staged[o, (y + ey) % NY, (x + ex) % NX] = value + f = np.stack([np.roll(staged[q], (D2Q9_EY[q], D2Q9_EX[q]), axis=(0, 1)) for q in range(9)]) + rho = f.sum(axis=0) + ux, uy = (CX[:, None, None] * f).sum(axis=0) / rho, (CY[:, None, None] * f).sum(axis=0) / rho + g = f + omega * (feq(rho, ux, uy) - f) + g[:, solid] = np.nan + return g + + def totals(g): g = g.astype(np.float64) return g.sum(), (CX[:, None, None] * g).sum(), (CY[:, None, None] * g).sum() @@ -80,6 +128,17 @@ class CpuContracts(unittest.TestCase): with self.subTest(cfg=cfg.collision, layout=getattr(lay, "periodic", lay)): with self.assertRaises(ValueError): PeriodicBlockStepper(cfg, lay, context=self.forbidden_context()) + moving = circle() + moving.state.vx = .01 + build_links = CircleGeometry.build_cut_links + for objects in ([moving], [circle((1.5, 11.5), 2.0)], [circle(), circle()]): + with self.assertRaises(ValueError): + PeriodicBlockStepper(config(), layout, context=self.forbidden_context(), objects=objects) + with patch.object(CircleGeometry, "build_cut_links", + lambda geometry, *a, **k: build_links(geometry, *a, **k)[1:]): + with self.assertRaisesRegex(ValueError, "cover every fluid pull"): + PeriodicBlockStepper(config(), layout, context=self.forbidden_context(), + objects=[circle()]) with patch.object(implementation, "MAX_GRID_DIM", 3): with self.assertRaisesRegex(ValueError, "launch limits"): PeriodicBlockStepper(config(), layout, context=self.forbidden_context()) @@ -107,6 +166,16 @@ class CpuContracts(unittest.TestCase): alloc.assert_not_called() upload.assert_not_called() + def test_static_circle_tables(self): + layout = BlockLayout2D.build(NY, NX, 12, 16) + tables = implementation._static_circles(layout, (circle(),)) + links = CircleGeometry(*CENTER, RADIUS).build_cut_links(NX, NY, domain_flags=dense_flags()) + self.assertEqual(len(tables["cell"]), len(links)) + np.testing.assert_array_equal(tables["q"], np.array([l.q for l in links], np.float32)) + solid_cells = layout.owner[tables["solid"]].astype(np.intp) + self.assertTrue((tables["flags"][solid_cells] & 2).all()) + self.assertEqual(int(np.count_nonzero(tables["flags"] & 2)), int(tables["solid"].sum())) + def fake_backend(self): backend = PeriodicBlockStepper.__new__(PeriodicBlockStepper) backend.layout = BlockLayout2D.build(NY, NX, 12, 16) @@ -117,6 +186,7 @@ class CpuContracts(unittest.TestCase): backend._stream = Mock() backend._module = Mock() backend._copy = backend._step = Mock() + backend.boundary, backend._curved, backend._tables, backend._links = None, None, {}, 0 backend._dst, backend._src = Mock(), Mock() backend._block, backend._grid, backend._copy_grid = (32, 1, 1), (1, 1, 1), (1, 1) backend._allocations = [Mock(), Mock()] @@ -211,6 +281,46 @@ class PeriodicBlockGPU(unittest.TestCase): *(a - b for a, b in zip(after, before))), flush=True) self.assertLessEqual(error, ORACLE_ATOL) + def test_static_circle_across_block_seams(self): + from CelerisLab.cuda.context import CudaContext + flags = dense_flags() + solid = ((flags & 2) != 0).reshape(NY, NX) + links = CircleGeometry(*CENTER, RADIUS).build_cut_links(NX, NY, domain_flags=flags) + for by, bx in SPLITS[1:]: # fixture must exercise receiver-local halo staging + layout = BlockLayout2D.build(NY, NX, by, bx) + mapped = layout.map_cut_links(links) + cells = int(np.prod(layout.local_shape)) + core = set(layout.owner.ravel().tolist()) + writes = {int(w) // (9 * cells) * cells + int(w) % cells for w in mapped.write} + self.assertTrue(writes - core, (by, bx)) + self.assertTrue({int(d) for d in mapped.donor if d != NO_DONOR} - core, (by, bx)) + g0 = initial_state() + results = [] + with CudaContext(0) as context: + for by, bx in SPLITS: + layout = BlockLayout2D.build(NY, NX, by, bx) + with PeriodicBlockStepper(config(), layout, context=context, objects=[circle()]) as backend: + for key, values, dtype in (("q", [l.q for l in links], np.float32), + ("dir", [l.dir for l in links], np.uint8), + ("fallback", [l.fallback for l in links], np.uint8)): + np.testing.assert_array_equal(backend.boundary[key], np.array(values, dtype)) + backend.initialize(g0) + backend.step(STEPS) + results.append(backend.download()) + for (by, bx), g in zip(SPLITS, results): + self.assertTrue(np.isfinite(g[:, ~solid]).all(), (by, bx)) + self.assertTrue(np.isnan(g[:, solid]).all(), (by, bx)) + np.testing.assert_array_equal(g, results[0], err_msg=f"split {(by, bx)}") + expected = circle_oracle(g0, STEPS, solid, links) + error = np.abs(results[0][:, ~solid].astype(np.float64) - expected[:, ~solid]).max() + start = g0.copy() + start[:, solid] = 0 + before, after = totals(start), totals(np.nan_to_num(results[0])) + print(f"circle max_abs_oracle_error={error:.3e} gate={ORACLE_ATOL:.0e} links={len(links)}", flush=True) + print("circle drift mass={:.3e} jx={:.3e} jy={:.3e}".format( + *(a - b for a, b in zip(after, before))), flush=True) + self.assertLessEqual(error, ORACLE_ATOL) + if __name__ == "__main__": unittest.main()