CelerisLab/README.md

CelerisLab

GPU-Accelerated Lattice Boltzmann Method (LBM) CFD Solver

CelerisLab is a high-performance computational fluid dynamics solver based on the Lattice Boltzmann Method, leveraging NVIDIA CUDA for GPU acceleration. It provides a Python API for scripting, real-time control loop integration, and scientific workflow automation.

Features

• GPU Acceleration: CUDA kernels for high-performance simulation (384x192 D2Q9: ~4400 MLUPS on V100)
• D2Q9 / D3Q19 Lattice: 2D and 3D lattice implementations
• Multiple Collision Models: SRT, TRT, and MRT operators; Smagorinsky LES subgrid model
• Dual Streaming Paths: Standard double-buffer pull and memory-efficient esoteric-pull (EsoPull). EsoPull is verified as numerically equivalent to double-buffer for D2Q9 curved-boundary MRT (Kan99b K2 validation).
• Curved Boundary Bouzidi: Immersed boundary support for complex geometries with wall velocity control
• Flexible Boundary Conditions: NEQ-extrapolation pressure outlet, parabolic/uniform velocity inlet, half-way bounce-back walls
• Rotating Body Control: Real-time setting of body rotation speeds via sim.set_body()
• Force / Torque / Sensor Readback: On-demand force, torque, and area-averaged sensor velocity
• Physics Validated: Strouhal numbers match Sah04 (confined cylinder) and Kan99b (rotating cylinder) references

Quick Start

Single cylinder

from CelerisLab import Simulation

sim = Simulation()
sim.add_body("circle", center=(50, 50), radius=10)
sim.initialize()

for step in range(10000):
    sim.run(1)

macro = sim.get_macroscopic()  # {"rho": ..., "ux": ..., "uy": ...}
force = sim.read_force(0)       # [fx, fy] on body 0
sim.close()

DRL control loop

from CelerisLab import Simulation

sim = Simulation()
sim.add_body("circle", center=(256, 128), radius=10)
sim.add_body("sensor", center=(300, 128), radius=10)
sim.initialize()

for episode in range(100):
    # Step the simulation (auto uploads action, downloads obs)
    sim.run(100)

    # Read individual body telemetry (primary API)
    data = sim.read_body(0)
    print(f"step={sim.stepper.step_count} "
          f"force=({data.force[0]:.4f},{data.force[1]:.4f}) "
          f"sensor=({data.sensor[0]:.4f},{data.sensor[1]:.4f})")

    # DRL policy inference (replace with your model)
    action_omega = 0.001 * (0.5 - data.force[0])

    # Set action (host-only, will be auto-uploaded next run)
    sim.set_body(0, omega=action_omega)

sim.close()

Async control (performance-oriented, custom stream)

import pycuda.driver as cuda

stream = cuda.Stream()
sim.set_body(0, omega=0.002)  # host-only
sim.run(100, stream=stream)   # action uploaded, steps run, obs downloaded on stream
# stream is synced inside run() -- obs is ready
data = sim.read_body(0)

Manual stream control (max overlap)

import pycuda.driver as cuda

stream = cuda.Stream()
# Skip transfers for the first batch, just enqueue kernels
sim.run(100, stream=stream, upload_act=False, sync_obs=False)

# ... other GPU work can overlap with the kernel launches ...

# Later: sync and read
stream.synchronize()
obs = sim.read_bodies(stream=stream)  # sync already done, just read pinned buffer

## Installation

### Prerequisites

- Python 3.8+
- NVIDIA GPU with CUDA Compute Capability 6.0+
- CUDA Toolkit 11.0+
- NVIDIA drivers

### Install from source

