commit 99c175042a0b9bf9c8af47a027c51a11d09dd28f Author: Frank14f <1515444314@qq.com> Date: Sun Feb 15 22:36:46 2026 +0800 Initial commit: CelerisLab v0.2.0 with src layout diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..f2ebeeb --- /dev/null +++ b/.gitignore @@ -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 diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..bcaf3b5 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/README.md b/README.md new file mode 100644 index 0000000..17751da --- /dev/null +++ b/README.md @@ -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 +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 diff --git a/configs/config_cuda.json b/configs/config_cuda.json new file mode 100644 index 0000000..3d4c10b --- /dev/null +++ b/configs/config_cuda.json @@ -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 +} \ No newline at end of file diff --git a/configs/config_flowfield.json b/configs/config_flowfield.json new file mode 100644 index 0000000..3f2f42d --- /dev/null +++ b/configs/config_flowfield.json @@ -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"] + } +} \ No newline at end of file diff --git a/configs/config_gym.json b/configs/config_gym.json new file mode 100644 index 0000000..544b7b4 --- /dev/null +++ b/configs/config_gym.json @@ -0,0 +1,3 @@ +{ + +} \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..b702f2c --- /dev/null +++ b/pyproject.toml @@ -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"] diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..3c67a73 --- /dev/null +++ b/setup.py @@ -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', + ], + }, +) \ No newline at end of file diff --git a/src/CelerisLab/__init__.py b/src/CelerisLab/__init__.py new file mode 100644 index 0000000..ba87619 --- /dev/null +++ b/src/CelerisLab/__init__.py @@ -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'] \ No newline at end of file diff --git a/src/CelerisLab/compiler.py b/src/CelerisLab/compiler.py new file mode 100644 index 0000000..c51cbe5 --- /dev/null +++ b/src/CelerisLab/compiler.py @@ -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) \ No newline at end of file diff --git a/src/CelerisLab/configs/config_cuda.json b/src/CelerisLab/configs/config_cuda.json new file mode 100644 index 0000000..3d4c10b --- /dev/null +++ b/src/CelerisLab/configs/config_cuda.json @@ -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 +} \ No newline at end of file diff --git a/src/CelerisLab/configs/config_flowfield.json b/src/CelerisLab/configs/config_flowfield.json new file mode 100644 index 0000000..3f2f42d --- /dev/null +++ b/src/CelerisLab/configs/config_flowfield.json @@ -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"] + } +} \ No newline at end of file diff --git a/src/CelerisLab/configs/config_gym.json b/src/CelerisLab/configs/config_gym.json new file mode 100644 index 0000000..544b7b4 --- /dev/null +++ b/src/CelerisLab/configs/config_gym.json @@ -0,0 +1,3 @@ +{ + +} \ No newline at end of file diff --git a/src/CelerisLab/driver.py b/src/CelerisLab/driver.py new file mode 100644 index 0000000..09832ad --- /dev/null +++ b/src/CelerisLab/driver.py @@ -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() \ No newline at end of file diff --git a/src/CelerisLab/kernels/D2Q9.cu b/src/CelerisLab/kernels/D2Q9.cu new file mode 100644 index 0000000..395d8ed --- /dev/null +++ b/src/CelerisLab/kernels/D2Q9.cu @@ -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; +} \ No newline at end of file diff --git a/src/CelerisLab/kernels/IO.cu b/src/CelerisLab/kernels/IO.cu new file mode 100644 index 0000000..e69de29 diff --git a/src/CelerisLab/kernels/const.h b/src/CelerisLab/kernels/const.h new file mode 100644 index 0000000..e21bcc3 --- /dev/null +++ b/src/CelerisLab/kernels/const.h @@ -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 \ No newline at end of file diff --git a/src/CelerisLab/kernels/kernel.cu b/src/CelerisLab/kernels/kernel.cu new file mode 100644 index 0000000..4bf36f8 --- /dev/null +++ b/src/CelerisLab/kernels/kernel.cu @@ -0,0 +1,222 @@ +// 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/kernels/macros.h b/src/CelerisLab/kernels/macros.h new file mode 100644 index 0000000..7a2466d --- /dev/null +++ b/src/CelerisLab/kernels/macros.h @@ -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 \ No newline at end of file diff --git a/src/CelerisLab/kernels/preproc.cu b/src/CelerisLab/kernels/preproc.cu new file mode 100644 index 0000000..c2b25c5 --- /dev/null +++ b/src/CelerisLab/kernels/preproc.cu @@ -0,0 +1,2 @@ +#include "macros.h" +#include "const.h" diff --git a/src/CelerisLab/preprocess.py b/src/CelerisLab/preprocess.py new file mode 100644 index 0000000..69aadc1 --- /dev/null +++ b/src/CelerisLab/preprocess.py @@ -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 \ No newline at end of file diff --git a/src/CelerisLab/utils.py b/src/CelerisLab/utils.py new file mode 100644 index 0000000..9410a29 --- /dev/null +++ b/src/CelerisLab/utils.py @@ -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 + )