From c33caf636d3d8fe02aa5a81631d3cec3c041516b Mon Sep 17 00:00:00 2001 From: Codex Date: Fri, 16 Oct 2026 23:28:29 -0700 Subject: [PATCH] Retain shifted density perturbation in opt-in blocked collision --- docs/shifted_blocked_candidate.md | 43 +++++++++++++++++++ src/CelerisLab/lbm/kernels/blocks.cu | 10 ++++- .../lbm/kernels/boundary/planar_velocity.cuh | 21 +++++++-- .../lbm/kernels/operators/macro.cuh | 17 ++++++-- src/CelerisLab/lbm/refinement.py | 8 ++-- tests/unit/cpp/shifted_block_operators.cpp | 36 ++++++++++++++++ tests/unit/test_shifted_block_operators.py | 32 ++++++++++++++ 7 files changed, 154 insertions(+), 13 deletions(-) create mode 100644 docs/shifted_blocked_candidate.md create mode 100644 tests/unit/cpp/shifted_block_operators.cpp create mode 100644 tests/unit/test_shifted_block_operators.py diff --git a/docs/shifted_blocked_candidate.md b/docs/shifted_blocked_candidate.md new file mode 100644 index 0000000..0b4e919 --- /dev/null +++ b/docs/shifted_blocked_candidate.md @@ -0,0 +1,43 @@ +# Opt-in shifted blocked candidate + +This candidate is based on fbe7377. The default remains unshifted. Enable with +`LBMConfig(ddf_shifting=True)` only on the supported routes below. + +| Route | Candidate support | +|---|---| +| PeriodicBlockStepper, FP32 D2Q9 SRT double-buffer, no LES | shifted or unshifted | +| FixedPatchStepper, every level blocked and same representation | shifted or unshifted | +| Passive stationary circular BFL boundary and raw exchange | supported candidate scope | +| Synchronized three-level child remap | private raw shifted snapshot, exact overlap/copy | +| Dense root/refinement or mixed representation | rejected for shifted | +| Slot/reservoir, surface actuation, velocity Jet | existing public shifted exclusion retained | +| Other collision/lattice/storage/streaming/LES combinations | existing exclusions retained | + +The D2Q9 equilibrium helper accepts a retained density offset only in shifted +compilation. BlockStep computes delta=sum(g) before adding the unit density; +velocity uses the existing macro calculation. The existing SRT collision and +all unshifted expressions remain. No global precision type or duplicate +collision implementation is introduced. Ordinary shifted Simulation callers +still use the existing rho−1 overload; this candidate does not change their +numerical method. + +Public initialize receives physical populations and performs shifted subtraction +in FP64 before the one FP32 storage cast. Public shifted download returns FP64 +physical populations with exact-rational D2Q9 background weights. Internal raw +snapshots retain FP32 g and a representation tag; runtime remap never decodes +and re-encodes overlap. Exact mass is interpreted as 1+sum(g), and momentum as +sum(c*g); analyses should use raw values for tiny perturbation budgets. The +FP32 BFL force constants have their own arithmetic rounding; no exact global +conservation claim follows from shifting. + +Planar helper formulas include their shifted reference-density offsets and have +actual-header host tests. Their public shifted route remains excluded until a +separate test-gated migration. Next velocity-Jet work must validate physical +wall density, incoming/outgoing populations, raw force/mass, ownership and +transfer interactions before relaxing that restriction. + +Evidence so far: selected-cell CUDA replay favored retained-delta shifting; +this is not full-solver validation. Host algebra/storage/remap tests and staged +GPU continuation, three-grid cylinder physics, raw budgets, Couette helper +stress and measured runtime are separate gates. Original failed experiments +remain evidence; no acceptance threshold is relaxed. diff --git a/src/CelerisLab/lbm/kernels/blocks.cu b/src/CelerisLab/lbm/kernels/blocks.cu index 739f1fb..f0704f6 100644 --- a/src/CelerisLab/lbm/kernels/blocks.cu +++ b/src/CelerisLab/lbm/kernels/blocks.cu @@ -18,8 +18,8 @@ #include "boundary/curved_boundary.cuh" #include "boundary/planar_velocity.cuh" -#if DIM != 2 || NQ != 9 || COLLISION_MODEL != 0 || STORE_PRECISION != 0 || USE_DDF_SHIFTING || USE_LES || STREAMING_MODEL != 0 -#error "Blocks require unshifted FP32 D2Q9 SRT double buffers without LES" +#if DIM != 2 || NQ != 9 || COLLISION_MODEL != 0 || STORE_PRECISION != 0 || USE_LES || STREAMING_MODEL != 0 +#error "Blocks require FP32 D2Q9 SRT double buffers without LES" #endif // Fixed-patch cell kinds, matching refinement.py ACTIVE/INNER_GHOST/OUTER_GHOST. @@ -143,7 +143,13 @@ extern "C" __global__ void BlockStep(const fpxx* in, fpxx* out, unsigned int hal } float rho, ux, uy, feq[NQ]; compute_rho_u(f, rho, ux, uy); +#if USE_DDF_SHIFTING + // Preserve the density perturbation before adding the unit reference. + const float drho = f[0]+f[1]+f[2]+f[3]+f[4]+f[5]+f[6]+f[7]+f[8]; + compute_feq(rho, ux, uy, feq, drho); +#else compute_feq(rho, ux, uy, feq); +#endif collide_srt_no_force(f, feq, d_params.omega); } } diff --git a/src/CelerisLab/lbm/kernels/boundary/planar_velocity.cuh b/src/CelerisLab/lbm/kernels/boundary/planar_velocity.cuh index 820950d..84e5d05 100644 --- a/src/CelerisLab/lbm/kernels/boundary/planar_velocity.cuh +++ b/src/CelerisLab/lbm/kernels/boundary/planar_velocity.cuh @@ -1,6 +1,6 @@ // On-site D2Q9 Zou-He velocity closure, paired lattice ordering. // Input is poststream f at the arrival time. Only unknown populations change; -// the caller owns collision/store. No force, shifting or off-lattice walls. +// the caller owns collision/store. No force or off-lattice walls. Shifted storage uses unit reference density. #ifndef CELERIS_BOUNDARY_PLANAR_VELOCITY_CUH #define CELERIS_BOUNDARY_PLANAR_VELOCITY_CUH #include @@ -28,7 +28,11 @@ __device__ __forceinline__ float bottom_wall_normal_velocity(unsigned int x) __device__ __forceinline__ void reconstruct_bottom_velocity(float* f, float u, float v) { - const float rho = (f[0] + f[1] + f[2] + 2.0f*(f[4] + f[6] + f[7])) / (1.0f-v); + const float rho = ( +#if USE_DDF_SHIFTING + 1.0f + +#endif + f[0] + f[1] + f[2] + 2.0f*(f[4] + f[6] + f[7])) / (1.0f-v); f[3] = f[4] + (2.0f/3.0f)*rho*v; f[5] = f[6] + 0.5f*(f[2]-f[1]) + 0.5f*rho*u + (1.0f/6.0f)*rho*v; f[8] = f[7] + 0.5f*(f[1]-f[2]) - 0.5f*rho*u + (1.0f/6.0f)*rho*v; @@ -36,7 +40,11 @@ __device__ __forceinline__ void reconstruct_bottom_velocity(float* f, float u, f __device__ __forceinline__ void reconstruct_top_velocity(float* f, float u, float v) { - const float rho = (f[0] + f[1] + f[2] + 2.0f*(f[3] + f[5] + f[8])) / (1.0f+v); + const float rho = ( +#if USE_DDF_SHIFTING + 1.0f + +#endif + f[0] + f[1] + f[2] + 2.0f*(f[3] + f[5] + f[8])) / (1.0f+v); f[4] = f[3] - (2.0f/3.0f)*rho*v; f[6] = f[5] + 0.5f*(f[1]-f[2]) - 0.5f*rho*u - (1.0f/6.0f)*rho*v; f[7] = f[8] + 0.5f*(f[2]-f[1]) + 0.5f*rho*u - (1.0f/6.0f)*rho*v; @@ -46,7 +54,12 @@ __device__ __forceinline__ void reconstruct_top_velocity(float* f, float u, floa // J_y=rho*v follows from known poststream populations (ordinary D2Q9). __device__ __forceinline__ void reconstruct_top_reservoir(float* f, float rho) { - const float jy = f[0]+f[1]+f[2]+2.0f*(f[3]+f[5]+f[8])-rho; +#if USE_DDF_SHIFTING + const float density_offset = rho - 1.0f; +#else + const float density_offset = rho; +#endif + const float jy = f[0]+f[1]+f[2]+2.0f*(f[3]+f[5]+f[8])-density_offset; f[4] = f[3] - (2.0f/3.0f)*jy; f[6] = f[5] + 0.5f*(f[1]-f[2]) - (1.0f/6.0f)*jy; f[7] = f[8] + 0.5f*(f[2]-f[1]) - (1.0f/6.0f)*jy; diff --git a/src/CelerisLab/lbm/kernels/operators/macro.cuh b/src/CelerisLab/lbm/kernels/operators/macro.cuh index 5fba6d0..f4d43a6 100644 --- a/src/CelerisLab/lbm/kernels/operators/macro.cuh +++ b/src/CelerisLab/lbm/kernels/operators/macro.cuh @@ -64,11 +64,12 @@ __device__ __forceinline__ void compute_rho_u(const float* __restrict__ f, #if LATTICE_MODEL == LATTICE_D2Q9 __device__ __forceinline__ void compute_feq(float rho, float ux, float uy, - float* __restrict__ feq) -{ + float* __restrict__ feq #if USE_DDF_SHIFTING - const float rhom1 = rho - 1.0f; // perturbation density + , float rhom1 #endif + ) +{ const float ux3 = 3.0f * ux; const float uy3 = 3.0f * uy; const float c3 = -1.5f * (ux * ux + uy * uy); // = -|u|²/(2cs²) @@ -112,6 +113,16 @@ __device__ __forceinline__ void compute_feq(float rho, float ux, float uy, #undef FEQ_CORE } +#if USE_DDF_SHIFTING +// Existing callers recover the offset from rho; blocked shifted collision +// supplies sum(g) directly to the overload above, retaining small perturbations. +__device__ __forceinline__ void compute_feq(float rho, float ux, float uy, + float* __restrict__ feq) +{ + compute_feq(rho, ux, uy, feq, rho - 1.0f); +} +#endif + #elif LATTICE_MODEL == LATTICE_D3Q19 __device__ __forceinline__ void compute_feq(float rho, float ux, float uy, float uz, diff --git a/src/CelerisLab/lbm/refinement.py b/src/CelerisLab/lbm/refinement.py index 8fe30b2..78b6a31 100644 --- a/src/CelerisLab/lbm/refinement.py +++ b/src/CelerisLab/lbm/refinement.py @@ -52,7 +52,7 @@ class FixedPatch2D: def from_config(cls, cfg: LBMConfig, bounds) -> "FixedPatch2D": """Prepare maps without changing cfg or allocating coarse-sized tables. - The initial backend contract is D2Q9/SRT/FP32/unshifted/double_buffer, + The initial backend contract is D2Q9/SRT/FP32/double_buffer, without LES. Two coarse cells of margin keep all halo donors in bounds; this does NOT certify that those donors are fluid or periodic. """ @@ -61,8 +61,8 @@ class FixedPatch2D: if (not checked.is_d2q9 or checked.collision != "SRT" or checked.streaming != "double_buffer" or checked.store_precision != "FP32" - or checked.ddf_shifting or checked.les_enabled): - raise ValueError("FixedPatch2D requires D2Q9/SRT/FP32/unshifted/double_buffer without LES") + or checked.les_enabled): + raise ValueError("FixedPatch2D requires D2Q9/SRT/FP32/double_buffer without LES") bounds = tuple(bounds) if len(bounds) != 4 or any(isinstance(v, (bool, np.bool_)) or not isinstance(v, Integral) for v in bounds): @@ -452,7 +452,7 @@ def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=(), checked.validate() if (patch.coarse_shape != (cfg.ny, cfg.nx) or not checked.is_d2q9 or checked.collision != "SRT" or checked.streaming != "double_buffer" - or checked.store_precision != "FP32" or checked.ddf_shifting or checked.les_enabled + or checked.store_precision != "FP32" or checked.les_enabled or patch.tau_coarse != 3*checked.viscosity+.5): raise ValueError("Fine boundary requires the matching fixed-patch configuration") from ..body.geometry.circle import CircleGeometry diff --git a/tests/unit/cpp/shifted_block_operators.cpp b/tests/unit/cpp/shifted_block_operators.cpp new file mode 100644 index 0000000..b4d2768 --- /dev/null +++ b/tests/unit/cpp/shifted_block_operators.cpp @@ -0,0 +1,36 @@ +#include +#include +#include +#define __device__ +#define __constant__ +#define __forceinline__ inline +#define LATTICE_D2Q9 0 +#define LATTICE_MODEL 0 +#define NQ 9 +#define DIM 2 +#define USE_DDF_SHIFTING 1 +#define WALL_NORMAL_VELOCITY 0.f +#define WALL_SLOT_WIDTH 0 +#define NX 16 +#include "core/descriptors.cuh" +#include "operators/macro.cuh" +#include "operators/collision_srt.cuh" +#include "boundary/planar_velocity.cuh" +void emit(const float* f) {for(int q=0;q<9;q++) printf("%.9g ",f[q]);puts("");} +int main() { + for(float delta : {0.f, 0x1p-30f, -0x1p-30f, .0002f}) + for(float u : {0.f,.005f}) { + float eq[9]; compute_feq(1.f+delta,u,-.2f*u,eq,delta); + printf("E %.9g %.9g %.9g\n",delta,u,-.2f*u);emit(eq); + } + for(int mode : {0,1,2}) for(float rho : {1.f,1.04f}) { + float f[9];const float u=.007f,v=.003f; + compute_feq(rho,u,v,f,rho-1.f); + for(int q=0;q<9;q++)f[q]+=.00001f*(q-4); + printf("B %d %.9g %.9g %.9g\n",mode,rho,u,v);emit(f); + if(mode==0)reconstruct_bottom_velocity(f,u,v); + else if(mode==1)reconstruct_top_velocity(f,u,v); + else reconstruct_top_reservoir(f,rho); + emit(f); + } +} diff --git a/tests/unit/test_shifted_block_operators.py b/tests/unit/test_shifted_block_operators.py new file mode 100644 index 0000000..17ff295 --- /dev/null +++ b/tests/unit/test_shifted_block_operators.py @@ -0,0 +1,32 @@ +"""Actual-header CPU algebra, not CUDA arithmetic/performance evidence.""" +import numpy as np +from .test_numerical_operators import compile_run + + +def test_retained_density_and_shifted_planar_headers(tmp_path): + rows=compile_run(tmp_path,'shifted_block_operators.cpp',{}) + cx=np.array([0,1,-1,0,0,1,-1,1,-1]);cy=np.array([0,0,0,1,-1,1,-1,-1,1]) + w=np.array([4/9]+[1/9]*4+[1/36]*4) + it=iter(rows) + for header in it: + kind,*numbers=header.split() + a=np.fromstring(next(it),sep=' ',dtype=np.float32).astype(float) + if kind=='E': + delta,u,v=np.array(numbers,dtype=np.float32).astype(float) + rho=float(np.float32(1+delta));cu=3*(cx*u+cy*v) + expected=w*(delta+rho*(cu+.5*cu*cu-1.5*(u*u+v*v))) + np.testing.assert_allclose(a,expected,rtol=3e-6,atol=1e-10) + if u==0: + np.testing.assert_allclose(a.sum(),delta,rtol=5e-8,atol=0) + if delta: assert np.any(a!=0) + else: + mode=int(numbers[0]);rho,u,v=np.array(numbers[1:],dtype=np.float32).astype(float) + b=np.fromstring(next(it),sep=' ',dtype=np.float32).astype(float) + known=cy!=(1 if mode==0 else -1) + np.testing.assert_array_equal(a[known],b[known]) + physical_rho=1+b.sum() + if mode==2: + assert abs(physical_rho-rho)<2e-8 + assert abs(b@cx)<2e-9 + else: + np.testing.assert_allclose([b@cx,b@cy],physical_rho*np.array([u,v]),atol=3e-9,rtol=0)