diff --git a/euler/Cargo.toml b/euler/Cargo.toml index a8ead29..e3b6f4d 100644 --- a/euler/Cargo.toml +++ b/euler/Cargo.toml @@ -16,6 +16,7 @@ sbp = { path = "../sbp" } arrayvec = "0.6.0" serde = { version = "1.0.115", default-features = false, optional = true, features = ["derive"] } integrate = { path = "../utils/integrate" } +once_cell = "1.7.2" [dev-dependencies] criterion = "0.3.2" diff --git a/euler/src/eval.rs b/euler/src/eval.rs index 74f6140..5e4c7cc 100644 --- a/euler/src/eval.rs +++ b/euler/src/eval.rs @@ -85,10 +85,11 @@ impl<'a, D: Dimension, BP: EvaluatorPressure> Evaluator eva.u(t, x.view(), y.view(), rho.view(), rhou.view_mut()); eva.v(t, x.view(), y.view(), rho.view(), rhov.view_mut()); eva.p(t, x, y, rho.view(), rhou.view(), rhov.view(), e.view_mut()); + let gamma = *GAMMA.get().expect("GAMMA is not defined"); azip!((rho in &rho, u in &rhou, v in &rhov, e in &mut e) { let p = *e; - *e = p / (GAMMA - 1.0) + rho * (u*u + v*v) / 2.0; + *e = p / (gamma - 1.0) + rho * (u*u + v*v) / 2.0; }); azip!((rho in &rho, rhou in &mut rhou) *rhou *= rho); diff --git a/euler/src/lib.rs b/euler/src/lib.rs index 7acda40..b1b6242 100644 --- a/euler/src/lib.rs +++ b/euler/src/lib.rs @@ -1,6 +1,7 @@ pub use arrayvec::ArrayVec; use ndarray::azip; use ndarray::prelude::*; +use once_cell::sync::OnceCell; use sbp::grid::{Grid, Metrics}; use sbp::operators::{InterpolationOperator, SbpOperator2d, UpwindOperator2d}; use sbp::utils::Direction; @@ -11,7 +12,7 @@ use eval::Evaluator; mod vortex; pub use vortex::{VortexParameters, Vortice}; -pub const GAMMA: Float = 1.4; +pub static GAMMA: OnceCell = OnceCell::new(); // A collection of buffers that allows one to efficiently // move to the next state @@ -566,6 +567,7 @@ fn upwind_dissipation( metrics: &Metrics, tmp: (&mut Field, &mut Field), ) { + let gamma = *GAMMA.get().expect("GAMMA is not defined"); for (((FieldValue { rho, rhou, rhov, e }, tmp0), tmp1), metric) in y .iter() .zip(tmp.0.iter_mut()) @@ -580,9 +582,9 @@ fn upwind_dissipation( let uhat = metric.detj_dxi_dx * u + metric.detj_dxi_dy * v; let vhat = metric.detj_deta_dx * u + metric.detj_deta_dy * v; - let p = pressure(GAMMA, rho, rhou, rhov, e); + let p = pressure(gamma, rho, rhou, rhov, e); assert!(p > 0.0); - let c = (GAMMA * p / rho).sqrt(); + let c = (gamma * p / rho).sqrt(); // The accurate hypot is very slow, and the accuracy is // not that important in this case @@ -625,9 +627,11 @@ fn fluxes(k: (&mut Field, &mut Field), y: &Field, metrics: &Metrics, wb: &mut Fi let rhov = y.rhov(); let e = y.e(); + let gamma = *GAMMA.get().expect("GAMMA is not defined"); + let mut p = wb.rho_mut(); azip!((p in &mut p, &rho in &rho, &rhou in &rhou, &rhov in &rhov, &e in &e) { - *p = pressure(GAMMA, rho, rhou, rhov, e) + *p = pressure(gamma, rho, rhou, rhov, e) }); let (mut c0, c1, mut c2, c3) = k.0.components_mut(); @@ -1059,6 +1063,7 @@ fn SAT_characteristic( assert_eq!(y.shape(), z.shape()); assert_eq!(y.shape()[0], 4); assert_eq!(y.shape()[1], detj.shape()[0]); + let gamma = *GAMMA.get().expect("GAMMA is not defined"); for (((((mut k, y), z), detj), detj_d_dx), detj_d_dy) in k .axis_iter_mut(ndarray::Axis(1)) @@ -1086,19 +1091,19 @@ fn SAT_characteristic( 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 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 alpha = rho / (sbp::consts::SQRT_2 * c); - let phi2_c2 = (phi2 + c * c) / (GAMMA - 1.0); + let phi2_c2 = (phi2 + c * c) / (gamma - 1.0); #[rustfmt::skip] let T = [ [ 1.0, 0.0, alpha, alpha], [ u, ky, alpha*(u + kx * c), alpha*(u - kx * c)], [ v, -kx, alpha*(v + ky * c), alpha*(v - ky * c)], - [phi2 / (GAMMA - 1.0), rho*(ky * u - kx * v), alpha*(phi2_c2 + c * theta), alpha*(phi2_c2 - c * theta)], + [phi2 / (gamma - 1.0), rho*(ky * u - kx * v), alpha*(phi2_c2 + c * theta), alpha*(phi2_c2 - c * theta)], ]; let U = kx_ * u + ky_ * v; let L = [ @@ -1110,10 +1115,10 @@ fn SAT_characteristic( let beta = 1.0 / (2.0 * c * c); #[rustfmt::skip] let TI = [ - [ 1.0 - phi2 / (c * c), (GAMMA - 1.0) * u / (c * c), (GAMMA - 1.0) * v / (c * c), -(GAMMA - 1.0) / (c * c)], + [ 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)/rho, ky/rho, -kx/rho, 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)], + [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], rhou - z[1], rhov - z[2], e - z[3]]; diff --git a/euler/src/vortex.rs b/euler/src/vortex.rs index c4262e2..ea5714c 100644 --- a/euler/src/vortex.rs +++ b/euler/src/vortex.rs @@ -33,13 +33,14 @@ impl eval::Evaluator for VortexParameters { rhov: ArrayViewMut, e: ArrayViewMut, ) { + let gamma = *GAMMA.get().expect("GAMMA is not defined"); let m = self.mach; - let p_inf = 1.0 / (GAMMA * m * m); + let p_inf = 1.0 / (gamma * m * m); let rho_inf: Float = 1.0; let u_inf: Float = 1.0; let v_inf: Float = 0.0; - let e_inf = p_inf / (GAMMA - 1.0) + rho_inf * (u_inf.powi(2) + v_inf.powi(2)) / 2.0; + let e_inf = p_inf / (gamma - 1.0) + rho_inf * (u_inf.powi(2) + v_inf.powi(2)) / 2.0; azip!((rho in rho, rhou in rhou, @@ -69,9 +70,9 @@ impl eval::Evaluator for VortexParameters { let dy = y - vortice.y0; let f = (1.0 - (dx*dx + dy*dy))/(rstar*rstar); - *rho = Float::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)); + *rho = Float::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 = Float::powf(*rho, GAMMA)*p_inf; + let p = Float::powf(*rho, gamma)*p_inf; let u = 1.0 - eps*dy/(2.0*PI*p_inf.sqrt()*rstar*rstar)*(f/2.0).exp(); let v = eps*dx/(2.0*PI*p_inf.sqrt()*rstar*rstar)*(f/2.0).exp(); @@ -80,7 +81,7 @@ impl eval::Evaluator for VortexParameters { assert!(*rho > 0.0); *rhou = *rho*u; *rhov = *rho*v; - *e = p/(GAMMA - 1.0) + *rho*(u*u + v*v)/2.0; + *e = p/(gamma - 1.0) + *rho*(u*u + v*v)/2.0; } } @@ -94,8 +95,8 @@ impl eval::Evaluator for VortexParameters { let dy = y - vortice.y0; let f = (1.0 - (dx*dx + dy*dy))/(rstar*rstar); - let rho_vortice = Float::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)); - let p = Float::powf(rho_vortice, GAMMA)*p_inf; + let rho_vortice = Float::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)); + let p = Float::powf(rho_vortice, gamma)*p_inf; let u = 1.0 - eps*dy/(2.0*PI*p_inf.sqrt()*rstar*rstar)*(f/2.0).exp(); let v = eps*dx/(2.0*PI*p_inf.sqrt()*rstar*rstar)*(f/2.0).exp(); @@ -106,7 +107,7 @@ impl eval::Evaluator for VortexParameters { assert!(*rho > 0.0); *rhou += rho_vortice*u - rho_inf*u_inf; *rhov += rho_vortice*v - rho_inf*v_inf; - *e += (p/(GAMMA - 1.0) + rho_vortice*(u*u + v*v)/2.0) - e_inf; + *e += (p/(gamma - 1.0) + rho_vortice*(u*u + v*v)/2.0) - e_inf; } }); } diff --git a/multigrid/src/eval.rs b/multigrid/src/eval.rs index 28be21f..7fe4c38 100644 --- a/multigrid/src/eval.rs +++ b/multigrid/src/eval.rs @@ -250,8 +250,8 @@ fn append_default_context() { pub fn default_context() -> HashMapContext { let mut context = math_consts_context! {}.unwrap(); - - context.set_value("GAMMA".into(), GAMMA.into()).unwrap(); + let gamma = *GAMMA.get().expect("GAMMA is not defined"); + context.set_value("GAMMA".into(), gamma.into()).unwrap(); context .set_function( diff --git a/multigrid/src/parsing.rs b/multigrid/src/parsing.rs index ffd8e9d..d9d518e 100644 --- a/multigrid/src/parsing.rs +++ b/multigrid/src/parsing.rs @@ -283,6 +283,10 @@ pub enum BoundaryConditions { NotNeeded, } +fn default_gamma() -> Float { + 1.4 +} + #[derive(Clone, Debug, Serialize, Deserialize)] /// Input configuration (json) pub struct Configuration { @@ -291,6 +295,8 @@ pub struct Configuration { pub initial_conditions: InputInitialConditions, #[serde(default)] pub boundary_conditions: InputBoundaryConditions, + #[serde(default = "default_gamma")] + pub gamma: Float, } pub struct RuntimeConfiguration { @@ -305,6 +311,8 @@ pub struct RuntimeConfiguration { impl Configuration { pub fn into_runtime(mut self) -> RuntimeConfiguration { + let gamma = self.gamma; + let _ = euler::GAMMA.set(gamma); let default = self.grids.shift_remove("default").unwrap_or_default(); let names = self.grids.keys().cloned().collect(); diff --git a/webfront/src/euler.rs b/webfront/src/euler.rs index 0cd926f..7d00db2 100644 --- a/webfront/src/euler.rs +++ b/webfront/src/euler.rs @@ -14,6 +14,7 @@ impl EulerUniverse { impl EulerUniverse { #[wasm_bindgen(constructor)] pub fn new_with_slice(height: usize, width: usize, x: &[f32], y: &[f32]) -> Self { + let _ = euler::GAMMA.set(1.4); 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(euler::System::new(x, y, operators::Upwind4))