--- a/src/CelerisLab/lbm/kernels/refinement.cu +++ b/src/CelerisLab/lbm/kernels/refinement.cu @@ -7,10 +7,17 @@ #include "core/descriptors.cuh" #include "core/layout.cuh" #include "core/params.cuh" +#include "core/flags.cuh" +#include "core/obs.cuh" #include "operators/macro.cuh" #include "operators/collision_srt.cuh" #include "streaming/pull_double_buffer.cuh" #include "boundary/planar_velocity.cuh" +#include "boundary/curved_boundary.cuh" + +extern "C" { +#include "step/curved_double.cu" +} #if DIM != 2 || NQ != 9 || COLLISION_MODEL != 0 || STORE_PRECISION != 0 || USE_DDF_SHIFTING || USE_LES || STREAMING_MODEL != 0 #error "Refinement requires unshifted FP32 D2Q9 SRT double buffers without LES" @@ -48,14 +55,15 @@ #endif } -extern "C" __global__ void PatchFineStep(const float* in, float* out) { +extern "C" __global__ void PatchFineStep(const float* in, float* out, const uint16_t* flags) { unsigned int x, y; unsigned long k; index_from_thread(x, y, k); if (x >= NX || y >= NY) return; float f[NQ]; // Outer halo is invalid after EACH hop, including when reusing the buffer // that held explosion data. Inner halo only streams; it never collides. - if (x < 2 || y < 2 || x >= NX - 2 || y >= NY - 2) { + if (x < 2 || y < 2 || x >= NX - 2 || y >= NY - 2 + || (flags != nullptr && is_solid(flags[k]))) { for (int q = 0; q < NQ; ++q) f[q] = NAN; } else { unsigned long j[NQ]; --- a/src/CelerisLab/lbm/kernels/step/aux_kernels.cu +++ b/src/CelerisLab/lbm/kernels/step/aux_kernels.cu @@ -24,61 +24,7 @@ #include "../streaming/esopull_semantic_helpers.cuh" #include "../boundary/esopull_curved_helpers.cuh" -__device__ __forceinline__ float action_omega(const float* action, int body_id) -{ - if (action == nullptr || body_id < 0) return 0.0f; - unsigned int base = (unsigned int)body_id * (unsigned int)OBS_BODY_SLOT_FLOATS; - return action[base + (unsigned int)(OBS_BODY_SLOT_FLOATS - 1)]; -} - -__global__ void CurvedBoundaryKernel( - fpxx* fi, - const unsigned int* cl_fluid_idx, - const unsigned char* cl_dir, - const float* cl_q, - const float* cl_rx, - const float* cl_ry, - const float* cl_rz, - const unsigned char* cl_fallback_class, - const int* cl_body_id, - const float* action, - float* obs, - unsigned int n_curved, - const float* surface_velocity = nullptr, float* link_exchange = nullptr) -{ - unsigned int tid = threadIdx.x + blockIdx.x * blockDim.x; - if (tid >= n_curved) return; - - const unsigned long k_f = (unsigned long)cl_fluid_idx[tid]; - const unsigned int dir = (unsigned int)cl_dir[tid]; - const float q = cl_q[tid]; - const int bid = cl_body_id[tid]; - const float rx = cl_rx[tid]; - const float ry = cl_ry[tid]; - const unsigned char fallback_class = cl_fallback_class[tid]; - const float omega = action_omega(action, bid); - -#if DIM == 2 - float Uw = -omega * ry; - float Vw = omega * rx; - if (surface_velocity != nullptr) { - // Independent prescribed fluid velocity on fixed geometry. - const float u = surface_velocity[2*tid], v = surface_velocity[2*tid+1]; - if (u != 0.0f || v != 0.0f) { Uw += u; Vw += v; } - } - (void)cl_rz; - apply_bouzidi_link(dir, q, k_f, fi, Uw, Vw, rx, ry, fallback_class, obs, bid, - link_exchange == nullptr ? nullptr : link_exchange + 3*tid); -#elif DIM == 3 - const float rz = cl_rz[tid]; - const float Uw = -omega * ry; - const float Vw = omega * rx; - const float Ww = 0.0f; - apply_bouzidi_link(dir, q, k_f, fi, Uw, Vw, Ww, - rx, ry, rz, fallback_class, obs, bid); -#endif -} -// End of CurvedBoundaryKernel +#include "curved_double.cu" // --------------------------------------------------------------------------- // PrepareEsoPullCurvedKernel -- read-only curved prep for EsoPull --- /dev/null +++ b/src/CelerisLab/lbm/kernels/step/curved_double.cu @@ -0,0 +1,61 @@ +// Shared compact double-buffer curved boundary; caller supplies core/BFL headers. +#ifndef CELERIS_CURVED_DOUBLE_CU +#define CELERIS_CURVED_DOUBLE_CU + +__device__ __forceinline__ float action_omega(const float* action, int body_id) +{ + if (action == nullptr || body_id < 0) return 0.0f; + unsigned int base = (unsigned int)body_id * (unsigned int)OBS_BODY_SLOT_FLOATS; + return action[base + (unsigned int)(OBS_BODY_SLOT_FLOATS - 1)]; +} + +__global__ void CurvedBoundaryKernel( + fpxx* fi, + const unsigned int* cl_fluid_idx, + const unsigned char* cl_dir, + const float* cl_q, + const float* cl_rx, + const float* cl_ry, + const float* cl_rz, + const unsigned char* cl_fallback_class, + const int* cl_body_id, + const float* action, + float* obs, + unsigned int n_curved, + const float* surface_velocity = nullptr, float* link_exchange = nullptr) +{ + unsigned int tid = threadIdx.x + blockIdx.x * blockDim.x; + if (tid >= n_curved) return; + + const unsigned long k_f = (unsigned long)cl_fluid_idx[tid]; + const unsigned int dir = (unsigned int)cl_dir[tid]; + const float q = cl_q[tid]; + const int bid = cl_body_id[tid]; + const float rx = cl_rx[tid]; + const float ry = cl_ry[tid]; + const unsigned char fallback_class = cl_fallback_class[tid]; + const float omega = action_omega(action, bid); + +#if DIM == 2 + float Uw = -omega * ry; + float Vw = omega * rx; + if (surface_velocity != nullptr) { + // Independent prescribed fluid velocity on fixed geometry. + const float u = surface_velocity[2*tid], v = surface_velocity[2*tid+1]; + if (u != 0.0f || v != 0.0f) { Uw += u; Vw += v; } + } + (void)cl_rz; + apply_bouzidi_link(dir, q, k_f, fi, Uw, Vw, rx, ry, fallback_class, obs, bid, + link_exchange == nullptr ? nullptr : link_exchange + 3*tid); +#elif DIM == 3 + const float rz = cl_rz[tid]; + const float Uw = -omega * ry; + const float Vw = omega * rx; + const float Ww = 0.0f; + apply_bouzidi_link(dir, q, k_f, fi, Uw, Vw, Ww, + rx, ry, rz, fallback_class, obs, bid); +#endif +} +// End of CurvedBoundaryKernel + +#endif --- a/src/CelerisLab/lbm/refinement.py +++ b/src/CelerisLab/lbm/refinement.py @@ -1,9 +1,9 @@ -"""CPU preparation for one fixed D2Q9 2:1 patch; no coupled stepper yet. +"""CPU preparation for one fixed D2Q9 2:1 patch and optional interior circles. Cell indices use y*nx+x; population slots use q*(nx*ny)+cell, matching core/layout.cuh. Bounds are half-open (x0, y0, x1, y1) in coarse cells. The four fine halo cells implement the pinned CPUREF volumetric geometry. -Topology does not encode physical boundary conditions or obstacle safety. +Base topology is independent of boundaries; prepare_fine_boundary checks them. """ from dataclasses import dataclass, replace from numbers import Integral @@ -11,7 +11,7 @@ import numpy as np from ..config import LBMConfig -from .descriptors import D2Q9_EX, D2Q9_EY +from .descriptors import D2Q9_EX, D2Q9_EY, FLUID, SOLID ACTIVE, INNER_GHOST, OUTER_GHOST = 0, 1, 2 HALO = 4 @@ -136,3 +136,157 @@ raise ValueError("step counts must be nonnegative integers") if fine_steps != 2 * coarse_steps: raise ValueError("levels are not synchronized: fine_steps must equal 2*coarse_steps") + + +@dataclass(frozen=True, eq=False) +class FineBoundary2D: + """Fixed local projection plus auditable, finite transfer dependency sets. + + transfer_reads rows are (hop, direction, cell), traced from coalescence at + hop 2 back to exploded G at hop 0. No previous coarse interval is traversed. + Cell footprints also include the whole explosion write/donor tables. + """ + objects: tuple + flags: np.ndarray + curved_host: tuple + boundary_cells: np.ndarray + transfer_fine_cells: np.ndarray + transfer_coarse_cells: np.ndarray + transfer_reads: np.ndarray + + +def _transfer_footprint(patch): + ny, nx = patch.fine_shape + cells = nx * ny + pending = [(2, *divmod(int(slot), cells)) for slot in patch.coalesce_fine.ravel()] + reads = set() + while pending: + hop, q, k = pending.pop() + if (hop, q, k) in reads: + continue + reads.add((hop, q, k)) + if hop == 0: + continue + y, x = divmod(k, nx) + kind = patch.fine_kind[y, x] + if kind == OUTER_GHOST: + raise ValueError("Transfer reads an invalid post-hop outer ghost") + # A collided output depends on all nine arrivals; ghosts only stream. + for direction in (range(9) if kind == ACTIVE else (q,)): + sx, sy = x - D2Q9_EX[direction], y - D2Q9_EY[direction] + if not (0 <= sx < nx and 0 <= sy < ny): + raise ValueError("Transfer dependency leaves the fine allocation") + pending.append((hop - 1, direction, sy * nx + sx)) + fine = set(map(int, patch.ghost_fine)) | {k for _, _, k in reads} + coarse = set(map(int, patch.ghost_coarse)) | set(map( + int, patch.coalesce_coarse % np.prod(patch.coarse_shape))) + return fine, coarse, _readonly(sorted(reads), np.int64) + + +def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False): + """Project fixed circles and reject boundary/transfer overlap before CUDA. + + Input centers/radii are in coarse physical coordinates (dx=dt=1). + SurfaceActuation retains its lattice-velocity meaning: input coarse units + are converted into independent fine-local objects with dx=dt=.5. + Original objects are not retained or mutated. No runtime geometry updates. + """ + if not objects: + return None + checked = replace(cfg) + checked.validate() + if (patch.coarse_shape != (cfg.ny, cfg.nx) or not checked.is_d2q9 + or checked.collision != "SRT" or checked.streaming != "double_buffer" + or checked.store_precision != "FP32" or checked.ddf_shifting or checked.les_enabled + or patch.tau_coarse != 3*checked.viscosity+.5): + raise ValueError("Fine boundary requires the matching fixed-patch configuration") + from ..body.geometry.circle import CircleGeometry + from ..body.objects import SimObject + from ..body.coupling.soa_packer import pack_cut_links_to_soa + + ny, nx = patch.fine_shape + fine_cfg = replace(cfg, nx=nx, ny=ny, viscosity=2 * cfg.viscosity) + flags = np.full(nx * ny, FLUID, np.uint16) + active = patch.fine_kind.ravel() == ACTIVE + local, links = [], [] + x0, y0, x1, y1 = patch.bounds + dx = dt = .5 + for bid, obj in enumerate(objects): + if (type(obj.geometry) is not CircleGeometry or obj.is_sensor or obj.is_force_region + or any(getattr(obj.state, name) != 0 for name in ("vx", "vy", "omega")) + or any(vars(obj.control).values())): + raise ValueError("Fine boundaries require stationary circles without other roles") + cx, cy = obj.geometry.center + radius = obj.geometry._radius + if (not np.isfinite([cx, cy, radius]).all() or radius <= 0 + or cx-radius < x0 or cx+radius > x1 or cy-radius < y0 or cy+radius > y1): + raise ValueError("Entire fine circle must lie inside the physical patch") + center = tuple((v-o)/dx-.5 for v, o in zip((cx, cy), patch.fine_origin)) + spec = obj.surface_actuation + if spec is not None: + spec = replace(spec, normal_velocity=spec.normal_velocity * dt / dx) + item = SimObject(bid, geometry=CircleGeometry(*center, radius/dx), + center=center, radius=radius/dx, surface_actuation=spec) + item.validate_surface_actuation(fine_cfg) + mask = item.get_flag_mask(nx, ny) + solid = (mask & SOLID) != 0 + if not solid.any() or np.any(solid & ~active) or np.any(flags[solid] & SOLID): + raise ValueError("Fine solids must be resolved, disjoint and entirely active") + flags[solid] = mask[solid] + local.append(item) + for item in local: + cut = item.geometry.build_cut_links(nx, ny, domain_flags=flags) + if not cut: + raise ValueError("Fine circle has no resolved cut links") + if item.surface_actuation is not None: + _, aperture = item.surface_actuation.project( + np.array([link.rx for link in cut], np.float32), + np.array([link.ry for link in cut], np.float32)) + if not aperture.any(): + raise ValueError("Surface aperture hits no resolved cut links") + for link in cut: + link.body_id = item.obj_id + links.extend(cut) + + boundary, writers, receivers = set(), set(), set() + for link in links: + r, q = link.fluid_idx, link.dir + y, x = divmod(r, nx) + sx, sy = x + D2Q9_EX[q], y + D2Q9_EY[q] + s = sy * nx + sx + if not (active[r] and flags[r] & FLUID and active[s] and flags[s] & SOLID): + raise ValueError("Cut-link receiver/source must be active fluid/solid") + writer = (s, q+1 if q % 2 else q-1) + if writer in writers: + raise ValueError("Duplicate fine boundary writer") + writers.add(writer) + receivers.add((r, q)) + boundary.update((r, s)) + if link.q < .5: + dx_, dy_ = x - D2Q9_EX[q], y - D2Q9_EY[q] + if not (0 <= dx_ < nx and 0 <= dy_ < ny and active[dy_*nx+dx_]): + raise ValueError("q<.5 donor crosses the fine active interface; fallback cannot hide it") + boundary.add(dy_*nx+dx_) + # Reverse coverage: every fluid pull across a solid neighbor needs a writer. + required = set() + for r in np.flatnonzero(active & ((flags & FLUID) != 0)): + for q in range(1, 9): + if flags[r + D2Q9_EX[q] + nx*D2Q9_EY[q]] & SOLID: + required.add((int(r), q)) + if required != receivers: + raise ValueError("Fine solid pull coverage does not match the cut links") + fine, coarse, reads = _transfer_footprint(patch) + if fine & (boundary | set(map(int, np.flatnonzero(flags & SOLID)))): + raise ValueError("Boundary stencil overlaps the two-hop transfer footprint") + for k in coarse: + y, x = divmod(k, cfg.nx) + if x0 <= x < x1 and y0 <= y < y1: + raise ValueError("Transfer coarse donor/target is not an uncovered bulk cell") + # Coarse reconstruction targets y=1,ny-2 pull from y=0..2,ny-3..ny-1. + if wall_aware and (y <= 2 or y >= cfg.ny-3): + raise ValueError("Transfer footprint touches the coarse wall reconstruction stencil") + return FineBoundary2D(tuple(local), _readonly(flags, np.uint16), + tuple(_readonly(a, a.dtype) for a in pack_cut_links_to_soa(links)), + _readonly(sorted(boundary), np.uint32), + _readonly(sorted(fine), np.uint32), + _readonly(sorted(coarse), np.uint32), reads) --- a/src/CelerisLab/lbm/refinement_stepper.py +++ b/src/CelerisLab/lbm/refinement_stepper.py @@ -1,5 +1,5 @@ # CelerisLab/lbm/refinement_stepper.py -"""Opt-in fixed patch backend; no Simulation/objects adapter. +"""Opt-in fixed patch backend with optional stationary fine-interior circles. The caller owns the CudaContext and must close this backend before that context. Only synchronized postcollision_g leaf populations are public. No checkpoint, @@ -14,7 +14,8 @@ import pycuda.driver as cuda from ..cuda import compiler_v2 as compiler -from .refinement import ACTIVE, HALO, FixedPatch2D +from .refinement import ACTIVE, HALO, FixedPatch2D, prepare_fine_boundary +from .descriptors import SOLID @dataclass @@ -33,7 +34,11 @@ ``coarse_boundary`` is ``periodic`` (default) or ``slot_reservoir``: a bottom on-site slot and top fixed-density reservoir from ``cfg``. Periodic mode does not apply the wall/slot boundary fields in ``cfg``. - Physical walls remain coarse; fine halo donors must be ordinary fluid. + Outer walls remain coarse; fine halo donors must be ordinary fluid. + ``fine_objects`` accepts fixed circles in coarse physical coordinates; + velocities are in coarse lattice units (dx=dt=1). Independent local copies + are projected at fine cell centers. Their complete stencil must be active + and disjoint from the actual two-hop transfer footprint. ``initialize(coarse_g, fine_g)`` accepts physical Q-major (9,ny,nx) arrays already in postcollision phase. It does not collide or build channel flags. @@ -41,7 +46,7 @@ On a device/launch failure the instance becomes unusable: close and recreate it, rather than pretending the partially advanced fine buffer rolled back. """ - def __init__(self, cfg, bounds, *, context, coarse_boundary="periodic"): + def __init__(self, cfg, bounds, *, context, coarse_boundary="periodic", fine_objects=()): if HALO != 4: raise RuntimeError("refinement.cu requires a four-cell fine halo") if coarse_boundary not in ("periodic", "slot_reservoir"): @@ -55,6 +60,11 @@ "velocity=0 and 2 <= wall_slot_width < nx") if self.patch.bounds[1] < 4 or self.patch.bounds[3] > cfg.ny - 4: raise ValueError("slot_reservoir requires y0 >= 4 and y1 <= ny-4") + self.fine_boundary = prepare_fine_boundary( + cfg, self.patch, tuple(fine_objects), wall_aware=self._wall_aware) + self._fine_leaf = self.patch.fine_kind == ACTIVE + if self.fine_boundary is not None: + self._fine_leaf &= (self.fine_boundary.flags.reshape(self.patch.fine_shape) & SOLID) == 0 self.ctx = context self._lock = RLock() self._status = "new" @@ -62,6 +72,7 @@ self._allocations = [] self._levels = [] self._stream = None + self._fine_flags = self._fine_curved = self._curved = None self._block = (cfg.threads_per_block, 1, 1) fine_cfg = replace(cfg, nx=self.patch.fine_shape[1], ny=self.patch.fine_shape[0], viscosity=2 * cfg.viscosity) @@ -76,7 +87,9 @@ config.validate() if not np.isclose(config.omega, 1 / tau, rtol=1e-14, atol=0): raise RuntimeError("level omega/tau mismatch") - build = compiler.generate_config(config, n_objects=0) + n_objects = (len(self.fine_boundary.objects) + if name == "PatchFineStep" and self.fine_boundary is not None else 0) + build = compiler.generate_config(config, n_objects=n_objects) try: ptx = compiler.compile_kernel(self.ctx.sm_arch, source="refinement.cu", build_dir=build) @@ -86,7 +99,7 @@ ptr, size = module.get_global("d_params") if size != struct.calcsize("=ffffI"): raise RuntimeError("LBMParams host/device layout mismatch") - cuda.memcpy_htod(ptr, struct.pack("=ffffI", config.omega, 0, 0, 0, 0)) + cuda.memcpy_htod(ptr, struct.pack("=ffffI", config.omega, 0, 0, 0, n_objects)) shape = (9, config.ny, config.nx) size = int(np.prod(shape)) * 4 self._levels.append(_Level(shape, module, self._allocate(size), @@ -100,6 +113,16 @@ self._targets = self._upload(p.coalesce_coarse[order]) self._children = self._upload(p.coalesce_fine[order]) self._explode = self._levels[1].module.get_function("PatchExplode") + if self.fine_boundary is not None: + from .curved_links import CurvedLinkSoA + boundary = self.fine_boundary + self._fine_flags = self._upload(boundary.flags) + # SoA owns its allocations; never also register them in _allocations. + self._fine_curved = CurvedLinkSoA() + self._fine_curved.assign_host(*boundary.curved_host) + self._fine_curved.assign_surface(boundary.objects, fine_cfg, boundary.flags) + self._fine_curved.upload() + self._curved = self._levels[1].module.get_function("CurvedBoundaryKernel") except BaseException: try: self.close() @@ -131,7 +154,7 @@ coarse[:, y0:y1, x0:x1] = np.nan if self._wall_aware: coarse[:, (0, -1), :] = np.nan - fine[:, self.patch.fine_kind != ACTIVE] = np.nan + fine[:, ~self._fine_leaf] = np.nan def initialize(self, coarse_g, fine_g): """Upload an explicit G_0; invalid non-leaves are discarded.""" @@ -146,7 +169,7 @@ leaf[y0:y1, x0:x1] = False if self._wall_aware: leaf[[0, -1]] = False - for array, mask in zip(arrays, (leaf, self.patch.fine_kind == ACTIVE)): + for array, mask in zip(arrays, (leaf, self._fine_leaf)): if not np.isfinite(array[:, mask]).all() or (array[:, mask].sum(axis=0) <= 0).any(): raise ValueError("Leaf populations must be finite with positive density") self._mask_nonleaves(*arrays) @@ -154,11 +177,32 @@ for level, array in zip(self._levels, arrays): cuda.memcpy_htod(level.current, array) cuda.memcpy_htod(level.spare, array) + if self._fine_curved is not None and self._fine_curved.link_exchange_gpu is not None: + cuda.memset_d8(self._fine_curved.link_exchange_gpu, 0, self._fine_curved.count * 12) self._step_count = 0 self._status = "ready" except BaseException: self._status = "failed" raise + + def _apply_fine_boundary(self, buffer, **kw): + """Stage solid source slots; exchange buffer holds this fine substep only.""" + cl = self._fine_curved + if cl is None: + if self.fine_boundary is not None: + raise RuntimeError("Fine boundary device buffers are missing") + return + if self._curved is None or cl.fluid_idx_gpu is None: + raise RuntimeError("Fine boundary device buffers are missing") + driven = any(obj.surface_actuation is not None for obj in self.fine_boundary.objects) + if driven and (cl.surface_velocity_gpu is None or cl.link_exchange_gpu is None): + raise RuntimeError("Surface boundary device buffers are missing") + self._curved(buffer, cl.fluid_idx_gpu, cl.dir_gpu, cl.q_gpu, + cl.rx_gpu, cl.ry_gpu, cl.rz_gpu, cl.fallback_class_gpu, cl.body_id_gpu, + np.uintp(0), np.uintp(0), np.uint32(cl.count), + cl.surface_velocity_gpu if cl.surface_velocity_gpu is not None else np.uintp(0), + cl.link_exchange_gpu if cl.link_exchange_gpu is not None else np.uintp(0), + grid=((cl.count + self._block[0] - 1) // self._block[0], 1), **kw) def step(self, n=1): """Advance G→G: one explosion, two fine collisions, one coarse collision.""" @@ -176,8 +220,13 @@ np.uint32(len(p.ghost_fine)), np.uint64(np.prod(p.coarse_shape)), grid=((len(p.ghost_fine) + self._block[0] - 1) // self._block[0], 1), **kw) - f.kernel(f.current, f.spare, grid=f.grid, **kw) - f.kernel(f.spare, f.current, grid=f.grid, **kw) + flags = self._fine_flags if self._fine_flags is not None else np.uintp(0) + if self.fine_boundary is not None and self._fine_flags is None: + raise RuntimeError("Fine boundary flags are missing") + self._apply_fine_boundary(f.current, **kw) + f.kernel(f.current, f.spare, flags, grid=f.grid, **kw) + self._apply_fine_boundary(f.spare, **kw) + f.kernel(f.spare, f.current, flags, grid=f.grid, **kw) c.kernel(c.current, c.spare, f.current, self._targets, self._children, np.uint32(len(p.coalesce_coarse)), *(np.uint32(v) for v in p.bounds), np.uint32(self._wall_aware), @@ -220,6 +269,13 @@ error = exc self._levels.clear() self._explode = None + if getattr(self, "_fine_curved", None) is not None: + try: + self._fine_curved.free() + except Exception as exc: + error = error or exc + self._fine_curved = None + self._curved = self._fine_flags = None for pointer in self._allocations: try: pointer.free() --- a/src/CelerisLab/lbm/stepper.py +++ b/src/CelerisLab/lbm/stepper.py @@ -147,13 +147,17 @@ ) else: assert cl.fluid_idx_gpu is not None + if cl.surface_velocity is not None and ( + cl.surface_velocity_gpu is None or cl.link_exchange_gpu is None): + raise RuntimeError("Surface boundary device buffers are missing") self.curved_fn( f.ddf_gpu, cl.fluid_idx_gpu, cl.dir_gpu, cl.q_gpu, cl.rx_gpu, cl.ry_gpu, cl.rz_gpu, cl.fallback_class_gpu, cl.body_id_gpu, action_gpu, obs_gpu, np.uint32(f.n_curved), - cl.surface_velocity_gpu or np.uintp(0), cl.link_exchange_gpu or np.uintp(0), + cl.surface_velocity_gpu if cl.surface_velocity_gpu is not None else np.uintp(0), + cl.link_exchange_gpu if cl.link_exchange_gpu is not None else np.uintp(0), block=(tpb, 1, 1), grid=grid_c, **launch_kw, ) --- /dev/null +++ b/tests/unit/test_fine_surface.py @@ -0,0 +1,215 @@ +"""CPU-only fine circle projection and transfer/boundary isolation contracts.""" +from dataclasses import asdict, replace + +import numpy as np +import pytest + +from CelerisLab.body.geometry.circle import CircleGeometry +from CelerisLab.body.objects import SimObject, SurfaceActuation +from CelerisLab.config import LBMConfig +from CelerisLab.lbm import refinement +from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY, FLUID, SOLID +from CelerisLab.lbm.refinement import ACTIVE, INNER_GHOST, FixedPatch2D + + +def setup_patch(): + cfg = LBMConfig(nx=24, ny=24, viscosity=.05) + cfg.validate() + return cfg, FixedPatch2D.from_config(cfg, (4, 4, 20, 20)) + + +def circle(center=(12.25, 12.25), radius=2.15, speed=.002): + return SimObject(0, geometry=CircleGeometry(*center, radius), center=center, + radius=radius, surface_actuation=SurfaceActuation(.7, .85, speed)) + + +@pytest.mark.parametrize('speed', [.002, -.002]) +def test_fine_projection_uses_cell_center_bridge_and_preserves_source(speed): + cfg, patch = setup_patch() + source = circle(center=(12.25, 12.25 if speed > 0 else 11.75), speed=speed) + expected_center = (20., 20. if speed > 0 else 19.) + assert patch.fine_origin == (2., 2.) + before = (source.center, source.radius, source.geometry.center, + source.geometry._radius, asdict(source.state), asdict(source.surface_actuation), + vars(cfg).copy()) + plan = refinement.prepare_fine_boundary(cfg, patch, [source]) + local, = plan.objects + assert local is not source and local.geometry is not source.geometry + np.testing.assert_allclose(local.center, expected_center, rtol=0, atol=1e-14) + np.testing.assert_allclose(local.geometry.center, expected_center, rtol=0, atol=1e-14) + assert local.radius == pytest.approx(4.3) + assert local.geometry._radius == pytest.approx(4.3) + assert local.surface_actuation.normal_velocity == speed # dt/dx is unchanged. + assert (source.center, source.radius, source.geometry.center, + source.geometry._radius, asdict(source.state), asdict(source.surface_actuation), + vars(cfg).copy()) == before + _, _, _, _, rx, ry, _, _ = plan.curved_host + np.testing.assert_allclose(np.hypot(rx, ry), 4.3, rtol=0, atol=6e-7) + # These world points lie on the ORIGINAL physical circle, including the +.5. + world_x = patch.fine_origin[0]+.5*(local.center[0]+rx.astype(float)+.5) + world_y = patch.fine_origin[1]+.5*(local.center[1]+ry.astype(float)+.5) + np.testing.assert_allclose(np.hypot(world_x-source.center[0], world_y-source.center[1]), + source.radius, rtol=0, atol=3e-7) + + +def test_fine_boundary_receivers_sources_and_reverse_coverage(): + cfg, patch = setup_patch() + plan = refinement.prepare_fine_boundary(cfg, patch, [circle()]) + ny, nx = patch.fine_shape + flags = plan.flags + assert flags.shape == (nx*ny,) and flags.dtype == np.uint16 + solid = (flags & SOLID) != 0 + assert solid.any() + assert np.all(patch.fine_kind.ravel()[solid] == ACTIVE) + receivers, directions, q, _, _, _, _, fallback = plan.curved_host + cx, cy = np.asarray(D2Q9_EX), np.asarray(D2Q9_EY) + sources = receivers.astype(np.int64)+cx[directions]+nx*cy[directions] + assert np.all((flags[receivers] & FLUID) != 0) + assert not np.any(flags[receivers] & SOLID) + assert np.all(solid[sources]) + assert np.all(patch.fine_kind.ravel()[receivers] == ACTIVE) + assert np.all(patch.fine_kind.ravel()[sources] == ACTIVE) + # A source population slot has a unique writer, not just a unique cell. + assert len(set(zip(sources.tolist(), directions.tolist()))) == len(receivers) + expected = set() + for cell in np.flatnonzero((patch.fine_kind.ravel() == ACTIVE) & ~solid): + y, x = divmod(int(cell), nx) + for direction in range(1, 9): + sx, sy = x+int(cx[direction]), y+int(cy[direction]) + if 0 <= sx < nx and 0 <= sy < ny and solid[sy*nx+sx]: + expected.add((int(cell), direction)) + assert set(zip(receivers.tolist(), directions.tolist())) == expected + donors = receivers[q < .5].astype(np.int64)-cx[directions[q < .5]]-nx*cy[directions[q < .5]] + assert np.all(patch.fine_kind.ravel()[donors] == ACTIVE) + assert not np.any(solid[donors]) + assert not fallback.any() + stencil = set(receivers.tolist()) | set(sources.tolist()) | set(donors.tolist()) + assert stencil == set(plan.boundary_cells.tolist()) + assert not set(plan.boundary_cells.tolist()) & set(plan.transfer_fine_cells.tolist()) + assert not set(np.flatnonzero(solid).tolist()) & set(plan.transfer_fine_cells.tolist()) + + +def test_transfer_dependencies_follow_two_pulls_and_expand_active_collision(): + cfg, patch = setup_patch() + plan = refinement.prepare_fine_boundary(cfg, patch, [circle()]) + ny, nx = patch.fine_shape + cells = nx*ny + # Closed form for exactly two hops, independent of the production traversal: + # the final inner ghost only pulls; its predecessor may have collided. + offsets = np.asarray(D2Q9_EX)+nx*np.asarray(D2Q9_EY) + visited = set() + for slot in patch.coalesce_fine.ravel(): + direction, target = divmod(int(slot), cells) + assert patch.fine_kind.ravel()[target] == INNER_GHOST + middle = target-int(offsets[direction]) + visited.update(((2, direction, target), (1, direction, middle))) + if patch.fine_kind.ravel()[middle] == ACTIVE: + visited.update((0, q, middle-int(offsets[q])) for q in range(9)) + else: + assert patch.fine_kind.ravel()[middle] == INNER_GHOST + visited.add((0, direction, middle-int(offsets[direction]))) + actual = {tuple(map(int, row)) for row in plan.transfer_reads} + assert actual == visited + assert len(plan.transfer_reads) == len(visited) + expected_fine = set(patch.ghost_fine.tolist()) | {cell for _, _, cell in visited} + expected_coarse = set(patch.ghost_coarse.tolist()) | set( + (patch.coalesce_coarse % (cfg.nx*cfg.ny)).tolist()) + assert set(plan.transfer_fine_cells.tolist()) == expected_fine + assert set(plan.transfer_coarse_cells.tolist()) == expected_coarse + for coarse_slot, children in zip(patch.coalesce_coarse, patch.coalesce_fine): + coarse_q, coarse_cell = divmod(int(coarse_slot), cfg.nx*cfg.ny) + coarse_y, coarse_x = divmod(coarse_cell, cfg.nx) + fine_q, fine_cell = np.divmod(children, cells) + np.testing.assert_array_equal(fine_q, coarse_q) + fine_y, fine_x = np.divmod(fine_cell, nx) + physical_centers = np.column_stack(( + patch.fine_origin[0]+.5*(fine_x+.5), + patch.fine_origin[1]+.5*(fine_y+.5))) + np.testing.assert_array_equal(physical_centers.mean(axis=0), + (coarse_x+.5, coarse_y+.5)) + # This case must exercise the collision fan-in, not only ghost straight lines. + active_middle = [(q, cell) for stage, q, cell in visited + if stage == 1 and patch.fine_kind.ravel()[cell] == ACTIVE] + assert active_middle + _, sample = active_middle[0] + for incoming in range(9): + source = sample-int(D2Q9_EX[incoming])-nx*int(D2Q9_EY[incoming]) + assert (0, incoming, source) in actual + + +def test_fine_plan_arrays_are_readonly_and_contiguous(): + cfg, patch = setup_patch() + plan = refinement.prepare_fine_boundary(cfg, patch, [circle()]) + assert len(plan.curved_host) == 8 + for array in (plan.flags, *plan.curved_host, plan.boundary_cells, + plan.transfer_fine_cells, plan.transfer_coarse_cells, plan.transfer_reads): + assert array.flags.c_contiguous and not array.flags.writeable + assert array.size + with pytest.raises(ValueError): + array.flat[0] = array.flat[0] + + +def test_empty_objects_keep_existing_refinement_path(): + cfg, patch = setup_patch() + assert refinement.prepare_fine_boundary(cfg, patch, []) is None + + +def test_fine_rejects_cross_level_low_q_donor_without_fallback(): + cfg, patch = setup_patch() + source = circle(center=(5.05, 12.25), radius=.7) + # In fine coordinates C=(5.6,20), R=1.4. Receiver (4,20), dir +x + # has q=.2, but its required reverse donor (3,20) is not active. + local = CircleGeometry(5.6, 20., 1.4) + links = local.build_cut_links(patch.fine_shape[1], patch.fine_shape[0]) + target, = [link for link in links if link.fluid_idx == 20*40+4 and link.dir == 1] + assert target.q == pytest.approx(.2) + assert patch.fine_kind[20, 3] != ACTIVE + with pytest.raises(ValueError): + refinement.prepare_fine_boundary(cfg, patch, [source]) + + +def test_fine_rejects_active_stencil_that_still_touches_transfer(): + cfg, patch = setup_patch() + source = circle(center=(6.05, 12.25), radius=.7) + # Move the cross-level example two fine cells inward: all boundary reads + # are now active, but the q=.2 reverse donor is still a transfer dependency. + local = CircleGeometry(7.6, 20., 1.4) + nx = patch.fine_shape[1] + solid = (local.build_flag_mask(nx, patch.fine_shape[0]) & SOLID) != 0 + assert np.all(patch.fine_kind.ravel()[solid] == ACTIVE) + links = local.build_cut_links(nx, patch.fine_shape[0]) + stencil = set() + for link in links: + offset = int(D2Q9_EX[link.dir])+nx*int(D2Q9_EY[link.dir]) + stencil.update((link.fluid_idx, link.fluid_idx+offset)) + if link.q < .5: + stencil.add(link.fluid_idx-offset) + assert np.all(patch.fine_kind.ravel()[list(stencil)] == ACTIVE) + with pytest.raises(ValueError, match='transfer footprint'): + refinement.prepare_fine_boundary(cfg, patch, [source]) + + +@pytest.mark.parametrize('change', [dict(nx=25), dict(collision='TRT'), dict(viscosity=.06)]) +def test_fine_rejects_mismatched_patch_configuration(change): + cfg, patch = setup_patch() + with pytest.raises(ValueError, match='matching'): + refinement.prepare_fine_boundary(replace(cfg, **change), patch, [circle()]) + + +def test_fine_rejects_nonzero_control_even_without_motion(): + cfg, patch = setup_patch() + obj = circle() + obj.control.ax = .001 + with pytest.raises(ValueError, match='stationary'): + refinement.prepare_fine_boundary(cfg, patch, [obj]) + + +@pytest.mark.parametrize('objects', [ + [circle(center=(4.25, 12.25), radius=1.)], + [circle(center=(12.5, 12.5), radius=.05)], + [circle(), circle(center=(13., 12.25))], +]) +def test_fine_rejects_ghost_solids_unresolved_circle_and_overlap(objects): + cfg, patch = setup_patch() + with pytest.raises(ValueError): + refinement.prepare_fine_boundary(cfg, patch, objects) --- /dev/null +++ b/tests/validation/multigrid/test_fine_surface_gpu.py @@ -0,0 +1,330 @@ +"""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, BOUNDS, context=context, fine_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.fine_boundary = SimpleNamespace(objects=[object_spec(.002)]) + 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, BOUNDS, context=context, fine_objects=[source]) + context.activate.assert_not_called() + allocate.assert_not_called()