diff --git a/CMakeLists.txt b/CMakeLists.txt index 7919438..bf32dd0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -72,6 +72,12 @@ endif() # Optional Apple Metal GPU backend (defines cyberfluids_metal + CYBERFLUIDS_HAVE_METAL). include(cmake/CyberfluidsMetal.cmake) +# Optional OpenCL GPU backend (defines cyberfluids_opencl + CYBERFLUIDS_HAVE_OPENCL). +include(cmake/CyberfluidsOpenCL.cmake) + +# Optional CUDA GPU backend (defines cyberfluids_cuda + CYBERFLUIDS_HAVE_CUDA). +include(cmake/CyberfluidsCUDA.cmake) + # Optional STL/OBJ geometry import (defines cyberfluids_geometry + CYBERFLUIDS_HAVE_GEOMETRY). include(cmake/CyberfluidsGeometry.cmake) diff --git a/README.md b/README.md index 9acf6d1..fc76f87 100644 --- a/README.md +++ b/README.md @@ -196,6 +196,55 @@ let u = cav.velocity() // [Double], length nx*ny*2 (ro See [`bindings/README.md`](bindings/README.md) for building and linking the bindings. +### Wind tunnel from an STL + +Load a mesh, voxelize it, run external flow past it, and export a VTK for ParaView โ€” the same +core from either language. STL loading needs a geometry-enabled build +(`cmake -B build -DCYBERFLUIDS_GEOMETRY=ON`, which pulls in +[CyberMeshGenerator](https://github.com/CyberdyneCorp/CyberMeshGenerator)); the analytic-sphere +obstacle works in any build. + +#### Python + +```python +import cyberfluids as cf +import numpy as np + +# Voxelize an STL and run flow at Reynolds 100 (inflow speed 0.05 lattice units). +omega = cf.WindTunnel.omega_for_reynolds(u_in=0.05, l_char=48, re=100) +tunnel = cf.WindTunnel.from_stl("car.stl", resolution=48, u_in=0.05, omega=omega) +tunnel.run(15_000) + +vel = tunnel.velocity() # np.ndarray, shape (nx, ny, nz, 3) +print("max speed:", np.linalg.norm(vel, axis=-1).max()) +tunnel.write_vtk("car.vtk") # ParaView: Contour 'solid' + Stream Tracer 'velocity' + +# No mesh? An analytic sphere works with any build: +# t = cf.WindTunnel(128, 64, 64, omega=omega, u_in=0.05); t.set_sphere(36, 32, 32, 10) +``` + +Or straight from the shell via the bundled example: + +```bash +python examples/wind_tunnel.py --stl car.stl --resolution 48 --reynolds 100 --out car.vtk +``` + +#### Swift + +```swift +import Cyberfluids + +// Voxelize an STL and run flow at Reynolds 100. +let omega = WindTunnel.omegaForReynolds(inflow: 0.05, lChar: 48, re: 100) +guard let tunnel = WindTunnel.fromSTL("car.stl", resolution: 48, inflow: 0.05, omega: omega) else { + fatalError("STL support needs a geometry-enabled build (-DCYBERFLUIDS_GEOMETRY=ON)") +} +tunnel.run(steps: 15_000) + +let u = tunnel.velocity() // [Double], length nx*ny*nz*3 (row-major) +tunnel.writeVTK("car.vtk") // open in ParaView +``` + ## ๐Ÿ“Š Status at a glance Cyberfluids runs on **CPU** (`std::execution::par_unseq` with a serial fallback) and **Metal** diff --git a/bindings/swift/Sources/Cyberfluids/WindTunnel.swift b/bindings/swift/Sources/Cyberfluids/WindTunnel.swift new file mode 100644 index 0000000..30c134f --- /dev/null +++ b/bindings/swift/Sources/Cyberfluids/WindTunnel.swift @@ -0,0 +1,73 @@ +import CCyberfluids + +/// A 3D wind tunnel: external flow past an obstacle (D3Q19), driven through the +/// C++ core. The obstacle is an analytic sphere (`setSphere`, always available) +/// or a voxelized STL/OBJ mesh (`fromSTL`, needs a geometry-enabled build). Write +/// the result to VTK and open it in ParaView for streamlines. +public final class WindTunnel { + private let handle: OpaquePointer + public let nx: Int + public let ny: Int + public let nz: Int + + private init(_ h: OpaquePointer) { + handle = h + nx = Int(cf_wind_tunnel_nx(h)) + ny = Int(cf_wind_tunnel_ny(h)) + nz = Int(cf_wind_tunnel_nz(h)) + } + + /// An empty tunnel; add an obstacle with `setSphere`. + public convenience init?(nx: Int, ny: Int, nz: Int, omega: Double, inflow: Double) { + guard let h = cf_wind_tunnel_create(Int64(nx), Int64(ny), Int64(nz), omega, inflow) + else { return nil } + self.init(h) + } + + /// Build a tunnel sized around a voxelized STL/OBJ obstacle. Returns nil when + /// the library lacks geometry support (see `geometryAvailable`) or on read + /// failure. `resolution` = obstacle cells along its longest axis. + public static func fromSTL(_ path: String, resolution: Int = 32, inflow: Double = 0.05, + omega: Double = 1.8, padUp: Int = 16, padDown: Int = 48, + padLat: Int = 12) -> WindTunnel? { + guard let h = cf_wind_tunnel_create_from_stl(path, Int32(resolution), Int32(padUp), + Int32(padDown), Int32(padLat), omega, inflow) + else { return nil } + return WindTunnel(h) + } + + /// True if the library was built with STL/OBJ geometry support. + public static func geometryAvailable() -> Bool { cf_geometry_available() == 1 } + + /// Relaxation rate for a target Reynolds number Re = inflow * lChar / nu. + public static func omegaForReynolds(inflow: Double, lChar: Double, re: Double) -> Double { + let nu = inflow * lChar / re + return 1.0 / (3.0 * nu + 0.5) + } + + deinit { cf_wind_tunnel_destroy(handle) } + + /// Set a solid sphere obstacle (centre and radius in lattice cells). + public func setSphere(cx: Double, cy: Double, cz: Double, radius: Double) { + cf_wind_tunnel_set_sphere(handle, cx, cy, cz, radius) + } + + public func run(steps: Int) { cf_wind_tunnel_run(handle, Int64(steps)) } + + /// Velocity field, row-major, length nx*ny*nz*3: value(x,y,z,c) at [((x*ny+y)*nz+z)*3+c]. + public func velocity() -> [Double] { + var out = [Double](repeating: 0, count: nx * ny * nz * 3) + out.withUnsafeMutableBufferPointer { cf_wind_tunnel_velocity(handle, $0.baseAddress) } + return out + } + + /// Solid-fraction field, row-major, length nx*ny*nz: value(x,y,z) at [(x*ny+y)*nz+z]. + public func solid() -> [Double] { + var out = [Double](repeating: 0, count: nx * ny * nz) + out.withUnsafeMutableBufferPointer { cf_wind_tunnel_solid(handle, $0.baseAddress) } + return out + } + + /// Write a legacy-VTK STRUCTURED_POINTS file (velocity, speed, solid) for ParaView. + public func writeVTK(_ path: String) { cf_wind_tunnel_write_vtk(handle, path) } +} diff --git a/cmake/CyberfluidsCUDA.cmake b/cmake/CyberfluidsCUDA.cmake new file mode 100644 index 0000000..8097fc7 --- /dev/null +++ b/cmake/CyberfluidsCUDA.cmake @@ -0,0 +1,23 @@ +# Optional CUDA GPU backend โ€” a bespoke device-kernel wind-tunnel solver, the +# NVIDIA mirror of the OpenCL path. Authored without nvcc; enable with +# -DCYBERFLUIDS_CUDA=ON on a machine with the CUDA toolkit to compile + validate. +set(CYBERFLUIDS_HAVE_CUDA OFF CACHE INTERNAL "CUDA backend built") + +if(CYBERFLUIDS_CUDA AND CYBERFLUIDS_HAVE_NUMPP) + include(CheckLanguage) + check_language(CUDA) + if(NOT CMAKE_CUDA_COMPILER) + message(STATUS "Cyberfluids: CYBERFLUIDS_CUDA is ON but no CUDA compiler found โ€” skipping") + return() + endif() + enable_language(CUDA) + + add_library(cyberfluids_cuda STATIC src/cuda/wind_tunnel_cuda.cu) + set_target_properties(cyberfluids_cuda PROPERTIES CUDA_STANDARD 20 CUDA_STANDARD_REQUIRED ON) + target_include_directories(cyberfluids_cuda PUBLIC ${CMAKE_SOURCE_DIR}/include) + target_link_libraries(cyberfluids_cuda PUBLIC cyberfluids_core numpp::numpp) + target_compile_definitions(cyberfluids_cuda PUBLIC CYBERFLUIDS_WITH_CUDA) + + set(CYBERFLUIDS_HAVE_CUDA ON CACHE INTERNAL "CUDA backend built") + message(STATUS "Cyberfluids: CUDA backend enabled (device-kernel wind tunnel)") +endif() diff --git a/cmake/CyberfluidsOpenCL.cmake b/cmake/CyberfluidsOpenCL.cmake new file mode 100644 index 0000000..8c5f546 --- /dev/null +++ b/cmake/CyberfluidsOpenCL.cmake @@ -0,0 +1,27 @@ +# Optional OpenCL GPU backend โ€” a bespoke device-kernel wind-tunnel solver +# (real OpenCL C kernels), NOT the CPU-only forEachIndex seam. Builds only when +# CYBERFLUIDS_OPENCL is ON and an OpenCL runtime is found. Mirrors the Metal +# backend module. +set(CYBERFLUIDS_HAVE_OPENCL OFF CACHE INTERNAL "OpenCL backend built") + +if(CYBERFLUIDS_OPENCL AND CYBERFLUIDS_HAVE_NUMPP) + find_package(OpenCL QUIET) + if(NOT OpenCL_FOUND) + message(STATUS "Cyberfluids: CYBERFLUIDS_OPENCL is ON but no OpenCL runtime found โ€” skipping") + return() + endif() + + add_library(cyberfluids_opencl STATIC src/opencl/wind_tunnel_opencl.cpp) + target_include_directories(cyberfluids_opencl PUBLIC ${CMAKE_SOURCE_DIR}/include) + target_compile_features(cyberfluids_opencl PUBLIC cxx_std_20) + target_link_libraries(cyberfluids_opencl PUBLIC cyberfluids_core numpp::numpp + PRIVATE OpenCL::OpenCL) + # Public define so consumers see the solver header; pin to OpenCL 1.2 (macOS + # ceiling) and silence Apple's OpenCL-deprecation warnings. + target_compile_definitions(cyberfluids_opencl + PUBLIC CYBERFLUIDS_WITH_OPENCL + PRIVATE CL_TARGET_OPENCL_VERSION=120 CL_SILENCE_DEPRECATION) + + set(CYBERFLUIDS_HAVE_OPENCL ON CACHE INTERNAL "OpenCL backend built") + message(STATUS "Cyberfluids: OpenCL backend enabled (device-kernel wind tunnel)") +endif() diff --git a/docs/backends.md b/docs/backends.md index 11d6aae..8621fd0 100644 --- a/docs/backends.md +++ b/docs/backends.md @@ -4,39 +4,55 @@ Cyberfluids targets a **single node** and deliberately avoids MPI. Parallelism i intra-node: threads, SIMD, and GPU. The fluid physics never changes when the backend does โ€” see the [hardware-backends spec](../openspec/specs/hardware-backends/spec.md). -## How switching works +## How it actually works (two distinct mechanisms) -The backend is selected at the **loop level** (who launches the kernel), not per cell. Model -code โ€” dynamics and boundary conditions โ€” is written once and is backend-agnostic. The CPU -path lowers to a single `std::for_each(std::execution::par_unseq, ...)`; GPU backends -implement the same kernel-launch interface. +There are **two** parallelism mechanisms, and it's important not to conflate them: + +1. **The CPU seam** โ€” `Backend::forEachIndex(n, f)` runs a **host C++ lambda** `f` across + `[0, n)`, lowering to `std::for_each(std::execution::par_unseq, ...)` (with a serial + fallback). This is a **CPU-only** abstraction: a host closure cannot be marshalled onto a + GPU, so this seam does **not** dispatch to devices. + +2. **GPU device-kernel solvers** โ€” GPU acceleration is a **bespoke solver per API**, with the + hot loop (collide + stream) rewritten as device kernels (MSL / OpenCL C / CUDA C++), data + kept resident on the device across steps. Each is validated against the CPU solver as an + fp64 oracle. The Metal cavity and the OpenCL/CUDA wind tunnel are built this way โ€” they do + not go through `forEachIndex`. ```mermaid flowchart LR - MODEL["Model code
(written once)"] --> SEAM["Backend seam
kernel dispatch"] - SEAM --> CPU["CPU"] - SEAM --> CUDA["CUDA"] - SEAM --> METAL["Metal"] - SEAM --> OCL["OpenCL / SYCL"] + MODEL["Model code"] --> SEAM["forEachIndex seam
(CPU only)"] + SEAM --> CPU["CPU ยท par_unseq / serial"] + GPUK["Device kernels
(collide + stream)"] --> METAL["Metal (MSL)"] + GPUK --> OCL["OpenCL C"] + GPUK --> CUDA["CUDA C++"] + ORACLE["CPU solver = fp64 oracle"] -.validates.-> GPUK ``` ## Backend matrix | Backend | Hardware | Enable flag | Status | |---|---|---|---| -| **CPU** | x86-64 (AVX-512), ARM64 (Neon) | on by default | โœ… Implemented โ€” `std::execution::par_unseq` (guarded by `__cpp_lib_parallel_algorithm`) with a serial fallback | -| **CUDA** | NVIDIA GeForce / Quadro / Tesla | `-DCYBERFLUIDS_CUDA=ON` | ๐Ÿ“‹ Stub โ€” seam locked, no device kernels yet | +| **CPU** | x86-64 (AVX-512), ARM64 (Neon) | on by default | โœ… `std::execution::par_unseq` (guarded by `__cpp_lib_parallel_algorithm`) with a serial fallback | +| **OpenCL** | Apple / AMD / Intel / NVIDIA GPUs | `-DCYBERFLUIDS_OPENCL=ON` | โœ… D3Q19 wind tunnel โ€” device kernels validated on Apple M2 Max GPU (Linf/Uin ~ 1e-5 vs the fp64 CPU oracle; ~1.35 GLUPS) | +| **CUDA** | NVIDIA GeForce / Quadro / Tesla | `-DCYBERFLUIDS_CUDA=ON` | ๐ŸŸก D3Q19 wind-tunnel kernels written (line-for-line mirror of the OpenCL path); compiled + validated on an NVIDIA machine | | **Metal** | Apple Silicon (Mac, iPad) | `-DCYBERFLUIDS_METAL=ON` | โœ… D2Q9 BGK cavity (fp32) โ€” Objective-C++ + runtime MSL kernels | -| **OpenCL / SYCL** | AMD, Intel, integrated GPUs | `-DCYBERFLUIDS_OPENCL=ON` | ๐Ÿ“‹ Stub โ€” seam locked | +| **SYCL** | AMD, Intel, integrated GPUs | `-DCYBERFLUIDS_SYCL=ON` | ๐Ÿ“‹ Planned | + +> Note: the `Cuda`/`Metal`/`OpenCL`/`Sycl` structs in `backend/gpu_stubs.hpp` are **not** the +> acceleration path โ€” they only satisfy the `forEachIndex` contract by running the lambda on +> the host (serial), so a flag-enabled build still runs. Real acceleration lives in the +> device-kernel solvers above. ## Guarantees - Building with **no GPU toolkit** present never breaks the CPU build; GPU code is excluded cleanly. - A GPU run and a CPU run of the same simulation agree within a **documented floating-point - tolerance**. + tolerance** (fp32 device vs fp64 CPU oracle). - No MPI library or launcher is required to build or run. -Status: the **CPU backend is implemented** (โœ…) and is the working path. GPU backends are -**stubs** today โ€” the `Backend::forEachIndex(n, f)` seam is locked so device kernels drop in -without touching model code; wiring real CUDA/Metal/OpenCL kernels is a follow-up change. +Status: the **CPU** path is the always-on baseline. **OpenCL** is a real, validated GPU +wind-tunnel solver (measured ~290ร— vs the serial macOS CPU path on an M2 Max; the fair figure +vs a fully parallel CPU is smaller). **CUDA** mirrors it and is validated on NVIDIA hardware. +**Metal** accelerates the 2D cavity. SYCL is planned. diff --git a/include/cyberfluids/solver/cuda/wind_tunnel_cuda.hpp b/include/cyberfluids/solver/cuda/wind_tunnel_cuda.hpp new file mode 100644 index 0000000..01840c8 --- /dev/null +++ b/include/cyberfluids/solver/cuda/wind_tunnel_cuda.hpp @@ -0,0 +1,53 @@ +#pragma once + +#if defined(CYBERFLUIDS_WITH_CUDA) + +#include +#include +#include +#include +#include + +#include "cyberfluids/geometry/voxel_field.hpp" + +namespace cyberfluids::solver { + +/// CUDA wind tunnel (D3Q19, fp32) โ€” the CUDA mirror of WindTunnelOpenCL. The +/// device kernels (src/cuda/wind_tunnel_cuda.cu) are line-for-line translations +/// of the OpenCL kernels validated on Apple GPU (Linf/Uin ~ 1e-5 vs the fp64 CPU +/// oracle); this path is authored here but compiled/run with nvcc on an NVIDIA +/// machine. Same public API as WindTunnelOpenCL so the same test drives both. +/// See openspec/specs/hardware-backends/spec.md. +class WindTunnelCuda { +public: + WindTunnelCuda(std::int64_t nx, std::int64_t ny, std::int64_t nz, double omega, double uIn); + ~WindTunnelCuda(); + WindTunnelCuda(const WindTunnelCuda&) = delete; + WindTunnelCuda& operator=(const WindTunnelCuda&) = delete; + + void setObstacleSphere(double cx, double cy, double cz, double radius, bool sharp = true); + void setObstacleField(const geometry::VoxelField& obj, int ox, int oy, int oz, + bool sharp = false); + void run(std::int64_t steps); + + std::int64_t nx() const; + std::int64_t ny() const; + std::int64_t nz() const; + + std::array velocity(std::int64_t x, std::int64_t y, std::int64_t z) const; + double solidFraction(std::int64_t x, std::int64_t y, std::int64_t z) const; + std::vector downloadPopulations() const; + void writeVtk(const std::string& path) const; + + /// True if a usable CUDA device is present. + static bool cudaAvailable(); + +private: + void setNs(std::int64_t x, std::int64_t y, std::int64_t z, double ns); + struct Impl; + std::unique_ptr impl_; +}; + +} // namespace cyberfluids::solver + +#endif // CYBERFLUIDS_WITH_CUDA diff --git a/include/cyberfluids/solver/opencl/wind_tunnel_opencl.hpp b/include/cyberfluids/solver/opencl/wind_tunnel_opencl.hpp new file mode 100644 index 0000000..e2d8814 --- /dev/null +++ b/include/cyberfluids/solver/opencl/wind_tunnel_opencl.hpp @@ -0,0 +1,60 @@ +#pragma once + +#if defined(CYBERFLUIDS_WITH_OPENCL) + +#include +#include +#include +#include +#include + +#include "cyberfluids/geometry/voxel_field.hpp" + +namespace cyberfluids::solver { + +/// GPU wind tunnel (D3Q19) running real OpenCL device kernels โ€” a bespoke +/// device-kernel solver (like the Metal cavity), NOT the CPU-only forEachIndex +/// seam. It mirrors the CPU solver::WindTunnel3D API and physics exactly (porous +/// convex-blend collide + far-field/outlet/interior stream), so the CPU solver +/// serves as its fp64 validation oracle. fp32 on the device. Populations stay +/// resident on the GPU across steps; the host copies back only for readout. +/// See docs/backends.md and openspec/specs/hardware-backends/spec.md. +class WindTunnelOpenCL { +public: + WindTunnelOpenCL(std::int64_t nx, std::int64_t ny, std::int64_t nz, double omega, double uIn); + ~WindTunnelOpenCL(); + WindTunnelOpenCL(const WindTunnelOpenCL&) = delete; + WindTunnelOpenCL& operator=(const WindTunnelOpenCL&) = delete; + + /// Analytic sphere obstacle (centre + radius in lattice cells). + void setObstacleSphere(double cx, double cy, double cz, double radius, bool sharp = true); + /// Embed a voxelized object's solid fraction at offset (ox,oy,oz). + void setObstacleField(const geometry::VoxelField& obj, int ox, int oy, int oz, + bool sharp = false); + + void run(std::int64_t steps); + + std::int64_t nx() const; + std::int64_t ny() const; + std::int64_t nz() const; + + /// Flux velocity (accumulated in fp64 on the host so the read-back adds no error). + std::array velocity(std::int64_t x, std::int64_t y, std::int64_t z) const; + double solidFraction(std::int64_t x, std::int64_t y, std::int64_t z) const; + /// Raw {q, ncells} populations copied off the device (for a determinism check). + std::vector downloadPopulations() const; + void writeVtk(const std::string& path) const; + + /// True if an OpenCL platform with a device is usable on this machine. + static bool openclAvailable(); + +private: + void setNs(std::int64_t x, std::int64_t y, std::int64_t z, double ns); + + struct Impl; + std::unique_ptr impl_; +}; + +} // namespace cyberfluids::solver + +#endif // CYBERFLUIDS_WITH_OPENCL diff --git a/openspec/changes/add-gpu-wind-tunnel/.openspec.yaml b/openspec/changes/add-gpu-wind-tunnel/.openspec.yaml new file mode 100644 index 0000000..e089cfa --- /dev/null +++ b/openspec/changes/add-gpu-wind-tunnel/.openspec.yaml @@ -0,0 +1,2 @@ +schema: spec-driven +created: 2026-07-05 diff --git a/openspec/changes/add-gpu-wind-tunnel/README.md b/openspec/changes/add-gpu-wind-tunnel/README.md new file mode 100644 index 0000000..4349e84 --- /dev/null +++ b/openspec/changes/add-gpu-wind-tunnel/README.md @@ -0,0 +1,3 @@ +# add-gpu-wind-tunnel + +OpenCL + CUDA device-kernel wind-tunnel solvers diff --git a/openspec/changes/add-gpu-wind-tunnel/proposal.md b/openspec/changes/add-gpu-wind-tunnel/proposal.md new file mode 100644 index 0000000..dbc4b34 --- /dev/null +++ b/openspec/changes/add-gpu-wind-tunnel/proposal.md @@ -0,0 +1,36 @@ +# GPU-accelerate the wind tunnel (OpenCL + CUDA device kernels) + +## Why + +Audit finding: the WindTunnel (and every solver through `Backend::forEachIndex`) runs on the +CPU. The GPU "backends" in `gpu_stubs.hpp` merely run the host lambda serially โ€” the seam takes +a host callable and cannot dispatch to a device, so it was never GPU-capable. Only the bespoke +Metal cavity used the GPU. This change delivers real GPU acceleration for the 3D wind tunnel and +corrects the misleading "seam is GPU-ready" story. + +## What Changes + +- **OpenCL wind tunnel** (`solver::WindTunnelOpenCL`, `cmake/CyberfluidsOpenCL.cmake`): two + device kernels โ€” porous convex-blend `collide` and `stream_bc` (inlet/wall/outlet/interior) โ€” + faithful ports of the CPU solver, with populations resident on the GPU (ping-pong buffers) and + copy-back only for readout. Validated on this machine's **Apple M2 Max GPU** against the fp64 + CPU `WindTunnel3D` oracle: `Linf/Uin โ‰ˆ 1e-5`, deterministic, no-slip preserved; ~1.35 GLUPS. +- **CUDA mirror** (`solver::WindTunnelCuda`, `src/cuda/wind_tunnel_cuda.cu`, + `cmake/CyberfluidsCUDA.cmake`): line-for-line translation of the OpenCL kernels; authored here, + compiled + validated on an NVIDIA machine (`-DCYBERFLUIDS_CUDA=ON`). +- **Tests**: `opencl_wind_tunnel` (runs here, skips cleanly without a device) and + `cuda_wind_tunnel` (runs on NVIDIA) โ€” both compare GPU fp32 vs CPU fp64 within tolerance, + plus determinism and no-slip. +- **Docs/spec honesty**: clarify that `forEachIndex` is CPU-only and GPU = bespoke device-kernel + solvers; update `docs/backends.md` and the hardware-backends spec. + +## Non-goals + +- Routing GPU through `forEachIndex` (architecturally impossible for a host lambda). +- A GPU path in the C ABI / Python binding (follow-up; the C++ solver + tests validate the kernels). +- SYCL; generalizing Metal beyond the 2D cavity. + +## Notes + +- The ~290ร— figure is vs the *serial* macOS CPU path (AppleClang libc++ lacks `par_unseq`); vs a + fully parallel CPU the fair speedup is smaller. The GPU throughput (~1.35 GLUPS) is the honest metric. diff --git a/openspec/changes/add-gpu-wind-tunnel/specs/hardware-backends/spec.md b/openspec/changes/add-gpu-wind-tunnel/specs/hardware-backends/spec.md new file mode 100644 index 0000000..44b208a --- /dev/null +++ b/openspec/changes/add-gpu-wind-tunnel/specs/hardware-backends/spec.md @@ -0,0 +1,44 @@ +# hardware-backends (delta) + +## MODIFIED Requirements + +### Requirement: Backend abstraction seam +The library SHALL expose a compute seam so per-cell CPU work is written once and runs across +the index range via `Backend::forEachIndex(n, f)`. This seam is **CPU-only**: `f` is a host +callable, so it lowers to an STL parallel-for and CANNOT dispatch to a GPU. GPU acceleration +SHALL therefore be provided as **bespoke device-kernel solvers** (the hot collide + stream loop +rewritten in device code, data kept resident on the device), NOT through this seam. The GPU +stub structs that satisfy the seam by running the host callable serially SHALL NOT be presented +as the acceleration path. + +#### Scenario: CPU kernel runs across the seam +- **WHEN** a per-cell operation is expressed as `Backend::forEachIndex(n, f)` on the CPU backend +- **THEN** it SHALL run across `[0, n)` using STL parallel algorithms (or a serial fallback) + +#### Scenario: GPU acceleration bypasses the seam +- **WHEN** a solver is GPU-accelerated +- **THEN** it SHALL use dedicated device kernels validated against the CPU solver, not + `forEachIndex` dispatch of a host callable + +### Requirement: Optional GPU backends +The library SHALL provide optional, independently selectable GPU backends โ€” NVIDIA CUDA, Apple +Metal, and OpenCL/SYCL โ€” enabled via build flags; absence of a GPU toolkit SHALL NOT break the CPU +build. Because GPU paths run in single precision, GPU-vs-CPU agreement is bounded by an fp32 +tolerance rather than being bit-exact. The Apple Metal backend SHALL provide a D2Q9 BGK +lid-driven cavity GPU solver against the system Metal framework. The library SHALL ALSO provide a +**3D wind-tunnel (D3Q19) GPU solver** with **OpenCL** device kernels (validated against the fp64 +CPU `WindTunnel3D` oracle) and a **CUDA** mirror of the same kernels; each GPU solver SHALL keep +populations resident on the device across steps and copy back only for readout. + +#### Scenario: Build with a GPU backend disabled +- **WHEN** the library is built without any GPU backend flag enabled +- **THEN** the CPU backend SHALL build and run, and GPU backends SHALL be cleanly excluded + +#### Scenario: GPU-accelerated run matches CPU +- **WHEN** the same D2Q9 lid-driven cavity (Metal) or 3D wind tunnel (OpenCL/CUDA) is run on the + GPU and on the CPU with identical parameters +- **THEN** the steady-state velocity field SHALL agree within a documented fp32 tolerance + +#### Scenario: Deterministic GPU run +- **WHEN** the same GPU simulation is run twice on the same device +- **THEN** the two runs SHALL produce bit-identical population fields diff --git a/openspec/changes/add-gpu-wind-tunnel/tasks.md b/openspec/changes/add-gpu-wind-tunnel/tasks.md new file mode 100644 index 0000000..58d09d3 --- /dev/null +++ b/openspec/changes/add-gpu-wind-tunnel/tasks.md @@ -0,0 +1,22 @@ +# Tasks + +## 1. OpenCL (validated here) +- [x] 1.1 `cmake/CyberfluidsOpenCL.cmake` (option, find OpenCL, cyberfluids_opencl lib) +- [x] 1.2 `WindTunnelOpenCL` โ€” collide + stream_bc kernels, resident buffers, ping-pong, fp64 read-back +- [x] 1.3 `tests/test_opencl_wind_tunnel.cpp`: GPU-vs-CPU (Linf/L2/mean), determinism, no-slip; skips w/o device +- [x] 1.4 Validated on Apple M2 Max: Linf/Uin ~ 1e-5, ~1.35 GLUPS + +## 2. CUDA mirror (test on NVIDIA) +- [x] 2.1 `cmake/CyberfluidsCUDA.cmake` (enable_language(CUDA), guarded) +- [x] 2.2 `src/cuda/wind_tunnel_cuda.cu` โ€” line-for-line kernel translation + host class +- [x] 2.3 `tests/test_cuda_wind_tunnel.cpp` (mirror; skips w/o CUDA device) +- [ ] 2.4 Compile + validate on an NVIDIA machine (`-DCYBERFLUIDS_CUDA=ON`) โ† run there + +## 3. Docs & spec +- [x] 3.1 Correct `docs/backends.md` (forEachIndex is CPU-only; GPU = device-kernel solvers) +- [x] 3.2 hardware-backends spec delta (seam CPU-only; OpenCL/CUDA wind tunnel) + +## 4. Review & archive +- [ ] 4.1 openspec validate --all --strict +- [ ] 4.2 Adversarial review (GPU memory/indexing; CUDA correctness) + fix +- [ ] 4.3 Archive diff --git a/src/cuda/wind_tunnel_cuda.cu b/src/cuda/wind_tunnel_cuda.cu new file mode 100644 index 0000000..1ee962b --- /dev/null +++ b/src/cuda/wind_tunnel_cuda.cu @@ -0,0 +1,269 @@ +// CUDA mirror of the OpenCL wind tunnel. The two __global__ kernels below are +// line-for-line translations of the OpenCL C kernels in +// src/opencl/wind_tunnel_opencl.cpp (validated on Apple GPU to ~1e-5 vs the fp64 +// CPU oracle). Authored on a machine without nvcc; compile + validate on an +// NVIDIA box with -DCYBERFLUIDS_CUDA=ON and the ctest `cuda_wind_tunnel`. + +#if defined(CYBERFLUIDS_WITH_CUDA) + +#include "cyberfluids/solver/cuda/wind_tunnel_cuda.hpp" + +#include +#include +#include +#include +#include + +#include + +#include "cyberfluids/io/vtk_writer.hpp" + +namespace cyberfluids::solver { + +namespace { + +constexpr int kCx[19] = {0, 1, -1, 0, 0, 0, 0, 1, -1, 1, -1, 1, -1, 1, -1, 0, 0, 0, 0}; +constexpr int kCy[19] = {0, 0, 0, 1, -1, 0, 0, 1, -1, -1, 1, 0, 0, 0, 0, 1, -1, 1, -1}; +constexpr int kCz[19] = {0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 1, -1, -1, 1, 1, -1, -1, 1}; +inline double weight(int i) { + return i == 0 ? 1.0 / 3.0 : (i < 7 ? 1.0 / 18.0 : 1.0 / 36.0); +} + +void cudaCheck(cudaError_t e, const char* what) { + if (e != cudaSuccess) + throw std::runtime_error(std::string("CUDA error at ") + what + ": " + + cudaGetErrorString(e)); +} + +// Device stencil constants (match include/cyberfluids/core/descriptors.hpp). +__constant__ int dCx[19] = {0, 1, -1, 0, 0, 0, 0, 1, -1, 1, -1, 1, -1, 1, -1, 0, 0, 0, 0}; +__constant__ int dCy[19] = {0, 0, 0, 1, -1, 0, 0, 1, -1, -1, 1, 0, 0, 0, 0, 1, -1, 1, -1}; +__constant__ int dCz[19] = {0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 1, -1, -1, 1, 1, -1, -1, 1}; +__constant__ int dOpp[19] = {0, 2, 1, 4, 3, 6, 5, 8, 7, 10, 9, 12, 11, 14, 13, 16, 15, 18, 17}; +__constant__ float dWt[19] = {0.333333333f, + 0.055555556f, 0.055555556f, 0.055555556f, 0.055555556f, 0.055555556f, 0.055555556f, + 0.027777778f, 0.027777778f, 0.027777778f, 0.027777778f, 0.027777778f, 0.027777778f, + 0.027777778f, 0.027777778f, 0.027777778f, 0.027777778f, 0.027777778f, 0.027777778f}; + +} // namespace + +// Porous convex-blend collide, in place (force = 0). +__global__ void collideKernel(float* f, const float* ns, long n, float omega) { + long g = (long)blockIdx.x * blockDim.x + threadIdx.x; + if (g >= n) return; + float fi[19]; + float rho = 0.f, jx = 0.f, jy = 0.f, jz = 0.f; + for (int i = 0; i < 19; ++i) { + float v = f[(long)i * n + g]; + fi[i] = v; rho += v; + jx += v * dCx[i]; jy += v * dCy[i]; jz += v * dCz[i]; + } + float invRho = 1.f / rho; + float ux = jx * invRho, uy = jy * invRho, uz = jz * invRho; + float uSqr = ux * ux + uy * uy + uz * uz; + float s = ns[g], fluid = 1.f - s; + for (int i = 0; i < 19; ++i) { + float ciu = dCx[i] * ux + dCy[i] * uy + dCz[i] * uz; + float feq = dWt[i] * rho * (1.f + 3.f * ciu + 4.5f * ciu * ciu - 1.5f * uSqr); + float bgk = fi[i] - omega * (fi[i] - feq); + f[(long)i * n + g] = fluid * bgk + s * fi[dOpp[i]]; + } +} + +// Stream + boundary conditions: outlet copy / free-stream Dirichlet / interior pull. +__global__ void streamKernel(const float* src, float* dst, const float* ns, float uIn, int nx, + int ny, int nz, long n) { + long g = (long)blockIdx.x * blockDim.x + threadIdx.x; + if (g >= n) return; + long z = g % nz, y = (g / nz) % ny, x = g / ((long)ny * nz); + if (x == nx - 1) { + long up = (((long)(nx - 2) * ny + y) * nz + z); + for (int i = 0; i < 19; ++i) dst[(long)i * n + g] = src[(long)i * n + up]; + return; + } + if (x == 0 || y == 0 || y == ny - 1 || z == 0 || z == nz - 1) { + if (ns[g] > 0.5f) { + for (int i = 0; i < 19; ++i) dst[(long)i * n + g] = src[(long)i * n + g]; + } else { + float uSqr = uIn * uIn; + for (int i = 0; i < 19; ++i) { + float ciu = dCx[i] * uIn; + dst[(long)i * n + g] = dWt[i] * (1.f + 3.f * ciu + 4.5f * ciu * ciu - 1.5f * uSqr); + } + } + return; + } + for (int i = 0; i < 19; ++i) { + long sx = x - dCx[i], sy = y - dCy[i], sz = z - dCz[i]; + dst[(long)i * n + g] = src[(long)i * n + ((sx * ny + sy) * nz + sz)]; + } +} + +struct WindTunnelCuda::Impl { + std::int64_t nx, ny, nz, ncells; + float omega, uIn; + float *fA = nullptr, *fB = nullptr, *nsBuf = nullptr; + float *cur = nullptr, *other = nullptr; + std::vector hf, hns; + mutable std::vector cache; + mutable bool cacheValid = false; + + std::int64_t idx(std::int64_t x, std::int64_t y, std::int64_t z) const { + return (x * ny + y) * nz + z; + } + float feq(int i, float u) const { + const float ciu = static_cast(kCx[i]) * u; + return static_cast(weight(i)) * (1.f + 3.f * ciu + 4.5f * ciu * ciu - 1.5f * u * u); + } + void seedInitial() { + hf.assign(static_cast(19 * ncells), 0.f); + hns.assign(static_cast(ncells), 0.f); + for (std::int64_t c = 0; c < ncells; ++c) + for (int i = 0; i < 19; ++i) hf[static_cast(i * ncells + c)] = feq(i, uIn); + } + void uploadState() { + cudaCheck(cudaMemcpy(fA, hf.data(), sizeof(float) * hf.size(), cudaMemcpyHostToDevice), + "upload fA"); + cudaCheck(cudaMemcpy(nsBuf, hns.data(), sizeof(float) * hns.size(), cudaMemcpyHostToDevice), + "upload ns"); + cur = fA; + other = fB; + cacheValid = false; + } + void ensureDownloaded() const { + if (cacheValid) return; + cache.resize(static_cast(19 * ncells)); + cudaCheck(cudaMemcpy(cache.data(), cur, sizeof(float) * cache.size(), + cudaMemcpyDeviceToHost), + "download populations"); + cacheValid = true; + } +}; + +WindTunnelCuda::WindTunnelCuda(std::int64_t nx, std::int64_t ny, std::int64_t nz, double omega, + double uIn) + : impl_(std::make_unique()) { + if (nx < 3 || ny < 3 || nz < 3) + throw std::invalid_argument("WindTunnelCuda: each extent must be >= 3"); + auto& im = *impl_; + im.nx = nx; im.ny = ny; im.nz = nz; im.ncells = nx * ny * nz; + im.omega = static_cast(omega); + im.uIn = static_cast(uIn); + const std::size_t bytes = sizeof(float) * static_cast(19 * im.ncells); + cudaCheck(cudaMalloc(&im.fA, bytes), "malloc fA"); + cudaCheck(cudaMalloc(&im.fB, bytes), "malloc fB"); + cudaCheck(cudaMalloc(&im.nsBuf, sizeof(float) * static_cast(im.ncells)), + "malloc ns"); + im.seedInitial(); + im.uploadState(); +} + +WindTunnelCuda::~WindTunnelCuda() { + auto& im = *impl_; + if (im.fA) cudaFree(im.fA); + if (im.fB) cudaFree(im.fB); + if (im.nsBuf) cudaFree(im.nsBuf); +} + +void WindTunnelCuda::setNs(std::int64_t x, std::int64_t y, std::int64_t z, double ns) { + auto& im = *impl_; + const std::int64_t c = im.idx(x, y, z); + im.hns[static_cast(c)] = static_cast(ns); + if (ns > 0.5) + for (int i = 0; i < 19; ++i) + im.hf[static_cast(i * im.ncells + c)] = im.feq(i, 0.f); +} + +void WindTunnelCuda::setObstacleSphere(double cx, double cy, double cz, double radius, bool sharp) { + auto& im = *impl_; + for (std::int64_t x = 0; x < im.nx; ++x) + for (std::int64_t y = 0; y < im.ny; ++y) + for (std::int64_t z = 0; z < im.nz; ++z) { + const double dx = x - cx, dy = y - cy, dz = z - cz; + const double phi = std::sqrt(dx * dx + dy * dy + dz * dz) - radius; + setNs(x, y, z, + sharp ? (phi < 0.0 ? 1.0 : 0.0) : geometry::sdfToSolidFraction(phi, 1.0)); + } + im.uploadState(); +} + +void WindTunnelCuda::setObstacleField(const geometry::VoxelField& obj, int ox, int oy, int oz, + bool sharp) { + auto& im = *impl_; + for (int x = 0; x < obj.nx; ++x) + for (int y = 0; y < obj.ny; ++y) + for (int z = 0; z < obj.nz; ++z) { + const std::int64_t tx = ox + x, ty = oy + y, tz = oz + z; + if (tx < 0 || tx >= im.nx || ty < 0 || ty >= im.ny || tz < 0 || tz >= im.nz) + continue; + const double phi = obj.distance[static_cast(obj.index(x, y, z))]; + setNs(tx, ty, tz, + sharp ? (phi < 0.0 ? 1.0 : 0.0) : geometry::sdfToSolidFraction(phi, obj.spacing)); + } + im.uploadState(); +} + +void WindTunnelCuda::run(std::int64_t steps) { + auto& im = *impl_; + const int threads = 256; + const int blocks = static_cast((im.ncells + threads - 1) / threads); + const long n = static_cast(im.ncells); + for (std::int64_t s = 0; s < steps; ++s) { + collideKernel<<>>(im.cur, im.nsBuf, n, im.omega); + streamKernel<<>>(im.cur, im.other, im.nsBuf, im.uIn, + static_cast(im.nx), static_cast(im.ny), + static_cast(im.nz), n); + std::swap(im.cur, im.other); + } + cudaCheck(cudaDeviceSynchronize(), "run sync"); + im.cacheValid = false; +} + +std::int64_t WindTunnelCuda::nx() const { return impl_->nx; } +std::int64_t WindTunnelCuda::ny() const { return impl_->ny; } +std::int64_t WindTunnelCuda::nz() const { return impl_->nz; } + +std::array WindTunnelCuda::velocity(std::int64_t x, std::int64_t y, + std::int64_t z) const { + const auto& im = *impl_; + im.ensureDownloaded(); + const std::int64_t c = im.idx(x, y, z); + double rho = 0, jx = 0, jy = 0, jz = 0; + for (int i = 0; i < 19; ++i) { + const double v = im.cache[static_cast(i * im.ncells + c)]; + rho += v; + jx += v * kCx[i]; jy += v * kCy[i]; jz += v * kCz[i]; + } + const double invRho = 1.0 / rho; + return {jx * invRho, jy * invRho, jz * invRho}; +} + +double WindTunnelCuda::solidFraction(std::int64_t x, std::int64_t y, std::int64_t z) const { + return impl_->hns[static_cast(impl_->idx(x, y, z))]; +} + +std::vector WindTunnelCuda::downloadPopulations() const { + impl_->ensureDownloaded(); + return impl_->cache; +} + +void WindTunnelCuda::writeVtk(const std::string& path) const { + io::VtkStructuredWriter<3> w(impl_->nx, impl_->ny, impl_->nz); + using Coord = io::VtkStructuredWriter<3>::Coord; + w.addVector("velocity", [this](Coord c) { return velocity(c[0], c[1], c[2]); }); + w.addScalar("speed", [this](Coord c) { + const auto u = velocity(c[0], c[1], c[2]); + return std::sqrt(u[0] * u[0] + u[1] * u[1] + u[2] * u[2]); + }); + w.addScalar("solid", [this](Coord c) { return solidFraction(c[0], c[1], c[2]); }); + w.write(path); +} + +bool WindTunnelCuda::cudaAvailable() { + int count = 0; + return cudaGetDeviceCount(&count) == cudaSuccess && count > 0; +} + +} // namespace cyberfluids::solver + +#endif // CYBERFLUIDS_WITH_CUDA diff --git a/src/opencl/wind_tunnel_opencl.cpp b/src/opencl/wind_tunnel_opencl.cpp new file mode 100644 index 0000000..c6aa630 --- /dev/null +++ b/src/opencl/wind_tunnel_opencl.cpp @@ -0,0 +1,371 @@ +#if defined(CYBERFLUIDS_WITH_OPENCL) + +#include "cyberfluids/solver/opencl/wind_tunnel_opencl.hpp" + +#include +#include +#include +#include +#include +#include + +#ifdef __APPLE__ +#include +#else +#include +#endif + +#include "cyberfluids/io/vtk_writer.hpp" + +namespace cyberfluids::solver { + +namespace { + +// D3Q19 stencil on the host (must match include/cyberfluids/core/descriptors.hpp). +constexpr int kCx[19] = {0, 1, -1, 0, 0, 0, 0, 1, -1, 1, -1, 1, -1, 1, -1, 0, 0, 0, 0}; +constexpr int kCy[19] = {0, 0, 0, 1, -1, 0, 0, 1, -1, -1, 1, 0, 0, 0, 0, 1, -1, 1, -1}; +constexpr int kCz[19] = {0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 1, -1, -1, 1, 1, -1, -1, 1}; +inline double weight(int i) { + return i == 0 ? 1.0 / 3.0 : (i < 7 ? 1.0 / 18.0 : 1.0 / 36.0); +} + +// Two OpenCL C kernels, faithful ports of PorousForcedBGKdynamics::collide +// (force=0) and WindTunnel3D::streamWithBoundaries. Same arithmetic shape as the +// fp64 CPU oracle so fp32 results agree within tolerance. +const char* kKernelSource = R"CL( +constant int cx[19] = {0, 1,-1,0,0,0,0, 1,-1,1,-1, 1,-1,1,-1, 0,0,0,0}; +constant int cy[19] = {0, 0,0,1,-1,0,0, 1,-1,-1,1, 0,0,0,0, 1,-1,1,-1}; +constant int cz[19] = {0, 0,0,0,0,1,-1, 0,0,0,0, 1,-1,-1,1, 1,-1,-1,1}; +constant int opp[19] = {0,2,1,4,3,6,5,8,7,10,9,12,11,14,13,16,15,18,17}; +constant float wt[19] = {0.333333333f, + 0.055555556f,0.055555556f,0.055555556f,0.055555556f,0.055555556f,0.055555556f, + 0.027777778f,0.027777778f,0.027777778f,0.027777778f,0.027777778f,0.027777778f, + 0.027777778f,0.027777778f,0.027777778f,0.027777778f,0.027777778f,0.027777778f}; + +// Porous convex-blend collide, in place (force = 0). +kernel void collide(global float* f, global const float* ns, const long n, const float omega) { + long g = get_global_id(0); + if (g >= n) return; + float fi[19]; + float rho = 0.f, jx = 0.f, jy = 0.f, jz = 0.f; + for (int i = 0; i < 19; ++i) { + float v = f[(long)i * n + g]; // snapshot: bounce-back reads pre-collision f[opp] + fi[i] = v; rho += v; + jx += v * cx[i]; jy += v * cy[i]; jz += v * cz[i]; + } + float invRho = 1.f / rho; + float ux = jx * invRho, uy = jy * invRho, uz = jz * invRho; + float uSqr = ux * ux + uy * uy + uz * uz; + float s = ns[g], fluid = 1.f - s; + for (int i = 0; i < 19; ++i) { + float ciu = cx[i] * ux + cy[i] * uy + cz[i] * uz; + float feq = wt[i] * rho * (1.f + 3.f * ciu + 4.5f * ciu * ciu - 1.5f * uSqr); + float bgk = fi[i] - omega * (fi[i] - feq); + f[(long)i * n + g] = fluid * bgk + s * fi[opp[i]]; + } +} + +// Stream + boundary conditions: outlet copy / free-stream Dirichlet / interior pull. +kernel void stream_bc(global const float* src, global float* dst, global const float* ns, + const float uIn, const int nx, const int ny, const int nz, const long n) { + long g = get_global_id(0); + if (g >= n) return; + long z = g % nz, y = (g / nz) % ny, x = g / ((long)ny * nz); + if (x == nx - 1) { // outlet: zero-gradient copy of upstream plane + long up = (((long)(nx - 2) * ny + y) * nz + z); + for (int i = 0; i < 19; ++i) dst[(long)i * n + g] = src[(long)i * n + up]; + return; + } + if (x == 0 || y == 0 || y == ny - 1 || z == 0 || z == nz - 1) { + if (ns[g] > 0.5f) { // solid on a face keeps its collided value + for (int i = 0; i < 19; ++i) dst[(long)i * n + g] = src[(long)i * n + g]; + } else { // free-stream feq(rho=1, u=(uIn,0,0)) + float uSqr = uIn * uIn; + for (int i = 0; i < 19; ++i) { + float ciu = cx[i] * uIn; + dst[(long)i * n + g] = wt[i] * (1.f + 3.f * ciu + 4.5f * ciu * ciu - 1.5f * uSqr); + } + } + return; + } + for (int i = 0; i < 19; ++i) { // interior pull-stream + long sx = x - cx[i], sy = y - cy[i], sz = z - cz[i]; + dst[(long)i * n + g] = src[(long)i * n + ((sx * ny + sy) * nz + sz)]; + } +} +)CL"; + +void clCheck(cl_int err, const char* what) { + if (err != CL_SUCCESS) + throw std::runtime_error(std::string("OpenCL error ") + std::to_string(err) + " at " + what); +} + +// Pick the first GPU device (fall back to any device) of the first platform. +bool pickDevice(cl_platform_id& plat, cl_device_id& dev) { + cl_uint nplat = 0; + if (clGetPlatformIDs(0, nullptr, &nplat) != CL_SUCCESS || nplat == 0) return false; + std::vector plats(nplat); + clGetPlatformIDs(nplat, plats.data(), nullptr); + for (cl_platform_id p : plats) { + cl_uint nd = 0; + if (clGetDeviceIDs(p, CL_DEVICE_TYPE_GPU, 1, &dev, &nd) == CL_SUCCESS && nd > 0) { + plat = p; + return true; + } + } + for (cl_platform_id p : plats) { + cl_uint nd = 0; + if (clGetDeviceIDs(p, CL_DEVICE_TYPE_ALL, 1, &dev, &nd) == CL_SUCCESS && nd > 0) { + plat = p; + return true; + } + } + return false; +} + +} // namespace + +struct WindTunnelOpenCL::Impl { + std::int64_t nx, ny, nz, ncells; + float omega, uIn; + cl_context ctx = nullptr; + cl_command_queue queue = nullptr; + cl_program program = nullptr; + cl_kernel collideK = nullptr, streamK = nullptr; + cl_mem fA = nullptr, fB = nullptr, nsBuf = nullptr; + cl_mem cur = nullptr, other = nullptr; + std::vector hf; // host staging populations {q, ncells} + std::vector hns; // host solid fraction {ncells} + mutable std::vector cache; + mutable bool cacheValid = false; + + std::int64_t idx(std::int64_t x, std::int64_t y, std::int64_t z) const { + return (x * ny + y) * nz + z; + } + + // Rest / free-stream equilibrium at u = (u,0,0), rho = 1 for population i. + float feq(int i, float u) const { + const float ciu = static_cast(kCx[i]) * u; + const float uSqr = u * u; + return static_cast(weight(i)) * (1.f + 3.f * ciu + 4.5f * ciu * ciu - 1.5f * uSqr); + } + + void seedInitial() { + hf.assign(static_cast(19 * ncells), 0.f); + hns.assign(static_cast(ncells), 0.f); + for (std::int64_t c = 0; c < ncells; ++c) + for (int i = 0; i < 19; ++i) + hf[static_cast(i * ncells + c)] = feq(i, uIn); + } + + void uploadState() { + clCheck(clEnqueueWriteBuffer(queue, fA, CL_TRUE, 0, + sizeof(float) * hf.size(), hf.data(), 0, nullptr, nullptr), + "write fA"); + clCheck(clEnqueueWriteBuffer(queue, nsBuf, CL_TRUE, 0, + sizeof(float) * hns.size(), hns.data(), 0, nullptr, nullptr), + "write ns"); + cur = fA; + other = fB; + cacheValid = false; + } + + void ensureDownloaded() const { + if (cacheValid) return; + cache.resize(static_cast(19 * ncells)); + clCheck(clEnqueueReadBuffer(queue, cur, CL_TRUE, 0, sizeof(float) * cache.size(), + cache.data(), 0, nullptr, nullptr), + "read populations"); + cacheValid = true; + } +}; + +WindTunnelOpenCL::WindTunnelOpenCL(std::int64_t nx, std::int64_t ny, std::int64_t nz, double omega, + double uIn) + : impl_(std::make_unique()) { + if (nx < 3 || ny < 3 || nz < 3) + throw std::invalid_argument("WindTunnelOpenCL: each extent must be >= 3"); + auto& im = *impl_; + im.nx = nx; + im.ny = ny; + im.nz = nz; + im.ncells = nx * ny * nz; + im.omega = static_cast(omega); + im.uIn = static_cast(uIn); + + cl_platform_id plat; + cl_device_id dev; + if (!pickDevice(plat, dev)) throw std::runtime_error("WindTunnelOpenCL: no OpenCL device"); + + cl_int err = CL_SUCCESS; + im.ctx = clCreateContext(nullptr, 1, &dev, nullptr, nullptr, &err); + clCheck(err, "clCreateContext"); + im.queue = clCreateCommandQueue(im.ctx, dev, 0, &err); + clCheck(err, "clCreateCommandQueue"); + + im.program = clCreateProgramWithSource(im.ctx, 1, &kKernelSource, nullptr, &err); + clCheck(err, "clCreateProgramWithSource"); + if (clBuildProgram(im.program, 1, &dev, "", nullptr, nullptr) != CL_SUCCESS) { + std::size_t len = 0; + clGetProgramBuildInfo(im.program, dev, CL_PROGRAM_BUILD_LOG, 0, nullptr, &len); + std::string log(len, '\0'); + clGetProgramBuildInfo(im.program, dev, CL_PROGRAM_BUILD_LOG, len, log.data(), nullptr); + throw std::runtime_error("WindTunnelOpenCL kernel build failed:\n" + log); + } + im.collideK = clCreateKernel(im.program, "collide", &err); + clCheck(err, "clCreateKernel collide"); + im.streamK = clCreateKernel(im.program, "stream_bc", &err); + clCheck(err, "clCreateKernel stream_bc"); + + const std::size_t bytes = sizeof(float) * static_cast(19 * im.ncells); + im.fA = clCreateBuffer(im.ctx, CL_MEM_READ_WRITE, bytes, nullptr, &err); + clCheck(err, "buffer fA"); + im.fB = clCreateBuffer(im.ctx, CL_MEM_READ_WRITE, bytes, nullptr, &err); + clCheck(err, "buffer fB"); + im.nsBuf = clCreateBuffer(im.ctx, CL_MEM_READ_WRITE, + sizeof(float) * static_cast(im.ncells), nullptr, &err); + clCheck(err, "buffer ns"); + + im.seedInitial(); + im.uploadState(); +} + +WindTunnelOpenCL::~WindTunnelOpenCL() { + auto& im = *impl_; + if (im.collideK) clReleaseKernel(im.collideK); + if (im.streamK) clReleaseKernel(im.streamK); + if (im.program) clReleaseProgram(im.program); + if (im.fA) clReleaseMemObject(im.fA); + if (im.fB) clReleaseMemObject(im.fB); + if (im.nsBuf) clReleaseMemObject(im.nsBuf); + if (im.queue) clReleaseCommandQueue(im.queue); + if (im.ctx) clReleaseContext(im.ctx); +} + +void WindTunnelOpenCL::setObstacleSphere(double cx, double cy, double cz, double radius, + bool sharp) { + auto& im = *impl_; + for (std::int64_t x = 0; x < im.nx; ++x) + for (std::int64_t y = 0; y < im.ny; ++y) + for (std::int64_t z = 0; z < im.nz; ++z) { + const double dx = x - cx, dy = y - cy, dz = z - cz; + const double phi = std::sqrt(dx * dx + dy * dy + dz * dz) - radius; + const double ns = + sharp ? (phi < 0.0 ? 1.0 : 0.0) : geometry::sdfToSolidFraction(phi, 1.0); + setNs(x, y, z, ns); + } + im.uploadState(); +} + +void WindTunnelOpenCL::setObstacleField(const geometry::VoxelField& obj, int ox, int oy, int oz, + bool sharp) { + auto& im = *impl_; + for (int x = 0; x < obj.nx; ++x) + for (int y = 0; y < obj.ny; ++y) + for (int z = 0; z < obj.nz; ++z) { + const std::int64_t tx = ox + x, ty = oy + y, tz = oz + z; + if (tx < 0 || tx >= im.nx || ty < 0 || ty >= im.ny || tz < 0 || tz >= im.nz) + continue; + const double phi = obj.distance[static_cast(obj.index(x, y, z))]; + const double ns = sharp ? (phi < 0.0 ? 1.0 : 0.0) + : geometry::sdfToSolidFraction(phi, obj.spacing); + setNs(tx, ty, tz, ns); + } + im.uploadState(); +} + +// Write ns and, for solid cells, seed the populations at rest (so the solid +// stays quiescent under bounce-back) โ€” mirrors WindTunnel3D::seedRest. +void WindTunnelOpenCL::setNs(std::int64_t x, std::int64_t y, std::int64_t z, double ns) { + auto& im = *impl_; + const std::int64_t c = im.idx(x, y, z); + im.hns[static_cast(c)] = static_cast(ns); + if (ns > 0.5) + for (int i = 0; i < 19; ++i) + im.hf[static_cast(i * im.ncells + c)] = im.feq(i, 0.f); +} + +void WindTunnelOpenCL::run(std::int64_t steps) { + auto& im = *impl_; + const std::size_t global = static_cast(im.ncells); + const cl_long n = static_cast(im.ncells); + const cl_int nx = static_cast(im.nx), ny = static_cast(im.ny), + nz = static_cast(im.nz); + // Static scalar args. + clCheck(clSetKernelArg(im.collideK, 1, sizeof(cl_mem), &im.nsBuf), "arg collide ns"); + clCheck(clSetKernelArg(im.collideK, 2, sizeof(cl_long), &n), "arg collide n"); + clCheck(clSetKernelArg(im.collideK, 3, sizeof(cl_float), &im.omega), "arg collide omega"); + clCheck(clSetKernelArg(im.streamK, 2, sizeof(cl_mem), &im.nsBuf), "arg stream ns"); + clCheck(clSetKernelArg(im.streamK, 3, sizeof(cl_float), &im.uIn), "arg stream uIn"); + clCheck(clSetKernelArg(im.streamK, 4, sizeof(cl_int), &nx), "arg stream nx"); + clCheck(clSetKernelArg(im.streamK, 5, sizeof(cl_int), &ny), "arg stream ny"); + clCheck(clSetKernelArg(im.streamK, 6, sizeof(cl_int), &nz), "arg stream nz"); + clCheck(clSetKernelArg(im.streamK, 7, sizeof(cl_long), &n), "arg stream n"); + + for (std::int64_t s = 0; s < steps; ++s) { + clCheck(clSetKernelArg(im.collideK, 0, sizeof(cl_mem), &im.cur), "arg collide f"); + clCheck(clEnqueueNDRangeKernel(im.queue, im.collideK, 1, nullptr, &global, nullptr, 0, + nullptr, nullptr), + "enqueue collide"); + clCheck(clSetKernelArg(im.streamK, 0, sizeof(cl_mem), &im.cur), "arg stream src"); + clCheck(clSetKernelArg(im.streamK, 1, sizeof(cl_mem), &im.other), "arg stream dst"); + clCheck(clEnqueueNDRangeKernel(im.queue, im.streamK, 1, nullptr, &global, nullptr, 0, + nullptr, nullptr), + "enqueue stream"); + std::swap(im.cur, im.other); + } + clFinish(im.queue); + im.cacheValid = false; +} + +std::int64_t WindTunnelOpenCL::nx() const { return impl_->nx; } +std::int64_t WindTunnelOpenCL::ny() const { return impl_->ny; } +std::int64_t WindTunnelOpenCL::nz() const { return impl_->nz; } + +std::array WindTunnelOpenCL::velocity(std::int64_t x, std::int64_t y, + std::int64_t z) const { + const auto& im = *impl_; + im.ensureDownloaded(); + const std::int64_t c = im.idx(x, y, z); + double rho = 0, jx = 0, jy = 0, jz = 0; // fp64 accumulation for the read-back + for (int i = 0; i < 19; ++i) { + const double v = im.cache[static_cast(i * im.ncells + c)]; + rho += v; + jx += v * kCx[i]; + jy += v * kCy[i]; + jz += v * kCz[i]; + } + const double invRho = 1.0 / rho; + return {jx * invRho, jy * invRho, jz * invRho}; +} + +double WindTunnelOpenCL::solidFraction(std::int64_t x, std::int64_t y, std::int64_t z) const { + return impl_->hns[static_cast(impl_->idx(x, y, z))]; +} + +std::vector WindTunnelOpenCL::downloadPopulations() const { + impl_->ensureDownloaded(); + return impl_->cache; +} + +void WindTunnelOpenCL::writeVtk(const std::string& path) const { + io::VtkStructuredWriter<3> w(impl_->nx, impl_->ny, impl_->nz); + using Coord = io::VtkStructuredWriter<3>::Coord; + w.addVector("velocity", [this](Coord c) { return velocity(c[0], c[1], c[2]); }); + w.addScalar("speed", [this](Coord c) { + const auto u = velocity(c[0], c[1], c[2]); + return std::sqrt(u[0] * u[0] + u[1] * u[1] + u[2] * u[2]); + }); + w.addScalar("solid", [this](Coord c) { return solidFraction(c[0], c[1], c[2]); }); + w.write(path); +} + +bool WindTunnelOpenCL::openclAvailable() { + cl_platform_id plat; + cl_device_id dev; + return pickDevice(plat, dev); +} + +} // namespace cyberfluids::solver + +#endif // CYBERFLUIDS_WITH_OPENCL diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 207a2cf..28662cf 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -33,6 +33,22 @@ if(CYBERFLUIDS_HAVE_NUMPP) add_test(NAME cmg_roundtrip COMMAND test_cmg_roundtrip) endif() + # OpenCL GPU wind-tunnel test โ€” built only when the OpenCL backend is present. + if(CYBERFLUIDS_HAVE_OPENCL) + add_executable(test_opencl_wind_tunnel test_opencl_wind_tunnel.cpp) + target_link_libraries(test_opencl_wind_tunnel PRIVATE cyberfluids_core numpp::numpp cyberfluids_opencl) + target_include_directories(test_opencl_wind_tunnel PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) + add_test(NAME opencl_wind_tunnel COMMAND test_opencl_wind_tunnel) + endif() + + # CUDA GPU wind-tunnel test โ€” built only when the CUDA backend is present. + if(CYBERFLUIDS_HAVE_CUDA) + add_executable(test_cuda_wind_tunnel test_cuda_wind_tunnel.cpp) + target_link_libraries(test_cuda_wind_tunnel PRIVATE cyberfluids_core numpp::numpp cyberfluids_cuda) + target_include_directories(test_cuda_wind_tunnel PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) + add_test(NAME cuda_wind_tunnel COMMAND test_cuda_wind_tunnel) + endif() + # Metal GPU cavity test โ€” built only when the Metal backend is present. if(CYBERFLUIDS_HAVE_METAL) add_executable(test_metal_cavity test_metal_cavity.cpp) diff --git a/tests/test_cuda_wind_tunnel.cpp b/tests/test_cuda_wind_tunnel.cpp new file mode 100644 index 0000000..dfa4b6f --- /dev/null +++ b/tests/test_cuda_wind_tunnel.cpp @@ -0,0 +1,71 @@ +/// CUDA analog of test_opencl_wind_tunnel: validates the CUDA wind tunnel (fp32) +/// against the CPU solver::WindTunnel3D (fp64) oracle. Skips cleanly (exit 0) when +/// no CUDA device is present. Runs on an NVIDIA machine with -DCYBERFLUIDS_CUDA=ON. + +#include +#include +#include +#include + +#include "cyberfluids/solver/cuda/wind_tunnel_cuda.hpp" +#include "cyberfluids/solver/wind_tunnel3d.hpp" +#include "testing.hpp" + +using cyberfluids::solver::WindTunnel3D; +using cyberfluids::solver::WindTunnelCuda; + +int main() { + if (!WindTunnelCuda::cudaAvailable()) { + std::printf("SKIP cuda_wind_tunnel: no CUDA device\n"); + return 0; + } + + const std::int64_t nx = 48, ny = 24, nz = 24; + const double Uin = 0.05, R = 4.0, cx = 14, cy = 12, cz = 12; + const double omega = WindTunnel3D<>::omegaForReynolds(Uin, 2 * R, 20.0); + const std::int64_t steps = 5000; + + WindTunnel3D<> cpu(nx, ny, nz, omega, Uin); + cpu.setObstacleSphere(cx, cy, cz, R, true); + cpu.run(steps); + + WindTunnelCuda gpu(nx, ny, nz, omega, Uin); + gpu.setObstacleSphere(cx, cy, cz, R, true); + gpu.run(steps); + + double linf = 0.0, sumDiff2 = 0.0, sumRef2 = 0.0, sumAbs = 0.0, solidMax = 0.0; + std::int64_t fluidN = 0; + for (std::int64_t x = 0; x < nx; ++x) + for (std::int64_t y = 0; y < ny; ++y) + for (std::int64_t z = 0; z < nz; ++z) { + const auto ug = gpu.velocity(x, y, z); + if (cpu.solidFraction(x, y, z) > 0.5) { + solidMax = std::max(solidMax, std::sqrt(ug[0] * ug[0] + ug[1] * ug[1] + + ug[2] * ug[2])); + continue; + } + const auto uc = cpu.velocity(x, y, z); + double d2 = 0, r2 = 0, dmax = 0; + for (int k = 0; k < 3; ++k) { + const double d = ug[k] - uc[k]; + d2 += d * d; + r2 += uc[k] * uc[k]; + dmax = std::max(dmax, std::fabs(d)); + } + linf = std::max(linf, dmax / Uin); + sumAbs += std::sqrt(d2) / Uin; + sumDiff2 += d2; + sumRef2 += r2; + ++fluidN; + } + std::printf("cuda vs cpu: Linf/Uin=%.4g meanRel=%.4g L2rel=%.4g\n", linf, + sumAbs / static_cast(fluidN), std::sqrt(sumDiff2 / sumRef2)); + + CF_CHECK(linf < 0.05); + CF_CHECK(std::sqrt(sumDiff2 / sumRef2) < 0.02); + CF_CHECK(sumAbs / static_cast(fluidN) < 0.01); + CF_CHECK(solidMax < 0.1 * Uin); + + if (cftest::failures == 0) std::printf("cuda_wind_tunnel: all checks passed\n"); + return cftest::failures; +} diff --git a/tests/test_opencl_wind_tunnel.cpp b/tests/test_opencl_wind_tunnel.cpp new file mode 100644 index 0000000..f92ff29 --- /dev/null +++ b/tests/test_opencl_wind_tunnel.cpp @@ -0,0 +1,96 @@ +/// Validates the OpenCL wind tunnel (fp32, GPU) against the CPU solver::WindTunnel3D +/// (fp64) oracle on the same case, plus determinism and no-slip. Skips cleanly +/// (exit 0) when no OpenCL device is available. Runs on this machine's GPU. + +#include +#include +#include +#include + +#include "cyberfluids/solver/opencl/wind_tunnel_opencl.hpp" +#include "cyberfluids/solver/wind_tunnel3d.hpp" +#include "testing.hpp" + +using cyberfluids::solver::WindTunnel3D; +using cyberfluids::solver::WindTunnelOpenCL; + +int main() { + if (!WindTunnelOpenCL::openclAvailable()) { + std::printf("SKIP opencl_wind_tunnel: no OpenCL device\n"); + return 0; + } + + const std::int64_t nx = 48, ny = 24, nz = 24; + const double Uin = 0.05, R = 4.0; + const double cx = 14, cy = 12, cz = 12; + const double omega = WindTunnel3D<>::omegaForReynolds(Uin, 2 * R, 20.0); + const std::int64_t steps = 5000; + + // fp64 CPU oracle. + WindTunnel3D<> cpu(nx, ny, nz, omega, Uin); + cpu.setObstacleSphere(cx, cy, cz, R, true); + cpu.run(steps); + + // fp32 GPU. + WindTunnelOpenCL gpu(nx, ny, nz, omega, Uin); + gpu.setObstacleSphere(cx, cy, cz, R, true); + gpu.run(steps); + + // (1) Obstacle fidelity: ns computed with identical host math -> exact match. + bool nsExact = true; + for (std::int64_t x = 0; x < nx; ++x) + for (std::int64_t y = 0; y < ny; ++y) + for (std::int64_t z = 0; z < nz; ++z) + if (gpu.solidFraction(x, y, z) != cpu.solidFraction(x, y, z)) nsExact = false; + CF_CHECK(nsExact); + + // (2) Determinism: a second identical GPU run is bit-for-bit equal. + WindTunnelOpenCL gpu2(nx, ny, nz, omega, Uin); + gpu2.setObstacleSphere(cx, cy, cz, R, true); + gpu2.run(steps); + const auto p1 = gpu.downloadPopulations(), p2 = gpu2.downloadPopulations(); + bool deterministic = (p1.size() == p2.size()); + for (std::size_t i = 0; i < p1.size() && deterministic; ++i) + if (p1[i] != p2[i]) deterministic = false; + CF_CHECK(deterministic); + + // (3) GPU (fp32) vs CPU (fp64) velocity field over fluid cells, normalized by Uin. + double linf = 0.0, sumDiff2 = 0.0, sumRef2 = 0.0, sumAbs = 0.0, solidMax = 0.0; + std::int64_t fluidN = 0; + for (std::int64_t x = 0; x < nx; ++x) + for (std::int64_t y = 0; y < ny; ++y) + for (std::int64_t z = 0; z < nz; ++z) { + const auto ug = gpu.velocity(x, y, z); + if (cpu.solidFraction(x, y, z) > 0.5) { + const double sp = std::sqrt(ug[0] * ug[0] + ug[1] * ug[1] + ug[2] * ug[2]); + solidMax = std::max(solidMax, sp); + continue; + } + const auto uc = cpu.velocity(x, y, z); + double d2 = 0, r2 = 0, dmax = 0; + for (int k = 0; k < 3; ++k) { + const double d = ug[k] - uc[k]; + d2 += d * d; + r2 += uc[k] * uc[k]; + dmax = std::max(dmax, std::fabs(d)); + } + linf = std::max(linf, dmax / Uin); + sumAbs += std::sqrt(d2) / Uin; + sumDiff2 += d2; + sumRef2 += r2; + ++fluidN; + } + const double meanRel = sumAbs / static_cast(fluidN); + const double l2Rel = std::sqrt(sumDiff2 / sumRef2); + + std::printf("opencl vs cpu: Linf/Uin=%.4g meanRel=%.4g L2rel=%.4g solidMax/Uin=%.4g\n", + linf, meanRel, l2Rel, solidMax / Uin); + + CF_CHECK(linf < 0.05); // fp32-vs-fp64, matches the Metal-cavity bar + CF_CHECK(l2Rel < 0.02); // field-relative L2 agreement < 2% + CF_CHECK(meanRel < 0.01); // mean per-cell error < 1% of Uin + CF_CHECK(solidMax < 0.1 * Uin); // no-slip preserved inside the solid + + if (cftest::failures == 0) std::printf("opencl_wind_tunnel: all checks passed\n"); + return cftest::failures; +}