diff --git a/docs/fixed_three_level_multipatch.md b/docs/fixed_three_level_multipatch.md new file mode 100644 index 0000000..b8190dd --- /dev/null +++ b/docs/fixed_three_level_multipatch.md @@ -0,0 +1,76 @@ +# Fixed three-level separated patches + +`FixedPatchStepper` supports exactly three static, blocked levels with separated +rectangle sets on either adjacent pair. This is not general AMR: there is no +per-region timestep, moving geometry, dynamic refinement or L1 body support. +All regions share one current/spare buffer pair per level. The existing schedule +advances L0 once, L1 twice and L2 four times per base step, with aggregate transfers. +No CUDA kernel changes are involved. + +## Explicit coordinate frames + +Keep the existing scalar `child=(x0,y0,x1,y1)` for a single scalar parent patch. +Its coordinates retain their original meaning: half-open integer indices in the +allocated L1 window, including its ghost band. + +The new keyword `child_rectangles_in_parent=[...]` is mutually exclusive with +`child`. It accepts a scalar or a set parent and interprets every rectangle in +the actual allocated L1 array, including ghost bands and holes. It never infers +region-local coordinates. The physical root cell edge corresponding to an L1 +index `(x,y)` is `parent.fine_origin + (x,y)/2`: + +- A scalar L1 parent has a shifted window origin. +- A set L1 parent has origin `(0,0)` and public shape `(2*ny,2*nx)`. +- A child set has public shape twice the allocated L1 shape. Its root origin is + the L1 allocated origin, even when the parent is a scalar window. + +For a 16x16 root, scalar parent `(4,4,12,12)` and scalar child `(8,8,16,16)` +are physically equivalent to parent set `[(4,4,12,12)]` with child set +`[(12,12,20,20)]`. A child set on the *scalar* parent instead uses +`[(8,8,16,16)]`. Reusing the same numbers after changing the parent frame can +move the physical region. + +Both adjacent pairs require a Chebyshev gap of at least four in their own +parent-cell units. The whole child union must preserve inherited ACTIVE +coverage, uncovered ACTIVE donors/receivers, surviving pulls, and the upstream +restriction footprint. Sparse holes remain unowned; inherited INNER remains +stream-only and OUTER remains invalid. Block sizes must align with every +support edge, including its four-cell ghost band. + +## Example + +```python +run = FixedPatchStepper( + cfg, [(4, 4, 20, 28), (24, 4, 44, 28)], context=ctx, + base_block=(4, 4), fine_block=(4, 4), child_block=(4, 4), + child_rectangles_in_parent=[(16, 16, 32, 40), (56, 16, 80, 40)], +) +``` + +This example uses a 48x32 root. `initialize` and `download` use three Q-major +dense arrays. Public arrays and masks still scale with the full allocated +extent; sparse device buffers do not remove that host-memory limit. + +Stationary whole-circle `child_objects` and `child_jets` use root edge +coordinates and a single L2 boundary object. Body IDs are physical object-list +positions, not region IDs or caller `obj_id`. Complete boundary stencils must +avoid the transfer footprint. Change a child Jet with +`set_jet_amplitude(index, value, level=2)`; L1 bodies/Jets remain unsupported. + +## Verification + +CPU geometry, table, pre-CUDA rejection, schedule and cleanup tests: + +```sh +PYTHONPATH=src python -m unittest discover -s tests/unit -p 'test_refinement*.py' -v +``` + +Explicit GPU gates are in +`tests/validation/multigrid/test_refinement_three_level_multipatch_gpu.py`. +They cover physical-frame equivalence, twins, asymmetric refinement, multiple +children, reversed order, independently varied aligned layouts, public/raw +masks and two L2 bodies with static/zero/driven/switched Jets. Snapshots are +0/1/2/8. CPU mocks do not execute CUDA and do not establish numerical accuracy. +GPU equivalence gates supplement the separately frozen independent dense +three-level numerical replay; passing either alone is not a conservation, +performance or physical-accuracy claim. diff --git a/src/CelerisLab/lbm/refinement.py b/src/CelerisLab/lbm/refinement.py index 288fffa..00aa6ac 100644 --- a/src/CelerisLab/lbm/refinement.py +++ b/src/CelerisLab/lbm/refinement.py @@ -318,11 +318,35 @@ def nest_patch(cfg, parent, child_bounds): would read a covered or OUTER cell without a child-supplied arrival; and child coverage of L1's own two-hop restriction footprint toward L0. """ + child = FixedPatch2D.from_config(_middle_config(cfg, parent), child_bounds) + return child, _nested_kind(parent, child) + + +def nest_patch_set(cfg, parent, child_rectangles_in_parent): + """Separated child union in the parent's allocated L1 array coordinates. + + The root edge coordinate of L1 index (x, y) is + ``parent.fine_origin + (x, y)/2``. For a set parent its origin is zero; + a scalar parent's window includes its ghost band and has a shifted origin. + The child set allocates a dense public extent twice that L1 array, while + its blocked storage omits holes. Both adjacent pairs retain gap >= 4 in + their respective parent-cell units. Validate support over the whole union. + """ + child = FixedPatchSet2D.from_config( + _middle_config(cfg, parent), child_rectangles_in_parent) + return child, _nested_kind(parent, child) + + +def _middle_config(cfg, parent): + fy, fx = parent.fine_shape + return replace(cfg, nx=fx, ny=fy, viscosity=2 * cfg.viscosity, y_wall_bc="bounce_back", + top_wall_bc="velocity", wall_slot_width=0, wall_normal_velocity=0.0, + velocity=0.0) + + +def _nested_kind(parent, child): + """Validate composed transfer support and preserve inherited ghost kinds.""" fy, fx = parent.fine_shape - middle = replace(cfg, nx=fx, ny=fy, viscosity=2 * cfg.viscosity, y_wall_bc="bounce_back", - top_wall_bc="velocity", wall_slot_width=0, wall_normal_velocity=0.0, - velocity=0.0) - child = FixedPatch2D.from_config(middle, child_bounds) kind = parent.fine_kind covered = child.covered active = (kind == ACTIVE) & ~covered @@ -346,7 +370,7 @@ def nest_patch(cfg, parent, child_bounds): raise ValueError("child covers the parent's restriction footprint") composed = kind.copy() composed[covered] = OUTER_GHOST - return child, _readonly(composed, np.uint8) + return _readonly(composed, np.uint8) def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=(), diff --git a/src/CelerisLab/lbm/refinement_stepper.py b/src/CelerisLab/lbm/refinement_stepper.py index 0fd2fcd..f0cc2f2 100644 --- a/src/CelerisLab/lbm/refinement_stepper.py +++ b/src/CelerisLab/lbm/refinement_stepper.py @@ -19,7 +19,7 @@ from .curved_links import surface_link_velocity from ..body.jets import jet_link_shapes, jet_link_velocity from .block_stepper import MAX_GRID_DIM from .refinement import (ACTIVE, HALO, OUTER_GHOST, FixedPatch2D, FixedPatchSet2D, - nest_patch, prepare_fine_boundary) + nest_patch, nest_patch_set, prepare_fine_boundary) from .descriptors import SOLID, D2Q9_EX, D2Q9_EY @@ -61,23 +61,34 @@ class _Pair: parent_cells: int # parent block window cells (explode donor stride) -def blocked_base_tables(patch, layout, order, wall_aware=False): +def blocked_base_tables(patch, layout, order, wall_aware=False, *, composed_kind=None): """Parent donor cells and receiver-local arrival slots, in transfer row order. Stage arrivals only after child reads and parent halo gather, immediately before parent pull. Source *owner* slots are wrong at block seams. + Nested pairs must pass composed_kind to retain inherited INNER/OUTER and + sparse holes; donors and receivers must be owned, uncovered ACTIVE cells. """ owner = layout.owner.ravel() q, receiver = np.divmod(patch.coalesce_coarse[order], np.uint64(owner.size)) + if composed_kind is None: + dense_kind = np.where(patch.covered, OUTER_GHOST, ACTIVE).astype(np.uint8) + if wall_aware: + dense_kind[[0, -1]] = OUTER_GHOST + else: + dense_kind = np.asarray(composed_kind) + if dense_kind.shape != patch.coarse_shape: + raise ValueError("composed parent kind has the wrong shape") + donors = patch.ghost_coarse.astype(np.intp) + receivers = receiver.astype(np.intp) + legal = layout.owned.ravel() & (dense_kind.ravel() == ACTIVE) & ~patch.covered.ravel() + if not legal[donors].all() or not legal[receivers].all(): + raise ValueError("child donors/receivers need owned uncovered parent ACTIVE cells") offset = (np.asarray(D2Q9_EY)[q.astype(np.intp)] * layout.local_shape[1] + np.asarray(D2Q9_EX)[q.astype(np.intp)]) - source = owner[receiver.astype(np.intp)].astype(np.int64) - offset - kind = np.full(layout.nblocks * int(np.prod(layout.local_shape)), OUTER_GHOST, np.uint8) - leaves = ~patch.covered - if wall_aware: - leaves[[0, -1]] = False - kind[owner[leaves.ravel()].astype(np.intp)] = ACTIVE - return (owner[patch.ghost_coarse].astype(np.uint64), layout.slot(source, q), kind) + source = owner[receivers].astype(np.int64) - offset + return (owner[donors].astype(np.uint64), layout.slot(source, q), + blocked_kind(layout, dense_kind)) @dataclass @@ -171,6 +182,16 @@ class FixedPatchStepper: shape, NaN outside the patches' ACTIVE cells. The step sequence is the blocked one above, with every patch's transfers in one launch each. + Three-level separated unions: ``child_rectangles_in_parent=[...]`` is + mutually exclusive with scalar ``child`` and accepts either parent form. + Coordinates index the allocated L1 array, including ghost bands and holes: + root edge = patch.fine_origin + L1 index / 2. A parent set has origin zero; + a scalar parent retains its shifted window origin. The child set's public + extent is twice the allocated L1 shape, even though storage is sparse. + Each pair requires gap >= 4 in its own parent-cell units. All regions share + one buffer pair per level and the fixed 1:2:4 schedule; L2-only bodies/Jets + below apply unchanged, with object IDs assigned over the whole level. + Three-level single chain: ``child`` bounds use L1 window cells; all three levels must be blocked. ``child_objects`` and ``child_jets`` use the L0 edge frame, with L2 dx=dt=.25 and unchanged lattice velocity. L1 bodies @@ -181,11 +202,14 @@ class FixedPatchStepper: """ def __init__(self, cfg, bounds, *, context, coarse_boundary="periodic", fine_objects=(), fine_block=None, base_block=None, jets=(), child=None, child_block=None, - child_objects=(), child_jets=()): + child_objects=(), child_jets=(), child_rectangles_in_parent=None): if HALO != 4: raise RuntimeError("refinement.cu requires a four-cell fine halo") if coarse_boundary not in ("periodic", "slot_reservoir"): raise ValueError("coarse_boundary must be periodic or slot_reservoir") + if child is not None and child_rectangles_in_parent is not None: + raise ValueError("child and child_rectangles_in_parent are mutually exclusive") + nested = child is not None or child_rectangles_in_parent is not None bounds = tuple(bounds) multipatch = len(bounds) > 0 and not isinstance(bounds[0], (Integral, np.integer)) if multipatch: @@ -205,7 +229,7 @@ class FixedPatchStepper: raise ValueError("slot_reservoir requires y0 >= 4 and y1 <= ny-4") fine_objects, jets = tuple(fine_objects), tuple(jets) child_objects, child_jets = tuple(child_objects), tuple(child_jets) - if child is None and (child_objects or child_jets): + if not nested and (child_objects or child_jets): raise ValueError("child_objects/child_jets require a nested child") if base_block is not None and fine_block is None: raise ValueError("base_block requires fine_block") @@ -230,17 +254,23 @@ class FixedPatchStepper: raise ValueError("fine blocks exceed launch limits") self.child_patch = self.child_layout = None self._middle_kind = None - if child is not None: - if multipatch or fine_block is None or base_block is None or child_block is None: - raise ValueError("a nested child needs one patch plus fine_block, base_block " - "and child_block") + if nested: + if child is not None and multipatch: + raise ValueError("scalar child requires one parent patch; use child_rectangles_in_parent") + if fine_block is None or base_block is None or child_block is None: + raise ValueError("a nested child needs fine_block, base_block and child_block") if tuple(fine_objects) or tuple(jets): raise ValueError("a nested child supports bodies/Jets only on L2; use child_objects/child_jets") if len(child_block) != 2: raise ValueError("child_block must be (by, bx)") - self.child_patch, self._middle_kind = nest_patch(cfg, self.patch, child) - self.child_layout = BlockLayout2D.build(*self.child_patch.fine_shape, *child_block, - periodic=(False, False)) + if child_rectangles_in_parent is not None: + self.child_patch, self._middle_kind = nest_patch_set( + cfg, self.patch, child_rectangles_in_parent) + self.child_layout = self.child_patch.fine_layout(child_block) + else: + self.child_patch, self._middle_kind = nest_patch(cfg, self.patch, child) + self.child_layout = BlockLayout2D.build(*self.child_patch.fine_shape, *child_block, + periodic=(False, False)) if (self.child_layout.nblocks > MAX_GRID_DIM or self.child_layout.local_shape[0] > MAX_GRID_DIM): raise ValueError("child blocks exceed launch limits") @@ -365,7 +395,8 @@ class FixedPatchStepper: g, cp = self._levels[2], self.child_patch order = np.argsort(cp.coalesce_coarse) ghosts, children, g_kind = blocked_fine_tables(cp, g.layout, order) - donors, targets, _ = blocked_base_tables(cp, f.layout, order) + donors, targets, _ = blocked_base_tables( + cp, f.layout, order, composed_kind=self._middle_kind) g.kind = self._upload(g_kind) g.halo_copy = g.module.get_function("BlockHaloCopy") self._child = _Pair(g.module.get_function("BlockExplode"), diff --git a/tests/unit/test_refinement_three_level_multipatch.py b/tests/unit/test_refinement_three_level_multipatch.py new file mode 100644 index 0000000..e5068b4 --- /dev/null +++ b/tests/unit/test_refinement_three_level_multipatch.py @@ -0,0 +1,314 @@ +"""CPU-only union geometry and production host-control-flow tests. + +The AST loader omits CUDA imports and supplies explicit mocks; it executes +unchanged production definitions, never CUDA kernels or a numerical emulator. +""" +import ast +from contextlib import nullcontext +from dataclasses import replace +from pathlib import Path +import types +import unittest +from unittest.mock import Mock, patch + +import numpy as np +from CelerisLab.config import LBMConfig +from CelerisLab.lbm import refinement as r + + +def host_module(): + path=Path(r.__file__).with_name('refinement_stepper.py') + tree=ast.parse(path.read_text()) + tree.body=[node for node in tree.body if not ( + isinstance(node,ast.Import) and any(a.name=='pycuda.driver' for a in node.names) + or isinstance(node,ast.ImportFrom) and node.module in ( + 'cuda','block_stepper','curved_links'))] + module=types.ModuleType('CelerisLab.lbm._test_host_refinement') + module.__package__='CelerisLab.lbm' + from CelerisLab.body.jets import surface_link_velocity + module.__dict__.update(cuda=Mock(),compiler=Mock(),MAX_GRID_DIM=65535, + surface_link_velocity=surface_link_velocity) + exec(compile(tree,str(path),'exec'),module.__dict__) + return module + + +host=host_module() +CFG = LBMConfig(nx=48, ny=32, viscosity=.1) +PARENTS = ((4,4,20,28), (24,4,44,28)) # exact gap=4 coarse cells +CHILDREN = ((16,16,32,40), (56,16,80,40)) # parent-allocated L1 cells + + +class Gate(unittest.TestCase): + def setUp(self): + self.parent = r.FixedPatchSet2D.from_config(CFG, PARENTS) + + def test_one_element_equivalence(self): + cfg = LBMConfig(nx=16,ny=16,viscosity=.1) + local = r.FixedPatch2D.from_config(cfg,(4,4,12,12)) + global_parent = r.FixedPatchSet2D.from_config(cfg,[(4,4,12,12)]) + a, ak = r.nest_patch(cfg,local,(8,8,16,16)) + b, bk = r.nest_patch_set(cfg,global_parent,[(12,12,20,20)]) + local_set, local_kind = r.nest_patch_set(cfg,local,[(8,8,16,16)]) + np.testing.assert_array_equal(ak,local_kind) + np.testing.assert_array_equal(a.fine_kind,local_set.patches[0].fine_kind) + np.testing.assert_array_equal(a.ghost_coarse,local_set.ghost_coarse) + np.testing.assert_array_equal(a.coalesce_coarse,local_set.coalesce_coarse) + np.testing.assert_array_equal(ak, bk[4:28,4:28]) + # L2 node-index offsets in global L2 grid: local parent origin=2 L0. + ay,ax=np.nonzero(a.fine_kind == r.ACTIVE) + by,bx=np.nonzero(b.fine_kind == r.ACTIVE) + np.testing.assert_array_equal(np.stack([ay+20,ax+20],-1),np.stack([by,bx],-1)) + np.testing.assert_array_equal(a.fine_kind,b.fine_kind[20:44,20:44]) + # Compare complete per-region maps after undoing set addressing. + bp=b.patches[0] + np.testing.assert_array_equal(a.coalesce_fine,bp.coalesce_fine) + aq,ac=np.divmod(a.coalesce_coarse,24*24) + bq,bc=np.divmod(bp.coalesce_coarse,32*32) + np.testing.assert_array_equal(aq,bq) + np.testing.assert_array_equal(np.stack(np.divmod(ac,24),-1)+4,np.stack(np.divmod(bc,32),-1)) + ay,ax=np.divmod(a.ghost_coarse,24) + by,bx=np.divmod(bp.ghost_coarse,32) + np.testing.assert_array_equal(np.stack([ay+4,ax+4],-1),np.stack([by,bx],-1)) + + def test_twin_tile_chains_and_order(self): + children=((16,16,32,40),(56,16,72,40)) + child,kind=r.nest_patch_set(CFG,self.parent,children) + other,ok=r.nest_patch_set(CFG,r.FixedPatchSet2D.from_config(CFG,PARENTS[::-1]),children[::-1]) + for name in ('fine_kind','ghost_fine','ghost_coarse','coalesce_coarse','coalesce_fine'): + np.testing.assert_array_equal(getattr(child,name),getattr(other,name)) + np.testing.assert_array_equal(kind,ok) + for p in child.patches: + self.assertEqual(p.fine_shape,child.patches[0].fine_shape) + np.testing.assert_array_equal(child.patches[0].fine_kind,child.patches[1].fine_kind) + self.assertEqual(child.tau_coarse,self.parent.tau_fine) + self.assertAlmostEqual(child.tau_fine,1.7) + + def test_asymmetric_no_child_region(self): + child,kind=r.nest_patch_set(CFG,self.parent,[CHILDREN[0]]) + np.testing.assert_array_equal(kind[:,48:],self.parent.fine_kind[:,48:]) + self.assertEqual(int(child.covered.sum()),16*24) + + def test_two_children_one_parent(self): + child,kind=r.nest_patch_set(CFG,self.parent,[(16,16,32,24),(16,28,32,40)]) + self.assertEqual(len(child.patches),2) + self.assertEqual(int(child.covered.sum()),16*20) + self.assertEqual(int((kind==r.ACTIVE).sum()),int((self.parent.fine_kind==r.ACTIVE).sum())-16*20) + + def test_holes_donors_receivers_footprint_and_invalid_pull(self): + child,_=r.nest_patch_set(CFG,self.parent,CHILDREN) + cases=[('inherited ACTIVE', lambda: r.nest_patch_set(CFG,self.parent,[(40,16,48,40)])), + ('donors/receivers', lambda: r.nest_patch_set(CFG,self.parent,[(8,16,32,40)]))] + for text,call in cases: + with self.subTest(text=text),self.assertRaisesRegex(ValueError,text):call() + # Synthetic upstream dependency inside child coverage isolates the final guard. + upstream=self.parent.coalesce_fine.copy() + upstream[0,:]=20*self.parent.fine_shape[1]+20 # q=0, deep inherited ACTIVE + with self.assertRaisesRegex(ValueError,'restriction footprint'): + r._nested_kind(replace(self.parent,coalesce_fine=upstream),child) + hole=int(np.flatnonzero(~self.parent.fine_layout((4,4)).owned.ravel())[0]) + donors=child.ghost_coarse.copy();donors[0]=hole + with self.assertRaisesRegex(ValueError,'donors/receivers'): + r._nested_kind(self.parent,replace(child,ghost_coarse=donors)) + receivers=child.coalesce_coarse.copy();receivers[0]=hole + with self.assertRaisesRegex(ValueError,'donors/receivers'): + r._nested_kind(self.parent,replace(child,coalesce_coarse=receivers)) + # Remove all supplying rows; covered-to-uncovered pulls must be rejected. + with self.assertRaisesRegex(ValueError,'surviving parent ACTIVE pull'): + r._nested_kind(self.parent,replace(child,coalesce_coarse=child.coalesce_coarse[:0],coalesce_fine=child.coalesce_fine[:0])) + + def test_existing_l2_body_jet_preparation_on_union(self): + from CelerisLab.body.geometry.circle import CircleGeometry + from CelerisLab.body.objects import SimObject + from CelerisLab.body.jets import Jet, Tube, Parallel, jet_link_shapes + child,_=r.nest_patch_set(CFG,self.parent,CHILDREN) + centers=((12.,14.),(34.,14.)) + objects=[SimObject(99,geometry=CircleGeometry(*c,.8),center=c,radius=.8) for c in centers] + jets=[Jet(Tube((c[0]-1.5,c[1]),0.,3.,.3),Parallel(0.),amplitude=.002) for c in centers] + middle=replace(CFG,nx=self.parent.fine_shape[1],ny=self.parent.fine_shape[0],viscosity=.2) + boundary=r.prepare_fine_boundary(middle,child,objects,jets=jets, + parent_origin=self.parent.fine_origin,parent_spacing=.5) + ids=boundary.curved_host[3] + np.testing.assert_array_equal(np.unique(ids),[0,1]) + self.assertEqual(int((ids==0).sum()),int((ids==1).sum())) + np.testing.assert_allclose(boundary.curved_hits[ids==1]-boundary.curved_hits[ids==0], + np.broadcast_to([22.,0.],boundary.curved_hits[ids==0].shape),atol=1e-14,rtol=0) + _,jid=jet_link_shapes(jets,boundary.curved_hits) + for i in (0,1): + self.assertTrue((ids[jid==i]==i).all()) + self.assertTrue((jid==i).any()) + + def test_sparse_parent_tables_preserve_inherited_kinds(self): + child,kind=r.nest_patch_set(CFG,self.parent,CHILDREN) + layout=self.parent.fine_layout((4,4));order=np.argsort(child.coalesce_coarse) + donors,targets,bkind=host.blocked_base_tables(child,layout,order,composed_kind=kind) + self.assertTrue((donors < layout.hole).all()) + self.assertTrue((targets < np.prod(layout.blocked_shape)).all()) + np.testing.assert_array_equal(bkind[layout.owner[layout.owned]],kind[layout.owned]) + inherited=self.parent.fine_kind != r.ACTIVE + np.testing.assert_array_equal(kind[inherited],self.parent.fine_kind[inherited]) + self.assertTrue((kind[child.covered] == r.OUTER_GHOST).all()) + l2=child.fine_layout((4,4)) + ghosts,children,ckind=host.blocked_fine_tables(child,l2,order) + self.assertTrue((ghosts20: + raise RuntimeError('table upload failed') + module.cuda.memcpy_htod.side_effect=upload + objects,jets=bodies_and_jets() if bodies else ((),()) + try: + run=module.FixedPatchStepper(CFG,PARENTS,context=context, + child_objects=objects,child_jets=jets, + **dict(KW,child_rectangles_in_parent=children)) + except RuntimeError: + for allocation in allocations: + allocation.free.assert_called_once_with() + raise + return run,module,allocations + + +class ConstructorAndSchedule(unittest.TestCase): + def test_invalid_inputs_fail_before_context_or_allocations(self): + context=Mock() + cases=[dict(child=(8,8,16,16)),dict(child_rectangles_in_parent=[]), + dict(child_rectangles_in_parent=[(16,16,32)]), + dict(child_rectangles_in_parent=[(16.5,16,32,40)]), + dict(child_rectangles_in_parent=[(16,16,32,24),(16,27,32,40)]), + dict(child_rectangles_in_parent=[(40,16,48,40)]), + dict(child_rectangles_in_parent=[(8,16,32,40)]), + dict(child_block=None),dict(fine_block=None),dict(base_block=None), + dict(child_block=(5,4)),dict(fine_block=(5,4)), + dict(child_block=(4,)),dict(child_rectangles_in_parent=(16,16,32,40))] + for changes in cases: + with self.subTest(changes=changes),self.assertRaises((ValueError,TypeError)): + host.FixedPatchStepper(CFG,PARENTS,context=context,**dict(KW,**changes)) + with self.assertRaisesRegex(ValueError,'closer'): + host.FixedPatchStepper(CFG,[(4,4,20,28),(23,4,44,28)],context=context,**KW) + objects,jets=bodies_and_jets() + with self.assertRaisesRegex(ValueError,'only on L2'): + host.FixedPatchStepper(CFG,PARENTS,context=context,fine_objects=objects,**KW) + from CelerisLab.body.jets import Jet,Tube,Parallel + spanning=Jet(Tube((10.5,14.),0.,25.,.3),Parallel(0.),amplitude=0.) + with self.assertRaisesRegex(ValueError,'one jet'): + host.FixedPatchStepper(CFG,PARENTS,context=context,child_objects=objects, + child_jets=[spanning],**KW) + context.activate.assert_not_called() + host.cuda.mem_alloc.assert_not_called() + + def test_aggregate_schedule_independent_of_region_count(self): + for children,bodies in ((CHILDREN,True),(CHILDREN[:1],False), + (((16,16,32,24),(16,28,32,40)),False)): + run,_,allocations=mock_backend(children=children,bodies=bodies) + self.assertEqual(len(run._levels),3) + c,f,g=run._levels + self.assertEqual(len({id(level.current) for level in run._levels} | + {id(level.spare) for level in run._levels}),6) + events=[] + initial_c=c.current + def record(name): + return lambda source,*a,**kw: events.append((name,source)) + run._explode.side_effect=record('explode01') + run._child.explode.side_effect=record('explode12') + for i,level in enumerate(run._levels): + level.kernel.side_effect=record(f'step{i}') + if bodies: + g.boundary.kernel.side_effect=record('bfl2') + run._child.coalesce.side_effect=lambda src,dst,*a,**kw: events.append(('restrict21',dst)) + run._coalesce.side_effect=lambda src,dst,*a,**kw: events.append(('restrict10',dst)) + run._status='ready' + expected=[('explode01',initial_c)] + for middle in (f.current,f.spare): + expected.append(('explode12',middle)) + for leaf in (g.current,g.spare): + expected.append(('halo2',leaf)) + if bodies: + expected.append(('bfl2',leaf)) + expected.append(('step2',leaf)) + expected.extend([('halo1',middle),('restrict21',middle),('step1',middle)]) + expected.extend([('halo0',initial_c),('restrict10',initial_c),('step0',initial_c)]) + def gather(level,buffer,**kw): + events.append((f'halo{next(i for i,l in enumerate(run._levels) if l is level)}',buffer)) + with patch.object(run,'_gather_halo',side_effect=gather): + run.step() + self.assertEqual(events,expected) + self.assertEqual([l.kernel.call_count for l in run._levels],[1,2,4]) + # A stale-feed implementation must fail the source-identity guard. + feeds=[] + advance=run._advance_child + def checked_advance(buffer, **kw): + self.assertIs(buffer,(f.current,f.spare)[len(feeds)]) + feeds.append(buffer) + advance(buffer,**kw) + with patch.object(run,'_advance_child',side_effect=lambda buffer,**kw: + checked_advance(f.current,**kw)): + with self.assertRaises(AssertionError): + run.step() + # A per-region shared-parent advance must fail the 1:2:4 guard. + run._status='ready' + counts=[0,0,0] + def count(index, multiple=1): + def step(*a,**kw): + counts[index]+=multiple + return step + c.kernel.side_effect=count(0) + f.kernel.side_effect=count(1,2) # mutant: twice per shared L1 hop + g.kernel.side_effect=count(2) + run.step() + with self.assertRaises(AssertionError): + self.assertEqual(counts,[1,2,4]) + run.close();run.close() + for allocation in allocations: + allocation.free.assert_called_once_with() + + def test_partial_construction_frees_every_allocation(self): + for options in (dict(fail_compile=True),dict(fail_upload=True)): + with self.subTest(options=options),self.assertRaisesRegex(RuntimeError,'failed'): + mock_backend(**options) + + +if __name__=='__main__': + unittest.main() diff --git a/tests/validation/multigrid/test_refinement_three_level_multipatch_gpu.py b/tests/validation/multigrid/test_refinement_three_level_multipatch_gpu.py new file mode 100644 index 0000000..5b21c66 --- /dev/null +++ b/tests/validation/multigrid/test_refinement_three_level_multipatch_gpu.py @@ -0,0 +1,178 @@ +"""Frozen GPU gates for the bounded three-level separated-union backend. + +Run explicitly on an available GPU; CPU geometry/mocks live in tests/unit. +Snapshots 0/1/2/8: finite public leaves, exact NaN masks, raw invalid cores; +reversed order and aligned layout changes must be bit-identical. The scalar +window comparison uses root physical coordinates, never equal array indices. +These equivalence checks supplement, not replace, independent dense CPU replay. +""" +from dataclasses import replace +import unittest + +import numpy as np + +from CelerisLab.config import LBMConfig +from CelerisLab.lbm.descriptors import SOLID +from CelerisLab.lbm.equilibrium import compute_feq_d2q9 +from CelerisLab.lbm.refinement import OUTER_GHOST +from CelerisLab.lbm.refinement_stepper import FixedPatchStepper +from test_refinement_blocked_base import raw, exact + +CFG = LBMConfig(nx=48, ny=32, viscosity=.1) +PARENTS = ((4, 4, 20, 28), (24, 4, 44, 28)) +CHILDREN = ((16, 16, 32, 40), (56, 16, 80, 40)) +TOPOLOGIES = { + 'twins': ((16, 16, 32, 40), (56, 16, 72, 40)), + 'asymmetric': (CHILDREN[0],), + 'multiple_children': ((16, 16, 32, 24), (16, 28, 32, 40)), +} +# 8-wide blocks require support edges divisible by 8, i.e. rectangle edges +# equal 2 modulo 4 in each pair's parent frame. The default fixture is 4-only. +ALIGNED_CFG = LBMConfig(nx=56, ny=40, viscosity=.1) +ALIGNED_PARENTS = ((6, 6, 22, 34), (26, 6, 50, 34)) +ALIGNED_CHILDREN = ((18, 18, 34, 50), (58, 18, 90, 50)) +STOPS = (0, 1, 2, 8) + + +def config(cfg, mode): + return cfg if mode == 'periodic' else replace( + cfg, velocity=0., wall_normal_velocity=.001, wall_slot_width=4, + y_wall_bc='on_site_velocity', top_wall_bc='reservoir') + + +def bodies(): + from CelerisLab.body.geometry.circle import CircleGeometry + from CelerisLab.body.objects import SimObject + return [SimObject(99, geometry=CircleGeometry(x, 14., .8), center=(x, 14.), radius=.8) + for x in (12., 34.)] + + +def jets(speed): + from CelerisLab.body.jets import Jet, Tube, Parallel + return [Jet(Tube((x-1.5, 14.), 0., 3., .3), Parallel(0.), amplitude=speed) + for x in (12., 34.)] + + +def frames(run): + middle_origin = np.asarray(run.patch.fine_origin) + return ((np.zeros(2), 1.), (middle_origin, .5), + (middle_origin + .5*np.asarray(run.child_patch.fine_origin), .25)) + + +def physical_state(shape, origin, spacing): + """Same smooth non-equilibrium G at the same physical cell centers.""" + y, x = np.indices(shape, dtype=np.float64) + x, y = origin[0] + spacing*(x+.5), origin[1] + spacing*(y+.5) + equilibrium = np.asarray(compute_feq_d2q9(1., .013, -.007), np.float64)[:, None, None] + q = np.arange(9)[:, None, None] + return (equilibrium + 1e-4*np.sin(.17*x + .23*y + q)).astype(np.float32) + + +def physical_leaves(array, origin, spacing): + y, x = np.nonzero(np.isfinite(array[0])) + coordinates = np.stack((origin[1]+spacing*(y+.5), origin[0]+spacing*(x+.5)), axis=-1) + return coordinates, array[:, y, x] + + +class ThreeLevelMultipatchGPU(unittest.TestCase): + def run_case(self, ctx, cfg=CFG, parents=PARENTS, children=CHILDREN, + mode='periodic', blocks=((4, 4), (4, 4)), objects=(), drive=(), switch=False, + scalar_child=None): + kwargs = dict(child_rectangles_in_parent=children) if scalar_child is None else dict(child=scalar_child) + with FixedPatchStepper(config(cfg, mode), parents, context=ctx, coarse_boundary=mode, + base_block=(4, 4), fine_block=blocks[0], child_block=blocks[1], + child_objects=objects, child_jets=drive, **kwargs) as run: + frame = frames(run) + run.initialize(*(physical_state(level.shape[1:], *f) + for level, f in zip(run._levels, frame))) + snapshots = [] + for stop in STOPS: + run.step(stop-run.step_count) + if switch and stop == 2: + for index in range(len(drive)): + run.set_jet_amplitude(index, -.002, level=2) + public, stored = run.download(), raw(run) + masks = [~run.patch.covered, run._fine_leaf, run._child_leaf] + if mode != 'periodic': + masks[0][[0, -1]] = False + invalid = [run.patch.covered.copy(), run._middle_kind == OUTER_GHOST, + run.child_patch.fine_kind == OUTER_GHOST] + if mode != 'periodic': + invalid[0][[0, -1]] = True + if run.child_boundary is not None: + invalid[2] |= (run.child_boundary.flags.reshape(run.child_patch.fine_shape) & SOLID) != 0 + for level, (g, mask, dense, bad) in enumerate(zip(public, masks, stored, invalid)): + label = f'{mode} stop={stop} level={level}' + self.assertTrue(np.isfinite(g[:, mask]).all(), label) + self.assertTrue(np.isnan(g[:, ~mask]).all(), label) + self.assertTrue(np.isnan(dense[:, bad]).all(), label) + snapshots.append((public, stored, frame)) + return snapshots + + def compare(self, actual, expected): + for stop, a, b in zip(STOPS, actual, expected): + for kind in (0, 1): + for level, (x, y) in enumerate(zip(a[kind], b[kind])): + exact(self, x, y, f'stop={stop} kind={kind} level={level}') + + def test_topologies_reversed_order_and_masks(self): + from CelerisLab.cuda.context import CudaContext + with CudaContext(0) as ctx: + for mode in ('periodic', 'slot_reservoir'): + for name, children in TOPOLOGIES.items(): + with self.subTest(mode=mode, topology=name): + expected = self.run_case(ctx, children=children, mode=mode) + actual = self.run_case(ctx, parents=PARENTS[::-1], children=children[::-1], mode=mode) + self.compare(actual, expected) + + def test_independently_varied_aligned_layouts(self): + from CelerisLab.cuda.context import CudaContext + with CudaContext(0) as ctx: + for mode in ('periodic', 'slot_reservoir'): + kwargs = dict(cfg=ALIGNED_CFG, parents=ALIGNED_PARENTS, + children=ALIGNED_CHILDREN, mode=mode) + expected = self.run_case(ctx, **kwargs) + for middle in ((4, 4), (8, 4), (4, 8)): + for leaf in ((4, 4), (8, 4), (4, 8)): + self.compare(self.run_case(ctx, blocks=(middle, leaf), **kwargs), expected) + + def test_one_element_equivalent_physical_frames(self): + from CelerisLab.cuda.context import CudaContext + cfg = LBMConfig(nx=16, ny=16, viscosity=.1) + parent = (4, 4, 12, 12) + with CudaContext(0) as ctx: + for mode in ('periodic', 'slot_reservoir'): + old = self.run_case(ctx, cfg=cfg, parents=parent, mode=mode, scalar_child=(8, 8, 16, 16)) + # Test both set-on-scalar and set-on-set allocated frames. + for parents, child in ((parent, (8, 8, 16, 16)), ((parent,), (12, 12, 20, 20))): + new = self.run_case(ctx, cfg=cfg, parents=parents, children=[child], mode=mode) + for a, b in zip(new, old): + for level in range(3): + ax, av = physical_leaves(a[0][level], *a[2][level]) + bx, bv = physical_leaves(b[0][level], *b[2][level]) + np.testing.assert_array_equal(ax, bx) + exact(self, av, bv, f'physical level={level}') + + def test_two_l2_bodies_zero_driven_switched_order_and_layout(self): + from CelerisLab.cuda.context import CudaContext + with CudaContext(0) as ctx: + for mode in ('periodic', 'slot_reservoir'): + kwargs = dict(cfg=ALIGNED_CFG, parents=ALIGNED_PARENTS, + children=ALIGNED_CHILDREN, mode=mode, objects=bodies()) + static = self.run_case(ctx, **kwargs) + zero = self.run_case(ctx, drive=jets(0.), **kwargs) + self.compare(zero, static) + driven = self.run_case(ctx, drive=jets(.002), **kwargs) + switched = self.run_case(ctx, drive=jets(.002), switch=True, **kwargs) + self.assertFalse(np.array_equal(driven[-1][0][2], static[-1][0][2], equal_nan=True)) + self.assertFalse(np.array_equal(driven[-1][0][2], switched[-1][0][2], equal_nan=True)) + for drive, switch, expected in (((), False, static), (jets(0.), False, zero), + (jets(.002), False, driven), (jets(.002), True, switched)): + for middle, leaf in (((8, 4), (4, 8)), ((4, 8), (8, 4))): + actual = self.run_case(ctx, blocks=(middle, leaf), drive=drive, switch=switch, + **dict(kwargs, parents=ALIGNED_PARENTS[::-1], children=ALIGNED_CHILDREN[::-1])) + self.compare(actual, expected) + + +if __name__ == '__main__': + unittest.main()