Finish riemann cleaning

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

@@ -1,9 +1,10 @@
 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_headers += exact_riemann_star_state.H
+CEXE_headers += exact_riemann_sample.H
+CEXE_headers += exact_riemann_shock.H
+CEXE_headers += exact_riemann_rarefaction.H
 CEXE_sources += extern_parameters.cpp
 CEXE_headers += extern_parameters.H
 

Util/exact_riemann/exact_riemann.H

@@ -7,8 +7,8 @@
 
 #include <castro_params.H>
 
-#include <riemann_sample.H>
-#include <riemann_star_state.H>
+#include <exact_riemann_sample.H>
+#include <exact_riemann_star_state.H>
 
 using namespace amrex::literals;
 
@@ -50,10 +50,13 @@ exact_riemann() {
 
     amrex::Real ustar, pstar, W_l, W_r;
 
-    riemann_star_state(problem::rho_l, problem::u_l, problem::p_l, xn,
-                       problem::rho_r, problem::u_r, problem::p_r, xn,
-                       ustar, pstar, W_l, W_r);
+    std::cout << "solving for star state: " << std::endl;
+
+    auto niter = riemann_star_state(problem::rho_l, problem::u_l, problem::p_l, xn,
+                                    problem::rho_r, problem::u_r, problem::p_r, xn,
+                                    ustar, pstar, W_l, W_r);
 
+    std::cout << "found star state, niter = " << niter << "; pstar, ustar = " << pstar << " " << ustar << std::endl;
 
     // find the solution as a function of xi = x/t
     std::ofstream ofile;

Util/exact_riemann/riemann_rarefaction.H → Util/exact_riemann/exact_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 → Util/exact_riemann/exact_riemann_sample.H

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

Util/exact_riemann/riemann_shock.H → Util/exact_riemann/exact_riemann_shock.H

@@ -40,18 +40,18 @@ W_s_shock(const amrex::Real W_s, const amrex::Real pstar,
     amrex::Real dedrho_p = eos_state.dedr - eos_state.dedT * eos_state.dpdr / eos_state.dpdT;
 
     fprime = 2.0_rt * W_s * (eos_state.e - e_s) -
-        2.0_rt * dedrho_p * (pstar - p_s) * std::pow(rhostar_s, 2) / W_s;
+        2.0_rt * dedrho_p * (pstar - p_s) * amrex::Math::powi<2>(rhostar_s) / W_s;
 
 }
 
 
 AMREX_INLINE
-void
+bool
 newton_shock(amrex::Real& W_s, const amrex::Real pstar,
              const amrex::Real rho_s, const amrex::Real p_s, const amrex::Real e_s,
              const amrex::Real* xn,
              amrex::Real* rhostar_hist, amrex::Real* Ws_hist,
-             eos_rep_t& eos_state, bool& converged) {
+             eos_rep_t& eos_state) {
 
     // Newton iterations -- we are zeroing the energy R-H jump condition
     // W^2 [e] = 1/2 [p^2]
@@ -60,7 +60,7 @@ newton_shock(amrex::Real& W_s, const amrex::Real pstar,
     //
     // and then compute f'
 
-    converged = false;
+    bool converged = false;
 
     int iter = 1;
     while (! converged && iter < castro::riemann_shock_maxiter) {
@@ -76,7 +76,7 @@ newton_shock(amrex::Real& W_s, const amrex::Real pstar,
             converged = true;
         }
 
-        W_s = amrex::min(2.0_rt * W_s, amrex::max(0.5_rt * W_s, W_s + dW));
+        W_s = std::clamp(W_s + dW, 0.5_rt * W_s, 2.0_rt * W_s);
 
         // store some history
 
@@ -85,6 +85,8 @@ newton_shock(amrex::Real& W_s, const amrex::Real pstar,
 
         iter++;
     }
+
+    return converged;
 }
 
 AMREX_INLINE
@@ -134,7 +136,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) *
@@ -154,10 +156,9 @@ shock(const amrex::Real pstar,
 
     // newton
 
-    bool converged;
-    newton_shock(W_s, pstar, rho_s, p_s, e_s, xn,
-                 rhostar_hist, Ws_hist,
-                 eos_state, converged);
+    auto converged = newton_shock(W_s, pstar, rho_s, p_s, e_s, xn,
+                                  rhostar_hist, Ws_hist,
+                                  eos_state);
 
 
     // now did we converge?
@@ -185,7 +186,7 @@ shock(const amrex::Real pstar,
     amrex::Real p_e = eos_state.dpdT / eos_state.dedT;
     amrex::Real p_rho = eos_state.dpdr - eos_state.dpdT * eos_state.dedr / eos_state.dedT;
 
-    amrex::Real p_tau = -std::pow(rhostar_s, 2) * p_rho;
+    amrex::Real p_tau = -amrex::Math::powi<2>(rhostar_s) * p_rho;
 
     amrex::Real dW2dpstar = (C*C - W_s*W_s) * W_s * W_s /
         ((0.5_rt * (pstar + p_s) * p_e - p_tau) * (pstar - p_s));

Util/exact_riemann/riemann_star_state.H → Util/exact_riemann/exact_riemann_star_state.H

@@ -3,14 +3,14 @@
 
 #include <AMReX_REAL.H>
 
-#include <riemann_support.H>
-#include <riemann_shock.H>
-#include <riemann_rarefaction.H>
+#include <riemann_constants.H>
+#include <exact_riemann_shock.H>
+#include <exact_riemann_rarefaction.H>
 
 using namespace amrex::literals;
 
 AMREX_INLINE
-void
+int
 riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex::Real p_l, const amrex::Real* xn_l,
                    const amrex::Real rho_r, const amrex::Real u_r, const amrex::Real p_r, const amrex::Real* xn_r,
                    amrex::Real& ustar, amrex::Real& pstar, amrex::Real& W_l, amrex::Real& W_r) {
@@ -74,13 +74,11 @@ riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex::
         pstar = ((W_r * p_l + W_l * p_r) + W_l * W_r * (u_l - u_r)) / (W_l + W_r);
     }
 
-    pstar = amrex::max(pstar, smallp);
+    pstar = std::max(pstar, smallp);
 
 
     // find the exact pstar and ustar
 
-    std::cout << "solving for star state: " << std::endl;
-
     // this procedure follows directly from Colella & Glaz 1985, section 1
     // the basic idea is that we want to find the pstar that satisfies:
     //
@@ -108,7 +106,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
@@ -151,7 +149,7 @@ riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex::
         }
 
         // get ready for the next iteration
-        pstar = pstar_new;
+        pstar = std::clamp(pstar_new, 0.5 * pstar, 2.0 * pstar);
 
         iter++;
     }
@@ -162,6 +160,7 @@ riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex::
 
     ustar = 0.5_rt * (ustar_l + ustar_r);
 
-    std::cout << "found pstar, ustar: " <<  pstar << " " << ustar << std::endl;
+    return iter;
+
 }
 #endif

Util/exact_riemann/inputs.bug8

@@ -16,3 +16,4 @@ problem.t     = 0.0008
 
 problem.npts  = 512
 
+castro.riemann_shock_maxiter = 100