Compare commits

...
Author SHA1 Message Date
ClaudeandClaude Opus 5.5 ae53fade13 Contract: name the in-EDGE half-cell projection; repair DFG runner current-checkout Jet path
- No NODE<->EDGE input conversion still leaves the root-edge -> fine-local
  half-cell projection (prepare_fine_boundary, fine_hits_to_edge_frame); say so.
- run_refined_dfg.py (jet_scenario != none, current checkout only): pass
  level=1 to FixedPatchStepper.set_jet_amplitude on reversal; use
  blocked_coarse_tables there. jet_scenario=none keeps the pinned source's
  blocked_base_tables; manifests and results unchanged.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-28 02:45:30 +00:00
ClaudeandClaude Opus 5.5 358bed1962 Naming/contract delta: _root_cfg, _device_boundaries, blocked_coarse_tables; NODE vs EDGE frame
- Coarse/fine name the two members of any pair (Hierarchy docstring); only the
  misleading blocked_base_tables -> blocked_coarse_tables. prepare_fine_boundary
  and blocked_fine_tables keep their pair-relative names.
- Depth-1 static circles use the integer-NODE frame (unchanged behavior); the
  root EDGE frame applies to refined levels only. Contract text corrected.
- _kernel_configs documents why it stays distinct from Hierarchy.configs.
- Drop the stale re-export comment; two tests import table builders from
  refinement_tables. Mechanical caller renames for this delta only.
