Frank14f/CelerisLab
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

refinement: remove the dense refinement route; level_blocks is required (D2)

Browse Source 
Deleted: kernels/refinement.cu and step/patch_step.cu (PatchExplode,
PatchFineStep, PatchCoarseStep, the CELERIS_REFINEMENT_AUDIT counters);
in FixedPatchStepper the dense branches of _build_level/_build_pair/_hop,
_apply_fine_boundary, the fine-level CurvedLinkSoA wiring (_fine_curved,
_fine_flags, _curved, _link_tables, its set_jet_amplitude/initialize/close
paths) and the now-vacuous 'requires all-blocked' guards (shifted,
projected, channel, raw exchange, several L1 rectangles, depth >= 3).

level_blocks is required: one tile (by, bx) of positive integers per level,
else ValueError before any device work naming the whole-level tile
((ny, nx) at the root, the fine window shape for one rectangle) as always
valid. No automatic tile choice. Removed route: None in level_blocks
(dense or mixed dense/blocked storage); parity of D/D, D/B and B/B is
recorded at 1573942 and 2469af7.

Kept: the uniform-grid solver (uniform.cu, Simulation, CurvedLinkSoA,
CurvedBoundaryKernel), FixedPatch2D, pinned historical runners.
Docs: compiler entries, BlockExplode comment, lifecycle, supported_settings.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Claude
2026-10-27 15:09:13 +00:00
co-authored by Claude Opus 5.5
parent 104f74632f
commit e18f7abe2b
7 changed files with 97 additions and 328 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/lifecycle.md
+4 -3
View File
@@ -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
docs/supported_settings.md
+3 -2
View File
@@ -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`.
src/CelerisLab/cuda/compiler.py
+1 -2
View File
@@ -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
src/CelerisLab/lbm/kernels/refinement.cu
-40
View File
@@ -1,40 +0,0 @@
// Fixed, unforced internal 2:1 patch. Separate compilation entry;
// default OneStep/channel kernels are deliberately not included.
#include <stdint.h>
#include <math.h>
#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"
src/CelerisLab/lbm/kernels/step/block_transfer.cu
+3 -3
View File
@@ -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)
src/CelerisLab/lbm/kernels/step/patch_step.cu
-95
View File
@@ -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);
}
src/CelerisLab/lbm/refinement_stepper.py
+86 -183
View File
@@ -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()