feat: add matched SRT channel closure for blocked AMR root

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")

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"

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,

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);
+}

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):