No numerical code or config derivation changed.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-28 02:34:07 +00:00
18 changed files with 109 additions and 74 deletions
+19 -4
View File
@@ -8,8 +8,16 @@ pending a GPU run).
- Level `l = 0 .. depth-1` (0 = root). Pair `hierarchy.pairs[l-1]` builds level `l`.
- Array shapes are `(ny, nx)`; rectangles, origins and centers are `(x, y)`.
- FixedPatchStepper bodies/Jets use the root EDGE frame (root cell `i` = `[i, i+1)`);
uniform `Simulation` objects use the integer-NODE frame. Nothing converts between them.
- Coarse/fine name the two members (parent/child) of ANY pair, not absolute L0/L1:
`fine_shape`, `fine_origin`, `prepare_fine_boundary`, `blocked_fine_tables`,
`blocked_coarse_tables`.
- Refined levels (l >= 1) of FixedPatchStepper place bodies/Jets in the root EDGE frame
(root cell `i` = `[i, i+1)`). Uniform `Simulation` objects AND the static circles of a
single-level (depth 1) FixedPatchStepper use the integer-NODE frame. Nothing converts
between them. This concerns input placement only: within the EDGE frame, refined
levels still project root-edge coordinates to fine-local cell indices with the
half-cell shift (`prepare_fine_boundary`: `(v - o)/dx - .5`; hits back via `fine_hits_to_edge_frame`
`origin + dx*(hit + .5)`), using `edge_frame`/`fine_origin` below.
- `Hierarchy.edge_frame(l)` -> `((ox, oy), h)`: cell `(i, j)` of level `l`'s allocated
array (ghost band and holes included) covers `[o + h*i, o + h*(i+1))`. A scalar
FixedPatch2D window carries the half-coarse-cell ghost shift `bounds - HALO/2` in
@@ -26,8 +34,11 @@ pending a GPU run).
| `base_layout`, `fine_layout`, `child_layout` (stepper properties) | `layouts[0]`, `layouts[1]`, `layouts[2]` (`FixedPatchSet2D.fine_layout(block)` is unchanged) |
| `fine_boundary`, `child_boundary` | `boundaries[1]`, `boundaries[2]` (host geometry; `None` = no bodies/Jets) |
| `_jets` | `jets[1]` |
| `_boundaries[l]` (dict from L1, level keys) | `_level_boundaries[l]` (list from L0, `None` at L0 unless single-level static circles) |
| `_fine_tables`, `_actuated`, `_jet_cache` | `_level_boundaries[1].tables` / `.actuated` / `.jet_cache` |
| `_boundaries[l]` (dict from L1, level keys) | `_device_boundaries[l]` (list from L0, `None` at L0 unless single-level static circles) |
| `_level_boundaries` (private, b4eacf7) | `_device_boundaries` |
| `_fine_tables`, `_actuated`, `_jet_cache` | `_device_boundaries[1].tables` / `.actuated` / `.jet_cache` |
| `_remap_cfg` (private) | `_root_cfg` |
| `refinement_tables.blocked_base_tables` (also via `refinement_stepper`) | `refinement_tables.blocked_coarse_tables` (import from `refinement_tables`; the stepper no longer re-exports) |
| `_fine_leaf`, `_child_leaf` | `leaf[1]`, `leaf[2]` (copies) |
| `_middle_kind` | `hierarchy.kinds[1]` |
| `_boundary_hosts` (private) | `_boundaries` (private list behind `boundaries`) |
@@ -56,6 +67,10 @@ the replaced level relative to its parent).
## Not done / open
- Two per-level config derivations stay distinct on purpose: `Hierarchy.configs`
(`_middle_config`, wall fields reset on every refined level, geometry validation) vs.
`FixedPatchStepper._kernel_configs` (reset only on level 1 of a slot_reservoir root;
generates each level's config.h). Folding them needs a config.h/PTX diff gate.
- No single-`Simulation` facade over FixedPatchStepper: the two own different state models
(uniform node frame, body sync, forcing, checkpoint vs. immutable blocked hierarchy, edge
frame, remap-as-new-instance). A facade would only forward calls; no concrete need found.
+9 -2
View File
@@ -26,8 +26,15 @@ class Hierarchy:
Index conventions (the single place they are defined):
* level l = 0 (root) .. depth-1; pair index p = l-1 builds level l from level p;
* array shapes are (ny, nx); rectangles, origins and centers are (x, y);
* the root EDGE frame puts root cell i at [i, i+1), centers at i+.5 (bodies and
Jets of FixedPatchStepper); uniform Simulation objects use the integer-NODE frame;
* a pair's two members are its COARSE (parent) and FINE (child) level, for every
pair, not absolute L0/L1 (fine_shape, fine_origin, blocked_coarse_tables, ...);
* refined levels (l >= 1) of FixedPatchStepper place bodies and Jets in the root
EDGE frame: root cell i covers [i, i+1), center i+.5. Uniform Simulation objects
and the static circles of a single-level (depth 1) FixedPatchStepper use the
integer-NODE frame (CircleGeometry.build_flag_mask at integer (x, y)). No
NODE<->EDGE input conversion is applied anywhere; inside the EDGE frame, refined
boundaries still project root-edge coordinates to fine-local cell indices with
the half-cell shift ((v - o)/h - .5, see prepare_fine_boundary);
* cell (i, j) of level l's allocated array (ghost band and holes included) covers
[o + h*i, o + h*(i+1)) with (o, h) = edge_frame(l); level 1 is built from the root
config directly, deeper levels from their parent frame (see edge_frame).
+27 -22
View File
@@ -22,7 +22,7 @@ from .channel import channel_masks, validate_channel_config, validate_channel_tr
from ..body.jets import jet_link_velocity
from .hierarchy import build_hierarchy
from .refinement import ACTIVE, HALO, FixedPatch2D, FixedPatchSet2D, prepare_fine_boundary
from .refinement_tables import (blocked_base_tables, blocked_fine_tables, # noqa: F401 (re-export)
from .refinement_tables import (blocked_coarse_tables, blocked_fine_tables,
blocked_kind, projected_explosion_tables, _prepare_blocked_boundary,
_BlockedBoundary,
BLOCKED_BOUNDARY_DEVICE_COLUMNS)
@@ -66,8 +66,12 @@ class _Level:
def _kernel_configs(cfg, pairs, wall_aware):
"""Per-level kernel LBMConfig (window shape, viscosity doubling per level).
Geometry configs live in Hierarchy.configs; level 1 of a slot_reservoir root
replaces the root wall fields by the internal bounce-back/velocity ghost setup.
Deliberately distinct from Hierarchy.configs (refinement._middle_config), which
resets the wall fields on EVERY refined level for geometry validation. Here only
level 1 of a slot_reservoir root resets them (internal bounce-back/velocity ghost
setup); other levels inherit their parent's fields. These configs generate each
level's compiled config.h, so folding the two derivations is a numerical-parity
change (config.h/PTX diff gate), not a style one.
"""
configs = [replace(cfg)]
for index, pair in enumerate(pairs):
@@ -119,7 +123,8 @@ class FixedPatchStepper:
supported).
``refinements=[]`` is one blocked level: periodic SRT root without LES,
Jets or exchange observation; its ``level_objects[0]`` are stationary
passive circles (static BFL).
passive circles (static BFL) in the integer-NODE frame, like Simulation
objects (see Hierarchy for frames).
``coarse_boundary`` is ``periodic`` (default) or ``slot_reservoir``:
a bottom on-site slot and top fixed-density reservoir from ``cfg``.
@@ -131,7 +136,7 @@ class FixedPatchStepper:
Periodic mode does not apply the wall/slot boundary fields in ``cfg``.
Outer walls remain coarse; fine halo donors must be ordinary fluid.
Bodies are fixed circles in root (L0) physical coordinates; velocities are
Bodies of refined levels are fixed circles in root (L0) physical coordinates; velocities are
in root lattice units (dx=dt=1). Independent local copies are projected at
fine cell centers; their complete stencil must be active and disjoint from
the actual two-hop transfer footprint. These use the edge/cell-center frame
@@ -258,7 +263,7 @@ class FixedPatchStepper:
raise ValueError(f"level {level} blocks exceed launch limits")
self._layouts.append(layout)
if self._channel:
base = blocked_base_tables(pairs[0], self._layouts[0], np.argsort(pairs[0].coalesce_coarse))
base = blocked_coarse_tables(pairs[0], self._layouts[0], np.argsort(pairs[0].coalesce_coarse))
validate_channel_transfers(pairs[0], self._layouts[0], base[1])
self._level_cfgs = _kernel_configs(cfg, pairs, self._wall_aware)
# Host boundary geometry per level in the root edge frame (Hierarchy.edge_frame).
@@ -294,7 +299,7 @@ class FixedPatchStepper:
if depth == 1 and level_objects[0]:
tables = static_circle_tables(self._layouts[0], level_objects[0])
root = _BlockedBoundary(None, tables, (), False)
self._level_boundaries = [root] + [
self._device_boundaries = [root] + [
_prepare_blocked_boundary(self._boundaries[level], self._layouts[level],
self._level_cfgs[level], level_jets[level],
_shifted_blocked=self._ddf_shifting)
@@ -310,7 +315,7 @@ class FixedPatchStepper:
if self._channel:
self._leaf[0] &= channel_masks((cfg.ny, cfg.nx))[0]
# Private constructor inputs for remap_finest.
self._remap_cfg = replace(cfg)
self._root_cfg = replace(cfg)
self._level_blocks = level_blocks
self.ctx = context
self._lock = RLock()
@@ -366,7 +371,7 @@ class FixedPatchStepper:
for index, layout in enumerate(self._layouts):
total += 2 * 4 * int(np.prod(layout.blocked_shape))
total += layout.halo_dst.nbytes + layout.halo_src.nbytes
b = self._level_boundaries[index]
b = self._device_boundaries[index]
if b is not None:
total += sum(np.asarray(v).nbytes for k, v in b.tables.items()
if k == "flags" or k in BLOCKED_BOUNDARY_DEVICE_COLUMNS)
@@ -389,19 +394,19 @@ class FixedPatchStepper:
grid = ((config.nx + self._block[0] - 1) // self._block[0], config.ny, layout.nblocks)
level = _Level(shape, module, self._allocate(size), self._allocate(size),
module.get_function(name), grid, layout)
level.boundary = self._level_boundaries[index]
level.boundary = self._device_boundaries[index]
level.halo_copy = module.get_function("BlockHaloCopy")
if len(layout.halo_dst):
level.halo_dst = self._upload(layout.halo_dst)
level.halo_src = self._upload(layout.halo_src)
if index == 0 and self._channel:
cfg = self._remap_cfg
cfg = self._root_cfg
ptr, size = module.get_global("d_block_channel")
if size != struct.calcsize("=IIf"):
raise RuntimeError("Block channel host/device layout mismatch")
cuda.memcpy_htod(ptr, struct.pack("=IIf", cfg.nx, cfg.ny, cfg.velocity))
if index == 0 and self._wall_aware:
cfg = self._remap_cfg
cfg = self._root_cfg
ptr, size = module.get_global("d_block_wall")
if size != struct.calcsize("=IIIf"):
raise RuntimeError("Block wall host/device layout mismatch")
@@ -415,11 +420,11 @@ class FixedPatchStepper:
order = np.argsort(pair.coalesce_coarse)
ghosts, children, child_kind = blocked_fine_tables(pair, child.layout, order)
if index == 0:
donors, targets, parent_kind = blocked_base_tables(pair, parent.layout, order,
donors, targets, parent_kind = blocked_coarse_tables(pair, parent.layout, order,
self._wall_aware)
else:
# Nested pairs keep inherited INNER/OUTER and sparse holes of the parent.
donors, targets, parent_kind = blocked_base_tables(
donors, targets, parent_kind = blocked_coarse_tables(
pair, parent.layout, order, composed_kind=self.hierarchy.kinds[index])
result = _Pair(child.module.get_function("BlockExplode"),
parent.module.get_function("BlockCoalesce"),
@@ -552,12 +557,12 @@ class FixedPatchStepper:
if (self._wall_aware or self._channel or self.hierarchy.depth != 3
or not isinstance(pairs[1], FixedPatchSet2D)
or self._boundaries[1] is not None
or any(b is not None for b in self._level_boundaries[:2])):
or any(b is not None for b in self._device_boundaries[:2])):
raise ValueError("remap_finest requires a periodic three-level hierarchy with a "
"level-2 set and no level-0/1 bodies")
cfg = self._remap_cfg
cfg = self._root_cfg
objects = tuple(objects)
actual = self._level_boundaries[2]
actual = self._device_boundaries[2]
jets = () if actual is None else actual.jets
old = RemapTopology.from_config(cfg, pairs[0], pairs[1].bounds,
child_objects=objects, child_jets=jets,
@@ -628,7 +633,7 @@ class FixedPatchStepper:
if (isinstance(level, bool) or not isinstance(level, Integral)
or not 1 <= level < self.hierarchy.depth):
raise ValueError("Jet level must name an existing refined level (1 .. depth-1)")
boundary = self._level_boundaries[level]
boundary = self._device_boundaries[level]
jets = list(self._level_jets[level])
jets[index] = jets[index].with_amplitude(amplitude) # validates first
velocity = None
@@ -735,7 +740,7 @@ class FixedPatchStepper:
self._exchange_valid = False
for _ in range(n):
if self._observe_raw_exchange: # stream-ordered: holds the last interval only
for b in self._level_boundaries:
for b in self._device_boundaries:
if b is not None:
cuda.memset_d8_async(b.interval_integral, 0,
len(b.tables["cell"]) * 32, self._stream)
@@ -758,7 +763,7 @@ class FixedPatchStepper:
raise RuntimeError("window exchange observation is disabled")
self._window_valid = False
try:
for b in self._level_boundaries:
for b in self._device_boundaries:
if b is not None and b.window_integral is not None:
cuda.memset_d8(b.window_integral, 0, len(b.tables["cell"]) * 32)
except BaseException:
@@ -771,7 +776,7 @@ class FixedPatchStepper:
if (isinstance(level, bool) or not isinstance(level, Integral)
or not 0 <= level < len(self._levels)):
raise ValueError("level must name an existing level")
b = self._level_boundaries[level]
b = self._device_boundaries[level]
if (b is None or isinstance(body_id, bool) or not isinstance(body_id, Integral)
or not 0 <= body_id < len(b.host.objects)):
raise ValueError("body_id must name a body on this level")
@@ -868,7 +873,7 @@ class FixedPatchStepper:
self._stream.synchronize()
except Exception as exc:
error = exc
for b in getattr(self, "_level_boundaries", ()):
for b in getattr(self, "_device_boundaries", ()):
if b is not None:
b.flags = b.device = b.kernel = None
b.last_raw = b.interval_integral = b.accumulate = None
+1 -1
View File
@@ -36,7 +36,7 @@ def blocked_kind(layout, dense_kind):
kind[owner[owned]] = np.asarray(dense_kind).ravel()[owned]
return kind
def blocked_base_tables(patch, layout, order, wall_aware=False, *, composed_kind=None):
def blocked_coarse_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
+2 -2
View File
@@ -7,7 +7,7 @@ from CelerisLab.config import LBMConfig
from CelerisLab.lbm.blocks import BlockLayout2D
from CelerisLab.lbm.channel import channel_masks, validate_channel_config, validate_channel_transfers
from CelerisLab.lbm.refinement import FixedPatch2D
from CelerisLab.lbm.refinement_tables import blocked_base_tables
from CelerisLab.lbm.refinement_tables import blocked_coarse_tables
from .test_numerical_operators import compile_run
@@ -31,7 +31,7 @@ def transfer(bounds, block=(4, 6)):
cfg = LBMConfig(nx=24, ny=20, viscosity=.1)
pair = FixedPatch2D.from_config(cfg, bounds)
layout = BlockLayout2D.build(20, 24, *block)
targets = blocked_base_tables(pair, layout, np.argsort(pair.coalesce_coarse))[1]
targets = blocked_coarse_tables(pair, layout, np.argsort(pair.coalesce_coarse))[1]
return pair, layout, targets
+5 -5
View File
@@ -102,7 +102,7 @@ def test_simulation_before_initialize_commits_host_only(monkeypatch):
def stepper_state(run):
boundary = run._level_boundaries[2]
boundary = run._device_boundaries[2]
return boundary.tables["surface_velocity"].copy(), boundary.jets, run._level_jets[2]
@@ -116,8 +116,8 @@ def test_stepper_invalid_input_writes_nothing_and_stays_ready():
with pytest.raises(error):
run.set_jet_amplitude(*args, **kw)
module.cuda.memcpy_htod.assert_not_called()
np.testing.assert_array_equal(run._level_boundaries[2].tables["surface_velocity"], velocity)
assert (run._level_boundaries[2].jets, run._level_jets[2]) == (jets, level_jets)
np.testing.assert_array_equal(run._device_boundaries[2].tables["surface_velocity"], velocity)
assert (run._device_boundaries[2].jets, run._level_jets[2]) == (jets, level_jets)
assert run._status == "ready"
run.set_jet_amplitude(0, .017, level=2)
assert run._level_jets[2][0].amplitude == .017
@@ -136,8 +136,8 @@ def test_stepper_device_failure_keeps_host_and_blocks_reuse(where):
run._stream.synchronize.side_effect = RuntimeError("upload failed")
with pytest.raises(RuntimeError, match="upload failed"):
run.set_jet_amplitude(0, .017, level=2)
np.testing.assert_array_equal(run._level_boundaries[2].tables["surface_velocity"], velocity)
assert (run._level_boundaries[2].jets, run._level_jets[2]) == (jets, level_jets)
np.testing.assert_array_equal(run._device_boundaries[2].tables["surface_velocity"], velocity)
assert (run._device_boundaries[2].jets, run._level_jets[2]) == (jets, level_jets)
assert run._status == "failed"
module.cuda.memcpy_htod.side_effect = None
run._stream.synchronize.side_effect = None
+2 -2
View File
@@ -85,7 +85,7 @@ class RawExchange(unittest.TestCase):
def test_four_level_schedule_publication_and_fixed_storage(self):
run, host, calls, arrays = backend()
b = run._level_boundaries[3]
b = run._device_boundaries[3]
rows = len(b.tables['cell'])
self.assertEqual(b.last_raw.size, rows * 16)
self.assertEqual(b.interval_integral.size, rows * 32)
@@ -128,7 +128,7 @@ class RawExchange(unittest.TestCase):
def test_disabled_has_nulls_and_no_observation_allocation_or_launch(self):
run, host, calls, _ = backend(False)
b = run._level_boundaries[3]
b = run._device_boundaries[3]
self.assertIs(b.last_raw, b.interval_integral)
self.assertIs(b.interval_integral, None)
run.step()
@@ -141,7 +141,7 @@ class Gate(unittest.TestCase):
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)
donors,targets,bkind=host.blocked_coarse_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])
@@ -158,10 +158,10 @@ class Gate(unittest.TestCase):
hole=int(np.flatnonzero(~layout.owned.ravel())[0])
bad=child.ghost_coarse.copy();bad[0]=hole
with self.assertRaisesRegex(ValueError,'owned uncovered'):
host.blocked_base_tables(replace(child,ghost_coarse=bad),layout,order,composed_kind=kind)
host.blocked_coarse_tables(replace(child,ghost_coarse=bad),layout,order,composed_kind=kind)
bad=child.coalesce_coarse.copy();bad[0]=hole
with self.assertRaisesRegex(ValueError,'owned uncovered'):
host.blocked_base_tables(replace(child,coalesce_coarse=bad),layout,order,composed_kind=kind)
host.blocked_coarse_tables(replace(child,coalesce_coarse=bad),layout,order,composed_kind=kind)
def bodies_and_jets(speed=.002):
+4 -4
View File
@@ -67,7 +67,7 @@ class ShiftedJets(unittest.TestCase):
def test_switch_upload_only_then_raw_remap_keeps_live_velocity(self):
run,module,allocations,storage,objects=ready();clone=None
try:
boundary=run._level_boundaries[2]
boundary=run._device_boundaries[2]
before=run._download_raw()
original=boundary.tables['surface_velocity'].copy()
module.cuda.memcpy_htod.reset_mock()
@@ -79,7 +79,7 @@ class ShiftedJets(unittest.TestCase):
for a,b in zip(before.values,run._download_raw().values):self.assertEqual(a.tobytes(),b.tobytes())
clone=run.remap_child(B,child_objects=objects)
self.assertEqual(clone._level_jets[2][0].amplitude,-jet().amplitude)
np.testing.assert_array_equal(clone._level_boundaries[2].tables['surface_velocity'],
np.testing.assert_array_equal(clone._device_boundaries[2].tables['surface_velocity'],
boundary.tables['surface_velocity'])
self.assertEqual(clone.step_count,7)
saved=run._download_raw()
@@ -92,12 +92,12 @@ class ShiftedJets(unittest.TestCase):
def test_invalid_amplitude_does_not_upload_or_commit(self):
run,module,_,_,_=ready()
try:
before=run._level_boundaries[2].tables['surface_velocity'].copy()
before=run._device_boundaries[2].tables['surface_velocity'].copy()
module.cuda.memcpy_htod.reset_mock()
for amplitude in (float('nan'),float('inf'),1.):
with self.assertRaises(ValueError):run.set_jet_amplitude(0,amplitude,level=2)
module.cuda.memcpy_htod.assert_not_called()
np.testing.assert_array_equal(before,run._level_boundaries[2].tables['surface_velocity'])
np.testing.assert_array_equal(before,run._device_boundaries[2].tables['surface_velocity'])
finally:run.close()
+1 -1
View File
@@ -73,7 +73,7 @@ class SingleLevel(unittest.TestCase):
def test_circle_solids_are_not_state(self):
run, _, _ = backend((circle(),))
solid = run._level_boundaries[0].tables["solid"]
solid = run._device_boundaries[0].tables["solid"]
self.assertTrue(solid.any())
np.testing.assert_array_equal(run._state[0], ~solid)
@@ -240,7 +240,7 @@ class SingleLevelGPU(unittest.TestCase):
with CudaContext(0) as context:
for split in SPLITS:
with single(context=context, block=split, objects=[circle()]) as backend:
tables = backend._level_boundaries[0].tables
tables = backend._device_boundaries[0].tables
for key, values, dtype in (("q", [l.q for l in links], np.float32),
("dir", [l.dir for l in links], np.uint8),
("fallback", [l.fallback for l in links], np.uint8)):
@@ -59,7 +59,7 @@ class ExchangeWindowCPU(unittest.TestCase):
stack.enter_context(patch.object(cuda, "memcpy_dtoh",
side_effect=lambda out, src: out.fill(0)))
backend = build(m)
b = backend._level_boundaries[1]
b = backend._device_boundaries[1]
self.assertIsNotNone(b.window_integral)
with self.assertRaises(RuntimeError):
backend.start_exchange_window() # not initialized
@@ -96,7 +96,7 @@ class ExchangeWindowCPU(unittest.TestCase):
with mocked_backend() as m, ExitStack() as stack:
stack.enter_context(patch.object(cuda, "memset_d8_async"))
backend = build(m, window=False)
self.assertIsNone(backend._level_boundaries[1].window_integral)
self.assertIsNone(backend._device_boundaries[1].window_integral)
backend.initialize(*[np.full(l.shape, 1 / 9, np.float32) for l in backend._levels])
backend.step(1)
self.assertEqual(len(accumulations(m)), 2)
@@ -120,7 +120,7 @@ class ExchangeWindowGPU(unittest.TestCase):
windowed.step(2)
plain.step(2)
windowed.start_exchange_window()
b = windowed._level_boundaries[1]
b = windowed._device_boundaries[1]
count = len(b.tables["cell"])
factors = np.array([.25, .25, .125, .25]) # h*dt, h*h*dt, h*h for level-1 h=.5
scale = np.zeros(4)
@@ -212,7 +212,7 @@ class ExchangeWindowShiftedMRTGPU(unittest.TestCase):
run.step(2)
if run is not plain:
run.start_exchange_window()
b = windowed._level_boundaries[1]
b = windowed._device_boundaries[1]
count = len(b.tables["cell"])
record = dict(scale=np.zeros(4), substeps=0)
audit_substeps(windowed, b, count, record)
@@ -228,7 +228,7 @@ class ExchangeWindowShiftedMRTGPU(unittest.TestCase):
same_bits(self, a, e)
rows = [np.empty((count, 4), np.float64) for _ in range(2)]
cuda.memcpy_dtoh(rows[0], b.window_integral)
cuda.memcpy_dtoh(rows[1], split._level_boundaries[1].window_integral)
cuda.memcpy_dtoh(rows[1], split._device_boundaries[1].window_integral)
np.testing.assert_array_equal(rows[0].view(np.uint64), rows[1].view(np.uint64))
out = windowed.read_exchange_window(0)
@@ -105,7 +105,7 @@ class RawExchangeGPU(unittest.TestCase):
arrays = [state(l.shape) for l in measured._levels]
measured.initialize(*arrays)
plain.initialize(*arrays)
b = measured._level_boundaries[1]
b = measured._device_boundaries[1]
count = len(b.tables['cell'])
kernel = measured._levels[1].kernel
samples = []
@@ -158,7 +158,7 @@ class RawExchangeGPU(unittest.TestCase):
FixedPatchStepper(cfg, context=ctx, **settings) as plain:
arrays = [state(l.shape, ux=0.) for l in run._levels]
run.initialize(*arrays); plain.initialize(*arrays)
b = run._level_boundaries[3]
b = run._device_boundaries[3]
count = len(b.tables['cell'])
kernel = run._levels[3].kernel
samples, independent_samples = [], []
@@ -213,7 +213,7 @@ class RawExchangeGPU(unittest.TestCase):
FixedPatchStepper(CFG, context=ctx, **fine((16, 20)), observe_raw_exchange=True) as amr:
amr.initialize(*(state(l.shape) for l in amr._levels))
uniform.initialize(state((9, 48, 64)))
b = amr._level_boundaries[1]
b = amr._device_boundaries[1]
count = len(b.tables['cell'])
links = obj.geometry.build_cut_links(64, 48)
np.testing.assert_array_equal([v.dir for v in links], b.tables['dir'])
@@ -221,7 +221,7 @@ class RawExchangeGPU(unittest.TestCase):
# Same fine-grid force arms after translation to the same center: the
# single level's static BFL gets the AMR Jet column and force arms; the raw
# output pointer is substituted at launch (no stepper-side accumulation).
root = uniform._level_boundaries[0]
root = uniform._device_boundaries[0]
velocity = uniform._upload(b.tables['surface_velocity'])
root.device.update(surface_velocity=velocity, rx=uniform._upload(b.tables['rx']),
ry=uniform._upload(b.tables['ry']))
@@ -24,7 +24,8 @@ from CelerisLab.lbm import refinement_stepper
from CelerisLab.lbm.blocks import BlockLayout2D
from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY
from CelerisLab.lbm.refinement import ACTIVE, OUTER_GHOST, FixedPatch2D
from CelerisLab.lbm.refinement_stepper import FixedPatchStepper, blocked_base_tables
from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
from CelerisLab.lbm.refinement_tables import blocked_coarse_tables
from test_refinement_blocked_fine import NX, NY, BOUNDS, config, state, circle, actuation
BASE_BLOCKS = ((24, 32), (12, 8), (3, 4), (1, 1))
@@ -108,7 +109,7 @@ class BlockedBaseCPU(unittest.TestCase):
leaves[[0, -1]] = False
for by, bx in BASE_BLOCKS:
layout = BlockLayout2D.build(NY, NX, by, bx, periodic=(True, mode == "periodic"))
donors, targets, kind = blocked_base_tables(p, layout, order, mode != "periodic")
donors, targets, kind = blocked_coarse_tables(p, layout, order, mode != "periodic")
ly, lx = layout.local_shape
L = ly * lx
expected, source_cells = [], []
@@ -138,7 +139,7 @@ class BlockedBaseCPU(unittest.TestCase):
p = FixedPatch2D.from_config(config("periodic"), BOUNDS)
order = np.argsort(p.coalesce_coarse)
layout = BlockLayout2D.build(NY, NX, 3, 4)
donors, targets, _ = blocked_base_tables(p, layout, order)
donors, targets, _ = blocked_coarse_tables(p, layout, order)
dense = np.arange(9 * NY * NX, dtype=np.float32).reshape(9, NY, NX)
initial = layout.to_blocked(dense)
arrivals = np.arange(len(targets), dtype=np.float32) + 10000
@@ -178,7 +179,7 @@ class BlockedBaseCPU(unittest.TestCase):
# With one parent block, all eight grazing sources are owner cores.
# Each is read for four fine ghost children: late explosion corrupts32.
single = BlockLayout2D.build(NY, NX, NY, NX)
donors, targets, _ = blocked_base_tables(p, single, order)
donors, targets, _ = blocked_coarse_tables(p, single, order)
initial = single.to_blocked(dense)
staged = initial.copy()
staged.ravel()[targets.astype(np.intp)] = arrivals
@@ -21,7 +21,8 @@ from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY
from CelerisLab.lbm.descriptors import SOLID
from CelerisLab.lbm.refinement import (ACTIVE, INNER_GHOST, OUTER_GHOST, FixedPatch2D,
prepare_fine_boundary)
from CelerisLab.lbm.refinement_stepper import FixedPatchStepper, blocked_fine_tables
from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
from CelerisLab.lbm.refinement_tables import blocked_fine_tables
NX, NY, BOUNDS, STEPS = 32, 24, (8, 6, 24, 18), 8
FINE_BLOCKS = ((32, 40), (16, 20), (16, 10)) # fine window 32x40: 1x1, 2x2, 4x2
@@ -186,8 +187,8 @@ class BlockedFineCPU(unittest.TestCase):
backend, calls = self.mocked_backend([circle(spec)])
c, f = backend._levels
coarse, cur, spare = c.current, f.current, f.spare
flags, kind = backend._level_boundaries[1].flags, backend._levels[1].kind
velocity = backend._level_boundaries[1].device.get("surface_velocity")
flags, kind = backend._device_boundaries[1].flags, backend._levels[1].kind
velocity = backend._device_boundaries[1].device.get("surface_velocity")
self.assertEqual(velocity is None, spec is None)
backend.step(1)
backend.close()
@@ -213,7 +214,7 @@ class BlockedFineCPU(unittest.TestCase):
for speed in (0.0, .002): # any spec, even zero amplitude
backend, calls = self.mocked_backend([circle(actuation(speed))])
self.assertIsNotNone(backend._fine_tables["surface_velocity"])
del backend._level_boundaries[1].device["surface_velocity"]
del backend._device_boundaries[1].device["surface_velocity"]
with self.assertRaisesRegex(RuntimeError, "Surface velocity device buffer is missing"):
backend.step(1)
names = [name for name, _ in calls]
@@ -107,7 +107,7 @@ class JetStepperCPU(unittest.TestCase):
def test_blocked_amplitude_update_rewrites_only_velocity_column(self):
backend, htod = mocked([circle()], [jet()], FINE_BLOCKS[1])
before = backend._fine_tables["surface_velocity"].copy()
device = backend._level_boundaries[1].device["surface_velocity"]
device = backend._device_boundaries[1].device["surface_velocity"]
with patch.object(cuda, "memcpy_htod", htod):
htod.reset_mock()
backend.set_jet_amplitude(0, -SPEED)
@@ -102,7 +102,7 @@ class Ledger:
def step(self):
run=self.run;self.samples.clear();run.step()
b=run._level_boundaries[self.index];h=run._exchange_spacing[self.index]
b=run._device_boundaries[self.index];h=run._exchange_spacing[self.index]
self.test.assertEqual(len(self.samples),2**self.index)
expected=np.zeros_like(self.samples[0],dtype=np.float64)
weights=np.array([h*h,h*h,h*h*h,h*h])
@@ -140,7 +140,7 @@ class ShiftedJetsGPU(unittest.TestCase):
level_objects=[(),(obj,)],level_jets=[(),(whole,)],observe_raw_exchange=True) as run:
run._initialize_raw(initial(run,rho,u=.013,v=0.))
expected=run._download_raw();ledger=Ledger(self,run,1)
self.assertTrue(run._level_boundaries[1].tables['aperture'].all())
self.assertTrue(run._device_boundaries[1].tables['aperture'].all())
for _ in range(4):ledger.step()
error=max(float(abs(a[:,leaf]-b[:,leaf]).max()) for a,b,leaf in
zip(run._download_raw().values,expected.values,run._leaf))
@@ -201,7 +201,7 @@ class ShiftedJetsGPU(unittest.TestCase):
self.exact_raw(before,old._download_raw())
self.exact_raw(_StoredPopulations(expected.g,True),new._download_raw())
self.assertEqual(new._level_jets[2][0].amplitude,-.001)
np.testing.assert_array_equal(new._level_boundaries[2].tables['surface_velocity'],old._level_boundaries[2].tables['surface_velocity'])
np.testing.assert_array_equal(new._device_boundaries[2].tables['surface_velocity'],old._device_boundaries[2].tables['surface_velocity'])
newledger=Ledger(self,new,2);ledger.step();newledger.step()
self.assertEqual((old.step_count,new.step_count),(count+1,count+1))
ledger.report('L2 source');newledger.report('L2 remapped')
+12 -6
View File
@@ -81,7 +81,7 @@ def public_jets(center, radius, scenario):
)
def jet_interval(handle, scenario, begin, n, diameter):
def jet_interval(handle, scenario, begin, n, diameter, level=None):
"""Apply a reversal only at a completed interval boundary, before advancing."""
if scenario != "none":
boundaries = [round(t * diameter / SPEED) for t in (250, 260, 270, 280, 290, 300, 305)]
@@ -92,8 +92,10 @@ def jet_interval(handle, scenario, begin, n, diameter):
if begin < switch:
return min(n, switch - begin), None
if begin == switch:
handle.set_jet_amplitude(0, -JET_PEAK)
handle.set_jet_amplitude(1, JET_PEAK)
# FixedPatchStepper needs an explicit Jet level; Simulation has one level.
kw = {} if level is None else {"level": level}
handle.set_jet_amplitude(0, -JET_PEAK, **kw)
handle.set_jet_amplitude(1, JET_PEAK, **kw)
return n, dict(step=begin, tstar=300.0, amplitudes=[-JET_PEAK, JET_PEAK])
return n, None
@@ -256,7 +258,7 @@ def host_geometry(source, h, collision="SRT", jet_scenario="none", placement="wi
validate_channel_config,
validate_channel_transfers,
)
from CelerisLab.lbm.refinement_tables import blocked_base_tables
from CelerisLab.lbm import refinement_tables
from CelerisLab.body.geometry.circle import CircleGeometry
from CelerisLab.body.objects import SimObject
@@ -268,7 +270,11 @@ def host_geometry(source, h, collision="SRT", jet_scenario="none", placement="wi
pair = FixedPatch2D.from_config(cfg, PLACEMENTS[placement])
layout = BlockLayout2D.build(cfg.ny, cfg.nx, *BASE_BLOCK)
validate_channel_config(cfg, [BASE_BLOCK, FINE_BLOCK])
targets = blocked_base_tables(pair, layout, np.argsort(pair.coalesce_coarse))[1]
# Historical pinned source (jet_scenario none) names it blocked_base_tables;
# the current checkout (public Jet scenarios) names it blocked_coarse_tables.
coarse_tables = (refinement_tables.blocked_base_tables if jet_scenario == "none"
else refinement_tables.blocked_coarse_tables)
targets = coarse_tables(pair, layout, np.argsort(pair.coalesce_coarse))[1]
validate_channel_transfers(pair, layout, targets)
BlockLayout2D.build(*pair.fine_shape, *FINE_BLOCK, periodic=(False, False))
body = SimObject(
@@ -448,7 +454,7 @@ def run_case(a, h, periodic, case):
)
n = min(max(1, round(width * p["diameter"] / SPEED)), max_steps - begin)
n, change = jet_interval(backend if refined else sim, a.jet_scenario,
begin, n, p["diameter"])
begin, n, p["diameter"], level=1 if refined else None)
if change is not None:
amplitude_changes.append(change)
impulse = np.zeros(2)