feat: add experimental projected explosion per refinement pair

src/CelerisLab/lbm/kernels/blocks.cu

@@ -63,6 +63,39 @@ extern "C" __global__ void BlockExplode(
         to[index_f(t % TOTAL_CELLS, q)] = from[q * coarse_cells + d % coarse_cells];
 }
 
+// Experimental Schornbaum--Ruede (arXiv:1508.07982, Algorithm 2)
+// and Chen et al. (doi:10.1016/j.physa.2005.09.036) population interpolation.
+// Independent implementation of the mathematical stencil; no upstream source copied.
+// Stencil cells are owner cores; a disabled axis repeats its center twice.
+extern "C" __global__ void BlockExplodeProjected(
+    const fpxx* coarse, fpxx* fine, const uint64_t* target, const uint64_t* donor,
+    unsigned int count, unsigned long coarse_cells,
+    const uint64_t* neighbors, const float* offsets)
+{
+    unsigned int row = blockIdx.x * blockDim.x + threadIdx.x;
+    if (row >= count) return;
+    const uint64_t t = target[row], d = donor[row];
+    fpxx* to = fine + (t / TOTAL_CELLS) * NQ * TOTAL_CELLS;
+    for (unsigned int q = 0; q < NQ; ++q) {
+        float v[4];
+        for (unsigned int a = 0; a < 4; ++a) {
+            const uint64_t n = neighbors[4ul * row + a];
+            v[a] = coarse[(n / coarse_cells) * NQ * coarse_cells
+                          + q * coarse_cells + n % coarse_cells];
+        }
+        float gx = .5f * (v[1] - v[0]), gy = .5f * (v[3] - v[2]);
+        const float ex = d_cx[q], ey = d_cy[q];
+        if (q != 0) {
+            const float along = (ex * gx + ey * gy) / (ex * ex + ey * ey);
+            gx -= ex * along;
+            gy -= ey * along;
+        }
+        to[index_f(t % TOTAL_CELLS, q)] =
+            coarse[(d / coarse_cells) * NQ * coarse_cells + q * coarse_cells + d % coarse_cells]
+            + offsets[2ul * row] * gx + offsets[2ul * row + 1] * gy;
+    }
+}
+
 // 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(

src/CelerisLab/lbm/refinement_stepper.py

@@ -22,7 +22,7 @@ from .block_stepper import MAX_GRID_DIM
 from .hierarchy import build_hierarchy, is_rectangle
 from .refinement import ACTIVE, HALO, FixedPatch2D, FixedPatchSet2D, prepare_fine_boundary
 from .refinement_tables import (blocked_base_tables, blocked_fine_tables,  # noqa: F401 (re-export)
-                                blocked_kind, _prepare_blocked_boundary,
+                                blocked_kind, projected_explosion_tables, _prepare_blocked_boundary,
                                 BLOCKED_BOUNDARY_DEVICE_COLUMNS)
 from .descriptors import SOLID
 from .raw_exchange import summarize_exchange
@@ -41,6 +41,8 @@ class _Pair:
     n_ghosts: int
     n_rows: int
     parent_cells: int    # parent block window cells (explode donor stride)
+    neighbors: object = None  # projected-only owner-core cells; constant allocates neither
+    offsets: object = None
 
 
 @dataclass
@@ -206,6 +208,12 @@ class FixedPatchStepper:
     halo -> BFL -> step four times per base step. Runtime child Jet changes
     require ``set_jet_amplitude(index, amplitude, level=2)``.
 
+    ``explosion="projected_linear"`` selects experimental transverse linear
+    reconstruction on each adjacent pair. Default "constant" is unchanged.
+    Projected mode requires periodic root, all-blocked SRT levels without
+    bodies/Jets or LES; this backend supplies no forcing. Geometry is generic,
+    but union/multilevel GPU and matched-physics qualification remain pending.
+
     Level lists (any depth): ``refinements[l]`` gives level l+1 in level-l
     integer cells (scalar window or rectangle list, see lbm.hierarchy);
     ``level_blocks``/``level_objects``/``level_jets`` hold one entry per level,
@@ -219,7 +227,7 @@ class FixedPatchStepper:
                  fine_block=None, base_block=None, jets=(), child=None, child_block=None,
                  child_objects=(), child_jets=(), child_rectangles_in_parent=None,
                  refinements=None, level_blocks=None, level_objects=None, level_jets=None,
-                 observe_raw_exchange=False, observe_window_exchange=False):
+                 observe_raw_exchange=False, observe_window_exchange=False, explosion="constant"):
         if HALO != 4:
             raise RuntimeError("refinement.cu requires a four-cell fine halo")
         if coarse_boundary not in ("periodic", "slot_reservoir", "channel"):
@@ -248,6 +256,15 @@ class FixedPatchStepper:
         level_blocks = [None if b is None else tuple(b) for b in level_blocks]
         _check_levels(refinements, level_blocks, level_objects, level_jets)
         depth = len(refinements) + 1
+        if explosion not in ("constant", "projected_linear"):
+            raise ValueError("explosion must be constant or projected_linear")
+        if explosion == "projected_linear" and (
+                coarse_boundary != "periodic" or any(b is None for b in level_blocks)
+                or any(level_objects) or any(level_jets)
+                or cfg.collision != "SRT" or cfg.les_enabled):
+            raise ValueError("projected explosion requires periodic all-blocked SRT levels "
+                             "without bodies/Jets or LES")
+        self._explosion = explosion
         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")
