feat: opt-in MRT collision for blocked/refined stepping

Additive compile-time COLLISION_MODEL == 2 branches in BlockStep,
BlockChannelStep and patch_collide repeat the uniform collide_dispatch
MRT sequence (zero forcing, omega clamp, collide_mrt_no_force). SRT
statements, including the shifted retained-delta/direct-momentum path,
are unchanged under the alternate branch. Guards accept SRT or MRT only;
TRT, LES, non-FP32, non-double-buffer and shifted dense/mixed refinement
stay rejected. Transfers, subcycling, BFL, raw ledgers, channel closures,
shifted-Jet SRT guard and PeriodicBlockStepper are untouched. Existing
MRT operator, shifted representation and fixed kinetic rates are reused.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

src/CelerisLab/lbm/channel.py

@@ -1,4 +1,4 @@
-"""Host contract for the blocked root's uniform-SRT channel boundary.
+"""Host contract for the blocked root's uniform SRT/MRT 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.
@@ -29,7 +29,7 @@ def validate_channel_config(cfg, blocks):
 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
+    Explosion donors and restriction receivers must use ordinary bulk collision
     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

src/CelerisLab/lbm/kernels/blocks.cu

@@ -14,12 +14,15 @@
 #include "core/obs.cuh"
 #include "operators/macro.cuh"
 #include "operators/collision_srt.cuh"
+#if COLLISION_MODEL == 2
+#include "operators/collision_mrt.cuh"
+#endif
 #include "streaming/pull_double_buffer.cuh"
 #include "boundary/curved_boundary.cuh"
 #include "boundary/planar_velocity.cuh"
 
-#if DIM != 2 || NQ != 9 || COLLISION_MODEL != 0 || STORE_PRECISION != 0 || USE_LES || STREAMING_MODEL != 0
-#error "Blocks require FP32 D2Q9 SRT double buffers without LES"
+#if DIM != 2 || NQ != 9 || (COLLISION_MODEL != 0 && COLLISION_MODEL != 2) || STORE_PRECISION != 0 || USE_LES || STREAMING_MODEL != 0
+#error "Blocks require FP32 D2Q9 SRT or MRT double buffers without LES"
 #endif
 
 // Fixed-patch cell kinds, matching refinement.py ACTIVE/INNER_GHOST/OUTER_GHOST.
@@ -141,6 +144,13 @@ extern "C" __global__ void BlockStep(const fpxx* in, fpxx* out, unsigned int hal
                         gx, d_block_wall.nx, d_block_wall.slot_width, d_block_wall.amplitude));
                 else if (gy == d_block_wall.ny - 2) reconstruct_top_reservoir(f, RHO);
             }
+#if COLLISION_MODEL == 2
+            // Uniform MRT sequence (collide_dispatch, zero forcing); SRT path below unchanged.
+            float rho, ux, uy;
+            compute_rho_u(f, rho, ux, uy);
+            collide_mrt_no_force(f, rho, ux, uy,
+                fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, d_params.omega)));
+#else
             float rho, ux, uy, feq[NQ];
             compute_rho_u(f, rho, ux, uy);
 #if USE_DDF_SHIFTING
@@ -153,6 +163,7 @@ extern "C" __global__ void BlockStep(const fpxx* in, fpxx* out, unsigned int hal
             compute_feq(rho, ux, uy, feq);
 #endif
             collide_srt_no_force(f, feq, d_params.omega);
+#endif
         }
     }
     stream_pull_store(k, f, out + base);

src/CelerisLab/lbm/kernels/refinement.cu

@@ -11,6 +11,9 @@
 #include "core/obs.cuh"
 #include "operators/macro.cuh"
 #include "operators/collision_srt.cuh"
+#if COLLISION_MODEL == 2
+#include "operators/collision_mrt.cuh"
+#endif
 #include "streaming/pull_double_buffer.cuh"
 #include "boundary/planar_velocity.cuh"
 #include "boundary/curved_boundary.cuh"
