update gresho problem to take Mach number as input (#2963)

this is easier than setting p0

Exec/hydro_tests/gresho_vortex/_prob_params

@@ -1,11 +1,15 @@
-p0      real     1.0_rt       y
+# initial Mach number
+M0      real     0.1_rt       y
 
+# ambient density
 rho0    real     1.0_rt       y
 
+# reference pressure -- this is set automatically
+p0      real     -1.0_rt
+
+# reference timescale
 t_r     real     1.0_rt       y
 
 x_r     real     0.0_rt
 
 q_r     real     0.0_rt
-
-nsub    integer   4           y

Exec/hydro_tests/gresho_vortex/inputs.2d

@@ -79,16 +79,8 @@ amr.derive_plot_vars=ALL
 # PROBLEM PARAMETERS
 problem.rho0 = 1.0
 
-# p0 = rho0 u_phi**2/(gamma * M**2) - 1/2  -- specify M to get p0
-
-# M ~ 0.001
-problem.p0 = 599999.5
-
-# M ~ 0.01
-# problem.p0 = 5999.5
-
-# M ~ 0.1
-# problem.p0 = 59.5
+# Mach number of the flow
+problem.M0 = 0.1
 
 # EOS
 eos.eos_assume_neutral = 1

Exec/hydro_tests/gresho_vortex/problem_initialize.H

@@ -27,6 +27,17 @@ void problem_initialize ()
     problem::x_r = probhi[0] - problo[0];
     problem::q_r = 0.4_rt * M_PI * problem::x_r / problem::t_r;
 
+    // pressure peaks at r = 1/5 from the center
+    // where p = p0 + 25/2 r^2, so peak pressure is p = p0 + 1/2
+    //
+    // Mach number is |u_phi| / sqrt(gamma p / rho), so we get solve
+    // for p0.  Note that the peak velocity is q,
+    //
+    // p0 = rho q^2 / (gamma M^2) - 1/2
+
+    problem::p0 = problem::rho0 * std::pow(problem::q_r, 2.0) /
+        (eos_rp::eos_gamma * std::pow(problem::M0, 2.0)) - 0.5_rt;
+
 }
 
 #endif

Exec/hydro_tests/gresho_vortex/problem_initialize_state_data.H

@@ -22,71 +22,46 @@ void problem_initialize_state_data (int i, int j, int k,
     const Real* dx = geomdata.CellSize();
 
     Real x = problo[0] + dx[0] * (static_cast<Real>(i) + 0.5_rt);
-    Real xl = problo[0] + dx[0] * (static_cast<Real>(i));
 
     Real y = 0.0;
-    Real yl = 0.0;
-
 #if AMREX_SPACEDIM >= 2
     y = problo[1] + dx[1] * (static_cast<Real>(j) + 0.5_rt);
-    yl = problo[1] + dx[1] * (static_cast<Real>(j));
 #endif
 
     Real z = 0.0;
-    Real zl = 0.0;
-
 #if AMREX_SPACEDIM == 3
     z = problo[2] + dx[2] * (static_cast<Real>(k) + 0.5_rt);
-    zl = problo[2] + dx[2] * (static_cast<Real>(k));
 #endif
 
     Real reint = 0.0;
     Real u_tot = 0.0;
 
-    for (int kk = 0; kk < problem::nsub; kk++) {
-        Real zz = zl + dx[2] * (static_cast<Real>(kk) + 0.5_rt) / problem::nsub;
-
-        for (int jj = 0; jj < problem::nsub; jj++) {
-            Real yy = yl + dx[1] * (static_cast<Real>(jj) + 0.5_rt) / problem::nsub;
-
-            for (int ii = 0; ii < problem::nsub; ii++) {
-                Real xx = xl + dx[0] * (static_cast<Real>(ii) + 0.5_rt) / problem::nsub;
-
-                Real r = std::sqrt((xx - problem::center[0]) * (xx - problem::center[0]) +
-                                   (yy - problem::center[1]) * (yy - problem::center[1]));
-
-                Real u_phi;
-                Real p;
+    Real r = std::sqrt(amrex::Math::powi<2>(x - problem::center[0]) +
+                       amrex::Math::powi<2>(y - problem::center[1]));
 
-                if (r < 0.2_rt) {
-                    u_phi = 5.0_rt * r;
-                    p = problem::p0 + 12.5_rt * r * r;
+    Real u_phi;
+    Real p;
 
-                } else if (r < 0.4_rt) {
-                    u_phi = 2.0_rt - 5.0_rt * r;
-                    p = problem::p0 + 12.5_rt * r * r + 4.0_rt *
-                        (1.0_rt - 5.0_rt * r - std::log(0.2_rt) + std::log(r));
+    if (r < 0.2_rt) {
+        u_phi = 5.0_rt * r;
+        p = problem::p0 + 12.5_rt * r * r;
 
-                } else {
-                    u_phi = 0.0_rt;
-                    p = problem::p0 - 2.0_rt + 4.0_rt * std::log(2.0_rt);
-                }
+    } else if (r < 0.4_rt) {
+        u_phi = 2.0_rt - 5.0_rt * r;
+        p = problem::p0 + 12.5_rt * r * r + 4.0_rt *
+            (1.0_rt - 5.0_rt * r - std::log(0.2_rt) + std::log(r));
 
-                u_tot += u_phi;
-                reint += p/(eos_rp::eos_gamma - 1.0_rt);
-            }
-        }
+    } else {
+        u_phi = 0.0_rt;
+        p = problem::p0 - 2.0_rt + 4.0_rt * std::log(2.0_rt);
     }
 
-    Real u_phi = u_tot / (problem::nsub * problem::nsub * problem::nsub);
-    reint = reint / (problem::nsub * problem::nsub * problem::nsub);
+    reint = p / (eos_rp::eos_gamma - 1.0_rt);
 
     state(i,j,k,URHO) = problem::rho0;
 
     // phi unit vector: \hat{\phi} = -sin(phi) \hat{x} + cos(phi) \hat{y}
     // with cos(phi) = x/r; sin(phi) = y/r
-    Real r = std::sqrt((x - problem::center[0]) * (x - problem::center[0]) +
-                       (y - problem::center[1]) * (y - problem::center[1]));
 
     //  -sin(phi) = y/r
     state(i,j,k,UMX) = -problem::rho0 * problem::q_r * u_phi * ((y - problem::center[1]) / r);