From 0fc841d58168830aeb4fb84a5ff64bf0cdce1941 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 7 Jun 2024 12:50:17 -0400 Subject: [PATCH] a bit of cleaning --- Util/exact_riemann/riemann_rarefaction.H | 77 +++++++++++++++--------- 1 file changed, 48 insertions(+), 29 deletions(-) diff --git a/Util/exact_riemann/riemann_rarefaction.H b/Util/exact_riemann/riemann_rarefaction.H index e1b1e41eaa..717aba4873 100644 --- a/Util/exact_riemann/riemann_rarefaction.H +++ b/Util/exact_riemann/riemann_rarefaction.H @@ -95,9 +95,12 @@ riemann_invariant_rhs2(const amrex::Real u, const amrex::Real tau, const amrex:: AMREX_INLINE std::pair -single_step(const amrex::Real pstar0, const amrex::Real dp, - const amrex::Real u0, const amrex::Real tau0, - const amrex::Real* xn, const int iwave, amrex::Real& T) { +single_step_p(const amrex::Real pstar0, const amrex::Real dp, + const amrex::Real u0, const amrex::Real tau0, + const amrex::Real* xn, const int iwave, amrex::Real& T) { + + // this takes a single step of the Riemann invariant system where + // p is the independent variable amrex::Real dtaudp1, dudp1; riemann_invariant_rhs(pstar0, tau0, u0, xn, iwave, T, dtaudp1, dudp1); @@ -121,6 +124,37 @@ single_step(const amrex::Real pstar0, const amrex::Real dp, } +AMREX_INLINE +std::pair +single_step_u(const amrex::Real u0, const amrex::Real du, + const amrex::Real tau0, const amrex::Real p0, + const amrex::Real* xn, const int iwave, amrex::Real& T) { + + // this takes a single step of the Riemann invariant system where + // u is the independent variable + + amrex::Real dtaudu1, dpdu1; + riemann_invariant_rhs2(u0, tau0, p0, xn, iwave, T, dtaudu1, dpdu1); + + amrex::Real dtaudu2, dpdu2; + riemann_invariant_rhs2(u0+0.5*du, tau0+0.5*du*dtaudu1, p0+0.5*du*dpdu1, + xn, iwave, T, dtaudu2, dpdu2); + + amrex::Real dtaudu3, dpdu3; + riemann_invariant_rhs2(u0+0.5*du, tau0+0.5*du*dtaudu2, p0+0.5*du*dpdu2, + xn, iwave, T, dtaudu3, dpdu3); + + amrex::Real dtaudu4, dpdu4; + riemann_invariant_rhs2(u0+du, tau0+du*dtaudu3, p0+du*dpdu3, + xn, iwave, T, dtaudu4, dpdu4); + + amrex::Real p = p0 + (1.0_rt/6.0_rt) * du * (dpdu1 + 2.0_rt * dpdu2 + 2.0_rt * dpdu3 + dpdu4); + amrex::Real tau = tau0 + (1.0_rt/6.0_rt) * du * (dtaudu1 + 2.0_rt * dtaudu2 + 2.0_rt * dtaudu3 + dtaudu4); + + return {tau, p}; + +} + AMREX_INLINE void rarefaction(const amrex::Real pstar, @@ -128,9 +162,6 @@ rarefaction(const amrex::Real pstar, const amrex::Real* xn, const int iwave, amrex::Real& Z_s, amrex::Real& W_s, amrex::Real& rhostar) { - const double S1{0.9}; - const double S2{4.0}; - // Compute Z_s = C(p*, rho*) for a rarefaction connecting the // state to the star region by integrating the Riemann invariant // from p_s to pstar. @@ -140,8 +171,13 @@ rarefaction(const amrex::Real pstar, // dtau/dp = -1/C**2 // du/dp = +/- 1/C (+ for 1-wave, - for 3-wave) + const double S1{0.9}; + const double S2{4.0}; + constexpr amrex::Real tol = 1.e-5; + // initial conditions + amrex::Real tau = 1.0_rt / rho_s; amrex::Real u = u_s; amrex::Real p = p_s; @@ -168,6 +204,8 @@ rarefaction(const amrex::Real pstar, amrex::Real rel_error = std::numeric_limits::max(); + // take a step + while (rel_error > tol) { dp = dp_new; if (p + dp < pstar) { @@ -177,11 +215,11 @@ rarefaction(const amrex::Real pstar, // take 2 half steps amrex::Real u_tmp; amrex::Real tau_tmp; - std::tie(u_tmp, tau_tmp) = single_step(p, 0.5*dp, u0, tau0, xn, iwave, T); - std::tie(u_new, tau_new) = single_step(p + 0.5*dp, 0.5*dp, u_tmp, tau_tmp, xn, iwave, T); + std::tie(u_tmp, tau_tmp) = single_step_p(p, 0.5*dp, u0, tau0, xn, iwave, T); + std::tie(u_new, tau_new) = single_step_p(p + 0.5*dp, 0.5*dp, u_tmp, tau_tmp, xn, iwave, T); // now take a single step to cover dp - auto [u_single, tau_single] = single_step(p, dp, u0, tau0, xn, iwave, T); + auto [u_single, tau_single] = single_step_p(p, dp, u0, tau0, xn, iwave, T); // estimate the relative error amrex::Real u_err = std::abs(u_new - u_single); @@ -192,10 +230,6 @@ rarefaction(const amrex::Real pstar, rel_error = std::max(u_err, tau_err); - //if (ntry > 10) { - // std::exit(-1); - //} - // adaptive step estimate double dp_est = S1 * dp * std::pow(std::abs(tol/rel_error), 0.2); @@ -292,29 +326,14 @@ rarefaction_to_u(const amrex::Real rho_s, const amrex::Real u_s, const amrex::Re std::cout << "integrating from u: " << u << " " << ustop << " " << xi << " " << c << std::endl; - amrex::Real du2 = 0.5_rt * du; - amrex::Real T = eos_state.T; while (! finished) { // do 4th-order RT - amrex::Real dtaudu1, dpdu1; - riemann_invariant_rhs2(u, tau, p, xn, iwave, T, dtaudu1, dpdu1); - - amrex::Real dtaudu2, dpdu2; - riemann_invariant_rhs2(u+du2, tau+du2*dtaudu1, p+du2*dpdu1, xn, iwave, T, dtaudu2, dpdu2); - - amrex::Real dtaudu3, dpdu3; - riemann_invariant_rhs2(u+du2, tau+du2*dtaudu2, p+du2*dpdu2, xn, iwave, T, dtaudu3, dpdu3); - - amrex::Real dtaudu4, dpdu4; - riemann_invariant_rhs2(u+du, tau+du*dtaudu3, p+du*dpdu3, xn, iwave, T, dtaudu4, dpdu4); - + std::tie(tau, p) = single_step_u(u, du, tau, p, xn, iwave, T); u += du; - p += (1.0_rt/6.0_rt) * du * (dpdu1 + 2.0_rt * dpdu2 + 2.0_rt * dpdu3 + dpdu4); - tau += (1.0_rt/6.0_rt) * du * (dtaudu1 + 2.0_rt * dtaudu2 + 2.0_rt * dtaudu3 + dtaudu4); // compute c