Compare commits

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

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

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

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

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

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-28 10:55:02 +00:00
ClaudeandClaude Opus 5.5 65339ab54e test_runtime_child_remap: inject remap body rejects through the boundary record
The removed _boundaries list made the two body-reject patches no-ops.
Inject them via _device_boundaries (level-2 record with foreign host
geometry; non-None level-1 record). set_jet_amplitude docstring names
the public jets view.

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

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

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

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-28 04:45:48 +00:00
Frank14f 7b652b2423 test: construct shear remap fixture through per-level API 2026-10-28 10:21:14 +08:00
Frank14f 7482d95a42 Update channel CPU shim include and retain frozen projected driver 2026-10-28 09:57:15 +08:00
Frank14f 55067652e4 Migrate current callers to level-indexed refinement API 2026-10-28 09:53:20 +08:00
49 changed files with 435 additions and 340 deletions
+7 -7
View File
@@ -117,8 +117,8 @@ with CudaContext(0) as ctx, FixedPatchStepper(cfg, context=ctx, refinements=[(8,
level_objects=[(), (body,)]) 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)),
np.broadcast_to(rest, (9, *run.patch.fine_shape)))
run.initialize(np.broadcast_to(rest, (9, *run.hierarchy.pairs[0].coarse_shape)),
np.broadcast_to(rest, (9, *run.hierarchy.pairs[0].fine_shape)))
run.step(128) # coarse steps; fine level takes two per step
coarse_g, fine_g = run.download() # NaN on covered, ghost and solid cells
```
@@ -143,7 +143,7 @@ Scope and limits:
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
`run.hierarchy.pairs[0].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
@@ -159,8 +159,8 @@ Scope and limits:
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:
run.initialize(np.broadcast_to(rest, (9, *run.patch.coarse_shape)),
np.broadcast_to(rest, (9, *run.patch.fine_shape))) # (9, 32, 64)
run.initialize(np.broadcast_to(rest, (9, *run.hierarchy.pairs[0].coarse_shape)),
np.broadcast_to(rest, (9, *run.hierarchy.pairs[0].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
@@ -193,10 +193,10 @@ Scope and limits:
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:
shapes = (run.patch.coarse_shape, run.patch.fine_shape, run.child_patch.fine_shape)
shapes = (run.hierarchy.pairs[0].coarse_shape, run.hierarchy.pairs[0].fine_shape, run.hierarchy.pairs[1].fine_shape)
run.initialize(*(np.broadcast_to(rest, (9, *s)) for s in shapes)) # 16², 24², 24²
run.step(4)
run.set_jet_amplitude(0, -.002, level=2) # default level=1 still addresses L1
run.set_jet_amplitude(0, -.002, level=2) # explicit level selects the driven level
run.step(4)
coarse_g, l1_g, l2_g = run.download() # solids/covered/nonleaves are NaN
```
+49 -13
View File
@@ -1,8 +1,11 @@
# src coherence stage (from main 914906e)
Owner direction 2026-10-28: unified style over historical compatibility, inside `src`.
No numerical change is intended; kernel text was moved verbatim (PTX/bit parity is
pending a GPU run).
No numerical change is intended. Stages after 5f69eba: 65339ab (one boundary/Jet record,
host preparation helpers, pair accessors), c5a3cdb (`block_collide` families), 442979b
(exchange metadata on the boundary record). Each was checked by parent/candidate
config.h/PTX comparison and host/GPU equality runs; c5a3cdb changed the PTX of the two
channel-root SRT modules only (omega load/clamp moved, outputs exact).
## Conventions (defined once, in `lbm.hierarchy.Hierarchy`)
@@ -28,20 +31,25 @@ pending a GPU run).
## FixedPatchStepper migration map
Historical: left column = names at or after 914906e that no longer exist; right column =
the current name.
| removed / changed | replacement |
|---|---|
| `patch`, `child_patch` | `hierarchy.pairs[0]`, `hierarchy.pairs[1]` |
| `base_layout`, `fine_layout`, `child_layout` (stepper properties) | `layouts[0]`, `layouts[1]`, `layouts[2]` (`FixedPatchSet2D.fine_layout(block)` is unchanged) |
| `fine_boundary`, `child_boundary` | `boundaries[1]`, `boundaries[2]` (host geometry; `None` = no bodies/Jets) |
| `_jets` | `jets[1]` |
| `_boundaries[l]` (dict from L1, level keys) | `_device_boundaries[l]` (list from L0, `None` at L0 unless single-level static circles) |
| `_level_boundaries` (private, b4eacf7) | `_device_boundaries` |
| `_fine_tables`, `_actuated`, `_jet_cache` | `_device_boundaries[1].tables` / `.actuated` / `.jet_cache` |
| `_boundaries[l]` (dict from L1, level keys), `_level_boundaries` (b4eacf7), `_device_boundaries` (65339ab) | `_boundaries[l]`: the one combined per-level `_BlockedBoundary` record (list from L0, `None` at L0 unless single-level static circles): host geometry `.host`, `.tables`, `.jets`, `.jet_cache`, exchange `.spacing`/`.centers`, device handles |
| `_fine_tables`, `_actuated`, `_jet_cache` | `_boundaries[1].tables` / `.actuated` / `.jet_cache` |
| `_level_jets` (private) | `_boundaries[l].jets` (public `jets[l]` is derived) |
| `_exchange_spacing[l]`, `_exchange_centers[l]` (private) | `_boundaries[l].spacing` (h_l, root units) / `.centers` (root edge frame); only levels with a boundary report exchange |
| `_level_blocks` (private) | `_layouts[l].block` |
| `_remap_cfg` (private) | `_root_cfg` |
| `refinement_tables.blocked_base_tables` (also via `refinement_stepper`) | `refinement_tables.blocked_coarse_tables` (import from `refinement_tables`; the stepper no longer re-exports) |
| `_fine_leaf`, `_child_leaf` | `leaf[1]`, `leaf[2]` (copies) |
| `_middle_kind` | `hierarchy.kinds[1]` |
| `_boundary_hosts` (private) | `_boundaries` (private list behind `boundaries`) |
| `_boundary_hosts` (private) | public `boundaries[l]` = `_boundaries[l].host` (host geometry only; the private record is the combined one) |
| `remap_child(rects, child_objects=...)` | `remap_finest(rects, objects=...)` |
| `set_jet_amplitude(i, a)` (implicit level 1) | `set_jet_amplitude(i, a, level=l)` (required) |
| `read_raw_exchange(b)` / `read_exchange_window(b)` (implicit level 1) | `read_raw_exchange(b, level=l)` / `read_exchange_window(b, level=l)` (required) |
@@ -56,13 +64,35 @@ the replaced level relative to its parent).
## Responsibilities
- Construction: `_check_levels` (inputs), `_kernel_configs` (per-level kernel configs),
`Hierarchy.edge_frame` (geometry frames), `_build_level` / `_build_pair` / `_upload_boundary`
(device). Observation reads share `_observed_boundary` + `_summarize`.
- Construction: `_check_levels` (inputs); host preparation, numpy only, before any device
work: `_level_layouts`, `_kernel_configs` (per-level kernel configs), `_host_boundaries`,
`_prepare_blocked_boundary`, `_fluid_leaves`; `Hierarchy.edge_frame` (geometry frames);
then `_planned_bytes()` (memory preflight) and `_build_level` / `_build_pair` /
`_upload_boundary` (device). A construction failure closes the partial instance.
Observation reads share `_observed_boundary` + `_summarize`.
- Pair accessors: `rectangles` (a scalar FixedPatch2D gives `(bounds,)`) and
`fine_level_offset` (`(0, 0)` for a set) on both pair types; the types and their
storage frames stay distinct.
- Per-level state, by lifecycle:
- Fixed plan structure (built once in `__init__`, read by `_planned_bytes()` before any
allocation, readable after `close()` through the public views): `_layouts`,
`_level_cfgs`, `_leaf`, `_state`, and in each `_boundaries[l]` record the host geometry
`.host`, link tables, `.jet_cache`, `.spacing`, `.centers`.
- Mutable boundary data in the same record: `.jets` and `.tables["surface_velocity"]`,
changed only by `set_jet_amplitude` (host prepare -> device write -> commit; an uncertain
write fails the instance); device handles (`.flags`, `.device`, `.kernel`, ledgers) are
set by `_upload_boundary` and cleared by `close()`.
- Runtime: `_levels` (`_Level`) holds compiled modules, kernels, the two population
buffers (swapped per step) and transfer/kind handles. It is created only after the
memory preflight and cleared on `close()`; `.layout`/`.boundary` reference the plan for
`_hop`. Folding the plan into `_Level` would mix both lifecycles in one record.
- Blocked kernels: `operators/block_collision.cuh block_collide` is the one collision extension
point of `BlockStep` and `BlockChannelStep`. MRT shares one path; SRT keeps its two explicit
arithmetic variants (level: unclamped, no-force, shifted direct-momentum feq; channel: uniform
OneStep arithmetic with omega clamp), selected by `BLOCK_ROOT_BC`.
point of `BlockStep` and `BlockChannelStep`. It computes rho/u once and reads no runtime
force (`d_params.fx/fy`). Two compile-time families: uniform-equivalent (MRT on any level,
SRT on the channel root `BLOCK_ROOT_BC == 2`): clamped omega, `collide_mrt_no_force` or
`collide_srt` with its zero Guo term; and level SRT (`BLOCK_ROOT_BC` 0/1): unclamped omega,
`collide_srt_no_force`, with the shifted direct-momentum feq under `USE_DDF_SHIFTING`.
The uniform path keeps `collide_dispatch` (runtime force, LES, TRT).
- `explosion="projected_linear"` stays an optional capability; `"constant"` stays default.
## Not done / open
@@ -72,6 +102,12 @@ the replaced level relative to its parent).
`FixedPatchStepper._kernel_configs` (reset only on level 1 of a slot_reservoir root;
generates each level's config.h). Folding them needs a config.h/PTX diff gate.
- No single-`Simulation` facade over FixedPatchStepper: the two own different state models
(uniform node frame, body sync, forcing, checkpoint vs. immutable blocked hierarchy, edge
(uniform node frame, body sync, forcing, checkpoint vs. fixed blocked hierarchy, edge
frame, remap-as-new-instance). A facade would only forward calls; no concrete need found.
- Deferred, each needing its own evidence first: sharing the channel-root cell/donor selector
with OneStep (the closure functions are already shared; a shared selector changes uniform
signatures and needs a uniform PTX gate); merging FixedPatch2D into FixedPatchSet2D
(different origin, ghost band, layout size and addressing; needs a geometry/parity table and
would need tiles that are not grid-aligned); one `clamp_omega` helper for uniform and
blocked paths (it lives in the hash-pinned LES header `turbulence_smag.cuh`).
- Tests, examples and validation callers are migrated separately (pinned runners stay frozen).
+1 -1
View File
@@ -130,7 +130,7 @@ coverage and two-hop BFL transfer-footprint checks still apply. This is manual
remapping, not adaptive body tracking or runtime body editing.
Supply all original root edge-frame definitions in their construction order.
Do not pass `current.child_boundary.objects`: those are projected fine-local copies.
Do not pass `current.boundaries[2].objects`: those are projected fine-local copies.
There is no retained mutable original-object registry. Before reading device
state, remapping reprojects the supplied definition onto the old topology and
compares exact ordered circle centers/radii, flags, physical cut links, body
+3 -3
View File
@@ -40,13 +40,13 @@ def config(mode):
def leaf_masks(run, mode):
"""Known fluid leaves: uncovered coarse cells (minus wall rows in slot mode)
and fine ACTIVE cells outside the body."""
p = run.patch
p = run.hierarchy.pairs[0]
x0, y0, x1, y1 = p.bounds
coarse = np.ones(p.coarse_shape, bool)
coarse[y0:y1, x0:x1] = False
if mode == "slot_reservoir":
coarse[[0, -1]] = False
solid = (run.fine_boundary.flags.reshape(p.fine_shape) & SOLID) != 0
solid = (run.boundaries[1].flags.reshape(p.fine_shape) & SOLID) != 0
return coarse, (p.fine_kind == ACTIVE) & ~solid
@@ -94,7 +94,7 @@ def main():
cfg, context=context, coarse_boundary=args.mode, refinements=[BOUNDS],
level_blocks=[tuple(args.base_block), tuple(args.fine_block)],
level_objects=[(), (body,)]) as run:
p = run.patch
p = run.hierarchy.pairs[0]
# Rest equilibrium is already postcollision G, so no initial collision is
# needed. Non-leaf entries (covered, ghost, solid) are discarded by initialize.
rest = np.asarray(compute_feq_d2q9(1.0, 0.0, 0.0), np.float32)[:, None, None]
+1 -1
View File
@@ -153,7 +153,7 @@ def main():
with CudaContext(0) as context, FixedPatchStepper(cfg, context=context,
coarse_boundary='slot_reservoir', refinements=[BOUNDS],
level_blocks=LEVEL_BLOCKS) as run:
p = run.patch
p = run.hierarchy.pairs[0]
masks = [np.ones(p.coarse_shape, bool), p.fine_kind == ACTIVE]
masks[0][[0, -1]] = False
masks[0][4:28, 16:32] = False
+1 -1
View File
@@ -50,7 +50,7 @@ def main():
store_precision="FP32", ddf_shifting=False, les_enabled=False)
with CudaContext(0) as context, FixedPatchStepper(cfg, context=context, refinements=[(4, 4, 12, 12)],
level_blocks=[(16, 16), (24, 24)]) as run:
p = run.patch
p = run.hierarchy.pairs[0]
coarse = shear_equilibrium(p.coarse_shape, (0, 0), 1.0, cfg.ny)
fine = shear_equilibrium(p.fine_shape, p.fine_origin, .5, cfg.ny)
coarse_leaf = np.ones(p.coarse_shape, dtype=bool)
+3 -1
View File
@@ -18,7 +18,9 @@ all NQ populations from an owner core cell into another block's halo cell.
Periodicity is explicit per axis. Across a non-periodic exterior, halo cells
have no pair and stay invalid (NaN after packing): no wrap, no implied BC.
Preparation only: no CUDA, no boundary conditions, no refinement levels.
Host preparation plus thin launch wrappers: no CUDA import or allocation
(launch_* only call kernels and device tables the caller owns), no refinement
levels; level_config and static_circle_tables prepare one blocked level.
The one-cell same-level halo here is unrelated to the four-cell fine halo of
the fixed refinement patch (refinement.py), whose semantics are unchanged.
+1 -1
View File
@@ -122,7 +122,7 @@ def leaf_partition(h):
for level, kind in enumerate(h.kinds):
if level:
pair = h.pairs[level - 1]
py, px = (pair.fine_level_offset if isinstance(pair, FixedPatch2D) else (0, 0))
py, px = pair.fine_level_offset
oy, ox = 2 * oy + py, 2 * ox + px
y, x = np.nonzero(kind == ACTIVE)
np.add.at(base, ((y + oy) >> level, (x + ox) >> level), 0.25 ** level)
@@ -1,10 +1,12 @@
// CelerisLab – operators/block_collision.cuh
// The one collision extension point of the blocked steps (BlockStep, BlockChannelStep).
// Compile-time selection: MRT shares one path; SRT keeps two explicit arithmetic
// variants, unchanged from their original kernels:
// level (BLOCK_ROOT_BC 0/1): unclamped omega, no forcing term; shifted storage
// builds feq from the retained drho and direct momentum (numerical variant);
// channel (BLOCK_ROOT_BC 2): matches uniform OneStep SRT arithmetic (omega clamp,
// zero Guo forcing term, compute_feq).
// Unforced by construction: FixedPatchStepper supplies no forcing, so no runtime
// d_params.fx/fy branch is read. Two compile-time arithmetic families:
// uniform-equivalent (MRT on every level; SRT on the channel root, BLOCK_ROOT_BC 2):
// collide_dispatch arithmetic with zero force and no LES: clamped omega;
// SRT keeps the zero Guo forcing term of collide_srt;
// level SRT (BLOCK_ROOT_BC 0/1): unclamped omega, no forcing term; shifted storage
// builds feq from the retained drho and direct momentum (numerical variant).
#ifndef CELERIS_OPERATORS_BLOCK_COLLISION_CUH
#define CELERIS_OPERATORS_BLOCK_COLLISION_CUH
@@ -16,17 +18,19 @@
#endif
__device__ __forceinline__ void block_collide(float* f) {
#if COLLISION_MODEL == 2
collide_mrt_unforced(f);
#elif BLOCK_ROOT_BC == 2
float rho, ux, uy, feq[NQ], forcing[NQ] = {};
float rho, ux, uy;
compute_rho_u(f, rho, ux, uy);
compute_feq(rho, ux, uy, feq);
#if COLLISION_MODEL == 2 || BLOCK_ROOT_BC == 2
const float omega = fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, d_params.omega));
#if COLLISION_MODEL == 2
collide_mrt_no_force(f, rho, ux, uy, omega);
#else
float feq[NQ], forcing[NQ] = {};
compute_feq(rho, ux, uy, feq);
collide_srt(f, feq, forcing, omega);
#endif
#else
float rho, ux, uy, feq[NQ];
compute_rho_u(f, rho, ux, uy);
float feq[NQ];
#if USE_DDF_SHIFTING
// Preserve the density perturbation before adding the unit reference.
const float drho = f[0]+f[1]+f[2]+f[3]+f[4]+f[5]+f[6]+f[7]+f[8];
@@ -114,16 +114,6 @@ __device__ __forceinline__ void collide_mrt_no_force(float* __restrict__ g,
collide_mrt(g, rho, ux, uy, Fin, omega);
}
// Unforced MRT with the clamped device omega, as in collide_dispatch with zero
// forcing. Shared by the blocked, channel-root and dense-patch collisions.
__device__ __forceinline__ void collide_mrt_unforced(float* __restrict__ g)
{
float rho, ux, uy;
compute_rho_u(g, rho, ux, uy);
collide_mrt_no_force(g, rho, ux, uy,
fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, d_params.omega)));
}
#elif LATTICE_MODEL == LATTICE_D3Q19
__device__ __forceinline__ void collide_mrt(float* __restrict__ g,
+12
View File
@@ -120,6 +120,11 @@ class FixedPatch2D:
(an index offset, not an (x, y) origin)."""
return 2 * self.bounds[1] - HALO, 2 * self.bounds[0] - HALO
@property
def rectangles(self) -> tuple:
"""The one parent rectangle, as a sequence like FixedPatchSet2D.rectangles."""
return (self.bounds,)
@property
def covered(self) -> np.ndarray:
"""Dense coarse mask of cells replaced by the patch (not leaves)."""
@@ -225,11 +230,18 @@ class FixedPatchSet2D:
def bounds(self) -> tuple:
return tuple(p.bounds for p in self.patches)
rectangles = bounds
@property
def fine_origin(self) -> tuple[float, float]:
"""Lower edge of the level-wide fine grid in coarse physical units."""
return 0.0, 0.0
@property
def fine_level_offset(self) -> tuple[int, int]:
"""The level-wide grid itself: zero (y, x) index offset."""
return 0, 0
@property
def covered(self) -> np.ndarray:
mask = np.zeros(self.coarse_shape, bool)
+109 -84
View File
@@ -111,6 +111,57 @@ def _check_levels(refinements, blocks, objects, jets):
raise ValueError("bodies/Jets are supported only on the finest level of a nested hierarchy")
# Host preparation (numpy only, before any device work); each returns one entry per level.
def _level_layouts(cfg, pairs, level_blocks, wall_aware):
"""Blocked storage layout of every level, checked against the launch limits."""
layouts = []
for level, block in enumerate(level_blocks):
if level == 0:
layout = BlockLayout2D.build(cfg.ny, cfg.nx, *block, periodic=(True, not wall_aware))
elif isinstance(pairs[level - 1], FixedPatchSet2D):
layout = pairs[level - 1].fine_layout(block)
else:
layout = BlockLayout2D.build(*pairs[level - 1].fine_shape, *block,
periodic=(False, False))
if layout.nblocks > MAX_GRID_DIM or layout.local_shape[0] > MAX_GRID_DIM:
raise ValueError(f"level {level} blocks exceed launch limits")
layouts.append(layout)
return layouts
def _host_boundaries(hierarchy, level_cfgs, level_objects, level_jets, wall_aware, shifted):
"""Host boundary geometry of every refined level in the root edge frame (None at L0)."""
hosts = [None]
for level in range(1, hierarchy.depth):
pair = hierarchy.pairs[level - 1]
if level == 1: # the root config carries the wall-aware boundary fields
host = prepare_fine_boundary(level_cfgs[0], pair, level_objects[1],
wall_aware=wall_aware, jets=level_jets[1],
_shifted_blocked=shifted)
else:
origin, spacing = hierarchy.edge_frame(level - 1)
host = prepare_fine_boundary(level_cfgs[level - 1], pair, level_objects[level],
jets=level_jets[level], parent_origin=origin,
parent_spacing=spacing, _shifted_blocked=shifted)
hosts.append(host)
return hosts
def _fluid_leaves(hierarchy, boundaries):
"""ACTIVE cells of every level that no static circle or body makes solid."""
leaves = [kind == ACTIVE for kind in hierarchy.kinds]
for level, b in enumerate(boundaries):
if b is None:
continue
if b.host is None: # single level: static root circles
leaves[level] &= ~b.tables["solid"]
else:
fine_shape = hierarchy.pairs[level - 1].fine_shape
leaves[level] &= (b.host.flags.reshape(fine_shape) & SOLID) == 0
return leaves
class FixedPatchStepper:
"""D2Q9 SRT/MRT FP32 fixed 2:1 refinement hierarchy, built from level lists.
@@ -246,40 +297,15 @@ class FixedPatchStepper:
self.hierarchy = build_hierarchy(cfg, refinements, wall_aware=self._wall_aware)
pairs = self.hierarchy.pairs
if self._wall_aware:
rects = (pairs[0].bounds,) if isinstance(pairs[0], FixedPatch2D) else pairs[0].bounds
if any(r[1] < 4 or r[3] > cfg.ny - 4 for r in rects):
if any(r[1] < 4 or r[3] > cfg.ny - 4 for r in pairs[0].rectangles):
raise ValueError("slot_reservoir requires y0 >= 4 and y1 <= ny-4")
self._layouts = []
for level, block in enumerate(level_blocks):
if level == 0:
layout = BlockLayout2D.build(cfg.ny, cfg.nx, *block,
periodic=(True, not self._wall_aware))
elif isinstance(pairs[level - 1], FixedPatchSet2D):
layout = pairs[level - 1].fine_layout(block)
else:
layout = BlockLayout2D.build(*pairs[level - 1].fine_shape, *block,
periodic=(False, False))
if layout.nblocks > MAX_GRID_DIM or layout.local_shape[0] > MAX_GRID_DIM:
raise ValueError(f"level {level} blocks exceed launch limits")
self._layouts.append(layout)
self._layouts = _level_layouts(cfg, pairs, level_blocks, self._wall_aware)
if self._channel:
base = blocked_coarse_tables(pairs[0], self._layouts[0], np.argsort(pairs[0].coalesce_coarse))
validate_channel_transfers(pairs[0], self._layouts[0], base[1])
self._level_cfgs = _kernel_configs(cfg, pairs, self._wall_aware)
# Host boundary geometry per level in the root edge frame (Hierarchy.edge_frame).
self._boundaries = [None]
for level in range(1, depth):
pair = pairs[level - 1]
if level == 1: # the root config carries the wall-aware boundary fields
host = prepare_fine_boundary(cfg, pair, level_objects[1],
wall_aware=self._wall_aware, jets=level_jets[1],
_shifted_blocked=self._ddf_shifting)
else:
origin, spacing = self.hierarchy.edge_frame(level - 1)
host = prepare_fine_boundary(self._level_cfgs[level - 1], pair, level_objects[level],
jets=level_jets[level], parent_origin=origin,
parent_spacing=spacing, _shifted_blocked=self._ddf_shifting)
self._boundaries.append(host)
hosts = _host_boundaries(self.hierarchy, self._level_cfgs, level_objects, level_jets,
self._wall_aware, self._ddf_shifting)
if not isinstance(observe_raw_exchange, bool):
raise ValueError("observe_raw_exchange must be bool")
self._observe_raw_exchange = observe_raw_exchange
@@ -291,32 +317,25 @@ class FixedPatchStepper:
self._window_valid = False
self._window_start = None
self._exchange_valid = False
self._exchange_spacing = tuple(self.hierarchy.edge_frame(l)[1] for l in range(depth))
self._exchange_centers = tuple(tuple(tuple(o.geometry.center) for o in objects)
for objects in level_objects)
self._level_jets = level_jets
root = None
if depth == 1 and level_objects[0]:
tables = static_circle_tables(self._layouts[0], level_objects[0])
root = _BlockedBoundary(None, tables, (), False)
self._device_boundaries = [root] + [
_prepare_blocked_boundary(self._boundaries[level], self._layouts[level],
# The one boundary/Jet record per level; host geometry and Jets are views of it.
self._boundaries = [root] + [
_prepare_blocked_boundary(hosts[level], self._layouts[level],
self._level_cfgs[level], level_jets[level],
_shifted_blocked=self._ddf_shifting)
for level in range(1, depth)]
self._leaf = [kind == ACTIVE for kind in self.hierarchy.kinds]
if root is not None:
self._leaf[0] &= ~root.tables["solid"]
for level in range(1, depth):
host = self._boundaries[level]
if host is not None:
self._leaf[level] &= (host.flags.reshape(pairs[level - 1].fine_shape) & SOLID) == 0
for level, b in enumerate(self._boundaries):
if b is not None: # exchange scaling (dt_l = h_l) and root-edge-frame moment centers
b.spacing = self.hierarchy.edge_frame(level)[1]
b.centers = tuple(tuple(o.geometry.center) for o in level_objects[level])
self._leaf = _fluid_leaves(self.hierarchy, self._boundaries)
self._state = [mask.copy() for mask in self._leaf]
if self._channel:
self._leaf[0] &= channel_masks((cfg.ny, cfg.nx))[0]
# Private constructor inputs for remap_finest.
self._root_cfg = replace(cfg)
self._level_blocks = level_blocks
self._root_cfg = replace(cfg) # private constructor input for remap_finest
self.ctx = context
self._lock = RLock()
self._status = "new"
@@ -371,7 +390,7 @@ class FixedPatchStepper:
for index, layout in enumerate(self._layouts):
total += 2 * 4 * int(np.prod(layout.blocked_shape))
total += layout.halo_dst.nbytes + layout.halo_src.nbytes
b = self._device_boundaries[index]
b = self._boundaries[index]
if b is not None:
total += sum(np.asarray(v).nbytes for k, v in b.tables.items()
if k == "flags" or k in BLOCKED_BOUNDARY_DEVICE_COLUMNS)
@@ -394,7 +413,7 @@ class FixedPatchStepper:
grid = ((config.nx + self._block[0] - 1) // self._block[0], config.ny, layout.nblocks)
level = _Level(shape, module, self._allocate(size), self._allocate(size),
module.get_function(name), grid, layout)
level.boundary = self._device_boundaries[index]
level.boundary = self._boundaries[index]
level.halo_copy = module.get_function("BlockHaloCopy")
if len(layout.halo_dst):
level.halo_dst = self._upload(layout.halo_dst)
@@ -556,19 +575,17 @@ class FixedPatchStepper:
pairs = self.hierarchy.pairs
if (self._wall_aware or self._channel or self.hierarchy.depth != 3
or not isinstance(pairs[1], FixedPatchSet2D)
or self._boundaries[1] is not None
or any(b is not None for b in self._device_boundaries[:2])):
or any(b is not None for b in self._boundaries[:2])):
raise ValueError("remap_finest requires a periodic three-level hierarchy with a "
"level-2 set and no level-0/1 bodies")
cfg = self._root_cfg
objects = tuple(objects)
actual = self._device_boundaries[2]
actual = self._boundaries[2]
jets = () if actual is None else actual.jets
old = RemapTopology.from_config(cfg, pairs[0], pairs[1].bounds,
child_objects=objects, child_jets=jets,
_shifted_blocked=self._ddf_shifting)
if (not same_physical_boundary(old.boundary, self._boundaries[2])
or (self._boundaries[2] is not None and actual is None)):
if not same_physical_boundary(old.boundary, None if actual is None else actual.host):
raise ValueError("objects must reproduce the original root-edge-frame level-2 boundary")
if actual is not None:
velocity = actual.tables.get("surface_velocity")
@@ -585,7 +602,7 @@ class FixedPatchStepper:
or not np.array_equal(old.surface_velocity, new.surface_velocity)):
raise ValueError("remap requires the identical pinned physical boundary and BFL support")
# Check sparse block alignment/launch size before any device read/allocation.
layout = new.child.fine_layout(self._level_blocks[2])
layout = new.child.fine_layout(self._layouts[2].block)
if layout.nblocks > MAX_GRID_DIM or layout.local_shape[0] > MAX_GRID_DIM:
raise ValueError("level 2 blocks exceed launch limits")
with self.ctx.activate():
@@ -603,7 +620,7 @@ class FixedPatchStepper:
phase="postcollision_g", **options)
candidate = type(self)(
cfg, context=self.ctx, refinements=[pairs[0].bounds, new.child.bounds],
level_blocks=list(self._level_blocks),
level_blocks=[layout.block for layout in self._layouts],
level_objects=[(), (), objects], level_jets=[(), (), jets],
observe_raw_exchange=self._observe_raw_exchange,
observe_window_exchange=self._observe_window, explosion=self._explosion)
@@ -622,7 +639,7 @@ class FixedPatchStepper:
return candidate
def set_jet_amplitude(self, index, amplitude, *, level):
"""Replace Jet ``index`` of ``level_jets[level]`` between steps; no geometry/table rebuild.
"""Replace Jet ``index`` of level ``level`` (``jets[level]``) between steps; no table rebuild.
Invalid input raises before any change. A device/sync error marks the
instance failed (as step/initialize do); host state commits only after
@@ -633,27 +650,24 @@ class FixedPatchStepper:
if (isinstance(level, bool) or not isinstance(level, Integral)
or not 1 <= level < self.hierarchy.depth):
raise ValueError("Jet level must name an existing refined level (1 .. depth-1)")
boundary = self._device_boundaries[level]
jets = list(self._level_jets[level])
boundary = self._boundaries[level]
# Host preparation: validate and build the new column without touching state.
jets = list(() if boundary is None else boundary.jets)
jets[index] = jets[index].with_amplitude(amplitude) # validates first
velocity = None
if boundary is not None and boundary.jet_cache is not None:
shape, jet_id = boundary.jet_cache
velocity = boundary.tables["surface_velocity"].copy()
driven = jet_id >= 0
velocity[driven] = jet_link_velocity(jets, shape, jet_id)[driven]
shape, jet_id = boundary.jet_cache
velocity = boundary.tables["surface_velocity"].copy()
driven = jet_id >= 0
velocity[driven] = jet_link_velocity(jets, shape, jet_id)[driven]
# Device action: an uncertain write fails the instance closed.
try:
if velocity is not None:
self._stream.synchronize() # no in-flight step reads the column
cuda.memcpy_htod(boundary.device["surface_velocity"], velocity)
self._stream.synchronize() # no in-flight step reads the column
cuda.memcpy_htod(boundary.device["surface_velocity"], velocity)
except BaseException:
self._status = "failed"
raise
if velocity is not None:
boundary.tables["surface_velocity"] = velocity
if boundary is not None:
boundary.jets = tuple(jets)
self._level_jets[level] = tuple(jets)
# Commit: the one boundary record holds both the column and the Jets.
boundary.tables["surface_velocity"] = velocity
boundary.jets = tuple(jets)
def _blocked_boundary(self, level, buffer, **kw):
b = level.boundary
@@ -696,7 +710,7 @@ class FixedPatchStepper:
b = level.boundary
if b is not None and b.last_raw is not None:
count = len(b.tables["cell"])
h = self._exchange_spacing[index] # acoustic scaling: dt_l = h_l
h = b.spacing # acoustic scaling: dt_l = h_l
b.accumulate(b.last_raw, b.interval_integral, np.uint32(count),
np.float64(h), np.float64(h),
grid=((count + self._block[0] - 1) // self._block[0], 1), **kw)
@@ -740,7 +754,7 @@ class FixedPatchStepper:
self._exchange_valid = False
for _ in range(n):
if self._observe_raw_exchange: # stream-ordered: holds the last interval only
for b in self._device_boundaries:
for b in self._boundaries:
if b is not None:
cuda.memset_d8_async(b.interval_integral, 0,
len(b.tables["cell"]) * 32, self._stream)
@@ -763,7 +777,7 @@ class FixedPatchStepper:
raise RuntimeError("window exchange observation is disabled")
self._window_valid = False
try:
for b in self._device_boundaries:
for b in self._boundaries:
if b is not None and b.window_integral is not None:
cuda.memset_d8(b.window_integral, 0, len(b.tables["cell"]) * 32)
except BaseException:
@@ -776,7 +790,7 @@ class FixedPatchStepper:
if (isinstance(level, bool) or not isinstance(level, Integral)
or not 0 <= level < len(self._levels)):
raise ValueError("level must name an existing level")
b = self._device_boundaries[level]
b = self._boundaries[level]
if (b is None or isinstance(body_id, bool) or not isinstance(body_id, Integral)
or not 0 <= body_id < len(b.host.objects)):
raise ValueError("body_id must name a body on this level")
@@ -784,9 +798,9 @@ class FixedPatchStepper:
def _summarize(self, b, rows, level, body_id, window, **kw):
return summarize_exchange(rows, b.tables["body_id"], b.tables["aperture"],
body_id=body_id, level=level, h=self._exchange_spacing[level],
body_id=body_id, level=level, h=b.spacing,
time_end=self._step_count, window=window,
reference_center=self._exchange_centers[level][body_id], **kw)
reference_center=b.centers[body_id], **kw)
def read_exchange_window(self, body_id, *, level):
"""Non-destructive read of [window start, current step count] in base-lattice units."""
@@ -873,7 +887,7 @@ class FixedPatchStepper:
self._stream.synchronize()
except Exception as exc:
error = exc
for b in getattr(self, "_device_boundaries", ()):
for b in getattr(self, "_boundaries", ()):
if b is not None:
b.flags = b.device = b.kernel = None
b.last_raw = b.interval_integral = b.accumulate = None
@@ -901,9 +915,20 @@ class FixedPatchStepper:
raise
# Read-only per-level state, indexed by level (0 = root). Pairs: hierarchy.pairs[l-1].
layouts = property(lambda self: tuple(self._layouts))
boundaries = property(lambda self: tuple(self._boundaries),
doc="Host boundary geometry per level (None: no bodies/Jets).")
jets = property(lambda self: tuple(self._level_jets))
leaf = property(lambda self: tuple(m.copy() for m in self._leaf),
doc="Per-level (ny, nx) masks of public (non-NaN) download cells.")
@property
def layouts(self):
return tuple(self._layouts)
@property
def boundaries(self):
"""Host boundary geometry per level (None: no bodies/Jets on a refined level)."""
return tuple(None if b is None else b.host for b in self._boundaries)
@property
def jets(self):
return tuple(() if b is None else b.jets for b in self._boundaries)
@property
def leaf(self):
"""Per-level (ny, nx) masks of public (non-NaN) download cells."""
return tuple(m.copy() for m in self._leaf)
+2
View File
@@ -130,6 +130,8 @@ class _BlockedBoundary:
interval_integral: object = None
accumulate: object = None
window_integral: object = None # optional observation-window FP64 ledger
spacing: float = None # level spacing h_l in root units (exchange scaling)
centers: tuple = () # body reference centers in the root edge frame
def _prepare_blocked_boundary(boundary, layout, cfg, jets, *, _shifted_blocked=False):
+1
View File
@@ -42,6 +42,7 @@ inline float __fdiv_rn(float a,float b){return a/b;}
#include "step/one_step_double.cu"
#else
#include "core/cell_kind.cuh"
#include "operators/block_collision.cuh"
#include "step/block_channel.cu"
#endif
int main(){
+7 -7
View File
@@ -102,8 +102,8 @@ def test_simulation_before_initialize_commits_host_only(monkeypatch):
def stepper_state(run):
boundary = run._device_boundaries[2]
return boundary.tables["surface_velocity"].copy(), boundary.jets, run._level_jets[2]
boundary = run._boundaries[2]
return boundary.tables["surface_velocity"].copy(), boundary.jets, run.jets[2]
def test_stepper_invalid_input_writes_nothing_and_stays_ready():
@@ -116,11 +116,11 @@ def test_stepper_invalid_input_writes_nothing_and_stays_ready():
with pytest.raises(error):
run.set_jet_amplitude(*args, **kw)
module.cuda.memcpy_htod.assert_not_called()
np.testing.assert_array_equal(run._device_boundaries[2].tables["surface_velocity"], velocity)
assert (run._device_boundaries[2].jets, run._level_jets[2]) == (jets, level_jets)
np.testing.assert_array_equal(run._boundaries[2].tables["surface_velocity"], velocity)
assert (run._boundaries[2].jets, run.jets[2]) == (jets, level_jets)
assert run._status == "ready"
run.set_jet_amplitude(0, .017, level=2)
assert run._level_jets[2][0].amplitude == .017
assert run.jets[2][0].amplitude == .017
finally:
run.close()
@@ -136,8 +136,8 @@ def test_stepper_device_failure_keeps_host_and_blocks_reuse(where):
run._stream.synchronize.side_effect = RuntimeError("upload failed")
with pytest.raises(RuntimeError, match="upload failed"):
run.set_jet_amplitude(0, .017, level=2)
np.testing.assert_array_equal(run._device_boundaries[2].tables["surface_velocity"], velocity)
assert (run._device_boundaries[2].jets, run._level_jets[2]) == (jets, level_jets)
np.testing.assert_array_equal(run._boundaries[2].tables["surface_velocity"], velocity)
assert (run._boundaries[2].jets, run.jets[2]) == (jets, level_jets)
assert run._status == "failed"
module.cuda.memcpy_htod.side_effect = None
run._stream.synchronize.side_effect = None
+4 -4
View File
@@ -94,7 +94,7 @@ class MultiplePinnedRemap(unittest.TestCase):
run.set_jet_amplitude(0, .017, level=2)
run.set_jet_amplitude(1, -.011, level=2)
before = run.download()
clone = run.remap_child((PINS+MOVED+ADDED)[::-1], child_objects=bodies)
clone = run.remap_finest((PINS+MOVED+ADDED)[::-1], objects=bodies)
try:
self.assertEqual(tuple(j.amplitude for j in clone._boundaries[2].jets), (.017, -.011))
np.testing.assert_array_equal(clone._boundaries[2].tables['surface_velocity'],
@@ -116,7 +116,7 @@ class MultiplePinnedRemap(unittest.TestCase):
start = len(allocations)
module.cuda.memcpy_dtoh.reset_mock()
with self.assertRaisesRegex(ValueError, 'boundary'):
run.remap_child(PINS+FLUID, child_objects=changed)
run.remap_finest(PINS+FLUID, objects=changed)
self.assertEqual(len(allocations), start)
module.cuda.memcpy_dtoh.assert_not_called()
with self.assertRaisesRegex(ValueError, 'pinned'):
@@ -128,7 +128,7 @@ class MultiplePinnedRemap(unittest.TestCase):
cases = [(PINS+FLUID, bodies[:1]), (PINS+FLUID, bodies[::-1]),
(PINS+FLUID, (bodies[0], body(center=(32.1, 10.)))),
(PINS+FLUID, (bodies[0], body(SurfaceActuation(0., .7, 0.), (32., 10.)))),
(PINS+FLUID, run.child_boundary.objects)]
(PINS+FLUID, run.boundaries[2].objects)]
for index in range(2):
changed = list(PINS+FLUID)
x0, y0, x1, y1 = changed[index]
@@ -141,7 +141,7 @@ class MultiplePinnedRemap(unittest.TestCase):
run._stream.synchronize.reset_mock()
for rects, supplied in cases:
with self.subTest(rects=rects), self.assertRaises(ValueError):
run.remap_child(rects, child_objects=supplied)
run.remap_finest(rects, objects=supplied)
self.assertEqual(len(allocations), start)
module.cuda.memcpy_dtoh.assert_not_called()
run._stream.synchronize.assert_not_called()
+8 -8
View File
@@ -61,11 +61,11 @@ class PinnedBodyRemap(unittest.TestCase):
run,module,allocations,storage,objects=ready()
run.set_jet_amplitude(0,.017,level=2)
before=run.download()
clone=run.remap_child(B,child_objects=objects)
clone=run.remap_finest(B,objects=objects)
self.assertEqual(clone._boundaries[2].jets[0].amplitude,.017)
np.testing.assert_array_equal(clone._boundaries[2].tables['surface_velocity'],
run._boundaries[2].tables['surface_velocity'])
self.assertTrue(same_physical_boundary(run.child_boundary,clone.child_boundary))
self.assertTrue(same_physical_boundary(run.boundaries[2],clone.boundaries[2]))
clone.step()
self.assertEqual((run.step_count,clone.step_count),(7,8))
for got,want in zip(run.download(),before): np.testing.assert_array_equal(got,want)
@@ -73,7 +73,7 @@ class PinnedBodyRemap(unittest.TestCase):
def test_static_body_without_actuation(self):
run,module,allocations,storage,objects=ready(jets=())
clone=run.remap_child(A,child_objects=objects)
clone=run.remap_finest(A,objects=objects)
self.assertFalse(clone._boundaries[2].actuated)
self.assertNotIn('surface_velocity',clone._boundaries[2].tables)
for a,b in zip(run.download(),clone.download()): np.testing.assert_array_equal(a,b)
@@ -82,13 +82,13 @@ class PinnedBodyRemap(unittest.TestCase):
def test_reject_before_device_reads_or_allocation(self):
run,module,allocations,storage,objects=ready()
start=len(allocations)
cases=[(B,()),(B,run.child_boundary.objects),(B,(body(center=(12.1,14.)),)),
cases=[(B,()),(B,run.boundaries[2].objects),(B,(body(center=(12.1,14.)),)),
(((24,16,36,40),B[1]),objects),(((16,28,32,40),B[1]),objects),
((B[1],),objects)]
module.cuda.memcpy_dtoh.reset_mock()
for rects,objs in cases:
with self.subTest(rects=rects),self.assertRaises(ValueError):
run.remap_child(rects,child_objects=objs)
run.remap_finest(rects,objects=objs)
self.assertEqual(len(allocations),start)
module.cuda.memcpy_dtoh.assert_not_called()
self.assertEqual((run._status,run.step_count),('ready',7))
@@ -100,9 +100,9 @@ class PinnedBodyRemap(unittest.TestCase):
(body(SurfaceActuation(0.,.7,speed)),),())
for spec in (SurfaceActuation(0.,.7,.003),SurfaceActuation(2.,.7,speed)):
module.cuda.memcpy_dtoh.reset_mock()
with self.assertRaises(ValueError): run.remap_child(B,child_objects=(body(spec),))
with self.assertRaises(ValueError): run.remap_finest(B,objects=(body(spec),))
module.cuda.memcpy_dtoh.assert_not_called()
clone=run.remap_child(B,child_objects=objects)
clone=run.remap_finest(B,objects=objects)
clone.close();run.close()
def test_allocation_failure_keeps_old_ready_and_count(self):
@@ -117,7 +117,7 @@ class PinnedBodyRemap(unittest.TestCase):
return allocate(size)
with patch.object(module.cuda,'mem_alloc',side_effect=fail):
with self.assertRaisesRegex(RuntimeError,'allocation'):
run.remap_child(B,child_objects=objects)
run.remap_finest(B,objects=objects)
self.assertEqual((run._status,run.step_count),('ready',7))
for got,want in zip(run.download(),before): np.testing.assert_array_equal(got,want)
for a in saved: a.free.assert_not_called()
+2 -2
View File
@@ -164,7 +164,7 @@ class ProjectedPairDispatch(unittest.TestCase):
run, host, _, _, storage = backend(refinements=FOUR[:2], explosion=explosion)
run._step_count = 7
old = run.download()
clone = run.remap_child(FOUR[1])
clone = run.remap_finest(FOUR[1])
self.assertEqual(clone._explosion, explosion)
self.assertEqual(clone.step_count, 7)
for pair in clone._pairs:
@@ -182,7 +182,7 @@ class ProjectedPairDispatch(unittest.TestCase):
refinements=FOUR[:2], explosion='projected_linear')
run._step_count = 3
before = run.download()
clone = run.remap_child(changed)
clone = run.remap_finest(changed)
self.assertEqual(clone._explosion, 'projected_linear')
self.assertEqual(clone.hierarchy.pairs[1].bounds, tuple(changed))
self.assertNotEqual(clone._pairs[1].n_ghosts, run._pairs[1].n_ghosts)
+2 -2
View File
@@ -85,7 +85,7 @@ class RawExchange(unittest.TestCase):
def test_four_level_schedule_publication_and_fixed_storage(self):
run, host, calls, arrays = backend()
b = run._device_boundaries[3]
b = run._boundaries[3]
rows = len(b.tables['cell'])
self.assertEqual(b.last_raw.size, rows * 16)
self.assertEqual(b.interval_integral.size, rows * 32)
@@ -128,7 +128,7 @@ class RawExchange(unittest.TestCase):
def test_disabled_has_nulls_and_no_observation_allocation_or_launch(self):
run, host, calls, _ = backend(False)
b = run._device_boundaries[3]
b = run._boundaries[3]
self.assertIs(b.last_raw, b.interval_integral)
self.assertIs(b.interval_integral, None)
run.step()
+5
View File
@@ -524,6 +524,11 @@ def test_public_jet_reversal_boundary_and_windows():
assert r.jet_interval(handle, "reverse", switch + 80, 80, diameter) == (80, None)
assert r.jet_interval(handle, "fixed", switch, 80, diameter) == (80, None)
assert not calls
refined_calls = []
refined = SimpleNamespace(set_jet_amplitude=lambda i, a, *, level: refined_calls.append((i, a, level)))
n, event = r.jet_interval(refined, "reverse", switch, 80, diameter, level=1)
assert n == 80 and event["step"] == switch and event["tstar"] == 300
assert refined_calls == [(0, -.001, 1), (1, .001, 1)]
rows = [dict(tstar=t, tstar_force=t - .5, cd=2.0, cl=-.1) for t in range(251, 306)]
windows = r.jet_diagnostics(rows, "reverse")
assert len(windows) == 6 and all(w["complete"] for w in windows)
+2 -2
View File
@@ -177,7 +177,7 @@ class FixedPatchSetTests(unittest.TestCase):
self.assertEqual((fy, fx), (32, 48))
rows = cells = 0
for p in s.patches:
y0, x0 = p.fine_level_origin
y0, x0 = p.fine_level_offset
ly, lx = p.fine_shape
np.testing.assert_array_equal(s.fine_kind[y0:y0 + ly, x0:x0 + lx], p.fine_kind)
gy, gx = np.divmod(p.ghost_fine.astype(np.int64), lx)
@@ -224,7 +224,7 @@ class FixedPatchSetTests(unittest.TestCase):
layout = s.fine_layout((4, 4))
support = np.zeros(s.fine_shape, bool)
for p in s.patches:
y0, x0 = p.fine_level_origin
y0, x0 = p.fine_level_offset
support[y0:y0 + p.fine_shape[0], x0:x0 + p.fine_shape[1]] = True
np.testing.assert_array_equal(layout.owned, support)
self.assertTrue(layout.owned.ravel()[s.ghost_fine.astype(np.intp)].all())
+16 -9
View File
@@ -1,5 +1,6 @@
"""Host guards and source contract for MRT reuse in refined/blocked stepping (no CUDA)."""
from pathlib import Path
import re
import unittest
from CelerisLab.config import LBMConfig
@@ -31,20 +32,26 @@ class RefinementMRTHost(unittest.TestCase):
"compute_feq(rho, ux, uy, feq);",
"collide_srt_no_force(f, feq, d_params.omega);"):
self.assertIn(text, block_collision)
self.assertIn("collide_srt(f, feq, forcing, omega);", channel)
self.assertIn("collide_srt(f, feq, forcing, omega);", block_collision)
self.assertIn("#if COLLISION_MODEL == 2 || BLOCK_ROOT_BC == 2", block_collision)
self.assertEqual(channel.count("block_collide(f);"), 1)
# Executable text only: the header comment explains the absent force read.
code = re.sub(r"//[^\n]*", "", channel + "\n" + block_collision)
self.assertIn("d_params.omega", code)
for text in ("collide_mrt_unforced", "collide_dispatch(", "d_params.fx", "d_params.fy"):
self.assertNotIn(text, code)
for source in (blocks,):
self.assertIn("(COLLISION_MODEL != 0 && COLLISION_MODEL != 2)", source)
self.assertIn("USE_LES", source)
self.assertIn("STORE_PRECISION != 0", source)
self.assertIn('#if COLLISION_MODEL == 2\n#include "operators/collision_mrt.cuh"\n#endif', source)
# One shared unforced MRT body; each blocked/patch path calls it exactly once.
mrt = (KERNELS / "operators/collision_mrt.cuh").read_text()
self.assertEqual(mrt.count("void collide_mrt_unforced("), 1)
for source in (blocks + "\n" + block_collision, channel):
self.assertEqual(source.count("collide_mrt_unforced(f);"), 1)
self.assertNotIn("collide_mrt_no_force(", source)
self.assertNotIn("compute_rho_u(f, rho, ux, uy);\n collide_mrt", source)
self.assertNotRegex(source, r"compute_rho_u\(f, rho, ux, uy\);\s*collide_mrt")
# One unforced MRT call, in block_collide only; the steps never collide directly.
self.assertEqual(block_collision.count("collide_mrt_no_force(f, rho, ux, uy, omega);"), 1)
self.assertEqual(block_collision.count("compute_rho_u(f, rho, ux, uy);"), 1)
self.assertEqual(block_collision.count("fmaxf(OMEGA_COLLISION_MIN, d_params.omega)"), 1)
for source in (blocks, channel):
self.assertNotIn("collide_mrt", source)
self.assertNotIn("collide_srt", source)
if __name__ == "__main__": unittest.main()
+12 -10
View File
@@ -1,5 +1,6 @@
"""CPU mocks execute production remap lifecycle; no GPU numerical claim."""
from dataclasses import replace
from types import SimpleNamespace
import unittest
from unittest.mock import patch
@@ -29,7 +30,7 @@ class RuntimeChildRemap(unittest.TestCase):
saved = {key: value.copy() for key, value in storage.items()}
old_allocations = list(allocations)
module.cuda.memcpy_dtoh.reset_mock()
clone = run.remap_child(A)
clone = run.remap_finest(A)
self.assertEqual(module.cuda.memcpy_dtoh.call_count, 3)
for a, b in zip(clone.download(), before):
np.testing.assert_array_equal(a, b)
@@ -48,7 +49,7 @@ class RuntimeChildRemap(unittest.TestCase):
run, module, allocations, storage = ready_backend()
before = run.download()
expected = remap(topology(A), topology(B), before).g
clone = run.remap_child(B)
clone = run.remap_finest(B)
for level, want, got in zip(clone._levels, expected, clone.download()):
np.testing.assert_array_equal(got, want)
for ptr in (level.current, level.spare):
@@ -65,17 +66,18 @@ class RuntimeChildRemap(unittest.TestCase):
module.cuda.memcpy_dtoh.reset_mock()
for status in ('new', 'failed', 'closed', 'advancing'):
run._status = status
with self.assertRaises(RuntimeError): run.remap_child(B)
with self.assertRaises(RuntimeError): run.remap_finest(B)
run._status = 'ready'
for bad in ([], [(40,16,48,40)], [(16,16,32,24),(16,27,32,55)]):
with self.assertRaises(ValueError): run.remap_child(bad)
hosts, pairs = run._boundary_hosts, run.hierarchy.pairs
with self.assertRaises(ValueError): run.remap_finest(bad)
records, pairs = run._boundaries, run.hierarchy.pairs
for attribute, value in (('_wall_aware', True),
('_boundary_hosts', [*hosts[:2], object()]), # child boundary
('_boundary_hosts', [hosts[0], object(), hosts[2]]), # fine boundary
# level-2 record whose host geometry differs from objects=()
('_boundaries', [*records[:2], SimpleNamespace(host=object(), jets=())]),
('_boundaries', [records[0], object(), records[2]]), # level-1 body
('hierarchy', replace(run.hierarchy, pairs=(pairs[0], None)))):
with patch.object(run, attribute, value):
with self.assertRaises(ValueError): run.remap_child(B)
with self.assertRaises(ValueError): run.remap_finest(B)
self.assertEqual(len(allocations), start)
module.cuda.memcpy_dtoh.assert_not_called()
run.close()
@@ -102,7 +104,7 @@ class RuntimeChildRemap(unittest.TestCase):
initialize(candidate, *arrays)
operation = patch.object(type(run), 'initialize', fail)
with operation, self.assertRaisesRegex(RuntimeError, 'candidate'):
run.remap_child(B)
run.remap_finest(B)
self.assertEqual((run._status, run.step_count), ('ready', 7))
for a, b in zip(run.download(), before): np.testing.assert_array_equal(a, b)
for allocation in old_allocations: allocation.free.assert_not_called()
@@ -116,7 +118,7 @@ class RuntimeChildRemap(unittest.TestCase):
target = run._stream.synchronize if operation == 'sync' else module.cuda.memcpy_dtoh
target.side_effect = RuntimeError('old device failure')
with self.assertRaisesRegex(RuntimeError, 'old device'):
run.remap_child(B)
run.remap_finest(B)
self.assertEqual(run._status, 'failed')
self.assertEqual(len(allocations), start)
target.side_effect = None
+7 -7
View File
@@ -67,7 +67,7 @@ class ShiftedJets(unittest.TestCase):
def test_switch_upload_only_then_raw_remap_keeps_live_velocity(self):
run,module,allocations,storage,objects=ready();clone=None
try:
boundary=run._device_boundaries[2]
boundary=run._boundaries[2]
before=run._download_raw()
original=boundary.tables['surface_velocity'].copy()
module.cuda.memcpy_htod.reset_mock()
@@ -77,13 +77,13 @@ class ShiftedJets(unittest.TestCase):
np.testing.assert_array_equal(boundary.tables['surface_velocity'],-original)
self.assertEqual(run.step_count,7)
for a,b in zip(before.values,run._download_raw().values):self.assertEqual(a.tobytes(),b.tobytes())
clone=run.remap_child(B,child_objects=objects)
self.assertEqual(clone._level_jets[2][0].amplitude,-jet().amplitude)
np.testing.assert_array_equal(clone._device_boundaries[2].tables['surface_velocity'],
clone=run.remap_finest(B,objects=objects)
self.assertEqual(clone.jets[2][0].amplitude,-jet().amplitude)
np.testing.assert_array_equal(clone._boundaries[2].tables['surface_velocity'],
boundary.tables['surface_velocity'])
self.assertEqual(clone.step_count,7)
saved=run._download_raw()
with self.assertRaises(ValueError):run.remap_child(B,child_objects=(body(center=(12.1,14.)),))
with self.assertRaises(ValueError):run.remap_finest(B,objects=(body(center=(12.1,14.)),))
for a,b in zip(saved.values,run._download_raw().values):self.assertEqual(a.tobytes(),b.tobytes())
finally:
if clone is not None:clone.close()
@@ -92,12 +92,12 @@ class ShiftedJets(unittest.TestCase):
def test_invalid_amplitude_does_not_upload_or_commit(self):
run,module,_,_,_=ready()
try:
before=run._device_boundaries[2].tables['surface_velocity'].copy()
before=run._boundaries[2].tables['surface_velocity'].copy()
module.cuda.memcpy_htod.reset_mock()
for amplitude in (float('nan'),float('inf'),1.):
with self.assertRaises(ValueError):run.set_jet_amplitude(0,amplitude,level=2)
module.cuda.memcpy_htod.assert_not_called()
np.testing.assert_array_equal(before,run._device_boundaries[2].tables['surface_velocity'])
np.testing.assert_array_equal(before,run._boundaries[2].tables['surface_velocity'])
finally:run.close()
+2 -2
View File
@@ -95,7 +95,7 @@ class ShiftedStorage(unittest.TestCase):
before = run._download_raw()
expected = remap(topology(A),topology(B),before.values)
with patch.object(run,'download',side_effect=AssertionError('public decode used')):
clone = run.remap_child(B)
clone = run.remap_finest(B)
self.assertEqual((run.step_count,clone.step_count),(7,7))
self.assertFalse(clone._exchange_valid)
for level, want, got in zip(clone._levels,expected.g,clone._download_raw().values):
@@ -114,7 +114,7 @@ class ShiftedStorage(unittest.TestCase):
physical = tuple(np.broadcast_to(_WEIGHTS, level.shape).copy() for level in run._levels)
for a in physical: a[1] += 2**-35
run.initialize(*physical)
for a,mask in zip(run._download_raw().values,run._leaf):
for a,mask in zip(run._download_raw().values,run.leaf):
self.assertTrue((a[1,mask] == np.float32(2**-35)).all())
self.assertEqual(run.step_count,0)
module.cuda.memcpy_htod.reset_mock()
+1 -1
View File
@@ -73,7 +73,7 @@ class SingleLevel(unittest.TestCase):
def test_circle_solids_are_not_state(self):
run, _, _ = backend((circle(),))
solid = run._device_boundaries[0].tables["solid"]
solid = run._boundaries[0].tables["solid"]
self.assertTrue(solid.any())
np.testing.assert_array_equal(run._state[0], ~solid)
+2 -2
View File
@@ -35,10 +35,10 @@ class UnionBodySupport(unittest.TestCase):
run,module,allocations,storage,objects=ready()
try:
run.set_jet_amplitude(0,.017,level=2)
clone=run.remap_child(((16,16,24,40),(24,16,32,40),A[1]),child_objects=objects)
clone=run.remap_finest(((16,16,24,40),(24,16,32,40),A[1]),objects=objects)
try:
self.assertEqual(clone._boundaries[2].jets[0].amplitude,.017)
self.assertTrue(same_physical_boundary(run.child_boundary,clone.child_boundary))
self.assertTrue(same_physical_boundary(run.boundaries[2],clone.boundaries[2]))
for a,b in zip(run.download(),clone.download()):np.testing.assert_array_equal(a,b)
finally:clone.close()
finally:run.close()
+1 -1
View File
@@ -40,7 +40,7 @@ class Equivalence(unittest.TestCase):
table_sets(FixedPatchSet2D.from_config(CFG, [(8, 8, 20, 22)])))
def test_separated_rectangles_reproduce_current_tables(self):
rects = [(4, 4, 10, 10), (16, 16, 26, 26)] # gap 6 >= MIN_PATCH_GAP
rects = [(4, 4, 10, 10), (16, 16, 26, 26)] # separated rectangles; union dependencies determine legality
self.assertEqual(table_sets(union(*rects)),
table_sets(FixedPatchSet2D.from_config(CFG, rects)))
@@ -240,7 +240,7 @@ class SingleLevelGPU(unittest.TestCase):
with CudaContext(0) as context:
for split in SPLITS:
with single(context=context, block=split, objects=[circle()]) as backend:
tables = backend._device_boundaries[0].tables
tables = backend._boundaries[0].tables
for key, values, dtype in (("q", [l.q for l in links], np.float32),
("dir", [l.dir for l in links], np.uint8),
("fallback", [l.fallback for l in links], np.uint8)):
@@ -113,7 +113,7 @@ class ChannelRootGPU(unittest.TestCase):
arrays = [state(*level.shape[1:]) for level in backend._levels]
backend.initialize(*arrays)
backend.step(4)
if bodies: backend.set_jet_amplitude(0, -.001)
if bodies: backend.set_jet_amplitude(0, -.001, level=1)
backend.step(4)
raw = backend._download_raw()
physical = backend.download()
@@ -128,7 +128,7 @@ class ChannelRootGPU(unittest.TestCase):
else:
for a,b in zip(reference, raw.values): np.testing.assert_array_equal(a,b)
with self.assertRaisesRegex(ValueError, "periodic"):
backend.remap_child([])
backend.remap_finest([])
if __name__ == "__main__": unittest.main()
@@ -59,44 +59,44 @@ class ExchangeWindowCPU(unittest.TestCase):
stack.enter_context(patch.object(cuda, "memcpy_dtoh",
side_effect=lambda out, src: out.fill(0)))
backend = build(m)
b = backend._device_boundaries[1]
b = backend._boundaries[1]
self.assertIsNotNone(b.window_integral)
with self.assertRaises(RuntimeError):
backend.start_exchange_window() # not initialized
arrays = [np.full(level.shape, 1 / 9, np.float32) for level in backend._levels]
backend.initialize(*arrays)
with self.assertRaisesRegex(RuntimeError, "no valid exchange window"):
backend.read_exchange_window(0)
backend.read_exchange_window(0, level=1)
backend.step(1)
self.assertEqual(len(accumulations(m)), 2) # closed window: interval ledger only
backend.start_exchange_window()
with self.assertRaisesRegex(RuntimeError, "at least one completed"):
backend.read_exchange_window(0)
backend.read_exchange_window(0, level=1)
m["calls"].clear()
backend.step(2)
calls = accumulations(m)
self.assertEqual(len(calls), 8) # two fine hops x two ledgers x two base steps
self.assertEqual({id(c[1][1]) for c in calls},
{id(b.interval_integral), id(b.window_integral)})
out = backend.read_exchange_window(0)
out = backend.read_exchange_window(0, level=1)
self.assertEqual((out["time_start"], out["time_end"], out["sample_count"]), (1., 3, 4))
self.assertEqual(backend.read_exchange_window(0)["time_end"], 3) # non-destructive
self.assertEqual(backend.read_exchange_window(0, level=1)["time_end"], 3) # non-destructive
backend.initialize(*arrays) # invalidates
with self.assertRaisesRegex(RuntimeError, "no valid exchange window"):
backend.read_exchange_window(0)
backend.read_exchange_window(0, level=1)
backend.start_exchange_window()
with patch.object(backend, "_hop", side_effect=RuntimeError("injected")):
with self.assertRaises(RuntimeError):
backend.step(1)
self.assertFalse(backend._window_valid)
with self.assertRaises(RuntimeError):
backend.read_exchange_window(0) # failed status
backend.read_exchange_window(0, level=1) # failed status
def test_disabled_default_has_no_window_cost(self):
with mocked_backend() as m, ExitStack() as stack:
stack.enter_context(patch.object(cuda, "memset_d8_async"))
backend = build(m, window=False)
self.assertIsNone(backend._device_boundaries[1].window_integral)
self.assertIsNone(backend._boundaries[1].window_integral)
backend.initialize(*[np.full(l.shape, 1 / 9, np.float32) for l in backend._levels])
backend.step(1)
self.assertEqual(len(accumulations(m)), 2)
@@ -120,7 +120,7 @@ class ExchangeWindowGPU(unittest.TestCase):
windowed.step(2)
plain.step(2)
windowed.start_exchange_window()
b = windowed._device_boundaries[1]
b = windowed._boundaries[1]
count = len(b.tables["cell"])
factors = np.array([.25, .25, .125, .25]) # h*dt, h*h*dt, h*h for level-1 h=.5
scale = np.zeros(4)
@@ -140,12 +140,12 @@ class ExchangeWindowGPU(unittest.TestCase):
summed = {part: [0., 0., 0., 0.] for part in ("total", "wall", "aperture")}
steps = 0
for speed in (0., .002, .002, -.001, 0., 0.): # amplitude updates stay in-window
windowed.set_jet_amplitude(0, speed)
plain.set_jet_amplitude(0, speed)
windowed.set_jet_amplitude(0, speed, level=1)
plain.set_jet_amplitude(0, speed, level=1)
windowed.step(1)
plain.step(1)
steps += 1
obs = windowed.read_raw_exchange(0)
obs = windowed.read_raw_exchange(0, level=1)
for part in summed:
o = obs[part]
summed[part] = [summed[part][0] + o["impulse"][0], summed[part][1] + o["impulse"][1],
@@ -153,7 +153,7 @@ class ExchangeWindowGPU(unittest.TestCase):
summed[part][3] + o["mass_into_fluid"]]
for a, e in zip(windowed._download_raw().values, plain._download_raw().values):
same_bits(self, a, e)
out = windowed.read_exchange_window(0)
out = windowed.read_exchange_window(0, level=1)
self.assertEqual((out["time_end"] - out["time_start"], out["sample_count"]),
(steps, 2 * steps))
self.assertEqual(substeps, 2 * steps)
@@ -165,11 +165,11 @@ class ExchangeWindowGPU(unittest.TestCase):
diff = np.abs(got - np.array(summed[part]))
self.assertTrue((diff <= bound).all(), f"{part}: {diff} > {bound}")
np.testing.assert_array_equal(o["force"], o["impulse"] / steps)
again = windowed.read_exchange_window(0)
again = windowed.read_exchange_window(0, level=1)
np.testing.assert_array_equal(again["total"]["impulse"], out["total"]["impulse"])
windowed.initialize(*arrays)
with self.assertRaises(RuntimeError):
windowed.read_exchange_window(0)
windowed.read_exchange_window(0, level=1)
def audit_substeps(stepper, b, count, record):
@@ -212,7 +212,7 @@ class ExchangeWindowShiftedMRTGPU(unittest.TestCase):
run.step(2)
if run is not plain:
run.start_exchange_window()
b = windowed._device_boundaries[1]
b = windowed._boundaries[1]
count = len(b.tables["cell"])
record = dict(scale=np.zeros(4), substeps=0)
audit_substeps(windowed, b, count, record)
@@ -220,7 +220,7 @@ class ExchangeWindowShiftedMRTGPU(unittest.TestCase):
for _ in range(steps):
windowed.step(1)
plain.step(1)
summed += components(windowed.read_raw_exchange(0))
summed += components(windowed.read_raw_exchange(0, level=1))
for a, e in zip(windowed._download_raw().values, plain._download_raw().values):
same_bits(self, a, e)
split.step(steps) # one call: same ledger arithmetic as five step(1) calls
@@ -228,31 +228,31 @@ class ExchangeWindowShiftedMRTGPU(unittest.TestCase):
same_bits(self, a, e)
rows = [np.empty((count, 4), np.float64) for _ in range(2)]
cuda.memcpy_dtoh(rows[0], b.window_integral)
cuda.memcpy_dtoh(rows[1], split._device_boundaries[1].window_integral)
cuda.memcpy_dtoh(rows[1], split._boundaries[1].window_integral)
np.testing.assert_array_equal(rows[0].view(np.uint64), rows[1].view(np.uint64))
out = windowed.read_exchange_window(0)
out = windowed.read_exchange_window(0, level=1)
self.assertEqual(record["substeps"], 2 * steps)
self.assertEqual((out["time_end"] - out["time_start"], out["sample_count"]),
(steps, 2 * steps))
bound = 4 * (count * record["substeps"]) * EPS * record["scale"] + 1e-300
diff = np.abs(components(out) - summed)
self.assertTrue((diff <= bound).all(), f"{diff} > {bound}")
np.testing.assert_array_equal(split.read_exchange_window(0)["total"]["impulse"],
np.testing.assert_array_equal(split.read_exchange_window(0, level=1)["total"]["impulse"],
out["total"]["impulse"])
windowed.step(0) # no interval: window and its duration unchanged
again = windowed.read_exchange_window(0)
again = windowed.read_exchange_window(0, level=1)
self.assertEqual(again["time_end"], out["time_end"])
np.testing.assert_array_equal(again["total"]["impulse"], out["total"]["impulse"])
windowed.start_exchange_window() # reset after a window
with self.assertRaisesRegex(RuntimeError, "at least one completed"):
windowed.read_exchange_window(0)
windowed.read_exchange_window(0, level=1)
record.update(scale=np.zeros(4), substeps=0)
windowed.step(1)
single = components(windowed.read_raw_exchange(0))
reset = windowed.read_exchange_window(0)
single = components(windowed.read_raw_exchange(0, level=1))
reset = windowed.read_exchange_window(0, level=1)
self.assertEqual(reset["time_end"] - reset["time_start"], 1)
bound = 4 * (count * record["substeps"]) * EPS * record["scale"] + 1e-300
self.assertTrue((np.abs(components(reset) - single) <= bound).all())
@@ -81,7 +81,7 @@ def run_blocked_patch(speed):
level_objects=[(), (source,)],
observe_raw_exchange=True) as backend:
c, f = backend._levels
flags = backend.fine_boundary.flags.reshape(40, 40)
flags = backend.boundaries[1].flags.reshape(40, 40)
active = np.zeros((40, 40), bool)
active[4:-4, 4:-4] = True
active &= ((flags & FLUID) != 0) & ((flags & SOLID) == 0)
@@ -101,7 +101,7 @@ def run_blocked_patch(speed):
with patch.object(f, "kernel", fine):
backend.step(2)
assert len(traces) == 4 and backend.step_count == 2
mass = backend.read_raw_exchange(0)["total"]["mass_into_fluid"]
mass = backend.read_raw_exchange(0, level=1)["total"]["mass_into_fluid"]
return traces, roi, backend.download(), float(mass)
@@ -56,7 +56,7 @@ def make(context):
def topology(run, rectangles, objects):
return RemapTopology.from_config(
CFG, run.patch, rectangles, child_objects=objects,
CFG, run.hierarchy.pairs[0], rectangles, child_objects=objects,
child_jets=run._boundaries[2].jets,
)
@@ -75,7 +75,7 @@ class MultiplePinnedRemapGPU(unittest.TestCase):
# Each jet selects a different body. Updating one must preserve the
# other body's current per-link velocities, including its driven links.
boundary = run._boundaries[2]
body_ids = run.child_boundary.curved_host[3]
body_ids = run.boundaries[2].curved_host[3]
for index, amplitude in enumerate(amplitudes):
before = boundary.tables['surface_velocity'].copy()
run.set_jet_amplitude(index, amplitude, level=2)
@@ -87,8 +87,8 @@ class MultiplePinnedRemapGPU(unittest.TestCase):
def check_boundary(self, old, new, mesh, amplitudes):
self.assertEqual(mesh.pinned_rectangles, PINS)
self.assertTrue(same_physical_boundary(old.child_boundary, new.child_boundary))
self.assertTrue(same_physical_boundary(new.child_boundary, mesh.boundary))
self.assertTrue(same_physical_boundary(old.boundaries[2], new.boundaries[2]))
self.assertTrue(same_physical_boundary(new.boundaries[2], mesh.boundary))
# Compare physical cut-link rows, never layout-dependent blocked addresses.
for run in (old, new):
boundary = run._boundaries[2]
@@ -107,7 +107,7 @@ class MultiplePinnedRemapGPU(unittest.TestCase):
before, count = old.download(), old.step_count
mesh = topology(old, START, objects)
# Input order is irrelevant even with several pins and fluid regions.
new = owned.enter_context(old.remap_child(START[::-1], child_objects=objects))
new = owned.enter_context(old.remap_finest(START[::-1], objects=objects))
self.assertIsNot(old, new)
self.assertEqual((old.step_count, new.step_count), (count, count))
self.check_boundary(old, new, mesh, amplitudes)
@@ -133,7 +133,7 @@ class MultiplePinnedRemapGPU(unittest.TestCase):
previous = current
self.change_amplitudes(previous, amplitudes)
before, count = previous.download(), previous.step_count
old_mesh = topology(previous, previous.child_patch.bounds, objects)
old_mesh = topology(previous, previous.hierarchy.pairs[1].bounds, objects)
new_mesh = topology(previous, rectangles, objects)
expected = remap_synchronized(
old_mesh, new_mesh, before, base_count=count,
@@ -141,7 +141,7 @@ class MultiplePinnedRemapGPU(unittest.TestCase):
phase='postcollision_g',
).g
current = owned.enter_context(
previous.remap_child(rectangles, child_objects=objects))
previous.remap_finest(rectangles, objects=objects))
self.assertIsNot(previous, current)
self.assertEqual((previous.step_count, current.step_count), (count, count))
self.check_boundary(previous, current, new_mesh, amplitudes)
@@ -36,8 +36,8 @@ class PinnedBodyRemapGPU(unittest.TestCase):
def check_boundary(self, before, after):
for name in ('flags','boundary_cells','curved_hits'):
exact(self,getattr(after.child_boundary,name),getattr(before.child_boundary,name),name)
for a,b in zip(after.child_boundary.curved_host,before.child_boundary.curved_host):
exact(self,getattr(after.boundaries[2],name),getattr(before.boundaries[2],name),name)
for a,b in zip(after.boundaries[2].curved_host,before.boundaries[2].curved_host):
exact(self,a,b,'physical cut links')
exact(self,after._boundaries[2].tables['surface_velocity'],
before._boundaries[2].tables['surface_velocity'],'current per-link velocities')
@@ -49,7 +49,7 @@ class PinnedBodyRemapGPU(unittest.TestCase):
old.initialize(*initial(old,False));old.step(2)
old.set_jet_amplitude(0,.017,level=2)
before=old.download()
new=owned.enter_context(old.remap_child(A,child_objects=objects))
new=owned.enter_context(old.remap_finest(A,objects=objects))
self.check_boundary(old,new);self.check_reset(new,before)
old.step();new.step()
for a,b in zip(new.download(),old.download()): exact(self,a,b,'no-op next step bitwise')
@@ -64,13 +64,13 @@ class PinnedBodyRemapGPU(unittest.TestCase):
for rects in (B,A):
previous=current;before=previous.download();count=previous.step_count
jets=previous._boundaries[2].jets
old=RemapTopology.from_config(CFG,previous.patch,previous.child_patch.bounds,
old=RemapTopology.from_config(CFG,previous.hierarchy.pairs[0],previous.hierarchy.pairs[1].bounds,
child_objects=objects,child_jets=jets)
new=RemapTopology.from_config(CFG,previous.patch,rects,
new=RemapTopology.from_config(CFG,previous.hierarchy.pairs[0],rects,
child_objects=objects,child_jets=jets)
expected=remap_synchronized(old,new,before,base_count=count,
level_steps=(count,2*count,4*count),phase='postcollision_g').g
current=owned.enter_context(previous.remap_child(rects,child_objects=objects))
current=owned.enter_context(previous.remap_finest(rects,objects=objects))
self.check_boundary(previous,current);self.check_reset(current,expected)
for a,b,om,nm in zip(current.download(),before,old.leaf_masks,new.leaf_masks):
exact(self,a[:,om & nm],b[:,om & nm],'physical-fluid overlap')
@@ -61,7 +61,7 @@ def census(run, index, buffer):
layout, b = level.layout, level.boundary
g = download(buffer, layout.blocked_shape)
dense = layout.to_dense(g)
leaf = run._leaf[index]
leaf = run.leaf[index]
ny, nx = leaf.shape
solid = (b.host.flags.reshape(ny, nx) & SOLID) != 0
y, x = np.nonzero(leaf)
@@ -105,7 +105,7 @@ class RawExchangeGPU(unittest.TestCase):
arrays = [state(l.shape) for l in measured._levels]
measured.initialize(*arrays)
plain.initialize(*arrays)
b = measured._device_boundaries[1]
b = measured._boundaries[1]
count = len(b.tables['cell'])
kernel = measured._levels[1].kernel
samples = []
@@ -119,7 +119,7 @@ class RawExchangeGPU(unittest.TestCase):
measured._stream.synchronize()
after = measured._levels[1].layout.to_dense(
download(dst, measured._levels[1].layout.blocked_shape))
change = totals(after, measured._leaf[1])-before
change = totals(after, measured.leaf[1])-before
expected = external + np.array([-raw[:, 0].astype(float).sum(),
-raw[:, 1].astype(float).sum(),
raw[:, 3].astype(float).sum()])
@@ -127,8 +127,8 @@ class RawExchangeGPU(unittest.TestCase):
samples.append(raw.copy())
measured._levels[1].kernel = audited_step
for speed in (0., .002, -.001, 0.):
measured.set_jet_amplitude(0, speed)
plain.set_jet_amplitude(0, speed)
measured.set_jet_amplitude(0, speed, level=1)
plain.set_jet_amplitude(0, speed, level=1)
samples.clear()
measured.step()
plain.step()
@@ -139,11 +139,11 @@ class RawExchangeGPU(unittest.TestCase):
np.testing.assert_array_equal(integral, expected)
last = download(b.last_raw, (count, 4))
np.testing.assert_array_equal(last, samples[-1])
obs = measured.read_raw_exchange(0)
obs = measured.read_raw_exchange(0, level=1)
assert obs['sample_count'] == 2
assert b.tables['aperture'].any()
measured.step(0)
np.testing.assert_array_equal(measured.read_raw_exchange(0)['total']['impulse'],
np.testing.assert_array_equal(measured.read_raw_exchange(0, level=1)['total']['impulse'],
obs['total']['impulse'])
def test_depth_four_interval_storage_and_sampling(self):
@@ -158,7 +158,7 @@ class RawExchangeGPU(unittest.TestCase):
FixedPatchStepper(cfg, context=ctx, **settings) as plain:
arrays = [state(l.shape, ux=0.) for l in run._levels]
run.initialize(*arrays); plain.initialize(*arrays)
b = run._device_boundaries[3]
b = run._boundaries[3]
count = len(b.tables['cell'])
kernel = run._levels[3].kernel
samples, independent_samples = [], []
@@ -171,7 +171,7 @@ class RawExchangeGPU(unittest.TestCase):
run._stream.synchronize()
after = run._levels[3].layout.to_dense(
download(dst, run._levels[3].layout.blocked_shape))
change = totals(after, run._leaf[3])-before
change = totals(after, run.leaf[3])-before
expected = external + np.array([-independent[:, 0].sum(),
-independent[:, 1].sum(), independent[:, 3].sum()])
np.testing.assert_allclose(change, expected, rtol=2e-4, atol=8e-5)
@@ -213,7 +213,7 @@ class RawExchangeGPU(unittest.TestCase):
FixedPatchStepper(CFG, context=ctx, **fine((16, 20)), observe_raw_exchange=True) as amr:
amr.initialize(*(state(l.shape) for l in amr._levels))
uniform.initialize(state((9, 48, 64)))
b = amr._device_boundaries[1]
b = amr._boundaries[1]
count = len(b.tables['cell'])
links = obj.geometry.build_cut_links(64, 48)
np.testing.assert_array_equal([v.dir for v in links], b.tables['dir'])
@@ -221,7 +221,7 @@ class RawExchangeGPU(unittest.TestCase):
# Same fine-grid force arms after translation to the same center: the
# single level's static BFL gets the AMR Jet column and force arms; the raw
# output pointer is substituted at launch (no stepper-side accumulation).
root = uniform._device_boundaries[0]
root = uniform._boundaries[0]
velocity = uniform._upload(b.tables['surface_velocity'])
root.device.update(surface_velocity=velocity, rx=uniform._upload(b.tables['rx']),
ry=uniform._upload(b.tables['ry']))
@@ -234,7 +234,7 @@ class RawExchangeGPU(unittest.TestCase):
samples.append(download(raw, (count, 4)))
root.kernel = measured_boundary
for speed in (0., .002, -.001, 0.):
amr.set_jet_amplitude(0, speed)
amr.set_jet_amplitude(0, speed, level=1)
cuda.memcpy_htod(velocity, b.tables['surface_velocity'])
samples.clear()
amr.step(); uniform.step(2)
@@ -205,7 +205,7 @@ class BlockedBaseCPU(unittest.TestCase):
level_blocks=[(3, 4), FINE_BLOCK],
level_objects=[(), (circle(actuation(.002)),)])
try:
backend.initialize(state(NY, NX, 1), state(*backend.patch.fine_shape, 2))
backend.initialize(state(NY, NX, 1), state(*backend.hierarchy.pairs[0].fine_shape, 2))
coarse, fine = backend._levels
current, spare, fine_current = coarse.current, coarse.spare, fine.current
backend.step(2)
@@ -256,7 +256,7 @@ class BlockedBaseCPU(unittest.TestCase):
with mocked_backend() as m:
backend = FixedPatchStepper(config("periodic"), context=m["context"], refinements=[BOUNDS],
level_blocks=[(3, 4), FINE_BLOCK])
backend.initialize(state(NY, NX, 1), state(*backend.patch.fine_shape, 2))
backend.initialize(state(NY, NX, 1), state(*backend.hierarchy.pairs[0].fine_shape, 2))
m["functions"]["BlockCoalesce"].side_effect = RuntimeError("injected parent stage failure")
current = backend._levels[0].current
with self.assertRaisesRegex(RuntimeError, "injected parent stage"):
@@ -126,7 +126,7 @@ class BlockedFineCPU(unittest.TestCase):
patch.object(cuda, "mem_alloc", side_effect=alloc), patch.object(cuda, "memcpy_htod"):
backend = FixedPatchStepper(config("periodic"), context=context, refinements=[BOUNDS],
level_blocks=[(12, 8), FINE_BLOCKS[0]])
self.assertEqual(len(backend.fine_layout.halo_dst), 0)
self.assertEqual(len(backend.layouts[1].halo_dst), 0)
self.assertIsNone(backend._levels[1].halo_dst)
backend._gather_halo(backend._levels[1], Mock()) # no launch without pairs
module.get_function.return_value.assert_not_called()
@@ -179,7 +179,7 @@ class BlockedFineCPU(unittest.TestCase):
backend = FixedPatchStepper(config("periodic"), context=context, refinements=[BOUNDS],
level_blocks=[(12, 8), FINE_BLOCKS[1]],
level_objects=[(), tuple(objects)])
backend.initialize(state(NY, NX, 1), state(*backend.patch.fine_shape, 2))
backend.initialize(state(NY, NX, 1), state(*backend.hierarchy.pairs[0].fine_shape, 2))
return backend, calls
def test_same_buffer_halo_bfl_step_launch_order(self):
@@ -187,8 +187,8 @@ class BlockedFineCPU(unittest.TestCase):
backend, calls = self.mocked_backend([circle(spec)])
c, f = backend._levels
coarse, cur, spare = c.current, f.current, f.spare
flags, kind = backend._device_boundaries[1].flags, backend._levels[1].kind
velocity = backend._device_boundaries[1].device.get("surface_velocity")
flags, kind = backend._boundaries[1].flags, backend._levels[1].kind
velocity = backend._boundaries[1].device.get("surface_velocity")
self.assertEqual(velocity is None, spec is None)
backend.step(1)
backend.close()
@@ -213,8 +213,8 @@ class BlockedFineCPU(unittest.TestCase):
def test_actuation_without_velocity_buffer_fails_before_bfl(self):
for speed in (0.0, .002): # any spec, even zero amplitude
backend, calls = self.mocked_backend([circle(actuation(speed))])
self.assertIsNotNone(backend._fine_tables["surface_velocity"])
del backend._device_boundaries[1].device["surface_velocity"]
self.assertIsNotNone(backend._boundaries[1].tables["surface_velocity"])
del backend._boundaries[1].device["surface_velocity"]
with self.assertRaisesRegex(RuntimeError, "Surface velocity device buffer is missing"):
backend.step(1)
names = [name for name, _ in calls]
@@ -304,7 +304,7 @@ class NormalBuildGPU(unittest.TestCase):
if not (x0 <= sx < x1 and y0 <= sy < y1):
corners.add(target)
self.assertEqual(len(corners), 8)
self.assertEqual(set(backend.patch.coalesce_coarse.tolist()), targets)
self.assertEqual(set(backend.hierarchy.pairs[0].coalesce_coarse.tolist()), targets)
backend.initialize(*initial)
before = discrete_mass(initial)
for step in range(1, 9):
@@ -383,8 +383,8 @@ class NormalBuildGPU(unittest.TestCase):
expected = scalar_interval(*[a.astype(np.float64) for a in initial], bounds)
with CudaContext(0) as context, FixedPatchStepper(cfg, context=context, refinements=[bounds],
level_blocks=[(13, 17), (14, 16)]) as backend:
self.assertEqual(backend.patch.fine_shape, (14, 16))
p = backend.patch
self.assertEqual(backend.hierarchy.pairs[0].fine_shape, (14, 16))
p = backend.hierarchy.pairs[0]
allocated = 2 * 4 * sum(int(np.prod(level.layout.blocked_shape)) for level in backend._levels) \
+ sum(a.nbytes for a in (p.ghost_fine, p.ghost_coarse, p.coalesce_coarse, p.coalesce_fine))
self.assertLess(allocated, 128*1024)
@@ -91,8 +91,8 @@ class CoupledGPU(unittest.TestCase):
ddf_shifting=shifted)
with FixedPatchStepper(cfg, context=context, refinements=[(4, 4, 12, 12)],
**storage) as backend:
self.assertAlmostEqual(backend.patch.tau_coarse, spec.tau_coarse)
self.assertAlmostEqual(backend.patch.tau_fine, spec.tau_fine)
self.assertAlmostEqual(backend.hierarchy.pairs[0].tau_coarse, spec.tau_coarse)
self.assertAlmostEqual(backend.hierarchy.pairs[0].tau_fine, spec.tau_fine)
self.assertEqual(tuple(level.shape for level in backend._levels),
((9, 16, 16), (9, 24, 24)))
# Physical FP64 g is encoded by initialize before FP32 rounding.
@@ -101,7 +101,7 @@ class CoupledGPU(unittest.TestCase):
backend.step(1)
actual = backend.download()
self.assertEqual(backend.step_count, 1)
backend.patch.check_synchronization(
backend.hierarchy.pairs[0].check_synchronization(
1, next_state.tick, phase="postcollision_g")
for got, want, arrival in zip(actual, expected, arrivals):
np.testing.assert_array_equal(np.isnan(got), np.isnan(want))
@@ -165,13 +165,13 @@ class CoupledGPU(unittest.TestCase):
with self.assertRaisesRegex(RuntimeError, "failed"):
operation()
raw = level_storage(fine)
self.assertTrue(np.isfinite(raw.ravel()[backend.patch.coalesce_fine]).all())
self.assertTrue(np.isfinite(raw.ravel()[backend.hierarchy.pairs[0].coalesce_fine]).all())
self.assertFalse(context.is_closed)
print("partial failure: two fine hops completed, no coarse step published; cleanup PASS", flush=True)
def eight_steps(self, backend, spec, maps):
"""Rest and seeded nonequilibrium, 8 steps against the scalar reference."""
p = backend.patch
p = backend.hierarchy.pairs[0]
for case in ("rest", "nonequilibrium_seed1"):
with self.subTest(case=case):
state = ref.make_state(spec, velocity=(0, 0) if case == "rest" else (.03, -.02))
@@ -62,8 +62,8 @@ class JetStepperCPU(unittest.TestCase):
C = (16.0, 12.0)
whole = Jet(Sector(C, .1, math.pi, 3.), Radial(C, inward=False), amplitude=0.)
backend, _ = mocked([circle()], [whole], FINE_BLOCKS[1])
_, jet_id = backend._jet_cache
self.assertGreater(len(np.unique(backend._fine_tables["block"][jet_id == 0])), 1)
_, jet_id = backend._boundaries[1].jet_cache
self.assertGreater(len(np.unique(backend._boundaries[1].tables["block"][jet_id == 0])), 1)
backend.close()
def test_fine_ownership_uses_remapped_physical_owner(self):
@@ -76,10 +76,10 @@ class JetStepperCPU(unittest.TestCase):
source = SimObject(99, geometry=CircleGeometry(*C, 2.15), center=C, radius=2.15)
whole = Jet(Sector(C, .1, math.pi, 3.), Radial(C, inward=False), amplitude=0.)
backend, _ = mocked([source], [whole], FINE_BLOCKS[1])
_, jet_id = backend._jet_cache
owners = backend.fine_boundary.curved_host[3][jet_id == 0]
_, jet_id = backend._boundaries[1].jet_cache
owners = backend.boundaries[1].curved_host[3][jet_id == 0]
self.assertEqual(np.unique(owners).tolist(), [0])
self.assertGreater(len(np.unique(backend._fine_tables["block"][jet_id == 0])), 1)
self.assertGreater(len(np.unique(backend._boundaries[1].tables["block"][jet_id == 0])), 1)
self.assertEqual(source.obj_id, 99)
backend.close()
cfg = config("periodic")
@@ -97,8 +97,8 @@ class JetStepperCPU(unittest.TestCase):
def test_blocked_jet_velocity_matches_legacy_rows(self):
legacy, _ = mocked([circle(actuation(SPEED))], (), FINE_BLOCKS[1])
driven, _ = mocked([circle()], [jet()], FINE_BLOCKS[1])
a = legacy._fine_tables["surface_velocity"]
b = driven._fine_tables["surface_velocity"]
a = legacy._boundaries[1].tables["surface_velocity"]
b = driven._boundaries[1].tables["surface_velocity"]
diff = np.abs(a - b).max()
print(f"fine jet vs legacy velocity rows: max |dU| = {diff:.3e} (exact={diff == 0})")
self.assertLessEqual(diff, 4 * np.finfo(np.float32).eps * SPEED)
@@ -106,20 +106,20 @@ class JetStepperCPU(unittest.TestCase):
def test_blocked_amplitude_update_rewrites_only_velocity_column(self):
backend, htod = mocked([circle()], [jet()], FINE_BLOCKS[1])
before = backend._fine_tables["surface_velocity"].copy()
device = backend._device_boundaries[1].device["surface_velocity"]
before = backend._boundaries[1].tables["surface_velocity"].copy()
device = backend._boundaries[1].device["surface_velocity"]
with patch.object(cuda, "memcpy_htod", htod):
htod.reset_mock()
backend.set_jet_amplitude(0, -SPEED)
backend.set_jet_amplitude(0, -SPEED, level=1)
(ptr, sent), _ = htod.call_args
self.assertEqual(htod.call_count, 1)
self.assertIs(ptr, device)
np.testing.assert_allclose(sent, -before, rtol=1e-6, atol=0)
np.testing.assert_array_equal(backend._fine_tables["surface_velocity"], sent)
np.testing.assert_array_equal(backend._boundaries[1].tables["surface_velocity"], sent)
with self.assertRaises(ValueError):
backend.set_jet_amplitude(0, 1.5)
backend.set_jet_amplitude(0, 1.5, level=1)
self.assertEqual(htod.call_count, 1)
self.assertEqual(backend._jets[0].amplitude, -SPEED)
self.assertEqual(backend.jets[1][0].amplitude, -SPEED)
backend.close()
def test_copy_failure_marks_failed_and_keeps_host_amplitude(self):
@@ -128,24 +128,24 @@ class JetStepperCPU(unittest.TestCase):
status = backend._status
for bad, error in ((1.5, ValueError), (float("nan"), ValueError)):
with self.assertRaises(error):
backend.set_jet_amplitude(0, bad)
backend.set_jet_amplitude(0, bad, level=1)
self.assertEqual(backend._status, status) # invalid input stays usable
with self.assertRaises(IndexError):
backend.set_jet_amplitude(1, -SPEED)
backend.set_jet_amplitude(1, -SPEED, level=1)
self.assertEqual(backend._status, status)
before = backend._fine_tables["surface_velocity"].copy()
before = backend._boundaries[1].tables["surface_velocity"].copy()
boom = Mock(side_effect=RuntimeError("copy failed"))
with patch.object(cuda, "memcpy_htod", boom), patch.object(cuda, "memcpy_htod_async", boom), \
patch.object(cuda, "pagelocked_empty", side_effect=lambda shape, dtype: np.empty(shape, dtype)):
with self.assertRaisesRegex(RuntimeError, "copy failed"):
backend.set_jet_amplitude(0, -SPEED)
backend.set_jet_amplitude(0, -SPEED, level=1)
self.assertEqual(backend._status, "failed", fine_block)
self.assertEqual(backend._jets[0].amplitude, SPEED)
np.testing.assert_array_equal(backend._fine_tables["surface_velocity"], before)
self.assertEqual(backend.jets[1][0].amplitude, SPEED)
np.testing.assert_array_equal(backend._boundaries[1].tables["surface_velocity"], before)
with self.assertRaisesRegex(RuntimeError, "failed"):
backend.step(1)
with self.assertRaisesRegex(RuntimeError, "failed"):
backend.set_jet_amplitude(0, -SPEED)
backend.set_jet_amplitude(0, -SPEED, level=1)
backend.close()
def test_jet_overlapping_legacy_aperture_rejected_before_cuda(self):
@@ -170,7 +170,7 @@ def run(mode, fine_block, *, objects, jets, switch=None):
backend.step(STEPS)
else:
backend.step(STEPS // 2)
backend.set_jet_amplitude(0, switch)
backend.set_jet_amplitude(0, switch, level=1)
backend.step(STEPS - STEPS // 2)
return backend.download(), raw(backend)
@@ -36,14 +36,14 @@ 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)
for p, w in zip(patch_set.patches, windows):
y0, x0 = p.fine_level_origin
y0, x0 = p.fine_level_offset
out[:, y0:y0 + p.fine_shape[0], x0:x0 + p.fine_shape[1]] = w
return out
def window(patch_set, index, fine):
p = patch_set.patches[index]
y0, x0 = p.fine_level_origin
y0, x0 = p.fine_level_offset
return fine[:, y0:y0 + p.fine_shape[0], x0:x0 + p.fine_shape[1]]
@@ -175,7 +175,7 @@ class MultipatchBodiesGPU(unittest.TestCase):
for step in self.STEPS:
if switch is not None and step == 8:
for index, amplitude in switch:
backend.set_jet_amplitude(index, amplitude)
backend.set_jet_amplitude(index, amplitude, level=1)
backend.step(step - backend.step_count)
out.append(backend.download() + raw(backend))
return out
@@ -128,8 +128,10 @@ class ThreeLevelCPU(unittest.TestCase):
run = mocked_child()
b = run._levels[2].boundary
before = b.tables["surface_velocity"].copy()
with self.assertRaises(TypeError):
run.set_jet_amplitude(0, -.002) # level is required
with self.assertRaises(IndexError):
run.set_jet_amplitude(0, -.002) # default stays L1
run.set_jet_amplitude(0, -.002, level=1) # level 1 has no Jets
for level in (0, 3, True):
with self.assertRaises(ValueError):
run.set_jet_amplitude(0, -.002, level=level)
@@ -138,7 +140,7 @@ class ThreeLevelCPU(unittest.TestCase):
self.assertEqual(upload.call_count, 1)
self.assertIs(upload.call_args.args[0], b.device["surface_velocity"])
np.testing.assert_array_equal(b.tables["surface_velocity"], -before)
self.assertEqual(run._jets, ())
self.assertEqual(run.jets[1], ())
with patch.object(cuda, "memcpy_htod", side_effect=RuntimeError("copy failed")):
with self.assertRaisesRegex(RuntimeError, "copy failed"):
run.set_jet_amplitude(0, .002, level=2)
@@ -150,17 +152,17 @@ class ThreeLevelCPU(unittest.TestCase):
def test_child_public_masks_include_solids(self):
from CelerisLab.lbm.descriptors import SOLID
run = mocked_child()
solid = run.child_boundary.flags.reshape(run.child_patch.fine_shape) & SOLID != 0
solid = run.boundaries[2].flags.reshape(run.hierarchy.pairs[1].fine_shape) & SOLID != 0
self.assertTrue(solid.any())
self.assertFalse(run._child_leaf[solid].any())
self.assertFalse(run.leaf[2][solid].any())
arrays = [np.ones(level.shape, np.float32) for level in run._levels]
run._mask_nonleaves(*arrays)
self.assertTrue(np.isnan(arrays[2][:, solid]).all())
with patch.object(cuda, "memcpy_htod") as upload:
run.initialize(*arrays)
sent = upload.call_args_list[-1].args[1]
dense = run.child_layout.to_dense(sent)
self.assertTrue(np.isnan(dense[:, ~run._child_leaf]).all())
dense = run.layouts[2].to_dense(sent)
self.assertTrue(np.isnan(dense[:, ~run.leaf[2]]).all())
run.close()
@@ -169,7 +171,7 @@ class ThreeLevelGPU(unittest.TestCase):
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:
shapes = ((16, 16), run.patch.fine_shape, run.child_patch.fine_shape)
shapes = ((16, 16), run.hierarchy.pairs[0].fine_shape, run.hierarchy.pairs[1].fine_shape)
run.initialize(*(np.broadcast_to(feq, (9, *s)).copy() for s in shapes))
run.step(8)
self.assertEqual(run.step_count, 8)
@@ -182,7 +184,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:
p1, p2 = run.patch, run.child_patch
p1, p2 = run.hierarchy.pairs[0], run.hierarchy.pairs[1]
run.initialize(state(16, 16, 1), state(*p1.fine_shape, 2), state(*p2.fine_shape, 3))
f = run._levels[1]
feeds, original = [], run._advance
@@ -195,7 +197,7 @@ class ThreeLevelGPU(unittest.TestCase):
run.step(1)
self.assertEqual(feeds, [cur, spare]) # current G, then the new G of hop 1
run.step(7)
middle = run._middle_kind
middle = run.hierarchy.kinds[1]
for g, leaf in zip(run.download(), (~p1.covered, middle == ACTIVE, p2.fine_kind == ACTIVE)):
self.assertTrue(np.isfinite(g[:, leaf]).all())
self.assertTrue(np.isnan(g[:, ~leaf]).all())
@@ -227,15 +229,15 @@ class ThreeLevelBodyGPU(unittest.TestCase):
if stop == 4 and switch:
run.set_jet_amplitude(0, -.002, level=2)
public, stored = run.download(), raw(run)
leaves = [~run.patch.covered, run._fine_leaf, run._child_leaf]
leaves = [~run.hierarchy.pairs[0].covered, run.leaf[1], run.leaf[2]]
if mode != "periodic":
leaves[0][[0, -1]] = False
for g, mask in zip(public, leaves):
self.assertTrue(np.isfinite(g[:, mask]).all())
self.assertTrue(np.isnan(g[:, ~mask]).all())
invalid = [run.patch.covered.copy(), run._middle_kind == OUTER_GHOST,
(run.child_patch.fine_kind == OUTER_GHOST) |
((run.child_boundary.flags.reshape(run.child_patch.fine_shape) & SOLID) != 0)]
invalid = [run.hierarchy.pairs[0].covered.copy(), run.hierarchy.kinds[1] == OUTER_GHOST,
(run.hierarchy.pairs[1].fine_kind == OUTER_GHOST) |
((run.boundaries[2].flags.reshape(run.hierarchy.pairs[1].fine_shape) & SOLID) != 0)]
if mode != "periodic":
invalid[0][[0, -1]] = True
for g, mask in zip(stored, invalid):
@@ -54,9 +54,9 @@ def jets(speed):
def frames(run):
middle_origin = np.asarray(run.patch.fine_origin)
middle_origin = np.asarray(run.hierarchy.pairs[0].fine_origin)
return ((np.zeros(2), 1.), (middle_origin, .5),
(middle_origin + .5*np.asarray(run.child_patch.fine_origin), .25))
(middle_origin + .5*np.asarray(run.hierarchy.pairs[1].fine_origin), .25))
def physical_state(shape, origin, spacing):
@@ -92,15 +92,15 @@ class ThreeLevelMultipatchGPU(unittest.TestCase):
for index in range(len(drive)):
run.set_jet_amplitude(index, -.002, level=2)
public, stored = run.download(), raw(run)
masks = [~run.patch.covered, run._fine_leaf, run._child_leaf]
masks = [~run.hierarchy.pairs[0].covered, run.leaf[1], run.leaf[2]]
if mode != 'periodic':
masks[0][[0, -1]] = False
invalid = [run.patch.covered.copy(), run._middle_kind == OUTER_GHOST,
run.child_patch.fine_kind == OUTER_GHOST]
invalid = [run.hierarchy.pairs[0].covered.copy(), run.hierarchy.kinds[1] == OUTER_GHOST,
run.hierarchy.pairs[1].fine_kind == OUTER_GHOST]
if mode != 'periodic':
invalid[0][[0, -1]] = True
if run.child_boundary is not None:
invalid[2] |= (run.child_boundary.flags.reshape(run.child_patch.fine_shape) & SOLID) != 0
if run.boundaries[2] is not None:
invalid[2] |= (run.boundaries[2].flags.reshape(run.hierarchy.pairs[1].fine_shape) & SOLID) != 0
for level, (g, mask, dense, bad) in enumerate(zip(public, masks, stored, invalid)):
label = f'{mode} stop={stop} level={level}'
self.assertTrue(np.isfinite(g[:, mask]).all(), label)
@@ -39,7 +39,7 @@ def initial(run, uniform):
class RuntimeChildRemapGPU(unittest.TestCase):
def check_public(self, run):
for array, leaf in zip(run.download(), (~run.patch.covered, run._fine_leaf, run._child_leaf)):
for array, leaf in zip(run.download(), (~run.hierarchy.pairs[0].covered, run.leaf[1], run.leaf[2])):
self.assertTrue(np.isfinite(array[:, leaf]).all())
self.assertTrue(np.isnan(array[:, ~leaf]).all())
@@ -57,7 +57,7 @@ class RuntimeChildRemapGPU(unittest.TestCase):
old.initialize(*initial(old, False))
old.step(2)
before = old.download()
new = owned.enter_context(old.remap_child(A))
new = owned.enter_context(old.remap_finest(A))
self.assertEqual((old.step_count, new.step_count), (2, 2))
self.check_reset(new, before)
for got, want in zip(old.download(), before): exact(self, got, want, 'old unchanged')
@@ -75,12 +75,12 @@ class RuntimeChildRemapGPU(unittest.TestCase):
current.step(2)
for rectangles in (B, A):
before, count = current.download(), current.step_count
old_mesh = RemapTopology.from_config(CFG, current.patch, current.child_patch.bounds)
new_mesh = RemapTopology.from_config(CFG, current.patch, rectangles)
old_mesh = RemapTopology.from_config(CFG, current.hierarchy.pairs[0], current.hierarchy.pairs[1].bounds)
new_mesh = RemapTopology.from_config(CFG, current.hierarchy.pairs[0], rectangles)
expected = remap_synchronized(old_mesh, new_mesh, before, base_count=count,
level_steps=(count, 2*count, 4*count), phase='postcollision_g').g
previous = current
current = owned.enter_context(previous.remap_child(rectangles))
current = owned.enter_context(previous.remap_finest(rectangles))
self.assertEqual(current.step_count, count)
self.check_reset(current, expected)
for got, want, old_g, old_leaf, new_leaf in zip(
@@ -8,14 +8,17 @@ from contextlib import ExitStack
import unittest
import numpy as np
from CelerisLab.cuda.context import CudaContext
from research.shear_refinement.example import (CHILDREN, make, initialize, observe,
topology, objects)
from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
from research.shear_refinement.example import (CFG, PARENTS, CHILDREN, initialize,
observe, topology, objects)
class ShearRefinementGPU(unittest.TestCase):
def test_measured_growth_matches_explicit_manual_continuation(self):
with CudaContext(0) as ctx, ExitStack() as owned:
old = owned.enter_context(make(ctx))
old = owned.enter_context(FixedPatchStepper(
CFG, context=ctx, refinements=[PARENTS, CHILDREN],
level_blocks=[(4, 4)] * 3, level_objects=[(), (), objects()]))
initialize(old)
old.step(2)
before = old.download()
@@ -25,10 +28,10 @@ class ShearRefinementGPU(unittest.TestCase):
self.assertGreater(report['requested'], 0)
self.assertEqual(report['corridor_unknown'], 0)
self.assertEqual(report['corridor_unavailable'], 0)
auto = owned.enter_context(old.remap_child(chosen, child_objects=objects()))
explicit = owned.enter_context(old.remap_child(manual, child_objects=objects()))
auto = owned.enter_context(old.remap_finest(chosen, objects=objects()))
explicit = owned.enter_context(old.remap_finest(manual, objects=objects()))
for run in (auto, explicit):
self.assertEqual(run.child_patch.bounds[0], CHILDREN[0])
self.assertEqual(run.hierarchy.pairs[1].bounds[0], CHILDREN[0])
self.assertEqual(run.step_count, 2)
run.step(2)
for a, b, leaf in zip(auto.download(), explicit.download(), topology(manual).leaf_masks):
@@ -58,7 +58,7 @@ def census(run,index,buffer):
"""Independent physical link/face census of actual staged receiver windows."""
level=run._levels[index];layout=level.layout;b=level.boundary
stored=download(buffer,layout.blocked_shape);dense=layout.to_dense(stored)
leaf=run._leaf[index];ny,nx=leaf.shape
leaf=run.leaf[index];ny,nx=leaf.shape
solid=(b.host.flags.reshape(ny,nx)&SOLID)!=0
y,x=np.nonzero(leaf);owner=layout.owner[y,x].astype(np.intp)
cells=int(np.prod(layout.local_shape));block,cell=np.divmod(owner,cells)
@@ -93,7 +93,7 @@ class Ledger:
self.max_row=max(self.max_row,float(abs(raw.astype(float)-independent).max()))
self.original(src,dst,*args,**kwargs);run._stream.synchronize()
after=level.layout.to_dense(download(dst,level.layout.blocked_shape))
change=raw_totals(after,run._leaf[index])-before
change=raw_totals(after,run.leaf[index])-before
expected=external+np.array([-independent[:,0].sum(),-independent[:,1].sum(),independent[:,3].sum()])
np.testing.assert_allclose(change,expected,rtol=2e-4,atol=8e-5)
self.max_balance=max(self.max_balance,float(abs(change-expected).max()))
@@ -102,7 +102,7 @@ class Ledger:
def step(self):
run=self.run;self.samples.clear();run.step()
b=run._device_boundaries[self.index];h=run._exchange_spacing[self.index]
b=run._boundaries[self.index];h=b.spacing
self.test.assertEqual(len(self.samples),2**self.index)
expected=np.zeros_like(self.samples[0],dtype=np.float64)
weights=np.array([h*h,h*h,h*h*h,h*h])
@@ -113,7 +113,7 @@ class Ledger:
for name in ('impulse','angular_impulse','mass_into_fluid'):
np.testing.assert_allclose(obs['total'][name],np.asarray(obs['wall'][name])+obs['aperture'][name],rtol=0,atol=2e-13)
self.test.assertEqual(obs['sample_count'],2**self.index)
for a,leaf in zip(run._download_raw().values,run._leaf):
for a,leaf in zip(run._download_raw().values,run.leaf):
self.test.assertTrue(np.isfinite(a[:,leaf]).all())
self.test.assertTrue((1+a[:,leaf].sum(axis=0,dtype=np.float64)>0).all())
return obs
@@ -140,10 +140,10 @@ class ShiftedJetsGPU(unittest.TestCase):
level_objects=[(),(obj,)],level_jets=[(),(whole,)],observe_raw_exchange=True) as run:
run._initialize_raw(initial(run,rho,u=.013,v=0.))
expected=run._download_raw();ledger=Ledger(self,run,1)
self.assertTrue(run._device_boundaries[1].tables['aperture'].all())
self.assertTrue(run._boundaries[1].tables['aperture'].all())
for _ in range(4):ledger.step()
error=max(float(abs(a[:,leaf]-b[:,leaf]).max()) for a,b,leaf in
zip(run._download_raw().values,expected.values,run._leaf))
zip(run._download_raw().values,expected.values,run.leaf))
self.assertLessEqual(error,1e-6) # inherited uniform-state cap
print(f'matched whole-circle rho={rho}: max_abs_raw_state={error:.9g}',flush=True)
ledger.report(f'matched whole-circle rho={rho}')
@@ -161,7 +161,7 @@ class ShiftedJetsGPU(unittest.TestCase):
ledger=Ledger(self,measured,1)
for speed in SPEEDS:
before=measured._download_raw();count=measured.step_count
measured.set_jet_amplitude(0,speed);plain.set_jet_amplitude(0,speed)
measured.set_jet_amplitude(0,speed, level=1);plain.set_jet_amplitude(0,speed, level=1)
self.exact_raw(before,measured._download_raw());self.assertEqual(measured.step_count,count)
ledger.step();plain.step()
self.exact_raw(measured._download_raw(),plain._download_raw())
@@ -193,15 +193,15 @@ class ShiftedJetsGPU(unittest.TestCase):
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
topology=lambda rects:RemapTopology.from_config(cfg,old.patch,rects,child_objects=(obj,),
child_jets=old._level_jets[2],_shifted_blocked=True)
topology=lambda rects:RemapTopology.from_config(cfg,old.hierarchy.pairs[0],rects,child_objects=(obj,),
child_jets=old.jets[2],_shifted_blocked=True)
expected=_remap_stored(topology(a),topology(b),before,shifted=True,base_count=count,
level_steps=(count,2*count,4*count),phase='postcollision_g')
new=owned.enter_context(old.remap_child(b,child_objects=(obj,)))
new=owned.enter_context(old.remap_finest(b,objects=(obj,)))
self.exact_raw(before,old._download_raw())
self.exact_raw(_StoredPopulations(expected.g,True),new._download_raw())
self.assertEqual(new._level_jets[2][0].amplitude,-.001)
np.testing.assert_array_equal(new._device_boundaries[2].tables['surface_velocity'],old._device_boundaries[2].tables['surface_velocity'])
self.assertEqual(new.jets[2][0].amplitude,-.001)
np.testing.assert_array_equal(new._boundaries[2].tables['surface_velocity'],old._boundaries[2].tables['surface_velocity'])
newledger=Ledger(self,new,2);ledger.step();newledger.step()
self.assertEqual((old.step_count,new.step_count),(count+1,count+1))
ledger.report('L2 source');newledger.report('L2 remapped')
@@ -61,7 +61,7 @@ class ShiftedRemapGPU(unittest.TestCase):
def check_finite(self,run):
raw=run._download_raw()
self.assertTrue(raw.shifted)
for a,leaf in zip(raw.values,run._leaf):
for a,leaf in zip(raw.values,run.leaf):
self.assertTrue(np.isfinite(a[:,leaf]).all())
self.assertTrue(np.isnan(a[:,~leaf]).all())
self.assertTrue((1+a[:,leaf].sum(axis=0,dtype=np.float64)>0).all())
@@ -71,7 +71,7 @@ class ShiftedRemapGPU(unittest.TestCase):
old=owned.enter_context(make(ctx))
old._initialize_raw(initial(old));old.step(2)
before=old._download_raw()
new=owned.enter_context(old.remap_child(A))
new=owned.enter_context(old.remap_finest(A))
self.check_reset(new,before.values)
self.snapshot_equal(old._download_raw(),before,'old unchanged')
self.snapshot_equal(new._download_raw(),before,'identity raw copy')
@@ -85,11 +85,11 @@ class ShiftedRemapGPU(unittest.TestCase):
old=owned.enter_context(make(ctx))
old._initialize_raw(initial(old));old.step(2)
before=old._download_raw();count=old.step_count
old_mesh=RemapTopology.from_config(CFG,old.patch,A)
new_mesh=RemapTopology.from_config(CFG,old.patch,B)
old_mesh=RemapTopology.from_config(CFG,old.hierarchy.pairs[0],A)
new_mesh=RemapTopology.from_config(CFG,old.hierarchy.pairs[0],B)
expected=_remap_stored(old_mesh,new_mesh,before,shifted=True,
base_count=count,level_steps=(count,2*count,4*count),phase='postcollision_g')
new=owned.enter_context(old.remap_child(B))
new=owned.enter_context(old.remap_finest(B))
self.check_reset(new,expected.g)
self.snapshot_equal(new._download_raw(),_StoredPopulations(expected.g,True),'host raw fill')
for a,b,om,nm in zip(before.values,new._download_raw().values,old_mesh.leaf_masks,new_mesh.leaf_masks):
@@ -103,7 +103,7 @@ class ShiftedRemapGPU(unittest.TestCase):
with CudaContext(0) as ctx,make(ctx) as run:
run._initialize_raw(initial(run));run.step(2)
before=run._download_raw()
with self.assertRaises(ValueError):run.remap_child([])
with self.assertRaises(ValueError):run.remap_finest([])
self.snapshot_equal(run._download_raw(),before,'rejection preserves source')
self.assertEqual(run.step_count,2)
run.step();self.check_finite(run)
+8 -4
View File
@@ -270,10 +270,14 @@ def host_geometry(source, h, collision="SRT", jet_scenario="none", placement="wi
pair = FixedPatch2D.from_config(cfg, PLACEMENTS[placement])
layout = BlockLayout2D.build(cfg.ny, cfg.nx, *BASE_BLOCK)
validate_channel_config(cfg, [BASE_BLOCK, FINE_BLOCK])
# Historical pinned source (jet_scenario none) names it blocked_base_tables;
# the current checkout (public Jet scenarios) names it blocked_coarse_tables.
coarse_tables = (refinement_tables.blocked_base_tables if jet_scenario == "none"
else refinement_tables.blocked_coarse_tables)
# Select by the imported source's own API, not by scenario: the pinned historical
# source (verify_source) defines blocked_base_tables, the current checkout (host
# tests, public Jet scenarios) blocked_coarse_tables. Exactly one must exist.
names = [n for n in ("blocked_base_tables", "blocked_coarse_tables")
if hasattr(refinement_tables, n)]
if len(names) != 1:
raise RuntimeError(f"Unexpected coarse transfer table API: {names}")
coarse_tables = getattr(refinement_tables, names[0])
targets = coarse_tables(pair, layout, np.argsort(pair.coalesce_coarse))[1]
validate_channel_transfers(pair, layout, targets)
BlockLayout2D.build(*pair.fine_shape, *FINE_BLOCK, periodic=(False, False))
+1 -1
View File
@@ -73,7 +73,7 @@ def run_variant(backend, variant, n_intervals, deadline):
else:
for _ in range(INTERVAL):
backend.step(1)
totals.append(components(backend.read_exchange_window(0)))
totals.append(components(backend.read_exchange_window(0, level=1)))
t.check_deadline(deadline)
return np.array(totals), (np.array(scales) if scales else None)