Compare commits

..
Author SHA1 Message Date
ClaudeandClaude Opus 5.5 1cbdc479cf test_block_step_sync: locate block_stepper without importing it
The module-level import of CelerisLab.lbm.block_stepper loaded
pycuda.driver into the shared test process, unlike the other CPU host
tests (which load mocked copies via AST). find_spec gives the source path
without the import. Test content and assertions unchanged.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-21 01:39:52 +00:00
36 changed files with 360 additions and 322 deletions

No files matched your search

+12 -13
View File
@@ -112,9 +112,8 @@ from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
cfg = LBMConfig(nx=32, ny=24, viscosity=.1) cfg = LBMConfig(nx=32, ny=24, viscosity=.1)
body = SimObject(0, geometry=CircleGeometry(16.0, 12.0, 2.15), center=(16.0, 12.0), radius=2.15, body = SimObject(0, geometry=CircleGeometry(16.0, 12.0, 2.15), center=(16.0, 12.0), radius=2.15,
surface_actuation=SurfaceActuation(0.7, 0.85, 0.002)) # angle, half-width, speed surface_actuation=SurfaceActuation(0.7, 0.85, 0.002)) # angle, half-width, speed
with CudaContext(0) as ctx, FixedPatchStepper(cfg, context=ctx, refinements=[(8, 6, 24, 18)], with CudaContext(0) as ctx, FixedPatchStepper(cfg, (8, 6, 24, 18), context=ctx,
level_blocks=[None, (16, 20)], fine_objects=[body], fine_block=(16, 20)) as run:
level_objects=[(), (body,)]) as run:
rest = np.asarray(compute_feq_d2q9(1.0, 0.0, 0.0), np.float32)[:, None, None] rest = np.asarray(compute_feq_d2q9(1.0, 0.0, 0.0), np.float32)[:, None, None]
# (9, ny, nx) postcollision populations per level; non-leaf cells are discarded. # (9, ny, nx) postcollision populations per level; non-leaf cells are discarded.
run.initialize(np.broadcast_to(rest, (9, *run.patch.coarse_shape)), run.initialize(np.broadcast_to(rest, (9, *run.patch.coarse_shape)),
@@ -134,19 +133,19 @@ Scope and limits:
- Bodies must be stationary circles whose whole boundary stencil lies in the fine - Bodies must be stationary circles whose whole boundary stencil lies in the fine
active region, clear of the coarse/fine transfer cells; violations are rejected active region, clear of the coarse/fine transfer cells; violations are rejected
before any GPU allocation. before any GPU allocation.
- `level_blocks=[None, (by, bx)]` stores the fine level in fixed-size blocks; `by`/`bx` must - `fine_block=(by, bx)` stores the fine level in fixed-size blocks; `by`/`bx` must
divide the WHOLE fine allocation, i.e. `2*(y1-y0)+8` by `2*(x1-x0)+8` including the divide the WHOLE fine allocation, i.e. `2*(y1-y0)+8` by `2*(x1-x0)+8` including the
four-cell ghost band. It gives results identical to dense storage and is four-cell ghost band. It gives results identical to dense storage and is
infrastructure for future block AMR, not a speedup. infrastructure for future block AMR, not a speedup.
- Several regions: pass a sequence of `(x0, y0, x1, y1)` rectangles as - Several regions: pass a sequence of `(x0, y0, x1, y1)` rectangles as
`refinements[0]`, with both levels blocked (`level_blocks`). Transfers are built once `bounds`, with both `fine_block` and `base_block`. Transfers are built once
from their union; adjacent, overlapping and concave unions are accepted when from their union; adjacent, overlapping and concave unions are accepted when
their actual dependencies close. There is no fixed rectangle-gap rule. their actual dependencies close. There is no fixed rectangle-gap rule.
The fine level is then one level-wide grid `(2*ny, 2*nx)`: The fine level is then one level-wide grid `(2*ny, 2*nx)`:
`run.patch.fine_shape` gives it, `initialize`/`download` use it, and only the `run.patch.fine_shape` gives it, `initialize`/`download` use it, and only the
grid-aligned tiles intersecting the union plus its four-cell ghost band are stored. grid-aligned tiles intersecting the union plus its four-cell ghost band are stored.
Block sizes must divide the level-wide grid; rectangle edges need not align Block sizes must divide the level-wide grid; rectangle edges need not align
with tiles. Order of rectangles does not matter. L1 bodies/Jets (`level_objects[1]`/`level_jets[1]`) work with tiles. Order of rectangles does not matter. `fine_objects`/`jets` work
as for one patch: each circle (same coarse cell-edge frame) must lie wholly as for one patch: each circle (same coarse cell-edge frame) must lie wholly
inside the union's active region with inside the union's active region with
its boundary stencil clear of transfer cells; a jet may not select links of two its boundary stencil clear of transfer cells; a jet may not select links of two
@@ -157,8 +156,8 @@ Scope and limits:
cfg = LBMConfig(nx=32, ny=16, viscosity=.1) cfg = LBMConfig(nx=32, ny=16, viscosity=.1)
rest = np.asarray(compute_feq_d2q9(1.0, 0.0, 0.0), np.float32)[:, None, None] rest = np.asarray(compute_feq_d2q9(1.0, 0.0, 0.0), np.float32)[:, None, None]
with CudaContext(0) as ctx, FixedPatchStepper( with CudaContext(0) as ctx, FixedPatchStepper(
cfg, context=ctx, refinements=[((4, 4, 12, 12), (16, 4, 24, 12))], cfg, ((4, 4, 12, 12), (16, 4, 24, 12)), context=ctx,
level_blocks=[(4, 4), (4, 4)]) as run: fine_block=(4, 4), base_block=(4, 4)) as run:
run.initialize(np.broadcast_to(rest, (9, *run.patch.coarse_shape)), run.initialize(np.broadcast_to(rest, (9, *run.patch.coarse_shape)),
np.broadcast_to(rest, (9, *run.patch.fine_shape))) # (9, 32, 64) np.broadcast_to(rest, (9, *run.patch.fine_shape))) # (9, 32, 64)
run.step(32) run.step(32)
@@ -167,12 +166,12 @@ Scope and limits:
``` ```
- Legacy three-level scalar chain: add `child=(x0, y0, x1, y1)` in - Legacy three-level scalar chain: add `child=(x0, y0, x1, y1)` in
**L1 window cells** (the patch's fine grid including its 4-cell ghost band) and **L1 window cells** (the patch's fine grid including its 4-cell ghost band) and
all three `level_blocks` entries. L2 must sit inside L1's `child_block`, with `fine_block` and `base_block`. L2 must sit inside L1's
ACTIVE cells, away from L1's ghost band and transfer footprint (checked before ACTIVE cells, away from L1's ghost band and transfer footprint (checked before
any GPU allocation). `initialize` and `download` take/return three arrays; each any GPU allocation). `initialize` and `download` take/return three arrays; each
level is NaN on its non-ACTIVE (and child-covered) cells. One coarse step = 2 L1 level is NaN on its non-ACTIVE (and child-covered) cells. One coarse step = 2 L1
steps = 4 L2 steps. One L1 patch only. Optional stationary circles and Jets steps = 4 L2 steps. One L1 patch only. Optional stationary circles and Jets
live on L2 via `level_objects[2]` and `level_jets[2]`; L1 bodies/Jets retain live on L2 via `child_objects` and `child_jets`; `fine_objects`/`jets` retain
their L1 meaning and must be empty in this mode. All object and Jet geometry their L1 meaning and must be empty in this mode. All object and Jet geometry
uses the **L0 edge frame**, not the child bounds' L1 window frame. L2 has uses the **L0 edge frame**, not the child bounds' L1 window frame. L2 has
dx=dt=0.25, so prescribed lattice velocity is unchanged. Each whole body and dx=dt=0.25, so prescribed lattice velocity is unchanged. Each whole body and
@@ -189,9 +188,9 @@ Scope and limits:
jet = Jet(Tube((6.5, 8.), 0., 3., .3), Parallel(0.), amplitude=.002) jet = Jet(Tube((6.5, 8.), 0., 3., .3), Parallel(0.), amplitude=.002)
rest = np.asarray(compute_feq_d2q9(1.0, 0.0, 0.0), np.float32)[:, None, None] rest = np.asarray(compute_feq_d2q9(1.0, 0.0, 0.0), np.float32)[:, None, None]
with CudaContext(0) as ctx, FixedPatchStepper( with CudaContext(0) as ctx, FixedPatchStepper(
cfg, context=ctx, refinements=[(4, 4, 12, 12), (8, 8, 16, 16)], cfg, (4, 4, 12, 12), context=ctx, fine_block=(4, 4), base_block=(4, 4),
level_blocks=[(4, 4)] * 3, level_objects=[(), (), (body,)], child=(8, 8, 16, 16), child_block=(4, 4),
level_jets=[(), (), (jet,)]) as run: child_objects=[body], child_jets=[jet]) as run:
shapes = (run.patch.coarse_shape, run.patch.fine_shape, run.child_patch.fine_shape) shapes = (run.patch.coarse_shape, run.patch.fine_shape, run.child_patch.fine_shape)
run.initialize(*(np.broadcast_to(rest, (9, *s)) for s in shapes)) # 16², 24², 24² run.initialize(*(np.broadcast_to(rest, (9, *s)) for s in shapes)) # 16², 24², 24²
run.step(4) run.step(4)
+1 -1
View File
@@ -1,6 +1,6 @@
# Blocked AMR root channel increment # Blocked AMR root channel increment
Opt in with `FixedPatchStepper(..., coarse_boundary="channel", refinements=[bounds], level_blocks=[(by,bx), (fy,fx)])`. Every level must be blocked. The supported boundary is the existing uniform D2Q9/SRT/FP32 double-buffer channel: local Zou-He west inlet, pressure-NEQ east outlet, stationary halfway bounce-back y walls. Uniform and parabolic inlet profiles and shifted/unshifted storage are supported; the existing configuration guards still reject MRT/LES and incompatible storage. This does not change the external DFG MRT reference. Opt in with `FixedPatchStepper(..., coarse_boundary="channel", base_block=(by,bx), fine_block=(fy,fx))`, or the existing level-list form. Every level must be blocked. The supported boundary is the existing uniform D2Q9/SRT/FP32 double-buffer channel: local Zou-He west inlet, pressure-NEQ east outlet, stationary halfway bounce-back y walls. Uniform and parabolic inlet profiles and shifted/unshifted storage are supported; the existing configuration guards still reject MRT/LES and incompatible storage. This does not change the external DFG MRT reference.
`step/block_channel.cu` is a small root-only kernel. Fine kernels, transfer kernels and recursive hop scheduling are unchanged. Global coordinates and the global profile height/velocity are separate from block-local storage strides. Corner ownership and collision phases match `OneStep`: `step/block_channel.cu` is a small root-only kernel. Fine kernels, transfer kernels and recursive hop scheduling are unchanged. Global coordinates and the global profile height/velocity are separate from block-local storage strides. Corner ownership and collision phases match `OneStep`:
+6 -6
View File
@@ -16,8 +16,8 @@ Keep the existing scalar `child=(x0,y0,x1,y1)` for a single scalar parent patch.
Its coordinates retain their original meaning: half-open integer indices in the Its coordinates retain their original meaning: half-open integer indices in the
allocated L1 window, including its ghost band. allocated L1 window, including its ghost band.
A rectangle list as the second `refinements` entry (formerly the keyword The new keyword `child_rectangles_in_parent=[...]` is mutually exclusive with
`child_rectangles_in_parent=[...]`) accepts a scalar or a set parent and interprets every rectangle in `child`. It accepts a scalar or a set parent and interprets every rectangle in
the actual allocated L1 array, including ghost bands and holes. It never infers the actual allocated L1 array, including ghost bands and holes. It never infers
region-local coordinates. The physical root cell edge corresponding to an L1 region-local coordinates. The physical root cell edge corresponding to an L1
index `(x,y)` is `parent.fine_origin + (x,y)/2`: index `(x,y)` is `parent.fine_origin + (x,y)/2`:
@@ -47,9 +47,9 @@ tile edges.
```python ```python
run = FixedPatchStepper( run = FixedPatchStepper(
cfg, context=ctx, cfg, [(4, 4, 20, 28), (24, 4, 44, 28)], context=ctx,
refinements=[[(4, 4, 20, 28), (24, 4, 44, 28)], [(16, 16, 32, 40), (56, 16, 80, 40)]], base_block=(4, 4), fine_block=(4, 4), child_block=(4, 4),
level_blocks=[(4, 4)] * 3, child_rectangles_in_parent=[(16, 16, 32, 40), (56, 16, 80, 40)],
) )
``` ```
@@ -57,7 +57,7 @@ This example uses a 48x32 root. `initialize` and `download` use three Q-major
dense arrays. Public arrays and masks still scale with the full allocated dense arrays. Public arrays and masks still scale with the full allocated
extent; sparse device buffers do not remove that host-memory limit. extent; sparse device buffers do not remove that host-memory limit.
Stationary whole-circle L2 bodies and Jets (`level_objects[2]`, `level_jets[2]`) use root edge Stationary whole-circle `child_objects` and `child_jets` use root edge
coordinates and a single L2 boundary object. Body IDs are physical object-list coordinates and a single L2 boundary object. Body IDs are physical object-list
positions, not region IDs or caller `obj_id`. Complete boundary stencils must positions, not region IDs or caller `obj_id`. Complete boundary stencils must
avoid the transfer footprint. Change a child Jet with avoid the transfer footprint. Change a child Jet with
+2 -2
View File
@@ -90,8 +90,8 @@ def main():
# The caller owns the CUDA context; the stepper must be closed first (the # The caller owns the CUDA context; the stepper must be closed first (the
# nested with-statement does that). # nested with-statement does that).
with CudaContext(0) as context, FixedPatchStepper( with CudaContext(0) as context, FixedPatchStepper(
cfg, context=context, coarse_boundary=args.mode, refinements=[BOUNDS], cfg, BOUNDS, context=context, coarse_boundary=args.mode,
level_blocks=[None, fine_block], level_objects=[(), (body,)]) as run: fine_objects=[body], fine_block=fine_block) as run:
p = run.patch p = run.patch
# Rest equilibrium is already postcollision G, so no initial collision is # Rest equilibrium is already postcollision G, so no initial collision is
# needed. Non-leaf entries (covered, ghost, solid) are discarded by initialize. # needed. Non-leaf entries (covered, ghost, solid) are discarded by initialize.
+2 -2
View File
@@ -148,8 +148,8 @@ def main():
samples=[], status='running', last_synchronized_coarse_step=0) samples=[], status='running', last_synchronized_coarse_step=0)
start, run = time.monotonic(), None start, run = time.monotonic(), None
try: try:
with CudaContext(0) as context, FixedPatchStepper(cfg, context=context, with CudaContext(0) as context, FixedPatchStepper(cfg, BOUNDS, context=context,
coarse_boundary='slot_reservoir', refinements=[BOUNDS]) as run: coarse_boundary='slot_reservoir') as run:
p = run.patch p = run.patch
masks = [np.ones(p.coarse_shape, bool), p.fine_kind == ACTIVE] masks = [np.ones(p.coarse_shape, bool), p.fine_kind == ACTIVE]
masks[0][[0, -1]] = False masks[0][[0, -1]] = False
+1 -1
View File
@@ -48,7 +48,7 @@ def main():
nx=16, ny=16, viscosity=.1, velocity=1e-3, rho=1.0, nx=16, ny=16, viscosity=.1, velocity=1e-3, rho=1.0,
lattice_model="D2Q9", collision="SRT", streaming="double_buffer", lattice_model="D2Q9", collision="SRT", streaming="double_buffer",
store_precision="FP32", ddf_shifting=False, les_enabled=False) store_precision="FP32", ddf_shifting=False, les_enabled=False)
with CudaContext(0) as context, PeriodicPatchStepper(cfg, context=context, refinements=[(4, 4, 12, 12)]) as run: with CudaContext(0) as context, PeriodicPatchStepper(cfg, (4, 4, 12, 12), context=context) as run:
p = run.patch p = run.patch
coarse = shear_equilibrium(p.coarse_shape, (0, 0), 1.0, cfg.ny) coarse = shear_equilibrium(p.coarse_shape, (0, 0), 1.0, cfg.ny)
fine = shear_equilibrium(p.fine_shape, p.fine_origin, .5, cfg.ny) fine = shear_equilibrium(p.fine_shape, p.fine_origin, .5, cfg.ny)
+135 -60
View File
@@ -65,8 +65,52 @@ def _level_name(level):
return ("base", "fine", "child")[level] if level < 3 else f"level {level}" return ("base", "fine", "child")[level] if level < 3 else f"level {level}"
def _legacy_levels(bounds, *, fine_block, base_block, child, child_block, child_rectangles_in_parent,
fine_objects, jets, child_objects, child_jets):
"""Normalise the original two/three-level keywords to per-level lists (original errors)."""
if child is not None and child_rectangles_in_parent is not None:
raise ValueError("child and child_rectangles_in_parent are mutually exclusive")
nested = child is not None or child_rectangles_in_parent is not None
bounds = tuple(bounds)
if not is_rectangle(bounds) and (fine_block is None or base_block is None):
raise ValueError("several patches require fine_block and base_block")
fine_objects, jets = tuple(fine_objects), tuple(jets)
child_objects, child_jets = tuple(child_objects), tuple(child_jets)
if not nested and (child_objects or child_jets):
raise ValueError("child_objects/child_jets require a nested child")
if base_block is not None and fine_block is None:
raise ValueError("base_block requires fine_block")
refinements, blocks = [bounds], [base_block, fine_block]
objects, level_jets = [(), fine_objects], [(), jets]
if nested:
if child is not None and not is_rectangle(bounds):
raise ValueError("scalar child requires one parent patch; use child_rectangles_in_parent")
if fine_block is None or base_block is None or child_block is None:
raise ValueError("a nested child needs fine_block, base_block and child_block")
if fine_objects or jets:
raise ValueError("a nested child supports bodies/Jets only on L2; use child_objects/child_jets")
refinements.append(tuple(child) if child is not None else tuple(child_rectangles_in_parent))
blocks.append(child_block)
objects.append(child_objects)
level_jets.append(child_jets)
return refinements, blocks, objects, level_jets
def _initial_levels(level_g, coarse_g, fine_g, child_g):
"""One tuple of per-level arrays from initialize()'s positional or legacy named form."""
if coarse_g is None and fine_g is None:
if child_g is not None and len(level_g) != 2:
raise TypeError("child_g needs exactly the coarse and fine arrays before it")
return tuple(level_g) + (() if child_g is None else (child_g,))
if level_g:
raise TypeError("pass level arrays positionally or as coarse_g/fine_g/child_g, not both")
if coarse_g is None or fine_g is None:
raise TypeError("coarse_g and fine_g are both required in the named form")
return (coarse_g, fine_g) + (() if child_g is None else (child_g,))
def _check_levels(refinements, blocks, objects, jets): def _check_levels(refinements, blocks, objects, jets):
"""Candidate limits of the level lists, checked before any device work.""" """Candidate limits shared by the list and the legacy entry points."""
depth = len(refinements) + 1 depth = len(refinements) + 1
if depth < 2: if depth < 2:
raise ValueError("at least one refinement level is required") raise ValueError("at least one refinement level is required")
@@ -90,14 +134,7 @@ def _check_levels(refinements, blocks, objects, jets):
class FixedPatchStepper: class FixedPatchStepper:
"""D2Q9 SRT/MRT FP32 fixed 2:1 refinement hierarchy, built from level lists. """D2Q9/SRT/FP32, one unforced internal 2:1 patch in an x-periodic box.
``refinements[l]`` gives level l+1 in level-l integer cells: one rectangle
(x0, y0, x1, y1) (FixedPatch2D, a fine window with its own origin) or a
sequence of rectangles (FixedPatchSet2D, a dependency-validated union on one
level-wide grid, public fine arrays (2*ny, 2*nx), NaN outside ACTIVE cells);
see lbm.hierarchy. ``level_blocks``, ``level_objects`` and ``level_jets``
hold one entry per level, L0 included (L0 bodies/Jets are not supported).
``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``.
@@ -108,53 +145,68 @@ class FixedPatchStepper:
actual channel closure footprint. Channel remapping is not supported. actual channel closure footprint. Channel remapping is not supported.
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.
``fine_objects`` accepts fixed circles in coarse physical coordinates;
velocities are in coarse lattice units (dx=dt=1). Independent local copies
are projected at fine cell centers. Their complete stencil must be active
and disjoint from the actual two-hop transfer footprint.
These use the edge/cell-center frame (coarse centers I+.5), whereas ordinary Simulation objects use the integer-node frame.
Bodies are fixed circles in root (L0) physical coordinates; velocities are ``initialize(coarse_g, fine_g)`` accepts physical Q-major (9,ny,nx) arrays
in root lattice units (dx=dt=1). Independent local copies are projected at already in postcollision phase. It does not collide or build channel flags.
fine cell centers; their complete stencil must be active and disjoint from Every ``step()`` publishes only after both levels complete on one stream.
the actual two-hop transfer footprint. These use the edge/cell-center frame On a device/launch failure the instance becomes unusable: close and recreate
(coarse centers I+.5), whereas ordinary Simulation objects use the it, rather than pretending the partially advanced fine buffer rolled back.
integer-node frame. Jets (body.jets.Jet) use the same frame and drive cut
links selected at their true hits; ``set_jet_amplitude(index, amplitude,
level=l)`` re-uploads the per-link velocity column only.
``initialize(*level_g)`` accepts one physical Q-major (9, ny, nx) array per
level, already in postcollision phase. It does not collide or build channel
flags. Every ``step()`` publishes only after all levels complete on one
stream. On a device/launch failure the instance becomes unusable: close and
recreate it, rather than pretending a partially advanced buffer rolled back.
A ``level_blocks`` entry (by, bx) stores that level in blocks.cu storage:
the ENTIRE level shape, including the four-cell patch ghost band, is
partitioned by a non-periodic BlockLayout2D with one-cell same-level halos.
Patch kinds stay per cell (OUTER invalid, INNER stream only, ACTIVE
collide); halos are gathered before each hop. Public arrays keep the dense
shapes; only the allocation differs. A blocked level requires blocked
refined levels below it; several L1 rectangles and depth >= 3 require all
levels blocked. Blocked bodies are stationary circles, optionally
surface-actuated: prepare_fine_boundary decides legality on the dense grid,
then cut links are staged receiver-local per hop: halo -> BFL -> step, same
buffer. A blocked parent gathers halos before receiver-local coalescence
staging and pull/collision; covered cells are poisoned on every step.
``fine_block=(by, bx)`` (opt-in) stores the fine level in
blocks.cu storage: the ENTIRE fine_shape, including the four-cell patch
ghost band, is partitioned by a non-periodic BlockLayout2D with one-cell
same-level halos. Patch kinds stay per cell (OUTER invalid, INNER stream
only, ACTIVE collide); halos are gathered before each fine hop. Public
arrays keep the dense (9, ny, nx) fine shape; only the allocation differs.
Blocked fine_objects are stationary circles, optionally surface-actuated:
prepare_fine_boundary decides legality on the dense fine grid, then cut
links are staged receiver-local per hop: halo -> BFL -> step, same buffer.
``base_block=(by, bx)`` also blocks the parent; it requires ``fine_block``.
Parent halos precede receiver-local coalescence staging and pull/collision;
covered cells are explicitly poisoned on every step. Public shapes stay dense.
``observe_raw_exchange=True`` records blocked physical cut-link exchange ``observe_raw_exchange=True`` records blocked physical cut-link exchange
without changing populations. read_raw_exchange reports the last complete without changing populations. read_raw_exchange reports the last complete
base interval or last substep in base-lattice units per unit depth. This base interval or last substep in base-lattice units per unit depth. This
is boundary momentum reaction, not a validated nozzle/stress load. Dense is boundary momentum reaction, not a validated nozzle/stress load. Dense
body levels are not supported by this diagnostic. body layouts are not supported by this new diagnostic in this candidate.
``observe_window_exchange=True`` (requires observe_raw_exchange) adds a second ``observe_window_exchange=True`` (requires observe_raw_exchange) adds a second
per-link FP64 ledger accumulated by the same kernel between per-link FP64 ledger accumulated by the same kernel between
start_exchange_window() and read_exchange_window(); disabled by default. start_exchange_window() and read_exchange_window(); disabled by default.
``jets`` (body.jets.Jet) use the same edge frame as ``fine_objects`` and
drive fine cut links selected at their true hits; ``set_jet_amplitude``
re-uploads the per-link velocity column only.
Nested levels: a child set's coordinates index the allocated parent array, Several patches: pass ``bounds`` as a sequence of (x0, y0, x1, y1)
including ghost bands and holes: root edge = parent fine_origin + index / 2 rectangles (FixedPatchSet2D, dependency-validated union). This requires
(zero origin for a parent set; a scalar parent keeps its shifted window fine_block and base_block; fine_objects/jets are projected once onto the
origin). A scalar child requires a scalar parent. Each pair requires gap >= 4 level-wide fine grid (each circle inside the union, body ids = object
in its own parent-cell units. All regions share one buffer pair per level positions, jet and writer checks level-wide), and
and the fixed schedule of 2**l hops per base step at level l. Bodies/Jets are supported only on the uses one level-wide fine grid (2*ny, 2*nx): public fine arrays have that
finest level of a nested hierarchy (L2 dx=dt=.25 in root units, unchanged shape, NaN outside the patches' ACTIVE cells. The step sequence is the
lattice velocity); every whole circle and its BFL stencil must fit strictly blocked one above, with unique union transfers in one launch each.
away from that level's transfer footprint.
Three-level separated unions: ``child_rectangles_in_parent=[...]`` is
mutually exclusive with scalar ``child`` and accepts either parent form.
Coordinates index the allocated L1 array, including ghost bands and holes:
root edge = patch.fine_origin + L1 index / 2. A parent set has origin zero;
a scalar parent retains its shifted window origin. The child set's public
extent is twice the allocated L1 shape, even though storage is sparse.
Each pair requires gap >= 4 in its own parent-cell units. All regions share
one buffer pair per level and the fixed 1:2:4 schedule; L2-only bodies/Jets
below apply unchanged, with object IDs assigned over the whole level.
Three-level single chain: ``child`` bounds use L1 window cells; all three
levels must be blocked. ``child_objects`` and ``child_jets`` use the L0
edge frame, with L2 dx=dt=.25 and unchanged lattice velocity. L1 bodies
and Jets remain unsupported in a chain. Every whole circle and its BFL
stencil must fit strictly away from the L2 transfer footprint. L2 runs
halo -> BFL -> step four times per base step. Runtime child Jet changes
require ``set_jet_amplitude(index, amplitude, level=2)``.
``explosion="projected_linear"`` selects experimental transverse linear ``explosion="projected_linear"`` selects experimental transverse linear
reconstruction on each adjacent pair. Default "constant" is unchanged. reconstruction on each adjacent pair. Default "constant" is unchanged.
@@ -162,21 +214,43 @@ class FixedPatchStepper:
bodies/Jets or LES; this backend supplies no forcing. Geometry is generic, bodies/Jets or LES; this backend supplies no forcing. Geometry is generic,
but union/multilevel GPU and matched-physics qualification remain pending. but union/multilevel GPU and matched-physics qualification remain pending.
One recursive hop advances every level: a child explodes from its parent's Level lists (any depth): ``refinements[l]`` gives level l+1 in level-l
current hop source, takes two hops, then the parent gathers halos, stages integer cells (scalar window or rectangle list, see lbm.hierarchy);
restriction, applies BFL and steps. ``level_blocks``/``level_objects``/``level_jets`` hold one entry per level,
L0 included. They are mutually exclusive with the legacy keywords above,
which are normalised to the same lists. One recursive hop advances every
level: a child explodes from its parent's current hop source, takes two
hops, then the parent gathers halos, stages restriction, applies BFL and
steps. Depth >= 3 must be fully blocked with bodies/Jets on the finest level.
""" """
def __init__(self, cfg, *, context, refinements, coarse_boundary="periodic", def __init__(self, cfg, bounds=None, *, context, coarse_boundary="periodic", fine_objects=(),
level_blocks=None, level_objects=None, level_jets=None, fine_block=None, base_block=None, jets=(), child=None, child_block=None,
child_objects=(), child_jets=(), child_rectangles_in_parent=None,
refinements=None, level_blocks=None, level_objects=None, level_jets=None,
observe_raw_exchange=False, observe_window_exchange=False, explosion="constant"): observe_raw_exchange=False, observe_window_exchange=False, explosion="constant"):
if HALO != 4: if HALO != 4:
raise RuntimeError("refinement.cu requires a four-cell fine halo") raise RuntimeError("refinement.cu requires a four-cell fine halo")
if coarse_boundary not in ("periodic", "slot_reservoir", "channel"): if coarse_boundary not in ("periodic", "slot_reservoir", "channel"):
raise ValueError("coarse_boundary must be periodic, slot_reservoir or channel") raise ValueError("coarse_boundary must be periodic, slot_reservoir or channel")
depth = len(refinements) + 1 if refinements is None:
level_blocks = [None] * depth if level_blocks is None else list(level_blocks) if bounds is None:
level_objects = [()] * depth if level_objects is None else list(level_objects) raise ValueError("bounds or refinements is required")
level_jets = [()] * depth if level_jets is None else list(level_jets) if level_blocks is not None or level_objects is not None or level_jets is not None:
raise ValueError("level_blocks/level_objects/level_jets require refinements")
refinements, level_blocks, level_objects, level_jets = _legacy_levels(
bounds, fine_block=fine_block, base_block=base_block, child=child,
child_block=child_block, child_rectangles_in_parent=child_rectangles_in_parent,
fine_objects=fine_objects, jets=jets, child_objects=child_objects,
child_jets=child_jets)
else:
if (bounds is not None or child is not None or child_rectangles_in_parent is not None
or fine_block is not None or base_block is not None or child_block is not None
or tuple(fine_objects) or tuple(jets) or tuple(child_objects) or tuple(child_jets)):
raise ValueError("refinements and the legacy patch keywords are mutually exclusive")
depth = len(refinements) + 1
level_blocks = [None] * depth if level_blocks is None else list(level_blocks)
level_objects = [()] * depth if level_objects is None else list(level_objects)
level_jets = [()] * depth if level_jets is None else list(level_jets)
level_objects = [tuple(o) for o in level_objects] level_objects = [tuple(o) for o in level_objects]
level_jets = [tuple(j) for j in level_jets] level_jets = [tuple(j) for j in level_jets]
level_blocks = [None if b is None else tuple(b) for b in level_blocks] level_blocks = [None if b is None else tuple(b) for b in level_blocks]
@@ -466,13 +540,13 @@ class FixedPatchStepper:
if self._ddf_shifting and any(layout is None for layout in self._layouts): if self._ddf_shifting and any(layout is None for layout in self._layouts):
raise ValueError("shifted populations require all-blocked levels") raise ValueError("shifted populations require all-blocked levels")
def initialize(self, *level_g): def initialize(self, *level_g, coarse_g=None, fine_g=None, child_g=None):
"""Upload physical postcollision G, encoding before any FP32 rounding. """Upload physical postcollision G, encoding before any FP32 rounding.
Shifted input should be physical FP64 to retain small perturbations. Shifted input should be physical FP64 to retain small perturbations.
Public initialization resets the step count; non-leaves are discarded. Public initialization resets the step count; non-leaves are discarded.
""" """
given = tuple(level_g) given = _initial_levels(level_g, coarse_g, fine_g, child_g)
with self._lock, self.ctx.activate(): with self._lock, self.ctx.activate():
self._require("new", "ready") self._require("new", "ready")
self._check_storage() self._check_storage()
@@ -596,9 +670,10 @@ class FixedPatchStepper:
old, new, source, base_count=count, level_steps=(count, 2*count, 4*count), old, new, source, base_count=count, level_steps=(count, 2*count, 4*count),
phase="postcollision_g", **options) phase="postcollision_g", **options)
candidate = type(self)( candidate = type(self)(
cfg, context=self.ctx, refinements=[pairs[0].bounds, new.child.bounds], cfg, pairs[0].bounds, context=self.ctx, base_block=base_block,
level_blocks=[base_block, fine_block, child_block], fine_block=fine_block, child_block=child_block,
level_objects=[(), (), child_objects], level_jets=[(), (), jets], child_rectangles_in_parent=new.child.bounds,
child_objects=child_objects, child_jets=jets,
observe_raw_exchange=self._observe_raw_exchange, observe_raw_exchange=self._observe_raw_exchange,
observe_window_exchange=self._observe_window, explosion=self._explosion) observe_window_exchange=self._observe_window, explosion=self._explosion)
try: try:
+4 -2
View File
@@ -1,6 +1,7 @@
"""Host contract of PeriodicBlockStepper.step(n): one stream wait per call (mocked driver).""" """Host contract of PeriodicBlockStepper.step(n): one stream wait per call (mocked driver)."""
import ast import ast
from contextlib import nullcontext from contextlib import nullcontext
import importlib.util
from pathlib import Path from pathlib import Path
import types import types
from unittest.mock import Mock from unittest.mock import Mock
@@ -8,14 +9,15 @@ import unittest
import numpy as np import numpy as np
import CelerisLab.lbm.block_stepper as b
from CelerisLab.config import LBMConfig from CelerisLab.config import LBMConfig
from CelerisLab.lbm.blocks import BlockLayout2D from CelerisLab.lbm.blocks import BlockLayout2D
def host_module(): def host_module():
"""block_stepper with pycuda/compiler replaced by mocks; no GPU context.""" """block_stepper with pycuda/compiler replaced by mocks; no GPU context."""
path = Path(b.__file__) # Locate the source without importing it: importing block_stepper loads pycuda.driver,
# which CPU-only import-boundary tests (test_build_contracts) must never see.
path = Path(importlib.util.find_spec("CelerisLab.lbm.block_stepper").origin)
tree = ast.parse(path.read_text()) tree = ast.parse(path.read_text())
tree.body = [node for node in tree.body if not ( tree.body = [node for node in tree.body if not (
isinstance(node, ast.Import) and any(a.name == "pycuda.driver" for a in node.names) isinstance(node, ast.Import) and any(a.name == "pycuda.driver" for a in node.names)
+2 -2
View File
@@ -80,8 +80,8 @@ def test_evolved_ghost_snapshot_and_physical_download():
host.cuda.memcpy_dtoh.side_effect = lambda a, ptr: np.copyto(a, storage[id(ptr)]) host.cuda.memcpy_dtoh.side_effect = lambda a, ptr: np.copyto(a, storage[id(ptr)])
context = Mock(activate=nullcontext, sm_arch="sm_70") context = Mock(activate=nullcontext, sm_arch="sm_70")
cfg = LBMConfig(nx=24, ny=20, viscosity=.1, velocity=.025) cfg = LBMConfig(nx=24, ny=20, viscosity=.1, velocity=.025)
backend = host.FixedPatchStepper(cfg, context=context, coarse_boundary="channel", backend = host.FixedPatchStepper(cfg, (6, 6, 18, 14), context=context,
refinements=[(6, 6, 18, 14)], level_blocks=[(4, 6), (4, 4)]) coarse_boundary="channel", base_block=(4, 6), fine_block=(4, 4))
try: try:
weights = np.array([4/9]+[1/9]*4+[1/36]*4)[:, None, None] weights = np.array([4/9]+[1/9]*4+[1/36]*4)[:, None, None]
arrays = [np.broadcast_to(weights, level.shape).copy() for level in backend._levels] arrays = [np.broadcast_to(weights, level.shape).copy() for level in backend._levels]
+2 -3
View File
@@ -157,9 +157,8 @@ class RawExchange(unittest.TestCase):
host = host_module() host = host_module()
ctx = Mock() ctx = Mock()
with self.assertRaisesRegex(ValueError, 'blocked body'): with self.assertRaisesRegex(ValueError, 'blocked body'):
host.FixedPatchStepper(LBMConfig(nx=32, ny=32, viscosity=.1), context=ctx, host.FixedPatchStepper(LBMConfig(nx=32, ny=32, viscosity=.1), (4, 4, 28, 28),
refinements=[(4, 4, 28, 28)], context=ctx, fine_objects=[body(radius=2.15)], observe_raw_exchange=True)
level_objects=[(), (body(radius=2.15),)], observe_raw_exchange=True)
ctx.activate.assert_not_called() ctx.activate.assert_not_called()
def test_all_depth_weights(self): def test_all_depth_weights(self):
@@ -173,10 +173,8 @@ def bodies_and_jets(speed=.002):
return bodies,jets return bodies,jets
def levels(parents=PARENTS,children=CHILDREN,blocks=((4,4),)*3,objects=((),(),()),jets=((),(),())): KW=dict(base_block=(4,4),fine_block=(4,4),child_block=(4,4),
"""Level-list constructor keywords for the L0 set, L1 set and L2 set used here.""" child_rectangles_in_parent=CHILDREN)
return dict(refinements=[parents,children],level_blocks=list(blocks),
level_objects=list(objects),level_jets=list(jets))
def mock_backend(*, children=CHILDREN, bodies=False, fail_compile=False, fail_upload=False): def mock_backend(*, children=CHILDREN, bodies=False, fail_compile=False, fail_upload=False):
@@ -203,8 +201,9 @@ def mock_backend(*, children=CHILDREN, bodies=False, fail_compile=False, fail_up
module.cuda.memcpy_htod.side_effect=upload module.cuda.memcpy_htod.side_effect=upload
objects,jets=bodies_and_jets() if bodies else ((),()) objects,jets=bodies_and_jets() if bodies else ((),())
try: try:
run=module.FixedPatchStepper(CFG,context=context, run=module.FixedPatchStepper(CFG,PARENTS,context=context,
**levels(children=children,objects=((),(),objects),jets=((),(),jets))) child_objects=objects,child_jets=jets,
**dict(KW,child_rectangles_in_parent=children))
except RuntimeError: except RuntimeError:
for allocation in allocations: for allocation in allocations:
allocation.free.assert_called_once_with() allocation.free.assert_called_once_with()
@@ -215,29 +214,28 @@ def mock_backend(*, children=CHILDREN, bodies=False, fail_compile=False, fail_up
class ConstructorAndSchedule(unittest.TestCase): class ConstructorAndSchedule(unittest.TestCase):
def test_invalid_inputs_fail_before_context_or_allocations(self): def test_invalid_inputs_fail_before_context_or_allocations(self):
context=Mock() context=Mock()
b=(4,4) cases=[dict(child=(8,8,16,16)),dict(child_rectangles_in_parent=[]),
cases=[dict(children=(8,8,16,16)),dict(children=[]), dict(child_rectangles_in_parent=[(16,16,32)]),
dict(children=[(16,16,32)]), dict(child_rectangles_in_parent=[(16.5,16,32,40)]),
dict(children=[(16.5,16,32,40)]), dict(child_rectangles_in_parent=[(16,16,32,24),(16,27,32,55)]),
dict(children=[(16,16,32,24),(16,27,32,55)]), dict(child_rectangles_in_parent=[(40,16,48,40)]),
dict(children=[(40,16,48,40)]), dict(child_rectangles_in_parent=[(8,16,32,40)]),
dict(children=[(8,16,32,40)]), dict(child_block=None),dict(fine_block=None),dict(base_block=None),
dict(blocks=(b,b,None)),dict(blocks=(b,None,b)),dict(blocks=(None,b,b)), dict(child_block=(5,4)),dict(fine_block=(5,4)),
dict(blocks=(b,b,(5,4))),dict(blocks=(b,(5,4),b)), dict(child_block=(4,)),dict(child_rectangles_in_parent=(16,16,32,40))]
dict(blocks=(b,b,(4,))),dict(children=(16,16,32,40))]
for changes in cases: for changes in cases:
with self.subTest(changes=changes),self.assertRaises((ValueError,TypeError)): with self.subTest(changes=changes),self.assertRaises((ValueError,TypeError)):
host.FixedPatchStepper(CFG,context=context,**levels(**changes)) host.FixedPatchStepper(CFG,PARENTS,context=context,**dict(KW,**changes))
with self.assertRaisesRegex(ValueError,'margin'): with self.assertRaisesRegex(ValueError,'margin'):
host.FixedPatchStepper(CFG,context=context,**levels(parents=[(1,4,20,28),(23,4,44,28)])) host.FixedPatchStepper(CFG,[(1,4,20,28),(23,4,44,28)],context=context,**KW)
objects,jets=bodies_and_jets() objects,jets=bodies_and_jets()
with self.assertRaisesRegex(ValueError,'only on the finest level'): with self.assertRaisesRegex(ValueError,'only on L2'):
host.FixedPatchStepper(CFG,context=context,**levels(objects=((),objects,()))) host.FixedPatchStepper(CFG,PARENTS,context=context,fine_objects=objects,**KW)
from CelerisLab.body.jets import Jet,Tube,Parallel from CelerisLab.body.jets import Jet,Tube,Parallel
spanning=Jet(Tube((10.5,14.),0.,25.,.3),Parallel(0.),amplitude=0.) spanning=Jet(Tube((10.5,14.),0.,25.,.3),Parallel(0.),amplitude=0.)
with self.assertRaisesRegex(ValueError,'one jet'): with self.assertRaisesRegex(ValueError,'one jet'):
host.FixedPatchStepper(CFG,context=context, host.FixedPatchStepper(CFG,PARENTS,context=context,child_objects=objects,
**levels(objects=((),(),objects),jets=((),(),[spanning]))) child_jets=[spanning],**KW)
context.activate.assert_not_called() context.activate.assert_not_called()
host.cuda.mem_alloc.assert_not_called() host.cuda.mem_alloc.assert_not_called()
@@ -315,23 +313,21 @@ class ConstructorAndSchedule(unittest.TestCase):
class RemovedLegacyAPI(unittest.TestCase): class InitialLevels(unittest.TestCase):
"""The pre-level-list constructor/initialize forms are gone and fail loudly.""" """initialize() entry adapter: positional per-level arrays or the legacy names, never both."""
def test_removed_constructor_forms_raise_type_error(self): def test_forms(self):
import inspect a, b, c = object(), object(), object()
context = Mock() f = host._initial_levels
with self.assertRaises(TypeError): # positional bounds self.assertEqual(f((a, b), None, None, None), (a, b))
host.FixedPatchStepper(CFG, PARENTS, context=context) self.assertEqual(f((a, b, c), None, None, None), (a, b, c))
removed = dict(bounds=PARENTS, fine_block=(4, 4), base_block=(4, 4), fine_objects=(), self.assertEqual(f((a, b), None, None, c), (a, b, c)) # old initialize(c, f, child_g=g)
jets=(), child=(8, 8, 16, 16), child_block=(4, 4), child_objects=(), self.assertEqual(f((), a, b, None), (a, b))
child_jets=(), child_rectangles_in_parent=CHILDREN) self.assertEqual(f((), a, b, c), (a, b, c))
for name, value in removed.items(): for args in (((a,), None, None, c), ((a, b, c), None, None, c), ((a,), a, b, None),
with self.subTest(keyword=name), self.assertRaises(TypeError): ((), a, None, None), ((), None, b, c)):
host.FixedPatchStepper(CFG, context=context, **levels(), **{name: value}) with self.subTest(args=args), self.assertRaises(TypeError):
context.activate.assert_not_called() f(*args)
self.assertEqual(list(inspect.signature(host.FixedPatchStepper.initialize).parameters),
["self", "level_g"]) # no coarse_g/fine_g/child_g
if __name__=='__main__': if __name__=='__main__':
+6 -5
View File
@@ -54,13 +54,14 @@ class ShiftedJets(unittest.TestCase):
(Jet(Tube((10.,14.),0.,26.,.3),Parallel(0.)),))] (Jet(Tube((10.,14.),0.,26.,.3),Parallel(0.)),))]
for objects,jets in cases: for objects,jets in cases:
with self.subTest(objects=len(objects),jets=len(jets)),self.assertRaises(ValueError): with self.subTest(objects=len(objects),jets=len(jets)),self.assertRaises(ValueError):
module.FixedPatchStepper(CFG,context=None,refinements=[backend.PARENTS,A], module.FixedPatchStepper(CFG,backend.PARENTS,context=None,
level_blocks=[(4,4)]*3,level_objects=[(),(),objects],level_jets=[(),(),jets]) base_block=(4,4),fine_block=(4,4),child_block=(4,4),
child_rectangles_in_parent=A,child_objects=objects,child_jets=jets)
for blocks in ((None,None,None),((4,4),None,(4,4))): for blocks in ((None,None,None),((4,4),None,(4,4))):
with self.assertRaises(ValueError): with self.assertRaises(ValueError):
module.FixedPatchStepper(CFG,context=None,refinements=[backend.PARENTS,A], module.FixedPatchStepper(CFG,backend.PARENTS,context=None,
level_blocks=list(blocks),level_objects=[(),(),(body(),)], base_block=blocks[0],fine_block=blocks[1],child_block=blocks[2],
level_jets=[(),(),(jet(),)]) child_rectangles_in_parent=A,child_objects=(body(),),child_jets=(jet(),))
module.cuda.mem_alloc.assert_not_called() module.cuda.mem_alloc.assert_not_called()
module.compiler.compile_kernel.assert_not_called() module.compiler.compile_kernel.assert_not_called()
+3 -4
View File
@@ -129,10 +129,9 @@ class ShiftedStorage(unittest.TestCase):
def test_dense_or_mixed_shifted_rejected_before_compilation(self): def test_dense_or_mixed_shifted_rejected_before_compilation(self):
module=backend.host_module() module=backend.host_module()
cfg=replace(backend.CFG,ddf_shifting=True) cfg=replace(backend.CFG,ddf_shifting=True)
for blocks in ([None,None],[(4,4),None],[None,(4,4)]): for kwargs in ({},dict(base_block=(4,4)),dict(fine_block=(4,4))):
with self.assertRaisesRegex(ValueError,'all-blocked|blocked refined levels below'): with self.assertRaisesRegex(ValueError,'all-blocked|base_block requires fine_block'):
module.FixedPatchStepper(cfg,context=None,refinements=[(4,4,20,28)], module.FixedPatchStepper(cfg,(4,4,20,28),context=None,**kwargs)
level_blocks=blocks)
module.compiler.compile_kernel.assert_not_called() module.compiler.compile_kernel.assert_not_called()
def test_uniform_block_codec_and_raw_snapshot(self): def test_uniform_block_codec_and_raw_snapshot(self):
@@ -106,10 +106,9 @@ class ChannelRootGPU(unittest.TestCase):
cfg = LBMConfig(nx=48, ny=32, viscosity=.1, velocity=.015, ddf_shifting=shifted) cfg = LBMConfig(nx=48, ny=32, viscosity=.1, velocity=.015, ddf_shifting=shifted)
reference = None reference = None
for block in ((32,48), (8,8), (4,1)): for block in ((32,48), (8,8), (4,1)):
with FixedPatchStepper(cfg, context=context, coarse_boundary="channel", with FixedPatchStepper(cfg, (8,6,40,26), context=context,
refinements=[(8,6,40,26)], level_blocks=[block, (8,8)], coarse_boundary="channel", base_block=block, fine_block=(8,8),
level_objects=[(), (body,) if bodies else ()], fine_objects=(body,) if bodies else (), jets=(jet,) if bodies else ()) as backend:
level_jets=[(), (jet,) if bodies else ()]) as backend:
arrays = [state(*level.shape[1:]) for level in backend._levels] arrays = [state(*level.shape[1:]) for level in backend._levels]
backend.initialize(*arrays) backend.initialize(*arrays)
backend.step(4) backend.step(4)
@@ -25,8 +25,8 @@ EPS = np.finfo(np.float64).eps
def build(m, window=True, raw=True): def build(m, window=True, raw=True):
return FixedPatchStepper(config("periodic"), context=m["context"], refinements=[BOUNDS], return FixedPatchStepper(config("periodic"), BOUNDS, context=m["context"], base_block=(3, 4),
level_blocks=[(3, 4), (16, 20)], level_objects=[(), (circle(),)], fine_block=(16, 20), fine_objects=[circle()],
observe_raw_exchange=raw, observe_window_exchange=window) observe_raw_exchange=raw, observe_window_exchange=window)
@@ -109,11 +109,10 @@ class ExchangeWindowGPU(unittest.TestCase):
from test_raw_exchange_gpu import BOUNDS as RB, CFG, circle as rc, jet, state from test_raw_exchange_gpu import BOUNDS as RB, CFG, circle as rc, jet, state
from CelerisLab.cuda.context import CudaContext from CelerisLab.cuda.context import CudaContext
with CudaContext(0) as ctx: with CudaContext(0) as ctx:
kw = dict(context=ctx, refinements=[RB], level_blocks=[None, (16, 20)], kw = dict(context=ctx, fine_block=(16, 20), fine_objects=[rc()], jets=[jet()],
level_objects=[(), (rc(),)], level_jets=[(), (jet(),)],
observe_raw_exchange=True) observe_raw_exchange=True)
with FixedPatchStepper(CFG, observe_window_exchange=True, **kw) as windowed, \ with FixedPatchStepper(CFG, RB, observe_window_exchange=True, **kw) as windowed, \
FixedPatchStepper(CFG, **kw) as plain: FixedPatchStepper(CFG, RB, **kw) as plain:
arrays = [state(level.shape) for level in windowed._levels] arrays = [state(level.shape) for level in windowed._levels]
windowed.initialize(*arrays) windowed.initialize(*arrays)
plain.initialize(*arrays) plain.initialize(*arrays)
@@ -201,11 +200,11 @@ class ExchangeWindowShiftedMRTGPU(unittest.TestCase):
from CelerisLab.cuda.context import CudaContext from CelerisLab.cuda.context import CudaContext
cfg = LBMConfig(nx=32, ny=24, viscosity=.1, collision="MRT", ddf_shifting=True) cfg = LBMConfig(nx=32, ny=24, viscosity=.1, collision="MRT", ddf_shifting=True)
with CudaContext(0) as ctx: with CudaContext(0) as ctx:
kw = dict(context=ctx, refinements=[RB], level_blocks=[(12, 8), (16, 20)], kw = dict(context=ctx, base_block=(12, 8), fine_block=(16, 20), fine_objects=[rc()],
level_objects=[(), (rc(),)], observe_raw_exchange=True) observe_raw_exchange=True)
with FixedPatchStepper(cfg, observe_window_exchange=True, **kw) as windowed, \ with FixedPatchStepper(cfg, RB, observe_window_exchange=True, **kw) as windowed, \
FixedPatchStepper(cfg, observe_window_exchange=True, **kw) as split, \ FixedPatchStepper(cfg, RB, observe_window_exchange=True, **kw) as split, \
FixedPatchStepper(cfg, **kw) as plain: FixedPatchStepper(cfg, RB, **kw) as plain:
arrays = [state(level.shape) for level in windowed._levels] arrays = [state(level.shape) for level in windowed._levels]
for run in (windowed, split, plain): for run in (windowed, split, plain):
run.initialize(*arrays) run.initialize(*arrays)
@@ -140,8 +140,7 @@ def run_patch(speed, *, audit=True):
build = (patch.object(compiler, "compile_kernel", audit_compiler(compiler.compile_kernel)) build = (patch.object(compiler, "compile_kernel", audit_compiler(compiler.compile_kernel))
if audit else nullcontext()) if audit else nullcontext())
with build: with build:
backend = implementation.FixedPatchStepper(cfg, context=context, refinements=[BOUNDS], backend = implementation.FixedPatchStepper(cfg, BOUNDS, context=context, fine_objects=[source])
level_objects=[(), (source,)])
with backend, context.activate(): with backend, context.activate():
c, f = backend._levels c, f = backend._levels
cl = backend._fine_curved cl = backend._fine_curved
@@ -326,7 +325,6 @@ def test_cross_interface_donor_rejects_before_cuda():
context.activate.side_effect = AssertionError("must reject before CUDA activation") context.activate.side_effect = AssertionError("must reject before CUDA activation")
with patch.object(cuda, "mem_alloc") as allocate: with patch.object(cuda, "mem_alloc") as allocate:
with pytest.raises(ValueError, match="donor crosses"): with pytest.raises(ValueError, match="donor crosses"):
implementation.FixedPatchStepper(cfg, context=context, refinements=[BOUNDS], implementation.FixedPatchStepper(cfg, BOUNDS, context=context, fine_objects=[source])
level_objects=[(), (source,)])
context.activate.assert_not_called() context.activate.assert_not_called()
allocate.assert_not_called() allocate.assert_not_called()
@@ -48,8 +48,9 @@ def make(context):
for (x, y), amplitude in zip(centers, (.003, -.002)) for (x, y), amplitude in zip(centers, (.003, -.002))
) )
run = FixedPatchStepper( run = FixedPatchStepper(
CFG, context=context, refinements=[PARENTS, START], level_blocks=[(4, 4)] * 3, CFG, PARENTS, context=context, base_block=(4, 4), fine_block=(4, 4),
level_objects=[(), (), objects], level_jets=[(), (), jets], child_block=(4, 4), child_rectangles_in_parent=START,
child_objects=objects, child_jets=jets,
) )
return run, objects return run, objects
@@ -25,8 +25,9 @@ B = ((16,16,32,40), (60,16,76,40))
def make(ctx): def make(ctx):
body = SimObject(71, geometry=CircleGeometry(12.,14.,.8),center=(12.,14.),radius=.8) body = SimObject(71, geometry=CircleGeometry(12.,14.,.8),center=(12.,14.),radius=.8)
jet = Jet(Tube((10.5,14.),0.,3.,.3),Parallel(0.),'parabolic',.003) jet = Jet(Tube((10.5,14.),0.,3.,.3),Parallel(0.),'parabolic',.003)
run = FixedPatchStepper(CFG,context=ctx,refinements=[PARENTS,A],level_blocks=[(4,4)]*3, run = FixedPatchStepper(CFG,PARENTS,context=ctx,base_block=(4,4),fine_block=(4,4),
level_objects=[(),(),(body,)],level_jets=[(),(),(jet,)]) child_block=(4,4),child_rectangles_in_parent=A,
child_objects=(body,),child_jets=(jet,))
return run,(body,) return run,(body,)
@@ -90,18 +90,14 @@ def census(run, index, buffer):
return totals(dense, leaf), external, independent return totals(dense, leaf), external, independent
def fine(block):
"""Two-level circle+jet construction on BOUNDS; block None keeps the dense fine level."""
return dict(refinements=[BOUNDS], level_blocks=[None, block],
level_objects=[(), (circle(),)], level_jets=[(), (jet(),)])
class RawExchangeGPU(unittest.TestCase): class RawExchangeGPU(unittest.TestCase):
def test_row_ledger_interval_jet_switch_and_no_population_effect(self): def test_row_ledger_interval_jet_switch_and_no_population_effect(self):
with CudaContext(0) as ctx: with CudaContext(0) as ctx:
for block in ((32, 40), (16, 20), (8, 10)): for block in ((32, 40), (16, 20), (8, 10)):
with FixedPatchStepper(CFG, context=ctx, **fine(block), observe_raw_exchange=True) as measured, \ with FixedPatchStepper(CFG, BOUNDS, context=ctx, fine_block=block,
FixedPatchStepper(CFG, context=ctx, **fine(block)) as plain: fine_objects=[circle()], jets=[jet()], observe_raw_exchange=True) as measured, \
FixedPatchStepper(CFG, BOUNDS, context=ctx, fine_block=block,
fine_objects=[circle()], jets=[jet()]) as plain:
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)
@@ -150,8 +146,9 @@ class RawExchangeGPU(unittest.TestCase):
"""Same fine resolution/initial data with dense storage: a numerical """Same fine resolution/initial data with dense storage: a numerical
reference for cut-link rows, not a uniform-whole-domain convergence study.""" reference for cut-link rows, not a uniform-whole-domain convergence study."""
with CudaContext(0) as ctx: with CudaContext(0) as ctx:
with FixedPatchStepper(CFG, context=ctx, **fine(None)) as dense, \ with FixedPatchStepper(CFG, BOUNDS, context=ctx, fine_objects=[circle()], jets=[jet()]) as dense, \
FixedPatchStepper(CFG, context=ctx, **fine((16, 20)), observe_raw_exchange=True) as blocked: FixedPatchStepper(CFG, BOUNDS, context=ctx, fine_block=(16, 20),
fine_objects=[circle()], jets=[jet()], observe_raw_exchange=True) as blocked:
arrays = [state(l.shape) for l in blocked._levels] arrays = [state(l.shape) for l in blocked._levels]
dense.initialize(*arrays) dense.initialize(*arrays)
blocked.initialize(*arrays) blocked.initialize(*arrays)
@@ -232,7 +229,8 @@ class RawExchangeGPU(unittest.TestCase):
layout = BlockLayout2D.build(48, 64, 8, 8) layout = BlockLayout2D.build(48, 64, 8, 8)
with CudaContext(0) as ctx: with CudaContext(0) as ctx:
with PeriodicBlockStepper(fine_cfg, layout, context=ctx, objects=[obj]) as uniform, \ with PeriodicBlockStepper(fine_cfg, layout, context=ctx, objects=[obj]) as uniform, \
FixedPatchStepper(CFG, context=ctx, **fine((16, 20)), observe_raw_exchange=True) as amr: FixedPatchStepper(CFG, BOUNDS, context=ctx, fine_block=(16, 20),
fine_objects=[circle()], jets=[jet()], 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._level_boundaries[1]
@@ -83,10 +83,11 @@ def mocked_backend(mode="periodic", fail_alloc=None):
class BlockedBaseCPU(unittest.TestCase): class BlockedBaseCPU(unittest.TestCase):
def test_reject_before_cuda(self): def test_reject_before_cuda(self):
context = Mock(activate=Mock(side_effect=AssertionError("CUDA reached"))) context = Mock(activate=Mock(side_effect=AssertionError("CUDA reached")))
for blocks in ([(12, 8), None], [(7, 8), FINE_BLOCK], [(12,), FINE_BLOCK]): for kwargs in ({"base_block": (12, 8)},
with self.subTest(blocks=blocks), self.assertRaises(ValueError): {"base_block": (7, 8), "fine_block": FINE_BLOCK},
FixedPatchStepper(config("periodic"), context=context, refinements=[BOUNDS], {"base_block": (12,), "fine_block": FINE_BLOCK}):
level_blocks=blocks) with self.subTest(kwargs=kwargs), self.assertRaises(ValueError):
FixedPatchStepper(config("periodic"), BOUNDS, context=context, **kwargs)
context.activate.assert_not_called() context.activate.assert_not_called()
def test_all_arrival_rows_donors_and_kind_masks(self): def test_all_arrival_rows_donors_and_kind_masks(self):
@@ -192,10 +193,9 @@ class BlockedBaseCPU(unittest.TestCase):
def test_launch_order_stream_wall_constants_and_buffer_swaps(self): def test_launch_order_stream_wall_constants_and_buffer_swaps(self):
with mocked_backend("slot_reservoir") as m: with mocked_backend("slot_reservoir") as m:
backend = FixedPatchStepper(config("slot_reservoir"), context=m["context"], backend = FixedPatchStepper(config("slot_reservoir"), BOUNDS, context=m["context"],
coarse_boundary="slot_reservoir", refinements=[BOUNDS], coarse_boundary="slot_reservoir", fine_block=FINE_BLOCK,
level_blocks=[(3, 4), FINE_BLOCK], base_block=(3, 4), fine_objects=[circle(actuation(.002))])
level_objects=[(), (circle(actuation(.002)),)])
try: try:
backend.initialize(state(NY, NX, 1), state(*backend.patch.fine_shape, 2)) backend.initialize(state(NY, NX, 1), state(*backend.patch.fine_shape, 2))
coarse, fine = backend._levels coarse, fine = backend._levels
@@ -236,16 +236,16 @@ class BlockedBaseCPU(unittest.TestCase):
for fail_at in (1, 6, 10): for fail_at in (1, 6, 10):
with self.subTest(fail_at=fail_at), mocked_backend(fail_alloc=fail_at) as m: with self.subTest(fail_at=fail_at), mocked_backend(fail_alloc=fail_at) as m:
with self.assertRaisesRegex(MemoryError, "injected allocation"): with self.assertRaisesRegex(MemoryError, "injected allocation"):
FixedPatchStepper(config("periodic"), context=m["context"], refinements=[BOUNDS], FixedPatchStepper(config("periodic"), BOUNDS, context=m["context"],
level_blocks=[(3, 4), FINE_BLOCK]) fine_block=FINE_BLOCK, base_block=(3, 4))
self.assertEqual(len(m["allocations"]), fail_at) self.assertEqual(len(m["allocations"]), fail_at)
for allocation in m["allocations"]: for allocation in m["allocations"]:
allocation.free.assert_called_once_with() allocation.free.assert_called_once_with()
def test_failed_parent_stage_cannot_publish_or_resume(self): def test_failed_parent_stage_cannot_publish_or_resume(self):
with mocked_backend() as m: with mocked_backend() as m:
backend = FixedPatchStepper(config("periodic"), context=m["context"], refinements=[BOUNDS], backend = FixedPatchStepper(config("periodic"), BOUNDS, context=m["context"],
level_blocks=[(3, 4), FINE_BLOCK]) fine_block=FINE_BLOCK, base_block=(3, 4))
backend.initialize(state(NY, NX, 1), state(*backend.patch.fine_shape, 2)) backend.initialize(state(NY, NX, 1), state(*backend.patch.fine_shape, 2))
m["functions"]["BlockCoalesce"].side_effect = RuntimeError("injected parent stage failure") m["functions"]["BlockCoalesce"].side_effect = RuntimeError("injected parent stage failure")
current = backend._levels[0].current current = backend._levels[0].current
@@ -289,10 +289,9 @@ class BlockedBaseGPU(unittest.TestCase):
for base_block in (None, *layouts): for base_block in (None, *layouts):
label = f"{mode} width={width} speed={speed} base={base_block}" label = f"{mode} width={width} speed={speed} base={base_block}"
snapshots = [] snapshots = []
with FixedPatchStepper(cfg, context=context, coarse_boundary=mode, with FixedPatchStepper(cfg, BOUNDS, context=context, coarse_boundary=mode,
refinements=[BOUNDS], level_blocks=[base_block, FINE_BLOCK], fine_block=FINE_BLOCK, base_block=base_block,
level_objects=[(), () if speed is None fine_objects=() if speed is None else [circle(actuation(speed))]) as backend:
else (circle(actuation(speed)),)]) as backend:
backend.initialize(*arrays) backend.initialize(*arrays)
for step in (0, 1, 2, 8): for step in (0, 1, 2, 8):
if step == 1 and base_block is not None: if step == 1 and base_block is not None:
@@ -105,10 +105,9 @@ class BlockedFineCPU(unittest.TestCase):
return Mock(activate=Mock(side_effect=AssertionError("CUDA reached")), sm_arch="sm_70") return Mock(activate=Mock(side_effect=AssertionError("CUDA reached")), sm_arch="sm_70")
def test_rejections_before_cuda(self): def test_rejections_before_cuda(self):
for kw in ({"level_blocks": [None, (7, 7)]}, {"level_blocks": [None, (16,)]}): for kw in ({"fine_block": (7, 7)}, {"fine_block": (16,)}):
with self.subTest(kw=kw), self.assertRaises(ValueError): with self.subTest(kw=kw), self.assertRaises(ValueError):
FixedPatchStepper(config("periodic"), context=self.forbidden_context(), FixedPatchStepper(config("periodic"), BOUNDS, context=self.forbidden_context(), **kw)
refinements=[BOUNDS], **kw)
def test_single_block_layout_allocates_no_empty_halo_tables(self): def test_single_block_layout_allocates_no_empty_halo_tables(self):
# A 1x1 non-periodic fine layout has zero halo pairs; CUDA rejects mem_alloc(0). # A 1x1 non-periodic fine layout has zero halo pairs; CUDA rejects mem_alloc(0).
@@ -126,8 +125,8 @@ class BlockedFineCPU(unittest.TestCase):
patch.object(compiler, "compile_kernel"), patch.object(compiler, "cleanup_build_dir"), \ patch.object(compiler, "compile_kernel"), patch.object(compiler, "cleanup_build_dir"), \
patch.object(compiler, "load_module", return_value=module), \ patch.object(compiler, "load_module", return_value=module), \
patch.object(cuda, "mem_alloc", side_effect=alloc), patch.object(cuda, "memcpy_htod"): patch.object(cuda, "mem_alloc", side_effect=alloc), patch.object(cuda, "memcpy_htod"):
backend = FixedPatchStepper(config("periodic"), context=context, refinements=[BOUNDS], backend = FixedPatchStepper(config("periodic"), BOUNDS, context=context,
level_blocks=[None, FINE_BLOCKS[0]]) fine_block=FINE_BLOCKS[0])
self.assertEqual(len(backend.fine_layout.halo_dst), 0) self.assertEqual(len(backend.fine_layout.halo_dst), 0)
self.assertIsNone(backend._levels[1].halo_dst) self.assertIsNone(backend._levels[1].halo_dst)
backend._gather_halo(backend._levels[1], Mock()) # no launch without pairs backend._gather_halo(backend._levels[1], Mock()) # no launch without pairs
@@ -178,9 +177,8 @@ class BlockedFineCPU(unittest.TestCase):
patch.object(cuda, "mem_alloc", side_effect=lambda size: Mock()), \ patch.object(cuda, "mem_alloc", side_effect=lambda size: Mock()), \
patch.object(cuda, "memcpy_htod"), \ patch.object(cuda, "memcpy_htod"), \
patch.object(curved_links, "CurvedLinkSoA", side_effect=AssertionError("dense SoA built")): patch.object(curved_links, "CurvedLinkSoA", side_effect=AssertionError("dense SoA built")):
backend = FixedPatchStepper(config("periodic"), context=context, refinements=[BOUNDS], backend = FixedPatchStepper(config("periodic"), BOUNDS, context=context,
level_blocks=[None, FINE_BLOCKS[1]], fine_objects=objects, fine_block=FINE_BLOCKS[1])
level_objects=[(), tuple(objects)])
backend.initialize(state(NY, NX, 1), state(*backend.patch.fine_shape, 2)) backend.initialize(state(NY, NX, 1), state(*backend.patch.fine_shape, 2))
return backend, calls return backend, calls
@@ -291,9 +289,8 @@ class BlockedFineGPU(unittest.TestCase):
results = [] results = []
with CudaContext(0) as context: with CudaContext(0) as context:
for fine_block in (None, *FINE_BLOCKS): for fine_block in (None, *FINE_BLOCKS):
with FixedPatchStepper(cfg, context=context, coarse_boundary=mode, with FixedPatchStepper(cfg, BOUNDS, context=context, coarse_boundary=mode,
refinements=[BOUNDS], fine_block=fine_block) as backend:
level_blocks=[None, fine_block]) as backend:
backend.initialize(g_coarse, g_fine) backend.initialize(g_coarse, g_fine)
backend.step(STEPS) backend.step(STEPS)
self.assertEqual(backend.step_count, STEPS) self.assertEqual(backend.step_count, STEPS)
@@ -321,9 +318,8 @@ class BlockedFineGPU(unittest.TestCase):
results = [] results = []
with CudaContext(0) as context: with CudaContext(0) as context:
for fine_block in (None, *FINE_BLOCKS): for fine_block in (None, *FINE_BLOCKS):
with FixedPatchStepper(cfg, context=context, coarse_boundary=mode, with FixedPatchStepper(cfg, BOUNDS, context=context, coarse_boundary=mode,
refinements=[BOUNDS], level_blocks=[None, fine_block], fine_objects=[circle()], fine_block=fine_block) as backend:
level_objects=[(), (circle(),)]) as backend:
backend.initialize(g_coarse, g_fine) backend.initialize(g_coarse, g_fine)
backend.step(STEPS) backend.step(STEPS)
results.append((backend.download(), raw(backend))) results.append((backend.download(), raw(backend)))
@@ -350,9 +346,9 @@ class BlockedFineGPU(unittest.TestCase):
results = [] results = []
with CudaContext(0) as context: with CudaContext(0) as context:
for fine_block in (None, *FINE_BLOCKS): for fine_block in (None, *FINE_BLOCKS):
with FixedPatchStepper(cfg, context=context, coarse_boundary=mode, with FixedPatchStepper(cfg, BOUNDS, context=context, coarse_boundary=mode,
refinements=[BOUNDS], level_blocks=[None, fine_block], fine_objects=[circle(actuation(speed))],
level_objects=[(), (circle(actuation(speed)),)]) as backend: fine_block=fine_block) as backend:
backend.initialize(g_coarse, g_fine) backend.initialize(g_coarse, g_fine)
backend.step(STEPS) backend.step(STEPS)
results.append((backend.download(), raw(backend))) results.append((backend.download(), raw(backend)))
@@ -94,8 +94,7 @@ class FaultGuardsCPU(unittest.TestCase):
("slot_reservoir", (12, 4, 20, 21))): ("slot_reservoir", (12, 4, 20, 21))):
with self.subTest(mode=mode, bounds=bounds): with self.subTest(mode=mode, bounds=bounds):
with self.assertRaises(ValueError): with self.assertRaises(ValueError):
implementation.FixedPatchStepper(cfg, context=context, coarse_boundary=mode, implementation.FixedPatchStepper(cfg, bounds, context=context, coarse_boundary=mode)
refinements=[bounds])
context.activate.assert_not_called() context.activate.assert_not_called()
# Both exact channel limits and the original periodic two-cell margin # Both exact channel limits and the original periodic two-cell margin
# must reach context activation, with no allocation or device required. # must reach context activation, with no allocation or device required.
@@ -103,8 +102,7 @@ class FaultGuardsCPU(unittest.TestCase):
("periodic", (2, 2, 6, 5))): ("periodic", (2, 2, 6, 5))):
with self.subTest(mode=mode, bounds=bounds): with self.subTest(mode=mode, bounds=bounds):
with self.assertRaisesRegex(RuntimeError, "validated before CUDA"): with self.assertRaisesRegex(RuntimeError, "validated before CUDA"):
implementation.FixedPatchStepper(cfg, context=context, coarse_boundary=mode, implementation.FixedPatchStepper(cfg, bounds, context=context, coarse_boundary=mode)
refinements=[bounds])
def fake_backend(self): def fake_backend(self):
backend = PeriodicPatchStepper.__new__(PeriodicPatchStepper) backend = PeriodicPatchStepper.__new__(PeriodicPatchStepper)
@@ -156,12 +154,12 @@ class FaultGuardsCPU(unittest.TestCase):
context = Mock(activate=nullcontext, sm_arch="sm_70") context = Mock(activate=nullcontext, sm_arch="sm_70")
with patch.object(implementation, "HALO", 2): with patch.object(implementation, "HALO", 2):
with self.assertRaisesRegex(RuntimeError, "four-cell"): with self.assertRaisesRegex(RuntimeError, "four-cell"):
PeriodicPatchStepper(cfg, context=context, refinements=[(4, 4, 12, 12)]) PeriodicPatchStepper(cfg, (4, 4, 12, 12), context=context)
p = implementation.FixedPatch2D.from_config(cfg, (4, 4, 12, 12)) p = implementation.FixedPatch2D.from_config(cfg, (4, 4, 12, 12))
with patch.object(cuda, "Stream"), patch.object(implementation.FixedPatch2D, "from_config", with patch.object(cuda, "Stream"), patch.object(implementation.FixedPatch2D, "from_config",
return_value=replace(p, tau_coarse=.9)): return_value=replace(p, tau_coarse=.9)):
with self.assertRaisesRegex(RuntimeError, "omega/tau"): with self.assertRaisesRegex(RuntimeError, "omega/tau"):
PeriodicPatchStepper(cfg, context=context, refinements=[p.bounds]) PeriodicPatchStepper(cfg, p.bounds, context=context)
module = Mock() module = Mock()
module.get_global.return_value = (Mock(), 16) module.get_global.return_value = (Mock(), 16)
with patch.object(cuda, "Stream"), patch.object(compiler, "generate_config"), \ with patch.object(cuda, "Stream"), patch.object(compiler, "generate_config"), \
@@ -169,7 +167,7 @@ class FaultGuardsCPU(unittest.TestCase):
patch.object(compiler, "load_module", return_value=module), \ patch.object(compiler, "load_module", return_value=module), \
patch.object(cuda, "memcpy_htod") as upload: patch.object(cuda, "memcpy_htod") as upload:
with self.assertRaisesRegex(RuntimeError, "layout mismatch"): with self.assertRaisesRegex(RuntimeError, "layout mismatch"):
PeriodicPatchStepper(cfg, context=context, refinements=[p.bounds]) PeriodicPatchStepper(cfg, p.bounds, context=context)
upload.assert_not_called() upload.assert_not_called()
@@ -243,7 +241,7 @@ class NormalBuildGPU(unittest.TestCase):
max_error = max_bc_error = 0.0 max_error = max_bc_error = 0.0
with CudaContext(0) as context: with CudaContext(0) as context:
with patch.object(compiler, "generate_config", side_effect=record_config): with patch.object(compiler, "generate_config", side_effect=record_config):
backend = implementation.FixedPatchStepper(cfg, context=context, refinements=[bounds], backend = implementation.FixedPatchStepper(cfg, bounds, context=context,
coarse_boundary="slot_reservoir") coarse_boundary="slot_reservoir")
with backend: with backend:
self.assertEqual(len(generated), 2) self.assertEqual(len(generated), 2)
@@ -313,8 +311,7 @@ class NormalBuildGPU(unittest.TestCase):
for a in (initial.coarse_f, initial.fine_f))) for a in (initial.coarse_f, initial.fine_f)))
max_pop = max_drift = max_moment = 0.0 max_pop = max_drift = max_moment = 0.0
with CudaContext(0) as context, PeriodicPatchStepper( with CudaContext(0) as context, PeriodicPatchStepper(
LBMConfig(nx=16, ny=16, viscosity=.1), context=context, LBMConfig(nx=16, ny=16, viscosity=.1), (4, 4, 12, 12), context=context) as backend:
refinements=[(4, 4, 12, 12)]) as backend:
backend.initialize(*postcollision(state, spec, maps)) backend.initialize(*postcollision(state, spec, maps))
before = totals(backend.download(), maps) before = totals(backend.download(), maps)
for step in range(1, 9): for step in range(1, 9):
@@ -345,7 +342,7 @@ class NormalBuildGPU(unittest.TestCase):
fine_mask[4:-4, 4:-4] = True fine_mask[4:-4, 4:-4] = True
initial[1][:, ~fine_mask] = np.nan initial[1][:, ~fine_mask] = np.nan
expected = scalar_interval(*[a.astype(np.float64) for a in initial], bounds) expected = scalar_interval(*[a.astype(np.float64) for a in initial], bounds)
with CudaContext(0) as context, PeriodicPatchStepper(cfg, context=context, refinements=[bounds]) as backend: with CudaContext(0) as context, PeriodicPatchStepper(cfg, bounds, context=context) as backend:
self.assertEqual(backend.patch.fine_shape, (14, 16)) self.assertEqual(backend.patch.fine_shape, (14, 16))
p = backend.patch p = backend.patch
allocated = 2 * 9 * (17*13+16*14) * 4 + sum(a.nbytes for a in allocated = 2 * 9 * (17*13+16*14) * 4 + sum(a.nbytes for a in
@@ -84,13 +84,13 @@ class CoupledGPU(unittest.TestCase):
cases = (("dense_unshifted", False, {}), cases = (("dense_unshifted", False, {}),
("blocked_shifted", True, ("blocked_shifted", True,
dict(level_blocks=[(4, 4), (6, 6)]))) dict(base_block=(4, 4), fine_block=(6, 6))))
with CudaContext(0) as context: with CudaContext(0) as context:
for label, shifted, storage in cases: for label, shifted, storage in cases:
with self.subTest(storage=label): with self.subTest(storage=label):
cfg = LBMConfig(nx=16, ny=16, viscosity=.0192, cfg = LBMConfig(nx=16, ny=16, viscosity=.0192,
ddf_shifting=shifted) ddf_shifting=shifted)
with FixedPatchStepper(cfg, context=context, refinements=[(4, 4, 12, 12)], with FixedPatchStepper(cfg, (4, 4, 12, 12), context=context,
**storage) as backend: **storage) as backend:
self.assertAlmostEqual(backend.patch.tau_coarse, spec.tau_coarse) self.assertAlmostEqual(backend.patch.tau_coarse, spec.tau_coarse)
self.assertAlmostEqual(backend.patch.tau_fine, spec.tau_fine) self.assertAlmostEqual(backend.patch.tau_fine, spec.tau_fine)
@@ -148,7 +148,7 @@ class CoupledGPU(unittest.TestCase):
maps = ref.build_maps(spec) maps = ref.build_maps(spec)
with CudaContext(0) as context: with CudaContext(0) as context:
with PeriodicPatchStepper(LBMConfig(nx=16, ny=16, viscosity=.1), with PeriodicPatchStepper(LBMConfig(nx=16, ny=16, viscosity=.1),
context=context, refinements=[(4, 4, 12, 12)]) as backend: (4, 4, 12, 12), context=context) as backend:
backend.initialize(*postcollision(ref.make_state(spec), spec, maps)) backend.initialize(*postcollision(ref.make_state(spec), spec, maps))
coarse, fine = backend._levels coarse, fine = backend._levels
@@ -187,7 +187,7 @@ class CoupledGPU(unittest.TestCase):
with CudaContext(0) as context: with CudaContext(0) as context:
with patch.object(compiler, "compile_kernel", audit_compile): with patch.object(compiler, "compile_kernel", audit_compile):
backend = PeriodicPatchStepper(cfg, context=context, refinements=[(4, 4, 12, 12)]) backend = PeriodicPatchStepper(cfg, (4, 4, 12, 12), context=context)
with backend: with backend:
# Production allocations: four DDF arrays and four compact tables. # Production allocations: four DDF arrays and four compact tables.
device_bytes = 2 * 9 * (16**2 + 24**2) * 4 device_bytes = 2 * 9 * (16**2 + 24**2) * 4
@@ -270,7 +270,7 @@ class CoupledGPU(unittest.TestCase):
backend.download() backend.download()
# Compile and execute the actual non-audited entry as well: one step. # Compile and execute the actual non-audited entry as well: one step.
with PeriodicPatchStepper(cfg, context=context, refinements=[(4, 4, 12, 12)]) as plain: with PeriodicPatchStepper(cfg, (4, 4, 12, 12), context=context) as plain:
rest = postcollision(ref.make_state(spec), spec, maps) rest = postcollision(ref.make_state(spec), spec, maps)
plain.initialize(*rest) plain.initialize(*rest)
plain.step() plain.step()
@@ -41,9 +41,8 @@ def mocked(objects, jets, fine_block):
patch.object(compiler, "load_module", return_value=module), \ patch.object(compiler, "load_module", return_value=module), \
patch.object(cuda, "mem_alloc", side_effect=lambda size: Mock()), \ patch.object(cuda, "mem_alloc", side_effect=lambda size: Mock()), \
patch.object(cuda, "memcpy_htod", htod): patch.object(cuda, "memcpy_htod", htod):
backend = FixedPatchStepper(config("periodic"), context=context, refinements=[BOUNDS], backend = FixedPatchStepper(config("periodic"), BOUNDS, context=context,
level_blocks=[None, fine_block], level_objects=[(), tuple(objects)], fine_objects=objects, fine_block=fine_block, jets=jets)
level_jets=[(), tuple(jets)])
return backend, htod return backend, htod
@@ -150,9 +149,8 @@ class JetStepperCPU(unittest.TestCase):
def test_jet_overlapping_legacy_aperture_rejected_before_cuda(self): def test_jet_overlapping_legacy_aperture_rejected_before_cuda(self):
with patch.object(cuda, "Stream", side_effect=AssertionError("CUDA touched")): with patch.object(cuda, "Stream", side_effect=AssertionError("CUDA touched")):
with self.assertRaisesRegex(ValueError, "SurfaceActuation"): with self.assertRaisesRegex(ValueError, "SurfaceActuation"):
FixedPatchStepper(config("periodic"), context=Mock(), refinements=[BOUNDS], FixedPatchStepper(config("periodic"), BOUNDS, context=Mock(),
level_objects=[(), (circle(actuation(SPEED)),)], fine_objects=[circle(actuation(SPEED))], jets=[jet()])
level_jets=[(), (jet(),)])
def run(mode, fine_block, *, objects, jets, switch=None): def run(mode, fine_block, *, objects, jets, switch=None):
@@ -160,9 +158,8 @@ def run(mode, fine_block, *, objects, jets, switch=None):
cfg = config(mode) cfg = config(mode)
patch_ = FixedPatch2D.from_config(cfg, BOUNDS) patch_ = FixedPatch2D.from_config(cfg, BOUNDS)
with CudaContext(0) as context, FixedPatchStepper( with CudaContext(0) as context, FixedPatchStepper(
cfg, context=context, coarse_boundary=mode, refinements=[BOUNDS], cfg, BOUNDS, context=context, coarse_boundary=mode, fine_objects=objects,
level_blocks=[None, fine_block], level_objects=[(), tuple(objects)], fine_block=fine_block, jets=jets) as backend:
level_jets=[(), tuple(jets)]) as backend:
backend.initialize(state(NY, NX, 1), state(*patch_.fine_shape, 2)) backend.initialize(state(NY, NX, 1), state(*patch_.fine_shape, 2))
if switch is None: if switch is None:
backend.step(STEPS) backend.step(STEPS)
@@ -42,9 +42,8 @@ class RefinementMRTGPU(unittest.TestCase):
collision="MRT", ddf_shifting=shifted) collision="MRT", ddf_shifting=shifted)
reference = None reference = None
for block in ((32, 48), (8, 8), (4, 1)): for block in ((32, 48), (8, 8), (4, 1)):
with FixedPatchStepper(cfg, context=context, coarse_boundary=boundary, with FixedPatchStepper(cfg, (8, 6, 40, 26), context=context, coarse_boundary=boundary,
refinements=[(8, 6, 40, 26)], level_blocks=[block, (8, 8)], base_block=block, fine_block=(8, 8), fine_objects=(body,)) as backend:
level_objects=[(), (body,)]) as backend:
backend.initialize(*[state(*level.shape[1:]) for level in backend._levels]) backend.initialize(*[state(*level.shape[1:]) for level in backend._levels])
backend.step(8) backend.step(8)
raw = backend._download_raw() raw = backend._download_raw()
@@ -64,8 +63,8 @@ class RefinementMRTGPU(unittest.TestCase):
results = {} results = {}
layouts = ((None, None), *((base, FINE_BLOCK) for base in (None, *BASE_BLOCKS))) layouts = ((None, None), *((base, FINE_BLOCK) for base in (None, *BASE_BLOCKS)))
for base_block, fine_block in layouts: for base_block, fine_block in layouts:
with FixedPatchStepper(cfg, context=context, refinements=[BOUNDS], with FixedPatchStepper(cfg, BOUNDS, context=context, fine_block=fine_block,
level_blocks=[base_block, fine_block]) as backend: base_block=base_block) as backend:
backend.initialize(*[fixture.state(*level.shape[1:], seed) for seed, level in enumerate(backend._levels)]) backend.initialize(*[fixture.state(*level.shape[1:], seed) for seed, level in enumerate(backend._levels)])
backend.step(4) backend.step(4)
results[base_block, fine_block] = raw(backend) results[base_block, fine_block] = raw(backend)
@@ -76,9 +75,9 @@ class RefinementMRTGPU(unittest.TestCase):
def test_shifted_all_blocked_three_level_mrt_smoke(self): def test_shifted_all_blocked_three_level_mrt_smoke(self):
"""Per-level shear rate omega_l = 1/(3*2**l*nu + 1/2); finite ACTIVE leaves.""" """Per-level shear rate omega_l = 1/(3*2**l*nu + 1/2); finite ACTIVE leaves."""
cfg = LBMConfig(nx=16, ny=16, viscosity=.1, collision="MRT", ddf_shifting=True) cfg = LBMConfig(nx=16, ny=16, viscosity=.1, collision="MRT", ddf_shifting=True)
kw = dict(refinements=[(4, 4, 12, 12), (8, 8, 16, 16)], level_blocks=[(4, 4)] * 3) kw = dict(fine_block=(4, 4), base_block=(4, 4), child=(8, 8, 16, 16), child_block=(4, 4))
with CudaContext(0) as context: with CudaContext(0) as context:
with FixedPatchStepper(cfg, context=context, **kw) as backend: with FixedPatchStepper(cfg, (4, 4, 12, 12), context=context, **kw) as backend:
self.assertEqual(len(backend._level_cfgs), 3) self.assertEqual(len(backend._level_cfgs), 3)
for level, level_cfg in enumerate(backend._level_cfgs): for level, level_cfg in enumerate(backend._level_cfgs):
checked = replace(level_cfg) checked = replace(level_cfg)
@@ -26,12 +26,6 @@ FACE = ((4, 4, 8, 12), (12, 4, 16, 12))
S1 = ((4, 4, 12, 12), (20, 4, 28, 12)) S1 = ((4, 4, 12, 12), (20, 4, 28, 12))
def levels(bounds, blocks=(None, None), objects=(), jets=(), **rest):
"""Two-level constructor keywords: one L1 rectangle or set, per-level blocks, L1 bodies/Jets."""
return dict(refinements=[bounds], level_blocks=list(blocks),
level_objects=[(), tuple(objects)], level_jets=[(), tuple(jets)], **rest)
def level_fine(patch_set, windows): def level_fine(patch_set, windows):
"""Level-wide fine array with each patch window filled from `windows`.""" """Level-wide fine array with each patch window filled from `windows`."""
out = np.full((9, *patch_set.fine_shape), np.nan, np.float32) out = np.full((9, *patch_set.fine_shape), np.nan, np.float32)
@@ -51,21 +45,21 @@ class MultipatchCPU(unittest.TestCase):
def test_rejections_before_cuda(self): def test_rejections_before_cuda(self):
cfg = config("periodic") cfg = config("periodic")
context = object() # never activated: every case must fail first context = object() # never activated: every case must fail first
for kw in (dict(), dict(blocks=(None, (8, 4))), dict(blocks=((4, 4), None)), for kw in (dict(), dict(fine_block=(8, 4)), dict(base_block=(4, 4)),
dict(blocks=((4, 4), (8, 4)), objects=[circle()]), dict(fine_block=(8, 4), base_block=(4, 4), fine_objects=[circle()]),
dict(blocks=((4, 4), (5, 4)))): dict(fine_block=(5, 4), base_block=(4, 4))):
with self.subTest(kw=kw), self.assertRaises(ValueError): with self.subTest(kw=kw), self.assertRaises(ValueError):
FixedPatchStepper(cfg, context=context, **levels(((4, 4, 8, 12), (14, 4, 18, 12)), **kw)) FixedPatchStepper(cfg, ((4, 4, 8, 12), (14, 4, 18, 12)), context=context, **kw)
with self.assertRaisesRegex(ValueError, "margin"): with self.assertRaisesRegex(ValueError, "margin"):
FixedPatchStepper(cfg, context=context, FixedPatchStepper(cfg, ((1, 4, 8, 12), (10, 4, 14, 12)), context=context,
**levels(((1, 4, 8, 12), (10, 4, 14, 12)), blocks=((4, 4), (8, 4)))) fine_block=(8, 4), base_block=(4, 4))
from CelerisLab.body.jets import Jet, Parallel, Tube from CelerisLab.body.jets import Jet, Parallel, Tube
cfg32 = LBMConfig(nx=32, ny=16, viscosity=.1) cfg32 = LBMConfig(nx=32, ny=16, viscosity=.1)
span = Jet(Tube((5., 8.), 0., 22., .6), Parallel(0.), amplitude=0.) span = Jet(Tube((5., 8.), 0., 22., .6), Parallel(0.), amplitude=0.)
with self.assertRaisesRegex(ValueError, "one jet"): # same obj_id cannot bypass it with self.assertRaisesRegex(ValueError, "one jet"): # same obj_id cannot bypass it
FixedPatchStepper(cfg32, context=context, FixedPatchStepper(cfg32, ((2, 2, 14, 14), (18, 2, 30, 14)), context=context,
**levels(((2, 2, 14, 14), (18, 2, 30, 14)), blocks=((4, 4), (4, 4)), fine_block=(4, 4), base_block=(4, 4),
objects=[body(8., 8.), body(24., 8.)], jets=[span])) fine_objects=[body(8., 8.), body(24., 8.)], jets=[span])
class MultipatchGPU(unittest.TestCase): class MultipatchGPU(unittest.TestCase):
@@ -73,7 +67,7 @@ class MultipatchGPU(unittest.TestCase):
def run_steps(self, context, cfg, bounds, coarse, fine, **kw): def run_steps(self, context, cfg, bounds, coarse, fine, **kw):
out = [] out = []
with FixedPatchStepper(cfg, context=context, **levels(bounds, **kw)) as backend: with FixedPatchStepper(cfg, bounds, context=context, **kw) as backend:
backend.initialize(coarse, fine) backend.initialize(coarse, fine)
for step in self.STEPS: for step in self.STEPS:
backend.step(step - backend.step_count) backend.step(step - backend.step_count)
@@ -88,7 +82,7 @@ class MultipatchGPU(unittest.TestCase):
p = FixedPatch2D.from_config(cfg, BOUNDS) p = FixedPatch2D.from_config(cfg, BOUNDS)
s = FixedPatchSet2D.from_config(cfg, (BOUNDS,)) s = FixedPatchSet2D.from_config(cfg, (BOUNDS,))
coarse, fine = state(NY, NX, 1), state(*p.fine_shape, 2) coarse, fine = state(NY, NX, 1), state(*p.fine_shape, 2)
kw = dict(coarse_boundary=mode, blocks=((3, 4), (8, 4))) kw = dict(coarse_boundary=mode, fine_block=(8, 4), base_block=(3, 4))
old = self.run_steps(context, cfg, BOUNDS, coarse, fine, **kw) old = self.run_steps(context, cfg, BOUNDS, coarse, fine, **kw)
new = self.run_steps(context, cfg, (BOUNDS,), coarse, level_fine(s, [fine]), **kw) new = self.run_steps(context, cfg, (BOUNDS,), coarse, level_fine(s, [fine]), **kw)
for step, a, b in zip(self.STEPS, old, new): for step, a, b in zip(self.STEPS, old, new):
@@ -105,7 +99,7 @@ class MultipatchGPU(unittest.TestCase):
one = FixedPatchSet2D.from_config(cfg16, (S1[0],)) one = FixedPatchSet2D.from_config(cfg16, (S1[0],))
two = FixedPatchSet2D.from_config(cfg32, S1) two = FixedPatchSet2D.from_config(cfg32, S1)
coarse, fine = state(16, 16, 3), state(*one.patches[0].fine_shape, 4) coarse, fine = state(16, 16, 3), state(*one.patches[0].fine_shape, 4)
kw = dict(blocks=((4, 4), (4, 4))) kw = dict(fine_block=(4, 4), base_block=(4, 4))
with CudaContext(0) as context: with CudaContext(0) as context:
ref = self.run_steps(context, cfg16, (S1[0],), coarse, level_fine(one, [fine]), **kw) ref = self.run_steps(context, cfg16, (S1[0],), coarse, level_fine(one, [fine]), **kw)
tiled = np.concatenate([coarse, coarse], axis=2) tiled = np.concatenate([coarse, coarse], axis=2)
@@ -124,7 +118,7 @@ class MultipatchGPU(unittest.TestCase):
s = FixedPatchSet2D.from_config(cfg, FACE) s = FixedPatchSet2D.from_config(cfg, FACE)
coarse = state(16, 24, 5) coarse = state(16, 24, 5)
fine = level_fine(s, [state(*p.fine_shape, 6 + i) for i, p in enumerate(s.patches)]) fine = level_fine(s, [state(*p.fine_shape, 6 + i) for i, p in enumerate(s.patches)])
kw = dict(blocks=((4, 4), (4, 4))) kw = dict(fine_block=(4, 4), base_block=(4, 4))
leaf = s.fine_kind == ACTIVE leaf = s.fine_kind == ACTIVE
support = s.fine_layout((4, 4)).owned support = s.fine_layout((4, 4)).owned
outer = support & (s.fine_kind == OUTER_GHOST) outer = support & (s.fine_kind == OUTER_GHOST)
@@ -170,7 +164,7 @@ class MultipatchBodiesGPU(unittest.TestCase):
def run_steps(self, context, cfg, bounds, coarse, fine, switch=None, **kw): def run_steps(self, context, cfg, bounds, coarse, fine, switch=None, **kw):
out = [] out = []
with FixedPatchStepper(cfg, context=context, **levels(bounds, **kw)) as backend: with FixedPatchStepper(cfg, bounds, context=context, **kw) as backend:
backend.initialize(coarse, fine) backend.initialize(coarse, fine)
for step in self.STEPS: for step in self.STEPS:
if switch is not None and step == 8: if switch is not None and step == 8:
@@ -188,7 +182,8 @@ class MultipatchBodiesGPU(unittest.TestCase):
p = FixedPatch2D.from_config(cfg, BOUNDS) p = FixedPatch2D.from_config(cfg, BOUNDS)
s = FixedPatchSet2D.from_config(cfg, (BOUNDS,)) s = FixedPatchSet2D.from_config(cfg, (BOUNDS,))
coarse, fine = state(NY, NX, 1), state(*p.fine_shape, 2) coarse, fine = state(NY, NX, 1), state(*p.fine_shape, 2)
kw = dict(coarse_boundary=mode, blocks=((3, 4), (8, 4)), objects=[circle()]) kw = dict(coarse_boundary=mode, fine_block=(8, 4), base_block=(3, 4),
fine_objects=[circle()])
old = self.run_steps(context, cfg, BOUNDS, coarse, fine, **kw) old = self.run_steps(context, cfg, BOUNDS, coarse, fine, **kw)
new = self.run_steps(context, cfg, (BOUNDS,), coarse, level_fine(s, [fine]), **kw) new = self.run_steps(context, cfg, (BOUNDS,), coarse, level_fine(s, [fine]), **kw)
for step, a, b in zip(self.STEPS, old, new): for step, a, b in zip(self.STEPS, old, new):
@@ -207,17 +202,17 @@ class MultipatchBodiesGPU(unittest.TestCase):
two = FixedPatchSet2D.from_config(cfg32, self.PATCHES) two = FixedPatchSet2D.from_config(cfg32, self.PATCHES)
coarse, fine = state(16, 16, 3), state(*one.patches[0].fine_shape, 4) coarse, fine = state(16, 16, 3), state(*one.patches[0].fine_shape, 4)
jet = lambda x: Jet(Tube((x, 8.), 0., 6., .6), Parallel(0.), amplitude=.002) jet = lambda x: Jet(Tube((x, 8.), 0., 6., .6), Parallel(0.), amplitude=.002)
kw = dict(blocks=((4, 4), (4, 4))) kw = dict(fine_block=(4, 4), base_block=(4, 4))
with CudaContext(0) as context: with CudaContext(0) as context:
for jets in ((), (jet(5.),)): for jets in ((), (jet(5.),)):
ref = self.run_steps(context, cfg16, self.PATCHES[:1], coarse, level_fine(one, [fine]), ref = self.run_steps(context, cfg16, self.PATCHES[:1], coarse, level_fine(one, [fine]),
switch=[(0, -.001)] if jets else None, switch=[(0, -.001)] if jets else None,
objects=[body(8., 8.)], jets=jets, **kw) fine_objects=[body(8., 8.)], jets=jets, **kw)
both = jets and (jet(5.), jet(21.)) both = jets and (jet(5.), jet(21.))
got = self.run_steps(context, cfg32, self.PATCHES, np.concatenate([coarse, coarse], 2), got = self.run_steps(context, cfg32, self.PATCHES, np.concatenate([coarse, coarse], 2),
level_fine(two, [fine, fine]), level_fine(two, [fine, fine]),
switch=[(0, -.001), (1, -.001)] if jets else None, switch=[(0, -.001), (1, -.001)] if jets else None,
objects=[body(8., 8., 5), body(24., 8., 5)], fine_objects=[body(8., 8., 5), body(24., 8., 5)],
jets=both or (), **kw) jets=both or (), **kw)
for step, a, b in zip(self.STEPS, ref, got): for step, a, b in zip(self.STEPS, ref, got):
label = f"jets={len(jets)} step {step}" label = f"jets={len(jets)} step {step}"
@@ -27,16 +27,7 @@ from test_refinement_blocked_base import raw
CFG = LBMConfig(nx=16, ny=16, viscosity=.1) CFG = LBMConfig(nx=16, ny=16, viscosity=.1)
L1, L2 = (4, 4, 12, 12), (8, 8, 16, 16) L1, L2 = (4, 4, 12, 12), (8, 8, 16, 16)
B = (4, 4) KW = dict(fine_block=(4, 4), base_block=(4, 4), child=L2, child_block=(4, 4))
def levels(l1=L1, l2=L2, blocks=(B, B, B), objects=((), (), ()), jets=((), (), ())):
"""Level-list keywords for the scalar L1 patch with one scalar L2 child."""
return dict(refinements=[l1, l2], level_blocks=list(blocks),
level_objects=list(objects), level_jets=list(jets))
KW = levels()
def child_body(spec=None): def child_body(spec=None):
@@ -61,29 +52,26 @@ def mocked_child():
patch.object(compiler, "load_module", return_value=module), \ patch.object(compiler, "load_module", return_value=module), \
patch.object(cuda, "mem_alloc", side_effect=lambda size: Mock()), \ patch.object(cuda, "mem_alloc", side_effect=lambda size: Mock()), \
patch.object(cuda, "memcpy_htod"): patch.object(cuda, "memcpy_htod"):
return FixedPatchStepper(CFG, context=context, return FixedPatchStepper(CFG, L1, context=context, child_objects=[child_body()],
**levels(objects=((), (), (child_body(),)), jets=((), (), (child_jet(),)))) child_jets=[child_jet()], **KW)
class ThreeLevelCPU(unittest.TestCase): class ThreeLevelCPU(unittest.TestCase):
def test_rejections_before_cuda(self): def test_rejections_before_cuda(self):
context = object() # never activated: every case must fail first context = object() # never activated: every case must fail first
bad = (levels(blocks=(B, B, None)), levels(blocks=(B, None, B)), bad = (dict(KW, child_block=None), dict(KW, fine_block=None),
levels(objects=((), (circle(),), ())), levels(l2=(4, 8, 16, 16)), dict(KW, fine_objects=[circle()]), dict(KW, child=(4, 8, 16, 16)),
levels(blocks=(B, B, (5, 4)))) dict(KW, child_block=(5, 4)))
for kw in bad: for kw in bad:
with self.subTest(kw=kw), self.assertRaises(ValueError): with self.subTest(kw=kw), self.assertRaises(ValueError):
FixedPatchStepper(CFG, context=context, **kw) FixedPatchStepper(CFG, L1, context=context, **kw)
# L2 bodies/Jets without an L2 level (one entry too many), and an L2 Jet without a body. for kw in (dict(child_objects=[child_body()]), dict(child_jets=[child_jet()]),
two_level = dict(refinements=[L1], level_blocks=[None, None]) dict(KW, child_jets=[child_jet()])):
for kw in (dict(two_level, level_objects=[(), (), (child_body(),)]),
dict(two_level, level_jets=[(), (), (child_jet(),)]),
levels(jets=((), (), (child_jet(),)))):
with self.assertRaises(ValueError): with self.assertRaises(ValueError):
FixedPatchStepper(CFG, context=context, **kw) FixedPatchStepper(CFG, L1, context=context, **kw)
with self.assertRaises(ValueError): # several L1 patches + scalar child: not yet with self.assertRaises(ValueError): # several L1 patches + child: not yet
FixedPatchStepper(LBMConfig(nx=32, ny=16, viscosity=.1), context=context, FixedPatchStepper(LBMConfig(nx=32, ny=16, viscosity=.1), ((4, 4, 12, 12), (20, 4, 28, 12)),
**levels(l1=((4, 4, 12, 12), (20, 4, 28, 12)))) context=context, **KW)
def test_four_halo_bfl_steps_and_fresh_middle_feeds(self): def test_four_halo_bfl_steps_and_fresh_middle_feeds(self):
@@ -168,7 +156,7 @@ class ThreeLevelGPU(unittest.TestCase):
def test_equilibrium_is_kept_on_every_level(self): def test_equilibrium_is_kept_on_every_level(self):
from CelerisLab.cuda.context import CudaContext from CelerisLab.cuda.context import CudaContext
feq = np.asarray(compute_feq_d2q9(1.0, .03, -.02), np.float32)[:, None, None] feq = np.asarray(compute_feq_d2q9(1.0, .03, -.02), np.float32)[:, None, None]
with CudaContext(0) as ctx, FixedPatchStepper(CFG, context=ctx, **KW) as run: with CudaContext(0) as ctx, FixedPatchStepper(CFG, L1, context=ctx, **KW) as run:
shapes = ((16, 16), run.patch.fine_shape, run.child_patch.fine_shape) shapes = ((16, 16), run.patch.fine_shape, run.child_patch.fine_shape)
run.initialize(*(np.broadcast_to(feq, (9, *s)).copy() for s in shapes)) run.initialize(*(np.broadcast_to(feq, (9, *s)).copy() for s in shapes))
run.step(8) run.step(8)
@@ -181,7 +169,7 @@ class ThreeLevelGPU(unittest.TestCase):
def test_nonequilibrium_masks_and_fresh_middle_feeds(self): def test_nonequilibrium_masks_and_fresh_middle_feeds(self):
from CelerisLab.cuda.context import CudaContext from CelerisLab.cuda.context import CudaContext
with CudaContext(0) as ctx, FixedPatchStepper(CFG, context=ctx, **KW) as run: with CudaContext(0) as ctx, FixedPatchStepper(CFG, L1, context=ctx, **KW) as run:
p1, p2 = run.patch, run.child_patch p1, p2 = run.patch, run.child_patch
run.initialize(state(16, 16, 1), state(*p1.fine_shape, 2), state(*p2.fine_shape, 3)) run.initialize(state(16, 16, 1), state(*p1.fine_shape, 2), state(*p2.fine_shape, 3))
f = run._levels[1] f = run._levels[1]
@@ -217,9 +205,8 @@ class ThreeLevelBodyGPU(unittest.TestCase):
CFG, velocity=0., wall_normal_velocity=.001, wall_slot_width=4, CFG, velocity=0., wall_normal_velocity=.001, wall_slot_width=4,
y_wall_bc="on_site_velocity", top_wall_bc="reservoir") y_wall_bc="on_site_velocity", top_wall_bc="reservoir")
with CudaContext(0) as ctx, FixedPatchStepper( with CudaContext(0) as ctx, FixedPatchStepper(
cfg, context=ctx, coarse_boundary=mode, cfg, L1, context=ctx, coarse_boundary=mode, child_objects=[child_body(spec)],
**levels(blocks=(B, B, block), objects=((), (), (child_body(spec),)), child_jets=jets, **dict(KW, child_block=block)) as run:
jets=((), (), tuple(jets)))) as run:
run.initialize(*(state(*level.shape[1:], i + 1) for i, level in enumerate(run._levels))) run.initialize(*(state(*level.shape[1:], i + 1) for i, level in enumerate(run._levels)))
snapshots = [] snapshots = []
for stop in (0, 1, 2, 4, 8): for stop in (0, 1, 2, 4, 8):
@@ -78,10 +78,10 @@ class ThreeLevelMultipatchGPU(unittest.TestCase):
def run_case(self, ctx, cfg=CFG, parents=PARENTS, children=CHILDREN, def run_case(self, ctx, cfg=CFG, parents=PARENTS, children=CHILDREN,
mode='periodic', blocks=((4, 4), (4, 4)), objects=(), drive=(), switch=False, mode='periodic', blocks=((4, 4), (4, 4)), objects=(), drive=(), switch=False,
scalar_child=None): scalar_child=None):
child = children if scalar_child is None else scalar_child kwargs = dict(child_rectangles_in_parent=children) if scalar_child is None else dict(child=scalar_child)
with FixedPatchStepper(config(cfg, mode), context=ctx, coarse_boundary=mode, with FixedPatchStepper(config(cfg, mode), parents, context=ctx, coarse_boundary=mode,
refinements=[parents, child], level_blocks=[(4, 4), blocks[0], blocks[1]], base_block=(4, 4), fine_block=blocks[0], child_block=blocks[1],
level_objects=[(), (), objects], level_jets=[(), (), drive]) as run: child_objects=objects, child_jets=drive, **kwargs) as run:
frame = frames(run) frame = frames(run)
run.initialize(*(physical_state(level.shape[1:], *f) run.initialize(*(physical_state(level.shape[1:], *f)
for level, f in zip(run._levels, frame))) for level, f in zip(run._levels, frame)))
@@ -26,8 +26,9 @@ FP64_ATOL = 3e-6
def make(ctx): def make(ctx):
return FixedPatchStepper(CFG, context=ctx, refinements=[PARENTS, A], return FixedPatchStepper(CFG, PARENTS, context=ctx, base_block=(4, 4),
level_blocks=[(4, 4)] * 3) fine_block=(4, 4), child_block=(4, 4),
child_rectangles_in_parent=A)
def initial(run, uniform): def initial(run, uniform):
@@ -136,8 +136,8 @@ class ShiftedJetsGPU(unittest.TestCase):
Parallel(0.),amplitude=.013) Parallel(0.),amplitude=.013)
for rho in (.99,1.01): for rho in (.99,1.01):
with self.subTest(rho=rho),CudaContext(0) as ctx: with self.subTest(rho=rho),CudaContext(0) as ctx:
with FixedPatchStepper(cfg,context=ctx,refinements=[(8,6,24,18)],level_blocks=[(4,4)]*2, with FixedPatchStepper(cfg,(8,6,24,18),context=ctx,base_block=(4,4),fine_block=(4,4),
level_objects=[(),(obj,)],level_jets=[(),(whole,)],observe_raw_exchange=True) as run: fine_objects=(obj,),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._level_boundaries[1].tables['aperture'].all())
@@ -153,8 +153,7 @@ class ShiftedJetsGPU(unittest.TestCase):
obj=circle((16.,12.),2.15);j=drive((16.,12.),2.15) obj=circle((16.,12.),2.15);j=drive((16.,12.),2.15)
for rho in (.99,1.01): for rho in (.99,1.01):
with self.subTest(rho=rho),CudaContext(0) as ctx,ExitStack() as owned: with self.subTest(rho=rho),CudaContext(0) as ctx,ExitStack() as owned:
settings=dict(context=ctx,refinements=[(8,6,24,18)],level_blocks=[(4,4)]*2, settings=dict(bounds=(8,6,24,18),context=ctx,base_block=(4,4),fine_block=(4,4),fine_objects=(obj,),jets=(j,))
level_objects=[(),(obj,)],level_jets=[(),(j,)])
measured=owned.enter_context(FixedPatchStepper(cfg,observe_raw_exchange=True,**settings)) measured=owned.enter_context(FixedPatchStepper(cfg,observe_raw_exchange=True,**settings))
plain=owned.enter_context(FixedPatchStepper(cfg,**settings)) plain=owned.enter_context(FixedPatchStepper(cfg,**settings))
measured._initialize_raw(initial(measured,rho));plain._initialize_raw(initial(plain,rho)) measured._initialize_raw(initial(measured,rho));plain._initialize_raw(initial(plain,rho))
@@ -187,9 +186,8 @@ class ShiftedJetsGPU(unittest.TestCase):
parents=((4,4,20,28),(24,4,44,28));a=((16,16,32,40),(56,16,72,40));b=(a[0],(60,16,76,40)) parents=((4,4,20,28),(24,4,44,28));a=((16,16,32,40),(56,16,72,40));b=(a[0],(60,16,76,40))
obj=circle((12.,14.),.8);j=drive((12.,14.),.8) obj=circle((12.,14.),.8);j=drive((12.,14.),.8)
with CudaContext(0) as ctx,ExitStack() as owned: with CudaContext(0) as ctx,ExitStack() as owned:
old=owned.enter_context(FixedPatchStepper(cfg,context=ctx,refinements=[parents,a], old=owned.enter_context(FixedPatchStepper(cfg,parents,context=ctx,base_block=(4,4),fine_block=(4,4),
level_blocks=[(4,4)]*3,level_objects=[(),(),(obj,)],level_jets=[(),(),(j,)], child_block=(4,4),child_rectangles_in_parent=a,child_objects=(obj,),child_jets=(j,),observe_raw_exchange=True))
observe_raw_exchange=True))
old._initialize_raw(initial(old,.99));ledger=Ledger(self,old,2) old._initialize_raw(initial(old,.99));ledger=Ledger(self,old,2)
for speed in (.002,-.001):old.set_jet_amplitude(0,speed,level=2);ledger.step() for speed in (.002,-.001):old.set_jet_amplitude(0,speed,level=2);ledger.step()
before=old._download_raw();count=old.step_count before=old._download_raw();count=old.step_count
@@ -22,8 +22,9 @@ B=((20,16,36,40),)
def make(ctx): def make(ctx):
return FixedPatchStepper(CFG,context=ctx,refinements=[PARENTS,A], return FixedPatchStepper(CFG,PARENTS,context=ctx,base_block=(4,4),
level_blocks=[(4,4)]*3) fine_block=(4,4),child_block=(4,4),
child_rectangles_in_parent=A)
def initial(run): def initial(run):
+2 -2
View File
@@ -126,8 +126,8 @@ def run_case(a, h, case):
cfg, pair, body, host = rr.host_geometry(a.source, h, a.collision) cfg, pair, body, host = rr.host_geometry(a.source, h, a.collision)
context = stack.enter_context(CudaContext(0)) context = stack.enter_context(CudaContext(0))
backend = stack.enter_context(FixedPatchStepper( backend = stack.enter_context(FixedPatchStepper(
cfg, context=context, coarse_boundary="channel", refinements=[rr.PLACEMENTS["wide"]], cfg, rr.PLACEMENTS["wide"], context=context, coarse_boundary="channel",
level_blocks=[rr.BASE_BLOCK, rr.FINE_BLOCK], level_objects=[(), (body,)], base_block=rr.BASE_BLOCK, fine_block=rr.FINE_BLOCK, fine_objects=(body,),
observe_raw_exchange=True, observe_window_exchange=True)) observe_raw_exchange=True, observe_window_exchange=True))
from CelerisLab.lbm.channel import channel_masks from CelerisLab.lbm.channel import channel_masks
from CelerisLab.lbm.refinement import ACTIVE from CelerisLab.lbm.refinement import ACTIVE
+5 -4
View File
@@ -279,11 +279,12 @@ def run_case(c, ctx, output, provenance, continuation=False):
PeriodicBlockStepper(fine_cfg, layout, context=ctx, objects=[fobj]) as uniform, PeriodicBlockStepper(fine_cfg, layout, context=ctx, objects=[fobj]) as uniform,
FixedPatchStepper( FixedPatchStepper(
cfg, cfg,
c["bounds"],
context=ctx, context=ctx,
refinements=[c["bounds"]], base_block=(8, 8),
level_blocks=[(8, 8), (8, 8)], fine_block=(8, 8),
level_objects=[(), (obj,)], fine_objects=[obj],
level_jets=[(), tuple(jets)], jets=jets,
observe_raw_exchange=True, observe_raw_exchange=True,
) as amr, ) as amr,
): ):
+5 -4
View File
@@ -383,12 +383,13 @@ def run_case(a, h, periodic, case):
backend = stack.enter_context( backend = stack.enter_context(
FixedPatchStepper( FixedPatchStepper(
cfg, cfg,
PLACEMENTS[a.placement],
context=context, context=context,
coarse_boundary="channel", coarse_boundary="channel",
refinements=[PLACEMENTS[a.placement]], base_block=BASE_BLOCK,
level_blocks=[BASE_BLOCK, FINE_BLOCK], fine_block=FINE_BLOCK,
level_objects=[(), (body,)], fine_objects=(body,),
level_jets=[(), public_jets((80.5, 81.0), 20.0, a.jet_scenario)], jets=public_jets((80.5, 81.0), 20.0, a.jet_scenario),
observe_raw_exchange=True, observe_raw_exchange=True,
**({"observe_window_exchange": True} if a.observation_mode == "window" else {}), **({"observe_window_exchange": True} if a.observation_mode == "window" else {}),
) )
+2 -3
View File
@@ -173,9 +173,8 @@ def time_refined(a, h, sched):
cfg, pair, body, host = r.host_geometry(a.source, h, "MRT") cfg, pair, body, host = r.host_geometry(a.source, h, "MRT")
context = stack.enter_context(CudaContext(0)) context = stack.enter_context(CudaContext(0))
backend = stack.enter_context(FixedPatchStepper( backend = stack.enter_context(FixedPatchStepper(
cfg, context=context, coarse_boundary="channel", refinements=[r.BOUNDS], cfg, r.BOUNDS, context=context, coarse_boundary="channel", base_block=r.BASE_BLOCK,
level_blocks=[r.BASE_BLOCK, r.FINE_BLOCK], level_objects=[(), (body,)], fine_block=r.FINE_BLOCK, fine_objects=(body,), observe_raw_exchange=True))
observe_raw_exchange=True))
masks = [channel_masks(pair.coarse_shape)[0] & ~pair.covered, # as run_refined_dfg masks = [channel_masks(pair.coarse_shape)[0] & ~pair.covered, # as run_refined_dfg
(pair.fine_kind == ACTIVE) & ((host.flags.reshape(pair.fine_shape) & SOLID) == 0)] (pair.fine_kind == ACTIVE) & ((host.flags.reshape(pair.fine_shape) & SOLID) == 0)]
fresh = (r.initial_state(pair.coarse_shape, (0.0, 0.0), 1.0, root=True), fresh = (r.initial_state(pair.coarse_shape, (0.0, 0.0), 1.0, root=True),
+2 -2
View File
@@ -130,8 +130,8 @@ def main():
cfg, pair, body, host = r.host_geometry(a.source, h, "MRT") cfg, pair, body, host = r.host_geometry(a.source, h, "MRT")
context = stack.enter_context(CudaContext(0)) context = stack.enter_context(CudaContext(0))
backend = stack.enter_context(FixedPatchStepper( backend = stack.enter_context(FixedPatchStepper(
cfg, context=context, coarse_boundary="channel", refinements=[r.BOUNDS], cfg, r.BOUNDS, context=context, coarse_boundary="channel", base_block=r.BASE_BLOCK,
level_blocks=[r.BASE_BLOCK, r.FINE_BLOCK], level_objects=[(), (body,)], observe_raw_exchange=True, fine_block=r.FINE_BLOCK, fine_objects=(body,), observe_raw_exchange=True,
observe_window_exchange=True)) observe_window_exchange=True))
masks = [channel_masks(pair.coarse_shape)[0] & ~pair.covered, masks = [channel_masks(pair.coarse_shape)[0] & ~pair.covered,
(pair.fine_kind == ACTIVE) & ((host.flags.reshape(pair.fine_shape) & SOLID) == 0)] (pair.fine_kind == ACTIVE) & ((host.flags.reshape(pair.fine_shape) & SOLID) == 0)]