diff --git a/README.md b/README.md
index e25d9f2..7de48be 100644
--- a/README.md
+++ b/README.md
@@ -90,6 +90,7 @@ sim.run(100, stream=stream, upload_act=False, sync_obs=False)
 # Later: sync and read
 stream.synchronize()
 obs = sim.read_bodies(stream=stream)  # sync already done, just read pinned buffer
+```
 
 ### Fixed 2:1 refinement patch (experimental backend)
 
@@ -99,10 +100,12 @@ normal surface velocity (blowing/suction). It is a standalone backend, not part
 of `Simulation`. Runnable demo: `examples/fixed_patch_circle.py`.
 
 ```python
+import numpy as np
 from CelerisLab.body.geometry.circle import CircleGeometry
 from CelerisLab.body.objects import SimObject, SurfaceActuation
 from CelerisLab.config import LBMConfig
 from CelerisLab.cuda.context import CudaContext
+from CelerisLab.lbm.equilibrium import compute_feq_d2q9
 from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
 
 cfg = LBMConfig(nx=32, ny=24, viscosity=.1)
@@ -110,7 +113,10 @@ body = SimObject(0, geometry=CircleGeometry(16.0, 12.0, 2.15), center=(16.0, 12.
                  surface_actuation=SurfaceActuation(0.7, 0.85, 0.002))  # angle, half-width, speed
 with CudaContext(0) as ctx, FixedPatchStepper(cfg, (8, 6, 24, 18), context=ctx,
                                               fine_objects=[body], fine_block=(16, 20)) as run:
-    run.initialize(coarse_g, fine_g)   # (9, ny, nx) postcollision populations per level
+    rest = np.asarray(compute_feq_d2q9(1.0, 0.0, 0.0), np.float32)[:, None, None]
+    # (9, ny, nx) postcollision populations per level; non-leaf cells are discarded.
+    run.initialize(np.broadcast_to(rest, (9, *run.patch.coarse_shape)),
+                   np.broadcast_to(rest, (9, *run.patch.fine_shape)))
     run.step(128)                      # coarse steps; fine level takes two per step
     coarse_g, fine_g = run.download()  # NaN on covered, ghost and solid cells
 ```
diff --git a/examples/fixed_patch_circle.py b/examples/fixed_patch_circle.py
index ee0f034..940106f 100644
--- a/examples/fixed_patch_circle.py
+++ b/examples/fixed_patch_circle.py
@@ -19,8 +19,9 @@ from CelerisLab.body.geometry.circle import CircleGeometry
 from CelerisLab.body.objects import SimObject, SurfaceActuation
 from CelerisLab.config import LBMConfig
 from CelerisLab.cuda.context import CudaContext
-from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY
+from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY, SOLID
 from CelerisLab.lbm.equilibrium import compute_feq_d2q9
+from CelerisLab.lbm.refinement import ACTIVE
 from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
 
 NX, NY, BOUNDS = 32, 24, (8, 6, 24, 18)  # coarse cells; patch is half-open (x0, y0, x1, y1)
@@ -35,15 +36,29 @@ def config(mode):
                      wall_slot_width=8, y_wall_bc="on_site_velocity", top_wall_bc="reservoir")
 
 
-def leaf_diagnostics(coarse, fine):
-    """Fluid leaves only (downloads are NaN on covered/ghost/solid cells).
+def leaf_masks(run, mode):
+    """Known fluid leaves: uncovered coarse cells (minus wall rows in slot mode)
+    and fine ACTIVE cells outside the body."""
+    p = run.patch
+    x0, y0, x1, y1 = p.bounds
+    coarse = np.ones(p.coarse_shape, bool)
+    coarse[y0:y1, x0:x1] = False
+    if mode == "slot_reservoir":
+        coarse[[0, -1]] = False
+    solid = (run.fine_boundary.flags.reshape(p.fine_shape) & SOLID) != 0
+    return coarse, (p.fine_kind == ACTIVE) & ~solid
+
+
+def leaf_diagnostics(coarse, fine, masks):
+    """Fluid leaves only; a non-finite population or rho <= 0 there is an error.
 
     Mass counts a fine cell as 1/4 of a coarse cell (2D, 2:1)."""
     out = {}
-    for name, g, area in (("coarse", coarse, 1.0), ("fine", fine, .25)):
-        leaf = np.isfinite(g).all(axis=0)
+    for name, g, leaf, area in (("coarse", coarse, masks[0], 1.0), ("fine", fine, masks[1], .25)):
         f = g[:, leaf].astype(np.float64)
         rho = f.sum(axis=0)
+        if not np.isfinite(f).all() or (rho <= 0).any():
+            raise RuntimeError(f"non-finite population or non-positive density on a {name} fluid leaf")
         ux = (np.array(D2Q9_EX)[:, None] * f).sum(axis=0) / rho
         uy = (np.array(D2Q9_EY)[:, None] * f).sum(axis=0) / rho
         out[name] = (area * rho.sum(), rho.min(), rho.max(), np.hypot(ux, uy).max())
@@ -83,11 +98,12 @@ def main():
         storage = f"blocked {fine_block}" if fine_block else "dense"
         print(f"mode={args.mode} fine storage={storage} speed={args.speed:+g} "
               f"(fine allocation {p.fine_shape[0]}x{p.fine_shape[1]}, fine origin {p.fine_origin})")
-        start = leaf_diagnostics(*run.download())
+        masks = leaf_masks(run, args.mode)
+        start = leaf_diagnostics(*run.download(), masks)
         for _ in range(0, args.steps, 32):
             run.step(min(32, args.steps - run.step_count))
             coarse, fine = run.download()  # synchronized postcollision G, NaN on non-leaves
-            d = leaf_diagnostics(coarse, fine)
+            d = leaf_diagnostics(coarse, fine, masks)
             mass = d["coarse"][0] + d["fine"][0]
             rho_lo = min(d["coarse"][1], d["fine"][1])
             rho_hi = max(d["coarse"][2], d["fine"][2])