Compare commits

..
Author SHA1 Message Date
ClaudeandClaude Opus 5.5 6d046f91c1 docs: sync src_coherence with the landed slices
Stage history and verification, historical label on the migration map,
host-preparation/device order, pair accessors, per-level state split into
fixed plan structure / mutable boundary data / runtime _Level, the two
block_collide families, and the deferred items with their evidence needs.
Doc-only.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-28 12:13:06 +00:00
ClaudeandClaude Opus 5.5 3361e97d5f docs: src_coherence migration map for the combined _boundaries record
Final private names (_device_boundaries/_level_jets/_exchange_*/_level_blocks
-> _boundaries record fields / layout.block), public boundaries = host view,
and the plan-versus-device (_Level) lifecycle. Doc-only.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-28 12:09:26 +00:00
ClaudeandClaude Opus 5.5 442979b4db FixedPatchStepper: exchange metadata on the boundary record; _boundaries name
- _exchange_spacing/_exchange_centers (two more per-level lists, only read
  for levels that have a boundary) become _BlockedBoundary.spacing/centers.
- _level_blocks dropped: remap_finest reads each layout's own .block.
- _device_boundaries -> _boundaries: the record holds host geometry, tables,
  Jets and device handles, so 'device' was misleading (the old host-only
  _boundaries list no longer exists). Tests migrated mechanically.
Remaining per-level lists (_layouts, _level_cfgs, _leaf, _state, _boundaries)
are the immutable construction plan, built before the memory preflight and
kept after close; _levels holds device handles only and is cleared on close.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-28 11:41:13 +00:00
ClaudeandClaude Opus 5.5 c5a3cdb757 test_refinement_mrt: check force reads in executable text, not comments
The block_collision.cuh header comment names the absent d_params.fx/fy
read, which tripped the source-text assertion. Strip // comments first
and also require the real d_params.omega read. Test-only; no kernel source
change.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-28 11:27:35 +00:00
ClaudeandClaude Opus 5.5 a4fc46d849 block_collide: one unforced uniform-equivalent branch; drop collide_mrt_unforced
The blocked collision had three arms (MRT wrapper, channel SRT, level SRT)
and a one-caller MRT wrapper in collision_mrt.cuh. Now block_collide computes
rho/u once and has two compile-time families:
- uniform-equivalent (MRT any level, SRT channel root): collide_dispatch
  arithmetic with zero force, no LES: one clamped omega, collide_mrt_no_force
  or collide_srt with its zero Guo term (unchanged);
- level SRT (BLOCK_ROOT_BC 0/1): unclamped, collide_srt_no_force, shifted
  momentum feq exception kept visible.
