Example/docs closure fixes

- README: close the Manual stream control code fence before the new
  section; snippet now builds rest-equilibrium inputs (runnable as shown).
- Example: diagnostics use the known fluid-leaf masks (uncovered coarse
  minus slot-mode wall rows; fine ACTIVE minus solids) and raise on any
  non-finite population or rho <= 0 there instead of filtering by
  finiteness. No solver changes.

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

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