Compare commits

...
Author SHA1 Message Date
ClaudeandClaude Opus 5.5 decd419e94 docs: cite only the recorded MRT parity pass (K1 7813334)
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-21 07:08:31 +00:00
ClaudeandClaude Opus 5.5 1c47cb8570 docs: supported settings review fixes (exact parity test IDs, divisibility guard)
Names the exact SRT/MRT parity tests and their sub-case scope (SRT: one
of 18 sub-cases, unidentified; MRT: 12 sub-cases, default outlet settings,
with the runs that recorded the pass), cites blocks.py:134 for tile
divisibility, labels refinement.py:457 as the fine-boundary re-check, and
says 'measured at 36ffc8b1' for the Taylor-Green observation.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-21 07:05:34 +00:00
ClaudeandClaude Opus 5.5 8e9fc72e19 docs: curved-wall observation states only the recorded matrix outcome
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-21 06:59:14 +00:00
ClaudeandClaude Opus 5.5 d2cd1114ed docs: supported settings for the 2D Jet + multilevel candidate
One page listing supported settings, each with the guard enforcing it
(file:line at c95a7ba) and the existing tests exercising it. Measured
observations are labelled by their exact tested case (SRT strict-parity
1/4320 values ~2 ULP with unresolved cause; V3 TG constant explosion
aggregate; curved-wall r1 pilot; V4 NOT_STEADY), not stated as general
bounds. Documentation only.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-21 06:59:02 +00:00
ClaudeandClaude Opus 5.5 c95a7baae4 test_projected_pairs: regenerate pinned build sources from NUMERICAL_REVISION kernels
The pinned BUILD_SOURCES are the a135252 qualification build inputs.
generate_config hashes the whole kernel snapshot, so after K1/K2 the
working-tree kernels (legitimately) no longer reproduce them. The test
now archives src/CelerisLab/lbm/kernels at NUMERICAL_REVISION and points
compiler_v2._KERNEL_DIR there; config headers are still rendered by
today's compiler. BUILD_SOURCES and the equality assertion are unchanged;
each size is a subTest; a missing NUMERICAL_REVISION fails (no skip).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-21 06:26:45 +00:00
ClaudeandClaude Opus 5.5 f4372a84c0 Merge a454f023 (claude/jet-shapes-once) for integration review
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-21 06:03:36 +00:00
ClaudeandClaude Opus 5.5 e48498d089 Merge 580fcf65 (grok/p1-docs-wording) for integration review
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-21 06:03:36 +00:00
ClaudeandClaude Opus 5.5 531739454e Merge 1cbdc479 (claude/step-single-sync) for integration review
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-21 06:03:36 +00:00
ClaudeandClaude Opus 5.5 b93de2ce24 Merge 1fa498c7 (claude/k2-block-coords) for integration review
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-21 06:03:36 +00:00
ClaudeandClaude Opus 5.5 a454f02355 Jets: compute jet link shapes once per boundary preparation
surface_link_velocity already runs jet_link_shapes; CurvedLinkSoA.
assign_surface and _prepare_blocked_boundary then ran it again with the
same jets and hits to fill the host jet cache. A private _surface_links
returns (velocity, aperture, shape, jet_id) from the one call; both
callers use it. surface_link_velocity keeps its two-value public
contract (now a thin wrapper). Same values, same validation order, same
cache contents; host geometry preparation only.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-21 05:57:11 +00:00
Grok 580fcf6525 docs(P1): align wording with level-list API only
Drop legacy constructor prose (child=/mutual exclusion), retarget
README and fixed_three_level_multipatch examples to refinements /
level_blocks, and fix stepper error/doc strings that still named
fine_block/base_block or a fixed gap >= 4 (unions use dependency
closure). Wording only; no behavior change.
2026-10-21 11:07:29 +08:00
ClaudeandClaude Opus 5.5 1cbdc479cf test_block_step_sync: locate block_stepper without importing it
The module-level import of CelerisLab.lbm.block_stepper loaded
pycuda.driver into the shared test process, unlike the other CPU host
tests (which load mocked copies via AST). find_spec gives the source path
without the import. Test content and assertions unchanged.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-21 01:39:52 +00:00
ClaudeandClaude Opus 5.5 1fa498c7d7 kernels: shared block-window to global coordinate helper
BlockStep (slot-reservoir wall rows) and BlockChannelStep computed the
global cell of a block-window cell with the same four lines.
block_global_xy (core/layout.cuh, 2D, __forceinline__) now holds them;
both kernels call it with their own domain width. Same integer
operations: no numerical, ABI or scheduling change. The slot-reservoir
kernel layout is unchanged (no split).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-21 01:24:27 +00:00
ClaudeandClaude Opus 5.5 781333407e kernels: one shared unforced MRT collision for blocked, channel and patch paths
BlockStep, BlockChannelStep and refinement.cu patch_collide each carried
the same compute_rho_u -> omega clamp -> collide_mrt_no_force sequence.
collide_mrt_unforced (collision_mrt.cuh, __forceinline__) now holds it
once; the three sites call it. Same operations in the same order: no
numerical, ABI or SRT change. The host source contract now requires one
helper definition and exactly one call per file, and forbids the old
inline sequence.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-21 01:15:47 +00:00
15 changed files with 220 additions and 89 deletions

No files matched your search

