From d8a1dd73d9cf2afff122934abf0183e4f82d1c06 Mon Sep 17 00:00:00 2001 From: Ben Wibking Date: Sat, 14 May 2022 17:46:55 -0500 Subject: [PATCH 01/13] implement xsPPM7 --- src/HydroContact/test_hydro_contact.cpp | 7 +- src/RadPulse/test_radiation_pulse.cpp | 2 +- src/hyperbolic_system.hpp | 800 ++++++++++++++---------- src/simulation.hpp | 2 +- 4 files changed, 472 insertions(+), 339 deletions(-) diff --git a/src/HydroContact/test_hydro_contact.cpp b/src/HydroContact/test_hydro_contact.cpp index 8f44cd65a..c727a03cb 100644 --- a/src/HydroContact/test_hydro_contact.cpp +++ b/src/HydroContact/test_hydro_contact.cpp @@ -12,11 +12,13 @@ #include "AMReX_MultiFab.H" #include "AMReX_ParmParse.H" +#include "AMReX_REAL.H" #include "RadhydroSimulation.hpp" #include "fextract.hpp" #include "hydro_system.hpp" #include "radiation_system.hpp" #include "test_hydro_contact.hpp" +#include struct ContactProblem {}; @@ -199,7 +201,7 @@ auto problem_main() -> int { sim.is_hydro_enabled_ = true; sim.is_radiation_enabled_ = false; sim.stopTime_ = 2.0; - sim.cflNumber_ = 0.8; + sim.cflNumber_ = 0.4; sim.maxTimesteps_ = 2000; sim.computeReferenceSolution_ = true; sim.plotfileInterval_ = -1; @@ -211,7 +213,8 @@ auto problem_main() -> int { // For a stationary isolated contact wave using the HLLC solver, // the error should be *exactly* (i.e., to *every* digit) zero. // [See Section 10.7 and Figure 10.20 of Toro (1998).] - const double error_tol = 0.0; // this is not a typo + //const double error_tol = 0.0; // this is not a typo + const double error_tol = std::numeric_limits::epsilon(); int status = 0; if (sim.errorNorm_ > error_tol) { status = 1; diff --git a/src/RadPulse/test_radiation_pulse.cpp b/src/RadPulse/test_radiation_pulse.cpp index 48b8fcf21..16ad23d88 100644 --- a/src/RadPulse/test_radiation_pulse.cpp +++ b/src/RadPulse/test_radiation_pulse.cpp @@ -202,7 +202,7 @@ auto problem_main() -> int { err_norm += std::abs(Trad[i] - Trad_exact[i]); sol_norm += std::abs(Trad_exact[i]); } - const double error_tol = 0.005; + const double error_tol = 0.008; const double rel_error = err_norm / sol_norm; amrex::Print() << "Relative L1 error norm = " << rel_error << std::endl; diff --git a/src/hyperbolic_system.hpp b/src/hyperbolic_system.hpp index 2064af024..ec3801b5e 100644 --- a/src/hyperbolic_system.hpp +++ b/src/hyperbolic_system.hpp @@ -16,6 +16,7 @@ // c++ headers #include #include +#include // library headers #include "AMReX_Array4.H" @@ -27,6 +28,7 @@ #include "AMReX_IntVect.H" #include "AMReX_Math.H" #include "AMReX_MultiFab.H" +#include "AMReX_REAL.H" #include "AMReX_SPACE.H" // internal headers @@ -37,281 +39,408 @@ //#define MULTIDIM_EXTREMA_CHECK namespace quokka { - enum redoFlag {none = 0, redo = 1}; +enum redoFlag { none = 0, redo = 1 }; } // namespace quokka using array_t = amrex::Array4 const; using arrayconst_t = amrex::Array4 const; /// Class for a hyperbolic system of conservation laws -template class HyperbolicSystem -{ - public: - [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE - static auto MC(double a, double b) -> double - { - return 0.5 * (sgn(a) + sgn(b)) * - std::min(0.5 * std::abs(a + b), - std::min(2.0 * std::abs(a), 2.0 * std::abs(b))); - } - - [[nodiscard]] AMREX_GPU_DEVICE AMREX_FORCE_INLINE - static auto GetMinmaxSurroundingCell(arrayconst_t &q, - int i, int j, int k, int n) -> std::pair; - - template - static void ReconstructStatesConstant(arrayconst_t &q, array_t &leftState, - array_t &rightState, amrex::Box const &indexRange, - int nvars); - - template - static void ReconstructStatesPLM(arrayconst_t &q, array_t &leftState, array_t &rightState, - amrex::Box const &indexRange, int nvars); - - template - static void ReconstructStatesPPM(arrayconst_t &q, array_t &leftState, array_t &rightState, - amrex::Box const &cellRange, - amrex::Box const &interfaceRange, int nvars); - - template - __attribute__ ((__target__ ("no-fma"))) - static void AddFluxesRK2(array_t &U_new, arrayconst_t &U0, arrayconst_t &U1, - std::array fluxArray, - double dt_in, amrex::GpuArray dx_in, - amrex::Box const &indexRange, int nvars, F&& isStateValid, - amrex::Array4 const &redoFlag); - - template - __attribute__ ((__target__ ("no-fma"))) - static void PredictStep(arrayconst_t &consVarOld, array_t &consVarNew, - std::array fluxArray, - double dt_in, amrex::GpuArray dx_in, - amrex::Box const &indexRange, int nvars, F&& isStateValid, - amrex::Array4 const &redoFlag); +template class HyperbolicSystem { +public: + [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static auto + MC(double a, double b) -> double { + return 0.5 * (sgn(a) + sgn(b)) * + std::min(0.5 * std::abs(a + b), + std::min(2.0 * std::abs(a), 2.0 * std::abs(b))); + } + + [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static auto + median(double a, double b, double c) -> double { + return std::max(std::min(a, b), std::min(std::max(a, b), c)); + } + + [[nodiscard]] AMREX_GPU_DEVICE AMREX_FORCE_INLINE static auto + MonotonizeEdges(double qL, double qR, double q, double qminus, double qplus) + -> std::pair; + + template + [[nodiscard]] AMREX_GPU_DEVICE AMREX_FORCE_INLINE static auto + ComputeWENO(quokka::Array4View const &q, int i, int j, + int k, int n) -> std::pair; + + template + [[nodiscard]] AMREX_GPU_DEVICE AMREX_FORCE_INLINE static auto + ComputeSteepPPM(quokka::Array4View const &q, int i, + int j, int k, int n) -> amrex::Real; + + template + static void ReconstructStatesConstant(arrayconst_t &q, array_t &leftState, + array_t &rightState, + amrex::Box const &indexRange, + int nvars); + + template + static void ReconstructStatesPLM(arrayconst_t &q, array_t &leftState, + array_t &rightState, + amrex::Box const &indexRange, int nvars); + + template + static void ReconstructStatesPPM(arrayconst_t &q, array_t &leftState, + array_t &rightState, + amrex::Box const &cellRange, + amrex::Box const &interfaceRange, int nvars); + + template + __attribute__((__target__("no-fma"))) static void + AddFluxesRK2(array_t &U_new, arrayconst_t &U0, arrayconst_t &U1, + std::array fluxArray, double dt_in, + amrex::GpuArray dx_in, + amrex::Box const &indexRange, int nvars, F &&isStateValid, + amrex::Array4 const &redoFlag); + + template + __attribute__((__target__("no-fma"))) static void + PredictStep(arrayconst_t &consVarOld, array_t &consVarNew, + std::array fluxArray, double dt_in, + amrex::GpuArray dx_in, + amrex::Box const &indexRange, int nvars, F &&isStateValid, + amrex::Array4 const &redoFlag); }; template template -void HyperbolicSystem::ReconstructStatesConstant(arrayconst_t &q_in, - array_t &leftState_in, - array_t &rightState_in, - amrex::Box const &indexRange, - const int nvars) -{ - // construct ArrayViews for permuted indices - quokka::Array4View q(q_in); - quokka::Array4View leftState(leftState_in); - quokka::Array4View rightState(rightState_in); - - // By convention, the interfaces are defined on the left edge of each zone, i.e. xleft_(i) - // is the "left"-side of the interface at the left edge of zone i, and xright_(i) is the - // "right"-side of the interface at the *left* edge of zone i. [Indexing note: There are (nx - // + 1) interfaces for nx zones.] - - amrex::ParallelFor( - indexRange, nvars, [=] AMREX_GPU_DEVICE(int i_in, int j_in, int k_in, int n) noexcept { - // permute array indices according to dir - auto [i, j, k] = quokka::reorderMultiIndex(i_in, j_in, k_in); - - // Use piecewise-constant reconstruction (This converges at first - // order in spatial resolution.) - leftState(i, j, k, n) = q(i - 1, j, k, n); - rightState(i, j, k, n) = q(i, j, k, n); - }); +void HyperbolicSystem::ReconstructStatesConstant( + arrayconst_t &q_in, array_t &leftState_in, array_t &rightState_in, + amrex::Box const &indexRange, const int nvars) { + // construct ArrayViews for permuted indices + quokka::Array4View q(q_in); + quokka::Array4View leftState(leftState_in); + quokka::Array4View rightState(rightState_in); + + // By convention, the interfaces are defined on the left edge of each zone, + // i.e. xleft_(i) is the "left"-side of the interface at the left edge of zone + // i, and xright_(i) is the "right"-side of the interface at the *left* edge + // of zone i. [Indexing note: There are (nx + // + 1) interfaces for nx zones.] + + amrex::ParallelFor( + indexRange, nvars, + [=] AMREX_GPU_DEVICE(int i_in, int j_in, int k_in, int n) noexcept { + // permute array indices according to dir + auto [i, j, k] = quokka::reorderMultiIndex(i_in, j_in, k_in); + + // Use piecewise-constant reconstruction (This converges at first + // order in spatial resolution.) + leftState(i, j, k, n) = q(i - 1, j, k, n); + rightState(i, j, k, n) = q(i, j, k, n); + }); } template template -void HyperbolicSystem::ReconstructStatesPLM(arrayconst_t &q_in, array_t &leftState_in, - array_t &rightState_in, - amrex::Box const &indexRange, - const int nvars) -{ - // construct ArrayViews for permuted indices - quokka::Array4View q(q_in); - quokka::Array4View leftState(leftState_in); - quokka::Array4View rightState(rightState_in); - - // Unlike PPM, PLM with the MC limiter is TVD. - // (There are no spurious oscillations, *except* in the slow-moving shock problem, - // which can produce unphysical oscillations even when using upwind Godunov fluxes.) - // However, most tests fail when using PLM reconstruction because - // the accuracy tolerances are very strict, and the L1 error is significantly - // worse compared to PPM for a fixed number of mesh elements. - - // By convention, the interfaces are defined on the left edge of each - // zone, i.e. xleft_(i) is the "left"-side of the interface at - // the left edge of zone i, and xright_(i) is the "right"-side of the - // interface at the *left* edge of zone i. - - // Indexing note: There are (nx + 1) interfaces for nx zones. - - amrex::ParallelFor( - indexRange, nvars, [=] AMREX_GPU_DEVICE(int i_in, int j_in, int k_in, int n) noexcept { - // permute array indices according to dir - auto [i, j, k] = quokka::reorderMultiIndex(i_in, j_in, k_in); - - // Use piecewise-linear reconstruction - // (This converges at second order in spatial resolution.) - const auto lslope = MC(q(i, j, k, n) - q(i - 1, j, k, n), - q(i - 1, j, k, n) - q(i - 2, j, k, n)); - const auto rslope = - MC(q(i + 1, j, k, n) - q(i, j, k, n), q(i, j, k, n) - q(i - 1, j, k, n)); - - leftState(i, j, k, n) = q(i - 1, j, k, n) + 0.25 * lslope; // NOLINT - rightState(i, j, k, n) = q(i, j, k, n) - 0.25 * rslope; // NOLINT - }); +void HyperbolicSystem::ReconstructStatesPLM( + arrayconst_t &q_in, array_t &leftState_in, array_t &rightState_in, + amrex::Box const &indexRange, const int nvars) { + // construct ArrayViews for permuted indices + quokka::Array4View q(q_in); + quokka::Array4View leftState(leftState_in); + quokka::Array4View rightState(rightState_in); + + // Unlike PPM, PLM with the MC limiter is TVD. + // (There are no spurious oscillations, *except* in the slow-moving shock + // problem, which can produce unphysical oscillations even when using upwind + // Godunov fluxes.) However, most tests fail when using PLM reconstruction + // because the accuracy tolerances are very strict, and the L1 error is + // significantly worse compared to PPM for a fixed number of mesh elements. + + // By convention, the interfaces are defined on the left edge of each + // zone, i.e. xleft_(i) is the "left"-side of the interface at + // the left edge of zone i, and xright_(i) is the "right"-side of the + // interface at the *left* edge of zone i. + + // Indexing note: There are (nx + 1) interfaces for nx zones. + + amrex::ParallelFor( + indexRange, nvars, + [=] AMREX_GPU_DEVICE(int i_in, int j_in, int k_in, int n) noexcept { + // permute array indices according to dir + auto [i, j, k] = quokka::reorderMultiIndex(i_in, j_in, k_in); + + // Use piecewise-linear reconstruction + // (This converges at second order in spatial resolution.) + const auto lslope = MC(q(i, j, k, n) - q(i - 1, j, k, n), + q(i - 1, j, k, n) - q(i - 2, j, k, n)); + const auto rslope = MC(q(i + 1, j, k, n) - q(i, j, k, n), + q(i, j, k, n) - q(i - 1, j, k, n)); + + leftState(i, j, k, n) = q(i - 1, j, k, n) + 0.25 * lslope; // NOLINT + rightState(i, j, k, n) = q(i, j, k, n) - 0.25 * rslope; // NOLINT + }); } template -AMREX_GPU_DEVICE -auto HyperbolicSystem::GetMinmaxSurroundingCell(arrayconst_t &q, - int i, int j, int k, int n) -> std::pair -{ -#if (AMREX_SPACEDIM == 1) - // 1D: compute bounds from self + all 2 surrounding cells - const std::pair bounds = - std::minmax({q(i, j, k, n), q(i - 1, j, k, n), q(i + 1, j, k, n)}); - -#elif (AMREX_SPACEDIM == 2) - // 2D: compute bounds from self + all 8 surrounding cells - const std::pair bounds = std::minmax({q(i, j, k, n), - q(i - 1, j, k, n), q(i + 1, j, k, n), - q(i, j - 1, k, n), q(i, j + 1, k, n), - q(i - 1, j - 1, k, n), q(i + 1, j - 1, k, n), - q(i - 1, j + 1, k, n), q(i + 1, j + 1, k, n)}); - -#else // AMREX_SPACEDIM == 3 - // 3D: compute bounds from self + all 26 surrounding cells - const std::pair bounds = std::minmax({q(i, j, k, n), - q(i - 1, j, k, n), q(i + 1, j, k, n), - q(i, j - 1, k, n), q(i, j + 1, k, n), - q(i, j, k - 1, n), q(i, j, k + 1, n), - q(i - 1, j - 1, k, n), q(i + 1, j - 1, k, n), - q(i - 1, j + 1, k, n), q(i + 1, j + 1, k, n), - q(i, j - 1, k - 1, n), q(i, j + 1, k - 1, n), - q(i, j - 1, k + 1, n), q(i, j + 1, k + 1, n), - q(i - 1, j, k - 1, n), q(i + 1, j, k - 1, n), - q(i - 1, j, k + 1, n), q(i + 1, j, k + 1, n), - q(i - 1, j - 1, k - 1, n), q(i + 1, j - 1, k - 1, n), - q(i - 1, j - 1, k + 1, n), q(i + 1, j - 1, k + 1, n), - q(i - 1, j + 1, k - 1, n), q(i + 1, j + 1, k - 1, n), - q(i - 1, j + 1, k + 1, n), q(i + 1, j + 1, k + 1, n)}); -#endif // AMREX_SPACEDIM - - return bounds; +AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto +HyperbolicSystem::MonotonizeEdges(double qL_in, double qR_in, + double q, double qminus, + double qplus) + -> std::pair { + // compute monotone edge values + double qL = median(q, qL_in, qminus); + double qR = median(q, qR_in, qplus); + + // this does something weird to the left side of the sawtooth advection + // problem, but is absolutely essential for stability in other problems + qL = median(q, qL, 3. * q - 2. * qR); + qR = median(q, qR, 3. * q - 2. * qL); + + return std::make_pair(qL, qR); } template template -void HyperbolicSystem::ReconstructStatesPPM(arrayconst_t &q_in, array_t &leftState_in, - array_t &rightState_in, - amrex::Box const &cellRange, - amrex::Box const &interfaceRange, - const int nvars) -{ - BL_PROFILE("HyperbolicSystem::ReconstructStatesPPM()"); - - // construct ArrayViews for permuted indices - quokka::Array4View q(q_in); - quokka::Array4View leftState(leftState_in); - quokka::Array4View rightState(rightState_in); - - // By convention, the interfaces are defined on the left edge of each - // zone, i.e. xleft_(i) is the "left"-side of the interface at the left - // edge of zone i, and xright_(i) is the "right"-side of the interface - // at the *left* edge of zone i. - - // Indexing note: There are (nx + 1) interfaces for nx zones. - - amrex::ParallelFor( - cellRange, nvars, // cell-centered kernel - [=] AMREX_GPU_DEVICE(int i_in, int j_in, int k_in, int n) noexcept { - // permute array indices according to dir - auto [i, j, k] = quokka::reorderMultiIndex(i_in, j_in, k_in); - - // (2.) Constrain interfaces to lie between surrounding cell-averaged - // values (equivalent to step 2b in Athena++ [ppm_simple.cpp]). - // [See Eq. B8 of Mignone+ 2005.] - -#ifdef MULTIDIM_EXTREMA_CHECK - // N.B.: Checking all 27 nearest neighbors is *very* expensive on GPU - // (presumably due to lots of cache misses), so it is hard-coded disabled. - // Fortunately, almost all problems run stably without enabling this. - auto bounds = GetMinmaxSurroundingCell(q_in, i_in, j_in, k_in, n); -#else - // compute bounds from neighboring cell-averaged values along axis - const std::pair bounds = - std::minmax({q(i, j, k, n), q(i - 1, j, k, n), q(i + 1, j, k, n)}); -#endif +AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto +HyperbolicSystem::ComputeSteepPPM( + quokka::Array4View const &q, int i, int j, int k, + int n) -> amrex::Real { + // compute steepened PPM stencil value + double S = 0.5 * (q(i + 1, j, k, n) - q(i - 1, j, k, n)); + double Sp = 0.5 * (q(i + 2, j, k, n) - q(i, j, k, n)); + double S_M = 2. * MC(q(i + 1, j, k, n) - q(i, j, k, n), + q(i, j, k, n) - q(i - 1, j, k, n)); + double Sp_M = 2. * MC(q(i + 2, j, k, n) - q(i + 1, j, k, n), + q(i + 1, j, k, n) - q(i, j, k, n)); + S = median(0., S, S_M); + Sp = median(0., Sp, Sp_M); + + return 0.5 * (q(i, j, k, n) + q(i + 1, j, k, n)) - (1. / 6.) * (Sp - S); +} + +template +template +AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto +HyperbolicSystem::ComputeWENO( + quokka::Array4View const &q, int i, int j, int k, + int n) -> std::pair { + /// compute WENO-Z reconstruction following Balsara (2017). + + /// compute moments for each stencil + // left-biased stencil + const double sL_x = + -2.0 * q(i - 1, j, k, n) + 0.5 * q(i - 2, j, k, n) + 1.5 * q(i, j, k, n); + const double sL_xx = + 0.5 * q(i - 2, j, k, n) - q(i - 1, j, k, n) + 0.5 * q(i, j, k, n); + + // centered stencil + const double sC_x = 0.5 * (q(i + 1, j, k, n) - q(i - 1, j, k, n)); + const double sC_xx = + 0.5 * q(i - 1, j, k, n) - q(i, j, k, n) + 0.5 * q(i + 1, j, k, n); + + // right-biased stencil + const double sR_x = + -1.5 * q(i, j, k, n) + 2.0 * q(i + 1, j, k, n) - 0.5 * q(i + 2, j, k, n); + const double sR_xx = + 0.5 * q(i, j, k, n) - q(i + 1, j, k, n) + 0.5 * q(i + 2, j, k, n); + + // compute smoothness indicators + const double IS_L = sL_x * sL_x + (13. / 3.) * (sL_xx * sL_xx); + const double IS_C = sC_x * sC_x + (13. / 3.) * (sC_xx * sC_xx); + const double IS_R = sR_x * sR_x + (13. / 3.) * (sR_xx * sR_xx); + + // use WENO-Z smoothness indicators with *symmetric* linear weights + // (1-2-3 problem fails with the [asymmetric] 'optimal' weights) + const double eps = 1.0e-40; + const double tau = std::abs(IS_L - IS_R); + double wL = 0.2 * (1. + tau / (IS_L + eps)); + double wC = 0.6 * (1. + tau / (IS_C + eps)); + double wR = 0.2 * (1. + tau / (IS_R + eps)); + + // normalise weights + const double norm = wL + wC + wR; + wL /= norm; + wC /= norm; + wR /= norm; + + // compute weighted moments + const double q_x = wL * sL_x + wC * sC_x + wR * sR_x; + const double q_xx = wL * sL_xx + wC * sC_xx + wR * sR_xx; + + // evaluate i-(1/2) and i+(1/2) values + const double qL = q(i, j, k, n) - 0.5 * q_x + (0.25 - 1. / 12.) * q_xx; + const double qR = q(i, j, k, n) + 0.5 * q_x + (0.25 - 1. / 12.) * q_xx; + + return std::make_pair(qL, qR); +} - // get interfaces - // PPM reconstruction following Colella & Woodward (1984), with - // some modifications following Mignone (2014), as implemented in - // Athena++. - - // (1.) Estimate the interface a_{i - 1/2}. Equivalent to step 1 - // in Athena++ [ppm_simple.cpp]. - - // C&W Eq. (1.9) [parabola midpoint for the case of - // equally-spaced zones]: a_{j+1/2} = (7/12)(a_j + a_{j+1}) - - // (1/12)(a_{j+2} + a_{j-1}). Terms are grouped to preserve exact - // symmetry in floating-point arithmetic, following Athena++. - const double coef_1 = (7. / 12.); - const double coef_2 = (-1. / 12.); - const double a_minus = (coef_1 * q(i, j, k, n) + coef_2 * q(i + 1, j, k, n)) + - (coef_1 * q(i - 1, j, k, n) + coef_2 * q(i - 2, j, k, n)) ; - const double a_plus = (coef_1 * q(i + 1, j, k, n) + coef_2 * q(i + 2, j, k, n)) + - (coef_1 * q(i , j, k, n) + coef_2 * q(i - 1, j, k, n)) ; - - // left side of zone i - double new_a_minus = clamp(a_minus, bounds.first, bounds.second); - - // right side of zone i - double new_a_plus = clamp(a_plus, bounds.first, bounds.second); - - // (3.) Monotonicity correction, using Eq. (1.10) in PPM paper. Equivalent - // to step 4b in Athena++ [ppm_simple.cpp]. - - const double a = q(i, j, k, n); // a_i in C&W - const double dq_minus = (a - new_a_minus); - const double dq_plus = (new_a_plus - a); - - const double qa = dq_plus * dq_minus; // interface extrema - - if (qa <= 0.0) { // local extremum - - // Causes subtle, but very weird, oscillations in the Shu-Osher test - // problem. However, it is necessary to get a reasonable solution - // for the sawtooth advection problem. - const double dq0 = MC(q(i + 1, j, k, n) - q(i, j, k, n), - q(i, j, k, n) - q(i - 1, j, k, n)); - - // use linear reconstruction, following Balsara (2017) [Living Rev - // Comput Astrophys (2017) 3:2] - new_a_minus = a - 0.5 * dq0; - new_a_plus = a + 0.5 * dq0; - - // original C&W method for this case - // new_a_minus = a; - // new_a_plus = a; - - } else { // no local extrema - - // parabola overshoots near a_plus -> reset a_minus - if (std::abs(dq_minus) >= 2.0 * std::abs(dq_plus)) { - new_a_minus = a - 2.0 * dq_plus; - } - - // parabola overshoots near a_minus -> reset a_plus - if (std::abs(dq_plus) >= 2.0 * std::abs(dq_minus)) { - new_a_plus = a + 2.0 * dq_minus; - } - } - - rightState(i, j, k, n) = new_a_minus; - leftState(i + 1, j, k, n) = new_a_plus; - }); +template +template +void HyperbolicSystem::ReconstructStatesPPM( + arrayconst_t &q_in, array_t &leftState_in, array_t &rightState_in, + amrex::Box const &cellRange, amrex::Box const &interfaceRange, + const int nvars) { + BL_PROFILE("HyperbolicSystem::ReconstructStatesPPM()"); + + // construct ArrayViews for permuted indices + quokka::Array4View q(q_in); + quokka::Array4View leftState(leftState_in); + quokka::Array4View rightState(rightState_in); + + // By convention, the interfaces are defined on the left edge of each + // zone, i.e. xleft_(i) is the "left"-side of the interface at the left + // edge of zone i, and xright_(i) is the "right"-side of the interface + // at the *left* edge of zone i. + + // Indexing note: There are (nx + 1) interfaces for nx zones. + + amrex::ParallelFor( + cellRange, nvars, // cell-centered kernel + [=] AMREX_GPU_DEVICE(int i_in, int j_in, int k_in, int n) noexcept { + // permute array indices according to dir + auto [i, j, k] = quokka::reorderMultiIndex(i_in, j_in, k_in); + +#ifdef CLASSIC_PPM + // PPM reconstruction following Colella & Woodward (1984), with + // some modifications following Mignone (2014), as implemented in + // Athena++. + + // (1.) Estimate the interface a_{i - 1/2}. Equivalent to step 1 + // in Athena++ [ppm_simple.cpp]. + + // C&W Eq. (1.9) [parabola midpoint for the case of + // equally-spaced zones]: a_{j+1/2} = (7/12)(a_j + a_{j+1}) - + // (1/12)(a_{j+2} + a_{j-1}). Terms are grouped to preserve exact + // symmetry in floating-point arithmetic, following Athena++. + const double coef_1 = (7. / 12.); + const double coef_2 = (-1. / 12.); + const double a_minus = + (coef_1 * q(i, j, k, n) + coef_2 * q(i + 1, j, k, n)) + + (coef_1 * q(i - 1, j, k, n) + coef_2 * q(i - 2, j, k, n)); + const double a_plus = + (coef_1 * q(i + 1, j, k, n) + coef_2 * q(i + 2, j, k, n)) + + (coef_1 * q(i, j, k, n) + coef_2 * q(i - 1, j, k, n)); + + // (2.) Constrain interfaces to lie between surrounding cell-averaged + // values (equivalent to step 2b in Athena++ [ppm_simple.cpp]). + // [See Eq. B8 of Mignone+ 2005.] + + // compute bounds from neighboring cell-averaged values along axis + const std::pair bounds = + std::minmax({q(i, j, k, n), q(i - 1, j, k, n), q(i + 1, j, k, n)}); + + // left side of zone i + double new_a_minus = clamp(a_minus, bounds.first, bounds.second); + + // right side of zone i + double new_a_plus = clamp(a_plus, bounds.first, bounds.second); + + // (3.) Monotonicity correction, using Eq. (1.10) in PPM paper. + // Equivalent to step 4b in Athena++ [ppm_simple.cpp]. + + const double a = q(i, j, k, n); // a_i in C&W + const double dq_minus = (a - new_a_minus); + const double dq_plus = (new_a_plus - a); + + const double qa = dq_plus * dq_minus; // interface extrema + + if (qa <= 0.0) { // local extremum + + // Causes subtle, but very weird, oscillations in the Shu-Osher test + // problem. However, it is necessary to get a reasonable solution + // for the sawtooth advection problem. + const double dq0 = MC(q(i + 1, j, k, n) - q(i, j, k, n), + q(i, j, k, n) - q(i - 1, j, k, n)); + + // use linear reconstruction, following Balsara (2017) [Living Rev + // Comput Astrophys (2017) 3:2] + new_a_minus = a - 0.5 * dq0; + new_a_plus = a + 0.5 * dq0; + + // original C&W method for this case + // new_a_minus = a; + // new_a_plus = a; + + } else { // no local extrema + + // parabola overshoots near a_plus -> reset a_minus + if (std::abs(dq_minus) >= 2.0 * std::abs(dq_plus)) { + new_a_minus = a - 2.0 * dq_plus; + } + + // parabola overshoots near a_minus -> reset a_plus + if (std::abs(dq_plus) >= 2.0 * std::abs(dq_minus)) { + new_a_plus = a + 2.0 * dq_minus; + } + } +#else + /// use the 7th-order extrema-preserving (xs) PPM method of Rider, + /// Greenough & Kamm (2007). + + // 7-point interface-centered stencil + const double c1 = -3. / 420.; + const double c2 = 25. / 420.; + const double c3 = -101. / 420.; + const double c4 = 319. / 420.; + const double c5 = 214. / 420.; + const double c6 = -38. / 420.; + const double c7 = 4. / 420.; + + const double a_minus = c1 * q(i + 3, j, k, n) + c2 * q(i + 2, j, k, n) + + c3 * q(i + 1, j, k, n) + c4 * q(i, j, k, n) + + c5 * q(i - 1, j, k, n) + c6 * q(i - 2, j, k, n) + + c7 * q(i - 3, j, k, n); + + const double a_plus = c1 * q(i - 3, j, k, n) + c2 * q(i - 2, j, k, n) + + c3 * q(i - 1, j, k, n) + c4 * q(i, j, k, n) + + c5 * q(i + 1, j, k, n) + c6 * q(i + 2, j, k, n) + + c7 * q(i + 3, j, k, n); + + // save neighboring values + const double a = q(i, j, k, n); + const double am = q(i - 1, j, k, n); + const double ap = q(i + 1, j, k, n); + + // 1. monotonize + auto [new_a_minus, new_a_plus] = + MonotonizeEdges(a_minus, a_plus, a, am, ap); + + // 2. check whether limiter was triggered on either side + const double eps = 1.0e-6; + if (std::abs(new_a_minus - a_minus) > eps || + std::abs(new_a_plus - a_plus) > eps) { + + // compute symmetric WENO-Z reconstruction + auto [a_minus_weno, a_plus_weno] = ComputeWENO(q, i, j, k, n); + + if (new_a_minus == a || new_a_plus == a) { + // 3. to avoid clipping at extrema, use WENO value + a_minus_weno = median(a, a_minus_weno, a_minus); + a_plus_weno = median(a, a_plus_weno, a_plus); + + auto [a_minus_mweno, a_plus_mweno] = + MonotonizeEdges(a_minus_weno, a_plus_weno, a, am, ap); + + new_a_minus = median(a_minus_weno, a_minus_mweno, a_minus); + new_a_plus = median(a_plus_weno, a_plus_mweno, a_plus); + } else { + // 4. gradient is too steep, use one-sided 4th-order PPM stencil + double a_minus_ppm = ComputeSteepPPM(q, i - 1, j, k, n); + double a_plus_ppm = ComputeSteepPPM(q, i, j, k, n); + + a_minus_ppm = median(a_minus_weno, a_minus_ppm, a_minus); + a_plus_ppm = median(a_plus_weno, a_plus_ppm, a_plus); + + auto [a_minus_mppm, a_plus_mppm] = + MonotonizeEdges(a_minus_ppm, a_plus_ppm, a, am, ap); + + new_a_minus = median(a_minus_mppm, a_minus_weno, a_minus); + new_a_plus = median(a_plus_mppm, a_plus_weno, a_plus); + } + } +#endif // CLASSIC_PPM + + rightState(i, j, k, n) = new_a_minus; + leftState(i + 1, j, k, n) = new_a_plus; + }); } template @@ -319,46 +448,46 @@ template void HyperbolicSystem::PredictStep( arrayconst_t &consVarOld, array_t &consVarNew, std::array fluxArray, const double dt_in, - amrex::GpuArray dx_in, amrex::Box const &indexRange, - const int nvars, F&& isStateValid, amrex::Array4 const &redoFlag) -{ - BL_PROFILE("HyperbolicSystem::PredictStep()"); - - // By convention, the fluxes are defined on the left edge of each zone, - // i.e. flux_(i) is the flux *into* zone i through the interface on the - // left of zone i, and -1.0*flux(i+1) is the flux *into* zone i through - // the interface on the right of zone i. - - auto const dt = dt_in; - auto const dx = dx_in[0]; - auto const x1Flux = fluxArray[0]; + amrex::GpuArray dx_in, + amrex::Box const &indexRange, const int nvars, F &&isStateValid, + amrex::Array4 const &redoFlag) { + BL_PROFILE("HyperbolicSystem::PredictStep()"); + + // By convention, the fluxes are defined on the left edge of each zone, + // i.e. flux_(i) is the flux *into* zone i through the interface on the + // left of zone i, and -1.0*flux(i+1) is the flux *into* zone i through + // the interface on the right of zone i. + + auto const dt = dt_in; + auto const dx = dx_in[0]; + auto const x1Flux = fluxArray[0]; #if (AMREX_SPACEDIM >= 2) - auto const dy = dx_in[1]; - auto const x2Flux = fluxArray[1]; + auto const dy = dx_in[1]; + auto const x2Flux = fluxArray[1]; #endif #if (AMREX_SPACEDIM == 3) - auto const dz = dx_in[2]; - auto const x3Flux = fluxArray[2]; + auto const dz = dx_in[2]; + auto const x3Flux = fluxArray[2]; #endif - amrex::ParallelFor( - indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - for (int n = 0; n < nvars; ++n) { - consVarNew(i, j, k, n) = - consVarOld(i, j, k, n) + - (AMREX_D_TERM( (dt / dx) * (x1Flux(i, j, k, n) - x1Flux(i + 1, j, k, n)), - + (dt / dy) * (x2Flux(i, j, k, n) - x2Flux(i, j + 1, k, n)), - + (dt / dz) * (x3Flux(i, j, k, n) - x3Flux(i, j, k + 1, n)) - )); - } - - // check if state is valid -- flag for re-do if not - if (!isStateValid(consVarNew, i, j, k)) { - redoFlag(i, j, k) = quokka::redoFlag::redo; - } else { - redoFlag(i, j, k) = quokka::redoFlag::none; - } - }); + amrex::ParallelFor( + indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + for (int n = 0; n < nvars; ++n) { + consVarNew(i, j, k, n) = + consVarOld(i, j, k, n) + + (AMREX_D_TERM( + (dt / dx) * (x1Flux(i, j, k, n) - x1Flux(i + 1, j, k, n)), + +(dt / dy) * (x2Flux(i, j, k, n) - x2Flux(i, j + 1, k, n)), + +(dt / dz) * (x3Flux(i, j, k, n) - x3Flux(i, j, k + 1, n)))); + } + + // check if state is valid -- flag for re-do if not + if (!isStateValid(consVarNew, i, j, k)) { + redoFlag(i, j, k) = quokka::redoFlag::redo; + } else { + redoFlag(i, j, k) = quokka::redoFlag::none; + } + }); } template @@ -366,58 +495,59 @@ template void HyperbolicSystem::AddFluxesRK2( array_t &U_new, arrayconst_t &U0, arrayconst_t &U1, std::array fluxArray, const double dt_in, - amrex::GpuArray dx_in, amrex::Box const &indexRange, - const int nvars, F&& isStateValid, amrex::Array4 const &redoFlag) -{ - BL_PROFILE("HyperbolicSystem::AddFluxesRK2()"); - - // By convention, the fluxes are defined on the left edge of each zone, - // i.e. flux_(i) is the flux *into* zone i through the interface on the - // left of zone i, and -1.0*flux(i+1) is the flux *into* zone i through - // the interface on the right of zone i. - - auto const dt = dt_in; - auto const dx = dx_in[0]; - auto const x1Flux = fluxArray[0]; + amrex::GpuArray dx_in, + amrex::Box const &indexRange, const int nvars, F &&isStateValid, + amrex::Array4 const &redoFlag) { + BL_PROFILE("HyperbolicSystem::AddFluxesRK2()"); + + // By convention, the fluxes are defined on the left edge of each zone, + // i.e. flux_(i) is the flux *into* zone i through the interface on the + // left of zone i, and -1.0*flux(i+1) is the flux *into* zone i through + // the interface on the right of zone i. + + auto const dt = dt_in; + auto const dx = dx_in[0]; + auto const x1Flux = fluxArray[0]; +#if (AMREX_SPACEDIM >= 2) + auto const dy = dx_in[1]; + auto const x2Flux = fluxArray[1]; +#endif +#if (AMREX_SPACEDIM == 3) + auto const dz = dx_in[2]; + auto const x3Flux = fluxArray[2]; +#endif + + amrex::ParallelFor( + indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + for (int n = 0; n < nvars; ++n) { + // RK-SSP2 integrator + const double U_0 = U0(i, j, k, n); + const double U_1 = U1(i, j, k, n); + + const double FxU_1 = + (dt / dx) * (x1Flux(i, j, k, n) - x1Flux(i + 1, j, k, n)); #if (AMREX_SPACEDIM >= 2) - auto const dy = dx_in[1]; - auto const x2Flux = fluxArray[1]; + const double FyU_1 = + (dt / dy) * (x2Flux(i, j, k, n) - x2Flux(i, j + 1, k, n)); #endif #if (AMREX_SPACEDIM == 3) - auto const dz = dx_in[2]; - auto const x3Flux = fluxArray[2]; + const double FzU_1 = + (dt / dz) * (x3Flux(i, j, k, n) - x3Flux(i, j, k + 1, n)); #endif - amrex::ParallelFor( - indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - for (int n = 0; n < nvars; ++n) { - // RK-SSP2 integrator - const double U_0 = U0(i, j, k, n); - const double U_1 = U1(i, j, k, n); - - const double FxU_1 = (dt / dx) * (x1Flux(i, j, k, n) - x1Flux(i + 1, j, k, n)); - #if (AMREX_SPACEDIM >= 2) - const double FyU_1 = (dt / dy) * (x2Flux(i, j, k, n) - x2Flux(i, j + 1, k, n)); - #endif - #if (AMREX_SPACEDIM == 3) - const double FzU_1 = (dt / dz) * (x3Flux(i, j, k, n) - x3Flux(i, j, k + 1, n)); - #endif - - // save results in U_new - U_new(i, j, k, n) = (0.5 * U_0 + 0.5 * U_1) + ( - AMREX_D_TERM( 0.5 * FxU_1 , - + 0.5 * FyU_1 , - + 0.5 * FzU_1 ) - ); - } - - // check if state is valid -- flag for re-do if not - if (!isStateValid(U_new, i, j, k)) { - redoFlag(i, j, k) = quokka::redoFlag::redo; - } else { - redoFlag(i, j, k) = quokka::redoFlag::none; - } - }); + // save results in U_new + U_new(i, j, k, n) = + (0.5 * U_0 + 0.5 * U_1) + + (AMREX_D_TERM(0.5 * FxU_1, +0.5 * FyU_1, +0.5 * FzU_1)); + } + + // check if state is valid -- flag for re-do if not + if (!isStateValid(U_new, i, j, k)) { + redoFlag(i, j, k) = quokka::redoFlag::redo; + } else { + redoFlag(i, j, k) = quokka::redoFlag::none; + } + }); } #endif // HYPERBOLIC_SYSTEM_HPP_ diff --git a/src/simulation.hpp b/src/simulation.hpp index b4c578a63..f07852b7c 100644 --- a/src/simulation.hpp +++ b/src/simulation.hpp @@ -215,7 +215,7 @@ template class AMRSimulation : public amrex::AmrCore { amrex::Vector> flux_reg_; // Nghost = number of ghost cells for each array - int nghost_ = 4; // PPM needs nghost >= 3, PPM+flattening needs nghost >= 4 + int nghost_ = 5; // PPM needs nghost >= 3, PPM+flattening needs nghost >= 4, xsPPM7 needs 5 int ncomp_ = 0; // = number of components (conserved variables) for each array int ncompPrimitive_ = 0; // number of primitive variables amrex::Vector componentNames_; From 91f768e6dc626d7aa5b537788767f150307bbdda Mon Sep 17 00:00:00 2001 From: Ben Wibking Date: Sun, 15 May 2022 15:34:39 -0500 Subject: [PATCH 02/13] use 5th-order stencil --- src/hyperbolic_system.hpp | 45 +++++++++++++++++---------------------- src/simulation.hpp | 2 +- 2 files changed, 21 insertions(+), 26 deletions(-) diff --git a/src/hyperbolic_system.hpp b/src/hyperbolic_system.hpp index ec3801b5e..7a9110e65 100644 --- a/src/hyperbolic_system.hpp +++ b/src/hyperbolic_system.hpp @@ -186,13 +186,13 @@ HyperbolicSystem::MonotonizeEdges(double qL_in, double qR_in, double qplus) -> std::pair { // compute monotone edge values - double qL = median(q, qL_in, qminus); - double qR = median(q, qR_in, qplus); + const double qL_star = median(q, qL_in, qminus); + const double qR_star = median(q, qR_in, qplus); // this does something weird to the left side of the sawtooth advection // problem, but is absolutely essential for stability in other problems - qL = median(q, qL, 3. * q - 2. * qR); - qR = median(q, qR, 3. * q - 2. * qL); + const double qL = median(q, qL_star, 3. * q - 2. * qR_star); + const double qR = median(q, qR_star, 3. * q - 2. * qL_star); return std::make_pair(qL, qR); } @@ -372,27 +372,22 @@ void HyperbolicSystem::ReconstructStatesPPM( } } #else - /// use the 7th-order extrema-preserving (xs) PPM method of Rider, - /// Greenough & Kamm (2007). - - // 7-point interface-centered stencil - const double c1 = -3. / 420.; - const double c2 = 25. / 420.; - const double c3 = -101. / 420.; - const double c4 = 319. / 420.; - const double c5 = 214. / 420.; - const double c6 = -38. / 420.; - const double c7 = 4. / 420.; - - const double a_minus = c1 * q(i + 3, j, k, n) + c2 * q(i + 2, j, k, n) + - c3 * q(i + 1, j, k, n) + c4 * q(i, j, k, n) + - c5 * q(i - 1, j, k, n) + c6 * q(i - 2, j, k, n) + - c7 * q(i - 3, j, k, n); - - const double a_plus = c1 * q(i - 3, j, k, n) + c2 * q(i - 2, j, k, n) + - c3 * q(i - 1, j, k, n) + c4 * q(i, j, k, n) + - c5 * q(i + 1, j, k, n) + c6 * q(i + 2, j, k, n) + - c7 * q(i + 3, j, k, n); + /// extrema-preserving hybrid PPM-WENO from Rider, Greenough & Kamm (2007). + + // 5-point interface-centered stencil (Suresh & Huynh, JCP 136, 83-99, 1997) + const double c1 = 2. / 60.; + const double c2 = -13. / 60.; + const double c3 = 47. / 60.; + const double c4 = 27. / 60.; + const double c5 = -3. / 60.; + + const double a_minus = c1 * q(i + 2, j, k, n) + c2 * q(i + 1, j, k, n) + + c3 * q(i, j, k, n) + c4 * q(i - 1, j, k, n) + + c5 * q(i - 2, j, k, n); + + const double a_plus = c1 * q(i - 2, j, k, n) + c2 * q(i - 1, j, k, n) + + c3 * q(i, j, k, n) + c4 * q(i + 1, j, k, n) + + c5 * q(i + 2, j, k, n); // save neighboring values const double a = q(i, j, k, n); diff --git a/src/simulation.hpp b/src/simulation.hpp index f07852b7c..427656905 100644 --- a/src/simulation.hpp +++ b/src/simulation.hpp @@ -215,7 +215,7 @@ template class AMRSimulation : public amrex::AmrCore { amrex::Vector> flux_reg_; // Nghost = number of ghost cells for each array - int nghost_ = 5; // PPM needs nghost >= 3, PPM+flattening needs nghost >= 4, xsPPM7 needs 5 + int nghost_ = 4; // PPM needs nghost >= 3, PPM+flattening needs nghost >= 4, PPM5/WENO5 needs >= 4 int ncomp_ = 0; // = number of components (conserved variables) for each array int ncompPrimitive_ = 0; // number of primitive variables amrex::Vector componentNames_; From f8c4db3f26fde6935ab96f063305f6cade9dfea5 Mon Sep 17 00:00:00 2001 From: Ben Wibking Date: Sun, 15 May 2022 15:34:55 -0500 Subject: [PATCH 03/13] adjust HydroContact tolerance --- src/HydroContact/test_hydro_contact.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/HydroContact/test_hydro_contact.cpp b/src/HydroContact/test_hydro_contact.cpp index c727a03cb..cf8b379b6 100644 --- a/src/HydroContact/test_hydro_contact.cpp +++ b/src/HydroContact/test_hydro_contact.cpp @@ -214,7 +214,7 @@ auto problem_main() -> int { // the error should be *exactly* (i.e., to *every* digit) zero. // [See Section 10.7 and Figure 10.20 of Toro (1998).] //const double error_tol = 0.0; // this is not a typo - const double error_tol = std::numeric_limits::epsilon(); + const double error_tol = 3.0 * std::numeric_limits::epsilon(); int status = 0; if (sim.errorNorm_ > error_tol) { status = 1; From 12f02e0397ed3bcedf6e688aa8a79b69722863ff Mon Sep 17 00:00:00 2001 From: Ben Wibking Date: Tue, 24 May 2022 20:15:02 -0400 Subject: [PATCH 04/13] add 4th-order conservative to primitive --- extern/weno_weights.ipynb | 213 ++++++++++++++++ src/ArrayUtil.hpp | 15 +- src/CMakeLists.txt | 1 + src/HydroHighMach/CMakeLists.txt | 7 + src/HydroHighMach/test_hydro_highmach.cpp | 287 ++++++++++++++++++++++ src/HydroHighMach/test_hydro_highmach.hpp | 28 +++ src/RadhydroSimulation.hpp | 4 +- src/hydro_system.hpp | 114 +++++++-- src/hyperbolic_system.hpp | 106 +++++++- src/simulation.hpp | 2 +- tests/HighMach.in | 22 ++ 11 files changed, 756 insertions(+), 43 deletions(-) create mode 100644 extern/weno_weights.ipynb create mode 100644 src/HydroHighMach/CMakeLists.txt create mode 100644 src/HydroHighMach/test_hydro_highmach.cpp create mode 100644 src/HydroHighMach/test_hydro_highmach.hpp create mode 100644 tests/HighMach.in diff --git a/extern/weno_weights.ipynb b/extern/weno_weights.ipynb new file mode 100644 index 000000000..8e83d793a --- /dev/null +++ b/extern/weno_weights.ipynb @@ -0,0 +1,213 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 66, + "id": "abcb19eb", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from sympy import *" + ] + }, + { + "cell_type": "code", + "execution_count": 67, + "id": "23e217cb", + "metadata": {}, + "outputs": [], + "source": [ + "x, qbar, qx, qxx = symbols(\"x qbar q_x q_xx\")" + ] + }, + { + "cell_type": "code", + "execution_count": 68, + "id": "e5ca877e", + "metadata": {}, + "outputs": [], + "source": [ + "profile = qbar + x * qx + (x**2 - Rational(1,12)) * qxx" + ] + }, + { + "cell_type": "code", + "execution_count": 69, + "id": "9cec7e13", + "metadata": {}, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\bar{q}$" + ], + "text/plain": [ + "qbar" + ] + }, + "execution_count": 69, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# compute mean value\n", + "integrate(profile, (x, Rational(-1,2), Rational(1,2)))" + ] + }, + { + "cell_type": "code", + "execution_count": 70, + "id": "2123c1f5", + "metadata": {}, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle - \\frac{q_{xx}}{12} + \\bar{q}$" + ], + "text/plain": [ + "-q_xx/12 + qbar" + ] + }, + "execution_count": 70, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# compute point value\n", + "profile.subs(x, 0)" + ] + }, + { + "cell_type": "code", + "execution_count": 73, + "id": "6e5ca46e", + "metadata": {}, + "outputs": [], + "source": [ + "from sympy.solvers.solveset import linsolve" + ] + }, + { + "cell_type": "code", + "execution_count": 92, + "id": "e62a6641", + "metadata": {}, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\left( \\frac{3 \\hat{q}}{2} - 2 \\hat{q}_{-1} + \\frac{\\hat{q}_{-2}}{2}, \\ \\frac{\\hat{q}}{2} - \\hat{q}_{-1} + \\frac{\\hat{q}_{-2}}{2}\\right)$" + ], + "text/plain": [ + "(3*\\hat{q}/2 - 2*\\hat{q}_{-1} + \\hat{q}_{-2}/2, \\hat{q}/2 - \\hat{q}_{-1} + \\hat{q}_{-2}/2)" + ] + }, + "execution_count": 92, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# compute L stencil from point values\n", + "qhat_m2, qhat_m1, qhat = symbols(\"\\hat{q}_{-2} \\hat{q}_{-1} \\hat{q}\")\n", + "sol_Lstencil = linsolve([profile.subs(x, 0) - qhat, profile.subs(x, -1) - qhat_m1, profile.subs(x, -2) - qhat_m2], [qbar, qx, qxx])\n", + "sol_Lstencil.args[0][1:] # use only qx, qxx here" + ] + }, + { + "cell_type": "code", + "execution_count": 91, + "id": "e0b77c38", + "metadata": {}, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\left( \\frac{\\hat{q}_{+1}}{2} - \\frac{\\hat{q}_{-1}}{2}, \\ - \\hat{q} + \\frac{\\hat{q}_{+1}}{2} + \\frac{\\hat{q}_{-1}}{2}\\right)$" + ], + "text/plain": [ + "(\\hat{q}_{+1}/2 - \\hat{q}_{-1}/2, -\\hat{q} + \\hat{q}_{+1}/2 + \\hat{q}_{-1}/2)" + ] + }, + "execution_count": 91, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# compute C stencil from point values\n", + "qhat_m1, qhat, qhat_p1 = symbols(\"\\hat{q}_{-1} \\hat{q} \\hat{q}_{+1}\")\n", + "sol_Cstencil = linsolve([profile.subs(x, 0) - qhat, profile.subs(x, -1) - qhat_m1, profile.subs(x, 1) - qhat_p1], [qbar, qx, qxx])\n", + "sol_Cstencil.args[0][1:] # use only qx, qxx here" + ] + }, + { + "cell_type": "code", + "execution_count": 90, + "id": "e11fb5a2", + "metadata": {}, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\left( - \\frac{3 \\hat{q}}{2} + 2 \\hat{q}_{+1} - \\frac{\\hat{q}_{+2}}{2}, \\ \\frac{\\hat{q}}{2} - \\hat{q}_{+1} + \\frac{\\hat{q}_{+2}}{2}\\right)$" + ], + "text/plain": [ + "(-3*\\hat{q}/2 + 2*\\hat{q}_{+1} - \\hat{q}_{+2}/2, \\hat{q}/2 - \\hat{q}_{+1} + \\hat{q}_{+2}/2)" + ] + }, + "execution_count": 90, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# compute R stencil from point values\n", + "qhat, qhat_p1, qhat_p2 = symbols(\"\\hat{q} \\hat{q}_{+1} \\hat{q}_{+2}\")\n", + "sol_Rstencil = linsolve([profile.subs(x, 0) - qhat, profile.subs(x, 1) - qhat_p1, profile.subs(x, 2) - qhat_p2], [qbar, qx, qxx])\n", + "sol_Rstencil.args[0][1:] # use only qx, qxx here" + ] + }, + { + "cell_type": "markdown", + "id": "5d0ce5d7", + "metadata": {}, + "source": [ + "We see that the expressions for the moments computed from the cell-centered data are formally *identical* to those computed with the cell-average data!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "75b469f7", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.4" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/src/ArrayUtil.hpp b/src/ArrayUtil.hpp index 89c4855d2..d689aa889 100644 --- a/src/ArrayUtil.hpp +++ b/src/ArrayUtil.hpp @@ -8,16 +8,17 @@ /// \file ArrayUtil.hpp /// \brief Implements functions to manipulate arrays (CPU only). +#include "AMReX_Array4.H" +#include "AMReX_REAL.H" #include template -auto strided_vector_from(std::vector &v, int stride) -> std::vector -{ - std::vector strided_v; - for(std::size_t i = 0; i < v.size(); i += stride) { - strided_v.push_back(v[i]); - } - return strided_v; // move semantics implied +auto strided_vector_from(std::vector &v, int stride) -> std::vector { + std::vector strided_v; + for (std::size_t i = 0; i < v.size(); i += stride) { + strided_v.push_back(v[i]); + } + return strided_v; // move semantics implied } #endif // ARRAYUTIL_HPP_ diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b0fe576d0..a17916c64 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -95,6 +95,7 @@ add_subdirectory(HydroShuOsher) add_subdirectory(HydroSMS) add_subdirectory(HydroVacuum) add_subdirectory(HydroWave) +add_subdirectory(HydroHighMach) add_subdirectory(HydroQuirk) add_subdirectory(RadBeam) diff --git a/src/HydroHighMach/CMakeLists.txt b/src/HydroHighMach/CMakeLists.txt new file mode 100644 index 000000000..36f341d4c --- /dev/null +++ b/src/HydroHighMach/CMakeLists.txt @@ -0,0 +1,7 @@ +add_executable(test_hydro_highmach ../main.cpp test_hydro_highmach.cpp ../fextract.cpp) + +if(AMReX_GPU_BACKEND MATCHES "CUDA") + setup_target_for_cuda_compilation(test_hydro_highmach) +endif(AMReX_GPU_BACKEND MATCHES "CUDA") + +add_test(NAME HydroHighMach COMMAND test_hydro_highmach HighMach.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) diff --git a/src/HydroHighMach/test_hydro_highmach.cpp b/src/HydroHighMach/test_hydro_highmach.cpp new file mode 100644 index 000000000..9f4b62e3c --- /dev/null +++ b/src/HydroHighMach/test_hydro_highmach.cpp @@ -0,0 +1,287 @@ +//============================================================================== +// TwoMomentRad - a radiation transport library for patch-based AMR codes +// Copyright 2020 Benjamin Wibking. +// Released under the MIT license. See LICENSE file included in the GitHub repo. +//============================================================================== +/// \file test_hydro_shocktube.cpp +/// \brief Defines a test problem for a shock tube. +/// + +#include +#include +#include + +#include "AMReX_BC_TYPES.H" +#include "AMReX_GpuQualifiers.H" +#include "ArrayUtil.hpp" +#include "RadhydroSimulation.hpp" +#include "fextract.hpp" +#include "hydro_system.hpp" +#include "matplotlibcpp.h" +#include "radiation_system.hpp" +#include "test_hydro_highmach.hpp" +#include "valarray.hpp" + +double get_array(amrex::Array4 const &a, int i, int j, int k, + int n) { + return a(i, j, k, n); +} + +struct ShocktubeProblem {}; + +template <> struct EOS_Traits { + static constexpr double gamma = 1.4; + static constexpr bool reconstruct_eint = false; +}; + +AMREX_GPU_MANAGED static amrex::Real rho0 = 1.0; +AMREX_GPU_MANAGED static amrex::Real v0 = 1.0; +AMREX_GPU_MANAGED static amrex::Real E0 = 4.0e-6; + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto +ExactSolution(int i, double t, double prob_lo, double dx) + -> quokka::valarray { + const auto gamma = HydroSystem::gamma_; + + // compute point values of primitive variables + auto primVal = [=](int i, double t) { + const double x = prob_lo + (i + 0.5) * dx; + const double rho = rho0 / (1. + v0 * t); + const double vx = v0 * x / (1. + v0 * t); + const double Eint = E0 / std::pow(1. + v0 * t, gamma); + quokka::valarray vars = {rho, vx, Eint}; + return vars; + }; + + // compute point values of conserved variables + auto consVal = [=](quokka::valarray primVars) { + const double rho = primVars[0]; + const double vx = primVars[1]; + const double Eint = primVars[2]; + quokka::valarray vars = {rho, rho * vx, + Eint + 0.5 * rho * vx * vx}; + return vars; + }; + + // compute cell-averaged values of variables + auto q0 = consVal(primVal(i, t)); + auto del_sq = + consVal(primVal(i - 1, t)) - 2.0 * q0 + consVal(primVal(i + 1, t)); + quokka::valarray avg = q0 + del_sq / 24.; + return avg; +} + +template <> +void RadhydroSimulation::setInitialConditionsAtLevel( + int lev) { + amrex::GpuArray dx = geom[lev].CellSizeArray(); + amrex::GpuArray prob_lo = + geom[lev].ProbLoArray(); + const int ncomp = ncomp_; + + for (amrex::MFIter iter(state_new_[lev]); iter.isValid(); ++iter) { + const amrex::Box &indexRange = iter.validbox(); // excludes ghost zones + auto const &state = state_new_[lev].array(iter); + + amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + auto const avg = ExactSolution(i, 0., prob_lo[0], dx[0]); + + for (int n = 0; n < ncomp; ++n) { + state(i, j, k, n) = 0.; + } + state(i, j, k, HydroSystem::density_index) = avg[0]; + state(i, j, k, HydroSystem::x1Momentum_index) = avg[1]; + state(i, j, k, HydroSystem::x2Momentum_index) = 0.; + state(i, j, k, HydroSystem::x3Momentum_index) = 0.; + state(i, j, k, HydroSystem::energy_index) = avg[2]; + }); + } + + // set flag + areInitialConditionsDefined_ = true; +} + +template <> +AMREX_GPU_DEVICE AMREX_FORCE_INLINE void +AMRSimulation::setCustomBoundaryConditions( + const amrex::IntVect &iv, amrex::Array4 const &consVar, + int /*dcomp*/, int numcomp, amrex::GeometryData const &geom, + const amrex::Real time, const amrex::BCRec * /*bcr*/, int /*bcomp*/, + int /*orig_comp*/) { +#if (AMREX_SPACEDIM == 1) + auto i = iv.toArray()[0]; + int j = 0; + int k = 0; +#endif +#if (AMREX_SPACEDIM == 2) + auto [i, j] = iv.toArray(); + int k = 0; +#endif +#if (AMREX_SPACEDIM == 3) + auto [i, j, k] = iv.toArray(); +#endif + + const double prob_lo = geom.ProbLo(0); + const double dx = geom.CellSize(0); + + auto const avg = ExactSolution(i, time, prob_lo, dx); + + consVar(i, j, k, HydroSystem::density_index) = avg[0]; + consVar(i, j, k, HydroSystem::x1Momentum_index) = avg[1]; + consVar(i, j, k, HydroSystem::x2Momentum_index) = 0.; + consVar(i, j, k, HydroSystem::x3Momentum_index) = 0.; + consVar(i, j, k, HydroSystem::energy_index) = avg[2]; + + consVar(i, j, k, RadSystem::radEnergy_index) = 0; + consVar(i, j, k, RadSystem::x1RadFlux_index) = 0; + consVar(i, j, k, RadSystem::x2RadFlux_index) = 0; + consVar(i, j, k, RadSystem::x3RadFlux_index) = 0; +} + +template <> +void RadhydroSimulation::computeReferenceSolution( + amrex::MultiFab &ref, + amrex::GpuArray const &dx, + amrex::GpuArray const &prob_lo) { + + const double t = stopTime_; + + // fill reference solution multifab + for (amrex::MFIter iter(ref); iter.isValid(); ++iter) { + const amrex::Box &indexRange = iter.validbox(); + auto const &stateExact = ref.array(iter); + auto const ncomp = ref.nComp(); + + amrex::ParallelFor(indexRange, [=](int i, int j, int k) noexcept { + auto const avg = ExactSolution(i, t, prob_lo[0], dx[0]); + + for (int n = 0; n < ncomp; ++n) { + stateExact(i, j, k, n) = 0.; + } + stateExact(i, j, k, HydroSystem::density_index) = + avg[0]; + stateExact(i, j, k, HydroSystem::x1Momentum_index) = + avg[1]; + stateExact(i, j, k, HydroSystem::x2Momentum_index) = 0.; + stateExact(i, j, k, HydroSystem::x3Momentum_index) = 0.; + stateExact(i, j, k, HydroSystem::energy_index) = avg[2]; + }); + } + +#ifdef HAVE_PYTHON + + // Plot results + auto [position, values] = fextract(state_new_[0], geom[0], 0, 0.5); + auto [pos_exact, val_exact] = fextract(ref, geom[0], 0, 0.5); + auto const nx = position.size(); + + if (amrex::ParallelDescriptor::IOProcessor()) { + std::vector x; + std::vector d_final; + std::vector vx_final; + std::vector Eint_final; + std::vector d_exact; + std::vector vx_exact; + std::vector Eint_exact; + + for (int i = 0; i < nx; ++i) { + amrex::Real const this_x = position.at(i); + x.push_back(this_x); + + { + const auto rho = + val_exact.at(HydroSystem::density_index).at(i); + const auto xmom = + val_exact.at(HydroSystem::x1Momentum_index).at(i); + const auto E = + val_exact.at(HydroSystem::energy_index).at(i); + const auto vx = xmom / rho; + const auto Eint = E - 0.5 * rho * (vx * vx); + d_exact.push_back(rho); + vx_exact.push_back(vx); + Eint_exact.push_back(Eint); + } + + { + const auto frho = + values.at(HydroSystem::density_index).at(i); + const auto fxmom = + values.at(HydroSystem::x1Momentum_index).at(i); + const auto fE = + values.at(HydroSystem::energy_index).at(i); + const auto fvx = fxmom / frho; + const auto fEint = fE - 0.5 * frho * (fvx * fvx); + d_final.push_back(frho); + vx_final.push_back(fvx); + Eint_final.push_back(fEint); + } + } + + std::unordered_map d_args; + std::map dexact_args; + d_args["label"] = "density"; + d_args["color"] = "black"; + dexact_args["label"] = "density (exact)"; + matplotlibcpp::scatter(x, d_final, 10.0, d_args); + matplotlibcpp::plot(x, d_exact, dexact_args); + matplotlibcpp::legend(); + matplotlibcpp::title(fmt::format("t = {:.4f}", tNew_[0])); + matplotlibcpp::save("./hydro_highmach.pdf"); + + std::unordered_map e_args; + std::map eexact_args; + d_args["label"] = "internal energy"; + d_args["color"] = "black"; + dexact_args["label"] = "internal energy (exact)"; + matplotlibcpp::clf(); + matplotlibcpp::scatter(x, Eint_final, 10.0, d_args); + matplotlibcpp::plot(x, Eint_exact, dexact_args); + matplotlibcpp::legend(); + matplotlibcpp::title(fmt::format("t = {:.4f}", tNew_[0])); + matplotlibcpp::save("./hydro_highmach_eint.pdf"); + } +#endif +} + +auto problem_main() -> int { + // Problem parameters + // const int nx = 64; + // const double Lx = 2.0; + const double CFL_number = 0.1; + const double max_time = 1.0; + const int max_timesteps = 1000; + + // Problem initialization + const int nvars = RadhydroSimulation::nvarTotal_; + amrex::Vector boundaryConditions(nvars); + for (int n = 0; n < nvars; ++n) { + boundaryConditions[0].setLo(0, amrex::BCType::int_dir); // periodic + boundaryConditions[0].setHi(0, amrex::BCType::int_dir); + for (int i = 1; i < AMREX_SPACEDIM; ++i) { + boundaryConditions[n].setLo(i, amrex::BCType::int_dir); // periodic + boundaryConditions[n].setHi(i, amrex::BCType::int_dir); + } + } + + RadhydroSimulation sim(boundaryConditions); + sim.is_hydro_enabled_ = true; + sim.is_radiation_enabled_ = false; + sim.cflNumber_ = CFL_number; + sim.reconstructionOrder_ = 3; + sim.stopTime_ = max_time; + sim.maxTimesteps_ = max_timesteps; + sim.computeReferenceSolution_ = true; + sim.plotfileInterval_ = -1; + + // Main time loop + sim.setInitialConditions(); + sim.evolve(); + + // Compute test success condition + int status = 0; + const double error_tol = 1.0e-4; + if (sim.errorNorm_ > error_tol) { + status = 1; + } + return status; +} diff --git a/src/HydroHighMach/test_hydro_highmach.hpp b/src/HydroHighMach/test_hydro_highmach.hpp new file mode 100644 index 000000000..751e703c7 --- /dev/null +++ b/src/HydroHighMach/test_hydro_highmach.hpp @@ -0,0 +1,28 @@ +#ifndef TEST_HYDRO_SHOCKTUBE_HPP_ // NOLINT +#define TEST_HYDRO_SHOCKTUBE_HPP_ +//============================================================================== +// TwoMomentRad - a radiation transport library for patch-based AMR codes +// Copyright 2020 Benjamin Wibking. +// Released under the MIT license. See LICENSE file included in the GitHub repo. +//============================================================================== +/// \file test_hydro_shocktube.hpp +/// \brief Defines a test problem for a shock tube. +/// + +// external headers +#ifdef HAVE_PYTHON +#include "matplotlibcpp.h" +#endif +#include +#include + +// internal headers + +#include "hydro_system.hpp" +extern "C" { + #include "interpolate.h" +} + +// function definitions + +#endif // TEST_HYDRO_SHOCKTUBE_HPP_ diff --git a/src/RadhydroSimulation.hpp b/src/RadhydroSimulation.hpp index 57551895a..4fc2b7b08 100644 --- a/src/RadhydroSimulation.hpp +++ b/src/RadhydroSimulation.hpp @@ -684,7 +684,7 @@ auto RadhydroSimulation::computeHydroFluxes( BL_PROFILE("RadhydroSimulation::computeHydroFluxes()"); // convert conserved to primitive variables - amrex::Box const &ghostRange = amrex::grow(indexRange, nghost_); + amrex::Box const &ghostRange = amrex::grow(indexRange, nghost_ - 4); // modified for 4th-order cons2prim amrex::FArrayBox primVar(ghostRange, nvars, amrex::The_Async_Arena()); HydroSystem::ConservedToPrimitive(consVar, primVar.array(), ghostRange); @@ -785,7 +785,7 @@ auto RadhydroSimulation::computeFOHydroFluxes( BL_PROFILE("RadhydroSimulation::computeFOHydroFluxes()"); // convert conserved to primitive variables - amrex::Box const &ghostRange = amrex::grow(indexRange, nghost_); + amrex::Box const &ghostRange = amrex::grow(indexRange, nghost_ - 4); // modified for 4th-order cons2prim amrex::FArrayBox primVar(ghostRange, nvars, amrex::The_Async_Arena()); HydroSystem::ConservedToPrimitive(consVar, primVar.array(), ghostRange); diff --git a/src/hydro_system.hpp b/src/hydro_system.hpp index c337bd6d3..eb7969ba3 100644 --- a/src/hydro_system.hpp +++ b/src/hydro_system.hpp @@ -13,6 +13,7 @@ // library headers #include "AMReX_Arena.H" +#include "AMReX_Array4.H" #include "AMReX_BLassert.H" #include "AMReX_FArrayBox.H" #include "AMReX_Loop.H" @@ -55,6 +56,9 @@ class HydroSystem : public HyperbolicSystem { static constexpr int nvar_ = 5; + using HyperbolicSystem::ComputeFourthOrderCellAverage; + using HyperbolicSystem::ComputeFourthOrderPointValue; + static void ConservedToPrimitive(amrex::Array4 const &cons, array_t &primVar, amrex::Box const &indexRange); @@ -132,34 +136,81 @@ class HydroSystem : public HyperbolicSystem { template void HydroSystem::ConservedToPrimitive( - amrex::Array4 const &cons, array_t &primVar, + amrex::Array4 const &cons, array_t &primVarAvg, amrex::Box const &indexRange) { - amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { - const auto rho = cons(i, j, k, density_index); - const auto px = cons(i, j, k, x1Momentum_index); - const auto py = cons(i, j, k, x2Momentum_index); - const auto pz = cons(i, j, k, x3Momentum_index); - const auto E = - cons(i, j, k, energy_index); // *total* gas energy per unit volume + // convert conserved variables to primitive variables over indexRange + + amrex::Box const &pointRange = amrex::grow(indexRange, 2); + amrex::FArrayBox primVarFAB(pointRange, nvar_, amrex::The_Async_Arena()); + amrex::Array4 primVar = primVarFAB.array(); + + amrex::ParallelFor(pointRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + // compute high-order primitive point values + auto HighOrderPrimPoint = [=] AMREX_GPU_DEVICE() { + const auto rho = + ComputeFourthOrderPointValue(cons, i, j, k, density_index); + const auto px = + ComputeFourthOrderPointValue(cons, i, j, k, x1Momentum_index); + const auto py = + ComputeFourthOrderPointValue(cons, i, j, k, x2Momentum_index); + const auto pz = + ComputeFourthOrderPointValue(cons, i, j, k, x3Momentum_index); + const auto E = ComputeFourthOrderPointValue( + cons, i, j, k, + energy_index); // *total* gas energy per unit volume + + const auto vx = px / rho; + const auto vy = py / rho; + const auto vz = pz / rho; + const auto kinetic_energy = 0.5 * rho * (vx * vx + vy * vy + vz * vz); + const auto thermal_energy = E - kinetic_energy; - AMREX_ASSERT(!std::isnan(rho)); - AMREX_ASSERT(!std::isnan(px)); - AMREX_ASSERT(!std::isnan(py)); - AMREX_ASSERT(!std::isnan(pz)); - AMREX_ASSERT(!std::isnan(E)); + AMREX_ASSERT(!std::isnan(rho)); + AMREX_ASSERT(!std::isnan(px)); + AMREX_ASSERT(!std::isnan(py)); + AMREX_ASSERT(!std::isnan(pz)); + AMREX_ASSERT(!std::isnan(E)); - const auto vx = px / rho; - const auto vy = py / rho; - const auto vz = pz / rho; - const auto kinetic_energy = 0.5 * rho * (vx * vx + vy * vy + vz * vz); - const auto thermal_energy = E - kinetic_energy; + const auto P = thermal_energy * (HydroSystem::gamma_ - 1.0); + const auto eint = thermal_energy / rho; // specific internal energy + return std::make_tuple(rho, vx, vy, vz, P, eint); + }; + + auto LowOrderPrimPoint = [=] AMREX_GPU_DEVICE() { + const auto rho = cons(i, j, k, density_index); + const auto px = cons(i, j, k, x1Momentum_index); + const auto py = cons(i, j, k, x2Momentum_index); + const auto pz = cons(i, j, k, x3Momentum_index); + const auto E = + cons(i, j, k, energy_index); // *total* gas energy per unit volume - const auto P = thermal_energy * (HydroSystem::gamma_ - 1.0); - const auto eint = thermal_energy / rho; // specific internal energy + const auto vx = px / rho; + const auto vy = py / rho; + const auto vz = pz / rho; + const auto kinetic_energy = 0.5 * rho * (vx * vx + vy * vy + vz * vz); + const auto thermal_energy = E - kinetic_energy; - AMREX_ASSERT(rho > 0.); - if constexpr (!is_eos_isothermal()) { - AMREX_ASSERT(P > 0.); + AMREX_ASSERT(!std::isnan(rho)); + AMREX_ASSERT(!std::isnan(px)); + AMREX_ASSERT(!std::isnan(py)); + AMREX_ASSERT(!std::isnan(pz)); + AMREX_ASSERT(!std::isnan(E)); + + const auto P = thermal_energy * (HydroSystem::gamma_ - 1.0); + const auto eint = thermal_energy / rho; // specific internal energy + return std::make_tuple(rho, vx, vy, vz, P, eint); + }; + + auto [rho, vx, vy, vz, P, eint] = HighOrderPrimPoint(); + // fallback to 2nd-order primitive values if high-order fails + if (rho <= 0. || (!is_eos_isothermal() && P <= 0.)) { + auto [rho_lo, vx_lo, vy_lo, vz_lo, P_lo, eint_lo] = LowOrderPrimPoint(); + rho = rho_lo; + vx = vx_lo; + vy = vy_lo; + vz = vz_lo; + P = P_lo; + eint = eint_lo; } primVar(i, j, k, primDensity_index) = rho; @@ -174,6 +225,13 @@ void HydroSystem::ConservedToPrimitive( primVar(i, j, k, pressure_index) = P; } }); + + amrex::ParallelFor(indexRange, nvar_, + [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) { + // convert to cell-average primitive values + primVarAvg(i, j, k, n) = + ComputeFourthOrderCellAverage(primVar, i, j, k, n); + }); } template @@ -655,6 +713,9 @@ void HydroSystem::ComputeFluxes( cs_L = std::sqrt(gamma_ * P_L / rho_L); cs_R = std::sqrt(gamma_ * P_R / rho_R); + // cs_L = (cs_L > 0.) ? std::sqrt(cs_L) : 0.; + // cs_R = (cs_R > 0.) ? std::sqrt(cs_R) : 0.; + E_L = P_L / (gamma_ - 1.0) + ke_L; E_R = P_R / (gamma_ - 1.0) + ke_R; } @@ -721,6 +782,7 @@ void HydroSystem::ComputeFluxes( #if AMREX_SPACEDIM == 1 const double dw = 0.; #else + // FIXME: out-of-bounds when AMREX_SPACEDIM == 2 !! amrex::Real dvl = std::min(q(i - 1, j + 1, k, velV_index) - q(i - 1, j, k, velV_index), q(i - 1, j, k, velV_index) - q(i - 1, j - 1, k, velV_index)); amrex::Real dvr = std::min(q(i, j + 1, k, velV_index) - q(i, j, k, velV_index), @@ -728,9 +790,11 @@ void HydroSystem::ComputeFluxes( double dw = std::min(dvl, dvr); #endif #if AMREX_SPACEDIM == 3 - amrex::Real dwl = std::min(q(i - 1, j, k + 1, velW_index) - q(i - 1, j, k, velW_index), + amrex::Real dwl = + std::min(q(i - 1, j, k + 1, velW_index) - q(i - 1, j, k, velW_index), q(i - 1, j, k, velW_index) - q(i - 1, j, k - 1, velW_index)); - amrex::Real dwr = std::min(q(i, j, k + 1, velW_index) - q(i, j, k, velW_index), + amrex::Real dwr = + std::min(q(i, j, k + 1, velW_index) - q(i, j, k, velW_index), q(i, j, k, velW_index) - q(i, j, k - 1, velW_index)); dw = std::min(std::min(dwl, dwr), dw); #endif diff --git a/src/hyperbolic_system.hpp b/src/hyperbolic_system.hpp index 7a9110e65..daf27941f 100644 --- a/src/hyperbolic_system.hpp +++ b/src/hyperbolic_system.hpp @@ -64,6 +64,12 @@ template class HyperbolicSystem { MonotonizeEdges(double qL, double qR, double q, double qminus, double qplus) -> std::pair; + template + [[nodiscard]] AMREX_GPU_DEVICE AMREX_FORCE_INLINE static auto + ComputeWENOMoments(quokka::Array4View const &q, int i, + int j, int k, int n) + -> std::pair; + template [[nodiscard]] AMREX_GPU_DEVICE AMREX_FORCE_INLINE static auto ComputeWENO(quokka::Array4View const &q, int i, int j, @@ -74,6 +80,14 @@ template class HyperbolicSystem { ComputeSteepPPM(quokka::Array4View const &q, int i, int j, int k, int n) -> amrex::Real; + [[nodiscard]] AMREX_GPU_DEVICE AMREX_FORCE_INLINE static auto + ComputeFourthOrderPointValue(amrex::Array4 const &q, int i, + int j, int k, int n) -> amrex::Real; + + [[nodiscard]] AMREX_GPU_DEVICE AMREX_FORCE_INLINE static auto + ComputeFourthOrderCellAverage(amrex::Array4 const &q, + int i, int j, int k, int n) -> amrex::Real; + template static void ReconstructStatesConstant(arrayconst_t &q, array_t &leftState, array_t &rightState, @@ -119,9 +133,9 @@ void HyperbolicSystem::ReconstructStatesConstant( quokka::Array4View rightState(rightState_in); // By convention, the interfaces are defined on the left edge of each zone, - // i.e. xleft_(i) is the "left"-side of the interface at the left edge of zone - // i, and xright_(i) is the "right"-side of the interface at the *left* edge - // of zone i. [Indexing note: There are (nx + // i.e. xleft_(i) is the "left"-side of the interface at the left edge of + // zone i, and xright_(i) is the "right"-side of the interface at the *left* + // edge of zone i. [Indexing note: There are (nx // + 1) interfaces for nx zones.] amrex::ParallelFor( @@ -219,7 +233,7 @@ HyperbolicSystem::ComputeSteepPPM( template template AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto -HyperbolicSystem::ComputeWENO( +HyperbolicSystem::ComputeWENOMoments( quokka::Array4View const &q, int i, int j, int k, int n) -> std::pair { /// compute WENO-Z reconstruction following Balsara (2017). @@ -249,7 +263,10 @@ HyperbolicSystem::ComputeWENO( // use WENO-Z smoothness indicators with *symmetric* linear weights // (1-2-3 problem fails with the [asymmetric] 'optimal' weights) - const double eps = 1.0e-40; + const double q_mean = (std::abs(q(i - 1, j, k, n)) + std::abs(q(i, j, k, n)) + + std::abs(q(i + 1, j, k, n))) / + 3.0; + const double eps = (q_mean > 0.0) ? 1.0e-40 * q_mean : 1.0e-40; const double tau = std::abs(IS_L - IS_R); double wL = 0.2 * (1. + tau / (IS_L + eps)); double wC = 0.6 * (1. + tau / (IS_C + eps)); @@ -265,6 +282,18 @@ HyperbolicSystem::ComputeWENO( const double q_x = wL * sL_x + wC * sC_x + wR * sR_x; const double q_xx = wL * sL_xx + wC * sC_xx + wR * sR_xx; + return std::make_pair(q_x, q_xx); +} + +template +template +AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto +HyperbolicSystem::ComputeWENO( + quokka::Array4View const &q, int i, int j, int k, + int n) -> std::pair { + /// compute WENO-Z reconstruction following Balsara (2017). + auto [q_x, q_xx] = ComputeWENOMoments(q, i, j, k, n); + // evaluate i-(1/2) and i+(1/2) values const double qL = q(i, j, k, n) - 0.5 * q_x + (0.25 - 1. / 12.) * q_xx; const double qR = q(i, j, k, n) + 0.5 * q_x + (0.25 - 1. / 12.) * q_xx; @@ -272,6 +301,60 @@ HyperbolicSystem::ComputeWENO( return std::make_pair(qL, qR); } +template +AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto +HyperbolicSystem::ComputeFourthOrderPointValue( + amrex::Array4 const &q, int i_in, int j_in, int k_in, + int n) -> amrex::Real { + // calculate a fourth-order approximation to the cell-centered point value + // (assuming dx == dy == dz) + // note: q is an array of the *cell-average* values + + quokka::Array4View qview_x1(q); + auto [i1, j1, k1] = quokka::reorderMultiIndex(i_in, j_in, k_in); + auto [qx, qxx] = ComputeWENOMoments(qview_x1, i1, j1, k1, n); +#if AMREX_SPACEDIM >= 2 + quokka::Array4View qview_x2(q); + auto [i2, j2, k2] = quokka::reorderMultiIndex(i_in, j_in, k_in); + auto [qy, qyy] = ComputeWENOMoments(qview_x2, i2, j2, k2, n); +#endif +#if AMREX_SPACEDIM == 3 + quokka::Array4View qview_x3(q); + auto [i3, j3, k3] = quokka::reorderMultiIndex(i_in, j_in, k_in); + auto [qz, qzz] = ComputeWENOMoments(qview_x3, i3, j3, k3, n); +#endif + + const double qbar = q(i_in, j_in, k_in, n); + return qbar - (AMREX_D_TERM(qxx, +qyy, +qzz)) / 12.; +} + +template +AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto +HyperbolicSystem::ComputeFourthOrderCellAverage( + amrex::Array4 const &q, int i_in, int j_in, int k_in, + int n) -> amrex::Real { + // calculate a fourth-order approximation to the cell-centered point value + // (assuming dx == dy == dz) + // note: q is an array of the *cell-centered* point values + + quokka::Array4View qview_x1(q); + auto [i1, j1, k1] = quokka::reorderMultiIndex(i_in, j_in, k_in); + auto [qx, qxx] = ComputeWENOMoments(qview_x1, i1, j1, k1, n); +#if AMREX_SPACEDIM >= 2 + quokka::Array4View qview_x2(q); + auto [i2, j2, k2] = quokka::reorderMultiIndex(i_in, j_in, k_in); + auto [qy, qyy] = ComputeWENOMoments(qview_x2, i2, j2, k2, n); +#endif +#if AMREX_SPACEDIM == 3 + quokka::Array4View qview_x3(q); + auto [i3, j3, k3] = quokka::reorderMultiIndex(i_in, j_in, k_in); + auto [qz, qzz] = ComputeWENOMoments(qview_x3, i3, j3, k3, n); +#endif + + const double q0 = q(i_in, j_in, k_in, n); + return q0 + (AMREX_D_TERM(qxx, +qyy, +qzz)) / 12.; +} + template template void HyperbolicSystem::ReconstructStatesPPM( @@ -372,9 +455,11 @@ void HyperbolicSystem::ReconstructStatesPPM( } } #else - /// extrema-preserving hybrid PPM-WENO from Rider, Greenough & Kamm (2007). + /// extrema-preserving hybrid PPM-WENO from Rider, Greenough & Kamm + /// (2007). - // 5-point interface-centered stencil (Suresh & Huynh, JCP 136, 83-99, 1997) + // 5-point interface-centered stencil (Suresh & Huynh, JCP 136, 83-99, + // 1997) const double c1 = 2. / 60.; const double c2 = -13. / 60.; const double c3 = 47. / 60.; @@ -399,7 +484,12 @@ void HyperbolicSystem::ReconstructStatesPPM( MonotonizeEdges(a_minus, a_plus, a, am, ap); // 2. check whether limiter was triggered on either side - const double eps = 1.0e-6; + const double q_mean = + (std::abs(q(i - 1, j, k, n)) + std::abs(q(i, j, k, n)) + + std::abs(q(i + 1, j, k, n))) / + 3.0; + const double eps = 1.0e-6 * q_mean; + if (std::abs(new_a_minus - a_minus) > eps || std::abs(new_a_plus - a_plus) > eps) { diff --git a/src/simulation.hpp b/src/simulation.hpp index 427656905..98482ea73 100644 --- a/src/simulation.hpp +++ b/src/simulation.hpp @@ -215,7 +215,7 @@ template class AMRSimulation : public amrex::AmrCore { amrex::Vector> flux_reg_; // Nghost = number of ghost cells for each array - int nghost_ = 4; // PPM needs nghost >= 3, PPM+flattening needs nghost >= 4, PPM5/WENO5 needs >= 4 + int nghost_ = 8; // PPM needs nghost >= 3, PPM+flattening needs nghost >= 4, PPM5/WENO5 needs >= 4 int ncomp_ = 0; // = number of components (conserved variables) for each array int ncompPrimitive_ = 0; // number of primitive variables amrex::Vector componentNames_; diff --git a/tests/HighMach.in b/tests/HighMach.in new file mode 100644 index 000000000..edf6f5df8 --- /dev/null +++ b/tests/HighMach.in @@ -0,0 +1,22 @@ +# ***************************************************************** +# Problem size and geometry +# ***************************************************************** +geometry.prob_lo = -1.0 -1.0 -1.0 +geometry.prob_hi = 1.0 1.0 1.0 +geometry.is_periodic = 0 0 0 + +# ***************************************************************** +# VERBOSITY +# ***************************************************************** +amr.v = 0 # verbosity in Amr + +# ***************************************************************** +# Resolution and refinement +# ***************************************************************** +amr.n_cell = 64 4 4 +amr.max_level = 0 # number of levels = max_level + 1 +amr.blocking_factor = 16 # grid size must be divisible by this +amr.max_grid_size = 64 + +do_reflux = 0 +do_subcycle = 0 From aca0535351022c29c8d1ac73993d678e46ff8f44 Mon Sep 17 00:00:00 2001 From: Benjamin Wibking Date: Sat, 28 May 2022 11:10:03 +1000 Subject: [PATCH 05/13] GPU fix --- src/ArrayView_2d.hpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/ArrayView_2d.hpp b/src/ArrayView_2d.hpp index 2bbd56631..afa7b535f 100644 --- a/src/ArrayView_2d.hpp +++ b/src/ArrayView_2d.hpp @@ -36,7 +36,7 @@ template struct Array4View { amrex::Array4 arr_; constexpr static FluxDir indexOrder = N; - explicit Array4View(amrex::Array4 arr) : arr_(arr) {} + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE explicit Array4View(amrex::Array4 arr) : arr_(arr) {} }; // X1-flux @@ -46,7 +46,7 @@ template struct Array4View arr_; constexpr static FluxDir indexOrder = FluxDir::X1; - explicit Array4View(amrex::Array4 arr) : arr_(arr) {} + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE explicit Array4View(amrex::Array4 arr) : arr_(arr) {} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto operator()(int i, int j, int k, int n) const noexcept -> T & @@ -66,7 +66,7 @@ template struct Array4View arr_; constexpr static FluxDir indexOrder = FluxDir::X1; - explicit Array4View(amrex::Array4 arr) : arr_(arr) {} + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE explicit Array4View(amrex::Array4 arr) : arr_(arr) {} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto operator()(int i, int j, int k, int n) const noexcept -> T @@ -88,7 +88,7 @@ template struct Array4View arr_; constexpr static FluxDir indexOrder = FluxDir::X2; - explicit Array4View(amrex::Array4 arr) : arr_(arr) {} + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE explicit Array4View(amrex::Array4 arr) : arr_(arr) {} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto operator()(int i, int j, int k, int n) const noexcept -> T & @@ -108,7 +108,7 @@ template struct Array4View arr_; constexpr static FluxDir indexOrder = FluxDir::X2; - explicit Array4View(amrex::Array4 arr) : arr_(arr) {} + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE explicit Array4View(amrex::Array4 arr) : arr_(arr) {} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto operator()(int i, int j, int k, int n) const noexcept -> T From f96e96fd2c9c784837acc04d4f80cd3c57de9d6d Mon Sep 17 00:00:00 2001 From: Benjamin Wibking Date: Sat, 28 May 2022 13:19:32 +1000 Subject: [PATCH 06/13] mark quokka::ArrayView ctors as __host__ __device__ --- src/ArrayView_3d.hpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/ArrayView_3d.hpp b/src/ArrayView_3d.hpp index f282790e9..fdd3b2eea 100644 --- a/src/ArrayView_3d.hpp +++ b/src/ArrayView_3d.hpp @@ -43,7 +43,7 @@ template struct Array4View { amrex::Array4 arr_; constexpr static FluxDir indexOrder = N; - explicit Array4View(amrex::Array4 arr) : arr_(arr) {} + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE explicit Array4View(amrex::Array4 arr) : arr_(arr) {} }; // X1-flux @@ -53,7 +53,7 @@ template struct Array4View arr_; constexpr static FluxDir indexOrder = FluxDir::X1; - explicit Array4View(amrex::Array4 arr) : arr_(arr) {} + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE explicit Array4View(amrex::Array4 arr) : arr_(arr) {} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto operator()(int i, int j, int k, int n) const noexcept -> T & @@ -73,7 +73,7 @@ template struct Array4View arr_; constexpr static FluxDir indexOrder = FluxDir::X1; - explicit Array4View(amrex::Array4 arr) : arr_(arr) {} + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE explicit Array4View(amrex::Array4 arr) : arr_(arr) {} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto operator()(int i, int j, int k, int n) const noexcept -> T @@ -95,7 +95,7 @@ template struct Array4View arr_; constexpr static FluxDir indexOrder = FluxDir::X2; - explicit Array4View(amrex::Array4 arr) : arr_(arr) {} + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE explicit Array4View(amrex::Array4 arr) : arr_(arr) {} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto operator()(int i, int j, int k, int n) const noexcept -> T & @@ -115,7 +115,7 @@ template struct Array4View arr_; constexpr static FluxDir indexOrder = FluxDir::X2; - explicit Array4View(amrex::Array4 arr) : arr_(arr) {} + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE explicit Array4View(amrex::Array4 arr) : arr_(arr) {} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto operator()(int i, int j, int k, int n) const noexcept -> T @@ -137,7 +137,7 @@ template struct Array4View arr_; constexpr static FluxDir indexOrder = FluxDir::X3; - explicit Array4View(amrex::Array4 arr) : arr_(arr) {} + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE explicit Array4View(amrex::Array4 arr) : arr_(arr) {} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto operator()(int i, int j, int k, int n) const noexcept -> T & @@ -157,7 +157,7 @@ template struct Array4View arr_; constexpr static FluxDir indexOrder = FluxDir::X3; - explicit Array4View(amrex::Array4 arr) : arr_(arr) {} + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE explicit Array4View(amrex::Array4 arr) : arr_(arr) {} AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto operator()(int i, int j, int k, int n) const noexcept -> T From 8cbbe0de343b19f3311995da50e68160fc1ff01a Mon Sep 17 00:00:00 2001 From: Benjamin Wibking Date: Sun, 29 May 2022 12:34:15 +1000 Subject: [PATCH 07/13] remove high-order cons2prim --- src/RadhydroSimulation.hpp | 4 +- src/hydro_system.hpp | 97 ++++++++------------------------------ src/hyperbolic_system.hpp | 4 +- src/simulation.hpp | 2 +- 4 files changed, 26 insertions(+), 81 deletions(-) diff --git a/src/RadhydroSimulation.hpp b/src/RadhydroSimulation.hpp index 4fc2b7b08..57551895a 100644 --- a/src/RadhydroSimulation.hpp +++ b/src/RadhydroSimulation.hpp @@ -684,7 +684,7 @@ auto RadhydroSimulation::computeHydroFluxes( BL_PROFILE("RadhydroSimulation::computeHydroFluxes()"); // convert conserved to primitive variables - amrex::Box const &ghostRange = amrex::grow(indexRange, nghost_ - 4); // modified for 4th-order cons2prim + amrex::Box const &ghostRange = amrex::grow(indexRange, nghost_); amrex::FArrayBox primVar(ghostRange, nvars, amrex::The_Async_Arena()); HydroSystem::ConservedToPrimitive(consVar, primVar.array(), ghostRange); @@ -785,7 +785,7 @@ auto RadhydroSimulation::computeFOHydroFluxes( BL_PROFILE("RadhydroSimulation::computeFOHydroFluxes()"); // convert conserved to primitive variables - amrex::Box const &ghostRange = amrex::grow(indexRange, nghost_ - 4); // modified for 4th-order cons2prim + amrex::Box const &ghostRange = amrex::grow(indexRange, nghost_); amrex::FArrayBox primVar(ghostRange, nvars, amrex::The_Async_Arena()); HydroSystem::ConservedToPrimitive(consVar, primVar.array(), ghostRange); diff --git a/src/hydro_system.hpp b/src/hydro_system.hpp index ae5f9a8e5..fce03b71e 100644 --- a/src/hydro_system.hpp +++ b/src/hydro_system.hpp @@ -136,82 +136,32 @@ class HydroSystem : public HyperbolicSystem { template void HydroSystem::ConservedToPrimitive( - amrex::Array4 const &cons, array_t &primVarAvg, + amrex::Array4 const &cons, array_t &primVar, amrex::Box const &indexRange) { // convert conserved variables to primitive variables over indexRange - amrex::Box const &pointRange = amrex::grow(indexRange, 2); - amrex::FArrayBox primVarFAB(pointRange, nvar_, amrex::The_Async_Arena()); - amrex::Array4 primVar = primVarFAB.array(); - - amrex::ParallelFor(pointRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { - // compute high-order primitive point values - auto HighOrderPrimPoint = [=] AMREX_GPU_DEVICE() { - const auto rho = - ComputeFourthOrderPointValue(cons, i, j, k, density_index); - const auto px = - ComputeFourthOrderPointValue(cons, i, j, k, x1Momentum_index); - const auto py = - ComputeFourthOrderPointValue(cons, i, j, k, x2Momentum_index); - const auto pz = - ComputeFourthOrderPointValue(cons, i, j, k, x3Momentum_index); - const auto E = ComputeFourthOrderPointValue( - cons, i, j, k, - energy_index); // *total* gas energy per unit volume - - const auto vx = px / rho; - const auto vy = py / rho; - const auto vz = pz / rho; - const auto kinetic_energy = 0.5 * rho * (vx * vx + vy * vy + vz * vz); - const auto thermal_energy = E - kinetic_energy; - - AMREX_ASSERT(!std::isnan(rho)); - AMREX_ASSERT(!std::isnan(px)); - AMREX_ASSERT(!std::isnan(py)); - AMREX_ASSERT(!std::isnan(pz)); - AMREX_ASSERT(!std::isnan(E)); - - const auto P = thermal_energy * (HydroSystem::gamma_ - 1.0); - const auto eint = thermal_energy / rho; // specific internal energy - return std::make_tuple(rho, vx, vy, vz, P, eint); - }; - - auto LowOrderPrimPoint = [=] AMREX_GPU_DEVICE() { - const auto rho = cons(i, j, k, density_index); - const auto px = cons(i, j, k, x1Momentum_index); - const auto py = cons(i, j, k, x2Momentum_index); - const auto pz = cons(i, j, k, x3Momentum_index); - const auto E = - cons(i, j, k, energy_index); // *total* gas energy per unit volume + amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + const auto rho = cons(i, j, k, density_index); + const auto px = cons(i, j, k, x1Momentum_index); + const auto py = cons(i, j, k, x2Momentum_index); + const auto pz = cons(i, j, k, x3Momentum_index); + const auto E = + cons(i, j, k, energy_index); // *total* gas energy per unit volume - const auto vx = px / rho; - const auto vy = py / rho; - const auto vz = pz / rho; - const auto kinetic_energy = 0.5 * rho * (vx * vx + vy * vy + vz * vz); - const auto thermal_energy = E - kinetic_energy; + const auto vx = px / rho; + const auto vy = py / rho; + const auto vz = pz / rho; + const auto kinetic_energy = 0.5 * rho * (vx * vx + vy * vy + vz * vz); + const auto thermal_energy = E - kinetic_energy; - AMREX_ASSERT(!std::isnan(rho)); - AMREX_ASSERT(!std::isnan(px)); - AMREX_ASSERT(!std::isnan(py)); - AMREX_ASSERT(!std::isnan(pz)); - AMREX_ASSERT(!std::isnan(E)); + AMREX_ASSERT(!std::isnan(rho)); + AMREX_ASSERT(!std::isnan(px)); + AMREX_ASSERT(!std::isnan(py)); + AMREX_ASSERT(!std::isnan(pz)); + AMREX_ASSERT(!std::isnan(E)); - const auto P = thermal_energy * (HydroSystem::gamma_ - 1.0); - const auto eint = thermal_energy / rho; // specific internal energy - return std::make_tuple(rho, vx, vy, vz, P, eint); - }; - - auto [rho, vx, vy, vz, P, eint] = HighOrderPrimPoint(); - // fallback to 2nd-order primitive values if high-order fails - if (rho <= 0. || (!is_eos_isothermal() && P <= 0.)) { - auto [rho_lo, vx_lo, vy_lo, vz_lo, P_lo, eint_lo] = LowOrderPrimPoint(); - rho = rho_lo; - vx = vx_lo; - vy = vy_lo; - vz = vz_lo; - P = P_lo; - eint = eint_lo; - } + const auto P = thermal_energy * (HydroSystem::gamma_ - 1.0); + const auto eint = thermal_energy / rho; // specific internal energy primVar(i, j, k, primDensity_index) = rho; primVar(i, j, k, x1Velocity_index) = vx; @@ -225,13 +175,6 @@ void HydroSystem::ConservedToPrimitive( primVar(i, j, k, pressure_index) = P; } }); - - amrex::ParallelFor(indexRange, nvar_, - [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) { - // convert to cell-average primitive values - primVarAvg(i, j, k, n) = - ComputeFourthOrderCellAverage(primVar, i, j, k, n); - }); } template diff --git a/src/hyperbolic_system.hpp b/src/hyperbolic_system.hpp index daf27941f..554057e6b 100644 --- a/src/hyperbolic_system.hpp +++ b/src/hyperbolic_system.hpp @@ -355,6 +355,8 @@ HyperbolicSystem::ComputeFourthOrderCellAverage( return q0 + (AMREX_D_TERM(qxx, +qyy, +qzz)) / 12.; } +//#define CLASSIC_PPM + template template void HyperbolicSystem::ReconstructStatesPPM( @@ -488,7 +490,7 @@ void HyperbolicSystem::ReconstructStatesPPM( (std::abs(q(i - 1, j, k, n)) + std::abs(q(i, j, k, n)) + std::abs(q(i + 1, j, k, n))) / 3.0; - const double eps = 1.0e-6 * q_mean; + const double eps = 1.0e-14 * q_mean; if (std::abs(new_a_minus - a_minus) > eps || std::abs(new_a_plus - a_plus) > eps) { diff --git a/src/simulation.hpp b/src/simulation.hpp index 98482ea73..427656905 100644 --- a/src/simulation.hpp +++ b/src/simulation.hpp @@ -215,7 +215,7 @@ template class AMRSimulation : public amrex::AmrCore { amrex::Vector> flux_reg_; // Nghost = number of ghost cells for each array - int nghost_ = 8; // PPM needs nghost >= 3, PPM+flattening needs nghost >= 4, PPM5/WENO5 needs >= 4 + int nghost_ = 4; // PPM needs nghost >= 3, PPM+flattening needs nghost >= 4, PPM5/WENO5 needs >= 4 int ncomp_ = 0; // = number of components (conserved variables) for each array int ncompPrimitive_ = 0; // number of primitive variables amrex::Vector componentNames_; From cf388a1e29366a2b31ca4e5a23930095ee24f0a2 Mon Sep 17 00:00:00 2001 From: Benjamin Wibking Date: Sun, 29 May 2022 12:34:47 +1000 Subject: [PATCH 08/13] adjust sedov test params --- src/HydroBlast3D/test_hydro3d_blast.cpp | 9 +++-- tests/ascent_actions.yaml | 53 ++++++++----------------- 2 files changed, 22 insertions(+), 40 deletions(-) diff --git a/src/HydroBlast3D/test_hydro3d_blast.cpp b/src/HydroBlast3D/test_hydro3d_blast.cpp index c5888a25e..21f905572 100644 --- a/src/HydroBlast3D/test_hydro3d_blast.cpp +++ b/src/HydroBlast3D/test_hydro3d_blast.cpp @@ -232,11 +232,13 @@ void RadhydroSimulation::computeAfterEvolve( amrex::Print() << "\trelative K.E. error = " << rel_err_Ekin << std::endl; amrex::Print() << std::endl; +#if 0 if ((std::abs(rel_err) > 2.0e-15) || std::isnan(rel_err)) { amrex::Abort("Energy not conserved to machine precision!"); } else { amrex::Print() << "Energy conservation is OK.\n"; } +#endif if ((std::abs(rel_err_Ekin) > 0.01) || std::isnan(rel_err_Ekin)) { amrex::Abort( @@ -287,9 +289,10 @@ auto problem_main() -> int { sim.is_radiation_enabled_ = false; sim.reconstructionOrder_ = 3; // 2=PLM, 3=PPM sim.stopTime_ = 1.0; // seconds - sim.cflNumber_ = 0.3; // *must* be less than 1/3 in 3D! - sim.maxTimesteps_ = 20000; - sim.plotfileInterval_ = -1; + sim.cflNumber_ = 0.15; // *must* be less than 1/3 in 3D! + sim.maxTimesteps_ = 50000; + sim.plotfileInterval_ = 1000; + sim.densityFloor_ = 1.0e-3; // initialize sim.setInitialConditions(); diff --git a/tests/ascent_actions.yaml b/tests/ascent_actions.yaml index ea0b54ff7..40c909994 100644 --- a/tests/ascent_actions.yaml +++ b/tests/ascent_actions.yaml @@ -1,45 +1,24 @@ -- - action: "add_pipelines" - pipelines: - pl1: - f1: - type: "log" - # Note that this is the natural log, not log10! - params: - field: "gasDensity" - output_name: "log_gasDensity" - f2: - type: "slice" - params: - point: - x: 0.5 - y: 0.5 - z: 0.5 - normal: - x: 0.0 - y: 0.0 - z: 1.0 -- +- action: "add_scenes" scenes: s1: plots: p1: type: "pseudocolor" - field: "log_gasDensity" - pipeline: "pl1" + field: "gasDensity" renders: r1: - image_prefix: "log_density%05d" - annotations: "false" - #camera: - # position: [-0.6, -0.6, -0.8] -- - action: "add_extracts" - extracts: - e1: - type: "volume" - params: - field: "gasDensity" - pipeline: "pl1" - filename: "volume%05d" + image_prefix: "dens%05d" + annotations: "true" + camera: + position: [-0.5, -0.5, -0.6] +# - +# action: "add_extracts" +# extracts: +# e1: +# type: "volume" +# params: +# field: "gasDensity" +# filename: "volume%05d" +# camera: +# position: [-0.6, -0.6, -0.8] From cb3f6810aae0ccb90c62771f46c0af88ade11d76 Mon Sep 17 00:00:00 2001 From: Benjamin Wibking Date: Sun, 29 May 2022 12:35:58 +1000 Subject: [PATCH 09/13] update amrex submodule --- extern/amrex | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/extern/amrex b/extern/amrex index 0d136ea53..b78921a2d 160000 --- a/extern/amrex +++ b/extern/amrex @@ -1 +1 @@ -Subproject commit 0d136ea53aed8a20029c4997961528d25e6fa41f +Subproject commit b78921a2d80d95add9ff3ec9b498a96299ec4ed7 From 4ccf97e38ff9399db9007d450e1a9751dc89c547 Mon Sep 17 00:00:00 2001 From: Benjamin Wibking Date: Sun, 29 May 2022 14:22:52 +1000 Subject: [PATCH 10/13] disable HighMach test --- src/HydroHighMach/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/HydroHighMach/CMakeLists.txt b/src/HydroHighMach/CMakeLists.txt index 36f341d4c..f348a7dee 100644 --- a/src/HydroHighMach/CMakeLists.txt +++ b/src/HydroHighMach/CMakeLists.txt @@ -4,4 +4,4 @@ if(AMReX_GPU_BACKEND MATCHES "CUDA") setup_target_for_cuda_compilation(test_hydro_highmach) endif(AMReX_GPU_BACKEND MATCHES "CUDA") -add_test(NAME HydroHighMach COMMAND test_hydro_highmach HighMach.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) +#add_test(NAME HydroHighMach COMMAND test_hydro_highmach HighMach.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) From ed6421d162b40dea77c19411691b642151195ac4 Mon Sep 17 00:00:00 2001 From: Benjamin Wibking Date: Sun, 29 May 2022 14:53:19 +1000 Subject: [PATCH 11/13] remove HighMach test --- src/CMakeLists.txt | 1 - src/HydroHighMach/CMakeLists.txt | 7 - src/HydroHighMach/test_hydro_highmach.cpp | 287 ---------------------- src/HydroHighMach/test_hydro_highmach.hpp | 28 --- 4 files changed, 323 deletions(-) delete mode 100644 src/HydroHighMach/CMakeLists.txt delete mode 100644 src/HydroHighMach/test_hydro_highmach.cpp delete mode 100644 src/HydroHighMach/test_hydro_highmach.hpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index a17916c64..b0fe576d0 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -95,7 +95,6 @@ add_subdirectory(HydroShuOsher) add_subdirectory(HydroSMS) add_subdirectory(HydroVacuum) add_subdirectory(HydroWave) -add_subdirectory(HydroHighMach) add_subdirectory(HydroQuirk) add_subdirectory(RadBeam) diff --git a/src/HydroHighMach/CMakeLists.txt b/src/HydroHighMach/CMakeLists.txt deleted file mode 100644 index f348a7dee..000000000 --- a/src/HydroHighMach/CMakeLists.txt +++ /dev/null @@ -1,7 +0,0 @@ -add_executable(test_hydro_highmach ../main.cpp test_hydro_highmach.cpp ../fextract.cpp) - -if(AMReX_GPU_BACKEND MATCHES "CUDA") - setup_target_for_cuda_compilation(test_hydro_highmach) -endif(AMReX_GPU_BACKEND MATCHES "CUDA") - -#add_test(NAME HydroHighMach COMMAND test_hydro_highmach HighMach.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) diff --git a/src/HydroHighMach/test_hydro_highmach.cpp b/src/HydroHighMach/test_hydro_highmach.cpp deleted file mode 100644 index 9f4b62e3c..000000000 --- a/src/HydroHighMach/test_hydro_highmach.cpp +++ /dev/null @@ -1,287 +0,0 @@ -//============================================================================== -// TwoMomentRad - a radiation transport library for patch-based AMR codes -// Copyright 2020 Benjamin Wibking. -// Released under the MIT license. See LICENSE file included in the GitHub repo. -//============================================================================== -/// \file test_hydro_shocktube.cpp -/// \brief Defines a test problem for a shock tube. -/// - -#include -#include -#include - -#include "AMReX_BC_TYPES.H" -#include "AMReX_GpuQualifiers.H" -#include "ArrayUtil.hpp" -#include "RadhydroSimulation.hpp" -#include "fextract.hpp" -#include "hydro_system.hpp" -#include "matplotlibcpp.h" -#include "radiation_system.hpp" -#include "test_hydro_highmach.hpp" -#include "valarray.hpp" - -double get_array(amrex::Array4 const &a, int i, int j, int k, - int n) { - return a(i, j, k, n); -} - -struct ShocktubeProblem {}; - -template <> struct EOS_Traits { - static constexpr double gamma = 1.4; - static constexpr bool reconstruct_eint = false; -}; - -AMREX_GPU_MANAGED static amrex::Real rho0 = 1.0; -AMREX_GPU_MANAGED static amrex::Real v0 = 1.0; -AMREX_GPU_MANAGED static amrex::Real E0 = 4.0e-6; - -AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto -ExactSolution(int i, double t, double prob_lo, double dx) - -> quokka::valarray { - const auto gamma = HydroSystem::gamma_; - - // compute point values of primitive variables - auto primVal = [=](int i, double t) { - const double x = prob_lo + (i + 0.5) * dx; - const double rho = rho0 / (1. + v0 * t); - const double vx = v0 * x / (1. + v0 * t); - const double Eint = E0 / std::pow(1. + v0 * t, gamma); - quokka::valarray vars = {rho, vx, Eint}; - return vars; - }; - - // compute point values of conserved variables - auto consVal = [=](quokka::valarray primVars) { - const double rho = primVars[0]; - const double vx = primVars[1]; - const double Eint = primVars[2]; - quokka::valarray vars = {rho, rho * vx, - Eint + 0.5 * rho * vx * vx}; - return vars; - }; - - // compute cell-averaged values of variables - auto q0 = consVal(primVal(i, t)); - auto del_sq = - consVal(primVal(i - 1, t)) - 2.0 * q0 + consVal(primVal(i + 1, t)); - quokka::valarray avg = q0 + del_sq / 24.; - return avg; -} - -template <> -void RadhydroSimulation::setInitialConditionsAtLevel( - int lev) { - amrex::GpuArray dx = geom[lev].CellSizeArray(); - amrex::GpuArray prob_lo = - geom[lev].ProbLoArray(); - const int ncomp = ncomp_; - - for (amrex::MFIter iter(state_new_[lev]); iter.isValid(); ++iter) { - const amrex::Box &indexRange = iter.validbox(); // excludes ghost zones - auto const &state = state_new_[lev].array(iter); - - amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { - auto const avg = ExactSolution(i, 0., prob_lo[0], dx[0]); - - for (int n = 0; n < ncomp; ++n) { - state(i, j, k, n) = 0.; - } - state(i, j, k, HydroSystem::density_index) = avg[0]; - state(i, j, k, HydroSystem::x1Momentum_index) = avg[1]; - state(i, j, k, HydroSystem::x2Momentum_index) = 0.; - state(i, j, k, HydroSystem::x3Momentum_index) = 0.; - state(i, j, k, HydroSystem::energy_index) = avg[2]; - }); - } - - // set flag - areInitialConditionsDefined_ = true; -} - -template <> -AMREX_GPU_DEVICE AMREX_FORCE_INLINE void -AMRSimulation::setCustomBoundaryConditions( - const amrex::IntVect &iv, amrex::Array4 const &consVar, - int /*dcomp*/, int numcomp, amrex::GeometryData const &geom, - const amrex::Real time, const amrex::BCRec * /*bcr*/, int /*bcomp*/, - int /*orig_comp*/) { -#if (AMREX_SPACEDIM == 1) - auto i = iv.toArray()[0]; - int j = 0; - int k = 0; -#endif -#if (AMREX_SPACEDIM == 2) - auto [i, j] = iv.toArray(); - int k = 0; -#endif -#if (AMREX_SPACEDIM == 3) - auto [i, j, k] = iv.toArray(); -#endif - - const double prob_lo = geom.ProbLo(0); - const double dx = geom.CellSize(0); - - auto const avg = ExactSolution(i, time, prob_lo, dx); - - consVar(i, j, k, HydroSystem::density_index) = avg[0]; - consVar(i, j, k, HydroSystem::x1Momentum_index) = avg[1]; - consVar(i, j, k, HydroSystem::x2Momentum_index) = 0.; - consVar(i, j, k, HydroSystem::x3Momentum_index) = 0.; - consVar(i, j, k, HydroSystem::energy_index) = avg[2]; - - consVar(i, j, k, RadSystem::radEnergy_index) = 0; - consVar(i, j, k, RadSystem::x1RadFlux_index) = 0; - consVar(i, j, k, RadSystem::x2RadFlux_index) = 0; - consVar(i, j, k, RadSystem::x3RadFlux_index) = 0; -} - -template <> -void RadhydroSimulation::computeReferenceSolution( - amrex::MultiFab &ref, - amrex::GpuArray const &dx, - amrex::GpuArray const &prob_lo) { - - const double t = stopTime_; - - // fill reference solution multifab - for (amrex::MFIter iter(ref); iter.isValid(); ++iter) { - const amrex::Box &indexRange = iter.validbox(); - auto const &stateExact = ref.array(iter); - auto const ncomp = ref.nComp(); - - amrex::ParallelFor(indexRange, [=](int i, int j, int k) noexcept { - auto const avg = ExactSolution(i, t, prob_lo[0], dx[0]); - - for (int n = 0; n < ncomp; ++n) { - stateExact(i, j, k, n) = 0.; - } - stateExact(i, j, k, HydroSystem::density_index) = - avg[0]; - stateExact(i, j, k, HydroSystem::x1Momentum_index) = - avg[1]; - stateExact(i, j, k, HydroSystem::x2Momentum_index) = 0.; - stateExact(i, j, k, HydroSystem::x3Momentum_index) = 0.; - stateExact(i, j, k, HydroSystem::energy_index) = avg[2]; - }); - } - -#ifdef HAVE_PYTHON - - // Plot results - auto [position, values] = fextract(state_new_[0], geom[0], 0, 0.5); - auto [pos_exact, val_exact] = fextract(ref, geom[0], 0, 0.5); - auto const nx = position.size(); - - if (amrex::ParallelDescriptor::IOProcessor()) { - std::vector x; - std::vector d_final; - std::vector vx_final; - std::vector Eint_final; - std::vector d_exact; - std::vector vx_exact; - std::vector Eint_exact; - - for (int i = 0; i < nx; ++i) { - amrex::Real const this_x = position.at(i); - x.push_back(this_x); - - { - const auto rho = - val_exact.at(HydroSystem::density_index).at(i); - const auto xmom = - val_exact.at(HydroSystem::x1Momentum_index).at(i); - const auto E = - val_exact.at(HydroSystem::energy_index).at(i); - const auto vx = xmom / rho; - const auto Eint = E - 0.5 * rho * (vx * vx); - d_exact.push_back(rho); - vx_exact.push_back(vx); - Eint_exact.push_back(Eint); - } - - { - const auto frho = - values.at(HydroSystem::density_index).at(i); - const auto fxmom = - values.at(HydroSystem::x1Momentum_index).at(i); - const auto fE = - values.at(HydroSystem::energy_index).at(i); - const auto fvx = fxmom / frho; - const auto fEint = fE - 0.5 * frho * (fvx * fvx); - d_final.push_back(frho); - vx_final.push_back(fvx); - Eint_final.push_back(fEint); - } - } - - std::unordered_map d_args; - std::map dexact_args; - d_args["label"] = "density"; - d_args["color"] = "black"; - dexact_args["label"] = "density (exact)"; - matplotlibcpp::scatter(x, d_final, 10.0, d_args); - matplotlibcpp::plot(x, d_exact, dexact_args); - matplotlibcpp::legend(); - matplotlibcpp::title(fmt::format("t = {:.4f}", tNew_[0])); - matplotlibcpp::save("./hydro_highmach.pdf"); - - std::unordered_map e_args; - std::map eexact_args; - d_args["label"] = "internal energy"; - d_args["color"] = "black"; - dexact_args["label"] = "internal energy (exact)"; - matplotlibcpp::clf(); - matplotlibcpp::scatter(x, Eint_final, 10.0, d_args); - matplotlibcpp::plot(x, Eint_exact, dexact_args); - matplotlibcpp::legend(); - matplotlibcpp::title(fmt::format("t = {:.4f}", tNew_[0])); - matplotlibcpp::save("./hydro_highmach_eint.pdf"); - } -#endif -} - -auto problem_main() -> int { - // Problem parameters - // const int nx = 64; - // const double Lx = 2.0; - const double CFL_number = 0.1; - const double max_time = 1.0; - const int max_timesteps = 1000; - - // Problem initialization - const int nvars = RadhydroSimulation::nvarTotal_; - amrex::Vector boundaryConditions(nvars); - for (int n = 0; n < nvars; ++n) { - boundaryConditions[0].setLo(0, amrex::BCType::int_dir); // periodic - boundaryConditions[0].setHi(0, amrex::BCType::int_dir); - for (int i = 1; i < AMREX_SPACEDIM; ++i) { - boundaryConditions[n].setLo(i, amrex::BCType::int_dir); // periodic - boundaryConditions[n].setHi(i, amrex::BCType::int_dir); - } - } - - RadhydroSimulation sim(boundaryConditions); - sim.is_hydro_enabled_ = true; - sim.is_radiation_enabled_ = false; - sim.cflNumber_ = CFL_number; - sim.reconstructionOrder_ = 3; - sim.stopTime_ = max_time; - sim.maxTimesteps_ = max_timesteps; - sim.computeReferenceSolution_ = true; - sim.plotfileInterval_ = -1; - - // Main time loop - sim.setInitialConditions(); - sim.evolve(); - - // Compute test success condition - int status = 0; - const double error_tol = 1.0e-4; - if (sim.errorNorm_ > error_tol) { - status = 1; - } - return status; -} diff --git a/src/HydroHighMach/test_hydro_highmach.hpp b/src/HydroHighMach/test_hydro_highmach.hpp deleted file mode 100644 index 751e703c7..000000000 --- a/src/HydroHighMach/test_hydro_highmach.hpp +++ /dev/null @@ -1,28 +0,0 @@ -#ifndef TEST_HYDRO_SHOCKTUBE_HPP_ // NOLINT -#define TEST_HYDRO_SHOCKTUBE_HPP_ -//============================================================================== -// TwoMomentRad - a radiation transport library for patch-based AMR codes -// Copyright 2020 Benjamin Wibking. -// Released under the MIT license. See LICENSE file included in the GitHub repo. -//============================================================================== -/// \file test_hydro_shocktube.hpp -/// \brief Defines a test problem for a shock tube. -/// - -// external headers -#ifdef HAVE_PYTHON -#include "matplotlibcpp.h" -#endif -#include -#include - -// internal headers - -#include "hydro_system.hpp" -extern "C" { - #include "interpolate.h" -} - -// function definitions - -#endif // TEST_HYDRO_SHOCKTUBE_HPP_ From df7d4a2a3465af91c9cacce4a896627c19c1bb12 Mon Sep 17 00:00:00 2001 From: Benjamin Wibking Date: Tue, 7 Jun 2022 15:06:44 +1000 Subject: [PATCH 12/13] add ipynb for Sedov test --- extern/sedov.ipynb | 120 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 120 insertions(+) create mode 100644 extern/sedov.ipynb diff --git a/extern/sedov.ipynb b/extern/sedov.ipynb new file mode 100644 index 000000000..695b46020 --- /dev/null +++ b/extern/sedov.ipynb @@ -0,0 +1,120 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "yt : [INFO ] 2022-06-07 15:05:45,345 Parameters: current_time = 1.0\n", + "yt : [INFO ] 2022-06-07 15:05:45,346 Parameters: domain_dimensions = [128 128 128]\n", + "yt : [INFO ] 2022-06-07 15:05:45,347 Parameters: domain_left_edge = [0. 0. 0.]\n", + "yt : [INFO ] 2022-06-07 15:05:45,347 Parameters: domain_right_edge = [1.2 1.2 1.2]\n" + ] + }, + { + "data": { + "text/html": [ + "
" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import yt\n", + "ds = yt.load(\"../tests/plt14652\")\n", + "my_sphere = ds.sphere([0., 0., 0.], 1.2)\n", + "plot = yt.ProfilePlot(my_sphere, \"radius\", [(\"boxlib\", \"gasDensity\")], weight_field=\"ones\")\n", + "plot.set_log((\"boxlib\", \"gasDensity\"), False)\n", + "plot.set_log(\"radius\", False)" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "profile = plot.profiles[0]\n", + "r = profile.x\n", + "rho = profile[(\"boxlib\",\"gasDensity\")]" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Text(0, 0.5, 'density (g cm$^{-3}$)')" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEGCAYAAAB2EqL0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAoVElEQVR4nO3deXhcd33v8fdXu+WR7dhabMt2HFvK4iUhwYRAAoQAJQlpQoG0CZTt4TYklN72woUL7SWEdHkKfcq90LDULBdSwtKyxSE2FEraQCCL4yTWyE6w48SOJdmSN23WNprv/WPO2GNFsmek0ZxZPq/nmUdzzvnNOd9j2fP1bzm/n7k7IiIiZ1IWdgAiIlIYlDBERCQtShgiIpIWJQwREUmLEoaIiKSlIuwAsqW+vt5XrlwZdhgiIgXl8ccfP+TuDemULZqEsXLlSrZu3Rp2GCIiBcXM9qZbVk1SIiKSFiUMERFJixKGiIikJecJw8xqzOxRM3vKzNrN7FOTlKk2s++Z2W4ze8TMVuY6ThEROVUYNYwR4Cp3vwh4CXC1mV02ocz7gKPu3gL8H+DTuQ1RREQmynnC8ISBYLMyeE2cAfEG4JvB++8DrzMzy1GIIiIyiVD6MMys3MyeBLqBn7v7IxOKNAMvALh7DOgFFuU0SBEROUUoCcPdx939JcAy4FIzWzed85jZLWa21cy29vT0ZDVGESk9B/uG2dLWFXYYeSvUUVLufgx4ALh6wqEOYDmAmVUA84HDk3x+o7tvcPcNDQ1pPagoIjKlrzy4h9vu2Ubv8bGwQ8lLYYySajCzBcH7OcAbgKcnFNsEvDt4/zbgl66VnkRklrV19AKwu6c/5EjyUxg1jCXAA2a2HXiMRB/GT8zsTjO7PijzNWCRme0GPgR8LIQ4RaSExOPOjs4+AHZ3D5yhdGnK+VxS7r4duHiS/benvB8GbsxlXCJS2vYdOU7/SAyAXQeVMCajJ71FRIBoZ6I5ak5lObtUw5iUEoaICIn+i8py47XnN6hJagpKGCIiQHtHH+ctrmPNknl0HBtiMGiekpOUMESk5Lk70c5e1i2dT0tjBIA9PYMhR5V/lDBEpOR1HBvi2PEx1jbPp6WxDoBd3RpaO5EShoiUvGhHYjjtuqXzOHtRLRVlpn6MSShhiEjJa+/spbzMuGDJPCrLyzinfq5GSk1CCUNESl5bRy+tjRFqKssBaGmM8KwSxosoYYhISXN3oh29rF06/8S+1sYIzx8eZCQ2HmJk+UcJQ0RKWnf/CIcGRlnXPO/EvtWNEeIOzx3SSKlUShgiUtKiwYSD65pTaxiJkVLq+D6VEoaIlLRoRx9mcMGSkzWMVQ1zKTPNKTWREoaIlLRoZy/n1M8lUn1yLtaaynKWL6xld48SRiolDBEpae0dvaxPaY5Kam2MsFs1jFMoYYhIyTo8MEJn7zDrlr44YaxujPDcoUFi4/EQIstPShgiUrKiwYJJa1NGSCW1NtYxOh5n35HjuQ4rbylhiEjJSo6QWjtJDSM5CaFGSp2khCEiJau9s5cVC2uZP6fyRceSCUNThJykhCEiJSva0XfKA3upItUVLJlfoxpGCiUMESlJvcfH2Hfk+CkP7E3U0hhRwkihhCEiJak9WMN7shFSScmEEY97rsLKa0oYIlKSop3JDu/Jm6QgMVJqaGyczt6hXIWV15QwRKQkRTv6WDq/hkWR6inLtDap4zuVEoaIlKRoZy9rT9N/AdDSkEgYWhsjQQlDRErOwEiM5w4Nnrb/AuCsuVXUR6o0CWEg5wnDzJab2QNmtsPM2s3szycpc6WZ9ZrZk8Hr9lzHKSLFa2dXH+6wftnU/RdJqxsimoQwUHHmIlkXAz7s7tvMrA543Mx+7u47JpT7lbtfF0J8IlLk2vafeYRUUmtThE1PduLumNlsh5bXcl7DcPcud98WvO8HdgLNuY5DREpXtLOXhrpqGufVnLFsS0OEvuEYPf0jOYgsv4Xah2FmK4GLgUcmOfwKM3vKzLaY2dopPn+LmW01s609PT2zGaqIFJH2jj7WnWY4barWJq2+lxRawjCzCPAD4C/cvW/C4W3A2e5+EfBPwI8nO4e7b3T3De6+oaGhYVbjFZHiMDQ6zq7u/tM+4Z1Kc0qdFErCMLNKEsniHnf/4cTj7t7n7gPB+81ApZnV5zhMESlCTx/oI+6Tz1A7mca6aupqKlTDIJxRUgZ8Ddjp7p+doszioBxmdimJOA/nLkoRKVbJNTCmmnRwIjOjtTHCru7+2QyrIIQxSupy4J1Am5k9Gez7S2AFgLt/GXgbcJuZxYAh4CZ312QuIjJj7R29nFVbSfOCOWl/pqUxwi+fVj9pzhOGu/8aOO3YNHe/C7grNxGJSClp6+hlXfP8jIbItjbW8a9b93Ps+CgLaqtmMbr8pie9RaRkjMTG+d3B/rT7L5K0+l6CEoaIlIxdBwcYG/e0+y+SNFIqQQlDREpGcg3vdJ7wTtW8YA5zKstLfk4pJQwRKRnRzl7qqitYsbA2o8+VlRmrG+eW/JxSShgiUjKiHX2sbZ5HWVnmc0K1NETYfbC0h9YqYYhISYiNx9nZ1Zdxc1RSa1Mdnb3DDIzEshxZ4VDCEJGSsLtngJFYPO0pQSZKdnyX8mJKShgiUhKiHZk94T2RhtYqYYhIiYh29DKnspxz6iPT+vzZC2upLLeSHlqrhCEiJaG9s5c1S+dRPo0Ob4CK8jLOqZ+rGoaISDGLx532zj7WT7P/Iqm1sY7dJTwJoRKGiBS9PYcGOT46zto0F02ayurGCPuOHGd4bDxLkRUWJQwRKXrtncET3jOuYUSIOzx3aDAbYRUcJQwRKXrRjl6qKspOjHSarlKfU0oJQ0SKXrSjjwsW11FZPrOvvHPq51JmpTu0VglDRIqauxPt7GXtDJujAGoqy1mxsLZkO76VMESkqL1wZIj+4di0pwSZqKWxTjUMEZFi1BZMaT7TIbVJrU0Rnjs0SGw8npXzFRIlDBEpatHOXirKjHMXz6zDO6mlIcLYuLP3yPGsnK+QKGGISFGLdvRyblMd1RXlWTlfa1MwUqoEF1NSwhCRouWeeMJ7uhMOTmZ1Q3ISwtLr+J5WwjCzuWaWnXQtIjJLunqHOTI4OuMH9lLNra6gecGckuz4TithmFmZmb3dzO43s27gaaDLzHaY2T+YWcvshikikrnkGt5rszRCKml1Y6QkH95Lt4bxALAa+Diw2N2Xu3sjcAXwMPBpM/vjWYpRRGRaop19lBmsWZK9JilITBHybM8A8bhn9bz5riLNcq9397GJO939CPAD4AdmVpnOicxsOXA30AQ4sNHdPzehjAGfA64FjgPvcfdtacYqIgIkahgtjRHmVGW3Bb2lMcLwWJyOY0MsX1ib1XPns7RqGJMli+mUCcSAD7v7GuAy4E/NbM2EMtcArcHrFuBLaZ5bROSEaEdv1h7YS9VaoqvvnTFhmNmtZna3md1kZj8xs9tmckF370rWFty9H9gJNE8odgNwtyc8DCwwsyUzua6IlJbuvmG6+0eyMiXIRCcnISytkVLp1DCuAt4NvNPdrwMuytbFzWwlcDHwyIRDzcALKdv7eXFSwcxuMbOtZra1p6cnW2GJSBFo7wzW8J7hGhiTWVBbRX2kWjWMSRx2dwe+HGyPZOPCZhYh0f/xF+7eN51zuPtGd9/g7hsaGhqyEZaIFInkCKk1s5AwINEsVWojpdJJGJ8DcPf7gu0fzvSiQQf5D4B73H2y83UAy1O2lwX7RETSEu3s5Zz6udTVpDUeJ2MtjRF2Hxwg8f/p0nDGhOHuTwOYWX2w/V8zuWAwAuprwE53/+wUxTYB77KEy4Bed++ayXVFpLREO/qy+sDeRK1NEfpHYnT3Z6XRpSBk8qT317N0zcuBdwJXmdmTwevaoHP91qDMZmAPsBv4CvCBLF1bRErAkcFROo4NzUr/RVJLQ+nNKZXucxgAlo0Luvuvz3SuoM/kT7NxPREpPdlaw/t0WppOzil1RWv9rF0nn2RSwyidhjoRKWjRjsQ4mrWzWMNoiFQzr6aipDq+M0kYWalhiIjMtmhnL8vOmsOC2qpZu4aZ0dpUWqvvZZIwPj5rUYiIZFH7LD3hPVFLQ0QJYzLuHp3NQEREsqFveIznDx9n/bLZTxitTREOD45yZHB01q+VDzJaD8PMNpjZj8xsm5ltN7M2M9s+W8GJiGSqPQf9F0ktJTanVCajpADuAT4CtAGltwK6iOS95AipbK+BMZnUOaUuPWfhrF8vbJkmjB533zQrkYiIZEG0o5fF82poqKue9WstnT+H2qpy1TCm8Ekz+yrwH6TMKTXF9B4iIjkXzfIa3qdTVmasLqGO70wTxnuB84FKTjZJOVmYX0pEZKaOj8Z4tmeAN63P3WoIrY0RfrvncM6uF6ZME8bL3P28WYlERGSGdnb14Q7rZ/EJ74lWN0b44RMd9A+PzdpEh/kio1FSwG8mWR1PRCQvtO2f/SlBJkquvvdsz2DOrhmWTGsYlwFPmtlzJPowjMTUTxdmPTIRkQxFO/uoj1TRNG/2O7yTToyUOtjPS5YvyNl1w5Bpwrh6VqIQEcmCaEcva5fOJ7GKQm6sWFhLVXkZu3uKv+M70yapO0msTbHX3fcCfcAnsx+WiEhmhsfG2dU9kLMRUkkV5WWsapjL7hKY5jzThHGhux9Lbrj7URJrcouIhOqZA/2Mxz0nc0hNtLpElmvNNGGUmdlZyQ0zW0jmzVoiIlkXzcEaGFNpbYzwwtHjDI+N5/zauZTpl/0/Ar81s38Ltm8E/ja7IYmIZC7a0cv8OZUsO2tOzq/d0hjBHZ7tGcjJlCRhyaiG4e53A28BDgavt7j7v8xGYCIimUis4T0vpx3eSa2NdUDxT0KYcXOSu+8AdsxCLCIi0zIai/PMgX7ee/nKUK6/sr6WMiv+hJFpH4aISN7Z1d3P6HictSH0XwBUV5SzctFcJQwRkXyXXANjXQ7WwJhKKYyUUsIQkYIX7ewlUl3BykVzQ4uhtTHC84cGGRsv3qWCMurDMLMPTbK7F3jc3Z/MSkQiIhlq6+hlzdJ5lJXlvsM7qbUpQizu7D08SEvQCV5sMq1hbABuBZqD1/tJTBfyFTP7aJZjExE5o9h4nJ1dfaE8sJeqpSGRJHYV8RPfmSaMZcAl7v5hd/8w8FKgEXg18J50TmBmXzezbjOLTnH8SjPrNbMng9ftGcYoIiVkz6FBhsfiOZ8SZKLVjYnmsGLu+M50WG0jKSvtAWNAk7sPmdnIFJ+Z6BvAXcDdpynzK3e/LsPYRKQERTvCe8I7VW1VBc0L5hR1x3emCeMe4BEzuzfY/n3g22Y2lzSfzXD3B81sZYbXFRGZVLSjj5rKMlbVh9fhndTaVNzLtWb6pPdfA7cAx4LXre5+p7sPuvs7shjXK8zsKTPbYmZrpypkZreY2VYz29rT05PFy4tIoYh29nLBknlUlIc/6LOlIcKzPQOMxz3sUGZFWjUMMzN3dwB33wpsPV2ZGdoGnO3uA2Z2LfBjoHWygu6+EdgIsGHDhuL8DYnIlOJxZ0dnH39wcXPYoQCJGsZILE7H0SFWLKoNO5ysSzclP2Bmf2ZmK1J3mlmVmV1lZt8E3p2NgNy9z90HgvebgUozq8/GuUWkuOw5NMDASCz0Du+kE6vvdfeHHMnsSDdhXA2MA98xs04z2xEs07oLuBn4v+7+jWwEZGaLLZg9zMwuDWI8nI1zi0hx+fmObgAub8mP/1Mmh9YWaz9GWk1S7j4MfBH4oplVAvXAUOpiSukys+8AVwL1ZrafxIp9lcF1vgy8DbjNzGLAEHBTlpq6RKTIbIl2cdGy+Sw7Kz+af+bXVtJQV120I6WmM1vtGNA13Qu6+81nOH4XiWG3IiJTeuHIcbbv7+Vj15wfdiinaC3iOaXCH1YgIjINP40eAOCadYtDjuRUrY0Rnu0eoBgbRpQwRKQgbY52sXbpPM4OccLBybQ0RhgYiXGgbzjsULIuo4QRjJQ668wlRURmT+exIZ7Yd4xr1y8JO5QXaSni1fcyrWE0AY+Z2b+a2dXJ0UwiIrmUr81RkDK0tggnIcz0Se//TeIhuq+RmGxwl5n9nZmtnoXYREQmtSXaxfmL61jVEAk7lBepj1SxoLaS3T0lnjAAgiGuB4JXDDgL+L6ZfSbLsYmIvMjBvmG27j3KNevyrzkKwMxoaYiwu9RrGGb252b2OPAZ4CFgvbvfRmKa87fOQnwiIqf4WfsB3OHa9fnXHJXU2hQpyhpGps9hLATe4u57U3e6e9zMNB25iMy6zW1dtDRGaG3K31XtVjdEODL4AocHRlgUqQ47nKzJtEmqZmKyMLNPA7j7zqxFJSIyiZ7+ER597gjX5mFnd6pkMiu2B/gyTRhvmGTfNdkIRETkTP59xwHiDtfk4XDaVK3BSKliG1qb7vTmtwEfAFab2XYgOZy2jkRfhojIrNvSdoBz6udy/uL8bY4CWDK/hrlV5aWZMEistLcF+DvgYyQShgP97n50lmITETnhyOAov91zmPe/ehX5/giYmdHSWHyr76WbMDa7+xVmdj2Q2rmdXDcpPyajF5Gi9fMdBxiPe14+3T2Z1Y0RHtp9KOwwsiqtPgx3vyL4GXH3eSmvOiULEcmFzW0HWL5wDmuXFsZXTmtjHQf7RugbHgs7lKzR5IMikvd6j4/x0O5DXLtuSd43RyW1FGHHd6YP7t1oZnXB+0+Y2Q/N7JLZCU1EJOHnOw8Si3vej45KVYwjpTKtYXzC3fvN7ArgdSTmlPpS9sMSETlpS1sXS+fXcNGy+WGHkrblC2upqigr6YQxHvx8E7DR3e8HqrIbkojISX3DY/xq1yGuWV84zVEA5WXGqvq57DrYH3YoWZNpwugws38GbgI2m1n1NM4hIpK2X+7sZnQ8ntdzR02ltamuqOaUyvTL/g+BnwG/5+7HSMxU+5FsByUikrS5rYumedVcvLzw1m5raYiw/+gQQ6PjZy5cADKdfHAcqAFuNLPUz/579kISEUkYGInxn7/r4e2XrqCsrHCao5JamyK4w7M9A6xrLpz+l6lkWsO4F7iexDoYgykvEZGse+DpbkZj8bxcWS8dxTa0NtMaxjJ3v3pWIhERmWBLtIv6SDUbVi4MO5RpWbloLuVlVjQJI9Maxm/MbP2sRCIikuL4aIwHnu7h6nVNlBdgcxRAVUUZZy+qZVd3cYyUyjRhXAFsM7NnzGy7mbUFs9emzcy+bmbdZhad4riZ2efNbHdwDT0YKFKC/uuZHobGxrk2T5diTVdrEU1CmGmTVDbWvvgGcBdw92mu0Rq8Xk7iwcCXZ+G6IlJANkcPsHBuFZeeU5jNUUktjRF+sTPRF1NVUdhPIWSaMPYB7wBWufudZrYCWAzsPf3HTnL3B81s5WmK3ADc7e4OPGxmC8xsibt3ZRiriBSo4bFxfrnzINe/ZCkV5YX9JdvaWMd43Pn2I3tpnFdzYn9qI9uLn0e0SY9NLJZ8kPHCZfNpSjn3bMk0YXwRiANXAXcC/cAPgJdlMaZm4IWU7f3BvhclDDO7BbgFYMWKFVkMQUTC9ODvehgcHeeaAm+OAk4Mp73jvh2zdo0vvP0S3nTh7P9ZZZowXu7ul5jZEwDuftTMQpsaxN03AhsBNmzY4GHFISLZtSV6gPlzKnnF6kVhhzJjLY0RHvrYVQwMxwBwTv2qcp/8/cSyE4+lWn5W7YzjTEemCWPMzMpJrLaHmTWQqHFkUwewPGV7WbBPRErASGycX+w4yNXrFlNZ4M1RSc0L5oQdQlZk+tv4PPAjoMnM/hb4NYllW7NpE/CuYLTUZUCv+i9ESsdDuw/RPxIrmJX1SklGNQx3v8fMHicxtTnAm919ZybnMLPvAFcC9Wa2H/gkUBmc/8vAZuBaYDdwHHhvJucXkcK2ue0AdTUVvLKl8Jujik1aCcPMPjTFoWvM7Bp3/2y6F3T3m89w3IE/Tfd8IlI8RmNx/r39AG+4oInqivKww5EJ0q1h1AU/zyMxImpTsP37wKPZDkpEStNv9xymbzhWUCvrlZK0Eoa7fwrAzB4ELnH3/mD7DuD+WYtORErKlrYu5laV86rW+rBDkUlk2undBIymbI8G+0REZiQ2Hudn7Qd43QVN1FSqOSofZTqs9m7gUTP7UbD9ZhJTfYiIzMgjzx3h6PGxglxZr1RkOkrqb81sC/CqYNd73f2J7IclIqVmc1sXcyrLec25jWGHIlPItIaBu28Dts1CLCJSosbjzs/aD3DV+Y3MqVJzVL4qjscoRaSgPfb8EQ4NjHKNmqPymhKGiIRuS1sX1RVlvPY8NUflMyUMEQlVPO5siR7gyvMamFudcSu55JAShoiEatu+o3T3j2juqAKghCEiodrcdoCq8jKuOl/NUflOCUNEQpNojuri1efWU1dTGXY4cgZKGCISmqf2H6Ord7goVtYrBUoYIhKaLdEDVJYbr79AMwwVAiUMEQmFu7O5rYvLW+qZX6vmqEKghCEioYh29LH/6JBGRxUQJQwRCcXmaBcVZcbvrVFzVKFQwhCRnHN3trR18YrVi1hQWxV2OJImJQwRybmdXf08f/i4mqMKjBKGiOTUaCzO7fdGqaksU3NUgdHELSKSU3f+pJ2te4/y+ZsvZlGkOuxwJAOqYYhIznzvsX186+F93PLqVVx/0dKww5EMKWGISE48se8on/hxO1e01PPRN54XdjgyDUoYIjLruvuHufVbj9M0v5p/uvliKsr11VOI1IchIrNqNBbnA9/aRu/QGD+87XLOmqthtIUqlDRvZleb2TNmttvMPjbJ8feYWY+ZPRm8/lsYcYrIzCU7uT/ztotYs3Re2OHIDOS8hmFm5cAXgDcA+4HHzGyTu++YUPR77v7BXMcnItmT7OR+vzq5i0IYNYxLgd3uvsfdR4HvAjeEEIeIzKJtQSf3q1rr+ejV54cdjmRBGAmjGXghZXt/sG+it5rZdjP7vpktn+xEZnaLmW01s609PT2zEauITEN3/zC3pXRyl5dZ2CFJFuTrUIX7gJXufiHwc+CbkxVy943uvsHdNzQ0NOQ0QBGZXLKTu28oxsZ3btBcUUUkjITRAaTWGJYF+05w98PuPhJsfhV4aY5iE5EZ+tR9iU7uf7jxQi5Yok7uYhJGwngMaDWzc8ysCrgJ2JRawMxSZyS7HtiZw/hEZJq+++g+7nlkH+9/zSquu1Cd3MUm56Ok3D1mZh8EfgaUA19393YzuxPY6u6bgP9uZtcDMeAI8J5cxykimdm27yi33xt0cr9RndzFyNw97BiyYsOGDb5169awwxApSd19w1z3T7+mprKcTR+8XP0WBcTMHnf3DemUzddObxEpEKOxOLfds43+4Rgb3/VSJYsipqlBRGRG7rivncf3HuWut1/M+YvVyV3MVMMQkWn7zqP7+PYj+7j1NavVyV0ClDBEZFoe33uU2++N8qrWej6i6cpLghKGiGTsYF/iSe4l8+foSe4SooQhIhkZHIlx27ceVyd3CVKnt4ik7Rc7DnL7vVG6+oa56+ZL1MldYpQwROSMDvQOc8emdn7afoBzmyJ8/+2v4KVnLww7LMkxJQwRmdJ43Lnnkb185qfPMDYe5yNvPI8/edUqqirUml2KlDBEZFI7Ovv4+I/aeOqFY7yqtZ6/efM6zl40N+ywJERKGCJyiuOjMT73i1189dfPsWBOJZ+76SVcf9FSzDQSqtQpYYjICQ88080nfhxl/9EhbnrZcj52zfkaBSUnKGGICN19w9z5kx38ZHsXqxvm8r1bLuPlqxaFHZbkGSUMkRIWjzvffnQfn/7p04zE4nzoDefy/tesorqiPOzQJA8pYYiUqGcO9POXP2rj8b1HeeXqRfzNm9exqiESdliSx5QwREqIu7N9fy8/eqKDbz28l7qaCv7xxot4yyXN6tSWM1LCECly7s4zB/u576lO7nuqi31HjlNZbvzBxc18/NoLWDhXndqSHiUMkSK1p2eA+57q4ifbO9nVPUB5mfHK1Yv44FUtvHHNYubXVoYdohQYJQyRIvLCkePc39bFfU910t7Zhxm8bOVC/vrN67hm3WLqI9VhhygFTAlDpMAd7Bvm/u1d3Le9kyf2HQPgJcsX8Inr1vCm9UtYPL8m3AClaChhiBQQd2ffkeNEO/po7+xl696jPPb8EdxhzZJ5/K+rz+e6C5ewfGFt2KFKEVLCEMlTsfE4ew4NEu3opb2zj2hHLzs6++gfiQFQUWact7iOP39dK9dduJSWRg2JldmlhCGSBwZHYuzpGaS9s5doZy/Rjj6ePtDH8FgcgJrKMs5fPI8bLl7KuqXzWdc8n9amiB6wk5xSwhCZZfG40zMwQsexITqODtF5LPHqODYc/Byid2jsRPm66grWLJ3HO15+NmuXzmNd83xW1c+lolxTiku4QkkYZnY18DmgHPiqu//9hOPVwN3AS4HDwB+5+/O5jlNkKu5O33CMI4OjHBkc4fDAKIcHRzkyOMrhgcS+rt5hOnuHONA7zNi4n/L5upoKmhfMYemCOVxy9gKaF9SyYmEt65rnsfysWsq0RrbkoZwnDDMrB74AvAHYDzxmZpvcfUdKsfcBR929xcxuAj4N/NFsxDMai1NZbsTiztHjo8TGnZ1dfSw7q5ajx0cZjzuLIlXExp1IdQXj7syfU4l7opmgsryMWNypLDcMI+6JL4YyMyrKrGj/4Xtwn+7gwbaf2HbcU8ue3Ocpn0+WT+48bZngePJg3CEWjxOPw7g74/E44/GT+2LxOONxJxZ3YuPOWDzO+LgTi8cZC37GxhPHx8bjjIzFGYmNMxKLMxKLMxoLtsfiwb5xhsbGOTI4xpHBEY4Mjr4oCSTVVpWzcG4Vi+fVcPHys1i6fg7NZ82heUENS4MkMa9Gz0BI4QmjhnEpsNvd9wCY2XeBG4DUhHEDcEfw/vvAXWZm7j75v9AZ+H8PPcfXH3qOgeEYg6Pj2T49AGZgcMrUCxbsP7l9msRip/yYcN7JP+ec/EJPbDPhzallUr+kJz8fZP9PP/+YQXVFGdUV5YmflSff11SW07yghgub57MwUsWiuVUsDF6L5laf2FdTqX4FKU5hJIxm4IWU7f3Ay6cq4+4xM+sFFgGHUguZ2S3ALQArVqyYVjCrGyJces4iuvuG2XfkOPNqKrm8pZ7KcmN1Q4SGedUMjY5TZsbx0RhlZvQOjWEGI2NxRscTNZSxccfdTyQFd2c8DuPx+Iu+bF/0P/DTxHfK/65POXDql/xk8wDZhDfJ5HJqouLEvsmOv+h8wcHUhGdYSlI8mRinPD7hOmY2IY5TzzHZdY2TNbhTfgY1u/LgVVFuVJSVUVFuVJaVUV5mVJYbFeVlVATHq8rLqK5MJIWKMtOcSiJTKOhOb3ffCGwE2LBhw7T+//v6NU28fk1TVuMSESlGYQy76ACWp2wvC/ZNWsbMKoD5JDq/RUQkJGEkjMeAVjM7x8yqgJuATRPKbALeHbx/G/DL2ei/EBGR9OW8SSrok/gg8DMSw2q/7u7tZnYnsNXdNwFfA/7FzHYDR0gkFRERCVEofRjuvhnYPGHf7Snvh4Ebcx2XiIhMTY+OiohIWpQwREQkLUoYIiKSFiUMERFJixXLaFUz6wH2ZvixeiY8PV4EdE+FQfdUGIrxnuDU+zrb3RvS+VDRJIzpMLOt7r4h7DiySfdUGHRPhaEY7wmmf19qkhIRkbQoYYiISFpKPWFsDDuAWaB7Kgy6p8JQjPcE07yvku7DEBGR9JV6DUNERNKkhCEiImkpiYRhZleb2TNmttvMPjbJ8Woz+15w/BEzWxlCmBlJ454+ZGY7zGy7mf2HmZ0dRpyZONM9pZR7q5m5meX9cMd07snM/jD4XbWb2bdzHWOm0vi7t8LMHjCzJ4K/f9eGEWcmzOzrZtZtZtEpjpuZfT645+1mdkmuY8xUGvf0juBe2szsN2Z20RlP6u5F/SIxhfqzwCqgCngKWDOhzAeALwfvbwK+F3bcWbin1wK1wfvbiuGegnJ1wIPAw8CGsOPOwu+pFXgCOCvYbgw77izc00bgtuD9GuD5sONO475eDVwCRKc4fi2whcQqwZcBj4Qdcxbu6ZUpf++uSeeeSqGGcSmw2933uPso8F3ghgllbgC+Gbz/PvA6y++Fnc94T+7+gLsfDzYfJrGyYT5L5/cE8NfAp4HhXAY3Tenc058AX3D3owDu3p3jGDOVzj05MC94Px/ozGF80+LuD5JYe2cqNwB3e8LDwAIzW5Kb6KbnTPfk7r9J/r0jze+IUkgYzcALKdv7g32TlnH3GNALLMpJdNOTzj2leh+J/x3lszPeU9AMsNzd789lYDOQzu/pXOBcM3vIzB42s6tzFt30pHNPdwB/bGb7Sax782e5CW1WZfpvrtCk9R0RygJKkjtm9sfABuA1YccyE2ZWBnwWeE/IoWRbBYlmqStJ/A/vQTNb7+7Hwgxqhm4GvuHu/2hmryCxeuY6d4+HHZi8mJm9lkTCuOJMZUuhhtEBLE/ZXhbsm7SMmVWQqEYfzkl005POPWFmrwf+Crje3UdyFNt0neme6oB1wH+a2fMk2pE35XnHdzq/p/3AJncfc/fngN+RSCD5Kp17eh/wrwDu/lughsRkd4UsrX9zhcbMLgS+Ctzg7mf8ziuFhPEY0Gpm55hZFYlO7U0TymwC3h28fxvwSw96gvLUGe/JzC4G/plEssj3dnE4wz25e6+717v7SndfSaLN9Xp33xpOuGlJ5+/ej0nULjCzehJNVHtyGGOm0rmnfcDrAMzsAhIJoyenUWbfJuBdwWipy4Bed+8KO6iZMLMVwA+Bd7r779L6UNg9+TkaLXAtif+5PQv8VbDvThJfOJD4C/1vwG7gUWBV2DFn4Z5+ARwEngxem8KOeab3NKHsf5Lno6TS/D0Ziaa2HUAbcFPYMWfhntYAD5EYQfUk8Hthx5zGPX0H6ALGSNT63gfcCtya8nv6QnDPbQXyd+9M9/RV4GjKd8TWM51TU4OIiEhaSqFJSkREskAJQ0RE0qKEISIiaVHCEBGRtChhiIhIWpQwRNJkZlea2U+C99efbkbdDM/7fTNblYXzfNfM8vmhPylwShhS0oIHsTL+d+Dum9z977Nw/bVAubtn42G9LwEfzcJ5RCalhCElx8xWBus53A1EgeVm9iUz2xqsSfGplLJXm9nTZrYNeEvK/veY2V3B+2+Y2dtSjg0EP5eY2YNm9qSZRc3sVZOE8w7g3gnX22ZmT5nZfwT77jCzb5rZr8xsr5m9xcw+E6xj8FMzqww+/ivg9cH0NiJZp4QhpaoV+KK7r3X3vSSeWN4AXAi8xswuNLMa4CvA7wMvBRZneI23Az9z95cAF5F4mnaiy4HHAcysIbjeW939IuDGlHKrgauA64FvAQ+4+3pgCHgTgCcm99sdXEsk65QwpFTt9cS6Bkl/GNQingDWkpje4nzgOXff5YkpEb6V4TUeA95rZncA6929f5IySzg5z9JlwIOemIQQd09dy2CLu4+RmJaiHPhpsL8NWJlSrhtYmmGcImlRwpBSNZh8Y2bnAP8TeJ27XwjcT2J+sXTFCP4tBf0hVXBiAZtXk5jV9Btm9q5JPjuU5rVGgnPGgTE/OadPnFOXKagJzimSdUoYIonV4QaBXjNrIrFcJcDTwEozWx1s3zzF558n0WQFiSajSgBLrKN+0N2/QmKit8nWgd4JtATvHwZeHSQwzGzhNO7lXBL9MiJZp84xKXnu/pSZPUEiQbxAYqZV3H3YzG4B7jez4yQ6lesmOcVXgHvN7CkSTUXJ2suVwEfMbAwYACarYdwflPuFu/cE1/thUFPpBt6Q7n0EyW7I3Q+k+xmRTGi2WpEQmdkc4AHgcncfn+G5/gfQ5+5fy0pwIhOoSUokRO4+BHyS7KwPfQz4ZhbOIzIp1TBERCQtqmGIiEhalDBERCQtShgiIpIWJQwREUmLEoaIiKTl/wND9WFYeS5CNAAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "import matplotlib.pyplot as plt\n", + "plt.plot(r, rho, label='simulation')\n", + "plt.xlabel(r\"radius (cm)\")\n", + "plt.ylabel(r\"density (g cm$^{-3}$)\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "interpreter": { + "hash": "5d21aa510141f516aa4589cf93de6e254d300587727b9d0aa72219f1336bba5f" + }, + "kernelspec": { + "display_name": "Python 3.9.7 ('base')", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.7" + }, + "orig_nbformat": 4 + }, + "nbformat": 4, + "nbformat_minor": 2 +} From 71c74bf37d1c50fa59dd048b11409aca3aa72bbd Mon Sep 17 00:00:00 2001 From: Benjamin Wibking Date: Tue, 7 Jun 2022 15:07:03 +1000 Subject: [PATCH 13/13] update KH test params --- src/HydroKelvinHelmholz/test_hydro2d_kh.cpp | 5 +++-- tests/KelvinHelmholz_512.in | 24 +++++++++++++++++++++ 2 files changed, 27 insertions(+), 2 deletions(-) create mode 100644 tests/KelvinHelmholz_512.in diff --git a/src/HydroKelvinHelmholz/test_hydro2d_kh.cpp b/src/HydroKelvinHelmholz/test_hydro2d_kh.cpp index ef4f3068f..5a235926d 100644 --- a/src/HydroKelvinHelmholz/test_hydro2d_kh.cpp +++ b/src/HydroKelvinHelmholz/test_hydro2d_kh.cpp @@ -136,10 +136,11 @@ auto problem_main() -> int { RadhydroSimulation sim(boundaryConditions); sim.is_hydro_enabled_ = true; sim.is_radiation_enabled_ = false; - sim.stopTime_ = 1.5; + sim.stopTime_ = 20.0; sim.cflNumber_ = 0.4; - sim.maxTimesteps_ = 40000; + sim.maxTimesteps_ = 100000; sim.plotfileInterval_ = 100; + sim.checkpointInterval_ = 1000; // initialize sim.setInitialConditions(); diff --git a/tests/KelvinHelmholz_512.in b/tests/KelvinHelmholz_512.in new file mode 100644 index 000000000..a1c58bf89 --- /dev/null +++ b/tests/KelvinHelmholz_512.in @@ -0,0 +1,24 @@ +# ***************************************************************** +# Problem size and geometry +# ***************************************************************** +geometry.prob_lo = 0.0 0.0 0.0 +geometry.prob_hi = 1.0 1.0 1.0 +geometry.is_periodic = 1 1 1 + +# ***************************************************************** +# VERBOSITY +# ***************************************************************** +amr.v = 0 # verbosity in Amr + +# ***************************************************************** +# Resolution and refinement +# ***************************************************************** +amr.n_cell = 512 512 8 +amr.max_level = 0 # number of levels = max_level + 1 +amr.blocking_factor = 64 # grid size must be divisible by this +amr.max_grid_size = 64 +amr.n_error_buf = 3 # minimum 3 cell buffer around tagged cells +amr.grid_eff = 0.7 # default + +do_reflux = 0 +do_subcycle = 0