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 90e742a..2737f0e 100644
--- a/crates/specialized/rtx-cfd/src/solvers/incompressible/piso_gpu.rs
+++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/piso_gpu.rs
@@ -191,7 +191,7 @@ impl PisoGpuSolver {
     async fn gpu_momentum_prediction(&mut self, dt: f32) -> CfdResult<()> {
         let buffers = self
             .gpu_buffers
-            .as_ref()
+            .as_mut()
             .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?;
 
         let config = self.cpu_solver.config();
@@ -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)
         };
@@ -416,7 +416,7 @@ impl PisoGpuSolver {
         {
             let buffers = self
                 .gpu_buffers
-                .as_ref()
+                .as_mut()
                 .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?;
             let u_correction = buffers.u_correction.clone();
             let v_correction = buffers.v_correction.clone();
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 fd649ac..d2b2162 100644
--- a/crates/specialized/rtx-cfd/src/solvers/incompressible/simple_gpu.rs
+++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/simple_gpu.rs
@@ -183,7 +183,7 @@ impl SimpleGpuSolver {
     async fn gpu_momentum_prediction(&mut self, dt: f32) -> CfdResult<()> {
         let buffers = self
             .gpu_buffers
-            .as_ref()
+            .as_mut()
             .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?;
 
         let config = self.cpu_solver.config();
@@ -250,7 +250,7 @@ impl SimpleGpuSolver {
     async fn gpu_pressure_correction(&mut self) -> CfdResult<f32> {
         let buffers = self
             .gpu_buffers
-            .as_ref()
+            .as_mut()
             .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?;
 
         let config = self.cpu_solver.config();
@@ -297,10 +297,10 @@ impl SimpleGpuSolver {
     }
 
     /// GPU-accelerated pressure update step
-    async fn gpu_pressure_update(&self, pressure_relaxation: f32) -> CfdResult<()> {
+    async fn gpu_pressure_update(&mut 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 5de38cf..77300e6 100644
--- a/crates/specialized/rtx-cfd/src/solvers/lbm/d2q9_gpu.rs
+++ b/crates/specialized/rtx-cfd/src/solvers/lbm/d2q9_gpu.rs
@@ -256,11 +256,15 @@ impl D2Q9GpuSolver {
     fn gpu_collision_step(&mut self) -> CfdResult<()> {
         let buffers = self
             .gpu_buffers
-            .as_mut()
+            .as_ref()
             .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?;
 
         // First compute equilibrium distributions
         self.gpu_compute_equilibrium()?;
+        let buffers = self
+            .gpu_buffers
+            .as_mut()
+            .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?;
 
         // Then perform collision: f = f + omega * (f_eq - f)
         let module = self
@@ -476,7 +480,7 @@ impl D2Q9GpuSolver {
     }
 
     /// Complete GPU LBM time step
-    pub fn gpu_step(&self) -> CfdResult<()> {
+    pub fn gpu_step(&mut self) -> CfdResult<()> {
         self.gpu_collision_step()?;
         self.gpu_streaming_step()?;
         self.gpu_apply_bounce_back_boundaries()?;
@@ -575,7 +579,7 @@ impl D2Q9GpuSolver {
     }
 
     /// Apply boundary conditions
-    fn gpu_apply_boundaries(&self, _boundary_conditions: &BoundaryConditions) -> CfdResult<()> {
+    fn gpu_apply_boundaries(&mut self, _boundary_conditions: &BoundaryConditions) -> CfdResult<()> {
         // For now, just apply bounce-back
         self.gpu_apply_bounce_back_boundaries()?;
         Ok(())
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 c5a3fa7..7303630 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(&mut self) -> CfdResult<()> {
         let buffers = self
             .gpu_buffers
-            .as_mut()
+            .as_ref()
             .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?;
 
         let omega = 1.0 / self.params.tau;
@@ -375,6 +375,10 @@ __global__ void d3q19_bounce_back_boundaries(
 
         // Step 2: Compute equilibrium distributions
         self.gpu_compute_equilibrium()?;
+        let buffers = self
+            .gpu_buffers
+            .as_mut()
+            .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?;
 
         // Step 3: Perform BGK collision
         let module = self
@@ -606,7 +610,7 @@ __global__ void d3q19_bounce_back_boundaries(
     }
 
     /// Complete GPU LBM time step
-    pub fn gpu_step(&self) -> CfdResult<()> {
+    pub fn gpu_step(&mut self) -> CfdResult<()> {
         self.gpu_collision_step()?;
         self.gpu_streaming_step()?;
         self.gpu_apply_bounce_back_boundaries()?;
@@ -667,8 +671,7 @@ __global__ void d3q19_bounce_back_boundaries(
     }
 
     /// Get macroscopic variables at a specific point (copy from GPU)
-    pub fn gpu_macroscopic_variables_at(
-        &self,
+    pub fn gpu_macroscopic_variables_at(&mut self,
         x: usize,
         y: usize,
         z: usize,
@@ -680,6 +683,10 @@ __global__ void d3q19_bounce_back_boundaries(
 
         // Extract macroscopic variables first
         self.gpu_extract_macroscopic_variables()?;
+        let buffers = self
+            .gpu_buffers
+            .as_ref()
+            .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?;
 
         // Copy arrays from GPU
         let density_host = self.kernel_manager.copy_from_device(&buffers.density)?;
@@ -723,7 +730,7 @@ __global__ void d3q19_bounce_back_boundaries(
     }
 
     /// Calculate total mass on GPU
-    pub fn gpu_total_mass(&self) -> CfdResult<f64> {
+    pub fn gpu_total_mass(&mut self) -> CfdResult<f64> {
         let buffers = self
             .gpu_buffers
             .as_ref()
@@ -736,7 +743,7 @@ __global__ void d3q19_bounce_back_boundaries(
     }
 
     /// Calculate kinetic energy on GPU
-    pub fn gpu_kinetic_energy(&self) -> CfdResult<f64> {
+    pub fn gpu_kinetic_energy(&mut self) -> CfdResult<f64> {
         let buffers = self
             .gpu_buffers
             .as_ref()
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 bd67bb5..dc8e5d0 100644
--- a/crates/specialized/rtx-cfd/src/turbulence/k_epsilon_gpu.rs
+++ b/crates/specialized/rtx-cfd/src/turbulence/k_epsilon_gpu.rs
@@ -131,7 +131,7 @@ impl KEpsilonGpuModel {
     pub fn set_initial_conditions(&mut self, k_init: f32, epsilon_init: f32) -> CfdResult<()> {
         let buffers = self
             .gpu_buffers
-            .as_ref()
+            .as_mut()
             .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?;
 
         // Initialize k and epsilon with uniform values