From 0a6ce88d0ec1f8fe4dabba82c5c0c35e11900083 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 1 Jul 2024 07:56:06 -0400 Subject: [PATCH] more work on getting the Riemann solver to fully use Castro stuff --- Source/hydro/riemann_constants.H | 3 +++ Util/exact_riemann/Make.package | 1 - Util/exact_riemann/inputs.test1.helm | 1 + Util/exact_riemann/inputs.test2.helm | 2 ++ Util/exact_riemann/riemann_rarefaction.H | 14 +++++--------- Util/exact_riemann/riemann_sample.H | 2 +- Util/exact_riemann/riemann_shock.H | 2 +- Util/exact_riemann/riemann_star_state.H | 11 ++++++----- Util/exact_riemann/riemann_support.H | 14 -------------- 9 files changed, 19 insertions(+), 31 deletions(-) delete mode 100644 Util/exact_riemann/riemann_support.H diff --git a/Source/hydro/riemann_constants.H b/Source/hydro/riemann_constants.H index 65ca6641c8..ad8fa2a141 100644 --- a/Source/hydro/riemann_constants.H +++ b/Source/hydro/riemann_constants.H @@ -9,6 +9,9 @@ namespace riemann_constants { constexpr amrex::Real smlp1 = 1.e-10_rt; constexpr amrex::Real small = 1.e-8_rt; constexpr amrex::Real smallu = 1.e-12_rt; + constexpr amrex::Real riemann_integral_tol = 1.e-8_rt; + constexpr amrex::Real riemann_u_tol = 1.e-6_rt; + constexpr amrex::Real riemann_p_tol = 1.e-8_rt; constexpr int HISTORY_SIZE=40; constexpr int PSTAR_BISECT_FACTOR = 5; } diff --git a/Util/exact_riemann/Make.package b/Util/exact_riemann/Make.package index 2cc1befab7..4916ffc34c 100644 --- a/Util/exact_riemann/Make.package +++ b/Util/exact_riemann/Make.package @@ -3,7 +3,6 @@ CEXE_sources += main.cpp CEXE_headers += exact_riemann.H CEXE_headers += riemann_star_state.H CEXE_headers += riemann_sample.H -CEXE_headers += riemann_support.H CEXE_sources += extern_parameters.cpp CEXE_headers += extern_parameters.H diff --git a/Util/exact_riemann/inputs.test1.helm b/Util/exact_riemann/inputs.test1.helm index b36d2f0bab..2810e90046 100644 --- a/Util/exact_riemann/inputs.test1.helm +++ b/Util/exact_riemann/inputs.test1.helm @@ -19,3 +19,4 @@ problem.t = 0.0008 problem.npts = 128 castro.T_guess = 1.e5 +castro.riemann_shock_maxiter = 20 \ No newline at end of file diff --git a/Util/exact_riemann/inputs.test2.helm b/Util/exact_riemann/inputs.test2.helm index 4a96e631a3..79bcd27016 100644 --- a/Util/exact_riemann/inputs.test2.helm +++ b/Util/exact_riemann/inputs.test2.helm @@ -15,3 +15,5 @@ problem.xjump = 0.5e5 problem.t = 0.00008 problem.npts = 128 + +castro.riemann_shock_maxiter = 100 diff --git a/Util/exact_riemann/riemann_rarefaction.H b/Util/exact_riemann/riemann_rarefaction.H index fd2e266342..8cdefed510 100644 --- a/Util/exact_riemann/riemann_rarefaction.H +++ b/Util/exact_riemann/riemann_rarefaction.H @@ -174,8 +174,6 @@ rarefaction(const amrex::Real pstar, 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; @@ -206,7 +204,7 @@ rarefaction(const amrex::Real pstar, // take a step - while (rel_error > tol) { + while (rel_error > riemann_constants::riemann_integral_tol) { dp = dp_new; if (p + dp < pstar) { dp = pstar - p; @@ -232,7 +230,7 @@ rarefaction(const amrex::Real pstar, // adaptive step estimate - double dp_est = S1 * dp * std::pow(std::abs(tol/rel_error), 0.2); + double dp_est = S1 * dp * std::pow(std::abs(riemann_constants::riemann_integral_tol/rel_error), 0.2); dp_new = dp_est; //std::clamp(S1*dp_est, dp/S2, S2*dp); } @@ -295,8 +293,6 @@ rarefaction_to_u(const amrex::Real rho_s, const amrex::Real u_s, const amrex::Re const double S1{0.9}; - constexpr amrex::Real tol = 1.e-8; - // initial conditions amrex::Real tau = 1.0_rt / rho_s; @@ -343,7 +339,7 @@ rarefaction_to_u(const amrex::Real rho_s, const amrex::Real u_s, const amrex::Re // take a step - while (rel_error > tol) { + while (rel_error > riemann_constants::riemann_integral_tol) { du = du_new; // take 2 half steps @@ -368,7 +364,7 @@ rarefaction_to_u(const amrex::Real rho_s, const amrex::Real u_s, const amrex::Re amrex::Real du_est{}; if (rel_error > 0) { - du_est = S1 * du * std::pow(std::abs(tol / rel_error), 0.2); + du_est = S1 * du * std::pow(std::abs(riemann_constants::riemann_integral_tol / rel_error), 0.2); } else { du_est = 10.0 * du; } @@ -417,7 +413,7 @@ rarefaction_to_u(const amrex::Real rho_s, const amrex::Real u_s, const amrex::Re } } - if (std::abs(du_new) < riemann_solver::tol * std::abs(u) + riemann_solver::tol) { + if (std::abs(du_new) < riemann_constants::riemann_u_tol * std::abs(u) + riemann_constants::riemann_u_tol) { finished = true; } diff --git a/Util/exact_riemann/riemann_sample.H b/Util/exact_riemann/riemann_sample.H index 8f0350fa4c..d123fd046b 100644 --- a/Util/exact_riemann/riemann_sample.H +++ b/Util/exact_riemann/riemann_sample.H @@ -3,7 +3,7 @@ #include -#include +#include #include using namespace amrex::literals; diff --git a/Util/exact_riemann/riemann_shock.H b/Util/exact_riemann/riemann_shock.H index a6bb04eb9e..ff21ea5932 100644 --- a/Util/exact_riemann/riemann_shock.H +++ b/Util/exact_riemann/riemann_shock.H @@ -134,7 +134,7 @@ shock(const amrex::Real pstar, // just doesn't work. In this case, we use the ideas from CG, Eq. 35, and // take W = sqrt(gamma p rho) - if (pstar - p_s < riemann_solver::tol_p * p_s) { + if (pstar - p_s < riemann_constants::riemann_p_tol * p_s) { W_s = std::sqrt(eos_state.gam1 * p_s * rho_s); } else { W_s = std::sqrt((pstar - p_s) * diff --git a/Util/exact_riemann/riemann_star_state.H b/Util/exact_riemann/riemann_star_state.H index 2448247129..2eb83221d5 100644 --- a/Util/exact_riemann/riemann_star_state.H +++ b/Util/exact_riemann/riemann_star_state.H @@ -3,7 +3,7 @@ #include -#include +#include #include #include @@ -108,7 +108,7 @@ riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex:: bool converged = false; int iter = 1; - while (! converged && iter < riemann_solver::max_iters) { + while (! converged && iter < castro::riemann_shock_maxiter) { // compute Z_l, W_l and Z_r, W_r -- the form of these depend // on whether the wave is a shock or a rarefaction @@ -146,13 +146,14 @@ riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex:: amrex::Real err1 = std::abs(ustar_r - ustar_l); amrex::Real err2 = pstar_new - pstar; - if (err1 < tol * amrex::max(std::abs(ustar_l), std::abs(ustar_r)) && - err2 < tol * pstar) { + std::cout << "iter = " << iter << " " << err1 / std::max(std::abs(ustar_l), std::abs(ustar_r)) << " " << err2 / pstar << std::endl; + + if (err2 < tol * pstar) { converged = true; } // get ready for the next iteration - pstar = pstar_new; + pstar = std::clamp(pstar_new, 0.5 * pstar, 2.0 * pstar); iter++; } diff --git a/Util/exact_riemann/riemann_support.H b/Util/exact_riemann/riemann_support.H deleted file mode 100644 index 56fd6b8aaf..0000000000 --- a/Util/exact_riemann/riemann_support.H +++ /dev/null @@ -1,14 +0,0 @@ -#ifndef RIEMANN_SUPPORT_H -#define RIEMANN_SUPPORT_H - -#include - -namespace riemann_solver { - const int max_iters = 100; - - const amrex::Real smallrho = 1.e-5_rt; - const amrex::Real tol = 1.e-6_rt; - const amrex::Real tol_p = 1.e-6_rt; -} - -#endif