@@ -19,8 +22,8 @@ extern "C" {
 #include "step/curved_double.cu"
 }
 
-#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"
+#if DIM != 2 || NQ != 9 || (COLLISION_MODEL != 0 && COLLISION_MODEL != 2) || STORE_PRECISION != 0 || USE_DDF_SHIFTING || USE_LES || STREAMING_MODEL != 0
+#error "Refinement requires unshifted FP32 D2Q9 SRT or MRT double buffers without LES"
 #endif
 
 // Enabled only in the validation build; no diagnostic arrays in normal builds.
@@ -33,10 +36,17 @@ __device__ unsigned int incoming_counts[NQ * TOTAL_CELLS];
 #endif
 
 __device__ __forceinline__ void patch_collide(float* f, unsigned long k) {
+#if COLLISION_MODEL == 2
+    float rho, ux, uy;
+    compute_rho_u(f, rho, ux, uy);
+    collide_mrt_no_force(f, rho, ux, uy,
+        fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, d_params.omega)));
+#else
     float rho, ux, uy, feq[NQ];
     compute_rho_u(f, rho, ux, uy);
     compute_feq(rho, ux, uy, feq);
     collide_srt_no_force(f, feq, d_params.omega);
+#endif
 #ifdef CELERIS_REFINEMENT_AUDIT
     ++collision_counts[k];
 #endif

src/CelerisLab/lbm/kernels/step/block_channel.cu

@@ -54,12 +54,20 @@ extern "C" __global__ void BlockChannelStep(
             f[7] = load_ddf(in + base, index_f(k, 8));
         }
         if (!wall && !outlet) {
+#if COLLISION_MODEL == 2
+            // Match uniform MRT arithmetic (collide_dispatch, zero forcing).
+            float rho, ux, uy;
+            compute_rho_u(f, rho, ux, uy);
+            collide_mrt_no_force(f, rho, ux, uy,
+                fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, d_params.omega)));
+#else
             // 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);
+#endif
         }
     }
     stream_pull_store(k, f, out + base);

src/CelerisLab/lbm/refinement.py

@@ -52,17 +52,17 @@ class FixedPatch2D:
     def from_config(cls, cfg: LBMConfig, bounds) -> "FixedPatch2D":
         """Prepare maps without changing cfg or allocating coarse-sized tables.
 
-        The initial backend contract is D2Q9/SRT/FP32/double_buffer,
+        The backend contract is D2Q9/SRT or MRT/FP32/double_buffer,
         without LES. Two coarse cells of margin keep all halo donors in bounds;
         this does NOT certify that those donors are fluid or periodic.
         """
         checked = replace(cfg)
         checked.validate()
-        if (not checked.is_d2q9 or checked.collision != "SRT"
+        if (not checked.is_d2q9 or checked.collision not in ("SRT", "MRT")
                 or checked.streaming != "double_buffer"
                 or checked.store_precision != "FP32"
                 or checked.les_enabled):
-            raise ValueError("FixedPatch2D requires D2Q9/SRT/FP32/double_buffer without LES")
+            raise ValueError("FixedPatch2D requires D2Q9/SRT or MRT/FP32/double_buffer without LES")
         bounds = tuple(bounds)
         if len(bounds) != 4 or any(isinstance(v, (bool, np.bool_))
                                    or not isinstance(v, Integral) for v in bounds):
@@ -451,7 +451,7 @@ def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=(),
     checked = replace(cfg)
     checked.validate()
     if (patch.coarse_shape != (cfg.ny, cfg.nx) or not checked.is_d2q9
-            or checked.collision != "SRT" or checked.streaming != "double_buffer"
+            or checked.collision not in ("SRT", "MRT") or checked.streaming != "double_buffer"
             or checked.store_precision != "FP32" or checked.les_enabled
             or patch.tau_coarse != 3*checked.viscosity+.5):
         raise ValueError("Fine boundary requires the matching fixed-patch configuration")