diff --git a/euler/src/lib.rs b/euler/src/lib.rs index 9de5928..13e605b 100644 --- a/euler/src/lib.rs +++ b/euler/src/lib.rs @@ -610,18 +610,15 @@ 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 - .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()) - .zip(metrics.detj_deta_dy().iter()) + 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_dxi_dx().iter()) + .zip(metrics.detj_dxi_dy().iter()) + .zip(metrics.detj_deta_dx().iter()) + .zip(metrics.detj_deta_dy().iter()) { let rho = y[0]; assert!(rho > 0.0); @@ -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)], ];