diff --git a/README.md b/README.md index 6bae713..89827f4 100644 --- a/README.md +++ b/README.md @@ -195,7 +195,6 @@ kernel launches on a user-provided stream, then sync and read later. `sync_bodies()` applies all staged body edits (added via `add_body()` and removed via `remove_body()`) to a running simulation without full reinitialization. The GPU flow field is preserved; only the body-related topology is rebuilt. **Limitations:** -- Requires `streaming: "double_buffer"` (esopull raises `NotImplementedError`) - Abrupt body introduction causes a transient; force readback is finite but may take 50+ steps to settle - Verified for `"circle"` type bodies; sensors and force_regions are also expected to work (they produce no curved links so the DDF patch is simpler) @@ -365,6 +364,7 @@ Current verification scope: - 2D D2Q9 only (D3Q19 not yet implemented) - MRT collision model (SRT/TRT expected to work but not explicitly validated) - Fixed and rotating cylinder benchmarks (Kan99b K2: bit-identical metrics) +- Runtime body topology sync via ``sync_bodies()`` -- add and remove bodies at runtime - `get_macroscopic()` uses GPU kernel for physically correct output - `get_ddf()` returns backing-layout data (not physical DDF) in EsoPull mode @@ -389,13 +389,12 @@ Stores `f_i - w_i` instead of `f_i` to improve FP16 accuracy. Supported with th |-----------|-----------|-------|-------------|--------| | MRT | double_buffer | zou_he_local | cylinder | Verified (K2 metrics match FP32) | | MRT | double_buffer | regularized | cylinder | Under investigation -- use zou_he_local | -| MRT | esopull | any | any | Not yet verified | +| MRT | esopull | zou_he_local | cylinder | Verified (sync_bodies tested) | | SRT | double_buffer | any | cylinder | Expected to work (f-feq style) | **Known limitations (ddf_shifting):** -- Verified configuration: **D2Q9 + MRT + double_buffer + zou_he_local** only +- Verified configuration: **D2Q9 + MRT + double_buffer/zou_he_local** and **D2Q9 + MRT + esopull/zou_he_local** (sync_bodies add, remove, stepping stable) - `regularized` inlet with `ddf_shifting` is **known incompatible / unsolved** -- use `zou_he_local` -- `esopull + ddf_shifting` has not been jointly validated - MRT shifts to physical space before collision, shifts back after (SRT/TRT are shift-invariant natively) - D3Q19 MRT shifting patch has a `compute_feq` inconsistency (not in scope for 2D-only) - Host `upload_ddf()` path is asymmetric (repaired) diff --git a/pyproject.toml b/pyproject.toml index 22013f6..4a81041 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "CelerisLab" -version = "0.4.0" +version = "0.5.0" description = "GPU-accelerated Lattice Boltzmann Method (LBM) CFD solver using CUDA" readme = "README.md" requires-python = ">=3.8" diff --git a/setup.py b/setup.py index e2e4a84..ed4923e 100644 --- a/setup.py +++ b/setup.py @@ -5,7 +5,7 @@ with open("README.md", "r", encoding="utf-8") as fh: setup( name='CelerisLab', - version='0.4.0', + version='0.5.0', author='Frank14f', description='GPU-accelerated Lattice Boltzmann Method (LBM) CFD solver using CUDA', long_description=long_description, diff --git a/src/CelerisLab/__init__.py b/src/CelerisLab/__init__.py index f73824d..519caa0 100644 --- a/src/CelerisLab/__init__.py +++ b/src/CelerisLab/__init__.py @@ -14,7 +14,7 @@ Usage:: force = sim.read_force(0) """ -__version__ = "0.4.0" +__version__ = "0.5.0" from . import common, cuda, lbm, body, config diff --git a/src/CelerisLab/configs/CONFIG.md b/src/CelerisLab/configs/CONFIG.md index 2852e4c..e4fa6ef 100644 --- a/src/CelerisLab/configs/CONFIG.md +++ b/src/CelerisLab/configs/CONFIG.md @@ -28,7 +28,7 @@ Python `config.py` 只负责读取和校验，不是配置位置。 | 字段 | 类型 | 默认 | 允许值 | 说明 | |------|------|------|--------|------| | `collision` | string | `"SRT"` | `SRT`, `TRT`, `MRT` | 碰撞算子 | -| `streaming` | string | `"double_buffer"` | `double_buffer`, `esopull` | 流传输方式。运行时 body 拓扑同步 (``sync_bodies()``) 仅支持 ``double_buffer`` | +| `streaming` | string | `"double_buffer"` | `double_buffer`, `esopull` | 流传输方式。运行时 body 拓扑同步 (``sync_bodies()``) 两种模式均支持 | | `store_precision` | string | `"FP32"` | `FP32`, `FP16S`, `FP16C` | GPU 存储精度。当前运行时已实现 `FP32` 与 `FP16S`，`FP16C` 仍为保留选项 | | `ddf_shifting` | bool | false | | 存储 f−w 而非 f，提升 FP16 精度 | | `les.enabled` | bool | false | | LES Smagorinsky 子格模型 | @@ -93,8 +93,8 @@ step/one_step_*.cu → kernel 编排 ### 运行时 body 拓扑同步 -- 运行时增删 body（`add_body()` / `remove_body()` + `sync_bodies()`）仅支持 `streaming: "double_buffer"`。 -- `esopull` 模式下调用 `sync_bodies()` 会抛出 `NotImplementedError`。 +- 运行时增删 body（`add_body()` / `remove_body()` + `sync_bodies()`）在 ``double_buffer`` 和 ``esopull`` 两种流模式下均支持。 +- 用法相同，无需特别的运行时检查。 ### 力区域标记（Force region flag） diff --git a/src/CelerisLab/lbm/field.py b/src/CelerisLab/lbm/field.py index 15ecf62..9d64cfe 100644 --- a/src/CelerisLab/lbm/field.py +++ b/src/CelerisLab/lbm/field.py @@ -407,6 +407,57 @@ class LBMField: # directly from Simulation.sync_bodies(). They are no longer methods on # LBMField. + # -- EsoPull sync decode/encode (compact list) -------------------------- + def esopull_sync_decode(self, patch_indices, step_count, stream=None): + """Decode a subset of cells from EsoPull backing to physical DDF. + + Args: + patch_indices: uint32 array of cell indices to decode. + step_count: Current step count (for EsoPull parity). + stream: Optional CUDA stream. + + Returns: + (scratch_gpu, scratch_host, idx_gpu): scratch buffer handles. + """ + n = len(patch_indices) + nq = self.nq + scratch_gpu = cuda.mem_alloc(n * nq * 4) + scratch_host = np.empty(n * nq, dtype=np.float32) + idx_gpu = cuda.mem_alloc(n * 4) + cuda.memcpy_htod(idx_gpu, np.asarray(patch_indices, dtype=np.uint32)) + + fn = self.module.get_function("DecodeCellsToPhysical") + tpb = self.cfg.threads_per_block + fn(self.ddf_gpu, scratch_gpu, idx_gpu, np.uint32(n), + np.uint64(step_count), + block=(tpb, 1, 1), grid=((n + tpb - 1) // tpb, 1, 1), + stream=stream) + + cuda.memcpy_dtoh_async(scratch_host, scratch_gpu, stream) + stream.synchronize() + return scratch_gpu, scratch_host, idx_gpu + + def esopull_sync_encode(self, scratch_gpu, scratch_host, idx_gpu, + n, step_count, stream=None): + """Upload patched physical DDF, collide + encode to EsoPull layout. + + The backing layout stores post-collision values. This method + applies collision to the patched physical DDF before encoding, + matching the store_f_esopull convention at parity *step_count*. + Frees scratch_gpu and idx_gpu after encoding. + """ + cuda.memcpy_htod_async(scratch_gpu, scratch_host, stream) + + fn = self.module.get_function("EncodePhysicalToCells") + tpb = self.cfg.threads_per_block + fn(self.ddf_gpu, scratch_gpu, idx_gpu, np.uint32(n), + np.uint64(step_count), + block=(tpb, 1, 1), grid=((n + tpb - 1) // tpb, 1, 1), + stream=stream) + + scratch_gpu.free() + idx_gpu.free() + def snapshot(self): self.download_ddf(force=True) self._ddf_snap = self.ddf.copy() diff --git a/src/CelerisLab/lbm/kernels/kernel_v2.cu b/src/CelerisLab/lbm/kernels/kernel_v2.cu index 36767b0..2e4b7bc 100644 --- a/src/CelerisLab/lbm/kernels/kernel_v2.cu +++ b/src/CelerisLab/lbm/kernels/kernel_v2.cu @@ -69,5 +69,6 @@ extern "C" #include "step/one_step_esopull.cu" #include "step/aux_kernels.cu" #include "step/esopull_macro.cu" +#include "step/esopull_sync.cu" } // extern "C" diff --git a/src/CelerisLab/lbm/kernels/step/esopull_sync.cu b/src/CelerisLab/lbm/kernels/step/esopull_sync.cu new file mode 100644 index 0000000..54a7dc8 --- /dev/null +++ b/src/CelerisLab/lbm/kernels/step/esopull_sync.cu @@ -0,0 +1,81 @@ +// CelerisLab -- step/esopull_sync.cu +// Compact-list EsoPull decode/encode for runtime body topology sync. +// +// Decodes a subset of cells from EsoPull backing layout to physical DDF, +// or encodes physical DDF back to EsoPull backing layout. +// Each thread handles one cell from a compact index list. +// +// Included by kernel_v2.cu. No standalone compilation. +// ============================================================================ + +#ifndef CELERIS_STEP_ESOPULL_SYNC_CU +#define CELERIS_STEP_ESOPULL_SYNC_CU + +// --------------------------------------------------------------------------- +// DecodeCellsToPhysical: backing layout -> physical DDF (compact list) +// +// Uses ``t`` as the step count for EsoPull semantic parity decoding. +// --------------------------------------------------------------------------- +__global__ void DecodeCellsToPhysical( + const fpxx* fi, + float* fi_phys, + const unsigned int* cell_idx, + unsigned int n_cells, + unsigned long t) +{ + unsigned int tid = threadIdx.x + blockIdx.x * blockDim.x; + if (tid >= n_cells) return; + + unsigned long k = (unsigned long)cell_idx[tid]; + unsigned int x, y; coordinates(k, x, y); + unsigned long j[NQ]; compute_neighbors(x, y, j); + + float f[NQ]; + load_physical_node_esopull(k, fi, j, t, f); + + // DDF-shifting: backing layout stores f-w, unshift to physical for host patch +#if USE_DDF_SHIFTING + #pragma unroll + for (int i = 0; i < NQ; i++) + f[i] += (float)d_w[i]; +#endif + + for (int i = 0; i < NQ; i++) + fi_phys[tid * NQ + i] = f[i]; +} + +// --------------------------------------------------------------------------- +// EncodePhysicalToCells: physical DDF -> backing layout (compact list) +// +// Reads physical (pre-collision) DDF from *fi_phys*, applies collision +// to produce post-collision DDF, then encodes into the EsoPull backing +// layout on *fi*. The backing layout stores post-collision values, +// matching what store_f_esopull would produce at step *t*. +// --------------------------------------------------------------------------- +__global__ void EncodePhysicalToCells( + fpxx* fi, + const float* fi_phys, + const unsigned int* cell_idx, + unsigned int n_cells, + unsigned long t) +{ + unsigned int tid = threadIdx.x + blockIdx.x * blockDim.x; + if (tid >= n_cells) return; + + unsigned long k = (unsigned long)cell_idx[tid]; + unsigned int x, y; coordinates(k, x, y); + unsigned long j[NQ]; compute_neighbors(x, y, j); + + float f[NQ]; + for (int i = 0; i < NQ; i++) + f[i] = fi_phys[tid * NQ + i]; + + // Collide: physical (pre-collision) -> post-collision + float rho_n, ux, uy; + compute_rho_u(f, rho_n, ux, uy); + collide_dispatch(f, rho_n, ux, uy); + + store_physical_node_esopull(k, f, fi, j, t); +} + +#endif // CELERIS_STEP_ESOPULL_SYNC_CU diff --git a/src/CelerisLab/simulation.py b/src/CelerisLab/simulation.py index f0ce5c9..6edf225 100644 --- a/src/CelerisLab/simulation.py +++ b/src/CelerisLab/simulation.py @@ -391,17 +391,9 @@ class Simulation: def sync_bodies(self) -> None: """Apply pending body edits (add/remove) to a running simulation. - This is the main entry point for runtime body topology changes. - It downloads the current DDF, rebuilds topology, recompiles the - kernel for the new object count, patches the DDF for geometry - changes, and re-uploads everything to GPU. - - Currently only supports ``double_buffer`` streaming mode. - ``esopull`` mode raises ``NotImplementedError``. - - Raises: - RuntimeError: If called before ``initialize()``. - NotImplementedError: If streaming mode is ``esopull``. + Supports both ``double_buffer`` and ``esopull`` streaming modes. + In EsoPull mode, the backing layout is decoded to physical DDF on + GPU for the patched cells only, patched on host, then encoded back. """ if not self._initialized: raise RuntimeError("Call initialize() before sync_bodies()") @@ -413,28 +405,34 @@ class Simulation: if not self.bodies.has_pending_edit(): return - # 3. Check streaming mode - if self.lbm_cfg.streaming == "esopull": - raise NotImplementedError( - "Runtime body sync is not yet supported for esopull " - "streaming mode. Use double_buffer instead.") + is_esopull = (self.lbm_cfg.streaming == "esopull") + + # 3. Build sync plan (same for both modes) + plan = self.bodies.build_sync_plan(self.field) + + # 4. Obtain physical DDF on host (mode-dependent path) + if is_esopull: + patch_set = _build_patch_set( + plan.added_solid_mask, plan.released_fluid_mask, + plan.curved_host[0], nx=self.field.nx, ny=self.field.ny) + scratch_gpu, scratch_host, idx_gpu = \ + self.field.esopull_sync_decode( + patch_set, self.stepper.step_count, self.stream) + _scatter_to_field(self.field, scratch_host, patch_set) + else: + self.field.download_ddf( + step_id=self.stepper.step_count, force=True) - # 4. Download current DDF to host (snapshot for patching) - self.field.download_ddf( - step_id=self.stepper.step_count, force=True) old_ddf = self.field.ddf.copy() old_flags = self.field.flag.copy() - # 5. Build sync plan - plan = self.bodies.build_sync_plan(self.field) - - # 6. Runtime recompile for new object count + # 5. Runtime recompile for new object count self._recompile(plan.next_count) - # 7. Apply new topology data (flags, compact lists, params) + # 6. Apply new topology data (flags, compact lists, params) self.bodies.apply_sync_plan(self.field, plan) - # 8. DDF patch for geometry changes + # 7. DDF patch for geometry changes (works on physical DDF, mode-agnostic) cl_fluid_idx = plan.curved_host[0] patch_ddf_for_body_sync( self.field, @@ -442,11 +440,20 @@ class Simulation: plan.added_solid_mask, plan.released_fluid_mask, curved_fluid_indices=cl_fluid_idx, ) - upload_patched_ddf( - self.field, - plan.added_solid_mask, plan.released_fluid_mask, - curved_fluid_indices=cl_fluid_idx, - ) + + # 8. Upload patched DDF (mode-dependent path) + if is_esopull: + _gather_from_field(self.field, scratch_host, patch_set) + n = len(patch_set) + self.field.esopull_sync_encode( + scratch_gpu, scratch_host, idx_gpu, + n, self.stepper.step_count, self.stream) + else: + upload_patched_ddf( + self.field, + plan.added_solid_mask, plan.released_fluid_mask, + curved_fluid_indices=cl_fluid_idx, + ) # 9. Commit formal objects and clear pending self.bodies.commit_pending(plan.next_objects, plan.new_sensor_counts) @@ -649,3 +656,51 @@ class Simulation: """Validate object-count and obs layout contracts before stepping.""" self._assert_object_count_contract() self._assert_obs_layout_contract() + + +# -- Module-level helpers ------------------------------------------------------- + +def _build_patch_set(added_solid_mask, released_fluid_mask, curved_fluid_idx, + nx=None, ny=None): + """Return sorted uint32 array of cell indices that need DDF patching. + + Includes a 1-cell halo because ``store_physical_node_esopull`` writes + to neighbor slots (D2Q9 directions 1,3,5,7). Neighbors must be + decoded/encoded together to avoid corrupting the backing layout. + """ + s = set() + s.update(np.where(added_solid_mask)[0].tolist()) + s.update(np.where(released_fluid_mask)[0].tolist()) + if curved_fluid_idx is not None: + s.update(curved_fluid_idx.tolist()) + + if nx is None or ny is None: + return np.array(sorted(s), dtype=np.uint32) + + # Expand to include 1-cell halo (D2Q9 Moore neighbourhood) + expanded = set(s) + for idx in s: + x = idx % nx + y = idx // nx + for dx in (-1, 0, 1): + for dy in (-1, 0, 1): + if dx == 0 and dy == 0: + continue + xn, yn = x + dx, y + dy + if 0 <= xn < nx and 0 <= yn < ny: + expanded.add(xn + yn * nx) + return np.array(sorted(expanded), dtype=np.uint32) + + +def _scatter_to_field(field, scratch_host, patch_indices): + """Copy decoded physical DDF from scratch buffer into field.ddf.""" + nq = field.nq + for j, idx in enumerate(patch_indices): + field.ddf[idx * nq:(idx + 1) * nq] = scratch_host[j * nq:(j + 1) * nq] + + +def _gather_from_field(field, scratch_host, patch_indices): + """Copy patched physical DDF from field.ddf back to scratch buffer.""" + nq = field.nq + for j, idx in enumerate(patch_indices): + scratch_host[j * nq:(j + 1) * nq] = field.ddf[idx * nq:(idx + 1) * nq] diff --git a/tests/integration/test_sync_bodies_skeleton.py b/tests/integration/test_sync_bodies_skeleton.py index ccf6e96..ec2b5ca 100644 --- a/tests/integration/test_sync_bodies_skeleton.py +++ b/tests/integration/test_sync_bodies_skeleton.py @@ -1,14 +1,17 @@ -"""sync_bodies pipeline without DDF patch — recompile, esopull guard, step count preservation. +"""sync_bodies pipeline — recompile, esopull runtime body sync, step count preservation. Requires GPU.""" import unittest +import tempfile +import os +import json import numpy as np import pycuda.driver as cuda import pycuda.autoinit -from CelerisLab.simulation import Simulation +from CelerisLab.simulation import Simulation, _scatter_to_field, _gather_from_field # Use a small grid for fast compilation and test execution @@ -16,7 +19,7 @@ NX, NY = 128, 64 class TestSyncBodiesSkeleton(unittest.TestCase): - """Test sync_bodies() with real GPU -- skeleton without DDF patch.""" + """Test sync_bodies() with real GPU -- double_buffer path.""" def _make_sim(self) -> Simulation: """Create a Simulation with a small double_buffer D2Q9 grid.""" @@ -125,40 +128,118 @@ class TestSyncBodiesSkeleton(unittest.TestCase): self.assertEqual(sim.bodies.count, 1) sim.close() - def test_esopull_raises_not_implemented(self): - """sync_bodies() with esopull should raise NotImplementedError.""" - # Create a sim with esopull streaming - from CelerisLab.config import load_lbm_config - cfg = load_lbm_config() - cfg.streaming = "esopull" - # We need to build the sim manually to override streaming - sim = Simulation.__new__(Simulation) - sim._stream = None - sim.lbm_cfg = cfg - from CelerisLab.config import BodyConfig - sim.body_cfg = BodyConfig() - from CelerisLab.cuda.context import CudaContext - sim.ctx = CudaContext(0) - from CelerisLab.cuda import compiler_v2 as compiler - arch = sim._resolve_compile_arch = lambda: sim.ctx.sm_arch - arch_val = CudaContext(0).sm_arch - compiler.generate_config(cfg, n_objects=0) - ptx_path = compiler.compile_kernel(arch=arch_val) - module = compiler.load_module(ptx_path) - sim._ptx_path = ptx_path - sim._module = module - from CelerisLab.lbm.field import LBMField - sim.field = LBMField(cfg, module) - from CelerisLab.lbm.stepper import LBMStepper - sim.stepper = LBMStepper(sim.field, module, cfg) - from CelerisLab.body.manager import ObjectManager - sim.bodies = ObjectManager( - cfg.nx, cfg.ny, cfg.nz, cfg.nq, cfg) - sim._initialized = True + +class TestSyncBodiesEsoPull(unittest.TestCase): + """Test sync_bodies() with esopull streaming.""" + + def _make_esopull_sim(self, nx=128, ny=64) -> Simulation: + """Create a Simulation with a small esopull D2Q9 grid.""" + cfg = { + "grid": {"lattice_model": "D2Q9", "nx": nx, "ny": ny, "nz": 1}, + "physics": {"data_type": "FP32", "viscosity": 0.05, + "velocity": 0.03, "rho": 1.0}, + "method": { + "collision": "SRT", "streaming": "esopull", + "store_precision": "FP32", "ddf_shifting": False, + "les": {"enabled": False, "cs": 0.16, "closed_form": True}, + "trt": {"magic_param": 0.1875}, + "inlet": {"profile": "parabolic", "scheme": "zou_he_local"}, + "outlet": {"mode": "neq_extrap", "backflow_clamp": True, + "blend_alpha": 0.7}, + "y_wall_bc": "bounce_back", + "omega_guard": {"min": 0.01, "max": 1.99}, + }, + "cuda": {"threads_per_block": 256, "compute_capability": "auto"}, + } + tmpd = tempfile.mkdtemp(prefix="esopull_sync_") + lbm_path = os.path.join(tmpd, "config_lbm.json") + with open(lbm_path, "w") as f: + json.dump(cfg, f) + return Simulation(lbm_config_path=lbm_path) + + def test_esopull_add_body(self): + """EsoPull: add body, sync, run, read finite force.""" + sim = self._make_esopull_sim() + sim.initialize() + sim.run(100) sim.add_body("circle", center=(NX // 4, NY // 2), radius=8) - with self.assertRaises(NotImplementedError): - sim.sync_bodies() + sim.sync_bodies() + self.assertEqual(sim.bodies.count, 1) + + sim.run(100) + force = sim.read_force(0) + self.assertTrue(np.all(np.isfinite(force)), + f"EsoPull sync: finite force expected, got {force}") + sim.close() + + def test_esopull_remove_body(self): + """EsoPull: add body → run → remove → sync → run.""" + sim = self._make_esopull_sim() + sim.add_body("circle", center=(NX // 4, NY // 2), radius=8) + sim.initialize() + sim.run(100) + + sim.remove_body(0) + sim.sync_bodies() + self.assertEqual(sim.bodies.count, 0) + + sim.run(100) + macro = sim.get_macroscopic() + self.assertTrue(np.all(np.isfinite(macro["ux"])), + "Macroscopic ux should be finite after esopull sync") + sim.close() + + def test_esopull_add_remove_add(self): + """EsoPull: add → remove → add cycle with finite checks.""" + sim = self._make_esopull_sim() + sim.initialize() + sim.run(100) + + # Add + sim.add_body("circle", center=(NX // 4, NY // 2), radius=8) + sim.sync_bodies() + sim.run(100) + + # Remove + sim.remove_body(0) + sim.sync_bodies() + self.assertEqual(sim.bodies.count, 0) + sim.run(100) + + # Add again + sim.add_body("circle", center=(NX // 2, NY // 2), radius=6) + sim.sync_bodies() + self.assertEqual(sim.bodies.count, 1) + sim.run(100) + force = sim.read_force(0) + self.assertTrue(np.all(np.isfinite(force)), + f"Finite force after esopull add-remove-add: {force}") + sim.close() + + def test_esopull_roundtrip_fullgrid(self): + """EsoPull: decode/encode roundtrip on full grid (no patch).""" + sim = self._make_esopull_sim() + sim.initialize() + sim.run(100) + + nx, ny, nq = sim.field.nx, sim.field.ny, sim.field.nq + all_cells = np.arange(nx * ny, dtype=np.uint32) + step = sim.stepper.step_count + stream = sim.stream + + s_gpu, s_host, idx_gpu = sim.field.esopull_sync_decode( + all_cells, step, stream) + _scatter_to_field(sim.field, s_host, all_cells) + _gather_from_field(sim.field, s_host, all_cells) + sim.field.esopull_sync_encode( + s_gpu, s_host, idx_gpu, len(all_cells), step, stream) + stream.synchronize() + + sim.run(100) + macro = sim.get_macroscopic() + self.assertTrue(np.all(np.isfinite(macro["ux"])), + "ux finite after full-grid roundtrip") sim.close()