From ae6b17c67d3d4001510cc54493fddded8c18ec28 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 26 Sep 2024 15:36:46 -0400 Subject: [PATCH] update gresho problem to take Mach number as input (#2963) this is easier than setting p0 --- Exec/hydro_tests/gresho_vortex/_prob_params | 10 +++- Exec/hydro_tests/gresho_vortex/inputs.2d | 12 +--- .../gresho_vortex/problem_initialize.H | 11 ++++ .../problem_initialize_state_data.H | 55 +++++-------------- 4 files changed, 35 insertions(+), 53 deletions(-) diff --git a/Exec/hydro_tests/gresho_vortex/_prob_params b/Exec/hydro_tests/gresho_vortex/_prob_params index 765bdc4577..15a8bbc082 100644 --- a/Exec/hydro_tests/gresho_vortex/_prob_params +++ b/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 diff --git a/Exec/hydro_tests/gresho_vortex/inputs.2d b/Exec/hydro_tests/gresho_vortex/inputs.2d index 665d0c9c11..8cb28b4bc1 100644 --- a/Exec/hydro_tests/gresho_vortex/inputs.2d +++ b/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 diff --git a/Exec/hydro_tests/gresho_vortex/problem_initialize.H b/Exec/hydro_tests/gresho_vortex/problem_initialize.H index 2dc9f62aae..6f21f8da81 100644 --- a/Exec/hydro_tests/gresho_vortex/problem_initialize.H +++ b/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 diff --git a/Exec/hydro_tests/gresho_vortex/problem_initialize_state_data.H b/Exec/hydro_tests/gresho_vortex/problem_initialize_state_data.H index f87fb240af..1cfa18348b 100644 --- a/Exec/hydro_tests/gresho_vortex/problem_initialize_state_data.H +++ b/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(i) + 0.5_rt); - Real xl = problo[0] + dx[0] * (static_cast(i)); Real y = 0.0; - Real yl = 0.0; - #if AMREX_SPACEDIM >= 2 y = problo[1] + dx[1] * (static_cast(j) + 0.5_rt); - yl = problo[1] + dx[1] * (static_cast(j)); #endif Real z = 0.0; - Real zl = 0.0; - #if AMREX_SPACEDIM == 3 z = problo[2] + dx[2] * (static_cast(k) + 0.5_rt); - zl = problo[2] + dx[2] * (static_cast(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(kk) + 0.5_rt) / problem::nsub; - - for (int jj = 0; jj < problem::nsub; jj++) { - Real yy = yl + dx[1] * (static_cast(jj) + 0.5_rt) / problem::nsub; - - for (int ii = 0; ii < problem::nsub; ii++) { - Real xx = xl + dx[0] * (static_cast(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);