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

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

### DRL control loop

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

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

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

```python
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.patch.coarse_shape)),
                   np.broadcast_to(rest, (9, *run.patch.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.patch.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.

  ```python
  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.patch.coarse_shape)),
                     np.broadcast_to(rest, (9, *run.patch.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.

  ```python
  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.patch.coarse_shape, run.patch.fine_shape, run.child_patch.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)  # default level=1 still addresses L1
      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](docs/validation/2d_scope.md).
  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

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

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

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

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

```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 | 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](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

See [docs/tests_overview.md](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)

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

### Async control (performance-oriented)

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

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

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

## Validation campaigns

- [Periodic Taylor–Green AMR](docs/validation/taylor_green.md): analytic periodic-decay pilot and frozen baseline-relative assessment.
- [Matched refined DFG](docs/validation/refined_dfg.md): frozen SRT uniform/refined comparison, requiring channel GPU gates before execution.
- [DFG 2D-2 Re100](docs/validation/dfg_re100.md): frozen periodicity and official-reference comparisons, with a declared noncanonical outlet.
- [Uniform LES shear decay](docs/validation/uniform_les_flow.md): independent continuum-oracle flow diagnostic, baseline and isolated corrected source variants.
- [Smagorinsky helper comparison](docs/validation/les_helpers.md): pinned baseline versus an output-local two-coefficient candidate, without production edits.
- [Jet-circle control volume](docs/validation/jet_circle_cv.md): 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

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