Compare commits

..
Author SHA1 Message Date
ClaudeandClaude Opus 5.5 106aa620b9 test: annulus mass ledger definition with synthetic seam/edge checks
Host-only stdlib definition of the per-hop identity
dM + Q_exact = C + H with Q_exact from pre-staging owner populations and
H = sum over every annulus cell/direction of (pulled - ideal owner input),
split into seam (internal) and edge parts. Synthetic tests: ideal pull,
collision source, interior-seam and edge perturbations, the double count
of staged flux + full H, a late owner snapshot hiding staging, and a
wrong pull direction. No solver change, no GPU.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-20 18:32:42 +00:00
40 changed files with 578 additions and 665 deletions
+36 -38
View File
@@ -112,9 +112,8 @@ from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
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,
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)],
level_blocks=[None, (16, 20)],
level_objects=[(), (body,)]) as run:
with CudaContext(0) as ctx, FixedPatchStepper(cfg, (8, 6, 24, 18), context=ctx,
fine_objects=[body], fine_block=(16, 20)) as run:
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.
run.initialize(np.broadcast_to(rest, (9, *run.patch.coarse_shape)),
@@ -125,8 +124,8 @@ with CudaContext(0) as ctx, FixedPatchStepper(cfg, context=ctx, refinements=[(8,
Scope and limits:
- The examples below use D2Q9, SRT, FP32, double buffer, unshifted, no LES or
forcing. Fixed 2:1 hierarchies use level lists of arbitrary depth; two- and
three-level examples below use that same constructor.
forcing. Fixed 2:1 hierarchies can use level lists of arbitrary depth; the
legacy two/three-level forms are shown below.
- `coarse_boundary="periodic"` (default) or `"slot_reservoir"` (bottom slot, top
reservoir; patch needs `y0 >= 4`, `y1 <= ny-4`). Periodic mode ignores the wall/slot config fields.
- Object centers/radii use the coarse **cell-edge** frame (cell `I` spans `[I, I+1)`),
@@ -134,51 +133,50 @@ Scope and limits:
- 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
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
four-cell ghost band. It gives results identical to dense storage and is
infrastructure for future block AMR, not a speedup.
- 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
their actual dependencies close. There is no fixed rectangle-gap rule.
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
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
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
inside the union's active region with
its boundary stencil clear of transfer cells; a jet may not select links of two
different bodies, even across regions. Wall-touching patches remain unsupported.
For nesting, add further `refinements` entries as shown below.
For nesting, use the level lists or legacy child forms below.
```python
cfg = LBMConfig(nx=32, ny=16, viscosity=.1)
rest = np.asarray(compute_feq_d2q9(1.0, 0.0, 0.0), np.float32)[:, None, None]
with CudaContext(0) as ctx, FixedPatchStepper(
cfg, context=ctx, refinements=[((4, 4, 12, 12), (16, 4, 24, 12))],
level_blocks=[(4, 4), (4, 4)]) as run:
cfg, ((4, 4, 12, 12), (16, 4, 24, 12)), context=ctx,
fine_block=(4, 4), base_block=(4, 4)) as run:
run.initialize(np.broadcast_to(rest, (9, *run.patch.coarse_shape)),
np.broadcast_to(rest, (9, *run.patch.fine_shape))) # (9, 32, 64)
run.step(32)
coarse_g, fine_g = run.download() # coarse: NaN on every covered cell;
# fine: NaN on every non-ACTIVE cell
```
- Three-level nested windows: pass two scalar rectangles in `refinements` —
L1 in L0 cells, L2 in **L1 window cells** (the parent's fine grid including its
4-cell ghost band) — plus all three `level_blocks` entries. L2 must sit inside
L1's ACTIVE cells, away from L1's ghost band and transfer footprint (checked
before 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 steps = 4 L2 steps. One L1 window only in this scalar
chain. Optional stationary circles and Jets live on L2 via `level_objects[2]`
and `level_jets[2]`; L1 bodies/Jets retain 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 dx=dt=0.25, so prescribed lattice
velocity is unchanged. Each whole body and its BFL stencil must be inside L2
ACTIVE cells, outside the transfer footprint. A Jet may select one physical
body across multiple blocks, never two bodies.
- 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
`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
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
steps = 4 L2 steps. One L1 patch only. Optional stationary circles and Jets
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
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
its BFL stencil must be inside L2 ACTIVE cells, outside the transfer footprint.
A Jet may select one physical body across multiple blocks, never two bodies.
```python
from CelerisLab.body.geometry.circle import CircleGeometry
@@ -190,9 +188,9 @@ Scope and limits:
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]
with CudaContext(0) as ctx, FixedPatchStepper(
cfg, context=ctx, refinements=[(4, 4, 12, 12), (8, 8, 16, 16)],
level_blocks=[(4, 4)] * 3, level_objects=[(), (), (body,)],
level_jets=[(), (), (jet,)]) as run:
cfg, (4, 4, 12, 12), context=ctx, fine_block=(4, 4), base_block=(4, 4),
child=(8, 8, 16, 16), child_block=(4, 4),
child_objects=[body], child_jets=[jet]) as run:
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.step(4)
@@ -200,16 +198,16 @@ Scope and limits:
run.step(4)
coarse_g, l1_g, l2_g = run.download() # solids/covered/nonleaves are NaN
```
- Level lists (the only constructor): `refinements[l]` describes level `l+1`
in allocated level-`l` integer cells; use a rectangle list for a union, or a
scalar rectangle for a single window (a scalar child requires a scalar
parent). `level_blocks`, `level_objects` and `level_jets` each have one entry
per level, including L0. Three or more levels require all-blocked storage and
bodies/Jets only on the finest level. Every adjacent pair must preserve
inherited ACTIVE coverage, valid donors/receivers, surviving pulls and the
upstream transfer footprint. Level `l` takes `2**l` steps per base interval.
This describes implemented topology and scheduling, not accuracy or
convergence qualification.
- Level lists: `refinements[l]` describes level `l+1` in allocated level-`l`
integer cells; use a rectangle list for a union, or a scalar rectangle for a
single window (a scalar child requires a scalar parent). `level_blocks`,
`level_objects` and `level_jets` each have one entry per level, including L0.
These lists are mutually exclusive with the legacy patch keywords. Three or
more levels require all-blocked storage and bodies/Jets only on the finest
level. Every adjacent pair must preserve inherited ACTIVE coverage, valid
donors/receivers, surviving pulls and the upstream transfer footprint.
Level `l` takes `2**l` steps per base interval. This describes implemented
topology and scheduling, not accuracy or convergence qualification.
- Blocked body levels can enable `observe_raw_exchange=True` and call
`read_raw_exchange(body_id, level=...)` after a completed base interval.
It reduces unique cut-link observations into wall/aperture/total force,
+1 -1
View File
@@ -1,6 +1,6 @@
# 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`:
+15 -15
View File
@@ -1,8 +1,8 @@
# Fixed three-level rectangle unions
This page describes three-level `FixedPatchStepper` setups built from the
`refinements` level list, with rectangle unions on either adjacent pair. The
same constructor supports arbitrary depth; three or more levels require
This page describes the legacy three-level entry point of `FixedPatchStepper`,
with rectangle unions on either adjacent pair. The `refinements` level-list
entry point also supports arbitrary depth; three or more levels require
all-blocked storage and bodies/Jets only on the finest level. This is not
general AMR: there is no per-region timestep, moving geometry, automatic
dynamic refinement or L1 body support in the three-level form.
@@ -12,15 +12,15 @@ No CUDA kernel changes are involved.
## Explicit coordinate frames
For a single scalar parent patch, pass a scalar rectangle as `refinements[1]`:
`(x0,y0,x1,y1)` in half-open integer indices of the allocated L1 window,
including its ghost band.
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
allocated L1 window, including its ghost band.
A rectangle list as `refinements[1]` 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 region-local coordinates. The physical root
cell edge corresponding to an L1 index `(x,y)` is
`parent.fine_origin + (x,y)/2`:
The new keyword `child_rectangles_in_parent=[...]` is mutually exclusive with
`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
region-local coordinates. The physical root cell edge corresponding to an L1
index `(x,y)` is `parent.fine_origin + (x,y)/2`:
- A scalar L1 parent has a shifted window origin.
- A set L1 parent has origin `(0,0)` and public shape `(2*ny,2*nx)`.
@@ -47,9 +47,9 @@ tile edges.
```python
run = FixedPatchStepper(
cfg, context=ctx,
refinements=[[(4, 4, 20, 28), (24, 4, 44, 28)], [(16, 16, 32, 40), (56, 16, 80, 40)]],
level_blocks=[(4, 4)] * 3,
cfg, [(4, 4, 20, 28), (24, 4, 44, 28)], context=ctx,
base_block=(4, 4), fine_block=(4, 4), child_block=(4, 4),
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
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
positions, not region IDs or caller `obj_id`. Complete boundary stencils must
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
# nested with-statement does that).
with CudaContext(0) as context, FixedPatchStepper(
cfg, context=context, coarse_boundary=args.mode, refinements=[BOUNDS],
level_blocks=[None, fine_block], level_objects=[(), (body,)]) as run:
cfg, BOUNDS, context=context, coarse_boundary=args.mode,
fine_objects=[body], fine_block=fine_block) as run:
p = run.patch
# Rest equilibrium is already postcollision G, so no initial collision is
# 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)
start, run = time.monotonic(), None
try:
with CudaContext(0) as context, FixedPatchStepper(cfg, context=context,
coarse_boundary='slot_reservoir', refinements=[BOUNDS]) as run:
with CudaContext(0) as context, FixedPatchStepper(cfg, BOUNDS, context=context,
coarse_boundary='slot_reservoir') as run:
p = run.patch
masks = [np.ones(p.coarse_shape, bool), p.fine_kind == ACTIVE]
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,
lattice_model="D2Q9", collision="SRT", streaming="double_buffer",
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
coarse = shear_equilibrium(p.coarse_shape, (0, 0), 1.0, cfg.ny)
fine = shear_equilibrium(p.fine_shape, p.fine_origin, .5, cfg.ny)
+4 -14
View File
@@ -204,21 +204,11 @@ class PeriodicBlockStepper:
raise
def step(self, n=1):
"""Advance n steps: halo refresh, then core pull/collide, per step.
All n steps are enqueued in order on the one stream and the host waits
once before returning (same launches and bits as n step(1) calls).
Asynchronous device errors surface at that single wait, so they are not
attributed to the step that caused them; call step(1) repeatedly to
localize one. On any failure the instance is failed and step_count keeps
its entry value.
"""
"""Advance n steps: halo refresh, then core pull/collide, per step."""
if isinstance(n, bool) or not isinstance(n, Integral) or n < 0:
raise ValueError("n must be a nonnegative integer")
with self._lock, self.ctx.activate():
self._require("ready")
if n == 0:
return
kw = dict(block=self._block, stream=self._stream)
try:
self._status = "advancing"
@@ -234,9 +224,9 @@ class PeriodicBlockStepper:
self._step(self._current, self._spare, np.uint32(self.layout.halo),
self._tables.get("flags", np.uintp(0)), np.uintp(0), # no kind
grid=self._grid, **kw)
self._current, self._spare = self._spare, self._current # host handles only
self._stream.synchronize()
self._step_count += n
self._stream.synchronize()
self._current, self._spare = self._spare, self._current
self._step_count += 1
self._status = "ready"
except BaseException:
self._status = "failed"
+144 -85
View File
@@ -65,8 +65,52 @@ def _level_name(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):
"""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
if depth < 2:
raise ValueError("at least one refinement level is required")
@@ -75,13 +119,13 @@ def _check_levels(refinements, blocks, objects, jets):
raise ValueError(f"{name} needs one entry per level ({depth})")
for level, block in enumerate(blocks):
if block is not None and len(block) != 2:
raise ValueError(f"level_blocks[{level}] must be (by, bx)")
raise ValueError(f"{_level_name(level)}_block must be (by, bx)")
if objects[0] or jets[0]:
raise ValueError("base-level bodies/Jets are not supported")
if any(blocks[level - 1] is not None and blocks[level] is None for level in range(1, depth)):
raise ValueError("a blocked level requires blocked refined levels below it")
if not is_rectangle(refinements[0]) and (blocks[0] is None or blocks[1] is None):
raise ValueError("several L1 patches require level_blocks[0] and level_blocks[1]")
raise ValueError("several patches require fine_block and base_block")
if depth > 2:
if any(block is None for block in blocks):
raise ValueError("three or more levels must all be blocked")
@@ -90,14 +134,7 @@ def _check_levels(refinements, blocks, objects, jets):
class FixedPatchStepper:
"""D2Q9 SRT/MRT FP32 fixed 2:1 refinement hierarchy, built from level lists.
``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).
"""D2Q9/SRT/FP32, one unforced internal 2:1 patch in an x-periodic box.
``coarse_boundary`` is ``periodic`` (default) or ``slot_reservoir``:
a bottom on-site slot and top fixed-density reservoir from ``cfg``.
@@ -108,55 +145,68 @@ class FixedPatchStepper:
actual channel closure footprint. Channel remapping is not supported.
Periodic mode does not apply the wall/slot boundary fields in ``cfg``.
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
in root lattice units (dx=dt=1). Independent local copies are projected at
fine cell centers; their complete stencil must be active and disjoint from
the actual two-hop transfer footprint. These use the edge/cell-center frame
(coarse centers I+.5), whereas ordinary Simulation objects use the
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.
``initialize(coarse_g, fine_g)`` accepts physical Q-major (9,ny,nx) arrays
already in postcollision phase. It does not collide or build channel flags.
Every ``step()`` publishes only after both levels complete on one stream.
On a device/launch failure the instance becomes unusable: close and recreate
it, rather than pretending the partially advanced fine buffer rolled back.
``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
without changing populations. read_raw_exchange reports the last complete
base interval or last substep in base-lattice units per unit depth. This
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
per-link FP64 ledger accumulated by the same kernel between
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,
including ghost bands and holes: root edge = parent fine_origin + index / 2
(zero origin for a parent set; a scalar parent keeps its shifted window
origin). A scalar child requires a scalar parent. Rectangle unions impose no
fixed gap rule: adjacent, overlapping and concave unions are accepted when
their actual transfer dependencies close (FixedPatchSet2D /
nest_patch_set). All regions share one buffer pair per level and the fixed
schedule of 2**l hops per base step at level l. Bodies/Jets are supported only on the
finest level of a nested hierarchy (L2 dx=dt=.25 in root units, unchanged
lattice velocity); every whole circle and its BFL stencil must fit strictly
away from that level's transfer footprint.
Several patches: pass ``bounds`` as a sequence of (x0, y0, x1, y1)
rectangles (FixedPatchSet2D, dependency-validated union). This requires
fine_block and base_block; fine_objects/jets are projected once onto the
level-wide fine grid (each circle inside the union, body ids = object
positions, jet and writer checks level-wide), and
uses one level-wide fine grid (2*ny, 2*nx): public fine arrays have that
shape, NaN outside the patches' ACTIVE cells. The step sequence is the
blocked one above, with unique union transfers in one launch each.
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
reconstruction on each adjacent pair. Default "constant" is unchanged.
@@ -164,21 +214,43 @@ class FixedPatchStepper:
bodies/Jets or LES; this backend supplies no forcing. Geometry is generic,
but union/multilevel GPU and matched-physics qualification remain pending.
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.
Level lists (any depth): ``refinements[l]`` gives level l+1 in level-l
integer cells (scalar window or rectangle list, see lbm.hierarchy);
``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",
level_blocks=None, level_objects=None, level_jets=None,
def __init__(self, cfg, bounds=None, *, context, coarse_boundary="periodic", fine_objects=(),
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"):
if HALO != 4:
raise RuntimeError("refinement.cu requires a four-cell fine halo")
if coarse_boundary not in ("periodic", "slot_reservoir", "channel"):
raise ValueError("coarse_boundary must be periodic, slot_reservoir or channel")
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)
if refinements is None:
if bounds is None:
raise ValueError("bounds or refinements is required")
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_jets = [tuple(j) for j in level_jets]
level_blocks = [None if b is None else tuple(b) for b in level_blocks]
@@ -468,13 +540,13 @@ class FixedPatchStepper:
if self._ddf_shifting and any(layout is None for layout in self._layouts):
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.
Shifted input should be physical FP64 to retain small perturbations.
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():
self._require("new", "ready")
self._check_storage()
@@ -598,9 +670,10 @@ class FixedPatchStepper:
old, new, source, base_count=count, level_steps=(count, 2*count, 4*count),
phase="postcollision_g", **options)
candidate = type(self)(
cfg, context=self.ctx, refinements=[pairs[0].bounds, new.child.bounds],
level_blocks=[base_block, fine_block, child_block],
level_objects=[(), (), child_objects], level_jets=[(), (), jets],
cfg, pairs[0].bounds, context=self.ctx, base_block=base_block,
fine_block=fine_block, child_block=child_block,
child_rectangles_in_parent=new.child.bounds,
child_objects=child_objects, child_jets=jets,
observe_raw_exchange=self._observe_raw_exchange,
observe_window_exchange=self._observe_window, explosion=self._explosion)
try:
@@ -749,41 +822,27 @@ class FixedPatchStepper:
grid=((pairs + self._block[0] - 1) // self._block[0], 1), **kw)
def step(self, n=1):
"""Advance G→G by n base intervals: level l takes 2**l hops per interval.
All n intervals are enqueued in order on the one stream, and the host waits
once, before returning; step(0) launches and waits for nothing. On return
the device work is complete, step_count and an open window cover all n
intervals, and the raw-exchange ledger holds the final interval only.
Asynchronous device errors surface at that single wait, so they are not
attributed to the interval that caused them; call step(1) repeatedly to
localize one (same launches, same bits). On any failure the instance is
failed, the exchange ledger and window are invalid, and step_count keeps
its value from entry to the failing call.
"""
"""Advance G→G by n base intervals: level l takes 2**l hops per interval."""
if isinstance(n, bool) or not isinstance(n, Integral) or n < 0:
raise ValueError("n must be a nonnegative integer")
with self._lock, self.ctx.activate():
self._require("ready")
if n == 0:
return
root = self._levels[0]
kw = dict(block=self._block, stream=self._stream)
try:
self._status = "advancing"
self._exchange_valid = False
for _ in range(n):
if self._observe_raw_exchange: # stream-ordered: holds the last interval only
if self._observe_raw_exchange:
self._exchange_valid = False
for b in self._level_boundaries:
if b is not None:
cuda.memset_d8_async(b.interval_integral, 0,
len(b.tables["cell"]) * 32, self._stream)
self._hop(0, root.current, root.spare, **kw)
root.current, root.spare = root.spare, root.current # host handles only
self._stream.synchronize()
self._step_count += n
self._exchange_valid = self._observe_raw_exchange
self._stream.synchronize()
root.current, root.spare = root.spare, root.current
self._step_count += 1
self._exchange_valid = self._observe_raw_exchange
self._status = "ready"
except BaseException:
self._window_valid = False
+101
View File
@@ -0,0 +1,101 @@
"""Synthetic checks of the annulus mass ledger. Host-only, stdlib; not solver evidence."""
import random
import sys
from pathlib import Path
import unittest
sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "validation"))
import annulus_ledger as al # noqa: E402
NX = NY = 24
# Annulus 4 <= r <= 8 around (12, 12); dyadic values keep every sum exact.
A = {(x, y) for x in range(NX) for y in range(NY) if 16 <= (x - 12) ** 2 + (y - 12) ** 2 <= 64}
DELTA = 0.375
def owner(seed=1):
rng = random.Random(seed)
return [[[rng.randint(1, 1024) / 1024 for _ in range(NX)] for _ in range(NY)] for _ in range(9)]
def collide(P, source=None):
"""Mass-conserving stand-in: rotate populations within each cell (+ optional source)."""
D = [[[P[(q + 1) % 9][y][x] for x in range(NX)] for y in range(NY)] for q in range(9)]
if source:
(x, y), s = source
D[0][y][x] += s
return D
def link(interior):
"""(x, y, q) with x in A and x - c_q inside (seam) or outside (edge) A."""
for (x, y) in sorted(A):
for q in range(1, 9):
src = ((x - al.CX[q]) % NX, (y - al.CY[q]) % NY)
if (src in A) == interior:
return x, y, q
raise AssertionError
class Ledger(unittest.TestCase):
def setUp(self):
self.O = owner()
self.P = al.pull(self.O, NX, NY)
def run_ledger(self, P, D):
return al.ledger(self.O, P, D, A, NX, NY)
def test_ideal_pull_has_no_h_and_closes(self):
r = self.run_ledger(self.P, collide(self.P))
self.assertEqual((r["H_seam"], r["H_edge"], r["C"], r["residual"]), (0, 0, 0, 0))
self.assertEqual(r["dM"], -r["Q_exact"])
self.assertNotEqual(r["Q_exact"], 0)
def test_collision_source_is_c(self):
x, y = sorted(A)[3]
r = self.run_ledger(self.P, collide(self.P, ((x, y), DELTA)))
self.assertEqual((r["C"], r["H"], r["residual"]), (DELTA, 0, 0))
def test_interior_seam_perturbation(self):
x, y, q = link(interior=True)
self.P[q][y][x] += DELTA
r = self.run_ledger(self.P, collide(self.P))
self.assertEqual((r["H_seam"], r["H_edge"], r["residual"]), (DELTA, 0, 0))
# Invisible to any flux through the edge of A; only H explains the mass change.
self.assertEqual(al.staged_flux(self.O, self.P, A, NX, NY), r["Q_exact"])
self.assertEqual(r["dM"] + r["Q_exact"], DELTA)
def test_edge_perturbation_and_staged_flux_double_count(self):
x, y, q = link(interior=False)
self.P[q][y][x] += DELTA
r = self.run_ledger(self.P, collide(self.P))
self.assertEqual((r["H_seam"], r["H_edge"], r["residual"]), (0, DELTA, 0))
staged = al.staged_flux(self.O, self.P, A, NX, NY)
self.assertEqual(staged, r["Q_exact"] - DELTA)
# Staged flux + full H counts the staging change twice.
self.assertEqual(r["dM"] + staged - r["C"] - r["H"], -DELTA)
# Staged flux alone closes only without the edge part of H.
self.assertEqual(r["dM"] + staged - r["C"] - r["H_seam"], 0)
def test_owner_snapshot_after_staging_hides_the_change(self):
# Snapshot taken after staging wrote the same slot into the owner image.
x, y, q = link(interior=False)
sx, sy = (x - al.CX[q]) % NX, (y - al.CY[q]) % NY
self.P[q][y][x] += DELTA
late = [[row[:] for row in plane] for plane in self.O]
late[q][sy][sx] += DELTA
r = al.ledger(late, self.P, collide(self.P), A, NX, NY)
self.assertEqual((r["H"], r["residual"]), (0, 0))
def test_wrong_pull_direction_is_flagged(self):
# Indexing error (push instead of pull): H is nonzero with no physical cause.
bad = [[[self.O[q][(y + al.CY[q]) % NY][(x + al.CX[q]) % NX] for x in range(NX)]
for y in range(NY)] for q in range(9)]
r = self.run_ledger(bad, collide(bad))
self.assertNotEqual(r["H"], 0)
self.assertEqual(r["residual"], 0)
if __name__ == "__main__":
unittest.main()
-102
View File
@@ -1,102 +0,0 @@
"""Host contract of PeriodicBlockStepper.step(n): one stream wait per call (mocked driver)."""
import ast
from contextlib import nullcontext
from pathlib import Path
import types
from unittest.mock import Mock
import unittest
import numpy as np
import CelerisLab.lbm.block_stepper as b
from CelerisLab.config import LBMConfig
from CelerisLab.lbm.blocks import BlockLayout2D
def host_module():
"""block_stepper with pycuda/compiler replaced by mocks; no GPU context."""
path = Path(b.__file__)
tree = ast.parse(path.read_text())
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)
or isinstance(node, ast.ImportFrom) and node.module == "cuda")]
module = types.ModuleType("CelerisLab.lbm._test_host_block_stepper")
module.__package__ = "CelerisLab.lbm"
module.__dict__.update(cuda=Mock(), compiler=Mock())
exec(compile(tree, str(path), "exec"), module.__dict__)
return module
def backend():
host, log = host_module(), []
host.cuda.mem_get_info.return_value = (1 << 40, 1 << 40)
host.cuda.mem_alloc.side_effect = lambda size: Mock(size=size) # distinct buffers
module = Mock()
module.get_global.return_value = (Mock(), 20)
module.get_function.side_effect = lambda name: Mock(
side_effect=lambda *a, **k: log.append((name, a[:2])))
host.compiler.load_module.return_value = module
cfg = LBMConfig(nx=16, ny=16, viscosity=.1)
run = host.PeriodicBlockStepper(cfg, BlockLayout2D.build(16, 16, 8, 8),
context=Mock(activate=nullcontext, sm_arch="sm_70"))
run._stream.synchronize.side_effect = lambda: log.append(("sync", ()))
run.initialize(np.full((9, 16, 16), 1 / 9))
log.clear()
return run, log
class BlockStepSync(unittest.TestCase):
def test_step_n_waits_once_and_matches_n_single_steps(self):
many, log_many = backend()
many.step(3)
single, log_single = backend()
for _ in range(3):
single.step(1)
names = [name for name, _ in log_many]
self.assertEqual(names.count("sync"), 1)
self.assertEqual(names[-1], "sync")
self.assertEqual(names.count("BlockStep"), 3)
# Same launch order and buffer roles as three step(1) calls.
strip = lambda log: [(n, tuple(id(x) for x in a)) for n, a in log if n != "sync"]
roles = lambda run, log: [(n, tuple({id(run._current): "c", id(run._spare): "s"}.get(x, "-")
for x in ids)) for n, ids in strip(log)]
self.assertEqual([n for n, _ in strip(log_many)], [n for n, _ in strip(log_single)])
self.assertEqual(roles(many, log_many), roles(single, log_single))
self.assertIsNot(many._current, many._spare)
self.assertEqual((many.step_count, single.step_count), (3, 3))
def test_step_zero_launches_and_waits_for_nothing(self):
run, log = backend()
run.step(0)
self.assertEqual((log, run.step_count), ([], 0))
def test_failure_keeps_entry_step_count(self):
run, log = backend()
run.step(1)
run._stream.synchronize.side_effect = RuntimeError("device fault")
with self.assertRaises(RuntimeError):
run.step(2)
self.assertEqual((run.step_count, run._status), (1, "failed"))
def test_launch_failure_mid_batch_fails_closed_without_waiting(self):
run, log = backend()
calls = []
step = run._step
def failing(*args, **kwargs):
calls.append(1)
if len(calls) == 2:
raise RuntimeError("launch failure")
return step(*args, **kwargs)
run._step = failing
with self.assertRaises(RuntimeError):
run.step(3)
self.assertNotIn("sync", [name for name, _ in log])
self.assertEqual((run.step_count, run._status), (0, "failed"))
with self.assertRaises(RuntimeError):
run.download()
if __name__ == "__main__":
unittest.main()
+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)])
context = Mock(activate=nullcontext, sm_arch="sm_70")
cfg = LBMConfig(nx=24, ny=20, viscosity=.1, velocity=.025)
backend = host.FixedPatchStepper(cfg, context=context, coarse_boundary="channel",
refinements=[(6, 6, 18, 14)], level_blocks=[(4, 6), (4, 4)])
backend = host.FixedPatchStepper(cfg, (6, 6, 18, 14), context=context,
coarse_boundary="channel", base_block=(4, 6), fine_block=(4, 4))
try:
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]
+2 -3
View File
@@ -157,9 +157,8 @@ class RawExchange(unittest.TestCase):
host = host_module()
ctx = Mock()
with self.assertRaisesRegex(ValueError, 'blocked body'):
host.FixedPatchStepper(LBMConfig(nx=32, ny=32, viscosity=.1), context=ctx,
refinements=[(4, 4, 28, 28)],
level_objects=[(), (body(radius=2.15),)], observe_raw_exchange=True)
host.FixedPatchStepper(LBMConfig(nx=32, ny=32, viscosity=.1), (4, 4, 28, 28),
context=ctx, fine_objects=[body(radius=2.15)], observe_raw_exchange=True)
ctx.activate.assert_not_called()
def test_all_depth_weights(self):
@@ -173,10 +173,8 @@ def bodies_and_jets(speed=.002):
return bodies,jets
def levels(parents=PARENTS,children=CHILDREN,blocks=((4,4),)*3,objects=((),(),()),jets=((),(),())):
"""Level-list constructor keywords for the L0 set, L1 set and L2 set used here."""
return dict(refinements=[parents,children],level_blocks=list(blocks),
level_objects=list(objects),level_jets=list(jets))
KW=dict(base_block=(4,4),fine_block=(4,4),child_block=(4,4),
child_rectangles_in_parent=CHILDREN)
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
objects,jets=bodies_and_jets() if bodies else ((),())
try:
run=module.FixedPatchStepper(CFG,context=context,
**levels(children=children,objects=((),(),objects),jets=((),(),jets)))
run=module.FixedPatchStepper(CFG,PARENTS,context=context,
child_objects=objects,child_jets=jets,
**dict(KW,child_rectangles_in_parent=children))
except RuntimeError:
for allocation in allocations:
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):
def test_invalid_inputs_fail_before_context_or_allocations(self):
context=Mock()
b=(4,4)
cases=[dict(children=(8,8,16,16)),dict(children=[]),
dict(children=[(16,16,32)]),
dict(children=[(16.5,16,32,40)]),
dict(children=[(16,16,32,24),(16,27,32,55)]),
dict(children=[(40,16,48,40)]),
dict(children=[(8,16,32,40)]),
dict(blocks=(b,b,None)),dict(blocks=(b,None,b)),dict(blocks=(None,b,b)),
dict(blocks=(b,b,(5,4))),dict(blocks=(b,(5,4),b)),
dict(blocks=(b,b,(4,))),dict(children=(16,16,32,40))]
cases=[dict(child=(8,8,16,16)),dict(child_rectangles_in_parent=[]),
dict(child_rectangles_in_parent=[(16,16,32)]),
dict(child_rectangles_in_parent=[(16.5,16,32,40)]),
dict(child_rectangles_in_parent=[(16,16,32,24),(16,27,32,55)]),
dict(child_rectangles_in_parent=[(40,16,48,40)]),
dict(child_rectangles_in_parent=[(8,16,32,40)]),
dict(child_block=None),dict(fine_block=None),dict(base_block=None),
dict(child_block=(5,4)),dict(fine_block=(5,4)),
dict(child_block=(4,)),dict(child_rectangles_in_parent=(16,16,32,40))]
for changes in cases:
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'):
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()
with self.assertRaisesRegex(ValueError,'only on the finest level'):
host.FixedPatchStepper(CFG,context=context,**levels(objects=((),objects,())))
with self.assertRaisesRegex(ValueError,'only on L2'):
host.FixedPatchStepper(CFG,PARENTS,context=context,fine_objects=objects,**KW)
from CelerisLab.body.jets import Jet,Tube,Parallel
spanning=Jet(Tube((10.5,14.),0.,25.,.3),Parallel(0.),amplitude=0.)
with self.assertRaisesRegex(ValueError,'one jet'):
host.FixedPatchStepper(CFG,context=context,
**levels(objects=((),(),objects),jets=((),(),[spanning])))
host.FixedPatchStepper(CFG,PARENTS,context=context,child_objects=objects,
child_jets=[spanning],**KW)
context.activate.assert_not_called()
host.cuda.mem_alloc.assert_not_called()
@@ -315,23 +313,21 @@ class ConstructorAndSchedule(unittest.TestCase):
class RemovedLegacyAPI(unittest.TestCase):
"""The pre-level-list constructor/initialize forms are gone and fail loudly."""
class InitialLevels(unittest.TestCase):
"""initialize() entry adapter: positional per-level arrays or the legacy names, never both."""
def test_removed_constructor_forms_raise_type_error(self):
import inspect
context = Mock()
with self.assertRaises(TypeError): # positional bounds
host.FixedPatchStepper(CFG, PARENTS, context=context)
removed = dict(bounds=PARENTS, fine_block=(4, 4), base_block=(4, 4), fine_objects=(),
jets=(), child=(8, 8, 16, 16), child_block=(4, 4), child_objects=(),
child_jets=(), child_rectangles_in_parent=CHILDREN)
for name, value in removed.items():
with self.subTest(keyword=name), self.assertRaises(TypeError):
host.FixedPatchStepper(CFG, context=context, **levels(), **{name: value})
context.activate.assert_not_called()
self.assertEqual(list(inspect.signature(host.FixedPatchStepper.initialize).parameters),
["self", "level_g"]) # no coarse_g/fine_g/child_g
def test_forms(self):
a, b, c = object(), object(), object()
f = host._initial_levels
self.assertEqual(f((a, b), None, None, None), (a, b))
self.assertEqual(f((a, b, c), None, None, None), (a, b, c))
self.assertEqual(f((a, b), None, None, c), (a, b, c)) # old initialize(c, f, child_g=g)
self.assertEqual(f((), a, b, None), (a, b))
self.assertEqual(f((), a, b, c), (a, b, c))
for args in (((a,), None, None, c), ((a, b, c), None, None, c), ((a,), a, b, None),
((), a, None, None), ((), None, b, c)):
with self.subTest(args=args), self.assertRaises(TypeError):
f(*args)
if __name__=='__main__':
+6 -5
View File
@@ -54,13 +54,14 @@ class ShiftedJets(unittest.TestCase):
(Jet(Tube((10.,14.),0.,26.,.3),Parallel(0.)),))]
for objects,jets in cases:
with self.subTest(objects=len(objects),jets=len(jets)),self.assertRaises(ValueError):
module.FixedPatchStepper(CFG,context=None,refinements=[backend.PARENTS,A],
level_blocks=[(4,4)]*3,level_objects=[(),(),objects],level_jets=[(),(),jets])
module.FixedPatchStepper(CFG,backend.PARENTS,context=None,
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))):
with self.assertRaises(ValueError):
module.FixedPatchStepper(CFG,context=None,refinements=[backend.PARENTS,A],
level_blocks=list(blocks),level_objects=[(),(),(body(),)],
level_jets=[(),(),(jet(),)])
module.FixedPatchStepper(CFG,backend.PARENTS,context=None,
base_block=blocks[0],fine_block=blocks[1],child_block=blocks[2],
child_rectangles_in_parent=A,child_objects=(body(),),child_jets=(jet(),))
module.cuda.mem_alloc.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):
module=backend.host_module()
cfg=replace(backend.CFG,ddf_shifting=True)
for blocks in ([None,None],[(4,4),None],[None,(4,4)]):
with self.assertRaisesRegex(ValueError,'all-blocked|blocked refined levels below'):
module.FixedPatchStepper(cfg,context=None,refinements=[(4,4,20,28)],
level_blocks=blocks)
for kwargs in ({},dict(base_block=(4,4)),dict(fine_block=(4,4))):
with self.assertRaisesRegex(ValueError,'all-blocked|base_block requires fine_block'):
module.FixedPatchStepper(cfg,(4,4,20,28),context=None,**kwargs)
module.compiler.compile_kernel.assert_not_called()
def test_uniform_block_codec_and_raw_snapshot(self):
-170
View File
@@ -1,170 +0,0 @@
"""Host contract of FixedPatchStepper.step(n): one stream wait per call (mocked driver).
All n base intervals are enqueued in order on the backend stream and the host
waits once before step returns. Launch order and buffer arguments must equal n
step(1) calls; completion metadata, raw-exchange and window ledgers, and failure
invalidation are published only after that wait.
"""
from contextlib import nullcontext
from unittest.mock import Mock
import unittest
import numpy as np
from CelerisLab.config import LBMConfig
from tests.unit.test_raw_exchange import body
from tests.unit.test_refinement_three_level_multipatch import host_module
BODY_LEVEL, HOPS = 3, 8 # body on L3: 2**3 hops per base interval
def backend(window=False):
"""Four-level host backend; every launch, memset and stream wait goes to one log."""
host = host_module()
log, allocations = [], []
def allocate(size):
allocations.append(Mock(size=size))
return allocations[-1]
module = Mock()
module.get_global.return_value = (Mock(), 20)
module.get_function.side_effect = lambda name: Mock(
side_effect=lambda *a, **k: log.append((name, a, k)))
host.compiler.load_module.return_value = module
host.cuda.mem_alloc.side_effect = allocate
host.cuda.memcpy_dtoh.side_effect = lambda array, pointer: array.fill(0)
host.cuda.memset_d8_async.side_effect = lambda *a: log.append(("memset", a, {}))
run = host.FixedPatchStepper(
LBMConfig(nx=48, ny=48, viscosity=.02), context=Mock(activate=nullcontext, sm_arch="sm_70"),
refinements=[[(6, 6, 26, 26)], [(20, 20, 44, 44)], [(56, 56, 72, 72)]],
level_blocks=[(4, 4)] * 4, level_objects=[(), (), (), (body(),)],
observe_raw_exchange=True, observe_window_exchange=window)
run._stream.synchronize.side_effect = lambda: log.append(("sync", (), {}))
run.initialize(*[np.full(level.shape, 1 / 9, np.float32) for level in run._levels])
log.clear()
return run, log, allocations
def signature(log, allocations):
"""Event names with device-pointer arguments replaced by allocation index."""
index = {id(a): i for i, a in enumerate(allocations)}
return [(name, tuple(index.get(id(a), "-") for a in args)) for name, args, _ in log]
def names(log):
return [name for name, _, _ in log]
class StepSync(unittest.TestCase):
def test_step1_waits_once_after_all_launches(self):
run, log, _ = backend()
run.step(1)
self.assertEqual(names(log).count("sync"), 1)
self.assertEqual(names(log)[-1], "sync")
self.assertEqual(names(log)[0], "memset")
self.assertEqual(names(log).count("BlockStep"), 1 + 2 + 4 + 8)
self.assertEqual(run.step_count, 1)
run.close()
def test_step_n_equals_n_step1_launches_with_one_wait(self):
single, single_log, single_alloc = backend()
for _ in range(3):
single.step(1)
batched, batched_log, batched_alloc = backend()
batched.step(3)
expected = [e for e in signature(single_log, single_alloc) if e[0] != "sync"]
got = signature(batched_log, batched_alloc)
self.assertEqual(got[-1][0], "sync")
self.assertEqual(names(batched_log).count("sync"), 1)
self.assertEqual(got[:-1], expected) # same order, same ping-pong buffers
self.assertEqual(names(batched_log).count("memset"), 3)
self.assertEqual(batched.step_count, single.step_count)
index = [{id(a): i for i, a in enumerate(alloc)} for alloc in (single_alloc, batched_alloc)]
self.assertEqual(index[1][id(batched._levels[0].current)],
index[0][id(single._levels[0].current)])
single.close()
batched.close()
def test_step0_launches_and_waits_for_nothing(self):
run, log, _ = backend()
run.step(2)
log.clear()
run.step(0)
self.assertEqual(log, [])
self.assertEqual(run.step_count, 2)
value = run.read_raw_exchange(0, level=BODY_LEVEL) # step(0) keeps the ledger valid
self.assertEqual((value["time_start"], value["time_end"]), (1, 2))
run.close()
def test_completion_metadata_covers_all_intervals(self):
run, _, _ = backend()
with self.assertRaisesRegex(RuntimeError, "completed"):
run.read_raw_exchange(0, level=BODY_LEVEL)
run.step(3)
value = run.read_raw_exchange(0, level=BODY_LEVEL)
self.assertEqual((value["time_start"], value["time_end"], value["sample_count"]),
(2, 3, HOPS)) # step(n) retains only the final base interval
run.step(2)
self.assertEqual(run.read_raw_exchange(0, level=BODY_LEVEL)["time_end"], 5)
run.close()
def test_window_accounting_across_batched_calls(self):
windows = []
for calls in ((1, 3), (4,), (1, 1, 1, 1)):
run, log, _ = backend(window=True)
run.step(1)
run.start_exchange_window()
log.clear()
for n in calls:
run.step(n)
self.assertEqual(names(log).count("BlockAccumulateExchange"), 2 * 4 * HOPS)
self.assertEqual(names(log).count("sync"), len(calls))
out = run.read_exchange_window(0, level=BODY_LEVEL)
windows.append((out["time_start"], out["time_end"], out["sample_count"]))
run.close()
self.assertEqual(windows, [(1., 5, 4 * HOPS)] * 3)
def test_failed_wait_publishes_nothing(self):
run, _, _ = backend(window=True)
run.step(2)
run.start_exchange_window()
run._stream.synchronize.side_effect = RuntimeError("device failure")
with self.assertRaisesRegex(RuntimeError, "device failure"):
run.step(3)
self.assertEqual(run.step_count, 2) # value at entry to the failing call
self.assertFalse(run._exchange_valid)
self.assertFalse(run._window_valid)
for read in (lambda: run.read_raw_exchange(0, level=BODY_LEVEL),
lambda: run.read_exchange_window(0, level=BODY_LEVEL), run.download,
lambda: run.step(1)):
with self.assertRaisesRegex(RuntimeError, "failed"):
read()
run._stream.synchronize.side_effect = None
run.close()
def test_launch_error_mid_batch_skips_wait_and_fails(self):
run, log, _ = backend()
run.step(1)
log.clear()
hop, calls = run._hop, []
def failing(index, *args, **kw):
if index == 0:
calls.append(1)
if len(calls) == 2:
raise RuntimeError("launch failure")
return hop(index, *args, **kw)
run._hop = failing
with self.assertRaisesRegex(RuntimeError, "launch failure"):
run.step(3)
self.assertNotIn("sync", names(log))
self.assertEqual(run.step_count, 1)
self.assertFalse(run._exchange_valid)
with self.assertRaisesRegex(RuntimeError, "failed"):
run.read_raw_exchange(0, level=BODY_LEVEL)
run.close()
if __name__ == "__main__":
unittest.main()
+75
View File
@@ -0,0 +1,75 @@
"""Per-hop mass ledger of a fluid cell set A (host-only, stdlib, D2Q9 pull scheme).
Diagnostic definition only; no solver code, no GPU. Arrays are indexed [q][y][x] on a
periodic nx x ny lattice; A is a set of (x, y) ACTIVE cells.
O owner populations before any staging (canonical owner slots, start of the hop)
P populations each cell actually pulls: P[q](x) is what x reads for direction q
(after halo copy, BFL/coalesce staging, ...)
D populations written by the cell after collision
With M0 = sum_{A,q} O and M1 = sum_{A,q} D (both owner snapshots of A):
M1 - M0 + Q_exact = C + H (exact)
C = sum_{x in A, q} (D_q(x) - P_q(x)) collision source
H = sum_{x in A, q} (P_q(x) - O_q(x - c_q)) actual minus ideal owner input,
split into H_seam (x - c_q in A: internal pulls, e.g. tile seams) and
H_edge (x - c_q not in A: inputs crossing the edge of A)
Q_exact = sum_{y in A, y+c_q not in A} O_q(y) - sum_{y not in A, y+c_q in A} O_q(y)
owner-based net outward flux (outgoing minus incoming)
Proof: sum_{A,q} O_q(x - c_q) - M0 = -Q_exact (relabel y = x - c_q per direction), and
sum_{A,q} D = C + H + sum_{A,q} O_q(x - c_q). Q_exact must use O, not staged values: a
flux built from staged inputs already contains H_edge, so adding the full H counts it
twice (see staged_flux). The edge of A must not cross cut links, otherwise O at the solid
end enters Q_exact.
"""
CX = (0, 1, -1, 0, 0, 1, -1, 1, -1)
CY = (0, 0, 0, 1, -1, 1, -1, -1, 1)
def _src(x, y, q, nx, ny):
return (x - CX[q]) % nx, (y - CY[q]) % ny
def pull(O, nx, ny):
"""Ideal pulled input: P[q](x) = O[q](x - c_q)."""
return [[[O[q][_src(x, y, q, nx, ny)[1]][_src(x, y, q, nx, ny)[0]]
for x in range(nx)] for y in range(ny)] for q in range(9)]
def ledger(O, P, D, A, nx, ny):
A = set(A)
out = dict(dM=0.0, C=0.0, H_seam=0.0, H_edge=0.0, Q_exact=0.0)
for (x, y) in A:
for q in range(9):
out["dM"] += D[q][y][x] - O[q][y][x]
out["C"] += D[q][y][x] - P[q][y][x]
sx, sy = _src(x, y, q, nx, ny)
key = "H_seam" if (sx, sy) in A else "H_edge"
out[key] += P[q][y][x] - O[q][sy][sx]
if (sx, sy) not in A:
out["Q_exact"] -= O[q][sy][sx] # incoming owner population
tx, ty = (x + CX[q]) % nx, (y + CY[q]) % ny
if (tx, ty) not in A:
out["Q_exact"] += O[q][y][x] # outgoing owner population
out["H"] = out["H_seam"] + out["H_edge"]
out["residual"] = out["dM"] + out["Q_exact"] - out["C"] - out["H"]
return out
def staged_flux(O, P, A, nx, ny):
"""Outward flux with incoming terms from the staged inputs P (the wrong combination)."""
A = set(A)
q_out = 0.0
for (x, y) in A:
for q in range(9):
sx, sy = _src(x, y, q, nx, ny)
if (sx, sy) not in A:
q_out -= P[q][y][x]
tx, ty = (x + CX[q]) % nx, (y + CY[q]) % ny
if (tx, ty) not in A:
q_out += O[q][y][x]
return q_out
@@ -106,10 +106,9 @@ class ChannelRootGPU(unittest.TestCase):
cfg = LBMConfig(nx=48, ny=32, viscosity=.1, velocity=.015, ddf_shifting=shifted)
reference = None
for block in ((32,48), (8,8), (4,1)):
with FixedPatchStepper(cfg, context=context, coarse_boundary="channel",
refinements=[(8,6,40,26)], level_blocks=[block, (8,8)],
level_objects=[(), (body,) if bodies else ()],
level_jets=[(), (jet,) if bodies else ()]) as backend:
with FixedPatchStepper(cfg, (8,6,40,26), context=context,
coarse_boundary="channel", base_block=block, fine_block=(8,8),
fine_objects=(body,) if bodies else (), jets=(jet,) if bodies else ()) as backend:
arrays = [state(*level.shape[1:]) for level in backend._levels]
backend.initialize(*arrays)
backend.step(4)
@@ -25,8 +25,8 @@ EPS = np.finfo(np.float64).eps
def build(m, window=True, raw=True):
return FixedPatchStepper(config("periodic"), context=m["context"], refinements=[BOUNDS],
level_blocks=[(3, 4), (16, 20)], level_objects=[(), (circle(),)],
return FixedPatchStepper(config("periodic"), BOUNDS, context=m["context"], base_block=(3, 4),
fine_block=(16, 20), fine_objects=[circle()],
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 CelerisLab.cuda.context import CudaContext
with CudaContext(0) as ctx:
kw = dict(context=ctx, refinements=[RB], level_blocks=[None, (16, 20)],
level_objects=[(), (rc(),)], level_jets=[(), (jet(),)],
kw = dict(context=ctx, fine_block=(16, 20), fine_objects=[rc()], jets=[jet()],
observe_raw_exchange=True)
with FixedPatchStepper(CFG, observe_window_exchange=True, **kw) as windowed, \
FixedPatchStepper(CFG, **kw) as plain:
with FixedPatchStepper(CFG, RB, observe_window_exchange=True, **kw) as windowed, \
FixedPatchStepper(CFG, RB, **kw) as plain:
arrays = [state(level.shape) for level in windowed._levels]
windowed.initialize(*arrays)
plain.initialize(*arrays)
@@ -201,11 +200,11 @@ class ExchangeWindowShiftedMRTGPU(unittest.TestCase):
from CelerisLab.cuda.context import CudaContext
cfg = LBMConfig(nx=32, ny=24, viscosity=.1, collision="MRT", ddf_shifting=True)
with CudaContext(0) as ctx:
kw = dict(context=ctx, refinements=[RB], level_blocks=[(12, 8), (16, 20)],
level_objects=[(), (rc(),)], observe_raw_exchange=True)
with FixedPatchStepper(cfg, observe_window_exchange=True, **kw) as windowed, \
FixedPatchStepper(cfg, observe_window_exchange=True, **kw) as split, \
FixedPatchStepper(cfg, **kw) as plain:
kw = dict(context=ctx, base_block=(12, 8), fine_block=(16, 20), fine_objects=[rc()],
observe_raw_exchange=True)
with FixedPatchStepper(cfg, RB, observe_window_exchange=True, **kw) as windowed, \
FixedPatchStepper(cfg, RB, observe_window_exchange=True, **kw) as split, \
FixedPatchStepper(cfg, RB, **kw) as plain:
arrays = [state(level.shape) for level in windowed._levels]
for run in (windowed, split, plain):
run.initialize(*arrays)
@@ -140,8 +140,7 @@ def run_patch(speed, *, audit=True):
build = (patch.object(compiler, "compile_kernel", audit_compiler(compiler.compile_kernel))
if audit else nullcontext())
with build:
backend = implementation.FixedPatchStepper(cfg, context=context, refinements=[BOUNDS],
level_objects=[(), (source,)])
backend = implementation.FixedPatchStepper(cfg, BOUNDS, context=context, fine_objects=[source])
with backend, context.activate():
c, f = backend._levels
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")
with patch.object(cuda, "mem_alloc") as allocate:
with pytest.raises(ValueError, match="donor crosses"):
implementation.FixedPatchStepper(cfg, context=context, refinements=[BOUNDS],
level_objects=[(), (source,)])
implementation.FixedPatchStepper(cfg, BOUNDS, context=context, fine_objects=[source])
context.activate.assert_not_called()
allocate.assert_not_called()
@@ -48,8 +48,9 @@ def make(context):
for (x, y), amplitude in zip(centers, (.003, -.002))
)
run = FixedPatchStepper(
CFG, context=context, refinements=[PARENTS, START], level_blocks=[(4, 4)] * 3,
level_objects=[(), (), objects], level_jets=[(), (), jets],
CFG, PARENTS, context=context, base_block=(4, 4), fine_block=(4, 4),
child_block=(4, 4), child_rectangles_in_parent=START,
child_objects=objects, child_jets=jets,
)
return run, objects
@@ -25,8 +25,9 @@ B = ((16,16,32,40), (60,16,76,40))
def make(ctx):
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)
run = FixedPatchStepper(CFG,context=ctx,refinements=[PARENTS,A],level_blocks=[(4,4)]*3,
level_objects=[(),(),(body,)],level_jets=[(),(),(jet,)])
run = FixedPatchStepper(CFG,PARENTS,context=ctx,base_block=(4,4),fine_block=(4,4),
child_block=(4,4),child_rectangles_in_parent=A,
child_objects=(body,),child_jets=(jet,))
return run,(body,)
@@ -90,18 +90,14 @@ def census(run, index, buffer):
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):
def test_row_ledger_interval_jet_switch_and_no_population_effect(self):
with CudaContext(0) as ctx:
for block in ((32, 40), (16, 20), (8, 10)):
with FixedPatchStepper(CFG, context=ctx, **fine(block), observe_raw_exchange=True) as measured, \
FixedPatchStepper(CFG, context=ctx, **fine(block)) as plain:
with FixedPatchStepper(CFG, BOUNDS, context=ctx, fine_block=block,
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]
measured.initialize(*arrays)
plain.initialize(*arrays)
@@ -150,8 +146,9 @@ class RawExchangeGPU(unittest.TestCase):
"""Same fine resolution/initial data with dense storage: a numerical
reference for cut-link rows, not a uniform-whole-domain convergence study."""
with CudaContext(0) as ctx:
with FixedPatchStepper(CFG, context=ctx, **fine(None)) as dense, \
FixedPatchStepper(CFG, context=ctx, **fine((16, 20)), observe_raw_exchange=True) as blocked:
with FixedPatchStepper(CFG, BOUNDS, context=ctx, fine_objects=[circle()], jets=[jet()]) as dense, \
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]
dense.initialize(*arrays)
blocked.initialize(*arrays)
@@ -232,7 +229,8 @@ class RawExchangeGPU(unittest.TestCase):
layout = BlockLayout2D.build(48, 64, 8, 8)
with CudaContext(0) as ctx:
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))
uniform.initialize(state((9, 48, 64)))
b = amr._level_boundaries[1]
@@ -83,10 +83,11 @@ def mocked_backend(mode="periodic", fail_alloc=None):
class BlockedBaseCPU(unittest.TestCase):
def test_reject_before_cuda(self):
context = Mock(activate=Mock(side_effect=AssertionError("CUDA reached")))
for blocks in ([(12, 8), None], [(7, 8), FINE_BLOCK], [(12,), FINE_BLOCK]):
with self.subTest(blocks=blocks), self.assertRaises(ValueError):
FixedPatchStepper(config("periodic"), context=context, refinements=[BOUNDS],
level_blocks=blocks)
for kwargs in ({"base_block": (12, 8)},
{"base_block": (7, 8), "fine_block": FINE_BLOCK},
{"base_block": (12,), "fine_block": FINE_BLOCK}):
with self.subTest(kwargs=kwargs), self.assertRaises(ValueError):
FixedPatchStepper(config("periodic"), BOUNDS, context=context, **kwargs)
context.activate.assert_not_called()
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):
with mocked_backend("slot_reservoir") as m:
backend = FixedPatchStepper(config("slot_reservoir"), context=m["context"],
coarse_boundary="slot_reservoir", refinements=[BOUNDS],
level_blocks=[(3, 4), FINE_BLOCK],
level_objects=[(), (circle(actuation(.002)),)])
backend = FixedPatchStepper(config("slot_reservoir"), BOUNDS, context=m["context"],
coarse_boundary="slot_reservoir", fine_block=FINE_BLOCK,
base_block=(3, 4), fine_objects=[circle(actuation(.002))])
try:
backend.initialize(state(NY, NX, 1), state(*backend.patch.fine_shape, 2))
coarse, fine = backend._levels
@@ -236,16 +236,16 @@ class BlockedBaseCPU(unittest.TestCase):
for fail_at in (1, 6, 10):
with self.subTest(fail_at=fail_at), mocked_backend(fail_alloc=fail_at) as m:
with self.assertRaisesRegex(MemoryError, "injected allocation"):
FixedPatchStepper(config("periodic"), context=m["context"], refinements=[BOUNDS],
level_blocks=[(3, 4), FINE_BLOCK])
FixedPatchStepper(config("periodic"), BOUNDS, context=m["context"],
fine_block=FINE_BLOCK, base_block=(3, 4))
self.assertEqual(len(m["allocations"]), fail_at)
for allocation in m["allocations"]:
allocation.free.assert_called_once_with()
def test_failed_parent_stage_cannot_publish_or_resume(self):
with mocked_backend() as m:
backend = FixedPatchStepper(config("periodic"), context=m["context"], refinements=[BOUNDS],
level_blocks=[(3, 4), FINE_BLOCK])
backend = FixedPatchStepper(config("periodic"), BOUNDS, context=m["context"],
fine_block=FINE_BLOCK, base_block=(3, 4))
backend.initialize(state(NY, NX, 1), state(*backend.patch.fine_shape, 2))
m["functions"]["BlockCoalesce"].side_effect = RuntimeError("injected parent stage failure")
current = backend._levels[0].current
@@ -289,10 +289,9 @@ class BlockedBaseGPU(unittest.TestCase):
for base_block in (None, *layouts):
label = f"{mode} width={width} speed={speed} base={base_block}"
snapshots = []
with FixedPatchStepper(cfg, context=context, coarse_boundary=mode,
refinements=[BOUNDS], level_blocks=[base_block, FINE_BLOCK],
level_objects=[(), () if speed is None
else (circle(actuation(speed)),)]) as backend:
with FixedPatchStepper(cfg, BOUNDS, context=context, coarse_boundary=mode,
fine_block=FINE_BLOCK, base_block=base_block,
fine_objects=() if speed is None else [circle(actuation(speed))]) as backend:
backend.initialize(*arrays)
for step in (0, 1, 2, 8):
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")
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):
FixedPatchStepper(config("periodic"), context=self.forbidden_context(),
refinements=[BOUNDS], **kw)
FixedPatchStepper(config("periodic"), BOUNDS, context=self.forbidden_context(), **kw)
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).
@@ -126,8 +125,8 @@ class BlockedFineCPU(unittest.TestCase):
patch.object(compiler, "compile_kernel"), patch.object(compiler, "cleanup_build_dir"), \
patch.object(compiler, "load_module", return_value=module), \
patch.object(cuda, "mem_alloc", side_effect=alloc), patch.object(cuda, "memcpy_htod"):
backend = FixedPatchStepper(config("periodic"), context=context, refinements=[BOUNDS],
level_blocks=[None, FINE_BLOCKS[0]])
backend = FixedPatchStepper(config("periodic"), BOUNDS, context=context,
fine_block=FINE_BLOCKS[0])
self.assertEqual(len(backend.fine_layout.halo_dst), 0)
self.assertIsNone(backend._levels[1].halo_dst)
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, "memcpy_htod"), \
patch.object(curved_links, "CurvedLinkSoA", side_effect=AssertionError("dense SoA built")):
backend = FixedPatchStepper(config("periodic"), context=context, refinements=[BOUNDS],
level_blocks=[None, FINE_BLOCKS[1]],
level_objects=[(), tuple(objects)])
backend = FixedPatchStepper(config("periodic"), BOUNDS, context=context,
fine_objects=objects, fine_block=FINE_BLOCKS[1])
backend.initialize(state(NY, NX, 1), state(*backend.patch.fine_shape, 2))
return backend, calls
@@ -291,9 +289,8 @@ class BlockedFineGPU(unittest.TestCase):
results = []
with CudaContext(0) as context:
for fine_block in (None, *FINE_BLOCKS):
with FixedPatchStepper(cfg, context=context, coarse_boundary=mode,
refinements=[BOUNDS],
level_blocks=[None, fine_block]) as backend:
with FixedPatchStepper(cfg, BOUNDS, context=context, coarse_boundary=mode,
fine_block=fine_block) as backend:
backend.initialize(g_coarse, g_fine)
backend.step(STEPS)
self.assertEqual(backend.step_count, STEPS)
@@ -321,9 +318,8 @@ class BlockedFineGPU(unittest.TestCase):
results = []
with CudaContext(0) as context:
for fine_block in (None, *FINE_BLOCKS):
with FixedPatchStepper(cfg, context=context, coarse_boundary=mode,
refinements=[BOUNDS], level_blocks=[None, fine_block],
level_objects=[(), (circle(),)]) as backend:
with FixedPatchStepper(cfg, BOUNDS, context=context, coarse_boundary=mode,
fine_objects=[circle()], fine_block=fine_block) as backend:
backend.initialize(g_coarse, g_fine)
backend.step(STEPS)
results.append((backend.download(), raw(backend)))
@@ -350,9 +346,9 @@ class BlockedFineGPU(unittest.TestCase):
results = []
with CudaContext(0) as context:
for fine_block in (None, *FINE_BLOCKS):
with FixedPatchStepper(cfg, context=context, coarse_boundary=mode,
refinements=[BOUNDS], level_blocks=[None, fine_block],
level_objects=[(), (circle(actuation(speed)),)]) as backend:
with FixedPatchStepper(cfg, BOUNDS, context=context, coarse_boundary=mode,
fine_objects=[circle(actuation(speed))],
fine_block=fine_block) as backend:
backend.initialize(g_coarse, g_fine)
backend.step(STEPS)
results.append((backend.download(), raw(backend)))
@@ -94,8 +94,7 @@ class FaultGuardsCPU(unittest.TestCase):
("slot_reservoir", (12, 4, 20, 21))):
with self.subTest(mode=mode, bounds=bounds):
with self.assertRaises(ValueError):
implementation.FixedPatchStepper(cfg, context=context, coarse_boundary=mode,
refinements=[bounds])
implementation.FixedPatchStepper(cfg, bounds, context=context, coarse_boundary=mode)
context.activate.assert_not_called()
# Both exact channel limits and the original periodic two-cell margin
# must reach context activation, with no allocation or device required.
@@ -103,8 +102,7 @@ class FaultGuardsCPU(unittest.TestCase):
("periodic", (2, 2, 6, 5))):
with self.subTest(mode=mode, bounds=bounds):
with self.assertRaisesRegex(RuntimeError, "validated before CUDA"):
implementation.FixedPatchStepper(cfg, context=context, coarse_boundary=mode,
refinements=[bounds])
implementation.FixedPatchStepper(cfg, bounds, context=context, coarse_boundary=mode)
def fake_backend(self):
backend = PeriodicPatchStepper.__new__(PeriodicPatchStepper)
@@ -156,12 +154,12 @@ class FaultGuardsCPU(unittest.TestCase):
context = Mock(activate=nullcontext, sm_arch="sm_70")
with patch.object(implementation, "HALO", 2):
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))
with patch.object(cuda, "Stream"), patch.object(implementation.FixedPatch2D, "from_config",
return_value=replace(p, tau_coarse=.9)):
with self.assertRaisesRegex(RuntimeError, "omega/tau"):
PeriodicPatchStepper(cfg, context=context, refinements=[p.bounds])
PeriodicPatchStepper(cfg, p.bounds, context=context)
module = Mock()
module.get_global.return_value = (Mock(), 16)
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(cuda, "memcpy_htod") as upload:
with self.assertRaisesRegex(RuntimeError, "layout mismatch"):
PeriodicPatchStepper(cfg, context=context, refinements=[p.bounds])
PeriodicPatchStepper(cfg, p.bounds, context=context)
upload.assert_not_called()
@@ -243,7 +241,7 @@ class NormalBuildGPU(unittest.TestCase):
max_error = max_bc_error = 0.0
with CudaContext(0) as context:
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")
with backend:
self.assertEqual(len(generated), 2)
@@ -313,8 +311,7 @@ class NormalBuildGPU(unittest.TestCase):
for a in (initial.coarse_f, initial.fine_f)))
max_pop = max_drift = max_moment = 0.0
with CudaContext(0) as context, PeriodicPatchStepper(
LBMConfig(nx=16, ny=16, viscosity=.1), context=context,
refinements=[(4, 4, 12, 12)]) as backend:
LBMConfig(nx=16, ny=16, viscosity=.1), (4, 4, 12, 12), context=context) as backend:
backend.initialize(*postcollision(state, spec, maps))
before = totals(backend.download(), maps)
for step in range(1, 9):
@@ -345,7 +342,7 @@ class NormalBuildGPU(unittest.TestCase):
fine_mask[4:-4, 4:-4] = True
initial[1][:, ~fine_mask] = np.nan
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))
p = backend.patch
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, {}),
("blocked_shifted", True,
dict(level_blocks=[(4, 4), (6, 6)])))
dict(base_block=(4, 4), fine_block=(6, 6))))
with CudaContext(0) as context:
for label, shifted, storage in cases:
with self.subTest(storage=label):
cfg = LBMConfig(nx=16, ny=16, viscosity=.0192,
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:
self.assertAlmostEqual(backend.patch.tau_coarse, spec.tau_coarse)
self.assertAlmostEqual(backend.patch.tau_fine, spec.tau_fine)
@@ -148,7 +148,7 @@ class CoupledGPU(unittest.TestCase):
maps = ref.build_maps(spec)
with CudaContext(0) as context:
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))
coarse, fine = backend._levels
@@ -187,7 +187,7 @@ class CoupledGPU(unittest.TestCase):
with CudaContext(0) as context:
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:
# Production allocations: four DDF arrays and four compact tables.
device_bytes = 2 * 9 * (16**2 + 24**2) * 4
@@ -270,7 +270,7 @@ class CoupledGPU(unittest.TestCase):
backend.download()
# 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)
plain.initialize(*rest)
plain.step()
@@ -41,9 +41,8 @@ def mocked(objects, jets, fine_block):
patch.object(compiler, "load_module", return_value=module), \
patch.object(cuda, "mem_alloc", side_effect=lambda size: Mock()), \
patch.object(cuda, "memcpy_htod", htod):
backend = FixedPatchStepper(config("periodic"), context=context, refinements=[BOUNDS],
level_blocks=[None, fine_block], level_objects=[(), tuple(objects)],
level_jets=[(), tuple(jets)])
backend = FixedPatchStepper(config("periodic"), BOUNDS, context=context,
fine_objects=objects, fine_block=fine_block, jets=jets)
return backend, htod
@@ -150,9 +149,8 @@ class JetStepperCPU(unittest.TestCase):
def test_jet_overlapping_legacy_aperture_rejected_before_cuda(self):
with patch.object(cuda, "Stream", side_effect=AssertionError("CUDA touched")):
with self.assertRaisesRegex(ValueError, "SurfaceActuation"):
FixedPatchStepper(config("periodic"), context=Mock(), refinements=[BOUNDS],
level_objects=[(), (circle(actuation(SPEED)),)],
level_jets=[(), (jet(),)])
FixedPatchStepper(config("periodic"), BOUNDS, context=Mock(),
fine_objects=[circle(actuation(SPEED))], jets=[jet()])
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)
patch_ = FixedPatch2D.from_config(cfg, BOUNDS)
with CudaContext(0) as context, FixedPatchStepper(
cfg, context=context, coarse_boundary=mode, refinements=[BOUNDS],
level_blocks=[None, fine_block], level_objects=[(), tuple(objects)],
level_jets=[(), tuple(jets)]) as backend:
cfg, BOUNDS, context=context, coarse_boundary=mode, fine_objects=objects,
fine_block=fine_block, jets=jets) as backend:
backend.initialize(state(NY, NX, 1), state(*patch_.fine_shape, 2))
if switch is None:
backend.step(STEPS)
@@ -42,9 +42,8 @@ class RefinementMRTGPU(unittest.TestCase):
collision="MRT", ddf_shifting=shifted)
reference = None
for block in ((32, 48), (8, 8), (4, 1)):
with FixedPatchStepper(cfg, context=context, coarse_boundary=boundary,
refinements=[(8, 6, 40, 26)], level_blocks=[block, (8, 8)],
level_objects=[(), (body,)]) as backend:
with FixedPatchStepper(cfg, (8, 6, 40, 26), context=context, coarse_boundary=boundary,
base_block=block, fine_block=(8, 8), fine_objects=(body,)) as backend:
backend.initialize(*[state(*level.shape[1:]) for level in backend._levels])
backend.step(8)
raw = backend._download_raw()
@@ -64,8 +63,8 @@ class RefinementMRTGPU(unittest.TestCase):
results = {}
layouts = ((None, None), *((base, FINE_BLOCK) for base in (None, *BASE_BLOCKS)))
for base_block, fine_block in layouts:
with FixedPatchStepper(cfg, context=context, refinements=[BOUNDS],
level_blocks=[base_block, fine_block]) as backend:
with FixedPatchStepper(cfg, BOUNDS, context=context, fine_block=fine_block,
base_block=base_block) as backend:
backend.initialize(*[fixture.state(*level.shape[1:], seed) for seed, level in enumerate(backend._levels)])
backend.step(4)
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):
"""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)
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 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)
for level, level_cfg in enumerate(backend._level_cfgs):
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))
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):
"""Level-wide fine array with each patch window filled from `windows`."""
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):
cfg = config("periodic")
context = object() # never activated: every case must fail first
for kw in (dict(), dict(blocks=(None, (8, 4))), dict(blocks=((4, 4), None)),
dict(blocks=((4, 4), (8, 4)), objects=[circle()]),
dict(blocks=((4, 4), (5, 4)))):
for kw in (dict(), dict(fine_block=(8, 4)), dict(base_block=(4, 4)),
dict(fine_block=(8, 4), base_block=(4, 4), fine_objects=[circle()]),
dict(fine_block=(5, 4), base_block=(4, 4))):
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"):
FixedPatchStepper(cfg, context=context,
**levels(((1, 4, 8, 12), (10, 4, 14, 12)), blocks=((4, 4), (8, 4))))
FixedPatchStepper(cfg, ((1, 4, 8, 12), (10, 4, 14, 12)), context=context,
fine_block=(8, 4), base_block=(4, 4))
from CelerisLab.body.jets import Jet, Parallel, Tube
cfg32 = LBMConfig(nx=32, ny=16, viscosity=.1)
span = Jet(Tube((5., 8.), 0., 22., .6), Parallel(0.), amplitude=0.)
with self.assertRaisesRegex(ValueError, "one jet"): # same obj_id cannot bypass it
FixedPatchStepper(cfg32, context=context,
**levels(((2, 2, 14, 14), (18, 2, 30, 14)), blocks=((4, 4), (4, 4)),
objects=[body(8., 8.), body(24., 8.)], jets=[span]))
FixedPatchStepper(cfg32, ((2, 2, 14, 14), (18, 2, 30, 14)), context=context,
fine_block=(4, 4), base_block=(4, 4),
fine_objects=[body(8., 8.), body(24., 8.)], jets=[span])
class MultipatchGPU(unittest.TestCase):
@@ -73,7 +67,7 @@ class MultipatchGPU(unittest.TestCase):
def run_steps(self, context, cfg, bounds, coarse, fine, **kw):
out = []
with FixedPatchStepper(cfg, context=context, **levels(bounds, **kw)) as backend:
with FixedPatchStepper(cfg, bounds, context=context, **kw) as backend:
backend.initialize(coarse, fine)
for step in self.STEPS:
backend.step(step - backend.step_count)
@@ -88,7 +82,7 @@ class MultipatchGPU(unittest.TestCase):
p = FixedPatch2D.from_config(cfg, BOUNDS)
s = FixedPatchSet2D.from_config(cfg, (BOUNDS,))
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)
new = self.run_steps(context, cfg, (BOUNDS,), coarse, level_fine(s, [fine]), **kw)
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],))
two = FixedPatchSet2D.from_config(cfg32, S1)
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:
ref = self.run_steps(context, cfg16, (S1[0],), coarse, level_fine(one, [fine]), **kw)
tiled = np.concatenate([coarse, coarse], axis=2)
@@ -124,7 +118,7 @@ class MultipatchGPU(unittest.TestCase):
s = FixedPatchSet2D.from_config(cfg, FACE)
coarse = state(16, 24, 5)
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
support = s.fine_layout((4, 4)).owned
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):
out = []
with FixedPatchStepper(cfg, context=context, **levels(bounds, **kw)) as backend:
with FixedPatchStepper(cfg, bounds, context=context, **kw) as backend:
backend.initialize(coarse, fine)
for step in self.STEPS:
if switch is not None and step == 8:
@@ -188,7 +182,8 @@ class MultipatchBodiesGPU(unittest.TestCase):
p = FixedPatch2D.from_config(cfg, BOUNDS)
s = FixedPatchSet2D.from_config(cfg, (BOUNDS,))
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)
new = self.run_steps(context, cfg, (BOUNDS,), coarse, level_fine(s, [fine]), **kw)
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)
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)
kw = dict(blocks=((4, 4), (4, 4)))
kw = dict(fine_block=(4, 4), base_block=(4, 4))
with CudaContext(0) as context:
for jets in ((), (jet(5.),)):
ref = self.run_steps(context, cfg16, self.PATCHES[:1], coarse, level_fine(one, [fine]),
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.))
got = self.run_steps(context, cfg32, self.PATCHES, np.concatenate([coarse, coarse], 2),
level_fine(two, [fine, fine]),
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)
for step, a, b in zip(self.STEPS, ref, got):
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)
L1, L2 = (4, 4, 12, 12), (8, 8, 16, 16)
B = (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()
KW = dict(fine_block=(4, 4), base_block=(4, 4), child=L2, child_block=(4, 4))
def child_body(spec=None):
@@ -61,29 +52,26 @@ def mocked_child():
patch.object(compiler, "load_module", return_value=module), \
patch.object(cuda, "mem_alloc", side_effect=lambda size: Mock()), \
patch.object(cuda, "memcpy_htod"):
return FixedPatchStepper(CFG, context=context,
**levels(objects=((), (), (child_body(),)), jets=((), (), (child_jet(),))))
return FixedPatchStepper(CFG, L1, context=context, child_objects=[child_body()],
child_jets=[child_jet()], **KW)
class ThreeLevelCPU(unittest.TestCase):
def test_rejections_before_cuda(self):
context = object() # never activated: every case must fail first
bad = (levels(blocks=(B, B, None)), levels(blocks=(B, None, B)),
levels(objects=((), (circle(),), ())), levels(l2=(4, 8, 16, 16)),
levels(blocks=(B, B, (5, 4))))
bad = (dict(KW, child_block=None), dict(KW, fine_block=None),
dict(KW, fine_objects=[circle()]), dict(KW, child=(4, 8, 16, 16)),
dict(KW, child_block=(5, 4)))
for kw in bad:
with self.subTest(kw=kw), self.assertRaises(ValueError):
FixedPatchStepper(CFG, context=context, **kw)
# L2 bodies/Jets without an L2 level (one entry too many), and an L2 Jet without a body.
two_level = dict(refinements=[L1], level_blocks=[None, None])
for kw in (dict(two_level, level_objects=[(), (), (child_body(),)]),
dict(two_level, level_jets=[(), (), (child_jet(),)]),
levels(jets=((), (), (child_jet(),)))):
FixedPatchStepper(CFG, L1, context=context, **kw)
for kw in (dict(child_objects=[child_body()]), dict(child_jets=[child_jet()]),
dict(KW, child_jets=[child_jet()])):
with self.assertRaises(ValueError):
FixedPatchStepper(CFG, context=context, **kw)
with self.assertRaises(ValueError): # several L1 patches + scalar child: not yet
FixedPatchStepper(LBMConfig(nx=32, ny=16, viscosity=.1), context=context,
**levels(l1=((4, 4, 12, 12), (20, 4, 28, 12))))
FixedPatchStepper(CFG, L1, context=context, **kw)
with self.assertRaises(ValueError): # several L1 patches + child: not yet
FixedPatchStepper(LBMConfig(nx=32, ny=16, viscosity=.1), ((4, 4, 12, 12), (20, 4, 28, 12)),
context=context, **KW)
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):
from CelerisLab.cuda.context import CudaContext
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)
run.initialize(*(np.broadcast_to(feq, (9, *s)).copy() for s in shapes))
run.step(8)
@@ -181,7 +169,7 @@ class ThreeLevelGPU(unittest.TestCase):
def test_nonequilibrium_masks_and_fresh_middle_feeds(self):
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
run.initialize(state(16, 16, 1), state(*p1.fine_shape, 2), state(*p2.fine_shape, 3))
f = run._levels[1]
@@ -217,9 +205,8 @@ class ThreeLevelBodyGPU(unittest.TestCase):
CFG, velocity=0., wall_normal_velocity=.001, wall_slot_width=4,
y_wall_bc="on_site_velocity", top_wall_bc="reservoir")
with CudaContext(0) as ctx, FixedPatchStepper(
cfg, context=ctx, coarse_boundary=mode,
**levels(blocks=(B, B, block), objects=((), (), (child_body(spec),)),
jets=((), (), tuple(jets)))) as run:
cfg, L1, context=ctx, coarse_boundary=mode, child_objects=[child_body(spec)],
child_jets=jets, **dict(KW, child_block=block)) as run:
run.initialize(*(state(*level.shape[1:], i + 1) for i, level in enumerate(run._levels)))
snapshots = []
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,
mode='periodic', blocks=((4, 4), (4, 4)), objects=(), drive=(), switch=False,
scalar_child=None):
child = children if scalar_child is None else scalar_child
with FixedPatchStepper(config(cfg, mode), context=ctx, coarse_boundary=mode,
refinements=[parents, child], level_blocks=[(4, 4), blocks[0], blocks[1]],
level_objects=[(), (), objects], level_jets=[(), (), drive]) as run:
kwargs = dict(child_rectangles_in_parent=children) if scalar_child is None else dict(child=scalar_child)
with FixedPatchStepper(config(cfg, mode), parents, context=ctx, coarse_boundary=mode,
base_block=(4, 4), fine_block=blocks[0], child_block=blocks[1],
child_objects=objects, child_jets=drive, **kwargs) as run:
frame = frames(run)
run.initialize(*(physical_state(level.shape[1:], *f)
for level, f in zip(run._levels, frame)))
@@ -26,8 +26,9 @@ FP64_ATOL = 3e-6
def make(ctx):
return FixedPatchStepper(CFG, context=ctx, refinements=[PARENTS, A],
level_blocks=[(4, 4)] * 3)
return FixedPatchStepper(CFG, PARENTS, context=ctx, base_block=(4, 4),
fine_block=(4, 4), child_block=(4, 4),
child_rectangles_in_parent=A)
def initial(run, uniform):
@@ -136,8 +136,8 @@ class ShiftedJetsGPU(unittest.TestCase):
Parallel(0.),amplitude=.013)
for rho in (.99,1.01):
with self.subTest(rho=rho),CudaContext(0) as ctx:
with FixedPatchStepper(cfg,context=ctx,refinements=[(8,6,24,18)],level_blocks=[(4,4)]*2,
level_objects=[(),(obj,)],level_jets=[(),(whole,)],observe_raw_exchange=True) as run:
with FixedPatchStepper(cfg,(8,6,24,18),context=ctx,base_block=(4,4),fine_block=(4,4),
fine_objects=(obj,),jets=(whole,),observe_raw_exchange=True) as run:
run._initialize_raw(initial(run,rho,u=.013,v=0.))
expected=run._download_raw();ledger=Ledger(self,run,1)
self.assertTrue(run._level_boundaries[1].tables['aperture'].all())
@@ -153,8 +153,7 @@ class ShiftedJetsGPU(unittest.TestCase):
obj=circle((16.,12.),2.15);j=drive((16.,12.),2.15)
for rho in (.99,1.01):
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,
level_objects=[(),(obj,)],level_jets=[(),(j,)])
settings=dict(bounds=(8,6,24,18),context=ctx,base_block=(4,4),fine_block=(4,4),fine_objects=(obj,),jets=(j,))
measured=owned.enter_context(FixedPatchStepper(cfg,observe_raw_exchange=True,**settings))
plain=owned.enter_context(FixedPatchStepper(cfg,**settings))
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))
obj=circle((12.,14.),.8);j=drive((12.,14.),.8)
with CudaContext(0) as ctx,ExitStack() as owned:
old=owned.enter_context(FixedPatchStepper(cfg,context=ctx,refinements=[parents,a],
level_blocks=[(4,4)]*3,level_objects=[(),(),(obj,)],level_jets=[(),(),(j,)],
observe_raw_exchange=True))
old=owned.enter_context(FixedPatchStepper(cfg,parents,context=ctx,base_block=(4,4),fine_block=(4,4),
child_block=(4,4),child_rectangles_in_parent=a,child_objects=(obj,),child_jets=(j,),observe_raw_exchange=True))
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()
before=old._download_raw();count=old.step_count
@@ -22,8 +22,9 @@ B=((20,16,36,40),)
def make(ctx):
return FixedPatchStepper(CFG,context=ctx,refinements=[PARENTS,A],
level_blocks=[(4,4)]*3)
return FixedPatchStepper(CFG,PARENTS,context=ctx,base_block=(4,4),
fine_block=(4,4),child_block=(4,4),
child_rectangles_in_parent=A)
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)
context = stack.enter_context(CudaContext(0))
backend = stack.enter_context(FixedPatchStepper(
cfg, context=context, coarse_boundary="channel", refinements=[rr.PLACEMENTS["wide"]],
level_blocks=[rr.BASE_BLOCK, rr.FINE_BLOCK], level_objects=[(), (body,)],
cfg, rr.PLACEMENTS["wide"], context=context, coarse_boundary="channel",
base_block=rr.BASE_BLOCK, fine_block=rr.FINE_BLOCK, fine_objects=(body,),
observe_raw_exchange=True, observe_window_exchange=True))
from CelerisLab.lbm.channel import channel_masks
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,
FixedPatchStepper(
cfg,
c["bounds"],
context=ctx,
refinements=[c["bounds"]],
level_blocks=[(8, 8), (8, 8)],
level_objects=[(), (obj,)],
level_jets=[(), tuple(jets)],
base_block=(8, 8),
fine_block=(8, 8),
fine_objects=[obj],
jets=jets,
observe_raw_exchange=True,
) as amr,
):
+5 -4
View File
@@ -383,12 +383,13 @@ def run_case(a, h, periodic, case):
backend = stack.enter_context(
FixedPatchStepper(
cfg,
PLACEMENTS[a.placement],
context=context,
coarse_boundary="channel",
refinements=[PLACEMENTS[a.placement]],
level_blocks=[BASE_BLOCK, FINE_BLOCK],
level_objects=[(), (body,)],
level_jets=[(), public_jets((80.5, 81.0), 20.0, a.jet_scenario)],
base_block=BASE_BLOCK,
fine_block=FINE_BLOCK,
fine_objects=(body,),
jets=public_jets((80.5, 81.0), 20.0, a.jet_scenario),
observe_raw_exchange=True,
**({"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")
context = stack.enter_context(CudaContext(0))
backend = stack.enter_context(FixedPatchStepper(
cfg, context=context, coarse_boundary="channel", refinements=[r.BOUNDS],
level_blocks=[r.BASE_BLOCK, r.FINE_BLOCK], level_objects=[(), (body,)],
observe_raw_exchange=True))
cfg, r.BOUNDS, context=context, coarse_boundary="channel", base_block=r.BASE_BLOCK,
fine_block=r.FINE_BLOCK, fine_objects=(body,), observe_raw_exchange=True))
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)]
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")
context = stack.enter_context(CudaContext(0))
backend = stack.enter_context(FixedPatchStepper(
cfg, context=context, coarse_boundary="channel", refinements=[r.BOUNDS],
level_blocks=[r.BASE_BLOCK, r.FINE_BLOCK], level_objects=[(), (body,)], observe_raw_exchange=True,
cfg, r.BOUNDS, context=context, coarse_boundary="channel", base_block=r.BASE_BLOCK,
fine_block=r.FINE_BLOCK, fine_objects=(body,), observe_raw_exchange=True,
observe_window_exchange=True))
masks = [channel_masks(pair.coarse_shape)[0] & ~pair.covered,
(pair.fine_kind == ACTIVE) & ((host.flags.reshape(pair.fine_shape) & SOLID) == 0)]