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

Fixed 2:1 refinement patch (experimental backend)

FixedPatchStepper advances one fixed 2:1 refinement patch inside a coarse grid, optionally with stationary circles in the fine region and a prescribed normal surface velocity (blowing/suction). It is a standalone backend, not part of Simulation. Runnable demo (blocked fine storage; omit --fine-block for dense): PYTHONPATH="$PWD/src" python examples/fixed_patch_circle.py --fine-block 16 20.

import numpy as np
from CelerisLab.body.geometry.circle import CircleGeometry
from CelerisLab.body.objects import SimObject, SurfaceActuation
from CelerisLab.config import LBMConfig
from CelerisLab.cuda.context import CudaContext
from CelerisLab.lbm.equilibrium import compute_feq_d2q9
from CelerisLab.lbm.refinement_stepper import FixedPatchStepper

cfg = LBMConfig(nx=32, ny=24, viscosity=.1)
body = SimObject(0, geometry=CircleGeometry(16.0, 12.0, 2.15), center=(16.0, 12.0), radius=2.15,
                 surface_actuation=SurfaceActuation(0.7, 0.85, 0.002))  # angle, half-width, speed
with CudaContext(0) as ctx, FixedPatchStepper(cfg, context=ctx, refinements=[(8, 6, 24, 18)],
                                              level_blocks=[None, (16, 20)],
                                              level_objects=[(), (body,)]) as run:
    rest = np.asarray(compute_feq_d2q9(1.0, 0.0, 0.0), np.float32)[:, None, None]
    # (9, ny, nx) postcollision populations per level; non-leaf cells are discarded.
    run.initialize(np.broadcast_to(rest, (9, *run.hierarchy.pairs[0].coarse_shape)),
                   np.broadcast_to(rest, (9, *run.hierarchy.pairs[0].fine_shape)))
    run.step(128)                      # coarse steps; fine level takes two per step
    coarse_g, fine_g = run.download()  # NaN on covered, ghost and solid cells

Scope and limits:

• The examples below use D2Q9, SRT, FP32, double buffer, unshifted, no LES or forcing. Fixed 2:1 hierarchies use level lists of arbitrary depth; two- and three-level examples below use that same constructor.

• coarse_boundary="periodic" (default) or "slot_reservoir" (bottom slot, top reservoir; patch needs y0 >= 4, y1 <= ny-4). Periodic mode ignores the wall/slot config fields.

• Object centers/radii use the coarse cell-edge frame (cell I spans [I, I+1)), unlike Simulation's integer-node frame. Surface speeds are in coarse lattice units.

• Bodies must be stationary circles whose whole boundary stencil lies in the fine active region, clear of the coarse/fine transfer cells; violations are rejected before any GPU allocation.

• level_blocks=[None, (by, bx)] stores the fine level in fixed-size blocks; by/bx must divide the WHOLE fine allocation, i.e. 2*(y1-y0)+8 by 2*(x1-x0)+8 including the four-cell ghost band. It gives results identical to dense storage and is infrastructure for future block AMR, not a speedup.

