Add OpenCL + CUDA GPU wind-tunnel solvers

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)
 

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**

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) }
+}

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()

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()

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<br/>(written once)"] --> SEAM["Backend seam<br/>kernel dispatch"]
-    SEAM --> CPU["CPU"]
-    SEAM --> CUDA["CUDA"]
-    SEAM --> METAL["Metal"]
-    SEAM --> OCL["OpenCL / SYCL"]
+    MODEL["Model code"] --> SEAM["forEachIndex seam<br/>(CPU only)"]
+    SEAM --> CPU["CPU · par_unseq / serial"]
+    GPUK["Device kernels<br/>(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.

include/cyberfluids/solver/cuda/wind_tunnel_cuda.hpp

@@ -0,0 +1,53 @@
+#pragma once
+
+#if defined(CYBERFLUIDS_WITH_CUDA)
+
+#include <array>
+#include <cstdint>
+#include <memory>
+#include <string>
+#include <vector>
+
+#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<double, 3> 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<float> 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> impl_;
+};
+
+}  // namespace cyberfluids::solver
+
+#endif  // CYBERFLUIDS_WITH_CUDA

include/cyberfluids/solver/opencl/wind_tunnel_opencl.hpp

@@ -0,0 +1,60 @@
+#pragma once
+
+#if defined(CYBERFLUIDS_WITH_OPENCL)
+
+#include <array>
+#include <cstdint>
+#include <memory>
+#include <string>
+#include <vector>
+
+#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<double, 3> 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<float> 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> impl_;
+};
+
+}  // namespace cyberfluids::solver
+
+#endif  // CYBERFLUIDS_WITH_OPENCL

openspec/changes/add-gpu-wind-tunnel/.openspec.yaml

@@ -0,0 +1,2 @@
+schema: spec-driven
+created: 2026-07-05

openspec/changes/add-gpu-wind-tunnel/README.md

@@ -0,0 +1,3 @@
+# add-gpu-wind-tunnel
+
+OpenCL + CUDA device-kernel wind-tunnel solvers