From 94e49ff9b5a811f4bb29d03a3c8c59b9e4ccf5bc Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Wed, 30 Jun 2021 17:35:05 +0200 Subject: [PATCH] Simplify vortex using Evaluator trait --- euler/src/lib.rs | 23 ++----- euler/src/vortex.rs | 143 +++++++++++++++++++++----------------------- 2 files changed, 73 insertions(+), 93 deletions(-) diff --git a/euler/src/lib.rs b/euler/src/lib.rs index b7801ab..7acda40 100644 --- a/euler/src/lib.rs +++ b/euler/src/lib.rs @@ -7,9 +7,9 @@ use sbp::utils::Direction; use sbp::Float; pub mod eval; - +use eval::Evaluator; mod vortex; -pub use vortex::{vortex, VortexParameters, Vortice}; +pub use vortex::{VortexParameters, Vortice}; pub const GAMMA: Float = 1.4; @@ -320,23 +320,8 @@ impl Field { time: Float, vortex_param: &VortexParameters, ) { - assert_eq!(x.shape(), y.shape()); - assert_eq!(x.shape()[1], self.nx()); - assert_eq!(x.shape()[0], self.ny()); - let (rho, rhou, rhov, e) = self.components_mut(); - let n = rho.len(); - - vortex( - rho.into_shape((n,)).unwrap(), - rhou.into_shape((n,)).unwrap(), - rhov.into_shape((n,)).unwrap(), - e.into_shape((n,)).unwrap(), - x.into_shape((n,)).unwrap(), - y.into_shape((n,)).unwrap(), - time, - &vortex_param, - ) + vortex_param.evaluate(time, x, y, rho, rhou, rhov, e) } fn iter(&self) -> impl ExactSizeIterator + '_ { let n = self.nx() * self.ny(); @@ -949,7 +934,7 @@ fn vortexify( fiter.next().unwrap(), ); let (y, x) = yx; - vortex(rho, rhou, rhov, e, x, y, time, &vparams); + vparams.evaluate(time, x, y, rho, rhou, rhov, e) } #[allow(non_snake_case)] diff --git a/euler/src/vortex.rs b/euler/src/vortex.rs index 4087f16..c4262e2 100644 --- a/euler/src/vortex.rs +++ b/euler/src/vortex.rs @@ -20,52 +20,71 @@ pub struct VortexParameters { pub mach: Float, } -#[allow(clippy::too_many_arguments)] -#[allow(clippy::many_single_char_names)] -pub fn vortex( - rho: ArrayViewMut1, - rhou: ArrayViewMut1, - rhov: ArrayViewMut1, - e: ArrayViewMut1, - x: ArrayView1, - y: ArrayView1, - time: Float, - vortex_param: &VortexParameters, -) { - assert_eq!(rho.len(), rhou.len()); - assert_eq!(rho.len(), rhov.len()); - assert_eq!(rho.len(), e.len()); - assert_eq!(rho.len(), x.len()); - assert_eq!(rho.len(), y.len()); - assert_eq!(x.shape(), y.shape()); +impl eval::Evaluator for VortexParameters { + #[allow(clippy::too_many_arguments)] + #[allow(clippy::many_single_char_names)] + fn evaluate( + &self, + time: Float, + x: ArrayView, + y: ArrayView, + rho: ArrayViewMut, + rhou: ArrayViewMut, + rhov: ArrayViewMut, + e: ArrayViewMut, + ) { + let m = self.mach; + let p_inf = 1.0 / (GAMMA * m * m); - let m = vortex_param.mach; - 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 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; + azip!((rho in rho, + rhou in rhou, + rhov in rhov, + e in e, + x in x, + y in y) + { - azip!((rho in rho, - rhou in rhou, - rhov in rhov, - e in e, - x in x, - y in y) - { + let mut iterator = self.vortices.iter(); - let mut iterator = vortex_param.vortices.iter(); + match iterator.next() { + None => { + *rho = rho_inf; + *rhou = rho_inf*u_inf; + *rhou = rho_inf*v_inf; + *e = e_inf; + return; + }, + Some(vortice) => { + use sbp::consts::PI; - match iterator.next() { - None => { - *rho = rho_inf; - *rhou = rho_inf*u_inf; - *rhou = rho_inf*v_inf; - *e = e_inf; - return; - }, - Some(vortice) => { + let rstar = vortice.rstar; + let eps = vortice.eps; + + let dx = (x - vortice.x0) - time; + 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)); + assert!(*rho > 0.0); + 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(); + + assert!(p > 0.0); + + assert!(*rho > 0.0); + *rhou = *rho*u; + *rhov = *rho*v; + *e = p/(GAMMA - 1.0) + *rho*(u*u + v*v)/2.0; + } + } + + for vortice in iterator { use sbp::consts::PI; let rstar = vortice.rstar; @@ -75,44 +94,20 @@ pub fn vortex( 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)); - assert!(*rho > 0.0); - let p = Float::powf(*rho, 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(); + assert!(rho_vortice > 0.0); assert!(p > 0.0); + *rho += rho_vortice - rho_inf; assert!(*rho > 0.0); - *rhou = *rho*u; - *rhov = *rho*v; - *e = p/(GAMMA - 1.0) + *rho*(u*u + v*v)/2.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; } - } - - for vortice in iterator { - use sbp::consts::PI; - - let rstar = vortice.rstar; - let eps = vortice.eps; - - let dx = (x - vortice.x0) - time; - 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 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(); - - assert!(rho_vortice > 0.0); - assert!(p > 0.0); - - *rho += rho_vortice - rho_inf; - 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; - } - }); + }); + } }