• Several regions: pass a sequence of (x0, y0, x1, y1) rectangles as refinements[0], with both levels blocked (level_blocks). Transfers are built once from their union; adjacent, overlapping and concave unions are accepted when their actual dependencies close. There is no fixed rectangle-gap rule. The fine level is then one level-wide grid (2*ny, 2*nx): run.hierarchy.pairs[0].fine_shape gives it, initialize/download use it, and only the grid-aligned tiles intersecting the union plus its four-cell ghost band are stored. Block sizes must divide the level-wide grid; rectangle edges need not align with tiles. Order of rectangles does not matter. L1 bodies/Jets (level_objects[1]/level_jets[1]) work as for one patch: each circle (same coarse cell-edge frame) must lie wholly inside the union's active region with its boundary stencil clear of transfer cells; a jet may not select links of two different bodies, even across regions. Wall-touching patches remain unsupported. For nesting, add further refinements entries as shown below.

  cfg = LBMConfig(nx=32, ny=16, viscosity=.1)
  rest = np.asarray(compute_feq_d2q9(1.0, 0.0, 0.0), np.float32)[:, None, None]
  with CudaContext(0) as ctx, FixedPatchStepper(
          cfg, context=ctx, refinements=[((4, 4, 12, 12), (16, 4, 24, 12))],
          level_blocks=[(4, 4), (4, 4)]) as run:
      run.initialize(np.broadcast_to(rest, (9, *run.hierarchy.pairs[0].coarse_shape)),
                     np.broadcast_to(rest, (9, *run.hierarchy.pairs[0].fine_shape)))  # (9, 32, 64)
      run.step(32)
      coarse_g, fine_g = run.download()  # coarse: NaN on every covered cell;
                                         # fine: NaN on every non-ACTIVE cell
  

