Compare commits
7
Commits
| Author | SHA1 | Date | |
|---|---|---|---|
|
|
6d046f91c1 | ||
|
|
3361e97d5f | ||
|
|
442979b4db | ||
|
|
c5a3cdb757 | ||
|
|
a4fc46d849 | ||
|
|
65339ab54e | ||
|
|
ae766f7112 |
No files matched your search
@@ -1,8 +1,11 @@
|
|||||||
# src coherence stage (from main 914906e)
|
# src coherence stage (from main 914906e)
|
||||||
|
|
||||||
Owner direction 2026-10-28: unified style over historical compatibility, inside `src`.
|
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
|
No numerical change is intended. Stages after 5f69eba: 65339ab (one boundary/Jet record,
|
||||||
pending a GPU run).
|
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`)
|
## Conventions (defined once, in `lbm.hierarchy.Hierarchy`)
|
||||||
|
|
||||||
@@ -28,20 +31,25 @@ pending a GPU run).
|
|||||||
|
|
||||||
## FixedPatchStepper migration map
|
## FixedPatchStepper migration map
|
||||||
|
|
||||||
|
Historical: left column = names at or after 914906e that no longer exist; right column =
|
||||||
|
the current name.
|
||||||
|
|
||||||
| removed / changed | replacement |
|
| removed / changed | replacement |
|
||||||
|---|---|
|
|---|---|
|
||||||
| `patch`, `child_patch` | `hierarchy.pairs[0]`, `hierarchy.pairs[1]` |
|
| `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) |
|
| `base_layout`, `fine_layout`, `child_layout` (stepper properties) | `layouts[0]`, `layouts[1]`, `layouts[2]` (`FixedPatchSet2D.fine_layout(block)` is unchanged) |
|
||||||
| `fine_boundary`, `child_boundary` | `boundaries[1]`, `boundaries[2]` (host geometry; `None` = no bodies/Jets) |
|
| `fine_boundary`, `child_boundary` | `boundaries[1]`, `boundaries[2]` (host geometry; `None` = no bodies/Jets) |
|
||||||
| `_jets` | `jets[1]` |
|
| `_jets` | `jets[1]` |
|
||||||
| `_boundaries[l]` (dict from L1, level keys) | `_device_boundaries[l]` (list from L0, `None` at L0 unless single-level static circles) |
|
| `_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 |
|
||||||
| `_level_boundaries` (private, b4eacf7) | `_device_boundaries` |
|
| `_fine_tables`, `_actuated`, `_jet_cache` | `_boundaries[1].tables` / `.actuated` / `.jet_cache` |
|
||||||
| `_fine_tables`, `_actuated`, `_jet_cache` | `_device_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` |
|
| `_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) |
|
| `refinement_tables.blocked_base_tables` (also via `refinement_stepper`) | `refinement_tables.blocked_coarse_tables` (import from `refinement_tables`; the stepper no longer re-exports) |
|
||||||
| `_fine_leaf`, `_child_leaf` | `leaf[1]`, `leaf[2]` (copies) |
|
| `_fine_leaf`, `_child_leaf` | `leaf[1]`, `leaf[2]` (copies) |
|
||||||
| `_middle_kind` | `hierarchy.kinds[1]` |
|
| `_middle_kind` | `hierarchy.kinds[1]` |
|
||||||
| `_boundary_hosts` (private) | `_boundaries` (private list behind `boundaries`) |
|
| `_boundary_hosts` (private) | 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=...)` |
|
| `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) |
|
| `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) |
|
| `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
|
## Responsibilities
|
||||||
|
|
||||||
- Construction: `_check_levels` (inputs), `_kernel_configs` (per-level kernel configs),
|
- Construction: `_check_levels` (inputs); host preparation, numpy only, before any device
|
||||||
`Hierarchy.edge_frame` (geometry frames), `_build_level` / `_build_pair` / `_upload_boundary`
|
work: `_level_layouts`, `_kernel_configs` (per-level kernel configs), `_host_boundaries`,
|
||||||
(device). Observation reads share `_observed_boundary` + `_summarize`.
|
`_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
|
- 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
|
point of `BlockStep` and `BlockChannelStep`. It computes rho/u once and reads no runtime
|
||||||
arithmetic variants (level: unclamped, no-force, shifted direct-momentum feq; channel: uniform
|
force (`d_params.fx/fy`). Two compile-time families: uniform-equivalent (MRT on any level,
|
||||||
OneStep arithmetic with omega clamp), selected by `BLOCK_ROOT_BC`.
|
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.
|
- `explosion="projected_linear"` stays an optional capability; `"constant"` stays default.
|
||||||
|
|
||||||
## Not done / open
|
## 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;
|
`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.
|
generates each level's config.h). Folding them needs a config.h/PTX diff gate.
|
||||||
- No single-`Simulation` facade over FixedPatchStepper: the two own different state models
|
- No single-`Simulation` facade over FixedPatchStepper: the two own different state models
|
||||||
(uniform node frame, body sync, forcing, checkpoint vs. immutable blocked hierarchy, edge
|
(uniform node frame, body sync, forcing, checkpoint vs. fixed blocked hierarchy, edge
|
||||||
frame, remap-as-new-instance). A facade would only forward calls; no concrete need found.
|
frame, remap-as-new-instance). A facade would only forward calls; no concrete need found.
|
||||||
|
- 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).
|
- Tests, examples and validation callers are migrated separately (pinned runners stay frozen).
|
||||||
@@ -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
|
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.
|
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 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.
|
the fixed refinement patch (refinement.py), whose semantics are unchanged.
|
||||||
|
|
||||||
|
|||||||
@@ -122,7 +122,7 @@ def leaf_partition(h):
|
|||||||
for level, kind in enumerate(h.kinds):
|
for level, kind in enumerate(h.kinds):
|
||||||
if level:
|
if level:
|
||||||
pair = h.pairs[level - 1]
|
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
|
oy, ox = 2 * oy + py, 2 * ox + px
|
||||||
y, x = np.nonzero(kind == ACTIVE)
|
y, x = np.nonzero(kind == ACTIVE)
|
||||||
np.add.at(base, ((y + oy) >> level, (x + ox) >> level), 0.25 ** level)
|
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).
|
// The one collision extension point of the blocked steps (BlockStep, BlockChannelStep).
|
||||||
// Compile-time selection: MRT shares one path; SRT keeps two explicit arithmetic
|
// Unforced by construction: FixedPatchStepper supplies no forcing, so no runtime
|
||||||
// variants, unchanged from their original kernels:
|
// d_params.fx/fy branch is read. Two compile-time arithmetic families:
|
||||||
// level (BLOCK_ROOT_BC 0/1): unclamped omega, no forcing term; shifted storage
|
// uniform-equivalent (MRT on every level; SRT on the channel root, BLOCK_ROOT_BC 2):
|
||||||
// builds feq from the retained drho and direct momentum (numerical variant);
|
// collide_dispatch arithmetic with zero force and no LES: clamped omega;
|
||||||
// channel (BLOCK_ROOT_BC 2): matches uniform OneStep SRT arithmetic (omega clamp,
|
// SRT keeps the zero Guo forcing term of collide_srt;
|
||||||
// zero Guo forcing term, compute_feq).
|
// 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
|
#ifndef CELERIS_OPERATORS_BLOCK_COLLISION_CUH
|
||||||
#define CELERIS_OPERATORS_BLOCK_COLLISION_CUH
|
#define CELERIS_OPERATORS_BLOCK_COLLISION_CUH
|
||||||
|
|
||||||
@@ -16,17 +18,19 @@
|
|||||||
#endif
|
#endif
|
||||||
|
|
||||||
__device__ __forceinline__ void block_collide(float* f) {
|
__device__ __forceinline__ void block_collide(float* f) {
|
||||||
#if COLLISION_MODEL == 2
|
float rho, ux, uy;
|
||||||
collide_mrt_unforced(f);
|
|
||||||
#elif BLOCK_ROOT_BC == 2
|
|
||||||
float rho, ux, uy, feq[NQ], forcing[NQ] = {};
|
|
||||||
compute_rho_u(f, 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));
|
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);
|
collide_srt(f, feq, forcing, omega);
|
||||||
|
#endif
|
||||||
#else
|
#else
|
||||||
float rho, ux, uy, feq[NQ];
|
float feq[NQ];
|
||||||
compute_rho_u(f, rho, ux, uy);
|
|
||||||
#if USE_DDF_SHIFTING
|
#if USE_DDF_SHIFTING
|
||||||
// Preserve the density perturbation before adding the unit reference.
|
// 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];
|
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);
|
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
|
#elif LATTICE_MODEL == LATTICE_D3Q19
|
||||||
|
|
||||||
__device__ __forceinline__ void collide_mrt(float* __restrict__ g,
|
__device__ __forceinline__ void collide_mrt(float* __restrict__ g,
|
||||||
|
|||||||
@@ -120,6 +120,11 @@ class FixedPatch2D:
|
|||||||
(an index offset, not an (x, y) origin)."""
|
(an index offset, not an (x, y) origin)."""
|
||||||
return 2 * self.bounds[1] - HALO, 2 * self.bounds[0] - HALO
|
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
|
@property
|
||||||
def covered(self) -> np.ndarray:
|
def covered(self) -> np.ndarray:
|
||||||
"""Dense coarse mask of cells replaced by the patch (not leaves)."""
|
"""Dense coarse mask of cells replaced by the patch (not leaves)."""
|
||||||
@@ -225,11 +230,18 @@ class FixedPatchSet2D:
|
|||||||
def bounds(self) -> tuple:
|
def bounds(self) -> tuple:
|
||||||
return tuple(p.bounds for p in self.patches)
|
return tuple(p.bounds for p in self.patches)
|
||||||
|
|
||||||
|
rectangles = bounds
|
||||||
|
|
||||||
@property
|
@property
|
||||||
def fine_origin(self) -> tuple[float, float]:
|
def fine_origin(self) -> tuple[float, float]:
|
||||||
"""Lower edge of the level-wide fine grid in coarse physical units."""
|
"""Lower edge of the level-wide fine grid in coarse physical units."""
|
||||||
return 0.0, 0.0
|
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
|
@property
|
||||||
def covered(self) -> np.ndarray:
|
def covered(self) -> np.ndarray:
|
||||||
mask = np.zeros(self.coarse_shape, bool)
|
mask = np.zeros(self.coarse_shape, bool)
|
||||||
|
|||||||
@@ -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")
|
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:
|
class FixedPatchStepper:
|
||||||
"""D2Q9 SRT/MRT FP32 fixed 2:1 refinement hierarchy, built from level lists.
|
"""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)
|
self.hierarchy = build_hierarchy(cfg, refinements, wall_aware=self._wall_aware)
|
||||||
pairs = self.hierarchy.pairs
|
pairs = self.hierarchy.pairs
|
||||||
if self._wall_aware:
|
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 pairs[0].rectangles):
|
||||||
if any(r[1] < 4 or r[3] > cfg.ny - 4 for r in rects):
|
|
||||||
raise ValueError("slot_reservoir requires y0 >= 4 and y1 <= ny-4")
|
raise ValueError("slot_reservoir requires y0 >= 4 and y1 <= ny-4")
|
||||||
self._layouts = []
|
self._layouts = _level_layouts(cfg, pairs, level_blocks, self._wall_aware)
|
||||||
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)
|
|
||||||
if self._channel:
|
if self._channel:
|
||||||
base = blocked_coarse_tables(pairs[0], self._layouts[0], np.argsort(pairs[0].coalesce_coarse))
|
base = blocked_coarse_tables(pairs[0], self._layouts[0], np.argsort(pairs[0].coalesce_coarse))
|
||||||
validate_channel_transfers(pairs[0], self._layouts[0], base[1])
|
validate_channel_transfers(pairs[0], self._layouts[0], base[1])
|
||||||
self._level_cfgs = _kernel_configs(cfg, pairs, self._wall_aware)
|
self._level_cfgs = _kernel_configs(cfg, pairs, self._wall_aware)
|
||||||
# Host boundary geometry per level in the root edge frame (Hierarchy.edge_frame).
|
hosts = _host_boundaries(self.hierarchy, self._level_cfgs, level_objects, level_jets,
|
||||||
self._boundaries = [None]
|
self._wall_aware, self._ddf_shifting)
|
||||||
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)
|
|
||||||
if not isinstance(observe_raw_exchange, bool):
|
if not isinstance(observe_raw_exchange, bool):
|
||||||
raise ValueError("observe_raw_exchange must be bool")
|
raise ValueError("observe_raw_exchange must be bool")
|
||||||
self._observe_raw_exchange = observe_raw_exchange
|
self._observe_raw_exchange = observe_raw_exchange
|
||||||
@@ -291,32 +317,25 @@ class FixedPatchStepper:
|
|||||||
self._window_valid = False
|
self._window_valid = False
|
||||||
self._window_start = None
|
self._window_start = None
|
||||||
self._exchange_valid = False
|
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
|
root = None
|
||||||
if depth == 1 and level_objects[0]:
|
if depth == 1 and level_objects[0]:
|
||||||
tables = static_circle_tables(self._layouts[0], level_objects[0])
|
tables = static_circle_tables(self._layouts[0], level_objects[0])
|
||||||
root = _BlockedBoundary(None, tables, (), False)
|
root = _BlockedBoundary(None, tables, (), False)
|
||||||
self._device_boundaries = [root] + [
|
# The one boundary/Jet record per level; host geometry and Jets are views of it.
|
||||||
_prepare_blocked_boundary(self._boundaries[level], self._layouts[level],
|
self._boundaries = [root] + [
|
||||||
|
_prepare_blocked_boundary(hosts[level], self._layouts[level],
|
||||||
self._level_cfgs[level], level_jets[level],
|
self._level_cfgs[level], level_jets[level],
|
||||||
_shifted_blocked=self._ddf_shifting)
|
_shifted_blocked=self._ddf_shifting)
|
||||||
for level in range(1, depth)]
|
for level in range(1, depth)]
|
||||||
self._leaf = [kind == ACTIVE for kind in self.hierarchy.kinds]
|
for level, b in enumerate(self._boundaries):
|
||||||
if root is not None:
|
if b is not None: # exchange scaling (dt_l = h_l) and root-edge-frame moment centers
|
||||||
self._leaf[0] &= ~root.tables["solid"]
|
b.spacing = self.hierarchy.edge_frame(level)[1]
|
||||||
for level in range(1, depth):
|
b.centers = tuple(tuple(o.geometry.center) for o in level_objects[level])
|
||||||
host = self._boundaries[level]
|
self._leaf = _fluid_leaves(self.hierarchy, self._boundaries)
|
||||||
if host is not None:
|
|
||||||
self._leaf[level] &= (host.flags.reshape(pairs[level - 1].fine_shape) & SOLID) == 0
|
|
||||||
self._state = [mask.copy() for mask in self._leaf]
|
self._state = [mask.copy() for mask in self._leaf]
|
||||||
if self._channel:
|
if self._channel:
|
||||||
self._leaf[0] &= channel_masks((cfg.ny, cfg.nx))[0]
|
self._leaf[0] &= channel_masks((cfg.ny, cfg.nx))[0]
|
||||||
# Private constructor inputs for remap_finest.
|
self._root_cfg = replace(cfg) # private constructor input for remap_finest
|
||||||
self._root_cfg = replace(cfg)
|
|
||||||
self._level_blocks = level_blocks
|
|
||||||
self.ctx = context
|
self.ctx = context
|
||||||
self._lock = RLock()
|
self._lock = RLock()
|
||||||
self._status = "new"
|
self._status = "new"
|
||||||
@@ -371,7 +390,7 @@ class FixedPatchStepper:
|
|||||||
for index, layout in enumerate(self._layouts):
|
for index, layout in enumerate(self._layouts):
|
||||||
total += 2 * 4 * int(np.prod(layout.blocked_shape))
|
total += 2 * 4 * int(np.prod(layout.blocked_shape))
|
||||||
total += layout.halo_dst.nbytes + layout.halo_src.nbytes
|
total += layout.halo_dst.nbytes + layout.halo_src.nbytes
|
||||||
b = self._device_boundaries[index]
|
b = self._boundaries[index]
|
||||||
if b is not None:
|
if b is not None:
|
||||||
total += sum(np.asarray(v).nbytes for k, v in b.tables.items()
|
total += sum(np.asarray(v).nbytes for k, v in b.tables.items()
|
||||||
if k == "flags" or k in BLOCKED_BOUNDARY_DEVICE_COLUMNS)
|
if k == "flags" or k in BLOCKED_BOUNDARY_DEVICE_COLUMNS)
|
||||||
@@ -394,7 +413,7 @@ class FixedPatchStepper:
|
|||||||
grid = ((config.nx + self._block[0] - 1) // self._block[0], config.ny, layout.nblocks)
|
grid = ((config.nx + self._block[0] - 1) // self._block[0], config.ny, layout.nblocks)
|
||||||
level = _Level(shape, module, self._allocate(size), self._allocate(size),
|
level = _Level(shape, module, self._allocate(size), self._allocate(size),
|
||||||
module.get_function(name), grid, layout)
|
module.get_function(name), grid, layout)
|
||||||
level.boundary = self._device_boundaries[index]
|
level.boundary = self._boundaries[index]
|
||||||
level.halo_copy = module.get_function("BlockHaloCopy")
|
level.halo_copy = module.get_function("BlockHaloCopy")
|
||||||
if len(layout.halo_dst):
|
if len(layout.halo_dst):
|
||||||
level.halo_dst = self._upload(layout.halo_dst)
|
level.halo_dst = self._upload(layout.halo_dst)
|
||||||
@@ -556,19 +575,17 @@ class FixedPatchStepper:
|
|||||||
pairs = self.hierarchy.pairs
|
pairs = self.hierarchy.pairs
|
||||||
if (self._wall_aware or self._channel or self.hierarchy.depth != 3
|
if (self._wall_aware or self._channel or self.hierarchy.depth != 3
|
||||||
or not isinstance(pairs[1], FixedPatchSet2D)
|
or not isinstance(pairs[1], FixedPatchSet2D)
|
||||||
or self._boundaries[1] is not None
|
or any(b is not None for b in self._boundaries[:2])):
|
||||||
or any(b is not None for b in self._device_boundaries[:2])):
|
|
||||||
raise ValueError("remap_finest requires a periodic three-level hierarchy with a "
|
raise ValueError("remap_finest requires a periodic three-level hierarchy with a "
|
||||||
"level-2 set and no level-0/1 bodies")
|
"level-2 set and no level-0/1 bodies")
|
||||||
cfg = self._root_cfg
|
cfg = self._root_cfg
|
||||||
objects = tuple(objects)
|
objects = tuple(objects)
|
||||||
actual = self._device_boundaries[2]
|
actual = self._boundaries[2]
|
||||||
jets = () if actual is None else actual.jets
|
jets = () if actual is None else actual.jets
|
||||||
old = RemapTopology.from_config(cfg, pairs[0], pairs[1].bounds,
|
old = RemapTopology.from_config(cfg, pairs[0], pairs[1].bounds,
|
||||||
child_objects=objects, child_jets=jets,
|
child_objects=objects, child_jets=jets,
|
||||||
_shifted_blocked=self._ddf_shifting)
|
_shifted_blocked=self._ddf_shifting)
|
||||||
if (not same_physical_boundary(old.boundary, self._boundaries[2])
|
if not same_physical_boundary(old.boundary, None if actual is None else actual.host):
|
||||||
or (self._boundaries[2] is not None and actual is None)):
|
|
||||||
raise ValueError("objects must reproduce the original root-edge-frame level-2 boundary")
|
raise ValueError("objects must reproduce the original root-edge-frame level-2 boundary")
|
||||||
if actual is not None:
|
if actual is not None:
|
||||||
velocity = actual.tables.get("surface_velocity")
|
velocity = actual.tables.get("surface_velocity")
|
||||||
@@ -585,7 +602,7 @@ class FixedPatchStepper:
|
|||||||
or not np.array_equal(old.surface_velocity, new.surface_velocity)):
|
or not np.array_equal(old.surface_velocity, new.surface_velocity)):
|
||||||
raise ValueError("remap requires the identical pinned physical boundary and BFL support")
|
raise ValueError("remap requires the identical pinned physical boundary and BFL support")
|
||||||
# Check sparse block alignment/launch size before any device read/allocation.
|
# 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:
|
if layout.nblocks > MAX_GRID_DIM or layout.local_shape[0] > MAX_GRID_DIM:
|
||||||
raise ValueError("level 2 blocks exceed launch limits")
|
raise ValueError("level 2 blocks exceed launch limits")
|
||||||
with self.ctx.activate():
|
with self.ctx.activate():
|
||||||
@@ -603,7 +620,7 @@ class FixedPatchStepper:
|
|||||||
phase="postcollision_g", **options)
|
phase="postcollision_g", **options)
|
||||||
candidate = type(self)(
|
candidate = type(self)(
|
||||||
cfg, context=self.ctx, refinements=[pairs[0].bounds, new.child.bounds],
|
cfg, 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],
|
level_objects=[(), (), objects], level_jets=[(), (), jets],
|
||||||
observe_raw_exchange=self._observe_raw_exchange,
|
observe_raw_exchange=self._observe_raw_exchange,
|
||||||
observe_window_exchange=self._observe_window, explosion=self._explosion)
|
observe_window_exchange=self._observe_window, explosion=self._explosion)
|
||||||
@@ -622,7 +639,7 @@ class FixedPatchStepper:
|
|||||||
return candidate
|
return candidate
|
||||||
|
|
||||||
def set_jet_amplitude(self, index, amplitude, *, level):
|
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
|
Invalid input raises before any change. A device/sync error marks the
|
||||||
instance failed (as step/initialize do); host state commits only after
|
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)
|
if (isinstance(level, bool) or not isinstance(level, Integral)
|
||||||
or not 1 <= level < self.hierarchy.depth):
|
or not 1 <= level < self.hierarchy.depth):
|
||||||
raise ValueError("Jet level must name an existing refined level (1 .. depth-1)")
|
raise ValueError("Jet level must name an existing refined level (1 .. depth-1)")
|
||||||
boundary = self._device_boundaries[level]
|
boundary = self._boundaries[level]
|
||||||
jets = list(self._level_jets[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
|
jets[index] = jets[index].with_amplitude(amplitude) # validates first
|
||||||
velocity = None
|
shape, jet_id = boundary.jet_cache
|
||||||
if boundary is not None and boundary.jet_cache is not None:
|
velocity = boundary.tables["surface_velocity"].copy()
|
||||||
shape, jet_id = boundary.jet_cache
|
driven = jet_id >= 0
|
||||||
velocity = boundary.tables["surface_velocity"].copy()
|
velocity[driven] = jet_link_velocity(jets, shape, jet_id)[driven]
|
||||||
driven = jet_id >= 0
|
# Device action: an uncertain write fails the instance closed.
|
||||||
velocity[driven] = jet_link_velocity(jets, shape, jet_id)[driven]
|
|
||||||
try:
|
try:
|
||||||
if velocity is not None:
|
self._stream.synchronize() # no in-flight step reads the column
|
||||||
self._stream.synchronize() # no in-flight step reads the column
|
cuda.memcpy_htod(boundary.device["surface_velocity"], velocity)
|
||||||
cuda.memcpy_htod(boundary.device["surface_velocity"], velocity)
|
|
||||||
except BaseException:
|
except BaseException:
|
||||||
self._status = "failed"
|
self._status = "failed"
|
||||||
raise
|
raise
|
||||||
if velocity is not None:
|
# Commit: the one boundary record holds both the column and the Jets.
|
||||||
boundary.tables["surface_velocity"] = velocity
|
boundary.tables["surface_velocity"] = velocity
|
||||||
if boundary is not None:
|
boundary.jets = tuple(jets)
|
||||||
boundary.jets = tuple(jets)
|
|
||||||
self._level_jets[level] = tuple(jets)
|
|
||||||
|
|
||||||
def _blocked_boundary(self, level, buffer, **kw):
|
def _blocked_boundary(self, level, buffer, **kw):
|
||||||
b = level.boundary
|
b = level.boundary
|
||||||
@@ -696,7 +710,7 @@ class FixedPatchStepper:
|
|||||||
b = level.boundary
|
b = level.boundary
|
||||||
if b is not None and b.last_raw is not None:
|
if b is not None and b.last_raw is not None:
|
||||||
count = len(b.tables["cell"])
|
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),
|
b.accumulate(b.last_raw, b.interval_integral, np.uint32(count),
|
||||||
np.float64(h), np.float64(h),
|
np.float64(h), np.float64(h),
|
||||||
grid=((count + self._block[0] - 1) // self._block[0], 1), **kw)
|
grid=((count + self._block[0] - 1) // self._block[0], 1), **kw)
|
||||||
@@ -740,7 +754,7 @@ class FixedPatchStepper:
|
|||||||
self._exchange_valid = False
|
self._exchange_valid = False
|
||||||
for _ in range(n):
|
for _ in range(n):
|
||||||
if self._observe_raw_exchange: # stream-ordered: holds the last interval only
|
if self._observe_raw_exchange: # stream-ordered: holds the last interval only
|
||||||
for b in self._device_boundaries:
|
for b in self._boundaries:
|
||||||
if b is not None:
|
if b is not None:
|
||||||
cuda.memset_d8_async(b.interval_integral, 0,
|
cuda.memset_d8_async(b.interval_integral, 0,
|
||||||
len(b.tables["cell"]) * 32, self._stream)
|
len(b.tables["cell"]) * 32, self._stream)
|
||||||
@@ -763,7 +777,7 @@ class FixedPatchStepper:
|
|||||||
raise RuntimeError("window exchange observation is disabled")
|
raise RuntimeError("window exchange observation is disabled")
|
||||||
self._window_valid = False
|
self._window_valid = False
|
||||||
try:
|
try:
|
||||||
for b in self._device_boundaries:
|
for b in self._boundaries:
|
||||||
if b is not None and b.window_integral is not None:
|
if b is not None and b.window_integral is not None:
|
||||||
cuda.memset_d8(b.window_integral, 0, len(b.tables["cell"]) * 32)
|
cuda.memset_d8(b.window_integral, 0, len(b.tables["cell"]) * 32)
|
||||||
except BaseException:
|
except BaseException:
|
||||||
@@ -776,7 +790,7 @@ class FixedPatchStepper:
|
|||||||
if (isinstance(level, bool) or not isinstance(level, Integral)
|
if (isinstance(level, bool) or not isinstance(level, Integral)
|
||||||
or not 0 <= level < len(self._levels)):
|
or not 0 <= level < len(self._levels)):
|
||||||
raise ValueError("level must name an existing level")
|
raise ValueError("level must name an existing level")
|
||||||
b = self._device_boundaries[level]
|
b = self._boundaries[level]
|
||||||
if (b is None or isinstance(body_id, bool) or not isinstance(body_id, Integral)
|
if (b is None or isinstance(body_id, bool) or not isinstance(body_id, Integral)
|
||||||
or not 0 <= body_id < len(b.host.objects)):
|
or not 0 <= body_id < len(b.host.objects)):
|
||||||
raise ValueError("body_id must name a body on this level")
|
raise ValueError("body_id must name a body on this level")
|
||||||
@@ -784,9 +798,9 @@ class FixedPatchStepper:
|
|||||||
|
|
||||||
def _summarize(self, b, rows, level, body_id, window, **kw):
|
def _summarize(self, b, rows, level, body_id, window, **kw):
|
||||||
return summarize_exchange(rows, b.tables["body_id"], b.tables["aperture"],
|
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,
|
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):
|
def read_exchange_window(self, body_id, *, level):
|
||||||
"""Non-destructive read of [window start, current step count] in base-lattice units."""
|
"""Non-destructive read of [window start, current step count] in base-lattice units."""
|
||||||
@@ -873,7 +887,7 @@ class FixedPatchStepper:
|
|||||||
self._stream.synchronize()
|
self._stream.synchronize()
|
||||||
except Exception as exc:
|
except Exception as exc:
|
||||||
error = exc
|
error = exc
|
||||||
for b in getattr(self, "_device_boundaries", ()):
|
for b in getattr(self, "_boundaries", ()):
|
||||||
if b is not None:
|
if b is not None:
|
||||||
b.flags = b.device = b.kernel = None
|
b.flags = b.device = b.kernel = None
|
||||||
b.last_raw = b.interval_integral = b.accumulate = None
|
b.last_raw = b.interval_integral = b.accumulate = None
|
||||||
@@ -901,9 +915,20 @@ class FixedPatchStepper:
|
|||||||
raise
|
raise
|
||||||
|
|
||||||
# Read-only per-level state, indexed by level (0 = root). Pairs: hierarchy.pairs[l-1].
|
# Read-only per-level state, indexed by level (0 = root). Pairs: hierarchy.pairs[l-1].
|
||||||
layouts = property(lambda self: tuple(self._layouts))
|
@property
|
||||||
boundaries = property(lambda self: tuple(self._boundaries),
|
def layouts(self):
|
||||||
doc="Host boundary geometry per level (None: no bodies/Jets).")
|
return tuple(self._layouts)
|
||||||
jets = property(lambda self: tuple(self._level_jets))
|
|
||||||
leaf = property(lambda self: tuple(m.copy() for m in self._leaf),
|
@property
|
||||||
doc="Per-level (ny, nx) masks of public (non-NaN) download cells.")
|
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)
|
||||||
@@ -130,6 +130,8 @@ class _BlockedBoundary:
|
|||||||
interval_integral: object = None
|
interval_integral: object = None
|
||||||
accumulate: object = None
|
accumulate: object = None
|
||||||
window_integral: object = None # optional observation-window FP64 ledger
|
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):
|
def _prepare_blocked_boundary(boundary, layout, cfg, jets, *, _shifted_blocked=False):
|
||||||
|
|||||||
@@ -102,7 +102,7 @@ def test_simulation_before_initialize_commits_host_only(monkeypatch):
|
|||||||
|
|
||||||
|
|
||||||
def stepper_state(run):
|
def stepper_state(run):
|
||||||
boundary = run._device_boundaries[2]
|
boundary = run._boundaries[2]
|
||||||
return boundary.tables["surface_velocity"].copy(), boundary.jets, run.jets[2]
|
return boundary.tables["surface_velocity"].copy(), boundary.jets, run.jets[2]
|
||||||
|
|
||||||
|
|
||||||
@@ -116,8 +116,8 @@ def test_stepper_invalid_input_writes_nothing_and_stays_ready():
|
|||||||
with pytest.raises(error):
|
with pytest.raises(error):
|
||||||
run.set_jet_amplitude(*args, **kw)
|
run.set_jet_amplitude(*args, **kw)
|
||||||
module.cuda.memcpy_htod.assert_not_called()
|
module.cuda.memcpy_htod.assert_not_called()
|
||||||
np.testing.assert_array_equal(run._device_boundaries[2].tables["surface_velocity"], velocity)
|
np.testing.assert_array_equal(run._boundaries[2].tables["surface_velocity"], velocity)
|
||||||
assert (run._device_boundaries[2].jets, run.jets[2]) == (jets, level_jets)
|
assert (run._boundaries[2].jets, run.jets[2]) == (jets, level_jets)
|
||||||
assert run._status == "ready"
|
assert run._status == "ready"
|
||||||
run.set_jet_amplitude(0, .017, level=2)
|
run.set_jet_amplitude(0, .017, level=2)
|
||||||
assert run.jets[2][0].amplitude == .017
|
assert run.jets[2][0].amplitude == .017
|
||||||
@@ -136,8 +136,8 @@ def test_stepper_device_failure_keeps_host_and_blocks_reuse(where):
|
|||||||
run._stream.synchronize.side_effect = RuntimeError("upload failed")
|
run._stream.synchronize.side_effect = RuntimeError("upload failed")
|
||||||
with pytest.raises(RuntimeError, match="upload failed"):
|
with pytest.raises(RuntimeError, match="upload failed"):
|
||||||
run.set_jet_amplitude(0, .017, level=2)
|
run.set_jet_amplitude(0, .017, level=2)
|
||||||
np.testing.assert_array_equal(run._device_boundaries[2].tables["surface_velocity"], velocity)
|
np.testing.assert_array_equal(run._boundaries[2].tables["surface_velocity"], velocity)
|
||||||
assert (run._device_boundaries[2].jets, run.jets[2]) == (jets, level_jets)
|
assert (run._boundaries[2].jets, run.jets[2]) == (jets, level_jets)
|
||||||
assert run._status == "failed"
|
assert run._status == "failed"
|
||||||
module.cuda.memcpy_htod.side_effect = None
|
module.cuda.memcpy_htod.side_effect = None
|
||||||
run._stream.synchronize.side_effect = None
|
run._stream.synchronize.side_effect = None
|
||||||
|
|||||||
@@ -96,11 +96,11 @@ class MultiplePinnedRemap(unittest.TestCase):
|
|||||||
before = run.download()
|
before = run.download()
|
||||||
clone = run.remap_finest((PINS+MOVED+ADDED)[::-1], objects=bodies)
|
clone = run.remap_finest((PINS+MOVED+ADDED)[::-1], objects=bodies)
|
||||||
try:
|
try:
|
||||||
self.assertEqual(tuple(j.amplitude for j in clone._device_boundaries[2].jets), (.017, -.011))
|
self.assertEqual(tuple(j.amplitude for j in clone._boundaries[2].jets), (.017, -.011))
|
||||||
np.testing.assert_array_equal(clone._device_boundaries[2].tables['surface_velocity'],
|
np.testing.assert_array_equal(clone._boundaries[2].tables['surface_velocity'],
|
||||||
run._device_boundaries[2].tables['surface_velocity'])
|
run._boundaries[2].tables['surface_velocity'])
|
||||||
clone.set_jet_amplitude(1, .005, level=2)
|
clone.set_jet_amplitude(1, .005, level=2)
|
||||||
self.assertEqual(run._device_boundaries[2].jets[1].amplitude, -.011)
|
self.assertEqual(run._boundaries[2].jets[1].amplitude, -.011)
|
||||||
clone.step()
|
clone.step()
|
||||||
self.assertEqual((run.step_count, clone.step_count), (7, 8))
|
self.assertEqual((run.step_count, clone.step_count), (7, 8))
|
||||||
for got, want in zip(run.download(), before):
|
for got, want in zip(run.download(), before):
|
||||||
|
|||||||
@@ -62,9 +62,9 @@ class PinnedBodyRemap(unittest.TestCase):
|
|||||||
run.set_jet_amplitude(0,.017,level=2)
|
run.set_jet_amplitude(0,.017,level=2)
|
||||||
before=run.download()
|
before=run.download()
|
||||||
clone=run.remap_finest(B,objects=objects)
|
clone=run.remap_finest(B,objects=objects)
|
||||||
self.assertEqual(clone._device_boundaries[2].jets[0].amplitude,.017)
|
self.assertEqual(clone._boundaries[2].jets[0].amplitude,.017)
|
||||||
np.testing.assert_array_equal(clone._device_boundaries[2].tables['surface_velocity'],
|
np.testing.assert_array_equal(clone._boundaries[2].tables['surface_velocity'],
|
||||||
run._device_boundaries[2].tables['surface_velocity'])
|
run._boundaries[2].tables['surface_velocity'])
|
||||||
self.assertTrue(same_physical_boundary(run.boundaries[2],clone.boundaries[2]))
|
self.assertTrue(same_physical_boundary(run.boundaries[2],clone.boundaries[2]))
|
||||||
clone.step()
|
clone.step()
|
||||||
self.assertEqual((run.step_count,clone.step_count),(7,8))
|
self.assertEqual((run.step_count,clone.step_count),(7,8))
|
||||||
@@ -74,8 +74,8 @@ class PinnedBodyRemap(unittest.TestCase):
|
|||||||
def test_static_body_without_actuation(self):
|
def test_static_body_without_actuation(self):
|
||||||
run,module,allocations,storage,objects=ready(jets=())
|
run,module,allocations,storage,objects=ready(jets=())
|
||||||
clone=run.remap_finest(A,objects=objects)
|
clone=run.remap_finest(A,objects=objects)
|
||||||
self.assertFalse(clone._device_boundaries[2].actuated)
|
self.assertFalse(clone._boundaries[2].actuated)
|
||||||
self.assertNotIn('surface_velocity',clone._device_boundaries[2].tables)
|
self.assertNotIn('surface_velocity',clone._boundaries[2].tables)
|
||||||
for a,b in zip(run.download(),clone.download()): np.testing.assert_array_equal(a,b)
|
for a,b in zip(run.download(),clone.download()): np.testing.assert_array_equal(a,b)
|
||||||
clone.close();run.close()
|
clone.close();run.close()
|
||||||
|
|
||||||
|
|||||||
@@ -85,7 +85,7 @@ class RawExchange(unittest.TestCase):
|
|||||||
|
|
||||||
def test_four_level_schedule_publication_and_fixed_storage(self):
|
def test_four_level_schedule_publication_and_fixed_storage(self):
|
||||||
run, host, calls, arrays = backend()
|
run, host, calls, arrays = backend()
|
||||||
b = run._device_boundaries[3]
|
b = run._boundaries[3]
|
||||||
rows = len(b.tables['cell'])
|
rows = len(b.tables['cell'])
|
||||||
self.assertEqual(b.last_raw.size, rows * 16)
|
self.assertEqual(b.last_raw.size, rows * 16)
|
||||||
self.assertEqual(b.interval_integral.size, rows * 32)
|
self.assertEqual(b.interval_integral.size, rows * 32)
|
||||||
@@ -128,7 +128,7 @@ class RawExchange(unittest.TestCase):
|
|||||||
|
|
||||||
def test_disabled_has_nulls_and_no_observation_allocation_or_launch(self):
|
def test_disabled_has_nulls_and_no_observation_allocation_or_launch(self):
|
||||||
run, host, calls, _ = backend(False)
|
run, host, calls, _ = backend(False)
|
||||||
b = run._device_boundaries[3]
|
b = run._boundaries[3]
|
||||||
self.assertIs(b.last_raw, b.interval_integral)
|
self.assertIs(b.last_raw, b.interval_integral)
|
||||||
self.assertIs(b.interval_integral, None)
|
self.assertIs(b.interval_integral, None)
|
||||||
run.step()
|
run.step()
|
||||||
|
|||||||
@@ -1,5 +1,6 @@
|
|||||||
"""Host guards and source contract for MRT reuse in refined/blocked stepping (no CUDA)."""
|
"""Host guards and source contract for MRT reuse in refined/blocked stepping (no CUDA)."""
|
||||||
from pathlib import Path
|
from pathlib import Path
|
||||||
|
import re
|
||||||
import unittest
|
import unittest
|
||||||
|
|
||||||
from CelerisLab.config import LBMConfig
|
from CelerisLab.config import LBMConfig
|
||||||
@@ -32,22 +33,25 @@ class RefinementMRTHost(unittest.TestCase):
|
|||||||
"collide_srt_no_force(f, feq, d_params.omega);"):
|
"collide_srt_no_force(f, feq, d_params.omega);"):
|
||||||
self.assertIn(text, block_collision)
|
self.assertIn(text, block_collision)
|
||||||
self.assertIn("collide_srt(f, feq, forcing, omega);", block_collision)
|
self.assertIn("collide_srt(f, feq, forcing, omega);", block_collision)
|
||||||
self.assertIn("#elif BLOCK_ROOT_BC == 2", block_collision)
|
self.assertIn("#if COLLISION_MODEL == 2 || BLOCK_ROOT_BC == 2", block_collision)
|
||||||
self.assertEqual(channel.count("block_collide(f);"), 1)
|
self.assertEqual(channel.count("block_collide(f);"), 1)
|
||||||
self.assertNotIn("collide_mrt_unforced(f);", channel)
|
# Executable text only: the header comment explains the absent force read.
|
||||||
|
code = re.sub(r"//[^\n]*", "", channel + "\n" + block_collision)
|
||||||
|
self.assertIn("d_params.omega", code)
|
||||||
|
for text in ("collide_mrt_unforced", "collide_dispatch(", "d_params.fx", "d_params.fy"):
|
||||||
|
self.assertNotIn(text, code)
|
||||||
for source in (blocks,):
|
for source in (blocks,):
|
||||||
self.assertIn("(COLLISION_MODEL != 0 && COLLISION_MODEL != 2)", source)
|
self.assertIn("(COLLISION_MODEL != 0 && COLLISION_MODEL != 2)", source)
|
||||||
self.assertIn("USE_LES", source)
|
self.assertIn("USE_LES", source)
|
||||||
self.assertIn("STORE_PRECISION != 0", source)
|
self.assertIn("STORE_PRECISION != 0", source)
|
||||||
self.assertIn('#if COLLISION_MODEL == 2\n#include "operators/collision_mrt.cuh"\n#endif', 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.
|
# One unforced MRT call, in block_collide only; the steps never collide directly.
|
||||||
mrt = (KERNELS / "operators/collision_mrt.cuh").read_text()
|
self.assertEqual(block_collision.count("collide_mrt_no_force(f, rho, ux, uy, omega);"), 1)
|
||||||
self.assertEqual(mrt.count("void collide_mrt_unforced("), 1)
|
self.assertEqual(block_collision.count("compute_rho_u(f, rho, ux, uy);"), 1)
|
||||||
for source in (blocks + "\n" + block_collision, channel + "\n" + block_collision):
|
self.assertEqual(block_collision.count("fmaxf(OMEGA_COLLISION_MIN, d_params.omega)"), 1)
|
||||||
self.assertEqual(source.count("collide_mrt_unforced(f);"), 1)
|
for source in (blocks, channel):
|
||||||
self.assertNotIn("collide_mrt_no_force(", source)
|
self.assertNotIn("collide_mrt", source)
|
||||||
self.assertNotIn("compute_rho_u(f, rho, ux, uy);\n collide_mrt", source)
|
self.assertNotIn("collide_srt", source)
|
||||||
self.assertNotRegex(source, r"compute_rho_u\(f, rho, ux, uy\);\s*collide_mrt")
|
|
||||||
|
|
||||||
|
|
||||||
if __name__ == "__main__": unittest.main()
|
if __name__ == "__main__": unittest.main()
|
||||||
@@ -1,5 +1,6 @@
|
|||||||
"""CPU mocks execute production remap lifecycle; no GPU numerical claim."""
|
"""CPU mocks execute production remap lifecycle; no GPU numerical claim."""
|
||||||
from dataclasses import replace
|
from dataclasses import replace
|
||||||
|
from types import SimpleNamespace
|
||||||
import unittest
|
import unittest
|
||||||
from unittest.mock import patch
|
from unittest.mock import patch
|
||||||
|
|
||||||
@@ -69,10 +70,11 @@ class RuntimeChildRemap(unittest.TestCase):
|
|||||||
run._status = 'ready'
|
run._status = 'ready'
|
||||||
for bad in ([], [(40,16,48,40)], [(16,16,32,24),(16,27,32,55)]):
|
for bad in ([], [(40,16,48,40)], [(16,16,32,24),(16,27,32,55)]):
|
||||||
with self.assertRaises(ValueError): run.remap_finest(bad)
|
with self.assertRaises(ValueError): run.remap_finest(bad)
|
||||||
hosts, pairs = run._boundaries, run.hierarchy.pairs
|
records, pairs = run._boundaries, run.hierarchy.pairs
|
||||||
for attribute, value in (('_wall_aware', True),
|
for attribute, value in (('_wall_aware', True),
|
||||||
('_boundaries', [*hosts[:2], object()]), # child boundary
|
# level-2 record whose host geometry differs from objects=()
|
||||||
('_boundaries', [hosts[0], object(), hosts[2]]), # fine boundary
|
('_boundaries', [*records[:2], SimpleNamespace(host=object(), jets=())]),
|
||||||
|
('_boundaries', [records[0], object(), records[2]]), # level-1 body
|
||||||
('hierarchy', replace(run.hierarchy, pairs=(pairs[0], None)))):
|
('hierarchy', replace(run.hierarchy, pairs=(pairs[0], None)))):
|
||||||
with patch.object(run, attribute, value):
|
with patch.object(run, attribute, value):
|
||||||
with self.assertRaises(ValueError): run.remap_finest(B)
|
with self.assertRaises(ValueError): run.remap_finest(B)
|
||||||
|
|||||||
@@ -67,7 +67,7 @@ class ShiftedJets(unittest.TestCase):
|
|||||||
def test_switch_upload_only_then_raw_remap_keeps_live_velocity(self):
|
def test_switch_upload_only_then_raw_remap_keeps_live_velocity(self):
|
||||||
run,module,allocations,storage,objects=ready();clone=None
|
run,module,allocations,storage,objects=ready();clone=None
|
||||||
try:
|
try:
|
||||||
boundary=run._device_boundaries[2]
|
boundary=run._boundaries[2]
|
||||||
before=run._download_raw()
|
before=run._download_raw()
|
||||||
original=boundary.tables['surface_velocity'].copy()
|
original=boundary.tables['surface_velocity'].copy()
|
||||||
module.cuda.memcpy_htod.reset_mock()
|
module.cuda.memcpy_htod.reset_mock()
|
||||||
@@ -79,7 +79,7 @@ class ShiftedJets(unittest.TestCase):
|
|||||||
for a,b in zip(before.values,run._download_raw().values):self.assertEqual(a.tobytes(),b.tobytes())
|
for a,b in zip(before.values,run._download_raw().values):self.assertEqual(a.tobytes(),b.tobytes())
|
||||||
clone=run.remap_finest(B,objects=objects)
|
clone=run.remap_finest(B,objects=objects)
|
||||||
self.assertEqual(clone.jets[2][0].amplitude,-jet().amplitude)
|
self.assertEqual(clone.jets[2][0].amplitude,-jet().amplitude)
|
||||||
np.testing.assert_array_equal(clone._device_boundaries[2].tables['surface_velocity'],
|
np.testing.assert_array_equal(clone._boundaries[2].tables['surface_velocity'],
|
||||||
boundary.tables['surface_velocity'])
|
boundary.tables['surface_velocity'])
|
||||||
self.assertEqual(clone.step_count,7)
|
self.assertEqual(clone.step_count,7)
|
||||||
saved=run._download_raw()
|
saved=run._download_raw()
|
||||||
@@ -92,12 +92,12 @@ class ShiftedJets(unittest.TestCase):
|
|||||||
def test_invalid_amplitude_does_not_upload_or_commit(self):
|
def test_invalid_amplitude_does_not_upload_or_commit(self):
|
||||||
run,module,_,_,_=ready()
|
run,module,_,_,_=ready()
|
||||||
try:
|
try:
|
||||||
before=run._device_boundaries[2].tables['surface_velocity'].copy()
|
before=run._boundaries[2].tables['surface_velocity'].copy()
|
||||||
module.cuda.memcpy_htod.reset_mock()
|
module.cuda.memcpy_htod.reset_mock()
|
||||||
for amplitude in (float('nan'),float('inf'),1.):
|
for amplitude in (float('nan'),float('inf'),1.):
|
||||||
with self.assertRaises(ValueError):run.set_jet_amplitude(0,amplitude,level=2)
|
with self.assertRaises(ValueError):run.set_jet_amplitude(0,amplitude,level=2)
|
||||||
module.cuda.memcpy_htod.assert_not_called()
|
module.cuda.memcpy_htod.assert_not_called()
|
||||||
np.testing.assert_array_equal(before,run._device_boundaries[2].tables['surface_velocity'])
|
np.testing.assert_array_equal(before,run._boundaries[2].tables['surface_velocity'])
|
||||||
finally:run.close()
|
finally:run.close()
|
||||||
|
|
||||||
|
|
||||||
|
|||||||
@@ -73,7 +73,7 @@ class SingleLevel(unittest.TestCase):
|
|||||||
|
|
||||||
def test_circle_solids_are_not_state(self):
|
def test_circle_solids_are_not_state(self):
|
||||||
run, _, _ = backend((circle(),))
|
run, _, _ = backend((circle(),))
|
||||||
solid = run._device_boundaries[0].tables["solid"]
|
solid = run._boundaries[0].tables["solid"]
|
||||||
self.assertTrue(solid.any())
|
self.assertTrue(solid.any())
|
||||||
np.testing.assert_array_equal(run._state[0], ~solid)
|
np.testing.assert_array_equal(run._state[0], ~solid)
|
||||||
|
|
||||||
|
|||||||
@@ -37,7 +37,7 @@ class UnionBodySupport(unittest.TestCase):
|
|||||||
run.set_jet_amplitude(0,.017,level=2)
|
run.set_jet_amplitude(0,.017,level=2)
|
||||||
clone=run.remap_finest(((16,16,24,40),(24,16,32,40),A[1]),objects=objects)
|
clone=run.remap_finest(((16,16,24,40),(24,16,32,40),A[1]),objects=objects)
|
||||||
try:
|
try:
|
||||||
self.assertEqual(clone._device_boundaries[2].jets[0].amplitude,.017)
|
self.assertEqual(clone._boundaries[2].jets[0].amplitude,.017)
|
||||||
self.assertTrue(same_physical_boundary(run.boundaries[2],clone.boundaries[2]))
|
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)
|
for a,b in zip(run.download(),clone.download()):np.testing.assert_array_equal(a,b)
|
||||||
finally:clone.close()
|
finally:clone.close()
|
||||||
|
|||||||
@@ -240,7 +240,7 @@ class SingleLevelGPU(unittest.TestCase):
|
|||||||
with CudaContext(0) as context:
|
with CudaContext(0) as context:
|
||||||
for split in SPLITS:
|
for split in SPLITS:
|
||||||
with single(context=context, block=split, objects=[circle()]) as backend:
|
with single(context=context, block=split, objects=[circle()]) as backend:
|
||||||
tables = backend._device_boundaries[0].tables
|
tables = backend._boundaries[0].tables
|
||||||
for key, values, dtype in (("q", [l.q for l in links], np.float32),
|
for key, values, dtype in (("q", [l.q for l in links], np.float32),
|
||||||
("dir", [l.dir for l in links], np.uint8),
|
("dir", [l.dir for l in links], np.uint8),
|
||||||
("fallback", [l.fallback for l in links], np.uint8)):
|
("fallback", [l.fallback for l in links], np.uint8)):
|
||||||
|
|||||||
@@ -59,7 +59,7 @@ class ExchangeWindowCPU(unittest.TestCase):
|
|||||||
stack.enter_context(patch.object(cuda, "memcpy_dtoh",
|
stack.enter_context(patch.object(cuda, "memcpy_dtoh",
|
||||||
side_effect=lambda out, src: out.fill(0)))
|
side_effect=lambda out, src: out.fill(0)))
|
||||||
backend = build(m)
|
backend = build(m)
|
||||||
b = backend._device_boundaries[1]
|
b = backend._boundaries[1]
|
||||||
self.assertIsNotNone(b.window_integral)
|
self.assertIsNotNone(b.window_integral)
|
||||||
with self.assertRaises(RuntimeError):
|
with self.assertRaises(RuntimeError):
|
||||||
backend.start_exchange_window() # not initialized
|
backend.start_exchange_window() # not initialized
|
||||||
@@ -96,7 +96,7 @@ class ExchangeWindowCPU(unittest.TestCase):
|
|||||||
with mocked_backend() as m, ExitStack() as stack:
|
with mocked_backend() as m, ExitStack() as stack:
|
||||||
stack.enter_context(patch.object(cuda, "memset_d8_async"))
|
stack.enter_context(patch.object(cuda, "memset_d8_async"))
|
||||||
backend = build(m, window=False)
|
backend = build(m, window=False)
|
||||||
self.assertIsNone(backend._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.initialize(*[np.full(l.shape, 1 / 9, np.float32) for l in backend._levels])
|
||||||
backend.step(1)
|
backend.step(1)
|
||||||
self.assertEqual(len(accumulations(m)), 2)
|
self.assertEqual(len(accumulations(m)), 2)
|
||||||
@@ -120,7 +120,7 @@ class ExchangeWindowGPU(unittest.TestCase):
|
|||||||
windowed.step(2)
|
windowed.step(2)
|
||||||
plain.step(2)
|
plain.step(2)
|
||||||
windowed.start_exchange_window()
|
windowed.start_exchange_window()
|
||||||
b = windowed._device_boundaries[1]
|
b = windowed._boundaries[1]
|
||||||
count = len(b.tables["cell"])
|
count = len(b.tables["cell"])
|
||||||
factors = np.array([.25, .25, .125, .25]) # h*dt, h*h*dt, h*h for level-1 h=.5
|
factors = np.array([.25, .25, .125, .25]) # h*dt, h*h*dt, h*h for level-1 h=.5
|
||||||
scale = np.zeros(4)
|
scale = np.zeros(4)
|
||||||
@@ -212,7 +212,7 @@ class ExchangeWindowShiftedMRTGPU(unittest.TestCase):
|
|||||||
run.step(2)
|
run.step(2)
|
||||||
if run is not plain:
|
if run is not plain:
|
||||||
run.start_exchange_window()
|
run.start_exchange_window()
|
||||||
b = windowed._device_boundaries[1]
|
b = windowed._boundaries[1]
|
||||||
count = len(b.tables["cell"])
|
count = len(b.tables["cell"])
|
||||||
record = dict(scale=np.zeros(4), substeps=0)
|
record = dict(scale=np.zeros(4), substeps=0)
|
||||||
audit_substeps(windowed, b, count, record)
|
audit_substeps(windowed, b, count, record)
|
||||||
@@ -228,7 +228,7 @@ class ExchangeWindowShiftedMRTGPU(unittest.TestCase):
|
|||||||
same_bits(self, a, e)
|
same_bits(self, a, e)
|
||||||
rows = [np.empty((count, 4), np.float64) for _ in range(2)]
|
rows = [np.empty((count, 4), np.float64) for _ in range(2)]
|
||||||
cuda.memcpy_dtoh(rows[0], b.window_integral)
|
cuda.memcpy_dtoh(rows[0], b.window_integral)
|
||||||
cuda.memcpy_dtoh(rows[1], split._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))
|
np.testing.assert_array_equal(rows[0].view(np.uint64), rows[1].view(np.uint64))
|
||||||
|
|
||||||
out = windowed.read_exchange_window(0, level=1)
|
out = windowed.read_exchange_window(0, level=1)
|
||||||
|
|||||||
@@ -57,7 +57,7 @@ def make(context):
|
|||||||
def topology(run, rectangles, objects):
|
def topology(run, rectangles, objects):
|
||||||
return RemapTopology.from_config(
|
return RemapTopology.from_config(
|
||||||
CFG, run.hierarchy.pairs[0], rectangles, child_objects=objects,
|
CFG, run.hierarchy.pairs[0], rectangles, child_objects=objects,
|
||||||
child_jets=run._device_boundaries[2].jets,
|
child_jets=run._boundaries[2].jets,
|
||||||
)
|
)
|
||||||
|
|
||||||
|
|
||||||
@@ -74,7 +74,7 @@ class MultiplePinnedRemapGPU(unittest.TestCase):
|
|||||||
def change_amplitudes(self, run, amplitudes):
|
def change_amplitudes(self, run, amplitudes):
|
||||||
# Each jet selects a different body. Updating one must preserve the
|
# Each jet selects a different body. Updating one must preserve the
|
||||||
# other body's current per-link velocities, including its driven links.
|
# other body's current per-link velocities, including its driven links.
|
||||||
boundary = run._device_boundaries[2]
|
boundary = run._boundaries[2]
|
||||||
body_ids = run.boundaries[2].curved_host[3]
|
body_ids = run.boundaries[2].curved_host[3]
|
||||||
for index, amplitude in enumerate(amplitudes):
|
for index, amplitude in enumerate(amplitudes):
|
||||||
before = boundary.tables['surface_velocity'].copy()
|
before = boundary.tables['surface_velocity'].copy()
|
||||||
@@ -91,7 +91,7 @@ class MultiplePinnedRemapGPU(unittest.TestCase):
|
|||||||
self.assertTrue(same_physical_boundary(new.boundaries[2], mesh.boundary))
|
self.assertTrue(same_physical_boundary(new.boundaries[2], mesh.boundary))
|
||||||
# Compare physical cut-link rows, never layout-dependent blocked addresses.
|
# Compare physical cut-link rows, never layout-dependent blocked addresses.
|
||||||
for run in (old, new):
|
for run in (old, new):
|
||||||
boundary = run._device_boundaries[2]
|
boundary = run._boundaries[2]
|
||||||
self.assertEqual(tuple(jet.amplitude for jet in boundary.jets), amplitudes)
|
self.assertEqual(tuple(jet.amplitude for jet in boundary.jets), amplitudes)
|
||||||
exact(self, boundary.tables['surface_velocity'], mesh.surface_velocity,
|
exact(self, boundary.tables['surface_velocity'], mesh.surface_velocity,
|
||||||
'latest physical per-link velocities')
|
'latest physical per-link velocities')
|
||||||
|
|||||||
@@ -39,9 +39,9 @@ class PinnedBodyRemapGPU(unittest.TestCase):
|
|||||||
exact(self,getattr(after.boundaries[2],name),getattr(before.boundaries[2],name),name)
|
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):
|
for a,b in zip(after.boundaries[2].curved_host,before.boundaries[2].curved_host):
|
||||||
exact(self,a,b,'physical cut links')
|
exact(self,a,b,'physical cut links')
|
||||||
exact(self,after._device_boundaries[2].tables['surface_velocity'],
|
exact(self,after._boundaries[2].tables['surface_velocity'],
|
||||||
before._device_boundaries[2].tables['surface_velocity'],'current per-link velocities')
|
before._boundaries[2].tables['surface_velocity'],'current per-link velocities')
|
||||||
self.assertEqual(after._device_boundaries[2].jets[0].amplitude,.017)
|
self.assertEqual(after._boundaries[2].jets[0].amplitude,.017)
|
||||||
|
|
||||||
def test_noop_next_step_bitwise(self):
|
def test_noop_next_step_bitwise(self):
|
||||||
with CudaContext(0) as ctx, ExitStack() as owned:
|
with CudaContext(0) as ctx, ExitStack() as owned:
|
||||||
@@ -63,7 +63,7 @@ class PinnedBodyRemapGPU(unittest.TestCase):
|
|||||||
current.set_jet_amplitude(0,.017,level=2)
|
current.set_jet_amplitude(0,.017,level=2)
|
||||||
for rects in (B,A):
|
for rects in (B,A):
|
||||||
previous=current;before=previous.download();count=previous.step_count
|
previous=current;before=previous.download();count=previous.step_count
|
||||||
jets=previous._device_boundaries[2].jets
|
jets=previous._boundaries[2].jets
|
||||||
old=RemapTopology.from_config(CFG,previous.hierarchy.pairs[0],previous.hierarchy.pairs[1].bounds,
|
old=RemapTopology.from_config(CFG,previous.hierarchy.pairs[0],previous.hierarchy.pairs[1].bounds,
|
||||||
child_objects=objects,child_jets=jets)
|
child_objects=objects,child_jets=jets)
|
||||||
new=RemapTopology.from_config(CFG,previous.hierarchy.pairs[0],rects,
|
new=RemapTopology.from_config(CFG,previous.hierarchy.pairs[0],rects,
|
||||||
|
|||||||
@@ -105,7 +105,7 @@ class RawExchangeGPU(unittest.TestCase):
|
|||||||
arrays = [state(l.shape) for l in measured._levels]
|
arrays = [state(l.shape) for l in measured._levels]
|
||||||
measured.initialize(*arrays)
|
measured.initialize(*arrays)
|
||||||
plain.initialize(*arrays)
|
plain.initialize(*arrays)
|
||||||
b = measured._device_boundaries[1]
|
b = measured._boundaries[1]
|
||||||
count = len(b.tables['cell'])
|
count = len(b.tables['cell'])
|
||||||
kernel = measured._levels[1].kernel
|
kernel = measured._levels[1].kernel
|
||||||
samples = []
|
samples = []
|
||||||
@@ -158,7 +158,7 @@ class RawExchangeGPU(unittest.TestCase):
|
|||||||
FixedPatchStepper(cfg, context=ctx, **settings) as plain:
|
FixedPatchStepper(cfg, context=ctx, **settings) as plain:
|
||||||
arrays = [state(l.shape, ux=0.) for l in run._levels]
|
arrays = [state(l.shape, ux=0.) for l in run._levels]
|
||||||
run.initialize(*arrays); plain.initialize(*arrays)
|
run.initialize(*arrays); plain.initialize(*arrays)
|
||||||
b = run._device_boundaries[3]
|
b = run._boundaries[3]
|
||||||
count = len(b.tables['cell'])
|
count = len(b.tables['cell'])
|
||||||
kernel = run._levels[3].kernel
|
kernel = run._levels[3].kernel
|
||||||
samples, independent_samples = [], []
|
samples, independent_samples = [], []
|
||||||
@@ -213,7 +213,7 @@ class RawExchangeGPU(unittest.TestCase):
|
|||||||
FixedPatchStepper(CFG, context=ctx, **fine((16, 20)), observe_raw_exchange=True) as amr:
|
FixedPatchStepper(CFG, context=ctx, **fine((16, 20)), observe_raw_exchange=True) as amr:
|
||||||
amr.initialize(*(state(l.shape) for l in amr._levels))
|
amr.initialize(*(state(l.shape) for l in amr._levels))
|
||||||
uniform.initialize(state((9, 48, 64)))
|
uniform.initialize(state((9, 48, 64)))
|
||||||
b = amr._device_boundaries[1]
|
b = amr._boundaries[1]
|
||||||
count = len(b.tables['cell'])
|
count = len(b.tables['cell'])
|
||||||
links = obj.geometry.build_cut_links(64, 48)
|
links = obj.geometry.build_cut_links(64, 48)
|
||||||
np.testing.assert_array_equal([v.dir for v in links], b.tables['dir'])
|
np.testing.assert_array_equal([v.dir for v in links], b.tables['dir'])
|
||||||
@@ -221,7 +221,7 @@ class RawExchangeGPU(unittest.TestCase):
|
|||||||
# Same fine-grid force arms after translation to the same center: the
|
# Same fine-grid force arms after translation to the same center: the
|
||||||
# single level's static BFL gets the AMR Jet column and force arms; the raw
|
# single level's static BFL gets the AMR Jet column and force arms; the raw
|
||||||
# output pointer is substituted at launch (no stepper-side accumulation).
|
# output pointer is substituted at launch (no stepper-side accumulation).
|
||||||
root = uniform._device_boundaries[0]
|
root = uniform._boundaries[0]
|
||||||
velocity = uniform._upload(b.tables['surface_velocity'])
|
velocity = uniform._upload(b.tables['surface_velocity'])
|
||||||
root.device.update(surface_velocity=velocity, rx=uniform._upload(b.tables['rx']),
|
root.device.update(surface_velocity=velocity, rx=uniform._upload(b.tables['rx']),
|
||||||
ry=uniform._upload(b.tables['ry']))
|
ry=uniform._upload(b.tables['ry']))
|
||||||
|
|||||||
@@ -187,8 +187,8 @@ class BlockedFineCPU(unittest.TestCase):
|
|||||||
backend, calls = self.mocked_backend([circle(spec)])
|
backend, calls = self.mocked_backend([circle(spec)])
|
||||||
c, f = backend._levels
|
c, f = backend._levels
|
||||||
coarse, cur, spare = c.current, f.current, f.spare
|
coarse, cur, spare = c.current, f.current, f.spare
|
||||||
flags, kind = backend._device_boundaries[1].flags, backend._levels[1].kind
|
flags, kind = backend._boundaries[1].flags, backend._levels[1].kind
|
||||||
velocity = backend._device_boundaries[1].device.get("surface_velocity")
|
velocity = backend._boundaries[1].device.get("surface_velocity")
|
||||||
self.assertEqual(velocity is None, spec is None)
|
self.assertEqual(velocity is None, spec is None)
|
||||||
backend.step(1)
|
backend.step(1)
|
||||||
backend.close()
|
backend.close()
|
||||||
@@ -213,8 +213,8 @@ class BlockedFineCPU(unittest.TestCase):
|
|||||||
def test_actuation_without_velocity_buffer_fails_before_bfl(self):
|
def test_actuation_without_velocity_buffer_fails_before_bfl(self):
|
||||||
for speed in (0.0, .002): # any spec, even zero amplitude
|
for speed in (0.0, .002): # any spec, even zero amplitude
|
||||||
backend, calls = self.mocked_backend([circle(actuation(speed))])
|
backend, calls = self.mocked_backend([circle(actuation(speed))])
|
||||||
self.assertIsNotNone(backend._device_boundaries[1].tables["surface_velocity"])
|
self.assertIsNotNone(backend._boundaries[1].tables["surface_velocity"])
|
||||||
del backend._device_boundaries[1].device["surface_velocity"]
|
del backend._boundaries[1].device["surface_velocity"]
|
||||||
with self.assertRaisesRegex(RuntimeError, "Surface velocity device buffer is missing"):
|
with self.assertRaisesRegex(RuntimeError, "Surface velocity device buffer is missing"):
|
||||||
backend.step(1)
|
backend.step(1)
|
||||||
names = [name for name, _ in calls]
|
names = [name for name, _ in calls]
|
||||||
|
|||||||
@@ -62,8 +62,8 @@ class JetStepperCPU(unittest.TestCase):
|
|||||||
C = (16.0, 12.0)
|
C = (16.0, 12.0)
|
||||||
whole = Jet(Sector(C, .1, math.pi, 3.), Radial(C, inward=False), amplitude=0.)
|
whole = Jet(Sector(C, .1, math.pi, 3.), Radial(C, inward=False), amplitude=0.)
|
||||||
backend, _ = mocked([circle()], [whole], FINE_BLOCKS[1])
|
backend, _ = mocked([circle()], [whole], FINE_BLOCKS[1])
|
||||||
_, jet_id = backend._device_boundaries[1].jet_cache
|
_, jet_id = backend._boundaries[1].jet_cache
|
||||||
self.assertGreater(len(np.unique(backend._device_boundaries[1].tables["block"][jet_id == 0])), 1)
|
self.assertGreater(len(np.unique(backend._boundaries[1].tables["block"][jet_id == 0])), 1)
|
||||||
backend.close()
|
backend.close()
|
||||||
|
|
||||||
def test_fine_ownership_uses_remapped_physical_owner(self):
|
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)
|
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.)
|
whole = Jet(Sector(C, .1, math.pi, 3.), Radial(C, inward=False), amplitude=0.)
|
||||||
backend, _ = mocked([source], [whole], FINE_BLOCKS[1])
|
backend, _ = mocked([source], [whole], FINE_BLOCKS[1])
|
||||||
_, jet_id = backend._device_boundaries[1].jet_cache
|
_, jet_id = backend._boundaries[1].jet_cache
|
||||||
owners = backend.boundaries[1].curved_host[3][jet_id == 0]
|
owners = backend.boundaries[1].curved_host[3][jet_id == 0]
|
||||||
self.assertEqual(np.unique(owners).tolist(), [0])
|
self.assertEqual(np.unique(owners).tolist(), [0])
|
||||||
self.assertGreater(len(np.unique(backend._device_boundaries[1].tables["block"][jet_id == 0])), 1)
|
self.assertGreater(len(np.unique(backend._boundaries[1].tables["block"][jet_id == 0])), 1)
|
||||||
self.assertEqual(source.obj_id, 99)
|
self.assertEqual(source.obj_id, 99)
|
||||||
backend.close()
|
backend.close()
|
||||||
cfg = config("periodic")
|
cfg = config("periodic")
|
||||||
@@ -97,8 +97,8 @@ class JetStepperCPU(unittest.TestCase):
|
|||||||
def test_blocked_jet_velocity_matches_legacy_rows(self):
|
def test_blocked_jet_velocity_matches_legacy_rows(self):
|
||||||
legacy, _ = mocked([circle(actuation(SPEED))], (), FINE_BLOCKS[1])
|
legacy, _ = mocked([circle(actuation(SPEED))], (), FINE_BLOCKS[1])
|
||||||
driven, _ = mocked([circle()], [jet()], FINE_BLOCKS[1])
|
driven, _ = mocked([circle()], [jet()], FINE_BLOCKS[1])
|
||||||
a = legacy._device_boundaries[1].tables["surface_velocity"]
|
a = legacy._boundaries[1].tables["surface_velocity"]
|
||||||
b = driven._device_boundaries[1].tables["surface_velocity"]
|
b = driven._boundaries[1].tables["surface_velocity"]
|
||||||
diff = np.abs(a - b).max()
|
diff = np.abs(a - b).max()
|
||||||
print(f"fine jet vs legacy velocity rows: max |dU| = {diff:.3e} (exact={diff == 0})")
|
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)
|
self.assertLessEqual(diff, 4 * np.finfo(np.float32).eps * SPEED)
|
||||||
@@ -106,8 +106,8 @@ class JetStepperCPU(unittest.TestCase):
|
|||||||
|
|
||||||
def test_blocked_amplitude_update_rewrites_only_velocity_column(self):
|
def test_blocked_amplitude_update_rewrites_only_velocity_column(self):
|
||||||
backend, htod = mocked([circle()], [jet()], FINE_BLOCKS[1])
|
backend, htod = mocked([circle()], [jet()], FINE_BLOCKS[1])
|
||||||
before = backend._device_boundaries[1].tables["surface_velocity"].copy()
|
before = backend._boundaries[1].tables["surface_velocity"].copy()
|
||||||
device = backend._device_boundaries[1].device["surface_velocity"]
|
device = backend._boundaries[1].device["surface_velocity"]
|
||||||
with patch.object(cuda, "memcpy_htod", htod):
|
with patch.object(cuda, "memcpy_htod", htod):
|
||||||
htod.reset_mock()
|
htod.reset_mock()
|
||||||
backend.set_jet_amplitude(0, -SPEED, level=1)
|
backend.set_jet_amplitude(0, -SPEED, level=1)
|
||||||
@@ -115,7 +115,7 @@ class JetStepperCPU(unittest.TestCase):
|
|||||||
self.assertEqual(htod.call_count, 1)
|
self.assertEqual(htod.call_count, 1)
|
||||||
self.assertIs(ptr, device)
|
self.assertIs(ptr, device)
|
||||||
np.testing.assert_allclose(sent, -before, rtol=1e-6, atol=0)
|
np.testing.assert_allclose(sent, -before, rtol=1e-6, atol=0)
|
||||||
np.testing.assert_array_equal(backend._device_boundaries[1].tables["surface_velocity"], sent)
|
np.testing.assert_array_equal(backend._boundaries[1].tables["surface_velocity"], sent)
|
||||||
with self.assertRaises(ValueError):
|
with self.assertRaises(ValueError):
|
||||||
backend.set_jet_amplitude(0, 1.5, level=1)
|
backend.set_jet_amplitude(0, 1.5, level=1)
|
||||||
self.assertEqual(htod.call_count, 1)
|
self.assertEqual(htod.call_count, 1)
|
||||||
@@ -133,7 +133,7 @@ class JetStepperCPU(unittest.TestCase):
|
|||||||
with self.assertRaises(IndexError):
|
with self.assertRaises(IndexError):
|
||||||
backend.set_jet_amplitude(1, -SPEED, level=1)
|
backend.set_jet_amplitude(1, -SPEED, level=1)
|
||||||
self.assertEqual(backend._status, status)
|
self.assertEqual(backend._status, status)
|
||||||
before = backend._device_boundaries[1].tables["surface_velocity"].copy()
|
before = backend._boundaries[1].tables["surface_velocity"].copy()
|
||||||
boom = Mock(side_effect=RuntimeError("copy failed"))
|
boom = Mock(side_effect=RuntimeError("copy failed"))
|
||||||
with patch.object(cuda, "memcpy_htod", boom), patch.object(cuda, "memcpy_htod_async", boom), \
|
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)):
|
patch.object(cuda, "pagelocked_empty", side_effect=lambda shape, dtype: np.empty(shape, dtype)):
|
||||||
@@ -141,7 +141,7 @@ class JetStepperCPU(unittest.TestCase):
|
|||||||
backend.set_jet_amplitude(0, -SPEED, level=1)
|
backend.set_jet_amplitude(0, -SPEED, level=1)
|
||||||
self.assertEqual(backend._status, "failed", fine_block)
|
self.assertEqual(backend._status, "failed", fine_block)
|
||||||
self.assertEqual(backend.jets[1][0].amplitude, SPEED)
|
self.assertEqual(backend.jets[1][0].amplitude, SPEED)
|
||||||
np.testing.assert_array_equal(backend._device_boundaries[1].tables["surface_velocity"], before)
|
np.testing.assert_array_equal(backend._boundaries[1].tables["surface_velocity"], before)
|
||||||
with self.assertRaisesRegex(RuntimeError, "failed"):
|
with self.assertRaisesRegex(RuntimeError, "failed"):
|
||||||
backend.step(1)
|
backend.step(1)
|
||||||
with self.assertRaisesRegex(RuntimeError, "failed"):
|
with self.assertRaisesRegex(RuntimeError, "failed"):
|
||||||
|
|||||||
@@ -102,7 +102,7 @@ class Ledger:
|
|||||||
|
|
||||||
def step(self):
|
def step(self):
|
||||||
run=self.run;self.samples.clear();run.step()
|
run=self.run;self.samples.clear();run.step()
|
||||||
b=run._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)
|
self.test.assertEqual(len(self.samples),2**self.index)
|
||||||
expected=np.zeros_like(self.samples[0],dtype=np.float64)
|
expected=np.zeros_like(self.samples[0],dtype=np.float64)
|
||||||
weights=np.array([h*h,h*h,h*h*h,h*h])
|
weights=np.array([h*h,h*h,h*h*h,h*h])
|
||||||
@@ -140,7 +140,7 @@ class ShiftedJetsGPU(unittest.TestCase):
|
|||||||
level_objects=[(),(obj,)],level_jets=[(),(whole,)],observe_raw_exchange=True) as run:
|
level_objects=[(),(obj,)],level_jets=[(),(whole,)],observe_raw_exchange=True) as run:
|
||||||
run._initialize_raw(initial(run,rho,u=.013,v=0.))
|
run._initialize_raw(initial(run,rho,u=.013,v=0.))
|
||||||
expected=run._download_raw();ledger=Ledger(self,run,1)
|
expected=run._download_raw();ledger=Ledger(self,run,1)
|
||||||
self.assertTrue(run._device_boundaries[1].tables['aperture'].all())
|
self.assertTrue(run._boundaries[1].tables['aperture'].all())
|
||||||
for _ in range(4):ledger.step()
|
for _ in range(4):ledger.step()
|
||||||
error=max(float(abs(a[:,leaf]-b[:,leaf]).max()) for a,b,leaf in
|
error=max(float(abs(a[:,leaf]-b[:,leaf]).max()) for a,b,leaf in
|
||||||
zip(run._download_raw().values,expected.values,run.leaf))
|
zip(run._download_raw().values,expected.values,run.leaf))
|
||||||
@@ -201,7 +201,7 @@ class ShiftedJetsGPU(unittest.TestCase):
|
|||||||
self.exact_raw(before,old._download_raw())
|
self.exact_raw(before,old._download_raw())
|
||||||
self.exact_raw(_StoredPopulations(expected.g,True),new._download_raw())
|
self.exact_raw(_StoredPopulations(expected.g,True),new._download_raw())
|
||||||
self.assertEqual(new.jets[2][0].amplitude,-.001)
|
self.assertEqual(new.jets[2][0].amplitude,-.001)
|
||||||
np.testing.assert_array_equal(new._device_boundaries[2].tables['surface_velocity'],old._device_boundaries[2].tables['surface_velocity'])
|
np.testing.assert_array_equal(new._boundaries[2].tables['surface_velocity'],old._boundaries[2].tables['surface_velocity'])
|
||||||
newledger=Ledger(self,new,2);ledger.step();newledger.step()
|
newledger=Ledger(self,new,2);ledger.step();newledger.step()
|
||||||
self.assertEqual((old.step_count,new.step_count),(count+1,count+1))
|
self.assertEqual((old.step_count,new.step_count),(count+1,count+1))
|
||||||
ledger.report('L2 source');newledger.report('L2 remapped')
|
ledger.report('L2 source');newledger.report('L2 remapped')
|
||||||
|
|||||||
Reference in new issue
Block a user