diff --git a/docs/development/raw_exchange.md b/docs/development/raw_exchange.md new file mode 100644 index 0000000..169016b --- /dev/null +++ b/docs/development/raw_exchange.md @@ -0,0 +1,65 @@ +# Raw multilevel boundary exchange (opt-in candidate) + +Use `FixedPatchStepper(..., observe_raw_exchange=True)` and, after a successful +base step, `read_raw_exchange(body_id, level=l, window="last_base_step")`. +The alternative window is `"last_substep"`. Body IDs are object positions within +the named level, not the object's original user label. The torque reference is +that circle's geometry.center (the center used to construct BFL arms), reported +in the root edge frame. A differing SimObject.center does not change this reference. Current geometry +restrictions are unchanged: fixed complete bodies on the finest level, isolated +from the transfer footprint. New observation currently requires blocked body +storage; dense-body configurations raise before device allocation. Existing +dense diagnostics are unchanged. + +Each result contains `wall`, `aperture`, `total`, with `force`, `torque`, +`impulse`, `angular_impulse`, `mass_into_fluid`, and `mass_rate`. These are raw +boundary momentum reaction, not validated stress or whole-nozzle mechanical +load. No GI/advective correction or guessed thrust is added. Apertures retain +their geometric classification when their prescribed velocity is zero. + +All results are base-lattice units per unit out-of-plane depth. For local +spacing and timestep h=dt, raw `(jx,jy,tz,dm)` converts to impulse +`(h*dt*jx,h*dt*jy,h*h*dt*tz)` and mass `h*h*dm`. +Force/torque are impulse/angular impulse divided by the reported window length. +`dm=b-a` is positive into the exterior fluid, while j is the opposite of the +fluid's momentum change across that cut link. Interpolated BFL can introduce +mass drift; dm is an observation, not an assertion of impermeability. + +`last_base_step` means [t-1,t] even after step(n); `last_substep` means [t-h,t] +and is not a common instantaneous sample across levels. A Jet amplitude set +before step(1) applies to all its substeps. Reads are repeatable; step(0) leaves +them unchanged. Initialization and a remap candidate invalidate observation. +Only successful stream synchronization publishes a window; failed instances +reject reads rather than returning partially advanced measurements. + +Storage is one FP32 [N,4] raw array and one FP64 [N,4] integral per measured +level, plus row metadata. No substep history is stored. Same-stream sequencing +is halo/restriction -> BFL -> independent row accumulation -> step. There are +no accumulation atomics or body obs offsets (blocked n_objects remains zero). +Disabled observation allocates neither observation buffer and launches no +accumulation kernel. The standalone periodic block caller passes explicit +null optional arguments. + +## Gates + +CPU (standard library unittest, plus existing NumPy dependency): + + PYTHONPATH=src python -m unittest discover -s tests/unit -p test_raw_exchange.py -v + +GPU, explicitly on a CUDA host with the existing environment: + + PYTHONPATH=src python -m unittest discover -s tests/validation/multigrid -p test_raw_exchange_gpu.py -v + +The GPU suite includes an independent actual-population cut-link census and +ACTIVE-fluid momentum/mass ledger with non-body boundary flux, observation +on/off population equivalence, dense/blocked row comparison, three block +partitions, zero/step/reversed/disabled Jet amplitudes, depth-four eight-sample +accumulation, and a short-time whole-domain uniform-fine comparison with the +same physical geometry, viscosity and Jet velocity. The uniform reference +uses a test-only adapter to the standalone kernel's optional inputs. + +These are implementation checks. They do not establish physical Jet/stress +accuracy, long-time AMR error, or convergence order. Resolution/long-time +studies and independently validated continuum load benchmarks remain separate. +GPU gates were prepared but not run in the authoring environment (no CUDA, +PyCUDA or nvcc); run them before merging this candidate. diff --git a/src/CelerisLab/lbm/block_stepper.py b/src/CelerisLab/lbm/block_stepper.py index 384d396..44a1109 100644 --- a/src/CelerisLab/lbm/block_stepper.py +++ b/src/CelerisLab/lbm/block_stepper.py @@ -211,6 +211,7 @@ class PeriodicBlockStepper: t = self._tables self._curved(self._current, t["block"], t["cell"], t["dir"], t["q"], t["fallback"], np.uint32(self._links), np.uintp(0), # no actuation + np.uintp(0), np.uintp(0), np.uintp(0), # no observation grid=self._link_grid, **kw) self._step(self._current, self._spare, np.uint32(self.layout.halo), self._tables.get("flags", np.uintp(0)), np.uintp(0), # no kind diff --git a/src/CelerisLab/lbm/kernels/blocks.cu b/src/CelerisLab/lbm/kernels/blocks.cu index 55b1849..739f1fb 100644 --- a/src/CelerisLab/lbm/kernels/blocks.cu +++ b/src/CelerisLab/lbm/kernels/blocks.cu @@ -76,11 +76,12 @@ extern "C" __global__ void BlockCoalesce( // cell, so the existing BFL writes the solid-source window image that this // block's next pull reads (no halo copy between this kernel and BlockStep). // surface_velocity (nullable, 2 floats per link): prescribed fluid velocity on -// fixed geometry, same convention as CurvedBoundaryKernel. No exchange output. +// fixed geometry, same convention as CurvedBoundaryKernel. Optional raw rows +// (jx,jy,torque,dm) are unique physical links, never body-observation slots. extern "C" __global__ void BlockCurvedBoundary( fpxx* f, const uint32_t* block, const uint32_t* cell, const unsigned char* dir, const float* q, const unsigned char* fallback, unsigned int count, - const float* surface_velocity) + const float* surface_velocity, const float* rx, const float* ry, float* raw) { unsigned int row = blockIdx.x * blockDim.x + threadIdx.x; if (row >= count) return; @@ -90,7 +91,24 @@ extern "C" __global__ void BlockCurvedBoundary( if (u != 0.0f || v != 0.0f) { Uw = u; Vw = v; } } apply_bouzidi_link(dir[row], q[row], cell[row], f + (unsigned long)block[row] * NQ * TOTAL_CELLS, - Uw, Vw, 0.0f, 0.0f, fallback[row], nullptr, -1); + Uw, Vw, rx ? rx[row] : 0.0f, ry ? ry[row] : 0.0f, + fallback[row], nullptr, -1, raw ? raw + 4ul * row : nullptr, + raw ? raw + 4ul * row + 3 : nullptr); +} + +// Same-stream BFL -> this row-owned accumulation -> pull step. No atomics, +// no cross-row reduction and no history allocation proportional to substeps. +extern "C" __global__ void BlockAccumulateExchange( + const float* raw, double* integral, unsigned int count, double h, double dt) +{ + unsigned int row = blockIdx.x * blockDim.x + threadIdx.x; + if (row >= count) return; + const unsigned long k = 4ul * row; + const double impulse = h * dt; + integral[k] += impulse * (double)raw[k]; + integral[k + 1] += impulse * (double)raw[k + 1]; + integral[k + 2] += h * impulse * (double)raw[k + 2]; + integral[k + 3] += h * h * (double)raw[k + 3]; } // Pull/collide/store core cells only; halo cells of out are left untouched. diff --git a/src/CelerisLab/lbm/kernels/boundary/curved_boundary.cuh b/src/CelerisLab/lbm/kernels/boundary/curved_boundary.cuh index 6b566f8..9bfab5c 100644 --- a/src/CelerisLab/lbm/kernels/boundary/curved_boundary.cuh +++ b/src/CelerisLab/lbm/kernels/boundary/curved_boundary.cuh @@ -82,7 +82,8 @@ __device__ inline void apply_bouzidi_link( float Uw, float Vw, float rx, float ry, unsigned char fallback_class, - float* __restrict__ obs, int body_id, float* link_exchange = nullptr) + float* __restrict__ obs, int body_id, float* link_exchange = nullptr, + float* link_mass = nullptr) { const unsigned int dir_opp = (unsigned int)opp_dir((int)dir); @@ -120,6 +121,10 @@ __device__ inline void apply_bouzidi_link( // corrected incoming population into the adjacent fluid node. store_ddf(fi, index_f(k_s, dir_opp), f_reflected); + // Equal opposite weights cancel even with shifted DDFs. + if (link_mass != nullptr) + *link_mass = f_reflected - f_toward; + if (obs != nullptr || link_exchange != nullptr) { #if USE_DDF_SHIFTING const float w_i = (float)d_w[dir]; diff --git a/src/CelerisLab/lbm/raw_exchange.py b/src/CelerisLab/lbm/raw_exchange.py new file mode 100644 index 0000000..d19fbea --- /dev/null +++ b/src/CelerisLab/lbm/raw_exchange.py @@ -0,0 +1,36 @@ +"""Host-only reduction of raw cut-link observations; no stress/GI correction. + +Rows are unique physical links. All quantities are 2D, per unit out-of-plane +thickness; centers are in the root edge frame, torque is about the named body. +""" +import numpy as np + + +def summarize_exchange(rows, body_ids, aperture, *, body_id, level, h, + time_end, window, reference_center): + """Reduce FP64 interval integrals or convert FP32 last-substep raw rows. + + dt_l=h_l under acoustic 2:1 refinement. Local momentum j -> base force + h*j -> impulse h*dt*j; torque adds one h; mass is h**2*dm. + """ + if window not in ("last_base_step", "last_substep"): + raise ValueError("unknown exchange window") + values = np.array(rows, dtype=np.float64, copy=True) + duration = 1.0 if window == "last_base_step" else h + if window == "last_substep": + values *= np.array([h*h, h*h, h*h*h, h*h]) + selected = np.asarray(body_ids) == body_id + aperture = np.asarray(aperture, bool) + + def reduce(mask): + integral = values[mask].sum(axis=0, dtype=np.float64) + return dict(impulse=integral[:2].copy(), angular_impulse=float(integral[2]), + mass_into_fluid=float(integral[3]), force=integral[:2] / duration, + torque=float(integral[2] / duration), mass_rate=float(integral[3] / duration)) + + return dict(wall=reduce(selected & ~aperture), aperture=reduce(selected & aperture), + total=reduce(selected), body_id=int(body_id), level=int(level), + time_start=time_end-duration, time_end=time_end, dt=h, + sample_count=round(1/h) if window == "last_base_step" else 1, + window=window, units="base_lattice_per_unit_depth", + definition="raw_boundary_exchange", reference_center=tuple(reference_center), valid=True) diff --git a/src/CelerisLab/lbm/refinement_stepper.py b/src/CelerisLab/lbm/refinement_stepper.py index aeb40da..1c8155a 100644 --- a/src/CelerisLab/lbm/refinement_stepper.py +++ b/src/CelerisLab/lbm/refinement_stepper.py @@ -23,6 +23,7 @@ from .refinement import ACTIVE, HALO, FixedPatch2D, FixedPatchSet2D, prepare_fin from .refinement_tables import (blocked_base_tables, blocked_fine_tables, # noqa: F401 (re-export) blocked_kind, _prepare_blocked_boundary) from .descriptors import SOLID +from .raw_exchange import summarize_exchange @dataclass @@ -158,8 +159,11 @@ class FixedPatchStepper: ``base_block=(by, bx)`` also blocks the parent; it requires ``fine_block``. Parent halos precede receiver-local coalescence staging and pull/collision; covered cells are explicitly poisoned on every step. Public shapes stay dense. - Blocked mode applies prescribed surface velocity only; its per-link - exchange diagnostic is unavailable (rows keep physical cut-link order). + ``observe_raw_exchange=True`` records blocked physical cut-link exchange + without changing populations. read_raw_exchange reports the last complete + base interval or last substep in base-lattice units per unit depth. This + is boundary momentum reaction, not a validated nozzle/stress load. Dense + body layouts are not supported by this new diagnostic in this candidate. ``jets`` (body.jets.Jet) use the same edge frame as ``fine_objects`` and drive fine cut links selected at their true hits; ``set_jet_amplitude`` re-uploads the per-link velocity column only. @@ -203,7 +207,8 @@ class FixedPatchStepper: def __init__(self, cfg, bounds=None, *, context, coarse_boundary="periodic", fine_objects=(), fine_block=None, base_block=None, jets=(), child=None, child_block=None, child_objects=(), child_jets=(), child_rectangles_in_parent=None, - refinements=None, level_blocks=None, level_objects=None, level_jets=None): + refinements=None, level_blocks=None, level_objects=None, level_jets=None, + observe_raw_exchange=False): if HALO != 4: raise RuntimeError("refinement.cu requires a four-cell fine halo") if coarse_boundary not in ("periodic", "slot_reservoir"): @@ -272,6 +277,7 @@ class FixedPatchStepper: # Per-level boundary geometry in the root edge frame (level l spacing 2**-l). self._boundary_hosts = [None] origin, spacing = (0.0, 0.0), 1.0 + spacings = [spacing] for level in range(1, depth): pair = pairs[level - 1] if level == 1: @@ -284,6 +290,17 @@ class FixedPatchStepper: self._boundary_hosts.append(host) origin = tuple(o + spacing * f for o, f in zip(origin, pair.fine_origin)) spacing /= 2 + spacings.append(spacing) + if not isinstance(observe_raw_exchange, bool): + raise ValueError("observe_raw_exchange must be bool") + if observe_raw_exchange and any(b is not None and layout is None + for b, layout in zip(self._boundary_hosts, self._layouts)): + raise ValueError("raw exchange currently requires blocked body levels") + self._observe_raw_exchange = observe_raw_exchange + self._exchange_valid = False + self._exchange_spacing = tuple(spacings) + self._exchange_centers = tuple(tuple(tuple(o.geometry.center) for o in objects) + for objects in level_objects) self._level_jets = level_jets self._level_boundaries = [None] + [ _prepare_blocked_boundary(self._boundary_hosts[level], self._layouts[level], @@ -318,7 +335,15 @@ class FixedPatchStepper: b = level.boundary if b is not None: b.flags = self._upload(b.tables["flags"]) - b.device = {k: self._upload(v) for k, v in b.tables.items() if k != "flags"} + columns = ("block", "cell", "dir", "q", "fallback", "surface_velocity") + b.device = {k: self._upload(v) for k, v in b.tables.items() if k in columns} + if self._observe_raw_exchange: + b.device.update(rx=self._upload(b.tables["rx"]), + ry=self._upload(b.tables["ry"])) + count = len(b.tables["cell"]) + b.last_raw = self._allocate(count * 4 * 4) + b.interval_integral = self._allocate(count * 4 * 8) + b.accumulate = level.module.get_function("BlockAccumulateExchange") b.kernel = level.module.get_function("BlockCurvedBoundary") if self._level_boundaries[1] is not None: b = self._level_boundaries[1] @@ -470,6 +495,7 @@ class FixedPatchStepper: 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._exchange_valid = False self._status = "ready" except BaseException: self._status = "failed" @@ -544,7 +570,8 @@ class FixedPatchStepper: cfg, pairs[0].bounds, context=self.ctx, base_block=base_block, fine_block=fine_block, child_block=child_block, child_rectangles_in_parent=new.child.bounds, - child_objects=child_objects, child_jets=jets) + child_objects=child_objects, child_jets=jets, + observe_raw_exchange=self._observe_raw_exchange) try: candidate.initialize(*transferred.g) candidate._step_count = count @@ -625,6 +652,8 @@ class FixedPatchStepper: raise RuntimeError("Surface velocity device buffer is missing") b.kernel(buffer, t["block"], t["cell"], t["dir"], t["q"], t["fallback"], np.uint32(count), velocity if velocity is not None else np.uintp(0), + t.get("rx", np.uintp(0)), t.get("ry", np.uintp(0)), + b.last_raw if b.last_raw is not None else np.uintp(0), grid=((count + self._block[0] - 1) // self._block[0], 1), **kw) def _advance(self, index, parent_buffer, **kw): @@ -664,6 +693,13 @@ class FixedPatchStepper: np.uint32(pair.n_rows), grid=((pair.n_rows + self._block[0] - 1) // self._block[0], 1), **kw) self._blocked_boundary(level, src, **kw) + b = level.boundary + if b is not None and b.last_raw is not None: + count = len(b.tables["cell"]) + h = self._exchange_spacing[index] # acoustic scaling: dt_l = h_l + b.accumulate(b.last_raw, b.interval_integral, np.uint32(count), + np.float64(h), np.float64(h), + grid=((count + self._block[0] - 1) // self._block[0], 1), **kw) flags = np.uintp(0) if level.boundary is None else level.boundary.flags level.kernel(src, dst, np.uint32(level.layout.halo), flags, level.kind, grid=level.grid, **kw) @@ -684,15 +720,57 @@ class FixedPatchStepper: try: self._status = "advancing" for _ in range(n): + if self._observe_raw_exchange: + self._exchange_valid = False + for b in self._level_boundaries: + if b is not None: + cuda.memset_d8_async(b.interval_integral, 0, + len(b.tables["cell"]) * 32, self._stream) self._hop(0, root.current, root.spare, **kw) self._stream.synchronize() root.current, root.spare = root.spare, root.current self._step_count += 1 + self._exchange_valid = self._observe_raw_exchange self._status = "ready" except BaseException: self._status = "failed" raise + def read_raw_exchange(self, body_id, *, level=1, window="last_base_step"): + """Read raw boundary reaction about its geometric center in the root edge frame. + + Repeated reads and step(0) preserve the window. last_substep samples + [t-h,t], not a common instantaneous time across refinement levels. + step(n) retains only the final complete base interval [t-1,t]. + Wall/aperture splitting is geometric, including zero-amplitude Jets. + """ + with self._lock, self.ctx.activate(): + self._require("ready") + if not self._observe_raw_exchange: + raise RuntimeError("raw exchange observation is disabled") + if not self._exchange_valid: + raise RuntimeError("raw exchange requires a completed base interval") + if (isinstance(level, bool) or not isinstance(level, Integral) + or not 0 <= level < len(self._levels)): + raise ValueError("level must name an existing level") + b = self._level_boundaries[level] + if (b is None or isinstance(body_id, bool) or not isinstance(body_id, Integral) + or not 0 <= body_id < len(b.host.objects)): + raise ValueError("body_id must name a body on this level") + if window not in ("last_base_step", "last_substep"): + raise ValueError("window must be last_base_step or last_substep") + count = len(b.tables["cell"]) + rows = np.empty((count, 4), np.float64 if window == "last_base_step" else np.float32) + try: + cuda.memcpy_dtoh(rows, b.interval_integral if window == "last_base_step" else b.last_raw) + except BaseException: + self._status = "failed" + raise + return summarize_exchange(rows, b.tables["body_id"], b.tables["aperture"], + body_id=body_id, level=level, h=self._exchange_spacing[level], + time_end=self._step_count, window=window, + reference_center=self._exchange_centers[level][body_id]) + def download(self): """Return synchronized G arrays (one per level), NaN on all non-leaves; no host cache.""" with self._lock, self.ctx.activate(): @@ -727,6 +805,7 @@ class FixedPatchStepper: for b in getattr(self, "_level_boundaries", ()): if b is not None: b.flags = b.device = b.kernel = None + b.last_raw = b.interval_integral = b.accumulate = None self._levels.clear() self._pairs = [] if getattr(self, "_fine_curved", None) is not None: diff --git a/src/CelerisLab/lbm/refinement_tables.py b/src/CelerisLab/lbm/refinement_tables.py index 5c7a452..0f8d5f3 100644 --- a/src/CelerisLab/lbm/refinement_tables.py +++ b/src/CelerisLab/lbm/refinement_tables.py @@ -77,6 +77,9 @@ class _BlockedBoundary: flags: object = None device: object = None kernel: object = None + last_raw: object = None + interval_integral: object = None + accumulate: object = None def _prepare_blocked_boundary(boundary, layout, cfg, jets): @@ -84,11 +87,17 @@ def _prepare_blocked_boundary(boundary, layout, cfg, jets): return None fluid_idx, dirs, q, body_id, rx, ry, _, fallback = boundary.curved_host tables, _ = blocked_link_tables(layout, boundary.flags, fluid_idx, dirs, q, fallback) + # Physical row identity excludes body_id: duplicate owners are always wrong. + keys = np.column_stack((fluid_idx, dirs)) + if len(np.unique(keys, axis=0)) != len(keys): + raise ValueError("duplicate physical boundary cut-link") + tables.update(body_id=np.asarray(body_id, np.int32), rx=np.asarray(rx, np.float32), + ry=np.asarray(ry, np.float32), aperture=np.zeros(len(dirs), bool)) actuated = bool(jets) or any(o.surface_actuation is not None for o in boundary.objects) if actuated: cfg = replace(cfg) cfg.validate() - tables["surface_velocity"], _ = surface_link_velocity( + tables["surface_velocity"], tables["aperture"] = surface_link_velocity( boundary.objects, cfg, boundary.flags, fluid_idx, dirs, body_id, rx, ry, jets=jets, hits=boundary.curved_hits) return _BlockedBoundary(boundary, tables, jets, actuated, diff --git a/tests/unit/test_raw_exchange.py b/tests/unit/test_raw_exchange.py new file mode 100644 index 0000000..28f8e70 --- /dev/null +++ b/tests/unit/test_raw_exchange.py @@ -0,0 +1,171 @@ +"""CPU reduction and real host scheduling with CUDA doubles, not GPU physics.""" +from contextlib import nullcontext +from dataclasses import replace +from unittest.mock import Mock, patch + +import numpy as np +import unittest + +from CelerisLab.body.geometry.circle import CircleGeometry +from CelerisLab.body.objects import SimObject, SurfaceActuation +from CelerisLab.config import LBMConfig +from CelerisLab.lbm.blocks import BlockLayout2D +from CelerisLab.lbm.refinement import FixedPatch2D, prepare_fine_boundary +from CelerisLab.lbm.refinement_tables import _prepare_blocked_boundary +from CelerisLab.lbm.raw_exchange import summarize_exchange +from test_refinement_three_level_multipatch import host_module + + +def body(center=(16., 16.), radius=.5, speed=0.): + return SimObject(91, geometry=CircleGeometry(*center, radius), center=center, radius=radius, + surface_actuation=SurfaceActuation(.7, .85, speed)) + + +def check_weights_split_and_window(level): + h = 2. ** -level + rows = np.array([[2, -3, 5, 7], [-1, 4, -2, -6], [99, 99, 99, 99]], np.float32) + kw = dict(body_id=0, level=level, h=h, time_end=3, reference_center=(7, 8)) + last = summarize_exchange(rows, [0, 0, 1], [False, True, False], window='last_substep', **kw) + integral = rows.astype(float)*[h*h, h*h, h**3, h*h]/h + full = summarize_exchange(integral, [0, 0, 1], [False, True, False], window='last_base_step', **kw) + np.testing.assert_array_equal(last['total']['force'], np.array([1, 1])*h) + assert last['total']['torque'] == 3*h*h + assert last['total']['mass_into_fluid'] == h*h + for name in ('force', 'torque', 'mass_rate'): + np.testing.assert_array_equal(last['total'][name], full['total'][name]) + for name in ('force', 'impulse', 'torque', 'angular_impulse', 'mass_into_fluid'): + np.testing.assert_array_equal(last['total'][name], + np.asarray(last['wall'][name])+last['aperture'][name]) + assert last['time_start'] == 3-h and full['time_start'] == 2 + assert last['sample_count'] == 1 and full['sample_count'] == 2**level + assert full['reference_center'] == (7, 8) + + +def test_final_substep_is_not_interval_mean(): + kw = dict(body_id=0, level=1, h=.5, time_end=1, reference_center=(0, 0)) + last = summarize_exchange([[4, 0, 0, 0]], [0], [False], window='last_substep', **kw) + full = summarize_exchange([[1.5, 0, 0, 0]], [0], [False], window='last_base_step', **kw) + assert last['total']['force'][0] == 2 and full['total']['force'][0] == 1.5 + + +def test_host_metadata_zero_jet_and_rows(): + cfg = LBMConfig(nx=32, ny=32, viscosity=.1) + pair = FixedPatch2D.from_config(cfg, (4, 4, 28, 28)) + boundary = prepare_fine_boundary(cfg, pair, [body(radius=2.15)]) + local = replace(cfg, nx=pair.fine_shape[1], ny=pair.fine_shape[0], viscosity=.2) + local.validate() + layout = BlockLayout2D.build(*pair.fine_shape, 8, 8, periodic=(False, False)) + b = _prepare_blocked_boundary(boundary, layout, local, ()) + for key, idx in [('body_id', 3), ('rx', 4), ('ry', 5)]: + np.testing.assert_array_equal(b.tables[key], boundary.curved_host[idx]) + assert b.tables['aperture'].any() and (~b.tables['aperture']).any() + assert not b.tables['surface_velocity'].any() + # Classification does not disappear when prescribed velocity is zero. + assert len(b.tables['cell']) == len(boundary.curved_host[0]) + + +def backend(observe=True, obj=None): + host = host_module() + calls = [] + module = Mock() + module.get_global.return_value = (Mock(), 20) + module.get_function.side_effect = lambda name: Mock(side_effect=lambda *a, **k: calls.append((name, a, k))) + host.compiler.load_module.return_value = module + host.cuda.mem_alloc.side_effect = lambda n: Mock(size=n) + host.cuda.memcpy_dtoh.side_effect = lambda array, pointer: array.fill(0) + ctx = Mock(activate=nullcontext, sm_arch='sm_70') + run = host.FixedPatchStepper(LBMConfig(nx=48, ny=48, viscosity=.02), context=ctx, + refinements=[[(6, 6, 26, 26)], [(20, 20, 44, 44)], [(56, 56, 72, 72)]], + level_blocks=[(4, 4)]*4, level_objects=[(), (), (), (body() if obj is None else obj,)], + observe_raw_exchange=observe) + arrays = [np.full(l.shape, 1/9, np.float32) for l in run._levels] + run.initialize(*arrays) + return run, host, calls, arrays + + +def test_four_level_schedule_publication_and_fixed_storage(): + run, host, calls, arrays = backend() + b = run._level_boundaries[3] + rows = len(b.tables['cell']) + assert b.last_raw.size == rows*16 and b.interval_integral.size == rows*32 + with unittest.TestCase().assertRaisesRegex(RuntimeError, 'completed'): + run.read_raw_exchange(0, level=3) + run.step() + names = [name for name, _, _ in calls] + assert names.count('BlockAccumulateExchange') == 8 + for i, (name, args, kw) in enumerate(calls): + if name == 'BlockAccumulateExchange': + assert names[i-1] == 'BlockCurvedBoundary' and names[i+1] == 'BlockStep' + assert args[0] is b.last_raw and args[1] is b.interval_integral + assert args[3:] == (.125, .125) + assert kw['stream'] is run._stream + assert host.cuda.memset_d8_async.call_count == 1 + value = run.read_raw_exchange(0, level=3) + assert value['sample_count'] == 8 and value['time_start'] == 0 + run.step(0) + assert run.read_raw_exchange(0, level=3)['time_end'] == 1 + assert run.read_raw_exchange(0, level=3, window='last_substep')['time_start'] == .875 + run.step(2) + assert run.read_raw_exchange(0, level=3)['time_start'] == 2 + run.initialize(*arrays) + with unittest.TestCase().assertRaisesRegex(RuntimeError, 'completed'): + run.read_raw_exchange(0, level=3) + run.close() + + + +def test_torque_reference_uses_geometry_center_not_optional_object_center(): + for center in (None, (3., 4.)): + obj = SimObject(91, geometry=CircleGeometry(16., 16., .5), center=center, radius=.5) + run, _, _, _ = backend(obj=obj) + run.step() + result = run.read_raw_exchange(0, level=3) + assert result['reference_center'] == (16., 16.) + assert result['reference_center'] != tuple(obj.center) + run.close() + +def test_disabled_has_nulls_and_no_observation_allocation_or_launch(): + run, host, calls, _ = backend(False) + b = run._level_boundaries[3] + assert b.last_raw is b.interval_integral is None + run.step() + for name, args, _ in calls: + if name == 'BlockCurvedBoundary': + assert tuple(args[8:]) == (0, 0, 0) + assert all(n != 'BlockAccumulateExchange' for n, _, _ in calls) + host.cuda.memset_d8_async.assert_not_called() + with unittest.TestCase().assertRaisesRegex(RuntimeError, 'disabled'): + run.read_raw_exchange(0, level=3) + run.close() + + +def test_failed_sync_never_publishes_partial_interval(): + run, host, calls, _ = backend() + run.step() + run._stream.synchronize.side_effect = RuntimeError('device failure') + with unittest.TestCase().assertRaisesRegex(RuntimeError, 'device failure'): + run.step() + assert not run._exchange_valid + with unittest.TestCase().assertRaisesRegex(RuntimeError, 'failed'): + run.read_raw_exchange(0, level=3) + run._stream.synchronize.side_effect = None + run.close() + + +def test_dense_observation_rejected_before_device_access(): + host = host_module() + ctx = Mock() + with unittest.TestCase().assertRaisesRegex(ValueError, 'blocked body'): + host.FixedPatchStepper(LBMConfig(nx=32, ny=32, viscosity=.1), (4, 4, 28, 28), + context=ctx, fine_objects=[body(radius=2.15)], observe_raw_exchange=True) + ctx.activate.assert_not_called() + + +def test_all_depth_weights(): + for level in (1, 2, 3, 7): + check_weights_split_and_window(level) + + +def load_tests(loader, tests, pattern): + return unittest.TestSuite(unittest.FunctionTestCase(value) for name, value in globals().items() + if name.startswith("test_") and callable(value)) diff --git a/tests/validation/multigrid/test_raw_exchange_gpu.py b/tests/validation/multigrid/test_raw_exchange_gpu.py new file mode 100644 index 0000000..e08da43 --- /dev/null +++ b/tests/validation/multigrid/test_raw_exchange_gpu.py @@ -0,0 +1,264 @@ +"""GPU gates for raw exchange. Run explicitly; requires PyCUDA and CUDA. + +Independent census reads actual staged incoming/outgoing populations, including +non-body ACTIVE-boundary flux, then checks postcollision fluid momentum/mass. +This is a discrete implementation gate, NOT a physical nozzle-load validation. +""" +from dataclasses import replace +import unittest +import numpy as np +import pycuda.driver as cuda + +from CelerisLab.body.geometry.circle import CircleGeometry +from CelerisLab.body.objects import SimObject +from CelerisLab.body.jets import Jet, Tube, Parallel +from CelerisLab.config import LBMConfig +from CelerisLab.cuda.context import CudaContext +from CelerisLab.lbm.refinement_stepper import FixedPatchStepper +from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY, SOLID + +CX, CY = np.asarray(D2Q9_EX), np.asarray(D2Q9_EY) +OPP = np.array([np.flatnonzero((CX == -x) & (CY == -y))[0] for x, y in zip(CX, CY)]) +W = np.array([4/9] + [1/9]*4 + [1/36]*4, np.float32) +CFG = LBMConfig(nx=32, ny=24, viscosity=.1) +BOUNDS = (8, 6, 24, 18) + + +def circle(): + return SimObject(17, geometry=CircleGeometry(16., 12., 2.15), center=(16., 12.), radius=2.15) + + +def jet(): + # Asymmetric aperture avoids hiding exchange errors by paired-jet cancellation. + return Jet(Tube((16., 12.), 0., 3., .7), Parallel(0.), amplitude=0.) + + +def state(shape, ux=.013): + cu = CX[:, None, None]*ux + values = W[:, None, None]*(1+3*cu+4.5*cu*cu-1.5*ux*ux) + return np.broadcast_to(values, shape).astype(np.float32).copy() + + +def download(pointer, shape, dtype=np.float32): + result = np.empty(shape, dtype) + cuda.memcpy_dtoh(result, pointer) + return result + + +def totals(g, leaf): + values = g[:, leaf].astype(np.float64) + return np.array([(values*CX[:, None]).sum(), (values*CY[:, None]).sum(), values.sum()]) + + +def census(run, index, buffer): + """Independent incoming/outgoing face census on ACTIVE fluid cells. + + Reads the receiver's actual block window, not a freshly gathered owner image + (that would erase BFL/restriction staging). Classifies external flux solely + from geometry, without consulting raw observation arrays. + """ + level = run._levels[index] + layout, b = level.layout, level.boundary + g = download(buffer, layout.blocked_shape) + dense = layout.to_dense(g) + leaf = run._leaf[index] + ny, nx = leaf.shape + solid = (b.host.flags.reshape(ny, nx) & SOLID) != 0 + y, x = np.nonzero(leaf) + owner = layout.owner[y, x].astype(np.intp) + cells = int(np.prod(layout.local_shape)) + block, cell = np.divmod(owner, cells) + flat = g.reshape(layout.nblocks, 9, cells) + external = np.zeros(3, np.float64) + for q, (cx, cy) in enumerate(zip(CX, CY)): + sy, sx = y-cy, x-cx + incoming = flat[block, q, cell-cy*layout.local_shape[1]-cx].astype(np.float64) + entering = ~leaf[sy, sx] & ~solid[sy, sx] + external += np.array([cx, cy, 1])*incoming[entering].sum() + dy, dx = y+cy, x+cx + leaving = ~leaf[dy, dx] & ~solid[dy, dx] + outgoing = flat[block, q, cell].astype(np.float64) + external -= np.array([cx, cy, 1])*outgoing[leaving].sum() + receiver, direction, _, _, rx, ry, _, _ = b.host.curved_host + owner = layout.owner.ravel()[receiver].astype(np.intp) + block, cell = np.divmod(owner, cells) + a = flat[block, direction, cell].astype(np.float64) + source = cell+CY[direction]*layout.local_shape[1]+CX[direction] + incoming = flat[block, OPP[direction], source].astype(np.float64) + jx, jy = CX[direction]*(a+incoming), CY[direction]*(a+incoming) + independent = np.column_stack((jx, jy, rx*jy-ry*jx, incoming-a)) + return totals(dense, leaf), external, independent + + +class RawExchangeGPU(unittest.TestCase): + def test_row_ledger_interval_jet_switch_and_no_population_effect(self): + with CudaContext(0) as ctx: + for block in ((32, 40), (16, 20), (8, 10)): + with FixedPatchStepper(CFG, BOUNDS, context=ctx, fine_block=block, + fine_objects=[circle()], jets=[jet()], observe_raw_exchange=True) as measured, \ + FixedPatchStepper(CFG, BOUNDS, context=ctx, fine_block=block, + fine_objects=[circle()], jets=[jet()]) as plain: + arrays = [state(l.shape) for l in measured._levels] + measured.initialize(*arrays) + plain.initialize(*arrays) + b = measured._level_boundaries[1] + count = len(b.tables['cell']) + kernel = measured._levels[1].kernel + samples = [] + + def audited_step(src, dst, *args, **kw): + measured._stream.synchronize() + before, external, independent = census(measured, 1, src) + raw = download(b.last_raw, (count, 4)) + np.testing.assert_allclose(raw, independent, rtol=2e-6, atol=2e-7) + kernel(src, dst, *args, **kw) + measured._stream.synchronize() + after = measured._levels[1].layout.to_dense( + download(dst, measured._levels[1].layout.blocked_shape)) + change = totals(after, measured._leaf[1])-before + expected = external + np.array([-raw[:, 0].astype(float).sum(), + -raw[:, 1].astype(float).sum(), + raw[:, 3].astype(float).sum()]) + np.testing.assert_allclose(change, expected, rtol=2e-4, atol=8e-5) + samples.append(raw.copy()) + measured._levels[1].kernel = audited_step + for speed in (0., .002, -.001, 0.): + measured.set_jet_amplitude(0, speed) + plain.set_jet_amplitude(0, speed) + samples.clear() + measured.step() + plain.step() + for a, e in zip(measured.download(), plain.download()): + np.testing.assert_array_equal(a, e) + integral = download(b.interval_integral, (count, 4), np.float64) + expected = sum(s.astype(float)*[.25, .25, .125, .25] for s in samples) + np.testing.assert_array_equal(integral, expected) + last = download(b.last_raw, (count, 4)) + np.testing.assert_array_equal(last, samples[-1]) + obs = measured.read_raw_exchange(0) + assert obs['sample_count'] == 2 + assert b.tables['aperture'].any() + measured.step(0) + np.testing.assert_array_equal(measured.read_raw_exchange(0)['total']['impulse'], + 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, BOUNDS, context=ctx, fine_objects=[circle()], jets=[jet()]) as dense, \ + FixedPatchStepper(CFG, BOUNDS, context=ctx, fine_block=(16, 20), + fine_objects=[circle()], jets=[jet()], 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) + drive = Jet(Tube((16., 16.), 0., .75, .2), Parallel(0.), amplitude=0.) + settings = dict(refinements=[[(6, 6, 26, 26)], [(20, 20, 44, 44)], [(56, 56, 72, 72)]], + level_blocks=[(4, 4)]*4, level_objects=[(), (), (), (obj,)], + level_jets=[(), (), (), (drive,)]) + with CudaContext(0) as ctx: + with FixedPatchStepper(cfg, context=ctx, observe_raw_exchange=True, **settings) as run, \ + FixedPatchStepper(cfg, context=ctx, **settings) as plain: + arrays = [state(l.shape, ux=0.) for l in run._levels] + run.initialize(*arrays); plain.initialize(*arrays) + b = run._level_boundaries[3] + count = len(b.tables['cell']) + kernel = run._levels[3].kernel + samples, independent_samples = [], [] + def audited_step(src, dst, *args, **kw): + run._stream.synchronize() + before, external, independent = census(run, 3, src) + raw = download(b.last_raw, (count, 4)) + np.testing.assert_allclose(raw, independent, rtol=2e-6, atol=2e-7) + kernel(src, dst, *args, **kw) + run._stream.synchronize() + after = run._levels[3].layout.to_dense( + download(dst, run._levels[3].layout.blocked_shape)) + change = totals(after, run._leaf[3])-before + expected = external + np.array([-independent[:, 0].sum(), + -independent[:, 1].sum(), independent[:, 3].sum()]) + np.testing.assert_allclose(change, expected, rtol=2e-4, atol=8e-5) + samples.append(raw) + independent_samples.append(independent) + run._levels[3].kernel = audited_step + for speed in (0., .002, 0.): + run.set_jet_amplitude(0, speed, level=3) + plain.set_jet_amplitude(0, speed, level=3) + samples.clear(); independent_samples.clear() + run.step(); plain.step() + self.assertEqual(len(samples), 8) + weights = np.array([1/64, 1/64, 1/512, 1/64]) + integral = download(b.interval_integral, (count, 4), np.float64) + np.testing.assert_array_equal(integral, sum(v.astype(float)*weights for v in samples)) + np.testing.assert_allclose(integral, sum(v*weights for v in independent_samples), + rtol=2e-6, atol=2e-7) + np.testing.assert_array_equal(download(b.last_raw, (count, 4)), samples[-1]) + for actual, expected in zip(run.download(), plain.download()): + np.testing.assert_array_equal(actual, expected) + result = run.read_raw_exchange(0, level=3) + self.assertEqual(result['sample_count'], 8) + self.assertEqual(result['time_start'], run.step_count-1) + if run.step_count == 1: + np.testing.assert_allclose(result['total']['force'], 0, atol=2e-7) + self.assertAlmostEqual(result['total']['torque'], 0, delta=2e-7) + + def test_uniform_whole_domain_fine_short_time_reference(self): + """Same root geometry, nu and Jet speed on a whole-domain fine mesh. + Eight local hops are a short-time equivalence check, not convergence. + The test adapter only enables the standalone kernel's optional inputs. + """ + from CelerisLab.lbm.blocks import BlockLayout2D + from CelerisLab.lbm.block_stepper import PeriodicBlockStepper + fine_cfg = replace(CFG, nx=64, ny=48, viscosity=.2) + obj = SimObject(0, geometry=CircleGeometry(31.5, 23.5, 4.3), + center=(31.5, 23.5), radius=4.3) + layout = BlockLayout2D.build(48, 64, 8, 8) + with CudaContext(0) as ctx: + with PeriodicBlockStepper(fine_cfg, layout, context=ctx, objects=[obj]) as uniform, \ + FixedPatchStepper(CFG, BOUNDS, context=ctx, fine_block=(16, 20), + fine_objects=[circle()], jets=[jet()], observe_raw_exchange=True) as amr: + amr.initialize(*(state(l.shape) for l in amr._levels)) + uniform.initialize(state((9, 48, 64))) + b = amr._level_boundaries[1] + count = len(b.tables['cell']) + links = obj.geometry.build_cut_links(64, 48) + np.testing.assert_array_equal([v.dir for v in links], b.tables['dir']) + np.testing.assert_allclose([v.q for v in links], b.tables['q'], atol=2e-7, rtol=0) + # Same fine-grid force arms after translation to the same center. + rx = uniform._upload(b.tables['rx']) + ry = uniform._upload(b.tables['ry']) + velocity = uniform._upload(b.tables['surface_velocity']) + raw = uniform._allocate(count*16) + original = uniform._curved + samples = [] + def measured_boundary(*args, **kw): + original(*args[:7], velocity, rx, ry, raw, **kw) + uniform._stream.synchronize() + samples.append(download(raw, (count, 4))) + uniform._curved = measured_boundary + for speed in (0., .002, -.001, 0.): + amr.set_jet_amplitude(0, speed) + cuda.memcpy_htod(velocity, b.tables['surface_velocity']) + samples.clear() + amr.step(); uniform.step(2) + expected = sum(v.astype(float)*[.25, .25, .125, .25] for v in samples) + actual = download(b.interval_integral, (count, 4), np.float64) + np.testing.assert_allclose(actual, expected, rtol=1e-5, atol=2e-6) + + +if __name__ == '__main__': + unittest.main()