Compare commits

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

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

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-28 02:34:07 +00:00
18 changed files with 109 additions and 74 deletions
+19 -4
View File
@@ -8,8 +8,16 @@ pending a GPU run).
- Level `l = 0 .. depth-1` (0 = root). Pair `hierarchy.pairs[l-1]` builds level `l`. - 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.
+9 -2
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).
+27 -22
View File
@@ -22,7 +22,7 @@ from .channel import channel_masks, validate_channel_config, validate_channel_tr
from ..body.jets import jet_link_velocity from ..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):
@@ -119,7 +123,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 +136,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
@@ -258,7 +263,7 @@ class FixedPatchStepper:
raise ValueError(f"level {level} blocks exceed launch limits") raise ValueError(f"level {level} blocks exceed launch limits")
self._layouts.append(layout) 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). # Host boundary geometry per level in the root edge frame (Hierarchy.edge_frame).
@@ -294,7 +299,7 @@ class FixedPatchStepper:
if depth == 1 and level_objects[0]: 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] + [ self._device_boundaries = [root] + [
_prepare_blocked_boundary(self._boundaries[level], self._layouts[level], _prepare_blocked_boundary(self._boundaries[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)
@@ -310,7 +315,7 @@ class FixedPatchStepper:
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 +371,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 +394,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 +420,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"),
@@ -552,12 +557,12 @@ class FixedPatchStepper:
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 self._boundaries[1] is not None
or any(b is not None for b in self._level_boundaries[:2])): or any(b is not None for b in self._device_boundaries[:2])):
raise ValueError("remap_finest requires a periodic three-level hierarchy with a " 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,
@@ -628,7 +633,7 @@ 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]) jets = list(self._level_jets[level])
jets[index] = jets[index].with_amplitude(amplitude) # validates first jets[index] = jets[index].with_amplitude(amplitude) # validates first
velocity = None velocity = None
@@ -735,7 +740,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 +763,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 +776,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 +873,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
+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._level_jets[2] return boundary.tables["surface_velocity"].copy(), boundary.jets, run._level_jets[2]
@@ -116,8 +116,8 @@ def test_stepper_invalid_input_writes_nothing_and_stays_ready():
with pytest.raises(error): 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._level_jets[2]) == (jets, level_jets) assert (run._device_boundaries[2].jets, run._level_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._level_jets[2][0].amplitude == .017 assert run._level_jets[2][0].amplitude == .017
@@ -136,8 +136,8 @@ def test_stepper_device_failure_keeps_host_and_blocks_reuse(where):
run._stream.synchronize.side_effect = RuntimeError("upload failed") 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._level_jets[2]) == (jets, level_jets) assert (run._device_boundaries[2].jets, run._level_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
+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):
+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_child(B,child_objects=objects) clone=run.remap_child(B,child_objects=objects)
self.assertEqual(clone._level_jets[2][0].amplitude,-jet().amplitude) self.assertEqual(clone._level_jets[2][0].amplitude,-jet().amplitude)
np.testing.assert_array_equal(clone._level_boundaries[2].tables['surface_velocity'], np.testing.assert_array_equal(clone._device_boundaries[2].tables['surface_velocity'],
boundary.tables['surface_velocity']) 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)
@@ -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) out = windowed.read_exchange_window(0)
@@ -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()
@@ -213,7 +214,7 @@ class BlockedFineCPU(unittest.TestCase):
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._fine_tables["surface_velocity"]) self.assertIsNotNone(backend._fine_tables["surface_velocity"])
del backend._level_boundaries[1].device["surface_velocity"] del backend._device_boundaries[1].device["surface_velocity"]
with self.assertRaisesRegex(RuntimeError, "Surface velocity device buffer is missing"): 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]
@@ -107,7 +107,7 @@ 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._fine_tables["surface_velocity"].copy() before = backend._fine_tables["surface_velocity"].copy()
device = backend._level_boundaries[1].device["surface_velocity"] device = backend._device_boundaries[1].device["surface_velocity"]
with patch.object(cuda, "memcpy_htod", htod): with patch.object(cuda, "memcpy_htod", htod):
htod.reset_mock() htod.reset_mock()
backend.set_jet_amplitude(0, -SPEED) backend.set_jet_amplitude(0, -SPEED)
@@ -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._level_jets[2][0].amplitude,-.001) self.assertEqual(new._level_jets[2][0].amplitude,-.001)
np.testing.assert_array_equal(new._level_boundaries[2].tables['surface_velocity'],old._level_boundaries[2].tables['surface_velocity']) np.testing.assert_array_equal(new._device_boundaries[2].tables['surface_velocity'],old._device_boundaries[2].tables['surface_velocity'])
newledger=Ledger(self,new,2);ledger.step();newledger.step() 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')
+12 -6
View File
@@ -81,7 +81,7 @@ def public_jets(center, radius, scenario):
) )
def jet_interval(handle, scenario, begin, n, diameter): def jet_interval(handle, scenario, begin, n, diameter, level=None):
"""Apply a reversal only at a completed interval boundary, before advancing.""" """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,8 +92,10 @@ def jet_interval(handle, scenario, begin, n, diameter):
if begin < switch: if begin < switch:
return min(n, switch - begin), None return min(n, switch - begin), None
if begin == switch: if begin == switch:
handle.set_jet_amplitude(0, -JET_PEAK) # FixedPatchStepper needs an explicit Jet level; Simulation has one level.
handle.set_jet_amplitude(1, JET_PEAK) kw = {} if level is None else {"level": level}
handle.set_jet_amplitude(0, -JET_PEAK, **kw)
handle.set_jet_amplitude(1, JET_PEAK, **kw)
return n, dict(step=begin, tstar=300.0, amplitudes=[-JET_PEAK, JET_PEAK]) return n, dict(step=begin, tstar=300.0, amplitudes=[-JET_PEAK, JET_PEAK])
return n, None return n, None
@@ -256,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
@@ -268,7 +270,11 @@ 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] # Historical pinned source (jet_scenario none) names it blocked_base_tables;
# the current checkout (public Jet scenarios) names it blocked_coarse_tables.
coarse_tables = (refinement_tables.blocked_base_tables if jet_scenario == "none"
else refinement_tables.blocked_coarse_tables)
targets = coarse_tables(pair, layout, np.argsort(pair.coalesce_coarse))[1]
validate_channel_transfers(pair, layout, targets) 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(
@@ -448,7 +454,7 @@ def run_case(a, h, periodic, case):
) )
n = min(max(1, round(width * p["diameter"] / SPEED)), max_steps - begin) n = min(max(1, round(width * p["diameter"] / SPEED)), max_steps - begin)
n, change = jet_interval(backend if refined else sim, a.jet_scenario, n, change = jet_interval(backend if refined else sim, a.jet_scenario,
begin, n, p["diameter"]) begin, n, p["diameter"], level=1 if refined else None)
if change is not None: if change is not None:
amplitude_changes.append(change) amplitude_changes.append(change)
impulse = np.zeros(2) impulse = np.zeros(2)