more cleaning

Util/exact_riemann/exact_riemann.H

@@ -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/inputs.bug8

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

Util/exact_riemann/inputs.test1.helm

@@ -19,4 +19,3 @@ problem.t     = 0.0008
 problem.npts  = 128
 
 castro.T_guess = 1.e5
-castro.riemann_shock_maxiter = 20

Util/exact_riemann/inputs.test2.helm

@@ -15,5 +15,3 @@ problem.xjump = 0.5e5
 problem.t     = 0.00008
 
 problem.npts  = 128
-
-castro.riemann_shock_maxiter = 100

Util/exact_riemann/riemann_shock.H

@@ -46,12 +46,12 @@ W_s_shock(const amrex::Real W_s, const amrex::Real pstar,
 
 
 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
@@ -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?

Util/exact_riemann/riemann_star_state.H

@@ -10,7 +10,7 @@
 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) {
@@ -79,8 +79,6 @@ riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex::
 
     // 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:
     //
@@ -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