From 0de750bad1164206ec5ab5902c3bf7f3c6bfd091 Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Wed, 16 Sep 2026 07:28:42 -0500 Subject: [PATCH] =?UTF-8?q?rtx-cfd=20legacy=20GPU:=20kernel=20outputs=20we?= =?UTF-8?q?re=20device-buffer=20CLONES=20(every=20write=20lost)=20?= =?UTF-8?q?=E2=80=94=20the=20buffers=20are=20borrowed=20mutably=20now;=20t?= =?UTF-8?q?he=20Poisson=20kernels=20keep=20boundary=20Dirichlet=20values;?= =?UTF-8?q?=20tests:=20the=20advection=20pulse=20marched=2050=20steps,=20t?= =?UTF-8?q?he=20Jacobi=20budget=2020k,=20the=20reduction=20reference=20in?= =?UTF-8?q?=20f64?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../rtx-cfd/src/kernels/cuda/poisson.cu | 6 +-- .../src/solvers/incompressible/piso_gpu.rs | 12 +++--- .../src/solvers/incompressible/simple_gpu.rs | 6 +-- .../rtx-cfd/src/solvers/lbm/d2q9_gpu.rs | 38 +++++++++---------- .../rtx-cfd/src/solvers/lbm/d3q19_gpu.rs | 36 +++++++++--------- .../rtx-cfd/src/turbulence/k_epsilon_gpu.rs | 2 +- .../rtx-cfd/tests/gpu_reduction_tests.rs | 4 +- .../rtx-cfd/tests/kernels_tests.rs | 16 ++++---- 8 files changed, 58 insertions(+), 62 deletions(-) diff --git a/crates/specialized/rtx-cfd/src/kernels/cuda/poisson.cu b/crates/specialized/rtx-cfd/src/kernels/cuda/poisson.cu index 64d425d..2cc6937 100644 --- a/crates/specialized/rtx-cfd/src/kernels/cuda/poisson.cu +++ b/crates/specialized/rtx-cfd/src/kernels/cuda/poisson.cu @@ -25,9 +25,9 @@ extern "C" __global__ void poisson_jacobi_2d( int idx = j * nx + i; - // Boundary conditions (Dirichlet - zero on boundaries for now) + // Boundary cells keep the Dirichlet values they hold. if (i == 0 || i == nx - 1 || j == 0 || j == ny - 1) { - phi_new[idx] = 0.0f; + phi_new[idx] = phi[idx]; return; } @@ -66,7 +66,6 @@ extern "C" __global__ void poisson_gauss_seidel_2d( // Boundary conditions if (i == 0 || i == nx - 1 || j == 0 || j == ny - 1) { - phi[idx] = 0.0f; return; } @@ -105,7 +104,6 @@ extern "C" __global__ void poisson_sor_2d( // Boundary conditions if (i == 0 || i == nx - 1 || j == 0 || j == ny - 1) { - phi[idx] = 0.0f; return; } diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/piso_gpu.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/piso_gpu.rs index 6b2e648..4478fbe 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/piso_gpu.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/piso_gpu.rs @@ -147,7 +147,7 @@ impl PisoGpuSolver { fn copy_to_gpu(&self, flow_field: &FlowField) -> CfdResult<()> { let buffers = self .gpu_buffers - .as_ref() + .as_mut() .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?; // Convert nalgebra matrices to flat vectors @@ -171,7 +171,7 @@ impl PisoGpuSolver { fn copy_from_gpu(&self, flow_field: &mut FlowField) -> CfdResult<()> { let buffers = self .gpu_buffers - .as_ref() + .as_mut() .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?; // Copy from GPU @@ -271,7 +271,7 @@ impl PisoGpuSolver { let (nx, ny) = { let buffers = self .gpu_buffers - .as_ref() + .as_mut() .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?; (buffers.nx, buffers.ny) }; @@ -313,7 +313,7 @@ impl PisoGpuSolver { let residual_norm = { let buffers = self .gpu_buffers - .as_ref() + .as_mut() .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?; self.matrix_ops_kernel.vector_norm(&buffers.mass_source)? }; @@ -464,8 +464,8 @@ impl PisoGpuSolver { .stream() .launch_builder(&func) .arg(&buffers.pressure_correction) - .arg(&mut buffers.u_correction.clone()) - .arg(&mut buffers.v_correction.clone()) + .arg(&mut buffers.u_correction) + .arg(&mut buffers.v_correction) .arg(&correction_factor) .arg(&dx) .arg(&dy) diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/simple_gpu.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/simple_gpu.rs index b09cb0b..777eda7 100644 --- a/crates/specialized/rtx-cfd/src/solvers/incompressible/simple_gpu.rs +++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/simple_gpu.rs @@ -136,7 +136,7 @@ impl SimpleGpuSolver { fn copy_to_gpu(&self, flow_field: &FlowField) -> CfdResult<()> { let buffers = self .gpu_buffers - .as_ref() + .as_mut() .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?; // Convert nalgebra matrices to flat vectors @@ -163,7 +163,7 @@ impl SimpleGpuSolver { fn copy_from_gpu(&self, flow_field: &mut FlowField) -> CfdResult<()> { let buffers = self .gpu_buffers - .as_ref() + .as_mut() .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?; // Copy from GPU @@ -300,7 +300,7 @@ impl SimpleGpuSolver { async fn gpu_pressure_update(&self, pressure_relaxation: f32) -> CfdResult<()> { let buffers = self .gpu_buffers - .as_ref() + .as_mut() .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?; // p = p + α_p * p' (pressure update with relaxation) diff --git a/crates/specialized/rtx-cfd/src/solvers/lbm/d2q9_gpu.rs b/crates/specialized/rtx-cfd/src/solvers/lbm/d2q9_gpu.rs index 1052da7..4894240 100644 --- a/crates/specialized/rtx-cfd/src/solvers/lbm/d2q9_gpu.rs +++ b/crates/specialized/rtx-cfd/src/solvers/lbm/d2q9_gpu.rs @@ -256,7 +256,7 @@ impl D2Q9GpuSolver { fn gpu_collision_step(&self) -> CfdResult<()> { let buffers = self .gpu_buffers - .as_ref() + .as_mut() .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?; // First compute equilibrium distributions @@ -284,7 +284,7 @@ impl D2Q9GpuSolver { self.kernel_manager .stream() .launch_builder(&func) - .arg(&mut buffers.f.clone()) + .arg(&mut buffers.f) .arg(&buffers.f_eq) .arg(&(self.omega as f32)) .arg(&(buffers.nx as i32)) @@ -303,7 +303,7 @@ impl D2Q9GpuSolver { pub fn gpu_streaming_step(&self) -> CfdResult<()> { let buffers = self .gpu_buffers - .as_ref() + .as_mut() .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?; let module = self @@ -327,7 +327,7 @@ impl D2Q9GpuSolver { self.kernel_manager .stream() .launch_builder(&func) - .arg(&mut buffers.f_temp.clone()) + .arg(&mut buffers.f_temp) .arg(&buffers.f) .arg(&(buffers.nx as i32)) .arg(&(buffers.ny as i32)) @@ -347,7 +347,7 @@ impl D2Q9GpuSolver { fn gpu_compute_macroscopic(&self) -> CfdResult<()> { let buffers = self .gpu_buffers - .as_ref() + .as_mut() .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?; let module = self @@ -371,9 +371,9 @@ impl D2Q9GpuSolver { self.kernel_manager .stream() .launch_builder(&func) - .arg(&mut buffers.density.clone()) - .arg(&mut buffers.velocity_x.clone()) - .arg(&mut buffers.velocity_y.clone()) + .arg(&mut buffers.density) + .arg(&mut buffers.velocity_x) + .arg(&mut buffers.velocity_y) .arg(&buffers.f) .arg(&(buffers.nx as i32)) .arg(&(buffers.ny as i32)) @@ -394,7 +394,7 @@ impl D2Q9GpuSolver { fn gpu_compute_equilibrium(&self) -> CfdResult<()> { let buffers = self .gpu_buffers - .as_ref() + .as_mut() .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?; let module = self @@ -418,7 +418,7 @@ impl D2Q9GpuSolver { self.kernel_manager .stream() .launch_builder(&func) - .arg(&mut buffers.f_eq.clone()) + .arg(&mut buffers.f_eq) .arg(&buffers.density) .arg(&buffers.velocity_x) .arg(&buffers.velocity_y) @@ -438,7 +438,7 @@ impl D2Q9GpuSolver { pub fn gpu_apply_bounce_back_boundaries(&self) -> CfdResult<()> { let buffers = self .gpu_buffers - .as_ref() + .as_mut() .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?; let module = self @@ -462,7 +462,7 @@ impl D2Q9GpuSolver { self.kernel_manager .stream() .launch_builder(&func) - .arg(&mut buffers.f.clone()) + .arg(&mut buffers.f) .arg(&(buffers.nx as i32)) .arg(&(buffers.ny as i32)) .launch(config) @@ -492,7 +492,7 @@ impl D2Q9GpuSolver { ) -> CfdResult<()> { let buffers = self .gpu_buffers - .as_ref() + .as_mut() .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?; let module = self @@ -516,10 +516,10 @@ impl D2Q9GpuSolver { self.kernel_manager .stream() .launch_builder(&func) - .arg(&mut buffers.f.clone()) - .arg(&mut buffers.density.clone()) - .arg(&mut buffers.velocity_x.clone()) - .arg(&mut buffers.velocity_y.clone()) + .arg(&mut buffers.f) + .arg(&mut buffers.density) + .arg(&mut buffers.velocity_x) + .arg(&mut buffers.velocity_y) .arg(&(density as f32)) .arg(&(velocity.x as f32)) .arg(&(velocity.y as f32)) @@ -539,7 +539,7 @@ impl D2Q9GpuSolver { pub fn get_macroscopic_at(&self, i: usize, j: usize) -> CfdResult<(f64, Vector2)> { let buffers = self .gpu_buffers - .as_ref() + .as_mut() .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?; let idx = j * buffers.nx + i; @@ -585,7 +585,7 @@ impl D2Q9GpuSolver { pub fn update_flow_field(&self, flow_field: &mut FlowField) -> CfdResult<()> { let buffers = self .gpu_buffers - .as_ref() + .as_mut() .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?; let (nx, ny, _, _) = flow_field.grid_info(); diff --git a/crates/specialized/rtx-cfd/src/solvers/lbm/d3q19_gpu.rs b/crates/specialized/rtx-cfd/src/solvers/lbm/d3q19_gpu.rs index a202696..9e5987c 100644 --- a/crates/specialized/rtx-cfd/src/solvers/lbm/d3q19_gpu.rs +++ b/crates/specialized/rtx-cfd/src/solvers/lbm/d3q19_gpu.rs @@ -365,7 +365,7 @@ __global__ void d3q19_bounce_back_boundaries( pub fn gpu_collision_step(&self) -> CfdResult<()> { let buffers = self .gpu_buffers - .as_ref() + .as_mut() .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?; let omega = 1.0 / self.params.tau; @@ -399,7 +399,7 @@ __global__ void d3q19_bounce_back_boundaries( self.kernel_manager .stream() .launch_builder(&func) - .arg(&mut buffers.f.clone()) + .arg(&mut buffers.f) .arg(&buffers.f_eq) .arg(&(omega as f32)) .arg(&(buffers.nx as i32)) @@ -419,7 +419,7 @@ __global__ void d3q19_bounce_back_boundaries( pub fn gpu_streaming_step(&self) -> CfdResult<()> { let buffers = self .gpu_buffers - .as_ref() + .as_mut() .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?; let module = self @@ -444,7 +444,7 @@ __global__ void d3q19_bounce_back_boundaries( self.kernel_manager .stream() .launch_builder(&func) - .arg(&mut buffers.f_temp.clone()) + .arg(&mut buffers.f_temp) .arg(&buffers.f) .arg(&(buffers.nx as i32)) .arg(&(buffers.ny as i32)) @@ -465,7 +465,7 @@ __global__ void d3q19_bounce_back_boundaries( fn gpu_extract_macroscopic_variables(&self) -> CfdResult<()> { let buffers = self .gpu_buffers - .as_ref() + .as_mut() .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?; let module = self @@ -495,10 +495,10 @@ __global__ void d3q19_bounce_back_boundaries( self.kernel_manager .stream() .launch_builder(&func) - .arg(&mut buffers.density.clone()) - .arg(&mut buffers.velocity_x.clone()) - .arg(&mut buffers.velocity_y.clone()) - .arg(&mut buffers.velocity_z.clone()) + .arg(&mut buffers.density) + .arg(&mut buffers.velocity_x) + .arg(&mut buffers.velocity_y) + .arg(&mut buffers.velocity_z) .arg(&buffers.f) .arg(&(buffers.nx as i32)) .arg(&(buffers.ny as i32)) @@ -520,7 +520,7 @@ __global__ void d3q19_bounce_back_boundaries( fn gpu_compute_equilibrium(&self) -> CfdResult<()> { let buffers = self .gpu_buffers - .as_ref() + .as_mut() .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?; let module = self @@ -545,7 +545,7 @@ __global__ void d3q19_bounce_back_boundaries( self.kernel_manager .stream() .launch_builder(&func) - .arg(&mut buffers.f_eq.clone()) + .arg(&mut buffers.f_eq) .arg(&buffers.density) .arg(&buffers.velocity_x) .arg(&buffers.velocity_y) @@ -567,7 +567,7 @@ __global__ void d3q19_bounce_back_boundaries( pub fn gpu_apply_bounce_back_boundaries(&self) -> CfdResult<()> { let buffers = self .gpu_buffers - .as_ref() + .as_mut() .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?; let module = self @@ -591,7 +591,7 @@ __global__ void d3q19_bounce_back_boundaries( self.kernel_manager .stream() .launch_builder(&func) - .arg(&mut buffers.f.clone()) + .arg(&mut buffers.f) .arg(&(buffers.nx as i32)) .arg(&(buffers.ny as i32)) .arg(&(buffers.nz as i32)) @@ -621,7 +621,7 @@ __global__ void d3q19_bounce_back_boundaries( ) -> CfdResult<()> { let buffers = self .gpu_buffers - .as_ref() + .as_mut() .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?; let module = self @@ -648,7 +648,7 @@ __global__ void d3q19_bounce_back_boundaries( self.kernel_manager .stream() .launch_builder(&func) - .arg(&mut buffers.f.clone()) + .arg(&mut buffers.f) .arg(&(density as f32)) .arg(&(velocity.x as f32)) .arg(&(velocity.y as f32)) @@ -675,7 +675,7 @@ __global__ void d3q19_bounce_back_boundaries( ) -> CfdResult { let buffers = self .gpu_buffers - .as_ref() + .as_mut() .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?; // Extract macroscopic variables first @@ -726,7 +726,7 @@ __global__ void d3q19_bounce_back_boundaries( pub fn gpu_total_mass(&self) -> CfdResult { let buffers = self .gpu_buffers - .as_ref() + .as_mut() .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?; self.gpu_extract_macroscopic_variables()?; @@ -739,7 +739,7 @@ __global__ void d3q19_bounce_back_boundaries( pub fn gpu_kinetic_energy(&self) -> CfdResult { let buffers = self .gpu_buffers - .as_ref() + .as_mut() .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?; self.gpu_extract_macroscopic_variables()?; diff --git a/crates/specialized/rtx-cfd/src/turbulence/k_epsilon_gpu.rs b/crates/specialized/rtx-cfd/src/turbulence/k_epsilon_gpu.rs index dc8e5d0..f5e7aa4 100644 --- a/crates/specialized/rtx-cfd/src/turbulence/k_epsilon_gpu.rs +++ b/crates/specialized/rtx-cfd/src/turbulence/k_epsilon_gpu.rs @@ -429,7 +429,7 @@ impl KEpsilonGpuModel { pub fn get_eddy_viscosity(&self) -> CfdResult> { let buffers = self .gpu_buffers - .as_ref() + .as_mut() .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?; let n = buffers.nx * buffers.ny * buffers.nz; diff --git a/crates/specialized/rtx-cfd/tests/gpu_reduction_tests.rs b/crates/specialized/rtx-cfd/tests/gpu_reduction_tests.rs index cd02c08..03bd221 100644 --- a/crates/specialized/rtx-cfd/tests/gpu_reduction_tests.rs +++ b/crates/specialized/rtx-cfd/tests/gpu_reduction_tests.rs @@ -37,7 +37,7 @@ mod gpu_reduction_tests { // Create test data let test_data: Vec = (1..=100).map(|i| i as f32).collect(); - let expected_sum = test_data.iter().sum::(); + let expected_sum = test_data.iter().map(|&v| v as f64).sum::() as f32; // f64 reference: a sequential f32 sum of a million terms carries ~1e-3 of rounding // Copy to GPU let gpu_data = manager.copy_to_device(&test_data)?; @@ -100,7 +100,7 @@ mod gpu_reduction_tests { // Create large test array let n = 1_000_000; let test_data: Vec = (0..n).map(|i| (i % 100) as f32).collect(); - let expected_sum = test_data.iter().sum::(); + let expected_sum = test_data.iter().map(|&v| v as f64).sum::() as f32; // f64 reference: a sequential f32 sum of a million terms carries ~1e-3 of rounding // Copy to GPU let gpu_data = manager.copy_to_device(&test_data)?; diff --git a/crates/specialized/rtx-cfd/tests/kernels_tests.rs b/crates/specialized/rtx-cfd/tests/kernels_tests.rs index 8f536b3..a84975c 100644 --- a/crates/specialized/rtx-cfd/tests/kernels_tests.rs +++ b/crates/specialized/rtx-cfd/tests/kernels_tests.rs @@ -46,17 +46,15 @@ mod cuda_tests { d_phi = kernel_manager.copy_to_device(&phi)?; // Run advection kernel - advection_kernel.apply( - &d_phi, - &mut d_phi_new, - velocity as f32, - dt as f32, - dx as f32, - )?; + // CFL = v dt / dx ≈ 0.1: march 50 steps so the pulse moves several cells. + for _ in 0..50 { + advection_kernel.apply(&d_phi, &mut d_phi_new, velocity as f32, dt as f32, dx as f32)?; + std::mem::swap(&mut d_phi, &mut d_phi_new); + } // Copy result back let mut result = vec![0.0f32; nx]; - result = kernel_manager.copy_from_device(&d_phi_new)?; + result = kernel_manager.copy_from_device(&d_phi)?; // Verify that the pulse has moved (mass conservation) let initial_mass: f32 = phi.iter().sum(); @@ -206,7 +204,7 @@ mod cuda_tests { d_source = kernel_manager.copy_to_device(&source)?; // Solve Poisson equation - let max_iterations = 1000; + let max_iterations = 20_000; // Jacobi needs O(n²) sweeps on this grid let tolerance = 1e-6; let iterations = poisson_kernel.solve_2d( &mut d_phi,