Compare commits
9
Commits
| Author | SHA1 | Date | |
|---|---|---|---|
|
|
c5a3cdb757 | ||
|
|
a4fc46d849 | ||
|
|
65339ab54e | ||
|
|
ae766f7112 | ||
|
|
5f69ebaa90 | ||
|
|
9d0588dcad | ||
|
|
7b652b2423 | ||
|
|
7482d95a42 | ||
|
|
55067652e4 |
@@ -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
|
||||
```
|
||||
|
||||
@@ -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
|
||||
|
||||
@@ -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]
|
||||
|
||||
@@ -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
|
||||
|
||||
@@ -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)
|
||||
|
||||
@@ -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.
|
||||
|
||||
|
||||
@@ -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,
|
||||
|
||||
@@ -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)
|
||||
|
||||
@@ -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
|
||||
@@ -294,23 +320,17 @@ class FixedPatchStepper:
|
||||
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)
|
||||
# The one boundary/Jet record per level; host geometry and Jets are views of it.
|
||||
self._device_boundaries = [root] + [
|
||||
_prepare_blocked_boundary(self._boundaries[level], self._layouts[level],
|
||||
_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
|
||||
self._leaf = _fluid_leaves(self.hierarchy, self._device_boundaries)
|
||||
self._state = [mask.copy() for mask in self._leaf]
|
||||
if self._channel:
|
||||
self._leaf[0] &= channel_masks((cfg.ny, cfg.nx))[0]
|
||||
@@ -556,7 +576,6 @@ 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])):
|
||||
raise ValueError("remap_finest requires a periodic three-level hierarchy with a "
|
||||
"level-2 set and no level-0/1 bodies")
|
||||
@@ -567,8 +586,7 @@ class FixedPatchStepper:
|
||||
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")
|
||||
@@ -622,7 +640,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
|
||||
@@ -634,26 +652,23 @@ class FixedPatchStepper:
|
||||
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])
|
||||
# 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
|
||||
@@ -901,9 +916,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._device_boundaries)
|
||||
|
||||
@property
|
||||
def jets(self):
|
||||
return tuple(() if b is None else b.jets for b in self._device_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)
|
||||
|
||||
@@ -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(){
|
||||
|
||||
@@ -103,7 +103,7 @@ 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]
|
||||
return boundary.tables["surface_velocity"].copy(), boundary.jets, run.jets[2]
|
||||
|
||||
|
||||
def test_stepper_invalid_input_writes_nothing_and_stays_ready():
|
||||
@@ -117,10 +117,10 @@ def test_stepper_invalid_input_writes_nothing_and_stays_ready():
|
||||
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)
|
||||
assert (run._device_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()
|
||||
|
||||
@@ -137,7 +137,7 @@ def test_stepper_device_failure_keeps_host_and_blocks_reuse(where):
|
||||
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)
|
||||
assert (run._device_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
|
||||
|
||||
@@ -94,13 +94,13 @@ 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'],
|
||||
run._boundaries[2].tables['surface_velocity'])
|
||||
self.assertEqual(tuple(j.amplitude for j in clone._device_boundaries[2].jets), (.017, -.011))
|
||||
np.testing.assert_array_equal(clone._device_boundaries[2].tables['surface_velocity'],
|
||||
run._device_boundaries[2].tables['surface_velocity'])
|
||||
clone.set_jet_amplitude(1, .005, level=2)
|
||||
self.assertEqual(run._boundaries[2].jets[1].amplitude, -.011)
|
||||
self.assertEqual(run._device_boundaries[2].jets[1].amplitude, -.011)
|
||||
clone.step()
|
||||
self.assertEqual((run.step_count, clone.step_count), (7, 8))
|
||||
for got, want in zip(run.download(), before):
|
||||
@@ -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()
|
||||
|
||||
@@ -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)
|
||||
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))
|
||||
clone=run.remap_finest(B,objects=objects)
|
||||
self.assertEqual(clone._device_boundaries[2].jets[0].amplitude,.017)
|
||||
np.testing.assert_array_equal(clone._device_boundaries[2].tables['surface_velocity'],
|
||||
run._device_boundaries[2].tables['surface_velocity'])
|
||||
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,22 +73,22 @@ 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)
|
||||
self.assertFalse(clone._boundaries[2].actuated)
|
||||
self.assertNotIn('surface_velocity',clone._boundaries[2].tables)
|
||||
clone=run.remap_finest(A,objects=objects)
|
||||
self.assertFalse(clone._device_boundaries[2].actuated)
|
||||
self.assertNotIn('surface_velocity',clone._device_boundaries[2].tables)
|
||||
for a,b in zip(run.download(),clone.download()): np.testing.assert_array_equal(a,b)
|
||||
clone.close();run.close()
|
||||
|
||||
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()
|
||||
|
||||
@@ -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)
|
||||
|
||||
@@ -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)
|
||||
|
||||
@@ -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())
|
||||
|
||||
@@ -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()
|
||||
|
||||
@@ -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._device_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=()
|
||||
('_device_boundaries', [*records[:2], SimpleNamespace(host=object(), jets=())]),
|
||||
('_device_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
|
||||
|
||||
@@ -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)
|
||||
clone=run.remap_finest(B,objects=objects)
|
||||
self.assertEqual(clone.jets[2][0].amplitude,-jet().amplitude)
|
||||
np.testing.assert_array_equal(clone._device_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()
|
||||
|
||||
@@ -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()
|
||||
|
||||
@@ -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.assertEqual(clone._device_boundaries[2].jets[0].amplitude,.017)
|
||||
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()
|
||||
|
||||
@@ -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)))
|
||||
|
||||
|
||||
@@ -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()
|
||||
|
||||
@@ -66,31 +66,31 @@ class ExchangeWindowCPU(unittest.TestCase):
|
||||
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:
|
||||
@@ -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):
|
||||
@@ -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
|
||||
@@ -231,28 +231,28 @@ class ExchangeWindowShiftedMRTGPU(unittest.TestCase):
|
||||
cuda.memcpy_dtoh(rows[1], split._device_boundaries[1].window_integral)
|
||||
np.testing.assert_array_equal(rows[0].view(np.uint64), rows[1].view(np.uint64))
|
||||
|
||||
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,8 +56,8 @@ def make(context):
|
||||
|
||||
def topology(run, rectangles, objects):
|
||||
return RemapTopology.from_config(
|
||||
CFG, run.patch, rectangles, child_objects=objects,
|
||||
child_jets=run._boundaries[2].jets,
|
||||
CFG, run.hierarchy.pairs[0], rectangles, child_objects=objects,
|
||||
child_jets=run._device_boundaries[2].jets,
|
||||
)
|
||||
|
||||
|
||||
@@ -74,8 +74,8 @@ class MultiplePinnedRemapGPU(unittest.TestCase):
|
||||
def change_amplitudes(self, run, amplitudes):
|
||||
# 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]
|
||||
boundary = run._device_boundaries[2]
|
||||
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,11 +87,11 @@ 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]
|
||||
boundary = run._device_boundaries[2]
|
||||
self.assertEqual(tuple(jet.amplitude for jet in boundary.jets), amplitudes)
|
||||
exact(self, boundary.tables['surface_velocity'], mesh.surface_velocity,
|
||||
'latest physical per-link velocities')
|
||||
@@ -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,12 +36,12 @@ 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')
|
||||
self.assertEqual(after._boundaries[2].jets[0].amplitude,.017)
|
||||
exact(self,after._device_boundaries[2].tables['surface_velocity'],
|
||||
before._device_boundaries[2].tables['surface_velocity'],'current per-link velocities')
|
||||
self.assertEqual(after._device_boundaries[2].jets[0].amplitude,.017)
|
||||
|
||||
def test_noop_next_step_bitwise(self):
|
||||
with CudaContext(0) as ctx, ExitStack() as owned:
|
||||
@@ -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')
|
||||
@@ -63,14 +63,14 @@ class PinnedBodyRemapGPU(unittest.TestCase):
|
||||
current.set_jet_amplitude(0,.017,level=2)
|
||||
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,
|
||||
jets=previous._device_boundaries[2].jets
|
||||
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)
|
||||
@@ -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):
|
||||
@@ -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)
|
||||
@@ -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):
|
||||
@@ -213,7 +213,7 @@ 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"])
|
||||
self.assertIsNotNone(backend._device_boundaries[1].tables["surface_velocity"])
|
||||
del backend._device_boundaries[1].device["surface_velocity"]
|
||||
with self.assertRaisesRegex(RuntimeError, "Surface velocity device buffer is missing"):
|
||||
backend.step(1)
|
||||
|
||||
@@ -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._device_boundaries[1].jet_cache
|
||||
self.assertGreater(len(np.unique(backend._device_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._device_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._device_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._device_boundaries[1].tables["surface_velocity"]
|
||||
b = driven._device_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()
|
||||
before = backend._device_boundaries[1].tables["surface_velocity"].copy()
|
||||
device = backend._device_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._device_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._device_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._device_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()))
|
||||
@@ -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
|
||||
@@ -143,7 +143,7 @@ class ShiftedJetsGPU(unittest.TestCase):
|
||||
self.assertTrue(run._device_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,14 +193,14 @@ 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)
|
||||
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'])
|
||||
newledger=Ledger(self,new,2);ledger.step();newledger.step()
|
||||
self.assertEqual((old.step_count,new.step_count),(count+1,count+1))
|
||||
|
||||
@@ -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)
|
||||
|
||||
@@ -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))
|
||||
|
||||
@@ -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)
|
||||
|
||||
|
||||
Reference in New Issue
Block a user