360 lines
18 KiB
Diff
360 lines
18 KiB
Diff
--- a/src/CelerisLab/lbm/refinement_stepper.py
|
|
+++ b/src/CelerisLab/lbm/refinement_stepper.py
|
|
@@ -1,5 +1,5 @@
|
|
# CelerisLab/lbm/refinement_stepper.py
|
|
-"""Opt-in fixed periodic patch backend; no Simulation/channel/objects adapter.
|
|
+"""Opt-in fixed patch backend; no Simulation/objects adapter.
|
|
|
|
The caller owns the CudaContext and must close this backend before that context.
|
|
Only synchronized postcollision_g leaf populations are public. No checkpoint,
|
|
@@ -27,8 +27,12 @@
|
|
grid: tuple
|
|
|
|
|
|
-class PeriodicPatchStepper:
|
|
- """D2Q9/SRT/FP32, one unforced 2:1 patch in a doubly periodic coarse box.
|
|
+class FixedPatchStepper:
|
|
+ """D2Q9/SRT/FP32, one unforced internal 2:1 patch in an x-periodic box.
|
|
+
|
|
+ ``coarse_boundary`` is ``periodic`` (default) or ``slot_reservoir``:
|
|
+ a bottom on-site slot and top fixed-density reservoir from ``cfg``.
|
|
+ Physical walls remain coarse; fine halo donors must be ordinary fluid.
|
|
|
|
``initialize(coarse_g, fine_g)`` accepts physical Q-major (9,ny,nx) arrays
|
|
already in postcollision phase. It does not collide or build channel flags.
|
|
@@ -36,10 +40,20 @@
|
|
On a device/launch failure the instance becomes unusable: close and recreate
|
|
it, rather than pretending the partially advanced fine buffer rolled back.
|
|
"""
|
|
- def __init__(self, cfg, bounds, *, context):
|
|
+ def __init__(self, cfg, bounds, *, context, coarse_boundary="periodic"):
|
|
if HALO != 4:
|
|
raise RuntimeError("refinement.cu requires a four-cell fine halo")
|
|
+ if coarse_boundary not in ("periodic", "slot_reservoir"):
|
|
+ raise ValueError("coarse_boundary must be periodic or slot_reservoir")
|
|
self.patch = FixedPatch2D.from_config(cfg, bounds)
|
|
+ self._wall_aware = coarse_boundary == "slot_reservoir"
|
|
+ if self._wall_aware:
|
|
+ if (cfg.y_wall_bc != "on_site_velocity" or cfg.top_wall_bc != "reservoir"
|
|
+ or cfg.velocity != 0 or not 2 <= cfg.wall_slot_width < cfg.nx):
|
|
+ raise ValueError("slot_reservoir requires on_site_velocity, reservoir, "
|
|
+ "velocity=0 and 2 <= wall_slot_width < nx")
|
|
+ if self.patch.bounds[1] < 4 or self.patch.bounds[3] > cfg.ny - 4:
|
|
+ raise ValueError("slot_reservoir requires y0 >= 4 and y1 <= ny-4")
|
|
self.ctx = context
|
|
self._lock = RLock()
|
|
self._status = "new"
|
|
@@ -50,6 +64,9 @@
|
|
self._block = (cfg.threads_per_block, 1, 1)
|
|
fine_cfg = replace(cfg, nx=self.patch.fine_shape[1], ny=self.patch.fine_shape[0],
|
|
viscosity=2 * cfg.viscosity)
|
|
+ if self._wall_aware:
|
|
+ fine_cfg = replace(fine_cfg, y_wall_bc="bounce_back", top_wall_bc="velocity",
|
|
+ wall_slot_width=0, wall_normal_velocity=0.0, velocity=0.0)
|
|
with self.ctx.activate():
|
|
try:
|
|
self._stream = cuda.Stream()
|
|
@@ -102,7 +119,7 @@
|
|
|
|
def _require(self, *states):
|
|
if self._status not in states:
|
|
- raise RuntimeError(f"Periodic patch is {self._status}; expected {' or '.join(states)}")
|
|
+ raise RuntimeError(f"Fixed patch is {self._status}; expected {' or '.join(states)}")
|
|
|
|
@property
|
|
def step_count(self):
|
|
@@ -111,10 +128,12 @@
|
|
def _mask_nonleaves(self, coarse, fine):
|
|
x0, y0, x1, y1 = self.patch.bounds
|
|
coarse[:, y0:y1, x0:x1] = np.nan
|
|
+ if self._wall_aware:
|
|
+ coarse[:, (0, -1), :] = np.nan
|
|
fine[:, self.patch.fine_kind != ACTIVE] = np.nan
|
|
|
|
def initialize(self, coarse_g, fine_g):
|
|
- """Upload an explicit periodic G_0; invalid non-leaves are discarded."""
|
|
+ """Upload an explicit G_0; invalid non-leaves are discarded."""
|
|
with self._lock, self.ctx.activate():
|
|
self._require("new", "ready")
|
|
arrays = [np.array(a, dtype=np.float32, order="C", copy=True)
|
|
@@ -124,6 +143,8 @@
|
|
x0, y0, x1, y1 = self.patch.bounds
|
|
leaf = np.ones(self.patch.coarse_shape, dtype=bool)
|
|
leaf[y0:y1, x0:x1] = False
|
|
+ if self._wall_aware:
|
|
+ leaf[[0, -1]] = False
|
|
for array, mask in zip(arrays, (leaf, self.patch.fine_kind == ACTIVE)):
|
|
if not np.isfinite(array[:, mask]).all() or (array[:, mask].sum(axis=0) <= 0).any():
|
|
raise ValueError("Leaf populations must be finite with positive density")
|
|
@@ -158,7 +179,8 @@
|
|
f.kernel(f.spare, f.current, grid=f.grid, **kw)
|
|
c.kernel(c.current, c.spare, f.current, self._targets, self._children,
|
|
np.uint32(len(p.coalesce_coarse)),
|
|
- *(np.uint32(v) for v in p.bounds), grid=c.grid, **kw)
|
|
+ *(np.uint32(v) for v in p.bounds), np.uint32(self._wall_aware),
|
|
+ grid=c.grid, **kw)
|
|
self._stream.synchronize()
|
|
c.current, c.spare = c.spare, c.current
|
|
self._step_count += 1
|
|
@@ -216,3 +238,7 @@
|
|
except Exception:
|
|
if exc[0] is None:
|
|
raise
|
|
+
|
|
+
|
|
+# Preserve the original entry point and its periodic default.
|
|
+PeriodicPatchStepper = FixedPatchStepper
|
|
--- a/src/CelerisLab/lbm/kernels/refinement.cu
|
|
+++ b/src/CelerisLab/lbm/kernels/refinement.cu
|
|
@@ -1,4 +1,4 @@
|
|
-// Fixed, unforced, doubly periodic 2:1 patch. Separate compilation entry;
|
|
+// Fixed, unforced internal 2:1 patch. Separate compilation entry;
|
|
// default OneStep/channel kernels are deliberately not included.
|
|
#include <stdint.h>
|
|
#include <math.h>
|
|
@@ -10,6 +10,7 @@
|
|
#include "operators/macro.cuh"
|
|
#include "operators/collision_srt.cuh"
|
|
#include "streaming/pull_double_buffer.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 "Refinement requires unshifted FP32 D2Q9 SRT double buffers without LES"
|
|
@@ -69,13 +70,15 @@
|
|
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 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];
|
|
@@ -100,6 +103,10 @@
|
|
}
|
|
// 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);
|
|
}
|
|
--- a/tests/validation/multigrid/test_refinement_followup.py
|
|
+++ b/tests/validation/multigrid/test_refinement_followup.py
|
|
@@ -1,10 +1,11 @@
|
|
-"""Only review gaps: CPU fault guards, normal-build seed1/8 and rectangle/1.
|
|
+"""Bounded fixed-patch checks: fault guards, periodic geometry and coarse slot BC.
|
|
|
|
Use the existing bounded runner/environment; no audit-build rerun. The scalar
|
|
-rectangle oracle covers one interval, not a second general-purpose reference.
|
|
+oracle covers the fixed tau pair and two boundary modes used by these cases.
|
|
"""
|
|
from contextlib import nullcontext
|
|
from dataclasses import replace
|
|
+import json
|
|
from pathlib import Path
|
|
import resource
|
|
from threading import RLock
|
|
@@ -21,8 +22,8 @@
|
|
PeriodicPatchStepper)
|
|
|
|
|
|
-def scalar_interval(coarse, fine, bounds):
|
|
- """Independent G→G single interval, geometric links, no production tables."""
|
|
+def scalar_interval(coarse, fine, bounds, slot=None):
|
|
+ """Independent G→G interval; optional (amplitude, width) bottom slot/top rho=1."""
|
|
_, ny, nx = coarse.shape
|
|
_, fy, fx = fine.shape
|
|
x0, y0, x1, y1 = bounds
|
|
@@ -53,7 +54,7 @@
|
|
result = np.full(coarse.shape, np.nan)
|
|
for y in range(ny):
|
|
for x in range(nx):
|
|
- if x0 <= x < x1 and y0 <= y < y1:
|
|
+ if (x0 <= x < x1 and y0 <= y < y1) or (slot is not None and y in (0, ny-1)):
|
|
continue
|
|
values = []
|
|
for q in range(9):
|
|
@@ -64,11 +65,45 @@
|
|
else:
|
|
value = coarse[q, (y-cy[q]) % ny, (x-cx[q]) % nx]
|
|
values.append(value)
|
|
+ if slot is not None and y == 1:
|
|
+ amplitude, width = slot
|
|
+ s = 2*((x+.5)-nx/2)/width
|
|
+ v = amplitude * max(0.0, 1-s*s)**2
|
|
+ rho = (values[0]+values[1]+values[2]
|
|
+ + 2*(values[4]+values[6]+values[7]))/(1-v)
|
|
+ values[3] = values[4] + 2*rho*v/3
|
|
+ values[5] = values[6] + (values[2]-values[1])/2 + rho*v/6
|
|
+ values[8] = values[7] + (values[1]-values[2])/2 + rho*v/6
|
|
+ elif slot is not None and y == ny-2:
|
|
+ jy = values[0]+values[1]+values[2]+2*(values[3]+values[5]+values[8])-1
|
|
+ values[4] = values[3] - 2*jy/3
|
|
+ values[6] = values[5] + (values[1]-values[2])/2 - jy/6
|
|
+ values[7] = values[8] + (values[2]-values[1])/2 - jy/6
|
|
result[:, y, x] = collide(values, .8)
|
|
return result, current
|
|
|
|
|
|
class FaultGuardsCPU(unittest.TestCase):
|
|
+ def test_boundary_mode_and_margin_contract(self):
|
|
+ cfg = LBMConfig(nx=32, ny=24, viscosity=.1, velocity=0, wall_normal_velocity=.001,
|
|
+ wall_slot_width=8, y_wall_bc="on_site_velocity", top_wall_bc="reservoir")
|
|
+ context = Mock()
|
|
+ context.activate.side_effect = RuntimeError("validated before CUDA")
|
|
+ for mode, bounds in (("unknown", (12, 4, 20, 10)),
|
|
+ ("slot_reservoir", (12, 3, 20, 10)),
|
|
+ ("slot_reservoir", (12, 4, 20, 21))):
|
|
+ with self.subTest(mode=mode, bounds=bounds):
|
|
+ with self.assertRaises(ValueError):
|
|
+ implementation.FixedPatchStepper(cfg, bounds, context=context, coarse_boundary=mode)
|
|
+ context.activate.assert_not_called()
|
|
+ # Both exact channel limits and the original periodic two-cell margin
|
|
+ # must reach context activation, with no allocation or device required.
|
|
+ for mode, bounds in (("slot_reservoir", (12, 4, 20, 20)),
|
|
+ ("periodic", (2, 2, 6, 5))):
|
|
+ with self.subTest(mode=mode, bounds=bounds):
|
|
+ with self.assertRaisesRegex(RuntimeError, "validated before CUDA"):
|
|
+ implementation.FixedPatchStepper(cfg, bounds, context=context, coarse_boundary=mode)
|
|
+
|
|
def fake_backend(self):
|
|
backend = PeriodicPatchStepper.__new__(PeriodicPatchStepper)
|
|
backend.ctx = Mock(activate=nullcontext)
|
|
@@ -145,6 +180,126 @@
|
|
def tearDown(self):
|
|
rss = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / 1024
|
|
print(f"host_peak_rss={rss:.3f} MiB", flush=True)
|
|
+
|
|
+ def test_slot_reservoir_eight_steps(self):
|
|
+ nx, ny, bounds, amplitude, width = 32, 24, (12, 4, 20, 10), .001, 8
|
|
+ x0, y0, x1, y1 = bounds
|
|
+ fx, fy = 2*(x1-x0)+8, 2*(y1-y0)+8
|
|
+ cfg = LBMConfig(nx=nx, ny=ny, viscosity=.1, velocity=0, rho=1,
|
|
+ wall_normal_velocity=amplitude, wall_slot_width=width,
|
|
+ y_wall_bc="on_site_velocity", top_wall_bc="reservoir")
|
|
+ cx = np.array((0, 1, -1, 0, 0, 1, -1, 1, -1))
|
|
+ cy = np.array((0, 0, 0, 1, -1, 1, -1, -1, 1))
|
|
+ weights = np.array((4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36))
|
|
+
|
|
+ def slot_velocity(x):
|
|
+ return amplitude * np.maximum(0, 1-(2*(x-nx/2)/width)**2)**2
|
|
+
|
|
+ def equilibrium(shape, origin, dx):
|
|
+ ys = origin[1] + dx*(np.arange(shape[0])+.5)
|
|
+ xs = origin[0] + dx*(np.arange(shape[1])+.5)
|
|
+ v = slot_velocity(xs)[None, :] * (1-(ys[:, None]-1.5)/(ny-3))
|
|
+ cu = cy[:, None, None] * v
|
|
+ return np.array(weights[:, None, None]*(1+3*cu+4.5*cu*cu-1.5*v*v), np.float32)
|
|
+
|
|
+ coarse_leaf = np.ones((ny, nx), bool)
|
|
+ coarse_leaf[[0, -1]] = False
|
|
+ coarse_leaf[y0:y1, x0:x1] = False
|
|
+ fine_leaf = np.zeros((fy, fx), bool)
|
|
+ fine_leaf[4:-4, 4:-4] = True
|
|
+ initial = [equilibrium((ny, nx), (0, 0), 1),
|
|
+ equilibrium((fy, fx), (x0-2, y0-2), .5)]
|
|
+ initial[0][:, ~coarse_leaf] = np.nan
|
|
+ initial[1][:, ~fine_leaf] = np.nan
|
|
+ # A postcollision zero-mass/zero-momentum perturbation at a grazing
|
|
+ # source distinguishes closed-corner coupling from ordinary coarse pull.
|
|
+ initial[0][[5, 6], y0-1, x0] += 1e-5
|
|
+ initial[0][[7, 8], y0-1, x0] -= 1e-5
|
|
+ expected = tuple(a.astype(np.float64) for a in initial)
|
|
+
|
|
+ def discrete_mass(arrays):
|
|
+ return (arrays[0][:, coarse_leaf].sum(dtype=np.float64)
|
|
+ + .25*arrays[1][:, fine_leaf].sum(dtype=np.float64))
|
|
+
|
|
+ # Observe the actual generated headers, without changing generation or
|
|
+ # compilation. In particular the fine module must have slot BC disabled.
|
|
+ generated = []
|
|
+ generate = compiler.generate_config
|
|
+
|
|
+ def record_config(*args, **kwargs):
|
|
+ build = generate(*args, **kwargs)
|
|
+ metadata = json.loads((Path(build)/".celeris-build").read_text())
|
|
+ definitions = {}
|
|
+ for line in (Path(build)/"config/config_method.h").read_text().splitlines():
|
|
+ parts = line.split()
|
|
+ if len(parts) == 3 and parts[0] == "#define":
|
|
+ definitions[parts[1]] = parts[2]
|
|
+ generated.append((metadata["config"], definitions))
|
|
+ return build
|
|
+
|
|
+ max_error = max_bc_error = 0.0
|
|
+ with CudaContext(0) as context:
|
|
+ with patch.object(compiler, "generate_config", side_effect=record_config):
|
|
+ backend = implementation.FixedPatchStepper(cfg, bounds, context=context,
|
|
+ coarse_boundary="slot_reservoir")
|
|
+ with backend:
|
|
+ self.assertEqual(len(generated), 2)
|
|
+ fine_config, fine_defines = generated[1]
|
|
+ self.assertEqual((fine_config["nx"], fine_config["ny"]), (fx, fy))
|
|
+ for name in ("Y_WALL_BC", "TOP_WALL_BC", "WALL_SLOT_WIDTH",
|
|
+ "WALL_NORMAL_VELOCITY", "TOP_WALL_NORMAL_VELOCITY"):
|
|
+ self.assertEqual(float(fine_defines[name].removesuffix("f")), 0, name)
|
|
+ self.assertEqual(int(generated[0][1]["Y_WALL_BC"]), 2)
|
|
+ self.assertEqual(int(generated[0][1]["TOP_WALL_BC"]), 1)
|
|
+ # Independently derive all target links, including the eight
|
|
+ # whose coarse sources are uncovered, from physical midpoints.
|
|
+ targets, corners = set(), set()
|
|
+ for y, x in np.argwhere(coarse_leaf):
|
|
+ for q in range(1, 9):
|
|
+ if x0 <= x+.5-cx[q]/2 <= x1 and y0 <= y+.5-cy[q]/2 <= y1:
|
|
+ target = q*nx*ny+y*nx+x
|
|
+ targets.add(target)
|
|
+ sx, sy = x-cx[q], y-cy[q]
|
|
+ if not (x0 <= sx < x1 and y0 <= sy < y1):
|
|
+ corners.add(target)
|
|
+ self.assertEqual(len(corners), 8)
|
|
+ self.assertEqual(set(backend.patch.coalesce_coarse.tolist()), targets)
|
|
+ backend.initialize(*initial)
|
|
+ before = discrete_mass(initial)
|
|
+ for step in range(1, 9):
|
|
+ expected = scalar_interval(*expected, bounds, slot=(amplitude, width))
|
|
+ backend.step()
|
|
+ self.assertEqual(backend.step_count, step)
|
|
+ raw = [np.empty(a.shape, np.float32) for a in initial]
|
|
+ for level, array in zip(backend._levels, raw):
|
|
+ cuda.memcpy_dtoh(array, level.current)
|
|
+ # Read raw storage before masking: active NaNs cannot be
|
|
+ # hidden by the public download's non-leaf cleanup.
|
|
+ self.assertTrue(np.isfinite(raw[0][:, coarse_leaf]).all())
|
|
+ self.assertTrue(np.isfinite(raw[1][:, fine_leaf]).all())
|
|
+ self.assertTrue(np.isnan(raw[0][:, ~coarse_leaf]).all())
|
|
+ for got, want in zip(raw, expected):
|
|
+ np.testing.assert_array_equal(np.isnan(got), np.isnan(want))
|
|
+ max_error = max(max_error, float(np.nanmax(np.abs(got-want))))
|
|
+ self.assertLessEqual(max_error, POP_ATOL)
|
|
+ public = backend.download()
|
|
+ np.testing.assert_array_equal(public[0], raw[0])
|
|
+ np.testing.assert_array_equal(public[1][:, fine_leaf], raw[1][:, fine_leaf])
|
|
+ self.assertTrue(np.isnan(public[1][:, ~fine_leaf]).all())
|
|
+ walls = raw[0][:, [1, ny-2], :].astype(np.float64)
|
|
+ rho = walls.sum(axis=0)
|
|
+ ux = np.einsum("q,qyx->yx", cx, walls)/rho
|
|
+ vy = np.einsum("q,qyx->yx", cy, walls)/rho
|
|
+ max_bc_error = max(max_bc_error, float(np.abs(ux).max()),
|
|
+ float(np.abs(vy[0]-slot_velocity(np.arange(nx)+.5)).max()),
|
|
+ float(np.abs(rho[1]-1).max()))
|
|
+ self.assertLessEqual(max_bc_error, 2e-6)
|
|
+ change = discrete_mass(raw)-before
|
|
+ bottom_flux, top_flux = (rho*vy).sum(axis=1)
|
|
+ print(f"slot/reservoir steps=8 max_population={max_error:.9g} "
|
|
+ f"max_BC_error={max_bc_error:.9g} discrete_leaf_mass_change={change:+.9g} "
|
|
+ f"final_sum_rho_v_bottom={bottom_flux:.9g} top={top_flux:.9g}; "
|
|
+ "raw halos/corners/fine BC isolation PASS", flush=True)
|
|
|
|
def test_nonequilibrium_seed1_eight_steps(self):
|
|
spec = ref.make_spec(collision="srt")
|