Frank14f/CelerisLab
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Initial commit: CelerisLab v0.2.0 with src layout

Browse Source
This commit is contained in:
Frank14f 2026-02-15 22:36:46 +08:00
commit 99c175042a
22 changed files with 1733 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

68
.gitignore vendored Normal file
View File

@ -0,0 +1,68 @@
# Byte-compiled / optimized / DLL files
__pycache__/
*.py[cod]
*$py.class
*.so
# Distribution / packaging
.Python
build/
develop-eggs/
dist/
downloads/
eggs/
.eggs/
lib/
lib64/
parts/
sdist/
var/
wheels/
*.egg-info/
.installed.cfg
*.egg
# CUDA compilation outputs
*.ptx
*.cubin
# IDE
.vscode/
.idea/
*.swp
*.swo
*~
# Jupyter Notebook
.ipynb_checkpoints
# PyCharm
.idea/
# Unit test / coverage reports
htmlcov/
.tox/
.coverage
.coverage.*
.cache
.pytest_cache/
nosetests.xml
coverage.xml
*.cover
# Environments
.env
.venv
env/
venv/
ENV/
env.bak/
venv.bak/
# MacOS
.DS_Store
# Temporary files
*.tmp
*.bak
*.log

21
LICENSE Normal file
View File

@ -0,0 +1,21 @@
MIT License
Copyright (c) 2026 Frank14f
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.

224
README.md Normal file
View File

