some more commentw

Util/exact_riemann/riemann_rarefaction.H

@@ -131,9 +131,14 @@ rarefaction(const amrex::Real pstar,
     const double S1{0.9};
     const double S2{4.0};
 
-    // Compute Z_s = C for a rarefaction connecting the state to the star
-    // region by integrating the Riemann invariant from p_s to pstar.
-    // This means solving a system of ODEs.  We use RK4.
+    // 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.
+    //
+    // The system is:
+    //
+    //    dtau/dp = -1/C**2
+    //    du/dp = +/- 1/C    (+ for 1-wave, - for 3-wave)
 
     constexpr amrex::Real tol = 1.e-5;
 

Util/exact_riemann/riemann_star_state.H

@@ -17,7 +17,7 @@ riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex::
 
     const amrex::Real tol = 1.e-10_rt;
     const amrex::Real smallp = 1.e-8_rt;
-    const amrex::Real SMALL = 1.e-13_rt;
+    const amrex::Real SMALL = 1.e-8_rt;
 
 
     // get the initial sound speeds
@@ -79,7 +79,9 @@ 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: " << rho_l << " " << u_l << " " << p_l << " " << rho_r << " " << u_r << " " << p_r << std::endl;
+    std::cout << "solving for star state: " << std::endl;
+    std::cout << rho_l << " " << u_l << " " << p_l << std::endl;
+    std::cout << rho_r << " " << u_r << " " << p_r << std::endl;
 
 
     // this procedure follows directly from Colella & Glaz 1985, section 1
@@ -116,6 +118,7 @@ riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex::
 
         amrex::Real Z_l, Z_r;
 
+        std::cout << "iteration: " << iter << " " << pstar << " " << p_l << " " << p_r << std::endl;
         // left wave
 
         if (pstar - p_l > SMALL * p_l) {