diff --git a/Cargo.toml b/Cargo.toml index ba48a31..4694ead 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -30,3 +30,6 @@ criterion = "0.3.0" name = "system" path = "src/benches/bench.rs" harness = false + +[profile.bench] +debug = true diff --git a/main.js b/main.js index 35071c1..a37bf4a 100644 --- a/main.js +++ b/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; - 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 dt = 0.01; let fieldPtr; if (chosenField.value === 0) { diff --git a/src/benches/bench.rs b/src/benches/bench.rs index 0d03a19..4a0ec82 100644 --- a/src/benches/bench.rs +++ b/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(universe: &mut System, 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(universe: &mut EulerSystem, n: usize) { + for _ in 0..n { + universe.advance(0.01); + } +} + +fn advance_euler_system_upwind(universe: &mut EulerSystem, 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::::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::::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::::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); diff --git a/src/euler.rs b/src/euler.rs index f7f3a61..b1ee6a0 100644 --- a/src/euler.rs +++ b/src/euler.rs @@ -443,7 +443,7 @@ fn SAT_characteristics( 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(), diff --git a/src/lib.rs b/src/lib.rs index ac8c56f..49908ca 100644 --- a/src/lib.rs +++ b/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 - 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; - } - } + self.0.vortex(x0, y0) } pub fn advance(&mut self, dt: f32) { @@ -239,6 +195,54 @@ impl EulerSystem { ); 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 EulerSystem {