From 5b4610441c09a0671823065943a326972e041e70 Mon Sep 17 00:00:00 2001 From: TRON Date: Sat, 17 Oct 2026 23:46:50 -0700 Subject: [PATCH] feat: add matched SRT channel closure for blocked AMR root --- docs/development/root_channel_candidate.md | 22 +++ src/CelerisLab/lbm/channel.py | 54 ++++++++ src/CelerisLab/lbm/kernels/blocks.cu | 3 + .../lbm/kernels/boundary/inlet/common.cuh | 14 +- .../lbm/kernels/step/block_channel.cu | 66 +++++++++ src/CelerisLab/lbm/refinement_stepper.py | 45 ++++-- tests/unit/cpp/channel_root.cpp | 84 +++++++++++ tests/unit/test_channel_root.py | 99 +++++++++++++ .../multigrid/test_channel_root_gpu.py | 131 ++++++++++++++++++ 9 files changed, 504 insertions(+), 14 deletions(-) create mode 100644 docs/development/root_channel_candidate.md create mode 100644 src/CelerisLab/lbm/channel.py create mode 100644 src/CelerisLab/lbm/kernels/step/block_channel.cu create mode 100644 tests/unit/cpp/channel_root.cpp create mode 100644 tests/unit/test_channel_root.py create mode 100644 tests/validation/multigrid/test_channel_root_gpu.py diff --git a/docs/development/root_channel_candidate.md b/docs/development/root_channel_candidate.md new file mode 100644 index 0000000..c8a188f --- /dev/null +++ b/docs/development/root_channel_candidate.md @@ -0,0 +1,22 @@ +# Blocked AMR root channel increment + +Opt in with `FixedPatchStepper(..., coarse_boundary="channel", base_block=(by,bx), fine_block=(fy,fx))`, or the existing level-list form. Every level must be blocked. The supported boundary is the existing uniform D2Q9/SRT/FP32 double-buffer channel: local Zou-He west inlet, pressure-NEQ east outlet, stationary halfway bounce-back y walls. Uniform and parabolic inlet profiles and shifted/unshifted storage are supported; the existing configuration guards still reject MRT/LES and incompatible storage. This does not change the external DFG MRT reference. + +`step/block_channel.cu` is a small root-only kernel. Fine kernels, transfer kernels and recursive hop scheduling are unchanged. Global coordinates and the global profile height/velocity are separate from block-local storage strides. Corner ownership and collision phases match `OneStep`: + +- y-edge ghosts (including corners): pull, pair swap, no collision +- west ghosts at interior y: Zou-He reconstruction then SRT collision +- east ghosts at interior y: pressure NEQ from the immutable left input node, no collision +- physical fluid beside y walls: same-node opposite-input correction then collision + +The root uses the uniform SRT equilibrium arithmetic, including the existing shifted rho-minus-one path, to preserve matched-root semantics. Fine `BlockStep` retains its existing shifted direct-momentum/retained-density arithmetic. This increment is not a migration of fine arithmetic. + +Root halo topology stays wrapped like the uniform kernel, solely to reproduce its complete evolved ghost state as well as physical populations. It does not make the physical inlet/outlet periodic. Evolved root state includes boundary ghosts, while public fluid-leaf downloads exclude them. Initialization must provide finite positive-density root ghost states. Private raw snapshots preserve them exactly for same-representation continuation; public leaf-only downloads are not restart states. Channel remapping remains explicitly unsupported. + +`channel.py` checks the actual explosion donors and restriction receivers are uncovered bulk fluid and checks actual receiver-local restriction slots against pressure-donor and wall-correction reads. No arbitrary large padding or block-aligned patch restrictions are added. Keep the physical domain dimensions exact: choose legal root block divisors (including full width) instead of padding a benchmark domain. Stationary circles/Jets retain the existing finest-level and transfer-footprint rules. + +## Gates + +- `tests/unit/test_channel_root.py`: host semantics, minimal footprint, overlap rejection, raw-state/public-leaf separation; actual production uniform and blocked kernels under a CPU launch shim, five non-equilibrium steps, all corners, profiles, representations, y seams and a one-cell-wide block outlet donor. +- `tests/validation/multigrid/test_channel_root_gpu.py`: exact uniform-versus-blocked full-state GPU parity; AMR finite-state and block-split invariance with/without a stationary circle and velocity Jet; amplitude update and raw snapshot round-trip. +- GPU tests require explicit device admission by the caller. CPU-header algebra is not CUDA execution or external physical validation. The current increment does not establish DFG drag/lift/pressure acceptance or a performance speedup; those need matched SRT controls and a separate benchmark run. diff --git a/src/CelerisLab/lbm/channel.py b/src/CelerisLab/lbm/channel.py new file mode 100644 index 0000000..8915f0c --- /dev/null +++ b/src/CelerisLab/lbm/channel.py @@ -0,0 +1,54 @@ +"""Host contract for the blocked root's uniform-SRT channel boundary. + +Boundary ghosts are evolved state, but are not physical fluid leaves. Halo +ownership remains periodic to reproduce the uniform kernel's wrapped pulls. +No boundary closure is applied on a fine level. +""" +import numpy as np + + +def channel_masks(shape): + """Physical fluid and bulk-transfer masks in GLOBAL root coordinates.""" + ny, nx = shape + fluid = np.zeros(shape, bool) + fluid[1:ny - 1, 1:nx - 1] = True + bulk = fluid.copy() + bulk[[1, ny - 2]] = False + return fluid, bulk + + +def validate_channel_config(cfg, blocks): + if any(block is None for block in blocks): + raise ValueError("channel requires all-blocked levels") + if (cfg.inlet_scheme != "zou_he_local" or cfg.outlet_mode != "neq_extrap" + or cfg.y_wall_bc != "bounce_back" or cfg.top_wall_bc != "velocity" + or cfg.wall_slot_width != 0 or cfg.wall_normal_velocity != 0): + raise ValueError("channel requires Zou-He inlet, pressure NEQ outlet and stationary bounce-back y walls") + + +def validate_channel_transfers(pair, layout, targets): + """Reject actual exchange/closure overlap, not a padded geometric envelope. + + Explosion donors and restriction receivers must use ordinary bulk SRT + semantics. Restriction writes receiver-local source slots; these may not + overwrite the channel's extra immutable-input reads (pressure donor and + same-node halfway-wall corrections). The slot check deliberately respects + block seams, rather than treating different halo images as the same slot. + """ + ny, nx = pair.coarse_shape + _, bulk = channel_masks(pair.coarse_shape) + legal = bulk.ravel() & ~pair.covered.ravel() + receivers = pair.coalesce_coarse % (ny * nx) + if not legal[pair.ghost_coarse].all() or not legal[receivers].all(): + raise ValueError("channel transfer donors and receivers must be uncovered bulk fluid") + owner = layout.owner + reads = [] + # Outlet closes from the full previous-step state at its left neighbor. + donor = owner[1:-1, -1].astype(np.int64) - 1 + for q in range(9): + reads.append(layout.slot(donor, q)) + for y, directions in ((1, (4, 6, 7)), (ny - 2, (3, 5, 8))): + for q in directions: + reads.append(layout.slot(owner[y, 1:-1], q)) + if np.intersect1d(targets, np.concatenate(reads)).size: + raise ValueError("channel restriction writes overlap immutable boundary closure reads") diff --git a/src/CelerisLab/lbm/kernels/blocks.cu b/src/CelerisLab/lbm/kernels/blocks.cu index 8235658..add72db 100644 --- a/src/CelerisLab/lbm/kernels/blocks.cu +++ b/src/CelerisLab/lbm/kernels/blocks.cu @@ -157,3 +157,6 @@ extern "C" __global__ void BlockStep(const fpxx* in, fpxx* out, unsigned int hal } stream_pull_store(k, f, out + base); } + +// Root channel closure; fine and periodic levels keep BlockStep. +#include "step/block_channel.cu" diff --git a/src/CelerisLab/lbm/kernels/boundary/inlet/common.cuh b/src/CelerisLab/lbm/kernels/boundary/inlet/common.cuh index cd37f5d..fab9a69 100644 --- a/src/CelerisLab/lbm/kernels/boundary/inlet/common.cuh +++ b/src/CelerisLab/lbm/kernels/boundary/inlet/common.cuh @@ -27,18 +27,22 @@ #define INLET_REGULARIZED_NEQ_DAMP 0.50f #endif -__device__ __forceinline__ float inlet_target_u(float y_coord) { +__device__ __forceinline__ float inlet_target_u(float y_coord, unsigned int ny, float velocity) { #if INLET_PROFILE == 0 - return U0; + return velocity; #else - const float y_clamped = fminf((float)(NY - 2), fmaxf(1.0f, y_coord)); - const float H = fmaxf((float)(NY - 2), 1.0f); + const float y_clamped = fminf((float)(ny - 2), fmaxf(1.0f, y_coord)); + const float H = fmaxf((float)(ny - 2), 1.0f); const float eta = (y_clamped - 0.5f) / H; const float shape = fmaxf(0.0f, 4.0f * eta * (1.0f - eta)); - return U0 * 1.5f * shape; + return velocity * 1.5f * shape; #endif } +__device__ __forceinline__ float inlet_target_u(float y_coord) { + return inlet_target_u(y_coord, NY, U0); +} + #if LATTICE_MODEL == LATTICE_D2Q9 __device__ __forceinline__ float west_velocity_rho_closure_d2q9( const float* __restrict__ f, diff --git a/src/CelerisLab/lbm/kernels/step/block_channel.cu b/src/CelerisLab/lbm/kernels/step/block_channel.cu new file mode 100644 index 0000000..4e95674 --- /dev/null +++ b/src/CelerisLab/lbm/kernels/step/block_channel.cu @@ -0,0 +1,66 @@ +// Root-only blocked counterpart of OneStep's D2Q9/SRT channel closure. +// NX/NY are LOCAL strides; all physical boundary decisions use global coordinates. +#include "../boundary/inlet/common.cuh" +#include "../boundary/inlet/zou_he_local.cuh" +#include "../boundary/outlet/pressure_neq.cuh" + +struct BlockChannel { unsigned int nx, ny; float velocity; }; +extern "C" { __constant__ BlockChannel d_block_channel; } + +// Same signature as BlockStep so the AMR hop/transfer schedule is unchanged. +extern "C" __global__ void BlockChannelStep( + const fpxx* in, fpxx* out, unsigned int halo, + const uint16_t* flags, const unsigned char* kind) +{ + 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 int bx = NX - 2 * halo, by = NY - 2 * halo; + const unsigned int nbx = d_block_channel.nx / bx; + const unsigned int gx = (blockIdx.z % nbx) * bx + x - halo; + const unsigned int gy = (blockIdx.z / nbx) * by + y - halo; + const unsigned long base = (unsigned long)blockIdx.z * NQ * TOTAL_CELLS; + const unsigned long cell = (unsigned long)blockIdx.z * TOTAL_CELLS + k; + float f[NQ]; + if (kind != nullptr && kind[cell] == KIND_OUTER) { + for (int q = 0; q < NQ; ++q) f[q] = NAN; + } else { + unsigned long j[NQ]; + compute_neighbors(x, y, j); + stream_pull_load(k, f, in + base, j); + const bool wall = gy == 0 || gy == d_block_channel.ny - 1; + const bool outlet = gx == d_block_channel.nx - 1; + if (wall) { + // Wall owns all four corners. Ghosts pull and swap, never collide. + for (int q = 1; q < NQ; q += 2) { + const float t = f[q]; f[q] = f[q + 1]; f[q + 1] = t; + } + } else if (gx == 0) { + apply_zou_he_left_velocity_inlet_d2q9( + f, inlet_target_u((float)gy, d_block_channel.ny, d_block_channel.velocity), 0.0f); + } else if (outlet) { + // Immutable input donor, receiver-local (also valid across a block seam). + float donor[NQ]; + for (int q = 0; q < NQ; ++q) + donor[q] = load_ddf(in + base, index_f(k - 1, q)); + apply_pressure_outlet_d2q9(f, donor); + } else if (gy == 1) { + f[3] = load_ddf(in + base, index_f(k, 4)); + f[5] = load_ddf(in + base, index_f(k, 6)); + f[8] = load_ddf(in + base, index_f(k, 7)); + } else if (gy == d_block_channel.ny - 2) { + f[4] = load_ddf(in + base, index_f(k, 3)); + f[6] = load_ddf(in + base, index_f(k, 5)); + f[7] = load_ddf(in + base, index_f(k, 8)); + } + if (!wall && !outlet) { + // Match uniform SRT arithmetic, including its shifted equilibrium path. + float rho, ux, uy, feq[NQ], forcing[NQ] = {}; + compute_rho_u(f, rho, ux, uy); + compute_feq(rho, ux, uy, feq); + const float omega = fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, d_params.omega)); + collide_srt(f, feq, forcing, omega); + } + } + stream_pull_store(k, f, out + base); +} diff --git a/src/CelerisLab/lbm/refinement_stepper.py b/src/CelerisLab/lbm/refinement_stepper.py index 0843349..143c805 100644 --- a/src/CelerisLab/lbm/refinement_stepper.py +++ b/src/CelerisLab/lbm/refinement_stepper.py @@ -16,6 +16,7 @@ import pycuda.driver as cuda from ..cuda import compiler_v2 as compiler from .blocks import BlockLayout2D +from .channel import channel_masks, validate_channel_config, validate_channel_transfers from ..body.jets import jet_link_velocity from .block_stepper import MAX_GRID_DIM from .hierarchy import build_hierarchy, is_rectangle @@ -135,6 +136,11 @@ class FixedPatchStepper: ``coarse_boundary`` is ``periodic`` (default) or ``slot_reservoir``: a bottom on-site slot and top fixed-density reservoir from ``cfg``. + ``channel`` requires all-blocked levels and matches uniform SRT Zou-He + inlet, pressure-NEQ outlet and stationary bounce-back y walls. Boundary + ghosts evolve internally but public downloads exclude them; initialization + must supply their finite positive-density state. Transfers must avoid the + actual channel closure footprint. Channel remapping is not supported. Periodic mode does not apply the wall/slot boundary fields in ``cfg``. Outer walls remain coarse; fine halo donors must be ordinary fluid. ``fine_objects`` accepts fixed circles in coarse physical coordinates; @@ -213,8 +219,8 @@ class FixedPatchStepper: observe_raw_exchange=False): 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") + if coarse_boundary not in ("periodic", "slot_reservoir", "channel"): + raise ValueError("coarse_boundary must be periodic, slot_reservoir or channel") if refinements is None: if bounds is None: raise ValueError("bounds or refinements is required") @@ -242,6 +248,9 @@ class FixedPatchStepper: 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") + self._channel = coarse_boundary == "channel" + if self._channel: + validate_channel_config(cfg, level_blocks) self._wall_aware = coarse_boundary == "slot_reservoir" if self._wall_aware and (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): @@ -268,6 +277,9 @@ class FixedPatchStepper: if layout.nblocks > MAX_GRID_DIM or layout.local_shape[0] > MAX_GRID_DIM: raise ValueError(f"{_level_name(level)} blocks exceed launch limits") self._layouts.append(layout) + if self._channel: + base = blocked_base_tables(pairs[0], self._layouts[0], np.argsort(pairs[0].coalesce_coarse)) + validate_channel_transfers(pairs[0], self._layouts[0], base[1]) # Kernel configurations per level (geometry configs live in self.hierarchy). fine_cfg = replace(cfg, nx=pairs[0].fine_shape[1], ny=pairs[0].fine_shape[0], viscosity=2 * cfg.viscosity) @@ -318,6 +330,9 @@ class FixedPatchStepper: host = self._boundary_hosts[level] if host is not None: self._leaf[level] &= (host.flags.reshape(pairs[level - 1].fine_shape) & SOLID) == 0 + self._state = [mask.copy() for mask in self._leaf] + if self._channel: + self._leaf[0] &= channel_masks((cfg.ny, cfg.nx))[0] # Private constructor inputs for the bounded manual L2 continuation. self._remap_cfg = replace(cfg) self._level_blocks = level_blocks @@ -384,7 +399,7 @@ class FixedPatchStepper: 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" + name, source = ("BlockChannelStep" if index == 0 and self._channel else "BlockStep"), "blocks.cu" config.validate() if not np.isclose(config.omega, 1 / tau, rtol=1e-14, atol=0): raise RuntimeError("level omega/tau mismatch") @@ -412,6 +427,12 @@ class FixedPatchStepper: 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._channel: + cfg = self._remap_cfg + ptr, size = module.get_global("d_block_channel") + if size != struct.calcsize("=IIf"): + raise RuntimeError("Block channel host/device layout mismatch") + cuda.memcpy_htod(ptr, struct.pack("=IIf", cfg.nx, cfg.ny, cfg.velocity)) if index == 0 and self._wall_aware: cfg = self._remap_cfg ptr, size = module.get_global("d_block_wall") @@ -475,6 +496,10 @@ class FixedPatchStepper: for array, leaf in zip(arrays, self._leaf): array[:, ~leaf] = np.nan + def _mask_nonstate(self, *arrays): + for array, state in zip(arrays, self._state): + array[:, ~state] = np.nan + def _check_storage(self): if any(bool(c.ddf_shifting) != self._ddf_shifting for c in self._level_cfgs): raise ValueError("mixed population representations are unsupported") @@ -497,7 +522,7 @@ class FixedPatchStepper: arrays = [np.array(a, dtype=dtype, order="C", copy=True) for a in given] if any(a.shape != level.shape for a, level in zip(arrays, self._levels)): raise ValueError("Expected Q-major (9,ny,nx) arrays for every level") - for array, mask in zip(arrays, self._leaf): + for array, mask in zip(arrays, self._state): 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") encoded = tuple(_encode_physical(a, self._ddf_shifting) for a in arrays) @@ -518,12 +543,12 @@ class FixedPatchStepper: arrays = [np.array(a, order="C", copy=True) for a in snapshot.values] if any(a.shape != level.shape for a, level in zip(arrays, self._levels)): raise ValueError("Expected Q-major (9,ny,nx) arrays for every level") - for array, mask in zip(arrays, self._leaf): + for array, mask in zip(arrays, self._state): values = array[:, mask] if (not np.isfinite(values).all() or (values.sum(axis=0, dtype=np.float64) + int(self._ddf_shifting) <= 0).any()): raise ValueError("Leaf populations must be finite with positive density") - self._mask_nonleaves(*arrays) + self._mask_nonstate(*arrays) arrays = [a if level.layout is None else level.layout.to_blocked(a) for a, level in zip(arrays, self._levels)] try: @@ -562,7 +587,7 @@ class FixedPatchStepper: self._require("ready") self._check_storage() pairs = self.hierarchy.pairs - if (self._wall_aware or self.hierarchy.depth != 3 + if (self._wall_aware or self._channel or self.hierarchy.depth != 3 or not isinstance(pairs[1], FixedPatchSet2D) or self._boundary_hosts[1] is not None or any(b is not None for b in self._level_boundaries[:2])): @@ -820,7 +845,9 @@ class FixedPatchStepper: def download(self): """Return physical G (FP64 for shifted storage), NaN on non-leaves.""" snapshot = self._download_raw() - return tuple(_decode_stored(a, snapshot.shifted) for a in snapshot.values) + arrays = tuple(_decode_stored(a, snapshot.shifted) for a in snapshot.values) + self._mask_nonleaves(*arrays) + return arrays def _download_raw(self): """Current owned populations with a representation tag; no add/subtract w.""" @@ -837,7 +864,7 @@ class FixedPatchStepper: raise 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) + self._mask_nonstate(*arrays) return _StoredPopulations(tuple(arrays), self._ddf_shifting) def close(self): diff --git a/tests/unit/cpp/channel_root.cpp b/tests/unit/cpp/channel_root.cpp new file mode 100644 index 0000000..3d7cd21 --- /dev/null +++ b/tests/unit/cpp/channel_root.cpp @@ -0,0 +1,84 @@ +// Actual production step kernels under a CPU launch shim. No CUDA claim. +#include +#include +#include +#include +#define __device__ +#define __constant__ +#define __global__ +#define __forceinline__ inline +#define __launch_bounds__(...) +#define LATTICE_D2Q9 0 +#define LATTICE_D3Q19 1 +#define LATTICE_MODEL 0 +#define NQ 9 +#define DIM 2 +#define NZ 1 +#define NT 1 +#define COLLISION_MODEL 0 +#define USE_LES 0 +#define RHO 1.f +#define U0 .025f +#define OMEGA_COLLISION_MIN .01f +#define OMEGA_COLLISION_MAX 1.99f +#define Y_WALL_BC 0 +#define OUTLET_MODE 0 +struct Dim {unsigned int x=0,y=0,z=0;} threadIdx,blockIdx,blockDim; +inline float __fdiv_rn(float a,float b){return a/b;} +#include "core/precision.cuh" +#include "core/descriptors.cuh" +#include "core/layout.cuh" +#include "core/params.cuh" +#include "core/flags.cuh" +#include "operators/macro.cuh" +#include "operators/collision_srt.cuh" +#include "streaming/pull_double_buffer.cuh" +#if UNIFORM +#include "operators/forcing_guo.cuh" +#include "operators/helpers.cuh" +#include "boundary/bounce_back.cuh" +#include "boundary/inlet_outlet.cuh" +#include "step/one_step_double.cu" +#else +constexpr unsigned char KIND_OUTER=2; +#include "step/block_channel.cu" +#endif +int main(){ + constexpr int gx=24,gy=20,n=gx*gy; + std::vector state(9*n),next(9*n); + for(int q=0;q<9;q++)for(int y=0;y flags(n,FLAG_FLUID); + for(int y=0;y input(nbx*nby*9*TOTAL_CELLS),output(input.size()); + d_block_channel={gx,gy,U0}; +#endif + for(int step=0;step<5;step++){ +#if UNIFORM + for(int y=0;y