ClaudeandClaude Opus 5.5
d318911221
The package already exposed it as CelerisLab.cuda.compiler; the module now
has that name and every tip caller, doc and current-tree test imports it
directly (no alias).
Historical-revision tooling is not edited: the Taylor-Green harness
(run_case.py, hash-pinned check_projected_pairs.py, assess receipts) only
runs against a135252, whose package has cuda/compiler_v2.py, and the J
campaign runners and dated docs describe past trees. Its two current-tree
tests now import that package from git (historical_package) instead of the
tip, and the kernel-root test first checks the accepted case.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
334 lines
16 KiB
Python
334 lines
16 KiB
Python
"""Two coarse/four fine steps: isolated fine circle, 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.
|
|
"""
|
|
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
|
|
import pycuda.driver as cuda
|
|
import pytest
|
|
|
|
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
|
|
|
|
|
|
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)
|
|
|
|
|
|
def read(pointer, shape, dtype=np.float32):
|
|
values = np.empty(shape, dtype)
|
|
cuda.memcpy_dtoh(values, pointer)
|
|
return values
|
|
|
|
|
|
def rest(shape):
|
|
return np.broadcast_to(W[:, None, None], (9, *shape)).astype(np.float32).copy()
|
|
|
|
|
|
def surface(speed):
|
|
return None if speed is None else dict(angle=.7, half_width=.85, normal_velocity=speed)
|
|
|
|
|
|
def object_spec(speed):
|
|
return SimObject(0, CircleGeometry(12.25, 12.25, 2.15), center=(12.25, 12.25),
|
|
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")
|
|
cases = importlib.util.module_from_spec(spec)
|
|
spec.loader.exec_module(cases)
|
|
with cases.make_sim(nx=48, ny=48, viscosity=.1, velocity=0, rho=1) as sim:
|
|
sim.add_body("circle", center=(24, 24), radius=4.3, surface_actuation=surface(speed))
|
|
sim.initialize()
|
|
# Uniform 48^2 has the same dx=dt=.5 representation. Its distant channel
|
|
# boundaries cannot influence the circle ROI within these four hops.
|
|
traces = []
|
|
for _ in range(4):
|
|
sim.run(1)
|
|
traces.append(sim.get_ddf().reshape(9, 48, 48)[:, 4:44, 4:44][:, roi].copy())
|
|
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)
|
|
|
|
|
|
@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),
|
|
center=(5.05, 12.25), radius=.7, surface_actuation=surface(.002))
|
|
context = Mock()
|
|
context.activate.side_effect = AssertionError("must reject before CUDA activation")
|
|
with patch.object(cuda, "mem_alloc") as allocate:
|
|
with pytest.raises(ValueError, match="donor crosses"):
|
|
implementation.FixedPatchStepper(cfg, context=context, refinements=[BOUNDS],
|
|
level_objects=[(), (source,)])
|
|
context.activate.assert_not_called()
|
|
allocate.assert_not_called()
|