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

use symmetrisation from [0] 

[0]: Two-Dimensional Inlet Simulation Using a
    Diagonal Implicit Algorithm, Chausee & Pulliam, AIAA Journal
This commit is contained in:
Magnus Ulimoen 2020-09-15 18:12:27 +02:00
parent 5161dd8ab6
commit ac83a2cbca
1 changed files with 27 additions and 29 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
euler/src/lib.rs
View File

@ -610,14 +610,11 @@ fn upwind_dissipation(
let mut tmp0 = tmp.0.view_mut().into_shape((4, n)).unwrap();
let mut tmp1 = tmp.1.view_mut().into_shape((4, n)).unwrap();
for (
((((((y, mut tmp0), mut tmp1), detj), detj_dxi_dx), detj_dxi_dy), detj_deta_dx),
detj_deta_dy,
) in yview
for ((((((y, mut tmp0), mut tmp1), detj_dxi_dx), detj_dxi_dy), detj_deta_dx), detj_deta_dy) in
yview
.axis_iter(ndarray::Axis(1))
.zip(tmp0.axis_iter_mut(ndarray::Axis(1)))
.zip(tmp1.axis_iter_mut(ndarray::Axis(1)))
.zip(metrics.detj().iter())
.zip(metrics.detj_dxi_dx().iter())
.zip(metrics.detj_dxi_dy().iter())
.zip(metrics.detj_deta_dx().iter())
@ -632,27 +629,27 @@ fn upwind_dissipation(
let u = rhou / rho;
let v = rhov / rho;
let uhat = detj_dxi_dx / detj * u + detj_dxi_dy / detj * v;
let vhat = detj_deta_dx / detj * u + detj_deta_dy / detj * v;
let uhat = detj_dxi_dx * u + detj_dxi_dy * v;
let vhat = detj_deta_dx * u + detj_deta_dy * v;
let p = pressure(GAMMA, rho, rhou, rhov, e);
assert!(p > 0.0);
let c = (GAMMA * p / rho).sqrt();
let alpha_u = uhat.abs() + c;
let alpha_v = vhat.abs() + c;
let alpha_u = uhat.abs() + c * Float::hypot(*detj_dxi_dx, *detj_dxi_dy);
let alpha_v = vhat.abs() + c * Float::hypot(*detj_deta_dx, *detj_deta_dy);
tmp0[0] = alpha_u * rho * detj;
tmp1[0] = alpha_v * rho * detj;
tmp0[0] = alpha_u * rho;
tmp1[0] = alpha_v * rho;
tmp0[1] = alpha_u * rhou * detj;
tmp1[1] = alpha_v * rhou * detj;
tmp0[1] = alpha_u * rhou;
tmp1[1] = alpha_v * rhou;
tmp0[2] = alpha_u * rhov * detj;
tmp1[2] = alpha_v * rhov * detj;
tmp0[2] = alpha_u * rhov;
tmp1[2] = alpha_v * rhov;
tmp0[3] = alpha_u * e * detj;
tmp1[3] = alpha_v * e * detj;
tmp0[3] = alpha_u * e;
tmp1[3] = alpha_v * e;
}
op.dissxi(tmp.0.rho(), k.0.rho_mut());
@ -1099,15 +1096,16 @@ fn SAT_characteristic(
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);
#[rustfmt::skip]
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],
[ 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)],
];
let U = kx_ * u + ky_ * v;
let L = [
@ -1120,7 +1118,7 @@ fn SAT_characteristic(
#[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)],
[ -(ky * u - kx * v), ky, -kx, 0.0],
[ -(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)],
];