@ -0,0 +1,224 @@
# CelerisLab
**GPU-Accelerated Lattice Boltzmann Method (LBM) CFD Solver**
CelerisLab is a high-performance computational fluid dynamics (CFD) solver based on the Lattice Boltzmann Method, leveraging NVIDIA CUDA for GPU acceleration. It provides a Python interface for easy integration into scientific workflows while maintaining high computational efficiency through CUDA kernels.
## 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
- **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
## Installation
### Prerequisites
- Python 3.8 or higher
- NVIDIA GPU with CUDA Compute Capability 6.0 or higher
- CUDA Toolkit 11.0 or higher
- NVIDIA drivers
### Install from source
```bash
git clone <repository_url>
cd CelerisLab
pip install -e . # Installs from src/ directory
```
### Dependencies
- `pycuda>=2020.1`: CUDA Python bindings
- `numpy>=1.19.0`: Numerical computing
- `scipy>=1.5.0`: Scientific computing (special functions for vortex initialization)
## Quick Start
### Basic Flow Simulation
```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]}")
```
### Configuration
CelerisLab searches for configuration files in the following order:
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
}
```
**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"]
}
}
```
## 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
```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
)
```
### Multiple Sensors
```python
# Add sensors in a line downstream of obstacle
for i in range(5):
flow.add_sensor(position=(100 + i*10, 50, 0))
# Read all sensors at once
sensor_data = flow.read_sensor() # Returns array of shape (n_sensors, 3)
```
### Vortex Initialization
```python
# Initialize Lamb-Oseen vortex
flow.add_vortex(
center=(50, 50, 0),
circulation=1.0,
core_radius=10.0,
vortex_type='Lamb'
)
```
## Environment Variables
- `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)
## 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
## Citation
If you use CelerisLab in your research, please cite:
```bibtex
@software{celerislab2026,
author = {Frank14f},
title = {CelerisLab: GPU-Accelerated Lattice Boltzmann Method Solver},
year = {2026},
url = {https://github.com/frank14f/CelerisLab}
}
```
## License
MIT License - see LICENSE file for details
## 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

9
configs/config_cuda.json Normal file
View File

@ -0,0 +1,9 @@
{
"multi_gpu": false,
"gpu_connection": "NVLink",
"required_cuda_capability": "7.0",
"threads_per_block": 128,
"X_1U": 128,
"Y_1U": 32,
"Z_1U": 1
}

13
configs/config_flowfield.json Normal file
View File

@ -0,0 +1,13 @@
{
"data_type": "FP32",
"dimensionality": 2,
"lattice": 9,
"field_dim_in_U": [10, 16, 1],
"viscosity": 0.002,
"velocity": 0.01,
"boundary_conditions": {
"x": ["parabolic", "outflow"],
"y": ["noslip", "noslip"],
"z": ["none", "none"]
}
}

3
configs/config_gym.json Normal file
View File

@ -0,0 +1,3 @@
{
}

48
pyproject.toml Normal file
View File

@ -0,0 +1,48 @@
[build-system]
requires = ["setuptools>=61.0", "wheel"]
build-backend = "setuptools.build_meta"
[project]
name = "CelerisLab"
version = "0.2.0"
description = "GPU-accelerated Lattice Boltzmann Method (LBM) CFD solver using CUDA"
readme = "README.md"
requires-python = ">=3.8"
license = {text = "MIT"}
authors = [
{name = "Frank14f"}
]
keywords = ["cfd", "lattice-boltzmann", "cuda", "gpu", "fluid-dynamics", "lbm"]
classifiers = [
"Development Status :: 3 - Alpha",
"Intended Audience :: Science/Research",
"Topic :: Scientific/Engineering :: Physics",
"License :: OSI Approved :: MIT License",
"Programming Language :: Python :: 3",
"Programming Language :: Python :: 3.8",
"Programming Language :: Python :: 3.9",
"Programming Language :: Python :: 3.10",
"Programming Language :: Python :: 3.11",
]
dependencies = [
"pycuda>=2020.1",
"numpy>=1.19.0",
"scipy>=1.5.0",
]
[project.urls]
Homepage = "https://github.com/frank14f/CelerisLab"
Repository = "https://github.com/frank14f/CelerisLab.git"
[project.scripts]
celerislab = "CelerisLab.driver:main"
[tool.setuptools]
package-dir = {"" = "src"}
[tool.setuptools.packages.find]
where = ["src"]
[tool.setuptools.package-data]
CelerisLab = ["kernels/*.cu", "kernels/*.h", "configs/*.json"]

45
setup.py Normal file
View File

@ -0,0 +1,45 @@
from setuptools import setup, find_packages
with open("README.md", "r", encoding="utf-8") as fh:
long_description = fh.read()
setup(
name='CelerisLab',
version='0.2.0',
author='Frank14f',
description='GPU-accelerated Lattice Boltzmann Method (LBM) CFD solver using CUDA',
long_description=long_description,
long_description_content_type="text/markdown",
url='https://github.com/frank14f/CelerisLab',
packages=find_packages(where='src'),
package_dir={'': 'src'},
package_data={
'CelerisLab': [
'kernels/*.cu',
'kernels/*.h',
'configs/*.json',
],
},
install_requires=[
'pycuda>=2020.1',
'numpy>=1.19.0',
'scipy>=1.5.0',
],
python_requires='>=3.8',
classifiers=[
'Development Status :: 3 - Alpha',
'Intended Audience :: Science/Research',
'Topic :: Scientific/Engineering :: Physics',
'License :: OSI Approved :: MIT License',
'Programming Language :: Python :: 3',
'Programming Language :: Python :: 3.8',
'Programming Language :: Python :: 3.9',
'Programming Language :: Python :: 3.10',
'Programming Language :: Python :: 3.11',
],
entry_points={
'console_scripts': [
'celerislab=CelerisLab.driver:main',
],
},
)

24
src/CelerisLab/__init__.py Normal file
View File

@ -0,0 +1,24 @@
# CelerisLab/__init__.py
"""
CelerisLab: GPU-Accelerated Lattice Boltzmann Method CFD Solver
"""
__version__ = '0.2.0'
# Always import utils (no pycuda dependency)
from . import utils
# Try to import FlowField (requires pycuda)
try:
from .driver import FlowField
__all__ = ['FlowField', 'utils']
except ImportError as e:
# PyCUDA not available, only utils module will be accessible
import warnings
warnings.warn(
f"FlowField not available: {e}. "
"Install pycuda to use the full CelerisLab functionality. "
"Utils module is still accessible for configuration management.",
ImportWarning
)
__all__ = ['utils']

87
src/CelerisLab/compiler.py Normal file
View File

@ -0,0 +1,87 @@
# CelerisLab/kernels/compiler.py
import subprocess
import re
import os
from .utils import FlowFieldConfig, CudaConfig
def kernel_path(file_name: str) -> str:
current_dir = os.path.dirname(os.path.abspath(__file__))
return os.path.join(current_dir, "kernels", file_name)
def read_lines(file_path):
with open(file_path, "r") as file:
lines = file.readlines()
return lines
def write_lines(file_path, lines):
with open(file_path, "w") as file:
file.writelines(lines)
def modify_macro(lines, macro_name, new_value):
pattern = re.compile(rf"(#define\s+{macro_name}\s+)(\S+)")
for i, line in enumerate(lines):
if pattern.match(line):
lines[i] = pattern.sub(rf"\g<1>{new_value}", line)
break
return lines
def modify_const(lines, const_name, new_type, new_value):
pattern = re.compile(rf"(__constant__\s+)(\S+\s+{const_name}\s*=\s*)([^;]+)(;)")
for i, line in enumerate(lines):
if pattern.match(line):
lines[i] = pattern.sub(rf"\g<1>{new_type} {const_name} = {new_value}\4", line)
break
return lines
def compile_kernel():
subprocess.run(
[
"nvcc",
"-ptx",
kernel_path("kernel.cu"),
"-o",
kernel_path("kernel.ptx"),
]
)
def config_kernal(config_cuda: CudaConfig, config_field: FlowFieldConfig):
lines = read_lines(kernel_path("macros.h"))
lines = modify_macro(lines, "MULT_GPU", config_cuda.multi_gpu)
lines = modify_macro(lines, "NT", config_cuda.threads_per_block)
lines = modify_macro(lines, "X_1U", config_cuda.unit_dimensions[0])
lines = modify_macro(lines, "Y_1U", config_cuda.unit_dimensions[1])
lines = modify_macro(lines, "Z_1U", config_cuda.unit_dimensions[2])
if config_field.data_type == "FP32":
lb_type = "float"
else:
raise ValueError(f"Unsupported data type {config_field.data_type}.")
lines = modify_macro(lines, "LBtype", lb_type)
lines = modify_macro(lines, "UX", config_field.field_dim_in_U[0])
lines = modify_macro(lines, "UY", config_field.field_dim_in_U[1])
lines = modify_macro(lines, "UZ", config_field.field_dim_in_U[2])
lines = modify_macro(lines, "NX", config_field.field_dim_in_U[0] * config_cuda.unit_dimensions[0])
lines = modify_macro(lines, "NY", config_field.field_dim_in_U[1] * config_cuda.unit_dimensions[1])
lines = modify_macro(lines, "NZ", config_field.field_dim_in_U[2] * config_cuda.unit_dimensions[2])
lines = modify_macro(lines, "DIM", config_field.dimensionality)
lines = modify_macro(lines, "NQ", config_field.lattice)
lines = modify_macro(lines, "VIS", config_field.viscosity)
lines = modify_macro(lines, "U0", config_field.velocity)
write_lines(kernel_path("macros.h"), lines)
def config_object(n_obj: int):
lines = read_lines(kernel_path("macros.h"))
lines = modify_macro(lines, "N_OBJS", n_obj)
write_lines(kernel_path("macros.h"), lines)
def config_sensor(n_sen: int):
lines = read_lines(kernel_path("macros.h"))
lines = modify_macro(lines, "N_SENS", n_sen)
write_lines(kernel_path("macros.h"), lines)

9
src/CelerisLab/configs/config_cuda.json Normal file
View File

@ -0,0 +1,9 @@
{
"multi_gpu": false,
"gpu_connection": "NVLink",
"required_cuda_capability": "7.0",
"threads_per_block": 128,
"X_1U": 128,
"Y_1U": 32,
"Z_1U": 1
}

13
src/CelerisLab/configs/config_flowfield.json Normal file
View File

@ -0,0 +1,13 @@
{
"data_type": "FP32",
"dimensionality": 2,
"lattice": 9,
"field_dim_in_U": [10, 16, 1],
"viscosity": 0.002,
"velocity": 0.01,
"boundary_conditions": {
"x": ["parabolic", "outflow"],
"y": ["noslip", "noslip"],
"z": ["none", "none"]
}
}

3
src/CelerisLab/configs/config_gym.json Normal file
View File

@ -0,0 +1,3 @@
{
}

391
src/CelerisLab/driver.py Normal file
View File

@ -0,0 +1,391 @@
# CelerisLab/driver.py
import pycuda.driver as cuda
import numpy as np
import struct
from scipy.special import jv, expi
from typing import List, Tuple, Union, Optional
from . import utils
from . import preprocess as preproc
from . import compiler
FLUID = 0b00000001
SOLID = 0b00000010
GAS = 0b00000100
INTERFACE = 0b00001000
SENSOR = 0b00010000
V_TAYLOR = np.int32(1)
class FlowField:
def __init__(
self,
field_config: utils.FlowFieldConfig,
cuda_config: utils.CudaConfig,
device_id: Union[int, List[int]] = None,
):
self.field_config = field_config
self.cuda_config = cuda_config
cuda.init()
# Sanity checks
if cuda_config.multi_gpu:
if device_id is None or isinstance(device_id, int):
raise ValueError("Multi-GPU support requires a list of device IDs.")
# self.devices = [cuda.Device(id) for id in device_id]
raise NotImplementedError("Multi-GPU support is not implemented yet.")
else:
if isinstance(device_id, list):
if len(device_id) > 1:
raise ValueError(
"Single-GPU mode does not support multiple device IDs."
)
device_id = device_id[0]
elif device_id is None:
device_id = 0
utils.check_cuda_device_availability(device_id)
self.device = cuda.Device(device_id)
self.context = self.device.make_context()
utils.check_cuda_capability(field_config, cuda_config, device_id)
# Config kernel
compiler.config_kernal(cuda_config, field_config)
compiler.config_object(int(0))
# compiler.config_sensor(int(0))
# Set constants
if field_config.data_type == "FP32":
self.DATA_TYPE = np.float32
else:
raise ValueError(f"Unsupported data type {field_config.data_type}.")
self.FIELD_SHAPE = tuple(
size * unit
for size, unit in zip(
field_config.field_dim_in_U, cuda_config.unit_dimensions
)
)
self.FIELD_SIZE = np.prod(self.FIELD_SHAPE)
self.LATTICE = field_config.lattice
self.DIM = field_config.dimensionality
if field_config.lattice == 9 and field_config.dimensionality == 2:
self.E = 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)
self.OPP = np.array([0, 3, 4, 1, 2, 7, 8, 5, 6], dtype=np.int32)
self.WW = np.array(
[4 / 9, 1 / 9, 1 / 9, 1 / 9, 1 / 9, 1 / 36, 1 / 36, 1 / 36, 1 / 36],
dtype=self.DATA_TYPE,
)
else:
raise NotImplementedError(
f"Unsupported lattice type {field_config.lattice} in {field_config.dimensionality} dimensions."
)
# Compile kernel
compiler.compile_kernel()
self.ptx = cuda.module_from_file(compiler.kernel_path("kernel.ptx"))
self.step = self.ptx.get_function("OneStep")
initflow = self.ptx.get_function("InitTubeFlow")
# Initialize memory
self.ddf = np.zeros(self.FIELD_SIZE * self.LATTICE, dtype=self.DATA_TYPE)
self.ddf_save = np.zeros(self.FIELD_SIZE * self.LATTICE, dtype=self.DATA_TYPE)
self.flag = np.ones(self.FIELD_SIZE, dtype=np.uint8)
self.indx = np.zeros(self.FIELD_SIZE, dtype=np.int32)
self.delta_curve = np.zeros(0, dtype=self.DATA_TYPE)
self.vortex_config = np.zeros(7, dtype=float)
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(1)
self.vortex_gpu = cuda.mem_alloc(self.vortex_config.nbytes)
self.objects = {}
self.action = np.zeros(0, dtype=self.DATA_TYPE)
self.obs = np.zeros(0, dtype=self.DATA_TYPE)
initflow(
self.flag_gpu,
self.ddf_gpu,
block=(self.cuda_config.threads_per_block, 1, 1),
grid=(
int(self.FIELD_SHAPE[0] / self.cuda_config.threads_per_block),
int(self.FIELD_SHAPE[1]),
int(self.FIELD_SHAPE[2]),
),
)
cuda.memcpy_dtoh(self.flag, self.flag_gpu)
cuda.memcpy_dtoh(self.ddf, self.ddf_gpu)
def add_cylinder(self, center: Tuple[float, float, float], radius: float, id_obj: Optional[int] = None):
x_c, y_c, z_c = center
if (
x_c - radius <= 0
or x_c + radius >= self.FIELD_SHAPE[0] - 1
or y_c - radius <= 0
or y_c + radius >= self.FIELD_SHAPE[1] - 1
):
raise ValueError("Cylinder is out of bounds.")
index = self.delta_curve.size if self.delta_curve.size > 0 else 0
if self.DATA_TYPE == np.float32:
id_object = np.int32(len(self.objects))
# max_id = max(self.objects.keys())
else:
raise ValueError(f"Unsupported data type {self.DATA_TYPE}.")
for x in range(int(x_c - radius) - 1, int(x_c + radius) + 1):
for y in range(int(y_c - radius) - 1, int(y_c + radius) + 1):
if (x - x_c) ** 2 + (y - y_c) ** 2 < radius**2:
k = x + y * self.FIELD_SHAPE[0]
self.flag[k] = SOLID
delta_temp = np.zeros(11, dtype=self.DATA_TYPE)
delta_temp[0] = id_object.view(self.DATA_TYPE)
for i in range(self.LATTICE):
x_neb = x + self.E[i][0]
y_neb = y + self.E[i][1]
if (x_neb - x_c) ** 2 + (y_neb - y_c) ** 2 >= radius**2:
self.flag[k] |= INTERFACE
x_i, y_i = preproc.find_circle_intersection(
x, y, x_neb, y_neb, x_c, y_c, radius
)
d_neb = np.sqrt((x_i - x_neb) ** 2 + (y_i - y_neb) ** 2)
delta_temp[i] = d_neb / np.sqrt(
self.E[i][0] ** 2 + self.E[i][1] ** 2
)
if self.flag[k] & INTERFACE:
delta_temp[9] = (y_c - y) / radius
delta_temp[10] = (x - x_c) / radius
self.delta_curve = np.concatenate(
(self.delta_curve, delta_temp)
)
self.indx[k] = index
index += delta_temp.size
self.objects[id_object] = {
"type": "cylinder",
"center": center,
"radius": radius,
}
if hasattr(self, "delta_gpu"):
self.delta_gpu.free()
self.delta_gpu = cuda.mem_alloc(self.delta_curve.nbytes)
self.action = np.zeros(len(self.objects), dtype=self.DATA_TYPE)
if hasattr(self, "action_gpu"):
self.action_gpu.free()
self.action_gpu = cuda.mem_alloc(self.action.nbytes)
self.obs = np.zeros(len(self.objects) * self.DIM, dtype=self.DATA_TYPE)
if hasattr(self, "obs_gpu"):
self.obs_gpu.free()
self.obs_gpu = cuda.mem_alloc(self.obs.nbytes)
cuda.memcpy_htod(self.delta_gpu, self.delta_curve)
cuda.memcpy_htod(self.flag_gpu, self.flag)
cuda.memcpy_htod(self.indx_gpu, self.indx)
compiler.config_object(len(self.objects))
compiler.compile_kernel()
self.ptx = cuda.module_from_file(compiler.kernel_path("kernel.ptx"))
self.step = self.ptx.get_function("OneStep")
def add_sensor(self, center: Tuple[float, float, float], radius: float):
x_c, y_c, z_c = center
if (
x_c - radius <= 0
or x_c + radius >= self.FIELD_SHAPE[0] - 1
or y_c - radius <= 0
or y_c + radius >= self.FIELD_SHAPE[1] - 1
):
raise ValueError("Sensor is out of bounds.")
id_object = len(self.objects)
for x in range(int(x_c - radius) - 1, int(x_c + radius) + 1):
for y in range(int(y_c - radius) - 1, int(y_c + radius) + 1):
if (x - x_c) ** 2 + (y - y_c) ** 2 < radius**2:
k = x + y * self.FIELD_SHAPE[0]
self.flag[k] |= SENSOR
self.indx[k] = id_object
self.objects[id_object] = {
"type": "sensor",
"center": center,
}
self.action = np.zeros(len(self.objects), dtype=self.DATA_TYPE)
if hasattr(self, "action_gpu"):
self.action_gpu.free()
self.action_gpu = cuda.mem_alloc(self.action.nbytes)
self.obs = np.zeros(len(self.objects) * self.DIM, dtype=self.DATA_TYPE)
if hasattr(self, "force_gpu"):
self.obs_gpu.free()
self.obs_gpu = cuda.mem_alloc(self.obs.nbytes)
cuda.memcpy_htod(self.flag_gpu, self.flag)
cuda.memcpy_htod(self.indx_gpu, self.indx)
compiler.config_object(len(self.objects))
compiler.compile_kernel()
self.ptx = cuda.module_from_file(compiler.kernel_path("kernel.ptx"))
self.step = self.ptx.get_function("OneStep")
def add_vortex(self, center: Tuple[float, float, float], radius: float, strength: float, direction: float, type: str):
x_c, y_c, z_c = center
if (
x_c - radius <= 0
or x_c + radius >= self.FIELD_SHAPE[0] - 1
or y_c - radius <= 0
or y_c + radius >= self.FIELD_SHAPE[1] - 1
):
raise ValueError("Vortex is out of bounds.")
if type not in ["lamb", "oseen", "taylor"]:
raise ValueError("Vortex type" + type + " is not supported.")
x = np.linspace(-x_c, self.FIELD_SHAPE[0] - 1 - x_c, self.FIELD_SHAPE[0])
y = np.linspace(-y_c, self.FIELD_SHAPE[1] - 1 - y_c, self.FIELD_SHAPE[1])
X, Y = np.meshgrid(x, y)
r = np.sqrt(X**2 + Y**2)
nu = self.field_config.viscosity
theta = np.arctan2(Y, X)
psi = np.zeros_like(r)
if type == "lamb":
b = 3.831705970207512
n = b / radius
u0 = strength
inside = r <= radius
outside = r > radius
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
elif type == "oseen":
# 4 nu t = radius^2 / 4
kappa = 2 * np.pi * radius **2 * strength
u_vor = - kappa / (2 * np.pi * r) * (1 - np.exp(-4 * r**2 / radius**2)) * np.sin(theta)
v_vor = kappa / (2 * np.pi * r) * (1 - np.exp(-4 * r**2 / radius**2)) * np.cos(theta)
zeta = 4 * r**2 / radius**2
p_vor = -kappa**2 / 8 / np.pi**2 / r**2 * (-2 * zeta * (expi(-zeta) - expi(-2 * zeta)) + (1 - np.exp(-zeta))**2)
elif type == "taylor":
# 4 nu t = radius^2
M = strength * np.pi * radius**4 / 8 / nu
u_vor = - M * r * 4 * nu / radius**4 * np.exp(-r**2 / radius**2) * np.sin(theta)
v_vor = M * r * 4 * nu / radius**4 * np.exp(-r**2 / radius**2) * np.cos(theta)
p_vor = -4 * M**2 * nu**2 * np.exp(-2 * r**2 / radius**2) / np.pi**2 / radius**6
cuda.memcpy_dtoh(self.ddf, self.ddf_gpu)
ddf_temp = self.ddf.copy().reshape((self.LATTICE, self.FIELD_SHAPE[1], self.FIELD_SHAPE[0])).transpose(2, 1, 0)
u_ddf = ddf_temp[:, :, 1] + ddf_temp[:, :, 5] + ddf_temp[:, :, 8] - ddf_temp[:, :, 3] - ddf_temp[:, :, 6] - ddf_temp[:, :, 7]
v_ddf = ddf_temp[:, :, 2] + ddf_temp[:, :, 5] + ddf_temp[:, :, 6] - ddf_temp[:, :, 4] - ddf_temp[:, :, 7] - ddf_temp[:, :, 8]
p_ddf = np.sum(ddf_temp, axis=2) / 3
for i in range(self.FIELD_SHAPE[0]):
for j in range(self.FIELD_SHAPE[1]):
k = i + j * self.FIELD_SHAPE[0]
if (j == 0 or j == self.FIELD_SHAPE[1] - 1) or (i == 0 or i == self.FIELD_SHAPE[0] - 1):
continue
else:
for e in range(self.LATTICE):
u = u_ddf[i, j] + u_vor[j, i]
v = v_ddf[i, j] + v_vor[j, i]
p = p_ddf[i, j] + p_vor[j, i]
eu = self.E[e][0] * u + self.E[e][1] * v
u2 = u ** 2 + v ** 2
self.ddf[k + e * self.FIELD_SIZE] = self.WW[e] * (3 * p + 3 * eu + 4.5 * eu ** 2 - 1.5 * u2)
cuda.memcpy_htod(self.ddf_gpu, self.ddf)
# def add_vortex_gpu(self, center: Tuple[float, float, float], radius: float, strength: float, direction: float, type: str):
# x_c, y_c, z_c = center
# if (
# x_c - radius <= 0
# or x_c + radius >= self.FIELD_SHAPE[0] - 1
# or y_c - radius <= 0
# or y_c + radius >= self.FIELD_SHAPE[1] - 1
# ):
# raise ValueError("Vortex is out of bounds.")
# if type not in ["lamb", "oseen", "taylor"]:
# raise ValueError("Vortex type" + type + " is not supported.")
# add_vortex = self.ptx.get_function("AddVortex")
# self.vortex_config[0:3] = np.array(center, dtype=float)
# self.vortex_config[3] = radius
# self.vortex_config[4] = strength
# self.vortex_config[5] = direction
# if type == "taylor":
# self.vortex_config[6] =
def run(self, num_steps: int, action_target: np.ndarray):
if (
action_target.size != len(self.objects)
or action_target.dtype != self.DATA_TYPE
):
raise ValueError("action data type or size does not match the objects.")
elif len(self.objects) == 0:
raise ValueError("No objects have been added to the flow field.")
weight = 0.1
stream = cuda.Stream()
action_pinned = cuda.pagelocked_empty_like(self.action)
action_pinned[:] = self.action
obs_pinned = cuda.pagelocked_empty_like(self.obs)
self.obs[:] = 0
for i in range(num_steps):
action_pinned = (1 - weight) * action_pinned + weight * action_target
cuda.memcpy_htod_async(self.action_gpu, action_pinned, stream)
self.step(
self.flag_gpu,
self.ddf_gpu,
self.temp_gpu,
self.indx_gpu,
self.delta_gpu,
self.action_gpu,
self.obs_gpu,
block=(self.cuda_config.threads_per_block, 1, 1),
grid=(
int(self.FIELD_SHAPE[0] / self.cuda_config.threads_per_block),
int(self.FIELD_SHAPE[1]),
int(self.FIELD_SHAPE[2]),
),
stream=stream,
)
self.ddf_gpu, self.temp_gpu = self.temp_gpu, self.ddf_gpu
cuda.memcpy_dtoh_async(obs_pinned, self.obs_gpu, stream)
cuda.memset_d32_async(self.obs_gpu, 0, self.obs.size, stream)
self.obs += obs_pinned
stream.synchronize()
self.obs = (self.obs / num_steps).astype(self.DATA_TYPE)
def apply_ddf(self):
cuda.memcpy_htod(self.ddf_gpu, self.ddf)
def get_ddf(self):
cuda.memcpy_dtoh(self.ddf, self.ddf_gpu)
def save_ddf(self):
self.ddf_save = self.ddf.copy()
def restore_ddf(self):
self.ddf = self.ddf_save.copy()
def __del__(self):
self.context.pop()

101
src/CelerisLab/kernels/D2Q9.cu Normal file
View File

@ -0,0 +1,101 @@
#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;
}

0
src/CelerisLab/kernels/IO.cu Normal file
View File

10
src/CelerisLab/kernels/const.h Normal file
View File

@ -0,0 +1,10 @@
// 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

222
src/CelerisLab/kernels/kernel.cu Normal file
View File

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

37
src/CelerisLab/kernels/macros.h Normal file
View File

@ -0,0 +1,37 @@
// CelerisLab/kernels/macros.h
// cuda parameters
#define MULT_GPU False
#define NT 128
#define X_1U 128
#define Y_1U 32
#define Z_1U 1
// flow parameters
#define LBtype float
#define UX 10
#define UY 16
#define UZ 1
#define NX 1280
#define NY 512
#define NZ 1
#define DIM 2
#define NQ 9
#define VIS 0.004
#define RHO 1.0
#define U0 0.01
// constants
#define PI 3.141592653589793238
#define FLUID 0b00000001
#define SOLID 0b00000010
#define GAS 0b00000100
#define INTERFACE 0b00001000
#define SENSOR 0b00010000
// vortex type
#define V_TAYLOR 0b00000001
// variables
#define N_OBJS 7
// #define N_SENS 2

2
src/CelerisLab/kernels/preproc.cu Normal file
View File

@ -0,0 +1,2 @@
#include "macros.h"
#include "const.h"

40
src/CelerisLab/preprocess.py Normal file
View File

@ -0,0 +1,40 @@
# CelerisLab/preprocess.py
import math
import numpy as np
from typing import Tuple
FLUID = 0b00000001
SOLID = 0b00000010
GAS = 0b00000100
INTERFACE = 0b00001000
SENSOR = 0b00010000
def find_circle_intersection(x, y, x_neb, y_neb, xc, yc, r0):
dx, dy = x_neb - x, y_neb - y
a = dx ** 2 + dy ** 2
b = 2 * (dx * (x - xc) + dy * (y - yc))
c = (x - xc) ** 2 + (y - yc) ** 2 - r0 ** 2
det = b ** 2 - 4 * a * c
if det < 0:
return None
t1 = (-b + math.sqrt(det)) / (2 * a)
t2 = (-b - math.sqrt(det)) / (2 * a)
if 0 <= t1 <= 1:
return x + t1 * dx, y + t1 * dy
elif 0 <= t2 <= 1:
return x + t2 * dx, y + t2 * dy
else:
return None
def find_sensor_area(radius):
area = 0
for i in range(np.floor(-radius), np.ceil(radius)):
for j in range(np.floor(-radius), np.ceil(radius)):
if i ** 2 + j ** 2 <= radius ** 2:
area += 1
return area

363
src/CelerisLab/utils.py Normal file
View File

@ -0,0 +1,363 @@
# CelerisLab/utils.py
import pycuda.driver as cuda
import subprocess
import json
import os
from typing import NamedTuple, Optional, List, Tuple, Union
class CudaDeviceInfo(NamedTuple):
name: str
compute_capability: str
multiprocessors: int
total_global_memory: int
max_shared_memory_per_block: int
max_threads_per_block: int
max_blocks_per_multiprocessor: int
device_interconnect: Optional[str] = None
class FlowFieldConfig(NamedTuple):
data_type: str
dimensionality: int
lattice: int
field_dim_in_U: Tuple[int, int, int]
viscosity: float
velocity: float
boundary_conditions: Tuple[str, str, str, str, str, str]
class CudaConfig(NamedTuple):
multi_gpu: bool
gpu_connection: str
required_cuda_capability: str
threads_per_block: int
unit_dimensions: Tuple[int, int, int]
def check_cuda_device_availability(device_id=0):
if cuda.Device.count() == 0:
raise RuntimeError("No CUDA device is available.")
if device_id < 0 or device_id >= cuda.Device.count():
raise ValueError(
f"Invalid device_id {device_id}. Must be between 0 and {cuda.Device.count() - 1}."
)
try:
subprocess.check_output(["nvidia-smi", "--version"])
except subprocess.CalledProcessError:
raise RuntimeError("nvidia-smi is not available or not installed correctly.")
def query_cuda_device_info(device_id=0) -> CudaDeviceInfo:
check_cuda_device_availability(device_id)
try:
output = subprocess.check_output(
["nvidia-smi", "-q", "-d", "TOPOLOGY", "-i", str(device_id)], text=True
)
if "NVLink" in output:
device_interconnect = "NVLink"
elif "PCIe" in output:
device_interconnect = "PCIe"
else:
device_interconnect = "Unknown"
except Exception as e:
device_interconnect = None
device = cuda.Device(device_id)
return CudaDeviceInfo(
name=device.name(),
compute_capability=f"{device.compute_capability()[0]}.{device.compute_capability()[1]}",
multiprocessors=device.get_attribute(
cuda.device_attribute.MULTIPROCESSOR_COUNT
),
total_global_memory=device.total_memory(),
max_shared_memory_per_block=device.get_attribute(
cuda.device_attribute.MAX_SHARED_MEMORY_PER_BLOCK
),
max_threads_per_block=device.get_attribute(
cuda.device_attribute.MAX_THREADS_PER_BLOCK
),
max_blocks_per_multiprocessor=device.get_attribute(
cuda.device_attribute.MAX_BLOCKS_PER_MULTIPROCESSOR
),
device_interconnect=device_interconnect,
)
def find_config_file(config_filename: str, config_path: Optional[str] = None) -> str:
"""
Find configuration file by searching in multiple locations.
Search priority:
1. Provided config_path (if given)
2. Environment variable CELERISLAB_CONFIG_DIR
3. Current working directory ./configs/
4. Package installation location (relative to this utils.py file)
Args:
config_filename: Name of the config file (e.g., 'config_cuda.json')
config_path: Optional explicit path to config file
Returns:
Absolute path to the config file
Raises:
FileNotFoundError: If config file cannot be found in any location
"""
search_paths = []
# Priority 1: Explicit path provided
if config_path:
search_paths.append(config_path)
# Priority 2: Environment variable
env_config_dir = os.environ.get('CELERISLAB_CONFIG_DIR')
if env_config_dir:
search_paths.append(os.path.join(env_config_dir, config_filename))
# Priority 3: Current working directory
search_paths.append(os.path.join(os.getcwd(), 'configs', config_filename))
# Priority 4: Package installation location (relative to this utils.py)
package_root = os.path.dirname(os.path.abspath(__file__))
search_paths.append(os.path.join(package_root, 'configs', config_filename))
# Search for the file
for path in search_paths:
if os.path.isfile(path):
return os.path.abspath(path)
# File not found, provide helpful error message
error_msg = f"Configuration file '{config_filename}' not found. Searched in:\n"
for path in search_paths:
error_msg += f" - {path}\n"
error_msg += "\nTo fix this, you can:\n"
error_msg += " 1. Set CELERISLAB_CONFIG_DIR environment variable\n"
error_msg += " 2. Place config files in ./configs/ directory\n"
error_msg += " 3. Provide explicit config_path parameter"
raise FileNotFoundError(error_msg)
def load_flow_field_config(config_path: Optional[str] = None) -> FlowFieldConfig:
"""
Load flow field configuration from JSON file.
Args:
config_path: Optional path to config file. If None, searches in standard locations.
Can be relative path like 'configs/config_flowfield.json' or just filename.
Returns:
FlowFieldConfig object
"""
# Determine config filename and full path
if config_path:
# Check if it's just a filename or a path
if os.path.basename(config_path) == config_path:
# Just a filename, search for it
config_file = find_config_file(config_path, None)
else:
# It's a path, use it if exists, otherwise try to find the basename
if os.path.isfile(config_path):
config_file = config_path
else:
config_file = find_config_file(os.path.basename(config_path), None)
else:
# No path provided, search for default filename
config_file = find_config_file('config_flowfield.json', None)
try:
with open(config_file, "r") as file:
config = json.load(file)
required_keys = [
"data_type",
"dimensionality",
"lattice",
"field_dim_in_U",
"viscosity",
"boundary_conditions",
]
if not all(key in config for key in required_keys):
raise ValueError("Missing required configuration items.")
if config["data_type"] not in ["FP32", "FP64"]:
raise ValueError("Data type must be either FP32 or FP64.")
if config["dimensionality"] not in [2, 3]:
raise ValueError("Dimensionality must be either 2 or 3.")
if config["dimensionality"] == 2 and config["field_dim_in_U"][2] != 1:
raise ValueError(
"Field dimensions must be 1 in the third dimension for 2D simulations."
)
if config["lattice"] not in [9]:
raise ValueError("Lattice must be either 9 or 19.")
boundary_conditions = tuple(
condition
for key in ["x", "y", "z"]
for condition in config["boundary_conditions"].get(key, [])
)
if len(boundary_conditions) != 6:
raise ValueError("Boundary conditions must contain exactly six elements.")
return FlowFieldConfig(
data_type=config["data_type"],
dimensionality=config["dimensionality"],
lattice=config["lattice"],
field_dim_in_U=tuple(config["field_dim_in_U"]),
viscosity=config["viscosity"],
velocity=config["velocity"],
boundary_conditions=boundary_conditions,
)
except Exception as e:
raise RuntimeError(f"Failed to load or parse the flow field configuration: {e}")
def load_cuda_config(config_path: Optional[str] = None) -> CudaConfig:
"""
Load CUDA configuration from JSON file.
Args:
config_path: Optional path to config file. If None, searches in standard locations.
Can be relative path like 'configs/config_cuda.json' or just filename.
Returns:
CudaConfig object
"""
# Determine config filename and full path
if config_path:
# Check if it's just a filename or a path
if os.path.basename(config_path) == config_path:
# Just a filename, search for it
config_file = find_config_file(config_path, None)
else:
# It's a path, use it if exists, otherwise try to find the basename
if os.path.isfile(config_path):
config_file = config_path
else:
config_file = find_config_file(os.path.basename(config_path), None)
else:
# No path provided, search for default filename
config_file = find_config_file('config_cuda.json', None)
try:
with open(config_file, "r") as file:
config = json.load(file)
required_keys = [
"multi_gpu",
"gpu_connection",
"required_cuda_capability",
"threads_per_block",
"X_1U",
"Y_1U",
"Z_1U",
]
if not all(key in config for key in required_keys):
raise ValueError("Missing required configuration items.")
return CudaConfig(
multi_gpu=config["multi_gpu"],
gpu_connection=config["gpu_connection"],
required_cuda_capability=config["required_cuda_capability"],
threads_per_block=config["threads_per_block"],
unit_dimensions=(config["X_1U"], config["Y_1U"], config["Z_1U"]),
)
except Exception as e:
raise RuntimeError(f"Failed to load or parse the CUDA configuration: {e}")
def check_cuda_capability(
field_config: FlowFieldConfig,
cuda_config: CudaConfig,
device_id: Union[int, List[int]] = None,
):
SAFE_FACTOR = 0.8
if cuda_config.multi_gpu:
if device_id is None or isinstance(device_id, int):
raise ValueError("Multi-GPU support requires a list of device IDs.")
raise NotImplementedError("Multi-GPU support is not implemented yet.")
else:
if isinstance(device_id, list):
if len(device_id) > 1:
raise ValueError(
"Single-GPU mode does not support multiple device IDs."
)
device_id = device_id[0]
elif device_id is None:
device_id = 0
device_info = query_cuda_device_info(device_id)
if device_info.compute_capability != cuda_config.required_cuda_capability:
raise ValueError(
f"Device {device_info.name} has compute capability {device_info.compute_capability}, but {cuda_config.required_cuda_capability} is required."
)
field_size = sum(
size * unit
for size, unit in zip(
field_config.field_dim_in_U, cuda_config.unit_dimensions
)
)
if (
device_info.total_global_memory * SAFE_FACTOR
< calc_field_memory_consumption(
field_size,
field_config.dimensionality,
field_config.lattice,
field_config.data_type,
)
):
raise ValueError(
f"Device {device_info.name} does not have enough memory to store the flow field."
)
if (
device_info.max_threads_per_block * SAFE_FACTOR
< cuda_config.threads_per_block
):
raise ValueError(
f"Device {device_info.name} does not have enough threads per block to run the simulation."
)
block_size = cuda_config.threads_per_block
if (
device_info.max_shared_memory_per_block * SAFE_FACTOR
< 2
* calc_field_memory_consumption(
block_size,
field_config.dimensionality,
field_config.lattice,
field_config.data_type,
)
):
raise ValueError(
f"Device {device_info.name} does not have enough shared memory per block to run the simulation."
)
def calc_field_memory_consumption(
field_size: int, dimensionality: int, directions: int, data_type: str
) -> int:
if data_type == "FP32":
data_size = 4
elif data_type == "FP64":
data_size = 8
else:
raise ValueError(f"Unsupported data type {data_type}.")
return (
field_size * directions * data_size * 2
+ field_size * dimensionality * data_size
+ field_size
)