mulimoen
/
SummationByParts
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Simplify vortex using Evaluator trait

This commit is contained in:
Magnus Ulimoen 2021-06-30 17:35:05 +02:00
parent 92ad7bc580
commit 94e49ff9b5
2 changed files with 73 additions and 93 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
euler/src/lib.rs
View File

@ -7,9 +7,9 @@ use sbp::utils::Direction;
use sbp::Float; use sbp::Float;
pub mod eval; pub mod eval;
use eval::Evaluator;
mod vortex; mod vortex;
pub use vortex::{vortex, VortexParameters, Vortice}; pub use vortex::{VortexParameters, Vortice};
pub const GAMMA: Float = 1.4; pub const GAMMA: Float = 1.4;
@ -320,23 +320,8 @@ impl Field {
time: Float, time: Float,
vortex_param: &VortexParameters, 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 (rho, rhou, rhov, e) = self.components_mut();
let n = rho.len(); vortex_param.evaluate(time, x, y, rho, rhou, rhov, e)
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,
)
} }
fn iter(&self) -> impl ExactSizeIterator<Item = FieldValue> + '_ { fn iter(&self) -> impl ExactSizeIterator<Item = FieldValue> + '_ {
let n = self.nx() * self.ny(); let n = self.nx() * self.ny();
@ -949,7 +934,7 @@ fn vortexify(
fiter.next().unwrap(), fiter.next().unwrap(),
); );
let (y, x) = yx; 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)] #[allow(non_snake_case)]

143
euler/src/vortex.rs
View File

@ -20,52 +20,71 @@ pub struct VortexParameters {
pub mach: Float, pub mach: Float,
} }
#[allow(clippy::too_many_arguments)] impl<D: Dimension> eval::Evaluator<D> for VortexParameters {
#[allow(clippy::many_single_char_names)] #[allow(clippy::too_many_arguments)]
pub fn vortex( #[allow(clippy::many_single_char_names)]
rho: ArrayViewMut1<Float>, fn evaluate(
rhou: ArrayViewMut1<Float>, &self,
rhov: ArrayViewMut1<Float>, time: Float,
e: ArrayViewMut1<Float>, x: ArrayView<Float, D>,
x: ArrayView1<Float>, y: ArrayView<Float, D>,
y: ArrayView1<Float>, rho: ArrayViewMut<Float, D>,
time: Float, rhou: ArrayViewMut<Float, D>,
vortex_param: &VortexParameters, rhov: ArrayViewMut<Float, D>,
) { e: ArrayViewMut<Float, D>,
assert_eq!(rho.len(), rhou.len()); ) {
assert_eq!(rho.len(), rhov.len()); let m = self.mach;
assert_eq!(rho.len(), e.len()); let p_inf = 1.0 / (GAMMA * m * m);
assert_eq!(rho.len(), x.len());
assert_eq!(rho.len(), y.len());
assert_eq!(x.shape(), y.shape());
let m = vortex_param.mach; let rho_inf: Float = 1.0;
let p_inf = 1.0 / (GAMMA * m * m); 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; azip!((rho in rho,
let u_inf: Float = 1.0; rhou in rhou,
let v_inf: Float = 0.0; rhov in rhov,
let e_inf = p_inf / (GAMMA - 1.0) + rho_inf * (u_inf.powi(2) + v_inf.powi(2)) / 2.0; e in e,
x in x,
y in y)
{
azip!((rho in rho, let mut iterator = self.vortices.iter();
rhou in rhou,
rhov in rhov,
e in e,
x in x,
y in y)
{
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() { let rstar = vortice.rstar;
None => { let eps = vortice.eps;
*rho = rho_inf;
*rhou = rho_inf*u_inf; let dx = (x - vortice.x0) - time;
*rhou = rho_inf*v_inf; let dy = y - vortice.y0;
*e = e_inf; let f = (1.0 - (dx*dx + dy*dy))/(rstar*rstar);
return;
}, *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));
Some(vortice) => { 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; use sbp::consts::PI;
let rstar = vortice.rstar; let rstar = vortice.rstar;
@ -75,44 +94,20 @@ pub fn vortex(
let dy = y - vortice.y0; let dy = y - vortice.y0;
let f = (1.0 - (dx*dx + dy*dy))/(rstar*rstar); 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)); 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));
assert!(*rho > 0.0); let p = Float::powf(rho_vortice, 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 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(); 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); assert!(p > 0.0);
*rho += rho_vortice - rho_inf;
assert!(*rho > 0.0); assert!(*rho > 0.0);
*rhou = *rho*u; *rhou += rho_vortice*u - rho_inf*u_inf;
*rhov = *rho*v; *rhov += rho_vortice*v - rho_inf*v_inf;
*e = p/(GAMMA - 1.0) + *rho*(u*u + v*v)/2.0; *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;
}
});
} }