Compare commits

...
Author SHA1 Message Date
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
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
11 changed files with 88 additions and 47 deletions
+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] = {};
+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