test: retire dense-route tests, require tiles everywhere (D2)

Retired with the dense refinement route (results recorded at 1573942 and
2469af7 or named as accepted losses):
- test_refinement_blocked_base: BlockedBaseGPU exact parity (3 tests),
  DenseBlockedRouteParity (6 tests incl. MRT and two bodies),
  DenseBlockedRouteCPU (capability matrix, ledger comparator);
- test_refinement_blocked_fine: BlockedFineGPU dense-vs-blocked (3 tests);
- test_refinement_mrt_gpu: test_dense_versus_blocked_mrt_two_level;
- test_raw_exchange_gpu: test_dense_fine_storage_reference (private
  link_exchange cut-link oracle, accepted loss);
- test_refinement_gpu: test_bounded_equivalence_and_ownership (dense
  CELERIS_REFINEMENT_AUDIT ownership counters, accepted loss); its physics
  stays in test_bounded_equivalence_blocked; eight_steps loses the audit
  branch and the misleading 'ownership/counts PASS' label;
- test_fine_surface_gpu: test_fine_surface_four_hops (dense per-hop flux
  ledger, CurvedLinkSoA reflection oracle, audit counts: accepted losses;
  raw exchange covers the boundary part only, not full flux closure) and the
  dense SoA fault tests (missing surface buffer, close frees SoA); blocked
  equivalents remain in test_refinement_blocked_fine/_base;
