重构：分级设计，EsoPull

2026-04-17 21:50:38 +08:00 · 2026-04-17 21:50:38 +08:00 · 654e3290a4
commit 654e3290a4
parent b339f302d3
40 changed files with 2767 additions and 619 deletions
--- a/README.md
+++ b/README.md
@ -7,13 +7,14 @@ CelerisLab is a high-performance computational fluid dynamics (CFD) solver based
 ## Features

 - **GPU Acceleration**: CUDA-based kernels for high-performance simulations
- **D2Q9 Lattice**: 2D nine-velocity lattice implementation
- **MRT Collision Model**: Multiple-Relaxation-Time collision operator for improved stability
+- **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)
 - **Immersed Boundary Method (IBM)**: Support for complex geometries (cylinders, arbitrary shapes)
- **Flexible Boundary Conditions**: Periodic, velocity inlet, pressure outlet
- **Real-time Sensors**: Monitor flow properties at specific locations during simulation
- **Vortex Initialization**: Built-in support for Lamb, Oseen, and Taylor vortices
- **Dynamic Compilation**: Runtime CUDA kernel compilation with configurable parameters
+- **Flexible Boundary Conditions**: NEQ-extrapolation pressure outlet, parabolic/uniform velocity inlet, half-way bounce-back walls
+- **Layered Configuration**: Compile-time parameters organized into Global / Method / Case / Debug tiers
+- **High-Re Validated**: Tested up to Re=5000 (2D cylinder); MRT+LES and SRT+LES stable; TRT+LES stable with tuned Lambda and WMAX
+- **Python API**: High-level `Simulation` class for scripting and RL integration

 ## Installation

@ -40,162 +41,149 @@ pip install -e .  # Installs from src/ directory

 ## Quick Start

-### Basic Flow Simulation
+```python
+from CelerisLab import Simulation
+
+sim = Simulation("configs/config_lbm.json")
+sim.add_cylinder(center=(50, 50), radius=10)
+sim.initialize()
+
+for step in range(10000):
+    sim.run(1)
+
+macro = sim.get_macroscopic()  # {"rho": ..., "ux": ..., "uy": ...}
+sim.close()
+```
+
+Or as a context manager:

 ```python
-from CelerisLab import FlowField, utils
-
-# Load configurations
-config_cuda = utils.load_cuda_config()  # Uses default or CELERISLAB_CONFIG_DIR
-config_field = utils.load_flow_field_config()
-
-# Initialize flow field
-flow = FlowField(
-    config_cuda=config_cuda,
-    config_field=config_field,
-    device_id=0
-)
-
-# Add a cylinder obstacle
-flow.add_cylinder(
-    center=(50, 50, 0),
-    radius=10,
-    velocity=(0, 0, 0),
-    use_IBM=True
-)
-
-# Add sensors to monitor flow
-flow.add_sensor(position=(70, 50, 0))
-
-# Run simulation
-for step in range(10000):
-    flow.run(1)
-    
-    # Read sensor data every 100 steps
-    if step % 100 == 0:
-        sensor_data = flow.read_sensor()
-        print(f"Step {step}: Velocity = {sensor_data[0]}")
+with Simulation("configs/config_lbm.json") as sim:
+    sim.add_cylinder(center=(96, 64), radius=12)
+    sim.initialize()
+    sim.run(5000)
+    data = sim.get_macroscopic()
 ```

-### Configuration
+## Configuration

-CelerisLab searches for configuration files in the following order:
+### `configs/config_lbm.json`

-1. **Explicit path**: Passed to `load_*_config(config_path)`
-2. **Environment variable**: `CELERISLAB_CONFIG_DIR` environment variable
-3. **Current directory**: `./configs/` in current working directory
-4. **Package default**: Bundled `CelerisLab/configs/` directory
-
-#### Configuration Files
-
-**config_cuda.json**: CUDA execution parameters
 ```json
 {
-  "multi_gpu": false,
-  "gpu_connection": "NVLINK",
-  "required_cuda_capability": "6.0",
-  "threads_per_block": 256,
-  "X_1U": 16,
-  "Y_1U": 16,
-  "Z_1U": 1
+  "dim": 2,
+  "nq": 9,
+  "nx": 384,
+  "ny": 192,
+  "nz": 1,
+  "viscosity": 0.0005,
+  "velocity": 0.04,
+  "rho": 1.0,
+  "collision": "MRT",
+  "streaming": "double_buffer",
+  "les_enabled": true,
+  "les_cs": 0.16,
+  "trt_magic_param": 0.001,
+  "omega_max": 1.90,
+  "inlet_profile": "parabolic",
+  "outlet_mode": "neq_extrap",
+  "compute_capability": "auto",
+  "threads_per_block": 256
 }
 ```

-**config_flowfield.json**: Flow physics parameters
-```json
-{
-  "data_type": "FP32",
-  "dimensionality": 2,
-  "lattice": 9,
-  "field_dim_in_U": [100, 100, 1],
-  "viscosity": 0.01,
-  "velocity": 0.1,
-  "boundary_conditions": {
-    "x": ["periodic", "periodic"],
-    "y": ["periodic", "periodic"],
-    "z": ["periodic", "periodic"]
-  }
-}
-```
+### Parameter tiers
+
+| Tier | Headers | Examples |
+|---|---|---|
+| Global/Grid | `config_grid.h` | DIM, NQ, NX, NY, NZ |
+| Global/Physics | `config_physics.h` | VIS, RHO, U0, flag constants |
+| Method | `config_method.h` | COLLISION_MODEL, USE_LES, TRT_MAGIC_PARAM, OMEGA_COLLISION_MAX |
+| Case | `config_objects.h` | N_OBJS |
+
+Headers are auto-generated by the compiler from `LBMConfig`; do not edit manually.

 ## API Reference

-### FlowField Class
-
-Main interface for running LBM simulations.
-
-#### Constructor
-```python
-FlowField(config_cuda, config_field, device_id=0)
-```
-
-#### Methods
-
- `add_cylinder(center, radius, velocity, use_IBM=False)`: Add cylindrical obstacle
- `add_sensor(position)`: Add flow monitoring sensor
- `add_vortex(center, circulation, core_radius, vortex_type='Lamb')`: Initialize vortex
- `run(n_steps)`: Execute simulation steps
- `read_sensor()`: Read current sensor values
- `get_ddf()`: Get distribution function data
- `apply_ddf(ddf)`: Set distribution function data
-
-### Utility Functions
-
- `load_cuda_config(config_path=None)`: Load CUDA configuration
- `load_flow_field_config(config_path=None)`: Load flow field configuration
- `check_cuda_device_availability(device_id=0)`: Verify CUDA device
- `get_device_info(device_id=0)`: Query GPU properties
- `estimate_memory_consumption(config_field, num_objects, num_sensors)`: Calculate memory usage
-
-## Advanced Usage
-
-### Custom Geometry with IBM
+### `Simulation`

 ```python
-# IBM enables smooth treatment of curved boundaries
-flow.add_cylinder(
-    center=(grid_x//2, grid_y//2, 0),
-    radius=20,
-    velocity=(0.0, 0.0, 0.0),
-    use_IBM=True  # Enables immersed boundary method
-)
+sim = Simulation(lbm_config_path=None, body_config_path=None, device_id=0)
+sim.add_cylinder(center, radius) -> int
+sim.add_sensor(center, radius)  -> int
+sim.initialize()
+sim.run(steps)
+sim.step(n=1)
+sim.get_macroscopic() -> {"rho": ndarray, "ux": ndarray, "uy": ndarray}
+sim.get_ddf()    -> ndarray
+sim.get_flags()  -> ndarray
+sim.update_runtime_params(omega=..., u_inlet=...)
+sim.snapshot() / sim.restore()
+sim.close()
 ```

-### Multiple Sensors
+### Vortex initialization

 ```python
-# Add sensors in a line downstream of obstacle
-for i in range(5):
-    flow.add_sensor(position=(100 + i*10, 50, 0))
+from CelerisLab.lbm.initializers import add_vortex

-# Read all sensors at once
-sensor_data = flow.read_sensor()  # Returns array of shape (n_sensors, 3)
+# Superimpose a Lamb–Oseen vortex on an existing LBMField
+add_vortex(sim.field, center=(50, 50), radius=10.0, strength=1.0, vortex_type="lamb")
 ```

-### Vortex Initialization
+## Collision & LES Recommendations

-```python
-# Initialize Lamb-Oseen vortex
-flow.add_vortex(
-    center=(50, 50, 0),
-    circulation=1.0,
-    core_radius=10.0,
-    vortex_type='Lamb'
-)
+| 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_max=1.90`, `trt_magic_param=0.001` |
+
+## Project Layout
+
+```
+src/CelerisLab/
+  simulation.py          High-level API
+  config.py              LBMConfig / BodyConfig dataclasses
+  cuda/
+    compiler_v2.py       Config header generation + nvcc + PTX load
+    context.py           CUDA context lifecycle
+  lbm/
+    field.py             GPU memory management
+    stepper.py           Time-step driver
+    initializers.py      Vortex superposition
+    kernels/
+      kernel_v2.cu       Kernel entry (thin wrapper)
+      config/            Auto-generated config headers
+      core/              Descriptors, layout, flags, params
+      operators/         Collision, LES, forcing
+      boundary/          Inlet, outlet, wall, curved, IBM
+      streaming/         Double-buffer & esopull
+      step/              Step orchestration
+  body/
+    objects.py           SimObject / Cylinder / Sensor
+    manager.py           ObjectManager + GPU sync
+  common/
+    preprocess.py        Geometry utilities
+tests/
+  test_stability_matrix.py    13-case stability matrix (Re × collision × LES × streaming)
+  test_high_re_validation.py  High-Re directed validation (Re5000, 2D/3D, parameter sweep)
+output/
+  CelerisLab_stage1_architecture.md   Architecture specification (v3)
+  refactor_brief_stage1.md            Refactoring brief
+  high_re_audit_round1.md             8-round audit log
+legacy/                               Superseded code (FlowField, compiler v1, macros.h)
 ```

-## Environment Variables
+## Performance

- `CELERISLAB_CONFIG_DIR`: Directory containing configuration JSON files
- `OMP_NUM_THREADS`: OpenMP thread count (recommend setting to 1 for GPU workflows)
- `MKL_NUM_THREADS`: Intel MKL thread count (recommend setting to 1)
+Tested on Tesla V100-SXM2-16GB (CUDA 12.4):

-## Performance Tips
-
-1. **Grid Size**: Choose dimensions that are multiples of `unit_dimensions` in config_cuda.json
-2. **Thread Block Size**: 256 threads/block works well for most GPUs
-3. **Memory**: Estimate memory with `utils.estimate_memory_consumption()` before large runs
-4. **Single-threaded Python**: Set `OMP_NUM_THREADS=1` to avoid CPU interference with GPU
+| Config | Grid | MLUPS |
+|---|---|---|
+| Re100 MRT noLES | 384×192 | ~4200 |
+| Re100 EsoPull SRT | 384×192 | ~3900 |
+| Re3000 MRT+LES | 384×192 | ~4360 |

 ## Citation

@ -212,13 +200,4 @@ If you use CelerisLab in your research, please cite:

 ## License

