Start euler cfd

Cargo.lock

@@ -1175,6 +1175,18 @@ dependencies = [
  "with_builtin_macros",
 ]
 
+[[package]]
+name = "juno_cfd"
+version = "0.1.0"
+dependencies = [
+ "async-std",
+ "clap",
+ "hercules_rt",
+ "juno_build",
+ "nom 6.2.2",
+ "with_builtin_macros",
+]
+
 [[package]]
 name = "juno_concat"
 version = "0.1.0"

Cargo.toml

@@ -37,5 +37,6 @@ members = [
 	"juno_samples/median_window",
 	"juno_samples/rodinia/bfs",
 	"juno_samples/rodinia/backprop",
+	"juno_samples/rodinia/cfd",
 	"juno_samples/rodinia/srad",
 ]

juno_samples/rodinia/cfd/Cargo.toml

+[package]
+name = "juno_cfd"
+version = "0.1.0"
+authors = ["Aaron Councilman <aaronjc4@illinois.edu>"]
+edition = "2021"
+
+[[bin]]
+name = "juno_cfd"
+path = "src/main.rs"
+
+[features]
+cuda = ["juno_build/cuda", "hercules_rt/cuda"]
+
+[build-dependencies]
+juno_build = { path = "../../../juno_build" }
+
+[dependencies]
+juno_build = { path = "../../../juno_build" }
+hercules_rt = { path = "../../../hercules_rt" }
+async-std = "*"
+clap = { version = "*", features = ["derive"] }
+with_builtin_macros = "0.1.0"
+nom = "*"

juno_samples/rodinia/cfd/build.rs

+use juno_build::JunoCompiler;
+
+fn main() {
+    #[cfg(feature = "cuda")]
+    JunoCompiler::new()
+        .file_in_src("euler.jn")
+        .unwrap()
+        .schedule_in_src("gpu.sch")
+        .unwrap()
+        .build()
+        .unwrap();
+    #[cfg(not(feature = "cuda"))]
+    JunoCompiler::new()
+        .file_in_src("euler.jn")
+        .unwrap()
+        .build()
+        .unwrap();
+}

juno_samples/rodinia/cfd/src/euler.jn

