"""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_v2 as 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()