add bench for euler

Cargo.toml

@@ -30,3 +30,6 @@ criterion = "0.3.0"
 name = "system"
 path = "src/benches/bench.rs"
 harness = false
+
+[profile.bench]
+debug = true

main.js

@@ -211,20 +211,7 @@ import { EulerUniverse, Universe, default as init, set_panic_hook as setPanicHoo
     function drawMe(timeOfDraw) {
         gl.clear(gl.COLOR_BUFFER_BIT | gl.DEPTH_BUFFER_BIT);
 
-        let dt = (timeOfDraw - t) / TIMEFACTOR;
+        let dt = 0.01;
-        t = timeOfDraw;
-        if (firstDraw || dt <= 0.0) {
-            firstDraw = false;
-            dt = MAX_DT;
-        } else {
-            if (dt >= MAX_DT) {
-                warnTime += 1;
-                if (warnTime !== -2) {
-                    console.warn("Can not keep up with framerate");
-                }
-                dt = MAX_DT;
-            }
-        }
 
         let fieldPtr;
         if (chosenField.value === 0) {

src/benches/bench.rs

@@ -1,6 +1,6 @@
 use criterion::{black_box, criterion_group, criterion_main, Criterion};
 use maxwell::operators::{SbpOperator, Upwind4, UpwindOperator, SBP4};
-use maxwell::System;
+use maxwell::{EulerSystem, System};
 
 fn advance_system<SBP: SbpOperator>(universe: &mut System<SBP>, n: usize) {
     for _ in 0..n {
@@ -51,4 +51,56 @@ fn performance_benchmark(c: &mut Criterion) {
 }
 
 criterion_group!(benches, performance_benchmark);
-criterion_main!(benches);
+
+fn advance_euler_system<SBP: SbpOperator>(universe: &mut EulerSystem<SBP>, n: usize) {
+    for _ in 0..n {
+        universe.advance(0.01);
+    }
+}
+
+fn advance_euler_system_upwind<UO: UpwindOperator>(universe: &mut EulerSystem<UO>, n: usize) {
+    for _ in 0..n {
+        universe.advance_upwind(0.01);
+    }
+}
+
+fn euler_performance_benchmark(c: &mut Criterion) {
+    let mut group = c.benchmark_group("EulerSystem");
+    group.sample_size(25);
+
+    let w = 40;
+    let h = 26;
+    let x = ndarray::Array1::linspace(-10.0, 10.0, w);
+    let x = x.broadcast((h, w)).unwrap();
+    let y = ndarray::Array1::linspace(-10.0, 10.0, h);
+    let y = y.broadcast((w, h)).unwrap().reversed_axes();
+
+    let mut universe = EulerSystem::<Upwind4>::new(x.into_owned(), y.into_owned());
+    group.bench_function("advance_euler", |b| {
+        b.iter(|| {
+            universe.vortex(0.0, 0.0);
+            advance_euler_system(&mut universe, black_box(10))
+        })
+    });
+
+    let mut universe = EulerSystem::<Upwind4>::new(x.into_owned(), y.into_owned());
+    group.bench_function("advance_euler_upwind", |b| {
+        b.iter(|| {
+            universe.vortex(0.0, 0.0);
+            advance_euler_system_upwind(&mut universe, black_box(10))
+        })
+    });
+
+    let mut universe = EulerSystem::<SBP4>::new(x.into_owned(), y.into_owned());
+    group.bench_function("advance_euler_trad4", |b| {
+        b.iter(|| {
+            universe.vortex(0.0, 0.0);
+            advance_euler_system(&mut universe, black_box(10))
+        })
+    });
+
+    group.finish();
+}
+
+criterion_group!(euler_benches, euler_performance_benchmark);
+criterion_main!(benches, euler_benches);

src/euler.rs

@@ -443,7 +443,7 @@ fn SAT_characteristics<SBP: SbpOperator>(
         let hi = (k.ny() - 1) as f32 * SBP::h()[0];
         let sign = -1.0;
         let tau = 1.0;
-        let slice = s![y.nx() - 1, ..];
+        let slice = s![y.ny() - 1, ..];
         SAT_characteristic(
             k.north_mut(),
             y.north(),

src/lib.rs

@@ -146,51 +146,7 @@ impl EulerUniverse {
     }
 
     pub fn init(&mut self, x0: f32, y0: f32) {
-        // Should parametrise such that we have radius, drop in pressure at center, etc
+        self.0.vortex(x0, y0)
-        let rstar = 1.0;
-        let eps = 3.0;
-        #[allow(non_snake_case)]
-        let M = 0.5;
-
-        let p_inf = 1.0 / (euler::GAMMA * M * M);
-        let t = 0.0;
-
-        let nx = self.0.grid.nx();
-        let ny = self.0.grid.ny();
-
-        for j in 0..ny {
-            for i in 0..nx {
-                let x = self.0.grid.x[(j, i)];
-                let y = self.0.grid.y[(j, i)];
-
-                let dx = (x - x0) - t;
-                let dy = y - y0;
-                let f = (1.0 - (dx * dx + dy * dy)) / (rstar * rstar);
-
-                use euler::GAMMA;
-                use std::f32::consts::PI;
-                let u =
-                    1.0 - eps * dy / (2.0 * PI * p_inf.sqrt() * rstar * rstar) * (f / 2.0).exp();
-                let v =
-                    0.0 + eps * dx / (2.0 * PI * p_inf.sqrt() * rstar * rstar) * (f / 2.0).exp();
-                let rho = f32::powf(
-                    1.0 - eps * eps * (GAMMA - 1.0) * M * M
-                        / (8.0 * PI * PI * p_inf * rstar * rstar)
-                        * f.exp(),
-                    1.0 / (GAMMA - 1.0),
-                );
-                assert!(rho > 0.0);
-                let p = p_inf * rho.powf(GAMMA);
-                assert!(p > 0.0);
-                let e = p / (GAMMA - 1.0) + rho * (u * u + v * v) / 2.0;
-                assert!(e > 0.0);
-
-                self.0.sys.0[(0, j, i)] = rho;
-                self.0.sys.0[(1, j, i)] = rho * u;
-                self.0.sys.0[(2, j, i)] = rho * v;
-                self.0.sys.0[(3, j, i)] = e;
-            }
-        }
     }
 
     pub fn advance(&mut self, dt: f32) {
@@ -239,6 +195,54 @@ impl<SBP: operators::SbpOperator> EulerSystem<SBP> {
         );
         std::mem::swap(&mut self.sys.0, &mut self.sys.1);
     }
+
+    pub fn vortex(&mut self, x0: f32, y0: f32) {
+        // Should parametrise such that we have radius, drop in pressure at center, etc
+        let rstar = 1.0;
+        let eps = 3.0;
+        #[allow(non_snake_case)]
+        let M = 0.5;
+
+        let p_inf = 1.0 / (euler::GAMMA * M * M);
+        let t = 0.0;
+
+        let nx = self.grid.nx();
+        let ny = self.grid.ny();
+
+        for j in 0..ny {
+            for i in 0..nx {
+                let x = self.grid.x[(j, i)];
+                let y = self.grid.y[(j, i)];
+
+                let dx = (x - x0) - t;
+                let dy = y - y0;
+                let f = (1.0 - (dx * dx + dy * dy)) / (rstar * rstar);
+
+                use euler::GAMMA;
+                use std::f32::consts::PI;
+                let u =
+                    1.0 - eps * dy / (2.0 * PI * p_inf.sqrt() * rstar * rstar) * (f / 2.0).exp();
+                let v =
+                    0.0 + eps * dx / (2.0 * PI * p_inf.sqrt() * rstar * rstar) * (f / 2.0).exp();
+                let rho = f32::powf(
+                    1.0 - eps * eps * (GAMMA - 1.0) * M * M
+                        / (8.0 * PI * PI * p_inf * rstar * rstar)
+                        * f.exp(),
+                    1.0 / (GAMMA - 1.0),
+                );
+                assert!(rho > 0.0);
+                let p = p_inf * rho.powf(GAMMA);
+                assert!(p > 0.0);
+                let e = p / (GAMMA - 1.0) + rho * (u * u + v * v) / 2.0;
+                assert!(e > 0.0);
+
+                self.sys.0[(0, j, i)] = rho;
+                self.sys.0[(1, j, i)] = rho * u;
+                self.sys.0[(2, j, i)] = rho * v;
+                self.sys.0[(3, j, i)] = e;
+            }
+        }
+    }
 }
 
 impl<SBP: operators::UpwindOperator> EulerSystem<SBP> {