diff --git a/benches/euler.rs b/benches/euler.rs index 4c8c280..07fa2a4 100644 --- a/benches/euler.rs +++ b/benches/euler.rs @@ -1,14 +1,14 @@ use criterion::{black_box, criterion_group, criterion_main, Criterion}; +use sbp::euler::System; use sbp::operators::{SbpOperator, Upwind4, UpwindOperator, SBP4}; -use sbp::EulerSystem; -fn advance_system(universe: &mut EulerSystem, n: usize) { +fn advance_system(universe: &mut System, n: usize) { for _ in 0..n { universe.advance(1.0 / 40.0 * 0.2); } } -fn advance_system_upwind(universe: &mut EulerSystem, n: usize) { +fn advance_system_upwind(universe: &mut System, n: usize) { for _ in 0..n { universe.advance_upwind(1.0 / 40.0 * 0.2); } @@ -25,26 +25,26 @@ fn performance_benchmark(c: &mut Criterion) { 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()); + let mut universe = System::::new(x.into_owned(), y.into_owned()); group.bench_function("advance", |b| { b.iter(|| { - universe.vortex(0.0, 0.0); + universe.init_with_vortex(0.0, 0.0); advance_system(&mut universe, black_box(20)) }) }); - let mut universe = EulerSystem::::new(x.into_owned(), y.into_owned()); + let mut universe = System::::new(x.into_owned(), y.into_owned()); group.bench_function("advance_upwind", |b| { b.iter(|| { - universe.vortex(0.0, 0.0); + universe.init_with_vortex(0.0, 0.0); advance_system_upwind(&mut universe, black_box(20)) }) }); - let mut universe = EulerSystem::::new(x.into_owned(), y.into_owned()); + let mut universe = System::::new(x.into_owned(), y.into_owned()); group.bench_function("advance_trad4", |b| { b.iter(|| { - universe.vortex(0.0, 0.0); + universe.init_with_vortex(0.0, 0.0); advance_system(&mut universe, black_box(20)) }) }); diff --git a/src/euler.rs b/src/euler.rs index 0787a9b..da73219 100644 --- a/src/euler.rs +++ b/src/euler.rs @@ -5,6 +5,106 @@ use ndarray::{azip, Zip}; pub const GAMMA: f32 = 1.4; +// A collection of buffers that allows one to efficiently +// move to the next state +pub struct System { + sys: (Field, Field), + wb: WorkBuffers, + grid: Grid, +} + +impl System { + pub fn new(x: ndarray::Array2, y: ndarray::Array2) -> Self { + let grid = Grid::new(x, y).expect( + "Could not create grid. Different number of elements compared to width*height?", + ); + let nx = grid.nx(); + let ny = grid.ny(); + Self { + sys: (Field::new(ny, nx), Field::new(ny, nx)), + grid, + wb: WorkBuffers::new(ny, nx), + } + } + + pub fn advance(&mut self, dt: f32) { + advance( + RHS_trad, + &self.sys.0, + &mut self.sys.1, + dt, + &self.grid, + Some(&mut self.wb), + ); + std::mem::swap(&mut self.sys.0, &mut self.sys.1); + } + + pub fn init_with_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 / (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 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; + } + } + } + + pub fn field(&self) -> &Field { + &self.sys.0 + } +} + +impl System { + pub fn advance_upwind(&mut self, dt: f32) { + advance( + RHS_upwind, + &self.sys.0, + &mut self.sys.1, + dt, + &self.grid, + Some(&mut self.wb), + ); + std::mem::swap(&mut self.sys.0, &mut self.sys.1); + } +} + #[derive(Clone, Debug)] /// A 4 x ny x nx array pub struct Field(pub(crate) Array3); diff --git a/src/grid.rs b/src/grid.rs index c8a0008..e3f049c 100644 --- a/src/grid.rs +++ b/src/grid.rs @@ -1,4 +1,6 @@ use ndarray::Array2; + +#[derive(Debug, Clone)] pub struct Grid where SBP: super::operators::SbpOperator, diff --git a/src/lib.rs b/src/lib.rs index b344a0c..0a0b710 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,6 +1,6 @@ use wasm_bindgen::prelude::*; -mod euler; +pub mod euler; mod grid; mod maxwell; pub mod operators; @@ -121,18 +121,12 @@ fn gaussian(x: f32, x0: f32, y: f32, y0: f32) -> f32 { 1.0 / (2.0 * f32::consts::PI * sigma * sigma) * (-(x * x + y * y) / (2.0 * sigma * sigma)).exp() } -pub struct EulerSystem { - sys: (euler::Field, euler::Field), - wb: euler::WorkBuffers, - grid: Grid, -} - #[wasm_bindgen] -pub struct EulerUniverse(EulerSystem); +pub struct EulerUniverse(euler::System); impl EulerUniverse { pub fn new(x: ndarray::Array2, y: ndarray::Array2) -> Self { - Self(EulerSystem::new(x, y)) + Self(euler::System::new(x, y)) } } @@ -142,11 +136,11 @@ impl EulerUniverse { pub fn new_with_slice(height: usize, width: usize, x: &[f32], y: &[f32]) -> Self { let x = ndarray::Array2::from_shape_vec((height, width), x.to_vec()).unwrap(); let y = ndarray::Array2::from_shape_vec((height, width), y.to_vec()).unwrap(); - Self(EulerSystem::new(x, y)) + Self(euler::System::new(x, y)) } pub fn init(&mut self, x0: f32, y0: f32) { - self.0.vortex(x0, y0) + self.0.init_with_vortex(x0, y0) } pub fn advance(&mut self, dt: f32) { @@ -158,105 +152,16 @@ impl EulerUniverse { } pub fn get_rho_ptr(&self) -> *const u8 { - self.0.sys.0.rho().as_ptr() as *const u8 + self.0.field().rho().as_ptr() as *const u8 } pub fn get_rhou_ptr(&self) -> *const u8 { - self.0.sys.0.rhou().as_ptr() as *const u8 + self.0.field().rhou().as_ptr() as *const u8 } pub fn get_rhov_ptr(&self) -> *const u8 { - self.0.sys.0.rhov().as_ptr() as *const u8 + self.0.field().rhov().as_ptr() as *const u8 } pub fn get_e_ptr(&self) -> *const u8 { - self.0.sys.0.e().as_ptr() as *const u8 - } -} - -impl EulerSystem { - pub fn new(x: ndarray::Array2, y: ndarray::Array2) -> Self { - let grid = Grid::new(x, y).expect( - "Could not create grid. Different number of elements compared to width*height?", - ); - let nx = grid.nx(); - let ny = grid.ny(); - Self { - sys: (euler::Field::new(ny, nx), euler::Field::new(ny, nx)), - grid, - wb: euler::WorkBuffers::new(ny, nx), - } - } - - pub fn advance(&mut self, dt: f32) { - euler::advance( - euler::RHS_trad, - &self.sys.0, - &mut self.sys.1, - dt, - &self.grid, - Some(&mut self.wb), - ); - 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 { - pub fn advance_upwind(&mut self, dt: f32) { - euler::advance( - euler::RHS_upwind, - &self.sys.0, - &mut self.sys.1, - dt, - &self.grid, - Some(&mut self.wb), - ); - std::mem::swap(&mut self.sys.0, &mut self.sys.1); + self.0.field().e().as_ptr() as *const u8 } }