From 338aa9afe2b947e9280a2fbfb255690faec60587 Mon Sep 17 00:00:00 2001 From: aeronauty Date: Wed, 25 Mar 2026 18:23:44 +0100 Subject: [PATCH 01/12] feat: add WebGPU CFD solver (Phase 1 - Euler with explicit time-stepping) Port g2d CUDA CFD solver concepts to WebGPU compute shaders for browser-based structured-grid compressible flow analysis. This adds a new "CFD mode" alongside the existing panel method solver. New crate: rustfoil-cfd - O-type structured mesh generation with geometric radial stretching - Freestream initial conditions (Q = [rho, rho*u, rho*v, E, nu_tilde]) - Boundary condition type mapping (wall, farfield, interior) - Solver configuration structs matching GPU uniform layout WGSL compute shaders (8 shaders): - metrics: grid Jacobian and metric terms - primitives: conservative-to-primitive variable conversion - reconstruct_muscl: 2nd-order MUSCL with minmod limiter - roe_flux: Roe approximate Riemann solver with entropy fix - residual: flux divergence RHS computation - update: explicit forward Euler with positivity enforcement - bc: slip/no-slip wall + characteristic farfield BCs - forces_reduce: parallel reduction for Cl, Cd, Cm + L2 norm TypeScript orchestration (CfdSolver): - GPU buffer management for ~31MB structured grid data - 8-pass per-timestep compute dispatch pipeline - Async GPU-to-CPU readback for convergence monitoring UI integration: - CfdPanel with grid/flow/physics controls and convergence display - Zustand cfdStore for solver state management - Extended SolverMode type with 'cfd' option Co-Authored-By: Claude Opus 4.6 (1M context) --- Cargo.lock | 9 + Cargo.toml | 1 + crates/rustfoil-cfd/Cargo.toml | 11 + crates/rustfoil-cfd/src/boundary.rs | 68 ++++ crates/rustfoil-cfd/src/config.rs | 132 ++++++++ crates/rustfoil-cfd/src/init.rs | 88 +++++ crates/rustfoil-cfd/src/lib.rs | 15 + crates/rustfoil-cfd/src/mesh.rs | 248 ++++++++++++++ crates/rustfoil-wasm/Cargo.toml | 1 + crates/rustfoil-wasm/src/lib.rs | 106 ++++++ .../src/components/panels/CfdPanel.tsx | 313 ++++++++++++++++++ flexfoil-ui/src/lib/webgpu/cfd/CfdBuffers.ts | 95 ++++++ flexfoil-ui/src/lib/webgpu/cfd/CfdConfig.ts | 83 +++++ .../src/lib/webgpu/cfd/CfdPipelines.ts | 214 ++++++++++++ flexfoil-ui/src/lib/webgpu/cfd/CfdSolver.ts | 258 +++++++++++++++ flexfoil-ui/src/lib/webgpu/cfd/index.ts | 11 + .../src/lib/webgpu/shaders/cfd/bc.wgsl | 125 +++++++ .../lib/webgpu/shaders/cfd/cfd_common.wgsl | 42 +++ .../lib/webgpu/shaders/cfd/forces_reduce.wgsl | 115 +++++++ .../src/lib/webgpu/shaders/cfd/metrics.wgsl | 58 ++++ .../lib/webgpu/shaders/cfd/primitives.wgsl | 37 +++ .../webgpu/shaders/cfd/reconstruct_muscl.wgsl | 109 ++++++ .../src/lib/webgpu/shaders/cfd/residual.wgsl | 51 +++ .../src/lib/webgpu/shaders/cfd/roe_flux.wgsl | 185 +++++++++++ .../src/lib/webgpu/shaders/cfd/update.wgsl | 47 +++ flexfoil-ui/src/stores/cfdStore.ts | 136 ++++++++ flexfoil-ui/src/types/index.ts | 2 +- 27 files changed, 2559 insertions(+), 1 deletion(-) create mode 100644 crates/rustfoil-cfd/Cargo.toml create mode 100644 crates/rustfoil-cfd/src/boundary.rs create mode 100644 crates/rustfoil-cfd/src/config.rs create mode 100644 crates/rustfoil-cfd/src/init.rs create mode 100644 crates/rustfoil-cfd/src/lib.rs create mode 100644 crates/rustfoil-cfd/src/mesh.rs create mode 100644 flexfoil-ui/src/components/panels/CfdPanel.tsx create mode 100644 flexfoil-ui/src/lib/webgpu/cfd/CfdBuffers.ts create mode 100644 flexfoil-ui/src/lib/webgpu/cfd/CfdConfig.ts create mode 100644 flexfoil-ui/src/lib/webgpu/cfd/CfdPipelines.ts create mode 100644 flexfoil-ui/src/lib/webgpu/cfd/CfdSolver.ts create mode 100644 flexfoil-ui/src/lib/webgpu/cfd/index.ts create mode 100644 flexfoil-ui/src/lib/webgpu/shaders/cfd/bc.wgsl create mode 100644 flexfoil-ui/src/lib/webgpu/shaders/cfd/cfd_common.wgsl create mode 100644 flexfoil-ui/src/lib/webgpu/shaders/cfd/forces_reduce.wgsl create mode 100644 flexfoil-ui/src/lib/webgpu/shaders/cfd/metrics.wgsl create mode 100644 flexfoil-ui/src/lib/webgpu/shaders/cfd/primitives.wgsl create mode 100644 flexfoil-ui/src/lib/webgpu/shaders/cfd/reconstruct_muscl.wgsl create mode 100644 flexfoil-ui/src/lib/webgpu/shaders/cfd/residual.wgsl create mode 100644 flexfoil-ui/src/lib/webgpu/shaders/cfd/roe_flux.wgsl create mode 100644 flexfoil-ui/src/lib/webgpu/shaders/cfd/update.wgsl create mode 100644 flexfoil-ui/src/stores/cfdStore.ts diff --git a/Cargo.lock b/Cargo.lock index da6e7d2d..e81327b6 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -605,6 +605,14 @@ dependencies = [ "serde_json", ] +[[package]] +name = "rustfoil-cfd" +version = "0.1.0" +dependencies = [ + "rustfoil-core", + "serde", +] + [[package]] name = "rustfoil-cli" version = "0.1.0" @@ -701,6 +709,7 @@ version = "0.1.0" dependencies = [ "console_error_panic_hook", "js-sys", + "rustfoil-cfd", "rustfoil-core", "rustfoil-inviscid", "rustfoil-solver", diff --git a/Cargo.toml b/Cargo.toml index 0524c4b8..37dde92f 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -12,6 +12,7 @@ members = [ "crates/rustfoil-cli", "crates/rustfoil-inviscid", "crates/rustfoil-xfoil", "crates/rustfoil-python", + "crates/rustfoil-cfd", ] [workspace.package] diff --git a/crates/rustfoil-cfd/Cargo.toml b/crates/rustfoil-cfd/Cargo.toml new file mode 100644 index 00000000..e30d8aa9 --- /dev/null +++ b/crates/rustfoil-cfd/Cargo.toml @@ -0,0 +1,11 @@ +[package] +name = "rustfoil-cfd" +version.workspace = true +edition.workspace = true +license.workspace = true +authors.workspace = true +description = "Structured mesh CFD solver support for RustFoil (mesh gen, initial conditions, BCs)" + +[dependencies] +rustfoil-core = { path = "../rustfoil-core" } +serde = { version = "1.0", features = ["derive"] } diff --git a/crates/rustfoil-cfd/src/boundary.rs b/crates/rustfoil-cfd/src/boundary.rs new file mode 100644 index 00000000..112456c1 --- /dev/null +++ b/crates/rustfoil-cfd/src/boundary.rs @@ -0,0 +1,68 @@ +//! Boundary condition type mapping for the structured O-grid. +//! +//! In an O-grid topology: +//! - j=0: airfoil surface (wall) +//! - j=nj-1: far-field boundary +//! - i-direction: periodic (wraps around the airfoil) +//! - Wake cut: handled via periodicity in i + +/// Boundary condition types (matches WGSL constants). +#[repr(u32)] +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum BcType { + /// Interior cell (no special treatment) + Interior = 0, + /// Solid wall (slip for Euler, no-slip for NS/RANS) + Wall = 1, + /// Far-field (characteristic-based) + FarField = 2, + /// Wake cut (periodic connection) + WakeCut = 3, +} + +/// Generate boundary condition type array for the entire grid. +/// +/// Returns a flat u32 array of length ni*nj. +pub fn generate_bc_types(ni: u32, nj: u32) -> Vec { + let ni = ni as usize; + let nj = nj as usize; + let mut bc = vec![BcType::Interior as u32; ni * nj]; + + // j=0: wall boundary + for i in 0..ni { + bc[i] = BcType::Wall as u32; + } + + // j=nj-1: far-field boundary + for i in 0..ni { + bc[(nj - 1) * ni + i] = BcType::FarField as u32; + } + + bc +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_bc_types() { + let bc = generate_bc_types(8, 4); + assert_eq!(bc.len(), 32); + + // Wall at j=0 + for i in 0..8 { + assert_eq!(bc[i], BcType::Wall as u32); + } + // Interior at j=1,2 + for j in 1..3 { + for i in 0..8 { + assert_eq!(bc[j * 8 + i], BcType::Interior as u32); + } + } + // FarField at j=3 + for i in 0..8 { + assert_eq!(bc[3 * 8 + i], BcType::FarField as u32); + } + } +} diff --git a/crates/rustfoil-cfd/src/config.rs b/crates/rustfoil-cfd/src/config.rs new file mode 100644 index 00000000..ab9c47f2 --- /dev/null +++ b/crates/rustfoil-cfd/src/config.rs @@ -0,0 +1,132 @@ +use serde::{Deserialize, Serialize}; + +/// Physics mode for the CFD solver. +#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)] +pub enum PhysicsMode { + /// Inviscid compressible (Euler equations) + Euler = 0, + /// Laminar Navier-Stokes + LaminarNS = 1, + /// Reynolds-Averaged Navier-Stokes with Spalart-Allmaras + RansSA = 2, +} + +/// Spatial reconstruction scheme. +#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)] +pub enum ReconstructionMode { + /// 2nd-order MUSCL with minmod limiter + Muscl = 0, + /// 5th-order WENO5 + Weno5 = 1, +} + +/// Time-stepping scheme. +#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)] +pub enum TimeSteppingMode { + /// Explicit forward Euler (Phase 1) + ExplicitEuler = 0, + /// Diagonalized ADI (Phase 2) + Dadi = 1, +} + +/// CFD solver configuration. +/// +/// This struct is shared between Rust/WASM and TypeScript (via serde). +/// Fields are laid out to match the GPU uniform buffer `CfdParams`. +#[derive(Debug, Clone, Serialize, Deserialize)] +pub struct CfdConfig { + /// Circumferential grid points (around airfoil) + pub ni: u32, + /// Radial grid points (away from airfoil) + pub nj: u32, + /// Ratio of specific heats (default 1.4 for air) + pub gamma: f32, + /// CFL number for time stepping + pub cfl: f32, + /// Freestream Mach number + pub mach_inf: f32, + /// Angle of attack in radians + pub alpha: f32, + /// Reynolds number (used for NS/RANS) + pub reynolds: f32, + /// Prandtl number (default 0.72 for air) + pub prandtl: f32, + /// Physics mode + pub physics: PhysicsMode, + /// Reconstruction scheme + pub reconstruction: ReconstructionMode, + /// Time-stepping scheme + pub time_stepping: TimeSteppingMode, + /// Far-field distance in chord lengths + pub far_field: f32, + /// First cell wall-normal spacing (for viscous grids) + pub ds0: f32, +} + +impl Default for CfdConfig { + fn default() -> Self { + Self { + ni: 256, + nj: 128, + gamma: 1.4, + cfl: 0.5, + mach_inf: 0.5, + alpha: 0.0, + reynolds: 1e6, + prandtl: 0.72, + physics: PhysicsMode::Euler, + reconstruction: ReconstructionMode::Muscl, + time_stepping: TimeSteppingMode::ExplicitEuler, + far_field: 20.0, + ds0: 1e-4, + } + } +} + +/// GPU uniform buffer layout for CfdParams. +/// Must match the WGSL struct exactly (std140/std430 layout). +#[repr(C)] +#[derive(Debug, Clone, Copy)] +pub struct CfdParamsGpu { + pub ni: u32, + pub nj: u32, + pub gamma: f32, + pub cfl: f32, + pub mach_inf: f32, + pub alpha: f32, + pub reynolds: f32, + pub prandtl: f32, + pub dt: f32, + pub iteration: u32, + pub physics_mode: u32, + pub reconstruction: u32, +} + +impl CfdParamsGpu { + pub fn from_config(config: &CfdConfig, dt: f32, iteration: u32) -> Self { + Self { + ni: config.ni, + nj: config.nj, + gamma: config.gamma, + cfl: config.cfl, + mach_inf: config.mach_inf, + alpha: config.alpha, + reynolds: config.reynolds, + prandtl: config.prandtl, + dt, + iteration, + physics_mode: config.physics as u32, + reconstruction: config.reconstruction as u32, + } + } + + /// Return as raw bytes for GPU buffer upload. + pub fn as_bytes(&self) -> &[u8] { + unsafe { + std::slice::from_raw_parts( + self as *const Self as *const u8, + std::mem::size_of::(), + ) + } + } +} diff --git a/crates/rustfoil-cfd/src/init.rs b/crates/rustfoil-cfd/src/init.rs new file mode 100644 index 00000000..b79b391a --- /dev/null +++ b/crates/rustfoil-cfd/src/init.rs @@ -0,0 +1,88 @@ +//! Initial condition computation for the CFD solver. +//! +//! Generates freestream uniform initial conditions for the conservative +//! variables Q = [rho, rho*u, rho*v, E, nu_tilde]. + +use crate::config::CfdConfig; + +/// Compute freestream initial conditions for the entire grid. +/// +/// Returns a flat f32 array of length ni*nj*5 containing: +/// Q = [rho, rho*u, rho*v, E, nu_tilde] at each cell. +/// +/// Uses non-dimensionalization where: +/// - rho_inf = 1.0 +/// - |V_inf| = M_inf * a_inf = M_inf (with a_inf = 1.0) +/// - p_inf = 1.0 / gamma +pub fn compute_initial_conditions(config: &CfdConfig) -> Vec { + let ni = config.ni as usize; + let nj = config.nj as usize; + let gamma = config.gamma; + let mach = config.mach_inf; + let alpha = config.alpha; // radians + + // Freestream state (non-dimensional) + let rho_inf: f32 = 1.0; + let u_inf: f32 = mach * alpha.cos(); + let v_inf: f32 = mach * alpha.sin(); + let p_inf: f32 = 1.0 / gamma; + let e_inf: f32 = p_inf / (gamma - 1.0) + 0.5 * rho_inf * (u_inf * u_inf + v_inf * v_inf); + + // SA turbulent viscosity initial value + let nu_tilde_inf: f32 = if config.physics as u32 == 2 { + // RANS: initialize to ~3x molecular viscosity + 3.0 * mach / config.reynolds + } else { + 0.0 + }; + + let n_total = ni * nj * 5; + let mut q = vec![0.0f32; n_total]; + + for j in 0..nj { + for i in 0..ni { + let base = (j * ni + i) * 5; + q[base] = rho_inf; + q[base + 1] = rho_inf * u_inf; + q[base + 2] = rho_inf * v_inf; + q[base + 3] = e_inf; + q[base + 4] = nu_tilde_inf; + } + } + + q +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_initial_conditions_euler() { + let config = CfdConfig { + ni: 4, + nj: 4, + mach_inf: 0.5, + alpha: 0.0, + gamma: 1.4, + ..CfdConfig::default() + }; + + let q = compute_initial_conditions(&config); + assert_eq!(q.len(), 4 * 4 * 5); + + // Check first cell + let rho = q[0]; + let rhou = q[1]; + let rhov = q[2]; + let e = q[3]; + let nu = q[4]; + + assert!((rho - 1.0).abs() < 1e-6); + assert!((rhou - 0.5).abs() < 1e-6); // rho * M * cos(0) + assert!(rhov.abs() < 1e-6); // rho * M * sin(0) + assert!(nu.abs() < 1e-6); // Euler: no SA + // E = p/(gamma-1) + 0.5*rho*V^2 = (1/1.4)/0.4 + 0.5*0.25 = 1.7857 + 0.125 = 1.9107 + assert!((e - 1.9107).abs() < 0.01); + } +} diff --git a/crates/rustfoil-cfd/src/lib.rs b/crates/rustfoil-cfd/src/lib.rs new file mode 100644 index 00000000..2cd2b1df --- /dev/null +++ b/crates/rustfoil-cfd/src/lib.rs @@ -0,0 +1,15 @@ +//! RustFoil CFD - Structured mesh generation and setup for GPU-accelerated CFD. +//! +//! This crate provides CPU-side support for the WebGPU-based 2D CFD solver: +//! - O-type structured mesh generation around airfoils +//! - Solver configuration and parameter management +//! - Initial condition computation (freestream state) +//! - Boundary condition type mapping +//! +//! All heavy computation (flux evaluation, time stepping, turbulence model) +//! runs on WebGPU compute shaders. This crate prepares the data for GPU upload. + +pub mod boundary; +pub mod config; +pub mod init; +pub mod mesh; diff --git a/crates/rustfoil-cfd/src/mesh.rs b/crates/rustfoil-cfd/src/mesh.rs new file mode 100644 index 00000000..d86cb0bb --- /dev/null +++ b/crates/rustfoil-cfd/src/mesh.rs @@ -0,0 +1,248 @@ +//! Structured O-grid mesh generation for CFD. +//! +//! Generates a body-fitted O-topology mesh around an airfoil by: +//! 1. Distributing points along the airfoil surface (circumferential direction) +//! 2. Algebraically growing layers outward (radial direction) +//! 3. Optionally smoothing with Poisson iteration + +/// Result of mesh generation: flat arrays ready for GPU upload. +#[derive(Debug, Clone)] +pub struct MeshResult { + /// Node x-coordinates, length ni*nj (row-major: j*ni + i) + pub x: Vec, + /// Node y-coordinates, length ni*nj + pub y: Vec, + /// Grid dimensions + pub ni: u32, + pub nj: u32, +} + +/// Generate an O-type structured mesh around an airfoil. +/// +/// # Arguments +/// * `airfoil_x` - Airfoil x-coordinates (CCW from TE lower, around LE, to TE upper) +/// * `airfoil_y` - Airfoil y-coordinates +/// * `ni` - Number of circumferential points (should match airfoil point count or resample) +/// * `nj` - Number of radial layers +/// * `far_field` - Far-field radius in chord lengths +/// * `ds0` - First cell wall spacing (normalized by chord) +/// +/// # Returns +/// `MeshResult` with flat coordinate arrays +pub fn generate_o_mesh( + airfoil_x: &[f64], + airfoil_y: &[f64], + ni: u32, + nj: u32, + far_field: f64, + ds0: f64, +) -> MeshResult { + let ni = ni as usize; + let nj = nj as usize; + // Resample airfoil to ni points using linear interpolation along arc length + let (body_x, body_y) = resample_airfoil(airfoil_x, airfoil_y, ni); + + // Compute airfoil centroid for radial direction + let cx: f64 = body_x.iter().sum::() / ni as f64; + let cy: f64 = body_y.iter().sum::() / ni as f64; + + // Compute outward normal directions at each surface point + let mut nx = vec![0.0f64; ni]; + let mut ny = vec![0.0f64; ni]; + for i in 0..ni { + // Use centroid-to-point direction as approximate outward normal + let dx = body_x[i] - cx; + let dy = body_y[i] - cy; + let mag = (dx * dx + dy * dy).sqrt().max(1e-15); + nx[i] = dx / mag; + ny[i] = dy / mag; + } + + // Generate radial distribution with geometric stretching + let radial_dist = generate_radial_distribution(nj, ds0, far_field); + + // Build mesh: algebraic growth along normals + let mut x = vec![0.0f32; ni * nj]; + let mut y = vec![0.0f32; ni * nj]; + + for j in 0..nj { + let r = radial_dist[j]; + for i in 0..ni { + let idx = j * ni + i; + x[idx] = (body_x[i] + r * nx[i]) as f32; + y[idx] = (body_y[i] + r * ny[i]) as f32; + } + } + + // Laplacian smoothing passes for interior points + laplacian_smooth(&mut x, &mut y, ni, nj, 50); + + MeshResult { + x, + y, + ni: ni as u32, + nj: nj as u32, + } +} + +/// Resample airfoil coordinates to exactly `n` points using arc-length parameterization. +fn resample_airfoil(ax: &[f64], ay: &[f64], n: usize) -> (Vec, Vec) { + let m = ax.len(); + if m == n { + return (ax.to_vec(), ay.to_vec()); + } + + // Compute cumulative arc length + let mut s = vec![0.0f64; m]; + for i in 1..m { + let dx = ax[i] - ax[i - 1]; + let dy = ay[i] - ay[i - 1]; + s[i] = s[i - 1] + (dx * dx + dy * dy).sqrt(); + } + let total_s = s[m - 1]; + + // Interpolate at uniform arc-length intervals + let mut rx = vec![0.0f64; n]; + let mut ry = vec![0.0f64; n]; + let mut j = 0usize; + for i in 0..n { + let target = total_s * (i as f64) / ((n - 1) as f64); + // Advance to the segment containing target + while j < m - 2 && s[j + 1] < target { + j += 1; + } + let seg_len = s[j + 1] - s[j]; + let t = if seg_len > 1e-15 { + (target - s[j]) / seg_len + } else { + 0.0 + }; + rx[i] = ax[j] + t * (ax[j + 1] - ax[j]); + ry[i] = ay[j] + t * (ay[j + 1] - ay[j]); + } + + (rx, ry) +} + +/// Generate radial distance distribution with geometric stretching. +/// +/// The first cell has height `ds0`, and cells grow geometrically until +/// the total radial extent reaches `far_field`. +fn generate_radial_distribution(nj: usize, ds0: f64, far_field: f64) -> Vec { + let mut dist = vec![0.0f64; nj]; + // j=0 is the body surface + dist[0] = 0.0; + + if nj < 2 { + return dist; + } + + // Find geometric stretch ratio using bisection + // Total distance: ds0 * (ratio^(nj-1) - 1) / (ratio - 1) = far_field + let ratio = find_stretch_ratio(nj - 1, ds0, far_field); + + let mut ds = ds0; + for j in 1..nj { + dist[j] = dist[j - 1] + ds; + ds *= ratio; + } + + dist +} + +/// Find geometric stretch ratio via bisection. +fn find_stretch_ratio(n: usize, ds0: f64, total: f64) -> f64 { + let mut lo = 1.0001f64; + let mut hi = 2.0f64; + + // Ensure hi is large enough + while geometric_sum(n, ds0, hi) < total { + hi *= 2.0; + } + + for _ in 0..100 { + let mid = 0.5 * (lo + hi); + let sum = geometric_sum(n, ds0, mid); + if sum < total { + lo = mid; + } else { + hi = mid; + } + } + + 0.5 * (lo + hi) +} + +fn geometric_sum(n: usize, ds0: f64, ratio: f64) -> f64 { + if (ratio - 1.0).abs() < 1e-10 { + return ds0 * n as f64; + } + ds0 * (ratio.powi(n as i32) - 1.0) / (ratio - 1.0) +} + +/// Simple Laplacian smoothing for interior mesh points. +fn laplacian_smooth(x: &mut [f32], y: &mut [f32], ni: usize, nj: usize, iters: usize) { + let omega = 0.3f32; // Under-relaxation factor + + for _ in 0..iters { + for j in 1..nj - 1 { + for i in 0..ni { + let idx = j * ni + i; + // Periodic in i-direction (O-grid wraps around) + let ip = if i + 1 < ni { i + 1 } else { 0 }; + let im = if i > 0 { i - 1 } else { ni - 1 }; + + let avg_x = 0.25 + * (x[j * ni + ip] + x[j * ni + im] + x[(j + 1) * ni + i] + x[(j - 1) * ni + i]); + let avg_y = 0.25 + * (y[j * ni + ip] + y[j * ni + im] + y[(j + 1) * ni + i] + y[(j - 1) * ni + i]); + + x[idx] += omega * (avg_x - x[idx]); + y[idx] += omega * (avg_y - y[idx]); + } + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + use std::f64::consts::PI; + + #[test] + fn test_radial_distribution() { + let dist = generate_radial_distribution(64, 0.001, 20.0); + assert_eq!(dist.len(), 64); + assert!((dist[0]).abs() < 1e-10); + assert!((dist[1] - 0.001).abs() < 1e-6); + assert!((dist[63] - 20.0).abs() / 20.0 < 0.01); + } + + #[test] + fn test_resample_identity() { + let x = vec![0.0, 0.5, 1.0]; + let y = vec![0.0, 0.1, 0.0]; + let (rx, ry) = resample_airfoil(&x, &y, 3); + for i in 0..3 { + assert!((rx[i] - x[i]).abs() < 1e-10); + assert!((ry[i] - y[i]).abs() < 1e-10); + } + } + + #[test] + fn test_generate_o_mesh() { + // Simple circle as "airfoil" for testing + let n = 64; + let mut ax = vec![0.0; n]; + let mut ay = vec![0.0; n]; + for i in 0..n { + let theta = 2.0 * PI * (i as f64) / (n as f64); + ax[i] = 0.5 + 0.5 * theta.cos(); + ay[i] = 0.5 * theta.sin(); + } + + let mesh = generate_o_mesh(&ax, &ay, 64, 32, 10.0, 0.01); + assert_eq!(mesh.x.len(), 64 * 32); + assert_eq!(mesh.y.len(), 64 * 32); + } +} diff --git a/crates/rustfoil-wasm/Cargo.toml b/crates/rustfoil-wasm/Cargo.toml index 09dd5058..d358e742 100644 --- a/crates/rustfoil-wasm/Cargo.toml +++ b/crates/rustfoil-wasm/Cargo.toml @@ -14,6 +14,7 @@ rustfoil-core = { path = "../rustfoil-core" } rustfoil-solver = { path = "../rustfoil-solver" } rustfoil-inviscid = { path = "../rustfoil-inviscid" } rustfoil-xfoil = { path = "../rustfoil-xfoil" } +rustfoil-cfd = { path = "../rustfoil-cfd" } wasm-bindgen = "0.2" serde = { version = "1.0", features = ["derive"] } serde-wasm-bindgen = "0.6" diff --git a/crates/rustfoil-wasm/src/lib.rs b/crates/rustfoil-wasm/src/lib.rs index d1dc97e7..37757106 100644 --- a/crates/rustfoil-wasm/src/lib.rs +++ b/crates/rustfoil-wasm/src/lib.rs @@ -3271,6 +3271,112 @@ impl Default for RustFoil { } } +// ============================================================================ +// CFD Mesh Generation & Setup +// ============================================================================ + +/// Generate a structured O-type mesh around an airfoil for the CFD solver. +/// +/// # Arguments +/// * `coords_flat` - Flat array of [x0, y0, x1, y1, ...] airfoil coordinates +/// * `ni` - Circumferential grid points +/// * `nj` - Radial grid points (layers away from airfoil) +/// * `far_field` - Far-field distance in chord lengths +/// * `ds0` - First cell wall-normal spacing +/// +/// # Returns +/// JsValue containing { x: Float32Array, y: Float32Array, ni: u32, nj: u32 } +#[wasm_bindgen] +pub fn cfd_generate_mesh( + coords_flat: &[f64], + ni: u32, + nj: u32, + far_field: f64, + ds0: f64, +) -> JsValue { + let n_pts = coords_flat.len() / 2; + let mut ax = Vec::with_capacity(n_pts); + let mut ay = Vec::with_capacity(n_pts); + for i in 0..n_pts { + ax.push(coords_flat[i * 2]); + ay.push(coords_flat[i * 2 + 1]); + } + + let mesh = rustfoil_cfd::mesh::generate_o_mesh(&ax, &ay, ni, nj, far_field, ds0); + + #[derive(Serialize)] + struct MeshOutput { + x: Vec, + y: Vec, + ni: u32, + nj: u32, + } + + let output = MeshOutput { + x: mesh.x, + y: mesh.y, + ni: mesh.ni, + nj: mesh.nj, + }; + + serde_wasm_bindgen::to_value(&output).unwrap_or(JsValue::NULL) +} + +/// Compute freestream initial conditions for the CFD solver. +/// +/// # Arguments +/// * `ni` - Grid circumferential points +/// * `nj` - Grid radial points +/// * `mach` - Freestream Mach number +/// * `alpha_deg` - Angle of attack in degrees +/// * `gamma` - Ratio of specific heats (1.4 for air) +/// * `physics_mode` - 0=Euler, 1=LaminarNS, 2=RANS_SA +/// * `reynolds` - Reynolds number (for RANS initialization) +/// +/// # Returns +/// Float32Array of length ni*nj*5 with Q = [rho, rho*u, rho*v, E, nu_tilde] +#[wasm_bindgen] +pub fn cfd_initial_conditions( + ni: u32, + nj: u32, + mach: f32, + alpha_deg: f32, + gamma: f32, + physics_mode: u32, + reynolds: f32, +) -> Vec { + use rustfoil_cfd::config::{CfdConfig, PhysicsMode}; + + let physics = match physics_mode { + 0 => PhysicsMode::Euler, + 1 => PhysicsMode::LaminarNS, + 2 => PhysicsMode::RansSA, + _ => PhysicsMode::Euler, + }; + + let config = CfdConfig { + ni, + nj, + gamma, + mach_inf: mach, + alpha: alpha_deg.to_radians(), + physics, + reynolds, + ..CfdConfig::default() + }; + + rustfoil_cfd::init::compute_initial_conditions(&config) +} + +/// Generate boundary condition type array for the CFD grid. +/// +/// # Returns +/// Uint32Array of length ni*nj with BC types (0=interior, 1=wall, 2=farfield) +#[wasm_bindgen] +pub fn cfd_boundary_types(ni: u32, nj: u32) -> Vec { + rustfoil_cfd::boundary::generate_bc_types(ni, nj) +} + #[cfg(test)] mod tests { use super::*; diff --git a/flexfoil-ui/src/components/panels/CfdPanel.tsx b/flexfoil-ui/src/components/panels/CfdPanel.tsx new file mode 100644 index 00000000..1c0ba62c --- /dev/null +++ b/flexfoil-ui/src/components/panels/CfdPanel.tsx @@ -0,0 +1,313 @@ +/** + * CfdPanel - Controls for the WebGPU CFD solver. + * + * Provides configuration for grid size, physics mode, Mach/alpha/Re, + * start/stop/reset controls, and convergence history plot. + */ + +import { useState, useCallback, useRef, useEffect } from 'react'; +import { useCfdStore } from '../../stores/cfdStore'; +import { useAirfoilStore } from '../../stores/airfoilStore'; +import { CfdSolver } from '../../lib/webgpu/cfd/CfdSolver'; +import type { CfdSolverResult } from '../../lib/webgpu/cfd/CfdSolver'; +import type { PhysicsMode } from '../../lib/webgpu/cfd/CfdConfig'; + +export function CfdPanel() { + const { + config, + isRunning, + iteration, + convergenceHistory, + forceHistory, + cl, + cd, + cm, + residualL2, + meshGenerated, + setNi, + setNj, + setMach, + setAlphaDeg, + setReynolds, + setCfl, + setPhysics, + setRunning, + setMeshGenerated, + appendResult, + reset, + } = useCfdStore(); + + const coordinates = useAirfoilStore((s) => s.coordinates); + + const solverRef = useRef(null); + const runningRef = useRef(false); + const [error, setError] = useState(null); + + // Sync running ref with state + useEffect(() => { + runningRef.current = isRunning; + }, [isRunning]); + + const handleGenerateMesh = useCallback(async () => { + if (!coordinates || coordinates.length < 6) { + setError('Load an airfoil first'); + return; + } + + try { + setError(null); + // Import WASM functions dynamically + const wasm = await import('rustfoil-wasm'); + const meshResult = wasm.cfd_generate_mesh( + coordinates, + config.ni, + config.nj, + config.farField, + config.ds0, + ); + + if (!meshResult || !meshResult.x) { + setError('Mesh generation failed'); + return; + } + + // Get initial conditions + const physicsMode = config.physics === 'euler' ? 0 : config.physics === 'laminar_ns' ? 1 : 2; + const initialQ = wasm.cfd_initial_conditions( + config.ni, + config.nj, + config.machInf, + config.alphaDeg, + config.gamma, + physicsMode, + config.reynolds, + ); + + const bcTypes = wasm.cfd_boundary_types(config.ni, config.nj); + + // Check for WebGPU + if (!navigator.gpu) { + setError('WebGPU not available in this browser'); + return; + } + + const adapter = await navigator.gpu.requestAdapter({ powerPreference: 'high-performance' }); + if (!adapter) { + setError('No WebGPU adapter found'); + return; + } + + const device = await adapter.requestDevice({ + requiredLimits: { + maxStorageBufferBindingSize: adapter.limits.maxStorageBufferBindingSize, + maxComputeWorkgroupsPerDimension: adapter.limits.maxComputeWorkgroupsPerDimension, + }, + }); + + // Destroy previous solver if any + solverRef.current?.destroy(); + + solverRef.current = CfdSolver.create( + device, + config, + new Float32Array(meshResult.x), + new Float32Array(meshResult.y), + new Float32Array(initialQ), + new Uint32Array(bcTypes), + ); + + setMeshGenerated(true); + reset(); + setMeshGenerated(true); // re-set after reset clears it + } catch (e) { + setError(`Mesh generation error: ${e instanceof Error ? e.message : String(e)}`); + } + }, [coordinates, config, setMeshGenerated, reset]); + + const handleStart = useCallback(async () => { + if (!solverRef.current) { + setError('Generate mesh first'); + return; + } + + setError(null); + setRunning(true); + runningRef.current = true; + + const BATCH_SIZE = 20; + + try { + while (runningRef.current) { + const result: CfdSolverResult = await solverRef.current.step(BATCH_SIZE); + appendResult(result); + + // Check convergence (residual dropped 6 orders of magnitude) + if (result.residualL2 < 1e-10) { + setRunning(false); + break; + } + + // Yield to UI + await new Promise((r) => setTimeout(r, 0)); + } + } catch (e) { + setError(`Solver error: ${e instanceof Error ? e.message : String(e)}`); + setRunning(false); + } + }, [appendResult, setRunning]); + + const handleStop = useCallback(() => { + setRunning(false); + runningRef.current = false; + }, [setRunning]); + + const handleReset = useCallback(() => { + handleStop(); + reset(); + solverRef.current?.destroy(); + solverRef.current = null; + setMeshGenerated(false); + }, [handleStop, reset, setMeshGenerated]); + + // Cleanup on unmount + useEffect(() => { + return () => { + solverRef.current?.destroy(); + }; + }, []); + + return ( +
+

CFD Solver (WebGPU)

+ + {error && ( +
+ {error} +
+ )} + + {/* Grid Settings */} +
+ Grid +
+ + +
+
+ + {/* Flow Conditions */} +
+ Flow Conditions +
+ + + + +
+
+ + {/* Physics Mode */} +
+ Physics + +
+ + {/* Controls */} +
+ + {!isRunning ? ( + + ) : ( + + )} + +
+ + {/* Results */} +
+ Results +
+
Iteration: {iteration}
+
Cl: {cl.toFixed(6)}
+
Cd: {cd.toFixed(6)}
+
Cm: {cm.toFixed(6)}
+
Residual L2: {residualL2.toExponential(3)}
+
+
+ + {/* Convergence mini-plot (text-based for now) */} + {convergenceHistory.length > 0 && ( +
+ Convergence +
+ {convergenceHistory.slice(-10).map((p) => ( +
+ [{p.iteration.toString().padStart(6)}] L2={p.residualL2.toExponential(3)} +
+ ))} +
+
+ )} +
+ ); +} + +const inputStyle: React.CSSProperties = { + width: '70px', + background: '#333', + color: '#eee', + border: '1px solid #555', + borderRadius: '3px', + padding: '2px 4px', + fontSize: '12px', +}; + +const btnStyle: React.CSSProperties = { + flex: 1, + padding: '4px 8px', + background: '#555', + color: '#eee', + border: '1px solid #666', + borderRadius: '3px', + cursor: 'pointer', + fontSize: '12px', +}; diff --git a/flexfoil-ui/src/lib/webgpu/cfd/CfdBuffers.ts b/flexfoil-ui/src/lib/webgpu/cfd/CfdBuffers.ts new file mode 100644 index 00000000..f37aa6b8 --- /dev/null +++ b/flexfoil-ui/src/lib/webgpu/cfd/CfdBuffers.ts @@ -0,0 +1,95 @@ +/** + * GPU buffer management for the CFD solver. + * Allocates all storage/uniform buffers needed for the compute pipeline. + */ + +import { NVAR } from './CfdConfig'; + +export interface CfdBufferSet { + /** Conservative variables Q = [rho, rho*u, rho*v, E, nu_tilde] */ + Q: GPUBuffer; + /** Primitive variables W = [rho, u, v, p, nu_tilde] */ + W: GPUBuffer; + /** Mesh x-coordinates */ + meshX: GPUBuffer; + /** Mesh y-coordinates */ + meshY: GPUBuffer; + /** Grid metrics [xi_x, xi_y, eta_x, eta_y, J] */ + metrics: GPUBuffer; + /** Numerical flux in xi-direction */ + fluxXi: GPUBuffer; + /** Numerical flux in eta-direction */ + fluxEta: GPUBuffer; + /** Reconstructed left states at xi-faces */ + leftXi: GPUBuffer; + /** Reconstructed right states at xi-faces */ + rightXi: GPUBuffer; + /** Reconstructed left states at eta-faces */ + leftEta: GPUBuffer; + /** Reconstructed right states at eta-faces */ + rightEta: GPUBuffer; + /** RHS residual */ + residual: GPUBuffer; + /** Boundary condition types (u32 per cell) */ + bcType: GPUBuffer; + /** Solver parameters uniform buffer */ + params: GPUBuffer; + /** Force accumulator [Cl, Cd, Cm] */ + forceAccum: GPUBuffer; + /** Residual L2 norm */ + residualNorm: GPUBuffer; + /** Readback buffer (MAP_READ) for async GPU->CPU data transfer */ + readback: GPUBuffer; +} + +export function createCfdBuffers(device: GPUDevice, ni: number, nj: number): CfdBufferSet { + const cellCount = ni * nj; + const qSize = cellCount * NVAR * 4; // f32 per variable + const coordSize = cellCount * 4; + const metricsSize = cellCount * 5 * 4; // 5 metrics per cell + + const STORAGE_RW = GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST | GPUBufferUsage.COPY_SRC; + const STORAGE_R = GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST; + + return { + Q: device.createBuffer({ label: 'CFD Q', size: qSize, usage: STORAGE_RW }), + W: device.createBuffer({ label: 'CFD W', size: qSize, usage: STORAGE_RW }), + meshX: device.createBuffer({ label: 'CFD mesh_x', size: coordSize, usage: STORAGE_R }), + meshY: device.createBuffer({ label: 'CFD mesh_y', size: coordSize, usage: STORAGE_R }), + metrics: device.createBuffer({ label: 'CFD metrics', size: metricsSize, usage: STORAGE_RW }), + fluxXi: device.createBuffer({ label: 'CFD flux_xi', size: qSize, usage: STORAGE_RW }), + fluxEta: device.createBuffer({ label: 'CFD flux_eta', size: qSize, usage: STORAGE_RW }), + leftXi: device.createBuffer({ label: 'CFD left_xi', size: qSize, usage: STORAGE_RW }), + rightXi: device.createBuffer({ label: 'CFD right_xi', size: qSize, usage: STORAGE_RW }), + leftEta: device.createBuffer({ label: 'CFD left_eta', size: qSize, usage: STORAGE_RW }), + rightEta: device.createBuffer({ label: 'CFD right_eta', size: qSize, usage: STORAGE_RW }), + residual: device.createBuffer({ label: 'CFD residual', size: qSize, usage: STORAGE_RW }), + bcType: device.createBuffer({ label: 'CFD bc_type', size: cellCount * 4, usage: STORAGE_R }), + params: device.createBuffer({ + label: 'CFD params', + size: 48, // CfdParams struct size + usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST, + }), + forceAccum: device.createBuffer({ + label: 'CFD force_accum', + size: 12, // 3 x f32 (Cl, Cd, Cm) + usage: STORAGE_RW | GPUBufferUsage.COPY_SRC, + }), + residualNorm: device.createBuffer({ + label: 'CFD residual_norm', + size: 16, // 4 x f32 + usage: STORAGE_RW | GPUBufferUsage.COPY_SRC, + }), + readback: device.createBuffer({ + label: 'CFD readback', + size: 28, // 7 x f32 (3 forces + 4 residual norms) + usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST, + }), + }; +} + +export function destroyCfdBuffers(buffers: CfdBufferSet): void { + for (const buf of Object.values(buffers)) { + (buf as GPUBuffer).destroy(); + } +} diff --git a/flexfoil-ui/src/lib/webgpu/cfd/CfdConfig.ts b/flexfoil-ui/src/lib/webgpu/cfd/CfdConfig.ts new file mode 100644 index 00000000..f693cdbe --- /dev/null +++ b/flexfoil-ui/src/lib/webgpu/cfd/CfdConfig.ts @@ -0,0 +1,83 @@ +/** + * CFD solver configuration types. + * Must match the Rust CfdConfig struct and WGSL CfdParams struct. + */ + +export type PhysicsMode = 'euler' | 'laminar_ns' | 'rans_sa'; +export type ReconstructionMode = 'muscl' | 'weno5'; +export type TimeSteppingMode = 'explicit_euler' | 'dadi'; + +export interface CfdConfig { + ni: number; + nj: number; + gamma: number; + cfl: number; + machInf: number; + alphaDeg: number; + reynolds: number; + prandtl: number; + physics: PhysicsMode; + reconstruction: ReconstructionMode; + timeStepping: TimeSteppingMode; + farField: number; + ds0: number; +} + +export const DEFAULT_CFD_CONFIG: CfdConfig = { + ni: 256, + nj: 128, + gamma: 1.4, + cfl: 0.5, + machInf: 0.5, + alphaDeg: 2.0, + reynolds: 1e6, + prandtl: 0.72, + physics: 'euler', + reconstruction: 'muscl', + timeStepping: 'explicit_euler', + farField: 20.0, + ds0: 1e-4, +}; + +const PHYSICS_MAP: Record = { + euler: 0, + laminar_ns: 1, + rans_sa: 2, +}; + +const RECON_MAP: Record = { + muscl: 0, + weno5: 1, +}; + +/** + * Pack CfdConfig into a Float32Array/Uint32Array for GPU uniform upload. + * Layout matches the WGSL CfdParams struct (48 bytes = 12 x u32/f32). + */ +export function packCfdParams( + config: CfdConfig, + dt: number, + iteration: number, +): ArrayBuffer { + const buffer = new ArrayBuffer(48); + const u32 = new Uint32Array(buffer); + const f32 = new Float32Array(buffer); + + u32[0] = config.ni; + u32[1] = config.nj; + f32[2] = config.gamma; + f32[3] = config.cfl; + f32[4] = config.machInf; + f32[5] = (config.alphaDeg * Math.PI) / 180.0; + f32[6] = config.reynolds; + f32[7] = config.prandtl; + f32[8] = dt; + u32[9] = iteration; + u32[10] = PHYSICS_MAP[config.physics]; + u32[11] = RECON_MAP[config.reconstruction]; + + return buffer; +} + +/** Number of conservative variables per cell */ +export const NVAR = 5; diff --git a/flexfoil-ui/src/lib/webgpu/cfd/CfdPipelines.ts b/flexfoil-ui/src/lib/webgpu/cfd/CfdPipelines.ts new file mode 100644 index 00000000..ad54642c --- /dev/null +++ b/flexfoil-ui/src/lib/webgpu/cfd/CfdPipelines.ts @@ -0,0 +1,214 @@ +/** + * Compute pipeline creation for the CFD solver. + * Each shader is loaded from WGSL source, with cfd_common.wgsl prepended. + */ + +import cfdCommonSrc from '../shaders/cfd/cfd_common.wgsl?raw'; +import metricsSrc from '../shaders/cfd/metrics.wgsl?raw'; +import primitivesSrc from '../shaders/cfd/primitives.wgsl?raw'; +import reconstructMusclSrc from '../shaders/cfd/reconstruct_muscl.wgsl?raw'; +import roeFluxSrc from '../shaders/cfd/roe_flux.wgsl?raw'; +import residualSrc from '../shaders/cfd/residual.wgsl?raw'; +import updateSrc from '../shaders/cfd/update.wgsl?raw'; +import bcSrc from '../shaders/cfd/bc.wgsl?raw'; +import forcesReduceSrc from '../shaders/cfd/forces_reduce.wgsl?raw'; + +import type { CfdBufferSet } from './CfdBuffers'; + +export interface CfdPipelineSet { + metrics: GPUComputePipeline; + primitives: GPUComputePipeline; + reconstructXi: GPUComputePipeline; + reconstructEta: GPUComputePipeline; + roeFluxXi: GPUComputePipeline; + roeFluxEta: GPUComputePipeline; + residual: GPUComputePipeline; + update: GPUComputePipeline; + bc: GPUComputePipeline; + forcesReduce: GPUComputePipeline; +} + +export interface CfdBindGroupSet { + metrics: GPUBindGroup; + primitives: GPUBindGroup; + reconstruct: GPUBindGroup; + roeFlux: GPUBindGroup; + residual: GPUBindGroup; + update: GPUBindGroup; + bc: GPUBindGroup; + forcesReduce: GPUBindGroup; +} + +function makeShader(device: GPUDevice, source: string, label: string): GPUShaderModule { + return device.createShaderModule({ + label, + code: cfdCommonSrc + '\n' + source, + }); +} + +export function createCfdPipelines(device: GPUDevice): CfdPipelineSet { + const auto = 'auto'; + + const metricsModule = makeShader(device, metricsSrc, 'cfd_metrics'); + const primitivesModule = makeShader(device, primitivesSrc, 'cfd_primitives'); + const reconstructModule = makeShader(device, reconstructMusclSrc, 'cfd_reconstruct_muscl'); + const roeFluxModule = makeShader(device, roeFluxSrc, 'cfd_roe_flux'); + const residualModule = makeShader(device, residualSrc, 'cfd_residual'); + const updateModule = makeShader(device, updateSrc, 'cfd_update'); + const bcModule = makeShader(device, bcSrc, 'cfd_bc'); + const forcesModule = makeShader(device, forcesReduceSrc, 'cfd_forces_reduce'); + + return { + metrics: device.createComputePipeline({ + label: 'cfd_metrics_pipeline', + layout: auto, + compute: { module: metricsModule, entryPoint: 'compute_metrics' }, + }), + primitives: device.createComputePipeline({ + label: 'cfd_primitives_pipeline', + layout: auto, + compute: { module: primitivesModule, entryPoint: 'conservative_to_primitive' }, + }), + reconstructXi: device.createComputePipeline({ + label: 'cfd_reconstruct_xi_pipeline', + layout: auto, + compute: { module: reconstructModule, entryPoint: 'reconstruct_xi' }, + }), + reconstructEta: device.createComputePipeline({ + label: 'cfd_reconstruct_eta_pipeline', + layout: auto, + compute: { module: reconstructModule, entryPoint: 'reconstruct_eta' }, + }), + roeFluxXi: device.createComputePipeline({ + label: 'cfd_roe_flux_xi_pipeline', + layout: auto, + compute: { module: roeFluxModule, entryPoint: 'roe_flux_xi' }, + }), + roeFluxEta: device.createComputePipeline({ + label: 'cfd_roe_flux_eta_pipeline', + layout: auto, + compute: { module: roeFluxModule, entryPoint: 'roe_flux_eta' }, + }), + residual: device.createComputePipeline({ + label: 'cfd_residual_pipeline', + layout: auto, + compute: { module: residualModule, entryPoint: 'compute_residual' }, + }), + update: device.createComputePipeline({ + label: 'cfd_update_pipeline', + layout: auto, + compute: { module: updateModule, entryPoint: 'update_solution' }, + }), + bc: device.createComputePipeline({ + label: 'cfd_bc_pipeline', + layout: auto, + compute: { module: bcModule, entryPoint: 'apply_boundary_conditions' }, + }), + forcesReduce: device.createComputePipeline({ + label: 'cfd_forces_reduce_pipeline', + layout: auto, + compute: { module: forcesModule, entryPoint: 'reduce_wall_forces' }, + }), + }; +} + +export function createCfdBindGroups( + device: GPUDevice, + pipelines: CfdPipelineSet, + buffers: CfdBufferSet, +): CfdBindGroupSet { + return { + metrics: device.createBindGroup({ + label: 'cfd_metrics_bg', + layout: pipelines.metrics.getBindGroupLayout(0), + entries: [ + { binding: 0, resource: { buffer: buffers.params } }, + { binding: 1, resource: { buffer: buffers.meshX } }, + { binding: 2, resource: { buffer: buffers.meshY } }, + { binding: 3, resource: { buffer: buffers.metrics } }, + ], + }), + primitives: device.createBindGroup({ + label: 'cfd_primitives_bg', + layout: pipelines.primitives.getBindGroupLayout(0), + entries: [ + { binding: 0, resource: { buffer: buffers.params } }, + { binding: 1, resource: { buffer: buffers.Q } }, + { binding: 2, resource: { buffer: buffers.W } }, + ], + }), + reconstruct: device.createBindGroup({ + label: 'cfd_reconstruct_bg', + layout: pipelines.reconstructXi.getBindGroupLayout(0), + entries: [ + { binding: 0, resource: { buffer: buffers.params } }, + { binding: 1, resource: { buffer: buffers.W } }, + { binding: 2, resource: { buffer: buffers.leftXi } }, + { binding: 3, resource: { buffer: buffers.rightXi } }, + { binding: 4, resource: { buffer: buffers.leftEta } }, + { binding: 5, resource: { buffer: buffers.rightEta } }, + ], + }), + roeFlux: device.createBindGroup({ + label: 'cfd_roe_flux_bg', + layout: pipelines.roeFluxXi.getBindGroupLayout(0), + entries: [ + { binding: 0, resource: { buffer: buffers.params } }, + { binding: 1, resource: { buffer: buffers.metrics } }, + { binding: 2, resource: { buffer: buffers.leftXi } }, + { binding: 3, resource: { buffer: buffers.rightXi } }, + { binding: 4, resource: { buffer: buffers.leftEta } }, + { binding: 5, resource: { buffer: buffers.rightEta } }, + { binding: 6, resource: { buffer: buffers.fluxXi } }, + { binding: 7, resource: { buffer: buffers.fluxEta } }, + ], + }), + residual: device.createBindGroup({ + label: 'cfd_residual_bg', + layout: pipelines.residual.getBindGroupLayout(0), + entries: [ + { binding: 0, resource: { buffer: buffers.params } }, + { binding: 1, resource: { buffer: buffers.metrics } }, + { binding: 2, resource: { buffer: buffers.fluxXi } }, + { binding: 3, resource: { buffer: buffers.fluxEta } }, + { binding: 4, resource: { buffer: buffers.residual } }, + ], + }), + update: device.createBindGroup({ + label: 'cfd_update_bg', + layout: pipelines.update.getBindGroupLayout(0), + entries: [ + { binding: 0, resource: { buffer: buffers.params } }, + { binding: 1, resource: { buffer: buffers.residual } }, + { binding: 2, resource: { buffer: buffers.Q } }, + { binding: 3, resource: { buffer: buffers.bcType } }, + ], + }), + bc: device.createBindGroup({ + label: 'cfd_bc_bg', + layout: pipelines.bc.getBindGroupLayout(0), + entries: [ + { binding: 0, resource: { buffer: buffers.params } }, + { binding: 1, resource: { buffer: buffers.Q } }, + { binding: 2, resource: { buffer: buffers.metrics } }, + { binding: 3, resource: { buffer: buffers.bcType } }, + { binding: 4, resource: { buffer: buffers.meshX } }, + { binding: 5, resource: { buffer: buffers.meshY } }, + ], + }), + forcesReduce: device.createBindGroup({ + label: 'cfd_forces_bg', + layout: pipelines.forcesReduce.getBindGroupLayout(0), + entries: [ + { binding: 0, resource: { buffer: buffers.params } }, + { binding: 1, resource: { buffer: buffers.Q } }, + { binding: 2, resource: { buffer: buffers.metrics } }, + { binding: 3, resource: { buffer: buffers.meshX } }, + { binding: 4, resource: { buffer: buffers.meshY } }, + { binding: 5, resource: { buffer: buffers.residual } }, + { binding: 6, resource: { buffer: buffers.forceAccum } }, + { binding: 7, resource: { buffer: buffers.residualNorm } }, + ], + }), + }; +} diff --git a/flexfoil-ui/src/lib/webgpu/cfd/CfdSolver.ts b/flexfoil-ui/src/lib/webgpu/cfd/CfdSolver.ts new file mode 100644 index 00000000..483f07c1 --- /dev/null +++ b/flexfoil-ui/src/lib/webgpu/cfd/CfdSolver.ts @@ -0,0 +1,258 @@ +/** + * WebGPU CFD Solver - Main orchestrator. + * + * Manages the GPU compute pipeline for structured-grid CFD. + * Dispatches per-timestep compute passes and handles async readback. + */ + +import { CfdConfig, packCfdParams } from './CfdConfig'; +import { CfdBufferSet, createCfdBuffers, destroyCfdBuffers } from './CfdBuffers'; +import { + CfdPipelineSet, + CfdBindGroupSet, + createCfdPipelines, + createCfdBindGroups, +} from './CfdPipelines'; + +export interface CfdSolverResult { + cl: number; + cd: number; + cm: number; + residualL2: number; + iteration: number; +} + +export class CfdSolver { + private device: GPUDevice; + private config: CfdConfig; + private buffers: CfdBufferSet; + private pipelines: CfdPipelineSet; + private bindGroups: CfdBindGroupSet; + private iteration = 0; + private _destroyed = false; + + private constructor( + device: GPUDevice, + config: CfdConfig, + buffers: CfdBufferSet, + pipelines: CfdPipelineSet, + bindGroups: CfdBindGroupSet, + ) { + this.device = device; + this.config = config; + this.buffers = buffers; + this.pipelines = pipelines; + this.bindGroups = bindGroups; + } + + /** + * Initialize the CFD solver with mesh data and configuration. + * + * @param device - WebGPU device (shared with visualization renderer) + * @param config - Solver configuration + * @param meshX - Flat f32 array of mesh x-coordinates (ni*nj) + * @param meshY - Flat f32 array of mesh y-coordinates (ni*nj) + * @param initialQ - Flat f32 array of initial conservative variables (ni*nj*5) + * @param bcTypes - Flat u32 array of boundary condition types (ni*nj) + */ + static create( + device: GPUDevice, + config: CfdConfig, + meshX: Float32Array, + meshY: Float32Array, + initialQ: Float32Array, + bcTypes: Uint32Array, + ): CfdSolver { + const buffers = createCfdBuffers(device, config.ni, config.nj); + const pipelines = createCfdPipelines(device); + const bindGroups = createCfdBindGroups(device, pipelines, buffers); + + // Upload mesh data + device.queue.writeBuffer(buffers.meshX, 0, meshX); + device.queue.writeBuffer(buffers.meshY, 0, meshY); + device.queue.writeBuffer(buffers.Q, 0, initialQ); + device.queue.writeBuffer(buffers.bcType, 0, bcTypes); + + // Upload initial params + const dt = config.cfl * 0.01 / (1.0 + config.machInf); + const paramsData = packCfdParams(config, dt, 0); + device.queue.writeBuffer(buffers.params, 0, paramsData); + + const solver = new CfdSolver(device, config, buffers, pipelines, bindGroups); + + // Compute metrics (one-time) + solver.computeMetrics(); + + return solver; + } + + /** Compute grid metrics (runs once at mesh setup). */ + private computeMetrics(): void { + const { ni, nj } = this.config; + const encoder = this.device.createCommandEncoder({ label: 'cfd_metrics' }); + + const pass = encoder.beginComputePass({ label: 'metrics' }); + pass.setPipeline(this.pipelines.metrics); + pass.setBindGroup(0, this.bindGroups.metrics); + pass.dispatchWorkgroups(Math.ceil(ni / 16), Math.ceil(nj / 16)); + pass.end(); + + this.device.queue.submit([encoder.finish()]); + } + + /** + * Run N timesteps of the solver. + * Returns a promise that resolves with forces and residual after readback. + */ + async step(nSteps: number): Promise { + if (this._destroyed) throw new Error('Solver destroyed'); + + const { ni, nj } = this.config; + const dt = this.config.cfl * 0.01 / (1.0 + this.config.machInf); + const wgXi = Math.ceil(ni / 16); + const wgEta = Math.ceil(nj / 16); + + const encoder = this.device.createCommandEncoder({ label: `cfd_steps_${this.iteration}` }); + + for (let s = 0; s < nSteps; s++) { + this.iteration++; + + // Update params uniform + const paramsData = packCfdParams(this.config, dt, this.iteration); + this.device.queue.writeBuffer(this.buffers.params, 0, paramsData); + + // Pass 1: Q -> W (conservative to primitive) + { + const pass = encoder.beginComputePass({ label: 'primitives' }); + pass.setPipeline(this.pipelines.primitives); + pass.setBindGroup(0, this.bindGroups.primitives); + pass.dispatchWorkgroups(wgXi, wgEta); + pass.end(); + } + + // Pass 2: Reconstruct xi (one workgroup per j-line) + { + const pass = encoder.beginComputePass({ label: 'reconstruct_xi' }); + pass.setPipeline(this.pipelines.reconstructXi); + pass.setBindGroup(0, this.bindGroups.reconstruct); + pass.dispatchWorkgroups(nj); + pass.end(); + } + + // Pass 3: Roe flux xi + { + const pass = encoder.beginComputePass({ label: 'roe_flux_xi' }); + pass.setPipeline(this.pipelines.roeFluxXi); + pass.setBindGroup(0, this.bindGroups.roeFlux); + pass.dispatchWorkgroups(wgXi, wgEta); + pass.end(); + } + + // Pass 4: Reconstruct eta (one workgroup per i-column) + { + const pass = encoder.beginComputePass({ label: 'reconstruct_eta' }); + pass.setPipeline(this.pipelines.reconstructEta); + pass.setBindGroup(0, this.bindGroups.reconstruct); + pass.dispatchWorkgroups(ni); + pass.end(); + } + + // Pass 5: Roe flux eta + { + const pass = encoder.beginComputePass({ label: 'roe_flux_eta' }); + pass.setPipeline(this.pipelines.roeFluxEta); + pass.setBindGroup(0, this.bindGroups.roeFlux); + pass.dispatchWorkgroups(wgXi, wgEta); + pass.end(); + } + + // Pass 6: Compute residual + { + const pass = encoder.beginComputePass({ label: 'residual' }); + pass.setPipeline(this.pipelines.residual); + pass.setBindGroup(0, this.bindGroups.residual); + pass.dispatchWorkgroups(wgXi, wgEta); + pass.end(); + } + + // Pass 7: Update solution (explicit forward Euler) + { + const pass = encoder.beginComputePass({ label: 'update' }); + pass.setPipeline(this.pipelines.update); + pass.setBindGroup(0, this.bindGroups.update); + pass.dispatchWorkgroups(wgXi, wgEta); + pass.end(); + } + + // Pass 8: Apply boundary conditions + { + const pass = encoder.beginComputePass({ label: 'bc' }); + pass.setPipeline(this.pipelines.bc); + pass.setBindGroup(0, this.bindGroups.bc); + pass.dispatchWorkgroups(Math.ceil(ni / 256)); + pass.end(); + } + } + + // Final pass: compute forces and residual + { + const pass = encoder.beginComputePass({ label: 'forces' }); + pass.setPipeline(this.pipelines.forcesReduce); + pass.setBindGroup(0, this.bindGroups.forcesReduce); + pass.dispatchWorkgroups(1); // Single workgroup reduction + pass.end(); + } + + // Copy forces and residual to readback buffer + encoder.copyBufferToBuffer(this.buffers.forceAccum, 0, this.buffers.readback, 0, 12); + encoder.copyBufferToBuffer(this.buffers.residualNorm, 0, this.buffers.readback, 12, 4); + + this.device.queue.submit([encoder.finish()]); + + // Async readback + await this.buffers.readback.mapAsync(GPUMapMode.READ); + const data = new Float32Array(this.buffers.readback.getMappedRange().slice(0)); + this.buffers.readback.unmap(); + + return { + cl: data[0], + cd: data[1], + cm: data[2], + residualL2: data[3], + iteration: this.iteration, + }; + } + + /** Get the current iteration count. */ + getIteration(): number { + return this.iteration; + } + + /** Get the Q buffer for visualization (e.g., to create velocity texture). */ + getQBuffer(): GPUBuffer { + return this.buffers.Q; + } + + /** Get mesh coordinate buffers for mesh visualization. */ + getMeshBuffers(): { x: GPUBuffer; y: GPUBuffer } { + return { x: this.buffers.meshX, y: this.buffers.meshY }; + } + + /** Update configuration (e.g., CFL change mid-run). */ + updateConfig(config: Partial): void { + Object.assign(this.config, config); + } + + /** Reset solver to initial conditions. */ + reset(initialQ: Float32Array): void { + this.iteration = 0; + this.device.queue.writeBuffer(this.buffers.Q, 0, initialQ); + } + + /** Clean up all GPU resources. */ + destroy(): void { + if (this._destroyed) return; + this._destroyed = true; + destroyCfdBuffers(this.buffers); + } +} diff --git a/flexfoil-ui/src/lib/webgpu/cfd/index.ts b/flexfoil-ui/src/lib/webgpu/cfd/index.ts new file mode 100644 index 00000000..a70413a3 --- /dev/null +++ b/flexfoil-ui/src/lib/webgpu/cfd/index.ts @@ -0,0 +1,11 @@ +/** + * CFD Solver - WebGPU-accelerated 2D compressible flow solver. + * + * Solves Euler / Navier-Stokes / RANS-SA equations on structured O-grids + * using Roe approximate Riemann solver with MUSCL/WENO5 reconstruction. + */ + +export { CfdSolver } from './CfdSolver'; +export type { CfdSolverResult } from './CfdSolver'; +export { DEFAULT_CFD_CONFIG, packCfdParams, NVAR } from './CfdConfig'; +export type { CfdConfig, PhysicsMode, ReconstructionMode, TimeSteppingMode } from './CfdConfig'; diff --git a/flexfoil-ui/src/lib/webgpu/shaders/cfd/bc.wgsl b/flexfoil-ui/src/lib/webgpu/shaders/cfd/bc.wgsl new file mode 100644 index 00000000..2bd7b328 --- /dev/null +++ b/flexfoil-ui/src/lib/webgpu/shaders/cfd/bc.wgsl @@ -0,0 +1,125 @@ +// Boundary condition application for the CFD solver. +// +// Handles: +// - Wall BC (j=0): slip wall (Euler) or no-slip (NS/RANS) +// - Far-field BC (j=nj-1): characteristic-based +// - Periodicity in i-direction is implicit in the O-grid topology + +@group(0) @binding(0) var params: CfdParams; +@group(0) @binding(1) var Q: array; +@group(0) @binding(2) var metrics: array; +@group(0) @binding(3) var bc_type: array; +@group(0) @binding(4) var mesh_x: array; +@group(0) @binding(5) var mesh_y: array; + +@compute @workgroup_size(256, 1, 1) +fn apply_boundary_conditions(@builtin(global_invocation_id) gid: vec3) { + let i = gid.x; + let ni = params.ni; + let nj = params.nj; + + if (i >= ni) { + return; + } + + let gamma = params.gamma; + + // === Wall boundary (j=0) === + { + let j = 0u; + let cell = cell_idx(i, j, ni); + + // Get interior cell values (j=1) + let base_int = q_idx(i, 1u, 0u, ni); + let rho_int = Q[base_int + 0u]; + let rhou_int = Q[base_int + 1u]; + let rhov_int = Q[base_int + 2u]; + let E_int = Q[base_int + 3u]; + let nu_int = Q[base_int + 4u]; + + // Surface normal (eta-direction at wall, pointing outward from body) + let m_base = cell * 5u; + let eta_x = metrics[m_base + 2u]; + let eta_y = metrics[m_base + 3u]; + let mag = sqrt(eta_x * eta_x + eta_y * eta_y); + let nx = select(eta_x / mag, 0.0, mag < 1e-15); + let ny = select(eta_y / mag, 0.0, mag < 1e-15); + + let u_int = rhou_int / rho_int; + let v_int = rhov_int / rho_int; + + // Reflect normal velocity component + let Vn = u_int * nx + v_int * ny; + + var u_wall: f32; + var v_wall: f32; + + if (params.physics_mode == 0u) { + // Euler: slip wall (zero normal velocity, keep tangential) + u_wall = u_int - 2.0 * Vn * nx; + v_wall = v_int - 2.0 * Vn * ny; + } else { + // NS/RANS: no-slip wall (zero velocity) + u_wall = -u_int; // Ghost cell for zero at wall + v_wall = -v_int; + } + + let base_wall = q_idx(i, 0u, 0u, ni); + Q[base_wall + 0u] = rho_int; // Extrapolate density + Q[base_wall + 1u] = rho_int * u_wall; + Q[base_wall + 2u] = rho_int * v_wall; + // Extrapolate pressure, set energy accordingly + let p_int = (gamma - 1.0) * (E_int - 0.5 * rho_int * (u_int * u_int + v_int * v_int)); + Q[base_wall + 3u] = p_int / (gamma - 1.0) + 0.5 * rho_int * (u_wall * u_wall + v_wall * v_wall); + // SA: zero at wall + Q[base_wall + 4u] = select(nu_int, 0.0, params.physics_mode == 2u); + } + + // === Far-field boundary (j=nj-1) === + { + let j = nj - 1u; + + // Freestream state + let mach = params.mach_inf; + let alpha = params.alpha; + let rho_inf: f32 = 1.0; + let u_inf = mach * cos(alpha); + let v_inf = mach * sin(alpha); + let p_inf = 1.0 / gamma; + let E_inf = p_inf / (gamma - 1.0) + 0.5 * rho_inf * (u_inf * u_inf + v_inf * v_inf); + + // Simple characteristic-based: use freestream for inflow, extrapolate for outflow + let base_int = q_idx(i, j - 1u, 0u, ni); + let rho_int = Q[base_int + 0u]; + let u_int = Q[base_int + 1u] / rho_int; + let v_int = Q[base_int + 2u] / rho_int; + + // Radial direction at far-field + let m_base = cell_idx(i, j, ni) * 5u; + let eta_x = metrics[m_base + 2u]; + let eta_y = metrics[m_base + 3u]; + let mag = sqrt(eta_x * eta_x + eta_y * eta_y); + let nx = select(eta_x / mag, 0.0, mag < 1e-15); + let ny = select(eta_y / mag, 0.0, mag < 1e-15); + + // Check if inflow or outflow based on interior velocity + let Vn = u_int * nx + v_int * ny; + + let base_ff = q_idx(i, j, 0u, ni); + if (Vn >= 0.0) { + // Outflow: extrapolate from interior + Q[base_ff + 0u] = Q[base_int + 0u]; + Q[base_ff + 1u] = Q[base_int + 1u]; + Q[base_ff + 2u] = Q[base_int + 2u]; + Q[base_ff + 3u] = Q[base_int + 3u]; + Q[base_ff + 4u] = Q[base_int + 4u]; + } else { + // Inflow: set to freestream + Q[base_ff + 0u] = rho_inf; + Q[base_ff + 1u] = rho_inf * u_inf; + Q[base_ff + 2u] = rho_inf * v_inf; + Q[base_ff + 3u] = E_inf; + Q[base_ff + 4u] = 0.0; + } + } +} diff --git a/flexfoil-ui/src/lib/webgpu/shaders/cfd/cfd_common.wgsl b/flexfoil-ui/src/lib/webgpu/shaders/cfd/cfd_common.wgsl new file mode 100644 index 00000000..aacc3c17 --- /dev/null +++ b/flexfoil-ui/src/lib/webgpu/shaders/cfd/cfd_common.wgsl @@ -0,0 +1,42 @@ +// Common struct definitions for CFD compute shaders. +// This file is prepended to each shader at pipeline creation time. + +struct CfdParams { + ni: u32, + nj: u32, + gamma: f32, + cfl: f32, + mach_inf: f32, + alpha: f32, + reynolds: f32, + prandtl: f32, + dt: f32, + iteration: u32, + physics_mode: u32, // 0=Euler, 1=LaminarNS, 2=RANS_SA + reconstruction: u32, // 0=MUSCL, 1=WENO5 +}; + +// Boundary condition type constants +const BC_INTERIOR: u32 = 0u; +const BC_WALL: u32 = 1u; +const BC_FARFIELD: u32 = 2u; +const BC_WAKE: u32 = 3u; + +// Number of conservative variables per cell +const NVAR: u32 = 5u; + +// Helper: flat index for cell (i, j) in a ni x nj grid +fn cell_idx(i: u32, j: u32, ni: u32) -> u32 { + return j * ni + i; +} + +// Helper: flat index for variable `v` at cell (i, j) +fn q_idx(i: u32, j: u32, v: u32, ni: u32) -> u32 { + return (j * ni + i) * NVAR + v; +} + +// Helper: wrap i-index for periodic O-grid +fn wrap_i(i: i32, ni: u32) -> u32 { + let n = i32(ni); + return u32(((i % n) + n) % n); +} diff --git a/flexfoil-ui/src/lib/webgpu/shaders/cfd/forces_reduce.wgsl b/flexfoil-ui/src/lib/webgpu/shaders/cfd/forces_reduce.wgsl new file mode 100644 index 00000000..f0cefc09 --- /dev/null +++ b/flexfoil-ui/src/lib/webgpu/shaders/cfd/forces_reduce.wgsl @@ -0,0 +1,115 @@ +// Parallel reduction to compute aerodynamic force coefficients. +// Integrates pressure and (optionally) shear stress along the wall (j=0). +// +// Outputs: [Cl, Cd, Cm] in the force_accum buffer. +// Also computes L2 residual norm in residual_norm buffer. + +@group(0) @binding(0) var params: CfdParams; +@group(0) @binding(1) var Q: array; +@group(0) @binding(2) var metrics: array; +@group(0) @binding(3) var mesh_x: array; +@group(0) @binding(4) var mesh_y: array; +@group(0) @binding(5) var residual: array; +@group(0) @binding(6) var force_accum: array; // [Cl, Cd, Cm] +@group(0) @binding(7) var residual_norm: array; // [L2_rho, L2_rhou, L2_rhov, L2_E] + +var shared_cl: array; +var shared_cd: array; +var shared_cm: array; +var shared_res: array; + +@compute @workgroup_size(256, 1, 1) +fn reduce_wall_forces(@builtin(local_invocation_id) lid: vec3, + @builtin(workgroup_id) wid: vec3, + @builtin(num_workgroups) nwg: vec3) { + let tid = lid.x; + let ni = params.ni; + let nj = params.nj; + let gamma = params.gamma; + let mach = params.mach_inf; + let alpha = params.alpha; + + // Reference quantities for force coefficients + let p_inf = 1.0 / gamma; + let q_inf = 0.5 * gamma * p_inf * mach * mach; // Dynamic pressure + + // Each thread processes multiple wall cells + var local_cl: f32 = 0.0; + var local_cd: f32 = 0.0; + var local_cm: f32 = 0.0; + var local_res_sq: f32 = 0.0; + + var idx = tid; + while (idx < ni) { + // Surface pressure at wall (j=0) + let base = q_idx(idx, 0u, 0u, ni); + let rho = Q[base + 0u]; + let rhou = Q[base + 1u]; + let rhov = Q[base + 2u]; + let E = Q[base + 3u]; + let u = rhou / max(rho, 1e-10); + let v = rhov / max(rho, 1e-10); + let p = (gamma - 1.0) * (E - 0.5 * rho * (u * u + v * v)); + let cp = (p - p_inf) / q_inf; + + // Surface element: ds = |dx/di| at wall + let ip = select(idx + 1u, 0u, idx + 1u >= ni); + let im = select(idx - 1u, ni - 1u, idx == 0u); + let dx = 0.5 * (mesh_x[cell_idx(ip, 0u, ni)] - mesh_x[cell_idx(im, 0u, ni)]); + let dy = 0.5 * (mesh_y[cell_idx(ip, 0u, ni)] - mesh_y[cell_idx(im, 0u, ni)]); + + // Outward normal (rotate tangent 90 degrees: tangent (dx,dy) -> normal (-dy, dx)) + // Convention: positive Cl is lift (upward) + let fx = -p * (-dy); // Force in x per unit span (pressure * normal_x * ds) + let fy = -p * dx; // Force in y + + // Rotate to wind axes + let ca = cos(alpha); + let sa = sin(alpha); + let drag_contrib = fx * ca + fy * sa; + let lift_contrib = -fx * sa + fy * ca; + + // Moment about (0.25, 0) - quarter chord + let xc = mesh_x[cell_idx(idx, 0u, ni)] - 0.25; + let yc = mesh_y[cell_idx(idx, 0u, ni)]; + let moment_contrib = xc * fy - yc * fx; // Pitching moment + + local_cd += drag_contrib / q_inf; + local_cl += lift_contrib / q_inf; + local_cm += moment_contrib / q_inf; + + // Accumulate residual L2 norm (all interior cells in this thread's range) + for (var j: u32 = 1u; j < nj - 1u; j = j + 1u) { + let r0 = residual[q_idx(idx, j, 0u, ni)]; + local_res_sq += r0 * r0; + } + + idx += 256u; + } + + shared_cl[tid] = local_cl; + shared_cd[tid] = local_cd; + shared_cm[tid] = local_cm; + shared_res[tid] = local_res_sq; + workgroupBarrier(); + + // Parallel reduction + var stride: u32 = 128u; + while (stride > 0u) { + if (tid < stride) { + shared_cl[tid] += shared_cl[tid + stride]; + shared_cd[tid] += shared_cd[tid + stride]; + shared_cm[tid] += shared_cm[tid + stride]; + shared_res[tid] += shared_res[tid + stride]; + } + workgroupBarrier(); + stride = stride >> 1u; + } + + if (tid == 0u) { + force_accum[0] = shared_cl[0]; + force_accum[1] = shared_cd[0]; + force_accum[2] = shared_cm[0]; + residual_norm[0] = sqrt(shared_res[0] / f32(ni * (nj - 2u))); + } +} diff --git a/flexfoil-ui/src/lib/webgpu/shaders/cfd/metrics.wgsl b/flexfoil-ui/src/lib/webgpu/shaders/cfd/metrics.wgsl new file mode 100644 index 00000000..057a651e --- /dev/null +++ b/flexfoil-ui/src/lib/webgpu/shaders/cfd/metrics.wgsl @@ -0,0 +1,58 @@ +// Compute grid metrics from mesh node coordinates. +// Runs once at mesh setup (not per timestep). +// +// Computes: xi_x, xi_y, eta_x, eta_y, J (Jacobian) +// Using central differences on the (x,y) coordinates. + +@group(0) @binding(0) var params: CfdParams; +@group(0) @binding(1) var mesh_x: array; +@group(0) @binding(2) var mesh_y: array; +@group(0) @binding(3) var metrics: array; +// metrics layout: [xi_x, xi_y, eta_x, eta_y, J] per cell = 5 floats + +@compute @workgroup_size(16, 16, 1) +fn compute_metrics(@builtin(global_invocation_id) gid: vec3) { + let i = gid.x; + let j = gid.y; + let ni = params.ni; + let nj = params.nj; + + if (i >= ni || j >= nj) { + return; + } + + // Compute x_xi, y_xi (derivatives w.r.t. xi = i-direction) + let ip = wrap_i(i32(i) + 1, ni); + let im = wrap_i(i32(i) - 1, ni); + + let x_xi = 0.5 * (mesh_x[cell_idx(ip, j, ni)] - mesh_x[cell_idx(im, j, ni)]); + let y_xi = 0.5 * (mesh_y[cell_idx(ip, j, ni)] - mesh_y[cell_idx(im, j, ni)]); + + // Compute x_eta, y_eta (derivatives w.r.t. eta = j-direction) + var x_eta: f32; + var y_eta: f32; + if (j == 0u) { + x_eta = mesh_x[cell_idx(i, 1u, ni)] - mesh_x[cell_idx(i, 0u, ni)]; + y_eta = mesh_y[cell_idx(i, 1u, ni)] - mesh_y[cell_idx(i, 0u, ni)]; + } else if (j == nj - 1u) { + x_eta = mesh_x[cell_idx(i, nj - 1u, ni)] - mesh_x[cell_idx(i, nj - 2u, ni)]; + y_eta = mesh_y[cell_idx(i, nj - 1u, ni)] - mesh_y[cell_idx(i, nj - 2u, ni)]; + } else { + x_eta = 0.5 * (mesh_x[cell_idx(i, j + 1u, ni)] - mesh_x[cell_idx(i, j - 1u, ni)]); + y_eta = 0.5 * (mesh_y[cell_idx(i, j + 1u, ni)] - mesh_y[cell_idx(i, j - 1u, ni)]); + } + + // Jacobian of the transformation: J = x_xi * y_eta - x_eta * y_xi + let J = x_xi * y_eta - x_eta * y_xi; + let inv_J = select(1.0 / J, 1e10, abs(J) < 1e-15); + + // Metric terms: inverse mapping + // xi_x = y_eta / J, xi_y = -x_eta / J + // eta_x = -y_xi / J, eta_y = x_xi / J + let base = cell_idx(i, j, ni) * 5u; + metrics[base + 0u] = y_eta * inv_J; // xi_x + metrics[base + 1u] = -x_eta * inv_J; // xi_y + metrics[base + 2u] = -y_xi * inv_J; // eta_x + metrics[base + 3u] = x_xi * inv_J; // eta_y + metrics[base + 4u] = J; // Jacobian +} diff --git a/flexfoil-ui/src/lib/webgpu/shaders/cfd/primitives.wgsl b/flexfoil-ui/src/lib/webgpu/shaders/cfd/primitives.wgsl new file mode 100644 index 00000000..224c29a4 --- /dev/null +++ b/flexfoil-ui/src/lib/webgpu/shaders/cfd/primitives.wgsl @@ -0,0 +1,37 @@ +// Convert conservative variables Q to primitive variables W. +// Q = [rho, rho*u, rho*v, E, nu_tilde] +// W = [rho, u, v, p, nu_tilde] + +@group(0) @binding(0) var params: CfdParams; +@group(0) @binding(1) var Q: array; +@group(0) @binding(2) var W: array; + +@compute @workgroup_size(16, 16, 1) +fn conservative_to_primitive(@builtin(global_invocation_id) gid: vec3) { + let i = gid.x; + let j = gid.y; + let ni = params.ni; + let nj = params.nj; + + if (i >= ni || j >= nj) { + return; + } + + let base = q_idx(i, j, 0u, ni); + let rho = max(Q[base + 0u], 1e-10); + let rhou = Q[base + 1u]; + let rhov = Q[base + 2u]; + let E = Q[base + 3u]; + let nu_t = Q[base + 4u]; + + let inv_rho = 1.0 / rho; + let u = rhou * inv_rho; + let v = rhov * inv_rho; + let p = max((params.gamma - 1.0) * (E - 0.5 * rho * (u * u + v * v)), 1e-10); + + W[base + 0u] = rho; + W[base + 1u] = u; + W[base + 2u] = v; + W[base + 3u] = p; + W[base + 4u] = nu_t; +} diff --git a/flexfoil-ui/src/lib/webgpu/shaders/cfd/reconstruct_muscl.wgsl b/flexfoil-ui/src/lib/webgpu/shaders/cfd/reconstruct_muscl.wgsl new file mode 100644 index 00000000..7f6be08c --- /dev/null +++ b/flexfoil-ui/src/lib/webgpu/shaders/cfd/reconstruct_muscl.wgsl @@ -0,0 +1,109 @@ +// 2nd-order MUSCL reconstruction with minmod limiter. +// Reconstructs left/right states at cell faces for the Roe solver. +// +// Two entry points: reconstruct_xi (along i-direction) and reconstruct_eta (along j-direction). +// Each workgroup processes one grid line. + +@group(0) @binding(0) var params: CfdParams; +@group(0) @binding(1) var W: array; +@group(0) @binding(2) var left_xi: array; +@group(0) @binding(3) var right_xi: array; +@group(0) @binding(4) var left_eta: array; +@group(0) @binding(5) var right_eta: array; + +fn minmod(a: f32, b: f32) -> f32 { + if (a * b <= 0.0) { + return 0.0; + } + return select(b, a, abs(a) < abs(b)); +} + +// Get primitive variable v at wrapped cell (i, j) +fn get_w(i: i32, j: i32, v: u32, ni: u32, nj: u32) -> f32 { + let ii = wrap_i(i, ni); + let jj = clamp(u32(j), 0u, nj - 1u); + return W[q_idx(ii, jj, v, ni)]; +} + +// Reconstruct along xi (i-direction) for a fixed j-line. +// Face (i+1/2, j) has left state from cell (i) and right state from cell (i+1). +@compute @workgroup_size(1, 1, 1) +fn reconstruct_xi(@builtin(global_invocation_id) gid: vec3) { + let j = gid.x; + let ni = params.ni; + let nj = params.nj; + + if (j >= nj) { + return; + } + + // Process each face along this j-line + for (var i: u32 = 0u; i < ni; i = i + 1u) { + for (var v: u32 = 0u; v < NVAR; v = v + 1u) { + // Get stencil values + let wim1 = get_w(i32(i) - 1, i32(j), v, ni, nj); + let wi = get_w(i32(i), i32(j), v, ni, nj); + let wip1 = get_w(i32(i) + 1, i32(j), v, ni, nj); + let wip2 = get_w(i32(i) + 2, i32(j), v, ni, nj); + + // Slopes + let dL = wi - wim1; + let dR = wip1 - wi; + let dL2 = wip1 - wi; + let dR2 = wip2 - wip1; + + // MUSCL reconstruction with minmod + let slope_L = minmod(dL, dR); + let slope_R = minmod(dL2, dR2); + + // Left state at face i+1/2 (extrapolated from cell i) + let wL = wi + 0.5 * slope_L; + // Right state at face i+1/2 (extrapolated from cell i+1) + let wR = wip1 - 0.5 * slope_R; + + // Store: face index = j * ni + i (face between cell i and i+1) + let face_idx = (j * ni + i) * NVAR + v; + left_xi[face_idx] = wL; + right_xi[face_idx] = wR; + } + } +} + +// Reconstruct along eta (j-direction) for a fixed i-column. +// Face (i, j+1/2) has left state from cell (i,j) and right state from cell (i,j+1). +@compute @workgroup_size(1, 1, 1) +fn reconstruct_eta(@builtin(global_invocation_id) gid: vec3) { + let i = gid.x; + let ni = params.ni; + let nj = params.nj; + + if (i >= ni) { + return; + } + + // Process each face along this i-column + for (var j: u32 = 0u; j < nj - 1u; j = j + 1u) { + for (var v: u32 = 0u; v < NVAR; v = v + 1u) { + let wjm1 = get_w(i32(i), i32(j) - 1, v, ni, nj); + let wj = get_w(i32(i), i32(j), v, ni, nj); + let wjp1 = get_w(i32(i), i32(j) + 1, v, ni, nj); + let wjp2 = get_w(i32(i), i32(j) + 2, v, ni, nj); + + let dL = wj - wjm1; + let dR = wjp1 - wj; + let dL2 = wjp1 - wj; + let dR2 = wjp2 - wjp1; + + let slope_L = minmod(dL, dR); + let slope_R = minmod(dL2, dR2); + + let wL = wj + 0.5 * slope_L; + let wR = wjp1 - 0.5 * slope_R; + + // Face index: j * ni + i (face between cell j and j+1 at column i) + let face_idx = (j * ni + i) * NVAR + v; + left_eta[face_idx] = wL; + right_eta[face_idx] = wR; + } + } +} diff --git a/flexfoil-ui/src/lib/webgpu/shaders/cfd/residual.wgsl b/flexfoil-ui/src/lib/webgpu/shaders/cfd/residual.wgsl new file mode 100644 index 00000000..78643043 --- /dev/null +++ b/flexfoil-ui/src/lib/webgpu/shaders/cfd/residual.wgsl @@ -0,0 +1,51 @@ +// Compute the RHS residual: R(Q) = -(1/J) * (dF/dxi + dG/deta) +// For explicit time stepping: Q_new = Q_old + dt * R(Q) + +@group(0) @binding(0) var params: CfdParams; +@group(0) @binding(1) var metrics: array; +@group(0) @binding(2) var flux_xi: array; +@group(0) @binding(3) var flux_eta: array; +@group(0) @binding(4) var residual: array; + +@compute @workgroup_size(16, 16, 1) +fn compute_residual(@builtin(global_invocation_id) gid: vec3) { + let i = gid.x; + let j = gid.y; + let ni = params.ni; + let nj = params.nj; + + if (i >= ni || j >= nj) { + return; + } + + // Get Jacobian at this cell + let m_base = cell_idx(i, j, ni) * 5u; + let J = metrics[m_base + 4u]; + let inv_J = select(1.0 / J, 0.0, abs(J) < 1e-15); + + for (var v: u32 = 0u; v < NVAR; v = v + 1u) { + // xi-direction flux difference: F(i+1/2) - F(i-1/2) + // Face (i) stores flux at face i+1/2 + let face_ip = (j * ni + i) * NVAR + v; + let im = wrap_i(i32(i) - 1, ni); + let face_im = (j * ni + im) * NVAR + v; + + let dF_dxi = flux_xi[face_ip] - flux_xi[face_im]; + + // eta-direction flux difference: G(j+1/2) - G(j-1/2) + var dG_deta: f32 = 0.0; + if (j > 0u && j < nj - 1u) { + let face_jp = (j * ni + i) * NVAR + v; + let face_jm = ((j - 1u) * ni + i) * NVAR + v; + dG_deta = flux_eta[face_jp] - flux_eta[face_jm]; + } else if (j == 0u) { + // At wall: only outward flux + let face_jp = (0u * ni + i) * NVAR + v; + dG_deta = flux_eta[face_jp]; // G(1/2) - 0 (wall flux handled by BC) + } + + // Residual: R = -(1/J) * (dF/dxi + dG/deta) + let r_idx = q_idx(i, j, v, ni); + residual[r_idx] = -inv_J * (dF_dxi + dG_deta); + } +} diff --git a/flexfoil-ui/src/lib/webgpu/shaders/cfd/roe_flux.wgsl b/flexfoil-ui/src/lib/webgpu/shaders/cfd/roe_flux.wgsl new file mode 100644 index 00000000..ee2cec01 --- /dev/null +++ b/flexfoil-ui/src/lib/webgpu/shaders/cfd/roe_flux.wgsl @@ -0,0 +1,185 @@ +// Roe approximate Riemann solver. +// Computes numerical flux at each cell face from reconstructed left/right primitive states. +// +// Two entry points: roe_flux_xi (xi-direction faces) and roe_flux_eta (eta-direction faces). + +@group(0) @binding(0) var params: CfdParams; +@group(0) @binding(1) var metrics: array; +@group(0) @binding(2) var left_xi: array; +@group(0) @binding(3) var right_xi: array; +@group(0) @binding(4) var left_eta: array; +@group(0) @binding(5) var right_eta: array; +@group(0) @binding(6) var flux_xi: array; +@group(0) @binding(7) var flux_eta: array; + +const EPSROE: f32 = 0.1; + +// Compute Roe flux in a given direction (nx, ny = metric-weighted face normal). +// W_L, W_R are primitive states [rho, u, v, p, nu_tilde]. +// Returns flux vector [F0, F1, F2, F3, F4]. +fn roe_flux_1d( + rhoL: f32, uL: f32, vL: f32, pL: f32, nuL: f32, + rhoR: f32, uR: f32, vR: f32, pR: f32, nuR: f32, + nx: f32, ny: f32, gamma: f32 +) -> array { + // Roe-averaged quantities + let sqrtL = sqrt(max(rhoL, 1e-10)); + let sqrtR = sqrt(max(rhoR, 1e-10)); + let inv_sum = 1.0 / (sqrtL + sqrtR); + + let u_roe = (sqrtL * uL + sqrtR * uR) * inv_sum; + let v_roe = (sqrtL * vL + sqrtR * vR) * inv_sum; + + let hL = (gamma * pL / ((gamma - 1.0) * rhoL)) + 0.5 * (uL * uL + vL * vL); + let hR = (gamma * pR / ((gamma - 1.0) * rhoR)) + 0.5 * (uR * uR + vR * vR); + let h_roe = (sqrtL * hL + sqrtR * hR) * inv_sum; + + let q2_roe = u_roe * u_roe + v_roe * v_roe; + let c2_roe = max((gamma - 1.0) * (h_roe - 0.5 * q2_roe), 1e-10); + let c_roe = sqrt(c2_roe); + + // Contravariant velocity + let face_mag = sqrt(nx * nx + ny * ny); + let face_inv = select(1.0 / face_mag, 0.0, face_mag < 1e-15); + let nx_hat = nx * face_inv; + let ny_hat = ny * face_inv; + + let V_roe = u_roe * nx_hat + v_roe * ny_hat; + + // Eigenvalues with entropy fix + let lam1 = abs(V_roe - c_roe * face_mag); + let lam2 = abs(V_roe); + let lam3 = abs(V_roe + c_roe * face_mag); + + let eps = EPSROE * c_roe * face_mag; + let l1 = select(lam1, (lam1 * lam1 + eps * eps) / (2.0 * eps), lam1 < eps); + let l2 = select(lam2, (lam2 * lam2 + eps * eps) / (2.0 * eps), lam2 < eps); + let l3 = select(lam3, (lam3 * lam3 + eps * eps) / (2.0 * eps), lam3 < eps); + + // Jump in primitive variables + let drho = rhoR - rhoL; + let du = uR - uL; + let dv = vR - vL; + let dp = pR - pL; + let dnu = nuR - nuL; + + let rho_roe = sqrtL * sqrtR; + let dVn = du * nx_hat + dv * ny_hat; + + // Wave strengths + let w1 = (dp - rho_roe * c_roe * face_mag * dVn) / (2.0 * c2_roe); + let w2 = drho - dp / c2_roe; + let w3 = (dp + rho_roe * c_roe * face_mag * dVn) / (2.0 * c2_roe); + + // Physical flux from left and right states + let VnL = uL * nx + vL * ny; + let VnR = uR * nx + vR * ny; + + let eL = pL / (gamma - 1.0) + 0.5 * rhoL * (uL * uL + vL * vL); + let eR = pR / (gamma - 1.0) + 0.5 * rhoR * (uR * uR + vR * vR); + + // Left flux + var fL: array; + fL[0] = rhoL * VnL; + fL[1] = rhoL * uL * VnL + pL * nx; + fL[2] = rhoL * vL * VnL + pL * ny; + fL[3] = (eL + pL) * VnL; + fL[4] = rhoL * nuL * VnL; + + // Right flux + var fR: array; + fR[0] = rhoR * VnR; + fR[1] = rhoR * uR * VnR + pR * nx; + fR[2] = rhoR * vR * VnR + pR * ny; + fR[3] = (eR + pR) * VnR; + fR[4] = rhoR * nuR * VnR; + + // Roe dissipation terms + // Wave 1 (u - c) + var d: array; + d[0] = l1 * w1; + d[1] = l1 * w1 * (u_roe - c_roe * nx_hat); + d[2] = l1 * w1 * (v_roe - c_roe * ny_hat); + d[3] = l1 * w1 * (h_roe - c_roe * face_mag * V_roe); + d[4] = 0.0; + + // Wave 2 (u) - entropy wave + shear + d[0] += l2 * w2; + d[1] += l2 * (w2 * u_roe + rho_roe * (du - dVn * nx_hat)); + d[2] += l2 * (w2 * v_roe + rho_roe * (dv - dVn * ny_hat)); + d[3] += l2 * (w2 * 0.5 * q2_roe + rho_roe * (u_roe * du + v_roe * dv - V_roe * dVn)); + d[4] += l2 * (drho * nuR + rho_roe * dnu); // SA variable: passive advection + + // Wave 3 (u + c) + d[0] += l3 * w3; + d[1] += l3 * w3 * (u_roe + c_roe * nx_hat); + d[2] += l3 * w3 * (v_roe + c_roe * ny_hat); + d[3] += l3 * w3 * (h_roe + c_roe * face_mag * V_roe); + + // Roe flux = 0.5 * (F_L + F_R - dissipation) + var flux: array; + for (var k: u32 = 0u; k < 5u; k = k + 1u) { + flux[k] = 0.5 * (fL[k] + fR[k] - d[k]); + } + + return flux; +} + +@compute @workgroup_size(16, 16, 1) +fn roe_flux_xi(@builtin(global_invocation_id) gid: vec3) { + let i = gid.x; + let j = gid.y; + let ni = params.ni; + let nj = params.nj; + + if (i >= ni || j >= nj) { + return; + } + + // Get left/right primitive states at face (i+1/2, j) + let face = (j * ni + i) * NVAR; + let rhoL = left_xi[face + 0u]; let uL = left_xi[face + 1u]; + let vL = left_xi[face + 2u]; let pL = left_xi[face + 3u]; let nuL = left_xi[face + 4u]; + let rhoR = right_xi[face + 0u]; let uR = right_xi[face + 1u]; + let vR = right_xi[face + 2u]; let pR = right_xi[face + 3u]; let nuR = right_xi[face + 4u]; + + // Face normal in xi-direction: (xi_x, xi_y) from metrics at cell (i, j) + let m_base = cell_idx(i, j, ni) * 5u; + let nx = metrics[m_base + 0u]; // xi_x + let ny = metrics[m_base + 1u]; // xi_y + + let flux = roe_flux_1d(rhoL, uL, vL, pL, nuL, rhoR, uR, vR, pR, nuR, nx, ny, params.gamma); + + for (var k: u32 = 0u; k < 5u; k = k + 1u) { + flux_xi[face + k] = flux[k]; + } +} + +@compute @workgroup_size(16, 16, 1) +fn roe_flux_eta(@builtin(global_invocation_id) gid: vec3) { + let i = gid.x; + let j = gid.y; + let ni = params.ni; + let nj = params.nj; + + if (i >= ni || j >= nj - 1u) { + return; + } + + let face = (j * ni + i) * NVAR; + let rhoL = left_eta[face + 0u]; let uL = left_eta[face + 1u]; + let vL = left_eta[face + 2u]; let pL = left_eta[face + 3u]; let nuL = left_eta[face + 4u]; + let rhoR = right_eta[face + 0u]; let uR = right_eta[face + 1u]; + let vR = right_eta[face + 2u]; let pR = right_eta[face + 3u]; let nuR = right_eta[face + 4u]; + + // Face normal in eta-direction: (eta_x, eta_y) from metrics at cell (i, j) + let m_base = cell_idx(i, j, ni) * 5u; + let nx = metrics[m_base + 2u]; // eta_x + let ny = metrics[m_base + 3u]; // eta_y + + let flux = roe_flux_1d(rhoL, uL, vL, pL, nuL, rhoR, uR, vR, pR, nuR, nx, ny, params.gamma); + + for (var k: u32 = 0u; k < 5u; k = k + 1u) { + flux_eta[face + k] = flux[k]; + } +} diff --git a/flexfoil-ui/src/lib/webgpu/shaders/cfd/update.wgsl b/flexfoil-ui/src/lib/webgpu/shaders/cfd/update.wgsl new file mode 100644 index 00000000..419ffd93 --- /dev/null +++ b/flexfoil-ui/src/lib/webgpu/shaders/cfd/update.wgsl @@ -0,0 +1,47 @@ +// Update solution: Q_new = Q_old + dt * residual +// For explicit forward Euler time stepping. + +@group(0) @binding(0) var params: CfdParams; +@group(0) @binding(1) var residual: array; +@group(0) @binding(2) var Q: array; +@group(0) @binding(3) var bc_type: array; + +@compute @workgroup_size(16, 16, 1) +fn update_solution(@builtin(global_invocation_id) gid: vec3) { + let i = gid.x; + let j = gid.y; + let ni = params.ni; + let nj = params.nj; + + if (i >= ni || j >= nj) { + return; + } + + // Only update interior cells (boundaries handled by BC shader) + let cell = cell_idx(i, j, ni); + let bc = bc_type[cell]; + if (bc != BC_INTERIOR) { + return; + } + + let dt = params.dt; + for (var v: u32 = 0u; v < NVAR; v = v + 1u) { + let idx = q_idx(i, j, v, ni); + Q[idx] = Q[idx] + dt * residual[idx]; + } + + // Enforce positivity of density and pressure + let base = q_idx(i, j, 0u, ni); + Q[base + 0u] = max(Q[base + 0u], 1e-10); // rho > 0 + + // Recompute pressure to ensure positivity + let rho = Q[base + 0u]; + let u = Q[base + 1u] / rho; + let v2 = Q[base + 2u] / rho; + let E = Q[base + 3u]; + let p = (params.gamma - 1.0) * (E - 0.5 * rho * (u * u + v2 * v2)); + if (p < 1e-10) { + // Reset energy to maintain minimum pressure + Q[base + 3u] = 1e-10 / (params.gamma - 1.0) + 0.5 * rho * (u * u + v2 * v2); + } +} diff --git a/flexfoil-ui/src/stores/cfdStore.ts b/flexfoil-ui/src/stores/cfdStore.ts new file mode 100644 index 00000000..928663b8 --- /dev/null +++ b/flexfoil-ui/src/stores/cfdStore.ts @@ -0,0 +1,136 @@ +/** + * Zustand store for CFD solver state management. + * + * Manages grid settings, physics mode, solver runtime state, + * convergence history, and force coefficient history. + */ + +import { create } from 'zustand'; +import type { + CfdConfig, + PhysicsMode, + ReconstructionMode, + TimeSteppingMode, +} from '../lib/webgpu/cfd/CfdConfig'; +import { DEFAULT_CFD_CONFIG } from '../lib/webgpu/cfd/CfdConfig'; + +export interface ConvergencePoint { + iteration: number; + residualL2: number; +} + +export interface ForcePoint { + iteration: number; + cl: number; + cd: number; + cm: number; +} + +interface CfdState { + // Configuration + config: CfdConfig; + + // Runtime state + isRunning: boolean; + iteration: number; + converged: boolean; + meshGenerated: boolean; + + // History for plots + convergenceHistory: ConvergencePoint[]; + forceHistory: ForcePoint[]; + + // Latest values + cl: number; + cd: number; + cm: number; + residualL2: number; +} + +interface CfdActions { + // Config setters + setNi: (ni: number) => void; + setNj: (nj: number) => void; + setMach: (mach: number) => void; + setAlphaDeg: (alpha: number) => void; + setReynolds: (re: number) => void; + setCfl: (cfl: number) => void; + setPhysics: (mode: PhysicsMode) => void; + setReconstruction: (mode: ReconstructionMode) => void; + setTimeStepping: (mode: TimeSteppingMode) => void; + setFarField: (dist: number) => void; + setDs0: (ds0: number) => void; + + // Runtime actions + setRunning: (running: boolean) => void; + setMeshGenerated: (generated: boolean) => void; + appendResult: (result: { + iteration: number; + cl: number; + cd: number; + cm: number; + residualL2: number; + }) => void; + setConverged: (converged: boolean) => void; + reset: () => void; +} + +type CfdStore = CfdState & CfdActions; + +const initialState: CfdState = { + config: { ...DEFAULT_CFD_CONFIG }, + isRunning: false, + iteration: 0, + converged: false, + meshGenerated: false, + convergenceHistory: [], + forceHistory: [], + cl: 0, + cd: 0, + cm: 0, + residualL2: 1, +}; + +export const useCfdStore = create((set) => ({ + ...initialState, + + setNi: (ni) => set((s) => ({ config: { ...s.config, ni } })), + setNj: (nj) => set((s) => ({ config: { ...s.config, nj } })), + setMach: (machInf) => set((s) => ({ config: { ...s.config, machInf } })), + setAlphaDeg: (alphaDeg) => set((s) => ({ config: { ...s.config, alphaDeg } })), + setReynolds: (reynolds) => set((s) => ({ config: { ...s.config, reynolds } })), + setCfl: (cfl) => set((s) => ({ config: { ...s.config, cfl } })), + setPhysics: (physics) => set((s) => ({ config: { ...s.config, physics } })), + setReconstruction: (reconstruction) => set((s) => ({ config: { ...s.config, reconstruction } })), + setTimeStepping: (timeStepping) => set((s) => ({ config: { ...s.config, timeStepping } })), + setFarField: (farField) => set((s) => ({ config: { ...s.config, farField } })), + setDs0: (ds0) => set((s) => ({ config: { ...s.config, ds0 } })), + + setRunning: (isRunning) => set({ isRunning }), + setMeshGenerated: (meshGenerated) => set({ meshGenerated }), + setConverged: (converged) => set({ converged }), + + appendResult: (result) => + set((s) => ({ + iteration: result.iteration, + cl: result.cl, + cd: result.cd, + cm: result.cm, + residualL2: result.residualL2, + convergenceHistory: [ + ...s.convergenceHistory, + { iteration: result.iteration, residualL2: result.residualL2 }, + ], + forceHistory: [ + ...s.forceHistory, + { + iteration: result.iteration, + cl: result.cl, + cd: result.cd, + cm: result.cm, + }, + ], + })), + + reset: () => set(initialState), +})); diff --git a/flexfoil-ui/src/types/index.ts b/flexfoil-ui/src/types/index.ts index 1ba33a00..c8ca6ec3 100644 --- a/flexfoil-ui/src/types/index.ts +++ b/flexfoil-ui/src/types/index.ts @@ -18,7 +18,7 @@ export interface AirfoilPoint extends Point { /** Control modes for airfoil manipulation */ export type ControlMode = 'parameters' | 'camber-spline' | 'thickness-spline' | 'inverse-design' | 'geometry-design'; -export type SolverMode = 'viscous' | 'inviscid'; +export type SolverMode = 'viscous' | 'inviscid' | 'cfd'; export type RunMode = 'alpha' | 'cl'; export type AxisVariable = 'alpha' | 'cl' | 'cd' | 'cm' | 'ld' | 'reynolds' | 'mach' | 'ncrit' | 'flapDeflection' | 'flapHingeX'; From 839e9382955f2693b6e26d9d8afdfdc69d0b0020 Mon Sep 17 00:00:00 2001 From: aeronauty Date: Wed, 25 Mar 2026 18:25:37 +0100 Subject: [PATCH 02/12] test: add config tests for rustfoil-cfd (GPU params, serde, sizing) Adds 6 new tests for CfdConfig/CfdParamsGpu: - Default config values - GPU uniform struct is exactly 48 bytes (matches WGSL) - from_config correctly maps physics/reconstruction enums - Byte serialization produces correct little-endian layout - PhysicsMode enum values match WGSL constants - Serde JSON roundtrip Total: 11 tests passing in rustfoil-cfd. CI already runs `cargo test --workspace` which covers this crate. Co-Authored-By: Claude Opus 4.6 (1M context) --- Cargo.lock | 1 + crates/rustfoil-cfd/Cargo.toml | 3 ++ crates/rustfoil-cfd/src/config.rs | 67 +++++++++++++++++++++++++++++++ 3 files changed, 71 insertions(+) diff --git a/Cargo.lock b/Cargo.lock index e81327b6..0c45e9ce 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -611,6 +611,7 @@ version = "0.1.0" dependencies = [ "rustfoil-core", "serde", + "serde_json", ] [[package]] diff --git a/crates/rustfoil-cfd/Cargo.toml b/crates/rustfoil-cfd/Cargo.toml index e30d8aa9..05f8fb14 100644 --- a/crates/rustfoil-cfd/Cargo.toml +++ b/crates/rustfoil-cfd/Cargo.toml @@ -9,3 +9,6 @@ description = "Structured mesh CFD solver support for RustFoil (mesh gen, initia [dependencies] rustfoil-core = { path = "../rustfoil-core" } serde = { version = "1.0", features = ["derive"] } + +[dev-dependencies] +serde_json = "1.0" diff --git a/crates/rustfoil-cfd/src/config.rs b/crates/rustfoil-cfd/src/config.rs index ab9c47f2..600110eb 100644 --- a/crates/rustfoil-cfd/src/config.rs +++ b/crates/rustfoil-cfd/src/config.rs @@ -130,3 +130,70 @@ impl CfdParamsGpu { } } } + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_default_config() { + let cfg = CfdConfig::default(); + assert_eq!(cfg.ni, 256); + assert_eq!(cfg.nj, 128); + assert!((cfg.gamma - 1.4).abs() < 1e-6); + assert!((cfg.prandtl - 0.72).abs() < 1e-6); + } + + #[test] + fn test_gpu_params_size() { + // Must be 48 bytes (12 x 4-byte fields) to match WGSL struct + assert_eq!(std::mem::size_of::(), 48); + } + + #[test] + fn test_gpu_params_from_config() { + let cfg = CfdConfig { + ni: 64, + nj: 32, + mach_inf: 0.5, + alpha: 0.035, // ~2 degrees + physics: PhysicsMode::Euler, + reconstruction: ReconstructionMode::Muscl, + ..CfdConfig::default() + }; + let params = CfdParamsGpu::from_config(&cfg, 0.001, 42); + assert_eq!(params.ni, 64); + assert_eq!(params.nj, 32); + assert_eq!(params.iteration, 42); + assert_eq!(params.physics_mode, 0); // Euler + assert_eq!(params.reconstruction, 0); // MUSCL + assert!((params.dt - 0.001).abs() < 1e-8); + } + + #[test] + fn test_gpu_params_bytes() { + let cfg = CfdConfig::default(); + let params = CfdParamsGpu::from_config(&cfg, 0.01, 0); + let bytes = params.as_bytes(); + assert_eq!(bytes.len(), 48); + // First 4 bytes should be ni=256 as u32 little-endian + let ni = u32::from_le_bytes([bytes[0], bytes[1], bytes[2], bytes[3]]); + assert_eq!(ni, 256); + } + + #[test] + fn test_physics_mode_values() { + assert_eq!(PhysicsMode::Euler as u32, 0); + assert_eq!(PhysicsMode::LaminarNS as u32, 1); + assert_eq!(PhysicsMode::RansSA as u32, 2); + } + + #[test] + fn test_serde_roundtrip() { + let cfg = CfdConfig::default(); + let json = serde_json::to_string(&cfg).unwrap(); + let cfg2: CfdConfig = serde_json::from_str(&json).unwrap(); + assert_eq!(cfg2.ni, cfg.ni); + assert_eq!(cfg2.physics, cfg.physics); + } +} From b871add960504c04ab78e26505bfdae184100199 Mon Sep 17 00:00:00 2001 From: aeronauty Date: Wed, 25 Mar 2026 22:05:25 +0100 Subject: [PATCH 03/12] fix: resolve TypeScript errors and wire CFD panel into UI - Add CFD to PANELS array (layoutConfig.ts) and factory (DockingLayout.tsx) - Add CFD to MobileLayout tabs and component map - Fix CfdSolver.ts: use type-only imports for verbatimModuleSyntax compliance - Fix CfdSolver.ts: cast Float32Array/Uint32Array for GPUAllowSharedBufferSource - Fix CfdPanel.tsx: use static imports from rustfoil-wasm, flatten coordinates to Float64Array, remove unused forceHistory destructure - Fix AirfoilCanvas.tsx: guard SolverMode='cfd' when calling panel method functions (computeStreamlines, computePsiGrid, computeDividingStreamline) - Add cfd_* function declarations to rustfoil-wasm.d.ts TypeScript compiles cleanly (tsc -b --noEmit passes with 0 errors). Co-Authored-By: Claude Opus 4.6 (1M context) --- flexfoil-ui/src/components/AirfoilCanvas.tsx | 9 ++++++--- flexfoil-ui/src/components/DockingLayout.tsx | 7 +++++++ flexfoil-ui/src/components/MobileLayout.tsx | 3 +++ .../src/components/panels/CfdPanel.tsx | 18 +++++++++++------- flexfoil-ui/src/layoutConfig.ts | 1 + flexfoil-ui/src/lib/webgpu/cfd/CfdSolver.ts | 19 ++++++++++--------- flexfoil-ui/src/rustfoil-wasm.d.ts | 3 +++ 7 files changed, 41 insertions(+), 19 deletions(-) diff --git a/flexfoil-ui/src/components/AirfoilCanvas.tsx b/flexfoil-ui/src/components/AirfoilCanvas.tsx index aaa585d0..c233be46 100644 --- a/flexfoil-ui/src/components/AirfoilCanvas.tsx +++ b/flexfoil-ui/src/components/AirfoilCanvas.tsx @@ -945,7 +945,8 @@ export function AirfoilCanvas() { } try { - const result = computeStreamlines(panels, displayAlpha, reynolds, adaptiveStreamlineCount, streamlineBounds, mach, ncrit, maxIterations, solverMode); + const panelSolverMode = solverMode === 'cfd' ? undefined : solverMode; + const result = computeStreamlines(panels, displayAlpha, reynolds, adaptiveStreamlineCount, streamlineBounds, mach, ncrit, maxIterations, panelSolverMode); if (result.success) { setStreamlines(result.streamlines); } @@ -976,7 +977,8 @@ export function AirfoilCanvas() { Math.min(200, Math.round(xRange * 40)), Math.min(120, Math.round(yRange * 40)) ]; - const result = computePsiGrid(panels, displayAlpha, reynolds, bounds, resolution, mach, ncrit, maxIterations, solverMode); + const panelSolverMode = solverMode === 'cfd' ? undefined : solverMode; + const result = computePsiGrid(panels, displayAlpha, reynolds, bounds, resolution, mach, ncrit, maxIterations, panelSolverMode); if (!result.success) { console.error('Psi grid computation failed:', result.error); @@ -1029,6 +1031,7 @@ export function AirfoilCanvas() { } let psi0Lines: [number, number][][] = []; + const dividingPanelMode = solverMode === 'cfd' ? undefined : solverMode; const dividingResult = computeDividingStreamline( panels, displayAlpha, @@ -1037,7 +1040,7 @@ export function AirfoilCanvas() { mach, ncrit, maxIterations, - solverMode + dividingPanelMode ); if (dividingResult.success && dividingResult.streamline.length >= 2) { diff --git a/flexfoil-ui/src/components/DockingLayout.tsx b/flexfoil-ui/src/components/DockingLayout.tsx index 84f7cde1..9c20c59b 100644 --- a/flexfoil-ui/src/components/DockingLayout.tsx +++ b/flexfoil-ui/src/components/DockingLayout.tsx @@ -18,6 +18,7 @@ import { DataExplorerPanel } from './panels/DataExplorerPanel'; import { PlotBuilderPanel } from './panels/PlotBuilderPanel'; import { CaseLogsPanel } from './panels/CaseLogsPanel'; import { DistributionPanel } from './panels/DistributionPanel'; +import { CfdPanel } from './panels/CfdPanel'; import { MenuBar } from './MenuBar'; import { CommandPalette } from './CommandPalette'; import { FeedbackWidget } from './FeedbackWidget'; @@ -390,6 +391,12 @@ function DesktopLayout({ wasmStatus }: DockingLayoutProps) { ); + case 'cfd': + return ( +
+ +
+ ); default: return
Unknown panel: {component}
; } diff --git a/flexfoil-ui/src/components/MobileLayout.tsx b/flexfoil-ui/src/components/MobileLayout.tsx index bb72e4d9..b6eac9df 100644 --- a/flexfoil-ui/src/components/MobileLayout.tsx +++ b/flexfoil-ui/src/components/MobileLayout.tsx @@ -11,6 +11,7 @@ import { DataExplorerPanel } from './panels/DataExplorerPanel'; import { PlotBuilderPanel } from './panels/PlotBuilderPanel'; import { CaseLogsPanel } from './panels/CaseLogsPanel'; import { DistributionPanel } from './panels/DistributionPanel'; +import { CfdPanel } from './panels/CfdPanel'; import { DarkModeToggle } from './DarkModeToggle'; import { LayoutContext } from '../contexts/LayoutContext'; import { useOnboarding } from '../onboarding'; @@ -30,6 +31,7 @@ const TABS: { id: PanelId; label: string }[] = [ { id: 'plot-builder', label: 'Plots' }, { id: 'case-logs', label: 'Logs' }, { id: 'distributions', label: 'Distrib.' }, + { id: 'cfd', label: 'CFD' }, ]; const PANEL_COMPONENTS: Record = { @@ -45,6 +47,7 @@ const PANEL_COMPONENTS: Record = { 'plot-builder': PlotBuilderPanel, 'case-logs': CaseLogsPanel, 'distributions': DistributionPanel, + 'cfd': CfdPanel, }; interface MobileLayoutProps { diff --git a/flexfoil-ui/src/components/panels/CfdPanel.tsx b/flexfoil-ui/src/components/panels/CfdPanel.tsx index 1c0ba62c..346b2c55 100644 --- a/flexfoil-ui/src/components/panels/CfdPanel.tsx +++ b/flexfoil-ui/src/components/panels/CfdPanel.tsx @@ -11,6 +11,7 @@ import { useAirfoilStore } from '../../stores/airfoilStore'; import { CfdSolver } from '../../lib/webgpu/cfd/CfdSolver'; import type { CfdSolverResult } from '../../lib/webgpu/cfd/CfdSolver'; import type { PhysicsMode } from '../../lib/webgpu/cfd/CfdConfig'; +import { cfd_generate_mesh, cfd_initial_conditions, cfd_boundary_types } from 'rustfoil-wasm'; export function CfdPanel() { const { @@ -18,7 +19,6 @@ export function CfdPanel() { isRunning, iteration, convergenceHistory, - forceHistory, cl, cd, cm, @@ -56,10 +56,14 @@ export function CfdPanel() { try { setError(null); - // Import WASM functions dynamically - const wasm = await import('rustfoil-wasm'); - const meshResult = wasm.cfd_generate_mesh( - coordinates, + // Flatten [x,y] pairs into a flat Float64Array for WASM + const flat = new Float64Array(coordinates.length * 2); + for (let i = 0; i < coordinates.length; i++) { + flat[i * 2] = coordinates[i].x; + flat[i * 2 + 1] = coordinates[i].y; + } + const meshResult = cfd_generate_mesh( + flat, config.ni, config.nj, config.farField, @@ -73,7 +77,7 @@ export function CfdPanel() { // Get initial conditions const physicsMode = config.physics === 'euler' ? 0 : config.physics === 'laminar_ns' ? 1 : 2; - const initialQ = wasm.cfd_initial_conditions( + const initialQ = cfd_initial_conditions( config.ni, config.nj, config.machInf, @@ -83,7 +87,7 @@ export function CfdPanel() { config.reynolds, ); - const bcTypes = wasm.cfd_boundary_types(config.ni, config.nj); + const bcTypes = cfd_boundary_types(config.ni, config.nj); // Check for WebGPU if (!navigator.gpu) { diff --git a/flexfoil-ui/src/layoutConfig.ts b/flexfoil-ui/src/layoutConfig.ts index ef702e57..ba2863ce 100644 --- a/flexfoil-ui/src/layoutConfig.ts +++ b/flexfoil-ui/src/layoutConfig.ts @@ -13,6 +13,7 @@ export const PANELS = [ { id: 'plot-builder', name: 'Plot Builder', component: 'plot-builder' }, { id: 'case-logs', name: 'Case Logs', component: 'case-logs' }, { id: 'distributions', name: 'Distributions', component: 'distributions' }, + { id: 'cfd', name: 'CFD', component: 'cfd' }, ] as const; export type PanelId = (typeof PANELS)[number]['id']; diff --git a/flexfoil-ui/src/lib/webgpu/cfd/CfdSolver.ts b/flexfoil-ui/src/lib/webgpu/cfd/CfdSolver.ts index 483f07c1..faf25a8f 100644 --- a/flexfoil-ui/src/lib/webgpu/cfd/CfdSolver.ts +++ b/flexfoil-ui/src/lib/webgpu/cfd/CfdSolver.ts @@ -5,11 +5,12 @@ * Dispatches per-timestep compute passes and handles async readback. */ -import { CfdConfig, packCfdParams } from './CfdConfig'; -import { CfdBufferSet, createCfdBuffers, destroyCfdBuffers } from './CfdBuffers'; +import type { CfdConfig } from './CfdConfig'; +import { packCfdParams } from './CfdConfig'; +import type { CfdBufferSet } from './CfdBuffers'; +import { createCfdBuffers, destroyCfdBuffers } from './CfdBuffers'; +import type { CfdPipelineSet, CfdBindGroupSet } from './CfdPipelines'; import { - CfdPipelineSet, - CfdBindGroupSet, createCfdPipelines, createCfdBindGroups, } from './CfdPipelines'; @@ -68,10 +69,10 @@ export class CfdSolver { const bindGroups = createCfdBindGroups(device, pipelines, buffers); // Upload mesh data - device.queue.writeBuffer(buffers.meshX, 0, meshX); - device.queue.writeBuffer(buffers.meshY, 0, meshY); - device.queue.writeBuffer(buffers.Q, 0, initialQ); - device.queue.writeBuffer(buffers.bcType, 0, bcTypes); + device.queue.writeBuffer(buffers.meshX, 0, meshX as Float32Array); + device.queue.writeBuffer(buffers.meshY, 0, meshY as Float32Array); + device.queue.writeBuffer(buffers.Q, 0, initialQ as Float32Array); + device.queue.writeBuffer(buffers.bcType, 0, bcTypes as Uint32Array); // Upload initial params const dt = config.cfl * 0.01 / (1.0 + config.machInf); @@ -246,7 +247,7 @@ export class CfdSolver { /** Reset solver to initial conditions. */ reset(initialQ: Float32Array): void { this.iteration = 0; - this.device.queue.writeBuffer(this.buffers.Q, 0, initialQ); + this.device.queue.writeBuffer(this.buffers.Q, 0, initialQ as Float32Array); } /** Clean up all GPU resources. */ diff --git a/flexfoil-ui/src/rustfoil-wasm.d.ts b/flexfoil-ui/src/rustfoil-wasm.d.ts index 5955d48b..1b946323 100644 --- a/flexfoil-ui/src/rustfoil-wasm.d.ts +++ b/flexfoil-ui/src/rustfoil-wasm.d.ts @@ -39,4 +39,7 @@ declare module 'rustfoil-wasm' { export const gdes_scale_thickness: (...args: any[]) => any; export const gdes_scale_camber: (...args: any[]) => any; export const full_inverse_design_mdes: (...args: any[]) => any; + export const cfd_generate_mesh: (coords_flat: Float64Array, ni: number, nj: number, far_field: number, ds0: number) => any; + export const cfd_initial_conditions: (ni: number, nj: number, mach: number, alpha_deg: number, gamma: number, physics_mode: number, reynolds: number) => Float32Array; + export const cfd_boundary_types: (ni: number, nj: number) => Uint32Array; } From 31490e46b0e7abc035f03e7546cd91ecac7249f8 Mon Sep 17 00:00:00 2001 From: aeronauty Date: Wed, 25 Mar 2026 22:18:17 +0100 Subject: [PATCH 04/12] fix: resolve WebGPU bind group layout mismatches and buffer sizing Root cause: layout 'auto' generates bind group layouts only for bindings actually referenced by each entry point, not all declared in the shader module. Fixes: - Split shared reconstruct/roeFlux bind groups into per-direction variants (reconstructXi/Eta, roeFluxXi/Eta) matching what each entry point actually uses - Add dummy reads in bc.wgsl for bc_type, mesh_x, mesh_y so the auto-layout includes all 6 declared bindings - Fix forceAccum buffer from 12 to 16 bytes (WebGPU minimum) - Fix readback buffer from 28 to 32 bytes to match - Fix copyBufferToBuffer sizes to copy full 16-byte aligned chunks - Restructure step() to one command encoder per timestep Co-Authored-By: Claude Opus 4.6 (1M context) --- flexfoil-ui/src/lib/webgpu/cfd/CfdBuffers.ts | 4 +-- .../src/lib/webgpu/cfd/CfdPipelines.ts | 34 +++++++++++++++---- flexfoil-ui/src/lib/webgpu/cfd/CfdSolver.ts | 33 ++++++++++-------- .../src/lib/webgpu/shaders/cfd/bc.wgsl | 9 +++++ 4 files changed, 57 insertions(+), 23 deletions(-) diff --git a/flexfoil-ui/src/lib/webgpu/cfd/CfdBuffers.ts b/flexfoil-ui/src/lib/webgpu/cfd/CfdBuffers.ts index f37aa6b8..346b59ce 100644 --- a/flexfoil-ui/src/lib/webgpu/cfd/CfdBuffers.ts +++ b/flexfoil-ui/src/lib/webgpu/cfd/CfdBuffers.ts @@ -72,7 +72,7 @@ export function createCfdBuffers(device: GPUDevice, ni: number, nj: number): Cfd }), forceAccum: device.createBuffer({ label: 'CFD force_accum', - size: 12, // 3 x f32 (Cl, Cd, Cm) + size: 16, // 4 x f32 (Cl, Cd, Cm, _pad) — must be >= 16 for storage usage: STORAGE_RW | GPUBufferUsage.COPY_SRC, }), residualNorm: device.createBuffer({ @@ -82,7 +82,7 @@ export function createCfdBuffers(device: GPUDevice, ni: number, nj: number): Cfd }), readback: device.createBuffer({ label: 'CFD readback', - size: 28, // 7 x f32 (3 forces + 4 residual norms) + size: 32, // 8 x f32 (4 force_accum + 4 residual_norm) usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST, }), }; diff --git a/flexfoil-ui/src/lib/webgpu/cfd/CfdPipelines.ts b/flexfoil-ui/src/lib/webgpu/cfd/CfdPipelines.ts index ad54642c..e9a3337d 100644 --- a/flexfoil-ui/src/lib/webgpu/cfd/CfdPipelines.ts +++ b/flexfoil-ui/src/lib/webgpu/cfd/CfdPipelines.ts @@ -31,8 +31,10 @@ export interface CfdPipelineSet { export interface CfdBindGroupSet { metrics: GPUBindGroup; primitives: GPUBindGroup; - reconstruct: GPUBindGroup; - roeFlux: GPUBindGroup; + reconstructXi: GPUBindGroup; + reconstructEta: GPUBindGroup; + roeFluxXi: GPUBindGroup; + roeFluxEta: GPUBindGroup; residual: GPUBindGroup; update: GPUBindGroup; bc: GPUBindGroup; @@ -137,29 +139,47 @@ export function createCfdBindGroups( { binding: 2, resource: { buffer: buffers.W } }, ], }), - reconstruct: device.createBindGroup({ - label: 'cfd_reconstruct_bg', + // Separate bind groups for xi/eta because auto-layout only includes + // bindings actually referenced by each entry point. + reconstructXi: device.createBindGroup({ + label: 'cfd_reconstruct_xi_bg', layout: pipelines.reconstructXi.getBindGroupLayout(0), entries: [ { binding: 0, resource: { buffer: buffers.params } }, { binding: 1, resource: { buffer: buffers.W } }, { binding: 2, resource: { buffer: buffers.leftXi } }, { binding: 3, resource: { buffer: buffers.rightXi } }, + ], + }), + reconstructEta: device.createBindGroup({ + label: 'cfd_reconstruct_eta_bg', + layout: pipelines.reconstructEta.getBindGroupLayout(0), + entries: [ + { binding: 0, resource: { buffer: buffers.params } }, + { binding: 1, resource: { buffer: buffers.W } }, { binding: 4, resource: { buffer: buffers.leftEta } }, { binding: 5, resource: { buffer: buffers.rightEta } }, ], }), - roeFlux: device.createBindGroup({ - label: 'cfd_roe_flux_bg', + roeFluxXi: device.createBindGroup({ + label: 'cfd_roe_flux_xi_bg', layout: pipelines.roeFluxXi.getBindGroupLayout(0), entries: [ { binding: 0, resource: { buffer: buffers.params } }, { binding: 1, resource: { buffer: buffers.metrics } }, { binding: 2, resource: { buffer: buffers.leftXi } }, { binding: 3, resource: { buffer: buffers.rightXi } }, + { binding: 6, resource: { buffer: buffers.fluxXi } }, + ], + }), + roeFluxEta: device.createBindGroup({ + label: 'cfd_roe_flux_eta_bg', + layout: pipelines.roeFluxEta.getBindGroupLayout(0), + entries: [ + { binding: 0, resource: { buffer: buffers.params } }, + { binding: 1, resource: { buffer: buffers.metrics } }, { binding: 4, resource: { buffer: buffers.leftEta } }, { binding: 5, resource: { buffer: buffers.rightEta } }, - { binding: 6, resource: { buffer: buffers.fluxXi } }, { binding: 7, resource: { buffer: buffers.fluxEta } }, ], }), diff --git a/flexfoil-ui/src/lib/webgpu/cfd/CfdSolver.ts b/flexfoil-ui/src/lib/webgpu/cfd/CfdSolver.ts index faf25a8f..e6b4621c 100644 --- a/flexfoil-ui/src/lib/webgpu/cfd/CfdSolver.ts +++ b/flexfoil-ui/src/lib/webgpu/cfd/CfdSolver.ts @@ -113,15 +113,16 @@ export class CfdSolver { const wgXi = Math.ceil(ni / 16); const wgEta = Math.ceil(nj / 16); - const encoder = this.device.createCommandEncoder({ label: `cfd_steps_${this.iteration}` }); - for (let s = 0; s < nSteps; s++) { this.iteration++; - // Update params uniform + // Update params uniform before encoding const paramsData = packCfdParams(this.config, dt, this.iteration); this.device.queue.writeBuffer(this.buffers.params, 0, paramsData); + // One encoder per step to avoid writeBuffer/encoder interleaving issues + const encoder = this.device.createCommandEncoder({ label: `cfd_step_${this.iteration}` }); + // Pass 1: Q -> W (conservative to primitive) { const pass = encoder.beginComputePass({ label: 'primitives' }); @@ -135,7 +136,7 @@ export class CfdSolver { { const pass = encoder.beginComputePass({ label: 'reconstruct_xi' }); pass.setPipeline(this.pipelines.reconstructXi); - pass.setBindGroup(0, this.bindGroups.reconstruct); + pass.setBindGroup(0, this.bindGroups.reconstructXi); pass.dispatchWorkgroups(nj); pass.end(); } @@ -144,7 +145,7 @@ export class CfdSolver { { const pass = encoder.beginComputePass({ label: 'roe_flux_xi' }); pass.setPipeline(this.pipelines.roeFluxXi); - pass.setBindGroup(0, this.bindGroups.roeFlux); + pass.setBindGroup(0, this.bindGroups.roeFluxXi); pass.dispatchWorkgroups(wgXi, wgEta); pass.end(); } @@ -153,7 +154,7 @@ export class CfdSolver { { const pass = encoder.beginComputePass({ label: 'reconstruct_eta' }); pass.setPipeline(this.pipelines.reconstructEta); - pass.setBindGroup(0, this.bindGroups.reconstruct); + pass.setBindGroup(0, this.bindGroups.reconstructEta); pass.dispatchWorkgroups(ni); pass.end(); } @@ -162,7 +163,7 @@ export class CfdSolver { { const pass = encoder.beginComputePass({ label: 'roe_flux_eta' }); pass.setPipeline(this.pipelines.roeFluxEta); - pass.setBindGroup(0, this.bindGroups.roeFlux); + pass.setBindGroup(0, this.bindGroups.roeFluxEta); pass.dispatchWorkgroups(wgXi, wgEta); pass.end(); } @@ -193,22 +194,26 @@ export class CfdSolver { pass.dispatchWorkgroups(Math.ceil(ni / 256)); pass.end(); } + + this.device.queue.submit([encoder.finish()]); } - // Final pass: compute forces and residual + // Final pass: compute forces and residual norm { + const encoder = this.device.createCommandEncoder({ label: 'cfd_forces' }); + const pass = encoder.beginComputePass({ label: 'forces' }); pass.setPipeline(this.pipelines.forcesReduce); pass.setBindGroup(0, this.bindGroups.forcesReduce); pass.dispatchWorkgroups(1); // Single workgroup reduction pass.end(); - } - // Copy forces and residual to readback buffer - encoder.copyBufferToBuffer(this.buffers.forceAccum, 0, this.buffers.readback, 0, 12); - encoder.copyBufferToBuffer(this.buffers.residualNorm, 0, this.buffers.readback, 12, 4); + // Copy forces and residual to readback buffer + encoder.copyBufferToBuffer(this.buffers.forceAccum, 0, this.buffers.readback, 0, 16); + encoder.copyBufferToBuffer(this.buffers.residualNorm, 0, this.buffers.readback, 16, 16); - this.device.queue.submit([encoder.finish()]); + this.device.queue.submit([encoder.finish()]); + } // Async readback await this.buffers.readback.mapAsync(GPUMapMode.READ); @@ -219,7 +224,7 @@ export class CfdSolver { cl: data[0], cd: data[1], cm: data[2], - residualL2: data[3], + residualL2: data[4], // offset 4: first f32 of residualNorm (at readback byte 16) iteration: this.iteration, }; } diff --git a/flexfoil-ui/src/lib/webgpu/shaders/cfd/bc.wgsl b/flexfoil-ui/src/lib/webgpu/shaders/cfd/bc.wgsl index 2bd7b328..02186adb 100644 --- a/flexfoil-ui/src/lib/webgpu/shaders/cfd/bc.wgsl +++ b/flexfoil-ui/src/lib/webgpu/shaders/cfd/bc.wgsl @@ -24,6 +24,15 @@ fn apply_boundary_conditions(@builtin(global_invocation_id) gid: vec3) { let gamma = params.gamma; + // Touch all bindings so the auto-layout includes them. + // The compiler cannot eliminate these since they are storage reads. + let _bc_touch = bc_type[0]; + let _mx_touch = mesh_x[0]; + let _my_touch = mesh_y[0]; + _ = _bc_touch; + _ = _mx_touch; + _ = _my_touch; + // === Wall boundary (j=0) === { let j = 0u; From f92f51a4c806a5284940a298973115b901e10f9e Mon Sep 17 00:00:00 2001 From: aeronauty Date: Wed, 25 Mar 2026 22:23:29 +0100 Subject: [PATCH 05/12] fix: add dummy metrics read in forces_reduce.wgsl for auto-layout The reduce_wall_forces entry point declared metrics at binding 2 but never read it, causing layout:auto to exclude binding 2. This caused CreateBindGroup to fail with binding index 2 not present in layout. Co-Authored-By: Claude Opus 4.6 (1M context) --- flexfoil-ui/src/lib/webgpu/shaders/cfd/forces_reduce.wgsl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/flexfoil-ui/src/lib/webgpu/shaders/cfd/forces_reduce.wgsl b/flexfoil-ui/src/lib/webgpu/shaders/cfd/forces_reduce.wgsl index f0cefc09..80760e41 100644 --- a/flexfoil-ui/src/lib/webgpu/shaders/cfd/forces_reduce.wgsl +++ b/flexfoil-ui/src/lib/webgpu/shaders/cfd/forces_reduce.wgsl @@ -29,6 +29,10 @@ fn reduce_wall_forces(@builtin(local_invocation_id) lid: vec3, let mach = params.mach_inf; let alpha = params.alpha; + // Touch all bindings so auto-layout includes them + let _m_touch = metrics[0]; + _ = _m_touch; + // Reference quantities for force coefficients let p_inf = 1.0 / gamma; let q_inf = 0.5 * gamma * p_inf * mach * mach; // Dynamic pressure From 23aebeded4073f3ad9311fb0ef3ad583a010bd14 Mon Sep 17 00:00:00 2001 From: aeronauty Date: Thu, 26 Mar 2026 07:48:23 +0100 Subject: [PATCH 06/12] test: add integration tests for full CFD pipeline Tests the complete mesh generation pipeline with NACA 0012: - Full pipeline: mesh gen, initial conditions, BC types, GPU params - Mesh quality: cell areas, no degenerate cells - RANS: SA turbulence initialization All pass, confirming Rust CFD code works correctly. Co-Authored-By: Claude Opus 4.6 (1M context) --- crates/rustfoil-cfd/tests/integration.rs | 270 +++++++++++++++++++++++ 1 file changed, 270 insertions(+) create mode 100644 crates/rustfoil-cfd/tests/integration.rs diff --git a/crates/rustfoil-cfd/tests/integration.rs b/crates/rustfoil-cfd/tests/integration.rs new file mode 100644 index 00000000..4258276b --- /dev/null +++ b/crates/rustfoil-cfd/tests/integration.rs @@ -0,0 +1,270 @@ +//! Integration tests for the full CFD mesh generation pipeline. +//! +//! Tests the complete workflow: NACA 0012 coordinates → O-mesh → initial conditions → BC types. + +use std::f64::consts::PI; + +use rustfoil_cfd::boundary::generate_bc_types; +use rustfoil_cfd::config::{CfdConfig, CfdParamsGpu, PhysicsMode}; +use rustfoil_cfd::init::compute_initial_conditions; +use rustfoil_cfd::mesh::generate_o_mesh; + +/// Generate NACA 0012 airfoil coordinates (closed trailing edge). +fn naca0012(n: usize) -> (Vec, Vec) { + let mut x = Vec::with_capacity(n); + let mut y = Vec::with_capacity(n); + + // Upper surface from TE to LE (i=0..n/2) + let half = n / 2; + for i in 0..half { + let theta = PI * (i as f64) / (half as f64); + let xc = 0.5 * (1.0 + theta.cos()); // cosine spacing: 1.0 → 0.0 + let t = 0.12; // max thickness ratio + let yt = 5.0 + * t + * (0.2969 * xc.sqrt() + - 0.1260 * xc + - 0.3516 * xc.powi(2) + + 0.2843 * xc.powi(3) + - 0.1015 * xc.powi(4)); + x.push(xc); + y.push(yt); + } + // Lower surface from LE to TE (i=n/2..n) + for i in 0..half { + let theta = PI * (i as f64) / (half as f64); + let xc = 0.5 * (1.0 - theta.cos()); // cosine spacing: 0.0 → 1.0 + let t = 0.12; + let yt = 5.0 + * t + * (0.2969 * xc.sqrt() + - 0.1260 * xc + - 0.3516 * xc.powi(2) + + 0.2843 * xc.powi(3) + - 0.1015 * xc.powi(4)); + x.push(xc); + y.push(-yt); + } + + (x, y) +} + +#[test] +fn test_full_pipeline_naca0012() { + // Step 1: Generate airfoil coordinates + let (ax, ay) = naca0012(128); + assert_eq!(ax.len(), 128); + assert_eq!(ay.len(), 128); + + // Verify airfoil bounds + let x_min = ax.iter().cloned().fold(f64::INFINITY, f64::min); + let x_max = ax.iter().cloned().fold(f64::NEG_INFINITY, f64::max); + assert!(x_min >= -0.01, "x_min should be near 0, got {x_min}"); + assert!(x_max <= 1.01, "x_max should be near 1, got {x_max}"); + + // Step 2: Generate O-mesh + let ni = 128u32; + let nj = 64u32; + let far_field = 15.0; + let ds0 = 0.001; + + let mesh = generate_o_mesh(&ax, &ay, ni, nj, far_field, ds0); + + assert_eq!(mesh.ni, ni); + assert_eq!(mesh.nj, nj); + assert_eq!(mesh.x.len(), (ni * nj) as usize); + assert_eq!(mesh.y.len(), (ni * nj) as usize); + + // Verify wall (j=0) points are near the airfoil + for i in 0..ni as usize { + let wx = mesh.x[i] as f64; + let wy = mesh.y[i] as f64; + // Wall points should be within the chord length range + assert!( + wx > -5.0 && wx < 5.0, + "Wall point x[{i}] = {wx} out of range" + ); + assert!( + wy > -5.0 && wy < 5.0, + "Wall point y[{i}] = {wy} out of range" + ); + } + + // Verify far-field (j=nj-1) points are far away + let jj = (nj - 1) as usize; + let mut max_dist = 0.0f64; + for i in 0..ni as usize { + let fx = mesh.x[jj * ni as usize + i] as f64; + let fy = mesh.y[jj * ni as usize + i] as f64; + let dist = (fx * fx + fy * fy).sqrt(); + if dist > max_dist { + max_dist = dist; + } + } + assert!( + max_dist > 5.0, + "Far-field points should be far from origin, max_dist={max_dist}" + ); + + // No NaN or Inf in mesh + for val in &mesh.x { + assert!(val.is_finite(), "mesh.x contains non-finite value"); + } + for val in &mesh.y { + assert!(val.is_finite(), "mesh.y contains non-finite value"); + } + + println!("Mesh generated: ni={}, nj={}", mesh.ni, mesh.nj); + println!(" Wall x range: [{:.4}, {:.4}]", + mesh.x[..ni as usize].iter().cloned().fold(f32::INFINITY, f32::min), + mesh.x[..ni as usize].iter().cloned().fold(f32::NEG_INFINITY, f32::max), + ); + println!(" Far-field max distance: {:.2}", max_dist); + + // Step 3: Compute initial conditions + let config = CfdConfig { + ni, + nj, + mach_inf: 0.5, + alpha: 2.0f32.to_radians(), + gamma: 1.4, + reynolds: 1e5, + physics: PhysicsMode::Euler, + ..CfdConfig::default() + }; + + let q = compute_initial_conditions(&config); + assert_eq!(q.len(), (ni * nj * 5) as usize); + + // Verify freestream state at an arbitrary cell + let cell = 100; + let rho = q[cell * 5]; + let rhou = q[cell * 5 + 1]; + let rhov = q[cell * 5 + 2]; + let e = q[cell * 5 + 3]; + let nu = q[cell * 5 + 4]; + + assert!((rho - 1.0).abs() < 1e-5, "rho = {rho}"); + assert!(rhou > 0.0, "rhou should be positive for alpha~2deg"); + assert!(rhov > 0.0, "rhov should be positive for alpha~2deg"); + assert!(e > 0.0, "energy should be positive"); + assert!(nu.abs() < 1e-10, "nu_tilde should be 0 for Euler"); + + // No NaN or Inf in Q + for val in &q { + assert!(val.is_finite(), "Q contains non-finite value"); + } + + println!(" Initial conditions: rho={rho:.4}, rhou={rhou:.4}, rhov={rhov:.4}, E={e:.4}"); + + // Step 4: Generate boundary condition types + let bc = generate_bc_types(ni, nj); + assert_eq!(bc.len(), (ni * nj) as usize); + + // Wall at j=0 + for i in 0..ni as usize { + assert_eq!(bc[i], 1, "bc[j=0, i={i}] should be Wall(1)"); + } + // Interior at j=1 + for i in 0..ni as usize { + assert_eq!(bc[ni as usize + i], 0, "bc[j=1, i={i}] should be Interior(0)"); + } + // FarField at j=nj-1 + for i in 0..ni as usize { + assert_eq!( + bc[(nj as usize - 1) * ni as usize + i], + 2, + "bc[j=nj-1, i={i}] should be FarField(2)" + ); + } + + println!(" BC types: wall={}, interior={}, farfield={}", + bc.iter().filter(|&&v| v == 1).count(), + bc.iter().filter(|&&v| v == 0).count(), + bc.iter().filter(|&&v| v == 2).count(), + ); + + // Step 5: Verify GPU params struct + let params = CfdParamsGpu::from_config(&config, 0.001, 0); + assert_eq!(params.ni, ni); + assert_eq!(params.nj, nj); + assert!((params.gamma - 1.4).abs() < 1e-6); + assert!((params.mach_inf - 0.5).abs() < 1e-6); + assert_eq!(params.physics_mode, 0); // Euler + let bytes = params.as_bytes(); + assert_eq!(bytes.len(), 48); + + println!(" GPU params: {} bytes", bytes.len()); + println!("\nFull CFD pipeline test PASSED ✓"); +} + +#[test] +fn test_mesh_quality_metrics() { + // Test that the mesh has reasonable aspect ratios and no degenerate cells + let (ax, ay) = naca0012(64); + let mesh = generate_o_mesh(&ax, &ay, 64, 32, 10.0, 0.005); + + let ni = 64usize; + let nj = 32usize; + let mut min_area = f64::INFINITY; + let mut max_area = 0.0f64; + let mut degenerate_count = 0; + + // Check cell areas (should all be positive for a valid mesh) + for j in 0..nj - 1 { + for i in 0..ni { + let ip = (i + 1) % ni; + + let x00 = mesh.x[j * ni + i] as f64; + let y00 = mesh.y[j * ni + i] as f64; + let x10 = mesh.x[j * ni + ip] as f64; + let y10 = mesh.y[j * ni + ip] as f64; + let x01 = mesh.x[(j + 1) * ni + i] as f64; + let y01 = mesh.y[(j + 1) * ni + i] as f64; + let x11 = mesh.x[(j + 1) * ni + ip] as f64; + let y11 = mesh.y[(j + 1) * ni + ip] as f64; + + // Cross product of diagonals gives 2x area + let area = 0.5 + * ((x11 - x00) * (y01 - y10) - (x01 - x10) * (y11 - y00)) + .abs(); + + if area < 1e-12 { + degenerate_count += 1; + } + if area < min_area { + min_area = area; + } + if area > max_area { + max_area = area; + } + } + } + + println!("Mesh quality: min_area={:.2e}, max_area={:.2e}, degenerate={}", min_area, max_area, degenerate_count); + assert!( + degenerate_count < (ni * nj / 10), + "Too many degenerate cells: {degenerate_count}" + ); + assert!(min_area > 0.0, "All cell areas should be positive"); +} + +#[test] +fn test_rans_initial_conditions() { + let config = CfdConfig { + ni: 32, + nj: 16, + mach_inf: 0.3, + alpha: 0.0, + reynolds: 1e6, + physics: PhysicsMode::RansSA, + ..CfdConfig::default() + }; + + let q = compute_initial_conditions(&config); + + // RANS should have non-zero nu_tilde + let nu = q[4]; // First cell, 5th variable + assert!(nu > 0.0, "RANS nu_tilde should be positive, got {nu}"); + println!("RANS nu_tilde_inf = {:.6e}", nu); +} From 17ea32bf68cfa44f5100ed3895a94e94ef51dae5 Mon Sep 17 00:00:00 2001 From: aeronauty Date: Thu, 26 Mar 2026 13:04:58 +0100 Subject: [PATCH 07/12] feat: add CFD mesh Python bindings, API, and /cfd visualization page Co-Authored-By: Claude Opus 4.6 (1M context) --- Cargo.lock | 1 + crates/rustfoil-python/Cargo.toml | 1 + crates/rustfoil-python/src/lib.rs | 91 +++++ .../flexfoil-python/src/flexfoil/server.py | 326 ++++++++++++++++++ 4 files changed, 419 insertions(+) diff --git a/Cargo.lock b/Cargo.lock index 0c45e9ce..6bdbd44a 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -672,6 +672,7 @@ version = "0.1.0" dependencies = [ "pyo3", "rayon", + "rustfoil-cfd", "rustfoil-core", "rustfoil-inviscid", "rustfoil-solver", diff --git a/crates/rustfoil-python/Cargo.toml b/crates/rustfoil-python/Cargo.toml index ddf88a70..2748ba9c 100644 --- a/crates/rustfoil-python/Cargo.toml +++ b/crates/rustfoil-python/Cargo.toml @@ -15,5 +15,6 @@ rustfoil-core = { path = "../rustfoil-core" } rustfoil-solver = { path = "../rustfoil-solver" } rustfoil-inviscid = { path = "../rustfoil-inviscid" } rustfoil-xfoil = { path = "../rustfoil-xfoil" } +rustfoil-cfd = { path = "../rustfoil-cfd" } pyo3 = { version = "0.23", features = ["extension-module"] } rayon = "1.10" diff --git a/crates/rustfoil-python/src/lib.rs b/crates/rustfoil-python/src/lib.rs index e0903f41..1d2b98a2 100644 --- a/crates/rustfoil-python/src/lib.rs +++ b/crates/rustfoil-python/src/lib.rs @@ -504,6 +504,94 @@ fn get_bl_distribution( Ok(d.into()) } +// ─── CFD mesh generation ─────────────────────────────────────────────── + +/// Generate an O-type structured mesh around an airfoil. +/// +/// Args: +/// coords_flat: Flat [x0,y0,x1,y1,...] airfoil coordinates +/// ni: Circumferential grid points +/// nj: Radial grid layers +/// far_field: Far-field distance in chord lengths +/// ds0: First cell wall-normal spacing +/// +/// Returns: +/// dict with keys: x (list[float]), y (list[float]), ni (int), nj (int) +#[pyfunction] +fn cfd_generate_mesh<'py>( + py: Python<'py>, + coords_flat: Vec, + ni: u32, + nj: u32, + far_field: f64, + ds0: f64, +) -> PyResult> { + let n_pts = coords_flat.len() / 2; + let mut ax = Vec::with_capacity(n_pts); + let mut ay = Vec::with_capacity(n_pts); + for i in 0..n_pts { + ax.push(coords_flat[i * 2]); + ay.push(coords_flat[i * 2 + 1]); + } + + let mesh = rustfoil_cfd::mesh::generate_o_mesh(&ax, &ay, ni, nj, far_field, ds0); + + let d = PyDict::new(py); + d.set_item("x", mesh.x.iter().map(|v| *v as f64).collect::>())?; + d.set_item("y", mesh.y.iter().map(|v| *v as f64).collect::>())?; + d.set_item("ni", mesh.ni)?; + d.set_item("nj", mesh.nj)?; + Ok(d) +} + +/// Compute freestream initial conditions for the CFD grid. +/// +/// Args: +/// ni, nj: Grid dimensions +/// mach: Freestream Mach number +/// alpha_deg: Angle of attack (degrees) +/// gamma: Ratio of specific heats +/// physics_mode: 0=Euler, 1=LaminarNS, 2=RANS_SA +/// reynolds: Reynolds number +/// +/// Returns: +/// list[float] of length ni*nj*5, Q=[rho, rho*u, rho*v, E, nu_tilde] +#[pyfunction] +fn cfd_initial_conditions( + ni: u32, + nj: u32, + mach: f32, + alpha_deg: f32, + gamma: f32, + physics_mode: u32, + reynolds: f32, +) -> Vec { + use rustfoil_cfd::config::{CfdConfig, PhysicsMode}; + let config = CfdConfig { + ni, + nj, + mach_inf: mach, + alpha: alpha_deg.to_radians(), + gamma, + reynolds, + physics: match physics_mode { + 1 => PhysicsMode::LaminarNS, + 2 => PhysicsMode::RansSA, + _ => PhysicsMode::Euler, + }, + ..CfdConfig::default() + }; + rustfoil_cfd::init::compute_initial_conditions(&config) +} + +/// Generate boundary condition type array for the O-grid. +/// +/// Returns: list[int] of length ni*nj (0=Interior, 1=Wall, 2=FarField) +#[pyfunction] +fn cfd_boundary_types(ni: u32, nj: u32) -> Vec { + rustfoil_cfd::boundary::generate_bc_types(ni, nj) +} + #[pymodule] fn _rustfoil(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_function(wrap_pyfunction!(analyze_faithful, m)?)?; @@ -515,5 +603,8 @@ fn _rustfoil(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_function(wrap_pyfunction!(deflect_flap, m)?)?; m.add_function(wrap_pyfunction!(parse_dat_file, m)?)?; m.add_function(wrap_pyfunction!(get_bl_distribution, m)?)?; + m.add_function(wrap_pyfunction!(cfd_generate_mesh, m)?)?; + m.add_function(wrap_pyfunction!(cfd_initial_conditions, m)?)?; + m.add_function(wrap_pyfunction!(cfd_boundary_types, m)?)?; Ok(()) } diff --git a/packages/flexfoil-python/src/flexfoil/server.py b/packages/flexfoil-python/src/flexfoil/server.py index 736846fe..fc5f1cfa 100644 --- a/packages/flexfoil-python/src/flexfoil/server.py +++ b/packages/flexfoil-python/src/flexfoil/server.py @@ -245,9 +245,77 @@ async def get_response(self, path: str, scope): raise +# ─── CFD endpoints ────────────────────────────────────────────────────── + +async def cfd_mesh(request: Request) -> JSONResponse: + """Generate a structured O-mesh around an airfoil.""" + from flexfoil._rustfoil import cfd_generate_mesh as _cfd_mesh + + body = await request.json() + coords = body.get("coordinates", []) + ni = int(body.get("ni", 128)) + nj = int(body.get("nj", 64)) + far_field = float(body.get("far_field", 15.0)) + ds0 = float(body.get("ds0", 0.001)) + + # Flatten [{x,y},...] to [x0,y0,x1,y1,...] if needed + if coords and isinstance(coords[0], dict): + flat = [] + for p in coords: + flat.append(p["x"]) + flat.append(p["y"]) + coords = flat + + result = _cfd_mesh(coords, ni, nj, far_field, ds0) + return JSONResponse({ + "x": result["x"], + "y": result["y"], + "ni": result["ni"], + "nj": result["nj"], + }) + + +async def cfd_naca_mesh(request: Request) -> JSONResponse: + """Generate a NACA 4-digit airfoil and mesh it in one call.""" + from flexfoil._rustfoil import generate_naca4, cfd_generate_mesh as _cfd_mesh + + body = await request.json() + code = int(body.get("naca", 12)) # NACA designation as integer + n_pts = int(body.get("n_points", 128)) + ni = int(body.get("ni", 128)) + nj = int(body.get("nj", 64)) + far_field = float(body.get("far_field", 15.0)) + ds0 = float(body.get("ds0", 0.001)) + + # Generate airfoil + foil = generate_naca4(code, n_pts) + coords_flat = [] + for x, y in foil: + coords_flat.append(x) + coords_flat.append(y) + + result = _cfd_mesh(coords_flat, ni, nj, far_field, ds0) + return JSONResponse({ + "x": result["x"], + "y": result["y"], + "ni": result["ni"], + "nj": result["nj"], + "airfoil": foil, + }) + + +async def cfd_page(request: Request) -> Response: + """Serve the CFD mesh visualization page.""" + html = _CFD_HTML + return Response(content=html, media_type="text/html") + + def _build_routes() -> list: routes = [ Route("/api/health", health), + Route("/api/cfd/mesh", cfd_mesh, methods=["POST"]), + Route("/api/cfd/naca-mesh", cfd_naca_mesh, methods=["POST"]), + Route("/cfd", cfd_page), Route("/api/runs", list_runs, methods=["GET"]), Route("/api/runs", create_run, methods=["POST"]), Route("/api/runs", delete_runs, methods=["DELETE"]), @@ -340,4 +408,262 @@ def _open(): print(f" API: {url}/api/health") print() + print(f" CFD Mesh: {url}/cfd") + print() + uvicorn.run(app, host=host, port=port, log_level="info") + + +# ─── Inline CFD mesh visualization frontend ───────────────────────────── + +_CFD_HTML = r""" + + + + +FlexFoil CFD Mesh + + + +
+

FlexFoil CFD Mesh

+ Structured O-grid mesh generation +
+
+
+
+ Airfoil + + +
+
+ Grid + + + + +
+
+ Display + + + +
+ +
Ready
+
+
+
+ +
Scroll to zoom, drag to pan
+
+
+ + + +""" From 809a4ff5b4186c336d03ed1ad161c123ce6db15a Mon Sep 17 00:00:00 2001 From: aeronauty Date: Thu, 26 Mar 2026 13:19:30 +0100 Subject: [PATCH 08/12] fix: improve CFD mesh visualization contrast - Background: #0a0a1a -> #111118 (slightly lighter) - Circumferential grid lines: 0.2 -> 0.55 alpha, teal - Radial grid lines: 0.15 -> 0.45 alpha, blue - Line width: 0.5 -> 0.7 - Airfoil wall: green -> white, lineWidth 2 -> 2.5 - Airfoil points: larger (3px -> 4px), brighter orange Co-Authored-By: Claude Opus 4.6 (1M context) --- .../flexfoil-python/src/flexfoil/server.py | 26 +++++++++---------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/packages/flexfoil-python/src/flexfoil/server.py b/packages/flexfoil-python/src/flexfoil/server.py index fc5f1cfa..9c0868ca 100644 --- a/packages/flexfoil-python/src/flexfoil/server.py +++ b/packages/flexfoil-python/src/flexfoil/server.py @@ -572,7 +572,7 @@ def _open(): ctx.setTransform(dpr, 0, 0, dpr, 0, 0); const W = canvas.clientWidth, H = canvas.clientHeight; - ctx.fillStyle = '#0a0a1a'; + ctx.fillStyle = '#111118'; ctx.fillRect(0, 0, W, H); if (!meshData) return; @@ -587,10 +587,9 @@ def _open(): const sy = (wy) => cy - (wy - viewY) * viewScale; if (showGrid) { - ctx.strokeStyle = 'rgba(0, 180, 140, 0.2)'; - ctx.lineWidth = 0.5; - - // Draw j-lines (circumferential) + // j-lines (circumferential) — teal + ctx.strokeStyle = 'rgba(0, 210, 170, 0.55)'; + ctx.lineWidth = 0.7; for (let j = 0; j < nj; j++) { ctx.beginPath(); for (let i = 0; i <= ni; i++) { @@ -602,8 +601,9 @@ def _open(): ctx.stroke(); } - // Draw i-lines (radial) - ctx.strokeStyle = 'rgba(0, 140, 180, 0.15)'; + // i-lines (radial) — blue + ctx.strokeStyle = 'rgba(60, 160, 220, 0.45)'; + ctx.lineWidth = 0.7; for (let i = 0; i < ni; i++) { ctx.beginPath(); for (let j = 0; j < nj; j++) { @@ -615,9 +615,9 @@ def _open(): } } - // Draw airfoil surface (j=0 wall boundary, thicker) - ctx.strokeStyle = '#00ff88'; - ctx.lineWidth = 2; + // Airfoil surface (j=0 wall boundary) + ctx.strokeStyle = '#ffffff'; + ctx.lineWidth = 2.5; ctx.beginPath(); for (let i = 0; i <= ni; i++) { const ii = i % ni; @@ -626,12 +626,12 @@ def _open(): } ctx.stroke(); - // Draw original airfoil points if available + // Original airfoil points if (showAirfoil && airfoil) { - ctx.fillStyle = '#ff6644'; + ctx.fillStyle = '#ff8855'; for (const [ax, ay] of airfoil) { const px = sx(ax), py = sy(ay); - ctx.fillRect(px - 1.5, py - 1.5, 3, 3); + ctx.fillRect(px - 2, py - 2, 4, 4); } } } From e34909f0298cb96f730b43742ed50a53e528b4f2 Mon Sep 17 00:00:00 2001 From: aeronauty Date: Thu, 26 Mar 2026 13:26:16 +0100 Subject: [PATCH 09/12] fix: set grid line alpha to 1.0 Co-Authored-By: Claude Opus 4.6 (1M context) --- packages/flexfoil-python/src/flexfoil/server.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/packages/flexfoil-python/src/flexfoil/server.py b/packages/flexfoil-python/src/flexfoil/server.py index 9c0868ca..f35a722d 100644 --- a/packages/flexfoil-python/src/flexfoil/server.py +++ b/packages/flexfoil-python/src/flexfoil/server.py @@ -588,7 +588,7 @@ def _open(): if (showGrid) { // j-lines (circumferential) — teal - ctx.strokeStyle = 'rgba(0, 210, 170, 0.55)'; + ctx.strokeStyle = 'rgba(0, 210, 170, 1.0)'; ctx.lineWidth = 0.7; for (let j = 0; j < nj; j++) { ctx.beginPath(); @@ -602,7 +602,7 @@ def _open(): } // i-lines (radial) — blue - ctx.strokeStyle = 'rgba(60, 160, 220, 0.45)'; + ctx.strokeStyle = 'rgba(60, 160, 220, 1.0)'; ctx.lineWidth = 0.7; for (let i = 0; i < ni; i++) { ctx.beginPath(); From 8a9f7c074c8a36707328c9b9e686d73e9e680313 Mon Sep 17 00:00:00 2001 From: aeronauty Date: Thu, 26 Mar 2026 13:44:41 +0100 Subject: [PATCH 10/12] fix: bright solid grid lines, black background Co-Authored-By: Claude Opus 4.6 (1M context) --- packages/flexfoil-python/src/flexfoil/server.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/packages/flexfoil-python/src/flexfoil/server.py b/packages/flexfoil-python/src/flexfoil/server.py index f35a722d..55106fbe 100644 --- a/packages/flexfoil-python/src/flexfoil/server.py +++ b/packages/flexfoil-python/src/flexfoil/server.py @@ -572,7 +572,7 @@ def _open(): ctx.setTransform(dpr, 0, 0, dpr, 0, 0); const W = canvas.clientWidth, H = canvas.clientHeight; - ctx.fillStyle = '#111118'; + ctx.fillStyle = '#000005'; ctx.fillRect(0, 0, W, H); if (!meshData) return; @@ -587,9 +587,9 @@ def _open(): const sy = (wy) => cy - (wy - viewY) * viewScale; if (showGrid) { - // j-lines (circumferential) — teal - ctx.strokeStyle = 'rgba(0, 210, 170, 1.0)'; - ctx.lineWidth = 0.7; + // j-lines (circumferential) — bright cyan + ctx.strokeStyle = '#00e8c0'; + ctx.lineWidth = 1; for (let j = 0; j < nj; j++) { ctx.beginPath(); for (let i = 0; i <= ni; i++) { @@ -601,9 +601,9 @@ def _open(): ctx.stroke(); } - // i-lines (radial) — blue - ctx.strokeStyle = 'rgba(60, 160, 220, 1.0)'; - ctx.lineWidth = 0.7; + // i-lines (radial) — bright sky blue + ctx.strokeStyle = '#50b0ff'; + ctx.lineWidth = 1; for (let i = 0; i < ni; i++) { ctx.beginPath(); for (let j = 0; j < nj; j++) { From 4685b615a59c4f00ef9cbf81117f0913112eeb7f Mon Sep 17 00:00:00 2001 From: aeronauty Date: Thu, 26 Mar 2026 13:46:51 +0100 Subject: [PATCH 11/12] fix: DPR-aware line width (1.5px min), brighter colors Co-Authored-By: Claude Opus 4.6 (1M context) --- packages/flexfoil-python/src/flexfoil/server.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/packages/flexfoil-python/src/flexfoil/server.py b/packages/flexfoil-python/src/flexfoil/server.py index 55106fbe..11eca0a9 100644 --- a/packages/flexfoil-python/src/flexfoil/server.py +++ b/packages/flexfoil-python/src/flexfoil/server.py @@ -587,9 +587,11 @@ def _open(): const sy = (wy) => cy - (wy - viewY) * viewScale; if (showGrid) { + const lw = Math.max(1.5, 3 / (window.devicePixelRatio || 1)); + // j-lines (circumferential) — bright cyan - ctx.strokeStyle = '#00e8c0'; - ctx.lineWidth = 1; + ctx.strokeStyle = '#00ffd0'; + ctx.lineWidth = lw; for (let j = 0; j < nj; j++) { ctx.beginPath(); for (let i = 0; i <= ni; i++) { @@ -602,8 +604,8 @@ def _open(): } // i-lines (radial) — bright sky blue - ctx.strokeStyle = '#50b0ff'; - ctx.lineWidth = 1; + ctx.strokeStyle = '#66ccff'; + ctx.lineWidth = lw; for (let i = 0; i < ni; i++) { ctx.beginPath(); for (let j = 0; j < nj; j++) { From bb0d4d5e4210b623b3f6eb9cf5b64376aa22cf2b Mon Sep 17 00:00:00 2001 From: aeronauty Date: Thu, 26 Mar 2026 13:51:42 +0100 Subject: [PATCH 12/12] fix: remove DPR canvas scaling to fix invisible grid lines The setTransform(dpr,...) approach was causing sub-pixel rendering issues. Switch to 1:1 pixel mapping (canvas.width = clientWidth) which ensures lines are always at least 1 physical pixel. Co-Authored-By: Claude Opus 4.6 (1M context) --- packages/flexfoil-python/src/flexfoil/server.py | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/packages/flexfoil-python/src/flexfoil/server.py b/packages/flexfoil-python/src/flexfoil/server.py index 11eca0a9..698a48fb 100644 --- a/packages/flexfoil-python/src/flexfoil/server.py +++ b/packages/flexfoil-python/src/flexfoil/server.py @@ -566,12 +566,12 @@ def _open(): function draw() { const canvas = $('canvas'); const ctx = canvas.getContext('2d'); - const dpr = window.devicePixelRatio || 1; - canvas.width = canvas.clientWidth * dpr; - canvas.height = canvas.clientHeight * dpr; - ctx.setTransform(dpr, 0, 0, dpr, 0, 0); + // Use 1:1 pixel mapping — no DPR scaling (avoids sub-pixel line issues) + const W = canvas.clientWidth; + const H = canvas.clientHeight; + canvas.width = W; + canvas.height = H; - const W = canvas.clientWidth, H = canvas.clientHeight; ctx.fillStyle = '#000005'; ctx.fillRect(0, 0, W, H); @@ -587,11 +587,9 @@ def _open(): const sy = (wy) => cy - (wy - viewY) * viewScale; if (showGrid) { - const lw = Math.max(1.5, 3 / (window.devicePixelRatio || 1)); - // j-lines (circumferential) — bright cyan ctx.strokeStyle = '#00ffd0'; - ctx.lineWidth = lw; + ctx.lineWidth = 1; for (let j = 0; j < nj; j++) { ctx.beginPath(); for (let i = 0; i <= ni; i++) { @@ -605,7 +603,7 @@ def _open(): // i-lines (radial) — bright sky blue ctx.strokeStyle = '#66ccff'; - ctx.lineWidth = lw; + ctx.lineWidth = 1; for (let i = 0; i < ni; i++) { ctx.beginPath(); for (let j = 0; j < nj; j++) {