diff --git a/docs/synchronized_leaf_remap.md b/docs/synchronized_leaf_remap.md index f4f4523..639cf53 100644 --- a/docs/synchronized_leaf_remap.md +++ b/docs/synchronized_leaf_remap.md @@ -113,22 +113,24 @@ Per-population conservation is not an advancement target under collision. ### Bounded pinned-body extension `replacement = current.remap_child(rectangles, child_objects=original_objects)` -now also accepts exactly one stationary Circle in one of exactly two L2 -rectangles. Its **complete containing rectangle must remain identical**; -only the separated fluid rectangle may change. Existing nesting, four-cell -gap, full-circle containment, reverse cut-link coverage and two-hop BFL +also accepts stationary Circles in separated L2 rectangles. Every body’s +**complete containing rectangle must remain identical**. Several bodies may +share one pinned rectangle; there need not be any fluid-only rectangles. +Separated fluid-only rectangles may be added, removed or moved. Existing +nesting, four-cell gap, full-circle containment, reverse cut-link coverage and two-hop BFL transfer-footprint checks still apply. This is manual remapping, not adaptive body tracking or runtime body editing. -Supply the original root edge-frame definition used at construction. Do not -pass `current.child_boundary.objects`: those are projected fine-local copies. +Supply all original root edge-frame definitions in their construction order. +Do not pass `current.child_boundary.objects`: those are projected fine-local copies. There is no retained mutable original-object registry. Before reading device state, remapping reprojects the supplied definition onto the old topology and -compares flags, physical cut links, body IDs, q, lever arms, hits, BFL support, +compares exact ordered circle centers/radii, flags, physical cut links, body +IDs, q, lever arms, hits, BFL support, legacy SurfaceActuation specification and current per-link velocity against the existing boundary. Omission, changed geometry/actuation and changed body -rectangle fail without device reads or allocations. Root object IDs are still -assigned by list position, as in the constructor. +rectangle fail without device reads or allocations; reordered bodies fail +too. Root object IDs are still assigned by list position, as in the constructor. Jets are taken from the live L2 boundary, including the most recent successful `set_jet_amplitude(..., level=2)`. There is no replacement-Jets argument or @@ -139,12 +141,20 @@ boundary and velocities must match exactly. typed `FineBoundary2D`; it derives L2 fluid leaves as ACTIVE minus SOLID. No caller mask is accepted. Solids, ghosts and holes are never donors and remain NaN, physical-fluid overlap is exact, and exposed L1 still restricts all four -old fluid children. The helper checks the old/new body pin and boundary too. +old fluid children. The helper checks the old/new body pins and boundary too. +Its immutable `pinned_rectangles` is a canonical, deduplicated tuple, empty for body-free layouts; rectangle input +order does not affect it. Scalar-parent window offsets remain part of the +coordinate conversion. Focused CPU mocks: `tests/unit/test_pinned_body_remap.py`. GPU declarations: `tests/validation/multigrid/test_pinned_body_remap_gpu.py`, covering no-op next-step bitwise identity, moved-fluid A→B→A, current-amplitude preservation, physical boundary/velocity equality, source independence and absolute counts. -These GPU cases have **not been run here**. They do not establish evolution +Multiple-body CPU coverage: `tests/unit/test_multiple_pinned_remap.py`. +`tests/validation/multigrid/test_multiple_pinned_remap_gpu.py` adds two pins, +two independently updated Jets, variable fluid-region counts and add/move/remove +immediate comparisons with CPU transfer. All nine solid-excluded area-weighted +population integrals, exact overlap and nonleaf NaNs are checked. +These new GPU cases have **not been run here**. They do not establish evolution accuracy. A 3e-6 FP32-vs-FP64 evolution gate requires a separate independent BFL reference; the body-free replay gate above does not establish that claim. diff --git a/src/CelerisLab/lbm/refinement_stepper.py b/src/CelerisLab/lbm/refinement_stepper.py index 846b5be..ab4c8cb 100644 --- a/src/CelerisLab/lbm/refinement_stepper.py +++ b/src/CelerisLab/lbm/refinement_stepper.py @@ -499,9 +499,9 @@ class FixedPatchStepper: Periodic three-level set children only; L0/L1 and block sizes stay fixed. A body run requires its original root-edge-frame ``child_objects``: - one stationary circle pinned in one of exactly two L2 rectangles. - The complete body rectangle stays identical; only the other fluid - rectangle may move. Current L2 Jets (including updated amplitudes) + stationary circles pinned in unchanged complete L2 rectangles. + Separated fluid-only rectangles may be added, moved or removed. + Current L2 Jets (including updated amplitudes) are preserved. Local projected boundary objects are not valid inputs. The lock and stream synchronization certify current postcollision G and the completed 1:2:4 counts internally. @@ -538,10 +538,10 @@ class FixedPatchStepper: raise ValueError("child_objects changed boundary ownership or surface velocity") new = RemapTopology.from_config(cfg, self.patch, rectangles, child_objects=child_objects, child_jets=jets) - if (old.pinned_rectangle != new.pinned_rectangle + if (old.pinned_rectangles != new.pinned_rectangles or 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 body rectangle and boundary") + raise ValueError("remap requires the identical pinned body rectangles and boundary") base_block, fine_block, child_block = self._remap_blocks # Check sparse block alignment/launch size before any device read/allocation. layout = new.child.fine_layout(child_block) diff --git a/src/CelerisLab/lbm/synchronized_remap.py b/src/CelerisLab/lbm/synchronized_remap.py index 204a697..3e10973 100644 --- a/src/CelerisLab/lbm/synchronized_remap.py +++ b/src/CelerisLab/lbm/synchronized_remap.py @@ -1,7 +1,7 @@ """CPU-only whole-volume G transfer, not a runtime adaptive solver. Periodic three-level grids with fixed L0/L1 topology are supported, including -one stationary circle pinned inside an unchanged L2 rectangle. The caller +stationary circles pinned inside unchanged L2 rectangles. The caller certifies that current owned leaves are quiescent, unshifted D2Q9 postcollision G at the supplied completed base interval. No device buffers, device boundary state, ghosts, collision or time advance @@ -25,12 +25,14 @@ class RemapTopology: Child rectangles are in the parent's allocated L1 coordinates, including its window offset/halo when the parent is a scalar patch. Dense masks are rebuilt from geometry, never inferred from finite values in input arrays. + ``pinned_rectangles`` canonically lists each body-containing rectangle once; + fluid-only layouts have no pins. """ parent: object child: FixedPatchSet2D leaf_masks: tuple boundary: FineBoundary2D | None = None - pinned_rectangle: tuple | None = None + pinned_rectangles: tuple = () surface_velocity: np.ndarray | None = None @classmethod @@ -43,24 +45,24 @@ class RemapTopology: parent = type(parent).from_config(cfg, parent.bounds) child, middle_kind = nest_patch_set(cfg, parent, child_rectangles) child_objects, child_jets = tuple(child_objects), tuple(child_jets) - if child_objects and (len(child_objects) != 1 or len(child.bounds) != 2): - raise ValueError("body remap requires one stationary circle and two L2 rectangles") fine_cfg = replace(cfg, nx=parent.fine_shape[1], ny=parent.fine_shape[0], viscosity=2 * cfg.viscosity) boundary = prepare_fine_boundary( fine_cfg, child, child_objects, jets=child_jets, parent_origin=parent.fine_origin, parent_spacing=.5) - pinned, velocity = None, None + pinned, velocity = (), None child_leaf = child.fine_kind == ACTIVE if boundary is not None: child_leaf &= (boundary.flags.reshape(child.fine_shape) & SOLID) == 0 # Whole circle containment, not merely the resolved solid cells. - obj = child_objects[0] - cx, cy = (np.asarray(obj.geometry.center) - parent.fine_origin) * 2 - radius = obj.geometry._radius * 2 - pinned = next(r for r in child.bounds if - r[0] <= cx-radius and cx+radius <= r[2] and - r[1] <= cy-radius and cy+radius <= r[3]) + containing = set() + for obj in child_objects: + cx, cy = (np.asarray(obj.geometry.center) - parent.fine_origin) * 2 + radius = obj.geometry._radius * 2 + containing.add(next(r for r in child.bounds if + r[0] <= cx-radius and cx+radius <= r[2] and + r[1] <= cy-radius and cy+radius <= r[3])) + pinned = tuple(sorted(containing)) velocity = boundary_velocity(cfg, child, boundary, child_jets) if velocity is not None: velocity.setflags(write=False) @@ -88,8 +90,11 @@ def same_physical_boundary(a, b): return a is b if not isinstance(a, FineBoundary2D) or not isinstance(b, FineBoundary2D): return False - if (tuple(o.surface_actuation for o in a.objects) - != tuple(o.surface_actuation for o in b.objects)): + # Exact ordered geometry also catches changes hidden by cut-table rounding. + def signatures(boundary): + return tuple((o.geometry.center, o.geometry._radius, o.surface_actuation) + for o in boundary.objects) + if signatures(a) != signatures(b): return False return all(np.array_equal(x, y) for x, y in zip( (a.flags, a.boundary_cells, a.curved_hits, *a.curved_host), @@ -162,10 +167,10 @@ def remap_synchronized(old, new, old_g, *, base_count, level_steps, phase, or not _same_parent(old.parent, new.parent) or old.child.tau_fine != new.child.tau_fine): raise ValueError("remap requires unchanged L0/L1 geometry and relaxation") - if (old.pinned_rectangle != new.pinned_rectangle + if (old.pinned_rectangles != new.pinned_rectangles or 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 body rectangle and physical boundary") + raise ValueError("remap requires the identical pinned body rectangles and physical boundary") if not isinstance(level_steps, (tuple, list)): raise ValueError("level_steps must be three completed integer counts") counts = (base_count, *level_steps) diff --git a/tests/unit/test_multiple_pinned_remap.py b/tests/unit/test_multiple_pinned_remap.py new file mode 100644 index 0000000..0ee275f --- /dev/null +++ b/tests/unit/test_multiple_pinned_remap.py @@ -0,0 +1,153 @@ +"""Multiple body pins and variable fluid sets; mocked runtime is host-only.""" +import unittest +from unittest.mock import patch + +import numpy as np + +from CelerisLab.config import LBMConfig +from CelerisLab.body.jets import Jet, Tube, Parallel +from CelerisLab.body.objects import SurfaceActuation +from CelerisLab.lbm.refinement import FixedPatch2D +from CelerisLab.lbm.synchronized_remap import RemapTopology +from tests.unit.test_pinned_body_remap import body, mesh +from tests.unit.test_refinement_three_level_multipatch import mock_backend +from tests.unit.test_synchronized_remap import fields, integrals, remap + +PINS = ((16, 16, 32, 24), (56, 16, 72, 24)) +FLUID = ((16, 32, 32, 40), (56, 32, 72, 40)) +ADDED = ((16, 44, 32, 48),) +MOVED = ((20, 32, 36, 40), (60, 32, 76, 40)) + + +def objects(): + return (body(center=(12., 10.)), body(center=(32., 10.))) + + +def ready(): + bodies = objects() + jets = tuple(Jet(Tube((x-1.5, 10.), 0., 3., .3), Parallel(0.), + 'parabolic', .003) for x in (12., 32.)) + with patch('tests.unit.test_refinement_three_level_multipatch.bodies_and_jets', + return_value=(bodies, jets)): + run, module, allocations = mock_backend(children=PINS+FLUID, 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(*fields(mesh(PINS+FLUID, bodies, jets), np.float32)) + run._step_count = 7 + return run, module, allocations, bodies + + +class MultiplePinnedRemap(unittest.TestCase): + def check_transfer(self, old, new, arrays, tolerance): + out = remap(old, new, arrays) + before = integrals(old, arrays) + np.testing.assert_allclose(integrals(new, out.g), before, rtol=0, + atol=tolerance*np.abs(before[:9]).sum()) + for a, b, om, nm in zip(arrays, out.g, old.leaf_masks, new.leaf_masks): + self.assertEqual(a[:, om & nm].tobytes(), b[:, om & nm].tobytes()) + self.assertTrue(np.isnan(b[:, ~nm]).all()) + return out.g + + def test_add_move_remove_and_body_only_all_populations(self): + for dtype, tolerance in ((np.float64, 2e-13), (np.float32, 2e-7)): + old = mesh(PINS+FLUID, objects()) + arrays = fields(old, dtype) + for rects in (PINS+FLUID+ADDED, PINS+MOVED+ADDED, PINS, PINS+FLUID): + new = mesh(rects, objects()) + self.assertEqual(new.pinned_rectangles, PINS) + arrays = self.check_transfer(old, new, arrays, tolerance) + old = new + + def test_shared_pin_order_invariance_and_body_free(self): + bodies = (body(center=(11., 14.)), body(center=(14., 14.))) + rects = ((16, 16, 32, 40), (56, 16, 72, 40)) + old, new = mesh(rects, bodies), mesh(rects[::-1], bodies) + self.assertEqual(old.pinned_rectangles, (rects[0],)) + self.assertEqual(new.pinned_rectangles, old.pinned_rectangles) + self.check_transfer(old, new, fields(old), 2e-13) + only = mesh((rects[0],), bodies) + self.check_transfer(old, only, fields(old), 2e-13) + self.assertEqual(mesh(rects, ()).pinned_rectangles, ()) + + def test_shifted_scalar_parent_coordinates(self): + cfg = LBMConfig(nx=48, ny=32, viscosity=.1) + parent = FixedPatch2D.from_config(cfg, (4, 4, 44, 28)) + self.assertEqual(parent.fine_origin, (2., 2.)) + bodies = (body(center=(14., 12.)), body(center=(34., 12.))) + def topology(rects): + return RemapTopology.from_config(cfg, parent, rects, child_objects=bodies) + old, new = topology(PINS+FLUID), topology(PINS+MOVED) + self.assertEqual(old.pinned_rectangles, PINS) + self.check_transfer(old, new, fields(old), 2e-13) + + def test_direct_helper_rejects_each_changed_pin_and_body_order(self): + old = mesh(PINS+FLUID, objects()) + candidates = [mesh(PINS+FLUID, objects()[::-1])] + for index in range(2): + changed = list(PINS+FLUID) + x0, y0, x1, y1 = changed[index] + changed[index] = (x0-4, y0, x1, y1) + candidates.append(mesh(changed, objects())) + for new in candidates: + # Invalid donors ensure boundary validation precedes array reads. + with self.assertRaisesRegex(ValueError, 'pinned'): + remap(old, new, None) + + def test_live_independent_jets_and_source_ownership(self): + run, module, allocations, bodies = ready() + run.set_jet_amplitude(0, .017, level=2) + run.set_jet_amplitude(1, -.011, level=2) + before = run.download() + clone = run.remap_child((PINS+MOVED+ADDED)[::-1], child_objects=bodies) + try: + self.assertEqual(tuple(j.amplitude for j in clone._boundaries[2].jets), (.017, -.011)) + np.testing.assert_array_equal(clone._boundaries[2].tables['surface_velocity'], + run._boundaries[2].tables['surface_velocity']) + clone.set_jet_amplitude(1, .005, level=2) + self.assertEqual(run._boundaries[2].jets[1].amplitude, -.011) + clone.step() + self.assertEqual((run.step_count, clone.step_count), (7, 8)) + for got, want in zip(run.download(), before): + np.testing.assert_array_equal(got, want) + finally: + clone.close() + run.close() + + def test_sub_ulp_cut_table_geometry_change_rejected(self): + run, module, allocations, bodies = ready() + changed = objects() + changed[1].geometry._radius = np.nextafter(.8, np.inf) + start = len(allocations) + module.cuda.memcpy_dtoh.reset_mock() + with self.assertRaisesRegex(ValueError, 'boundary'): + run.remap_child(PINS+FLUID, child_objects=changed) + self.assertEqual(len(allocations), start) + module.cuda.memcpy_dtoh.assert_not_called() + with self.assertRaisesRegex(ValueError, 'pinned'): + remap(mesh(PINS+FLUID, bodies), mesh(PINS+FLUID, changed), None) + run.close() + + def test_each_pin_and_body_identity_rejected_before_device_access(self): + run, module, allocations, bodies = ready() + cases = [(PINS+FLUID, bodies[:1]), (PINS+FLUID, bodies[::-1]), + (PINS+FLUID, (bodies[0], body(center=(32.1, 10.)))), + (PINS+FLUID, (bodies[0], body(SurfaceActuation(0., .7, 0.), (32., 10.)))), + (PINS+FLUID, run.child_boundary.objects)] + for index in range(2): + changed = list(PINS+FLUID) + x0, y0, x1, y1 = changed[index] + changed[index] = (x0-4, y0, x1, y1) + cases.append((changed, bodies)) + cases.append((tuple(r for i, r in enumerate(PINS+FLUID) if i != index), bodies)) + start = len(allocations) + module.cuda.memcpy_dtoh.reset_mock() + run._stream.synchronize.reset_mock() + for rects, supplied in cases: + with self.subTest(rects=rects), self.assertRaises(ValueError): + run.remap_child(rects, child_objects=supplied) + self.assertEqual(len(allocations), start) + module.cuda.memcpy_dtoh.assert_not_called() + run._stream.synchronize.assert_not_called() + self.assertEqual((run._status, run.step_count), ('ready', 7)) + run.close() diff --git a/tests/unit/test_pinned_body_remap.py b/tests/unit/test_pinned_body_remap.py index 9d8df56..c105b75 100644 --- a/tests/unit/test_pinned_body_remap.py +++ b/tests/unit/test_pinned_body_remap.py @@ -55,7 +55,7 @@ class PinnedBodyRemap(unittest.TestCase): self.assertTrue(same_physical_boundary(old.boundary,new.boundary)) back=remap(new,old,out.g) self.assertTrue(np.isnan(back.g[2][:,~old.leaf_masks[2]]).all()) - self.assertEqual(old.pinned_rectangle,A[0]) + self.assertEqual(old.pinned_rectangles,(A[0],)) def test_current_jets_and_source_independence(self): run,module,allocations,storage,objects=ready() diff --git a/tests/validation/multigrid/test_multiple_pinned_remap_gpu.py b/tests/validation/multigrid/test_multiple_pinned_remap_gpu.py new file mode 100644 index 0000000..f518dd1 --- /dev/null +++ b/tests/validation/multigrid/test_multiple_pinned_remap_gpu.py @@ -0,0 +1,175 @@ +"""Manual GPU gates for two pinned bodies and several fluid-only L2 regions. + +Run explicitly on an approved idle GPU. The CPU helper checks immediate +whole-population transfer at the same physical time, not independent evolution +accuracy. Continuation checks cover finite leaves and no-op bitwise identity; +they do not establish an independent BFL replay or a numerical error bound. +""" +from contextlib import ExitStack +import unittest + +import numpy as np + +from CelerisLab.body.geometry.circle import CircleGeometry +from CelerisLab.body.jets import Jet, Parallel, Tube +from CelerisLab.body.objects import SimObject +from CelerisLab.cuda.context import CudaContext +from CelerisLab.lbm.refinement_stepper import FixedPatchStepper +from CelerisLab.lbm.synchronized_remap import ( + RemapTopology, remap_synchronized, same_physical_boundary, +) +from test_refinement_blocked_base import exact +from test_refinement_three_level_multipatch_gpu import CFG, PARENTS +from test_runtime_child_remap_gpu import RuntimeChildRemapGPU, initial + + +# All bounds use the allocated L1 frame, with support aligned to 4x4 blocks. +PINS = ((16, 16, 32, 24), (56, 16, 72, 24)) +FLUID = ((16, 32, 32, 40), (56, 32, 72, 40)) +EXTRA = ((16, 44, 32, 48),) +MOVED = ((20, 32, 36, 40), (60, 32, 76, 40)) +START = PINS + FLUID +CHANGES = ( + ('add', START + EXTRA, (.017, -.011)), + ('move', PINS + MOVED + EXTRA, (-.009, .005)), + ('remove', PINS + MOVED, (.013, -.007)), +) + + +def make(context): + centers = ((12., 10.), (32., 10.)) + objects = tuple( + SimObject(71 + index, geometry=CircleGeometry(*center, .8), + center=center, radius=.8) + for index, center in enumerate(centers) + ) + jets = tuple( + Jet(Tube((x - 1.5, y), 0., 3., .3), Parallel(0.), 'parabolic', amplitude) + for (x, y), amplitude in zip(centers, (.003, -.002)) + ) + run = FixedPatchStepper( + CFG, PARENTS, context=context, base_block=(4, 4), fine_block=(4, 4), + child_block=(4, 4), child_rectangles_in_parent=START, + child_objects=objects, child_jets=jets, + ) + return run, objects + + +def topology(run, rectangles, objects): + return RemapTopology.from_config( + CFG, run.patch, rectangles, child_objects=objects, + child_jets=run._boundaries[2].jets, + ) + + +def population_integrals(mesh, arrays): + """Nine separate area integrals in the root frame, accumulated in FP64.""" + return sum(array[:, leaf].sum(axis=1, dtype=np.float64) / 4**level + for level, (array, leaf) in enumerate(zip(arrays, mesh.leaf_masks))) + + +class MultiplePinnedRemapGPU(unittest.TestCase): + check_public = RuntimeChildRemapGPU.check_public + check_reset = RuntimeChildRemapGPU.check_reset + + def change_amplitudes(self, run, amplitudes): + # Each jet selects a different body. Updating one must preserve the + # other body's current per-link velocities, including its driven links. + boundary = run._boundaries[2] + body_ids = run.child_boundary.curved_host[3] + for index, amplitude in enumerate(amplitudes): + before = boundary.tables['surface_velocity'].copy() + run.set_jet_amplitude(index, amplitude, level=2) + after = boundary.tables['surface_velocity'] + exact(self, after[body_ids != index], before[body_ids != index], + 'other body velocity unchanged') + self.assertTrue(np.any(after[body_ids == index] != before[body_ids == index])) + self.assertEqual(tuple(jet.amplitude for jet in boundary.jets), amplitudes) + + def check_boundary(self, old, new, mesh, amplitudes): + self.assertEqual(mesh.pinned_rectangles, PINS) + self.assertTrue(same_physical_boundary(old.child_boundary, new.child_boundary)) + self.assertTrue(same_physical_boundary(new.child_boundary, mesh.boundary)) + # Compare physical cut-link rows, never layout-dependent blocked addresses. + for run in (old, new): + boundary = run._boundaries[2] + self.assertEqual(tuple(jet.amplitude for jet in boundary.jets), amplitudes) + exact(self, boundary.tables['surface_velocity'], mesh.surface_velocity, + 'latest physical per-link velocities') + + def test_noop_next_step_bitwise(self): + with CudaContext(0) as context, ExitStack() as owned: + old, objects = make(context) + owned.enter_context(old) + old.initialize(*initial(old, False)) + old.step(2) + amplitudes = (.017, -.011) + self.change_amplitudes(old, amplitudes) + before, count = old.download(), old.step_count + mesh = topology(old, START, objects) + # Input order is irrelevant even with several pins and fluid regions. + new = owned.enter_context(old.remap_child(START[::-1], child_objects=objects)) + self.assertIsNot(old, new) + self.assertEqual((old.step_count, new.step_count), (count, count)) + self.check_boundary(old, new, mesh, amplitudes) + self.check_reset(new, before) + for got, want in zip(old.download(), before): + exact(self, got, want, 'no-op source unchanged') + old.step() + new.step() + self.assertEqual((old.step_count, new.step_count), (count + 1, count + 1)) + for got, want in zip(new.download(), old.download()): + exact(self, got, want, 'no-op next step bitwise') + self.check_public(old) + self.check_public(new) + + def test_add_move_remove_immediate_transfer_and_source_independence(self): + with CudaContext(0) as context, ExitStack() as owned: + current, objects = make(context) + owned.enter_context(current) + current.initialize(*initial(current, False)) + current.step(2) + for action, rectangles, amplitudes in CHANGES: + with self.subTest(action=action): + previous = current + self.change_amplitudes(previous, amplitudes) + before, count = previous.download(), previous.step_count + old_mesh = topology(previous, previous.child_patch.bounds, objects) + new_mesh = topology(previous, rectangles, objects) + expected = remap_synchronized( + old_mesh, new_mesh, before, base_count=count, + level_steps=(count, 2 * count, 4 * count), + phase='postcollision_g', + ).g + current = owned.enter_context( + previous.remap_child(rectangles, child_objects=objects)) + self.assertIsNot(previous, current) + self.assertEqual((previous.step_count, current.step_count), (count, count)) + self.check_boundary(previous, current, new_mesh, amplitudes) + self.check_reset(current, expected) + actual = current.download() + for got, want, source, old_leaf, new_leaf in zip( + actual, expected, before, old_mesh.leaf_masks, new_mesh.leaf_masks): + exact(self, got, want, 'immediate CPU transfer parity') + exact(self, got[:, old_leaf & new_leaf], source[:, old_leaf & new_leaf], + 'same-level physical fluid overlap') + self.assertTrue(np.isfinite(got[:, new_leaf]).all()) + self.assertTrue(np.isnan(got[:, ~new_leaf]).all()) + # Conservation is asserted only across remap, not collision. + np.testing.assert_allclose( + population_integrals(new_mesh, actual), + population_integrals(old_mesh, before), + rtol=2 * np.finfo(np.float32).eps, atol=0, + err_msg='all nine area-weighted population integrals', + ) + self.check_public(current) + current.step() + self.assertEqual((previous.step_count, current.step_count), (count, count + 1)) + for got, want in zip(previous.download(), before): + exact(self, got, want, 'source unchanged after candidate advances') + self.check_boundary(previous, current, new_mesh, amplitudes) + self.check_public(current) + + +if __name__ == '__main__': + unittest.main()