Compare commits

...
Author SHA1 Message Date
ClaudeandClaude Opus 5.5 65339ab54e test_runtime_child_remap: inject remap body rejects through the boundary record
The removed _boundaries list made the two body-reject patches no-ops.
Inject them via _device_boundaries (level-2 record with foreign host
geometry; non-None level-1 record). set_jet_amplitude docstring names
the public jets view.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-28 10:22:29 +00:00
ClaudeandClaude Opus 5.5 ae766f7112 FixedPatchStepper: one boundary/Jet record per level, host preparation helpers
- Jets and host boundary geometry are views of the per-level _BlockedBoundary
  record (drop duplicate _level_jets/_boundaries registration).
- Host preparation split into _level_layouts, _host_boundaries, _fluid_leaves
  (numpy only, before any device work); constructor order unchanged.
- set_jet_amplitude: explicit prepare -> device write -> commit; behaviour
  unchanged (invalid input raises before any change, uncertain write fails closed).
- Pair accessor: FixedPatch2D.rectangles, FixedPatchSet2D.rectangles and
  fine_level_offset=(0, 0); removes two incidental isinstance branches.
  Geometry/representation unchanged; scalar/set types stay distinct.
- @property read-only views; blocks.py responsibility docstring corrected.
- Migrate test_runtime_child_remap to the public boundaries view.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-28 10:09:25 +00:00
ClaudeandClaude Opus 5.5 5f69ebaa90 run_refined_dfg: pick the coarse table builder by the imported source's API
jet_scenario=none host tests run against the current checkout (no
blocked_base_tables), while historical runs import the pinned source (no
blocked_coarse_tables). Select whichever single name the source defines;
fail if neither or both. No src alias; pins and manifests unchanged.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-28 05:11:34 +00:00
ClaudeandClaude Opus 5.5 9d0588dcad Merge codex/src-coherence-callers (7b652b2) into claude/src-coherence (ae53fad)
Conflicts (6 files): tests took the caller-migration side, then the
358bed1 renames (_level_boundaries -> _device_boundaries); run_refined_dfg.py
kept ae53fad (same jet_interval level= repair plus blocked_coarse_tables on the
current-checkout path). Follow-up renames in 5 auto-merged tests;
blocked_fine_tables imported from refinement_tables. Frozen drivers untouched
(run_jet_circle_cv.py, check_projected_pairs.py remap_child, historical DFG path).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-28 04:45:48 +00:00
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
26 changed files with 250 additions and 170 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`. - 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)`. - 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)`); - Coarse/fine name the two members (parent/child) of ANY pair, not absolute L0/L1:
uniform `Simulation` objects use the integer-NODE frame. Nothing converts between them. `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 - `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 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 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) | | `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) | | `fine_boundary`, `child_boundary` | `boundaries[1]`, `boundaries[2]` (host geometry; `None` = no bodies/Jets) |
| `_jets` | `jets[1]` | | `_jets` | `jets[1]` |
| `_boundaries[l]` (dict from L1, level keys) | `_level_boundaries[l]` (list from L0, `None` at L0 unless single-level static circles) | | `_boundaries[l]` (dict from L1, level keys) | `_device_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` | | `_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) | | `_fine_leaf`, `_child_leaf` | `leaf[1]`, `leaf[2]` (copies) |
| `_middle_kind` | `hierarchy.kinds[1]` | | `_middle_kind` | `hierarchy.kinds[1]` |
| `_boundary_hosts` (private) | `_boundaries` (private list behind `boundaries`) | | `_boundary_hosts` (private) | `_boundaries` (private list behind `boundaries`) |
@@ -56,6 +67,10 @@ the replaced level relative to its parent).
## Not done / open ## 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 - No single-`Simulation` facade over FixedPatchStepper: the two own different state models
(uniform node frame, body sync, forcing, checkpoint vs. immutable blocked hierarchy, edge (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. frame, remap-as-new-instance). A facade would only forward calls; no concrete need found.
+3 -1
View File
@@ -18,7 +18,9 @@ all NQ populations from an owner core cell into another block's halo cell.
Periodicity is explicit per axis. Across a non-periodic exterior, halo cells Periodicity is explicit per axis. Across a non-periodic exterior, halo cells
have no pair and stay invalid (NaN after packing): no wrap, no implied BC. have no pair and stay invalid (NaN after packing): no wrap, no implied BC.
Preparation only: no CUDA, no boundary conditions, no refinement levels. Host preparation plus thin launch wrappers: no CUDA import or allocation
(launch_* only call kernels and device tables the caller owns), no refinement
levels; level_config and static_circle_tables prepare one blocked level.
The one-cell same-level halo here is unrelated to the four-cell fine halo of The one-cell same-level halo here is unrelated to the four-cell fine halo of
the fixed refinement patch (refinement.py), whose semantics are unchanged. the fixed refinement patch (refinement.py), whose semantics are unchanged.
+10 -3
View File
@@ -26,8 +26,15 @@ class Hierarchy:
Index conventions (the single place they are defined): 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; * 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); * 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 * a pair's two members are its COARSE (parent) and FINE (child) level, for every
Jets of FixedPatchStepper); uniform Simulation objects use the integer-NODE frame; 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 * 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 [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). config directly, deeper levels from their parent frame (see edge_frame).
@@ -115,7 +122,7 @@ def leaf_partition(h):
for level, kind in enumerate(h.kinds): for level, kind in enumerate(h.kinds):
if level: if level:
pair = h.pairs[level - 1] pair = h.pairs[level - 1]
py, px = (pair.fine_level_offset if isinstance(pair, FixedPatch2D) else (0, 0)) py, px = pair.fine_level_offset
oy, ox = 2 * oy + py, 2 * ox + px oy, ox = 2 * oy + py, 2 * ox + px
y, x = np.nonzero(kind == ACTIVE) y, x = np.nonzero(kind == ACTIVE)
np.add.at(base, ((y + oy) >> level, (x + ox) >> level), 0.25 ** level) np.add.at(base, ((y + oy) >> level, (x + ox) >> level), 0.25 ** level)
+12
View File
@@ -120,6 +120,11 @@ class FixedPatch2D:
(an index offset, not an (x, y) origin).""" (an index offset, not an (x, y) origin)."""
return 2 * self.bounds[1] - HALO, 2 * self.bounds[0] - HALO return 2 * self.bounds[1] - HALO, 2 * self.bounds[0] - HALO
@property
def rectangles(self) -> tuple:
"""The one parent rectangle, as a sequence like FixedPatchSet2D.rectangles."""
return (self.bounds,)
@property @property
def covered(self) -> np.ndarray: def covered(self) -> np.ndarray:
"""Dense coarse mask of cells replaced by the patch (not leaves).""" """Dense coarse mask of cells replaced by the patch (not leaves)."""
@@ -225,11 +230,18 @@ class FixedPatchSet2D:
def bounds(self) -> tuple: def bounds(self) -> tuple:
return tuple(p.bounds for p in self.patches) return tuple(p.bounds for p in self.patches)
rectangles = bounds
@property @property
def fine_origin(self) -> tuple[float, float]: def fine_origin(self) -> tuple[float, float]:
"""Lower edge of the level-wide fine grid in coarse physical units.""" """Lower edge of the level-wide fine grid in coarse physical units."""
return 0.0, 0.0 return 0.0, 0.0
@property
def fine_level_offset(self) -> tuple[int, int]:
"""The level-wide grid itself: zero (y, x) index offset."""
return 0, 0
@property @property
def covered(self) -> np.ndarray: def covered(self) -> np.ndarray:
mask = np.zeros(self.coarse_shape, bool) mask = np.zeros(self.coarse_shape, bool)
+116 -85
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 ..body.jets import jet_link_velocity
from .hierarchy import build_hierarchy from .hierarchy import build_hierarchy
from .refinement import ACTIVE, HALO, FixedPatch2D, FixedPatchSet2D, prepare_fine_boundary 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, blocked_kind, projected_explosion_tables, _prepare_blocked_boundary,
_BlockedBoundary, _BlockedBoundary,
BLOCKED_BOUNDARY_DEVICE_COLUMNS) BLOCKED_BOUNDARY_DEVICE_COLUMNS)
@@ -66,8 +66,12 @@ class _Level:
def _kernel_configs(cfg, pairs, wall_aware): def _kernel_configs(cfg, pairs, wall_aware):
"""Per-level kernel LBMConfig (window shape, viscosity doubling per level). """Per-level kernel LBMConfig (window shape, viscosity doubling per level).
Geometry configs live in Hierarchy.configs; level 1 of a slot_reservoir root Deliberately distinct from Hierarchy.configs (refinement._middle_config), which
replaces the root wall fields by the internal bounce-back/velocity ghost setup. 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)] configs = [replace(cfg)]
for index, pair in enumerate(pairs): for index, pair in enumerate(pairs):
@@ -107,6 +111,57 @@ def _check_levels(refinements, blocks, objects, jets):
raise ValueError("bodies/Jets are supported only on the finest level of a nested hierarchy") raise ValueError("bodies/Jets are supported only on the finest level of a nested hierarchy")
# Host preparation (numpy only, before any device work); each returns one entry per level.
def _level_layouts(cfg, pairs, level_blocks, wall_aware):
"""Blocked storage layout of every level, checked against the launch limits."""
layouts = []
for level, block in enumerate(level_blocks):
if level == 0:
layout = BlockLayout2D.build(cfg.ny, cfg.nx, *block, periodic=(True, not wall_aware))
elif isinstance(pairs[level - 1], FixedPatchSet2D):
layout = pairs[level - 1].fine_layout(block)
else:
layout = BlockLayout2D.build(*pairs[level - 1].fine_shape, *block,
periodic=(False, False))
if layout.nblocks > MAX_GRID_DIM or layout.local_shape[0] > MAX_GRID_DIM:
raise ValueError(f"level {level} blocks exceed launch limits")
layouts.append(layout)
return layouts
def _host_boundaries(hierarchy, level_cfgs, level_objects, level_jets, wall_aware, shifted):
"""Host boundary geometry of every refined level in the root edge frame (None at L0)."""
hosts = [None]
for level in range(1, hierarchy.depth):
pair = hierarchy.pairs[level - 1]
if level == 1: # the root config carries the wall-aware boundary fields
host = prepare_fine_boundary(level_cfgs[0], pair, level_objects[1],
wall_aware=wall_aware, jets=level_jets[1],
_shifted_blocked=shifted)
else:
origin, spacing = hierarchy.edge_frame(level - 1)
host = prepare_fine_boundary(level_cfgs[level - 1], pair, level_objects[level],
jets=level_jets[level], parent_origin=origin,
parent_spacing=spacing, _shifted_blocked=shifted)
hosts.append(host)
return hosts
def _fluid_leaves(hierarchy, boundaries):
"""ACTIVE cells of every level that no static circle or body makes solid."""
leaves = [kind == ACTIVE for kind in hierarchy.kinds]
for level, b in enumerate(boundaries):
if b is None:
continue
if b.host is None: # single level: static root circles
leaves[level] &= ~b.tables["solid"]
else:
fine_shape = hierarchy.pairs[level - 1].fine_shape
leaves[level] &= (b.host.flags.reshape(fine_shape) & SOLID) == 0
return leaves
class FixedPatchStepper: class FixedPatchStepper:
"""D2Q9 SRT/MRT FP32 fixed 2:1 refinement hierarchy, built from level lists. """D2Q9 SRT/MRT FP32 fixed 2:1 refinement hierarchy, built from level lists.
@@ -119,7 +174,8 @@ class FixedPatchStepper:
supported). supported).
``refinements=[]`` is one blocked level: periodic SRT root without LES, ``refinements=[]`` is one blocked level: periodic SRT root without LES,
Jets or exchange observation; its ``level_objects[0]`` are stationary 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``: ``coarse_boundary`` is ``periodic`` (default) or ``slot_reservoir``:
a bottom on-site slot and top fixed-density reservoir from ``cfg``. a bottom on-site slot and top fixed-density reservoir from ``cfg``.
@@ -131,7 +187,7 @@ class FixedPatchStepper:
Periodic mode does not apply the wall/slot boundary fields in ``cfg``. Periodic mode does not apply the wall/slot boundary fields in ``cfg``.
Outer walls remain coarse; fine halo donors must be ordinary fluid. 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 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 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 the actual two-hop transfer footprint. These use the edge/cell-center frame
@@ -241,40 +297,15 @@ class FixedPatchStepper:
self.hierarchy = build_hierarchy(cfg, refinements, wall_aware=self._wall_aware) self.hierarchy = build_hierarchy(cfg, refinements, wall_aware=self._wall_aware)
pairs = self.hierarchy.pairs pairs = self.hierarchy.pairs
if self._wall_aware: if self._wall_aware:
rects = (pairs[0].bounds,) if isinstance(pairs[0], FixedPatch2D) else pairs[0].bounds if any(r[1] < 4 or r[3] > cfg.ny - 4 for r in pairs[0].rectangles):
if any(r[1] < 4 or r[3] > cfg.ny - 4 for r in rects):
raise ValueError("slot_reservoir requires y0 >= 4 and y1 <= ny-4") raise ValueError("slot_reservoir requires y0 >= 4 and y1 <= ny-4")
self._layouts = [] self._layouts = _level_layouts(cfg, pairs, level_blocks, self._wall_aware)
for level, block in enumerate(level_blocks):
if level == 0:
layout = BlockLayout2D.build(cfg.ny, cfg.nx, *block,
periodic=(True, not self._wall_aware))
elif isinstance(pairs[level - 1], FixedPatchSet2D):
layout = pairs[level - 1].fine_layout(block)
else:
layout = BlockLayout2D.build(*pairs[level - 1].fine_shape, *block,
periodic=(False, False))
if layout.nblocks > MAX_GRID_DIM or layout.local_shape[0] > MAX_GRID_DIM:
raise ValueError(f"level {level} blocks exceed launch limits")
self._layouts.append(layout)
if self._channel: 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]) validate_channel_transfers(pairs[0], self._layouts[0], base[1])
self._level_cfgs = _kernel_configs(cfg, pairs, self._wall_aware) self._level_cfgs = _kernel_configs(cfg, pairs, self._wall_aware)
# Host boundary geometry per level in the root edge frame (Hierarchy.edge_frame). hosts = _host_boundaries(self.hierarchy, self._level_cfgs, level_objects, level_jets,
self._boundaries = [None] self._wall_aware, self._ddf_shifting)
for level in range(1, depth):
pair = pairs[level - 1]
if level == 1: # the root config carries the wall-aware boundary fields
host = prepare_fine_boundary(cfg, pair, level_objects[1],
wall_aware=self._wall_aware, jets=level_jets[1],
_shifted_blocked=self._ddf_shifting)
else:
origin, spacing = self.hierarchy.edge_frame(level - 1)
host = prepare_fine_boundary(self._level_cfgs[level - 1], pair, level_objects[level],
jets=level_jets[level], parent_origin=origin,
parent_spacing=spacing, _shifted_blocked=self._ddf_shifting)
self._boundaries.append(host)
if not isinstance(observe_raw_exchange, bool): if not isinstance(observe_raw_exchange, bool):
raise ValueError("observe_raw_exchange must be bool") raise ValueError("observe_raw_exchange must be bool")
self._observe_raw_exchange = observe_raw_exchange self._observe_raw_exchange = observe_raw_exchange
@@ -289,28 +320,22 @@ class FixedPatchStepper:
self._exchange_spacing = tuple(self.hierarchy.edge_frame(l)[1] for l in range(depth)) self._exchange_spacing = tuple(self.hierarchy.edge_frame(l)[1] for l in range(depth))
self._exchange_centers = tuple(tuple(tuple(o.geometry.center) for o in objects) self._exchange_centers = tuple(tuple(tuple(o.geometry.center) for o in objects)
for objects in level_objects) for objects in level_objects)
self._level_jets = level_jets
root = None root = None
if depth == 1 and level_objects[0]: if depth == 1 and level_objects[0]:
tables = static_circle_tables(self._layouts[0], level_objects[0]) tables = static_circle_tables(self._layouts[0], level_objects[0])
root = _BlockedBoundary(None, tables, (), False) root = _BlockedBoundary(None, tables, (), False)
self._level_boundaries = [root] + [ # The one boundary/Jet record per level; host geometry and Jets are views of it.
_prepare_blocked_boundary(self._boundaries[level], self._layouts[level], self._device_boundaries = [root] + [
_prepare_blocked_boundary(hosts[level], self._layouts[level],
self._level_cfgs[level], level_jets[level], self._level_cfgs[level], level_jets[level],
_shifted_blocked=self._ddf_shifting) _shifted_blocked=self._ddf_shifting)
for level in range(1, depth)] for level in range(1, depth)]
self._leaf = [kind == ACTIVE for kind in self.hierarchy.kinds] self._leaf = _fluid_leaves(self.hierarchy, self._device_boundaries)
if root is not None:
self._leaf[0] &= ~root.tables["solid"]
for level in range(1, depth):
host = self._boundaries[level]
if host is not None:
self._leaf[level] &= (host.flags.reshape(pairs[level - 1].fine_shape) & SOLID) == 0
self._state = [mask.copy() for mask in self._leaf] self._state = [mask.copy() for mask in self._leaf]
if self._channel: if self._channel:
self._leaf[0] &= channel_masks((cfg.ny, cfg.nx))[0] self._leaf[0] &= channel_masks((cfg.ny, cfg.nx))[0]
# Private constructor inputs for remap_finest. # Private constructor inputs for remap_finest.
self._remap_cfg = replace(cfg) self._root_cfg = replace(cfg)
self._level_blocks = level_blocks self._level_blocks = level_blocks
self.ctx = context self.ctx = context
self._lock = RLock() self._lock = RLock()
@@ -366,7 +391,7 @@ class FixedPatchStepper:
for index, layout in enumerate(self._layouts): for index, layout in enumerate(self._layouts):
total += 2 * 4 * int(np.prod(layout.blocked_shape)) total += 2 * 4 * int(np.prod(layout.blocked_shape))
total += layout.halo_dst.nbytes + layout.halo_src.nbytes total += layout.halo_dst.nbytes + layout.halo_src.nbytes
b = self._level_boundaries[index] b = self._device_boundaries[index]
if b is not None: if b is not None:
total += sum(np.asarray(v).nbytes for k, v in b.tables.items() total += sum(np.asarray(v).nbytes for k, v in b.tables.items()
if k == "flags" or k in BLOCKED_BOUNDARY_DEVICE_COLUMNS) if k == "flags" or k in BLOCKED_BOUNDARY_DEVICE_COLUMNS)
@@ -389,19 +414,19 @@ class FixedPatchStepper:
grid = ((config.nx + self._block[0] - 1) // self._block[0], config.ny, layout.nblocks) grid = ((config.nx + self._block[0] - 1) // self._block[0], config.ny, layout.nblocks)
level = _Level(shape, module, self._allocate(size), self._allocate(size), level = _Level(shape, module, self._allocate(size), self._allocate(size),
module.get_function(name), grid, layout) 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") level.halo_copy = module.get_function("BlockHaloCopy")
if len(layout.halo_dst): if len(layout.halo_dst):
level.halo_dst = self._upload(layout.halo_dst) level.halo_dst = self._upload(layout.halo_dst)
level.halo_src = self._upload(layout.halo_src) level.halo_src = self._upload(layout.halo_src)
if index == 0 and self._channel: if index == 0 and self._channel:
cfg = self._remap_cfg cfg = self._root_cfg
ptr, size = module.get_global("d_block_channel") ptr, size = module.get_global("d_block_channel")
if size != struct.calcsize("=IIf"): if size != struct.calcsize("=IIf"):
raise RuntimeError("Block channel host/device layout mismatch") raise RuntimeError("Block channel host/device layout mismatch")
cuda.memcpy_htod(ptr, struct.pack("=IIf", cfg.nx, cfg.ny, cfg.velocity)) cuda.memcpy_htod(ptr, struct.pack("=IIf", cfg.nx, cfg.ny, cfg.velocity))
if index == 0 and self._wall_aware: if index == 0 and self._wall_aware:
cfg = self._remap_cfg cfg = self._root_cfg
ptr, size = module.get_global("d_block_wall") ptr, size = module.get_global("d_block_wall")
if size != struct.calcsize("=IIIf"): if size != struct.calcsize("=IIIf"):
raise RuntimeError("Block wall host/device layout mismatch") raise RuntimeError("Block wall host/device layout mismatch")
@@ -415,11 +440,11 @@ class FixedPatchStepper:
order = np.argsort(pair.coalesce_coarse) order = np.argsort(pair.coalesce_coarse)
ghosts, children, child_kind = blocked_fine_tables(pair, child.layout, order) ghosts, children, child_kind = blocked_fine_tables(pair, child.layout, order)
if index == 0: 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) self._wall_aware)
else: else:
# Nested pairs keep inherited INNER/OUTER and sparse holes of the parent. # 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]) pair, parent.layout, order, composed_kind=self.hierarchy.kinds[index])
result = _Pair(child.module.get_function("BlockExplode"), result = _Pair(child.module.get_function("BlockExplode"),
parent.module.get_function("BlockCoalesce"), parent.module.get_function("BlockCoalesce"),
@@ -551,19 +576,17 @@ class FixedPatchStepper:
pairs = self.hierarchy.pairs pairs = self.hierarchy.pairs
if (self._wall_aware or self._channel or self.hierarchy.depth != 3 if (self._wall_aware or self._channel or self.hierarchy.depth != 3
or not isinstance(pairs[1], FixedPatchSet2D) or not isinstance(pairs[1], FixedPatchSet2D)
or self._boundaries[1] is not None or any(b is not None for b in self._device_boundaries[:2])):
or any(b is not None for b in self._level_boundaries[:2])):
raise ValueError("remap_finest requires a periodic three-level hierarchy with a " raise ValueError("remap_finest requires a periodic three-level hierarchy with a "
"level-2 set and no level-0/1 bodies") "level-2 set and no level-0/1 bodies")
cfg = self._remap_cfg cfg = self._root_cfg
objects = tuple(objects) objects = tuple(objects)
actual = self._level_boundaries[2] actual = self._device_boundaries[2]
jets = () if actual is None else actual.jets jets = () if actual is None else actual.jets
old = RemapTopology.from_config(cfg, pairs[0], pairs[1].bounds, old = RemapTopology.from_config(cfg, pairs[0], pairs[1].bounds,
child_objects=objects, child_jets=jets, child_objects=objects, child_jets=jets,
_shifted_blocked=self._ddf_shifting) _shifted_blocked=self._ddf_shifting)
if (not same_physical_boundary(old.boundary, self._boundaries[2]) if not same_physical_boundary(old.boundary, None if actual is None else actual.host):
or (self._boundaries[2] is not None and actual is None)):
raise ValueError("objects must reproduce the original root-edge-frame level-2 boundary") raise ValueError("objects must reproduce the original root-edge-frame level-2 boundary")
if actual is not None: if actual is not None:
velocity = actual.tables.get("surface_velocity") velocity = actual.tables.get("surface_velocity")
@@ -617,7 +640,7 @@ class FixedPatchStepper:
return candidate return candidate
def set_jet_amplitude(self, index, amplitude, *, level): def set_jet_amplitude(self, index, amplitude, *, level):
"""Replace Jet ``index`` of ``level_jets[level]`` between steps; no geometry/table rebuild. """Replace Jet ``index`` of level ``level`` (``jets[level]``) between steps; no table rebuild.
Invalid input raises before any change. A device/sync error marks the Invalid input raises before any change. A device/sync error marks the
instance failed (as step/initialize do); host state commits only after instance failed (as step/initialize do); host state commits only after
@@ -628,27 +651,24 @@ class FixedPatchStepper:
if (isinstance(level, bool) or not isinstance(level, Integral) if (isinstance(level, bool) or not isinstance(level, Integral)
or not 1 <= level < self.hierarchy.depth): or not 1 <= level < self.hierarchy.depth):
raise ValueError("Jet level must name an existing refined level (1 .. depth-1)") 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]) # Host preparation: validate and build the new column without touching state.
jets = list(() if boundary is None else boundary.jets)
jets[index] = jets[index].with_amplitude(amplitude) # validates first jets[index] = jets[index].with_amplitude(amplitude) # validates first
velocity = None shape, jet_id = boundary.jet_cache
if boundary is not None and boundary.jet_cache is not None: velocity = boundary.tables["surface_velocity"].copy()
shape, jet_id = boundary.jet_cache driven = jet_id >= 0
velocity = boundary.tables["surface_velocity"].copy() velocity[driven] = jet_link_velocity(jets, shape, jet_id)[driven]
driven = jet_id >= 0 # Device action: an uncertain write fails the instance closed.
velocity[driven] = jet_link_velocity(jets, shape, jet_id)[driven]
try: try:
if velocity is not None: self._stream.synchronize() # no in-flight step reads the column
self._stream.synchronize() # no in-flight step reads the column cuda.memcpy_htod(boundary.device["surface_velocity"], velocity)
cuda.memcpy_htod(boundary.device["surface_velocity"], velocity)
except BaseException: except BaseException:
self._status = "failed" self._status = "failed"
raise raise
if velocity is not None: # Commit: the one boundary record holds both the column and the Jets.
boundary.tables["surface_velocity"] = velocity boundary.tables["surface_velocity"] = velocity
if boundary is not None: boundary.jets = tuple(jets)
boundary.jets = tuple(jets)
self._level_jets[level] = tuple(jets)
def _blocked_boundary(self, level, buffer, **kw): def _blocked_boundary(self, level, buffer, **kw):
b = level.boundary b = level.boundary
@@ -735,7 +755,7 @@ class FixedPatchStepper:
self._exchange_valid = False self._exchange_valid = False
for _ in range(n): for _ in range(n):
if self._observe_raw_exchange: # stream-ordered: holds the last interval only 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: if b is not None:
cuda.memset_d8_async(b.interval_integral, 0, cuda.memset_d8_async(b.interval_integral, 0,
len(b.tables["cell"]) * 32, self._stream) len(b.tables["cell"]) * 32, self._stream)
@@ -758,7 +778,7 @@ class FixedPatchStepper:
raise RuntimeError("window exchange observation is disabled") raise RuntimeError("window exchange observation is disabled")
self._window_valid = False self._window_valid = False
try: 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: if b is not None and b.window_integral is not None:
cuda.memset_d8(b.window_integral, 0, len(b.tables["cell"]) * 32) cuda.memset_d8(b.window_integral, 0, len(b.tables["cell"]) * 32)
except BaseException: except BaseException:
@@ -771,7 +791,7 @@ class FixedPatchStepper:
if (isinstance(level, bool) or not isinstance(level, Integral) if (isinstance(level, bool) or not isinstance(level, Integral)
or not 0 <= level < len(self._levels)): or not 0 <= level < len(self._levels)):
raise ValueError("level must name an existing level") 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) if (b is None or isinstance(body_id, bool) or not isinstance(body_id, Integral)
or not 0 <= body_id < len(b.host.objects)): or not 0 <= body_id < len(b.host.objects)):
raise ValueError("body_id must name a body on this level") raise ValueError("body_id must name a body on this level")
@@ -868,7 +888,7 @@ class FixedPatchStepper:
self._stream.synchronize() self._stream.synchronize()
except Exception as exc: except Exception as exc:
error = exc error = exc
for b in getattr(self, "_level_boundaries", ()): for b in getattr(self, "_device_boundaries", ()):
if b is not None: if b is not None:
b.flags = b.device = b.kernel = None b.flags = b.device = b.kernel = None
b.last_raw = b.interval_integral = b.accumulate = None b.last_raw = b.interval_integral = b.accumulate = None
@@ -896,9 +916,20 @@ class FixedPatchStepper:
raise raise
# Read-only per-level state, indexed by level (0 = root). Pairs: hierarchy.pairs[l-1]. # Read-only per-level state, indexed by level (0 = root). Pairs: hierarchy.pairs[l-1].
layouts = property(lambda self: tuple(self._layouts)) @property
boundaries = property(lambda self: tuple(self._boundaries), def layouts(self):
doc="Host boundary geometry per level (None: no bodies/Jets).") return tuple(self._layouts)
jets = property(lambda self: tuple(self._level_jets))
leaf = property(lambda self: tuple(m.copy() for m in self._leaf), @property
doc="Per-level (ny, nx) masks of public (non-NaN) download cells.") def boundaries(self):
"""Host boundary geometry per level (None: no bodies/Jets on a refined level)."""
return tuple(None if b is None else b.host for b in self._device_boundaries)
@property
def jets(self):
return tuple(() if b is None else b.jets for b in self._device_boundaries)
@property
def leaf(self):
"""Per-level (ny, nx) masks of public (non-NaN) download cells."""
return tuple(m.copy() for m in self._leaf)
+1 -1
View File
@@ -36,7 +36,7 @@ def blocked_kind(layout, dense_kind):
kind[owner[owned]] = np.asarray(dense_kind).ravel()[owned] kind[owner[owned]] = np.asarray(dense_kind).ravel()[owned]
return kind 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. """Parent donor cells and receiver-local arrival slots, in transfer row order.
Stage arrivals only after child reads and parent halo gather, immediately 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.blocks import BlockLayout2D
from CelerisLab.lbm.channel import channel_masks, validate_channel_config, validate_channel_transfers from CelerisLab.lbm.channel import channel_masks, validate_channel_config, validate_channel_transfers
from CelerisLab.lbm.refinement import FixedPatch2D 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 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) cfg = LBMConfig(nx=24, ny=20, viscosity=.1)
pair = FixedPatch2D.from_config(cfg, bounds) pair = FixedPatch2D.from_config(cfg, bounds)
layout = BlockLayout2D.build(20, 24, *block) 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 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): def stepper_state(run):
boundary = run._level_boundaries[2] boundary = run._device_boundaries[2]
return boundary.tables["surface_velocity"].copy(), boundary.jets, run.jets[2] return boundary.tables["surface_velocity"].copy(), boundary.jets, run.jets[2]
@@ -116,8 +116,8 @@ def test_stepper_invalid_input_writes_nothing_and_stays_ready():
with pytest.raises(error): with pytest.raises(error):
run.set_jet_amplitude(*args, **kw) run.set_jet_amplitude(*args, **kw)
module.cuda.memcpy_htod.assert_not_called() module.cuda.memcpy_htod.assert_not_called()
np.testing.assert_array_equal(run._level_boundaries[2].tables["surface_velocity"], velocity) np.testing.assert_array_equal(run._device_boundaries[2].tables["surface_velocity"], velocity)
assert (run._level_boundaries[2].jets, run.jets[2]) == (jets, level_jets) assert (run._device_boundaries[2].jets, run.jets[2]) == (jets, level_jets)
assert run._status == "ready" assert run._status == "ready"
run.set_jet_amplitude(0, .017, level=2) run.set_jet_amplitude(0, .017, level=2)
assert run.jets[2][0].amplitude == .017 assert run.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") run._stream.synchronize.side_effect = RuntimeError("upload failed")
with pytest.raises(RuntimeError, match="upload failed"): with pytest.raises(RuntimeError, match="upload failed"):
run.set_jet_amplitude(0, .017, level=2) run.set_jet_amplitude(0, .017, level=2)
np.testing.assert_array_equal(run._level_boundaries[2].tables["surface_velocity"], velocity) np.testing.assert_array_equal(run._device_boundaries[2].tables["surface_velocity"], velocity)
assert (run._level_boundaries[2].jets, run.jets[2]) == (jets, level_jets) assert (run._device_boundaries[2].jets, run.jets[2]) == (jets, level_jets)
assert run._status == "failed" assert run._status == "failed"
module.cuda.memcpy_htod.side_effect = None module.cuda.memcpy_htod.side_effect = None
run._stream.synchronize.side_effect = None run._stream.synchronize.side_effect = None
+4 -4
View File
@@ -96,11 +96,11 @@ class MultiplePinnedRemap(unittest.TestCase):
before = run.download() before = run.download()
clone = run.remap_finest((PINS+MOVED+ADDED)[::-1], objects=bodies) clone = run.remap_finest((PINS+MOVED+ADDED)[::-1], objects=bodies)
try: try:
self.assertEqual(tuple(j.amplitude for j in clone._level_boundaries[2].jets), (.017, -.011)) self.assertEqual(tuple(j.amplitude for j in clone._device_boundaries[2].jets), (.017, -.011))
np.testing.assert_array_equal(clone._level_boundaries[2].tables['surface_velocity'], np.testing.assert_array_equal(clone._device_boundaries[2].tables['surface_velocity'],
run._level_boundaries[2].tables['surface_velocity']) run._device_boundaries[2].tables['surface_velocity'])
clone.set_jet_amplitude(1, .005, level=2) clone.set_jet_amplitude(1, .005, level=2)
self.assertEqual(run._level_boundaries[2].jets[1].amplitude, -.011) self.assertEqual(run._device_boundaries[2].jets[1].amplitude, -.011)
clone.step() clone.step()
self.assertEqual((run.step_count, clone.step_count), (7, 8)) self.assertEqual((run.step_count, clone.step_count), (7, 8))
for got, want in zip(run.download(), before): for got, want in zip(run.download(), before):
+5 -5
View File
@@ -62,9 +62,9 @@ class PinnedBodyRemap(unittest.TestCase):
run.set_jet_amplitude(0,.017,level=2) run.set_jet_amplitude(0,.017,level=2)
before=run.download() before=run.download()
clone=run.remap_finest(B,objects=objects) clone=run.remap_finest(B,objects=objects)
self.assertEqual(clone._level_boundaries[2].jets[0].amplitude,.017) self.assertEqual(clone._device_boundaries[2].jets[0].amplitude,.017)
np.testing.assert_array_equal(clone._level_boundaries[2].tables['surface_velocity'], np.testing.assert_array_equal(clone._device_boundaries[2].tables['surface_velocity'],
run._level_boundaries[2].tables['surface_velocity']) run._device_boundaries[2].tables['surface_velocity'])
self.assertTrue(same_physical_boundary(run.boundaries[2],clone.boundaries[2])) self.assertTrue(same_physical_boundary(run.boundaries[2],clone.boundaries[2]))
clone.step() clone.step()
self.assertEqual((run.step_count,clone.step_count),(7,8)) self.assertEqual((run.step_count,clone.step_count),(7,8))
@@ -74,8 +74,8 @@ class PinnedBodyRemap(unittest.TestCase):
def test_static_body_without_actuation(self): def test_static_body_without_actuation(self):
run,module,allocations,storage,objects=ready(jets=()) run,module,allocations,storage,objects=ready(jets=())
clone=run.remap_finest(A,objects=objects) clone=run.remap_finest(A,objects=objects)
self.assertFalse(clone._level_boundaries[2].actuated) self.assertFalse(clone._device_boundaries[2].actuated)
self.assertNotIn('surface_velocity',clone._level_boundaries[2].tables) self.assertNotIn('surface_velocity',clone._device_boundaries[2].tables)
for a,b in zip(run.download(),clone.download()): np.testing.assert_array_equal(a,b) for a,b in zip(run.download(),clone.download()): np.testing.assert_array_equal(a,b)
clone.close();run.close() clone.close();run.close()
+2 -2
View File
@@ -85,7 +85,7 @@ class RawExchange(unittest.TestCase):
def test_four_level_schedule_publication_and_fixed_storage(self): def test_four_level_schedule_publication_and_fixed_storage(self):
run, host, calls, arrays = backend() run, host, calls, arrays = backend()
b = run._level_boundaries[3] b = run._device_boundaries[3]
rows = len(b.tables['cell']) rows = len(b.tables['cell'])
self.assertEqual(b.last_raw.size, rows * 16) self.assertEqual(b.last_raw.size, rows * 16)
self.assertEqual(b.interval_integral.size, rows * 32) 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): def test_disabled_has_nulls_and_no_observation_allocation_or_launch(self):
run, host, calls, _ = backend(False) 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.last_raw, b.interval_integral)
self.assertIs(b.interval_integral, None) self.assertIs(b.interval_integral, None)
run.step() run.step()
@@ -141,7 +141,7 @@ class Gate(unittest.TestCase):
def test_sparse_parent_tables_preserve_inherited_kinds(self): def test_sparse_parent_tables_preserve_inherited_kinds(self):
child,kind=r.nest_patch_set(CFG,self.parent,CHILDREN) child,kind=r.nest_patch_set(CFG,self.parent,CHILDREN)
layout=self.parent.fine_layout((4,4));order=np.argsort(child.coalesce_coarse) 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((donors < layout.hole).all())
self.assertTrue((targets < np.prod(layout.blocked_shape)).all()) self.assertTrue((targets < np.prod(layout.blocked_shape)).all())
np.testing.assert_array_equal(bkind[layout.owner[layout.owned]],kind[layout.owned]) 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]) hole=int(np.flatnonzero(~layout.owned.ravel())[0])
bad=child.ghost_coarse.copy();bad[0]=hole bad=child.ghost_coarse.copy();bad[0]=hole
with self.assertRaisesRegex(ValueError,'owned uncovered'): 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 bad=child.coalesce_coarse.copy();bad[0]=hole
with self.assertRaisesRegex(ValueError,'owned uncovered'): 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): def bodies_and_jets(speed=.002):
+5 -3
View File
@@ -1,5 +1,6 @@
"""CPU mocks execute production remap lifecycle; no GPU numerical claim.""" """CPU mocks execute production remap lifecycle; no GPU numerical claim."""
from dataclasses import replace from dataclasses import replace
from types import SimpleNamespace
import unittest import unittest
from unittest.mock import patch from unittest.mock import patch
@@ -69,10 +70,11 @@ class RuntimeChildRemap(unittest.TestCase):
run._status = 'ready' run._status = 'ready'
for bad in ([], [(40,16,48,40)], [(16,16,32,24),(16,27,32,55)]): for bad in ([], [(40,16,48,40)], [(16,16,32,24),(16,27,32,55)]):
with self.assertRaises(ValueError): run.remap_finest(bad) with self.assertRaises(ValueError): run.remap_finest(bad)
hosts, pairs = run._boundaries, run.hierarchy.pairs records, pairs = run._device_boundaries, run.hierarchy.pairs
for attribute, value in (('_wall_aware', True), for attribute, value in (('_wall_aware', True),
('_boundaries', [*hosts[:2], object()]), # child boundary # level-2 record whose host geometry differs from objects=()
('_boundaries', [hosts[0], object(), hosts[2]]), # fine boundary ('_device_boundaries', [*records[:2], SimpleNamespace(host=object(), jets=())]),
('_device_boundaries', [records[0], object(), records[2]]), # level-1 body
('hierarchy', replace(run.hierarchy, pairs=(pairs[0], None)))): ('hierarchy', replace(run.hierarchy, pairs=(pairs[0], None)))):
with patch.object(run, attribute, value): with patch.object(run, attribute, value):
with self.assertRaises(ValueError): run.remap_finest(B) with self.assertRaises(ValueError): run.remap_finest(B)
+4 -4
View File
@@ -67,7 +67,7 @@ class ShiftedJets(unittest.TestCase):
def test_switch_upload_only_then_raw_remap_keeps_live_velocity(self): def test_switch_upload_only_then_raw_remap_keeps_live_velocity(self):
run,module,allocations,storage,objects=ready();clone=None run,module,allocations,storage,objects=ready();clone=None
try: try:
boundary=run._level_boundaries[2] boundary=run._device_boundaries[2]
before=run._download_raw() before=run._download_raw()
original=boundary.tables['surface_velocity'].copy() original=boundary.tables['surface_velocity'].copy()
module.cuda.memcpy_htod.reset_mock() 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()) for a,b in zip(before.values,run._download_raw().values):self.assertEqual(a.tobytes(),b.tobytes())
clone=run.remap_finest(B,objects=objects) clone=run.remap_finest(B,objects=objects)
self.assertEqual(clone.jets[2][0].amplitude,-jet().amplitude) self.assertEqual(clone.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']) boundary.tables['surface_velocity'])
self.assertEqual(clone.step_count,7) self.assertEqual(clone.step_count,7)
saved=run._download_raw() saved=run._download_raw()
@@ -92,12 +92,12 @@ class ShiftedJets(unittest.TestCase):
def test_invalid_amplitude_does_not_upload_or_commit(self): def test_invalid_amplitude_does_not_upload_or_commit(self):
run,module,_,_,_=ready() run,module,_,_,_=ready()
try: 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() module.cuda.memcpy_htod.reset_mock()
for amplitude in (float('nan'),float('inf'),1.): for amplitude in (float('nan'),float('inf'),1.):
with self.assertRaises(ValueError):run.set_jet_amplitude(0,amplitude,level=2) with self.assertRaises(ValueError):run.set_jet_amplitude(0,amplitude,level=2)
module.cuda.memcpy_htod.assert_not_called() 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() finally:run.close()
+1 -1
View File
@@ -73,7 +73,7 @@ class SingleLevel(unittest.TestCase):
def test_circle_solids_are_not_state(self): def test_circle_solids_are_not_state(self):
run, _, _ = backend((circle(),)) run, _, _ = backend((circle(),))
solid = run._level_boundaries[0].tables["solid"] solid = run._device_boundaries[0].tables["solid"]
self.assertTrue(solid.any()) self.assertTrue(solid.any())
np.testing.assert_array_equal(run._state[0], ~solid) np.testing.assert_array_equal(run._state[0], ~solid)
+1 -1
View File
@@ -37,7 +37,7 @@ class UnionBodySupport(unittest.TestCase):
run.set_jet_amplitude(0,.017,level=2) run.set_jet_amplitude(0,.017,level=2)
clone=run.remap_finest(((16,16,24,40),(24,16,32,40),A[1]),objects=objects) clone=run.remap_finest(((16,16,24,40),(24,16,32,40),A[1]),objects=objects)
try: try:
self.assertEqual(clone._level_boundaries[2].jets[0].amplitude,.017) self.assertEqual(clone._device_boundaries[2].jets[0].amplitude,.017)
self.assertTrue(same_physical_boundary(run.boundaries[2],clone.boundaries[2])) self.assertTrue(same_physical_boundary(run.boundaries[2],clone.boundaries[2]))
for a,b in zip(run.download(),clone.download()):np.testing.assert_array_equal(a,b) for a,b in zip(run.download(),clone.download()):np.testing.assert_array_equal(a,b)
finally:clone.close() finally:clone.close()
@@ -240,7 +240,7 @@ class SingleLevelGPU(unittest.TestCase):
with CudaContext(0) as context: with CudaContext(0) as context:
for split in SPLITS: for split in SPLITS:
with single(context=context, block=split, objects=[circle()]) as backend: 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), for key, values, dtype in (("q", [l.q for l in links], np.float32),
("dir", [l.dir for l in links], np.uint8), ("dir", [l.dir for l in links], np.uint8),
("fallback", [l.fallback 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", stack.enter_context(patch.object(cuda, "memcpy_dtoh",
side_effect=lambda out, src: out.fill(0))) side_effect=lambda out, src: out.fill(0)))
backend = build(m) backend = build(m)
b = backend._level_boundaries[1] b = backend._device_boundaries[1]
self.assertIsNotNone(b.window_integral) self.assertIsNotNone(b.window_integral)
with self.assertRaises(RuntimeError): with self.assertRaises(RuntimeError):
backend.start_exchange_window() # not initialized backend.start_exchange_window() # not initialized
@@ -96,7 +96,7 @@ class ExchangeWindowCPU(unittest.TestCase):
with mocked_backend() as m, ExitStack() as stack: with mocked_backend() as m, ExitStack() as stack:
stack.enter_context(patch.object(cuda, "memset_d8_async")) stack.enter_context(patch.object(cuda, "memset_d8_async"))
backend = build(m, window=False) 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.initialize(*[np.full(l.shape, 1 / 9, np.float32) for l in backend._levels])
backend.step(1) backend.step(1)
self.assertEqual(len(accumulations(m)), 2) self.assertEqual(len(accumulations(m)), 2)
@@ -120,7 +120,7 @@ class ExchangeWindowGPU(unittest.TestCase):
windowed.step(2) windowed.step(2)
plain.step(2) plain.step(2)
windowed.start_exchange_window() windowed.start_exchange_window()
b = windowed._level_boundaries[1] b = windowed._device_boundaries[1]
count = len(b.tables["cell"]) count = len(b.tables["cell"])
factors = np.array([.25, .25, .125, .25]) # h*dt, h*h*dt, h*h for level-1 h=.5 factors = np.array([.25, .25, .125, .25]) # h*dt, h*h*dt, h*h for level-1 h=.5
scale = np.zeros(4) scale = np.zeros(4)
@@ -212,7 +212,7 @@ class ExchangeWindowShiftedMRTGPU(unittest.TestCase):
run.step(2) run.step(2)
if run is not plain: if run is not plain:
run.start_exchange_window() run.start_exchange_window()
b = windowed._level_boundaries[1] b = windowed._device_boundaries[1]
count = len(b.tables["cell"]) count = len(b.tables["cell"])
record = dict(scale=np.zeros(4), substeps=0) record = dict(scale=np.zeros(4), substeps=0)
audit_substeps(windowed, b, count, record) audit_substeps(windowed, b, count, record)
@@ -228,7 +228,7 @@ class ExchangeWindowShiftedMRTGPU(unittest.TestCase):
same_bits(self, a, e) same_bits(self, a, e)
rows = [np.empty((count, 4), np.float64) for _ in range(2)] rows = [np.empty((count, 4), np.float64) for _ in range(2)]
cuda.memcpy_dtoh(rows[0], b.window_integral) 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)) np.testing.assert_array_equal(rows[0].view(np.uint64), rows[1].view(np.uint64))
out = windowed.read_exchange_window(0, level=1) out = windowed.read_exchange_window(0, level=1)
@@ -57,7 +57,7 @@ def make(context):
def topology(run, rectangles, objects): def topology(run, rectangles, objects):
return RemapTopology.from_config( return RemapTopology.from_config(
CFG, run.hierarchy.pairs[0], rectangles, child_objects=objects, CFG, run.hierarchy.pairs[0], rectangles, child_objects=objects,
child_jets=run._level_boundaries[2].jets, child_jets=run._device_boundaries[2].jets,
) )
@@ -74,7 +74,7 @@ class MultiplePinnedRemapGPU(unittest.TestCase):
def change_amplitudes(self, run, amplitudes): def change_amplitudes(self, run, amplitudes):
# Each jet selects a different body. Updating one must preserve the # Each jet selects a different body. Updating one must preserve the
# other body's current per-link velocities, including its driven links. # other body's current per-link velocities, including its driven links.
boundary = run._level_boundaries[2] boundary = run._device_boundaries[2]
body_ids = run.boundaries[2].curved_host[3] body_ids = run.boundaries[2].curved_host[3]
for index, amplitude in enumerate(amplitudes): for index, amplitude in enumerate(amplitudes):
before = boundary.tables['surface_velocity'].copy() before = boundary.tables['surface_velocity'].copy()
@@ -91,7 +91,7 @@ class MultiplePinnedRemapGPU(unittest.TestCase):
self.assertTrue(same_physical_boundary(new.boundaries[2], mesh.boundary)) self.assertTrue(same_physical_boundary(new.boundaries[2], mesh.boundary))
# Compare physical cut-link rows, never layout-dependent blocked addresses. # Compare physical cut-link rows, never layout-dependent blocked addresses.
for run in (old, new): for run in (old, new):
boundary = run._level_boundaries[2] boundary = run._device_boundaries[2]
self.assertEqual(tuple(jet.amplitude for jet in boundary.jets), amplitudes) self.assertEqual(tuple(jet.amplitude for jet in boundary.jets), amplitudes)
exact(self, boundary.tables['surface_velocity'], mesh.surface_velocity, exact(self, boundary.tables['surface_velocity'], mesh.surface_velocity,
'latest physical per-link velocities') 'latest physical per-link velocities')
@@ -39,9 +39,9 @@ class PinnedBodyRemapGPU(unittest.TestCase):
exact(self,getattr(after.boundaries[2],name),getattr(before.boundaries[2],name),name) exact(self,getattr(after.boundaries[2],name),getattr(before.boundaries[2],name),name)
for a,b in zip(after.boundaries[2].curved_host,before.boundaries[2].curved_host): for a,b in zip(after.boundaries[2].curved_host,before.boundaries[2].curved_host):
exact(self,a,b,'physical cut links') exact(self,a,b,'physical cut links')
exact(self,after._level_boundaries[2].tables['surface_velocity'], exact(self,after._device_boundaries[2].tables['surface_velocity'],
before._level_boundaries[2].tables['surface_velocity'],'current per-link velocities') before._device_boundaries[2].tables['surface_velocity'],'current per-link velocities')
self.assertEqual(after._level_boundaries[2].jets[0].amplitude,.017) self.assertEqual(after._device_boundaries[2].jets[0].amplitude,.017)
def test_noop_next_step_bitwise(self): def test_noop_next_step_bitwise(self):
with CudaContext(0) as ctx, ExitStack() as owned: with CudaContext(0) as ctx, ExitStack() as owned:
@@ -63,7 +63,7 @@ class PinnedBodyRemapGPU(unittest.TestCase):
current.set_jet_amplitude(0,.017,level=2) current.set_jet_amplitude(0,.017,level=2)
for rects in (B,A): for rects in (B,A):
previous=current;before=previous.download();count=previous.step_count previous=current;before=previous.download();count=previous.step_count
jets=previous._level_boundaries[2].jets jets=previous._device_boundaries[2].jets
old=RemapTopology.from_config(CFG,previous.hierarchy.pairs[0],previous.hierarchy.pairs[1].bounds, old=RemapTopology.from_config(CFG,previous.hierarchy.pairs[0],previous.hierarchy.pairs[1].bounds,
child_objects=objects,child_jets=jets) child_objects=objects,child_jets=jets)
new=RemapTopology.from_config(CFG,previous.hierarchy.pairs[0],rects, new=RemapTopology.from_config(CFG,previous.hierarchy.pairs[0],rects,
@@ -105,7 +105,7 @@ class RawExchangeGPU(unittest.TestCase):
arrays = [state(l.shape) for l in measured._levels] arrays = [state(l.shape) for l in measured._levels]
measured.initialize(*arrays) measured.initialize(*arrays)
plain.initialize(*arrays) plain.initialize(*arrays)
b = measured._level_boundaries[1] b = measured._device_boundaries[1]
count = len(b.tables['cell']) count = len(b.tables['cell'])
kernel = measured._levels[1].kernel kernel = measured._levels[1].kernel
samples = [] samples = []
@@ -158,7 +158,7 @@ class RawExchangeGPU(unittest.TestCase):
FixedPatchStepper(cfg, context=ctx, **settings) as plain: FixedPatchStepper(cfg, context=ctx, **settings) as plain:
arrays = [state(l.shape, ux=0.) for l in run._levels] arrays = [state(l.shape, ux=0.) for l in run._levels]
run.initialize(*arrays); plain.initialize(*arrays) run.initialize(*arrays); plain.initialize(*arrays)
b = run._level_boundaries[3] b = run._device_boundaries[3]
count = len(b.tables['cell']) count = len(b.tables['cell'])
kernel = run._levels[3].kernel kernel = run._levels[3].kernel
samples, independent_samples = [], [] samples, independent_samples = [], []
@@ -213,7 +213,7 @@ class RawExchangeGPU(unittest.TestCase):
FixedPatchStepper(CFG, context=ctx, **fine((16, 20)), observe_raw_exchange=True) as amr: FixedPatchStepper(CFG, context=ctx, **fine((16, 20)), observe_raw_exchange=True) as amr:
amr.initialize(*(state(l.shape) for l in amr._levels)) amr.initialize(*(state(l.shape) for l in amr._levels))
uniform.initialize(state((9, 48, 64))) uniform.initialize(state((9, 48, 64)))
b = amr._level_boundaries[1] b = amr._device_boundaries[1]
count = len(b.tables['cell']) count = len(b.tables['cell'])
links = obj.geometry.build_cut_links(64, 48) links = obj.geometry.build_cut_links(64, 48)
np.testing.assert_array_equal([v.dir for v in links], b.tables['dir']) 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 # 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 # 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). # 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']) velocity = uniform._upload(b.tables['surface_velocity'])
root.device.update(surface_velocity=velocity, rx=uniform._upload(b.tables['rx']), root.device.update(surface_velocity=velocity, rx=uniform._upload(b.tables['rx']),
ry=uniform._upload(b.tables['ry'])) 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.blocks import BlockLayout2D
from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY
from CelerisLab.lbm.refinement import ACTIVE, OUTER_GHOST, FixedPatch2D 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 from test_refinement_blocked_fine import NX, NY, BOUNDS, config, state, circle, actuation
BASE_BLOCKS = ((24, 32), (12, 8), (3, 4), (1, 1)) BASE_BLOCKS = ((24, 32), (12, 8), (3, 4), (1, 1))
@@ -108,7 +109,7 @@ class BlockedBaseCPU(unittest.TestCase):
leaves[[0, -1]] = False leaves[[0, -1]] = False
for by, bx in BASE_BLOCKS: for by, bx in BASE_BLOCKS:
layout = BlockLayout2D.build(NY, NX, by, bx, periodic=(True, mode == "periodic")) 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 ly, lx = layout.local_shape
L = ly * lx L = ly * lx
expected, source_cells = [], [] expected, source_cells = [], []
@@ -138,7 +139,7 @@ class BlockedBaseCPU(unittest.TestCase):
p = FixedPatch2D.from_config(config("periodic"), BOUNDS) p = FixedPatch2D.from_config(config("periodic"), BOUNDS)
order = np.argsort(p.coalesce_coarse) order = np.argsort(p.coalesce_coarse)
layout = BlockLayout2D.build(NY, NX, 3, 4) 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) dense = np.arange(9 * NY * NX, dtype=np.float32).reshape(9, NY, NX)
initial = layout.to_blocked(dense) initial = layout.to_blocked(dense)
arrivals = np.arange(len(targets), dtype=np.float32) + 10000 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. # With one parent block, all eight grazing sources are owner cores.
# Each is read for four fine ghost children: late explosion corrupts32. # Each is read for four fine ghost children: late explosion corrupts32.
single = BlockLayout2D.build(NY, NX, NY, NX) 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) initial = single.to_blocked(dense)
staged = initial.copy() staged = initial.copy()
staged.ravel()[targets.astype(np.intp)] = arrivals 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.descriptors import SOLID
from CelerisLab.lbm.refinement import (ACTIVE, INNER_GHOST, OUTER_GHOST, FixedPatch2D, from CelerisLab.lbm.refinement import (ACTIVE, INNER_GHOST, OUTER_GHOST, FixedPatch2D,
prepare_fine_boundary) 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 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 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)]) backend, calls = self.mocked_backend([circle(spec)])
c, f = backend._levels c, f = backend._levels
coarse, cur, spare = c.current, f.current, f.spare coarse, cur, spare = c.current, f.current, f.spare
flags, kind = backend._level_boundaries[1].flags, backend._levels[1].kind flags, kind = backend._device_boundaries[1].flags, backend._levels[1].kind
velocity = backend._level_boundaries[1].device.get("surface_velocity") velocity = backend._device_boundaries[1].device.get("surface_velocity")
self.assertEqual(velocity is None, spec is None) self.assertEqual(velocity is None, spec is None)
backend.step(1) backend.step(1)
backend.close() backend.close()
@@ -212,8 +213,8 @@ class BlockedFineCPU(unittest.TestCase):
def test_actuation_without_velocity_buffer_fails_before_bfl(self): def test_actuation_without_velocity_buffer_fails_before_bfl(self):
for speed in (0.0, .002): # any spec, even zero amplitude for speed in (0.0, .002): # any spec, even zero amplitude
backend, calls = self.mocked_backend([circle(actuation(speed))]) backend, calls = self.mocked_backend([circle(actuation(speed))])
self.assertIsNotNone(backend._level_boundaries[1].tables["surface_velocity"]) self.assertIsNotNone(backend._device_boundaries[1].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"): with self.assertRaisesRegex(RuntimeError, "Surface velocity device buffer is missing"):
backend.step(1) backend.step(1)
names = [name for name, _ in calls] names = [name for name, _ in calls]
@@ -62,8 +62,8 @@ class JetStepperCPU(unittest.TestCase):
C = (16.0, 12.0) C = (16.0, 12.0)
whole = Jet(Sector(C, .1, math.pi, 3.), Radial(C, inward=False), amplitude=0.) whole = Jet(Sector(C, .1, math.pi, 3.), Radial(C, inward=False), amplitude=0.)
backend, _ = mocked([circle()], [whole], FINE_BLOCKS[1]) backend, _ = mocked([circle()], [whole], FINE_BLOCKS[1])
_, jet_id = backend._level_boundaries[1].jet_cache _, jet_id = backend._device_boundaries[1].jet_cache
self.assertGreater(len(np.unique(backend._level_boundaries[1].tables["block"][jet_id == 0])), 1) self.assertGreater(len(np.unique(backend._device_boundaries[1].tables["block"][jet_id == 0])), 1)
backend.close() backend.close()
def test_fine_ownership_uses_remapped_physical_owner(self): def test_fine_ownership_uses_remapped_physical_owner(self):
@@ -76,10 +76,10 @@ class JetStepperCPU(unittest.TestCase):
source = SimObject(99, geometry=CircleGeometry(*C, 2.15), center=C, radius=2.15) source = SimObject(99, geometry=CircleGeometry(*C, 2.15), center=C, radius=2.15)
whole = Jet(Sector(C, .1, math.pi, 3.), Radial(C, inward=False), amplitude=0.) whole = Jet(Sector(C, .1, math.pi, 3.), Radial(C, inward=False), amplitude=0.)
backend, _ = mocked([source], [whole], FINE_BLOCKS[1]) backend, _ = mocked([source], [whole], FINE_BLOCKS[1])
_, jet_id = backend._level_boundaries[1].jet_cache _, jet_id = backend._device_boundaries[1].jet_cache
owners = backend.boundaries[1].curved_host[3][jet_id == 0] owners = backend.boundaries[1].curved_host[3][jet_id == 0]
self.assertEqual(np.unique(owners).tolist(), [0]) self.assertEqual(np.unique(owners).tolist(), [0])
self.assertGreater(len(np.unique(backend._level_boundaries[1].tables["block"][jet_id == 0])), 1) self.assertGreater(len(np.unique(backend._device_boundaries[1].tables["block"][jet_id == 0])), 1)
self.assertEqual(source.obj_id, 99) self.assertEqual(source.obj_id, 99)
backend.close() backend.close()
cfg = config("periodic") cfg = config("periodic")
@@ -97,8 +97,8 @@ class JetStepperCPU(unittest.TestCase):
def test_blocked_jet_velocity_matches_legacy_rows(self): def test_blocked_jet_velocity_matches_legacy_rows(self):
legacy, _ = mocked([circle(actuation(SPEED))], (), FINE_BLOCKS[1]) legacy, _ = mocked([circle(actuation(SPEED))], (), FINE_BLOCKS[1])
driven, _ = mocked([circle()], [jet()], FINE_BLOCKS[1]) driven, _ = mocked([circle()], [jet()], FINE_BLOCKS[1])
a = legacy._level_boundaries[1].tables["surface_velocity"] a = legacy._device_boundaries[1].tables["surface_velocity"]
b = driven._level_boundaries[1].tables["surface_velocity"] b = driven._device_boundaries[1].tables["surface_velocity"]
diff = np.abs(a - b).max() diff = np.abs(a - b).max()
print(f"fine jet vs legacy velocity rows: max |dU| = {diff:.3e} (exact={diff == 0})") print(f"fine jet vs legacy velocity rows: max |dU| = {diff:.3e} (exact={diff == 0})")
self.assertLessEqual(diff, 4 * np.finfo(np.float32).eps * SPEED) self.assertLessEqual(diff, 4 * np.finfo(np.float32).eps * SPEED)
@@ -106,8 +106,8 @@ class JetStepperCPU(unittest.TestCase):
def test_blocked_amplitude_update_rewrites_only_velocity_column(self): def test_blocked_amplitude_update_rewrites_only_velocity_column(self):
backend, htod = mocked([circle()], [jet()], FINE_BLOCKS[1]) backend, htod = mocked([circle()], [jet()], FINE_BLOCKS[1])
before = backend._level_boundaries[1].tables["surface_velocity"].copy() before = backend._device_boundaries[1].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): with patch.object(cuda, "memcpy_htod", htod):
htod.reset_mock() htod.reset_mock()
backend.set_jet_amplitude(0, -SPEED, level=1) backend.set_jet_amplitude(0, -SPEED, level=1)
@@ -115,7 +115,7 @@ class JetStepperCPU(unittest.TestCase):
self.assertEqual(htod.call_count, 1) self.assertEqual(htod.call_count, 1)
self.assertIs(ptr, device) self.assertIs(ptr, device)
np.testing.assert_allclose(sent, -before, rtol=1e-6, atol=0) np.testing.assert_allclose(sent, -before, rtol=1e-6, atol=0)
np.testing.assert_array_equal(backend._level_boundaries[1].tables["surface_velocity"], sent) np.testing.assert_array_equal(backend._device_boundaries[1].tables["surface_velocity"], sent)
with self.assertRaises(ValueError): with self.assertRaises(ValueError):
backend.set_jet_amplitude(0, 1.5, level=1) backend.set_jet_amplitude(0, 1.5, level=1)
self.assertEqual(htod.call_count, 1) self.assertEqual(htod.call_count, 1)
@@ -133,7 +133,7 @@ class JetStepperCPU(unittest.TestCase):
with self.assertRaises(IndexError): with self.assertRaises(IndexError):
backend.set_jet_amplitude(1, -SPEED, level=1) backend.set_jet_amplitude(1, -SPEED, level=1)
self.assertEqual(backend._status, status) self.assertEqual(backend._status, status)
before = backend._level_boundaries[1].tables["surface_velocity"].copy() before = backend._device_boundaries[1].tables["surface_velocity"].copy()
boom = Mock(side_effect=RuntimeError("copy failed")) boom = Mock(side_effect=RuntimeError("copy failed"))
with patch.object(cuda, "memcpy_htod", boom), patch.object(cuda, "memcpy_htod_async", boom), \ with patch.object(cuda, "memcpy_htod", boom), patch.object(cuda, "memcpy_htod_async", boom), \
patch.object(cuda, "pagelocked_empty", side_effect=lambda shape, dtype: np.empty(shape, dtype)): patch.object(cuda, "pagelocked_empty", side_effect=lambda shape, dtype: np.empty(shape, dtype)):
@@ -141,7 +141,7 @@ class JetStepperCPU(unittest.TestCase):
backend.set_jet_amplitude(0, -SPEED, level=1) backend.set_jet_amplitude(0, -SPEED, level=1)
self.assertEqual(backend._status, "failed", fine_block) self.assertEqual(backend._status, "failed", fine_block)
self.assertEqual(backend.jets[1][0].amplitude, SPEED) self.assertEqual(backend.jets[1][0].amplitude, SPEED)
np.testing.assert_array_equal(backend._level_boundaries[1].tables["surface_velocity"], before) np.testing.assert_array_equal(backend._device_boundaries[1].tables["surface_velocity"], before)
with self.assertRaisesRegex(RuntimeError, "failed"): with self.assertRaisesRegex(RuntimeError, "failed"):
backend.step(1) backend.step(1)
with self.assertRaisesRegex(RuntimeError, "failed"): with self.assertRaisesRegex(RuntimeError, "failed"):
@@ -102,7 +102,7 @@ class Ledger:
def step(self): def step(self):
run=self.run;self.samples.clear();run.step() 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) self.test.assertEqual(len(self.samples),2**self.index)
expected=np.zeros_like(self.samples[0],dtype=np.float64) expected=np.zeros_like(self.samples[0],dtype=np.float64)
weights=np.array([h*h,h*h,h*h*h,h*h]) 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: level_objects=[(),(obj,)],level_jets=[(),(whole,)],observe_raw_exchange=True) as run:
run._initialize_raw(initial(run,rho,u=.013,v=0.)) run._initialize_raw(initial(run,rho,u=.013,v=0.))
expected=run._download_raw();ledger=Ledger(self,run,1) 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() for _ in range(4):ledger.step()
error=max(float(abs(a[:,leaf]-b[:,leaf]).max()) for a,b,leaf in error=max(float(abs(a[:,leaf]-b[:,leaf]).max()) for a,b,leaf in
zip(run._download_raw().values,expected.values,run.leaf)) 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(before,old._download_raw())
self.exact_raw(_StoredPopulations(expected.g,True),new._download_raw()) self.exact_raw(_StoredPopulations(expected.g,True),new._download_raw())
self.assertEqual(new.jets[2][0].amplitude,-.001) self.assertEqual(new.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() newledger=Ledger(self,new,2);ledger.step();newledger.step()
self.assertEqual((old.step_count,new.step_count),(count+1,count+1)) self.assertEqual((old.step_count,new.step_count),(count+1,count+1))
ledger.report('L2 source');newledger.report('L2 remapped') ledger.report('L2 source');newledger.report('L2 remapped')
+15 -6
View File
@@ -81,7 +81,7 @@ def public_jets(center, radius, scenario):
) )
def jet_interval(handle, scenario, begin, n, diameter, *, level=None): def jet_interval(handle, scenario, begin, n, diameter, level=None):
"""Apply a reversal only at a completed interval boundary, before advancing.""" """Apply a reversal only at a completed interval boundary, before advancing."""
if scenario != "none": if scenario != "none":
boundaries = [round(t * diameter / SPEED) for t in (250, 260, 270, 280, 290, 300, 305)] boundaries = [round(t * diameter / SPEED) for t in (250, 260, 270, 280, 290, 300, 305)]
@@ -92,9 +92,10 @@ def jet_interval(handle, scenario, begin, n, diameter, *, level=None):
if begin < switch: if begin < switch:
return min(n, switch - begin), None return min(n, switch - begin), None
if begin == switch: if begin == switch:
target = {} if level is None else dict(level=level) # FixedPatchStepper needs an explicit Jet level; Simulation has one level.
handle.set_jet_amplitude(0, -JET_PEAK, **target) kw = {} if level is None else {"level": level}
handle.set_jet_amplitude(1, JET_PEAK, **target) 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, dict(step=begin, tstar=300.0, amplitudes=[-JET_PEAK, JET_PEAK])
return n, None return n, None
@@ -257,7 +258,7 @@ def host_geometry(source, h, collision="SRT", jet_scenario="none", placement="wi
validate_channel_config, validate_channel_config,
validate_channel_transfers, 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.geometry.circle import CircleGeometry
from CelerisLab.body.objects import SimObject from CelerisLab.body.objects import SimObject
@@ -269,7 +270,15 @@ def host_geometry(source, h, collision="SRT", jet_scenario="none", placement="wi
pair = FixedPatch2D.from_config(cfg, PLACEMENTS[placement]) pair = FixedPatch2D.from_config(cfg, PLACEMENTS[placement])
layout = BlockLayout2D.build(cfg.ny, cfg.nx, *BASE_BLOCK) layout = BlockLayout2D.build(cfg.ny, cfg.nx, *BASE_BLOCK)
validate_channel_config(cfg, [BASE_BLOCK, FINE_BLOCK]) validate_channel_config(cfg, [BASE_BLOCK, FINE_BLOCK])
targets = blocked_base_tables(pair, layout, np.argsort(pair.coalesce_coarse))[1] # Select by the imported source's own API, not by scenario: the pinned historical
# source (verify_source) defines blocked_base_tables, the current checkout (host
# tests, public Jet scenarios) blocked_coarse_tables. Exactly one must exist.
names = [n for n in ("blocked_base_tables", "blocked_coarse_tables")
if hasattr(refinement_tables, n)]
if len(names) != 1:
raise RuntimeError(f"Unexpected coarse transfer table API: {names}")
coarse_tables = getattr(refinement_tables, names[0])
targets = coarse_tables(pair, layout, np.argsort(pair.coalesce_coarse))[1]
validate_channel_transfers(pair, layout, targets) validate_channel_transfers(pair, layout, targets)
BlockLayout2D.build(*pair.fine_shape, *FINE_BLOCK, periodic=(False, False)) BlockLayout2D.build(*pair.fine_shape, *FINE_BLOCK, periodic=(False, False))
body = SimObject( body = SimObject(