diff --git a/docs/lifecycle.md b/docs/lifecycle.md index 839d554..56694db 100644 --- a/docs/lifecycle.md +++ b/docs/lifecycle.md @@ -5,7 +5,7 @@ it is; it adds no framework. Line references are to candidate 0f8fdbc3. | Class | Contents | How it changes | Owner rule | |---|---|---|---| -| 1. Module and its globals | Compile-time (visible to the compiler, fixed per module): `config.h` macros incl. `NX`/`NY`, `COLLISION_MODEL`, `BLOCK_ROOT_BC`, `WALL_*`. Runtime per-module globals (not visible to the compiler; written by host code): `__constant__` `d_params`, `d_w`, `d_cx`, `d_cy`, `d_block_channel`, `d_block_wall`; plus bound kernel handles | Globals are written once by the owner right after `load_module`. Uniform grid: `Simulation._recompile` (construction, body-topology sync) loads a new module and re-uploads; `update_runtime_params` rewrites `d_params` (omega, force, n_objects) between runs. Steppers: never rewritten after construction, so they are immutable for the instance's lifetime although they are runtime values. The same quantity can be a macro in one route and a global in another (dense slot root: `WALL_SLOT_WIDTH`/`WALL_NORMAL_VELOCITY` macros; blocked slot root: `d_block_wall`), which can change generated arithmetic. | Every owner loads its own module (`compiler.load_module` = `module_from_file`, no cache). A loaded module is never shared between live owners: globals are per module, so sharing would share mutable `d_params`. | +| 1. Module and its globals | Compile-time (visible to the compiler, fixed per module): `config.h` macros incl. `NX`/`NY`, `COLLISION_MODEL`, `BLOCK_ROOT_BC`, `WALL_*`. Runtime per-module globals (not visible to the compiler; written by host code): `__constant__` `d_params`, `d_w`, `d_cx`, `d_cy`, `d_block_channel`, `d_block_wall`; plus bound kernel handles | Globals are written once by the owner right after `load_module`. Uniform grid: `Simulation._recompile` (construction, body-topology sync) loads a new module and re-uploads; `update_runtime_params` rewrites `d_params` (omega, force, n_objects) between runs. Steppers: never rewritten after construction, so they are immutable for the instance's lifetime although they are runtime values. The blocked slot/reservoir root reads its slot from the `d_block_wall` constant (the former dense route used `WALL_*` macros; both were bit-identical in the parity gates). | Every owner loads its own module (`compiler.load_module` = `module_from_file`, no cache). A loaded module is never shared between live owners: globals are per module, so sharing would share mutable `d_params`. | | 2. Geometry, layout, tables | flags, curved-link tables, block layouts and halo pairs, transfer (explode/coalesce) tables, block kinds | Immutable for the life of an instance. Uniform grid: rebuilt by body sync. `FixedPatchStepper`: a geometry change is a new instance (`remap_child`), committed only after a successful transfer. | No in-place table edits. | | 3. Controls and dependent boundary columns | Jet amplitude -> per-link `surface_velocity`; body actions (uniform grid) | Between steps, without a table or module rebuild (`set_jet_amplitude`). | Host copy commits only after the device upload succeeded; a device error fails the instance. | | 4. Per-step state | population double buffers and their parity, raw-exchange ledgers, observations | Device-only during `step(n)` / `run(n)`; buffer swaps are host-handle swaps. | Valid on the host only after the single wait at the end of the call. | @@ -15,10 +15,11 @@ Storage routes are separate concerns: - Original uniform grid (`uniform.cu`, `Simulation`): whole-domain `NX`/`NY`. - Blocked levels (`blocked.cu`; every `FixedPatchStepper` level, including the single level of `refinements=[]`): `NX`/`NY` are one block window including halo. -- Dense refinement levels (`refinement.cu`, `level_blocks[l] is None`). +- Dense refinement levels (`refinement.cu`) were removed: `level_blocks` is + required and every refinement level is blocked. A single-rectangle refinement (`FixedPatch2D`) is stored in its own fine -window frame in both dense and blocked form (`refinement_stepper.py:222-224`, +window frame (`FixedPatchStepper.__init__`, `BlockLayout2D.build(*pair.fine_shape, ...)`). Only a rectangle set (`FixedPatchSet2D`) uses the level-wide frame `(2*ny, 2*nx)`, whose `fine_layout` requires that full grid to be divisible by the tile diff --git a/docs/supported_settings.md b/docs/supported_settings.md index 2a7a323..f75e301 100644 --- a/docs/supported_settings.md +++ b/docs/supported_settings.md @@ -10,12 +10,13 @@ measured; they are not general bounds. | Setting | Supported | Guard | |---|---|---| -| Dimension / lattice | 2D D2Q9 only (3D refinement is future work) | `lbm/refinement.py` `FixedPatch2D.from_config`; `#error` guards in `lbm/kernels/blocks.cu` and `lbm/kernels/refinement.cu` | +| Dimension / lattice | 2D D2Q9 only (3D refinement is future work) | `lbm/refinement.py` `FixedPatch2D.from_config`; `#error` guards in `lbm/kernels/blocked.cu` | | Precision / streaming | FP32 storage, double buffer | same guards | | Collision on refinement levels | SRT or MRT; TRT and LES rejected | `lbm/refinement.py` `FixedPatch2D.from_config` (re-checked for fine boundaries in `prepare_fine_boundary`) | -| Shifted (ddf_shifting) storage | only when every level is blocked; dense or mixed rejected | `lbm/refinement_stepper.py` `FixedPatchStepper.__init__`, `FixedPatchStepper._check_storage`; `#error` guard in `lbm/kernels/refinement.cu` (dense patch kernel unshifted only) | +| Shifted (ddf_shifting) storage | supported (every level is blocked); one representation on all levels | `lbm/refinement_stepper.py` `FixedPatchStepper._check_storage` | | Fine relaxation | 1/tau_fine must lie inside the configured omega guards | `lbm/refinement.py` `FixedPatch2D.from_config` | | Single blocked level `FixedPatchStepper(refinements=[], level_blocks=[(by, bx)], level_objects=[circles])` (replaces `PeriodicBlockStepper`, removed in 8e1feff) | doubly periodic, SRT only, no LES/Jets/projection/exchange observation, stationary passive circles; `download()` returns a 1-tuple | `lbm/refinement_stepper.py` `_check_levels`, `FixedPatchStepper.__init__`; `lbm/blocks.py` `static_circle_tables` | +| Refinement storage | blocked on every level: `level_blocks` is required, one tile `(by, bx)` per level dividing that level; a whole-level tile (`(ny, nx)` at the root, the fine window shape for one rectangle) always divides. Dense and mixed dense/blocked storage were removed (parity of D/D, D/B and B/B recorded at 1573942 and 2469af7) | `lbm/refinement_stepper.py` `_check_levels`, `FixedPatchStepper.__init__` | Tests: `tests/unit/test_refinement_mrt.py`, `tests/validation/multigrid/test_refinement_mrt_gpu.py`, `tests/unit/test_shifted_storage.py`, `tests/unit/test_shifted_block_operators.py`. diff --git a/src/CelerisLab/cuda/compiler.py b/src/CelerisLab/cuda/compiler.py index 8701296..fd20a29 100644 --- a/src/CelerisLab/cuda/compiler.py +++ b/src/CelerisLab/cuda/compiler.py @@ -7,8 +7,7 @@ build_module runs, in a private directory it removes afterwards: 1. build_dir = generate_config(lbm_cfg, n_objects) — isolated source snapshot 2. compile_kernel(arch, source=entry, build_dir=build_dir) — nvcc -ptx 3. load_module(ptx_path) — PyCUDA module from PTX -Entries: uniform.cu (uniform grid), blocked.cu (blocked level), refinement.cu -(dense refinement level). +Entries: uniform.cu (uniform grid), blocked.cu (every blocked level, single or refined). """ from __future__ import annotations diff --git a/src/CelerisLab/lbm/kernels/refinement.cu b/src/CelerisLab/lbm/kernels/refinement.cu deleted file mode 100644 index 5d6ac9d..0000000 --- a/src/CelerisLab/lbm/kernels/refinement.cu +++ /dev/null @@ -1,40 +0,0 @@ -// Fixed, unforced internal 2:1 patch. Separate compilation entry; -// default OneStep/channel kernels are deliberately not included. -#include -#include -#include "config.h" -#include "core/precision.cuh" -#include "core/descriptors.cuh" -#include "core/layout.cuh" -#include "core/params.cuh" -#include "core/flags.cuh" -#include "core/obs.cuh" -#include "core/cell_kind.cuh" -#include "operators/macro.cuh" -#include "operators/collision_srt.cuh" -#if COLLISION_MODEL == 2 -#include "operators/collision_mrt.cuh" -#endif -#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 && COLLISION_MODEL != 2) || STORE_PRECISION != 0 || USE_DDF_SHIFTING || USE_LES || STREAMING_MODEL != 0 -#error "Refinement requires unshifted FP32 D2Q9 SRT or MRT double buffers without LES" -#endif - -// Enabled only in the validation build; no diagnostic arrays in normal builds. -#ifdef CELERIS_REFINEMENT_AUDIT -extern "C" { -__device__ unsigned int collision_counts[TOTAL_CELLS]; -__device__ unsigned int explosion_counts[TOTAL_CELLS]; -__device__ unsigned int incoming_counts[NQ * TOTAL_CELLS]; -} -#endif - -// Kernels: patch collision helper, PatchExplode, PatchFineStep, PatchCoarseStep. -#include "step/patch_step.cu" diff --git a/src/CelerisLab/lbm/kernels/step/block_transfer.cu b/src/CelerisLab/lbm/kernels/step/block_transfer.cu index 9b14fed..f6bdfd8 100644 --- a/src/CelerisLab/lbm/kernels/step/block_transfer.cu +++ b/src/CelerisLab/lbm/kernels/step/block_transfer.cu @@ -16,8 +16,8 @@ extern "C" __global__ void BlockHaloCopy( } // Fixed-patch explosion into blocked fine storage: target = owner block cell -// id of a patch-ghost cell; donors are parent block cell ids. Dense parent -// storage is the one-window case with coarse_cells = full domain size. +// id of a patch-ghost cell; donors are parent block cell ids and coarse_cells +// is the parent's block window size. extern "C" __global__ void BlockExplode( const fpxx* coarse, fpxx* fine, const uint64_t* target, const uint64_t* donor, unsigned int count, unsigned long coarse_cells) @@ -66,7 +66,7 @@ extern "C" __global__ void BlockExplodeProjected( } // Targets are receiver-local SOURCE slots for the immediately following pull. -// Keep child summation order identical to PatchCoarseStep; no collision here. +// Child summation order ((c0 + c1) + c2) + c3 is part of the numerics; no collision here. extern "C" __global__ void BlockCoalesce( const fpxx* fine, fpxx* coarse, const uint64_t* target, const uint64_t* children, unsigned int count) diff --git a/src/CelerisLab/lbm/kernels/step/patch_step.cu b/src/CelerisLab/lbm/kernels/step/patch_step.cu deleted file mode 100644 index 1f64ea7..0000000 --- a/src/CelerisLab/lbm/kernels/step/patch_step.cu +++ /dev/null @@ -1,95 +0,0 @@ -// CelerisLab/lbm/kernels/step/patch_step.cu -// Dense fixed-patch collision, explosion and the fine/coarse steps. -// Part of the refinement.cu compilation entry; relies on its includes, guards -// and the optional CELERIS_REFINEMENT_AUDIT counters declared there. - -__device__ __forceinline__ void patch_collide(float* f, unsigned long k) { -#if COLLISION_MODEL == 2 - collide_mrt_unforced(f); -#else - float rho, ux, uy, feq[NQ]; - compute_rho_u(f, rho, ux, uy); - compute_feq(rho, ux, uy, feq); - collide_srt_no_force(f, feq, d_params.omega); -#endif -#ifdef CELERIS_REFINEMENT_AUDIT - ++collision_counts[k]; -#endif -} - -extern "C" __global__ void PatchExplode( - const float* coarse, float* fine, const uint32_t* targets, - const uint32_t* donors, unsigned int count, unsigned long coarse_cells) -{ - unsigned int row = blockIdx.x * blockDim.x + threadIdx.x; - if (row >= count) return; - for (unsigned int q = 0; q < NQ; ++q) - fine[index_f(targets[row], q)] = coarse[q * coarse_cells + donors[row]]; -#ifdef CELERIS_REFINEMENT_AUDIT - ++explosion_counts[targets[row]]; -#endif -} - -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 - || (flags != nullptr && is_solid(flags[k]))) { - fill_invalid(f); - } else { - unsigned long j[NQ]; - compute_neighbors(x, y, j); // These interior targets never wrap. - stream_pull_load(k, f, in, j); - if (x >= 4 && y >= 4 && x < NX - 4 && y < NY - 4) - patch_collide(f, k); - } - stream_pull_store(k, f, out); -} - -extern "C" __global__ void PatchCoarseStep( - const float* in, float* out, const float* fine, - const uint64_t* targets, const uint64_t* children, unsigned int count, - unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1, - unsigned int wall_aware) -{ - unsigned int x, y; unsigned long k; - index_from_thread(x, y, k); - if (x >= NX || y >= NY) return; - // Both buffers were poisoned at initialization; covered cells stay invalid. - if (x >= x0 && x < x1 && y >= y0 && y < y1) return; - if (wall_aware && (y == 0 || y == NY - 1)) return; - unsigned long j[NQ]; - compute_neighbors(x, y, j); - float f[NQ]; - stream_pull_load(k, f, in, j); - if (x >= x0 - 1 && x <= x1 && y >= y0 - 1 && y <= y1) { - for (unsigned int q = 1; q < NQ; ++q) { - uint64_t slot = index_f(k, q); - unsigned int lo = 0, hi = count; - while (lo < hi) { - unsigned int mid = lo + (hi - lo) / 2; - if (targets[mid] < slot) lo = mid + 1; - else hi = mid; - } - if (lo < count && targets[lo] == slot) { - const uint64_t* s = children + 4ul * lo; - f[q] = 0.25f * (fine[s[0]] + fine[s[1]] + fine[s[2]] + fine[s[3]]); -#ifdef CELERIS_REFINEMENT_AUDIT - ++incoming_counts[slot]; -#endif - } - } - } - // Includes the eight grazing-corner links with UNCOVERED coarse donors. - // All replacement arrivals must be in place before computing rho/u. - if (wall_aware) { - if (y == 1) reconstruct_bottom_velocity(f, 0.0f, bottom_wall_normal_velocity(x)); - else if (y == NY - 2) reconstruct_top_reservoir(f, RHO); - } - patch_collide(f, k); - stream_pull_store(k, f, out); -} diff --git a/src/CelerisLab/lbm/refinement_stepper.py b/src/CelerisLab/lbm/refinement_stepper.py index 1fb7624..fd53ef7 100644 --- a/src/CelerisLab/lbm/refinement_stepper.py +++ b/src/CelerisLab/lbm/refinement_stepper.py @@ -19,7 +19,7 @@ from .blocks import (BlockLayout2D, launch_curved_boundary, launch_halo_copy, le static_circle_tables) from .channel import channel_masks, validate_channel_config, validate_channel_transfers from ..body.jets import jet_link_velocity -from .hierarchy import build_hierarchy, is_rectangle +from .hierarchy import build_hierarchy from .refinement import ACTIVE, HALO, FixedPatch2D, FixedPatchSet2D, prepare_fine_boundary from .refinement_tables import (blocked_base_tables, blocked_fine_tables, # noqa: F401 (re-export) blocked_kind, projected_explosion_tables, _prepare_blocked_boundary, @@ -75,21 +75,19 @@ def _check_levels(refinements, blocks, objects, jets): if len(values) != depth: raise ValueError(f"{name} needs one entry per level ({depth})") for level, block in enumerate(blocks): - if block is not None and len(block) != 2: - raise ValueError(f"level_blocks[{level}] must be (by, bx)") + if (not isinstance(block, tuple) or len(block) != 2 + or any(isinstance(v, (bool, np.bool_)) or not isinstance(v, Integral) or v <= 0 + for v in block)): + raise ValueError(f"level_blocks[{level}] must be a tile (by, bx) of positive integers; " + "a whole-level tile ((ny, nx) at the root, the fine window shape " + "for one rectangle) always divides its level") if depth == 1: # one blocked level: the root carries the static circles - if blocks[0] is None or jets[0]: - raise ValueError("a single level must be blocked and has no Jets") + if jets[0]: + raise ValueError("a single level has no Jets") return if objects[0] or jets[0]: raise ValueError("base-level bodies/Jets are not supported") - if any(blocks[level - 1] is not None and blocks[level] is None for level in range(1, depth)): - raise ValueError("a blocked level requires blocked refined levels below it") - if not is_rectangle(refinements[0]) and (blocks[0] is None or blocks[1] is None): - raise ValueError("several L1 patches require level_blocks[0] and level_blocks[1]") if depth > 2: - if any(block is None for block in blocks): - raise ValueError("three or more levels must all be blocked") if any(objects[level] or jets[level] for level in range(1, depth - 1)): raise ValueError("bodies/Jets are supported only on the finest level of a nested hierarchy") @@ -101,15 +99,16 @@ class FixedPatchStepper: (x0, y0, x1, y1) (FixedPatch2D, a fine window with its own origin) or a sequence of rectangles (FixedPatchSet2D, a dependency-validated union on one level-wide grid, public fine arrays (2*ny, 2*nx), NaN outside ACTIVE cells); - see lbm.hierarchy. ``level_blocks``, ``level_objects`` and ``level_jets`` - hold one entry per level, L0 included (L0 bodies/Jets are not supported). + see lbm.hierarchy. ``level_blocks`` (required), ``level_objects`` and + ``level_jets`` hold one entry per level, L0 included (L0 bodies/Jets are not + supported). ``refinements=[]`` is one blocked level: periodic SRT root without LES, Jets or exchange observation; its ``level_objects[0]`` are stationary passive circles (static BFL). ``coarse_boundary`` is ``periodic`` (default) or ``slot_reservoir``: a bottom on-site slot and top fixed-density reservoir from ``cfg``. - ``channel`` requires all-blocked levels and matches uniform SRT Zou-He + ``channel`` matches uniform SRT Zou-He inlet, pressure-NEQ outlet and stationary bounce-back y walls. Boundary ghosts evolve internally but public downloads exclude them; initialization must supply their finite positive-density state. Transfers must avoid the @@ -132,14 +131,13 @@ class FixedPatchStepper: stream. On a device/launch failure the instance becomes unusable: close and recreate it, rather than pretending a partially advanced buffer rolled back. - A ``level_blocks`` entry (by, bx) stores that level in blocked.cu storage: - the ENTIRE level shape, including the four-cell patch ghost band, is - partitioned by a non-periodic BlockLayout2D with one-cell same-level halos. - Patch kinds stay per cell (OUTER invalid, INNER stream only, ACTIVE - collide); halos are gathered before each hop. Public arrays keep the dense - shapes; only the allocation differs. A blocked level requires blocked - refined levels below it; several L1 rectangles and depth >= 3 require all - levels blocked. Blocked bodies are stationary circles, optionally + Every level is stored in blocked.cu storage with its ``level_blocks`` tile + (by, bx): the ENTIRE level shape, including the four-cell patch ghost band, + is partitioned by a BlockLayout2D with one-cell same-level halos; the tile + must divide the level (a whole-level tile always does). Patch kinds stay + per cell (OUTER invalid, INNER stream only, ACTIVE collide); halos are + gathered before each hop. Public arrays are dense (9, ny, nx); only the + allocation is blocked. Bodies are stationary circles, optionally surface-actuated: prepare_fine_boundary decides legality on the dense grid, then cut links are staged receiver-local per hop: halo -> BFL -> step, same buffer. A blocked parent gathers halos before receiver-local coalescence @@ -148,8 +146,7 @@ class FixedPatchStepper: ``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 levels are not supported by this diagnostic. + is boundary momentum reaction, not a validated nozzle/stress load. ``observe_window_exchange=True`` (requires observe_raw_exchange) adds a second per-link FP64 ledger accumulated by the same kernel between start_exchange_window() and read_exchange_window(); disabled by default. @@ -180,30 +177,32 @@ class FixedPatchStepper: level_blocks=None, level_objects=None, level_jets=None, observe_raw_exchange=False, observe_window_exchange=False, explosion="constant"): if HALO != 4: - raise RuntimeError("refinement.cu requires a four-cell fine halo") + raise RuntimeError("the refinement transfers require a four-cell fine halo") if coarse_boundary not in ("periodic", "slot_reservoir", "channel"): raise ValueError("coarse_boundary must be periodic, slot_reservoir or channel") depth = len(refinements) + 1 - level_blocks = [None] * depth if level_blocks is None else list(level_blocks) + if level_blocks is None: + raise ValueError("level_blocks is required: one tile (by, bx) per level; a whole-level " + "tile ((ny, nx) at the root, the fine window shape for one rectangle) " + "always divides its level") + level_blocks = list(level_blocks) level_objects = [()] * depth if level_objects is None else list(level_objects) level_jets = [()] * depth if level_jets is None else list(level_jets) level_objects = [tuple(o) for o in level_objects] level_jets = [tuple(j) for j in level_jets] - level_blocks = [None if b is None else tuple(b) for b in level_blocks] + level_blocks = [tuple(b) if isinstance(b, (list, tuple)) else b for b in level_blocks] _check_levels(refinements, level_blocks, level_objects, level_jets) depth = len(refinements) + 1 if explosion not in ("constant", "projected_linear"): raise ValueError("explosion must be constant or projected_linear") if explosion == "projected_linear" and ( - coarse_boundary != "periodic" or any(b is None for b in level_blocks) + coarse_boundary != "periodic" or any(level_objects) or any(level_jets) or cfg.collision != "SRT" or cfg.les_enabled): - raise ValueError("projected explosion requires periodic all-blocked SRT levels " + raise ValueError("projected explosion requires periodic SRT levels " "without bodies/Jets or LES") self._explosion = explosion self._ddf_shifting = bool(cfg.ddf_shifting) - if self._ddf_shifting and any(b is None for b in level_blocks): - raise ValueError("shifted populations require all-blocked levels; dense/mixed is unsupported") if depth == 1 and (coarse_boundary != "periodic" or cfg.collision != "SRT" or cfg.les_enabled or explosion != "constant" or observe_raw_exchange or observe_window_exchange): @@ -225,18 +224,16 @@ class FixedPatchStepper: raise ValueError("slot_reservoir requires y0 >= 4 and y1 <= ny-4") self._layouts = [] for level, block in enumerate(level_blocks): - layout = None - if block is not None: - if level == 0: - layout = BlockLayout2D.build(cfg.ny, cfg.nx, *block, - periodic=(True, not self._wall_aware)) - elif isinstance(pairs[level - 1], FixedPatchSet2D): - layout = pairs[level - 1].fine_layout(block) - else: - layout = BlockLayout2D.build(*pairs[level - 1].fine_shape, *block, - periodic=(False, False)) - if layout.nblocks > MAX_GRID_DIM or layout.local_shape[0] > MAX_GRID_DIM: - raise ValueError(f"{_level_name(level)} blocks exceed launch limits") + if level == 0: + layout = BlockLayout2D.build(cfg.ny, cfg.nx, *block, + periodic=(True, not self._wall_aware)) + elif isinstance(pairs[level - 1], FixedPatchSet2D): + layout = pairs[level - 1].fine_layout(block) + else: + layout = BlockLayout2D.build(*pairs[level - 1].fine_shape, *block, + periodic=(False, False)) + if layout.nblocks > MAX_GRID_DIM or layout.local_shape[0] > MAX_GRID_DIM: + raise ValueError(f"{_level_name(level)} blocks exceed launch limits") self._layouts.append(layout) if self._channel: base = blocked_base_tables(pairs[0], self._layouts[0], np.argsort(pairs[0].coalesce_coarse)) @@ -274,9 +271,6 @@ class FixedPatchStepper: 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 if not isinstance(observe_window_exchange, bool): raise ValueError("observe_window_exchange must be bool") @@ -319,8 +313,6 @@ class FixedPatchStepper: self._allocations = [] self._levels, self._pairs = [], [] self._stream = None - self._fine_flags = self._fine_curved = self._curved = None - self._link_tables = {} self._block = (cfg.threads_per_block, 1, 1) with self.ctx.activate(): try: @@ -348,23 +340,6 @@ class FixedPatchStepper: b.window_integral = self._allocate(count * 4 * 8) b.accumulate = level.module.get_function("BlockAccumulateExchange") b.kernel = level.module.get_function("BlockCurvedBoundary") - if depth == 1: - pass - elif self._level_boundaries[1] is not None: - b = self._level_boundaries[1] - self._fine_flags, self._link_tables, self._curved = b.flags, b.device, b.kernel - elif self._boundary_hosts[1] is not None: - from .curved_links import CurvedLinkSoA - boundary = self._boundary_hosts[1] - 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, - jets=self._level_jets[1], - hits=boundary.curved_hits) - self._fine_curved.upload() - self._curved = self._levels[1].module.get_function("CurvedBoundaryKernel") except BaseException: try: self.close() @@ -380,11 +355,8 @@ class FixedPatchStepper: """ total = 0 for index, layout in enumerate(self._layouts): - cells = int(np.prod(layout.blocked_shape if layout is not None - else (9, *self.hierarchy.kinds[index].shape))) - total += 2 * 4 * cells - if layout is not None: - total += layout.halo_dst.nbytes + layout.halo_src.nbytes + total += 2 * 4 * int(np.prod(layout.blocked_shape)) + total += layout.halo_dst.nbytes + layout.halo_src.nbytes b = self._level_boundaries[index] if b is not None: total += sum(np.asarray(v).nbytes for k, v in b.tables.items() @@ -397,83 +369,62 @@ class FixedPatchStepper: pairs = self.hierarchy.pairs tau = (pairs[0].tau_coarse if pairs else 3 * config.viscosity + .5) if index == 0 \ else pairs[index - 1].tau_fine - name, source = ("PatchCoarseStep" if index == 0 else "PatchFineStep"), "refinement.cu" - if layout is not None: - config = level_config(config, layout, tau) - name, source = ("BlockChannelStep" if index == 0 and self._channel else "BlockStep"), "blocked.cu" - else: - config.validate() - if not np.isclose(config.omega, 1 / tau, rtol=1e-14, atol=0): - raise RuntimeError("level omega/tau mismatch") - host = self._boundary_hosts[index] - n_objects = len(host.objects) if name == "PatchFineStep" and host is not None else 0 - block_root_bc = 0 - if index == 0 and layout is not None: - block_root_bc = 2 if self._channel else 1 if self._wall_aware else 0 - module = compiler.build_module(config, source, arch=self.ctx.sm_arch, n_objects=n_objects, + config = level_config(config, layout, tau) + name = "BlockChannelStep" if index == 0 and self._channel else "BlockStep" + block_root_bc = (2 if self._channel else 1 if self._wall_aware else 0) if index == 0 else 0 + module = compiler.build_module(config, "blocked.cu", arch=self.ctx.sm_arch, block_root_bc=block_root_bc) - compiler.set_params(module, config.omega, n_objects=n_objects) + compiler.set_params(module, config.omega) shape = (9, *self.hierarchy.kinds[index].shape) - size = int(np.prod(shape if layout is None else layout.blocked_shape)) * 4 - grid = ((config.nx + self._block[0] - 1) // self._block[0], config.ny, - 1 if layout is None else layout.nblocks) + size = int(np.prod(layout.blocked_shape)) * 4 + grid = ((config.nx + self._block[0] - 1) // self._block[0], config.ny, layout.nblocks) level = _Level(shape, module, self._allocate(size), self._allocate(size), module.get_function(name), grid, layout) level.boundary = self._level_boundaries[index] - if layout is not None: - level.halo_copy = module.get_function("BlockHaloCopy") - if len(layout.halo_dst): - level.halo_dst = self._upload(layout.halo_dst) - level.halo_src = self._upload(layout.halo_src) - if index == 0 and self._channel: - cfg = self._remap_cfg - ptr, size = module.get_global("d_block_channel") - if size != struct.calcsize("=IIf"): - raise RuntimeError("Block channel host/device layout mismatch") - cuda.memcpy_htod(ptr, struct.pack("=IIf", cfg.nx, cfg.ny, cfg.velocity)) - if index == 0 and self._wall_aware: - cfg = self._remap_cfg - ptr, size = module.get_global("d_block_wall") - if size != struct.calcsize("=IIIf"): - raise RuntimeError("Block wall host/device layout mismatch") - cuda.memcpy_htod(ptr, struct.pack("=IIIf", cfg.nx, cfg.ny, - cfg.wall_slot_width, cfg.wall_normal_velocity)) + level.halo_copy = module.get_function("BlockHaloCopy") + if len(layout.halo_dst): + level.halo_dst = self._upload(layout.halo_dst) + level.halo_src = self._upload(layout.halo_src) + if index == 0 and self._channel: + cfg = self._remap_cfg + ptr, size = module.get_global("d_block_channel") + if size != struct.calcsize("=IIf"): + raise RuntimeError("Block channel host/device layout mismatch") + cuda.memcpy_htod(ptr, struct.pack("=IIf", cfg.nx, cfg.ny, cfg.velocity)) + if index == 0 and self._wall_aware: + cfg = self._remap_cfg + ptr, size = module.get_global("d_block_wall") + if size != struct.calcsize("=IIIf"): + raise RuntimeError("Block wall host/device layout mismatch") + cuda.memcpy_htod(ptr, struct.pack("=IIIf", cfg.nx, cfg.ny, + cfg.wall_slot_width, cfg.wall_normal_velocity)) return level def _build_pair(self, index, pair): """Upload the transfer tables of levels (index, index+1) and the child's block kinds.""" parent, child = self._levels[index], self._levels[index + 1] order = np.argsort(pair.coalesce_coarse) - fine = None if child.layout is None else blocked_fine_tables(pair, child.layout, order) - base = None - if parent.layout is not None: - if index == 0: - base = blocked_base_tables(pair, parent.layout, order, self._wall_aware) - else: - # Nested pairs keep inherited INNER/OUTER and sparse holes of the parent. - base = blocked_base_tables(pair, parent.layout, order, - composed_kind=self.hierarchy.kinds[index]) - ghosts = pair.ghost_fine if fine is None else fine[0] - donors = pair.ghost_coarse if base is None else base[0] - if fine is not None: - donors = donors.astype(np.uint64) - targets = pair.coalesce_coarse[order] if base is None else base[1] - children = pair.coalesce_fine[order] if fine is None else fine[1] - result = _Pair(child.module.get_function("PatchExplode" if fine is None else "BlockExplode"), - None if parent.layout is None else parent.module.get_function("BlockCoalesce"), - self._upload(ghosts), self._upload(donors), self._upload(targets), - self._upload(children), len(ghosts), len(targets), - int(np.prod(pair.coarse_shape if parent.layout is None - else parent.layout.local_shape))) + ghosts, children, child_kind = blocked_fine_tables(pair, child.layout, order) + if index == 0: + donors, targets, parent_kind = blocked_base_tables(pair, parent.layout, order, + self._wall_aware) + else: + # Nested pairs keep inherited INNER/OUTER and sparse holes of the parent. + donors, targets, parent_kind = blocked_base_tables( + pair, parent.layout, order, composed_kind=self.hierarchy.kinds[index]) + result = _Pair(child.module.get_function("BlockExplode"), + parent.module.get_function("BlockCoalesce"), + self._upload(ghosts), self._upload(donors.astype(np.uint64)), + self._upload(targets), self._upload(children), len(ghosts), len(targets), + int(np.prod(parent.layout.local_shape))) if self._explosion == "projected_linear": neighbors, offsets = projected_explosion_tables( pair, parent.layout, self.hierarchy.kinds[index]) result.explode = child.module.get_function("BlockExplodeProjected") result.neighbors, result.offsets = self._upload(neighbors), self._upload(offsets) - if child.layout is not None: - child.kind = self._upload(blocked_kind(child.layout, self.hierarchy.kinds[index + 1])) - if index == 0 and base is not None: - parent.kind = self._upload(base[2]) + child.kind = self._upload(blocked_kind(child.layout, self.hierarchy.kinds[index + 1])) + if index == 0: + parent.kind = self._upload(parent_kind) return result def _allocate(self, size): @@ -506,8 +457,6 @@ class FixedPatchStepper: def _check_storage(self): if any(bool(c.ddf_shifting) != self._ddf_shifting for c in self._level_cfgs): raise ValueError("mixed population representations are unsupported") - if self._ddf_shifting and any(layout is None for layout in self._layouts): - raise ValueError("shifted populations require all-blocked levels") def initialize(self, *level_g): """Upload physical postcollision G, encoding before any FP32 rounding. @@ -552,16 +501,13 @@ class FixedPatchStepper: or (values.sum(axis=0, dtype=np.float64) + int(self._ddf_shifting) <= 0).any()): raise ValueError("Leaf populations must be finite with positive density") self._mask_nonstate(*arrays) - arrays = [a if level.layout is None else level.layout.to_blocked(a) - for a, level in zip(arrays, self._levels)] + arrays = [level.layout.to_blocked(a) for a, level in zip(arrays, self._levels)] try: for level, array in zip(self._levels, arrays): cuda.memcpy_htod(level.current, array) # A single level starts its spare NaN so stale halos cannot hide. cuda.memcpy_htod(level.spare, np.full_like(array, np.nan) if len(self._levels) == 1 else 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._exchange_valid = False self._window_valid = False @@ -660,26 +606,6 @@ class FixedPatchStepper: raise return candidate - 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._boundary_hosts[1] 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 = bool(self._level_jets[1]) or any( - obj.surface_actuation is not None for obj in self._boundary_hosts[1].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 set_jet_amplitude(self, index, amplitude, *, level=1): """Replace one jet's amplitude between steps; no geometry/table rebuild. @@ -703,9 +629,7 @@ class FixedPatchStepper: driven = jet_id >= 0 velocity[driven] = jet_link_velocity(jets, shape, jet_id)[driven] try: - if level == 1 and self._fine_curved is not None: - self._fine_curved.set_jet_amplitude(index, amplitude, self._stream) - elif velocity is not None: + if velocity is not None: self._stream.synchronize() # no in-flight step reads the column cuda.memcpy_htod(boundary.device["surface_velocity"], velocity) except BaseException: @@ -747,19 +671,6 @@ class FixedPatchStepper: nested = index + 1 < len(self._levels) if nested: self._advance(index + 1, src, **kw) - if level.layout is None: - if index == 0: # dense root: PatchCoarseStep fuses restriction and collision - pair = self._pairs[0] - level.kernel(src, dst, self._levels[1].current, pair.targets, pair.children, - np.uint32(pair.n_rows), *(np.uint32(v) for v in self.hierarchy.pairs[0].bounds), - np.uint32(self._wall_aware), grid=level.grid, **kw) - return - if self._boundary_hosts[index] is not None and self._fine_flags is None: - raise RuntimeError("Fine boundary flags are missing") - self._apply_fine_boundary(src, **kw) - flags = self._fine_flags if self._fine_flags is not None else np.uintp(0) - level.kernel(src, dst, flags, grid=level.grid, **kw) - return # Child donors were read above; no halo gather may overwrite staged arrivals below. self._gather_halo(level, src, **kw) if nested: @@ -924,7 +835,7 @@ class FixedPatchStepper: with self._lock, self.ctx.activate(): self._require("ready") self._check_storage() - arrays = [np.empty(level.shape if level.layout is None else level.layout.blocked_shape, + arrays = [np.empty(level.layout.blocked_shape, np.float32) for level in self._levels] try: for level, array in zip(self._levels, arrays): @@ -932,7 +843,7 @@ class FixedPatchStepper: except BaseException: self._status = "failed" raise - arrays = [a if level.layout is None else level.layout.to_dense(a) + arrays = [level.layout.to_dense(a) for a, level in zip(arrays, self._levels)] self._mask_nonstate(*arrays) return _StoredPopulations(tuple(arrays), self._ddf_shifting) @@ -958,14 +869,6 @@ class FixedPatchStepper: b.window_integral = None self._levels.clear() self._pairs = [] - 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 - self._link_tables = {} for pointer in self._allocations: try: pointer.free()