-MIT License - see LICENSE file for details
-
-## Contributing
-
-Contributions are welcome! Please feel free to submit issues and pull requests.
-
-## Acknowledgments
-
- Built with PyCUDA by Andreas Klöckner
- Inspired by the palabos C++ LBM library
+MIT License — see LICENSE file for details.
--- a/legacy/README.md
+++ b/legacy/README.md
@ -0,0 +1,19 @@
+# Legacy Code Archive
+
+This directory contains code that has been superseded by the current architecture but is kept for reference.
+
+## Contents
+
+| File / Dir | Replaced By | Reason |
+|---|---|---|
+| `lbm_driver.py` | `src/CelerisLab/simulation.py` + `lbm/field.py` + `lbm/stepper.py` | Monolithic FlowField class. New Simulation API separates concerns: CudaContext / LBMField / LBMStepper / ObjectManager. |
+| `cuda_compiler_v1.py` | `src/CelerisLab/cuda/compiler_v2.py` | macros.h-based build system. New compiler writes typed config/*.h headers per architectural layer. |
+| `macros.h` | `src/CelerisLab/lbm/kernels/config/*.h` | Single flat macro file. Now split into config_grid.h / config_physics.h / config_method.h / config_objects.h matching the Global/Method/Case/Debug parameter hierarchy. |
+| `common_utils.py` | `src/CelerisLab/config.py` + `src/CelerisLab/cuda/context.py` | FlowFieldConfig / CudaConfig NamedTuples and their JSON loaders. Replaced by LBMConfig / BodyConfig dataclasses (config.py) and CudaContext (cuda/context.py). |
+| `lbm_configs/` | `src/CelerisLab/configs/` | Old JSON config format used by FlowField / compiler_v1. |
+
+## Notes
+
+- None of these files is imported by any active module.
+- `lbm_driver.py` (FlowField) depended on `cuda_compiler_v1.py` and `common_utils.py`; all three were removed from src together.
+- `macros.h` was the old single-file configuration for `kernel_v2.cu`; kernel_v2.cu now includes `config.h` which aggregates `config/*.h`.
--- a/src/CelerisLab/common/utils.py
+++ b/src/CelerisLab/common/utils.py
--- a/src/CelerisLab/cuda/compiler.py
+++ b/src/CelerisLab/cuda/compiler.py
--- a/src/CelerisLab/lbm/configs/config_cuda.json
+++ b/src/CelerisLab/lbm/configs/config_cuda.json
--- a/src/CelerisLab/lbm/configs/config_flowfield.json
+++ b/src/CelerisLab/lbm/configs/config_flowfield.json
--- a/src/CelerisLab/lbm/driver.py
+++ b/src/CelerisLab/lbm/driver.py
--- a/src/CelerisLab/lbm/kernels/macros.h
+++ b/src/CelerisLab/lbm/kernels/macros.h
@ -3,21 +3,21 @@
 // cuda parameters
 #define MULT_GPU False
 #define NT 128
-#define X_1U 256
-#define Y_1U 128
-#define Z_1U 32
+#define X_1U 384
+#define Y_1U 192
+#define Z_1U 1

 // flow parameters
 #define LBtype float
 #define UX 1
 #define UY 1
 #define UZ 1
-#define NX 256
-#define NY 128
-#define NZ 32
-#define DIM 3
-#define NQ 19
-#define VIS 0.0096000000
+#define NX 384
+#define NY 192
+#define NZ 1
+#define DIM 2
+#define NQ 9
+#define VIS 0.0144000000
 #define RHO 1.0
 #define U0 0.04

@ -69,7 +69,7 @@

 // Smagorinsky constant C_s
 #ifndef LES_CS
-#define LES_CS 0.16f
+#define LES_CS 0.160000f
 #endif

 // Inlet profile: 1=parabolic (channel), 0=uniform (external flow)
@ -85,7 +85,7 @@
 // Outlet blend factor for damped outlet mode (OUTLET_MODE=2):
 // f_out = a*(non-eq extrapolation) + (1-a)*(zero-gradient copy)
 #ifndef OUTLET_BLEND_ALPHA
-#define OUTLET_BLEND_ALPHA 0.70f
+#define OUTLET_BLEND_ALPHA 0.700f
 #endif

 // Outlet backflow clamp: 0=off, 1=force non-negative streamwise velocity at outlet target
@ -104,5 +104,5 @@

 // TRT magic parameter Lambda used to map omega+ -> omega-
 #ifndef TRT_MAGIC_PARAM
-#define TRT_MAGIC_PARAM 0.1875f
+#define TRT_MAGIC_PARAM 0.187500f
 #endif
--- a/src/CelerisLab/init.py
+++ b/src/CelerisLab/init.py
@ -2,43 +2,28 @@
 """
 CelerisLab: GPU-Accelerated Computational Physics Simulation Library

-A modular framework for GPU-accelerated physics simulations including:
- Lattice Boltzmann Method (LBM) for fluid dynamics
- Future: Finite Element Method (FEM), Smoothed Particle Hydrodynamics (SPH), etc.
+Usage::

-Usage:
-    # Direct import of main classes
-    from CelerisLab import FlowField
-    
-    # Module-specific imports
-    from CelerisLab.lbm import FlowField
-    
-    # Namespace style
-    import CelerisLab as cl
-    solver = cl.lbm.FlowField(...)
+    from CelerisLab import Simulation
+
+    sim = Simulation("configs/config_lbm.json")
+    sim.add_cylinder((100, 50), radius=10)
+    sim.initialize()
+    sim.run(1000)
+    data = sim.get_macroscopic()
+
+Legacy FlowField API is still importable from CelerisLab.lbm.driver.
 """

-__version__ = '0.2.0'
+__version__ = "0.3.0"

-# Import submodules
-from . import common
-from . import cuda
-from . import lbm
+from . import common, cuda, lbm, body, config

-# Import commonly used utilities for convenience
-from .common import utils
-
-# Attempt to import main classes for direct access
 try:
-    from .lbm import FlowField
-    __all__ = ['lbm', 'common', 'cuda', 'utils', 'FlowField', '__version__']
+    from .simulation import Simulation
+    __all__ = ["Simulation", "lbm", "body", "common", "cuda", "config",
+               "__version__"]
 except ImportError as e:
-    # PyCUDA not available, only submodules accessible
    import warnings
-    warnings.warn(
-        f"FlowField not available: {e}. "
-        "Install pycuda to use the full CelerisLab functionality. "
-        "Utils and other modules are still accessible.",
-        ImportWarning
-    )
-    __all__ = ['lbm', 'common', 'cuda', 'utils', '__version__']
+    warnings.warn(f"Simulation not available: {e}", ImportWarning)
+    __all__ = ["lbm", "body", "common", "cuda", "config", "__version__"]
--- a/src/CelerisLab/body/init.py
+++ b/src/CelerisLab/body/init.py
@ -0,0 +1,8 @@
+# CelerisLab/body/__init__.py
+"""
+Body / object management for immersed and rigid objects.
+"""
+from .objects import SimObject, Cylinder, Sensor
+from .manager import ObjectManager
+
+__all__ = ["SimObject", "Cylinder", "Sensor", "ObjectManager"]
--- a/src/CelerisLab/body/manager.py
+++ b/src/CelerisLab/body/manager.py
@ -0,0 +1,167 @@
+# CelerisLab/body/manager.py
+"""
+ObjectManager — batch management of SimObjects.
+
+Responsibilities:
+  - Add / remove / query objects
+  - Build merged flag, indx, delta arrays from all objects
+  - Allocate action / obs GPU buffers
+  - Sync geometry to GPU
+  - (future) detect collisions, exchange forces, update states
+"""
+
+import numpy as np
+import pycuda.driver as cuda
+from typing import Dict, List, Optional
+
+from .objects import SimObject, FLUID, OBSTACLE, SENSOR_FLAG, INTERFACE, SOLID
+
+
+class ObjectManager:
+    """Central registry for all simulation objects."""
+
+    def __init__(self, nx: int, ny: int, nq: int, dim: int = 2):
+        self.nx = nx
+        self.ny = ny
+        self.nq = nq
+        self.dim = dim
+        self._objects: Dict[int, SimObject] = {}
+        self._next_id = 0
+
+        # GPU buffers (allocated on first sync)
+        self.action_gpu: Optional[cuda.DeviceAllocation] = None
+        self.obs_gpu: Optional[cuda.DeviceAllocation] = None
+
+        # Host buffers
+        self.action = np.zeros(0, dtype=np.float32)
+        self.obs = np.zeros(0, dtype=np.float32)
+
+    # -- Object CRUD --------------------------------------------------------
+    def add(self, obj: SimObject) -> int:
+        """Register an object and return its id."""
+        obj.obj_id = self._next_id
+        self._objects[self._next_id] = obj
+        self._next_id += 1
+        self._resize_buffers()
+        return obj.obj_id
+
+    def remove(self, obj_id: int):
+        del self._objects[obj_id]
+        self._resize_buffers()
+
+    def get(self, obj_id: int) -> SimObject:
+        return self._objects[obj_id]
+
+    @property
+    def objects(self) -> List[SimObject]:
+        return list(self._objects.values())
+
+    @property
+    def count(self) -> int:
+        return len(self._objects)
+
+    # -- Buffer management ---------------------------------------------------
+    def _resize_buffers(self):
+        n = self.count
+        self.action = np.zeros(max(n, 1), dtype=np.float32)
+        self.obs = np.zeros(max(n * self.dim, 1), dtype=np.float32)
+
+    # -- Build merged arrays -------------------------------------------------
+    def build_flags(self, base_flags: np.ndarray) -> np.ndarray:
+        """Merge all object flag masks onto *base_flags* (modified in-place)."""
+        for obj in self._objects.values():
+            mask = obj.get_flag_mask(self.nx, self.ny)
+            # Set obstacle/sensor bits; also mark solid+interface for curved BC
+            for k in range(mask.size):
+                if mask[k] != 0:
+                    if mask[k] == OBSTACLE:
+                        base_flags[k] = SOLID | INTERFACE | OBSTACLE
+                    else:
+                        base_flags[k] |= mask[k]
+        return base_flags
+
+    def build_indx(self, base_indx: np.ndarray) -> np.ndarray:
+        """Merge per-cell object indices."""
+        for obj in self._objects.values():
+            indx = obj.get_indx_map(self.nx, self.ny)
+            nonzero = indx != 0
+            base_indx[nonzero] = indx[nonzero]
+        return base_indx
+
+    def build_delta(self) -> np.ndarray:
+        """Concatenate curved-boundary data from all objects."""
+        parts = []
+        for obj in self._objects.values():
+            d = obj.get_delta_curve(self.nx, self.ny, self.nq)
+            if d.size > 0:
+                parts.append(d)
+        if parts:
+            return np.concatenate(parts)
+        return np.zeros(0, dtype=np.float32)
+
+    # -- GPU sync ------------------------------------------------------------
+    def sync_to_gpu(self, field):
+        """Upload merged flags, indx, delta, and action/obs buffers."""
+        field.flag = self.build_flags(field.flag)
+        field.indx = self.build_indx(field.indx)
+        field.delta = self.build_delta()
+
+        field.upload_flags()
+        field.upload_delta()
+        field.update_params(n_objects=self.count)
+
+        # Alloc action/obs GPU
+        self._free_gpu()
+        self.action_gpu = cuda.mem_alloc(self.action.nbytes)
+        self.obs_gpu = cuda.mem_alloc(self.obs.nbytes)
+        cuda.memcpy_htod(self.action_gpu, self.action)
+        cuda.memcpy_htod(self.obs_gpu, self.obs)
+
+        # Impose rest equilibrium on newly non-fluid nodes
+        self._rest_nonfluid(field)
+
+    def _rest_nonfluid(self, field):
+        """Set DDF of non-fluid nodes to rest equilibrium w_i * rho0."""
+        field.download_ddf()
+        nq = field.nq
+        nx, ny = field.nx, field.ny
+        n = nx * ny
+        # D2Q9 weights
+        if nq == 9:
+            w = np.array(
+                [4/9,1/9,1/9,1/9,1/9,1/36,1/36,1/36,1/36],
+                dtype=np.float32,
+            )
+        elif nq == 19:
+            w = np.array(
+                [1/3]+[1/18]*6+[1/36]*12, dtype=np.float32,
+            )
+        else:
+            return
+        f = field.ddf.reshape(nq, -1)
+        nonfluid = (field.flag & FLUID) == 0
+        for i in range(nq):
+            f[i, nonfluid] = w[i] * field.cfg.rho
+        field.ddf = f.reshape(-1)
+        field.upload_ddf()
+
+    def _free_gpu(self):
+        if self.action_gpu is not None:
+            self.action_gpu.free()
+            self.action_gpu = None
+        if self.obs_gpu is not None:
+            self.obs_gpu.free()
+            self.obs_gpu = None
+
+    # -- Future placeholders -------------------------------------------------
+    def update_states(self, dt: float):
+        """Integrate object motion (placeholder)."""
+        pass
+
+    def detect_collisions(self):
+        """Multi-body collision detection (placeholder)."""
+        pass
+
+    def exchange_forces(self, obs: np.ndarray):
+        """Distribute fluid forces to objects (placeholder)."""
+        pass
--- a/src/CelerisLab/body/objects.py
+++ b/src/CelerisLab/body/objects.py
@ -0,0 +1,189 @@
+# CelerisLab/body/objects.py
+"""
+Lightweight flat object model for immersed / rigid bodies and sensors.
+
+Design:
+  - SimObject is a thin base with standard interface (state, control,
+    get_flag_mask, get_delta_curve).
+  - Concrete types (Cylinder, Sensor, …) override geometry methods.
+  - No deep inheritance tree.  Users can subclass SimObject directly
+    for custom shapes — just implement get_flag_mask / get_delta_curve.
+"""
+
+import math
+import numpy as np
+from dataclasses import dataclass, field
+from typing import Optional, Tuple
+
+# Reuse flag constants (must match config/config_physics.h)
+FLUID = 0x01
+SOLID = 0x02
+OBSTACLE = 0x20
+INTERFACE = 0x08
+SENSOR_FLAG = 0x10
+
+
+# ---------------------------------------------------------------------------
+# State / Control containers
+# ---------------------------------------------------------------------------
+@dataclass
+class ObjectState:
+    """Position, velocity, orientation (2D for now)."""
+    x: float = 0.0
+    y: float = 0.0
+    theta: float = 0.0
+    vx: float = 0.0
+    vy: float = 0.0
+    omega: float = 0.0
+
+
+@dataclass
+class ObjectControl:
+    """Control input (force / velocity / displacement target)."""
+    ax: float = 0.0
+    ay: float = 0.0
+    alpha: float = 0.0
+    mode: int = 0  # 0=force, 1=velocity, 2=displacement
+
+
+# ---------------------------------------------------------------------------
+# Base class
+# ---------------------------------------------------------------------------
+class SimObject:
+    """Base for all simulation objects."""
+
+    obj_type: str = "generic"
+
+    def __init__(self, obj_id: int, center: Tuple[float, ...],
+                 radius: float = 0.0):
+        self.obj_id = obj_id
+        self.center = center
+        self.radius = radius
+        self.state = ObjectState(x=center[0], y=center[1])
+        self.control = ObjectControl()
+
+    def get_flag_mask(self, nx: int, ny: int) -> np.ndarray:
+        """Return (n,) uint8 array with flag bits set for this object."""
+        raise NotImplementedError
+
+    def get_delta_curve(self, nx: int, ny: int, nq: int) -> np.ndarray:
+        """Return flat float32 array of curved-boundary interpolation data."""
+        return np.zeros(0, dtype=np.float32)
+
+    def get_indx_map(self, nx: int, ny: int) -> np.ndarray:
+        """Return (n,) int32 with per-cell object index or offset."""
+        return np.zeros(nx * ny, dtype=np.int32)
+
+
+# ---------------------------------------------------------------------------
+# Concrete types
+# ---------------------------------------------------------------------------
+class Cylinder(SimObject):
+    obj_type = "cylinder"
+
+    def get_flag_mask(self, nx: int, ny: int) -> np.ndarray:
+        n = nx * ny
+        mask = np.zeros(n, dtype=np.uint8)
+        xc, yc = self.state.x, self.state.y
+        r = self.radius
+        x0 = max(1, int(xc - r) - 1)
+        x1 = min(nx - 2, int(xc + r) + 1)
+        y0 = max(1, int(yc - r) - 1)
+        y1 = min(ny - 2, int(yc + r) + 1)
+        for x in range(x0, x1 + 1):
+            for y in range(y0, y1 + 1):
+                if (x - xc)**2 + (y - yc)**2 < r**2:
+                    k = x + y * nx
+                    mask[k] = OBSTACLE
+        return mask
+
+    def get_delta_curve(self, nx: int, ny: int, nq: int) -> np.ndarray:
+        """Compute curved-boundary interpolation coefficients.
+
+        Uses ray–circle intersection (Yu-Mei-Shyy 2003 scheme).
+        Returns packed array: [obj_id, delta_1..delta_NQ, Uw_x, Uw_y] per node.
+        """
+        from ..common.preprocess import find_circle_intersection
+        xc, yc = self.state.x, self.state.y
+        r = self.radius
+
+        # D2Q9 velocity vectors
+        ex = [0, 1, -1, 0, 0, 1, -1, 1, -1]
+        ey = [0, 0, 0, 1, -1, 1, -1, -1, 1]
+
+        entries = []
+        x0 = max(1, int(xc - r) - 2)
+        x1 = min(nx - 2, int(xc + r) + 2)
+        y0 = max(1, int(yc - r) - 2)
+        y1 = min(ny - 2, int(yc + r) + 2)
+
+        for x in range(x0, x1 + 1):
+            for y in range(y0, y1 + 1):
+                if (x - xc)**2 + (y - yc)**2 >= r**2:
+                    continue
+                # This node is inside the obstacle
+                deltas = np.zeros(nq, dtype=np.float32)
+                has_fluid_nb = False
+                for i in range(nq):
+                    xn, yn = x + ex[i], y + ey[i]
+                    if 0 <= xn < nx and 0 <= yn < ny:
+                        if (xn - xc)**2 + (yn - yc)**2 >= r**2:
+                            hit = find_circle_intersection(
+                                float(x), float(y),
+                                float(xn), float(yn),
+                                xc, yc, r,
+                            )
+                            if hit is not None:
+                                dx = xn - x
+                                dy = yn - y
+                                dist = math.sqrt(
+                                    (hit[0]-x)**2 + (hit[1]-y)**2
+                                )
+                                norm = math.sqrt(dx*dx + dy*dy)
+                                deltas[i] = dist / norm if norm > 0 else 0.5
+                                has_fluid_nb = True
+                if has_fluid_nb:
+                    # Pack: [obj_id_as_float, delta_0..delta_NQ-1, Uw_x, Uw_y]
+                    entry = np.zeros(1 + nq + 2, dtype=np.float32)
+                    entry[0] = np.float32(self.obj_id).view(np.float32)
+                    entry[1:1+nq] = deltas
+                    entry[1+nq] = 0.0   # Uw_x (set by action at runtime)
+                    entry[1+nq+1] = 0.0 # Uw_y
+                    entries.append(entry)
+
+        if entries:
+            return np.concatenate(entries)
+        return np.zeros(0, dtype=np.float32)
+
+
+class Sensor(SimObject):
+    obj_type = "sensor"
+
+    def get_flag_mask(self, nx: int, ny: int) -> np.ndarray:
+        n = nx * ny
+        mask = np.zeros(n, dtype=np.uint8)
+        xc, yc = self.state.x, self.state.y
+        r = self.radius
+        x0 = max(1, int(xc - r) - 1)
+        x1 = min(nx - 2, int(xc + r) + 1)
+        y0 = max(1, int(yc - r) - 1)
+        y1 = min(ny - 2, int(yc + r) + 1)
+        for x in range(x0, x1 + 1):
+            for y in range(y0, y1 + 1):
+                if (x - xc)**2 + (y - yc)**2 < r**2:
+                    mask[x + y * nx] = SENSOR_FLAG
+        return mask
+
+    def get_indx_map(self, nx: int, ny: int) -> np.ndarray:
+        indx = np.zeros(nx * ny, dtype=np.int32)
+        xc, yc = self.state.x, self.state.y
+        r = self.radius
+        x0 = max(1, int(xc - r) - 1)
+        x1 = min(nx - 2, int(xc + r) + 1)
+        y0 = max(1, int(yc - r) - 1)
+        y1 = min(ny - 2, int(yc + r) + 1)
+        for x in range(x0, x1 + 1):
+            for y in range(y0, y1 + 1):
+                if (x - xc)**2 + (y - yc)**2 < r**2:
+                    indx[x + y * nx] = self.obj_id
+        return indx
--- a/src/CelerisLab/common/init.py
+++ b/src/CelerisLab/common/init.py
@ -3,7 +3,6 @@
 Common utilities and preprocessing functions.
 """

-from . import utils
 from . import preprocess

-__all__ = ['utils', 'preprocess']
+__all__ = ['preprocess']
--- a/src/CelerisLab/config.py
+++ b/src/CelerisLab/config.py
@ -0,0 +1,176 @@
+# CelerisLab/config.py
+"""
+Unified configuration system.
+
+Two JSON files:
+  config_lbm.json  – grid, physics, method, cuda
+  config_body.json – object list
+
+LBMConfig / BodyConfig are plain dataclasses; they translate to
+macro dicts consumed by compiler.py for header generation.
+"""
+
+import json
+import os
+from dataclasses import dataclass, field
+from typing import Any, Dict, List, Optional
+
+# ---------------------------------------------------------------------------
+# String ↔ integer mappings (used in JSON and macro generation)
+# ---------------------------------------------------------------------------
+COLLISION_MAP = {"SRT": 0, "TRT": 1, "MRT": 2}
+STREAMING_MAP = {"double_buffer": 0, "esopull": 1}
+PRECISION_MAP = {"FP32": 0, "FP16S": 1, "FP16C": 2}
+INLET_MAP = {"uniform": 0, "parabolic": 1}
+OUTLET_MAP = {"neq_extrap": 0, "zero_gradient": 1, "blended": 2}
+DTYPE_MAP = {"FP32": "float", "FP64": "double"}
+
+
+# ---------------------------------------------------------------------------
+# LBM config
+# ---------------------------------------------------------------------------
+@dataclass
+class LBMConfig:
+    # Grid
+    dim: int = 2
+    nq: int = 9
+    nx: int = 512
+    ny: int = 256
+    nz: int = 1
+
+    # Physics
+    data_type: str = "FP32"
+    viscosity: float = 0.002
+    velocity: float = 0.03
+    rho: float = 1.0
+
+    # Method – collision / streaming / storage
+    collision: str = "SRT"
+    streaming: str = "double_buffer"
+    store_precision: str = "FP32"
+    ddf_shifting: bool = False
+    les_enabled: bool = False
+    les_cs: float = 0.16
+    trt_magic_param: float = 0.1875
+    inlet_profile: str = "parabolic"
+    outlet_mode: str = "neq_extrap"
+    outlet_blend_alpha: float = 0.7
+    outlet_backflow_clamp: bool = True
+    omega_min: float = 0.01
+    omega_max: float = 1.999
+
+    # CUDA
+    threads_per_block: int = 128
+    compute_capability: str = "auto"
+
+    # Derived (computed on validate)
+    omega: float = 0.0
+
+    def validate(self):
+        assert self.dim in (2, 3), f"dim must be 2 or 3, got {self.dim}"
+        assert self.nq in (9, 19), f"nq must be 9 or 19, got {self.nq}"
+        assert self.collision in COLLISION_MAP, f"Unknown collision: {self.collision}"
+        assert self.streaming in STREAMING_MAP, f"Unknown streaming: {self.streaming}"
+        assert self.data_type in DTYPE_MAP, f"Unknown data_type: {self.data_type}"
+        assert self.nx > 0 and self.ny > 0 and self.nz > 0
+        if self.dim == 2:
+            assert self.nz == 1, "nz must be 1 for 2D"
+            assert self.nq == 9, "D2Q9 required for dim=2"
+        else:
+            assert self.nq == 19, "D3Q19 required for dim=3"
+        self.omega = 1.0 / (3.0 * self.viscosity + 0.5)
+
+    def to_macros(self) -> Dict[str, Any]:
+        """Return flat dict of macro_name → value for header generation."""
+        self.validate()
+        return {
+            "NT": self.threads_per_block,
+            "MULT_GPU": 0,
+            "NX": self.nx,
+            "NY": self.ny,
+            "NZ": self.nz,
+            "DIM": self.dim,
+            "NQ": self.nq,
+            "LBtype": DTYPE_MAP[self.data_type],
+            "VIS": f"{self.viscosity:.10f}",
+            "RHO": f"{self.rho:.1f}",
+            "U0": f"{self.velocity}",
+            "COLLISION_MODEL": COLLISION_MAP[self.collision],
+            "STREAMING_MODEL": STREAMING_MAP[self.streaming],
+            "STORE_PRECISION": PRECISION_MAP[self.store_precision],
+            "USE_DDF_SHIFTING": int(self.ddf_shifting),
+            "USE_LES": int(self.les_enabled),
+            "LES_CS": f"{self.les_cs:.6f}f",
+            "INLET_PROFILE": INLET_MAP[self.inlet_profile],
+            "OUTLET_MODE": OUTLET_MAP[self.outlet_mode],
+            "OUTLET_BLEND_ALPHA": f"{self.outlet_blend_alpha:.3f}f",
+            "OUTLET_BACKFLOW_CLAMP": int(self.outlet_backflow_clamp),
+            "OMEGA_COLLISION_MIN": f"{self.omega_min:.2f}f",
+            "OMEGA_COLLISION_MAX": f"{self.omega_max:.3f}f",
+            "TRT_MAGIC_PARAM": f"{self.trt_magic_param:.6f}f",
+        }
+
+
+# ---------------------------------------------------------------------------
+# Body config
+# ---------------------------------------------------------------------------
+@dataclass
+class BodyConfig:
+    objects: List[Dict[str, Any]] = field(default_factory=list)
+
+
+# ---------------------------------------------------------------------------
+# Loaders
+# ---------------------------------------------------------------------------
+def _find_config(filename: str, explicit_path: Optional[str] = None) -> str:
+    """Search for config file in standard locations."""
+    candidates = []
+    if explicit_path:
+        candidates.append(explicit_path)
+    env_dir = os.environ.get("CELERISLAB_CONFIG_DIR")
+    if env_dir:
+        candidates.append(os.path.join(env_dir, filename))
+    candidates.append(os.path.join(os.getcwd(), "configs", filename))
+    pkg_root = os.path.dirname(os.path.abspath(__file__))
+    candidates.append(os.path.join(pkg_root, "configs", filename))
+    for p in candidates:
+        if os.path.isfile(p):
+            return os.path.abspath(p)
+    raise FileNotFoundError(
+        f"Config '{filename}' not found.  Searched:\n"
+        + "\n".join(f"  - {p}" for p in candidates)
+    )
+
+
+def load_lbm_config(path: Optional[str] = None) -> LBMConfig:
+    fpath = _find_config("config_lbm.json", path)
+    with open(fpath) as f:
+        d = json.load(f)
+    g, p, m, c = d["grid"], d["physics"], d["method"], d["cuda"]
+    cfg = LBMConfig(
+        dim=g["dim"], nq=g["nq"], nx=g["nx"], ny=g["ny"], nz=g["nz"],
+        data_type=p["data_type"], viscosity=p["viscosity"],
+        velocity=p["velocity"], rho=p["rho"],
+        collision=m["collision"], streaming=m["streaming"],
+        store_precision=m["store_precision"],
+        ddf_shifting=m["ddf_shifting"],
+        les_enabled=m["les"]["enabled"], les_cs=m["les"]["cs"],
+        trt_magic_param=m["trt"]["magic_param"],
+        inlet_profile=m["inlet"]["profile"],
+        outlet_mode=m["outlet"]["mode"],
+        outlet_blend_alpha=m["outlet"]["blend_alpha"],
+        outlet_backflow_clamp=m["outlet"]["backflow_clamp"],
+        omega_min=m["omega_guard"]["min"],
+        omega_max=m["omega_guard"]["max"],
+        threads_per_block=c["threads_per_block"],
+        compute_capability=c["compute_capability"],
+    )
+    cfg.validate()
+    return cfg
+
+
+def load_body_config(path: Optional[str] = None) -> BodyConfig:
+    fpath = _find_config("config_body.json", path)
+    with open(fpath) as f:
+        d = json.load(f)
+    return BodyConfig(objects=d.get("objects", []))
--- a/src/CelerisLab/configs/config_body.json
+++ b/src/CelerisLab/configs/config_body.json
@ -0,0 +1,3 @@
+{
+    "objects": []
+}
--- a/src/CelerisLab/configs/config_lbm.json
+++ b/src/CelerisLab/configs/config_lbm.json
@ -0,0 +1,34 @@
+{
+    "grid": {
+        "dim": 2,
+        "nq": 9,
+        "nx": 512,
+        "ny": 256,
+        "nz": 1
+    },
+    "physics": {
+        "data_type": "FP32",
+        "viscosity": 0.002,
+        "velocity": 0.03,
+        "rho": 1.0
+    },
+    "method": {
+        "collision": "SRT",
+        "streaming": "double_buffer",
+        "store_precision": "FP32",
+        "ddf_shifting": false,
+        "les": {"enabled": false, "cs": 0.16},
+        "trt": {"magic_param": 0.1875},
+        "inlet": {"profile": "parabolic"},
+        "outlet": {
+            "mode": "neq_extrap",
+            "backflow_clamp": true,
+            "blend_alpha": 0.7
+        },
+        "omega_guard": {"min": 0.01, "max": 1.999}
+    },
+    "cuda": {
+        "threads_per_block": 128,
+        "compute_capability": "auto"
+    }
+}
--- a/src/CelerisLab/cuda/init.py
+++ b/src/CelerisLab/cuda/init.py
@ -3,6 +3,6 @@
 CUDA compiler and kernel management utilities.
 """

-from . import compiler
+from . import compiler_v2 as compiler

 __all__ = ['compiler']
--- a/src/CelerisLab/cuda/compiler_v2.py
+++ b/src/CelerisLab/cuda/compiler_v2.py
@ -0,0 +1,172 @@
+# CelerisLab/cuda/compiler.py
+"""
+Kernel configuration, compilation, and module loading.
+
+Workflow:
+  1. generate_config(lbm_cfg, n_objects) — write config/*.h from LBMConfig
+  2. compile_kernel(arch)               — nvcc -ptx kernel_v2.cu
+  3. load_module(ptx_path)              — PyCUDA module from PTX
+"""
+
+import os
+import subprocess
+from typing import Optional
+
+import pycuda.driver as cuda
+
+from ..config import LBMConfig
+
+# ---------------------------------------------------------------------------
+# Paths
+# ---------------------------------------------------------------------------
+_KERNEL_DIR = os.path.join(
+    os.path.dirname(os.path.dirname(os.path.abspath(__file__))),
+    "lbm", "kernels",
+)
+
+
+def kernel_path(filename: str) -> str:
+    return os.path.join(_KERNEL_DIR, filename)
+
+
+# ---------------------------------------------------------------------------
+# Header generation (f-string templates, no external deps)
+# ---------------------------------------------------------------------------
+_HEADER = "// AUTO-GENERATED by CelerisLab compiler – DO NOT EDIT MANUALLY\n"
+
+
+def _write(path: str, content: str):
+    os.makedirs(os.path.dirname(path), exist_ok=True)
+    with open(path, "w") as f:
+        f.write(content)
+
+
+def generate_config(cfg: LBMConfig, n_objects: int = 0):
+    """Render all config/*.h files from *cfg*."""
+    m = cfg.to_macros()
+    cdir = os.path.join(_KERNEL_DIR, "config")
+
+    _write(os.path.join(cdir, "config_grid.h"), f"""{_HEADER}\
+// Layer: Global / Grid
+#ifndef CELERIS_CONFIG_GRID_H
+#define CELERIS_CONFIG_GRID_H
+
+#define NT   {m['NT']}
+#define MULT_GPU {m['MULT_GPU']}
+
+#define NX  {m['NX']}
+#define NY  {m['NY']}
+#define NZ  {m['NZ']}
+#define DIM {m['DIM']}
+#define NQ  {m['NQ']}
+
+#endif
+""")
+
+    _write(os.path.join(cdir, "config_physics.h"), f"""{_HEADER}\
+// Layer: Global / Physics
+#ifndef CELERIS_CONFIG_PHYSICS_H
+#define CELERIS_CONFIG_PHYSICS_H
+
+#define LBtype {m['LBtype']}
+#define VIS {m['VIS']}
+#define RHO {m['RHO']}
+#define U0  {m['U0']}
+
+#define PI 3.141592653589793238
+
+#define FLUID     0x01
+#define SOLID     0x02
+#define GAS       0x04
+#define INTERFACE 0x08
+#define SENSOR    0x10
+#define OBSTACLE  0x20
+
+#define V_TAYLOR  1
+
+#endif
+""")
+
+    _write(os.path.join(cdir, "config_method.h"), f"""{_HEADER}\
+// Layer: Method
+#ifndef CELERIS_CONFIG_METHOD_H
+#define CELERIS_CONFIG_METHOD_H
+
+#define COLLISION_MODEL  {m['COLLISION_MODEL']}
+#define STREAMING_MODEL  {m['STREAMING_MODEL']}
+#define STORE_PRECISION  {m['STORE_PRECISION']}
+#define USE_DDF_SHIFTING {m['USE_DDF_SHIFTING']}
+
+#define USE_LES {m['USE_LES']}
+#define LES_CS  {m['LES_CS']}
+
+#define INLET_PROFILE        {m['INLET_PROFILE']}
+#define OUTLET_MODE          {m['OUTLET_MODE']}
+#define OUTLET_BLEND_ALPHA   {m['OUTLET_BLEND_ALPHA']}
+#define OUTLET_BACKFLOW_CLAMP {m['OUTLET_BACKFLOW_CLAMP']}
+
+#define OMEGA_COLLISION_MIN {m['OMEGA_COLLISION_MIN']}
+#define OMEGA_COLLISION_MAX {m['OMEGA_COLLISION_MAX']}
+#define TRT_MAGIC_PARAM     {m['TRT_MAGIC_PARAM']}
+
+#endif
+""")
+
+    _write(os.path.join(cdir, "config_objects.h"), f"""{_HEADER}\
+// Layer: Case / Objects
+#ifndef CELERIS_CONFIG_OBJECTS_H
+#define CELERIS_CONFIG_OBJECTS_H
+
+#define N_OBJS {n_objects}
+
+#endif
+""")
+
+
+# ---------------------------------------------------------------------------
+# Compilation
+# ---------------------------------------------------------------------------
+def compile_kernel(arch: str = "sm_70",
+                   source: str = "kernel_v2.cu",
+                   output: Optional[str] = None) -> str:
+    """Compile *source* to PTX.  Returns path to the .ptx file."""
+    src = kernel_path(source)
+    if output is None:
+        output = source.replace(".cu", ".ptx")
+    dst = kernel_path(output)
+    subprocess.run(
+        ["nvcc", "-ptx", f"-arch={arch}", src, "-o", dst],
+        check=True,
+    )
+    return dst
+
+
+# ---------------------------------------------------------------------------
+# Module loading
+# ---------------------------------------------------------------------------
+def load_module(ptx_path: Optional[str] = None) -> cuda.Module:
+    """Load a compiled PTX as a PyCUDA module."""
+    if ptx_path is None:
+        ptx_path = kernel_path("kernel_v2.ptx")
+    return cuda.module_from_file(ptx_path)
+
+
+# ---------------------------------------------------------------------------
+# Back-compat helpers (used by test scripts)
+# ---------------------------------------------------------------------------
+
+def compile_kernel_v2(arch: str = "sm_70") -> str:
+    """Alias for compile_kernel() kept for test-script compatibility."""
+    return compile_kernel(arch=arch)
+
+
+def read_lines(path: str):
+    """Read a text file and return a list of lines (with newlines)."""
+    with open(path, "r") as f:
+        return f.readlines()
+
+
+def write_lines(path: str, lines):
+    """Write a list of lines back to a text file."""
+    with open(path, "w") as f:
+        f.writelines(lines)
--- a/src/CelerisLab/cuda/context.py
+++ b/src/CelerisLab/cuda/context.py
@ -0,0 +1,79 @@
+# CelerisLab/cuda/context.py
+"""
+CUDA device selection, context lifecycle, and device info queries.
+"""
+
+import subprocess
+from typing import Optional
+
+import pycuda.driver as cuda
+
+
+class CudaContext:
+    """Manages a single CUDA device context.
+
+    Usage::
+
+        ctx = CudaContext(device_id=0)
+        # ... use PyCUDA ...
+        ctx.close()
+
+    Or as a context-manager::
+
+        with CudaContext(0) as ctx:
+            ...
+    """
+
+    def __init__(self, device_id: int = 0):
+        cuda.init()
+        if device_id < 0 or device_id >= cuda.Device.count():
+            raise ValueError(
+                f"device_id {device_id} invalid; "
+                f"{cuda.Device.count()} device(s) available."
+            )
+        self.device_id = device_id
+        self.device = cuda.Device(device_id)
+        self._ctx = self.device.make_context()
+
+    # -- properties ----------------------------------------------------------
+    @property
+    def name(self) -> str:
+        return self.device.name()
+
+    @property
+    def compute_capability(self) -> str:
+        cc = self.device.compute_capability()
+        return f"{cc[0]}.{cc[1]}"
+
+    @property
+    def sm_arch(self) -> str:
+        """Return 'sm_XX' string for nvcc -arch."""
+        cc = self.device.compute_capability()
+        return f"sm_{cc[0]}{cc[1]}"
+
+    @property
+    def total_memory(self) -> int:
+        return self.device.total_memory()
+
+    # -- lifecycle -----------------------------------------------------------
+    def close(self):
+        if self._ctx is not None:
+            self._ctx.pop()
+            self._ctx = None
+
+    def __enter__(self):
+        return self
+
+    def __exit__(self, *exc):
+        self.close()
+
+    def __del__(self):
+        self.close()
+
+
+def detect_compute_capability(device_id: int = 0) -> str:
+    """Return 'sm_XX' for the given device without creating a persistent context."""
+    cuda.init()
+    dev = cuda.Device(device_id)
+    cc = dev.compute_capability()
+    return f"sm_{cc[0]}{cc[1]}"
--- a/src/CelerisLab/lbm/init.py
+++ b/src/CelerisLab/lbm/init.py
@ -4,13 +4,14 @@ Lattice Boltzmann Method (LBM) module for fluid simulation.
 """

 try:
-    from .driver import FlowField
-    __all__ = ['FlowField']
+    from .field import LBMField
+    from .stepper import LBMStepper
+    from .initializers import add_lamb_oseen, add_taylor_green
+    __all__ = ["LBMField", "LBMStepper", "add_lamb_oseen", "add_taylor_green"]
 except ImportError as e:
    import warnings
-    warnings.warn(
-        f"LBM module not fully available: {e}. "
-        "Install pycuda to use the full LBM functionality.",
-        ImportWarning
-    )
+    warnings.warn(f"LBM module not fully available: {e}", ImportWarning)
    __all__ = []
+
+# Legacy import — kept while test scripts still reference FlowField
+
--- a/src/CelerisLab/lbm/field.py
+++ b/src/CelerisLab/lbm/field.py
@ -0,0 +1,163 @@
+# CelerisLab/lbm/field.py
+"""
+LBMField — GPU memory for distribution functions, flags, and macroscopic data.
+
+Owns:
+  ddf_gpu / temp_gpu   – distribution-function double buffers
+  flag_gpu             – cell-type flags
+  indx_gpu             – per-cell object index (for curved BC / sensors)
+  delta_gpu            – curved-boundary interpolation data
+
+Provides:
+  upload_params()      – push LBMParams to __constant__ memory  (fixes d_params bug)
+  get_ddf() / set_ddf()
+  get_macroscopic()    – rho, ux, uy [, uz]
+  snapshot() / restore()
+"""
+
+import ctypes
+import struct
+import numpy as np
+import pycuda.driver as cuda
+
+from ..config import LBMConfig
+
+
+class LBMField:
+    """Allocate and manage GPU arrays for one LBM domain."""
+
+    def __init__(self, cfg: LBMConfig, module: cuda.Module):
+        self.cfg = cfg
+        self.module = module
+
+        self.nx, self.ny, self.nz = cfg.nx, cfg.ny, cfg.nz
+        self.nq = cfg.nq
+        self.dim = cfg.dim
+        self.n = self.nx * self.ny * self.nz
+        self.dtype = np.float32  # extend when FP64 supported
+
+        # Host arrays
+        self.ddf = np.zeros(self.n * self.nq, dtype=self.dtype)
+        self.flag = np.ones(self.n, dtype=np.uint8) * 0x01  # FLUID
+        self.indx = np.zeros(self.n, dtype=np.int32)
+        self.delta = np.zeros(0, dtype=self.dtype)
+
+        # GPU allocations
+        self.ddf_gpu = cuda.mem_alloc(self.ddf.nbytes)
+        self.temp_gpu = cuda.mem_alloc(self.ddf.nbytes)
+        self.flag_gpu = cuda.mem_alloc(self.flag.nbytes)
+        self.indx_gpu = cuda.mem_alloc(self.indx.nbytes)
+        self.delta_gpu = cuda.mem_alloc(max(4, self.delta.nbytes))
+
+        # Snapshot
+        self._ddf_snap: np.ndarray | None = None
+
+        # Upload d_params immediately
+        self._upload_params()
+
+    # -- d_params upload -----------------------------------------------------
+    def _upload_params(self):
+        """Pack LBMParams struct and upload to __constant__ d_params."""
+        cfg = self.cfg
+        # Must match struct LBMParams in core/params.cuh layout:
+        #   uint Nx,Ny,Nz; ulong N; float omega,omega_bulk;
+        #   float fx,fy,fz; float rho_ref,u_inlet; uint n_objects;
+        fmt = "IIILff fff ff I"
+        data = struct.pack(
+            fmt,
+            cfg.nx, cfg.ny, cfg.nz,
+            self.n,
+            cfg.omega, 0.0,          # omega, omega_bulk
+            0.0, 0.0, 0.0,           # fx, fy, fz
+            cfg.rho, cfg.velocity,    # rho_ref, u_inlet
+            0,                        # n_objects (updated later)
+        )
+        ptr, size = self.module.get_global("d_params")
+        cuda.memcpy_htod(ptr, data)
+
+    def update_params(self, **kwargs):
+        """Re-upload d_params after changing runtime-adjustable values.
+
+        Accepted keys: omega, omega_bulk, fx, fy, fz, rho_ref, u_inlet, n_objects.
+        """
+        cfg = self.cfg
+        omega = kwargs.get("omega", cfg.omega)
+        omega_bulk = kwargs.get("omega_bulk", 0.0)
+        fx = kwargs.get("fx", 0.0)
+        fy = kwargs.get("fy", 0.0)
+        fz = kwargs.get("fz", 0.0)
+        rho_ref = kwargs.get("rho_ref", cfg.rho)
+        u_inlet = kwargs.get("u_inlet", cfg.velocity)
+        n_objects = kwargs.get("n_objects", 0)
+
+        fmt = "IIILff fff ff I"
+        data = struct.pack(
+            fmt,
+            cfg.nx, cfg.ny, cfg.nz, self.n,
+            omega, omega_bulk,
+            fx, fy, fz,
+            rho_ref, u_inlet,
+            n_objects,
+        )
+        ptr, _ = self.module.get_global("d_params")
+        cuda.memcpy_htod(ptr, data)
+
+    # -- Host ↔ GPU transfers ------------------------------------------------
+    def upload_ddf(self):
+        cuda.memcpy_htod(self.ddf_gpu, self.ddf)
+        cuda.memcpy_htod(self.temp_gpu, self.ddf)
+
+    def download_ddf(self):
+        cuda.memcpy_dtoh(self.ddf, self.ddf_gpu)
+
+    def upload_flags(self):
+        cuda.memcpy_htod(self.flag_gpu, self.flag)
+        cuda.memcpy_htod(self.indx_gpu, self.indx)
+
+    def upload_delta(self):
+        if self.delta.nbytes > 0:
+            self.delta_gpu.free()
+            self.delta_gpu = cuda.mem_alloc(self.delta.nbytes)
+            cuda.memcpy_htod(self.delta_gpu, self.delta)
+
+    # -- Macroscopic field extraction ----------------------------------------
+    def get_macroscopic(self):
+        """Download DDF and compute rho, ux, uy [, uz] on host."""
+        self.download_ddf()
+        nq, n = self.nq, self.n
+        shape = (self.nz, self.ny, self.nx) if self.dim == 3 else (self.ny, self.nx)
+
+        if nq == 9:
+            f = self.ddf.reshape(9, self.ny, self.nx)
+            rho = np.sum(f, axis=0)
+            cx = np.array([0, 1, -1, 0, 0, 1, -1, 1, -1], dtype=np.float32)
+            cy = np.array([0, 0, 0, 1, -1, 1, -1, -1, 1], dtype=np.float32)
+            rho_safe = np.where(np.abs(rho) > 1e-12, rho, 1.0)
+            ux = sum(f[i] * cx[i] for i in range(9)) / rho_safe
+            uy = sum(f[i] * cy[i] for i in range(9)) / rho_safe
+            return {"rho": rho, "ux": ux, "uy": uy}
+
+        if nq == 19:
+            f = self.ddf.reshape(19, self.nz, self.ny, self.nx)
+            rho = np.sum(f, axis=0)
+            cx = np.array([0,1,-1,0,0,0,0,1,-1,1,-1,0,0,1,-1,1,-1,0,0])
+            cy = np.array([0,0,0,1,-1,0,0,1,-1,0,0,1,-1,-1,1,0,0,1,-1])
+            cz = np.array([0,0,0,0,0,1,-1,0,0,1,-1,1,-1,0,0,-1,1,-1,1])
+            rho_safe = np.where(np.abs(rho) > 1e-12, rho, 1.0)
+            ux = sum(f[i]*cx[i] for i in range(19)) / rho_safe
+            uy = sum(f[i]*cy[i] for i in range(19)) / rho_safe
+            uz = sum(f[i]*cz[i] for i in range(19)) / rho_safe
+            return {"rho": rho, "ux": ux, "uy": uy, "uz": uz}
+
+        raise ValueError(f"Unsupported nq={nq}")
+
+    # -- Snapshots -----------------------------------------------------------
+    def snapshot(self):
+        self.download_ddf()
+        self._ddf_snap = self.ddf.copy()
+
+    def restore(self):
+        if self._ddf_snap is None:
+            raise RuntimeError("No snapshot to restore")
+        self.ddf = self._ddf_snap.copy()
+        self.upload_ddf()
--- a/src/CelerisLab/lbm/initializers.py
+++ b/src/CelerisLab/lbm/initializers.py
@ -0,0 +1,134 @@
+# CelerisLab/lbm/initializers.py
+"""
+Flow-field initialization helpers (vortex superposition, etc.).
+
+All functions operate on an LBMField — they download DDF, modify on host,
+and re-upload.  Only 2D D2Q9 vortices implemented for now.
+"""
+
+import numpy as np
+from scipy.special import jv, expi
+from typing import Tuple
+
+import pycuda.driver as cuda
+
+# D2Q9 velocity vectors and weights
+_E9 = np.array(
+    [0, 0, 1, 0, 0, 1, -1, 0, 0, -1, 1, 1, -1, 1, -1, -1, 1, -1],
+    dtype=np.int32,
+).reshape(9, 2)
+_W9 = np.array(
+    [4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36],
+    dtype=np.float32,
+)
+
+
+def add_vortex(field, center: Tuple[float, float],
+               radius: float, strength: float,
+               vortex_type: str = "lamb"):
+    """Superimpose a vortex onto an existing 2D D2Q9 flow field.
+
+    Supported types: "lamb", "oseen", "taylor".
+    """
+    if field.cfg.nq != 9 or field.cfg.dim != 2:
+        raise NotImplementedError("Vortex init only for 2D D2Q9")
+    if vortex_type not in ("lamb", "oseen", "taylor"):
+        raise ValueError(f"Unknown vortex type: {vortex_type}")
+
+    nx, ny = field.nx, field.ny
+    x_c, y_c = center
+    nu = field.cfg.viscosity
+
+    x = np.linspace(-x_c, nx - 1 - x_c, nx)
+    y = np.linspace(-y_c, ny - 1 - y_c, ny)
+    X, Y = np.meshgrid(x, y)
+    r = np.sqrt(X**2 + Y**2)
+    theta = np.arctan2(Y, X)
+
+    if vortex_type == "lamb":
+        u_vor, v_vor, p_vor = _lamb_vortex(r, theta, radius, strength)
+    elif vortex_type == "oseen":
+        u_vor, v_vor, p_vor = _oseen_vortex(r, theta, radius, strength)
+    else:
+        u_vor, v_vor, p_vor = _taylor_vortex(r, theta, radius, strength, nu)
+
+    # Download current DDF
+    field.download_ddf()
+    f = field.ddf.reshape(9, ny, nx)
+
+    # Compute current macroscopic from DDF
+    rho_old = np.sum(f, axis=0)
+    ux_old = (f[1]+f[5]+f[8]-f[3]-f[6]-f[7])
+    uy_old = (f[2]+f[5]+f[6]-f[4]-f[7]-f[8])
+    p_old = rho_old / 3.0
+
+    # Superimpose
+    # Vortex arrays are (ny, nx); u_vor[j, i] adds to node (i, j)
+    u_new = ux_old + u_vor
+    v_new = uy_old + v_vor
+    p_new = p_old + p_vor
+
+    # Reconstruct DDF from new macroscopic (equilibrium)
+    for j in range(1, ny - 1):
+        for i in range(1, nx - 1):
+            u, v, p = float(u_new[j, i]), float(v_new[j, i]), float(p_new[j, i])
+            for e in range(9):
+                eu = _E9[e, 0] * u + _E9[e, 1] * v
+                u2 = u * u + v * v
+                f[e, j, i] = _W9[e] * (3*p + 3*eu + 4.5*eu*eu - 1.5*u2)
+
+    field.ddf = f.reshape(-1)
+    field.upload_ddf()
+
+
+# ---------------------------------------------------------------------------
+# Vortex models
+# ---------------------------------------------------------------------------
+def _lamb_vortex(r, theta, radius, strength):
+    b = 3.831705970207512
+    n = b / radius
+    u0 = strength
+    inside = r <= radius
+    outside = r > radius
+    psi = np.zeros_like(r)
+    psi[inside] = (
+        (2*u0 / n / jv(0, b) * jv(1, n*r[inside]) - u0*r[inside])
+        * np.sin(theta[inside])
+    )
+    psi[outside] = -u0 * radius**2 / r[outside] * np.sin(theta[outside])
+    u_vor = np.gradient(psi, axis=0)
+    v_vor = -np.gradient(psi, axis=1)
+    p_vor = (
+        -2*(np.gradient(v_vor, axis=1) - np.gradient(u_vor, axis=0))*psi
+        - (u_vor**2 + v_vor**2) / 2
+    )
+    return u_vor, v_vor, p_vor
+
+
+def _oseen_vortex(r, theta, radius, strength):
+    kappa = 2 * np.pi * radius**2 * strength
+    r_safe = np.where(r > 1e-12, r, 1e-12)
+    exp_term = 1 - np.exp(-4 * r**2 / radius**2)
+    u_vor = -kappa / (2*np.pi*r_safe) * exp_term * np.sin(theta)
+    v_vor = kappa / (2*np.pi*r_safe) * exp_term * np.cos(theta)
+    zeta = 4 * r**2 / radius**2
+    p_vor = (
+        -kappa**2 / (8*np.pi**2 * r_safe**2)
+        * (-2*zeta*(expi(-zeta) - expi(-2*zeta))
+           + (1 - np.exp(-zeta))**2)
+    )
+    return u_vor, v_vor, p_vor
+
+
+def _taylor_vortex(r, theta, radius, strength, nu):
+    M = strength * np.pi * radius**4 / (8 * nu)
+    coeff = M * 4 * nu / radius**4
+    exp_term = np.exp(-r**2 / radius**2)
+    u_vor = -coeff * r * exp_term * np.sin(theta)
+    v_vor = coeff * r * exp_term * np.cos(theta)
+    p_vor = (
+        -4 * M**2 * nu**2
+        * np.exp(-2 * r**2 / radius**2)
+        / (np.pi**2 * radius**6)
+    )
+    return u_vor, v_vor, p_vor
--- a/src/CelerisLab/lbm/kernels/D2Q9.cu
+++ b/src/CelerisLab/lbm/kernels/D2Q9.cu
@ -1,101 +0,0 @@
-#include "macros.h"
-#include "const.h"
-
-__device__ void Index_lattice(int &x, int &y, int &k) {
-	// Only for D2
-	x = threadIdx.x + NT * blockIdx.x;
-	y = blockIdx.y;
-	k = y * NX + x;
-}
-
-__device__ void CollisionKernel(LBtype* g, LBtype* m) {
-	// Only for D2Q9
-	LBtype p, u, v;
-    LBtype niu = 1.0 / (0.5 + 3 * VIS);
-
-	u = (g[1]+g[5]+g[8]-g[3]-g[6]-g[7])/RHO;
-	v = (g[2]+g[5]+g[6]-g[4]-g[7]-g[8])/RHO;
-	p = (g[0]+g[1]+g[2]+g[3]+g[4]+g[5]+g[6]+g[7]+g[8])/3.0;
-
-	m[0]=   g[0]  +g[1]  +g[2]  +g[3]  +g[4]  +g[5]  +g[6]  +g[7]  +g[8];
-	m[1]=-4*g[0]  -g[1]  -g[2]  -g[3]  -g[4]+2*g[5]+2*g[6]+2*g[7]+2*g[8];
-	m[2]= 4*g[0]-2*g[1]-2*g[2]-2*g[3]-2*g[4]  +g[5]  +g[6]  +g[7]  +g[8];
-	m[3]=          g[1]         -g[3]         +g[5]  -g[6]  -g[7]  +g[8];
-	m[4]=       -2*g[1]       +2*g[3]         +g[5]  -g[6]  -g[7]  +g[8];
-	m[5]=                 g[2]         -g[4]  +g[5]  +g[6]  -g[7]  -g[8];
-	m[6]=              -2*g[2]       +2*g[4]  +g[5]  +g[6]  -g[7]  -g[8];
-	m[7]=          g[1]  -g[2]  +g[3]  -g[4];
-	m[8]=                                      g[5]  -g[6]  +g[7]  -g[8];	
-	
-	m[0]=1.00*( 3*p 	            -m[0]);
-	m[1]=1.20*(-6*p +3*RHO*(u*u+v*v)-m[1]);
-	m[2]=1.20*( 3*p -3*RHO*(u*u+v*v)-m[2]);
-	m[3]=1.00*( RHO*u               -m[3]);
-	m[4]=1.20*(-RHO*u               -m[4]);
-	m[5]=1.00*( RHO*v               -m[5]);
-	m[6]=1.20*(-RHO*v               -m[6]);
-	m[7]= niu*( RHO*(u*u-v*v)       -m[7]);
-	m[8]= niu*( RHO*u*v             -m[8]);
-
-	g[0]=g[0]+(  m[0]  -m[1]  +m[2]                                   )/ 9.0;
-	g[1]=g[1]+(4*m[0]  -m[1]-2*m[2]+6*m[3]-6*m[4]              +9*m[7])/36.0;
-	g[2]=g[2]+(4*m[0]  -m[1]-2*m[2]              +6*m[5]-6*m[6]-9*m[7])/36.0;
-	g[3]=g[3]+(4*m[0]  -m[1]-2*m[2]-6*m[3]+6*m[4]              +9*m[7])/36.0;
-	g[4]=g[4]+(4*m[0]  -m[1]-2*m[2]              -6*m[5]+6*m[6]-9*m[7])/36.0;
-	g[5]=g[5]+(4*m[0]+2*m[1]  +m[2]+6*m[3]+3*m[4]+6*m[5]+3*m[6]+9*m[8])/36.0;
-	g[6]=g[6]+(4*m[0]+2*m[1]  +m[2]-6*m[3]-3*m[4]+6*m[5]+3*m[6]-9*m[8])/36.0;
-	g[7]=g[7]+(4*m[0]+2*m[1]  +m[2]-6*m[3]-3*m[4]-6*m[5]-3*m[6]+9*m[8])/36.0;
-	g[8]=g[8]+(4*m[0]+2*m[1]  +m[2]+6*m[3]+3*m[4]-6*m[5]-3*m[6]-9*m[8])/36.0;
-}
-
-__device__ void ParabolicInlet(LBtype* f, LBtype* f_neb, LBtype y) {
-	LBtype p, u, v, yy;
-	LBtype feq1, feq5, feq8, feqn1, feqn5, feqn8;
-
-	p=(f_neb[0]+f_neb[1]+f_neb[2]+f_neb[3]+f_neb[4]+f_neb[5]+f_neb[6]+f_neb[7]+f_neb[8])/3.0;
-	yy=(y-0.5*(NY-1))/(NY-2.0);
-	u=U0*1.5*(1-4*yy*yy);
-	v=0.0;
-
-	feq1=(2*p+RHO*(2*u*u+2*u    -v*v)  )/ 6.0;
-	feq5=(  p+RHO*(  u*u+3*u*v+u+v*v+v))/12.0;
-	feq8=(  p+RHO*(  u*u-3*u*v+u+v*v-v))/12.0;
-
-	u=(f_neb[1]+f_neb[5]+f_neb[8]-f_neb[3]-f_neb[6]-f_neb[7])/RHO;
-	v=(f_neb[2]+f_neb[5]+f_neb[6]-f_neb[4]-f_neb[7]-f_neb[8])/RHO;
-
-	feqn1=(2*p+RHO*(2*u*u+2*u    -v*v)  )/ 6.0;
-	feqn5=(  p+RHO*(  u*u+3*u*v+u+v*v+v))/12.0;
-	feqn8=(  p+RHO*(  u*u-3*u*v+u+v*v-v))/12.0;
-
-	f[1]=f_neb[1]-feqn1+feq1;
-	f[5]=f_neb[5]-feqn5+feq5;
-	f[8]=f_neb[8]-feqn8+feq8;
-}
-
-__device__ void PressureOutlet(LBtype* f, LBtype* f_neb, LBtype y) {
-	// Edit to Parabolic Outlet temporarily
-	LBtype p, u, v, yy;	
-	LBtype feq3, feq6, feq7, feqn3, feqn6, feqn7;
-
-	p=0.0;
-
-	yy=(y-0.5*(NY-1))/(NY-2.0);
-	u=U0*1.5*(1-4*yy*yy);
-	v=0.0;
-
-	feq3=(2*p-RHO*(-2*u*u+2*u    +v*v)  )/ 6.0;
-	feq6=(  p+RHO*(   u*u-3*u*v-u+v*v+v))/12.0;
-	feq7=(  p+RHO*(   u*u+3*u*v-u+v*v-v))/12.0;
-
-	u=(f_neb[1]+f_neb[5]+f_neb[8]-f_neb[3]-f_neb[6]-f_neb[7])/RHO;
-	v=(f_neb[2]+f_neb[5]+f_neb[6]-f_neb[4]-f_neb[7]-f_neb[8])/RHO;
-	// p=(f_neb[0]+f_neb[1]+f_neb[2]+f_neb[3]+f_neb[4]+f_neb[5]+f_neb[6]+f_neb[7]+f_neb[8])/3.0;
-	feqn3=(2*p-RHO*(-2*u*u+2*u    +v*v)  )/ 6.0;
-	feqn6=(  p+RHO*(   u*u-3*u*v-u+v*v+v))/12.0;
-	feqn7=(  p+RHO*(   u*u+3*u*v-u+v*v-v))/12.0;
-
-	f[3]=f_neb[3]-feqn3+feq3;
-	f[6]=f_neb[6]-feqn6+feq6;
-	f[7]=f_neb[7]-feqn7+feq7;
-}
--- a/src/CelerisLab/lbm/kernels/IO.cu
+++ b/src/CelerisLab/lbm/kernels/IO.cu
--- a/src/CelerisLab/lbm/kernels/boundary/inlet_outlet.cuh
+++ b/src/CelerisLab/lbm/kernels/boundary/inlet_outlet.cuh
@ -43,8 +43,13 @@ __device__ __forceinline__ float inlet_target_u(float y_coord) {
 #if INLET_PROFILE == 0
    return U0;
 #else
-    float yy = (y_coord - 0.5f * (NY - 1)) / (NY - 2.0f);
-    return U0 * 1.5f * (1.0f - 4.0f * yy * yy);
+    // Define profile on fluid-node band y in [1, NY-2] and clamp to non-negative
+    // to avoid near-corner backflow injection.
+    const float y_clamped = fminf((float)(NY - 2), fmaxf(1.0f, y_coord));
+    const float H = fmaxf((float)(NY - 2), 1.0f);
+    const float eta = (y_clamped - 0.5f) / H;  // maps first/last fluid rows near 0/1
+    const float shape = fmaxf(0.0f, 4.0f * eta * (1.0f - eta));
+    return U0 * 1.5f * shape;
 #endif
 }

@ -72,21 +77,18 @@ __device__ inline void apply_parabolic_inlet(float* __restrict__ f,
    // Target velocity (parabolic profile)
    float u_target = inlet_target_u(y_coord);
    float v_target = 0.0f;
+    const float rho_in = RHO;

    float feq_tar[9], feq_neb[9];
-    compute_feq(rho_neb, u_target, v_target, feq_tar);
+    compute_feq(rho_in, u_target, v_target, feq_tar);
    compute_feq(rho_neb, u_neb, v_neb, feq_neb);

 #if COLLISION_MODEL == 1
-    // TRT path: reconstruct full population set at inlet using damped donor
-    // non-equilibrium transport. This follows the high-Re stable family that
-    // replaces all boundary-node populations, reducing odd-mode contamination.
+    // TRT path: reconstruct unknown incoming populations only (cx>0 at west inlet).
    const float beta = INLET_TRT_NEQ_DAMP;
-    #pragma unroll
-    for (int i = 0; i < 9; i++) {
-        const float fneq = f_neb[i] - feq_neb[i];
-        f[i] = feq_tar[i] + beta * fneq;
-    }
+    f[1] = feq_tar[1] + beta * (f_neb[1] - feq_neb[1]);
+    f[5] = feq_tar[5] + beta * (f_neb[5] - feq_neb[5]);
+    f[7] = feq_tar[7] + beta * (f_neb[7] - feq_neb[7]);
 #else
    const float beta = 1.0f;
    f[1] = feq_tar[1] + beta * (f_neb[1] - feq_neb[1]);
@ -172,10 +174,11 @@ __device__ inline void apply_parabolic_inlet_3d(float* __restrict__ f,

    // Target velocity (parabolic in y, uniform in z)
    float u_tar = inlet_target_u(y_coord);
+    const float rho_in = RHO;

    // feq arrays
    float feq_tar[19], feq_neb[19];
-    compute_feq(rho_neb, u_tar, 0.0f, 0.0f, feq_tar);
+    compute_feq(rho_in, u_tar, 0.0f, 0.0f, feq_tar);
    compute_feq(rho_neb, un, vn, wn, feq_neb);

    // Reconstruct cx>0 directions: i = 1, 7, 9, 13, 15
@ -242,4 +245,4 @@ __device__ inline void apply_pressure_outlet_3d(float* __restrict__ f,

 #endif // NQ == 19

-#endif // CELERIS_BOUNDARY_INLET_OUTLET_CUH
+#endif // CELERIS_BOUNDARY_INLET_OUTLET_CUH
--- a/src/CelerisLab/lbm/kernels/config.h
+++ b/src/CelerisLab/lbm/kernels/config.h
@ -0,0 +1,18 @@
+// CelerisLab – config.h
+// ============================================================================
+// Top-level configuration aggregator.
+// Includes all layer-specific config headers generated by compiler.py.
+//
+// Users: read these files to see current parameters; do NOT edit them.
+//        Modify the JSON config and re-run the Python build step instead.
+// ============================================================================
+
+#ifndef CELERIS_CONFIG_H
+#define CELERIS_CONFIG_H
+
+#include "config/config_grid.h"
+#include "config/config_physics.h"
+#include "config/config_method.h"
+#include "config/config_objects.h"
+
+#endif // CELERIS_CONFIG_H
--- a/src/CelerisLab/lbm/kernels/config/config_grid.h
+++ b/src/CelerisLab/lbm/kernels/config/config_grid.h
@ -0,0 +1,11 @@
+// AUTO-GENERATED by test_stability_matrix.py
+#ifndef CELERIS_CONFIG_GRID_H
+#define CELERIS_CONFIG_GRID_H
+#define NT   128
+#define MULT_GPU 0
+#define NX  384
+#define NY  192
+#define NZ  1
+#define DIM 2
+#define NQ  9
+#endif
--- a/src/CelerisLab/lbm/kernels/config/config_method.h
+++ b/src/CelerisLab/lbm/kernels/config/config_method.h
@ -0,0 +1,17 @@
+// AUTO-GENERATED by test_stability_matrix.py
+#ifndef CELERIS_CONFIG_METHOD_H
+#define CELERIS_CONFIG_METHOD_H
+#define COLLISION_MODEL 0
+#define STREAMING_MODEL 1
+#define STORE_PRECISION 0
+#define USE_DDF_SHIFTING 0
+#define USE_LES 0
+#define LES_CS 0.160000f
+#define INLET_PROFILE 1
+#define OUTLET_MODE 0
+#define OUTLET_BLEND_ALPHA 0.700f
+#define OUTLET_BACKFLOW_CLAMP 1
+#define OMEGA_COLLISION_MIN 0.01f
+#define OMEGA_COLLISION_MAX 1.999f
+#define TRT_MAGIC_PARAM 0.187500f
+#endif
--- a/src/CelerisLab/lbm/kernels/config/config_objects.h
+++ b/src/CelerisLab/lbm/kernels/config/config_objects.h
@ -0,0 +1,5 @@
+// AUTO-GENERATED by test_stability_matrix.py
+#ifndef CELERIS_CONFIG_OBJECTS_H
+#define CELERIS_CONFIG_OBJECTS_H
+#define N_OBJS 0
+#endif
--- a/src/CelerisLab/lbm/kernels/config/config_physics.h
+++ b/src/CelerisLab/lbm/kernels/config/config_physics.h
@ -0,0 +1,16 @@
+// AUTO-GENERATED by test_stability_matrix.py
+#ifndef CELERIS_CONFIG_PHYSICS_H
+#define CELERIS_CONFIG_PHYSICS_H
+#define LBtype float
+#define VIS 0.0144000000
+#define RHO 1.0
+#define U0  0.04
+#define PI 3.141592653589793238
+#define FLUID       0x01
+#define SOLID       0x02
+#define GAS         0x04
+#define INTERFACE   0x08
+#define SENSOR      0x10
+#define OBSTACLE    0x20
+#define V_TAYLOR    1
+#endif
--- a/src/CelerisLab/lbm/kernels/const.h
+++ b/src/CelerisLab/lbm/kernels/const.h
@ -1,10 +0,0 @@
-// CelerisLab/kernels/const.h
-
-#ifndef CONST_H
-#define CONST_H
-
-__constant__ int e[9][2] = {{0, 0}, {1, 0}, {0, 1}, {-1, 0}, {0, -1}, {1, 1}, {-1, 1}, {-1, -1}, {1, -1}};
-__constant__ int opp[9] = {0, 3, 4, 1, 2, 7, 8, 5, 6};
-__constant__ float w[9] = {4/9., 1/9., 1/9., 1/9., 1/9., 1/36., 1/36., 1/36., 1/36.};
-
-#endif
--- a/src/CelerisLab/lbm/kernels/core/flags.cuh
+++ b/src/CelerisLab/lbm/kernels/core/flags.cuh
@ -47,11 +47,12 @@

 // ---------------------------------------------------------------------------
 // Legacy compatibility (current driver.py uses uint8 with these bits)
+// Migration target: uint16 flag field (see config/config_physics.h)
 // ---------------------------------------------------------------------------
 #define LEGACY_FLUID     0x01
 #define LEGACY_SOLID     0x02
 #define LEGACY_GAS       0x04
-#define LEGACY_OBSTACLE  0x04   // obstacle / immersed body (triggers BB at adjacent fluid)
+#define LEGACY_OBSTACLE  0x20   // obstacle / immersed body (0x20 avoids GAS collision)
 #define LEGACY_INTERFACE 0x08
 #define LEGACY_SENSOR    0x10

--- a/src/CelerisLab/lbm/kernels/kernel.cu
+++ b/src/CelerisLab/lbm/kernels/kernel.cu
@ -1,222 +0,0 @@
-// CelerisLab/kernels/kernel.cu
-
-#include <stdio.h>
-#include <stdint.h>
-#include <cuda.h>
-
-#include "macros.h"
-#include "const.h"
-#include "D2Q9.cu"
-
-extern "C"
-{
-	__global__ void OneStep(uint8_t *flag, LBtype *f, LBtype *f_temp, int32_t *indx, LBtype *delta, LBtype *action, LBtype *obs)
-	{
-		__shared__ LBtype f_share[NT * NQ];
-		__shared__ LBtype obs_share[(N_OBJS * DIM > 0) ? N_OBJS * DIM : 1];
-
-		int x, y, k;
-		LBtype g[NQ], m[NQ];
-		Index_lattice(x, y, k); // Only for D2
-		int totalCells = NX * NY;
-		int id = indx[k];
-
-		for (int i = 0; i < NQ; i++)
-		{
-			f_share[threadIdx.x + i * NT] = f[k + i * totalCells];
-		}
-		for (int i = threadIdx.x; i < N_OBJS * DIM; i+=NT)
-		{
-			obs_share[i] = 0;
-		}
-
-		__syncthreads();
-
-		for (int i = 0; i < NQ; i++)
-		{
-			g[i] = f_share[threadIdx.x + i * NT];
-		}
-
-		if (flag[k] & FLUID)
-		{
-			CollisionKernel(g, m);
-
-			for (int i = 0; i < NQ; i++)
-			{
-				f_share[threadIdx.x + i * NT] = g[i];
-			}
-		}
-		else if (flag[k] & SOLID)
-		{
-			if (x == 0)
-			{
-				for (int i = 0; i < NQ; i++)
-				{
-					m[i] = f_share[threadIdx.x + i * NT + 1];
-				}
-				ParabolicInlet(g, m, y);
-			}
-			else if (x == NX - 1)
-			{
-				for (int i = 0; i < NQ; i++)
-				{
-					m[i] = f_share[threadIdx.x + i * NT - 1];
-				}
-				PressureOutlet(g, m, y);
-			}
-
-			for (int i = 0; i < NQ; i++)
-			{
-				f_share[threadIdx.x + i * NT] = g[i];
-			}
-		}
-
-		__syncthreads();
-
-		for (int i = 0; i < NQ; i++)
-		{
-			int x_neb = x + e[i][0];
-			int y_neb = y + e[i][1];
-
-			if (y != 0 && y != NY - 1)
-			{
-				if ((y == 1 && y_neb == 0) || (y == NY - 2 && y_neb == NY - 1))
-				{
-					f_temp[k + opp[i] * totalCells] = f_share[threadIdx.x + i * NT];
-				}
-				else
-				{
-					int k_neb = ((y_neb * NX + x_neb) + totalCells) % totalCells;
-					f_temp[k_neb + i * totalCells] = f_share[threadIdx.x + i * NT];
-				}
-			}
-		}
-
-		__syncthreads();
-
-		if (flag[k] & SOLID && flag[k] & INTERFACE)
-		{
-			LBtype Uw, Vw;
-			int id_obj = *reinterpret_cast<int*>(&delta[id]);
-			Uw = action[id_obj] * delta[id + 9];
-			Vw = action[id_obj] * delta[id + 10];
-
-			int x_neb, y_neb, k_neb;
-			for (int i = 1; i < 9; i++)
-			{
-				x_neb = x + e[i][0];
-				y_neb = y + e[i][1];
-				k_neb = x_neb + y_neb * NX;
-				if (flag[k_neb] & FLUID)
-				{
-					LBtype q = delta[id + i];
-					int k_neb2 = (y + 2 * e[i][1]) * NX + (x + 2 * e[i][0]);
-					LBtype temp = 6 * w[i] * (e[i][0] * Uw + e[i][1] * Vw);
-					f_temp[k_neb + i * totalCells] = (q * f_temp[k + opp[i] * totalCells] \
-						+ (1 - q) * f_temp[k_neb + opp[i] * totalCells] \
-						+ q * f_temp[k_neb2 + i * totalCells] + temp) / (1 + q);
-					f_temp[k + i * totalCells] = temp * Uw;
-					k_neb2 = (y - e[i][1]) * NX + (x - e[i][0]);
-					f_temp[k_neb2 + i * totalCells] = temp * Vw;
-
-					temp = f_temp[k_neb + i * totalCells] + f_temp[k + opp[i] * totalCells];
-					k_neb2 = (y - e[i][1]) * NX + (x - e[i][0]);
-					atomicAdd(&obs_share[DIM * id_obj], -temp * e[i][0] + f_temp[k + i * totalCells]);
-					atomicAdd(&obs_share[DIM * id_obj + 1], -temp * e[i][1] + f_temp[k_neb2 + i * totalCells]);
-				}
-			}
-		}
-		if (flag[k] & SENSOR)
-		{
-			LBtype u, v;
-			u = (g[1]+g[5]+g[8]-g[3]-g[6]-g[7])/RHO;
-			v = (g[2]+g[5]+g[6]-g[4]-g[7]-g[8])/RHO;
-			atomicAdd(&obs_share[DIM * id], u);
-			atomicAdd(&obs_share[DIM * id + 1], v);
-		}
-
-		__syncthreads();
-
-		for (int i = threadIdx.x; i < N_OBJS * DIM; i+=NT)
-		{
-			atomicAdd(&obs[i], obs_share[i]);
-		}
-	}
-
-	__global__ void InitTubeFlow(uint8_t *flag, LBtype *f)
-	{
-		__shared__ LBtype f_share[NT * NQ];
-		__shared__ uint8_t flag_share[NT];
-		int x, y, k;
-		LBtype u;
-		Index_lattice(x, y, k);
-		int totalCells = NX * NY;
-
-		flag_share[threadIdx.x] = flag[k];
-		for (int i = 0; i < NQ; i++)
-		{
-			f_share[threadIdx.x + i * NT] = f[k + i * totalCells];
-		}
-
-		__syncthreads();
-
-		u = U0 * 1.5 * (1 - 4 * (y - 0.5 * (NY - 1)) * (y - 0.5 * (NY - 1)) / ((NY - 2) * (NY - 2)));
-		if (y == 0 || y == NY - 1 || x == 0 || x == NX - 1)
-		{
-			flag_share[threadIdx.x] = SOLID;
-			for (int i = 0; i < NQ; i++)
-			{
-				f_share[threadIdx.x + i * NT] = 0;
-			}
-		}
-		else
-		{
-			flag_share[threadIdx.x] = FLUID;
-			for (int i = 0; i < NQ; i++)
-			{
-				f_share[threadIdx.x + i * NT] = w[i] * RHO * (3 * e[i][0] * u + \
-					4.5 * e[i][0] * e[i][0] * u * u - 1.5 * u * u);
-			}
-		}
-
-		__syncthreads();
-
-		flag[k] = flag_share[threadIdx.x];
-		for (int i = 0; i < NQ; i++)
-		{
-			f[k + i * totalCells] = f_share[threadIdx.x + i * NT];
-		}
-	}
-
-	// __global__ void AddVortex(LBtype *f, int32_t *config)
-	// {
-	// 	__shared__ LBtype f_share[NT * NQ];
-	// 	int x, y, k;
-	// 	LBtype u, v, u_vor, v_vor;
-	// 	Index_lattice(x, y, k);
-	// 	int totalCells = NX * NY;
-
-	// 	for (int i = 0; i < NQ; i++)
-	// 	{
-	// 		f_share[threadIdx.x + i * NT] = f[k + i * totalCells];
-	// 	}
-
-	// 	__syncthreads();
-
-	// 	u = f_share[threadIdx.x + 1 * NT] - f_share[threadIdx.x + 3 * NT] + f_share[threadIdx.x + 5 * NT] - f_share[threadIdx.x + 6 * NT] - f_share[threadIdx.x + 7 * NT] + f_share[threadIdx.x + 8 * NT];
-	// 	v = f_share[threadIdx.x + 2 * NT] - f_share[threadIdx.x + 4 * NT] + f_share[threadIdx.x + 5 * NT] + f_share[threadIdx.x + 6 * NT] - f_share[threadIdx.x + 7 * NT] - f_share[threadIdx.x + 8 * NT];
-
-	// 	if type & V_TAYLOR
-	// 	{
-	// 		u_vor = -2 * PI * U0 * sin(2 * PI * x / NX) * sin(2 * PI * y / NY);
-	// 		v_vor = 2 * PI * U0 * cos(2 * PI * x / NX) * cos(2 * PI * y / NY);
-	// 	}
-	// 	else
-	// 	{
-	// 		u_vor = 0;
-	// 		v_vor = 0;
-	// 	}
-
-
-	// }
-}
--- a/src/CelerisLab/lbm/kernels/kernel_v2.cu
+++ b/src/CelerisLab/lbm/kernels/kernel_v2.cu
@ -20,7 +20,7 @@
 // ---------------------------------------------------------------------------
 // Layer 0: Configuration (compile-time)
 // ---------------------------------------------------------------------------
-#include "macros.h"
+#include "config.h"

 // ---------------------------------------------------------------------------
 // Layer 1: Core primitives
@ -138,6 +138,8 @@ __global__ void InitTubeFlow_v2(uint8_t* flag, fpxx* fi)

 // ----- Main step (double-buffer) -----
 // Signature compatible with driver.py: flag, fi_in, fi_out, indx, delta, action, obs
+// TODO(Phase5b): Extract shared __device__ body to eliminate duplication with
+// StreamCollideDouble in step/one_step_double.cu (~140 lines overlap).
 __global__ void OneStep(
    uint8_t*  flag,
    fpxx*     fi_in,
@ -329,4 +331,252 @@ __global__ void OneStep(
 #endif
 }

+// ---------------------------------------------------------------------------
+// Esoteric-Pull wrappers (exported as extern "C" for PyCUDA)
+// ---------------------------------------------------------------------------
+
+// Thin forwarder: StreamCollideEsoPull is already __global__ in
+// step/one_step_esopull.cu.  We just need to ensure its symbol is
+// inside this extern "C" block so PyCUDA can resolve it by name.
+// Because it is already defined as __global__ in a header included above,
+// we cannot redefine it here.  Instead we provide a dedicated wrapper.
+__global__ void EsoPullStep(
+    fpxx*          fi,
+    uint8_t*       flag,
+    int32_t*       indx,
+    float*         delta,
+    float*         action,
+    float*         obs,
+    unsigned long  t)
+{
+    // Delegate to StreamCollideEsoPull with NULL rho/u/force arrays
+    // (we don't need macroscopic output during stepping).
+    // We cannot call a __global__ from another __global__, so we inline
+    // the body directly.  However since StreamCollideEsoPull is included
+    // as a full kernel above, we use a device-function extraction approach.
+    // For simplicity and correctness, we re-invoke the thread mapping here.
+
+#if DIM == 2
+    unsigned int x, y;
+    unsigned long k;
+    index_from_thread(x, y, k);
+    if (x >= (unsigned int)NX || y >= (unsigned int)NY) return;
+#elif DIM == 3
+    unsigned int x, y, z;
+    unsigned long k;
+    index_from_thread(x, y, z, k);
+    if (x >= (unsigned int)NX || y >= (unsigned int)NY || z >= (unsigned int)NZ) return;
+#endif
+
+    uint8_t fl = flag[k];
+
+    unsigned long j[NQ];
+    compute_neighbors(k, j);
+
+    float f[NQ];
+    load_f_esopull(k, f, fi, j, t);
+
+    // Solid nodes: BB + store, then return (essential for EsoPull correctness)
+    if ((fl & LEGACY_SOLID) && !(fl & LEGACY_INTERFACE) && !(fl & LEGACY_SENSOR)) {
+        if (x != 0 && x != (unsigned int)(NX - 1)) {
+            #pragma unroll
+            for (int i = 1; i < NQ; i += 2) {
+                float tmp = f[i]; f[i] = f[i+1]; f[i+1] = tmp;
+            }
+            store_f_esopull(k, f, fi, j, t);
+            return;
+        }
+    }
+
+    float rho_n, ux, uy;
+#if NQ == 9
+    compute_rho_u(f, rho_n, ux, uy);
+
+    // Inlet / outlet
+    if (fl & LEGACY_SOLID) {
+        bool interior_y = (y > 0u) && (y < (unsigned int)(NY - 1));
+        if (x == 0 && interior_y) {
+            unsigned long k_neb = linear_index(x + 1u, y);
+            unsigned long j_neb[NQ];
+            compute_neighbors(k_neb, j_neb);
+            float f_neb[NQ];
+            load_f_esopull(k_neb, f_neb, fi, j_neb, t);
+            apply_parabolic_inlet(f, f_neb, (float)y);
+        }
+        else if (x == (unsigned int)(NX - 1) && interior_y) {
+            unsigned long k_neb = linear_index(x - 1u, y);
+            unsigned long j_neb[NQ];
+            compute_neighbors(k_neb, j_neb);
+            float f_neb[NQ];
+            load_f_esopull(k_neb, f_neb, fi, j_neb, t);
+            apply_pressure_outlet(f, f_neb, (float)y);
+        } else {
+            #pragma unroll
+            for (int i = 1; i < NQ; i += 2) {
+                float tmp = f[i]; f[i] = f[i+1]; f[i+1] = tmp;
+            }
+        }
+    }
+
+    if (fl & LEGACY_OBSTACLE) {
+        #pragma unroll
+        for (int i = 1; i < NQ; i += 2) {
+            float tmp = f[i]; f[i] = f[i+1]; f[i+1] = tmp;
+        }
+    }
+
+    // Collision (fluid only)
+    if (fl & LEGACY_FLUID) {
+        float feq[NQ], Fin[NQ];
+        compute_rho_u(f, rho_n, ux, uy);
+        compute_feq(rho_n, ux, uy, feq);
+        zero_forcing(Fin);
+        float omega_col = d_params.omega;
+    #if USE_LES
+        omega_col = compute_omega_smag(f, feq, rho_n, omega_col);
+    #endif
+        omega_col = fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, omega_col));
+
+#if COLLISION_MODEL == 0
+        collide_srt(f, feq, Fin, omega_col);
+#elif COLLISION_MODEL == 1
+        collide_trt(f, feq, Fin, omega_col);
+#elif COLLISION_MODEL == 2
+        collide_mrt(f, rho_n, ux, uy, Fin, omega_col);
+#endif
+    }
+#elif NQ == 19
+    float uz;
+    compute_rho_u(f, rho_n, ux, uy, uz);
+
+    if (fl & LEGACY_SOLID) {
+        bool interior_y = (y > 0u) && (y < (unsigned int)(NY - 1));
+        if (x == 0 && interior_y) {
+            unsigned long k_neb = linear_index(x + 1u, y, z);
+            unsigned long j_neb[NQ];
+            compute_neighbors(k_neb, j_neb);
+            float f_neb[NQ];
+            load_f_esopull(k_neb, f_neb, fi, j_neb, t);
+            apply_parabolic_inlet_3d(f, f_neb, (float)y);
+        }
+        else if (x == (unsigned int)(NX - 1) && interior_y) {
+            unsigned long k_neb = linear_index(x - 1u, y, z);
+            unsigned long j_neb[NQ];
+            compute_neighbors(k_neb, j_neb);
+            float f_neb[NQ];
+            load_f_esopull(k_neb, f_neb, fi, j_neb, t);
+            apply_pressure_outlet_3d(f, f_neb, (float)y);
+        } else {
+            #pragma unroll
+            for (int i = 1; i < NQ; i += 2) {
+                float tmp = f[i]; f[i] = f[i+1]; f[i+1] = tmp;
+            }
+        }
+    }
+
+    if (fl & LEGACY_OBSTACLE) {
+        #pragma unroll
+        for (int i = 1; i < NQ; i += 2) {
+            float tmp = f[i]; f[i] = f[i+1]; f[i+1] = tmp;
+        }
+    }
+
+    if (fl & LEGACY_FLUID) {
+        float feq[NQ], Fin[NQ];
+        compute_rho_u(f, rho_n, ux, uy, uz);
+        compute_feq(rho_n, ux, uy, uz, feq);
+        zero_forcing(Fin);
+        float omega_col = d_params.omega;
+#if USE_LES
+        omega_col = compute_omega_smag(f, feq, rho_n, omega_col);
+#endif
+        omega_col = fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, omega_col));
+#if COLLISION_MODEL == 0
+        collide_srt(f, feq, Fin, omega_col);
+#elif COLLISION_MODEL == 1
+        collide_trt(f, feq, Fin, omega_col);
+#elif COLLISION_MODEL == 2
+        collide_mrt(f, rho_n, ux, uy, uz, Fin, omega_col);
+#endif
+    }
+#endif
+
+    store_f_esopull(k, f, fi, j, t);
+
+    // Sensor
+    if (fl & LEGACY_SENSOR) {
+        int id_obj = indx[k];
+        atomicAdd(&obs[DIM * id_obj],     ux);
+        atomicAdd(&obs[DIM * id_obj + 1], uy);
+    }
+}
+
+// ----- Esoteric-Pull initialization -----
+// Writes equilibrium to ALL NQ slots directly (not through esoteric store).
+// This ensures both even/odd read patterns get valid data on the first step.
+// Reference: FluidX3D lbm.cpp — kernel_initialize writes all slots directly.
+__global__ void InitEsoPull(uint8_t* flag, fpxx* fi)
+{
+#if DIM == 2
+    unsigned int x, y;
+    unsigned long k;
+    index_from_thread(x, y, k);
+    if (x >= (unsigned int)NX || y >= (unsigned int)NY) return;
+
+    float feq[NQ];
+    if (y == 0 || y == NY - 1 || x == 0 || x == NX - 1) {
+        flag[k] = LEGACY_SOLID;
+        for (int i = 0; i < NQ; i++) {
+            feq[i] = d_w[i] * RHO;
+#if USE_DDF_SHIFTING
+            feq[i] -= d_w[i];
+#endif
+        }
+    } else {
+        flag[k] = (uint8_t)LEGACY_FLUID;
+        float u_init = inlet_target_u((float)y);
+        for (int i = 0; i < NQ; i++) {
+            float cu = (float)d_cx[i] * u_init;
+            feq[i] = d_w[i] * RHO * (1.0f + 3.0f*cu + 4.5f*cu*cu - 1.5f*u_init*u_init);
+#if USE_DDF_SHIFTING
+            feq[i] -= d_w[i];
+#endif
+        }
+    }
+    // Direct store to ALL slots (not esoteric pattern)
+    for (int i = 0; i < NQ; i++) {
+        store_ddf(fi, index_f(k, (unsigned int)i), feq[i]);
+    }
+#elif DIM == 3
+    unsigned int x, y, z;
+    unsigned long k;
+    index_from_thread(x, y, z, k);
+    if (x >= (unsigned int)NX || y >= (unsigned int)NY || z >= (unsigned int)NZ) return;
+
+    float feq[NQ];
+    if (y == 0 || y == NY - 1 || x == 0 || x == NX - 1) {
+        flag[k] = LEGACY_SOLID;
+        for (int i = 0; i < NQ; i++) {
+            feq[i] = d_w[i] * RHO;
+#if USE_DDF_SHIFTING
+            feq[i] -= d_w[i];
+#endif
+        }
+    } else {
+        flag[k] = (uint8_t)LEGACY_FLUID;
+        float u_init = inlet_target_u((float)y);
+        for (int i = 0; i < NQ; i++) {
+            float cu = (float)d_cx[i] * u_init;
+            feq[i] = d_w[i] * RHO * (1.0f + 3.0f*cu + 4.5f*cu*cu - 1.5f*u_init*u_init);
+#if USE_DDF_SHIFTING
+            feq[i] -= d_w[i];
+#endif
+        }
+    }
+    for (int i = 0; i < NQ; i++) {
+        store_ddf(fi, index_f(k, (unsigned int)i), feq[i]);
+    }
+#endif
+}
+
 }  // extern "C"
--- a/src/CelerisLab/lbm/kernels/preproc.cu
+++ b/src/CelerisLab/lbm/kernels/preproc.cu
@ -1,2 +0,0 @@
-#include "macros.h"
-#include "const.h"
--- a/src/CelerisLab/lbm/kernels/step/one_step_esopull.cu
+++ b/src/CelerisLab/lbm/kernels/step/one_step_esopull.cu
@ -44,12 +44,6 @@ __global__ void StreamCollideEsoPull(

    uint8_t fl = flag[k];

-    // Skip pure solid / gas
-    if ((fl & LEGACY_SOLID) && !(fl & LEGACY_INTERFACE) && !(fl & LEGACY_SENSOR)) {
-        // For inlet/outlet solid nodes, we still process below
-        if (x != 0 && x != (unsigned int)(NX - 1)) return;
-    }
-
    // ----- Neighbor indices -----
    unsigned long j[NQ];
    compute_neighbors(k, j);
@ -58,6 +52,22 @@ __global__ void StreamCollideEsoPull(
    float f[NQ];
    load_f_esopull(k, f, fi, j, t);

+    // ----- Solid / wall nodes: bounce-back then store (essential for EsoPull) -----
+    // Unlike double-buffer, EsoPull requires solid nodes to participate in
+    // store_f so that adjacent fluid nodes pull the correct reflected DDF.
+    // Reference: FluidX3D kernel.cpp apply_moving_boundaries()
+    if ((fl & LEGACY_SOLID) && !(fl & LEGACY_INTERFACE) && !(fl & LEGACY_SENSOR)) {
+        if (x != 0 && x != (unsigned int)(NX - 1)) {
+            // Pure wall / interior solid: bounce-back (swap pairs) then store
+            #pragma unroll
+            for (int i = 1; i < NQ; i += 2) {
+                float ttmp = f[i]; f[i] = f[i+1]; f[i+1] = ttmp;
+            }
+            store_f_esopull(k, f, fi, j, t);
+            return;
+        }
+    }
+
    // ----- Compute macroscopic quantities -----
    float rho_n, ux, uy;
 #if NQ == 9
@ -103,9 +113,9 @@ __global__ void StreamCollideEsoPull(
        }
    }

-    if (y == 1 || y == (unsigned int)(NY - 2)) {
-        apply_wall_bb_d2q9(y, f);
-    }
+    // Wall BB at y=1/NY-2 is NOT needed for EsoPull: solid wall nodes
+    // now do BB+store, so fluid nodes pull correctly reflected DDFs.
+
 #elif NQ == 19
    if (fl & LEGACY_SOLID) {
        bool interior_y = (y > 0u) && (y < (unsigned int)(NY - 1));
@ -139,9 +149,7 @@ __global__ void StreamCollideEsoPull(
        }
    }

-    if (y == 1 || y == (unsigned int)(NY - 2)) {
-        apply_wall_bb_d3q19_y(y, f);
-    }
+    // Wall BB at y=1/NY-2 NOT needed for EsoPull (solid nodes do BB+store).
 #endif

    // ----- Forcing -----
--- a/src/CelerisLab/lbm/stepper.py
+++ b/src/CelerisLab/lbm/stepper.py
@ -0,0 +1,79 @@
+# CelerisLab/lbm/stepper.py
+"""
+LBMStepper — time-advance driver for LBM kernels.
+
+Owns kernel function handles and manages the double-buffer swap.
+Does NOT return observations by default — callers use field.get_macroscopic()
+or body manager methods to pull data when they need it.
+"""
+
+import numpy as np
+import pycuda.driver as cuda
+
+
+class LBMStepper:
+    """Drive the LBM kernel forward in time."""
+
+    def __init__(self, field, module: cuda.Module, cfg):
+        self.field = field
+        self.module = module
+        self.cfg = cfg
+
+        # Kernel handles
+        self.step_fn = module.get_function("OneStep")
+        self.init_fn = module.get_function("InitTubeFlow_v2")
+
+        # Launch geometry
+        tpb = cfg.threads_per_block
+        self.block = (tpb, 1, 1)
+        if cfg.dim == 2:
+            self.grid = (cfg.nx // tpb, cfg.ny, 1)
+        else:
+            self.grid = (cfg.nx // tpb, cfg.ny, cfg.nz)
+
+        self._step_count = 0
+
+    # -- Initialization ------------------------------------------------------
+    def initialize(self):
+        """Run the init kernel to set up channel flow + flags."""
+        f = self.field
+        self.init_fn(
+            f.flag_gpu, f.ddf_gpu,
+            block=self.block, grid=self.grid,
+        )
+        # Copy init state to both buffers
+        cuda.memcpy_dtod(f.temp_gpu, f.ddf_gpu, f.ddf.nbytes)
+        # Sync host
+        cuda.memcpy_dtoh(f.flag, f.flag_gpu)
+        cuda.memcpy_dtoh(f.ddf, f.ddf_gpu)
+
+    # -- Stepping ------------------------------------------------------------
+    def step(self, n: int = 1, action_gpu=None, obs_gpu=None):
+        """Advance *n* time steps.
+
+        Optional action_gpu / obs_gpu are raw device pointers for
+        object interaction (passed through to the kernel).
+        """
+        f = self.field
+        # Provide dummy pointers if no objects
+        dummy = cuda.mem_alloc(4) if action_gpu is None else None
+        act = action_gpu or dummy
+        ob = obs_gpu or dummy
+
+        for _ in range(n):
+            self.step_fn(
+                f.flag_gpu, f.ddf_gpu, f.temp_gpu,
+                f.indx_gpu, f.delta_gpu,
+                act, ob,
+                block=self.block, grid=self.grid,
+            )
+            # Swap buffers
+            f.ddf_gpu, f.temp_gpu = f.temp_gpu, f.ddf_gpu
+            self._step_count += 1
+
+        if dummy is not None:
+            dummy.free()
+
+    @property
+    def step_count(self) -> int:
+        return self._step_count
--- a/src/CelerisLab/simulation.py
+++ b/src/CelerisLab/simulation.py
@ -0,0 +1,165 @@
+# CelerisLab/simulation.py
+"""
+Top-level orchestrator — assembles LBM field, stepper, body manager,
+and CUDA context into a single coherent simulation.
+
+Usage::
+
+    sim = Simulation("configs/config_lbm.json")
+    sim.add_cylinder((100, 50), radius=10)
+    sim.initialize()
+    sim.run(1000)
+    macro = sim.get_macroscopic()   # {"rho": ..., "ux": ..., "uy": ...}
+"""
+
+from typing import Dict, Optional, Tuple, Any
+
+import numpy as np
+import pycuda.driver as cuda
+
+from .config import LBMConfig, BodyConfig, load_lbm_config, load_body_config
+from .cuda.context import CudaContext
+from .cuda import compiler_v2 as compiler
+from .lbm.field import LBMField
+from .lbm.stepper import LBMStepper
+from .body.objects import Cylinder, Sensor, SimObject
+from .body.manager import ObjectManager
+
+
+class Simulation:
+    """High-level simulation handle.
+
+    LBM field/stepper and body manager live at the same level;
+    this class orchestrates them.
+    """
+
+    def __init__(self,
+                 lbm_config_path: Optional[str] = None,
+                 body_config_path: Optional[str] = None,
+                 device_id: int = 0):
+        # Load configs
+        self.lbm_cfg = load_lbm_config(lbm_config_path)
+        self.body_cfg = load_body_config(body_config_path)
+
+        # CUDA context
+        self.ctx = CudaContext(device_id)
+        arch = self.ctx.sm_arch
+        if self.lbm_cfg.compute_capability != "auto":
+            arch = f"sm_{''.join(self.lbm_cfg.compute_capability.split('.'))}"
+
+        # Compile kernel
+        compiler.generate_config(self.lbm_cfg, n_objects=0)
+        self._ptx_path = compiler.compile_kernel(arch=arch)
+        self._module = compiler.load_module(self._ptx_path)
+
+        # LBM field & stepper
+        self.field = LBMField(self.lbm_cfg, self._module)
+        self.stepper = LBMStepper(self.field, self._module, self.lbm_cfg)
+
+        # Body manager
+        self.bodies = ObjectManager(
+            self.lbm_cfg.nx, self.lbm_cfg.ny,
+            self.lbm_cfg.nq, self.lbm_cfg.dim,
+        )
+
+        self._initialized = False
+
+    # -- Object management ---------------------------------------------------
+    def add_cylinder(self, center: Tuple[float, float],
+                     radius: float) -> int:
+        obj = Cylinder(obj_id=-1, center=center, radius=radius)
+        return self.bodies.add(obj)
+
+    def add_sensor(self, center: Tuple[float, float],
+                   radius: float) -> int:
+        obj = Sensor(obj_id=-1, center=center, radius=radius)
+        return self.bodies.add(obj)
+
+    def add_object(self, obj: SimObject) -> int:
+        return self.bodies.add(obj)
+
+    # -- Compilation ---------------------------------------------------------
+    def recompile(self):
+        """Re-generate config headers and recompile kernel.
+
+        Call after changing compile-time parameters (collision model, etc.).
+        """
+        arch = self.ctx.sm_arch
+        compiler.generate_config(self.lbm_cfg, n_objects=self.bodies.count)
+        self._ptx_path = compiler.compile_kernel(arch=arch)
+        self._module = compiler.load_module(self._ptx_path)
+        # Reconnect field and stepper to new module
+        self.field.module = self._module
+        self.field._upload_params()
+        self.stepper = LBMStepper(
+            self.field, self._module, self.lbm_cfg,
+        )
+
+    # -- Initialization ------------------------------------------------------
+    def initialize(self):
+        """Initialize flow field and sync objects to GPU."""
+        # Recompile if objects were added after construction
+        if self.bodies.count > 0:
+            self.recompile()
+        self.stepper.initialize()
+        if self.bodies.count > 0:
+            self.bodies.sync_to_gpu(self.field)
+        self._initialized = True
+
+    # -- Stepping ------------------------------------------------------------
+    def run(self, steps: int):
+        """Advance simulation by *steps* time steps."""
+        if not self._initialized:
+            raise RuntimeError("Call initialize() first")
+        self.stepper.step(
+            steps,
+            action_gpu=self.bodies.action_gpu,
+            obs_gpu=self.bodies.obs_gpu,
+        )
+
+    def step(self, n: int = 1):
+        """Advance *n* steps (convenience for interactive use)."""
+        self.run(n)
+
+    # -- Data access ---------------------------------------------------------
+    def get_macroscopic(self) -> Dict[str, np.ndarray]:
+        """Download DDF and return rho, ux, uy [, uz]."""
+        return self.field.get_macroscopic()
+
+    def get_ddf(self) -> np.ndarray:
+        self.field.download_ddf()
+        return self.field.ddf.copy()
+
+    def get_flags(self) -> np.ndarray:
+        return self.field.flag.copy()
+
+    # -- Runtime parameter updates -------------------------------------------
+    def update_runtime_params(self, **kwargs):
+        """Update __constant__ d_params without recompiling.
+
+        Accepted: omega, omega_bulk, fx, fy, fz, rho_ref, u_inlet, n_objects.
+        """
+        self.field.update_params(**kwargs)
+
+    # -- Snapshots -----------------------------------------------------------
+    def snapshot(self):
+        self.field.snapshot()
+
+    def restore(self):
+        self.field.restore()
+
+    # -- Cleanup -------------------------------------------------------------
+    def close(self):
+        self.ctx.close()
+
+    def __enter__(self):
+        return self
+
+    def __exit__(self, *exc):
+        self.close()
+
+    def __del__(self):
+        try:
+            self.close()
+        except Exception:
+            pass
--- a/tests/test_high_re_validation.py
+++ b/tests/test_high_re_validation.py
@ -83,6 +83,19 @@ def lattice_weights(nq):
    raise ValueError(f"Unsupported nq={nq}")


+def inlet_target_profile_1d(ny, u0, inlet_profile):
+    if int(inlet_profile) == 0:
+        return np.full(ny, float(u0), dtype=np.float32)
+
+    # Mirror boundary/inlet_outlet.cuh::inlet_target_u for consistent diagnostics.
+    y = np.arange(ny, dtype=np.float32)
+    y_clamped = np.clip(y, 1.0, float(ny - 2))
+    H = max(float(ny - 2), 1.0)
+    eta = (y_clamped - 0.5) / H
+    shape = np.clip(4.0 * eta * (1.0 - eta), 0.0, None)
+    return (float(u0) * 1.5 * shape).astype(np.float32)
+
+
 def impose_rest_state_on_nonfluid(cfg, host_ddf):
    nq = cfg["nq"]
    nx, ny, nz = cfg["nx"], cfg["ny"], cfg["nz"]
@ -135,12 +148,7 @@ def compute_case_diagnostics(cfg, host_ddf):
        x_probe = 1
        line_mask = fluid[:, x_probe]
        line_u = ux[:, x_probe]
-        y = np.arange(ny, dtype=np.float32)
-        if int(cfg.get("inlet_profile", 0)) == 0:
-            target = np.full(ny, float(cfg["u0"]), dtype=np.float32)
-        else:
-            yy = (y - 0.5 * (ny - 1)) / (ny - 2.0)
-            target = float(cfg["u0"]) * 1.5 * (1.0 - 4.0 * yy * yy)
+        target = inlet_target_profile_1d(ny, cfg["u0"], cfg.get("inlet_profile", 0))

        if np.any(line_mask):
            diff = line_u[line_mask] - target[line_mask]
@ -166,13 +174,9 @@ def compute_case_diagnostics(cfg, host_ddf):
                col_u.append(float(np.mean(ux[1:ny - 1, xp][col_mask])))
                col_r.append(float(np.mean(rho[1:ny - 1, xp][col_mask])))

-        if int(cfg.get("inlet_profile", 0)) == 0:
-            u_target_mean = float(cfg["u0"])
-        else:
-            y_int = np.arange(1, ny - 1, dtype=np.float32)
-            yy_int = (y_int - 0.5 * (ny - 1)) / (ny - 2.0)
-            target_int = float(cfg["u0"]) * 1.5 * (1.0 - 4.0 * yy_int * yy_int)
-            u_target_mean = float(np.mean(target_int)) if target_int.size > 0 else float(cfg["u0"])
+        target_full = inlet_target_profile_1d(ny, cfg["u0"], cfg.get("inlet_profile", 0))
+        target_int = target_full[1:ny - 1]
+        u_target_mean = float(np.mean(target_int)) if target_int.size > 0 else float(cfg["u0"])

        if len(col_u) >= 4:
            col_u_arr = np.array(col_u, dtype=np.float64)
@ -273,13 +277,9 @@ def compute_case_diagnostics(cfg, host_ddf):
            col_u.append(float(np.mean(ux[:, 1:ny - 1, xp][col_mask])))
            col_r.append(float(np.mean(rho[:, 1:ny - 1, xp][col_mask])))

-    if int(cfg.get("inlet_profile", 0)) == 0:
-        u_target_mean = float(cfg["u0"])
-    else:
-        y_int = np.arange(1, ny - 1, dtype=np.float32)
-        yy_int = (y_int - 0.5 * (ny - 1)) / (ny - 2.0)
-        target_int = float(cfg["u0"]) * 1.5 * (1.0 - 4.0 * yy_int * yy_int)
-        u_target_mean = float(np.mean(target_int)) if target_int.size > 0 else float(cfg["u0"])
+    target_full = inlet_target_profile_1d(ny, cfg["u0"], cfg.get("inlet_profile", 0))
+    target_int = target_full[1:ny - 1]
+    u_target_mean = float(np.mean(target_int)) if target_int.size > 0 else float(cfg["u0"])

    if len(col_u) >= 4:
        col_u_arr = np.array(col_u, dtype=np.float64)
@ -545,42 +545,74 @@ def compute_vis_omega(reynolds, diameter, u0):
 def set_macros(nx, ny, nz, dim, nq, vis, u0, collision_model, use_les, les_cs,
               outlet_mode, outlet_backflow_clamp, outlet_blend_alpha,
               omega_collision_max, inlet_profile, trt_magic_param):
-    lines = compiler.read_lines(compiler.kernel_path("macros.h"))
-    defs = {
-        "MULT_GPU": "False",
-        "NT": 128,
-        "X_1U": nx,
-        "Y_1U": ny,
-        "Z_1U": nz,
-        "LBtype": "float",
-        "UX": 1,
-        "UY": 1,
-        "UZ": 1,
-        "NX": nx,
-        "NY": ny,
-        "NZ": nz,
-        "DIM": dim,
-        "NQ": nq,
-        "VIS": f"{vis:.10f}",
-        "RHO": "1.0",
-        "U0": u0,
-        "N_OBJS": 0,
-        "COLLISION_MODEL": collision_model,
-        "STREAMING_MODEL": 0,
-        "STORE_PRECISION": 0,
-        "USE_DDF_SHIFTING": 0,
-        "USE_LES": int(use_les),
-        "LES_CS": f"{les_cs:.6f}f",
-        "INLET_PROFILE": int(inlet_profile),
-        "OUTLET_MODE": int(outlet_mode),
-        "OUTLET_BACKFLOW_CLAMP": int(outlet_backflow_clamp),
-        "OUTLET_BLEND_ALPHA": f"{float(outlet_blend_alpha):.3f}f",
-        "OMEGA_COLLISION_MAX": f"{float(omega_collision_max):.3f}f",
-        "TRT_MAGIC_PARAM": f"{float(trt_magic_param):.6f}f",
-    }
-    for name, value in defs.items():
-        lines = compiler.modify_macro(lines, name, value)
-    compiler.write_lines(compiler.kernel_path("macros.h"), lines)
+    """Write kernel config headers (config/*.h) — kernel_v2.cu uses config.h, not macros.h."""
+    cfg_dir = os.path.join(compiler.kernel_path("config"), "")
+    os.makedirs(cfg_dir, exist_ok=True)
+
+    with open(compiler.kernel_path("config/config_grid.h"), "w") as f:
+        f.write(f"""\
+// AUTO-GENERATED by test_high_re_validation — DO NOT EDIT MANUALLY
+#ifndef CELERIS_CONFIG_GRID_H
+#define CELERIS_CONFIG_GRID_H
+#define NT   128
+#define MULT_GPU False
+#define NX  {nx}
+#define NY  {ny}
+#define NZ  {nz}
+#define DIM {dim}
+#define NQ  {nq}
+#endif
+""")
+
+    with open(compiler.kernel_path("config/config_physics.h"), "w") as f:
+        f.write(f"""\
+// AUTO-GENERATED by test_high_re_validation — DO NOT EDIT MANUALLY
+#ifndef CELERIS_CONFIG_PHYSICS_H
+#define CELERIS_CONFIG_PHYSICS_H
+#define LBtype float
+#define VIS {vis:.10f}
+#define RHO 1.0
+#define U0  {u0}
+#define PI  3.141592653589793238
+#define FLUID     0x01
+#define SOLID     0x02
+#define GAS       0x04
+#define INTERFACE 0x08
+#define SENSOR    0x10
+#define OBSTACLE  0x20
+#define V_TAYLOR  1
+#endif
+""")
+
+    with open(compiler.kernel_path("config/config_method.h"), "w") as f:
+        f.write(f"""\
+// AUTO-GENERATED by test_high_re_validation — DO NOT EDIT MANUALLY
+#ifndef CELERIS_CONFIG_METHOD_H
+#define CELERIS_CONFIG_METHOD_H
+#define COLLISION_MODEL  {collision_model}
+#define STREAMING_MODEL  0
+#define STORE_PRECISION  0
+#define USE_DDF_SHIFTING 0
+#define USE_LES {int(use_les)}
+#define LES_CS  {les_cs:.6f}f
+#define INLET_PROFILE        {int(inlet_profile)}
+#define OUTLET_MODE          {int(outlet_mode)}
+#define OUTLET_BLEND_ALPHA   {float(outlet_blend_alpha):.3f}f
+#define OUTLET_BACKFLOW_CLAMP {int(outlet_backflow_clamp)}
+#define OMEGA_COLLISION_MIN 0.01f
+#define OMEGA_COLLISION_MAX {float(omega_collision_max):.4f}f
+#define TRT_MAGIC_PARAM     {float(trt_magic_param):.6f}f
+#endif
+""")
+
+    with open(compiler.kernel_path("config/config_objects.h"), "w") as f:
+        f.write("""\
+// AUTO-GENERATED by test_high_re_validation — DO NOT EDIT MANUALLY
+#ifndef CELERIS_CONFIG_OBJECTS_H
+#define CELERIS_CONFIG_OBJECTS_H
+#define N_OBJS 0
+#endif
+""")


 def build_flags_2d(nx, ny, cx, cy, radius):
@ -921,8 +953,20 @@ def main():
    parser.add_argument("--matrix-steps3d", type=int, default=600)
    args = parser.parse_args()

-    macro_path = compiler.kernel_path("macros.h")
-    macro_backup = compiler.read_lines(macro_path)
+    # Backup config/*.h (kernel_v2.cu uses config.h, not macros.h)
+    cfg_files = [
+        compiler.kernel_path("config/config_grid.h"),
+        compiler.kernel_path("config/config_physics.h"),
+        compiler.kernel_path("config/config_method.h"),
+        compiler.kernel_path("config/config_objects.h"),
+    ]
+    cfg_backups = {}
+    for p in cfg_files:
+        try:
+            with open(p) as f:
+                cfg_backups[p] = f.read()
+        except FileNotFoundError:
+            cfg_backups[p] = None

    out_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "output")
    os.makedirs(out_dir, exist_ok=True)
@ -978,7 +1022,10 @@ def main():
        print(f"Pass rate: {n_pass}/{len(results)}")
        print(f"Saved: {out_json}")
    finally:
-        compiler.write_lines(macro_path, macro_backup)
+        for p, content in cfg_backups.items():
+            if content is not None:
+                with open(p, "w") as f:
+                    f.write(content)


 if __name__ == "__main__":
--- a/tests/test_stability_matrix.py
+++ b/tests/test_stability_matrix.py
@ -0,0 +1,755 @@
+#!/usr/bin/env python3
+"""
+Stability Matrix Test
+=====================
+Tests three collision models (SRT/TRT/MRT) at low and high Re (with/without LES),
+plus Esoteric-Pull streaming at low Re with SRT.
+
+Outputs:
+  - Flow-field images (velocity, vorticity, streamlines) for each case
+  - Diagnostic JSON with stability metrics
+  - EsoPull vs double-buffer comparison plots
+
+Usage:
+  python3 tests/test_stability_matrix.py [--device 0] [--steps 2000]
+"""
+
+import argparse
+import json
+import math
+import os
+import struct
+import sys
+import time
+
+sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "src"))
+
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import numpy as np
+import pycuda.driver as cuda
+
+from CelerisLab.cuda import compiler
+
+# ---------------------------------------------------------------------------
+# Constants
+# ---------------------------------------------------------------------------
+FLUID = 0x01
+SOLID = 0x02
+OBSTACLE = 0x20  # fixed: was 0x04
+
+COLLISION_NAMES = {0: "SRT", 1: "TRT", 2: "MRT"}
+
+
+# ---------------------------------------------------------------------------
+# Helpers
+# ---------------------------------------------------------------------------
+def compute_vis_omega(re, diameter, u0):
+    vis = u0 * diameter / re
+    omega = 1.0 / (3.0 * vis + 0.5)
+    return vis, omega
+
+
+def lattice_weights(nq):
+    if nq == 9:
+        return np.array([4/9] + [1/9]*4 + [1/36]*4, dtype=np.float32)
+    if nq == 19:
+        return np.array([1/3] + [1/18]*6 + [1/36]*12, dtype=np.float32)
+    raise ValueError(f"nq={nq}")
+
+
+def build_flags_2d(nx, ny, cx, cy, radius):
+    flag = np.ones(nx * ny, dtype=np.uint8) * FLUID
+    for y in range(ny):
+        for x in range(nx):
+            k = y * nx + x
+            if y == 0 or y == ny - 1 or x == 0 or x == nx - 1:
+                flag[k] = SOLID
+            elif (x - cx)**2 + (y - cy)**2 < radius**2:
+                flag[k] = OBSTACLE
+    return flag
+
+
+def set_macros(nx, ny, dim, nq, vis, u0, collision_model, use_les, streaming_model,
+               omega_collision_max=1.999, inlet_profile=1, trt_magic_param=0.1875,
+               les_cs=0.16):
+    """Write config/*.h files used by kernel_v2.cu."""
+    cfg_dir = os.path.join(os.path.dirname(compiler.kernel_path("config.h")), "config")
+
+    # config_grid.h
+    with open(os.path.join(cfg_dir, "config_grid.h"), "w") as f:
+        f.write(f"""\
+// AUTO-GENERATED by test_stability_matrix.py
+#ifndef CELERIS_CONFIG_GRID_H
+#define CELERIS_CONFIG_GRID_H
+#define NT   128
+#define MULT_GPU 0
+#define NX  {nx}
+#define NY  {ny}
+#define NZ  1
+#define DIM {dim}
+#define NQ  {nq}
+#endif
+""")
+
+    # config_physics.h
+    with open(os.path.join(cfg_dir, "config_physics.h"), "w") as f:
+        f.write(f"""\
+// AUTO-GENERATED by test_stability_matrix.py
+#ifndef CELERIS_CONFIG_PHYSICS_H
+#define CELERIS_CONFIG_PHYSICS_H
+#define LBtype float
+#define VIS {vis:.10f}
+#define RHO 1.0
+#define U0  {u0}
+#define PI 3.141592653589793238
+#define FLUID       0x01
+#define SOLID       0x02
+#define GAS         0x04
+#define INTERFACE   0x08
+#define SENSOR      0x10
+#define OBSTACLE    0x20
+#define V_TAYLOR    1
+#endif
+""")
+
+    # config_method.h
+    with open(os.path.join(cfg_dir, "config_method.h"), "w") as f:
+        f.write(f"""\
+// AUTO-GENERATED by test_stability_matrix.py
+#ifndef CELERIS_CONFIG_METHOD_H
+#define CELERIS_CONFIG_METHOD_H
+#define COLLISION_MODEL {collision_model}
+#define STREAMING_MODEL {streaming_model}
+#define STORE_PRECISION 0
+#define USE_DDF_SHIFTING 0
+#define USE_LES {int(use_les)}
+#define LES_CS {les_cs:.6f}f
+#define INLET_PROFILE {int(inlet_profile)}
+#define OUTLET_MODE 0
+#define OUTLET_BLEND_ALPHA 0.700f
+#define OUTLET_BACKFLOW_CLAMP 1
+#define OMEGA_COLLISION_MIN 0.01f
+#define OMEGA_COLLISION_MAX {float(omega_collision_max):.3f}f
+#define TRT_MAGIC_PARAM {float(trt_magic_param):.6f}f
+#endif
+""")
+
+    # config_objects.h
+    with open(os.path.join(cfg_dir, "config_objects.h"), "w") as f:
+        f.write("""\
+// AUTO-GENERATED by test_stability_matrix.py
+#ifndef CELERIS_CONFIG_OBJECTS_H
+#define CELERIS_CONFIG_OBJECTS_H
+#define N_OBJS 0
+#endif
+""")
+
+
+def pack_d_params(nx, ny, omega, u0):
+    """Pack LBMParams struct for __constant__ memory upload."""
+    return struct.pack(
+        "IIIQfffffffI",
+        nx, ny, 1,           # Nx, Ny, Nz
+        nx * ny,              # N
+        omega,                # omega
+        1.1,                  # omega_bulk
+        0.0, 0.0, 0.0,       # fx, fy, fz
+        1.0,                  # rho_ref
+        u0,                   # u_inlet
+        0,                    # n_objects
+    )
+
+
+def impose_rest_on_nonfluid(flag, host_ddf, nq, nx, ny):
+    w = lattice_weights(nq)
+    f = host_ddf.reshape(nq, ny, nx)
+    nonfluid = flag.reshape(ny, nx) != FLUID
+    for i in range(nq):
+        f[i, nonfluid] = w[i]
+    return host_ddf
+
+
+def compute_macros_2d(host_ddf, nq, nx, ny, flag):
+    """Compute rho, ux, uy from DDF."""
+    cx9 = [0, 1, -1, 0, 0, 1, -1, 1, -1]
+    cy9 = [0, 0, 0, 1, -1, 1, -1, -1, 1]
+    f = host_ddf.reshape(nq, ny, nx)
+    rho = np.sum(f, axis=0)
+    ux = np.zeros_like(rho)
+    uy = np.zeros_like(rho)
+    for i in range(nq):
+        ux += cx9[i] * f[i]
+        uy += cy9[i] * f[i]
+    rho_safe = np.where(np.abs(rho) > 1e-12, rho, 1.0)
+    ux /= rho_safe
+    uy /= rho_safe
+    return rho, ux, uy
+
+
+def diagnose(rho, ux, uy, flag, nx, ny):
+    """Compute stability diagnostics."""
+    fluid = flag.reshape(ny, nx) == FLUID
+    nan_count = int(np.isnan(rho).sum())
+    rho_min = float(np.nanmin(rho))
+    rho_max = float(np.nanmax(rho))
+    mass = float(np.nansum(rho[fluid]))
+    vel = np.sqrt(ux**2 + uy**2)
+
+    # Ma check
+    ma_max = float(np.nanmax(vel[fluid])) * math.sqrt(3.0) if np.any(fluid) else 0.0
+
+    # Vorticity RMS in wake region
+    vort = np.gradient(uy, axis=1) - np.gradient(ux, axis=0)
+    wake_mask = fluid & (np.arange(nx)[None, :] > nx // 3)
+    vort_rms = float(np.sqrt(np.nanmean(vort[wake_mask]**2))) if np.any(wake_mask) else 0.0
+
+    stable = nan_count == 0 and rho_min > 0.0 and rho_max < 2.0
+    return {
+        "nan_count": nan_count,
+        "rho_min": rho_min,
+        "rho_max": rho_max,
+        "mass": mass,
+        "ma_max": ma_max,
+        "vort_rms": vort_rms,
+        "stable": stable,
+    }
+
+
+def plot_flow(rho, ux, uy, flag, nx, ny, title, out_path):
+    """Plot velocity magnitude, vorticity, and streamlines."""
+    fluid_mask = flag.reshape(ny, nx) != FLUID
+    vel = np.sqrt(ux**2 + uy**2)
+    vel_m = np.ma.array(vel, mask=fluid_mask)
+    vort = np.gradient(uy, axis=1) - np.gradient(ux, axis=0)
+    vort_m = np.ma.array(vort, mask=fluid_mask)
+
+    fig, axes = plt.subplots(1, 3, figsize=(18, 5))
+
+    # Velocity magnitude
+    im0 = axes[0].imshow(vel_m, origin="lower", aspect="auto", cmap="turbo")
+    plt.colorbar(im0, ax=axes[0], label="|u|")
+    axes[0].set_title("Velocity Magnitude")
+
+    # Vorticity
+    vals = vort[~fluid_mask]
+    if vals.size > 0:
+        vmax = max(float(np.percentile(np.abs(vals), 99)), 1e-8)
+    else:
+        vmax = 1e-6
+    im1 = axes[1].imshow(vort_m, origin="lower", aspect="auto", cmap="RdBu_r",
+                          vmin=-vmax, vmax=vmax)
+    plt.colorbar(im1, ax=axes[1], label="vorticity")
+    axes[1].set_title("Vorticity")
+
+    # Streamlines
+    X, Y = np.meshgrid(np.arange(nx), np.arange(ny))
+    ux_s = np.ma.array(ux, mask=fluid_mask)
+    uy_s = np.ma.array(uy, mask=fluid_mask)
+    speed = np.ma.sqrt(ux_s**2 + uy_s**2)
+    axes[2].streamplot(X, Y, ux_s, uy_s, color=speed, cmap="viridis",
+                       density=2.0, linewidth=0.7)
+    axes[2].set_xlim(0, nx)
+    axes[2].set_ylim(0, ny)
+    axes[2].set_title("Streamlines")
+
+    fig.suptitle(title, fontsize=13)
+    fig.tight_layout()
+    fig.savefig(out_path, dpi=150)
+    plt.close(fig)
+    return out_path
+
+
+# ---------------------------------------------------------------------------
+# Case runner: double-buffer
+# ---------------------------------------------------------------------------
+def run_double_buffer(device_id, cfg, out_dir):
+    """Run a case with standard double-buffer streaming."""
+    nx, ny = cfg["nx"], cfg["ny"]
+    nq = cfg["nq"]
+    n = nx * ny
+
+    set_macros(nx, ny, cfg["dim"], nq, cfg["vis"], cfg["u0"],
+               cfg["collision_model"], cfg["use_les"], streaming_model=0,
+               omega_collision_max=cfg.get("omega_max", 1.999),
+               trt_magic_param=cfg.get("trt_magic", 0.1875))
+    compiler.compile_kernel_v2()
+
+    cuda.init()
+    dev = cuda.Device(device_id)
+    ctx = dev.make_context()
+    try:
+        mod = cuda.module_from_file(compiler.kernel_path("kernel_v2.ptx"))
+        init_fn = mod.get_function("InitTubeFlow_v2")
+        step_fn = mod.get_function("OneStep")
+
+        # Upload d_params
+        params_ptr, params_size = mod.get_global("d_params")
+        params_data = pack_d_params(nx, ny, cfg["omega"], cfg["u0"])
+        if len(params_data) < params_size:
+            params_data += b"\x00" * (params_size - len(params_data))
+        cuda.memcpy_htod(params_ptr, params_data)
+
+        fsize = n * nq * 4
+        d_fi = cuda.mem_alloc(fsize)
+        d_fi2 = cuda.mem_alloc(fsize)
+        d_flag = cuda.mem_alloc(n)
+        d_indx = cuda.mem_alloc(n * 4)
+        d_delta = cuda.mem_alloc(4)
+        d_action = cuda.mem_alloc(4)
+        d_obs = cuda.mem_alloc(4)
+        cuda.memset_d32(d_indx, 0, n)
+        cuda.memset_d32(d_delta, 0, 1)
+        cuda.memset_d32(d_action, 0, 1)
+        cuda.memset_d32(d_obs, 0, 1)
+
+        block = (128, 1, 1)
+        grid = ((nx + 127) // 128, ny, 1)
+
+        init_fn(d_flag, d_fi, block=block, grid=grid)
+        cuda.memcpy_htod(d_flag, cfg["flag"])
+
+        host0 = np.empty(n * nq, dtype=np.float32)
+        cuda.memcpy_dtoh(host0, d_fi)
+        host0 = impose_rest_on_nonfluid(cfg["flag"], host0, nq, nx, ny)
+        cuda.memcpy_htod(d_fi, host0)
+        cuda.memcpy_htod(d_fi2, host0)
+
+        steps = cfg["steps"]
+        report = max(steps // 5, 1)
+        t0 = time.time()
+        diverged_step = None
+        for s in range(steps):
+            step_fn(d_flag, d_fi, d_fi2, d_indx, d_delta, d_action, d_obs,
+                    block=block, grid=grid)
+            d_fi, d_fi2 = d_fi2, d_fi
+
+            if (s + 1) % report == 0:
+                cuda.Context.synchronize()
+                h = np.empty(n * nq, dtype=np.float32)
+                cuda.memcpy_dtoh(h, d_fi)
+                rho_c = h.reshape(nq, ny, nx).sum(axis=0)
+                nc = int(np.isnan(rho_c).sum())
+                center = float(rho_c[ny // 2, nx // 2])
+                print(f"    step {s+1:6d}: rho_center={center:.6f} nan={nc}")
+                if nc > 0:
+                    diverged_step = s + 1
+                    break
+
+        cuda.Context.synchronize()
+        elapsed = time.time() - t0
+        host = np.empty(n * nq, dtype=np.float32)
+        cuda.memcpy_dtoh(host, d_fi)
+
+        rho, ux, uy = compute_macros_2d(host, nq, nx, ny, cfg["flag"])
+        diag = diagnose(rho, ux, uy, cfg["flag"], nx, ny)
+        diag["elapsed"] = elapsed
+        diag["mlups"] = n * steps / elapsed / 1e6 if elapsed > 0 else 0
+        diag["diverged_step"] = diverged_step
+
+        tag = cfg["tag"]
+        plot_path = plot_flow(rho, ux, uy, cfg["flag"], nx, ny, tag,
+                              os.path.join(out_dir, f"{tag}.png"))
+        diag["plot"] = plot_path
+        return diag
+    finally:
+        ctx.pop()
+
+
+# ---------------------------------------------------------------------------
+# Case runner: Esoteric-Pull (single buffer)
+# ---------------------------------------------------------------------------
+def run_esopull(device_id, cfg, out_dir):
+    """Run a case with Esoteric-Pull single-buffer streaming."""
+    nx, ny = cfg["nx"], cfg["ny"]
+    nq = cfg["nq"]
+    n = nx * ny
+
+    set_macros(nx, ny, cfg["dim"], nq, cfg["vis"], cfg["u0"],
+               cfg["collision_model"], cfg["use_les"], streaming_model=1,
+               omega_collision_max=cfg.get("omega_max", 1.999),
+               trt_magic_param=cfg.get("trt_magic", 0.1875))
+    compiler.compile_kernel_v2()
+
+    cuda.init()
+    dev = cuda.Device(device_id)
+    ctx = dev.make_context()
+    try:
+        mod = cuda.module_from_file(compiler.kernel_path("kernel_v2.ptx"))
+        init_fn = mod.get_function("InitEsoPull")
+        step_fn = mod.get_function("EsoPullStep")
+
+        # Upload d_params
+        params_ptr, params_size = mod.get_global("d_params")
+        params_data = pack_d_params(nx, ny, cfg["omega"], cfg["u0"])
+        if len(params_data) < params_size:
+            params_data += b"\x00" * (params_size - len(params_data))
+        cuda.memcpy_htod(params_ptr, params_data)
+
+        fsize = n * nq * 4
+        d_fi = cuda.mem_alloc(fsize)
+        d_flag = cuda.mem_alloc(n)
+        d_indx = cuda.mem_alloc(n * 4)
+        d_delta = cuda.mem_alloc(4)
+        d_action = cuda.mem_alloc(4)
+        d_obs = cuda.mem_alloc(4)
+        cuda.memset_d32(d_indx, 0, n)
+        cuda.memset_d32(d_delta, 0, 1)
+        cuda.memset_d32(d_action, 0, 1)
+        cuda.memset_d32(d_obs, 0, 1)
+
+        block = (128, 1, 1)
+        grid = ((nx + 127) // 128, ny, 1)
+
+        init_fn(d_flag, d_fi, block=block, grid=grid)
+        cuda.memcpy_htod(d_flag, cfg["flag"])
+
+        # Note: for EsoPull, we don't impose_rest_on_nonfluid on the raw
+        # DDF because the data is stored in esoteric layout.  InitEsoPull
+        # already stores rest equilibrium for solid nodes.
+
+        steps = cfg["steps"]
+        report = max(steps // 5, 1)
+        t0 = time.time()
+        diverged_step = None
+
+        for s in range(steps):
+            t_val = np.uint64(s)  # timestep counter for load/store parity
+            step_fn(d_fi, d_flag, d_indx, d_delta, d_action, d_obs,
+                    t_val, block=block, grid=grid)
+
+            if (s + 1) % report == 0:
+                cuda.Context.synchronize()
+                # For diagnostics, download raw DDF and decode from esopull layout
+                h = np.empty(n * nq, dtype=np.float32)
+                cuda.memcpy_dtoh(h, d_fi)
+                # Esoteric layout: at this point the DDF is in post-store layout
+                # for timestep s.  To compute macros we need to "undo" the esoteric
+                # read pattern.  A simpler approach: compute rho = sum(fi) per node.
+                # Because sum is invariant under slot permutation, rho is correct.
+                # But ux/uy need correct direction assignment.
+                # For diagnostic, use a simple sum-based stability check.
+                f_arr = h.reshape(nq, ny, nx)
+                rho_c = f_arr.sum(axis=0)
+                nc = int(np.isnan(rho_c).sum())
+                center = float(rho_c[ny // 2, nx // 2])
+                print(f"    step {s+1:6d}: rho_center={center:.6f} nan={nc}")
+                if nc > 0:
+                    diverged_step = s + 1
+                    break
+
+        cuda.Context.synchronize()
+        elapsed = time.time() - t0
+
+        # For final macros, do one more step that also writes to rho/u arrays.
+        # But we don't have UpdateMacro for EsoPull yet.  Instead, use the
+        # approach: run a "read-only" macro computation from the esoteric layout.
+        # For correctness, we load from the proper esoteric positions on host.
+        h = np.empty(n * nq, dtype=np.float32)
+        cuda.memcpy_dtoh(h, d_fi)
+
+        rho, ux, uy = _decode_esopull_macros(h, nq, nx, ny, cfg["flag"], steps)
+        diag = diagnose(rho, ux, uy, cfg["flag"], nx, ny)
+        diag["elapsed"] = elapsed
+        diag["mlups"] = n * steps / elapsed / 1e6 if elapsed > 0 else 0
+        diag["diverged_step"] = diverged_step
+
+        tag = cfg["tag"]
+        plot_path = plot_flow(rho, ux, uy, cfg["flag"], nx, ny, tag,
+                              os.path.join(out_dir, f"{tag}.png"))
+        diag["plot"] = plot_path
+        return diag
+    finally:
+        ctx.pop()
+
+
+def _decode_esopull_macros(host_ddf, nq, nx, ny, flag, last_t):
+    """Decode macroscopic quantities from esoteric-pull layout on host.
+
+    After step t (0-based), the store was done at parity t.
+    The next load would use parity t+1.  To read correct DDFs we mimic
+    load_f_esopull at t_read = last_t (the parity of the *next* step to execute).
+    """
+    fi = host_ddf.reshape(nq, ny * nx)  # fi[direction, node]
+    t_read = last_t  # parity for the load that would happen next
+
+    cx9 = np.array([0, 1, -1, 0, 0, 1, -1, 1, -1], dtype=np.float32)
+    cy9 = np.array([0, 0, 0, 1, -1, 1, -1, -1, 1], dtype=np.float32)
+
+    # Compute neighbor table once
+    j_table = np.zeros((nq, ny * nx), dtype=np.int64)
+    for y in range(ny):
+        for x in range(nx):
+            k = y * nx + x
+            xp = (x + 1) % nx
+            xm = (x - 1) % nx
+            yp = (y + 1) % ny
+            ym = (y - 1) % ny
+            j_table[0, k] = k
+            j_table[1, k] = yp * nx + xp if nq > 1 else k  # placeholder
+            j_table[2, k] = ym * nx + xm if nq > 2 else k
+            # D2Q9 neighbors: j[i] = neighbor in direction c_i
+            if nq == 9:
+                j_table[1, k] = y * nx + xp   # +x
+                j_table[2, k] = y * nx + xm   # -x
+                j_table[3, k] = yp * nx + x   # +y
+                j_table[4, k] = ym * nx + x   # -y
+                j_table[5, k] = yp * nx + xp  # +x+y
+                j_table[6, k] = ym * nx + xm  # -x-y
+                j_table[7, k] = ym * nx + xp  # +x-y
+                j_table[8, k] = yp * nx + xm  # -x+y
+
+    n = nx * ny
+    f_decoded = np.zeros((nq, n), dtype=np.float32)
+    f_decoded[0] = fi[0]
+
+    for i in range(1, nq, 2):
+        if t_read & 1:
+            # Odd: f[i] from fi[n, i], f[i+1] from fi[j[i], i+1]
+            f_decoded[i] = fi[i]
+            f_decoded[i + 1] = fi[i + 1, j_table[i]]
+        else:
+            # Even: f[i] from fi[n, i+1], f[i+1] from fi[j[i], i]
+            f_decoded[i] = fi[i + 1]
+            f_decoded[i + 1] = fi[i, j_table[i]]
+
+    f_decoded = f_decoded.reshape(nq, ny, nx)
+    rho = f_decoded.sum(axis=0)
+    rho_safe = np.where(np.abs(rho) > 1e-12, rho, 1.0)
+    ux = np.zeros_like(rho)
+    uy = np.zeros_like(rho)
+    for i in range(nq):
+        ux += cx9[i] * f_decoded[i]
+        uy += cy9[i] * f_decoded[i]
+    ux /= rho_safe
+    uy /= rho_safe
+    return rho, ux, uy
+
+
+# ---------------------------------------------------------------------------
+# Case builders
+# ---------------------------------------------------------------------------
+def build_cases(steps_low, steps_high):
+    """Build the full test matrix."""
+    # Grid params (moderate size for fast testing)
+    nx, ny = 384, 192
+    cx_ob, cy_ob, radius = 96.0, 96.0, 18.0
+    u0 = 0.04
+
+    cases = []
+    for re_val, re_label, n_steps, use_les in [
+        (100.0,  "Re100",  steps_low,  False),
+        (100.0,  "Re100",  steps_low,  True),
+        (3000.0, "Re3000", steps_high, False),
+        (3000.0, "Re3000", steps_high, True),
+    ]:
+        for cm in (0, 1, 2):
+            vis, omega = compute_vis_omega(re_val, 2.0 * radius, u0)
+            les_tag = "LES" if use_les else "noLES"
+            cm_name = COLLISION_NAMES[cm]
+            tag = f"DB_{re_label}_{cm_name}_{les_tag}"
+            cases.append({
+                "tag": tag,
+                "nx": nx, "ny": ny,
+                "dim": 2, "nq": 9,
+                "cx": cx_ob, "cy": cy_ob, "radius": radius,
+                "flag": build_flags_2d(nx, ny, cx_ob, cy_ob, radius),
+                "u0": u0,
+                "vis": vis,
+                "omega": omega,
+                "collision_model": cm,
+                "use_les": use_les,
+                "steps": n_steps,
+                "streaming": "double_buffer",
+                "omega_max": 1.999,
+                "trt_magic": 0.1875,
+            })
+
+    # EsoPull case: low Re, SRT only
+    re_eso = 100.0
+    vis_eso, omega_eso = compute_vis_omega(re_eso, 2.0 * radius, u0)
+    cases.append({
+        "tag": "EsoPull_Re100_SRT_noLES",
+        "nx": nx, "ny": ny,
+        "dim": 2, "nq": 9,
+        "cx": cx_ob, "cy": cy_ob, "radius": radius,
+        "flag": build_flags_2d(nx, ny, cx_ob, cy_ob, radius),
+        "u0": u0,
+        "vis": vis_eso,
+        "omega": omega_eso,
+        "collision_model": 0,
+        "use_les": False,
+        "steps": steps_low,
+        "streaming": "esopull",
+        "omega_max": 1.999,
+        "trt_magic": 0.1875,
+    })
+
+    return cases
+
+
+# ---------------------------------------------------------------------------
+# Comparison plot: EsoPull vs DoubleBuffer
+# ---------------------------------------------------------------------------
+def plot_comparison(results, out_dir):
+    """Compare EsoPull and DoubleBuffer at matching Re/collision settings."""
+    eso_key = "EsoPull_Re100_SRT_noLES"
+    db_key = "DB_Re100_SRT_noLES"
+
+    eso = results.get(eso_key)
+    db = results.get(db_key)
+    if eso is None or db is None:
+        return None
+
+    fig, axes = plt.subplots(2, 3, figsize=(18, 10))
+    fig.suptitle("EsoPull vs DoubleBuffer — Re100 SRT noLES", fontsize=14)
+
+    labels = ["DoubleBuffer", "EsoPull"]
+    for row, (r, label) in enumerate([(db, labels[0]), (eso, labels[1])]):
+        vel_img = plt.imread(r["plot"]) if os.path.exists(r["plot"]) else None
+        if vel_img is not None:
+            axes[row, 0].imshow(vel_img)
+            axes[row, 0].set_title(f"{label}: flow field")
+            axes[row, 0].axis("off")
+        else:
+            axes[row, 0].text(0.5, 0.5, f"No image for {label}",
+                              ha="center", va="center", transform=axes[row, 0].transAxes)
+            axes[row, 0].set_title(label)
+
+        # Metrics bar chart
+        metrics = {
+            "rho_min": r.get("rho_min", 0),
+            "rho_max": r.get("rho_max", 0),
+            "ma_max": r.get("ma_max", 0),
+            "vort_rms": r.get("vort_rms", 0),
+        }
+        bars = list(metrics.keys())
+        vals = [float(metrics[b]) for b in bars]
+        axes[row, 1].barh(bars, vals, color=["steelblue", "salmon", "green", "purple"])
+        axes[row, 1].set_title(f"{label}: diagnostics")
+
+        # Stability text
+        text_lines = [
+            f"stable: {r.get('stable', '?')}",
+            f"nan_count: {r.get('nan_count', '?')}",
+            f"mass: {r.get('mass', 0):.2f}",
+            f"MLUPS: {r.get('mlups', 0):.1f}",
+            f"diverged_step: {r.get('diverged_step', 'None')}",
+        ]
+        axes[row, 2].text(0.1, 0.5, "\n".join(text_lines), fontsize=12,
+                          family="monospace", va="center",
+                          transform=axes[row, 2].transAxes)
+        axes[row, 2].set_title(f"{label}: summary")
+        axes[row, 2].axis("off")
+
+    fig.tight_layout()
+    cmp_path = os.path.join(out_dir, "esopull_vs_doublebuffer.png")
+    fig.savefig(cmp_path, dpi=150)
+    plt.close(fig)
+    return cmp_path
+
+
+# ---------------------------------------------------------------------------
+# Main
+# ---------------------------------------------------------------------------
+def main():
+    parser = argparse.ArgumentParser(description="Stability matrix test")
+    parser.add_argument("--device", type=int, default=0)
+    parser.add_argument("--steps-low", type=int, default=3000,
+                        help="Steps for low-Re cases")
+    parser.add_argument("--steps-high", type=int, default=6000,
+                        help="Steps for high-Re cases")
+    parser.add_argument("--only-esopull", action="store_true",
+                        help="Only run the EsoPull test")
+    args = parser.parse_args()
+
+    # Backup config/*.h files (kernel_v2.cu uses config.h, NOT macros.h)
+    cfg_dir = os.path.join(os.path.dirname(compiler.kernel_path("config.h")), "config")
+    config_files = ["config_grid.h", "config_physics.h", "config_method.h", "config_objects.h"]
+    config_backups = {}
+    for cf in config_files:
+        path = os.path.join(cfg_dir, cf)
+        with open(path, "r") as f:
+            config_backups[path] = f.read()
+
+    out_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)),
+                           "..", "output", "stability_matrix")
+    os.makedirs(out_dir, exist_ok=True)
+
+    cases = build_cases(args.steps_low, args.steps_high)
+    if args.only_esopull:
+        cases = [c for c in cases if c["streaming"] == "esopull"]
+
+    results = {}
+    try:
+        for i, cfg in enumerate(cases):
+            tag = cfg["tag"]
+            streaming = cfg["streaming"]
+            print(f"\n[{i+1}/{len(cases)}] {tag}")
+            print(f"  Re={cfg['u0']*2*cfg['radius']/cfg['vis']:.0f}, "
+                  f"omega={cfg['omega']:.4f}, "
+                  f"collision={COLLISION_NAMES[cfg['collision_model']]}, "
+                  f"LES={cfg['use_les']}, streaming={streaming}")
+
+            if streaming == "esopull":
+                diag = run_esopull(args.device, cfg, out_dir)
+            else:
+                diag = run_double_buffer(args.device, cfg, out_dir)
+
+            diag["tag"] = tag
+            diag["streaming"] = streaming
+            diag["collision"] = COLLISION_NAMES[cfg["collision_model"]]
+            diag["use_les"] = cfg["use_les"]
+            diag["re"] = cfg["u0"] * 2 * cfg["radius"] / cfg["vis"]
+            results[tag] = diag
+
+            status = "PASS" if diag["stable"] else "FAIL"
+            print(f"  => {status}: rho=[{diag['rho_min']:.4f}, {diag['rho_max']:.4f}], "
+                  f"nan={diag['nan_count']}, ma_max={diag['ma_max']:.4f}, "
+                  f"MLUPS={diag['mlups']:.1f}")
+
+        # Comparison plot
+        cmp_path = plot_comparison(results, out_dir)
+        if cmp_path:
+            print(f"\nComparison plot: {cmp_path}")
+
+        # Summary table
+        print("\n" + "=" * 100)
+        print(f"{'Tag':<35s} {'Stream':<8s} {'Col':<5s} {'LES':<5s} "
+              f"{'Re':>6s} {'Stable':>7s} {'rho_min':>9s} {'rho_max':>9s} "
+              f"{'Ma_max':>8s} {'MLUPS':>7s}")
+        print("-" * 100)
+        for tag, r in results.items():
+            print(f"{tag:<35s} {r['streaming']:<8s} {r['collision']:<5s} "
+                  f"{'Y' if r['use_les'] else 'N':<5s} "
+                  f"{r['re']:6.0f} {'PASS' if r['stable'] else 'FAIL':>7s} "
+                  f"{r['rho_min']:9.5f} {r['rho_max']:9.5f} "
+                  f"{r['ma_max']:8.5f} {r['mlups']:7.1f}")
+        print("=" * 100)
+
+        # Save JSON
+        json_path = os.path.join(out_dir, "stability_matrix_results.json")
+        json_results = {}
+        for k, v in results.items():
+            jr = {}
+            for rk, rv in v.items():
+                if isinstance(rv, (np.integer, np.floating)):
+                    jr[rk] = float(rv)
+                elif isinstance(rv, np.bool_):
+                    jr[rk] = bool(rv)
+                else:
+                    jr[rk] = rv
+            json_results[k] = jr
+        with open(json_path, "w") as f:
+            json.dump(json_results, f, indent=2)
+        print(f"\nResults saved: {json_path}")
+
+    finally:
+        for path, content in config_backups.items():
+            with open(path, "w") as f:
+                f.write(content)
+
+
+if __name__ == "__main__":
+    main()