From fbe73779c1767658357ee3e876289ea7ad4fc109 Mon Sep 17 00:00:00 2001 From: Codex Date: Fri, 16 Oct 2026 20:51:29 -0700 Subject: [PATCH] Build refinement transfers from union and pin physical body support --- src/CelerisLab/lbm/refinement.py | 162 +++++++++++------- src/CelerisLab/lbm/refinement_stepper.py | 15 +- src/CelerisLab/lbm/synchronized_remap.py | 22 +-- tests/unit/test_hierarchy.py | 4 +- tests/unit/test_multiple_pinned_remap.py | 8 +- tests/unit/test_pinned_body_remap.py | 8 +- tests/unit/test_refinement.py | 24 +-- .../test_refinement_three_level_multipatch.py | 6 +- tests/unit/test_runtime_child_remap.py | 2 +- tests/unit/test_synchronized_remap.py | 3 +- tests/unit/test_union_body_support.py | 44 +++++ tests/unit/test_union_refinement.py | 161 +++++++++++++++++ .../multigrid/test_refinement_multipatch.py | 4 +- 13 files changed, 357 insertions(+), 106 deletions(-) create mode 100644 tests/unit/test_union_body_support.py create mode 100644 tests/unit/test_union_refinement.py diff --git a/src/CelerisLab/lbm/refinement.py b/src/CelerisLab/lbm/refinement.py index 00aa6ac..8fe30b2 100644 --- a/src/CelerisLab/lbm/refinement.py +++ b/src/CelerisLab/lbm/refinement.py @@ -151,20 +151,17 @@ class FixedPatch2D: raise ValueError("levels are not synchronized: fine_steps must equal 2*coarse_steps") -MIN_PATCH_GAP = 4 # coarse cells: the full 4-fine-cell supports stay disjoint +MIN_PATCH_GAP = 4 # Legacy exported constant; union transfers impose no gap policy. @dataclass(frozen=True, eq=False) class FixedPatchSet2D: - """Several separated FixedPatch2D rectangles sharing one 2:1 fine level. + """Rectangles describing one 2:1 fine union, with unique transfer ownership. - The fine level is the whole (2*ny, 2*nx) grid in level coordinates; only - the patch supports (active region + HALO ghost band) are stored, by a - sparse BlockLayout2D from fine_layout(). Tables have FixedPatch2D meaning - and row/child order, re-addressed to level cells: ghost_fine are level - fine cells, coalesce_fine are q*N+cell level slots; coarse ids unchanged. - Patches are kept in canonical (y0, x0) order, so input order is irrelevant. - Active regions must be >= MIN_PATCH_GAP coarse cells apart (Chebyshev). + Input patches remain canonical metadata. Numerical kinds and tables are + built once from their union; adjacent, overlapping and concave unions are + accepted when their actual transfer dependencies close. Sparse storage owns + every whole tile intersecting the support, independent of rectangle edges. """ coarse_shape: tuple[int, int] patches: tuple @@ -182,42 +179,43 @@ class FixedPatchSet2D: if not rectangles: raise ValueError("at least one patch rectangle is required") patches = sorted((FixedPatch2D.from_config(cfg, r) for r in rectangles), - key=lambda p: (p.bounds[1], p.bounds[0])) - for i, a in enumerate(patches): - for b in patches[i + 1:]: - (ax0, ay0, ax1, ay1), (bx0, by0, bx1, by1) = a.bounds, b.bounds - if max(ay0 - by1, by0 - ay1, ax0 - bx1, bx0 - ax1) < MIN_PATCH_GAP: - raise ValueError(f"patches {a.bounds} and {b.bounds} are closer " - f"than {MIN_PATCH_GAP} coarse cells") + key=lambda p: (p.bounds[1], p.bounds[0], p.bounds[3], p.bounds[2])) ny, nx = patches[0].coarse_shape - fy, fx = 2 * ny, 2 * nx - if fy * fx * 9 > np.iinfo(np.uint64).max // 2: - raise ValueError("fine level exceeds 64-bit indexing") - kind = np.full((fy, fx), OUTER_GHOST, np.uint8) - ghost_fine, coalesce_fine = [], [] - for p in patches: - y0, x0 = p.fine_level_origin - ly, lx = p.fine_shape - kind[y0:y0 + ly, x0:x0 + lx] = p.fine_kind - gy, gx = np.divmod(p.ghost_fine.astype(np.int64), lx) - ghost_fine.append((gy + y0) * fx + gx + x0) - q, cell = np.divmod(p.coalesce_fine.astype(np.int64), ly * lx) - cy, cx = np.divmod(cell, lx) - coalesce_fine.append(q * fy * fx + (cy + y0) * fx + cx + x0) - ghost_coarse = np.concatenate([p.ghost_coarse for p in patches]) - coalesce_coarse = np.concatenate([p.coalesce_coarse for p in patches]) covered = np.zeros((ny, nx), bool) for p in patches: covered |= p.covered - cells = ny * nx - if (covered.ravel()[ghost_coarse].any() - or covered.ravel()[coalesce_coarse % np.uint64(cells)].any() - or len(np.unique(coalesce_coarse)) != len(coalesce_coarse)): - raise ValueError("patch transfers overlap another patch") - return cls((ny, nx), tuple(patches), patches[0].tau_coarse, patches[0].tau_fine, - _readonly(kind, np.uint8), _readonly(np.concatenate(ghost_fine), np.uint64), - _readonly(ghost_coarse, np.uint32), _readonly(coalesce_coarse, np.uint64), - _readonly(np.concatenate(coalesce_fine), np.uint64).reshape(-1, 4)) + fy, fx = 2 * ny, 2 * nx + if fy * fx * 9 > np.iinfo(np.uint64).max // 2: + raise ValueError("fine level exceeds 64-bit indexing") + active = np.kron(covered, np.ones((2, 2), bool)) + inner = _dilate(active, 2) + support = _dilate(active, HALO) + kind = np.full((fy, fx), OUTER_GHOST, np.uint8) + kind[inner] = INNER_GHOST + kind[active] = ACTIVE + gy, gx = np.nonzero(support & ~active) + ghost_fine = (gy * fx + gx).astype(np.uint64) + ghost_coarse = ((gy >> 1) * nx + (gx >> 1)).astype(np.uint32) + targets, children = [], [] + ry, rx = np.nonzero(_dilate(covered, 1) & ~covered) + for y, x in zip(ry.tolist(), rx.tolist()): + for q in range(1, 9): + if _closed_union_contains(covered, 2 * x + 1 - D2Q9_EX[q], 2 * y + 1 - D2Q9_EY[q]): + targets.append(q * ny * nx + y * nx + x) + children.append([q * fy * fx + (2 * y + dy) * fx + 2 * x + dx + for dy in (0, 1) for dx in (0, 1)]) + targets = np.asarray(targets, np.uint64) + children = np.asarray(children, np.uint64).reshape(-1, 4) + cells = np.asarray(children % np.uint64(fy * fx), np.intp) + if not (kind.ravel()[cells] == INNER_GHOST).all(): + raise ValueError("coalesce children must be INNER ghost cells") + pair = cls((ny, nx), tuple(patches), patches[0].tau_coarse, patches[0].tau_fine, + _readonly(kind, np.uint8), _readonly(ghost_fine, np.uint64), + _readonly(ghost_coarse, np.uint32), _readonly(targets, np.uint64), + _readonly(children, np.uint64)) + _transfer_footprint(pair) + _check_parent_pull_closure(pair) + return pair @property def fine_shape(self) -> tuple[int, int]: @@ -240,22 +238,61 @@ class FixedPatchSet2D: return mask def fine_layout(self, block): - """Sparse non-periodic layout storing exactly the patch supports.""" + """Sparse non-periodic tiles covering union support; domain divisibility remains.""" from .blocks import BlockLayout2D block = tuple(block) if len(block) != 2 or any(isinstance(v, (bool, np.bool_)) or not isinstance(v, Integral) or v <= 0 for v in block): raise ValueError("fine block must be two positive integers (by, bx)") by, bx = map(int, block) - origins = set() - for p in self.patches: - y0, x0 = p.fine_level_origin - ly, lx = p.fine_shape - if y0 % by or x0 % bx or ly % by or lx % bx: - raise ValueError(f"patch {p.bounds} support is not aligned to fine blocks {block}") - origins.update((y, x) for y in range(y0, y0 + ly, by) for x in range(x0, x0 + lx, bx)) + support = _dilate(self.fine_kind == ACTIVE, HALO) + ty, tx = np.nonzero(support) + origins = sorted({(int(y) // by * by, int(x) // bx * bx) for y, x in zip(ty, tx)}) return BlockLayout2D.build(*self.fine_shape, by, bx, periodic=(False, False), - origins=sorted(origins)) + origins=origins) + + +def _dilate(mask, steps): + """Chebyshev dilation with a zero (non-periodic) exterior.""" + out = np.asarray(mask, bool) + for _ in range(steps): + grown = out.copy() + grown[1:] |= out[:-1] + grown[:-1] |= out[1:] + wide = grown.copy() + wide[:, 1:] |= grown[:, :-1] + wide[:, :-1] |= grown[:, 1:] + out = wide + return out + + +def _closed_union_contains(covered, mx, my): + """Does doubled point (mx, my) lie in a closed covered square [2X, 2X+2] x [2Y, 2Y+2]?""" + xs = (mx // 2 - 1, mx // 2) if mx % 2 == 0 else ((mx - 1) // 2,) + ys = (my // 2 - 1, my // 2) if my % 2 == 0 else ((my - 1) // 2,) + ny, nx = covered.shape + return any(0 <= y < ny and 0 <= x < nx and covered[y, x] for y in ys for x in xs) + + +def _check_parent_pull_closure(pair): + """Every uncovered parent cell pulling q from a covered cell needs its (R, q) row; no row + may stage into an uncovered parent cell unless its link passes the closed union.""" + ny, nx = pair.coarse_shape + rows = set(map(int, pair.coalesce_coarse)) + if len(rows) != len(pair.coalesce_coarse): + raise ValueError("duplicate coalesce writer (receiver, q)") + covered = pair.covered + for y, x in zip(*np.nonzero(~covered)): + for q in range(1, 9): + sy, sx = y - D2Q9_EY[q], x - D2Q9_EX[q] + if 0 <= sy < ny and 0 <= sx < nx and covered[sy, sx]: + if q * ny * nx + y * nx + x not in rows: + raise ValueError("a surviving parent ACTIVE pull lacks its coalesce row") + for slot in rows: + q, cell = divmod(slot, ny * nx) + y, x = divmod(cell, nx) + if covered[y, x] or not _closed_union_contains(covered, 2*x + 1 - D2Q9_EX[q], 2*y + 1 - D2Q9_EY[q]): + raise ValueError("coalesce row outside the union boundary") @dataclass(frozen=True, eq=False) @@ -323,14 +360,14 @@ def nest_patch(cfg, parent, child_bounds): def nest_patch_set(cfg, parent, child_rectangles_in_parent): - """Separated child union in the parent's allocated L1 array coordinates. + """Child rectangle 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. + its blocked storage omits holes. Validate the actual adjacent-pair support over the whole union; + no rectangle-gap policy is imposed. """ child = FixedPatchSet2D.from_config( _middle_config(cfg, parent), child_rectangles_in_parent) @@ -373,6 +410,17 @@ def _nested_kind(parent, child): return _readonly(composed, np.uint8) +def _circle_in_coverage(covered, cx, cy, radius): + """Exact disk versus uncovered parent squares; no containing-rectangle rule.""" + ny, nx = covered.shape + if cx-radius < 0 or cy-radius < 0 or cx+radius > nx or cy+radius > ny: + return False + y, x = np.nonzero(~covered) + dx = np.maximum(np.maximum(x-cx, cx-(x+1)), 0.) + dy = np.maximum(np.maximum(y-cy, cy-(y+1)), 0.) + return not np.any(dx*dx + dy*dy < radius*radius) + + def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=(), parent_origin=(0., 0.), parent_spacing=1.): """Project fixed circles and reject boundary/transfer overlap before CUDA. @@ -386,7 +434,7 @@ def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=(), Original objects are not retained or mutated. No runtime geometry updates. ``patch`` may be a FixedPatchSet2D: objects are then projected once onto the level-wide fine grid (fine_origin (0, 0)); each circle must lie inside - one patch, body ids are the positions in ``objects`` (never patch ids), and + the union, body ids are the positions in ``objects`` (never patch ids), and jets and the (receiver, dir) writer checks run once over the whole level. ``jets`` use the same coarse cell-edge frame as the objects; their selection is checked against true fine hits here. Jet amplitudes need no @@ -416,7 +464,6 @@ def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=(), flags = np.full(nx * ny, FLUID, np.uint16) active = patch.fine_kind.ravel() == ACTIVE local, links = [], [] - rects = patch.bounds if isinstance(patch, FixedPatchSet2D) else (patch.bounds,) covered = patch.covered dx = dt = .5 * parent_spacing fine_origin = origin + parent_spacing * np.asarray(patch.fine_origin) @@ -430,9 +477,8 @@ def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=(), px, py = (np.asarray((cx, cy)) - origin) / parent_spacing pr = radius / parent_spacing if (not np.isfinite([cx, cy, radius]).all() or radius <= 0 - or not any(x0 <= px-pr and px+pr <= x1 and y0 <= py-pr - and py+pr <= y1 for x0, y0, x1, y1 in rects)): - raise ValueError("Entire fine circle must lie inside one physical patch") + or not _circle_in_coverage(covered, px, py, pr)): + raise ValueError("Entire fine circle must lie inside the physical patch union") center = tuple((v-o)/dx-.5 for v, o in zip((cx, cy), fine_origin)) spec = obj.surface_actuation if spec is not None: diff --git a/src/CelerisLab/lbm/refinement_stepper.py b/src/CelerisLab/lbm/refinement_stepper.py index 44e95a1..7e05e01 100644 --- a/src/CelerisLab/lbm/refinement_stepper.py +++ b/src/CelerisLab/lbm/refinement_stepper.py @@ -170,13 +170,13 @@ class FixedPatchStepper: re-uploads the per-link velocity column only. Several patches: pass ``bounds`` as a sequence of (x0, y0, x1, y1) - rectangles (FixedPatchSet2D, >= 4 coarse cells apart). This requires + rectangles (FixedPatchSet2D, dependency-validated union). This requires fine_block and base_block; fine_objects/jets are projected once onto the - level-wide fine grid (each circle inside one patch, body ids = object + level-wide fine grid (each circle inside the union, body ids = object positions, jet and writer checks level-wide), and uses one level-wide fine grid (2*ny, 2*nx): public fine arrays have that shape, NaN outside the patches' ACTIVE cells. The step sequence is the - blocked one above, with every patch's transfers in one launch each. + blocked one above, with unique union transfers in one launch each. Three-level separated unions: ``child_rectangles_in_parent=[...]`` is mutually exclusive with scalar ``child`` and accepts either parent form. @@ -507,8 +507,8 @@ 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``: - stationary circles pinned in unchanged complete L2 rectangles. - Separated fluid-only rectangles may be added, moved or removed. + stationary circles with identical physical boundary and BFL support. + Other fluid coverage may change, even in the same connected union. Current L2 Jets (including updated amplitudes) are preserved. Local projected boundary objects are not valid inputs. The lock and stream synchronization certify @@ -547,10 +547,9 @@ class FixedPatchStepper: 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) - if (old.pinned_rectangles != new.pinned_rectangles - or not same_physical_boundary(old.boundary, new.boundary) + 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 body rectangles and boundary") + raise ValueError("remap requires the identical pinned physical boundary and BFL support") base_block, fine_block, child_block = self._level_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 3e10973..f39d7d2 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 -stationary circles pinned inside unchanged L2 rectangles. The caller +stationary circles with unchanged physical boundary/BFL support. 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 @@ -20,13 +20,14 @@ from ..body.jets import surface_link_velocity @dataclass(frozen=True, eq=False) class RemapTopology: - """Validated fixed parent and separated L2 set; use ``from_config``. + """Validated fixed parent and L2 union; use ``from_config``. 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. + ``pinned_rectangles`` is legacy metadata for single-rectangle containment. + Eligibility pins physical boundary/stencil data, not entire rectangles or + connected fluid regions; each topology independently validates transfer closure. """ parent: object child: FixedPatchSet2D @@ -59,9 +60,11 @@ class RemapTopology: 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])) + match = 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]), None) + if match is not None: + containing.add(match) pinned = tuple(sorted(containing)) velocity = boundary_velocity(cfg, child, boundary, child_jets) if velocity is not None: @@ -167,10 +170,9 @@ 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_rectangles != new.pinned_rectangles - or not same_physical_boundary(old.boundary, new.boundary) + 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 body rectangles and physical boundary") + raise ValueError("remap requires the identical pinned physical boundary and BFL support") 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_hierarchy.py b/tests/unit/test_hierarchy.py index 843efca..930d756 100644 --- a/tests/unit/test_hierarchy.py +++ b/tests/unit/test_hierarchy.py @@ -50,11 +50,11 @@ class FourLevels(unittest.TestCase): cases = [([L1, [(10, 20, 30, 44)]], "inherited ACTIVE"), ([L1, L2[:1] + [(54, 54, 58, 58)]], "inherited ACTIVE"), ([L1, L2, [(40, 40, 60, 60)]], "donors/receivers"), - ([L1, L2, [(41, 41, 60, 60)]], "donors/receivers"), - ([L1, L2, [(56, 56, 64, 64), (66, 56, 72, 64)]], "closer than")] + ([L1, L2, [(41, 41, 60, 60)]], "donors/receivers")] for levels, reason in cases: with self.subTest(levels=levels), self.assertRaisesRegex(ValueError, reason): build_hierarchy(CFG, levels) + build_hierarchy(CFG, [L1, L2, [(56, 56, 64, 64), (66, 56, 72, 64)]]) # gap 2 closes build_hierarchy(CFG, [L1, L2, [(42, 42, 60, 60)]]) # minimal legal margin diff --git a/tests/unit/test_multiple_pinned_remap.py b/tests/unit/test_multiple_pinned_remap.py index 0ee275f..a68d771 100644 --- a/tests/unit/test_multiple_pinned_remap.py +++ b/tests/unit/test_multiple_pinned_remap.py @@ -81,14 +81,9 @@ class MultiplePinnedRemap(unittest.TestCase): 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): + def test_direct_helper_rejects_changed_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'): @@ -138,6 +133,7 @@ class MultiplePinnedRemap(unittest.TestCase): changed = list(PINS+FLUID) x0, y0, x1, y1 = changed[index] changed[index] = (x0-4, y0, x1, y1) + changed[index] = (x0+8, y0, x1, y1) # removes actual body support cases.append((changed, bodies)) cases.append((tuple(r for i, r in enumerate(PINS+FLUID) if i != index), bodies)) start = len(allocations) diff --git a/tests/unit/test_pinned_body_remap.py b/tests/unit/test_pinned_body_remap.py index c105b75..24c05f9 100644 --- a/tests/unit/test_pinned_body_remap.py +++ b/tests/unit/test_pinned_body_remap.py @@ -83,7 +83,7 @@ class PinnedBodyRemap(unittest.TestCase): run,module,allocations,storage,objects=ready() start=len(allocations) cases=[(B,()),(B,run.child_boundary.objects),(B,(body(center=(12.1,14.)),)), - (((20,16,36,40),B[1]),objects),(((16,20,32,40),B[1]),objects), + (((24,16,36,40),B[1]),objects),(((16,28,32,40),B[1]),objects), ((B[1],),objects)] module.cuda.memcpy_dtoh.reset_mock() for rects,objs in cases: @@ -124,7 +124,9 @@ class PinnedBodyRemap(unittest.TestCase): for a in allocations[len(saved):]: a.free.assert_called_once() run.close() - def test_direct_helper_rejects_changed_pinned_rectangle(self): + def test_direct_helper_allows_changed_rectangle_with_same_boundary(self): old=mesh(A,(body(),)) new=mesh(((12,16,32,40),A[1]),(body(),)) - with self.assertRaisesRegex(ValueError,'pinned'): remap(old,new,fields(old)) + out = remap(old,new,fields(old)) + self.assertTrue(same_physical_boundary(old.boundary,new.boundary)) + self.assertIs(out.topology,new) diff --git a/tests/unit/test_refinement.py b/tests/unit/test_refinement.py index 5dd80c2..48c84ee 100644 --- a/tests/unit/test_refinement.py +++ b/tests/unit/test_refinement.py @@ -181,14 +181,14 @@ class FixedPatchSetTests(unittest.TestCase): ly, lx = p.fine_shape np.testing.assert_array_equal(s.fine_kind[y0:y0 + ly, x0:x0 + lx], p.fine_kind) gy, gx = np.divmod(p.ghost_fine.astype(np.int64), lx) - np.testing.assert_array_equal(s.ghost_fine[cells:cells + len(gy)], - (gy + y0) * fx + gx + x0) - np.testing.assert_array_equal(s.ghost_coarse[cells:cells + len(gy)], p.ghost_coarse) + ghosts = dict(zip(map(int, s.ghost_fine), map(int, s.ghost_coarse))) + for f, c in zip((gy + y0) * fx + gx + x0, p.ghost_coarse): + self.assertEqual(ghosts[int(f)], int(c)) q, c = np.divmod(p.coalesce_fine.astype(np.int64), ly * lx) cy, cx = np.divmod(c, lx) - np.testing.assert_array_equal(s.coalesce_fine[rows:rows + len(q)], - q * fy * fx + (cy + y0) * fx + cx + x0) - np.testing.assert_array_equal(s.coalesce_coarse[rows:rows + len(q)], p.coalesce_coarse) + rows_by_target = dict(zip(map(int, s.coalesce_coarse), s.coalesce_fine)) + for target, children in zip(p.coalesce_coarse, q * fy * fx + (cy + y0) * fx + cx + x0): + np.testing.assert_array_equal(rows_by_target[int(target)], children) rows, cells = rows + len(q), cells + len(gy) self.assertEqual((rows, cells), (len(s.coalesce_coarse), len(s.ghost_fine))) self.assertEqual(int((s.fine_kind == ACTIVE).sum()), 4 * int(s.covered.sum())) @@ -207,10 +207,12 @@ class FixedPatchSetTests(unittest.TestCase): np.testing.assert_array_equal(s.coalesce_coarse, p.coalesce_coarse) np.testing.assert_array_equal(s.ghost_coarse, p.ghost_coarse) - def test_gap_and_overlap_rejected(self): - for gap in (0, 2, 3): - with self.subTest(gap=gap), self.assertRaisesRegex(ValueError, "closer"): - self.build(((4, 4, 8, 12), (8 + gap, 4, 12 + gap, 12))) + def test_gap_and_overlap_use_dependency_checks(self): + for gap in (-2, 0, 1, 2, 3): + with self.subTest(gap=gap): + pair = self.build(((4, 4, 8, 12), (8 + gap, 4, 12 + gap, 12))) + refinement._transfer_footprint(pair) + refinement._check_parent_pull_closure(pair) with self.assertRaises(ValueError): self.build(()) with self.assertRaises(ValueError): # FixedPatch2D margin still applies @@ -226,7 +228,7 @@ class FixedPatchSetTests(unittest.TestCase): support[y0:y0 + p.fine_shape[0], x0:x0 + p.fine_shape[1]] = True np.testing.assert_array_equal(layout.owned, support) self.assertTrue(layout.owned.ravel()[s.ghost_fine.astype(np.intp)].all()) - with self.assertRaisesRegex(ValueError, "aligned"): + with self.assertRaisesRegex(ValueError, "divide"): self.build(self.FACE).fine_layout((5, 4)) for block in ((0, 4), (4, 0), (-4, 4), (4.0, 4), (True, 4), (4,), (4, 4, 4), (np.bool_(1), 4)): with self.subTest(block=block), self.assertRaisesRegex(ValueError, "positive integers"): diff --git a/tests/unit/test_refinement_three_level_multipatch.py b/tests/unit/test_refinement_three_level_multipatch.py index 2268c5a..a66f17b 100644 --- a/tests/unit/test_refinement_three_level_multipatch.py +++ b/tests/unit/test_refinement_three_level_multipatch.py @@ -217,7 +217,7 @@ class ConstructorAndSchedule(unittest.TestCase): 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=[(16,16,32,24),(16,27,32,55)]), 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), @@ -226,8 +226,8 @@ class ConstructorAndSchedule(unittest.TestCase): 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) + with self.assertRaisesRegex(ValueError,'margin'): + host.FixedPatchStepper(CFG,[(1,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) diff --git a/tests/unit/test_runtime_child_remap.py b/tests/unit/test_runtime_child_remap.py index 384ebfc..4bd3492 100644 --- a/tests/unit/test_runtime_child_remap.py +++ b/tests/unit/test_runtime_child_remap.py @@ -67,7 +67,7 @@ class RuntimeChildRemap(unittest.TestCase): run._status = status with self.assertRaises(RuntimeError): run.remap_child(B) run._status = 'ready' - for bad in ([], [(40,16,48,40)], [(16,16,32,24),(16,27,32,40)]): + for bad in ([], [(40,16,48,40)], [(16,16,32,24),(16,27,32,55)]): with self.assertRaises(ValueError): run.remap_child(bad) hosts, pairs = run._boundary_hosts, run.hierarchy.pairs for attribute, value in (('_wall_aware', True), diff --git a/tests/unit/test_synchronized_remap.py b/tests/unit/test_synchronized_remap.py index 30a4ca5..826600f 100644 --- a/tests/unit/test_synchronized_remap.py +++ b/tests/unit/test_synchronized_remap.py @@ -175,8 +175,7 @@ class SynchronizedRemapTests(unittest.TestCase): self.assertFalse(np.array_equal(self.g[2][:, lost], returned.g[2][:, lost])) def test_invalid_geometry_counts_representation_frame(self): - for rectangles in ([], [(2, 2, 10, 10)], [(16, 16, 32, 40)] * 2, - [(16, 16, 33.5, 40)], [(16, 16, 32, 40), (33, 16, 36, 40)]): + for rectangles in ([], [(2, 2, 10, 10)], [(16, 16, 33.5, 40)]): with self.assertRaises(ValueError): topology(rectangles) parent = self.old.parent diff --git a/tests/unit/test_union_body_support.py b/tests/unit/test_union_body_support.py new file mode 100644 index 0000000..a58af07 --- /dev/null +++ b/tests/unit/test_union_body_support.py @@ -0,0 +1,44 @@ +"""Union seam bodies and actual support pinning; CPU host-only.""" +import unittest +import numpy as np +from CelerisLab.config import LBMConfig +from CelerisLab.lbm.refinement import FixedPatchSet2D,prepare_fine_boundary,_circle_in_coverage +from CelerisLab.lbm.synchronized_remap import same_physical_boundary +from tests.unit import test_refinement as refinement_tests +from tests.unit.test_pinned_body_remap import mesh,body,ready,A +from tests.unit.test_synchronized_remap import fields,remap,integrals + +class UnionBodySupport(unittest.TestCase): + def test_seam_circle_matches_single_rectangle(self): + cfg=LBMConfig(nx=32,ny=16,viscosity=.1) + obj=refinement_tests.PatchSetBoundaryTests.body(16.,8.) + def b(rects):return prepare_fine_boundary(cfg,FixedPatchSet2D.from_config(cfg,rects),[obj]) + self.assertTrue(same_physical_boundary(b([(4,2,28,14)]),b([(4,2,16,14),(16,2,28,14)]))) + with self.assertRaisesRegex(ValueError,'union'): + b([(4,2,16,14),(17,2,28,14)]) + + def test_continuous_disk_not_only_sampled_solids(self): + mask=np.ones((8,8),bool);mask[3,3]=False + self.assertFalse(_circle_in_coverage(mask,3.99,3.99,.02)) + self.assertTrue(_circle_in_coverage(mask,4.1,4.1,.1)) + + def test_repartition_and_connected_fluid_changes_do_not_pin_rectangles(self): + old=mesh(A,(body(),));initial=fields(old) + for rects in (((16,16,24,40),(24,16,32,40),A[1]),((12,16,32,40),A[1]),((20,16,36,40),A[1])): + new=mesh(rects,(body(),));self.assertTrue(same_physical_boundary(old.boundary,new.boundary)) + out=remap(old,new,initial) + np.testing.assert_allclose(integrals(old,initial),integrals(new,out.g),rtol=2e-13,atol=2e-12) + for a,b,om,nm in zip(initial,out.g,old.leaf_masks,new.leaf_masks): + np.testing.assert_array_equal(a[:,om&nm],b[:,om&nm]) + + def test_runtime_seam_repartition_preserves_boundary_and_live_jets(self): + run,module,allocations,storage,objects=ready() + try: + run.set_jet_amplitude(0,.017,level=2) + clone=run.remap_child(((16,16,24,40),(24,16,32,40),A[1]),child_objects=objects) + try: + self.assertEqual(clone._boundaries[2].jets[0].amplitude,.017) + self.assertTrue(same_physical_boundary(run.child_boundary,clone.child_boundary)) + for a,b in zip(run.download(),clone.download()):np.testing.assert_array_equal(a,b) + finally:clone.close() + finally:run.close() diff --git a/tests/unit/test_union_refinement.py b/tests/unit/test_union_refinement.py new file mode 100644 index 0000000..c38bdad --- /dev/null +++ b/tests/unit/test_union_refinement.py @@ -0,0 +1,161 @@ +"""CPU tests: union tables vs current FixedPatchSet2D; numpy + CelerisLab on PYTHONPATH. +Optional numerics: ADAPTER_DIR= [ADAPTER_FILE=adapter.py or +multilevel_numeric_adapter.py (default: whichever exists)].""" +from dataclasses import replace +import importlib.util +import os +import unittest + +import numpy as np +from CelerisLab.config import LBMConfig +from CelerisLab.lbm.refinement import (ACTIVE, OUTER_GHOST, FixedPatchSet2D, _middle_config, + _nested_kind, _transfer_footprint) +from CelerisLab.lbm.refinement_tables import blocked_kind + +from CelerisLab.lbm.refinement import (FixedPatchSet2D as UnionPatchSet2D, + _check_parent_pull_closure as check_parent_pull_closure, _dilate as dilate) + +def table_sets(pair): + return (pair.fine_kind.tobytes(), pair.covered.tobytes(), + frozenset(zip(map(int,pair.ghost_fine),map(int,pair.ghost_coarse))), + frozenset((int(t),tuple(map(int,c))) for t,c in zip(pair.coalesce_coarse,pair.coalesce_fine))) + +CFG = LBMConfig(nx=30, ny=30, viscosity=.1) +A, B, U = (8, 8, 16, 20), (16, 8, 24, 20), (8, 8, 24, 20) + + +def union(*rects, cfg=CFG): + return UnionPatchSet2D.from_config(cfg, rects) + + +class Equivalence(unittest.TestCase): + def test_abutting_boxes_equal_their_union_rectangle(self): + reference = table_sets(FixedPatchSet2D.from_config(CFG, [U])) + for rects in ((A, B), (B, A), (U,)): + with self.subTest(rects=rects): + self.assertEqual(table_sets(union(*rects)), reference) + + def test_vertical_abutting(self): + self.assertEqual(table_sets(union((8, 8, 20, 14), (8, 14, 20, 22))), + table_sets(FixedPatchSet2D.from_config(CFG, [(8, 8, 20, 22)]))) + + def test_separated_rectangles_reproduce_current_tables(self): + rects = [(4, 4, 10, 10), (16, 16, 26, 26)] # gap 6 >= MIN_PATCH_GAP + self.assertEqual(table_sets(union(*rects)), + table_sets(FixedPatchSet2D.from_config(CFG, rects))) + + def test_layout_tiles_are_free_of_rectangle_alignment(self): + pair = union(A, B) + support = dilate(pair.fine_kind == ACTIVE, 4) + reference = FixedPatchSet2D.from_config(CFG, [U]).fine_layout((4, 4)) + np.testing.assert_array_equal(pair.fine_layout((4, 4)).owned, reference.owned) + layout = pair.fine_layout((6, 10)) # not aligned to the rectangles + self.assertTrue(layout.owned[support].all()) + kind = blocked_kind(layout, pair.fine_kind) + owned = layout.owned + np.testing.assert_array_equal(kind[layout.owner[owned]], pair.fine_kind[owned]) + + +class NewlyLegal(unittest.TestCase): + def check(self, pair): + _transfer_footprint(pair) + check_parent_pull_closure(pair) + ghosts = pair.ghost_fine.astype(np.intp) + self.assertFalse((pair.fine_kind.ravel()[ghosts] == ACTIVE).any()) + self.assertFalse(pair.covered.ravel()[pair.ghost_coarse.astype(np.intp)].any()) + + def test_small_gaps_and_l_shape(self): + for rects in (((8, 8, 16, 20), (17, 8, 24, 20)), # gap 1 + ((8, 8, 16, 20), (18, 8, 24, 20)), # gap 2 + ((8, 8, 20, 14), (8, 14, 14, 22))): # L-shape + with self.subTest(rects=rects): + self.check(union(*rects)) + + def test_corner_contact_only_by_dependency_checks(self): + # Policy: no geometric ban; accepted iff the direction-specific checks pass. + try: + pair = union((8, 8, 16, 16), (16, 16, 24, 24)) + except ValueError as error: + self.assertRegex(str(error), "outer ghost|allocation|INNER|surviving|coalesce") + else: + self.check(pair) + + +class Negatives(unittest.TestCase): + def test_inputs_are_validated_before_conversion(self): + for bad in ((8.5, 8, 16, 20), ("8", 8, 16, 20), (True, 8, 16, 20), "abcd", (8, 8, 16)): + with self.subTest(bad=bad), self.assertRaises(ValueError): + union(bad, B) + self.assertEqual(union(tuple(np.int64(v) for v in A), B).bounds, (A, B)) + pair = union(A, B) + for block in ((4.0, 4), (True, 4), (0, 4), (4,), ("4", 4)): + with self.subTest(block=block), self.assertRaises(ValueError): + pair.fine_layout(block) + + def test_missing_donor_at_domain_margin(self): + with self.assertRaisesRegex(ValueError, "margin"): + union((1, 8, 10, 20), (10, 8, 20, 20)) + + def test_missing_or_duplicate_coalesce_row(self): + pair = union(A, B) + with self.assertRaisesRegex(ValueError, "lacks its coalesce row"): + check_parent_pull_closure(replace(pair, coalesce_coarse=pair.coalesce_coarse[1:], + coalesce_fine=pair.coalesce_fine[1:])) + doubled = replace(pair, coalesce_coarse=np.r_[pair.coalesce_coarse, pair.coalesce_coarse[:1]], + coalesce_fine=np.r_[pair.coalesce_fine, pair.coalesce_fine[:1]]) + with self.assertRaisesRegex(ValueError, "duplicate"): + check_parent_pull_closure(doubled) + + def test_direction_specific_dependency_on_outer_cell(self): + pair = union(A, B) + active = pair.fine_kind == ACTIVE + kind = pair.fine_kind.copy() + kind[dilate(active, 2) & ~dilate(active, 1)] = OUTER_GHOST # thin the INNER band + with self.assertRaisesRegex(ValueError, "outer ghost"): + _transfer_footprint(replace(pair, fine_kind=kind)) + + def test_nested_union_child_donors_must_be_parent_active(self): + parent = FixedPatchSet2D.from_config(CFG, [(6, 6, 24, 24)]) # L1 ACTIVE 12..48 + middle = _middle_config(CFG, parent) + with self.assertRaisesRegex(ValueError, "donors/receivers"): + _nested_kind(parent, UnionPatchSet2D.from_config(middle, [(12, 20, 30, 40), (30, 20, 40, 40)])) + child = UnionPatchSet2D.from_config(middle, [(16, 20, 30, 40), (30, 20, 40, 40)]) + reference = FixedPatchSet2D.from_config(middle, [(16, 20, 40, 40)]) + self.assertEqual(table_sets(child), table_sets(reference)) + np.testing.assert_array_equal(_nested_kind(parent, child), _nested_kind(parent, reference)) + + +@unittest.skipUnless(os.environ.get("ADAPTER_DIR"), "set ADAPTER_DIR for the CPU D2Q9 adapter run") +class AdapterNumerics(unittest.TestCase): + def test_eight_steps_bitwise(self): + folder = os.environ["ADAPTER_DIR"] + name = os.environ.get("ADAPTER_FILE") or next( + n for n in ("multilevel_numeric_adapter.py", "adapter.py") + if os.path.exists(os.path.join(folder, n))) + path = os.path.join(folder, name) + spec = importlib.util.spec_from_file_location("multilevel_numeric_adapter", path) + adapter = importlib.util.module_from_spec(spec) + spec.loader.exec_module(adapter) + results = [] + for pair in (union(A, B), FixedPatchSet2D.from_config(CFG, [U])): + kinds = [np.where(pair.covered, OUTER_GHOST, ACTIVE).astype(np.uint8), pair.fine_kind] + taus = [pair.tau_coarse, pair.tau_fine] + initial = [] + for level, kind in enumerate(kinds): + ny, nx = kind.shape + y = (np.arange(ny) + .5) / 2 ** level + ux = 1e-3 * np.sin(2 * np.pi * y / 30)[:, None] * np.ones((1, nx)) + cx = np.array([0, 1, -1, 0, 0, 1, -1, 1, -1], float)[:, None, None] + w = np.array([4/9] + [1/9]*4 + [1/36]*4)[:, None, None] + g = w * (1 + 3*cx*ux + 4.5*(cx*ux)**2 - 1.5*ux*ux) + g[:, kind != ACTIVE] = np.nan + initial.append(g) + run = adapter.make_scheduler(kinds, taus, [pair], initial) + run.step(8) + results.append([level.current.copy() for level in run.levels]) + for a, b in zip(*results): + np.testing.assert_array_equal(a, b) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/validation/multigrid/test_refinement_multipatch.py b/tests/validation/multigrid/test_refinement_multipatch.py index 4cf2815..3fe05ac 100644 --- a/tests/validation/multigrid/test_refinement_multipatch.py +++ b/tests/validation/multigrid/test_refinement_multipatch.py @@ -50,8 +50,8 @@ class MultipatchCPU(unittest.TestCase): dict(fine_block=(5, 4), base_block=(4, 4))): with self.subTest(kw=kw), self.assertRaises(ValueError): FixedPatchStepper(cfg, ((4, 4, 8, 12), (14, 4, 18, 12)), context=context, **kw) - with self.assertRaisesRegex(ValueError, "closer"): - FixedPatchStepper(cfg, ((4, 4, 8, 12), (10, 4, 14, 12)), context=context, + with self.assertRaisesRegex(ValueError, "margin"): + FixedPatchStepper(cfg, ((1, 4, 8, 12), (10, 4, 14, 12)), context=context, fine_block=(8, 4), base_block=(4, 4)) from CelerisLab.body.jets import Jet, Parallel, Tube cfg32 = LBMConfig(nx=32, ny=16, viscosity=.1)