Gate shifted velocity Jets to all-blocked callers

src/CelerisLab/body/jets.py

@@ -189,13 +189,14 @@ def edge_to_node_frame(points):
 
 
 def surface_link_velocity(objects, cfg, flags, fluid_idx, dir, body_id, rx, ry,
-                          jets=(), hits=None):
+                          jets=(), hits=None, *, _shifted_blocked=False):
     """Per-link prescribed surface velocity (n, 2) FP32 and aperture mask.
 
     Host-only and storage-neutral: rows follow the given physical cut-link
     order. Returns (None, None) when neither an object carries surface
     actuation nor a jet is given. Jets need true FP64 ``hits`` (n, 2) and may
     not share a link with each other or with a legacy aperture.
+    Shifted storage is privately enabled only by the all-blocked stepper.
     """
     enabled = [obj for obj in objects if obj.surface_actuation is not None]
     if not enabled and not jets:
@@ -219,7 +220,9 @@ def surface_link_velocity(objects, cfg, flags, fluid_idx, dir, body_id, rx, ry,
         if hits is None or len(hits) != len(dir):
             raise ValueError("Jets require true hit points row-aligned with the cut links")
         if (cfg is None or not cfg.is_d2q9 or cfg.store_precision != "FP32"
-                or cfg.streaming != "double_buffer" or cfg.ddf_shifting):
+                or cfg.streaming != "double_buffer"
+                or (cfg.ddf_shifting and not (_shifted_blocked is True
+                    and cfg.collision == "SRT" and not cfg.les_enabled))):
             raise ValueError("Jets require D2Q9/FP32/unshifted/double_buffer")
         shape, jet_id = jet_link_shapes(jets, hits)
         driven = jet_id >= 0

src/CelerisLab/lbm/refinement.py

@@ -422,7 +422,7 @@ def _circle_in_coverage(covered, cx, cy, radius):
 
 
 def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=(),
