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 4478fbe..90e742a 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_mut()
+            .as_ref()
             .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_mut()
+            .as_ref()
             .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?;
 
         // Copy from GPU
@@ -191,7 +191,7 @@ impl PisoGpuSolver {
     async fn gpu_momentum_prediction(&mut self, dt: f32) -> CfdResult<()> {
         let buffers = self
             .gpu_buffers
-            .as_mut()
+            .as_ref()
             .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_mut()
+                .as_ref()
                 .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?;
             (buffers.nx, buffers.ny)
         };
@@ -284,7 +284,7 @@ impl PisoGpuSolver {
         let iterations = {
             let buffers = self
                 .gpu_buffers
-                .as_mut()
+                .as_ref()
                 .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?;
             self.poisson_kernel.solve_2d(
                 &mut buffers.pressure_correction,
@@ -303,7 +303,7 @@ impl PisoGpuSolver {
         {
             let buffers = self
                 .gpu_buffers
-                .as_mut()
+                .as_ref()
                 .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?;
             self.matrix_ops_kernel
                 .axpy(1.0, &buffers.pressure_correction, &mut buffers.p)?;
@@ -313,7 +313,7 @@ impl PisoGpuSolver {
         let residual_norm = {
             let buffers = self
                 .gpu_buffers
-                .as_mut()
+                .as_ref()
                 .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?;
             self.matrix_ops_kernel.vector_norm(&buffers.mass_source)?
         };
@@ -332,7 +332,7 @@ impl PisoGpuSolver {
     ) -> CfdResult<()> {
         let buffers = self
             .gpu_buffers
-            .as_mut()
+            .as_ref()
             .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?;
 
         // Clone the references we need for velocity field
@@ -416,7 +416,7 @@ impl PisoGpuSolver {
         {
             let buffers = self
                 .gpu_buffers
-                .as_mut()
+                .as_ref()
                 .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 777eda7..fd649ac 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_mut()
+            .as_ref()
             .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_mut()
+            .as_ref()
             .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?;
 
         // Copy from GPU
@@ -183,7 +183,7 @@ impl SimpleGpuSolver {
     async fn gpu_momentum_prediction(&mut self, dt: f32) -> CfdResult<()> {
         let buffers = self
             .gpu_buffers
-            .as_mut()
+            .as_ref()
             .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_mut()
+            .as_ref()
             .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?;
 
         let config = self.cpu_solver.config();
@@ -282,7 +282,7 @@ impl SimpleGpuSolver {
     async fn gpu_velocity_correction(&mut self) -> CfdResult<()> {
         let buffers = self
             .gpu_buffers
-            .as_mut()
+            .as_ref()
             .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?;
 
         // This step would require a custom velocity correction kernel
@@ -300,7 +300,7 @@ impl SimpleGpuSolver {
     async fn gpu_pressure_update(&self, pressure_relaxation: f32) -> CfdResult<()> {
         let buffers = self
             .gpu_buffers
-            .as_mut()
+            .as_ref()
             .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 4894240..5de38cf 100644
--- a/crates/specialized/rtx-cfd/src/solvers/lbm/d2q9_gpu.rs
+++ b/crates/specialized/rtx-cfd/src/solvers/lbm/d2q9_gpu.rs
@@ -253,7 +253,7 @@ impl D2Q9GpuSolver {
     }
 
     /// GPU-accelerated collision step
-    fn gpu_collision_step(&self) -> CfdResult<()> {
+    fn gpu_collision_step(&mut self) -> CfdResult<()> {
         let buffers = self
             .gpu_buffers
             .as_mut()
@@ -300,7 +300,7 @@ impl D2Q9GpuSolver {
     }
 
     /// GPU-accelerated streaming step
-    pub fn gpu_streaming_step(&self) -> CfdResult<()> {
+    pub fn gpu_streaming_step(&mut self) -> CfdResult<()> {
         let buffers = self
             .gpu_buffers
             .as_mut()
@@ -344,7 +344,7 @@ impl D2Q9GpuSolver {
     }
 
     /// GPU-accelerated macroscopic variable computation
-    fn gpu_compute_macroscopic(&self) -> CfdResult<()> {
+    fn gpu_compute_macroscopic(&mut self) -> CfdResult<()> {
         let buffers = self
             .gpu_buffers
             .as_mut()
@@ -391,7 +391,7 @@ impl D2Q9GpuSolver {
     }
 
     /// Compute equilibrium distributions on GPU
-    fn gpu_compute_equilibrium(&self) -> CfdResult<()> {
+    fn gpu_compute_equilibrium(&mut self) -> CfdResult<()> {
         let buffers = self
             .gpu_buffers
             .as_mut()
@@ -435,7 +435,7 @@ impl D2Q9GpuSolver {
     }
 
     /// GPU-accelerated bounce-back boundary conditions
-    pub fn gpu_apply_bounce_back_boundaries(&self) -> CfdResult<()> {
+    pub fn gpu_apply_bounce_back_boundaries(&mut self) -> CfdResult<()> {
         let buffers = self
             .gpu_buffers
             .as_mut()
@@ -539,7 +539,7 @@ impl D2Q9GpuSolver {
     pub fn get_macroscopic_at(&self, i: usize, j: usize) -> CfdResult<(f64, Vector2<f64>)> {
         let buffers = self
             .gpu_buffers
-            .as_mut()
+            .as_ref()
             .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_mut()
+            .as_ref()
             .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 9e5987c..c5a3fa7 100644
--- a/crates/specialized/rtx-cfd/src/solvers/lbm/d3q19_gpu.rs
+++ b/crates/specialized/rtx-cfd/src/solvers/lbm/d3q19_gpu.rs
@@ -362,7 +362,7 @@ __global__ void d3q19_bounce_back_boundaries(
     }
 
     /// GPU-accelerated BGK collision step
-    pub fn gpu_collision_step(&self) -> CfdResult<()> {
+    pub fn gpu_collision_step(&mut self) -> CfdResult<()> {
         let buffers = self
             .gpu_buffers
             .as_mut()
@@ -416,7 +416,7 @@ __global__ void d3q19_bounce_back_boundaries(
     }
 
     /// GPU-accelerated streaming step
-    pub fn gpu_streaming_step(&self) -> CfdResult<()> {
+    pub fn gpu_streaming_step(&mut self) -> CfdResult<()> {
         let buffers = self
             .gpu_buffers
             .as_mut()
@@ -462,7 +462,7 @@ __global__ void d3q19_bounce_back_boundaries(
     }
 
     /// Extract macroscopic variables on GPU
-    fn gpu_extract_macroscopic_variables(&self) -> CfdResult<()> {
+    fn gpu_extract_macroscopic_variables(&mut self) -> CfdResult<()> {
         let buffers = self
             .gpu_buffers
             .as_mut()
@@ -517,7 +517,7 @@ __global__ void d3q19_bounce_back_boundaries(
     }
 
     /// Compute equilibrium distributions on GPU
-    fn gpu_compute_equilibrium(&self) -> CfdResult<()> {
+    fn gpu_compute_equilibrium(&mut self) -> CfdResult<()> {
         let buffers = self
             .gpu_buffers
             .as_mut()
@@ -564,7 +564,7 @@ __global__ void d3q19_bounce_back_boundaries(
     }
 
     /// GPU-accelerated bounce-back boundary conditions
-    pub fn gpu_apply_bounce_back_boundaries(&self) -> CfdResult<()> {
+    pub fn gpu_apply_bounce_back_boundaries(&mut self) -> CfdResult<()> {
         let buffers = self
             .gpu_buffers
             .as_mut()
@@ -675,7 +675,7 @@ __global__ void d3q19_bounce_back_boundaries(
     ) -> CfdResult<MacroscopicVariables3D> {
         let buffers = self
             .gpu_buffers
-            .as_mut()
+            .as_ref()
             .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<f64> {
         let buffers = self
             .gpu_buffers
-            .as_mut()
+            .as_ref()
             .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<f64> {
         let buffers = self
             .gpu_buffers
-            .as_mut()
+            .as_ref()
             .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 f5e7aa4..bd67bb5 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_mut()
+            .as_ref()
             .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?;
 
         // Initialize k and epsilon with uniform values
@@ -429,7 +429,7 @@ impl KEpsilonGpuModel {
     pub fn get_eddy_viscosity(&self) -> CfdResult<Vec<f32>> {
         let buffers = self
             .gpu_buffers
-            .as_mut()
+            .as_ref()
             .ok_or_else(|| CfdError::gpu_error("GPU buffers not initialized"))?;
 
         let n = buffers.nx * buffers.ny * buffers.nz;