diff --git a/README.md b/README.md index 23fb20a..3e9caf8 100644 --- a/README.md +++ b/README.md @@ -135,13 +135,16 @@ Scope and limits: divide the WHOLE fine allocation, i.e. `2*(y1-y0)+8` by `2*(x1-x0)+8` including the four-cell ghost band. It gives results identical to dense storage and is infrastructure for future block AMR, not a speedup. -- Several separated patches (body-free): pass a sequence of `(x0, y0, x1, y1)` +- Several separated patches: pass a sequence of `(x0, y0, x1, y1)` rectangles as `bounds`, at least 4 coarse cells apart, with both `fine_block` and `base_block`. The fine level is then one level-wide grid `(2*ny, 2*nx)`: `run.patch.fine_shape` gives it, `initialize`/`download` use it, and only the patch windows (active region plus 4-cell ghost band, origin `(2*y0-4, 2*x0-4)`) are stored; each window must align with `fine_block`. Order of rectangles does - not matter. No bodies/Jets, wall-touching patches or nesting yet. + not matter. `fine_objects`/`jets` work as for one patch: each circle (same + coarse cell-edge frame) must lie wholly inside one patch's active region with + its boundary stencil clear of transfer cells; a jet may not select links of two + different bodies, even across patches. No wall-touching patches or nesting yet. ```python cfg = LBMConfig(nx=32, ny=16, viscosity=.1) @@ -175,7 +178,7 @@ Scope and limits: coarse_g, l1_g, l2_g = run.download() ``` - Not available: per-link surface exchange/force readout in blocked mode, - checkpoints, moving bodies, bodies/Jets with several patches or three levels, + checkpoints, moving bodies, bodies/Jets with three levels, several patches per level with nesting, dynamic refinement, `Simulation` integration. Validation so far is storage/operator equivalence and short runs, not accuracy, convergence or long-time stability. diff --git a/src/CelerisLab/lbm/refinement.py b/src/CelerisLab/lbm/refinement.py index 43ea6c8..5a769a0 100644 --- a/src/CelerisLab/lbm/refinement.py +++ b/src/CelerisLab/lbm/refinement.py @@ -227,6 +227,11 @@ class FixedPatchSet2D: def bounds(self) -> tuple: return tuple(p.bounds for p in self.patches) + @property + def fine_origin(self) -> tuple[float, float]: + """Lower edge of the level-wide fine grid in coarse physical units.""" + return 0.0, 0.0 + @property def covered(self) -> np.ndarray: mask = np.zeros(self.coarse_shape, bool) @@ -351,6 +356,10 @@ def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=()): SurfaceActuation retains its lattice-velocity meaning: input coarse units are converted into independent fine-local objects with dx=dt=.5. 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 + 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 rescaling (lattice velocity scales by dt/dx = 1 under acoustic 2:1). @@ -375,7 +384,8 @@ 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 = [], [] - x0, y0, x1, y1 = patch.bounds + rects = patch.bounds if isinstance(patch, FixedPatchSet2D) else (patch.bounds,) + covered = patch.covered dx = dt = .5 for bid, obj in enumerate(objects): if (type(obj.geometry) is not CircleGeometry or obj.is_sensor or obj.is_force_region @@ -385,8 +395,9 @@ def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=()): cx, cy = obj.geometry.center radius = obj.geometry._radius if (not np.isfinite([cx, cy, radius]).all() or radius <= 0 - or cx-radius < x0 or cx+radius > x1 or cy-radius < y0 or cy+radius > y1): - raise ValueError("Entire fine circle must lie inside the physical patch") + or not any(x0 <= cx-radius and cx+radius <= x1 and y0 <= cy-radius + and cy+radius <= y1 for x0, y0, x1, y1 in rects)): + raise ValueError("Entire fine circle must lie inside one physical patch") center = tuple((v-o)/dx-.5 for v, o in zip((cx, cy), patch.fine_origin)) spec = obj.surface_actuation if spec is not None: @@ -446,7 +457,7 @@ def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=()): raise ValueError("Boundary stencil overlaps the two-hop transfer footprint") for k in coarse: y, x = divmod(k, cfg.nx) - if x0 <= x < x1 and y0 <= y < y1: + if covered[y, x]: raise ValueError("Transfer coarse donor/target is not an uncovered bulk cell") # Rows 0/ny-1 are never computed and rows 1/ny-2 are reconstructed after # coalescence; rows 2/ny-3 are ordinary bulk cells only read as donors. diff --git a/src/CelerisLab/lbm/refinement_stepper.py b/src/CelerisLab/lbm/refinement_stepper.py index ead846d..4385bac 100644 --- a/src/CelerisLab/lbm/refinement_stepper.py +++ b/src/CelerisLab/lbm/refinement_stepper.py @@ -134,7 +134,9 @@ class FixedPatchStepper: Several patches: pass ``bounds`` as a sequence of (x0, y0, x1, y1) rectangles (FixedPatchSet2D, >= 4 coarse cells apart). This requires - fine_block and base_block, is body-free (no fine_objects/jets yet), and + 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 + 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. @@ -150,8 +152,6 @@ class FixedPatchStepper: if multipatch: if fine_block is None or base_block is None: raise ValueError("several patches require fine_block and base_block") - if tuple(fine_objects) or tuple(jets): - raise ValueError("several patches do not support fine_objects or jets yet") self.patch = FixedPatchSet2D.from_config(cfg, bounds) else: self.patch = FixedPatch2D.from_config(cfg, bounds) diff --git a/tests/unit/test_refinement.py b/tests/unit/test_refinement.py index 0f2df56..a993e35 100644 --- a/tests/unit/test_refinement.py +++ b/tests/unit/test_refinement.py @@ -233,6 +233,56 @@ class FixedPatchSetTests(unittest.TestCase): self.build(self.FACE).fine_layout(block) +class PatchSetBoundaryTests(unittest.TestCase): + """fine_objects/jets on a level-wide patch set (16x32, patches 4 cells apart).""" + PATCHES = ((2, 2, 14, 14), (18, 2, 30, 14)) + + @staticmethod + def body(x, y, obj_id=0): + from CelerisLab.body.geometry.circle import CircleGeometry + from CelerisLab.body.objects import SimObject + return SimObject(obj_id, geometry=CircleGeometry(x, y, 2.15), center=(x, y), radius=2.15) + + def prepare(self, nx, rects, objects, jets=()): + cfg = LBMConfig(nx=nx, ny=16, viscosity=.1) + s = refinement.FixedPatchSet2D.from_config(cfg, rects) + return s, refinement.prepare_fine_boundary(cfg, s, objects, jets=jets) + + def test_two_bodies_are_the_one_body_case_translated(self): + _, one = self.prepare(16, self.PATCHES[:1], [self.body(8., 8.)]) + _, two = self.prepare(32, self.PATCHES, [self.body(8., 8., 7), self.body(24., 8., 7)]) + f1, f2 = one.flags.reshape(32, 32), two.flags.reshape(32, 64) + np.testing.assert_array_equal(f2[:, :32], f1) + np.testing.assert_array_equal(f2[:, 32:], f1) + n = len(one.curved_host[0]) + self.assertEqual(len(two.curved_host[0]), 2 * n) + body_id = two.curved_host[3] + self.assertEqual(sorted(set(body_id.tolist())), [0, 1]) # positions, not obj_id 7 + # Translation must reproduce identical link geometry (prerequisite of the + # bitwise GPU tiling gate): dirs, q and fallback equal, receivers shifted. + fidx, dirs, q, fallback = (two.curved_host[i] for i in (0, 1, 2, 7)) + for column in (dirs, q, fallback): + np.testing.assert_array_equal(column[:n], column[n:]) + y, x = np.divmod(fidx.astype(np.int64), 64) + np.testing.assert_array_equal(y[:n], y[n:]) + np.testing.assert_array_equal(x[:n] + 32, x[n:]) + + def test_body_outside_a_single_patch_rejected(self): + for x in (16., 13.): # in the gap / crossing the patch edge + with self.subTest(x=x), self.assertRaises(ValueError): + self.prepare(32, self.PATCHES, [self.body(x, 8.)]) + + def test_one_jet_over_two_bodies_rejected_even_at_zero_amplitude(self): + from CelerisLab.body.jets import Jet, Parallel, Tube + bodies = [self.body(8., 8.), self.body(24., 8.)] + span = Jet(Tube((5., 8.), 0., 22., .6), Parallel(0.), amplitude=0.) + with self.assertRaisesRegex(ValueError, "one jet"): + self.prepare(32, self.PATCHES, bodies, jets=[span]) + own = [Jet(Tube((5., 8.), 0., 6., .6), Parallel(0.), amplitude=.002), + Jet(Tube((21., 8.), 0., 6., .6), Parallel(0.), amplitude=.002)] + self.prepare(32, self.PATCHES, bodies, jets=own) # one jet per body: accepted + + class NestPatchTests(unittest.TestCase): """Three-level single chain: L2 bounds in L1 window cells (probe fixture).""" def setUp(self): diff --git a/tests/validation/multigrid/test_refinement_multipatch.py b/tests/validation/multigrid/test_refinement_multipatch.py index 3891ad3..4cf2815 100644 --- a/tests/validation/multigrid/test_refinement_multipatch.py +++ b/tests/validation/multigrid/test_refinement_multipatch.py @@ -53,6 +53,13 @@ class MultipatchCPU(unittest.TestCase): with self.assertRaisesRegex(ValueError, "closer"): FixedPatchStepper(cfg, ((4, 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) + span = Jet(Tube((5., 8.), 0., 22., .6), Parallel(0.), amplitude=0.) + with self.assertRaisesRegex(ValueError, "one jet"): # same obj_id cannot bypass it + FixedPatchStepper(cfg32, ((2, 2, 14, 14), (18, 2, 30, 14)), context=context, + fine_block=(4, 4), base_block=(4, 4), + fine_objects=[body(8., 8.), body(24., 8.)], jets=[span]) class MultipatchGPU(unittest.TestCase): @@ -144,5 +151,76 @@ class MultipatchGPU(unittest.TestCase): print(f"face gap4 step {step}: weighted leaf mass {mass:.9g} (diagnostic)", flush=True) +def body(x, y, obj_id=0): + from CelerisLab.body.geometry.circle import CircleGeometry + from CelerisLab.body.objects import SimObject + return SimObject(obj_id, geometry=CircleGeometry(x, y, 2.15), center=(x, y), radius=2.15) + + +class MultipatchBodiesGPU(unittest.TestCase): + """fine_objects/jets on several patches (frozen gates, bitwise).""" + PATCHES = ((2, 2, 14, 14), (18, 2, 30, 14)) + STEPS = (0, 1, 2, 8) + + def run_steps(self, context, cfg, bounds, coarse, fine, switch=None, **kw): + out = [] + with FixedPatchStepper(cfg, bounds, context=context, **kw) as backend: + backend.initialize(coarse, fine) + for step in self.STEPS: + if switch is not None and step == 8: + for index, amplitude in switch: + backend.set_jet_amplitude(index, amplitude) + backend.step(step - backend.step_count) + out.append(backend.download() + raw(backend)) + return out + + def test_one_body_matches_single_patch_path(self): + from CelerisLab.cuda.context import CudaContext + with CudaContext(0) as context: + for mode in ("periodic", "slot_reservoir"): + cfg = config(mode) + p = FixedPatch2D.from_config(cfg, BOUNDS) + s = FixedPatchSet2D.from_config(cfg, (BOUNDS,)) + coarse, fine = state(NY, NX, 1), state(*p.fine_shape, 2) + kw = dict(coarse_boundary=mode, fine_block=(8, 4), base_block=(3, 4), + fine_objects=[circle()]) + old = self.run_steps(context, cfg, BOUNDS, coarse, fine, **kw) + new = self.run_steps(context, cfg, (BOUNDS,), coarse, level_fine(s, [fine]), **kw) + for step, a, b in zip(self.STEPS, old, new): + label = f"{mode} step {step}" + exact(self, b[0], a[0], label + " coarse") + exact(self, window(s, 0, b[1]), a[1], label + " fine") + exact(self, window(s, 0, b[3]), a[3], label + " raw fine") + + def test_two_bodies_and_jets_equal_tiled_single(self): + """Translation copy across a block seam, with jets and an amplitude switch.""" + from CelerisLab.body.jets import Jet, Parallel, Tube + from CelerisLab.cuda.context import CudaContext + cfg16 = LBMConfig(nx=16, ny=16, viscosity=.1) + cfg32 = replace(cfg16, nx=32) + one = FixedPatchSet2D.from_config(cfg16, self.PATCHES[:1]) + two = FixedPatchSet2D.from_config(cfg32, self.PATCHES) + coarse, fine = state(16, 16, 3), state(*one.patches[0].fine_shape, 4) + jet = lambda x: Jet(Tube((x, 8.), 0., 6., .6), Parallel(0.), amplitude=.002) + kw = dict(fine_block=(4, 4), base_block=(4, 4)) + with CudaContext(0) as context: + for jets in ((), (jet(5.),)): + ref = self.run_steps(context, cfg16, self.PATCHES[:1], coarse, level_fine(one, [fine]), + switch=[(0, -.001)] if jets else None, + fine_objects=[body(8., 8.)], jets=jets, **kw) + both = jets and (jet(5.), jet(21.)) + got = self.run_steps(context, cfg32, self.PATCHES, np.concatenate([coarse, coarse], 2), + level_fine(two, [fine, fine]), + switch=[(0, -.001), (1, -.001)] if jets else None, + fine_objects=[body(8., 8., 5), body(24., 8., 5)], + jets=both or (), **kw) + for step, a, b in zip(self.STEPS, ref, got): + label = f"jets={len(jets)} step {step}" + exact(self, b[0], np.concatenate([a[0], a[0]], axis=2), label + " coarse") + for i in (0, 1): + exact(self, window(two, i, b[1]), window(one, 0, a[1]), f"{label} fine {i}") + exact(self, window(two, i, b[3]), window(one, 0, a[3]), f"{label} raw fine {i}") + + if __name__ == "__main__": unittest.main()