-                          parent_origin=(0., 0.), parent_spacing=1.):
+                          parent_origin=(0., 0.), parent_spacing=1., _shifted_blocked=False):
     """Project fixed circles and reject boundary/transfer overlap before CUDA.
 
     Input centers/radii and Jets are in the root edge frame. ``parent_origin``
@@ -555,7 +555,7 @@ def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=(),
         jet_cfg.validate()
         fluid_idx, dirs, _, body_id, rx, ry, _, _ = packed
         surface_link_velocity(local, jet_cfg, flags, fluid_idx, dirs, body_id, rx, ry,
-                              jets=jets, hits=hits)
+                              jets=jets, hits=hits, _shifted_blocked=_shifted_blocked)
     return FineBoundary2D(tuple(local), _readonly(flags, np.uint16),
                           tuple(_readonly(a, a.dtype) for a in packed),
                           _readonly(sorted(boundary), np.uint32),

src/CelerisLab/lbm/refinement_stepper.py

@@ -287,11 +287,12 @@ class FixedPatchStepper:
             pair = pairs[level - 1]
             if level == 1:
                 host = prepare_fine_boundary(cfg, pair, level_objects[1],
-                                             wall_aware=self._wall_aware, jets=level_jets[1])
+                                             wall_aware=self._wall_aware, jets=level_jets[1],
+                                             _shifted_blocked=self._ddf_shifting)
             else:
                 host = prepare_fine_boundary(self._level_cfgs[level - 1], pair, level_objects[level],
                                              jets=level_jets[level], parent_origin=origin,
-                                             parent_spacing=spacing)
+                                             parent_spacing=spacing, _shifted_blocked=self._ddf_shifting)
             self._boundary_hosts.append(host)
             origin = tuple(o + spacing * f for o, f in zip(origin, pair.fine_origin))
             spacing /= 2
@@ -309,7 +310,8 @@ class FixedPatchStepper:
         self._level_jets = level_jets
         self._level_boundaries = [None] + [
             _prepare_blocked_boundary(self._boundary_hosts[level], self._layouts[level],
-                                      self._level_cfgs[level], level_jets[level])
+                                      self._level_cfgs[level], level_jets[level],
+                                      _shifted_blocked=self._ddf_shifting)
             for level in range(1, depth)]
         self._leaf = [kind == ACTIVE for kind in self.hierarchy.kinds]
         for level in range(1, depth):
@@ -558,6 +560,7 @@ class FixedPatchStepper:
 
         with self._lock:
             self._require("ready")
+            self._check_storage()
             pairs = self.hierarchy.pairs
             if (self._wall_aware or self.hierarchy.depth != 3
                     or not isinstance(pairs[1], FixedPatchSet2D)
@@ -569,19 +572,22 @@ class FixedPatchStepper:
             actual = self._level_boundaries[2]
             jets = () if actual is None else actual.jets
             old = RemapTopology.from_config(cfg, pairs[0], pairs[1].bounds,
-                                           child_objects=child_objects, child_jets=jets)
+                                           child_objects=child_objects, child_jets=jets,
+                                           _shifted_blocked=self._ddf_shifting)
             if (not same_physical_boundary(old.boundary, self._boundary_hosts[2])
                     or (self._boundary_hosts[2] is not None and actual is None)):
                 raise ValueError("child_objects must reproduce the original root-edge-frame boundary")
             if actual is not None:
                 velocity = actual.tables.get("surface_velocity")
                 if velocity is None:
-                    velocity = boundary_velocity(cfg, pairs[1], actual.host, jets)
+                    velocity = boundary_velocity(cfg, pairs[1], actual.host, jets,
+                                                 _shifted_blocked=self._ddf_shifting)
                 if (not same_physical_boundary(old.boundary, actual.host)
                         or not np.array_equal(old.surface_velocity, velocity)):
                     raise ValueError("child_objects changed boundary ownership or surface velocity")
             new = RemapTopology.from_config(cfg, pairs[0], rectangles,
-                                           child_objects=child_objects, child_jets=jets)
+                                           child_objects=child_objects, child_jets=jets,
+                                           _shifted_blocked=self._ddf_shifting)
             if (not same_physical_boundary(old.boundary, new.boundary)
                     or not np.array_equal(old.surface_velocity, new.surface_velocity)):
                 raise ValueError("remap requires the identical pinned physical boundary and BFL support")

src/CelerisLab/lbm/refinement_tables.py

@@ -87,7 +87,7 @@ class _BlockedBoundary:
     accumulate: object = None
 
 
-def _prepare_blocked_boundary(boundary, layout, cfg, jets):
+def _prepare_blocked_boundary(boundary, layout, cfg, jets, *, _shifted_blocked=False):
     if boundary is None or layout is None:
         return None
     fluid_idx, dirs, q, body_id, rx, ry, _, fallback = boundary.curved_host
@@ -102,6 +102,6 @@ def _prepare_blocked_boundary(boundary, layout, cfg, jets):
         cfg.validate()
         tables["surface_velocity"], tables["aperture"] = surface_link_velocity(
             boundary.objects, cfg, boundary.flags, fluid_idx, dirs, body_id, rx, ry,
-            jets=jets, hits=boundary.curved_hits)
+            jets=jets, hits=boundary.curved_hits, _shifted_blocked=_shifted_blocked)
     return _BlockedBoundary(boundary, tables, jets, actuated,
                             jet_link_shapes(jets, boundary.curved_hits) if jets else None)

src/CelerisLab/lbm/synchronized_remap.py

@@ -39,7 +39,7 @@ class RemapTopology:
 
     @classmethod
     def from_config(cls, cfg, parent, child_rectangles, *, coarse_boundary="periodic", objects=(),
-                    child_objects=(), child_jets=()):
+                    child_objects=(), child_jets=(), _shifted_blocked=False):
         if coarse_boundary != "periodic" or tuple(objects) or cfg.nz != 1:
             raise ValueError("remap requires periodic 2D state; use child_objects for a pinned L2 body")
         if type(parent) not in (FixedPatch2D, FixedPatchSet2D):
@@ -51,7 +51,7 @@ class RemapTopology:
                            viscosity=2 * cfg.viscosity)
         boundary = prepare_fine_boundary(
             fine_cfg, child, child_objects, jets=child_jets,
-            parent_origin=parent.fine_origin, parent_spacing=.5)
+            parent_origin=parent.fine_origin, parent_spacing=.5, _shifted_blocked=_shifted_blocked)
         pinned, velocity = (), None
         child_leaf = child.fine_kind == ACTIVE
         if boundary is not None:
@@ -67,7 +67,8 @@ class RemapTopology:
                 if match is not None:
                     containing.add(match)
             pinned = tuple(sorted(containing))
-            velocity = boundary_velocity(cfg, child, boundary, child_jets)
+            velocity = boundary_velocity(cfg, child, boundary, child_jets,
+                                         _shifted_blocked=_shifted_blocked)
             if velocity is not None:
                 velocity.setflags(write=False)
         masks = (~parent.covered, middle_kind == ACTIVE, child_leaf)
@@ -76,7 +77,7 @@ class RemapTopology:
         return cls(parent, child, masks, boundary, pinned, velocity)
 
 
-def boundary_velocity(cfg, child, boundary, jets=()):
+def boundary_velocity(cfg, child, boundary, jets=(), *, _shifted_blocked=False):
     """Physical cut-link order; never blocked addresses or stale Jet snapshots."""
     local_cfg = replace(cfg, nx=child.fine_shape[1], ny=child.fine_shape[0],
                         viscosity=4 * cfg.viscosity)
@@ -84,7 +85,7 @@ def boundary_velocity(cfg, child, boundary, jets=()):
     fluid, dirs, _, ids, rx, ry, _, _ = boundary.curved_host
     velocity, _ = surface_link_velocity(
         boundary.objects, local_cfg, boundary.flags, fluid, dirs, ids, rx, ry,
-        jets=jets, hits=boundary.curved_hits)
+        jets=jets, hits=boundary.curved_hits, _shifted_blocked=_shifted_blocked)
     return velocity
 
 

tests/unit/test_shifted_jets.py

@@ -0,0 +1,105 @@
+"""Shifted velocity Jets: bounded opt-in and host lifecycle, no CUDA numerics."""
+from dataclasses import replace
+import unittest
+from unittest.mock import patch
+import numpy as np
+from CelerisLab.body.jets import Jet, Tube, Parallel
+from CelerisLab.body.objects import SurfaceActuation
+from CelerisLab.lbm.synchronized_remap import RemapTopology
+from CelerisLab.lbm.refinement import FixedPatchSet2D
+from CelerisLab.lbm.population_storage import _StoredPopulations
+from tests.unit import test_refinement_three_level_multipatch as backend
+from tests.unit.test_pinned_body_remap import A, B, body, jet, mesh
+from tests.unit.test_shifted_storage import shifted_fields
+
+CFG=replace(backend.CFG,ddf_shifting=True)
+
+
+def topology(objects, jets, *, enabled=False, cfg=CFG):
+    return RemapTopology.from_config(cfg,FixedPatchSet2D.from_config(cfg,backend.PARENTS),A,
+        child_objects=objects,child_jets=jets,_shifted_blocked=enabled)
+
+
+def ready():
+    objects,jets=(body(),),(jet(),)
+    with patch.object(backend,'CFG',CFG),patch.object(backend,'bodies_and_jets',return_value=(objects,jets)):
+        run,module,allocations=backend.mock_backend(children=A,bodies=True)
+    storage={}
+    module.cuda.memcpy_htod.side_effect=lambda p,a:storage.__setitem__(id(p),np.array(a,copy=True))
+    module.cuda.memcpy_dtoh.side_effect=lambda a,p:np.copyto(a,storage[id(p)])
+    run._initialize_raw(_StoredPopulations(shifted_fields(mesh(A,objects,jets)),True))
+    run._step_count=7
+    return run,module,allocations,storage,objects
+
+
+class ShiftedJets(unittest.TestCase):
+    def test_default_guard_and_storage_neutral_host_rows(self):
+        with self.assertRaisesRegex(ValueError,'unshifted'):
+            topology((body(),),(jet(),))
+        shifted=topology((body(),),(jet(),),enabled=True)
+        unshifted=topology((body(),),(jet(),),cfg=backend.CFG)
+        np.testing.assert_array_equal(shifted.surface_velocity,unshifted.surface_velocity)
+        self.assertEqual(shifted.surface_velocity.dtype,np.float32)
+        legacy=body(SurfaceActuation(.7,.85,.002))
+        with self.assertRaisesRegex(ValueError,'unshifted'):
+            topology((legacy,),(),enabled=True)
+
+    def test_rejections_before_device_allocation(self):
+        module=backend.host_module()
+        moving=body();moving.state.omega=.01
+        cases=[((moving,),(jet(),)),((body(),),(jet(),jet())),
+               ((body(),),(Jet(Tube((100.,100.),0.,1.,.1),Parallel(0.)),)),
+               ((body(center=(8.,14.)),),(jet(),)),
+               ((body(),body(center=(34.,14.))),
+                (Jet(Tube((10.,14.),0.,26.,.3),Parallel(0.)),))]
+        for objects,jets in cases:
+            with self.subTest(objects=len(objects),jets=len(jets)),self.assertRaises(ValueError):
+                module.FixedPatchStepper(CFG,backend.PARENTS,context=None,
+                    base_block=(4,4),fine_block=(4,4),child_block=(4,4),
+                    child_rectangles_in_parent=A,child_objects=objects,child_jets=jets)
+        for blocks in ((None,None,None),((4,4),None,(4,4))):
+            with self.assertRaises(ValueError):
+                module.FixedPatchStepper(CFG,backend.PARENTS,context=None,
+                    base_block=blocks[0],fine_block=blocks[1],child_block=blocks[2],
+                    child_rectangles_in_parent=A,child_objects=(body(),),child_jets=(jet(),))
+        module.cuda.mem_alloc.assert_not_called()
+        module.compiler.compile_kernel.assert_not_called()
+
+    def test_switch_upload_only_then_raw_remap_keeps_live_velocity(self):
+        run,module,allocations,storage,objects=ready();clone=None
+        try:
+            boundary=run._level_boundaries[2]
+            before=run._download_raw()
+            original=boundary.tables['surface_velocity'].copy()
+            module.cuda.memcpy_htod.reset_mock()
+            run.set_jet_amplitude(0,-jet().amplitude,level=2)
+            module.cuda.memcpy_htod.assert_called_once()
+            self.assertIs(module.cuda.memcpy_htod.call_args.args[0],boundary.device['surface_velocity'])
+            np.testing.assert_array_equal(boundary.tables['surface_velocity'],-original)
+            self.assertEqual(run.step_count,7)
+            for a,b in zip(before.values,run._download_raw().values):self.assertEqual(a.tobytes(),b.tobytes())
+            clone=run.remap_child(B,child_objects=objects)
+            self.assertEqual(clone._level_jets[2][0].amplitude,-jet().amplitude)
+            np.testing.assert_array_equal(clone._level_boundaries[2].tables['surface_velocity'],
+                                          boundary.tables['surface_velocity'])
+            self.assertEqual(clone.step_count,7)
+            saved=run._download_raw()
+            with self.assertRaises(ValueError):run.remap_child(B,child_objects=(body(center=(12.1,14.)),))
+            for a,b in zip(saved.values,run._download_raw().values):self.assertEqual(a.tobytes(),b.tobytes())
+        finally:
+            if clone is not None:clone.close()
+            run.close()
+
+    def test_invalid_amplitude_does_not_upload_or_commit(self):
+        run,module,_,_,_=ready()
+        try:
+            before=run._level_boundaries[2].tables['surface_velocity'].copy()
+            module.cuda.memcpy_htod.reset_mock()
+            for amplitude in (float('nan'),float('inf'),1.):
+                with self.assertRaises(ValueError):run.set_jet_amplitude(0,amplitude,level=2)
+            module.cuda.memcpy_htod.assert_not_called()
+            np.testing.assert_array_equal(before,run._level_boundaries[2].tables['surface_velocity'])
+        finally:run.close()
+
+
+if __name__=='__main__':unittest.main()

tests/validation/multigrid/test_shifted_jets_gpu.py

@@ -0,0 +1,208 @@
+"""Shifted velocity-Jet GPU gates. Run only on an authorized device.
+
+PYTHONPATH=src python tests/validation/multigrid/test_shifted_jets_gpu.py -v
+Uniform initial curved-body states rho=.99/1.01; amplitude 0,+.002,-.001,0;
+four-level substep ledger; three-level live-Jet raw remap. No pressure actuator.
+The census is an implementation-exchange ledger, not an independent physical
+traction/nozzle-thrust model. Frozen inherited caps: raw rows rtol=2e-6/atol=2e-7; fluid ledger
+rtol=2e-4/atol=8e-5. Observer on/off, interval accumulation and raw copies exact.
+Report actual residuals too: passing old loose caps alone is not an accuracy gain.
+"""
+from contextlib import ExitStack
+import unittest
+import numpy as np
+import pycuda.driver as cuda
+from CelerisLab.config import LBMConfig
+from CelerisLab.cuda.context import CudaContext
+from CelerisLab.body.geometry.circle import CircleGeometry
+from CelerisLab.body.objects import SimObject
+from CelerisLab.body.jets import Jet, Tube, Parallel
+from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY, SOLID
+from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
+from CelerisLab.lbm.population_storage import _StoredPopulations
+from CelerisLab.lbm.synchronized_remap import RemapTopology, _remap_stored
+
+CX,CY=np.asarray(D2Q9_EX),np.asarray(D2Q9_EY)
+OPP=np.array([np.flatnonzero((CX==-x)&(CY==-y))[0] for x,y in zip(CX,CY)])
+W=np.array([4/9,*([1/9]*4),*([1/36]*4)],np.float64)
+SPEEDS=(0.,.002,-.001,0.)
+
+
+def circle(center,radius):
+    return SimObject(17,geometry=CircleGeometry(*center,radius),center=center,radius=radius)
+
+
+def drive(center,radius):
+    return Jet(Tube(center,0.,radius*1.5,radius*.32),Parallel(0.),amplitude=0.)
+
+
+def initial(run,rho,u=.013,v=-.007):
+    cu=3*(CX*u+CY*v)
+    g=W*((rho-1)+rho*(cu+.5*cu*cu-1.5*(u*u+v*v)))
+    return _StoredPopulations(tuple(np.broadcast_to(g[:,None,None],l.shape).astype(np.float32).copy()
+                                    for l in run._levels),True)
+
+
+def download(pointer,shape,dtype=np.float32):
+    out=np.empty(shape,dtype);cuda.memcpy_dtoh(out,pointer);return out
+
+
+def raw_totals(a,leaf):
+    g=a[:,leaf].astype(np.float64)
+    # Only differences at the same fixed leaf mask are used, so mass's unit
+    # background cancels. Momentum never requires adding weights.
+    return np.array([(g*CX[:,None]).sum(),(g*CY[:,None]).sum(),g.sum()])
+
+
+def census(run,index,buffer):
+    """Independent physical link/face census of actual staged receiver windows."""
+    level=run._levels[index];layout=level.layout;b=level.boundary
+    stored=download(buffer,layout.blocked_shape);dense=layout.to_dense(stored)
+    leaf=run._leaf[index];ny,nx=leaf.shape
+    solid=(b.host.flags.reshape(ny,nx)&SOLID)!=0
+    y,x=np.nonzero(leaf);owner=layout.owner[y,x].astype(np.intp)
+    cells=int(np.prod(layout.local_shape));block,cell=np.divmod(owner,cells)
+    flat=stored.reshape(layout.nblocks,9,cells);external=np.zeros(3,np.float64)
+    for q,(cx,cy) in enumerate(zip(CX,CY)):
+        sy,sx=y-cy,x-cx;dy,dx=y+cy,x+cx
+        incoming=flat[block,q,cell-cy*layout.local_shape[1]-cx].astype(np.float64)+W[q]
+        outgoing=flat[block,q,cell].astype(np.float64)+W[q]
+        entering=~leaf[sy,sx]&~solid[sy,sx];leaving=~leaf[dy,dx]&~solid[dy,dx]
+        external+=np.array([cx,cy,1])*(incoming[entering].sum()-outgoing[leaving].sum())
+    receiver,direction,_,_,rx,ry,_,_=b.host.curved_host
+    owner=layout.owner.ravel()[receiver].astype(np.intp);block,cell=np.divmod(owner,cells)
+    outgoing=flat[block,direction,cell].astype(np.float64)
+    source=cell+CY[direction]*layout.local_shape[1]+CX[direction]
+    reflected=flat[block,OPP[direction],source].astype(np.float64)
+    physical_sum=outgoing+reflected+2*W[direction]
+    jx,jy=CX[direction]*physical_sum,CY[direction]*physical_sum
+    rows=np.column_stack((jx,jy,rx*jy-ry*jx,reflected-outgoing))
+    return raw_totals(dense,leaf),external,rows
+
+
+class Ledger:
+    def __init__(self,test,run,index):
+        self.test,self.run,self.index=test,run,index
+        self.samples=[];self.max_row=0.;self.max_balance=0.
+        level=run._levels[index];self.original=level.kernel;b=level.boundary
+        def audited(src,dst,*args,**kwargs):
+            run._stream.synchronize()
+            before,external,independent=census(run,index,src)
+            raw=download(b.last_raw,(len(b.tables['cell']),4))
+            np.testing.assert_allclose(raw,independent,rtol=2e-6,atol=2e-7)
+            self.max_row=max(self.max_row,float(abs(raw.astype(float)-independent).max()))
+            self.original(src,dst,*args,**kwargs);run._stream.synchronize()
+            after=level.layout.to_dense(download(dst,level.layout.blocked_shape))
+            change=raw_totals(after,run._leaf[index])-before
+            expected=external+np.array([-independent[:,0].sum(),-independent[:,1].sum(),independent[:,3].sum()])
+            np.testing.assert_allclose(change,expected,rtol=2e-4,atol=8e-5)
+            self.max_balance=max(self.max_balance,float(abs(change-expected).max()))
+            self.samples.append(raw.copy())
+        level.kernel=audited
+
+    def step(self):
+        run=self.run;self.samples.clear();run.step()
+        b=run._level_boundaries[self.index];h=run._exchange_spacing[self.index]
+        self.test.assertEqual(len(self.samples),2**self.index)
+        expected=np.zeros_like(self.samples[0],dtype=np.float64)
+        weights=np.array([h*h,h*h,h*h*h,h*h])
+        for raw in self.samples:expected+=raw.astype(np.float64)*weights
+        actual=download(b.interval_integral,expected.shape,np.float64)
+        np.testing.assert_array_equal(actual,expected)
+        obs=run.read_raw_exchange(0,level=self.index)
+        for name in ('impulse','angular_impulse','mass_into_fluid'):
+            np.testing.assert_allclose(obs['total'][name],np.asarray(obs['wall'][name])+obs['aperture'][name],rtol=0,atol=2e-13)
+        self.test.assertEqual(obs['sample_count'],2**self.index)
+        for a,leaf in zip(run._download_raw().values,run._leaf):
+            self.test.assertTrue(np.isfinite(a[:,leaf]).all())
+            self.test.assertTrue((1+a[:,leaf].sum(axis=0,dtype=np.float64)>0).all())
+        return obs
+
+    def report(self,label):
+        print(f'{label}: max_abs_raw_row={self.max_row:.9g}, max_abs_fluid_ledger={self.max_balance:.9g}',flush=True)
+
+
+class ShiftedJetsGPU(unittest.TestCase):
+    def exact_raw(self,a,b):
+        self.assertTrue(a.shifted and b.shifted)
+        for x,y in zip(a.values,b.values):self.assertEqual(x.tobytes(),y.tobytes())
+
+    def test_matching_whole_circle_parallel_velocity_preserves_uniform_state(self):
+        cfg=LBMConfig(nx=32,ny=24,viscosity=.1,ddf_shifting=True)
+        center=(16.,12.);radius=2.15;obj=circle(center,radius)
+        # Rectangle contains every circle hit; prescribed surface velocity equals
+        # the uniform fluid velocity. Geometry itself remains fixed.
+        whole=Jet(Tube((center[0]-radius-1,center[1]),0.,2*radius+2,radius+1),
+                  Parallel(0.),amplitude=.013)
+        for rho in (.99,1.01):
+            with self.subTest(rho=rho),CudaContext(0) as ctx:
+                with FixedPatchStepper(cfg,(8,6,24,18),context=ctx,base_block=(4,4),fine_block=(4,4),
+                        fine_objects=(obj,),jets=(whole,),observe_raw_exchange=True) as run:
+                    run._initialize_raw(initial(run,rho,u=.013,v=0.))
+                    expected=run._download_raw();ledger=Ledger(self,run,1)
+                    self.assertTrue(run._level_boundaries[1].tables['aperture'].all())
+                    for _ in range(4):ledger.step()
+                    error=max(float(abs(a[:,leaf]-b[:,leaf]).max()) for a,b,leaf in
+                              zip(run._download_raw().values,expected.values,run._leaf))
+                    self.assertLessEqual(error,1e-6)  # inherited uniform-state cap
+                    print(f'matched whole-circle rho={rho}: max_abs_raw_state={error:.9g}',flush=True)
+                    ledger.report(f'matched whole-circle rho={rho}')
+
+    def test_uniform_initial_curved_rho_both_signs_switch_and_observer(self):
+        cfg=LBMConfig(nx=32,ny=24,viscosity=.1,ddf_shifting=True)
+        obj=circle((16.,12.),2.15);j=drive((16.,12.),2.15)
+        for rho in (.99,1.01):
+            with self.subTest(rho=rho),CudaContext(0) as ctx,ExitStack() as owned:
+                settings=dict(bounds=(8,6,24,18),context=ctx,base_block=(4,4),fine_block=(4,4),fine_objects=(obj,),jets=(j,))
+                measured=owned.enter_context(FixedPatchStepper(cfg,observe_raw_exchange=True,**settings))
+                plain=owned.enter_context(FixedPatchStepper(cfg,**settings))
+                measured._initialize_raw(initial(measured,rho));plain._initialize_raw(initial(plain,rho))
+                ledger=Ledger(self,measured,1)
+                for speed in SPEEDS:
+                    before=measured._download_raw();count=measured.step_count
+                    measured.set_jet_amplitude(0,speed);plain.set_jet_amplitude(0,speed)
+                    self.exact_raw(before,measured._download_raw());self.assertEqual(measured.step_count,count)
+                    ledger.step();plain.step()
+                    self.exact_raw(measured._download_raw(),plain._download_raw())
+                ledger.report(f'L1 rho={rho}')
+
+    def test_four_level_interval_ledger(self):
+        cfg=LBMConfig(nx=48,ny=48,viscosity=.02,ddf_shifting=True)
+        obj=circle((16.,16.),.5);j=drive((16.,16.),.5)
+        with CudaContext(0) as ctx,ExitStack() as owned:
+            settings=dict(context=ctx,refinements=[[(6,6,26,26)],[(20,20,44,44)],[(56,56,72,72)]],
+                          level_blocks=[(4,4)]*4,level_objects=[(),(),(),(obj,)],level_jets=[(),(),(),(j,)])
+            measured=owned.enter_context(FixedPatchStepper(cfg,observe_raw_exchange=True,**settings))
+            plain=owned.enter_context(FixedPatchStepper(cfg,**settings))
+            measured._initialize_raw(initial(measured,1.01));plain._initialize_raw(initial(plain,1.01))
+            ledger=Ledger(self,measured,3)
+            for speed in SPEEDS:
+                measured.set_jet_amplitude(0,speed,level=3);plain.set_jet_amplitude(0,speed,level=3)
+                ledger.step();plain.step();self.exact_raw(measured._download_raw(),plain._download_raw())
+            ledger.report('L3 rho=1.01')
+
+    def test_live_jet_remap_then_both_ledgers_continue(self):
+        cfg=LBMConfig(nx=48,ny=32,viscosity=.1,ddf_shifting=True)
+        parents=((4,4,20,28),(24,4,44,28));a=((16,16,32,40),(56,16,72,40));b=(a[0],(60,16,76,40))
+        obj=circle((12.,14.),.8);j=drive((12.,14.),.8)
+        with CudaContext(0) as ctx,ExitStack() as owned:
+            old=owned.enter_context(FixedPatchStepper(cfg,parents,context=ctx,base_block=(4,4),fine_block=(4,4),
+                child_block=(4,4),child_rectangles_in_parent=a,child_objects=(obj,),child_jets=(j,),observe_raw_exchange=True))
+            old._initialize_raw(initial(old,.99));ledger=Ledger(self,old,2)
+            for speed in (.002,-.001):old.set_jet_amplitude(0,speed,level=2);ledger.step()
+            before=old._download_raw();count=old.step_count
+            topology=lambda rects:RemapTopology.from_config(cfg,old.patch,rects,child_objects=(obj,),
+                child_jets=old._level_jets[2],_shifted_blocked=True)
+            expected=_remap_stored(topology(a),topology(b),before,shifted=True,base_count=count,
+                level_steps=(count,2*count,4*count),phase='postcollision_g')
+            new=owned.enter_context(old.remap_child(b,child_objects=(obj,)))
+            self.exact_raw(before,old._download_raw())
+            self.exact_raw(_StoredPopulations(expected.g,True),new._download_raw())
+            self.assertEqual(new._level_jets[2][0].amplitude,-.001)
+            np.testing.assert_array_equal(new._level_boundaries[2].tables['surface_velocity'],old._level_boundaries[2].tables['surface_velocity'])
+            newledger=Ledger(self,new,2);ledger.step();newledger.step()
+            self.assertEqual((old.step_count,new.step_count),(count+1,count+1))
+            ledger.report('L2 source');newledger.report('L2 remapped')
+
+
+if __name__=='__main__':unittest.main()