+const NNB : usize = 4;
+
+type Normals<nelr: usize> = struct {
+  x: f32[NNB, nelr],
+  y: f32[NNB, nelr],
+  z: f32[NNB, nelr],
+};
+
+type Momentum<nelr: usize> = struct {
+  x: f32[nelr],
+  y: f32[nelr],
+  z: f32[nelr],
+};
+
+type Variables<nelr: usize> = struct {
+  density: f32[nelr],
+  momentum: Momentum::<nelr>,
+  energy: f32[nelr],
+};
+
+type float3 = struct { x: f32, y: f32, z: f32 };
+
+type Variable = struct {
+  density: f32,
+  momentum: float3,
+  energy: f32,
+};
+
+fn compute_velocity(density: f32, momentum: float3) -> float3 {
+  return float3 { x: momentum.x / density,
+                  y: momentum.y / density,
+                  z: momentum.z / density };
+}
+
+fn compute_speed_sqd(velocity: float3) -> f32 {
+  return velocity.x * velocity.x + velocity.y * velocity.y + velocity.z * velocity.z;
+}
+
+const GAMMA : f32 = 1.4;
+
+fn compute_pressure(density: f32, density_energy: f32, speed_sqd: f32) -> f32 {
+  return (GAMMA - 1.0) * (density_energy - 0.5 * density * speed_sqd);
+}
+
+fn compute_speed_of_sound(density: f32, pressure: f32) -> f32 {
+  return sqrt!(GAMMA * pressure / density);
+}
+
+fn compute_step_factor<nelr: usize>(variables: Variables::<nelr>, areas: f32[nelr]) -> f32[nelr] {
+  let step_factors : f32[nelr];
+
+  for i = 0 to nelr {
+    let density = variables.density[i];
+
+    let momentum : float3;
+    momentum.x = variables.momentum.x[i];
+    momentum.y = variables.momentum.y[i];
+    momentum.z = variables.momentum.z[i];
+
+    let density_energy = variables.energy[i];
+
+    let velocity       = compute_velocity(density, momentum);
+    let speed_sqd      = compute_speed_sqd(velocity);
+    let pressure       = compute_pressure(density, density_energy, speed_sqd);
+    let speed_of_sound = compute_speed_of_sound(density, pressure);
+
+    step_factors[i] = 0.5 / sqrt!(areas[i]) * (sqrt!(speed_sqd) + speed_of_sound);
+  }
+
+  return step_factors;
+}
+
+fn compute_flux_contribution(
+  density: f32,
+  momentum: float3,
+  density_energy: f32,
+  pressure: f32,
+  velocity: float3,
+) -> (float3, float3, float3, float3) {
+  let fc_momentum_x = float3 { x: velocity.x * momentum.x + pressure,
+                               y: velocity.x * momentum.y,
+                               z: velocity.x * momentum.z };
+  let fc_momentum_y = float3 { x: fc_momentum_x.y,
+                               y: velocity.y * momentum.y + pressure,
+                               z: velocity.y * momentum.z };
+  let fc_momentum_z = float3 { x: fc_momentum_x.z,
+                               y: fc_momentum_y.z,
+                               z: velocity.z * momentum.z + pressure };
+  
+  let de_p = density_energy + pressure;
+  let fc_density_energy = float3 { x: velocity.x * de_p,
+                                   y: velocity.y * de_p,
+                                   z: velocity.z * de_p };
+
+  return (fc_momentum_x, fc_momentum_y, fc_momentum_z, fc_density_energy);
+}
+
+fn compute_flux<nelr: usize>(
+  variables: Variables::<nelr>,
+  elements_surrounding_elements: i32[NNB, nelr],
+  normals: Normals::<nelr>,
+  ff_variable: Variable,
+  ff_flux_contribution_density_energy: float3,
+  ff_flux_contribution_momentum_x: float3,
+  ff_flux_contribution_momentum_y: float3,
+  ff_flux_contribution_momentum_z: float3,
+) -> Variables::<nelr> {
+  const smoothing_coefficient : f32 = 0.2;
+  let fluxes: Variables::<nelr>;
+
+  for i = 0 to nelr {
+    let density_i = variables.density[i];
+
+    let momentum_i = float3 { x: variables.momentum.x[i],
+                              y: variables.momentum.y[i],
+                              z: variables.momentum.z[i] };
+
+    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);
+    let speed_i          = sqrt!(speed_sqd_i);
+    let pressure_i       = compute_pressure(density_i, density_energy_i, speed_sqd_i);
+    let speed_of_sound_i = compute_speed_of_sound(density_i, pressure_i);
+
+    let flux_contribution_i_momentum_x : float3;
+    let flux_contribution_i_momentum_y : float3;
+    let flux_contribution_i_momentum_z : float3;
+    let flux_contribution_i_density_energy: f32;
+    let (flux_contribution_i_momentum_x, flux_contribution_i_momentum_y,
+         flux_contribution_i_momentum_z, flux_contribution_i_density_energy)
+      = compute_flux_contribution(density_i, momentum_i, density_energy_i, pressure_i, velocity_i);
+
+    let flux_i_density : f32 = 0;
+    let flux_i_momentum = float3 { x: 0.0, y: 0.0, z: 0.0 };
+    let flux_i_density_energy : f32 = 0.0;
+
+    for j = 0 to NNB {
+      let nb = elements_surrounding_elements[j, i];
+      let normal = float3 {
+        x: normals.x[j, i],
+        y: normals.y[j, i],
+        z: normals.z[j, i],
+      };
+      let normal_len = sqrt!(normal.x*normal.x + normal.y*normal.y + normal.z*normal.z);
+
+      if nb >= 0 { // a legitimate neighbor
+        let nb = nb as usize;
+        let density_nb = variables.density[nb];
+        let momentum_nb = float3 {
+          x: variables.momentum.x[nb],
+          y: variables.momentum.y[nb],
+          z: variables.momentum.z[nb],
+        };
+        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);
+        let speed_of_sound_nb = compute_speed_of_sound(density_nb, pressure_nb);
+
+        let (flux_contribution_nb_momentum_x, flux_contribution_nb_momentum_y,
+             flux_contribution_nb_momentum_z, flux_contribution_nb_density_energy)
+          = compute_flux_contribution(density_nb, momentum_nb, density_energy_nb, pressure_nb, velocity_nb);
+        
+        // artificial viscosity
+        let factor = -normal_len * smoothing_coefficient * 0.5
+                   * (speed_i + sqrt!(speed_sqd_nb) + speed_of_sound_i + speed_of_sound_nb);
+        flux_i_density += factor * (density_i - density_nb);
+        flux_i_density_energy += factor * (density_energy_i - density_energy_nb);
+        flux_i_momentum.x += factor * (momentum_i.x - momentum_nb.x);
+        flux_i_momentum.y += factor * (momentum_i.y - momentum_nb.y);
+        flux_i_momentum.z += factor * (momentum_i.z - momentum_nb.z);
+
+        // accumulate cell-centered fluxes
+        let factor = 0.5 * normal.x;
+        flux_i_density += factor * (momentum_nb.x + momentum_i.x);
+        flux_i_density_energy += factor * (flux_contribution_nb_density_energy.x + flux_contribution_i_density_energy.x);
+        flux_i_momentum.x += factor * (flux_contribution_nb_momentum_x.x + flux_contribution_i_momentum_x.x);
+        flux_i_momentum.y += factor * (flux_contribution_nb_momentum_y.x + flux_contribution_i_momentum_y.x);
+        flux_i_momentum.z += factor * (flux_contribution_nb_momentum_z.x + flux_contribution_i_momentum_z.x);
+
+        let factor = 0.5 * normal.y;
+        flux_i_density += factor * (momentum_nb.y + momentum_i.y);
+        flux_i_density_energy += factor * (flux_contribution_nb_density_energy.y + flux_contribution_i_density_energy.y);
+        flux_i_momentum.x += factor * (flux_contribution_nb_momentum_x.y + flux_contribution_i_momentum_x.y);
+        flux_i_momentum.y += factor * (flux_contribution_nb_momentum_y.y + flux_contribution_i_momentum_y.y);
+        flux_i_momentum.z += factor * (flux_contribution_nb_momentum_z.y + flux_contribution_i_momentum_z.y);
+
+        let factor = 0.5 * normal.z;
+        flux_i_density += factor * (momentum_nb.z + momentum_i.z);
+        flux_i_density_energy += factor * (flux_contribution_nb_density_energy.z + flux_contribution_i_density_energy.z);
+        flux_i_momentum.x += factor * (flux_contribution_nb_momentum_x.z + flux_contribution_i_momentum_x.z);
+        flux_i_momentum.y += factor * (flux_contribution_nb_momentum_y.z + flux_contribution_i_momentum_y.z);
+        flux_i_momentum.z += factor * (flux_contribution_nb_momentum_z.z + flux_contribution_i_momentum_z.z);
+      } else if nb == -1 { // a wing boundary
+        flux_i_momentum.x += normal.x * pressure_i;
+        flux_i_momentum.y += normal.y * pressure_i;
+        flux_i_momentum.z += normal.z * pressure_i;
+      } else if nb == -2 { // a far field boundary
+        let factor = 0.5 * normal.x;
+        flux_i_density += factor * (ff_variable.momentum.x + momentum_i.x);
+        flux_i_density_energy += factor * (ff_flux_contribution_density_energy.x + flux_contribution_i_density_energy.x);
+        flux_i_momentum.x += factor * (ff_flux_contribution_momentum_x.x + flux_contribution_i_momentum_x.x);
+        flux_i_momentum.y += factor * (ff_flux_contribution_momentum_y.x + flux_contribution_i_momentum_y.x);
+        flux_i_momentum.z += factor * (ff_flux_contribution_momentum_z.x + flux_contribution_i_momentum_z.x);
+
+        let factor = 0.5 * normal.y;
+        flux_i_density += factor * (ff_variable.momentum.y + momentum_i.y);
+        flux_i_density_energy += factor * (ff_flux_contribution_density_energy.y + flux_contribution_i_density_energy.y);
+        flux_i_momentum.x += factor * (ff_flux_contribution_momentum_x.y + flux_contribution_i_momentum_x.y);
+        flux_i_momentum.y += factor * (ff_flux_contribution_momentum_y.y + flux_contribution_i_momentum_y.y);
+        flux_i_momentum.z += factor * (ff_flux_contribution_momentum_z.y + flux_contribution_i_momentum_z.y);
+
+        let factor = 0.5 * normal.z;
+        flux_i_density += factor * (ff_variable.momentum.y + momentum_i.z);
+        flux_i_density_energy += factor * (ff_flux_contribution_density_energy.z + flux_contribution_i_density_energy.z);
+        flux_i_momentum.x += factor * (ff_flux_contribution_momentum_x.z + flux_contribution_i_momentum_x.z);
+        flux_i_momentum.y += factor * (ff_flux_contribution_momentum_y.z + flux_contribution_i_momentum_y.z);
+        flux_i_momentum.z += factor * (ff_flux_contribution_momentum_z.z + flux_contribution_i_momentum_z.z);
+      }
+    }
+
+    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;
+  }
+
+  return fluxes;
+}
+
+const RK : usize = 3;
+
+fn time_step<nelr: usize>(
+  j: usize,
+  old_variables: Variables::<nelr>,
+  step_factors: f32[nelr],
+  fluxes: Variables::<nelr>,
+) -> Variables::<nelr> {
+  let variables : Variables::<nelr>;
+
+  for i = 0 to nelr {
+    let factor = step_factors[i] / (RK + 1 - j) as f32;
+
+    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];
+  }
+  
+  return variables;
+}
+
+#[entry]
+fn euler<nelr: usize>(
+  iterations: usize,
+  variables: Variables::<nelr>,
+  areas: f32[nelr],
+  elements_surrounding_elements: i32[NNB, nelr],
+  normals: Normals::<nelr>,
+  ff_variable: Variable,
+  ff_flux_contribution_density_energy: float3,
+  ff_flux_contribution_momentum_x: float3,
+  ff_flux_contribution_momentum_y: float3,
+  ff_flux_contribution_momentum_z: float3,
+  step_factors: f32[nelr],
+ ) -> Variables::<nelr> {
+  for i = 0 to iterations {
+    let old_variables = variables;
+    let step_factors = compute_step_factor::<nelr>(variables, areas);
+
+    for j = 0 to RK {
+      let fluxes = compute_flux::<nelr>(variables, elements_surrounding_elements,
+                                        normals, ff_variable, ff_flux_contribution_density_energy,
+                                        ff_flux_contribution_momentum_x,
+                                        ff_flux_contribution_momentum_y,
+                                        ff_flux_contribution_momentum_z);
+      variables = time_step::<nelr>(j, old_variables, step_factors, fluxes);
+    }
+  }
+
+  return variables;
+}

juno_samples/rodinia/cfd/src/main.rs

+#![feature(concat_idents)]
+
+juno_build::juno!("euler");
+
+fn main() {
+    todo!()
+}