- test_raw_exchange.test_dense_observation_rejected_before_device_access.
New TileContractCPU: level_blocks omitted/None/wrong length/None entry/
wrong arity/non-integer/bool/non-positive/non-tuple/not dividing all raise
before CUDA; whole-level tiles reach CUDA. Rejection tests and fixtures pass
tiles; launch-order and kernel-text tests follow the blocked root;
examples/fixed_periodic_patch uses whole-level tiles (+ host tile test).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Claude
2026-10-27 15:12:20 +00:00
co-authored by Claude Opus 5.5
parent e18f7abe2b
commit 2b10b01c49
14 changed files with 84 additions and 777 deletions
+5
View File
@@ -46,6 +46,11 @@ class ExampleTiles(unittest.TestCase):
FixedPatchStepper(cfg, context=forbidden(), coarse_boundary="slot_reservoir",
refinements=[bounds], level_blocks=[tuple(b) for b in blocks])
def test_fixed_periodic_patch_tiles(self):
with self.assertRaisesRegex(AssertionError, "CUDA reached"):
FixedPatchStepper(LBMConfig(nx=16, ny=16, viscosity=.1, velocity=1e-3), context=forbidden(),
refinements=[(4, 4, 12, 12)], level_blocks=[(16, 16), (24, 24)])
def test_fixed_patch_circle_default_tiles(self):
bounds = literal("fixed_patch_circle.py", "BOUNDS")
nx, ny = literal("fixed_patch_circle.py", "NX"), literal("fixed_patch_circle.py", "NY")
+2 -4
View File
@@ -146,7 +146,6 @@ class ProjectedPairDispatch(unittest.TestCase):
for cfg, changes in [(CFG, {'explosion':'unknown'}),
(CFG, {'coarse_boundary':'slot_reservoir'}),
(CFG, {'coarse_boundary':'channel'}),
(CFG, {'level_blocks':[None,None], 'refinements':[(6,6,26,26)]}),
(CFG, {'level_objects':[(),(object(),)]}),
(CFG, {'level_jets':[(),(object(),)]}),
(replace(CFG,collision='MRT'), {}),
@@ -155,9 +154,8 @@ class ProjectedPairDispatch(unittest.TestCase):
host.FixedPatchStepper(cfg, **dict(options, **changes))
context.activate.assert_not_called()
host.cuda.mem_alloc.assert_not_called()
# The legacy constant path still allows MRT and dense storage.
for options in ({'cfg':replace(CFG,collision='MRT')},
{'level_blocks':[None,None], 'refinements':[(6,6,26,26)]}):
# The constant path still allows MRT.
for options in ({'cfg':replace(CFG,collision='MRT')},):
run, _, _, _, _ = backend(**options)
run.close()
-9
View File
@@ -153,15 +153,6 @@ class RawExchange(unittest.TestCase):
run._stream.synchronize.side_effect = None
run.close()
def test_dense_observation_rejected_before_device_access(self):
host = host_module()
ctx = Mock()
with self.assertRaisesRegex(ValueError, 'blocked body'):
host.FixedPatchStepper(LBMConfig(nx=32, ny=32, viscosity=.1), context=ctx,
refinements=[(4, 4, 28, 28)],
level_objects=[(), (body(radius=2.15),)], observe_raw_exchange=True)
ctx.activate.assert_not_called()
def test_all_depth_weights(self):
for level in (1, 2, 3, 7):
self._check_weights_split_and_window(level)
+2 -5
View File
@@ -27,23 +27,20 @@ class RefinementMRTHost(unittest.TestCase):
block_collision = (KERNELS / "operators/block_collision.cuh").read_text()
self.assertIn('#include "operators/block_collision.cuh"', blocks)
self.assertEqual(blocks.count("block_collide(f);"), 1)
refinement = (KERNELS / "refinement.cu").read_text() + (KERNELS / "step/patch_step.cu").read_text()
for text in ("compute_feq_from_momentum(drho, jx, jy, ux, uy, feq);",
"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_no_force(f, feq, d_params.omega);", refinement)
for source in (blocks, refinement):
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)
self.assertIn("USE_DDF_SHIFTING ||", refinement) # dense/mixed shifted stays rejected
# 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, refinement):
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)
+2 -2
View File
@@ -125,11 +125,11 @@ class ShiftedStorage(unittest.TestCase):
with self.assertRaisesRegex(ValueError,'mixed'): run._download_raw()
finally: run.close()
def test_dense_or_mixed_shifted_rejected_before_compilation(self):
def test_missing_tiles_rejected_before_compilation(self):
module=backend.host_module()
cfg=replace(backend.CFG,ddf_shifting=True)
for blocks in ([None,None],[(4,4),None],[None,(4,4)]):
with self.assertRaisesRegex(ValueError,'all-blocked|blocked refined levels below'):
with self.assertRaisesRegex(ValueError,'must be a tile'):
module.FixedPatchStepper(cfg,context=None,refinements=[(4,4,20,28)],
level_blocks=blocks)
module.compiler.build_module.assert_not_called()
@@ -1,13 +1,17 @@
"""Two coarse/four fine steps: isolated fine circle, no flux-accuracy claim.
"""Two coarse/four fine steps: isolated fine circle on blocked storage, no flux-accuracy claim.
Run only through the bounded GPU runner after a fresh resource check. The test
hooks read actual kernel inputs/outputs; they never change the stepping order.
Run only through the bounded GPU runner after a fresh resource check. The fine
step hook reads the actual output buffer; it never changes the stepping order.
Retired with the dense refinement route (accepted losses, see the D2 commit):
the dense-kernel per-hop flux ledger (fine/coarse transfer + boundary budget),
the CurvedLinkSoA reflection oracle and the CELERIS_REFINEMENT_AUDIT collision
counts. The blocked test keeps the four-hop comparison with the uniform solver,
the raw-exchange mass sign and passive == zero actuation; raw exchange covers
the boundary part only, not full population flux closure.
"""
import importlib.util
from contextlib import nullcontext
from pathlib import Path
from threading import RLock
from types import SimpleNamespace
from unittest.mock import Mock, patch
import numpy as np
@@ -18,7 +22,6 @@ import CelerisLab
from CelerisLab.body.geometry.circle import CircleGeometry
from CelerisLab.body.objects import SimObject
from CelerisLab.config import LBMConfig
from CelerisLab.cuda import compiler
from CelerisLab.cuda.context import CudaContext
from CelerisLab.lbm import refinement_stepper as implementation
from CelerisLab.lbm.descriptors import FLUID, SOLID
@@ -27,10 +30,8 @@ from CelerisLab.lbm.descriptors import FLUID, SOLID
ROOT = Path(__file__).resolve().parents[3]
CX = np.array([0, 1, -1, 0, 0, 1, -1, 1, -1])
CY = np.array([0, 0, 0, 1, -1, 1, -1, -1, 1])
OPP = np.array([0, 2, 1, 4, 3, 6, 5, 8, 7])
W = np.array([4/9]+[1/9]*4+[1/36]*4)
MOMENTS = np.stack((np.ones(9), CX, CY))
GEOMETRY = ("fluid_idx", "dir", "q", "body_id", "rx", "ry", "rz", "fallback_class")
BOUNDS = (4, 4, 20, 20)
@@ -53,197 +54,6 @@ def object_spec(speed):
radius=2.15, surface_actuation=surface(speed))
def audit_compiler(original):
def compile_private(*args, **kwargs):
entry = Path(kwargs["build_dir"])/"fine_surface_audit.cu"
entry.write_text('#define CELERIS_REFINEMENT_AUDIT\n#include "refinement.cu"\n')
kwargs["source"] = entry.name
return original(*args, **kwargs)
return compile_private
def reflection_oracle(before, cl, speed, *, verify_projection=True):
"""FP64 BFL from frozen spatial populations, independent arc projection."""
flat = before.astype(float).reshape(9, -1)
radius = np.hypot(cl.rx.astype(float), cl.ry.astype(float))
delta = (np.arctan2(cl.ry.astype(float), cl.rx.astype(float))-.7+np.pi) % (2*np.pi)-np.pi
amplitude = 0 if speed is None else speed
strength = amplitude*np.maximum(0, 1-(delta/.85)**2)**2
velocity = strength[:, None]*np.column_stack((cl.rx, cl.ry))/radius[:, None]
if verify_projection and cl.surface_velocity is not None:
np.testing.assert_allclose(cl.surface_velocity, velocity, atol=3e-9, rtol=0)
incoming = []
for i, (k, d, q, fallback) in enumerate(zip(cl.fluid_idx, cl.dir, cl.q, cl.fallback_class)):
k, d, q = int(k), int(d), float(q)
rho = flat[:, k].sum()
dot = CX[d]*velocity[i, 0]+CY[d]*velocity[i, 1]
if fallback:
value = flat[d, k]-6*W[d]*rho*dot
elif q < .5:
value = 2*q*flat[d, k]+(1-2*q)*flat[d, k-CX[d]-40*CY[d]]-6*W[d]*rho*dot
else:
value = flat[d, k]/(2*q)+(1-1/(2*q))*flat[OPP[d], k]-3*W[d]*rho*dot/q
incoming.append(value)
return np.array(incoming)
def fine_flux(before, staged, active, solid):
"""Population ledger across the fine leaf boundary: solid versus halo.
Internal fluid-to-fluid streaming cancels pairwise; solid and halo packets
are counted in their actual direction, before any collision.
"""
boundary, transfer = np.zeros(9), np.zeros(9)
magnitude = 0.0
for d, (cx, cy) in enumerate(zip(CX, CY)):
incoming = active & ~np.roll(active, (cy, cx), (0, 1))
outgoing = active & ~np.roll(active, (-cy, -cx), (0, 1))
from_solid = np.roll(solid, (cy, cx), (0, 1))
into_solid = np.roll(solid, (-cy, -cx), (0, 1))
received = np.roll(staged[d], (cy, cx), (0, 1)).astype(float)
sent = before[d].astype(float)
for vector, inc, out in ((boundary, incoming & from_solid, outgoing & into_solid),
(transfer, incoming & ~from_solid, outgoing & ~into_solid)):
assert np.isfinite(received[inc]).all() and np.isfinite(sent[out]).all()
vector[d] = received[inc].sum()-sent[out].sum()
magnitude += np.abs(received[inc]).sum()+np.abs(sent[out]).sum()
return boundary, transfer, magnitude
def coarse_flux(before, fine, active):
"""Derive every replacement from physical midpoints, not production tables."""
x0, y0, x1, y1 = BOUNDS
incoming = np.stack([np.roll(before[d], (CY[d], CX[d]), (0, 1)) for d in range(9)]).astype(float)
magnitude = 0.0
for y, x in np.argwhere(active):
for d in range(1, 9):
if x0 <= x+.5-CX[d]/2 <= x1 and y0 <= y+.5-CY[d]/2 <= y1:
a, b = 4+2*(x-x0), 4+2*(y-y0)
children = fine[d, b:b+2, a:a+2].astype(float)
assert children.shape == (2, 2) and np.isfinite(children).all()
incoming[d, y, x] = children.sum()/4
magnitude += np.abs(children).sum()/4
assert np.isfinite(incoming[:, active]).all()
# This retains the eight grazing-corner replacements of uncovered donors.
flux = (incoming[:, active]-before[:, active].astype(float)).sum(axis=1)
magnitude += np.abs(incoming[:, active]).sum()+np.abs(before[:, active]).sum(dtype=float)
return flux, magnitude
def run_patch(speed, *, audit=True):
cfg = LBMConfig(nx=24, ny=24, viscosity=.05, velocity=0, rho=1)
coarse_active = np.ones((24, 24), bool)
coarse_active[4:20, 4:20] = False
source = object_spec(speed)
state = vars(source.state).copy()
with CudaContext(0) as context:
build = (patch.object(compiler, "compile_kernel", audit_compiler(compiler.compile_kernel))
if audit else nullcontext())
with build:
backend = implementation.FixedPatchStepper(cfg, context=context, refinements=[BOUNDS],
level_objects=[(), (source,)])
with backend, context.activate():
c, f = backend._levels
cl = backend._fine_curved
geometry = {name: getattr(cl, name).copy() for name in GEOMETRY}
flags = backend.fine_boundary.flags.reshape(40, 40)
solid = (flags & SOLID) != 0
active = np.zeros((40, 40), bool)
active[4:-4, 4:-4] = True
active &= (flags & FLUID) != 0
active &= ~solid
local = backend.fine_boundary.objects[0]
assert local.center == (20., 20.) and local.radius == 4.3
assert source.center == (12.25, 12.25) and source.radius == 2.15
yy, xx = np.mgrid[:40, :40]
roi = active & ((xx-20)**2+(yy-20)**2 <= 9**2)
backend.initialize(rest((24, 24)), rest((40, 40)))
if audit:
for level in backend._levels:
for name in ("collision_counts", "explosion_counts", "incoming_counts"):
pointer, size = level.module.get_global(name)
cuda.memset_d8(pointer, 0, size)
initial = backend.download()
initial_pop = initial[0][:, coarse_active].sum(axis=1, dtype=float)
initial_pop += .25*initial[1][:, active].sum(axis=1, dtype=float)
boundary, fine_transfer, coarse_transfer = np.zeros(9), np.zeros(9), np.zeros(9)
budget_sum, reflection_error, drive_mass = 0., 0., 0.
traces = []
original_curved, original_fine, original_coarse = backend._curved, f.kernel, c.kernel
def curved(*args, **kwargs):
nonlocal boundary, fine_transfer, budget_sum, reflection_error, drive_mass
backend._stream.synchronize()
before = read(args[0], f.shape)
expected = reflection_oracle(before, cl, speed)
original_curved(*args, **kwargs)
backend._stream.synchronize()
staged = read(args[0], f.shape)
d, k = cl.dir.astype(int), cl.fluid_idx.astype(int)
values = staged.reshape(9, -1)[OPP[d], k+CX[d]+40*CY[d]]
passive = reflection_oracle(before, cl, None, verify_projection=False)
drive_mass += .25*float(np.sum(values-passive))
error = float(np.max(np.abs(values-expected)))
reflection_error = max(reflection_error, error)
assert error <= 1e-7
np.testing.assert_array_equal(staged[:, active], before[:, active])
b, transfer, magnitude = fine_flux(before, staged, active, solid)
boundary += .25*b
fine_transfer += .25*transfer
budget_sum += .25*magnitude
def fine(*args, **kwargs):
nonlocal budget_sum
backend._stream.synchronize()
staged = read(args[0], f.shape)
original_fine(*args, **kwargs)
backend._stream.synchronize()
advanced = read(args[1], f.shape)
assert np.isfinite(advanced[:, active]).all()
assert np.isnan(advanced[:, solid]).all()
# Include actual collision operands in the FP32 rounding budget.
budget_sum += .25*np.abs(staged[:, active]).sum(dtype=float)
traces.append(advanced[:, roi].copy())
def coarse(*args, **kwargs):
nonlocal coarse_transfer, budget_sum
backend._stream.synchronize()
before, fine_values = read(args[0], c.shape), read(args[2], f.shape)
transfer, magnitude = coarse_flux(before, fine_values, coarse_active)
coarse_transfer += transfer
budget_sum += magnitude
return original_coarse(*args, **kwargs)
with patch.object(backend, "_curved", curved), patch.object(f, "kernel", fine), \
patch.object(c, "kernel", coarse):
backend.step(2)
actual = backend.download()
population = actual[0][:, coarse_active].sum(axis=1, dtype=float)
population += .25*actual[1][:, active].sum(axis=1, dtype=float)
ledger = boundary+fine_transfer+coarse_transfer
error = MOMENTS@(population-initial_pop-ledger)
budget = 8*np.finfo(np.float32).eps*budget_sum
assert np.isfinite(budget_sum) and budget_sum > 0
assert np.max(np.abs(error)) <= budget
if speed:
assert speed*drive_mass > 0
assert len(traces) == 4 and backend.step_count == 2
if audit:
for level, mask, count in ((c, coarse_active, 2), (f, active, 4)):
pointer, size = level.module.get_global("collision_counts")
counts = read(pointer, mask.shape, np.uint32)
assert size == counts.nbytes
np.testing.assert_array_equal(counts, count*mask)
for name, values in geometry.items():
np.testing.assert_array_equal(getattr(cl, name), values)
assert vars(source.state) == state and source.geometry.center == source.center
print(f"speed={speed} audit={audit} BFL={reflection_error:.9g} ledger={error.tolist()} "
f"budget={budget:.9g} boundary={list(MOMENTS@boundary)} "
f"transfer={list(MOMENTS@(fine_transfer+coarse_transfer))} "
f"drive_mass={drive_mass:+.9g}; counts={'PASS' if audit else 'normal build'}", flush=True)
return traces, roi, actual
def uniform_reference(speed, roi):
# Reuse the existing explicit config helper without importing unrelated tests.
spec = importlib.util.spec_from_file_location("fine_surface_cases", ROOT/"tests/integration/_cases.py")
@@ -261,24 +71,6 @@ def uniform_reference(speed, roi):
return traces
@pytest.mark.parametrize("speed", [0., .002, -.002])
def test_fine_surface_four_hops(speed):
assert Path(CelerisLab.__file__).resolve() == ROOT/"src/CelerisLab/__init__.py"
assert Path(compiler._KERNEL_DIR).resolve() == ROOT/"src/CelerisLab/lbm/kernels"
assert Path(implementation.__file__).resolve() == ROOT/"src/CelerisLab/lbm/refinement_stepper.py"
actual, roi, final = run_patch(speed)
reference = uniform_reference(speed, roi)
for got, wanted in zip(actual, reference):
np.testing.assert_allclose(got, wanted, atol=3e-6, rtol=0)
np.testing.assert_allclose(MOMENTS@got.astype(float), MOMENTS@wanted.astype(float),
atol=2e-5, rtol=0)
if speed == 0:
passive, passive_roi, passive_final = run_patch(None, audit=False)
np.testing.assert_array_equal(passive_roi, roi)
for got, wanted in zip(actual+list(final), passive+list(passive_final)):
np.testing.assert_array_equal(got, wanted)
def run_blocked_patch(speed):
"""All-blocked run (whole-level tiles) of the same case; fine ROI after each fine hop."""
cfg = LBMConfig(nx=24, ny=24, viscosity=.05, velocity=0, rho=1)
@@ -334,46 +126,6 @@ def test_fine_surface_blocked_matches_uniform_four_hops(speed):
print(f"blocked speed={speed} four hops vs uniform PASS mass_into_fluid={mass:+.9g}", flush=True)
@pytest.mark.parametrize("missing", ["surface_velocity_gpu", "link_exchange_gpu"])
def test_missing_surface_device_buffer_rejects_before_launch(missing):
backend = implementation.FixedPatchStepper.__new__(implementation.FixedPatchStepper)
backend._boundary_hosts = [None, SimpleNamespace(objects=[object_spec(.002)])]
backend._level_jets = [(), ()]
backend._fine_curved = SimpleNamespace(fluid_idx_gpu=object(),
surface_velocity_gpu=object(), link_exchange_gpu=object())
setattr(backend._fine_curved, missing, None)
backend._curved = Mock()
with pytest.raises(RuntimeError, match="Surface boundary device buffers are missing"):
backend._apply_fine_boundary(object())
backend._curved.assert_not_called()
def test_close_frees_fine_boundary_after_stream_failure_once():
backend = implementation.FixedPatchStepper.__new__(implementation.FixedPatchStepper)
backend.ctx = SimpleNamespace(activate=nullcontext)
backend._lock = RLock()
backend._status = "ready"
stream, soa, first, second = Mock(), Mock(), Mock(), Mock()
stream.synchronize.side_effect = RuntimeError("stream synchronization failed")
backend._stream = stream
backend._fine_curved = soa
backend._fine_flags = first
backend._curved = Mock()
backend._levels = [object(), object()]
backend._allocations = [first, second]
with pytest.raises(RuntimeError, match="stream synchronization failed"):
backend.close()
assert backend._status == "closed"
assert backend._levels == [] and backend._allocations == []
assert backend._fine_curved is None and backend._fine_flags is None
assert backend._curved is None and backend._stream is None
backend.close()
stream.synchronize.assert_called_once_with()
soa.free.assert_called_once_with()
first.free.assert_called_once_with()
second.free.assert_called_once_with()
def test_cross_interface_donor_rejects_before_cuda():
cfg = LBMConfig(nx=24, ny=24, viscosity=.05)
source = SimObject(0, CircleGeometry(5.05, 12.25, .7),
@@ -92,9 +92,8 @@ def census(run, index, buffer):
def fine(block):
"""Two-level circle+jet construction on BOUNDS, all-blocked with a (12, 8) root tile;
block None is the dense route (dense reference test only, retired with the dense route)."""
return dict(refinements=[BOUNDS], level_blocks=[None if block is None else (12, 8), block],
"""Two-level circle+jet construction on BOUNDS, blocked with a (12, 8) root tile."""
return dict(refinements=[BOUNDS], level_blocks=[(12, 8), block],
level_objects=[(), (circle(),)], level_jets=[(), (jet(),)])
class RawExchangeGPU(unittest.TestCase):
@@ -147,26 +146,6 @@ class RawExchangeGPU(unittest.TestCase):
np.testing.assert_array_equal(measured.read_raw_exchange(0)['total']['impulse'],
obs['total']['impulse'])
def test_dense_fine_storage_reference(self):
"""Same fine resolution/initial data with dense storage: a numerical
reference for cut-link rows, not a uniform-whole-domain convergence study."""
with CudaContext(0) as ctx:
with FixedPatchStepper(CFG, context=ctx, **fine(None)) as dense, \
FixedPatchStepper(CFG, context=ctx, **fine((16, 20)), observe_raw_exchange=True) as blocked:
arrays = [state(l.shape) for l in blocked._levels]
dense.initialize(*arrays)
blocked.initialize(*arrays)
for speed in (0., .002, 0.):
dense.set_jet_amplitude(0, speed)
blocked.set_jet_amplitude(0, speed)
dense.step(); blocked.step()
for a, e in zip(blocked.download(), dense.download()):
np.testing.assert_array_equal(a, e)
count = dense._fine_curved.count
expected = download(dense._fine_curved.link_exchange_gpu, (count, 3))
actual = download(blocked._level_boundaries[1].last_raw, (count, 4))
np.testing.assert_array_equal(actual[:, :3], expected)
def test_depth_four_interval_storage_and_sampling(self):
cfg = LBMConfig(nx=48, ny=48, viscosity=.02)
obj = SimObject(0, geometry=CircleGeometry(16., 16., .5), center=(16., 16.), radius=.5)
@@ -3,36 +3,11 @@
Run CPU and GPU classes in separate processes. GPU callers must first check
resources/device availability; this file never selects or reserves a GPU.
Dense/blocked route parity (DenseBlockedRouteParity): the same configuration,
initial state, body/Jet schedule and tilings on three storage routes, only
level_blocks varying: D/D [None, None], D/B [None, FINE_BLOCK] (dense root,
the mixed pattern of current examples) and B/B [(12, 8), FINE_BLOCK]. Gate:
every public and owner array identical in FP32 bits and NaN masks at steps
0, 1, 2, 8 and after re-initialize; D/B and B/B raw-exchange ledgers equal.
No tolerance. A mismatch is reported with route pair/level/step and diagnosed
(boundary semantics vs arithmetic ordering), never widened.
Route capabilities (rejections checked before CUDA by test_route_capability_matrix):
combination D/D D/B B/B
periodic / slot_reservoir root, fine circle yes yes yes (parity gate)
actuated circle, velocity Jet + amplitude yes yes yes (parity gate)
raw/window exchange on a body level rejected yes yes
shifted storage rejected rejected yes
projected explosion rejected rejected yes
channel root rejected rejected yes
several L1 rectangles rejected rejected yes
depth >= 3 with a dense level rejected rejected yes
MRT, periodic root, circle/actuation/Jet yes yes yes (parity gate)
two fine bodies (Jet + passive, two actuated) yes yes yes (parity gate)
MRT with slot/reservoir root rejected on every route (on_site_velocity is SRT-only)
Slot/reservoir walls are compile-time WALL_* macros on a dense root and the
d_block_wall constant on a blocked root: BC-encoding parity, bit-identical by
this gate. Dense fine bodies use CurvedLinkSoA + CurvedBoundaryKernel
(n_objects > 0), blocked ones BlockCurvedBoundary; same links and physics,
compared through the population rows. MRT body-free is also covered by
test_refinement_mrt_gpu::test_dense_versus_blocked_mrt_two_level. Parked:
shifted and depth >= 3 (no dense route to compare).
Dense refinement storage was removed: every level is blocked and level_blocks
is required (TileContractCPU). The retired dense/blocked route parity gates
(BlockedBaseGPU exact parity, DenseBlockedRouteParity incl. MRT and two
bodies, the route capability matrix) passed and are recorded at 1573942 and
2469af7.
"""
from contextlib import ExitStack, contextmanager, nullcontext
from dataclasses import replace
@@ -51,7 +26,6 @@ from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY
from CelerisLab.lbm.refinement import ACTIVE, OUTER_GHOST, FixedPatch2D
from CelerisLab.lbm.refinement_stepper import FixedPatchStepper, blocked_base_tables
from test_refinement_blocked_fine import NX, NY, BOUNDS, config, state, circle, actuation
from test_refinement_jets import SPEED, jet
BASE_BLOCKS = ((24, 32), (12, 8), (3, 4), (1, 1))
FINE_BLOCK = (16, 20)
@@ -298,224 +272,33 @@ class BlockedBaseCPU(unittest.TestCase):
allocation.free.assert_called_once_with()
class BlockedBaseGPU(unittest.TestCase):
def test_bodyfree_exact_parity_and_raw_invalid_slots(self):
self.compare([(None, ((24, 32), (3, 4), (1, 1)), 8)])
class TileContractCPU(unittest.TestCase):
"""level_blocks is required: one positive (by, bx) tile per level, checked before CUDA."""
def test_actuated_circle_exact_parity_and_raw_invalid_slots(self):
self.compare([(.002, ((12, 8),), 8), (-.002, ((12, 8),), 8)])
def test_non_power_two_slot_width_exact_parity(self):
self.compare([(None, ((3, 4),), 7), (None, ((3, 4),), 10)], ("slot_reservoir",))
def compare(self, cases, modes=("periodic", "slot_reservoir")):
from CelerisLab.cuda.context import CudaContext
expected_source = Path(__file__).resolve().parents[3] / "src/CelerisLab/lbm/refinement_stepper.py"
self.assertEqual(Path(refinement_stepper.__file__).resolve(), expected_source)
with CudaContext(0) as context:
for mode in modes:
for speed, layouts, width in cases:
cfg = config(mode)
if mode == "slot_reservoir":
cfg = replace(cfg, wall_slot_width=width)
p = FixedPatch2D.from_config(cfg, BOUNDS)
arrays = state(NY, NX, 1), state(*p.fine_shape, 2)
reference = None
for base_block in (None, *layouts):
label = f"{mode} width={width} speed={speed} base={base_block}"
snapshots = []
with FixedPatchStepper(cfg, context=context, coarse_boundary=mode,
refinements=[BOUNDS], level_blocks=[base_block, FINE_BLOCK],
level_objects=[(), () if speed is None
else (circle(actuation(speed)),)]) as backend:
backend.initialize(*arrays)
for step in (0, 1, 2, 8):
if step == 1 and base_block is not None:
# Blocked staging dirties invalid input slots; its output
# must poison them. Dense parent does not stage there.
coarse = backend._levels[0]
spare = coarse.layout.to_blocked(np.full(coarse.shape, 42, np.float32))
cuda.memcpy_htod(coarse.spare, spare)
backend.step(step - backend.step_count)
public, owner = backend.download(), raw(backend)
self.assertEqual(backend.step_count, step)
self.assertTrue(np.isnan(owner[0][:, 6:18, 8:24]).all(), label)
if mode != "periodic":
self.assertTrue(np.isnan(owner[0][:, (0, -1), :]).all(), label)
snapshots.append(public + owner)
backend.initialize(*arrays)
self.assertEqual(backend.step_count, 0)
for actual, expected in zip(backend.download() + raw(backend), snapshots[0]):
exact(self, actual, expected, label + " reset")
if reference is None:
reference = snapshots
else:
for step, actual, expected in zip((0, 1, 2, 8), snapshots, reference):
for index, (a, b) in enumerate(zip(actual, expected)):
exact(self, a, b, f"{label} step={step} array={index}")
print(f"PASS {label}: full finite bits/NaN masks at 0,1,2,8 + reset", flush=True)
ROUTES = (("D/D", [None, None]), ("D/B", [None, FINE_BLOCK]), ("B/B", [(12, 8), FINE_BLOCK]))
# Two-body rows: patch (4, 4, 28, 20) gives a 40 x 56 fine window, tiled (8, 8).
WIDE_BOUNDS = (4, 4, 28, 20)
WIDE_ROUTES = (("D/D", [None, None]), ("D/B", [None, (8, 8)]), ("B/B", [(12, 8), (8, 8)]))
LEFT, RIGHT, RADIUS = (12.0, 12.0), (20.0, 12.0), 2.15
def body(ident, center, spec=None):
from CelerisLab.body.geometry.circle import CircleGeometry
from CelerisLab.body.objects import SimObject
return SimObject(ident, geometry=CircleGeometry(*center, RADIUS), center=center, radius=RADIUS,
surface_actuation=spec)
TWO_BODY_CASES = (
("jet + passive", lambda: [body(0, LEFT), body(1, RIGHT)],
lambda: [actuation(SPEED).as_jet(LEFT, RADIUS)], -SPEED),
("actuated +/-", lambda: [body(0, LEFT, actuation(SPEED)), body(1, RIGHT, actuation(-SPEED))],
lambda: [], None),
)
def same_summary(a, b):
"""Exact equality of nested exchange summaries (dicts, sequences, arrays, numbers)."""
if isinstance(a, dict):
return a.keys() == b.keys() and all(same_summary(a[k], b[k]) for k in a)
if isinstance(a, (list, tuple)):
return len(a) == len(b) and all(same_summary(x, y) for x, y in zip(a, b))
a, b = np.asarray(a), np.asarray(b)
if a.dtype.kind in "biufc" and b.dtype.kind in "biufc":
return a.shape == b.shape and np.array_equal(a, b, equal_nan=True)
return a.dtype == b.dtype and np.array_equal(a, b) # labels such as the window name
class DenseBlockedRouteCPU(unittest.TestCase):
def test_same_summary_compares_labels_and_every_number(self):
ledger = {"window": "last_base_step", "time_start": 3, "total": {"force": [1.0, -2.0],
"torque": 0.5}, "aperture": {"force": np.array([np.nan, 1.0])}}
self.assertTrue(same_summary(ledger, replace_value(ledger, ())))
for path, value in ((("window",), "last_substep"), (("time_start",), 4),
(("total", "force"), [1.0, -2.000001]), (("total", "torque"), 0.25),
(("aperture", "force"), np.array([0.0, 1.0]))):
with self.subTest(path=path):
self.assertFalse(same_summary(ledger, replace_value(ledger, path, value)))
self.assertFalse(same_summary(replace_value(ledger, path, value), ledger))
def test_route_capability_matrix(self):
def test_level_blocks_required_and_validated(self):
context = Mock(activate=Mock(side_effect=AssertionError("CUDA reached")))
periodic = config("periodic")
inner = (BOUNDS[0] + 2, BOUNDS[1] + 2, BOUNDS[2] - 2, BOUNDS[3] - 2)
second = (BOUNDS[0], BOUNDS[1], BOUNDS[0] + 2, BOUNDS[1] + 2)
cases = (
("raw exchange on dense body level", periodic, dict(
level_blocks=[None, None], level_objects=[(), (circle(),)], observe_raw_exchange=True)),
("shifted D/D", replace(periodic, ddf_shifting=True), dict(level_blocks=[None, None])),
("shifted D/B", replace(periodic, ddf_shifting=True), dict(level_blocks=[None, FINE_BLOCK])),
("projected D/D", periodic, dict(level_blocks=[None, None], explosion="projected_linear")),
("projected D/B", periodic, dict(level_blocks=[None, FINE_BLOCK], explosion="projected_linear")),
("channel D/B", periodic, dict(level_blocks=[None, FINE_BLOCK], coarse_boundary="channel")),
("several L1 rectangles", periodic, dict(refinements=[[BOUNDS, second]], level_blocks=[None, None])),
("depth 3 dense", periodic, dict(refinements=[BOUNDS, inner], level_blocks=[None, None, None])),
)
for label, cfg, kw in cases:
kw.setdefault("refinements", [BOUNDS])
cfg = config("periodic")
for label, kw in (("omitted", {}),
("None list", dict(level_blocks=None)),
("too few levels", dict(level_blocks=[(12, 8)])),
("too many levels", dict(level_blocks=[(12, 8), FINE_BLOCK, FINE_BLOCK])),
("None root", dict(level_blocks=[None, FINE_BLOCK])),
("None fine", dict(level_blocks=[(12, 8), None])),
("wrong arity", dict(level_blocks=[(12,), FINE_BLOCK])),
("non-integer", dict(level_blocks=[(12, 8.0), FINE_BLOCK])),
("bool", dict(level_blocks=[(True, 8), FINE_BLOCK])),
("non-positive", dict(level_blocks=[(0, 8), FINE_BLOCK])),
("non-tuple", dict(level_blocks=[12, FINE_BLOCK])),
("not dividing", dict(level_blocks=[(7, 8), FINE_BLOCK]))):
with self.subTest(label), self.assertRaises(ValueError):
FixedPatchStepper(cfg, context=context, **kw)
FixedPatchStepper(cfg, context=context, refinements=[BOUNDS], **kw)
with self.assertRaisesRegex(ValueError, "whole-level tile"):
FixedPatchStepper(cfg, context=context, refinements=[BOUNDS])
context.activate.assert_not_called()
def replace_value(tree, path, value=None):
"""Deep copy of a nested dict with the item at ``path`` replaced (no change for ())."""
if isinstance(tree, dict):
return {k: (value if path == (k,) else replace_value(v, path[1:], value) if path[:1] == (k,)
else replace_value(v, (), None)) for k, v in tree.items()}
return np.array(tree, copy=True) if isinstance(tree, np.ndarray) else tree
class DenseBlockedRouteParity(unittest.TestCase):
"""D/D vs D/B vs B/B, unshifted SRT, one circle or one velocity Jet; see module docstring."""
def run_routes(self, mode, objects=(), jets=(), switch=None, cfg=None, bounds=BOUNDS, routes=ROUTES):
from CelerisLab.cuda.context import CudaContext
cfg = cfg or config(mode)
p = FixedPatch2D.from_config(cfg, bounds)
arrays = state(NY, NX, 1), state(*p.fine_shape, 2)
bodies = sorted({obj.obj_id for obj in objects})
results = {}
with CudaContext(0) as context:
for name, blocks in routes:
observe = bool(objects) and blocks[1] is not None
snapshots, ledgers = [], []
with FixedPatchStepper(cfg, context=context, coarse_boundary=mode, refinements=[bounds],
level_blocks=blocks, level_objects=[(), tuple(objects)],
level_jets=[(), tuple(jets)],
observe_raw_exchange=observe) as backend:
backend.initialize(*arrays)
for step in (0, 1, 2, 8):
if step == 8 and switch is not None:
backend.set_jet_amplitude(0, switch)
backend.step(step - backend.step_count)
snapshots.append(backend.download() + raw(backend))
if observe and step:
ledgers.append(tuple(backend.read_raw_exchange(body, window=window)
for body in bodies
for window in ("last_base_step", "last_substep")))
backend.initialize(*arrays)
snapshots.append(backend.download() + raw(backend))
results[name] = snapshots, ledgers
return results
def check(self, label, results):
reference, _ = results["D/D"]
for name in ("D/B", "B/B"):
for step, actual, expected in zip((0, 1, 2, 8, "reset"), results[name][0], reference):
for index, (a, b) in enumerate(zip(actual, expected)):
exact(self, a, b, f"{label} {name} vs D/D step={step} array={index}")
if results["D/B"][1]:
for step, a, b in zip((1, 2, 8), results["D/B"][1], results["B/B"][1]):
self.assertTrue(same_summary(a, b), f"{label} raw exchange D/B vs B/B step={step}")
for name, (snapshots, _) in results.items():
masses = [float(np.nansum(level.astype(np.float64))) for level in snapshots[3][:2]]
print(f"{label} {name}: step 8 mass by level {masses}", flush=True)
print(f"PASS {label}: D/D = D/B = B/B bits/NaN masks at 0,1,2,8 + reset", flush=True)
def test_body_free(self):
for mode in ("periodic", "slot_reservoir"):
self.check(f"{mode} body-free", self.run_routes(mode))
def test_passive_circle(self):
for mode in ("periodic", "slot_reservoir"):
self.check(f"{mode} circle", self.run_routes(mode, objects=[circle()]))
def test_actuated_circle(self):
for mode in ("periodic", "slot_reservoir"):
for speed in (SPEED, -SPEED):
self.check(f"{mode} actuated {speed}",
self.run_routes(mode, objects=[circle(actuation(speed))]))
def test_velocity_jet_with_amplitude_switch(self):
for mode in ("periodic", "slot_reservoir"):
self.check(f"{mode} jet switch",
self.run_routes(mode, objects=[circle()], jets=[jet()], switch=-SPEED))
def test_mrt_bodies(self):
"""MRT, periodic root (on_site_velocity slot/reservoir is SRT-only on every route)."""
cfg = replace(config("periodic"), collision="MRT")
for label, objects, jets, switch in (("circle", [circle()], (), None),
("actuated +", [circle(actuation(SPEED))], (), None),
("actuated -", [circle(actuation(-SPEED))], (), None),
("jet switch", [circle()], [jet()], -SPEED)):
self.check(f"MRT periodic {label}",
self.run_routes("periodic", objects, jets, switch, cfg=cfg))
def test_two_bodies(self):
"""Two circles on a wider patch (same tiles on every route): Jet + passive, two actuated."""
for mode in ("periodic", "slot_reservoir"):
for label, objects, jets, switch in TWO_BODY_CASES:
self.check(f"{mode} two bodies {label}",
self.run_routes(mode, objects(), jets(), switch, bounds=WIDE_BOUNDS,
routes=WIDE_ROUTES))
# The documented whole-level tiles are always valid (construction reaches CUDA).
whole = [(NY, NX), FixedPatch2D.from_config(cfg, BOUNDS).fine_shape]
with self.assertRaisesRegex(AssertionError, "CUDA reached"):
FixedPatchStepper(cfg, context=context, refinements=[BOUNDS], level_blocks=whole)
if __name__ == "__main__":
@@ -1,9 +1,8 @@
"""Opt-in blocked fine storage for FixedPatchStepper: CPU tables plus GPU gate.
"""Blocked fine storage for FixedPatchStepper: CPU tables, staging and launch contracts.
Frozen GPU gate (same stream, 8 coarse steps, nonequilibrium start), in both
periodic and slot_reservoir modes: dense fine storage vs 1x1, 2x2 and 4x2
fine blocks give bit-identical masked downloads AND bit-identical raw coarse
and unpacked owner-core fine arrays, including patch-ghost values/NaN masks.
The former GPU gate (dense fine storage vs 1x1, 2x2 and 4x2 fine blocks,
bit-identical) was retired with the dense route; its result and the
dense/blocked route parity are recorded at 1573942 and 2469af7.
"""
from contextlib import nullcontext
import unittest
@@ -125,7 +124,7 @@ class BlockedFineCPU(unittest.TestCase):
patch.object(compiler, "load_module", return_value=module), \
patch.object(cuda, "mem_alloc", side_effect=alloc), patch.object(cuda, "memcpy_htod"):
backend = FixedPatchStepper(config("periodic"), context=context, refinements=[BOUNDS],
level_blocks=[None, FINE_BLOCKS[0]])
level_blocks=[(12, 8), FINE_BLOCKS[0]])
self.assertEqual(len(backend.fine_layout.halo_dst), 0)
self.assertIsNone(backend._levels[1].halo_dst)
backend._gather_halo(backend._levels[1], Mock()) # no launch without pairs
@@ -177,7 +176,7 @@ class BlockedFineCPU(unittest.TestCase):
patch.object(cuda, "memcpy_htod"), \
patch.object(curved_links, "CurvedLinkSoA", side_effect=AssertionError("dense SoA built")):
backend = FixedPatchStepper(config("periodic"), context=context, refinements=[BOUNDS],
level_blocks=[None, FINE_BLOCKS[1]],
level_blocks=[(12, 8), FINE_BLOCKS[1]],
level_objects=[(), tuple(objects)])
backend.initialize(state(NY, NX, 1), state(*backend.patch.fine_shape, 2))
return backend, calls
@@ -187,8 +186,8 @@ class BlockedFineCPU(unittest.TestCase):
backend, calls = self.mocked_backend([circle(spec)])
c, f = backend._levels
coarse, cur, spare = c.current, f.current, f.spare
flags, kind = backend._fine_flags, backend._levels[1].kind
velocity = backend._link_tables.get("surface_velocity")
flags, kind = backend._level_boundaries[1].flags, backend._levels[1].kind
velocity = backend._level_boundaries[1].device.get("surface_velocity")
self.assertEqual(velocity is None, spec is None)
backend.step(1)
backend.close()
@@ -196,8 +195,9 @@ class BlockedFineCPU(unittest.TestCase):
self.assertEqual(order, [("BlockExplode", coarse),
("BlockHaloCopy", cur), ("BlockCurvedBoundary", cur), ("BlockStep", cur),
("BlockHaloCopy", spare), ("BlockCurvedBoundary", spare),
("BlockStep", spare), ("PatchCoarseStep", coarse)])
steps = [a for name, a in calls if name == "BlockStep"]
("BlockStep", spare), ("BlockHaloCopy", coarse),
("BlockCoalesce", cur), ("BlockStep", coarse)])
steps = [a for name, a in calls if name == "BlockStep"][:2] # the two fine hops
self.assertEqual([a[1] for a in steps], [spare, cur])
for a in steps:
self.assertIs(a[3], flags) # blocked flags, not null
@@ -213,7 +213,7 @@ class BlockedFineCPU(unittest.TestCase):
for speed in (0.0, .002): # any spec, even zero amplitude
backend, calls = self.mocked_backend([circle(actuation(speed))])
self.assertIsNotNone(backend._fine_tables["surface_velocity"])
del backend._link_tables["surface_velocity"]
del backend._level_boundaries[1].device["surface_velocity"]
with self.assertRaisesRegex(RuntimeError, "Surface velocity device buffer is missing"):
backend.step(1)
names = [name for name, _ in calls]
@@ -279,93 +279,5 @@ class BlockedFineCPU(unittest.TestCase):
self.assertEqual(set(np.unique(kind).tolist()) - {ACTIVE, INNER_GHOST, OUTER_GHOST}, set())
class BlockedFineGPU(unittest.TestCase):
def test_dense_and_blocked_fine_storage_are_bit_identical(self):
from CelerisLab.cuda.context import CudaContext
for mode in ("periodic", "slot_reservoir"):
cfg = config(mode)
patch = FixedPatch2D.from_config(cfg, BOUNDS)
g_coarse, g_fine = state(NY, NX, 1), state(*patch.fine_shape, 2)
results = []
with CudaContext(0) as context:
for fine_block in (None, *FINE_BLOCKS):
with FixedPatchStepper(cfg, context=context, coarse_boundary=mode,
refinements=[BOUNDS],
level_blocks=[None, fine_block]) as backend:
backend.initialize(g_coarse, g_fine)
backend.step(STEPS)
self.assertEqual(backend.step_count, STEPS)
results.append((backend.download(), raw(backend)))
(dense_public, dense_raw), blocked = results[0], results[1:]
kind = patch.fine_kind
self.assertTrue(np.isnan(dense_raw[1][:, kind == OUTER_GHOST]).all(), mode)
# INNER populations may be NaN after hop 2 (pulled from poisoned OUTER);
# their exact pattern is pinned by the dense-vs-blocked comparison below.
self.assertTrue(np.isfinite(dense_raw[1][:, kind == ACTIVE]).all(), mode)
for fine_block, (public, raw_arrays) in zip(FINE_BLOCKS, blocked):
msg = f"{mode} fine_block={fine_block}"
for a, b in zip(public + tuple(raw_arrays), dense_public + tuple(dense_raw)):
np.testing.assert_array_equal(a, b, err_msg=msg) # values and NaN masks
print(f"{mode}: dense vs {len(blocked)} fine-block layouts bit-identical "
f"after {STEPS} coarse steps", flush=True)
def test_passive_circle_dense_and_blocked_fine_storage_are_bit_identical(self):
from CelerisLab.cuda.context import CudaContext
for mode in ("periodic", "slot_reservoir"):
cfg = config(mode)
patch, boundary = fine_circle(mode)
solid = (boundary.flags.reshape(patch.fine_shape) & SOLID) != 0
g_coarse, g_fine = state(NY, NX, 1), state(*patch.fine_shape, 2)
results = []
with CudaContext(0) as context:
for fine_block in (None, *FINE_BLOCKS):
with FixedPatchStepper(cfg, context=context, coarse_boundary=mode,
refinements=[BOUNDS], level_blocks=[None, fine_block],
level_objects=[(), (circle(),)]) as backend:
backend.initialize(g_coarse, g_fine)
backend.step(STEPS)
results.append((backend.download(), raw(backend)))
(dense_public, dense_raw), blocked = results[0], results[1:]
fine = dense_raw[1]
self.assertTrue(np.isfinite(fine[:, (patch.fine_kind == ACTIVE) & ~solid]).all(), mode)
self.assertTrue(np.isnan(fine[:, solid]).all(), mode)
self.assertTrue(np.isnan(fine[:, patch.fine_kind == OUTER_GHOST]).all(), mode)
for fine_block, (public, raw_arrays) in zip(FINE_BLOCKS, blocked):
msg = f"circle {mode} fine_block={fine_block}"
for a, b in zip(public + tuple(raw_arrays), dense_public + tuple(dense_raw)):
np.testing.assert_array_equal(a, b, err_msg=msg) # values and NaN masks
print(f"circle {mode}: dense vs {len(blocked)} fine-block layouts bit-identical "
f"after {STEPS} coarse steps", flush=True)
def test_actuated_circle_dense_and_blocked_fine_storage_are_bit_identical(self):
from CelerisLab.cuda.context import CudaContext
for mode in ("periodic", "slot_reservoir"):
for speed in (.002, -.002):
cfg = config(mode)
patch, boundary = fine_circle(mode, actuation(speed))
solid = (boundary.flags.reshape(patch.fine_shape) & SOLID) != 0
g_coarse, g_fine = state(NY, NX, 1), state(*patch.fine_shape, 2)
results = []
with CudaContext(0) as context:
for fine_block in (None, *FINE_BLOCKS):
with FixedPatchStepper(cfg, context=context, coarse_boundary=mode,
refinements=[BOUNDS], level_blocks=[None, fine_block],
level_objects=[(), (circle(actuation(speed)),)]) as backend:
backend.initialize(g_coarse, g_fine)
backend.step(STEPS)
results.append((backend.download(), raw(backend)))
(dense_public, dense_raw), blocked = results[0], results[1:]
fine = dense_raw[1]
self.assertTrue(np.isfinite(fine[:, (patch.fine_kind == ACTIVE) & ~solid]).all())
self.assertTrue(np.isnan(fine[:, solid]).all())
self.assertTrue(np.isnan(fine[:, patch.fine_kind == OUTER_GHOST]).all())
for fine_block, (public, raw_arrays) in zip(FINE_BLOCKS, blocked):
msg = f"actuated {speed} {mode} fine_block={fine_block}"
for a, b in zip(public + tuple(raw_arrays), dense_public + tuple(dense_raw)):
np.testing.assert_array_equal(a, b, err_msg=msg) # values and NaN masks
print(f"actuated {speed:+.3f} {mode}: dense vs {len(blocked)} fine-block layouts "
f"bit-identical after {STEPS} coarse steps", flush=True)
if __name__ == "__main__":
unittest.main()
@@ -53,24 +53,12 @@ def totals(arrays, maps):
def level_storage(level):
"""Raw device storage of one level as dense (9, ny, nx), dense or blocked."""
if level.layout is None:
raw = np.empty(level.shape, np.float32)
cuda.memcpy_dtoh(raw, level.current)
return raw
"""Raw device storage of one blocked level as dense (9, ny, nx)."""
blocked = np.empty(level.layout.blocked_shape, np.float32)
cuda.memcpy_dtoh(blocked, level.current)
return level.layout.to_dense(blocked)
def global_array(module, name, shape):
out = np.empty(shape, dtype=np.uint32)
ptr, size = module.get_global(name)
assert size == out.nbytes
cuda.memcpy_dtoh(out, ptr)
return out
class CoupledGPU(unittest.TestCase):
def test_actual_tg_rates_one_interval(self):
"""Algebraic rate/phase check, not a new physical Taylor--Green benchmark.
@@ -181,11 +169,8 @@ class CoupledGPU(unittest.TestCase):
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, audit):
"""Rest and seeded nonequilibrium, 8 steps against the scalar reference.
audit=True also checks the dense CELERIS_REFINEMENT_AUDIT ownership counters.
"""
def eight_steps(self, backend, spec, maps):
"""Rest and seeded nonequilibrium, 8 steps against the scalar reference."""
p = backend.patch
for case in ("rest", "nonequilibrium_seed1"):
with self.subTest(case=case):
@@ -196,11 +181,6 @@ class CoupledGPU(unittest.TestCase):
for a in (state.coarse_f, state.fine_f)))
initial = postcollision(state, spec, maps)
backend.initialize(*initial)
if audit: # dense audit build only
for level in backend._levels:
for name in ("collision_counts", "explosion_counts", "incoming_counts"):
ptr, size = level.module.get_global(name)
cuda.memset_d8(ptr, 0, size)
before = totals(backend.download(), maps)
max_error = max_drift = max_moment_error = max_corner = 0.0
for step in range(1, 9):
@@ -220,17 +200,7 @@ class CoupledGPU(unittest.TestCase):
self.assertLessEqual(drift, MOMENT_ATOL)
self.assertLessEqual(moment_error, MOMENT_ATOL)
max_drift, max_moment_error = max(max_drift, drift), max(max_moment_error, moment_error)
c, f = backend._levels
if audit:
np.testing.assert_array_equal(global_array(c.module, "collision_counts", (16, 16)),
step * maps.coarse_active)
np.testing.assert_array_equal(global_array(f.module, "collision_counts", (24, 24)),
2 * step * maps.fine_active)
np.testing.assert_array_equal(global_array(f.module, "explosion_counts", (24, 24)),
step * ~maps.fine_active)
owners = np.zeros((9, 16, 16), np.uint32)
owners.ravel()[p.coalesce_coarse] = step
np.testing.assert_array_equal(global_array(c.module, "incoming_counts", owners.shape), owners)
f = backend._levels[1]
raw = level_storage(f)
np.testing.assert_array_equal(np.isfinite(raw), diagnostic["hop2_valid"])
self.assertTrue(np.isfinite(raw.ravel()[p.coalesce_fine]).all())
@@ -245,7 +215,7 @@ class CoupledGPU(unittest.TestCase):
max_corner = max(max_corner, error)
print(f"{case}: steps=8 max_population={max_error:.9g} "
f"max_leaf_drift={max_drift:.9g} max_leaf_error={max_moment_error:.9g} "
f"max_corner_packet={max_corner:.9g}; ownership/counts PASS", flush=True)
f"max_corner_packet={max_corner:.9g}; PASS", flush=True)
def test_bounded_equivalence_blocked(self):
"""The same eight-step reference and launch-failure checks on whole-level tiles."""
@@ -255,7 +225,7 @@ class CoupledGPU(unittest.TestCase):
with CudaContext(0) as context:
with PeriodicPatchStepper(cfg, context=context, refinements=[(4, 4, 12, 12)],
level_blocks=[(16, 16), (24, 24)]) as backend:
self.eight_steps(backend, spec, maps, audit=False)
self.eight_steps(backend, spec, maps)
with patch.object(backend._pairs[0], "explode", side_effect=RuntimeError("injected launch failure")):
with self.assertRaisesRegex(RuntimeError, "injected"):
backend.step()
@@ -264,59 +234,5 @@ class CoupledGPU(unittest.TestCase):
operation()
self.assertFalse(context.is_closed)
def test_bounded_equivalence_and_ownership(self):
self.assertEqual(Path(CelerisLab.__file__).resolve(), ROOT / "src/CelerisLab/__init__.py")
self.assertEqual(Path(compiler._KERNEL_DIR).resolve(), ROOT / "src/CelerisLab/lbm/kernels")
spec = ref.make_spec(collision="srt")
maps = ref.build_maps(spec)
cfg = LBMConfig(nx=16, ny=16, viscosity=.1)
compile_original = compiler.compile_kernel
def audit_compile(*args, **kwargs):
# Only this private build's entry differs. Production source unchanged.
entry = Path(kwargs["build_dir"]) / "refinement_audit.cu"
entry.write_text('#define CELERIS_REFINEMENT_AUDIT\n#include "refinement.cu"\n')
kwargs["source"] = entry.name
return compile_original(*args, **kwargs)
with CudaContext(0) as context:
with patch.object(compiler, "compile_kernel", audit_compile):
backend = PeriodicPatchStepper(cfg, context=context, refinements=[(4, 4, 12, 12)])
with backend:
# Production allocations: four DDF arrays and four compact tables.
device_bytes = 2 * 9 * (16**2 + 24**2) * 4
p = backend.patch
device_bytes += sum(a.nbytes for a in (p.ghost_fine, p.ghost_coarse,
p.coalesce_coarse, p.coalesce_fine))
audit_bytes = sum((2 + 9) * np.prod(level.shape[1:]) * 4
for level in backend._levels)
self.assertLess(device_bytes + audit_bytes, 128 * 1024)
print(f"device_arrays={device_bytes} audit_arrays={audit_bytes} bytes", flush=True)
self.eight_steps(backend, spec, maps, audit=True)
# Fail before any additional time step; no stale synchronized export.
with patch.object(backend._pairs[0], "explode", side_effect=RuntimeError("injected launch failure")):
with self.assertRaisesRegex(RuntimeError, "injected"):
backend.step()
with self.assertRaisesRegex(RuntimeError, "failed"):
backend.download()
with self.assertRaisesRegex(RuntimeError, "failed"):
backend.step()
self.assertFalse(context.is_closed)
with self.assertRaisesRegex(RuntimeError, "closed"):
backend.download()
# Compile and execute the actual non-audited entry as well: one step.
with PeriodicPatchStepper(cfg, context=context, refinements=[(4, 4, 12, 12)]) as plain:
rest = postcollision(ref.make_state(spec), spec, maps)
plain.initialize(*rest)
plain.step()
for got, expected in zip(plain.download(), rest):
np.testing.assert_allclose(got, expected, rtol=0, atol=POP_ATOL, equal_nan=True)
print("normal-build rest: steps=1 PASS; source/kernel paths pinned", flush=True)
rss = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / 1024
print(f"host_peak_rss={rss:.3f} MiB", flush=True)
if __name__ == "__main__":
unittest.main(verbosity=2)
@@ -4,12 +4,12 @@ Jets use the same coarse cell-edge frame as fine_objects. The equivalent jet
of Astra's actuated fixture is SurfaceActuation(0.7, 0.85, U).as_jet((16, 12),
2.15). Frozen GPU gates (8 coarse steps, both coarse boundary modes):
G1 jet run: 1x1/2x2/4x2 fine blocks (blocked (12, 8) root) bit-identical to the
FINE_BLOCKS[0] reference (values + NaN masks); dense leg kept until dense retirement;
FINE_BLOCKS[0] reference (values + NaN masks);
G2 jet vs legacy actuation (blocked reference): max |dG| over finite fine ACTIVE fluid
populations <= 1e-6 (velocity rows differ only by FP32 rounding; reported);
G3 amplitude change after 4 steps (+U -> -U), 4 more steps: blocked layouts
bit-identical to the FINE_BLOCKS[0] reference (dense leg until retirement) and the
run differs from the unchanged run.
bit-identical to the FINE_BLOCKS[0] reference and the run differs from the
unchanged run.
"""
from contextlib import nullcontext
import unittest
@@ -29,8 +29,8 @@ SPEED = .002
def root_and(fine_block):
"""All-blocked levels with a (12, 8) root tile; fine_block None is the dense route."""
return [None if fine_block is None else (12, 8), fine_block]
"""Blocked levels with a (12, 8) root tile."""
return [(12, 8), fine_block]
def jet(speed=SPEED):
@@ -107,7 +107,7 @@ class JetStepperCPU(unittest.TestCase):
def test_blocked_amplitude_update_rewrites_only_velocity_column(self):
backend, htod = mocked([circle()], [jet()], FINE_BLOCKS[1])
before = backend._fine_tables["surface_velocity"].copy()
device = backend._link_tables["surface_velocity"]
device = backend._level_boundaries[1].device["surface_velocity"]
with patch.object(cuda, "memcpy_htod", htod):
htod.reset_mock()
backend.set_jet_amplitude(0, -SPEED)
@@ -123,7 +123,7 @@ class JetStepperCPU(unittest.TestCase):
backend.close()
def test_copy_failure_marks_failed_and_keeps_host_amplitude(self):
for fine_block in (None, FINE_BLOCKS[1]): # dense SoA and blocked tables
for fine_block in (FINE_BLOCKS[1],): # blocked tables
backend, _ = mocked([circle()], [jet()], fine_block)
status = backend._status
for bad, error in ((1.5, ValueError), (float("nan"), ValueError)):
@@ -133,9 +133,7 @@ class JetStepperCPU(unittest.TestCase):
with self.assertRaises(IndexError):
backend.set_jet_amplitude(1, -SPEED)
self.assertEqual(backend._status, status)
dense = backend._fine_curved
before = (dense.surface_velocity if dense is not None
else backend._fine_tables["surface_velocity"]).copy()
before = backend._fine_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)):
@@ -143,11 +141,7 @@ class JetStepperCPU(unittest.TestCase):
backend.set_jet_amplitude(0, -SPEED)
self.assertEqual(backend._status, "failed", fine_block)
self.assertEqual(backend._jets[0].amplitude, SPEED)
if dense is not None:
self.assertEqual(dense.jets[0].amplitude, SPEED)
np.testing.assert_array_equal(dense.surface_velocity, before)
else:
np.testing.assert_array_equal(backend._fine_tables["surface_velocity"], before)
np.testing.assert_array_equal(backend._fine_tables["surface_velocity"], before)
with self.assertRaisesRegex(RuntimeError, "failed"):
backend.step(1)
with self.assertRaisesRegex(RuntimeError, "failed"):
@@ -158,6 +152,7 @@ class JetStepperCPU(unittest.TestCase):
with patch.object(cuda, "Stream", side_effect=AssertionError("CUDA touched")):
with self.assertRaisesRegex(ValueError, "SurfaceActuation"):
FixedPatchStepper(config("periodic"), context=Mock(), refinements=[BOUNDS],
level_blocks=root_and(FINE_BLOCKS[1]),
level_objects=[(), (circle(actuation(SPEED)),)],
level_jets=[(), (jet(),)])
@@ -192,9 +187,6 @@ class JetStepperGPU(unittest.TestCase):
public, raw_arrays = run(mode, fine_block, objects=[circle()], jets=[jet()])
for a, b in zip(public + tuple(raw_arrays), reference[0] + tuple(reference[1])):
np.testing.assert_array_equal(a, b, err_msg=f"G1 {mode} {fine_block}")
dense = run(mode, None, objects=[circle()], jets=[jet()]) # dense leg (D2 retires)
for a, b in zip(dense[0] + tuple(dense[1]), reference[0] + tuple(reference[1])):
np.testing.assert_array_equal(a, b, err_msg=f"G1 {mode} dense")
legacy = run(mode, FINE_BLOCKS[0], objects=[circle(actuation(SPEED))], jets=())
diff = np.abs(reference[1][1][:, fluid] - legacy[1][1][:, fluid]).max()
print(f"G2 {mode}: jet vs legacy fine max |dG| = {diff:.3e}", flush=True)
@@ -205,12 +197,12 @@ class JetStepperGPU(unittest.TestCase):
reference = run(mode, FINE_BLOCKS[0], objects=[circle()], jets=[jet()], switch=-SPEED)
steady = run(mode, FINE_BLOCKS[0], objects=[circle()], jets=[jet()])
self.assertFalse(np.array_equal(reference[1][1], steady[1][1], equal_nan=True), mode)
for fine_block in (*FINE_BLOCKS[1:], None): # None: dense leg (D2 retires)
for fine_block in FINE_BLOCKS[1:]:
public, raw_arrays = run(mode, fine_block, objects=[circle()], jets=[jet()],
switch=-SPEED)
for a, b in zip(public + tuple(raw_arrays), reference[0] + tuple(reference[1])):
np.testing.assert_array_equal(a, b, err_msg=f"G3 {mode} {fine_block}")
print(f"G3 {mode}: amplitude switch, {len(FINE_BLOCKS) - 1} layouts + dense "
print(f"G3 {mode}: amplitude switch, {len(FINE_BLOCKS) - 1} layouts "
f"identical to {FINE_BLOCKS[0]}", flush=True)
@@ -54,25 +54,6 @@ class RefinementMRTGPU(unittest.TestCase):
else:
for a, b in zip(reference, raw.values): np.testing.assert_array_equal(a, b)
def test_dense_versus_blocked_mrt_two_level(self):
"""Unshifted MRT: all-dense PatchCoarseStep/PatchFineStep vs mixed/all-blocked storage."""
from test_refinement_blocked_base import BASE_BLOCKS, FINE_BLOCK
from test_refinement_blocked_fine import NX, NY, BOUNDS
import test_refinement_blocked_fine as fixture
cfg = LBMConfig(nx=NX, ny=NY, viscosity=.1, collision="MRT")
with CudaContext(0) as context:
results = {}
layouts = ((None, None), *((base, FINE_BLOCK) for base in (None, *BASE_BLOCKS)))
for base_block, fine_block in layouts:
with FixedPatchStepper(cfg, context=context, refinements=[BOUNDS],
level_blocks=[base_block, fine_block]) as backend:
backend.initialize(*[fixture.state(*level.shape[1:], seed) for seed, level in enumerate(backend._levels)])
backend.step(4)
results[base_block, fine_block] = raw(backend)
for layout in layouts[1:]:
for level, (a, b) in enumerate(zip(results[layout], results[None, None])):
exact(self, a, b, f"layout={layout} level={level}")
def test_shifted_all_blocked_three_level_mrt_smoke(self):
"""Per-level shear rate omega_l = 1/(3*2**l*nu + 1/2); finite ACTIVE leaves."""
cfg = LBMConfig(nx=16, ny=16, viscosity=.1, collision="MRT", ddf_shifting=True)
@@ -75,7 +75,7 @@ class ThreeLevelCPU(unittest.TestCase):
with self.subTest(kw=kw), self.assertRaises(ValueError):
FixedPatchStepper(CFG, context=context, **kw)
# L2 bodies/Jets without an L2 level (one entry too many), and an L2 Jet without a body.
two_level = dict(refinements=[L1], level_blocks=[None, None])
two_level = dict(refinements=[L1], level_blocks=[B, B])
for kw in (dict(two_level, level_objects=[(), (), (child_body(),)]),
dict(two_level, level_jets=[(), (), (child_jet(),)]),
levels(jets=((), (), (child_jet(),)))):