@@ -479,6 +496,11 @@ class FixedPatchStepper:
                        self._upload(children), len(ghosts), len(targets),
                        int(np.prod(pair.coarse_shape if parent.layout is None
                                    else parent.layout.local_shape)))
+        if self._explosion == "projected_linear":
+            neighbors, offsets = projected_explosion_tables(
+                pair, parent.layout, self.hierarchy.kinds[index])
+            result.explode = child.module.get_function("BlockExplodeProjected")
+            result.neighbors, result.offsets = self._upload(neighbors), self._upload(offsets)
         if child.layout is not None:
             child.kind = self._upload(blocked_kind(child.layout, self.hierarchy.kinds[index + 1]))
         if index == 0 and base is not None:
@@ -653,7 +675,7 @@ class FixedPatchStepper:
                 child_rectangles_in_parent=new.child.bounds,
                 child_objects=child_objects, child_jets=jets,
                 observe_raw_exchange=self._observe_raw_exchange,
-                observe_window_exchange=self._observe_window)
+                observe_window_exchange=self._observe_window, explosion=self._explosion)
             try:
                 if self._ddf_shifting:
                     candidate._initialize_raw(_StoredPopulations(transferred.g, transferred.shifted))
@@ -746,6 +768,7 @@ class FixedPatchStepper:
         pair, level = self._pairs[index - 1], self._levels[index]
         pair.explode(parent_buffer, level.current, pair.ghosts, pair.donors, np.uint32(pair.n_ghosts),
                      np.uint64(pair.parent_cells),
+                     *((pair.neighbors, pair.offsets) if pair.neighbors is not None else ()),
                      grid=((pair.n_ghosts + self._block[0] - 1) // self._block[0], 1), **kw)
         self._hop(index, level.current, level.spare, **kw)
         self._hop(index, level.spare, level.current, **kw)

src/CelerisLab/lbm/refinement_tables.py

@@ -66,6 +66,50 @@ def blocked_base_tables(patch, layout, order, wall_aware=False, *, composed_kind
             blocked_kind(layout, dense_kind))
 
 
+def projected_explosion_tables(pair, parent_layout, composed_parent_kind):
+    """Experimental transverse linear explosion, in the pair's existing row order.
+
+    Schornbaum & Rüde (2016), arXiv:1508.07982, Algorithm 2, following
+    Chen et al., doi:10.1016/j.physa.2005.09.036. Independently implemented
+    from the mathematics: centered gradients, perpendicular projection and
+    child offsets +/-1/4. The kernel applies the projection per population.
+
+    Only original explosion donors that are owned, uncovered ACTIVE cells
+    may supply gradients. Both axis neighbors must be legal; otherwise repeat
+    the center twice for zero slope. Owner-core addresses avoid stale halos.
+    The same rule applies to scalar windows, unions and every adjacent pair;
+    nested layouts retain their own per-axis (usually nonperiodic) topology.
+    """
+    ny, nx = pair.coarse_shape
+    kind = np.asarray(composed_parent_kind)
+    if parent_layout.shape != pair.coarse_shape or kind.shape != pair.coarse_shape:
+        raise ValueError("projected explosion parent layout/kind has the wrong shape")
+    donor = pair.ghost_coarse.astype(np.intp)
+    legal = np.zeros(ny * nx, bool)
+    legal[donor] = True
+    legal &= (parent_layout.owned.ravel() & ~pair.covered.ravel()
+              & (kind.ravel() == ACTIVE))
+    if not legal[donor].all():
+        raise ValueError("projected explosion donors need owned uncovered parent ACTIVE cells")
+    y, x = np.divmod(donor, nx)
+    neighbors = np.repeat(donor[:, None], 4, axis=1)
+    for axis, (dx, dy) in enumerate(((1, 0), (0, 1))):
+        minus = ((y - dy) % ny) * nx + (x - dx) % nx
+        plus = ((y + dy) % ny) * nx + (x + dx) % nx
+        valid = legal[minus] & legal[plus]
+        if not parent_layout.periodic[axis]:
+            coordinate, extent = (x, nx) if axis == 0 else (y, ny)
+            valid &= (coordinate > 0) & (coordinate < extent - 1)
+        neighbors[valid, 2 * axis] = minus[valid]
+        neighbors[valid, 2 * axis + 1] = plus[valid]
+    fy, fx = np.divmod(pair.ghost_fine, pair.fine_shape[1])
+    ox, oy = pair.fine_origin
+    offsets = np.column_stack((ox + (fx + .5) / 2 - (x + .5),
+                               oy + (fy + .5) / 2 - (y + .5)))
+    return (np.ascontiguousarray(parent_layout.owner.ravel()[neighbors], dtype=np.uint64),
+            np.ascontiguousarray(offsets, dtype=np.float32))
+
+
 # Required device columns plus optional actuation; observation arms are uploaded
 # separately only when enabled. body_id/aperture remain host reduction metadata.
 BLOCKED_BOUNDARY_DEVICE_COLUMNS = ("block", "cell", "dir", "q", "fallback", "surface_velocity")