• Three-level nested windows: pass two scalar rectangles in refinements — L1 in L0 cells, L2 in L1 window cells (the parent's fine grid including its 4-cell ghost band) — plus all three level_blocks entries. L2 must sit inside L1's ACTIVE cells, away from L1's ghost band and transfer footprint (checked before any GPU allocation). initialize and download take/return three arrays; each level is NaN on its non-ACTIVE (and child-covered) cells. One coarse step = 2 L1 steps = 4 L2 steps. One L1 window only in this scalar chain. Optional stationary circles and Jets live on L2 via level_objects[2] and level_jets[2]; L1 bodies/Jets retain their L1 meaning and must be empty in this mode. All object and Jet geometry uses the L0 edge frame, not the child bounds' L1 window frame. L2 has dx=dt=0.25, so prescribed lattice velocity is unchanged. Each whole body and its BFL stencil must be inside L2 ACTIVE cells, outside the transfer footprint. A Jet may select one physical body across multiple blocks, never two bodies.

  from CelerisLab.body.geometry.circle import CircleGeometry
  from CelerisLab.body.objects import SimObject
  from CelerisLab.body.jets import Jet, Tube, Parallel
  
  cfg = LBMConfig(nx=16, ny=16, viscosity=.1)
  body = SimObject(0, geometry=CircleGeometry(8., 8., .8), center=(8., 8.), radius=.8)
  jet = Jet(Tube((6.5, 8.), 0., 3., .3), Parallel(0.), amplitude=.002)
  rest = np.asarray(compute_feq_d2q9(1.0, 0.0, 0.0), np.float32)[:, None, None]
  with CudaContext(0) as ctx, FixedPatchStepper(
          cfg, context=ctx, refinements=[(4, 4, 12, 12), (8, 8, 16, 16)],
          level_blocks=[(4, 4)] * 3, level_objects=[(), (), (body,)],
          level_jets=[(), (), (jet,)]) as run:
      shapes = (run.hierarchy.pairs[0].coarse_shape, run.hierarchy.pairs[0].fine_shape, run.hierarchy.pairs[1].fine_shape)
      run.initialize(*(np.broadcast_to(rest, (9, *s)) for s in shapes))  # 16², 24², 24²
      run.step(4)
      run.set_jet_amplitude(0, -.002, level=2)  # explicit level selects the driven level
      run.step(4)
      coarse_g, l1_g, l2_g = run.download()  # solids/covered/nonleaves are NaN
  

• Level lists (the only constructor): refinements[l] describes level l+1 in allocated level-l integer cells; use a rectangle list for a union, or a scalar rectangle for a single window (a scalar child requires a scalar parent). level_blocks, level_objects and level_jets each have one entry per level, including L0. Three or more levels require all-blocked storage and bodies/Jets only on the finest level. Every adjacent pair must preserve inherited ACTIVE coverage, valid donors/receivers, surviving pulls and the upstream transfer footprint. Level l takes 2**l steps per base interval. This describes implemented topology and scheduling, not accuracy or convergence qualification.

• Blocked body levels can enable observe_raw_exchange=True and call read_raw_exchange(body_id, level=...) after a completed base interval. It reduces unique cut-link observations into wall/aperture/total force, torque and mass exchange in base-lattice units per unit depth, for the last base interval (default) or window="last_substep". observe_window_exchange=True, alongside observe_raw_exchange=True, enables start_exchange_window() / read_exchange_window(body_id, level=...) across completed base intervals. These are raw boundary-reaction diagnostics, not a general stress or nozzle-load accuracy claim; dense body layouts are unsupported.

• Not available: checkpoints, moving bodies, L1 bodies/Jets in a three-level chain, automatic dynamic refinement, Simulation integration. Qualification depends on the tested configuration; see the validation scope. CPU topology/schedule tests alone do not establish GPU numerical correctness, general convergence or long-time stability.

Installation

Prerequisites

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

Install from source

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 (v0.5.1, confirmed bit-identical to double-buffer):

• 2D D2Q9 only (D3Q19 not yet implemented)
• MRT collision with both regularized and zou_he_local inlets
• Fixed and rotating cylinder benchmarks:

Benchmark	D	EsoPull CD	Double-buffer CD	CD diff
Kan99b K2	20	1.101	1.137	<3.2%
Kan99b K2	30	1.082	1.146	<5.6%

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

Known limitation (v0.5.1): FP16S combined with Bouzidi curved boundaries produces ~30-40% CD error even with ddf_shifting enabled. This is inherent to Bouzidi's per-direction DDF reads — FP16 quantization noise (~3e-5 per value) is not averaged across directions. DDF shifting is essential for FP16S (keeps values near 0 where FP16 has best precision), but does not fully resolve the Bouzidi incompatibility.

Verified benchmarks:

• Sah04 S2: St error within 1.5% (channel + curved + inlet/outlet)
• Kan99b K2: ~30% CD error with ddf_shifting (FP16S quantization noise through Bouzidi)
• Kan99b K2 without ddf_shifting: >100% error (unusable)
• High-blockage cases (S4 beta=0.9): May diverge earlier than FP32

For force-critical applications with curved boundaries, use FP32 storage. FP16S is suitable for applications where the curved boundary is simple (large D) or force accuracy is secondary.

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	Verified (K2 metrics match FP32)
MRT	esopull	zou_he_local	cylinder	Verified (sync_bodies tested, bit-identical to double-buffer)
MRT	esopull	regularized	cylinder	Verified (bit-identical to double-buffer)
SRT	double_buffer	any	cylinder	Expected to work (f-feq style)

Known limitations (ddf_shifting):

• Verified configurations as of v0.5.1: D2Q9 + MRT + double_buffer + any_inlet and D2Q9 + MRT + esopull + any_inlet (sync_bodies add, remove, stepping stable)
• 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")

Validation campaigns

• Periodic Taylor–Green AMR: analytic periodic-decay pilot and frozen baseline-relative assessment.
• Matched refined DFG: frozen SRT uniform/refined comparison, requiring channel GPU gates before execution.
• DFG 2D-2 Re100: frozen periodicity and official-reference comparisons, with a declared noncanonical outlet.
• Uniform LES shear decay: independent continuum-oracle flow diagnostic, baseline and isolated corrected source variants.
• Smagorinsky helper comparison: pinned baseline versus an output-local two-coefficient candidate, without production edits.
• Jet-circle control volume: bounded AMR versus uniform-grid conservation diagnostic.

The runners are in tests/validation, their CPU tests in tests/unit, and compact reference/source guards in tests/validation/data. Current Re100/Jet source guards explicitly pin the reviewed AMR-channel experiment; historical d79 evidence remains tied to its original checkpoint commits. New runs record actual Git HEAD separately from source origin. These campaigns preserve pinned production code; successful CPU checks alone do not establish GPU or physical validation. Keep generated results outside tracked files.

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.