No runtime d_params.fx/fy is read. Uniform kernels unchanged (only the
unused inline wrapper is removed from collision_mrt.cuh). Source-text test
updated to the new single call sites.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-28 10:55:02 +00:00
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
25 changed files with 269 additions and 192 deletions
+49 -13
View File
@@ -1,8 +1,11 @@
# src coherence stage (from main 914906e)
Owner direction 2026-10-28: unified style over historical compatibility, inside `src`.
No numerical change is intended; kernel text was moved verbatim (PTX/bit parity is
pending a GPU run).
No numerical change is intended. Stages after 5f69eba: 65339ab (one boundary/Jet record,
host preparation helpers, pair accessors), c5a3cdb (`block_collide` families), 442979b
(exchange metadata on the boundary record). Each was checked by parent/candidate
config.h/PTX comparison and host/GPU equality runs; c5a3cdb changed the PTX of the two
channel-root SRT modules only (omega load/clamp moved, outputs exact).
## Conventions (defined once, in `lbm.hierarchy.Hierarchy`)
@@ -28,20 +31,25 @@ pending a GPU run).
## FixedPatchStepper migration map
Historical: left column = names at or after 914906e that no longer exist; right column =
the current name.
| removed / changed | replacement |
|---|---|
| `patch`, `child_patch` | `hierarchy.pairs[0]`, `hierarchy.pairs[1]` |
| `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) | `_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` |
| `_boundaries[l]` (dict from L1, level keys), `_level_boundaries` (b4eacf7), `_device_boundaries` (65339ab) | `_boundaries[l]`: the one combined per-level `_BlockedBoundary` record (list from L0, `None` at L0 unless single-level static circles): host geometry `.host`, `.tables`, `.jets`, `.jet_cache`, exchange `.spacing`/`.centers`, device handles |
| `_fine_tables`, `_actuated`, `_jet_cache` | `_boundaries[1].tables` / `.actuated` / `.jet_cache` |
| `_level_jets` (private) | `_boundaries[l].jets` (public `jets[l]` is derived) |
| `_exchange_spacing[l]`, `_exchange_centers[l]` (private) | `_boundaries[l].spacing` (h_l, root units) / `.centers` (root edge frame); only levels with a boundary report exchange |
| `_level_blocks` (private) | `_layouts[l].block` |
| `_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`) |
| `_boundary_hosts` (private) | public `boundaries[l]` = `_boundaries[l].host` (host geometry only; the private record is the combined one) |
| `remap_child(rects, child_objects=...)` | `remap_finest(rects, objects=...)` |
| `set_jet_amplitude(i, a)` (implicit level 1) | `set_jet_amplitude(i, a, level=l)` (required) |
| `read_raw_exchange(b)` / `read_exchange_window(b)` (implicit level 1) | `read_raw_exchange(b, level=l)` / `read_exchange_window(b, level=l)` (required) |
@@ -56,13 +64,35 @@ the replaced level relative to its parent).
## Responsibilities
- Construction: `_check_levels` (inputs), `_kernel_configs` (per-level kernel configs),
`Hierarchy.edge_frame` (geometry frames), `_build_level` / `_build_pair` / `_upload_boundary`
(device). Observation reads share `_observed_boundary` + `_summarize`.
- Construction: `_check_levels` (inputs); host preparation, numpy only, before any device
work: `_level_layouts`, `_kernel_configs` (per-level kernel configs), `_host_boundaries`,
`_prepare_blocked_boundary`, `_fluid_leaves`; `Hierarchy.edge_frame` (geometry frames);
then `_planned_bytes()` (memory preflight) and `_build_level` / `_build_pair` /
`_upload_boundary` (device). A construction failure closes the partial instance.
Observation reads share `_observed_boundary` + `_summarize`.
- Pair accessors: `rectangles` (a scalar FixedPatch2D gives `(bounds,)`) and
`fine_level_offset` (`(0, 0)` for a set) on both pair types; the types and their
storage frames stay distinct.
- Per-level state, by lifecycle:
- Fixed plan structure (built once in `__init__`, read by `_planned_bytes()` before any
allocation, readable after `close()` through the public views): `_layouts`,
`_level_cfgs`, `_leaf`, `_state`, and in each `_boundaries[l]` record the host geometry
`.host`, link tables, `.jet_cache`, `.spacing`, `.centers`.
- Mutable boundary data in the same record: `.jets` and `.tables["surface_velocity"]`,
changed only by `set_jet_amplitude` (host prepare -> device write -> commit; an uncertain
write fails the instance); device handles (`.flags`, `.device`, `.kernel`, ledgers) are
set by `_upload_boundary` and cleared by `close()`.
- Runtime: `_levels` (`_Level`) holds compiled modules, kernels, the two population
buffers (swapped per step) and transfer/kind handles. It is created only after the
memory preflight and cleared on `close()`; `.layout`/`.boundary` reference the plan for
`_hop`. Folding the plan into `_Level` would mix both lifecycles in one record.
- Blocked kernels: `operators/block_collision.cuh block_collide` is the one collision extension
point of `BlockStep` and `BlockChannelStep`. MRT shares one path; SRT keeps its two explicit
arithmetic variants (level: unclamped, no-force, shifted direct-momentum feq; channel: uniform
OneStep arithmetic with omega clamp), selected by `BLOCK_ROOT_BC`.
point of `BlockStep` and `BlockChannelStep`. It computes rho/u once and reads no runtime
force (`d_params.fx/fy`). Two compile-time families: uniform-equivalent (MRT on any level,
SRT on the channel root `BLOCK_ROOT_BC == 2`): clamped omega, `collide_mrt_no_force` or
`collide_srt` with its zero Guo term; and level SRT (`BLOCK_ROOT_BC` 0/1): unclamped omega,
`collide_srt_no_force`, with the shifted direct-momentum feq under `USE_DDF_SHIFTING`.
The uniform path keeps `collide_dispatch` (runtime force, LES, TRT).
- `explosion="projected_linear"` stays an optional capability; `"constant"` stays default.
## Not done / open
@@ -72,6 +102,12 @@ the replaced level relative to its parent).
`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
(uniform node frame, body sync, forcing, checkpoint vs. fixed blocked hierarchy, edge
frame, remap-as-new-instance). A facade would only forward calls; no concrete need found.
- Deferred, each needing its own evidence first: sharing the channel-root cell/donor selector
with OneStep (the closure functions are already shared; a shared selector changes uniform
signatures and needs a uniform PTX gate); merging FixedPatch2D into FixedPatchSet2D
(different origin, ghost band, layout size and addressing; needs a geometry/parity table and
would need tiles that are not grid-aligned); one `clamp_omega` helper for uniform and
blocked paths (it lives in the hash-pinned LES header `turbulence_smag.cuh`).
- Tests, examples and validation callers are migrated separately (pinned runners stay frozen).
+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
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 fixed refinement patch (refinement.py), whose semantics are unchanged.
+1 -1
View File
@@ -122,7 +122,7 @@ def leaf_partition(h):
for level, kind in enumerate(h.kinds):
if level:
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
y, x = np.nonzero(kind == ACTIVE)
np.add.at(base, ((y + oy) >> level, (x + ox) >> level), 0.25 ** level)
@@ -1,10 +1,12 @@
// CelerisLab – operators/block_collision.cuh
// The one collision extension point of the blocked steps (BlockStep, BlockChannelStep).
// Compile-time selection: MRT shares one path; SRT keeps two explicit arithmetic
// variants, unchanged from their original kernels:
// level (BLOCK_ROOT_BC 0/1): unclamped omega, no forcing term; shifted storage
// builds feq from the retained drho and direct momentum (numerical variant);
// channel (BLOCK_ROOT_BC 2): matches uniform OneStep SRT arithmetic (omega clamp,
// zero Guo forcing term, compute_feq).
// Unforced by construction: FixedPatchStepper supplies no forcing, so no runtime
// d_params.fx/fy branch is read. Two compile-time arithmetic families:
// uniform-equivalent (MRT on every level; SRT on the channel root, BLOCK_ROOT_BC 2):
// collide_dispatch arithmetic with zero force and no LES: clamped omega;
// SRT keeps the zero Guo forcing term of collide_srt;
// level SRT (BLOCK_ROOT_BC 0/1): unclamped omega, no forcing term; shifted storage
// builds feq from the retained drho and direct momentum (numerical variant).
#ifndef CELERIS_OPERATORS_BLOCK_COLLISION_CUH
#define CELERIS_OPERATORS_BLOCK_COLLISION_CUH
@@ -16,17 +18,19 @@
#endif
__device__ __forceinline__ void block_collide(float* f) {
#if COLLISION_MODEL == 2
collide_mrt_unforced(f);
#elif BLOCK_ROOT_BC == 2
float rho, ux, uy, feq[NQ], forcing[NQ] = {};
float rho, ux, uy;
compute_rho_u(f, rho, ux, uy);
compute_feq(rho, ux, uy, feq);
#if COLLISION_MODEL == 2 || BLOCK_ROOT_BC == 2
const float omega = fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, d_params.omega));
#if COLLISION_MODEL == 2
collide_mrt_no_force(f, rho, ux, uy, omega);
#else
float feq[NQ], forcing[NQ] = {};
compute_feq(rho, ux, uy, feq);
collide_srt(f, feq, forcing, omega);
#endif
#else
float rho, ux, uy, feq[NQ];
compute_rho_u(f, rho, ux, uy);
float feq[NQ];
#if USE_DDF_SHIFTING
// Preserve the density perturbation before adding the unit reference.
const float drho = f[0]+f[1]+f[2]+f[3]+f[4]+f[5]+f[6]+f[7]+f[8];
@@ -114,16 +114,6 @@ __device__ __forceinline__ void collide_mrt_no_force(float* __restrict__ g,
collide_mrt(g, rho, ux, uy, Fin, omega);
}
// Unforced MRT with the clamped device omega, as in collide_dispatch with zero
// forcing. Shared by the blocked, channel-root and dense-patch collisions.
__device__ __forceinline__ void collide_mrt_unforced(float* __restrict__ g)
{
float rho, ux, uy;
compute_rho_u(g, rho, ux, uy);
collide_mrt_no_force(g, rho, ux, uy,
fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, d_params.omega)));
}
#elif LATTICE_MODEL == LATTICE_D3Q19
__device__ __forceinline__ void collide_mrt(float* __restrict__ g,
+12
View File
@@ -120,6 +120,11 @@ class FixedPatch2D:
(an index offset, not an (x, y) origin)."""
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
def covered(self) -> np.ndarray:
"""Dense coarse mask of cells replaced by the patch (not leaves)."""
@@ -225,11 +230,18 @@ class FixedPatchSet2D:
def bounds(self) -> tuple:
return tuple(p.bounds for p in self.patches)
rectangles = bounds
@property
def fine_origin(self) -> tuple[float, float]:
"""Lower edge of the level-wide fine grid in coarse physical units."""
return 0.0, 0.0
@property
def fine_level_offset(self) -> tuple[int, int]:
"""The level-wide grid itself: zero (y, x) index offset."""
return 0, 0
@property
def covered(self) -> np.ndarray:
mask = np.zeros(self.coarse_shape, bool)
+109 -84
View File
@@ -111,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")
# 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:
"""D2Q9 SRT/MRT FP32 fixed 2:1 refinement hierarchy, built from level lists.
@@ -246,40 +297,15 @@ class FixedPatchStepper:
self.hierarchy = build_hierarchy(cfg, refinements, wall_aware=self._wall_aware)
pairs = self.hierarchy.pairs
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 rects):
if any(r[1] < 4 or r[3] > cfg.ny - 4 for r in pairs[0].rectangles):
raise ValueError("slot_reservoir requires y0 >= 4 and y1 <= ny-4")
self._layouts = []
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)
self._layouts = _level_layouts(cfg, pairs, level_blocks, self._wall_aware)
if self._channel:
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).
self._boundaries = [None]
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)
hosts = _host_boundaries(self.hierarchy, self._level_cfgs, level_objects, level_jets,
self._wall_aware, self._ddf_shifting)
if not isinstance(observe_raw_exchange, bool):
raise ValueError("observe_raw_exchange must be bool")
self._observe_raw_exchange = observe_raw_exchange
@@ -291,32 +317,25 @@ class FixedPatchStepper:
self._window_valid = False
self._window_start = None
self._exchange_valid = False
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)
for objects in level_objects)
self._level_jets = level_jets
root = None
if depth == 1 and level_objects[0]:
tables = static_circle_tables(self._layouts[0], level_objects[0])
root = _BlockedBoundary(None, tables, (), False)
self._device_boundaries = [root] + [
_prepare_blocked_boundary(self._boundaries[level], self._layouts[level],
# The one boundary/Jet record per level; host geometry and Jets are views of it.
self._boundaries = [root] + [
_prepare_blocked_boundary(hosts[level], self._layouts[level],
self._level_cfgs[level], level_jets[level],
_shifted_blocked=self._ddf_shifting)
for level in range(1, depth)]
self._leaf = [kind == ACTIVE for kind in self.hierarchy.kinds]
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
for level, b in enumerate(self._boundaries):
if b is not None: # exchange scaling (dt_l = h_l) and root-edge-frame moment centers
b.spacing = self.hierarchy.edge_frame(level)[1]
b.centers = tuple(tuple(o.geometry.center) for o in level_objects[level])
self._leaf = _fluid_leaves(self.hierarchy, self._boundaries)
self._state = [mask.copy() for mask in self._leaf]
if self._channel:
self._leaf[0] &= channel_masks((cfg.ny, cfg.nx))[0]
# Private constructor inputs for remap_finest.
self._root_cfg = replace(cfg)
self._level_blocks = level_blocks
self._root_cfg = replace(cfg) # private constructor input for remap_finest
self.ctx = context
self._lock = RLock()
self._status = "new"
@@ -371,7 +390,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._device_boundaries[index]
b = self._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)
@@ -394,7 +413,7 @@ 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._device_boundaries[index]
level.boundary = self._boundaries[index]
level.halo_copy = module.get_function("BlockHaloCopy")
if len(layout.halo_dst):
level.halo_dst = self._upload(layout.halo_dst)
@@ -556,19 +575,17 @@ class FixedPatchStepper:
pairs = self.hierarchy.pairs
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._device_boundaries[:2])):
or any(b is not None for b in self._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._root_cfg
objects = tuple(objects)
actual = self._device_boundaries[2]
actual = self._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,
_shifted_blocked=self._ddf_shifting)
if (not same_physical_boundary(old.boundary, self._boundaries[2])
or (self._boundaries[2] is not None and actual is None)):
if not same_physical_boundary(old.boundary, None if actual is None else actual.host):
raise ValueError("objects must reproduce the original root-edge-frame level-2 boundary")
if actual is not None:
velocity = actual.tables.get("surface_velocity")
@@ -585,7 +602,7 @@ class FixedPatchStepper:
or not np.array_equal(old.surface_velocity, new.surface_velocity)):
raise ValueError("remap requires the identical pinned physical boundary and BFL support")
# Check sparse block alignment/launch size before any device read/allocation.
layout = new.child.fine_layout(self._level_blocks[2])
layout = new.child.fine_layout(self._layouts[2].block)
if layout.nblocks > MAX_GRID_DIM or layout.local_shape[0] > MAX_GRID_DIM:
raise ValueError("level 2 blocks exceed launch limits")
with self.ctx.activate():
@@ -603,7 +620,7 @@ class FixedPatchStepper:
phase="postcollision_g", **options)
candidate = type(self)(
cfg, context=self.ctx, refinements=[pairs[0].bounds, new.child.bounds],
level_blocks=list(self._level_blocks),
level_blocks=[layout.block for layout in self._layouts],
level_objects=[(), (), objects], level_jets=[(), (), jets],
observe_raw_exchange=self._observe_raw_exchange,
observe_window_exchange=self._observe_window, explosion=self._explosion)
@@ -622,7 +639,7 @@ class FixedPatchStepper:
return candidate
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
instance failed (as step/initialize do); host state commits only after
@@ -633,27 +650,24 @@ 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._device_boundaries[level]
jets = list(self._level_jets[level])
boundary = self._boundaries[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
velocity = None
if boundary is not None and boundary.jet_cache is not None:
shape, jet_id = boundary.jet_cache
velocity = boundary.tables["surface_velocity"].copy()
driven = jet_id >= 0
velocity[driven] = jet_link_velocity(jets, shape, jet_id)[driven]
shape, jet_id = boundary.jet_cache
velocity = boundary.tables["surface_velocity"].copy()
driven = jet_id >= 0
velocity[driven] = jet_link_velocity(jets, shape, jet_id)[driven]
# Device action: an uncertain write fails the instance closed.
try:
if velocity is not None:
self._stream.synchronize() # no in-flight step reads the column
cuda.memcpy_htod(boundary.device["surface_velocity"], velocity)
self._stream.synchronize() # no in-flight step reads the column
cuda.memcpy_htod(boundary.device["surface_velocity"], velocity)
except BaseException:
self._status = "failed"
raise
if velocity is not None:
boundary.tables["surface_velocity"] = velocity
if boundary is not None:
boundary.jets = tuple(jets)
self._level_jets[level] = tuple(jets)
# Commit: the one boundary record holds both the column and the Jets.
boundary.tables["surface_velocity"] = velocity
boundary.jets = tuple(jets)
def _blocked_boundary(self, level, buffer, **kw):
b = level.boundary
@@ -696,7 +710,7 @@ class FixedPatchStepper:
b = level.boundary
if b is not None and b.last_raw is not None:
count = len(b.tables["cell"])
h = self._exchange_spacing[index] # acoustic scaling: dt_l = h_l
h = b.spacing # acoustic scaling: dt_l = h_l
b.accumulate(b.last_raw, b.interval_integral, np.uint32(count),
np.float64(h), np.float64(h),
grid=((count + self._block[0] - 1) // self._block[0], 1), **kw)
@@ -740,7 +754,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._device_boundaries:
for b in self._boundaries:
if b is not None:
cuda.memset_d8_async(b.interval_integral, 0,
len(b.tables["cell"]) * 32, self._stream)
@@ -763,7 +777,7 @@ class FixedPatchStepper:
raise RuntimeError("window exchange observation is disabled")
self._window_valid = False
try:
for b in self._device_boundaries:
for b in self._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:
@@ -776,7 +790,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._device_boundaries[level]
b = self._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")
@@ -784,9 +798,9 @@ class FixedPatchStepper:
def _summarize(self, b, rows, level, body_id, window, **kw):
return summarize_exchange(rows, b.tables["body_id"], b.tables["aperture"],
body_id=body_id, level=level, h=self._exchange_spacing[level],
body_id=body_id, level=level, h=b.spacing,
time_end=self._step_count, window=window,
reference_center=self._exchange_centers[level][body_id], **kw)
reference_center=b.centers[body_id], **kw)
def read_exchange_window(self, body_id, *, level):
"""Non-destructive read of [window start, current step count] in base-lattice units."""
@@ -873,7 +887,7 @@ class FixedPatchStepper:
self._stream.synchronize()
except Exception as exc:
error = exc
for b in getattr(self, "_device_boundaries", ()):
for b in getattr(self, "_boundaries", ()):
if b is not None:
b.flags = b.device = b.kernel = None
b.last_raw = b.interval_integral = b.accumulate = None
@@ -901,9 +915,20 @@ class FixedPatchStepper:
raise
# Read-only per-level state, indexed by level (0 = root). Pairs: hierarchy.pairs[l-1].
layouts = property(lambda self: tuple(self._layouts))
boundaries = property(lambda self: tuple(self._boundaries),
doc="Host boundary geometry per level (None: no bodies/Jets).")
jets = property(lambda self: tuple(self._level_jets))
leaf = property(lambda self: tuple(m.copy() for m in self._leaf),
doc="Per-level (ny, nx) masks of public (non-NaN) download cells.")
@property
def layouts(self):
return tuple(self._layouts)
@property
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._boundaries)
@property
def jets(self):
return tuple(() if b is None else b.jets for b in self._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)
+2
View File
@@ -130,6 +130,8 @@ class _BlockedBoundary:
interval_integral: object = None
accumulate: object = None
window_integral: object = None # optional observation-window FP64 ledger
spacing: float = None # level spacing h_l in root units (exchange scaling)
centers: tuple = () # body reference centers in the root edge frame
def _prepare_blocked_boundary(boundary, layout, cfg, jets, *, _shifted_blocked=False):
+5 -5
View File
@@ -102,7 +102,7 @@ def test_simulation_before_initialize_commits_host_only(monkeypatch):
def stepper_state(run):
boundary = run._device_boundaries[2]
boundary = run._boundaries[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):
run.set_jet_amplitude(*args, **kw)
module.cuda.memcpy_htod.assert_not_called()
np.testing.assert_array_equal(run._device_boundaries[2].tables["surface_velocity"], velocity)
assert (run._device_boundaries[2].jets, run.jets[2]) == (jets, level_jets)
np.testing.assert_array_equal(run._boundaries[2].tables["surface_velocity"], velocity)
assert (run._boundaries[2].jets, run.jets[2]) == (jets, level_jets)
assert run._status == "ready"
run.set_jet_amplitude(0, .017, level=2)
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")
with pytest.raises(RuntimeError, match="upload failed"):
run.set_jet_amplitude(0, .017, level=2)
np.testing.assert_array_equal(run._device_boundaries[2].tables["surface_velocity"], velocity)
assert (run._device_boundaries[2].jets, run.jets[2]) == (jets, level_jets)
np.testing.assert_array_equal(run._boundaries[2].tables["surface_velocity"], velocity)
assert (run._boundaries[2].jets, run.jets[2]) == (jets, level_jets)
assert run._status == "failed"
module.cuda.memcpy_htod.side_effect = None
run._stream.synchronize.side_effect = None
+4 -4
View File
@@ -96,11 +96,11 @@ class MultiplePinnedRemap(unittest.TestCase):
before = run.download()
clone = run.remap_finest((PINS+MOVED+ADDED)[::-1], objects=bodies)
try:
self.assertEqual(tuple(j.amplitude for j in clone._device_boundaries[2].jets), (.017, -.011))
np.testing.assert_array_equal(clone._device_boundaries[2].tables['surface_velocity'],
run._device_boundaries[2].tables['surface_velocity'])
self.assertEqual(tuple(j.amplitude for j in clone._boundaries[2].jets), (.017, -.011))
np.testing.assert_array_equal(clone._boundaries[2].tables['surface_velocity'],
run._boundaries[2].tables['surface_velocity'])
clone.set_jet_amplitude(1, .005, level=2)
self.assertEqual(run._device_boundaries[2].jets[1].amplitude, -.011)
self.assertEqual(run._boundaries[2].jets[1].amplitude, -.011)
clone.step()
self.assertEqual((run.step_count, clone.step_count), (7, 8))
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)
before=run.download()
clone=run.remap_finest(B,objects=objects)
self.assertEqual(clone._device_boundaries[2].jets[0].amplitude,.017)
np.testing.assert_array_equal(clone._device_boundaries[2].tables['surface_velocity'],
run._device_boundaries[2].tables['surface_velocity'])
self.assertEqual(clone._boundaries[2].jets[0].amplitude,.017)
np.testing.assert_array_equal(clone._boundaries[2].tables['surface_velocity'],
run._boundaries[2].tables['surface_velocity'])
self.assertTrue(same_physical_boundary(run.boundaries[2],clone.boundaries[2]))
clone.step()
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):
run,module,allocations,storage,objects=ready(jets=())
clone=run.remap_finest(A,objects=objects)
self.assertFalse(clone._device_boundaries[2].actuated)
self.assertNotIn('surface_velocity',clone._device_boundaries[2].tables)
self.assertFalse(clone._boundaries[2].actuated)
self.assertNotIn('surface_velocity',clone._boundaries[2].tables)
for a,b in zip(run.download(),clone.download()): np.testing.assert_array_equal(a,b)
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):
run, host, calls, arrays = backend()
b = run._device_boundaries[3]
b = run._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._device_boundaries[3]
b = run._boundaries[3]
self.assertIs(b.last_raw, b.interval_integral)
self.assertIs(b.interval_integral, None)
run.step()
+14 -10
View File
@@ -1,5 +1,6 @@
"""Host guards and source contract for MRT reuse in refined/blocked stepping (no CUDA)."""
from pathlib import Path
import re
import unittest
from CelerisLab.config import LBMConfig
@@ -32,22 +33,25 @@ class RefinementMRTHost(unittest.TestCase):
"collide_srt_no_force(f, feq, d_params.omega);"):
self.assertIn(text, block_collision)
self.assertIn("collide_srt(f, feq, forcing, omega);", block_collision)
self.assertIn("#elif BLOCK_ROOT_BC == 2", block_collision)
self.assertIn("#if COLLISION_MODEL == 2 || BLOCK_ROOT_BC == 2", block_collision)
self.assertEqual(channel.count("block_collide(f);"), 1)
self.assertNotIn("collide_mrt_unforced(f);", channel)
# Executable text only: the header comment explains the absent force read.
code = re.sub(r"//[^\n]*", "", channel + "\n" + block_collision)
self.assertIn("d_params.omega", code)
for text in ("collide_mrt_unforced", "collide_dispatch(", "d_params.fx", "d_params.fy"):
self.assertNotIn(text, code)
for source in (blocks,):
self.assertIn("(COLLISION_MODEL != 0 && COLLISION_MODEL != 2)", source)
self.assertIn("USE_LES", source)
self.assertIn("STORE_PRECISION != 0", source)
self.assertIn('#if COLLISION_MODEL == 2\n#include "operators/collision_mrt.cuh"\n#endif', source)
# One shared unforced MRT body; each blocked/patch path calls it exactly once.
mrt = (KERNELS / "operators/collision_mrt.cuh").read_text()
self.assertEqual(mrt.count("void collide_mrt_unforced("), 1)
for source in (blocks + "\n" + block_collision, channel + "\n" + block_collision):
self.assertEqual(source.count("collide_mrt_unforced(f);"), 1)
self.assertNotIn("collide_mrt_no_force(", source)
self.assertNotIn("compute_rho_u(f, rho, ux, uy);\n collide_mrt", source)
self.assertNotRegex(source, r"compute_rho_u\(f, rho, ux, uy\);\s*collide_mrt")
# One unforced MRT call, in block_collide only; the steps never collide directly.
self.assertEqual(block_collision.count("collide_mrt_no_force(f, rho, ux, uy, omega);"), 1)
self.assertEqual(block_collision.count("compute_rho_u(f, rho, ux, uy);"), 1)
self.assertEqual(block_collision.count("fmaxf(OMEGA_COLLISION_MIN, d_params.omega)"), 1)
for source in (blocks, channel):
self.assertNotIn("collide_mrt", source)
self.assertNotIn("collide_srt", source)
if __name__ == "__main__": unittest.main()
+5 -3
View File
@@ -1,5 +1,6 @@
"""CPU mocks execute production remap lifecycle; no GPU numerical claim."""
from dataclasses import replace
from types import SimpleNamespace
import unittest
from unittest.mock import patch
@@ -69,10 +70,11 @@ class RuntimeChildRemap(unittest.TestCase):
run._status = 'ready'
for bad in ([], [(40,16,48,40)], [(16,16,32,24),(16,27,32,55)]):
with self.assertRaises(ValueError): run.remap_finest(bad)
hosts, pairs = run._boundaries, run.hierarchy.pairs
records, pairs = run._boundaries, run.hierarchy.pairs
for attribute, value in (('_wall_aware', True),
('_boundaries', [*hosts[:2], object()]), # child boundary
('_boundaries', [hosts[0], object(), hosts[2]]), # fine boundary
# level-2 record whose host geometry differs from objects=()
('_boundaries', [*records[:2], SimpleNamespace(host=object(), jets=())]),
('_boundaries', [records[0], object(), records[2]]), # level-1 body
('hierarchy', replace(run.hierarchy, pairs=(pairs[0], None)))):
with patch.object(run, attribute, value):
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):
run,module,allocations,storage,objects=ready();clone=None
try:
boundary=run._device_boundaries[2]
boundary=run._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_finest(B,objects=objects)
self.assertEqual(clone.jets[2][0].amplitude,-jet().amplitude)
np.testing.assert_array_equal(clone._device_boundaries[2].tables['surface_velocity'],
np.testing.assert_array_equal(clone._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._device_boundaries[2].tables['surface_velocity'].copy()
before=run._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._device_boundaries[2].tables['surface_velocity'])
np.testing.assert_array_equal(before,run._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._device_boundaries[0].tables["solid"]
solid = run._boundaries[0].tables["solid"]
self.assertTrue(solid.any())
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)
clone=run.remap_finest(((16,16,24,40),(24,16,32,40),A[1]),objects=objects)
try:
self.assertEqual(clone._device_boundaries[2].jets[0].amplitude,.017)
self.assertEqual(clone._boundaries[2].jets[0].amplitude,.017)
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)
finally:clone.close()
@@ -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._device_boundaries[0].tables
tables = backend._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._device_boundaries[1]
b = backend._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._device_boundaries[1].window_integral)
self.assertIsNone(backend._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._device_boundaries[1]
b = windowed._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._device_boundaries[1]
b = windowed._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._device_boundaries[1].window_integral)
cuda.memcpy_dtoh(rows[1], split._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, level=1)
@@ -57,7 +57,7 @@ def make(context):
def topology(run, rectangles, objects):
return RemapTopology.from_config(
CFG, run.hierarchy.pairs[0], rectangles, child_objects=objects,
child_jets=run._device_boundaries[2].jets,
child_jets=run._boundaries[2].jets,
)
@@ -74,7 +74,7 @@ class MultiplePinnedRemapGPU(unittest.TestCase):
def change_amplitudes(self, run, amplitudes):
# Each jet selects a different body. Updating one must preserve the
# other body's current per-link velocities, including its driven links.
boundary = run._device_boundaries[2]
boundary = run._boundaries[2]
body_ids = run.boundaries[2].curved_host[3]
for index, amplitude in enumerate(amplitudes):
before = boundary.tables['surface_velocity'].copy()
@@ -91,7 +91,7 @@ class MultiplePinnedRemapGPU(unittest.TestCase):
self.assertTrue(same_physical_boundary(new.boundaries[2], mesh.boundary))
# Compare physical cut-link rows, never layout-dependent blocked addresses.
for run in (old, new):
boundary = run._device_boundaries[2]
boundary = run._boundaries[2]
self.assertEqual(tuple(jet.amplitude for jet in boundary.jets), amplitudes)
exact(self, boundary.tables['surface_velocity'], mesh.surface_velocity,
'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)
for a,b in zip(after.boundaries[2].curved_host,before.boundaries[2].curved_host):
exact(self,a,b,'physical cut links')
exact(self,after._device_boundaries[2].tables['surface_velocity'],
before._device_boundaries[2].tables['surface_velocity'],'current per-link velocities')
self.assertEqual(after._device_boundaries[2].jets[0].amplitude,.017)
exact(self,after._boundaries[2].tables['surface_velocity'],
before._boundaries[2].tables['surface_velocity'],'current per-link velocities')
self.assertEqual(after._boundaries[2].jets[0].amplitude,.017)
def test_noop_next_step_bitwise(self):
with CudaContext(0) as ctx, ExitStack() as owned:
@@ -63,7 +63,7 @@ class PinnedBodyRemapGPU(unittest.TestCase):
current.set_jet_amplitude(0,.017,level=2)
for rects in (B,A):
previous=current;before=previous.download();count=previous.step_count
jets=previous._device_boundaries[2].jets
jets=previous._boundaries[2].jets
old=RemapTopology.from_config(CFG,previous.hierarchy.pairs[0],previous.hierarchy.pairs[1].bounds,
child_objects=objects,child_jets=jets)
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]
measured.initialize(*arrays)
plain.initialize(*arrays)
b = measured._device_boundaries[1]
b = measured._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._device_boundaries[3]
b = run._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._device_boundaries[1]
b = amr._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._device_boundaries[0]
root = uniform._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']))
@@ -187,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._device_boundaries[1].flags, backend._levels[1].kind
velocity = backend._device_boundaries[1].device.get("surface_velocity")
flags, kind = backend._boundaries[1].flags, backend._levels[1].kind
velocity = backend._boundaries[1].device.get("surface_velocity")
self.assertEqual(velocity is None, spec is None)
backend.step(1)
backend.close()
@@ -213,8 +213,8 @@ class BlockedFineCPU(unittest.TestCase):
def test_actuation_without_velocity_buffer_fails_before_bfl(self):
for speed in (0.0, .002): # any spec, even zero amplitude
backend, calls = self.mocked_backend([circle(actuation(speed))])
self.assertIsNotNone(backend._device_boundaries[1].tables["surface_velocity"])
del backend._device_boundaries[1].device["surface_velocity"]
self.assertIsNotNone(backend._boundaries[1].tables["surface_velocity"])
del backend._boundaries[1].device["surface_velocity"]
with self.assertRaisesRegex(RuntimeError, "Surface velocity device buffer is missing"):
backend.step(1)
names = [name for name, _ in calls]
@@ -62,8 +62,8 @@ class JetStepperCPU(unittest.TestCase):
C = (16.0, 12.0)
whole = Jet(Sector(C, .1, math.pi, 3.), Radial(C, inward=False), amplitude=0.)
backend, _ = mocked([circle()], [whole], FINE_BLOCKS[1])
_, jet_id = backend._device_boundaries[1].jet_cache
self.assertGreater(len(np.unique(backend._device_boundaries[1].tables["block"][jet_id == 0])), 1)
_, jet_id = backend._boundaries[1].jet_cache
self.assertGreater(len(np.unique(backend._boundaries[1].tables["block"][jet_id == 0])), 1)
backend.close()
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)
whole = Jet(Sector(C, .1, math.pi, 3.), Radial(C, inward=False), amplitude=0.)
backend, _ = mocked([source], [whole], FINE_BLOCKS[1])
_, jet_id = backend._device_boundaries[1].jet_cache
_, jet_id = backend._boundaries[1].jet_cache
owners = backend.boundaries[1].curved_host[3][jet_id == 0]
self.assertEqual(np.unique(owners).tolist(), [0])
self.assertGreater(len(np.unique(backend._device_boundaries[1].tables["block"][jet_id == 0])), 1)
self.assertGreater(len(np.unique(backend._boundaries[1].tables["block"][jet_id == 0])), 1)
self.assertEqual(source.obj_id, 99)
backend.close()
cfg = config("periodic")
@@ -97,8 +97,8 @@ class JetStepperCPU(unittest.TestCase):
def test_blocked_jet_velocity_matches_legacy_rows(self):
legacy, _ = mocked([circle(actuation(SPEED))], (), FINE_BLOCKS[1])
driven, _ = mocked([circle()], [jet()], FINE_BLOCKS[1])
a = legacy._device_boundaries[1].tables["surface_velocity"]
b = driven._device_boundaries[1].tables["surface_velocity"]
a = legacy._boundaries[1].tables["surface_velocity"]
b = driven._boundaries[1].tables["surface_velocity"]
diff = np.abs(a - b).max()
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)
@@ -106,8 +106,8 @@ class JetStepperCPU(unittest.TestCase):
def test_blocked_amplitude_update_rewrites_only_velocity_column(self):
backend, htod = mocked([circle()], [jet()], FINE_BLOCKS[1])
before = backend._device_boundaries[1].tables["surface_velocity"].copy()
device = backend._device_boundaries[1].device["surface_velocity"]
before = backend._boundaries[1].tables["surface_velocity"].copy()
device = backend._boundaries[1].device["surface_velocity"]
with patch.object(cuda, "memcpy_htod", htod):
htod.reset_mock()
backend.set_jet_amplitude(0, -SPEED, level=1)
@@ -115,7 +115,7 @@ class JetStepperCPU(unittest.TestCase):
self.assertEqual(htod.call_count, 1)
self.assertIs(ptr, device)
np.testing.assert_allclose(sent, -before, rtol=1e-6, atol=0)
np.testing.assert_array_equal(backend._device_boundaries[1].tables["surface_velocity"], sent)
np.testing.assert_array_equal(backend._boundaries[1].tables["surface_velocity"], sent)
with self.assertRaises(ValueError):
backend.set_jet_amplitude(0, 1.5, level=1)
self.assertEqual(htod.call_count, 1)
@@ -133,7 +133,7 @@ class JetStepperCPU(unittest.TestCase):
with self.assertRaises(IndexError):
backend.set_jet_amplitude(1, -SPEED, level=1)
self.assertEqual(backend._status, status)
before = backend._device_boundaries[1].tables["surface_velocity"].copy()
before = backend._boundaries[1].tables["surface_velocity"].copy()
boom = Mock(side_effect=RuntimeError("copy failed"))
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)):
@@ -141,7 +141,7 @@ class JetStepperCPU(unittest.TestCase):
backend.set_jet_amplitude(0, -SPEED, level=1)
self.assertEqual(backend._status, "failed", fine_block)
self.assertEqual(backend.jets[1][0].amplitude, SPEED)
np.testing.assert_array_equal(backend._device_boundaries[1].tables["surface_velocity"], before)
np.testing.assert_array_equal(backend._boundaries[1].tables["surface_velocity"], before)
with self.assertRaisesRegex(RuntimeError, "failed"):
backend.step(1)
with self.assertRaisesRegex(RuntimeError, "failed"):
@@ -102,7 +102,7 @@ class Ledger:
def step(self):
run=self.run;self.samples.clear();run.step()
b=run._device_boundaries[self.index];h=run._exchange_spacing[self.index]
b=run._boundaries[self.index];h=b.spacing
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._device_boundaries[1].tables['aperture'].all())
self.assertTrue(run._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.jets[2][0].amplitude,-.001)
np.testing.assert_array_equal(new._device_boundaries[2].tables['surface_velocity'],old._device_boundaries[2].tables['surface_velocity'])
np.testing.assert_array_equal(new._boundaries[2].tables['surface_velocity'],old._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')