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

```python
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()
```

### Multi-body control loop

```python
from CelerisLab import Simulation

sim = Simulation()
# Three rotating cylinders
sim.add_body("circle", center=(1006, 150), radius=10)
sim.add_body("circle", center=(1015, 140), radius=10)
sim.add_body("circle", center=(1015, 160), radius=10)
# Downstream velocity sensor
sim.add_body("sensor", center=(1050, 150), radius=10)
sim.initialize()

for step in range(100):
    # Set body rotation speeds (implicit GPU upload)
    sim.set_body(0, omega=0.002)
    sim.set_body(1, omega=-0.001)
    sim.set_body(2, omega=0.001)

    # Advance 10 LBM steps
    sim.run(10)

    # Read telemetry
    fx, fy = sim.read_force(0)
    ux, uy = sim.read_sensor(3)
    print(f"force=({fx:.4f},{fy:.4f}) sensor=({ux:.4f},{uy:.4f})")

sim.close()
```

## 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

```python
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 |
| `sim.add_body(type="sensor", center=(x,y), radius=r)` | int body_id | Add a velocity sensor |
| `sim.add_cylinder(center, radius)` | int body_id | Backward-compat alias |
| `sim.add_sensor(center, radius)` | int body_id | Backward-compat alias |
| `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, checkpoint_interval=0)` | Run N LBM steps |
| `sim.set_body(id, omega=...)` | Set body rotation speed (implicit GPU upload, ~1 μs) |
| `sim.read_force(id)` -> ndarray | Force vector [fx, fy] (2D) |
| `sim.read_torque(id)` -> ndarray | Torque [tz] (2D) |
| `sim.read_sensor(id)` -> ndarray | Area-averaged velocity via GPU sensor kernel |
| `sim.set_force(id, fx=..., fy=...)` | Set force density on a force_region object (notice: see persistence note below) |

### force_region usage

```python
# 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()` writes the action buffer directly but does not update the object's state record. If `sync_to_gpu()` is called afterward, the force will be reset to zero. For the common usage pattern (initialize -> set_force -> run -> set_force -> run ...), this is not an issue. A future update will add proper force storage in the object state.

### 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 | None needed |
| `"force_region"` | None (zero mask) | **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)

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

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

```python
stream = cuda.Stream()
sim.bodies.zero_force_segment_async(stream)
sim.stepper.step(
    1,
    action_gpu=sim.bodies.action_gpu,
    obs_gpu=sim.bodies.obs_gpu,
    stream=stream,
)
stream.synchronize()
force = sim.read_force(0)
```

## 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

```json
{
  "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 | MLUPS |
|--------|-------|
| Re100 MRT noLES | ~4400 |

### 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](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:

```bash
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

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

## Project Layout

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

```python
sim.set_body(0, omega=0.002)
sim.run(10)
force = sim.read_force(0)
```

### Async control (performance-oriented)

```python
sim.set_body(0, omega=0.002)  # implicit H2D, ~1 μs
sim.stepper.step(10, ..., stream=sim.stream)
sim.bodies.download_obs_full_async(sim.stream)
sim.stream.synchronize()
force = sim.read_force(0)
```

## Vortex initialization

```python
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

```python
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

```bibtex
@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.