From cdcc0afded55bd30537366c3225e4d786e407c77 Mon Sep 17 00:00:00 2001
From: Aaron Councilman <aaronjc4@illinois.edu>
Date: Sat, 1 Mar 2025 13:43:19 -0600
Subject: [PATCH] Modify cfd to break up input products

---
 juno_samples/rodinia/cfd/src/lib.rs      | 229 +++++++++++++----------
 juno_samples/rodinia/cfd/src/rust_cfd.rs | 195 ++++++++++---------
 juno_samples/rodinia/cfd/src/setup.rs    | 140 ++++++++++----
 3 files changed, 335 insertions(+), 229 deletions(-)

diff --git a/juno_samples/rodinia/cfd/src/lib.rs b/juno_samples/rodinia/cfd/src/lib.rs
index e054a605..62ee59f4 100644
--- a/juno_samples/rodinia/cfd/src/lib.rs
+++ b/juno_samples/rodinia/cfd/src/lib.rs
@@ -24,144 +24,180 @@ pub struct CFDInputs {
 fn run_euler(
     nelr: usize,
     iterations: usize,
-    mut variables: AlignedSlice<f32>,
+    mut variables: Variables,
     areas: &[f32],
     elements_surrounding_elements: &[i32],
-    normals: &[f32],
-    ff_variable: &[f32],
+    normals: &Normals,
+    ff_variable: &Variable,
     ff_fc_density_energy: &Float3,
     ff_fc_momentum_x: &Float3,
     ff_fc_momentum_y: &Float3,
     ff_fc_momentum_z: &Float3,
-) -> Vec<f32> {
-    let mut variables = HerculesMutBox::from(variables.as_mut_slice());
+) -> Variables {
+    let mut v_density = HerculesMutBox::from(variables.density.as_mut_slice());
+    let mut v_momentum_x = HerculesMutBox::from(variables.momentum.x.as_mut_slice());
+    let mut v_momentum_y = HerculesMutBox::from(variables.momentum.y.as_mut_slice());
+    let mut v_momentum_z = HerculesMutBox::from(variables.momentum.z.as_mut_slice());
+    let mut v_energy = HerculesMutBox::from(variables.energy.as_mut_slice());
+
     let areas = HerculesImmBox::from(areas);
     let elements_surrounding_elements = HerculesImmBox::from(elements_surrounding_elements);
-    let normals = HerculesImmBox::from(normals);
-    let ff_variable = HerculesImmBox::from(ff_variable);
 
-    // TODO: Make hercules box handle structs, for now we'll copy into a vec
-    let ff_fc_density_energy = vec![
-        ff_fc_density_energy.x,
-        ff_fc_density_energy.y,
-        ff_fc_density_energy.z,
-    ];
-    let ff_fc_density_energy = HerculesImmBox::from(ff_fc_density_energy.as_slice());
-    let ff_fc_momentum_x = vec![ff_fc_momentum_x.x, ff_fc_momentum_x.y, ff_fc_momentum_x.z];
-    let ff_fc_momentum_x = HerculesImmBox::from(ff_fc_momentum_x.as_slice());
-    let ff_fc_momentum_y = vec![ff_fc_momentum_y.x, ff_fc_momentum_y.y, ff_fc_momentum_y.z];
-    let ff_fc_momentum_y = HerculesImmBox::from(ff_fc_momentum_y.as_slice());
-    let ff_fc_momentum_z = vec![ff_fc_momentum_z.x, ff_fc_momentum_z.y, ff_fc_momentum_z.z];
-    let ff_fc_momentum_z = HerculesImmBox::from(ff_fc_momentum_z.as_slice());
+    let normals_x = HerculesImmBox::from(normals.x.as_slice());
+    let normals_y = HerculesImmBox::from(normals.y.as_slice());
+    let normals_z = HerculesImmBox::from(normals.z.as_slice());
 
     let mut runner = runner!(euler);
-    let (density, x_m, y_m, z_m, energy) = async_std::task::block_on(async {
+    let (density, momentum_x, momentum_y, momentum_z, energy) = 
+    async_std::task::block_on(async {
         runner
             .run(
                 nelr as u64,
                 iterations as u64,
-                variables.to(),
+                v_density.to(),
+                v_momentum_x.to(),
+                v_momentum_y.to(),
+                v_momentum_z.to(),
+                v_energy.to(),
                 areas.to(),
                 elements_surrounding_elements.to(),
-                normals.to(),
-                ff_variable.to(),
-                ff_fc_density_energy.to(),
-                ff_fc_momentum_x.to(),
-                ff_fc_momentum_y.to(),
-                ff_fc_momentum_z.to(),
+                normals_x.to(),
+                normals_y.to(),
+                normals_z.to(),
+                ff_variable.density,
+                ff_variable.momentum.x,
+                ff_variable.momentum.y,
+                ff_variable.momentum.z,
+                ff_variable.energy,
+                ff_fc_density_energy.x,
+                ff_fc_density_energy.y,
+                ff_fc_density_energy.z,
+                ff_fc_momentum_x.x,
+                ff_fc_momentum_x.y,
+                ff_fc_momentum_x.z,
+                ff_fc_momentum_y.x,
+                ff_fc_momentum_y.y,
+                ff_fc_momentum_y.z,
+                ff_fc_momentum_z.x,
+                ff_fc_momentum_z.y,
+                ff_fc_momentum_z.z,
             )
             .await
     });
 
-    let total = vec![];
-    let density = HerculesMutBox::from(density).as_slice();
-    let x_m = HerculesMutBox::from(x_m).as_slice();
-    let y_m = HerculesMutBox::from(y_m).as_slice();
-    let z_m = HerculesMutBox::from(z_m).as_slice();
-    let energy = HerculesMutBox::from(energy).as_slice();
-    total.extend(density.into_iter().map(|id: &mut f32| *id));
-    total.extend(x_m.into_iter().map(|id: &mut f32| *id));
-    total.extend(y_m.into_iter().map(|id: &mut f32| *id));
-    total.extend(z_m.into_iter().map(|id: &mut f32| *id));
-    total.extend(energy.into_iter().map(|id: &mut f32| *id));
-    total
+    Variables {
+        density: AlignedSlice::from_slice(HerculesMutBox::from(density).as_slice()),
+        momentum: Momentum {
+            x: AlignedSlice::from_slice(HerculesMutBox::from(momentum_x).as_slice()),
+            y: AlignedSlice::from_slice(HerculesMutBox::from(momentum_y).as_slice()),
+            z: AlignedSlice::from_slice(HerculesMutBox::from(momentum_z).as_slice()),
+        },
+        energy: AlignedSlice::from_slice(HerculesMutBox::from(energy).as_slice()),
+    }
 }
 
 fn run_pre_euler(
     nelr: usize,
     iterations: usize,
-    mut variables: AlignedSlice<f32>,
+    mut variables: Variables,
     areas: &[f32],
     elements_surrounding_elements: &[i32],
-    normals: &[f32],
-    ff_variable: &[f32],
+    normals: &Normals,
+    ff_variable: &Variable,
     ff_fc_density_energy: &Float3,
     ff_fc_momentum_x: &Float3,
     ff_fc_momentum_y: &Float3,
     ff_fc_momentum_z: &Float3,
-) -> Vec<f32> {
-    let mut variables = HerculesMutBox::from(variables.as_mut_slice());
+) -> Variables {
+    let mut v_density = HerculesMutBox::from(variables.density.as_mut_slice());
+    let mut v_momentum_x = HerculesMutBox::from(variables.momentum.x.as_mut_slice());
+    let mut v_momentum_y = HerculesMutBox::from(variables.momentum.y.as_mut_slice());
+    let mut v_momentum_z = HerculesMutBox::from(variables.momentum.z.as_mut_slice());
+    let mut v_energy = HerculesMutBox::from(variables.energy.as_mut_slice());
+
     let areas = HerculesImmBox::from(areas);
     let elements_surrounding_elements = HerculesImmBox::from(elements_surrounding_elements);
-    let normals = HerculesImmBox::from(normals);
-    let ff_variable = HerculesImmBox::from(ff_variable);
 
-    // TODO: Make hercules box handle structs, for now we'll copy into a vec
-    let ff_fc_density_energy = vec![
-        ff_fc_density_energy.x,
-        ff_fc_density_energy.y,
-        ff_fc_density_energy.z,
-    ];
-    let ff_fc_density_energy = HerculesImmBox::from(ff_fc_density_energy.as_slice());
-    let ff_fc_momentum_x = vec![ff_fc_momentum_x.x, ff_fc_momentum_x.y, ff_fc_momentum_x.z];
-    let ff_fc_momentum_x = HerculesImmBox::from(ff_fc_momentum_x.as_slice());
-    let ff_fc_momentum_y = vec![ff_fc_momentum_y.x, ff_fc_momentum_y.y, ff_fc_momentum_y.z];
-    let ff_fc_momentum_y = HerculesImmBox::from(ff_fc_momentum_y.as_slice());
-    let ff_fc_momentum_z = vec![ff_fc_momentum_z.x, ff_fc_momentum_z.y, ff_fc_momentum_z.z];
-    let ff_fc_momentum_z = HerculesImmBox::from(ff_fc_momentum_z.as_slice());
+    let normals_x = HerculesImmBox::from(normals.x.as_slice());
+    let normals_y = HerculesImmBox::from(normals.y.as_slice());
+    let normals_z = HerculesImmBox::from(normals.z.as_slice());
 
     let mut runner = runner!(pre_euler);
-
-    let variables = variables.to();
-
-    let (density, x_m, y_m, z_m, energy) = async_std::task::block_on(async {
+    let (density, momentum_x, momentum_y, momentum_z, energy) = 
+    async_std::task::block_on(async {
         runner
             .run(
                 nelr as u64,
                 iterations as u64,
-                variables,
+                v_density.to(),
+                v_momentum_x.to(),
+                v_momentum_y.to(),
+                v_momentum_z.to(),
+                v_energy.to(),
                 areas.to(),
                 elements_surrounding_elements.to(),
-                normals.to(),
-                ff_variable.to(),
-                ff_fc_density_energy.to(),
-                ff_fc_momentum_x.to(),
-                ff_fc_momentum_y.to(),
-                ff_fc_momentum_z.to(),
+                normals_x.to(),
+                normals_y.to(),
+                normals_z.to(),
+                ff_variable.density,
+                ff_variable.momentum.x,
+                ff_variable.momentum.y,
+                ff_variable.momentum.z,
+                ff_variable.energy,
+                ff_fc_density_energy.x,
+                ff_fc_density_energy.y,
+                ff_fc_density_energy.z,
+                ff_fc_momentum_x.x,
+                ff_fc_momentum_x.y,
+                ff_fc_momentum_x.z,
+                ff_fc_momentum_y.x,
+                ff_fc_momentum_y.y,
+                ff_fc_momentum_y.z,
+                ff_fc_momentum_z.x,
+                ff_fc_momentum_z.y,
+                ff_fc_momentum_z.z,
             )
             .await
     });
 
-    let total = vec![];
-    let density = HerculesMutBox::from(density).as_slice();
-    let x_m = HerculesMutBox::from(x_m).as_slice();
-    let y_m = HerculesMutBox::from(y_m).as_slice();
-    let z_m = HerculesMutBox::from(z_m).as_slice();
-    let energy = HerculesMutBox::from(energy).as_slice();
-    total.extend(density.into_iter().map(|id: &mut f32| *id));
-    total.extend(x_m.into_iter().map(|id: &mut f32| *id));
-    total.extend(y_m.into_iter().map(|id: &mut f32| *id));
-    total.extend(z_m.into_iter().map(|id: &mut f32| *id));
-    total.extend(energy.into_iter().map(|id: &mut f32| *id));
-    total
+    Variables {
+        density: AlignedSlice::from_slice(HerculesMutBox::from(density).as_slice()),
+        momentum: Momentum {
+            x: AlignedSlice::from_slice(HerculesMutBox::from(momentum_x).as_slice()),
+            y: AlignedSlice::from_slice(HerculesMutBox::from(momentum_y).as_slice()),
+            z: AlignedSlice::from_slice(HerculesMutBox::from(momentum_z).as_slice()),
+        },
+        energy: AlignedSlice::from_slice(HerculesMutBox::from(energy).as_slice()),
+    }
 }
 
 fn compare_float(x: f32, y: f32) -> bool {
     (x - y).abs() < 1e-5
 }
 
-fn compare_floats(xs: &[f32], ys: &[f32]) -> bool {
-    xs.len() == ys.len() && xs.iter().zip(ys.iter()).all(|(x, y)| compare_float(*x, *y))
+fn compare_floats(xs: &Variables, ys: &Variables) -> bool {
+    let xs_density = xs.density.as_slice();
+    let xs_momentum_x = xs.momentum.x.as_slice();
+    let xs_momentum_y = xs.momentum.y.as_slice();
+    let xs_momentum_z = xs.momentum.z.as_slice();
+    let xs_energy = xs.energy.as_slice();
+
+    let ys_density = ys.density.as_slice();
+    let ys_momentum_x = ys.momentum.x.as_slice();
+    let ys_momentum_y = ys.momentum.y.as_slice();
+    let ys_momentum_z = ys.momentum.z.as_slice();
+    let ys_energy = ys.energy.as_slice();
+
+    xs_density.len() == ys_density.len()
+        && xs_density.iter().zip(ys_density.iter()).all(|(x, y)| compare_float(*x, *y))
+        && xs_momentum_x.len() == ys_momentum_x.len()
+        && xs_momentum_x.iter().zip(ys_momentum_x.iter()).all(|(x, y)| compare_float(*x, *y))
+        && xs_momentum_y.len() == ys_momentum_y.len()
+        && xs_momentum_y.iter().zip(ys_momentum_y.iter()).all(|(x, y)| compare_float(*x, *y))
+        && xs_momentum_z.len() == ys_momentum_z.len()
+        && xs_momentum_z.iter().zip(ys_momentum_z.iter()).all(|(x, y)| compare_float(*x, *y))
+        && xs_energy.len() == ys_energy.len()
+        && xs_energy.iter().zip(ys_energy.iter()).all(|(x, y)| compare_float(*x, *y))
 }
 
 pub fn cfd_harness(args: CFDInputs) {
@@ -172,8 +208,6 @@ pub fn cfd_harness(args: CFDInputs) {
         pre_euler,
     } = args;
 
-    assert!(block_size % 16 == 0, "Hercules expects all arrays to be 64-byte aligned, cfd uses structs of arrays that are annoying to deal with if the block_size is not a multiple of 16");
-
     let FarFieldConditions {
         ff_variable,
         ff_fc_momentum_x,
@@ -189,7 +223,7 @@ pub fn cfd_harness(args: CFDInputs) {
         normals,
     } = read_domain_geometry(data_file, block_size);
 
-    let variables = initialize_variables(nelr, ff_variable.as_slice());
+    let variables = initialize_variables(nelr, &ff_variable);
 
     let res_juno = if pre_euler {
         run_pre_euler(
@@ -198,8 +232,8 @@ pub fn cfd_harness(args: CFDInputs) {
             variables.clone(),
             areas.as_slice(),
             elements_surrounding_elements.as_slice(),
-            normals.as_slice(),
-            ff_variable.as_slice(),
+            &normals,
+            &ff_variable,
             &ff_fc_density_energy,
             &ff_fc_momentum_x,
             &ff_fc_momentum_y,
@@ -212,8 +246,8 @@ pub fn cfd_harness(args: CFDInputs) {
             variables.clone(),
             areas.as_slice(),
             elements_surrounding_elements.as_slice(),
-            normals.as_slice(),
-            ff_variable.as_slice(),
+            &normals,
+            &ff_variable,
             &ff_fc_density_energy,
             &ff_fc_momentum_x,
             &ff_fc_momentum_y,
@@ -227,8 +261,8 @@ pub fn cfd_harness(args: CFDInputs) {
             variables,
             areas.as_slice(),
             elements_surrounding_elements.as_slice(),
-            normals.as_slice(),
-            ff_variable.as_slice(),
+            &normals,
+            &ff_variable,
             &ff_fc_density_energy,
             &ff_fc_momentum_x,
             &ff_fc_momentum_y,
@@ -241,8 +275,8 @@ pub fn cfd_harness(args: CFDInputs) {
             variables,
             areas.as_slice(),
             elements_surrounding_elements.as_slice(),
-            normals.as_slice(),
-            ff_variable.as_slice(),
+            &normals,
+            &ff_variable,
             &ff_fc_density_energy,
             &ff_fc_momentum_x,
             &ff_fc_momentum_y,
@@ -250,8 +284,7 @@ pub fn cfd_harness(args: CFDInputs) {
         )
     };
 
-    if !compare_floats(&res_juno, res_rust.as_slice()) {
-        assert_eq!(res_juno.len(), res_rust.as_slice().len());
+    if !compare_floats(&res_juno, &res_rust) {
         panic!("Mismatch in results");
     }
 }
diff --git a/juno_samples/rodinia/cfd/src/rust_cfd.rs b/juno_samples/rodinia/cfd/src/rust_cfd.rs
index 42479df8..6d92c2b5 100644
--- a/juno_samples/rodinia/cfd/src/rust_cfd.rs
+++ b/juno_samples/rodinia/cfd/src/rust_cfd.rs
@@ -58,19 +58,19 @@ fn compute_speed_of_sound(density: f32, pressure: f32) -> f32 {
     (GAMMA * pressure / density).sqrt()
 }
 
-fn compute_step_factor(nelr: usize, variables: &[f32], areas: &[f32]) -> Vec<f32> {
+fn compute_step_factor(nelr: usize, variables: &Variables, areas: &[f32]) -> Vec<f32> {
     let mut step_factors = vec![0.0; nelr];
 
     for i in 0..nelr {
-        let density = variables[i + VAR_DENSITY * nelr];
+        let density = variables.density[i];
 
         let momentum = Float3 {
-            x: variables[i + (VAR_MOMENTUM + 0) * nelr],
-            y: variables[i + (VAR_MOMENTUM + 1) * nelr],
-            z: variables[i + (VAR_MOMENTUM + 2) * nelr],
+            x: variables.momentum.x[i],
+            y: variables.momentum.y[i],
+            z: variables.momentum.z[i],
         };
 
-        let density_energy = variables[i + VAR_DENSITY_ENERGY * nelr];
+        let density_energy = variables.energy[i];
 
         let velocity = compute_velocity(density, &momentum);
         let speed_sqd = compute_speed_sqd(&velocity);
@@ -85,28 +85,36 @@ fn compute_step_factor(nelr: usize, variables: &[f32], areas: &[f32]) -> Vec<f32
 
 fn compute_flux_euler(
     nelr: usize,
-    variables: &[f32],
+    variables: &Variables,
     elements_surrounding_elements: &[i32],
-    normals: &[f32],
-    ff_variable: &[f32],
+    normals: &Normals,
+    ff_variable: &Variable,
     ff_fc_momentum_x: &Float3,
     ff_fc_momentum_y: &Float3,
     ff_fc_momentum_z: &Float3,
     ff_fc_density_energy: &Float3,
-) -> Vec<f32> {
+) -> Variables {
     let smoothing_coefficient: f32 = 0.2;
 
-    let mut fluxes = vec![0.0; nelr * NVAR];
+    let mut fluxes = Variables {
+        density: AlignedSlice::of_len(nelr),
+        momentum: Momentum {
+            x: AlignedSlice::of_len(nelr),
+            y: AlignedSlice::of_len(nelr),
+            z: AlignedSlice::of_len(nelr),
+        },
+        energy: AlignedSlice::of_len(nelr),
+    };
 
     for i in 0..nelr {
-        let density_i = variables[i + VAR_DENSITY * nelr];
+        let density_i = variables.density[i];
         let momentum_i = Float3 {
-            x: variables[i + (VAR_MOMENTUM + 0) * nelr],
-            y: variables[i + (VAR_MOMENTUM + 1) * nelr],
-            z: variables[i + (VAR_MOMENTUM + 2) * nelr],
+            x: variables.momentum.x[i],
+            y: variables.momentum.y[i],
+            z: variables.momentum.z[i],
         };
 
-        let density_energy_i = variables[i + VAR_DENSITY_ENERGY * nelr];
+        let density_energy_i = variables.energy[i];
 
         let velocity_i = compute_velocity(density_i, &momentum_i);
         let speed_sqd_i = compute_speed_sqd(&velocity_i);
@@ -134,9 +142,9 @@ fn compute_flux_euler(
         for j in 0..NNB {
             let nb = elements_surrounding_elements[i + j * nelr];
             let normal = Float3 {
-                x: normals[i + (j + 0 * NNB) * nelr],
-                y: normals[i + (j + 1 * NNB) * nelr],
-                z: normals[i + (j + 2 * NNB) * nelr],
+                x: normals.x[i + j * nelr],
+                y: normals.y[i + j * nelr],
+                z: normals.z[i + j * nelr],
             };
             let normal_len =
                 (normal.x * normal.x + normal.y * normal.y + normal.z * normal.z).sqrt();
@@ -144,13 +152,13 @@ fn compute_flux_euler(
             if nb >= 0 {
                 // legitimate neighbor
                 let nb = nb as usize;
-                let density_nb = variables[nb + VAR_DENSITY * nelr];
+                let density_nb = variables.density[nb];
                 let momentum_nb = Float3 {
-                    x: variables[nb + (VAR_MOMENTUM + 0) * nelr],
-                    y: variables[nb + (VAR_MOMENTUM + 1) * nelr],
-                    z: variables[nb + (VAR_MOMENTUM + 2) * nelr],
+                    x: variables.momentum.x[nb],
+                    y: variables.momentum.y[nb],
+                    z: variables.momentum.z[nb],
                 };
-                let density_energy_nb = variables[nb + VAR_DENSITY_ENERGY * nelr];
+                let density_energy_nb = variables.energy[nb];
                 let velocity_nb = compute_velocity(density_nb, &momentum_nb);
                 let speed_sqd_nb = compute_speed_sqd(&velocity_nb);
                 let pressure_nb = compute_pressure(density_nb, density_energy_nb, speed_sqd_nb);
@@ -205,21 +213,21 @@ fn compute_flux_euler(
             } else if nb == -2 {
                 // a far field boundary
                 let factor = 0.5 * normal.x;
-                flux_i_density += factor * (ff_variable[VAR_MOMENTUM + 0] + momentum_i.x);
+                flux_i_density += factor * (ff_variable.momentum.x + momentum_i.x);
                 flux_i_density_energy += factor * (ff_fc_density_energy.x + fc_i_density_energy.x);
                 flux_i_momentum.x += factor * (ff_fc_momentum_x.x + fc_i_momentum_x.x);
                 flux_i_momentum.y += factor * (ff_fc_momentum_y.x + fc_i_momentum_y.x);
                 flux_i_momentum.z += factor * (ff_fc_momentum_z.x + fc_i_momentum_z.x);
 
                 let factor = 0.5 * normal.y;
-                flux_i_density += factor * (ff_variable[VAR_MOMENTUM + 1] + momentum_i.y);
+                flux_i_density += factor * (ff_variable.momentum.y + momentum_i.y);
                 flux_i_density_energy += factor * (ff_fc_density_energy.y + fc_i_density_energy.y);
                 flux_i_momentum.x += factor * (ff_fc_momentum_x.y + fc_i_momentum_x.y);
                 flux_i_momentum.y += factor * (ff_fc_momentum_y.y + fc_i_momentum_y.y);
                 flux_i_momentum.z += factor * (ff_fc_momentum_z.y + fc_i_momentum_z.y);
 
                 let factor = 0.5 * normal.z;
-                flux_i_density += factor * (ff_variable[VAR_MOMENTUM + 2] + momentum_i.z);
+                flux_i_density += factor * (ff_variable.momentum.z + momentum_i.z);
                 flux_i_density_energy += factor * (ff_fc_density_energy.z + fc_i_density_energy.z);
                 flux_i_momentum.x += factor * (ff_fc_momentum_x.z + fc_i_momentum_x.z);
                 flux_i_momentum.y += factor * (ff_fc_momentum_y.z + fc_i_momentum_y.z);
@@ -227,11 +235,11 @@ fn compute_flux_euler(
             }
         }
 
-        fluxes[i + VAR_DENSITY * nelr] = flux_i_density;
-        fluxes[i + (VAR_MOMENTUM + 0) * nelr] = flux_i_momentum.x;
-        fluxes[i + (VAR_MOMENTUM + 1) * nelr] = flux_i_momentum.y;
-        fluxes[i + (VAR_MOMENTUM + 2) * nelr] = flux_i_momentum.z;
-        fluxes[i + VAR_DENSITY_ENERGY * nelr] = flux_i_density_energy;
+        fluxes.density[i] = flux_i_density;
+        fluxes.momentum.x[i] = flux_i_momentum.x;
+        fluxes.momentum.y[i] = flux_i_momentum.y;
+        fluxes.momentum.z[i] = flux_i_momentum.z;
+        fluxes.energy[i] = flux_i_density_energy;
     }
 
     fluxes
@@ -239,7 +247,7 @@ fn compute_flux_euler(
 
 fn compute_flux_contributions(
     nelr: usize,
-    variables: &[f32],
+    variables: &Variables,
 ) -> (Vec<f32>, Vec<f32>, Vec<f32>, Vec<f32>) {
     let mut fc_momentum_x = vec![0.0; nelr * NDIM];
     let mut fc_momentum_y = vec![0.0; nelr * NDIM];
@@ -247,14 +255,14 @@ fn compute_flux_contributions(
     let mut fc_density_energy = vec![0.0; nelr * NDIM];
 
     for i in 0..nelr {
-        let density_i = variables[i + VAR_DENSITY * nelr];
+        let density_i = variables.density[i];
         let momentum_i = Float3 {
-            x: variables[i + (VAR_MOMENTUM + 0) * nelr],
-            y: variables[i + (VAR_MOMENTUM + 1) * nelr],
-            z: variables[i + (VAR_MOMENTUM + 2) * nelr],
+            x: variables.momentum.x[i],
+            y: variables.momentum.y[i],
+            z: variables.momentum.z[i],
         };
 
-        let density_energy_i = variables[i + VAR_DENSITY_ENERGY * nelr];
+        let density_energy_i = variables.energy[i];
 
         let velocity_i = compute_velocity(density_i, &momentum_i);
         let speed_sqd_i = compute_speed_sqd(&velocity_i);
@@ -298,32 +306,40 @@ fn compute_flux_contributions(
 
 fn compute_flux_pre_euler(
     nelr: usize,
-    variables: &[f32],
+    variables: &Variables,
     elements_surrounding_elements: &[i32],
-    normals: &[f32],
+    normals: &Normals,
     fc_momentum_x: &[f32],
     fc_momentum_y: &[f32],
     fc_momentum_z: &[f32],
     fc_density_energy: &[f32],
-    ff_variable: &[f32],
+    ff_variable: &Variable,
     ff_fc_momentum_x: &Float3,
     ff_fc_momentum_y: &Float3,
     ff_fc_momentum_z: &Float3,
     ff_fc_density_energy: &Float3,
-) -> Vec<f32> {
+) -> Variables {
     let smoothing_coefficient: f32 = 0.2;
 
-    let mut fluxes = vec![0.0; nelr * NVAR];
+    let mut fluxes = Variables {
+        density: AlignedSlice::of_len(nelr),
+        momentum: Momentum {
+            x: AlignedSlice::of_len(nelr),
+            y: AlignedSlice::of_len(nelr),
+            z: AlignedSlice::of_len(nelr),
+        },
+        energy: AlignedSlice::of_len(nelr),
+    };
 
     for i in 0..nelr {
-        let density_i = variables[i + VAR_DENSITY * nelr];
+        let density_i = variables.density[i];
         let momentum_i = Float3 {
-            x: variables[i + (VAR_MOMENTUM + 0) * nelr],
-            y: variables[i + (VAR_MOMENTUM + 1) * nelr],
-            z: variables[i + (VAR_MOMENTUM + 2) * nelr],
+            x: variables.momentum.x[i],
+            y: variables.momentum.y[i],
+            z: variables.momentum.z[i],
         };
 
-        let density_energy_i = variables[i + VAR_DENSITY_ENERGY * nelr];
+        let density_energy_i = variables.energy[i];
 
         let velocity_i = compute_velocity(density_i, &momentum_i);
         let speed_sqd_i = compute_speed_sqd(&velocity_i);
@@ -363,22 +379,22 @@ fn compute_flux_pre_euler(
         for j in 0..NNB {
             let nb = elements_surrounding_elements[i + j * nelr];
             let normal = Float3 {
-                x: normals[i + (j + 0 * NNB) * nelr],
-                y: normals[i + (j + 1 * NNB) * nelr],
-                z: normals[i + (j + 2 * NNB) * nelr],
+                x: normals.x[i + j * nelr],
+                y: normals.y[i + j * nelr],
+                z: normals.z[i + j * nelr],
             };
             let normal_len =
                 (normal.x * normal.x + normal.y * normal.y + normal.z * normal.z).sqrt();
 
             if nb >= 0 {
                 let nb = nb as usize;
-                let density_nb = variables[nb + VAR_DENSITY * nelr];
+                let density_nb = variables.density[nb];
                 let momentum_nb = Float3 {
-                    x: variables[nb + (VAR_MOMENTUM + 0) * nelr],
-                    y: variables[nb + (VAR_MOMENTUM + 1) * nelr],
-                    z: variables[nb + (VAR_MOMENTUM + 2) * nelr],
+                    x: variables.momentum.x[nb],
+                    y: variables.momentum.y[nb],
+                    z: variables.momentum.z[nb],
                 };
-                let density_energy_nb = variables[nb + VAR_DENSITY_ENERGY * nelr];
+                let density_energy_nb = variables.energy[nb];
                 let velocity_nb = compute_velocity(density_nb, &momentum_nb);
                 let speed_sqd_nb = compute_speed_sqd(&velocity_nb);
                 let pressure_nb = compute_pressure(density_nb, density_energy_nb, speed_sqd_nb);
@@ -443,21 +459,21 @@ fn compute_flux_pre_euler(
                 flux_i_momentum.z += normal.z * pressure_i;
             } else if nb == -2 {
                 let factor = 0.5 * normal.x;
-                flux_i_density += factor * (ff_variable[VAR_MOMENTUM + 0] + momentum_i.x);
+                flux_i_density += factor * (ff_variable.momentum.x + momentum_i.x);
                 flux_i_density_energy += factor * (ff_fc_density_energy.x + fc_i_density_energy.x);
                 flux_i_momentum.x += factor * (ff_fc_momentum_x.x + fc_i_momentum_x.x);
                 flux_i_momentum.y += factor * (ff_fc_momentum_y.x + fc_i_momentum_y.x);
                 flux_i_momentum.z += factor * (ff_fc_momentum_z.x + fc_i_momentum_z.x);
 
                 let factor = 0.5 * normal.y;
-                flux_i_density += factor * (ff_variable[VAR_MOMENTUM + 1] + momentum_i.y);
+                flux_i_density += factor * (ff_variable.momentum.y + momentum_i.y);
                 flux_i_density_energy += factor * (ff_fc_density_energy.y + fc_i_density_energy.y);
                 flux_i_momentum.x += factor * (ff_fc_momentum_x.y + fc_i_momentum_x.y);
                 flux_i_momentum.y += factor * (ff_fc_momentum_y.y + fc_i_momentum_y.y);
                 flux_i_momentum.z += factor * (ff_fc_momentum_z.y + fc_i_momentum_z.y);
 
                 let factor = 0.5 * normal.z;
-                flux_i_density += factor * (ff_variable[VAR_MOMENTUM + 2] + momentum_i.z);
+                flux_i_density += factor * (ff_variable.momentum.z + momentum_i.z);
                 flux_i_density_energy += factor * (ff_fc_density_energy.z + fc_i_density_energy.z);
                 flux_i_momentum.x += factor * (ff_fc_momentum_x.z + fc_i_momentum_x.z);
                 flux_i_momentum.y += factor * (ff_fc_momentum_y.z + fc_i_momentum_y.z);
@@ -465,11 +481,11 @@ fn compute_flux_pre_euler(
             }
         }
 
-        fluxes[i + VAR_DENSITY * nelr] = flux_i_density;
-        fluxes[i + (VAR_MOMENTUM + 0) * nelr] = flux_i_momentum.x;
-        fluxes[i + (VAR_MOMENTUM + 1) * nelr] = flux_i_momentum.y;
-        fluxes[i + (VAR_MOMENTUM + 2) * nelr] = flux_i_momentum.z;
-        fluxes[i + VAR_DENSITY_ENERGY * nelr] = flux_i_density_energy;
+        fluxes.density[i] = flux_i_density;
+        fluxes.momentum.x[i] = flux_i_momentum.x;
+        fluxes.momentum.y[i] = flux_i_momentum.y;
+        fluxes.momentum.z[i] = flux_i_momentum.z;
+        fluxes.energy[i] = flux_i_density_energy;
     }
 
     fluxes
@@ -478,47 +494,42 @@ fn compute_flux_pre_euler(
 fn time_step(
     j: usize,
     nelr: usize,
-    old_variables: &[f32],
-    variables: &mut AlignedSlice<f32>,
+    old_variables: &Variables,
+    variables: &mut Variables,
     step_factors: &[f32],
-    fluxes: &[f32],
+    fluxes: &Variables,
 ) {
     for i in 0..nelr {
         let factor = step_factors[i] / (RK + 1 - j) as f32;
-        variables[i + VAR_DENSITY * nelr] =
-            old_variables[i + VAR_DENSITY * nelr] + factor * fluxes[i + VAR_DENSITY * nelr];
-        variables[i + (VAR_MOMENTUM + 0) * nelr] = old_variables[i + (VAR_MOMENTUM + 0) * nelr]
-            + factor * fluxes[i + (VAR_MOMENTUM + 0) * nelr];
-        variables[i + (VAR_MOMENTUM + 1) * nelr] = old_variables[i + (VAR_MOMENTUM + 1) * nelr]
-            + factor * fluxes[i + (VAR_MOMENTUM + 1) * nelr];
-        variables[i + (VAR_MOMENTUM + 2) * nelr] = old_variables[i + (VAR_MOMENTUM + 2) * nelr]
-            + factor * fluxes[i + (VAR_MOMENTUM + 2) * nelr];
-        variables[i + VAR_DENSITY_ENERGY * nelr] = old_variables[i + VAR_DENSITY_ENERGY * nelr]
-            + factor * fluxes[i + VAR_DENSITY_ENERGY * nelr];
+        variables.density[i] = old_variables.density[i] + factor * fluxes.density[i];
+        variables.momentum.x[i] = old_variables.momentum.x[i] + factor * fluxes.momentum.x[i];
+        variables.momentum.y[i] = old_variables.momentum.y[i] + factor * fluxes.momentum.y[i];
+        variables.momentum.z[i] = old_variables.momentum.z[i] + factor * fluxes.momentum.z[i];
+        variables.energy[i] = old_variables.energy[i] + factor * fluxes.energy[i];
     }
 }
 
 pub fn euler(
     nelr: usize,
     iterations: usize,
-    mut variables: AlignedSlice<f32>,
+    mut variables: Variables,
     areas: &[f32],
     elements_surrounding_elements: &[i32],
-    normals: &[f32],
-    ff_variable: &[f32],
+    normals: &Normals,
+    ff_variable: &Variable,
     ff_fc_density_energy: &Float3,
     ff_fc_momentum_x: &Float3,
     ff_fc_momentum_y: &Float3,
     ff_fc_momentum_z: &Float3,
-) -> AlignedSlice<f32> {
+) -> Variables {
     for i in 0..iterations {
         let old_variables = variables.clone();
-        let step_factors = compute_step_factor(nelr, variables.as_slice(), areas);
+        let step_factors = compute_step_factor(nelr, &variables, areas);
 
         for j in 0..RK {
             let fluxes = compute_flux_euler(
                 nelr,
-                variables.as_slice(),
+                &variables,
                 elements_surrounding_elements,
                 normals,
                 ff_variable,
@@ -530,7 +541,7 @@ pub fn euler(
             time_step(
                 j,
                 nelr,
-                old_variables.as_slice(),
+                &old_variables,
                 &mut variables,
                 &step_factors,
                 &fluxes,
@@ -544,26 +555,26 @@ pub fn euler(
 pub fn pre_euler(
     nelr: usize,
     iterations: usize,
-    mut variables: AlignedSlice<f32>,
+    mut variables: Variables,
     areas: &[f32],
     elements_surrounding_elements: &[i32],
-    normals: &[f32],
-    ff_variable: &[f32],
+    normals: &Normals,
+    ff_variable: &Variable,
     ff_fc_density_energy: &Float3,
     ff_fc_momentum_x: &Float3,
     ff_fc_momentum_y: &Float3,
     ff_fc_momentum_z: &Float3,
-) -> AlignedSlice<f32> {
+) -> Variables {
     for i in 0..iterations {
         let old_variables = variables.clone();
-        let step_factors = compute_step_factor(nelr, variables.as_slice(), areas);
+        let step_factors = compute_step_factor(nelr, &variables, areas);
 
         for j in 0..RK {
             let (fc_momentum_x, fc_momentum_y, fc_momentum_z, fc_density_energy) =
-                compute_flux_contributions(nelr, variables.as_slice());
+                compute_flux_contributions(nelr, &variables);
             let fluxes = compute_flux_pre_euler(
                 nelr,
-                variables.as_slice(),
+                &variables,
                 elements_surrounding_elements,
                 normals,
                 &fc_momentum_x,
@@ -579,7 +590,7 @@ pub fn pre_euler(
             time_step(
                 j,
                 nelr,
-                old_variables.as_slice(),
+                &old_variables,
                 &mut variables,
                 &step_factors,
                 &fluxes,
diff --git a/juno_samples/rodinia/cfd/src/setup.rs b/juno_samples/rodinia/cfd/src/setup.rs
index b996f057..0da00251 100644
--- a/juno_samples/rodinia/cfd/src/setup.rs
+++ b/juno_samples/rodinia/cfd/src/setup.rs
@@ -1,10 +1,10 @@
+use std::fmt::{Debug, Error, Formatter};
 use std::fs::File;
 use std::io::Read;
+use std::ops::{Index, IndexMut};
 use std::path::Path;
 use std::str::FromStr;
 
-use std::ops::{Index, IndexMut};
-
 use nom::Parser;
 
 use crate::rust_cfd;
@@ -14,10 +14,6 @@ pub const ff_mach: f32 = 1.2;
 pub const NDIM: usize = 3;
 pub const NNB: usize = 4;
 pub const RK: usize = 3;
-pub const VAR_DENSITY: usize = 0;
-pub const VAR_MOMENTUM: usize = 1;
-pub const VAR_DENSITY_ENERGY: usize = VAR_MOMENTUM + NDIM;
-pub const NVAR: usize = VAR_DENSITY_ENERGY + 1;
 pub const deg_angle_of_attack: f32 = 0.0;
 
 #[repr(align(64))]
@@ -31,6 +27,15 @@ pub struct AlignedSlice<T> {
     length: usize,
 }
 
+impl<T> Debug for AlignedSlice<T>
+where
+    T: Debug
+{
+    fn fmt(&self, f: &mut Formatter<'_>) -> Result<(), Error> {
+        write!(f, "{:?}", self.as_slice())
+    }
+}
+
 impl<T> AlignedSlice<T>
 where
     T: Default,
@@ -53,6 +58,17 @@ where
     }
 }
 
+impl<T> AlignedSlice<T>
+where
+    T: Clone + Default,
+{
+    pub fn from_slice(slice: &[T]) -> Self {
+        let mut res = Self::of_len(slice.len());
+        res.as_mut_slice().clone_from_slice(slice);
+        res
+    }
+}
+
 impl<T> AlignedSlice<T> {
     pub fn as_slice(&self) -> &[T] {
         &self.slice[self.offset..self.offset + self.length]
@@ -90,15 +106,20 @@ where
     }
 }
 
-#[repr(C)]
 pub struct Float3 {
     pub x: f32,
     pub y: f32,
     pub z: f32,
 }
 
+pub struct Variable {
+    pub density: f32,
+    pub momentum: Float3,
+    pub energy: f32,
+}
+
 pub struct FarFieldConditions {
-    pub ff_variable: AlignedSlice<f32>,
+    pub ff_variable: Variable,
     pub ff_fc_momentum_x: Float3,
     pub ff_fc_momentum_y: Float3,
     pub ff_fc_momentum_z: Float3,
@@ -106,14 +127,14 @@ pub struct FarFieldConditions {
 }
 
 pub fn set_far_field_conditions() -> FarFieldConditions {
-    let mut ff_variable: AlignedSlice<f32> = AlignedSlice::of_len(NVAR);
+    let mut ff_variable = Variable { density: 0.0, momentum: Float3 { x: 0.0, y: 0.0, z: 0.0 }, energy: 0.0 };
 
     let angle_of_attack = std::f32::consts::PI / 180.0 * deg_angle_of_attack;
 
-    ff_variable[VAR_DENSITY] = 1.4;
+    ff_variable.density = 1.4;
 
     let ff_pressure = 1.0;
-    let ff_speed_of_sound = (GAMMA * ff_pressure / ff_variable[VAR_DENSITY]).sqrt();
+    let ff_speed_of_sound = (GAMMA * ff_pressure / ff_variable.density).sqrt();
     let ff_speed = ff_mach * ff_speed_of_sound;
 
     let ff_velocity = Float3 {
@@ -122,23 +143,23 @@ pub fn set_far_field_conditions() -> FarFieldConditions {
         z: 0.0,
     };
 
-    ff_variable[VAR_MOMENTUM + 0] = ff_variable[VAR_DENSITY] * ff_velocity.x;
-    ff_variable[VAR_MOMENTUM + 1] = ff_variable[VAR_DENSITY] * ff_velocity.y;
-    ff_variable[VAR_MOMENTUM + 2] = ff_variable[VAR_DENSITY] * ff_velocity.z;
+    ff_variable.momentum.x = ff_variable.density * ff_velocity.x;
+    ff_variable.momentum.y = ff_variable.density * ff_velocity.y;
+    ff_variable.momentum.z = ff_variable.density * ff_velocity.z;
 
-    ff_variable[VAR_DENSITY_ENERGY] =
-        ff_variable[VAR_DENSITY] * (0.5 * (ff_speed * ff_speed)) + (ff_pressure / (GAMMA - 1.0));
+    ff_variable.energy =
+        ff_variable.density * (0.5 * (ff_speed * ff_speed)) + (ff_pressure / (GAMMA - 1.0));
 
     let ff_momentum = Float3 {
-        x: ff_variable[VAR_MOMENTUM + 0],
-        y: ff_variable[VAR_MOMENTUM + 1],
-        z: ff_variable[VAR_MOMENTUM + 2],
+        x: ff_variable.momentum.x,
+        y: ff_variable.momentum.y,
+        z: ff_variable.momentum.z,
     };
     let (ff_fc_momentum_x, ff_fc_momentum_y, ff_fc_momentum_z, ff_fc_density_energy) =
         rust_cfd::compute_flux_contribution(
-            ff_variable[VAR_DENSITY],
+            ff_variable.density,
             &ff_momentum,
-            ff_variable[VAR_DENSITY_ENERGY],
+            ff_variable.energy,
             ff_pressure,
             &ff_velocity,
         );
@@ -152,11 +173,17 @@ pub fn set_far_field_conditions() -> FarFieldConditions {
     }
 }
 
+pub struct Normals {
+    pub x: AlignedSlice<f32>,
+    pub y: AlignedSlice<f32>,
+    pub z: AlignedSlice<f32>,
+}
+
 pub struct GeometryData {
     pub nelr: usize,
     pub areas: AlignedSlice<f32>,
     pub elements_surrounding_elements: AlignedSlice<i32>,
-    pub normals: AlignedSlice<f32>,
+    pub normals: Normals,
 }
 
 pub fn read_domain_geometry<P>(path: P, block_size: usize) -> GeometryData
@@ -183,7 +210,11 @@ fn cfd_parser<'a>(block_size: usize, text: &'a str) -> nom::IResult<&'a str, Geo
 
     let mut areas: AlignedSlice<f32> = AlignedSlice::of_len(nelr);
     let mut elements_surrounding_elements: AlignedSlice<i32> = AlignedSlice::of_len(nelr * NNB);
-    let mut normals: AlignedSlice<f32> = AlignedSlice::of_len(NDIM * NNB * nelr);
+    let mut normals = Normals {
+        x: AlignedSlice::of_len(NNB * nelr),
+        y: AlignedSlice::of_len(NNB * nelr),
+        z: AlignedSlice::of_len(NNB * nelr),
+    };
 
     // read in data
     let mut text = text;
@@ -209,14 +240,21 @@ fn cfd_parser<'a>(block_size: usize, text: &'a str) -> nom::IResult<&'a str, Geo
             let val = i32::from_str(val).unwrap();
             let val = if neg { -val } else { val };
             elements_surrounding_elements[i + j * nelr] = if val < 0 { -1 } else { val } - 1; // it's coming in with Fortran numbering
-
-            for k in 0..NDIM {
-                let t = nom::character::complete::multispace0(text)?.0;
-                let (t, val) = nom::number::complete::float(t)?;
-                text = t;
-
-                normals[i + (j + k * NNB) * nelr] = -val;
-            }
+                
+            let t = nom::character::complete::multispace0(text)?.0;
+            let (t, val) = nom::number::complete::float(t)?;
+            text = t;
+            normals.x[i + j * nelr] = -val;
+                
+            let t = nom::character::complete::multispace0(text)?.0;
+            let (t, val) = nom::number::complete::float(t)?;
+            text = t;
+            normals.y[i + j * nelr] = -val;
+                
+            let t = nom::character::complete::multispace0(text)?.0;
+            let (t, val) = nom::number::complete::float(t)?;
+            text = t;
+            normals.z[i + j * nelr] = -val;
         }
     }
 
@@ -227,9 +265,9 @@ fn cfd_parser<'a>(block_size: usize, text: &'a str) -> nom::IResult<&'a str, Geo
         for j in 0..NNB {
             elements_surrounding_elements[i + j * nelr] =
                 elements_surrounding_elements[last + j * nelr];
-            for k in 0..NDIM {
-                normals[i + (j + k * NNB) * nelr] = normals[last + (j + k * NNB) * nelr];
-            }
+            normals.x[i + j * nelr] = normals.x[last + j * nelr];
+            normals.y[i + j * nelr] = normals.y[last + j * nelr];
+            normals.z[i + j * nelr] = normals.z[last + j * nelr];
         }
     }
 
@@ -245,13 +283,37 @@ fn cfd_parser<'a>(block_size: usize, text: &'a str) -> nom::IResult<&'a str, Geo
     ))
 }
 
-pub fn initialize_variables(nelr: usize, ff_variable: &[f32]) -> AlignedSlice<f32> {
-    let mut variables = AlignedSlice::of_len(nelr * NVAR);
+#[derive(Clone)]
+pub struct Momentum {
+    pub x: AlignedSlice<f32>,
+    pub y: AlignedSlice<f32>,
+    pub z: AlignedSlice<f32>,
+}
+
+#[derive(Clone)]
+pub struct Variables {
+    pub density: AlignedSlice<f32>,
+    pub momentum: Momentum,
+    pub energy: AlignedSlice<f32>,
+}
+
+pub fn initialize_variables(nelr: usize, ff_variable: &Variable) -> Variables {
+    let mut variables = Variables {
+        density: AlignedSlice::of_len(nelr),
+        momentum: Momentum {
+            x: AlignedSlice::of_len(nelr),
+            y: AlignedSlice::of_len(nelr),
+            z: AlignedSlice::of_len(nelr),
+        },
+        energy: AlignedSlice::of_len(nelr),
+    };
 
     for i in 0..nelr {
-        for j in 0..NVAR {
-            variables[i + j * nelr] = ff_variable[j];
-        }
+        variables.density[i] = ff_variable.density;
+        variables.momentum.x[i] = ff_variable.momentum.x;
+        variables.momentum.y[i] = ff_variable.momentum.y;
+        variables.momentum.z[i] = ff_variable.momentum.z;
+        variables.energy[i] = ff_variable.energy;
     }
 
     variables
-- 
GitLab