```bash
git clone <repository_url>
cd CelerisLab
pip install -e .

Dependencies

• pycuda>=2020.1 — CUDA Python bindings
• numpy>=1.19.0 — numerical computing
• scipy>=1.5.0 — special functions for vortex initialization

API Reference

Simulation

sim = Simulation(
    lbm_config_path: Optional[str] = None,   # path to config JSON
    body_config_path: Optional[str] = None,   # path to body config JSON
    device_id: int = 0,                       # GPU device index
)

Body creation

Method	Returns	Description
sim.add_body(type="circle", center=(x,y), radius=r)	int body_id	Add a cylinder body (primary API)
sim.add_body(type="sensor", center=(x,y), radius=r)	int body_id	Add a velocity sensor
sim.add_body(type="force_region", center=(x,y), radius=r)	int body_id	Add a force application region
sim.add_cylinder(center, radius)	int body_id	Convenience wrapper (deprecated)
sim.add_sensor(center, radius)	int body_id	Convenience wrapper (deprecated)
sim.add_object(obj)	int body_id	Add pre-configured SimObject

Future geometry types (polygon, mesh) will use the same add_body() function with a different type parameter.

Runtime control

Method	Description
sim.initialize()	Recompile if needed, flow field + sync objects to GPU
sim.run(steps, *, upload_act=True, sync_obs=True, zero_obs=True, stream=None)	Run N LBM steps. See stream subsection below.
sim.set_body(id, omega=...)	Set body rotation speed (host array only, uploaded at next run())
sim.read_body(id, *, normalize=True) -> BodyTelemetry	Unified telemetry: {force, torque, sensor} from pinned buffer
sim.read_bodies() -> ndarray	Flat array of all bodies' telemetry (batch DRL read)
sim.read_force(id, *, normalize=True) -> ndarray	Force vector [fx, fy]
sim.read_torque(id, *, normalize=True) -> ndarray	Torque [tz]
sim.read_sensor(id, *, normalize=True) -> ndarray	Area-averaged velocity; time-normalised when normalize=True
sim.set_force(id, fx=..., fy=...)	Set force density on a force_region object

Action/obs transfer model: set_body() / set_force() are host-only — they modify the host action array without triggering GPU upload. The GPU buffer is automatically updated at the start of the next run() call when upload_act=True (the default).

Obs telemetry model: GPU kernels accumulate force, torque, and sensor readings into the obs_gpu buffer via atomicAdd. By default, run(zero_obs=True) clears the entire obs_gpu buffer (all three segments) and resets an internal step counter before stepping. After the step group, telemetry is downloaded to a pinned host buffer when sync_obs=True.

All three readback methods accept a normalize keyword:

• normalize=True (default): divides the raw GPU value by the accumulated step count, yielding a per-step average — the physically meaningful quantity for most use cases.
• normalize=False: returns the raw GPU-accumulated sum (no time division).

Sensor special handling: Area-normalisation (dividing by the number of sensor cells) is always applied internally in read_sensor(), regardless of the normalize flag. The normalize parameter only controls the additional time-normalisation step.

run() parameters:

• steps: Number of LBM steps.
• upload_act (default True): Upload host action array to action_gpu before stepping.
• sync_obs (default True): Download obs_gpu to host pinned buffer after stepping.
• zero_obs (default True): Zero all obs segments (force, torque, sensor) on GPU and reset the step accumulator before the step group. Set False to accumulate telemetry across multiple run() calls.
• stream (default None): CUDA stream for all operations. None uses an internal stream.
• checkpoint_interval (default 0): If >0, save an HDF5 checkpoint every N steps.

Use upload_act=False, sync_obs=False to skip all transfers and enqueue pure kernel launches on a user-provided stream, then sync and read later.

Runtime body topology sync

Method	Description
sim.remove_body(id)	Stage a body for removal (committed at next sync_bodies())
sim.sync_bodies()	Commit pending add/remove edits: recompile kernel, rebuild flags/compact lists, patch DDF, re-upload to GPU

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:

• 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)

# Add a body to an already-initialized simulation
sim = Simulation()
sim.initialize()
sim.run(500)
sim.add_body("circle", center=(256, 128), radius=10)
sim.sync_bodies()     # recompile + patch
sim.run(500)
force = sim.read_force(0)

# Remove the same body at runtime
sim.remove_body(0)
sim.sync_bodies()     # recompile + patch flags/DDF
sim.run(500)

Note: If run() is called without a preceding sync_bodies(), any staged edits are silently discarded.

force_region usage

# Create a circular force application region
fr_id = sim.add_body("force_region", center=(50, 50), radius=15)

# Set force density (lattice units, implicit GPU upload)
sim.set_force(fr_id, fx=0.001, fy=0.0)

# The region applies Guo forcing on each step. Zero force = no-op.
sim.set_force(fr_id, fx=0.0, fy=0.0)  # disable force

Persistence note: set_force() only updates the host action array. The GPU buffer is synced at the next run() call. If sync_to_gpu() is called manually before run(), the force will be reset to zero. For the common usage pattern (initialize -> set_force -> run -> set_force -> run ...), this is not an issue.

Comparison: body types

Type	Flag overlay	Produces cut-links	Readback	Runtime control
"circle"	OBSTACLE + BC_CURVED	Yes (Bouzidi)	force/torque	set_body(id, omega=...)
"sensor"	FLUID + SENSOR_FLAG	No	area-averaged velocity (always); optional per-step average	None needed
"force_region"	FLUID + FRC_REGION	No	None	set_force(id, fx=..., fy=...)

Data access

Method	Description
sim.get_macroscopic()	Download DDF, return dict with rho/ux/uy
sim.get_ddf()	Download raw DDF array
sim.get_flags()	Copy host-side flag array
sim.update_runtime_params(omega=..., fx=..., fy=...)	Update runtime constants without recompile

Checkpoint / Snapshot

Method	Description
sim.save_checkpoint(path) -> str	HDF5 checkpoint with full state
sim.load_checkpoint(path)	Restore from HDF5 (config must match)
sim.snapshot() / sim.restore()	In-memory field snapshot

Low-level access

Attribute	Description
sim.bodies	ObjectManager for direct GPU buffer access (action_gpu, obs_gpu)
sim.stream	Internal CUDA stream for async operations
sim.field	LBMField (GPU memory + curved/sensor SoA handles)
sim.stepper	LBMStepper for fine-grained step control

LBMStepper (advanced usage)

stepper.step(n=1, *, action_gpu, obs_gpu, stream=None)

When fine-grained control is needed (e.g., custom async patterns), step manually:

stream = cuda.Stream()
sim.bodies.zero_obs_async(stream)
sim.stepper.step(
    1,
    action_gpu=sim.bodies.action_gpu,
    obs_gpu=sim.bodies.obs_gpu,
    stream=stream,
)
stream.synchronize()
sim.bodies.increment_obs_steps(1)   # manually track steps for normalize
force = sim.read_force(0)           # normalize=True: divides by 1 step

Configuration

Config file location

Simulation() resolves config_lbm.json in this order:

1. Explicit path argument to Simulation(path)
2. $CELERISLAB_CONFIG_DIR/config_lbm.json
3. ./configs/config_lbm.json (current working directory)
4. The copy shipped inside the installed package

Config structure

{
  "grid": {
    "lattice_model": "D2Q9",
    "nx": 512, "ny": 256, "nz": 1
  },
  "physics": {
    "data_type": "FP32",
    "viscosity": 0.0035,
    "velocity": 0.03,
    "rho": 1.0
  },
  "method": {
    "collision": "SRT",
    "streaming": "double_buffer",
    "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",
      "trt_neq_damp": 0.5,
      "regularized_neq_damp": 0.5
    },
    "outlet": {
      "mode": "neq_extrap",
      "backflow_clamp": true,
      "blend_alpha": 0.7,
      "srt_neq_damp": 0.5
    },
    "y_wall_bc": "bounce_back",
    "omega_guard": { "min": 0.01, "max": 1.99 }
  },
  "cuda": {
    "threads_per_block": 256,
    "compute_capability": "auto"
  }
}

Full parameter documentation lives in src/CelerisLab/configs/CONFIG.md.

Performance

Benchmarks (V100, D2Q9, 384x192)

Config	Streaming	MLUPS
Re100 MRT noLES	double_buffer	~4400
Re100 MRT noLES	esopull	~4400

EsoPull streaming mode

EsoPull (Esoteric-Pull) is a single-buffer streaming scheme that uses half the memory of double-buffer. It is fully supported for 2D D2Q9 with curved boundaries, rotating cylinders, sensors, and force regions.

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

Enable via config: "streaming": "esopull"

FP16S store precision

Half-precision storage is supported for the DDF buffer. All computations are performed in FP32; only storage uses FP16 with a scaling factor.

Verified benchmarks:

• Sah04 S2: St error within 1.5% (channel + curved + inlet/outlet)
• Kan99b K2: Shows quantization sensitivity (St ~16% deviation from FP32 at Re=100)
• High-blockage cases (S4 beta=0.9): May diverge earlier than FP32

Enable via config: "store_precision": "FP16S"

ddf_shifting mode

Stores f_i - w_i instead of f_i to improve FP16 accuracy. Supported with the following verified combinations:

Collision	Streaming	Inlet	Curved body	Status
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	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 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
• 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)
• Checkpoint now enforces streaming and ddf_shifting match

Performance characteristics

The GPU is the primary runtime cost. Python overhead is minimal.

384x192 (validation grid): GPU kernel time is ~78 μs/step, of which OneStep is ~5.9 μs (MRT D2Q9). The remaining time is dominated by pycuda kernel launch overhead (~37 μs per launch).

3000x300 (production grid): Estimated GPU compute time is ~530 μs/step, with pycuda overhead fixed at ~111 μs, yielding ~83% GPU utilization.

sim.set_body() and sim.read_force() data transfers are negligible (~1 μs for 72 bytes).

For a detailed breakdown, see docs/performance_analysis.md.

Body Module Architecture

body/
  __init__.py          Package exports
  objects.py           SimObject container + ObjectState / ObjectControl
  manager.py           ObjectManager: GPU buffer lifecycle, sync, telemetry
  registry.py          BodyRegistry: pure add/remove/query
  action_smoother.py   ActionSmoother for control input ramping
  geometry/            Shape implementations (CircleGeometry, Geometry ABC)
  coupling/            Body-fluid coupling: SoA packing, force/torque
  preprocess/          Grid preprocessing: flag overlay, cut-link building

Module Boundaries

• body/ — geometry, rigid-body state, preprocessing, force/torque readback
• lbm/ — lattice Boltzmann kernels, field memory, stepper
• cuda/ — compilation pipeline, context lifecycle
• common/ — shared utilities (checkpoint, render, streakline pathline)

Geometry is separated from boundary methods. CircleGeometry produces geometry-agnostic CutLink records. The SoA packer (body/coupling/soa_packer.py) is the single point that knows the kernel memory layout. Adding a new shape (polygon, mesh) requires only a new Geometry subclass.

Validated Benchmarks

Benchmark	Description	Key metrics	Precision
Sah04 S1-S4	Confined stationary cylinder	Strouhal matching Sahin & Owens (2004)	St error < 5%
Kan99b K0-K5	Rotating cylinder in open domain	St, Cd, Cl matching Kang et al. (1999)	See tolerance table
Sensor accuracy	GPU sensor vs CPU flow-field average	Match to 1e-9	Verified

Run validation scripts:

conda run -n pycuda_3_10 python tests/validation/run_kan99b_rotating_cylinder.py
conda run -n pycuda_3_10 python tests/validation/run_sah04_st_matrix.py
conda run -n pycuda_3_10 python tests/validation/test_sensor_accuracy.py

Performance baseline

conda run -n pycuda_3_10 python tests/validation/run_perf_baseline.py \
  --lattice-model D2Q9 --nx 384 --ny 192 --collision MRT

Project Layout

See docs/tests_overview.md for a complete guide to the test suite.

src/CelerisLab/
  simulation.py          High-level API
  config.py              LBMConfig / BodyConfig dataclasses
  body/                  Object management, geometry, GPU sync
  cuda/                  CUDA context, compilation, PTX load
  lbm/                   Field, stepper, kernels (CUDA source)
  common/                Preprocess, checkpoint, render, streakline
tests/
  validation/            Regression runners (Kan99b, Sah04, sensor, perf)
  postproc/              Post-processing scripts (exp_ctrl_matrix, streakline)
  specs/                 Validation spec documents
  audit/                 Audit reports (archived, see docs/)
  output/                Test outputs (force CSV, vorticity PNG, checkpoints)
docs/
  performance_analysis.md   GPU/Python profiling report
  audit/                    Audit findings (round 1-2, kernel layer, body refactor notes)
  validation_specs/         Validation methodology documents
legacy/                   Superseded code (FlowField, compiler v1, macros.h)
ref/                      External reference implementations (FluidX3D)

Collision model recommendations

Use case	Recommended config
Low Re (<= 500)	SRT or TRT, LES off
Medium Re (500-2000)	MRT or SRT+LES
High Re (2000-5000)	MRT+LES (most robust); SRT+LES; TRT+LES with omega_guard.max in 1.90-1.99

Common control loop patterns

Sync control (simple)

sim.set_body(0, omega=0.002)
sim.run(10)
data = sim.read_body(0)

Async control (performance-oriented)

sim.set_body(0, omega=0.002)  # host-only, ~1 μs
sim.stepper.step(10, ..., stream=sim.stream)
sim.bodies.increment_obs_steps(10)             # track steps for normalize
sim.bodies.download_obs_full_async(sim.stream)
sim.stream.synchronize()
force = sim.read_force(0)                     # per-step average force

Use sim.run() for the common case -- it stores the step count automatically:

sim.set_body(0, omega=0.002)
sim.run(10, stream=sim.stream)
force = sim.read_force(0)                     # per-step average force

Vortex initialization

from CelerisLab.lbm.initializers import add_vortex
add_vortex(sim.field, center=(50, 50), radius=10.0, strength=1.0, vortex_type="lamb")

Streakline visualization

from CelerisLab.common.streakline import Streakline, ReleaseConfig, IntegratorConfig

streak = Streakline(release_points=..., nx=nx, ny=ny)
for step in range(steps):
    sim.run(1)
    if step % sample_every == 0:
        macro = sim.get_macroscopic()
        streak.observe(ux=macro["ux"], uy=macro["uy"], step=step)
streak.render("streakline.png")

Citation

@software{celerislab2026,
  author = {Frank14f},
  title = {CelerisLab: GPU-Accelerated Lattice Boltzmann Method Solver},
  year = {2026},
  url = {https://github.com/frank14f/CelerisLab}
}

License

MIT License — see LICENSE file for details.