+27 -26
View File
@@ -125,8 +125,8 @@ with CudaContext(0) as ctx, FixedPatchStepper(cfg, context=ctx, refinements=[(8,
Scope and limits: Scope and limits:
- The examples below use D2Q9, SRT, FP32, double buffer, unshifted, no LES or - The examples below use D2Q9, SRT, FP32, double buffer, unshifted, no LES or
forcing. Fixed 2:1 hierarchies can use level lists of arbitrary depth; the forcing. Fixed 2:1 hierarchies use level lists of arbitrary depth; two- and
legacy two/three-level forms are shown below. three-level examples below use that same constructor.
- `coarse_boundary="periodic"` (default) or `"slot_reservoir"` (bottom slot, top - `coarse_boundary="periodic"` (default) or `"slot_reservoir"` (bottom slot, top
reservoir; patch needs `y0 >= 4`, `y1 <= ny-4`). Periodic mode ignores the wall/slot config fields. reservoir; patch needs `y0 >= 4`, `y1 <= ny-4`). Periodic mode ignores the wall/slot config fields.
- Object centers/radii use the coarse **cell-edge** frame (cell `I` spans `[I, I+1)`), - Object centers/radii use the coarse **cell-edge** frame (cell `I` spans `[I, I+1)`),
@@ -151,7 +151,7 @@ Scope and limits:
inside the union's active region with inside the union's active region with
its boundary stencil clear of transfer cells; a jet may not select links of two its boundary stencil clear of transfer cells; a jet may not select links of two
different bodies, even across regions. Wall-touching patches remain unsupported. different bodies, even across regions. Wall-touching patches remain unsupported.
For nesting, use the level lists or legacy child forms below. For nesting, add further `refinements` entries as shown below.
```python ```python
cfg = LBMConfig(nx=32, ny=16, viscosity=.1) cfg = LBMConfig(nx=32, ny=16, viscosity=.1)
@@ -165,19 +165,20 @@ Scope and limits:
coarse_g, fine_g = run.download() # coarse: NaN on every covered cell; coarse_g, fine_g = run.download() # coarse: NaN on every covered cell;
# fine: NaN on every non-ACTIVE cell # fine: NaN on every non-ACTIVE cell
``` ```
- Legacy three-level scalar chain: add `child=(x0, y0, x1, y1)` in - Three-level nested windows: pass two scalar rectangles in `refinements` —
**L1 window cells** (the patch's fine grid including its 4-cell ghost band) and L1 in L0 cells, L2 in **L1 window cells** (the parent's fine grid including its
all three `level_blocks` entries. L2 must sit inside L1's 4-cell ghost band) — plus all three `level_blocks` entries. L2 must sit inside
ACTIVE cells, away from L1's ghost band and transfer footprint (checked before L1's ACTIVE cells, away from L1's ghost band and transfer footprint (checked
any GPU allocation). `initialize` and `download` take/return three arrays; each before any GPU allocation). `initialize` and `download` take/return three
level is NaN on its non-ACTIVE (and child-covered) cells. One coarse step = 2 L1 arrays; each level is NaN on its non-ACTIVE (and child-covered) cells. One
steps = 4 L2 steps. One L1 patch only. Optional stationary circles and Jets coarse step = 2 L1 steps = 4 L2 steps. One L1 window only in this scalar
live on L2 via `level_objects[2]` and `level_jets[2]`; L1 bodies/Jets retain chain. Optional stationary circles and Jets live on L2 via `level_objects[2]`
their L1 meaning and must be empty in this mode. All object and Jet geometry and `level_jets[2]`; L1 bodies/Jets retain their L1 meaning and must be empty
uses the **L0 edge frame**, not the child bounds' L1 window frame. L2 has in this mode. All object and Jet geometry uses the **L0 edge frame**, not the
dx=dt=0.25, so prescribed lattice velocity is unchanged. Each whole body and child bounds' L1 window frame. L2 has dx=dt=0.25, so prescribed lattice
its BFL stencil must be inside L2 ACTIVE cells, outside the transfer footprint. velocity is unchanged. Each whole body and its BFL stencil must be inside L2
A Jet may select one physical body across multiple blocks, never two bodies. ACTIVE cells, outside the transfer footprint. A Jet may select one physical
body across multiple blocks, never two bodies.
```python ```python
from CelerisLab.body.geometry.circle import CircleGeometry from CelerisLab.body.geometry.circle import CircleGeometry
@@ -199,16 +200,16 @@ Scope and limits:
run.step(4) run.step(4)
coarse_g, l1_g, l2_g = run.download() # solids/covered/nonleaves are NaN coarse_g, l1_g, l2_g = run.download() # solids/covered/nonleaves are NaN
``` ```
- Level lists: `refinements[l]` describes level `l+1` in allocated level-`l` - Level lists (the only constructor): `refinements[l]` describes level `l+1`
integer cells; use a rectangle list for a union, or a scalar rectangle for a in allocated level-`l` integer cells; use a rectangle list for a union, or a
single window (a scalar child requires a scalar parent). `level_blocks`, scalar rectangle for a single window (a scalar child requires a scalar
`level_objects` and `level_jets` each have one entry per level, including L0. parent). `level_blocks`, `level_objects` and `level_jets` each have one entry
These lists are mutually exclusive with the legacy patch keywords. Three or per level, including L0. Three or more levels require all-blocked storage and
more levels require all-blocked storage and bodies/Jets only on the finest bodies/Jets only on the finest level. Every adjacent pair must preserve
level. Every adjacent pair must preserve inherited ACTIVE coverage, valid inherited ACTIVE coverage, valid donors/receivers, surviving pulls and the
donors/receivers, surviving pulls and the upstream transfer footprint. upstream transfer footprint. Level `l` takes `2**l` steps per base interval.
Level `l` takes `2**l` steps per base interval. This describes implemented This describes implemented topology and scheduling, not accuracy or
topology and scheduling, not accuracy or convergence qualification. convergence qualification.
- Blocked body levels can enable `observe_raw_exchange=True` and call - Blocked body levels can enable `observe_raw_exchange=True` and call
`read_raw_exchange(body_id, level=...)` after a completed base interval. `read_raw_exchange(body_id, level=...)` after a completed base interval.
It reduces unique cut-link observations into wall/aperture/total force, It reduces unique cut-link observations into wall/aperture/total force,
+11 -11
View File
@@ -1,8 +1,8 @@
# Fixed three-level rectangle unions # Fixed three-level rectangle unions
This page describes the legacy three-level entry point of `FixedPatchStepper`, This page describes three-level `FixedPatchStepper` setups built from the
with rectangle unions on either adjacent pair. The `refinements` level-list `refinements` level list, with rectangle unions on either adjacent pair. The
entry point also supports arbitrary depth; three or more levels require same constructor supports arbitrary depth; three or more levels require
all-blocked storage and bodies/Jets only on the finest level. This is not all-blocked storage and bodies/Jets only on the finest level. This is not
general AMR: there is no per-region timestep, moving geometry, automatic general AMR: there is no per-region timestep, moving geometry, automatic
dynamic refinement or L1 body support in the three-level form. dynamic refinement or L1 body support in the three-level form.
@@ -12,15 +12,15 @@ No CUDA kernel changes are involved.
## Explicit coordinate frames ## Explicit coordinate frames
Keep the existing scalar `child=(x0,y0,x1,y1)` for a single scalar parent patch. For a single scalar parent patch, pass a scalar rectangle as `refinements[1]`:
Its coordinates retain their original meaning: half-open integer indices in the `(x0,y0,x1,y1)` in half-open integer indices of the allocated L1 window,
allocated L1 window, including its ghost band. including its ghost band.
A rectangle list as the second `refinements` entry (formerly the keyword A rectangle list as `refinements[1]` accepts a scalar or a set parent and
`child_rectangles_in_parent=[...]`) accepts a scalar or a set parent and interprets every rectangle in interprets every rectangle in the actual allocated L1 array, including ghost
the actual allocated L1 array, including ghost bands and holes. It never infers bands and holes. It never infers region-local coordinates. The physical root
region-local coordinates. The physical root cell edge corresponding to an L1 cell edge corresponding to an L1 index `(x,y)` is
index `(x,y)` is `parent.fine_origin + (x,y)/2`: `parent.fine_origin + (x,y)/2`:
- A scalar L1 parent has a shifted window origin. - A scalar L1 parent has a shifted window origin.
- A set L1 parent has origin `(0,0)` and public shape `(2*ny,2*nx)`. - A set L1 parent has origin `(0,0)` and public shape `(2*ny,2*nx)`.
+87
View File
@@ -0,0 +1,87 @@
# Supported settings: 2D Jet + fixed multilevel refinement (candidate c95a7ba)
Every entry names the guard that enforces it (file:line at c95a7ba) and, where one
exists, an existing test that exercises it. A setting not listed here is not
claimed. Observations at the end are tied to the exact case in which they were
measured; they are not general bounds.
## Lattice, storage, collision
| Setting | Supported | Guard |
|---|---|---|
| Dimension / lattice | 2D D2Q9 only (3D refinement is future work) | `lbm/refinement.py:65`, `lbm/kernels/blocks.cu:25`, `lbm/kernels/refinement.cu:26` |
| Precision / streaming | FP32 storage, double buffer | same guards |
| Collision on refinement levels | SRT or MRT; TRT and LES rejected | `lbm/refinement.py:65` (re-checked for fine boundaries at `lbm/refinement.py:457`) |
| Shifted (ddf_shifting) storage | only when every level is blocked; dense or mixed rejected | `lbm/refinement_stepper.py:198`, `:469`; `lbm/kernels/refinement.cu:26` (dense patch kernel unshifted only) |
| Fine relaxation | 1/tau_fine must lie inside the configured omega guards | `lbm/refinement.py:80` |
| Uniform blocked backend `PeriodicBlockStepper` | doubly periodic, SRT only, stationary passive circles | `lbm/block_stepper.py:84`, `:91`, `:45`-`:67` |
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`.
## Refinement hierarchy (`FixedPatchStepper(cfg, *, context, refinements, ...)`)
| Setting | Supported | Guard |
|---|---|---|
| Depth | at least one refinement level; one entry per level in `level_blocks`/`level_objects`/`level_jets` | `lbm/refinement_stepper.py:72`, `:75` |
| Single L1 rectangle | two-coarse-cell margin to the domain edge | `lbm/refinement.py:73` |
| Several L1 rectangles (union) | require `level_blocks[0]` and `level_blocks[1]`; transfers must close (dependency checks, no fixed gap rule) | `lbm/refinement_stepper.py:84`; `lbm/refinement.py:283`-`:295` |
| Nested levels | depth >= 3 must be all blocked; a blocked level requires blocked levels below; a scalar child needs a scalar parent | `lbm/refinement_stepper.py:82`, `:87`; `lbm/hierarchy.py:54` |
| Nested transfer support | child only over inherited ACTIVE cells, donors/receivers uncovered ACTIVE, no child over the parent restriction footprint | `lbm/refinement.py:391`-`:407` |
| Tile shape | positive `(by, bx)`; the level allocation must be divisible by the tile (no partial tiles yet) | `lbm/refinement_stepper.py:78`, `lbm/refinement.py:246` (shape); `lbm/blocks.py:134` (divisibility) |
| Coarse boundary | `periodic`, `slot_reservoir` or `channel` | `lbm/refinement_stepper.py:177` |
| `slot_reservoir` | on-site velocity bottom slot, reservoir top, velocity 0, 2 <= slot width < nx, patch rows 4..ny-4 | `lbm/refinement_stepper.py:205`, `:212`; `lbm/refinement.py:546` |
| `channel` | all levels blocked; Zou-He inlet, pressure-NEQ outlet, stationary bounce-back walls; transfers away from the closure | `lbm/channel.py:22`, `:26`, `:43`, `:54` |
| Explosion | `constant` (default); `projected_linear` only periodic, all-blocked SRT, no bodies/Jets, no LES (experimental) | `lbm/refinement_stepper.py:188`, `:193` |
| L2 runtime remap (`remap_child`) | periodic three-level L2 sets only, no L0/L1 bodies, identical pinned boundary | `lbm/refinement_stepper.py:557`, `:567`-`:581` |
Tests: `tests/unit/test_refinement_three_level_multipatch.py`, `tests/validation/multigrid/test_refinement_multipatch.py`,
`test_refinement_three_level.py`, `test_refinement_three_level_multipatch_gpu.py`, `test_refinement_blocked_base.py`,
`test_refinement_blocked_fine.py`, `test_channel_root_gpu.py`, `tests/unit/test_channel_root.py`,
`tests/unit/test_projected_pair_stepper.py`, `test_runtime_child_remap_gpu.py`, `test_pinned_body_remap_gpu.py`.
## Bodies and Jets on refined levels
| Setting | Supported | Guard |
|---|---|---|
| Bodies | stationary circles without sensor/force-region roles, entirely inside the patch union, fully resolved, disjoint, stencil away from the two-hop transfer footprint | `lbm/refinement.py:474`-`:542` |
| Level | not on L0; only on the finest level of a nested hierarchy | `lbm/refinement_stepper.py:80`, `:89` |
| Jet condition | velocity only; `pressure` is reserved and raises NotImplementedError | `body/jets.py:148`, `:150` |
| Jet geometry | each Jet selects resolved links of exactly one stationary body; no Jet/Jet or Jet/SurfaceActuation overlap; abs(amplitude) < 1 lattice | `body/jets.py:155`, `:175`, `:177`, `:255`, `:261`, `:263` |
| Jets with shifted storage | only in the all-blocked stepper with SRT and no LES | `body/jets.py:249` |
| SurfaceActuation | unshifted storage only | `body/objects.py:155` |
| Raw boundary exchange (force) | blocked body levels only; the window ledger requires `observe_raw_exchange` | `lbm/refinement_stepper.py:264`, `:269` |
| Checkpoint with Jets | refused | `simulation.py:687` |
Raw exchange is the cut-link momentum reaction (wall + aperture), not a validated
nozzle thrust or stress load. Jet amplitude is a lattice velocity, not a flow rate.
Tests: `tests/unit/test_jets.py`, `test_jets_review.py`, `test_shifted_jets.py`, `test_raw_exchange.py`,
`tests/validation/multigrid/test_refinement_jets.py`, `test_shifted_jets_gpu.py`, `test_raw_exchange_gpu.py`,
`test_exchange_window.py`, `test_fine_surface_gpu.py`.
## Observations from specific tested cases (not general bounds)
- SRT channel root vs uniform OneStep, strict-parity test
`tests/validation/multigrid/test_channel_root_gpu.py::ChannelRootGPU::test_uniform_matched_root`:
in one of its 18 sub-cases (2 storage modes x 3 profile/clamp/damp settings x 3 blocks; which
one is not recorded), 1 of the 4320 compared values differs, max abs 2.33e-10, max rel 1.77e-7
(about 2 ULP). The test remains failed; its cause and any general bound are unresolved.
- MRT root parity test
`tests/validation/multigrid/test_refinement_mrt_gpu.py::RefinementMRTGPU::test_uniform_matched_mrt_root`
(12 sub-cases: 2 storage modes x 2 profiles x 3 blocks, default outlet settings) passed bitwise
in the Horsea focused GPU run at K1 7813334 (2026-10-21).
- Periodic Taylor-Green Re 100, measured at 36ffc8b1 (V3 matched pair, N128 depth 2):
every per-sample screen passed; the constant explosion's aggregate worst velocity L2,
0.004418, exceeds its preregistered 0.002 budget, while projected_linear passed that budget.
- Curved-wall pressure (pilot harness, exact-q analytic annulus, tau 0.8): the coarse swirl r1
case exceeded its terminal density-gain cap (1.105e-4 vs 1e-4), so the frozen 9-case matrix
is NOT_ESTABLISHED. No validated fix exists. This is a pilot harness, not production code.
- DFG 2D-1 Re 20 (V4, SRT, t* 60, pressure-NEQ outlet and local Zou-He inlet): refined and uniform
runs were NOT_STEADY under the frozen steady rule, so no interval verdict was given.
## Not supported in this release
3D, dynamic/automatic AMR (tagging, time-varying hierarchies beyond `remap_child`),
moving bodies on refined levels, pressure Jets, Jet checkpoints, LES or TRT on refined
levels, partial (non-dividing) tiles.
+14 -3
View File
@@ -207,19 +207,30 @@ def surface_link_velocity(objects, cfg, flags, fluid_idx, dir, body_id, rx, ry,
jets=(), hits=None, *, _shifted_blocked=False): jets=(), hits=None, *, _shifted_blocked=False):
"""Per-link prescribed surface velocity (n, 2) FP32 and aperture mask. """Per-link prescribed surface velocity (n, 2) FP32 and aperture mask.
See _surface_links, which also returns the unit-amplitude jet shapes.
"""
return _surface_links(objects, cfg, flags, fluid_idx, dir, body_id, rx, ry, jets, hits,
_shifted_blocked=_shifted_blocked)[:2]
def _surface_links(objects, cfg, flags, fluid_idx, dir, body_id, rx, ry,
jets=(), hits=None, *, _shifted_blocked=False):
"""(velocity, aperture, jet shape, jet id); the last two are None without jets.
Host-only and storage-neutral: rows follow the given physical cut-link Host-only and storage-neutral: rows follow the given physical cut-link
order. Returns (None, None) when neither an object carries surface order. Returns all None when neither an object carries surface
actuation nor a jet is given. Jets need true FP64 ``hits`` (n, 2) and may actuation nor a jet is given. Jets need true FP64 ``hits`` (n, 2) and may
not share a link with each other or with a legacy aperture. not share a link with each other or with a legacy aperture.
Shifted storage is privately enabled only by the all-blocked stepper. Shifted storage is privately enabled only by the all-blocked stepper.
""" """
enabled = [obj for obj in objects if obj.surface_actuation is not None] enabled = [obj for obj in objects if obj.surface_actuation is not None]
if not enabled and not jets: if not enabled and not jets:
return None, None return None, None, None, None
from .preprocess.flag_overlay import validate_cut_link_ownership from .preprocess.flag_overlay import validate_cut_link_ownership
validate_cut_link_ownership(flags, fluid_idx, dir, cfg.nq) validate_cut_link_ownership(flags, fluid_idx, dir, cfg.nq)
surface_velocity = np.zeros((len(dir), 2), np.float32) surface_velocity = np.zeros((len(dir), 2), np.float32)
surface_aperture = np.zeros(len(dir), bool) surface_aperture = np.zeros(len(dir), bool)
shape = jet_id = None
for obj in enabled: for obj in enabled:
obj.validate_surface_actuation(cfg) obj.validate_surface_actuation(cfg)
selected = body_id == obj.obj_id selected = body_id == obj.obj_id
@@ -252,4 +263,4 @@ def surface_link_velocity(objects, cfg, flags, fluid_idx, dir, body_id, rx, ry,
raise ValueError("Jet overlaps a SurfaceActuation aperture") raise ValueError("Jet overlaps a SurfaceActuation aperture")
surface_velocity[driven] = jet_link_velocity(jets, shape, jet_id)[driven] surface_velocity[driven] = jet_link_velocity(jets, shape, jet_id)[driven]
surface_aperture |= driven surface_aperture |= driven
return surface_velocity, surface_aperture return surface_velocity, surface_aperture, shape, jet_id
+2 -4
View File
@@ -94,12 +94,10 @@ class CurvedLinkSoA:
""" """
self.surface_velocity = self.aperture = None self.surface_velocity = self.aperture = None
self.jets, self.jet_shape, self.jet_id = tuple(jets), None, None self.jets, self.jet_shape, self.jet_id = tuple(jets), None, None
self.surface_velocity, self.aperture = surface_link_velocity( from ..body.jets import _surface_links
self.surface_velocity, self.aperture, self.jet_shape, self.jet_id = _surface_links(
objects, cfg, flags, self.fluid_idx, self.dir, self.body_id, self.rx, self.ry, objects, cfg, flags, self.fluid_idx, self.dir, self.body_id, self.rx, self.ry,
jets=self.jets, hits=hits) jets=self.jets, hits=hits)
if self.jets:
from ..body.jets import jet_link_shapes
self.jet_shape, self.jet_id = jet_link_shapes(self.jets, hits)
def set_jet_amplitude(self, index, amplitude, stream=None): def set_jet_amplitude(self, index, amplitude, stream=None):
"""Rewrite jet rows from cached shapes; re-upload only the velocity column. """Rewrite jet rows from cached shapes; re-upload only the velocity column.
+3 -9
View File
@@ -168,21 +168,15 @@ extern "C" __global__ void BlockStep(const fpxx* in, fpxx* out, unsigned int hal
stream_pull_load(k, f, in + base, j); stream_pull_load(k, f, in + base, j);
if (cell_kind == KIND_ACTIVE) { if (cell_kind == KIND_ACTIVE) {
if (d_block_wall.nx != 0) { if (d_block_wall.nx != 0) {
const unsigned int bx = NX - 2 * halo, by = NY - 2 * halo; unsigned int gx, gy;
const unsigned int nbx = d_block_wall.nx / bx; block_global_xy(x, y, halo, d_block_wall.nx, gx, gy);
const unsigned int gx = (blockIdx.z % nbx) * bx + x - halo;
const unsigned int gy = (blockIdx.z / nbx) * by + y - halo;
if (gy == 1) if (gy == 1)
reconstruct_bottom_velocity(f, 0.0f, bottom_wall_normal_velocity( reconstruct_bottom_velocity(f, 0.0f, bottom_wall_normal_velocity(
gx, d_block_wall.nx, d_block_wall.slot_width, d_block_wall.amplitude)); gx, d_block_wall.nx, d_block_wall.slot_width, d_block_wall.amplitude));
else if (gy == d_block_wall.ny - 2) reconstruct_top_reservoir(f, RHO); else if (gy == d_block_wall.ny - 2) reconstruct_top_reservoir(f, RHO);
} }
#if COLLISION_MODEL == 2 #if COLLISION_MODEL == 2
// Uniform MRT sequence (collide_dispatch, zero forcing); SRT path below unchanged. collide_mrt_unforced(f); // SRT path below unchanged.
float rho, ux, uy;
compute_rho_u(f, rho, ux, uy);
collide_mrt_no_force(f, rho, ux, uy,
fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, d_params.omega)));
#else #else
float rho, ux, uy, feq[NQ]; float rho, ux, uy, feq[NQ];
compute_rho_u(f, rho, ux, uy); compute_rho_u(f, rho, ux, uy);
@@ -53,6 +53,18 @@ __device__ __forceinline__ void index_from_thread(unsigned int& x,
k = (unsigned long)y * NX + x; k = (unsigned long)y * NX + x;
} }
// Global cell of window cell (x, y) in block blockIdx.z of a row-major,
// full-domain block layout domain_nx cells wide. NX/NY are the local strides.
__device__ __forceinline__ void block_global_xy(unsigned int x, unsigned int y,
unsigned int halo, unsigned int domain_nx,
unsigned int& gx, unsigned int& gy)
{
const unsigned int bx = NX - 2 * halo, by = NY - 2 * halo;
const unsigned int nbx = domain_nx / bx;
gx = (blockIdx.z % nbx) * bx + x - halo;
gy = (blockIdx.z / nbx) * by + y - halo;
}
// --------------------------------------------------------------------------- // ---------------------------------------------------------------------------
// Neighbor list (periodic BCs via modular arithmetic) // Neighbor list (periodic BCs via modular arithmetic)
// --------------------------------------------------------------------------- // ---------------------------------------------------------------------------
@@ -114,6 +114,16 @@ __device__ __forceinline__ void collide_mrt_no_force(float* __restrict__ g,
collide_mrt(g, rho, ux, uy, Fin, omega); collide_mrt(g, rho, ux, uy, Fin, omega);
} }
// Unforced MRT with the clamped device omega, as in collide_dispatch with zero
// forcing. Shared by the blocked, channel-root and dense-patch collisions.
__device__ __forceinline__ void collide_mrt_unforced(float* __restrict__ g)
{
float rho, ux, uy;
compute_rho_u(g, rho, ux, uy);
collide_mrt_no_force(g, rho, ux, uy,
fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, d_params.omega)));
}
#elif LATTICE_MODEL == LATTICE_D3Q19 #elif LATTICE_MODEL == LATTICE_D3Q19
__device__ __forceinline__ void collide_mrt(float* __restrict__ g, __device__ __forceinline__ void collide_mrt(float* __restrict__ g,
+1 -4
View File
@@ -37,10 +37,7 @@ __device__ unsigned int incoming_counts[NQ * TOTAL_CELLS];
__device__ __forceinline__ void patch_collide(float* f, unsigned long k) { __device__ __forceinline__ void patch_collide(float* f, unsigned long k) {
#if COLLISION_MODEL == 2 #if COLLISION_MODEL == 2
float rho, ux, uy; collide_mrt_unforced(f);
compute_rho_u(f, rho, ux, uy);
collide_mrt_no_force(f, rho, ux, uy,
fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, d_params.omega)));
#else #else
float rho, ux, uy, feq[NQ]; float rho, ux, uy, feq[NQ];
compute_rho_u(f, rho, ux, uy); compute_rho_u(f, rho, ux, uy);
@@ -15,10 +15,8 @@ extern "C" __global__ void BlockChannelStep(
unsigned int x, y; unsigned long k; unsigned int x, y; unsigned long k;
index_from_thread(x, y, k); index_from_thread(x, y, k);
if (x < halo || y < halo || x >= NX - halo || y >= NY - halo) return; if (x < halo || y < halo || x >= NX - halo || y >= NY - halo) return;
const unsigned int bx = NX - 2 * halo, by = NY - 2 * halo; unsigned int gx, gy;
const unsigned int nbx = d_block_channel.nx / bx; block_global_xy(x, y, halo, d_block_channel.nx, gx, gy);
const unsigned int gx = (blockIdx.z % nbx) * bx + x - halo;
const unsigned int gy = (blockIdx.z / nbx) * by + y - halo;
const unsigned long base = (unsigned long)blockIdx.z * NQ * TOTAL_CELLS; const unsigned long base = (unsigned long)blockIdx.z * NQ * TOTAL_CELLS;
const unsigned long cell = (unsigned long)blockIdx.z * TOTAL_CELLS + k; const unsigned long cell = (unsigned long)blockIdx.z * TOTAL_CELLS + k;
float f[NQ]; float f[NQ];
@@ -55,11 +53,7 @@ extern "C" __global__ void BlockChannelStep(
} }
if (!wall && !outlet) { if (!wall && !outlet) {
#if COLLISION_MODEL == 2 #if COLLISION_MODEL == 2
// Match uniform MRT arithmetic (collide_dispatch, zero forcing). collide_mrt_unforced(f);
float rho, ux, uy;
compute_rho_u(f, rho, ux, uy);
collide_mrt_no_force(f, rho, ux, uy,
fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, d_params.omega)));
#else #else
// Match uniform SRT arithmetic, including its shifted equilibrium path. // Match uniform SRT arithmetic, including its shifted equilibrium path.
float rho, ux, uy, feq[NQ], forcing[NQ] = {}; float rho, ux, uy, feq[NQ], forcing[NQ] = {};
+7 -5
View File
@@ -75,13 +75,13 @@ def _check_levels(refinements, blocks, objects, jets):
raise ValueError(f"{name} needs one entry per level ({depth})") raise ValueError(f"{name} needs one entry per level ({depth})")
for level, block in enumerate(blocks): for level, block in enumerate(blocks):
if block is not None and len(block) != 2: if block is not None and len(block) != 2:
raise ValueError(f"{_level_name(level)}_block must be (by, bx)") raise ValueError(f"level_blocks[{level}] must be (by, bx)")
if objects[0] or jets[0]: if objects[0] or jets[0]:
raise ValueError("base-level bodies/Jets are not supported") 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)): 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") 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): if not is_rectangle(refinements[0]) and (blocks[0] is None or blocks[1] is None):
raise ValueError("several patches require fine_block and base_block") raise ValueError("several L1 patches require level_blocks[0] and level_blocks[1]")
if depth > 2: if depth > 2:
if any(block is None for block in blocks): if any(block is None for block in blocks):
raise ValueError("three or more levels must all be blocked") raise ValueError("three or more levels must all be blocked")
@@ -149,9 +149,11 @@ class FixedPatchStepper:
Nested levels: a child set's coordinates index the allocated parent array, Nested levels: a child set's coordinates index the allocated parent array,
including ghost bands and holes: root edge = parent fine_origin + index / 2 including ghost bands and holes: root edge = parent fine_origin + index / 2
(zero origin for a parent set; a scalar parent keeps its shifted window (zero origin for a parent set; a scalar parent keeps its shifted window
origin). A scalar child requires a scalar parent. Each pair requires gap >= 4 origin). A scalar child requires a scalar parent. Rectangle unions impose no
in its own parent-cell units. All regions share one buffer pair per level fixed gap rule: adjacent, overlapping and concave unions are accepted when
and the fixed schedule of 2**l hops per base step at level l. Bodies/Jets are supported only on the their actual transfer dependencies close (FixedPatchSet2D /
nest_patch_set). All regions share one buffer pair per level and the fixed
schedule of 2**l hops per base step at level l. Bodies/Jets are supported only on the
finest level of a nested hierarchy (L2 dx=dt=.25 in root units, unchanged finest level of a nested hierarchy (L2 dx=dt=.25 in root units, unchanged
lattice velocity); every whole circle and its BFL stencil must fit strictly lattice velocity); every whole circle and its BFL stencil must fit strictly
away from that level's transfer footprint. away from that level's transfer footprint.
+6 -4
View File
@@ -8,7 +8,7 @@ from dataclasses import dataclass, replace
import numpy as np import numpy as np
from .blocks import blocked_link_tables from .blocks import blocked_link_tables
from ..body.jets import jet_link_shapes, surface_link_velocity from ..body.jets import _surface_links
from .descriptors import D2Q9_EX, D2Q9_EY from .descriptors import D2Q9_EX, D2Q9_EY
from .refinement import ACTIVE, OUTER_GHOST from .refinement import ACTIVE, OUTER_GHOST
@@ -142,11 +142,13 @@ def _prepare_blocked_boundary(boundary, layout, cfg, jets, *, _shifted_blocked=F
tables.update(body_id=np.asarray(body_id, np.int32), rx=np.asarray(rx, np.float32), 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)) 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) actuated = bool(jets) or any(o.surface_actuation is not None for o in boundary.objects)
jet_cache = None
if actuated: if actuated:
cfg = replace(cfg) cfg = replace(cfg)
cfg.validate() cfg.validate()
tables["surface_velocity"], tables["aperture"] = surface_link_velocity( tables["surface_velocity"], tables["aperture"], shape, jet_id = _surface_links(
boundary.objects, cfg, boundary.flags, fluid_idx, dirs, body_id, rx, ry, boundary.objects, cfg, boundary.flags, fluid_idx, dirs, body_id, rx, ry,
jets=jets, hits=boundary.curved_hits, _shifted_blocked=_shifted_blocked) jets=jets, hits=boundary.curved_hits, _shifted_blocked=_shifted_blocked)
return _BlockedBoundary(boundary, tables, jets, actuated, if jets:
jet_link_shapes(jets, boundary.curved_hits) if jets else None) jet_cache = (shape, jet_id)
return _BlockedBoundary(boundary, tables, jets, actuated, jet_cache)
+4 -2
View File
@@ -1,6 +1,7 @@
"""Host contract of PeriodicBlockStepper.step(n): one stream wait per call (mocked driver).""" """Host contract of PeriodicBlockStepper.step(n): one stream wait per call (mocked driver)."""
import ast import ast
from contextlib import nullcontext from contextlib import nullcontext
import importlib.util
from pathlib import Path from pathlib import Path
import types import types
from unittest.mock import Mock from unittest.mock import Mock
@@ -8,14 +9,15 @@ import unittest
import numpy as np import numpy as np
import CelerisLab.lbm.block_stepper as b
from CelerisLab.config import LBMConfig from CelerisLab.config import LBMConfig
from CelerisLab.lbm.blocks import BlockLayout2D from CelerisLab.lbm.blocks import BlockLayout2D
def host_module(): def host_module():
"""block_stepper with pycuda/compiler replaced by mocks; no GPU context.""" """block_stepper with pycuda/compiler replaced by mocks; no GPU context."""
path = Path(b.__file__) # Locate the source without importing it: importing block_stepper loads pycuda.driver,
# which CPU-only import-boundary tests (test_build_contracts) must never see.
path = Path(importlib.util.find_spec("CelerisLab.lbm.block_stepper").origin)
tree = ast.parse(path.read_text()) tree = ast.parse(path.read_text())
tree.body = [node for node in tree.body if not ( tree.body = [node for node in tree.body if not (
isinstance(node, ast.Import) and any(a.name == "pycuda.driver" for a in node.names) isinstance(node, ast.Import) and any(a.name == "pycuda.driver" for a in node.names)
+6 -1
View File
@@ -37,8 +37,13 @@ class RefinementMRTHost(unittest.TestCase):
self.assertIn("STORE_PRECISION != 0", source) self.assertIn("STORE_PRECISION != 0", source)
self.assertIn('#if COLLISION_MODEL == 2\n#include "operators/collision_mrt.cuh"\n#endif', source) self.assertIn('#if COLLISION_MODEL == 2\n#include "operators/collision_mrt.cuh"\n#endif', source)
self.assertIn("USE_DDF_SHIFTING ||", refinement) # dense/mixed shifted stays rejected self.assertIn("USE_DDF_SHIFTING ||", refinement) # dense/mixed shifted stays rejected
# One shared unforced MRT body; each blocked/patch path calls it exactly once.
mrt = (KERNELS / "operators/collision_mrt.cuh").read_text()
self.assertEqual(mrt.count("void collide_mrt_unforced("), 1)
for source in (blocks, channel, refinement): for source in (blocks, channel, refinement):
self.assertEqual(source.count("collide_mrt_no_force("), 1) self.assertEqual(source.count("collide_mrt_unforced(f);"), 1)
self.assertNotIn("collide_mrt_no_force(", source)
self.assertNotIn("compute_rho_u(f, rho, ux, uy);\n collide_mrt", source)
if __name__ == "__main__": unittest.main() if __name__ == "__main__": unittest.main()
@@ -249,19 +249,35 @@ class ProjectedPairsTests(unittest.TestCase):
self.assertFalse(assess(incomplete)["comparisons"][3]["passed"]) self.assertFalse(assess(incomplete)["comparisons"][3]["passed"])
def test_pinned_build_sources_match_host_generated_inputs(self): def test_pinned_build_sources_match_host_generated_inputs(self):
"""The pinned fingerprints are what today's compiler generates from the qualified
NUMERICAL_REVISION kernels; the current kernels may legitimately differ."""
import tarfile
from CelerisLab.cuda import compiler_v2 from CelerisLab.cuda import compiler_v2
for size in (128, 256): kernels = "src/CelerisLab/lbm/kernels"
cfg = replace(run_case.config(size, 100), nx=10, ny=10) if subprocess.run(["git", "-C", str(REPO), "cat-file", "-e",
cfg.validate() f"{NUMERICAL_REVISION}^{{commit}}"], capture_output=True).returncode:
root = compiler_v2.generate_config(cfg) self.fail(f"{NUMERICAL_REVISION} is not in this clone; the pinned build "
try: "sources cannot be checked")
actual = json.loads((Path(root) / ".celeris-build").read_text()) archive = subprocess.run(["git", "-C", str(REPO), "archive", "--format=tar",
expected = build_metadata(size) NUMERICAL_REVISION, kernels],
self.assertEqual(actual, {key: expected[key] for key in capture_output=True, check=True).stdout
("source_sha256", "config", "n_objects")}) with tempfile.TemporaryDirectory() as snapshot:
finally: with tarfile.open(fileobj=io.BytesIO(archive)) as tar:
compiler_v2.cleanup_build_dir(root) tar.extractall(snapshot)
with patch.object(compiler_v2, "_KERNEL_DIR", str(Path(snapshot) / kernels)):
for size in (128, 256):
with self.subTest(size=size):
cfg = replace(run_case.config(size, 100), nx=10, ny=10)
cfg.validate()
root = compiler_v2.generate_config(cfg)
try:
actual = json.loads((Path(root) / ".celeris-build").read_text())
expected = build_metadata(size)
self.assertEqual(actual, {key: expected[key] for key in
("source_sha256", "config", "n_objects")})
finally:
compiler_v2.cleanup_build_dir(root)
def test_source_preflight_requires_clean_exact_revision_and_tree(self): def test_source_preflight_requires_clean_exact_revision_and_tree(self):
import oracle import oracle