From 4673fbfbf7756c47f39daf4851ab5e6476d3db87 Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Sat, 25 Jan 2020 20:40:55 +0100 Subject: [PATCH] checkpoint --- main.js | 30 +-- src/euler.rs | 589 +++++++++++++++++++++++++++++++++++++++++++++++++++ src/grid.rs | 24 ++- src/lib.rs | 132 +++++++++++- 4 files changed, 758 insertions(+), 17 deletions(-) create mode 100644 src/euler.rs diff --git a/main.js b/main.js index 5b079f6..f47830f 100644 --- a/main.js +++ b/main.js @@ -1,4 +1,4 @@ -import { Universe, default as init, set_panic_hook as setPanicHook } from "./maxwell.js"; +import { EulerUniverse, Universe, default as init, set_panic_hook as setPanicHook } from "./maxwell.js"; /** * Initialises and runs the Maxwell solver, @@ -49,7 +49,7 @@ import { Universe, default as init, set_panic_hook as setPanicHook } from "./max if (uChosenField == 0) { r = vField + 0.5; } else if (uChosenField == 1) { - g = (vField + 1.0)/2.0; + g = 2.5*(vField - 1.0) + 0.5; } else { b = vField + 0.5; } @@ -103,8 +103,8 @@ import { Universe, default as init, set_panic_hook as setPanicHook } from "./max for (let j = 0; j < height; j += 1) { for (let i = 0; i < width; i += 1) { const n = width*j + i; - x[n] = i / (width - 1.0); - y[n] = j / (height - 1.0); + x[n] = 10.0*(i / (width - 1.0) - 0.5); + y[n] = 20.0*(j / (height - 1.0)); if (DIAMOND) { @@ -116,7 +116,7 @@ import { Universe, default as init, set_panic_hook as setPanicHook } from "./max } } - const universe = new Universe(width, height, x, y); + const universe = new EulerUniverse(height, width, x, y); // Transfer x, y to cpu, prepare fBuffer @@ -201,7 +201,7 @@ import { Universe, default as init, set_panic_hook as setPanicHook } from "./max }; chosenField.cycle(); - universe.init(0.5, 0.5); + universe.init(0, 10); /** * Integrates and draws the next iteration @@ -227,14 +227,18 @@ import { Universe, default as init, set_panic_hook as setPanicHook } from "./max let fieldPtr; if (chosenField.value === 0) { - fieldPtr = universe.get_ex_ptr(); + fieldPtr = universe.get_rho_ptr(); } else if (chosenField.value === 1) { - fieldPtr = universe.get_hz_ptr(); + fieldPtr = universe.get_rhou_ptr(); + } else if (chosenField.value == 2) { + fieldPtr = universe.get_rhov_ptr(); } else { - fieldPtr = universe.get_ey_ptr(); - } + fieldPtr = universe.get_e_ptr(); + }; const field = new Float32Array(wasm.memory.buffer, fieldPtr, width*height); gl.bufferData(gl.ARRAY_BUFFER, field, gl.DYNAMIC_DRAW); + console.log(field.reduce((min, v) => v < min ? v : min)); + console.log(field.reduce((max, v) => v > max ? v : max)); { const offset = 0; @@ -243,8 +247,8 @@ import { Universe, default as init, set_panic_hook as setPanicHook } from "./max gl.drawElements(gl.TRIANGLES, vertexCount, type, offset); } - universe.advance_upwind(dt/2); - universe.advance_upwind(dt/2); + universe.advance(dt/2); + universe.advance(dt/2); window.requestAnimationFrame(drawMe); } @@ -267,7 +271,7 @@ import { Universe, default as init, set_panic_hook as setPanicHook } from "./max // Must adjust for bbox and transformations for x/y const mousex = event.clientX / window.innerWidth; const mousey = event.clientY / window.innerHeight; - universe.init(mousex, 1.0 - mousey); + universe.init(10*(mousex-0.5), 20.0*(1.0 - mousey)); }, {"passive": true}); resizeCanvas(); diff --git a/src/euler.rs b/src/euler.rs new file mode 100644 index 0000000..f5257e5 --- /dev/null +++ b/src/euler.rs @@ -0,0 +1,589 @@ +use super::operators::{SbpOperator, UpwindOperator}; +use super::Grid; +use ndarray::prelude::*; +use ndarray::{azip, Zip}; + +pub const GAMMA: f32 = 1.4; + +#[derive(Clone, Debug)] +/// A 4 x ny x nx array +pub struct Field(pub(crate) Array3); + +impl std::ops::Deref for Field { + type Target = Array3; + fn deref(&self) -> &Self::Target { + &self.0 + } +} + +impl std::ops::DerefMut for Field { + fn deref_mut(&mut self) -> &mut Self::Target { + &mut self.0 + } +} + +impl Field { + pub fn new(ny: usize, nx: usize) -> Self { + let field = Array3::zeros((4, ny, nx)); + + Self(field) + } + + pub fn nx(&self) -> usize { + self.0.shape()[2] + } + pub fn ny(&self) -> usize { + self.0.shape()[1] + } + + pub fn rho(&self) -> ArrayView2 { + self.slice(s![0, .., ..]) + } + pub fn rhou(&self) -> ArrayView2 { + self.slice(s![1, .., ..]) + } + pub fn rhov(&self) -> ArrayView2 { + self.slice(s![2, .., ..]) + } + pub fn e(&self) -> ArrayView2 { + self.slice(s![3, .., ..]) + } + + pub fn rho_mut(&mut self) -> ArrayViewMut2 { + self.slice_mut(s![0, .., ..]) + } + pub fn rhou_mut(&mut self) -> ArrayViewMut2 { + self.slice_mut(s![1, .., ..]) + } + pub fn rhov_mut(&mut self) -> ArrayViewMut2 { + self.slice_mut(s![2, .., ..]) + } + pub fn e_mut(&mut self) -> ArrayViewMut2 { + self.slice_mut(s![3, .., ..]) + } + + #[allow(unused)] + pub fn components( + &self, + ) -> ( + ArrayView2, + ArrayView2, + ArrayView2, + ArrayView2, + ) { + (self.rho(), self.rhou(), self.rhov(), self.e()) + } + #[allow(unused)] + pub fn components_mut( + &mut self, + ) -> ( + ArrayViewMut2, + ArrayViewMut2, + ArrayViewMut2, + ArrayViewMut2, + ) { + let mut iter = self.0.outer_iter_mut(); + + let rho = iter.next().unwrap(); + let rhou = iter.next().unwrap(); + let rhov = iter.next().unwrap(); + let e = iter.next().unwrap(); + assert_eq!(iter.next(), None); + + (rho, rhou, rhov, e) + } + + fn north(&self) -> ArrayView2 { + self.slice(s![.., self.ny() - 1, ..]) + } + fn south(&self) -> ArrayView2 { + self.slice(s![.., 0, ..]) + } + fn east(&self) -> ArrayView2 { + self.slice(s![.., .., self.nx() - 1]) + } + fn west(&self) -> ArrayView2 { + self.slice(s![.., .., 0]) + } + fn north_mut(&mut self) -> ArrayViewMut2 { + let ny = self.ny(); + self.slice_mut(s![.., ny - 1, ..]) + } + fn south_mut(&mut self) -> ArrayViewMut2 { + self.slice_mut(s![.., 0, ..]) + } + fn east_mut(&mut self) -> ArrayViewMut2 { + let nx = self.nx(); + self.slice_mut(s![.., .., nx - 1]) + } + fn west_mut(&mut self) -> ArrayViewMut2 { + self.slice_mut(s![.., .., 0]) + } +} + +#[allow(unused)] +pub(crate) fn advance_upwind( + prev: &Field, + fut: &mut Field, + dt: f32, + grid: &Grid, + work_buffers: Option<&mut WorkBuffers>, +) where + UO: UpwindOperator, +{ + assert_eq!(prev.0.shape(), fut.0.shape()); + + let mut wb: WorkBuffers; + let (y, k, tmp) = if let Some(x) = work_buffers { + (&mut x.y, &mut x.buf, &mut x.tmp) + } else { + wb = WorkBuffers::new(prev.nx(), prev.ny()); + (&mut wb.y, &mut wb.buf, &mut wb.tmp) + }; + + let boundaries = BoundaryTerms { + north: Boundary::This, + south: Boundary::This, + west: Boundary::This, + east: Boundary::This, + }; + + for i in 0..4 { + // y = y0 + c*kn + y.assign(&prev); + match i { + 0 => {} + 1 | 2 => { + y.scaled_add(1.0 / 2.0 * dt, &k[i - 1]); + } + 3 => { + y.scaled_add(dt, &k[i - 1]); + } + _ => { + unreachable!(); + } + }; + + RHS_upwind(&mut k[i], &y, grid, &boundaries, tmp); + } + + Zip::from(&mut fut.0) + .and(&prev.0) + .and(&*k[0]) + .and(&*k[1]) + .and(&*k[2]) + .and(&*k[3]) + .apply(|y1, &y0, &k1, &k2, &k3, &k4| *y1 = y0 + dt / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4)); +} + +/// Solving (Au)_x + (Bu)_y +/// with: +/// A B +/// [ 0, 0, 0] [ 0, 1, 0] +/// [ 0, 0, -1] [ 1, 0, 0] +/// [ 0, -1, 0] [ 0, 0, 0] +pub(crate) fn advance( + prev: &Field, + fut: &mut Field, + dt: f32, + grid: &Grid, + work_buffers: Option<&mut WorkBuffers>, +) where + SBP: SbpOperator, +{ + assert_eq!(prev.0.shape(), fut.0.shape()); + + let mut wb: WorkBuffers; + let (y, k, tmp) = if let Some(x) = work_buffers { + (&mut x.y, &mut x.buf, &mut x.tmp) + } else { + wb = WorkBuffers::new(prev.nx(), prev.ny()); + (&mut wb.y, &mut wb.buf, &mut wb.tmp) + }; + + let boundaries = BoundaryTerms { + north: Boundary::This, + south: Boundary::This, + west: Boundary::This, + east: Boundary::This, + }; + + for i in 0..4 { + // y = y0 + c*kn + y.assign(&prev); + match i { + 0 => {} + 1 | 2 => { + y.scaled_add(1.0 / 2.0 * dt, &k[i - 1]); + } + 3 => { + y.scaled_add(dt, &k[i - 1]); + } + _ => { + unreachable!(); + } + }; + + RHS(&mut k[i], &y, grid, &boundaries, tmp); + } + + Zip::from(&mut fut.0) + .and(&prev.0) + .and(&*k[0]) + .and(&*k[1]) + .and(&*k[2]) + .and(&*k[3]) + .apply(|y1, &y0, &k1, &k2, &k3, &k4| *y1 = y0 + dt / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4)); +} + +fn pressure(gamma: f32, rho: f32, rhou: f32, rhov: f32, e: f32) -> f32 { + (gamma - 1.0) * (e - (rhou * rhou + rhov * rhov) / (2.0 * rho)) +} + +#[allow(non_snake_case)] +/// This flux is rotated by the grid metrics +/// (Au)_x + (Bu)_y = 1/J [ +/// (J xi_x Au)_xi + (J eta_x Au)_eta +/// (J xi_y Bu)_xi + (J eta_y Bu)_eta +/// ] +/// where J is the grid determinant +/// +/// This is used both in fluxes and SAT terms +fn RHS( + k: &mut Field, + y: &Field, + grid: &Grid, + boundaries: &BoundaryTerms, + tmp: &mut (Field, Field, Field, Field), +) { + let ehat = &mut tmp.0; + let fhat = &mut tmp.1; + fluxes([ehat, fhat], y, grid); + let de = &mut tmp.2; + let df = &mut tmp.3; + + SBP::diffxi(ehat.rho(), de.rho_mut()); + SBP::diffxi(ehat.rhou(), de.rhou_mut()); + SBP::diffxi(ehat.rhov(), de.rhov_mut()); + SBP::diffxi(ehat.e(), de.e_mut()); + + SBP::diffeta(fhat.rho(), df.rho_mut()); + SBP::diffeta(fhat.rhou(), df.rhou_mut()); + SBP::diffeta(fhat.rhov(), df.rhov_mut()); + SBP::diffeta(fhat.e(), df.e_mut()); + + // And dissipation... + + ndarray::azip!((out in &mut k.0, + eflux in &de.0, + fflux in &df.0, + detj in &grid.detj.broadcast((4, y.ny(), y.nx())).unwrap()) { + *out = (-eflux - fflux)/detj + }); + + SAT_characteristics(k, y, grid, boundaries); +} + +#[allow(non_snake_case)] +#[allow(unused)] +fn RHS_upwind( + k: &mut Field, + y: &Field, + grid: &Grid, + boundaries: &BoundaryTerms, + tmp: &mut (Field, Field, Field, Field), +) { + // fluxes(k, y, grid, tmp); + // dissipation(k, y, grid, tmp); + + SAT_characteristics(k, y, grid, boundaries); + + azip!((k in &mut k.0, + &detj in &grid.detj.broadcast((3, y.ny(), y.nx())).unwrap()) { + *k /= detj; + }); +} + +fn fluxes(k: [&mut Field; 2], y: &Field, grid: &Grid) { + let j_dxi_dx = grid.detj_dxi_dx.view(); + let j_dxi_dy = grid.detj_dxi_dy.view(); + let j_deta_dx = grid.detj_deta_dx.view(); + let j_deta_dy = grid.detj_deta_dy.view(); + + let rho = y.rho(); + let rhou = y.rhou(); + let rhov = y.rhov(); + let e = y.e(); + + for j in 0..y.ny() { + for i in 0..y.nx() { + let rho = rho[(j, i)]; + let rhou = rhou[(j, i)]; + let rhov = rhov[(j, i)]; + let e = e[(j, i)]; + let p = pressure(GAMMA, rho, rhou, rhov, e); + + let ef = [ + rhou, + rhou * rhou / rho + p, + rhou * rhov / rho, + rhou * (e + p) / rho, + ]; + let ff = [ + rhov, + rhou * rhov / rho, + rhov * rhov / rho + p, + rhov * (e + p) / rho, + ]; + + for comp in 0..4 { + let eflux = j_dxi_dx[(j, i)] * ef[comp] + j_dxi_dy[(j, i)] * ff[comp]; + let fflux = j_deta_dx[(j, i)] * ef[comp] + j_deta_dy[(j, i)] * ff[comp]; + + k[0][(comp, j, i)] = eflux; + k[1][(comp, j, i)] = fflux; + } + } + } +} + +#[allow(unused)] +fn dissipation( + k: &mut Field, + y: &Field, + grid: &Grid, + tmp: &mut (Array2, Array2, Array2, Array2), +) { + todo!() +} + +#[derive(Clone, Debug)] +pub enum Boundary { + This, +} + +#[derive(Clone, Debug)] +pub struct BoundaryTerms { + pub north: Boundary, + pub south: Boundary, + pub east: Boundary, + pub west: Boundary, +} + +#[allow(non_snake_case)] +/// Boundary conditions (SAT) +fn SAT_characteristics( + k: &mut Field, + y: &Field, + grid: &Grid, + _boundaries: &BoundaryTerms, +) { + // North boundary + { + let hi = (k.ny() - 1) as f32 * SBP::h()[0]; + let sign = -1.0; + let tau = 1.0; + let slice = s![y.nx() - 1, ..]; + SAT_characteristic( + k.north_mut(), + y.north(), + y.south(), // Self South + hi, + sign, + tau, + grid.detj.slice(slice), + grid.detj_deta_dx.slice(slice), + grid.detj_deta_dy.slice(slice), + ); + } + // South boundary + { + let hi = (k.ny() - 1) as f32 * SBP::h()[0]; + let sign = 1.0; + let tau = -1.0; + let slice = s![0, ..]; + SAT_characteristic( + k.south_mut(), + y.south(), + y.north(), // Self North + hi, + sign, + tau, + grid.detj.slice(slice), + grid.detj_deta_dx.slice(slice), + grid.detj_deta_dy.slice(slice), + ); + } + // West Boundary + { + let hi = (k.nx() - 1) as f32 * SBP::h()[0]; + let sign = 1.0; + let tau = -1.0; + let slice = s![.., 0]; + println!("{:?}", slice); + SAT_characteristic( + k.west_mut(), + y.west(), + y.east(), // Self East + hi, + sign, + tau, + grid.detj.slice(slice), + grid.detj_dxi_dx.slice(slice), + grid.detj_dxi_dy.slice(slice), + ); + } + // East Boundary + { + let hi = (k.nx() - 1) as f32 * SBP::h()[0]; + let sign = -1.0; + let tau = 1.0; + let slice = s![.., y.nx() - 1]; + SAT_characteristic( + k.east_mut(), + y.east(), + y.west(), // Self West + hi, + sign, + tau, + grid.detj.slice(slice), + grid.detj_dxi_dx.slice(slice), + grid.detj_dxi_dy.slice(slice), + ); + } +} + +#[allow(non_snake_case)] +/// Boundary conditions (SAT) +fn SAT_characteristic( + mut k: ArrayViewMut2, + y: ArrayView2, + z: ArrayView2, // Size 4 x n (all components in line) + hi: f32, + sign: f32, + tau: f32, + detj: ArrayView1, + detj_d_dx: ArrayView1, + detj_d_dy: ArrayView1, +) { + assert_eq!(detj.shape(), detj_d_dx.shape()); + assert_eq!(detj.shape(), detj_d_dy.shape()); + assert_eq!(y.shape(), z.shape()); + assert_eq!(y.shape()[0], 4); + assert_eq!(y.shape()[1], detj.shape()[0]); + + for i in 0..z.shape()[1] { + let rho = y[(0, i)]; + let rhou = y[(1, i)]; + let rhov = y[(2, i)]; + let e = y[(3, i)]; + + let kx_ = detj_d_dx[i] / detj[i]; + let ky_ = detj_d_dy[i] / detj[i]; + + let (kx, ky) = { + let r = f32::hypot(kx_, ky_); + (kx_ / r, ky_ / r) + }; + + let u = rhou / rho; + let v = rhov / rho; + + let theta = kx * u + ky * v; + + let p = pressure(GAMMA, rho, rhou, rhov, e); + let c = (GAMMA * p / rho).sqrt(); + let phi2 = (GAMMA - 1.0) * (u * u + v * v) / 2.0; + + let phi2_c2 = (phi2 + c * c) / (GAMMA - 1.0); + + let T = [ + [1.0, 0.0, 1.0, 1.0], + [u, ky, u + kx * c, u - kx * c], + [v, -kx, v + ky * c, v - ky * c], + [ + phi2 / (GAMMA - 1.0), + ky * u - kx * v, + phi2_c2 + c * theta, + phi2_c2 - c * theta, + ], + ]; + let U = kx_ * u + ky_ * v; + let L = [ + U, + U, + U + c * f32::hypot(kx_, ky_), + U - c * f32::hypot(kx_, ky_), + ]; + let beta = 1.0 / (2.0 * c * c); + let TI = [ + [ + 1.0 - phi2 / (c * c), + (GAMMA - 1.0) * u / (c * c), + (GAMMA - 1.0) * v / (c * c), + -(GAMMA - 1.0) / (c * c), + ], + [-(ky * u - kx * v), ky, -kx, 0.0], + [ + beta * (phi2 - c * theta), + beta * (kx * c - (GAMMA - 1.0) * u), + beta * (ky * c - (GAMMA - 1.0) * v), + beta * (GAMMA - 1.0), + ], + [ + beta * (phi2 + c * theta), + -beta * (kx * c + (GAMMA - 1.0) * u), + -beta * (ky * c + (GAMMA - 1.0) * v), + beta * (GAMMA - 1.0), + ], + ]; + + let res = [ + rho - z[(0, i)], + rhou - z[(1, i)], + rhov - z[(2, i)], + e - z[(3, i)], + ]; + let mut TIres = [0.0; 4]; + for row in 0..4 { + for col in 0..4 { + TIres[row] += TI[row][col] * res[col]; + } + } + + // L + sign(abs(L)) * TIres + let mut LTIres = [0.0; 4]; + for row in 0..4 { + LTIres[row] = (L[row] + sign * L[row].abs()) * TIres[row]; + } + + // T*LTIres + let mut TLTIres = [0.0; 4]; + for row in 0..4 { + for col in 0..4 { + TLTIres[row] += T[row][col] * LTIres[col]; + } + } + + for comp in 0..4 { + k[(comp, i)] += hi * tau * TLTIres[comp]; + } + } +} + +pub struct WorkBuffers { + y: Field, + buf: [Field; 4], + tmp: (Field, Field, Field, Field), +} + +impl WorkBuffers { + pub fn new(nx: usize, ny: usize) -> Self { + let arr3 = Field::new(nx, ny); + Self { + y: arr3.clone(), + buf: [arr3.clone(), arr3.clone(), arr3.clone(), arr3.clone()], + tmp: (arr3.clone(), arr3.clone(), arr3.clone(), arr3), + } + } +} diff --git a/src/grid.rs b/src/grid.rs index eff39d1..c8a0008 100644 --- a/src/grid.rs +++ b/src/grid.rs @@ -16,9 +16,10 @@ where } impl Grid { - pub fn new(nx: usize, ny: usize, x: &[f32], y: &[f32]) -> Result { - let x = Array2::from_shape_vec((ny, nx), x.to_vec())?; - let y = Array2::from_shape_vec((ny, nx), y.to_vec())?; + pub fn new(x: Array2, y: Array2) -> Result { + assert_eq!(x.shape(), y.shape()); + let ny = x.shape()[0]; + let nx = x.shape()[1]; let mut dx_dxi = Array2::zeros((ny, nx)); SBP::diffxi(x.view(), dx_dxi.view_mut()); @@ -63,4 +64,21 @@ impl Grid { operator: std::marker::PhantomData, }) } + pub fn new_from_slice( + ny: usize, + nx: usize, + x: &[f32], + y: &[f32], + ) -> Result { + let x = Array2::from_shape_vec((ny, nx), x.to_vec())?; + let y = Array2::from_shape_vec((ny, nx), y.to_vec())?; + + Self::new(x, y) + } + pub fn nx(&self) -> usize { + self.x.shape()[1] + } + pub fn ny(&self) -> usize { + self.x.shape()[0] + } } diff --git a/src/lib.rs b/src/lib.rs index b60f2fe..835f555 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,5 +1,6 @@ use wasm_bindgen::prelude::*; +mod euler; mod grid; mod maxwell; pub mod operators; @@ -62,7 +63,7 @@ impl System { assert_eq!((width * height), x.len()); assert_eq!((width * height), y.len()); - let grid = Grid::new(width, height, x, y).expect( + let grid = Grid::new_from_slice(height, width, x, y).expect( "Could not create grid. Different number of elements compared to width*height?", ); Self { @@ -119,3 +120,132 @@ 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); + +impl EulerUniverse { + pub fn new(x: ndarray::Array2, y: ndarray::Array2) -> Self { + Self(EulerSystem::new(x, y)) + } +} + +#[wasm_bindgen] +impl EulerUniverse { + #[wasm_bindgen(constructor)] + 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)) + } + + pub fn init(&mut self, x0: f32, y0: f32) { + // Should parametrise such that we have radius, drop in pressure at center, etc + let rstar = 0.5; + let eps = 1.0; + #[allow(non_snake_case)] + let M = 0.1; + + 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) - t; + 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 = rho.powf(GAMMA) * p_inf; + 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) { + self.0.advance(dt) + } + + pub fn advance_upwind(&mut self, _dt: f32) { + todo!() + } + + pub fn get_rho_ptr(&self) -> *const u8 { + self.0.sys.0.rho().as_ptr() as *const u8 + } + pub fn get_rhou_ptr(&self) -> *const u8 { + self.0.sys.0.rhou().as_ptr() as *const u8 + } + pub fn get_rhov_ptr(&self) -> *const u8 { + self.0.sys.0.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( + &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); + } +} + +#[test] +fn start_and_advance_euler() { + let x = ndarray::Array2::from_shape_fn((20, 10), |(_j, i)| i as f32 / (10 - 1) as f32); + let y = ndarray::Array2::from_shape_fn((20, 10), |(j, _i)| j as f32 / (20 - 1) as f32); + let mut universe = EulerUniverse::new(x, y); + universe.init(0.5, 0.5); + universe.advance(0.01); +}