diff --git a/README.md b/README.md index 17751da..1ebe5a5 100644 --- 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. diff --git a/legacy/README.md b/legacy/README.md new file mode 100644 index 0000000..6fbca5d --- /dev/null +++ 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`. diff --git a/src/CelerisLab/common/utils.py b/legacy/common_utils.py similarity index 100% rename from src/CelerisLab/common/utils.py rename to legacy/common_utils.py diff --git a/src/CelerisLab/cuda/compiler.py b/legacy/cuda_compiler_v1.py similarity index 100% rename from src/CelerisLab/cuda/compiler.py rename to legacy/cuda_compiler_v1.py diff --git a/src/CelerisLab/lbm/configs/config_cuda.json b/legacy/lbm_configs/config_cuda.json similarity index 100% rename from src/CelerisLab/lbm/configs/config_cuda.json rename to legacy/lbm_configs/config_cuda.json diff --git a/src/CelerisLab/lbm/configs/config_flowfield.json b/legacy/lbm_configs/config_flowfield.json similarity index 100% rename from src/CelerisLab/lbm/configs/config_flowfield.json rename to legacy/lbm_configs/config_flowfield.json diff --git a/src/CelerisLab/lbm/driver.py b/legacy/lbm_driver.py similarity index 100% rename from src/CelerisLab/lbm/driver.py rename to legacy/lbm_driver.py diff --git a/src/CelerisLab/lbm/kernels/macros.h b/legacy/macros.h similarity index 90% rename from src/CelerisLab/lbm/kernels/macros.h rename to legacy/macros.h index 0ed1b25..7962b26 100644 --- a/src/CelerisLab/lbm/kernels/macros.h +++ b/legacy/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 \ No newline at end of file diff --git a/src/CelerisLab/__init__.py b/src/CelerisLab/__init__.py index e64195c..c3587de 100644 --- 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__"] diff --git a/src/CelerisLab/body/__init__.py b/src/CelerisLab/body/__init__.py new file mode 100644 index 0000000..1a2ba7d --- /dev/null +++ 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"] diff --git a/src/CelerisLab/body/manager.py b/src/CelerisLab/body/manager.py new file mode 100644 index 0000000..e45891b --- /dev/null +++ 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 diff --git a/src/CelerisLab/body/objects.py b/src/CelerisLab/body/objects.py new file mode 100644 index 0000000..5b23691 --- /dev/null +++ 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 diff --git a/src/CelerisLab/common/__init__.py b/src/CelerisLab/common/__init__.py index 7663ec6..0aea2aa 100644 --- 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'] diff --git a/src/CelerisLab/config.py b/src/CelerisLab/config.py new file mode 100644 index 0000000..ccb264e --- /dev/null +++ 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", [])) diff --git a/src/CelerisLab/configs/config_body.json b/src/CelerisLab/configs/config_body.json new file mode 100644 index 0000000..46d51d6 --- /dev/null +++ b/src/CelerisLab/configs/config_body.json @@ -0,0 +1,3 @@ +{ + "objects": [] +} diff --git a/src/CelerisLab/configs/config_lbm.json b/src/CelerisLab/configs/config_lbm.json new file mode 100644 index 0000000..842750e --- /dev/null +++ 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" + } +} diff --git a/src/CelerisLab/cuda/__init__.py b/src/CelerisLab/cuda/__init__.py index 8151d57..7e439dc 100644 --- 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'] diff --git a/src/CelerisLab/cuda/compiler_v2.py b/src/CelerisLab/cuda/compiler_v2.py new file mode 100644 index 0000000..f6a438a --- /dev/null +++ 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) diff --git a/src/CelerisLab/cuda/context.py b/src/CelerisLab/cuda/context.py new file mode 100644 index 0000000..199f3c3 --- /dev/null +++ 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]}" diff --git a/src/CelerisLab/lbm/__init__.py b/src/CelerisLab/lbm/__init__.py index 0a8f33b..5e94790 100644 --- 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 + diff --git a/src/CelerisLab/lbm/field.py b/src/CelerisLab/lbm/field.py new file mode 100644 index 0000000..e78e89a --- /dev/null +++ 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() diff --git a/src/CelerisLab/lbm/initializers.py b/src/CelerisLab/lbm/initializers.py new file mode 100644 index 0000000..3a9555f --- /dev/null +++ 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 diff --git a/src/CelerisLab/lbm/kernels/D2Q9.cu b/src/CelerisLab/lbm/kernels/D2Q9.cu deleted file mode 100644 index 395d8ed..0000000 --- a/src/CelerisLab/lbm/kernels/D2Q9.cu +++ /dev/null @@ -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; -} \ No newline at end of file diff --git a/src/CelerisLab/lbm/kernels/IO.cu b/src/CelerisLab/lbm/kernels/IO.cu deleted file mode 100644 index e69de29..0000000 diff --git a/src/CelerisLab/lbm/kernels/boundary/inlet_outlet.cuh b/src/CelerisLab/lbm/kernels/boundary/inlet_outlet.cuh index 3b1e284..6e47956 100644 --- 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 \ No newline at end of file diff --git a/src/CelerisLab/lbm/kernels/config.h b/src/CelerisLab/lbm/kernels/config.h new file mode 100644 index 0000000..52f27f6 --- /dev/null +++ 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 diff --git a/src/CelerisLab/lbm/kernels/config/config_grid.h b/src/CelerisLab/lbm/kernels/config/config_grid.h new file mode 100644 index 0000000..6d5aee7 --- /dev/null +++ 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 diff --git a/src/CelerisLab/lbm/kernels/config/config_method.h b/src/CelerisLab/lbm/kernels/config/config_method.h new file mode 100644 index 0000000..cc1bc1f --- /dev/null +++ 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 diff --git a/src/CelerisLab/lbm/kernels/config/config_objects.h b/src/CelerisLab/lbm/kernels/config/config_objects.h new file mode 100644 index 0000000..a416917 --- /dev/null +++ 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 diff --git a/src/CelerisLab/lbm/kernels/config/config_physics.h b/src/CelerisLab/lbm/kernels/config/config_physics.h new file mode 100644 index 0000000..99d3b46 --- /dev/null +++ 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 diff --git a/src/CelerisLab/lbm/kernels/const.h b/src/CelerisLab/lbm/kernels/const.h deleted file mode 100644 index e21bcc3..0000000 --- a/src/CelerisLab/lbm/kernels/const.h +++ /dev/null @@ -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 \ No newline at end of file diff --git a/src/CelerisLab/lbm/kernels/core/flags.cuh b/src/CelerisLab/lbm/kernels/core/flags.cuh index 7524651..8c135df 100644 --- 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 diff --git a/src/CelerisLab/lbm/kernels/kernel.cu b/src/CelerisLab/lbm/kernels/kernel.cu deleted file mode 100644 index 4bf36f8..0000000 --- a/src/CelerisLab/lbm/kernels/kernel.cu +++ /dev/null @@ -1,222 +0,0 @@ -// CelerisLab/kernels/kernel.cu - -#include -#include -#include - -#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(&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; - // } - - - // } -} \ No newline at end of file diff --git a/src/CelerisLab/lbm/kernels/kernel_v2.cu b/src/CelerisLab/lbm/kernels/kernel_v2.cu index 717e4b2..36c353f 100644 --- 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" diff --git a/src/CelerisLab/lbm/kernels/preproc.cu b/src/CelerisLab/lbm/kernels/preproc.cu deleted file mode 100644 index c2b25c5..0000000 --- a/src/CelerisLab/lbm/kernels/preproc.cu +++ /dev/null @@ -1,2 +0,0 @@ -#include "macros.h" -#include "const.h" diff --git a/src/CelerisLab/lbm/kernels/step/one_step_esopull.cu b/src/CelerisLab/lbm/kernels/step/one_step_esopull.cu index d93b094..1f44f53 100644 --- 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 ----- diff --git a/src/CelerisLab/lbm/stepper.py b/src/CelerisLab/lbm/stepper.py new file mode 100644 index 0000000..59843b5 --- /dev/null +++ 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 diff --git a/src/CelerisLab/simulation.py b/src/CelerisLab/simulation.py new file mode 100644 index 0000000..24676e2 --- /dev/null +++ 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 diff --git a/tests/test_high_re_validation.py b/tests/test_high_re_validation.py index 8920894..835a8de 100644 --- 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__": diff --git a/tests/test_stability_matrix.py b/tests/test_stability_matrix.py new file mode 100644 index 0000000..463b00d --- /dev/null +++ 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()