more work on getting the Riemann solver to fully use Castro stuff

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;
 }

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
 

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

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

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;
         }
 

Util/exact_riemann/riemann_sample.H

@@ -3,7 +3,7 @@
 
 #include <AMReX_REAL.H>
 
-#include <riemann_support.H>
+#include <riemann_constants.H>
 #include <riemann_rarefaction.H>
 
 using namespace amrex::literals;

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) *

Util/exact_riemann/riemann_star_state.H

@@ -3,7 +3,7 @@
 
 #include <AMReX_REAL.H>
 
-#include <riemann_support.H>
+#include <riemann_constants.H>
 #include <riemann_shock.H>
 #include <riemann_rarefaction.H>
 
@@ -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++;
     }