From 13eaa749323587becbbb62115144084150a173ab Mon Sep 17 00:00:00 2001 From: neco kriel Date: Fri, 3 Oct 2025 15:54:41 +1000 Subject: [PATCH 01/34] draft hydro version of test --- inputs/hydro_balsara_vortex.in | 35 ++++ .../HydroBalsaraVortex/CMakeLists.txt | 8 + .../test_hydro_balsara_vortex.cpp | 149 ++++++++++++++++++ 3 files changed, 192 insertions(+) create mode 100644 inputs/hydro_balsara_vortex.in create mode 100644 src/problems/HydroBalsaraVortex/CMakeLists.txt create mode 100644 src/problems/HydroBalsaraVortex/test_hydro_balsara_vortex.cpp diff --git a/inputs/hydro_balsara_vortex.in b/inputs/hydro_balsara_vortex.in new file mode 100644 index 0000000000..f29981bc96 --- /dev/null +++ b/inputs/hydro_balsara_vortex.in @@ -0,0 +1,35 @@ +# ***************************************************************** +# 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 = 1 + +# ***************************************************************** +# Resolution and refinement +# ***************************************************************** +amr.n_cell = 32 32 16 +amr.max_level = 0 +amr.max_grid_size = 32 +amr.blocking_factor_x = 32 +amr.blocking_factor_y = 16 +amr.blocking_factor_z = 16 + +do_reflux = 0 +do_subcycle = 0 + +plotfile_interval = -1 +plotfile_prefix = hydro_balsara_vortex_plt +checkpoint_interval = -1 +cfl = 0.3 +stop_time = 2.0 +max_timesteps = 10000 + +hydro.rk_integrator_order = 2 +hydro.reconstruction_order = 5 +hydro.use_dual_energy = 0 diff --git a/src/problems/HydroBalsaraVortex/CMakeLists.txt b/src/problems/HydroBalsaraVortex/CMakeLists.txt new file mode 100644 index 0000000000..009ff341f0 --- /dev/null +++ b/src/problems/HydroBalsaraVortex/CMakeLists.txt @@ -0,0 +1,8 @@ +if (AMReX_SPACEDIM EQUAL 3) + add_executable(test_hydro_balsara_vortex test_hydro_balsara_vortex.cpp ${QuokkaObjSources}) + if(AMReX_GPU_BACKEND MATCHES "CUDA") + setup_target_for_cuda_compilation(test_hydro_balsara_vortex) + endif() + + add_test(NAME HydroBalsaraVortex COMMAND test_hydro_balsara_vortex ../inputs/hydro_balsara_vortex.in ${QuokkaTestParams} WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) +endif() diff --git a/src/problems/HydroBalsaraVortex/test_hydro_balsara_vortex.cpp b/src/problems/HydroBalsaraVortex/test_hydro_balsara_vortex.cpp new file mode 100644 index 0000000000..3c412f6072 --- /dev/null +++ b/src/problems/HydroBalsaraVortex/test_hydro_balsara_vortex.cpp @@ -0,0 +1,149 @@ +//============================================================================== +// Copyright 2025 Neco Kriel. +// Released under the MIT license. See LICENSE file included in the GitHub repo. +//============================================================================== +/// \file test_balsara_vortex_hllc.cpp +/// \brief Hydrodynamic (HLLC) Balsara vortex (B = 0; constant-ρ, prescribed p). +/// + +#include +#include +#include + +#include "AMReX_Array.H" +#include "AMReX_Array4.H" +#include "AMReX_REAL.H" + +#include "QuokkaSimulation.hpp" +#include "grid.hpp" +#include "hydro/EOS.hpp" +#include "physics_info.hpp" +#include "util/BC.hpp" + +struct HydroBalsaraVortex {}; + +template <> struct quokka::EOS_Traits { + static constexpr double gamma = 5.0 / 3.0; + static constexpr double mean_molecular_weight = C::m_u; + static constexpr double boltzmann_constant = C::k_B; +}; + +template <> struct Physics_Traits { + static constexpr bool is_hydro_enabled = true; + static constexpr int numMassScalars = 0; + static constexpr int numPassiveScalars = numMassScalars + 0; + static constexpr bool is_self_gravity_enabled = false; + static constexpr bool is_radiation_enabled = false; + static constexpr bool is_mhd_enabled = false; + static constexpr int nGroups = 1; + static constexpr UnitSystem unit_system = UnitSystem::CGS; +}; + +// background state +constexpr double gamma_gas = quokka::EOS_Traits::gamma; +constexpr double bg_density = 1.0; +constexpr double bg_pressure = 1.0; +// vortex parameters +constexpr double vortex_speed = 5.0 / (2.0 * M_PI); +constexpr double vortex_x0 = 0.5; +constexpr double vortex_y0 = 0.5; +// drift is off by default +constexpr double vortex_drift_x1 = 0.0; +constexpr double vortex_drift_x2 = 0.0; +constexpr double vortex_drift_x3 = 0.0; + +AMREX_GPU_DEVICE +inline void computeVortexSolution( + int i, int j, int k, + amrex::Array4 const &state, + amrex::GpuArray const &dx, + amrex::GpuArray const &prob_lo) +{ + const amrex::Real x1_L = prob_lo[0] + i * dx[0]; + const amrex::Real x2_L = prob_lo[1] + j * dx[1]; + const amrex::Real x1_C = x1_L + 0.5 * dx[0]; + const amrex::Real x2_C = x2_L + 0.5 * dx[1]; + + const double rel_x1 = static_cast(x1_C) - vortex_x0; + const double rel_x2 = static_cast(x2_C) - vortex_y0; + const double radius_sq = rel_x1 * rel_x1 + rel_x2 * rel_x2; + + const double density = bg_density; + const double pressure = bg_pressure - 0.5 * vortex_speed * vortex_speed * gcem::exp(1.0 - radius_sq); + + const double delta_vel_x1 = -rel_x2 * vortex_speed * gcem::exp(0.5 * (1.0 - radius_sq)); + const double delta_vel_x2 = rel_x1 * vortex_speed * gcem::exp(0.5 * (1.0 - radius_sq)); + const double vel_x1 = vortex_drift_x1 + delta_vel_x1; + const double vel_x2 = vortex_drift_x2 + delta_vel_x2; + const double vel_x3 = vortex_drift_x3; + + const double mom_x1 = density * vel_x1; + const double mom_x2 = density * vel_x2; + const double mom_x3 = density * vel_x3; + + const double Eint = pressure / (gamma_gas - 1.0); + const double Ekin = 0.5 * density * (vel_x1 * vel_x1 + vel_x2 * vel_x2 + vel_x3 * vel_x3); + const double Etot = Eint + Ekin; + + state(i, j, k, HydroSystem::density_index) = density; + state(i, j, k, HydroSystem::x1Momentum_index) = mom_x1; + state(i, j, k, HydroSystem::x2Momentum_index) = mom_x2; + state(i, j, k, HydroSystem::x3Momentum_index) = mom_x3; + state(i, j, k, HydroSystem::energy_index) = Etot; + state(i, j, k, HydroSystem::internalEnergy_index) = Eint; +} + +template <> +void QuokkaSimulation::setInitialConditionsOnGrid(quokka::grid const &grid_elem) +{ + const amrex::GpuArray dx = grid_elem.dx_; + const amrex::GpuArray prob_lo = grid_elem.prob_lo_; + const amrex::Array4 &state_cc = grid_elem.array_; + const amrex::Box &indexRange = grid_elem.indexRange_; + const int ncomp_cc = Physics_Indices::nvarTotal_cc; + + amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + for (int icomp = 0; icomp < ncomp_cc; ++icomp) { + state_cc(i, j, k, icomp) = 0.0; + } + computeVortexSolution(i, j, k, state_cc, dx, prob_lo); + }); +} + +template <> +void QuokkaSimulation::computeReferenceSolution( + amrex::MultiFab &ref, + amrex::GpuArray const &dx, + amrex::GpuArray const &prob_lo) +{ + 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, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + for (int icomp = 0; icomp < ncomp; ++icomp) { + stateExact(i, j, k, icomp) = 0.0; + } + computeVortexSolution(i, j, k, stateExact, dx, prob_lo); + }); + } +} + +auto problem_main() -> int +{ + auto BCs_cc = quokka::BC(quokka::BCType::int_dir); // periodic + + QuokkaSimulation sim(BCs_cc); + sim.computeReferenceSolution_ = true; + + sim.setInitialConditions(); + sim.evolve(); + + int status = 0; + const double error_tol = 0.002; + if (sim.errorNorm_ > error_tol) { + status = 1; + } + return status; +} From ab38c06451cc22b1375af708bddd14be47b052db Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 3 Oct 2025 06:02:05 +0000 Subject: [PATCH 02/34] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../test_hydro_balsara_vortex.cpp | 23 ++++++++----------- 1 file changed, 9 insertions(+), 14 deletions(-) diff --git a/src/problems/HydroBalsaraVortex/test_hydro_balsara_vortex.cpp b/src/problems/HydroBalsaraVortex/test_hydro_balsara_vortex.cpp index 3c412f6072..e5fd0b6beb 100644 --- a/src/problems/HydroBalsaraVortex/test_hydro_balsara_vortex.cpp +++ b/src/problems/HydroBalsaraVortex/test_hydro_balsara_vortex.cpp @@ -20,7 +20,8 @@ #include "physics_info.hpp" #include "util/BC.hpp" -struct HydroBalsaraVortex {}; +struct HydroBalsaraVortex { +}; template <> struct quokka::EOS_Traits { static constexpr double gamma = 5.0 / 3.0; @@ -53,11 +54,8 @@ constexpr double vortex_drift_x2 = 0.0; constexpr double vortex_drift_x3 = 0.0; AMREX_GPU_DEVICE -inline void computeVortexSolution( - int i, int j, int k, - amrex::Array4 const &state, - amrex::GpuArray const &dx, - amrex::GpuArray const &prob_lo) +inline void computeVortexSolution(int i, int j, int k, amrex::Array4 const &state, amrex::GpuArray const &dx, + amrex::GpuArray const &prob_lo) { const amrex::Real x1_L = prob_lo[0] + i * dx[0]; const amrex::Real x2_L = prob_lo[1] + j * dx[1]; @@ -67,8 +65,8 @@ inline void computeVortexSolution( const double rel_x1 = static_cast(x1_C) - vortex_x0; const double rel_x2 = static_cast(x2_C) - vortex_y0; const double radius_sq = rel_x1 * rel_x1 + rel_x2 * rel_x2; - - const double density = bg_density; + + const double density = bg_density; const double pressure = bg_pressure - 0.5 * vortex_speed * vortex_speed * gcem::exp(1.0 - radius_sq); const double delta_vel_x1 = -rel_x2 * vortex_speed * gcem::exp(0.5 * (1.0 - radius_sq)); @@ -93,8 +91,7 @@ inline void computeVortexSolution( state(i, j, k, HydroSystem::internalEnergy_index) = Eint; } -template <> -void QuokkaSimulation::setInitialConditionsOnGrid(quokka::grid const &grid_elem) +template <> void QuokkaSimulation::setInitialConditionsOnGrid(quokka::grid const &grid_elem) { const amrex::GpuArray dx = grid_elem.dx_; const amrex::GpuArray prob_lo = grid_elem.prob_lo_; @@ -111,10 +108,8 @@ void QuokkaSimulation::setInitialConditionsOnGrid(quokka::gr } template <> -void QuokkaSimulation::computeReferenceSolution( - amrex::MultiFab &ref, - amrex::GpuArray const &dx, - amrex::GpuArray const &prob_lo) +void QuokkaSimulation::computeReferenceSolution(amrex::MultiFab &ref, amrex::GpuArray const &dx, + amrex::GpuArray const &prob_lo) { for (amrex::MFIter iter(ref); iter.isValid(); ++iter) { const amrex::Box &indexRange = iter.validbox(); From 568e193855dba299dc97111745057eff08491627 Mon Sep 17 00:00:00 2001 From: neco kriel Date: Fri, 3 Oct 2025 16:07:35 +1000 Subject: [PATCH 03/34] add hydro test --- src/problems/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/src/problems/CMakeLists.txt b/src/problems/CMakeLists.txt index 6338f80f9d..51832a4d2f 100644 --- a/src/problems/CMakeLists.txt +++ b/src/problems/CMakeLists.txt @@ -14,6 +14,7 @@ add_subdirectory(ParticleSF) add_subdirectory(HydroBlast2D) add_subdirectory(HydroBlast3D) add_subdirectory(HydroContact) +add_subdirectory(HydroBalsaraVortex) add_subdirectory(HydroKelvinHelmholz) add_subdirectory(HydroLeblanc) add_subdirectory(HydroRichtmeyerMeshkov) From a1bbe22c9f0f8937d69baa765eda41778b06ad39 Mon Sep 17 00:00:00 2001 From: neco kriel Date: Fri, 3 Oct 2025 16:58:11 +1000 Subject: [PATCH 04/34] run on same domain as Balsara, at higher Nres + for longer --- inputs/hydro_balsara_vortex.in | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/inputs/hydro_balsara_vortex.in b/inputs/hydro_balsara_vortex.in index f29981bc96..83d21437e4 100644 --- a/inputs/hydro_balsara_vortex.in +++ b/inputs/hydro_balsara_vortex.in @@ -1,9 +1,9 @@ # ***************************************************************** # 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 +geometry.prob_lo = -5.0 -5.0 -5.0 +geometry.prob_hi = 5.0 5.0 5.0 +geometry.is_periodic = 1 1 1 # ***************************************************************** # VERBOSITY @@ -13,21 +13,21 @@ amr.v = 1 # ***************************************************************** # Resolution and refinement # ***************************************************************** -amr.n_cell = 32 32 16 +amr.n_cell = 64 64 16 amr.max_level = 0 -amr.max_grid_size = 32 -amr.blocking_factor_x = 32 -amr.blocking_factor_y = 16 +amr.max_grid_size = 64 +amr.blocking_factor_x = 64 +amr.blocking_factor_y = 64 amr.blocking_factor_z = 16 do_reflux = 0 do_subcycle = 0 -plotfile_interval = -1 +plotfile_interval = 10 plotfile_prefix = hydro_balsara_vortex_plt checkpoint_interval = -1 cfl = 0.3 -stop_time = 2.0 +stop_time = 10.0 max_timesteps = 10000 hydro.rk_integrator_order = 2 From d6b09fb40decb87a10f980e96eb4a65e58ee3442 Mon Sep 17 00:00:00 2001 From: neco kriel Date: Fri, 3 Oct 2025 16:58:24 +1000 Subject: [PATCH 05/34] center problem --- .../test_hydro_balsara_vortex.cpp | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/src/problems/HydroBalsaraVortex/test_hydro_balsara_vortex.cpp b/src/problems/HydroBalsaraVortex/test_hydro_balsara_vortex.cpp index e5fd0b6beb..222177fde9 100644 --- a/src/problems/HydroBalsaraVortex/test_hydro_balsara_vortex.cpp +++ b/src/problems/HydroBalsaraVortex/test_hydro_balsara_vortex.cpp @@ -8,7 +8,6 @@ #include #include -#include #include "AMReX_Array.H" #include "AMReX_Array4.H" @@ -46,8 +45,8 @@ constexpr double bg_density = 1.0; constexpr double bg_pressure = 1.0; // vortex parameters constexpr double vortex_speed = 5.0 / (2.0 * M_PI); -constexpr double vortex_x0 = 0.5; -constexpr double vortex_y0 = 0.5; +constexpr double vortex_center_x1 = 0.0; +constexpr double vortex_center_x2 = 0.0; // drift is off by default constexpr double vortex_drift_x1 = 0.0; constexpr double vortex_drift_x2 = 0.0; @@ -62,15 +61,17 @@ inline void computeVortexSolution(int i, int j, int k, amrex::Array4(x1_C) - vortex_x0; - const double rel_x2 = static_cast(x2_C) - vortex_y0; - const double radius_sq = rel_x1 * rel_x1 + rel_x2 * rel_x2; + const double delta_x1_from_center = static_cast(x1_C) - vortex_center_x1; + const double delta_x2_from_center = static_cast(x2_C) - vortex_center_x2; + const double radius_sq = delta_x1_from_center * delta_x1_from_center + delta_x2_from_center * delta_x2_from_center; + const double radial_profile = std::exp(0.5 * (1.0 - radius_sq)); + const double radial_profile_sq = std::exp(1.0 - radius_sq); const double density = bg_density; - const double pressure = bg_pressure - 0.5 * vortex_speed * vortex_speed * gcem::exp(1.0 - radius_sq); + const double pressure = bg_pressure - 0.5 * density * vortex_speed * vortex_speed * radial_profile_sq; - const double delta_vel_x1 = -rel_x2 * vortex_speed * gcem::exp(0.5 * (1.0 - radius_sq)); - const double delta_vel_x2 = rel_x1 * vortex_speed * gcem::exp(0.5 * (1.0 - radius_sq)); + const double delta_vel_x1 = -delta_x2_from_center * vortex_speed * radial_profile; + const double delta_vel_x2 = delta_x1_from_center * vortex_speed * radial_profile; const double vel_x1 = vortex_drift_x1 + delta_vel_x1; const double vel_x2 = vortex_drift_x2 + delta_vel_x2; const double vel_x3 = vortex_drift_x3; From 8b1e3fe63daf1d5ba67f694b9d66c1925979ca80 Mon Sep 17 00:00:00 2001 From: neco kriel Date: Thu, 9 Oct 2025 09:27:59 +1100 Subject: [PATCH 06/34] centralise input file --- inputs/{hydro_balsara_vortex.in => balsara_vortex.in} | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) rename inputs/{hydro_balsara_vortex.in => balsara_vortex.in} (90%) diff --git a/inputs/hydro_balsara_vortex.in b/inputs/balsara_vortex.in similarity index 90% rename from inputs/hydro_balsara_vortex.in rename to inputs/balsara_vortex.in index 83d21437e4..46e8fd0bbe 100644 --- a/inputs/hydro_balsara_vortex.in +++ b/inputs/balsara_vortex.in @@ -24,12 +24,15 @@ do_reflux = 0 do_subcycle = 0 plotfile_interval = 10 -plotfile_prefix = hydro_balsara_vortex_plt +plotfile_prefix = balsara_vortex_plt checkpoint_interval = -1 cfl = 0.3 -stop_time = 10.0 max_timesteps = 10000 hydro.rk_integrator_order = 2 hydro.reconstruction_order = 5 hydro.use_dual_energy = 0 + +setup.vortex_Mach = 0.1 +setup.advection = 1 +setup.num_orbits = 3 From 331bcf88eb7f3e763e2550f1e4f92d71e0045c5f Mon Sep 17 00:00:00 2001 From: neco kriel Date: Thu, 9 Oct 2025 09:28:21 +1100 Subject: [PATCH 07/34] use const arrays --- src/hydro/HLLD.hpp | 47 +++++++++++++++++++++++----------------------- 1 file changed, 24 insertions(+), 23 deletions(-) diff --git a/src/hydro/HLLD.hpp b/src/hydro/HLLD.hpp index 5e804f551a..c5000437d9 100644 --- a/src/hydro/HLLD.hpp +++ b/src/hydro/HLLD.hpp @@ -40,16 +40,16 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState u_star_R{}; // frequently used term - double const bx_sq = SQUARE(bx); + const double bx_sq = SQUARE(bx); // compute L/R states for select conserved variables // (group transverse vector components for floating-point associativity symmetry) // magnetic pressure - double const pb_L = 0.5 * (bx_sq + (SQUARE(sL.by) + SQUARE(sL.bz))); - double const pb_R = 0.5 * (bx_sq + (SQUARE(sR.by) + SQUARE(sR.bz))); + const double pb_L = 0.5 * (bx_sq + (SQUARE(sL.by) + SQUARE(sL.bz))); + const double pb_R = 0.5 * (bx_sq + (SQUARE(sR.by) + SQUARE(sR.bz))); // kinetic energy - double const ke_L = 0.5 * sL.rho * (SQUARE(sL.u) + (SQUARE(sL.v) + SQUARE(sL.w))); - double const ke_R = 0.5 * sR.rho * (SQUARE(sR.u) + (SQUARE(sR.v) + SQUARE(sR.w))); + const double ke_L = 0.5 * sL.rho * (SQUARE(sL.u) + (SQUARE(sL.v) + SQUARE(sL.w))); + const double ke_R = 0.5 * sR.rho * (SQUARE(sR.u) + (SQUARE(sR.v) + SQUARE(sR.w))); // set left conserved states u_L.rho = sL.rho; u_L.mx = sL.u * u_L.rho; @@ -86,8 +86,8 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState compression @@ -131,10 +131,11 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState 0.0 ? 1.0 : -1.0); + const double rho_sum_inv = 1.0 / (rho_sqrt_L + rho_sqrt_R); + const double bx_sign = (bx > 0.0 ? 1.0 : -1.0); u_dstar_L.rho = u_star_L.rho; u_dstar_R.rho = u_star_R.rho; u_dstar_L.mx = u_star_L.mx; From 6c53409943068e04fee38c7d819d74f5a452af5b Mon Sep 17 00:00:00 2001 From: neco kriel Date: Thu, 9 Oct 2025 09:29:38 +1100 Subject: [PATCH 08/34] setup from input file --- .../test_hydro_balsara_vortex.cpp | 52 ++++++++++++++++--- 1 file changed, 44 insertions(+), 8 deletions(-) diff --git a/src/problems/HydroBalsaraVortex/test_hydro_balsara_vortex.cpp b/src/problems/HydroBalsaraVortex/test_hydro_balsara_vortex.cpp index 222177fde9..f91d4c1eb1 100644 --- a/src/problems/HydroBalsaraVortex/test_hydro_balsara_vortex.cpp +++ b/src/problems/HydroBalsaraVortex/test_hydro_balsara_vortex.cpp @@ -2,15 +2,19 @@ // Copyright 2025 Neco Kriel. // Released under the MIT license. See LICENSE file included in the GitHub repo. //============================================================================== -/// \file test_balsara_vortex_hllc.cpp -/// \brief Hydrodynamic (HLLC) Balsara vortex (B = 0; constant-ρ, prescribed p). +/// \file test_hydro_balsara_vortex.cpp +/// \brief hydro (HLLC) Balsara vortex (constant-ρ, prescribed p). /// #include #include +#include +#include #include "AMReX_Array.H" #include "AMReX_Array4.H" +#include "AMReX_Gpu.H" +#include "AMReX_ParmParse.H" #include "AMReX_REAL.H" #include "QuokkaSimulation.hpp" @@ -43,14 +47,15 @@ template <> struct Physics_Traits { constexpr double gamma_gas = quokka::EOS_Traits::gamma; constexpr double bg_density = 1.0; constexpr double bg_pressure = 1.0; +constexpr double sound_speed = gcem::sqrt(gamma_gas * bg_pressure / bg_density); // vortex parameters -constexpr double vortex_speed = 5.0 / (2.0 * M_PI); +AMREX_GPU_MANAGED double vortex_Mach = 0.01; // NOLINT +// domain extends over [-5, 5] by default constexpr double vortex_center_x1 = 0.0; constexpr double vortex_center_x2 = 0.0; // drift is off by default -constexpr double vortex_drift_x1 = 0.0; -constexpr double vortex_drift_x2 = 0.0; -constexpr double vortex_drift_x3 = 0.0; +AMREX_GPU_MANAGED double vortex_drift_x1 = 0.0; // NOLINT +AMREX_GPU_MANAGED double vortex_drift_x2 = 0.0; // NOLINT AMREX_GPU_DEVICE inline void computeVortexSolution(int i, int j, int k, amrex::Array4 const &state, amrex::GpuArray const &dx, @@ -65,16 +70,17 @@ inline void computeVortexSolution(int i, int j, int k, amrex::Array4(x2_C) - vortex_center_x2; const double radius_sq = delta_x1_from_center * delta_x1_from_center + delta_x2_from_center * delta_x2_from_center; const double radial_profile = std::exp(0.5 * (1.0 - radius_sq)); - const double radial_profile_sq = std::exp(1.0 - radius_sq); + const double radial_profile_sq = radial_profile * radial_profile; const double density = bg_density; + const double vortex_speed = vortex_Mach * sound_speed; const double pressure = bg_pressure - 0.5 * density * vortex_speed * vortex_speed * radial_profile_sq; const double delta_vel_x1 = -delta_x2_from_center * vortex_speed * radial_profile; const double delta_vel_x2 = delta_x1_from_center * vortex_speed * radial_profile; const double vel_x1 = vortex_drift_x1 + delta_vel_x1; const double vel_x2 = vortex_drift_x2 + delta_vel_x2; - const double vel_x3 = vortex_drift_x3; + const double vel_x3 = 0.0; const double mom_x1 = density * vel_x1; const double mom_x2 = density * vel_x2; @@ -128,9 +134,39 @@ void QuokkaSimulation::computeReferenceSolution(amrex::Multi auto problem_main() -> int { + amrex::ParmParse const hpp("setup"); + + int advection_int = 0; + int num_orbits = 1; + hpp.query("vortex_Mach", vortex_Mach); + hpp.query("advection", advection_int); + hpp.query("num_orbits", num_orbits); + const double vortex_speed = vortex_Mach * sound_speed; + const bool is_advection_enabled = (advection_int != 0); + auto BCs_cc = quokka::BC(quokka::BCType::int_dir); // periodic QuokkaSimulation sim(BCs_cc); + + double stop_time = 0.0; + if (is_advection_enabled) { + const double advection_speed = vortex_speed; + vortex_drift_x2 = vortex_drift_x1 = advection_speed / std::sqrt(2.0); + const double length_x1 = sim.geom[0].ProbLength(0); + const double length_x2 = sim.geom[0].ProbLength(1); + if (std::abs(length_x1 - length_x2) > 1e-12) { + amrex::Abort("The domain must be square for advection."); + } + const double advection_distance = std::sqrt(length_x1 * length_x1 + length_x2 * length_x2); + const double advection_duration = advection_distance / vortex_speed; + stop_time = static_cast(num_orbits) * advection_duration; + } else { + vortex_drift_x1 = vortex_drift_x2 = 0.0; + const double orbital_duration = 2 * M_PI / vortex_speed; + stop_time = static_cast(num_orbits) * orbital_duration; + } + + sim.stopTime_ = stop_time; sim.computeReferenceSolution_ = true; sim.setInitialConditions(); From fd4f04f417b24b96f27de5714b668cba4f999136 Mon Sep 17 00:00:00 2001 From: neco kriel Date: Thu, 9 Oct 2025 09:32:52 +1100 Subject: [PATCH 09/34] Merge development branch into feature branch --- src/problems/CMakeLists.txt | 1 + src/problems/HydroBalsaraVortex/CMakeLists.txt | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/problems/CMakeLists.txt b/src/problems/CMakeLists.txt index 51832a4d2f..b226a14c06 100644 --- a/src/problems/CMakeLists.txt +++ b/src/problems/CMakeLists.txt @@ -66,6 +66,7 @@ add_subdirectory(RadhydroBB) add_subdirectory(AlfvenWaveLinear) add_subdirectory(AlfvenWaveCircular) +add_subdirectory(MHDBalsaraVortex) add_subdirectory(CurrentSheet) add_subdirectory(OrszagTang) add_subdirectory(FastWave) diff --git a/src/problems/HydroBalsaraVortex/CMakeLists.txt b/src/problems/HydroBalsaraVortex/CMakeLists.txt index 009ff341f0..5f796d2dfc 100644 --- a/src/problems/HydroBalsaraVortex/CMakeLists.txt +++ b/src/problems/HydroBalsaraVortex/CMakeLists.txt @@ -4,5 +4,5 @@ if (AMReX_SPACEDIM EQUAL 3) setup_target_for_cuda_compilation(test_hydro_balsara_vortex) endif() - add_test(NAME HydroBalsaraVortex COMMAND test_hydro_balsara_vortex ../inputs/hydro_balsara_vortex.in ${QuokkaTestParams} WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) + add_test(NAME HydroBalsaraVortex COMMAND test_hydro_balsara_vortex ../inputs/balsara_vortex.in ${QuokkaTestParams} WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) endif() From da5d0206c9af92f7018ce04cab6d7818c1640190 Mon Sep 17 00:00:00 2001 From: neco kriel Date: Thu, 9 Oct 2025 09:34:22 +1100 Subject: [PATCH 10/34] merge development branch into feature branch --- src/problems/MHDBalsaraVortex/CMakeLists.txt | 8 + .../test_mhd_balsara_vortex.cpp | 204 ++++++++++++++++++ 2 files changed, 212 insertions(+) create mode 100644 src/problems/MHDBalsaraVortex/CMakeLists.txt create mode 100644 src/problems/MHDBalsaraVortex/test_mhd_balsara_vortex.cpp diff --git a/src/problems/MHDBalsaraVortex/CMakeLists.txt b/src/problems/MHDBalsaraVortex/CMakeLists.txt new file mode 100644 index 0000000000..6fdc38c8b6 --- /dev/null +++ b/src/problems/MHDBalsaraVortex/CMakeLists.txt @@ -0,0 +1,8 @@ +if (AMReX_SPACEDIM EQUAL 3) + add_executable(test_mhd_balsara_vortex test_mhd_balsara_vortex.cpp ${QuokkaObjSources}) + if(AMReX_GPU_BACKEND MATCHES "CUDA") + setup_target_for_cuda_compilation(test_mhd_balsara_vortex) + endif() + + add_test(NAME HydroBalsaraVortex COMMAND test_mhd_balsara_vortex ../inputs/balsara_vortex.in ${QuokkaTestParams} WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) +endif() diff --git a/src/problems/MHDBalsaraVortex/test_mhd_balsara_vortex.cpp b/src/problems/MHDBalsaraVortex/test_mhd_balsara_vortex.cpp new file mode 100644 index 0000000000..514cde24a7 --- /dev/null +++ b/src/problems/MHDBalsaraVortex/test_mhd_balsara_vortex.cpp @@ -0,0 +1,204 @@ +//============================================================================== +// Copyright 2025 Neco Kriel. +// Released under the MIT license. See LICENSE file included in the GitHub repo. +//============================================================================== +/// \file test_hydro_balsara_vortex.cpp +/// \brief hydro (HLLC) Balsara vortex (constant-ρ, prescribed p). +/// + +#include +#include +#include + +#include "AMReX_Array.H" +#include "AMReX_Array4.H" +#include "AMReX_Gpu.H" +#include "AMReX_ParmParse.H" +#include "AMReX_REAL.H" + +#include "QuokkaSimulation.hpp" +#include "grid.hpp" +#include "hydro/EOS.hpp" +#include "physics_info.hpp" +#include "util/BC.hpp" + +struct MHDBalsaraVortex { +}; + +template <> struct quokka::EOS_Traits { + static constexpr double gamma = 5.0 / 3.0; + static constexpr double mean_molecular_weight = C::m_u; + static constexpr double boltzmann_constant = C::k_B; +}; + +template <> struct Physics_Traits { + static constexpr bool is_hydro_enabled = true; + static constexpr int numMassScalars = 0; + static constexpr int numPassiveScalars = numMassScalars + 0; + static constexpr bool is_self_gravity_enabled = false; + static constexpr bool is_radiation_enabled = false; + static constexpr bool is_mhd_enabled = true; + static constexpr int nGroups = 1; + static constexpr UnitSystem unit_system = UnitSystem::CGS; +}; + +// background state +constexpr double gamma_gas = quokka::EOS_Traits::gamma; +constexpr double bg_density = 1.0; +constexpr double bg_pressure = 1.0; +constexpr double sound_speed = gcem::sqrt(gamma_gas * bg_pressure / bg_density); +// vortex parameters +AMREX_GPU_MANAGED double vortex_Mach = 0.01; // NOLINT +// domain extends over [-5, 5] by default +constexpr double vortex_center_x1 = 0.0; +constexpr double vortex_center_x2 = 0.0; +// drift is off by default +AMREX_GPU_MANAGED double vortex_drift_x1 = 0.0; // NOLINT +AMREX_GPU_MANAGED double vortex_drift_x2 = 0.0; // NOLINT + +AMREX_GPU_DEVICE +inline void computeVortexSolution(int i, int j, int k, amrex::Array4 const &state, amrex::GpuArray const &dx, + amrex::GpuArray const &prob_lo) +{ + const amrex::Real x1_L = prob_lo[0] + i * dx[0]; + const amrex::Real x2_L = prob_lo[1] + j * dx[1]; + const amrex::Real x1_C = x1_L + 0.5 * dx[0]; + const amrex::Real x2_C = x2_L + 0.5 * dx[1]; + + const double delta_x1_from_center = static_cast(x1_C) - vortex_center_x1; + const double delta_x2_from_center = static_cast(x2_C) - vortex_center_x2; + const double radius_sq = delta_x1_from_center * delta_x1_from_center + delta_x2_from_center * delta_x2_from_center; + const double radial_profile = std::exp(0.5 * (1.0 - radius_sq)); + const double radial_profile_sq = radial_profile * radial_profile; + + const double density = bg_density; + const double vortex_speed = vortex_Mach * sound_speed; + const double pressure = bg_pressure - 0.5 * density * vortex_speed * vortex_speed * radial_profile_sq; + + const double delta_vel_x1 = -delta_x2_from_center * vortex_speed * radial_profile; + const double delta_vel_x2 = delta_x1_from_center * vortex_speed * radial_profile; + const double vel_x1 = vortex_drift_x1 + delta_vel_x1; + const double vel_x2 = vortex_drift_x2 + delta_vel_x2; + const double vel_x3 = 0.0; + + const double mom_x1 = density * vel_x1; + const double mom_x2 = density * vel_x2; + const double mom_x3 = density * vel_x3; + + const double Eint = pressure / (gamma_gas - 1.0); + const double Ekin = 0.5 * density * (vel_x1 * vel_x1 + vel_x2 * vel_x2 + vel_x3 * vel_x3); + const double Etot = Eint + Ekin; + + state(i, j, k, HydroSystem::density_index) = density; + state(i, j, k, HydroSystem::x1Momentum_index) = mom_x1; + state(i, j, k, HydroSystem::x2Momentum_index) = mom_x2; + state(i, j, k, HydroSystem::x3Momentum_index) = mom_x3; + state(i, j, k, HydroSystem::energy_index) = Etot; + state(i, j, k, HydroSystem::internalEnergy_index) = Eint; +} + +template <> void QuokkaSimulation::setInitialConditionsOnGrid(quokka::grid const &grid_elem) +{ + const amrex::GpuArray dx = grid_elem.dx_; + const amrex::GpuArray prob_lo = grid_elem.prob_lo_; + const amrex::Array4 &state_cc = grid_elem.array_; + const amrex::Box &indexRange = grid_elem.indexRange_; + const int ncomp_cc = Physics_Indices::nvarTotal_cc; + + amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + for (int icomp = 0; icomp < ncomp_cc; ++icomp) { + state_cc(i, j, k, icomp) = 0.0; + } + computeVortexSolution(i, j, k, state_cc, dx, prob_lo); + }); +} + +template <> void QuokkaSimulation::setInitialConditionsOnGridFaceVars(quokka::grid const &grid_elem) +{ + // extract grid information + const amrex::Array4 &state_fc = grid_elem.array_; + const amrex::Box &indexRange = grid_elem.indexRange_; + + const int ncomp_fc = Physics_Indices::nvarPerDim_fc; + // loop over the grid and set the initial condition + amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + for (int n = 0; n < ncomp_fc; ++n) { + state_fc(i, j, k, n) = 0; // fill all b-field quantities with zeros + } + }); +} + +template <> +void QuokkaSimulation::computeReferenceSolution(amrex::MultiFab &ref, amrex::GpuArray const &dx, + amrex::GpuArray const &prob_lo) +{ + 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, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + for (int icomp = 0; icomp < ncomp; ++icomp) { + stateExact(i, j, k, icomp) = 0.0; + } + computeVortexSolution(i, j, k, stateExact, dx, prob_lo); + }); + } +} + +auto problem_main() -> int +{ + amrex::ParmParse const hpp("setup"); + + int advection_int = 0; + int num_orbits = 1; + hpp.query("vortex_Mach", vortex_Mach); + hpp.query("advection", advection_int); + hpp.query("num_orbits", num_orbits); + const double vortex_speed = vortex_Mach * sound_speed; + const bool is_advection_enabled = (advection_int != 0); + + auto BCs_cc = quokka::BC(quokka::BCType::int_dir); + + const int nvars_fc = Physics_Indices::nvarTotal_fc; + amrex::Vector BCs_fc(nvars_fc); + for (int icomp = 0; icomp < nvars_fc; ++icomp) { + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + BCs_fc[icomp].setLo(idim, amrex::BCType::int_dir); + BCs_fc[icomp].setHi(idim, amrex::BCType::int_dir); + } + } + + QuokkaSimulation sim(BCs_cc, BCs_fc); + + double stop_time = 0.0; + if (is_advection_enabled) { + const double advection_speed = vortex_speed; + vortex_drift_x2 = vortex_drift_x1 = advection_speed / std::sqrt(2.0); + const double length_x1 = sim.geom[0].ProbLength(0); + const double length_x2 = sim.geom[0].ProbLength(1); + if (std::abs(length_x1 - length_x2) > 1e-12) { + amrex::Abort("The domain must be square for advection."); + } + const double advection_distance = std::sqrt(length_x1 * length_x1 + length_x2 * length_x2); + const double advection_duration = advection_distance / vortex_speed; + stop_time = static_cast(num_orbits) * advection_duration; + } else { + vortex_drift_x1 = vortex_drift_x2 = 0.0; + const double orbital_duration = 2 * M_PI / vortex_speed; + stop_time = static_cast(num_orbits) * orbital_duration; + } + + sim.stopTime_ = stop_time; + sim.computeReferenceSolution_ = true; + + sim.setInitialConditions(); + sim.evolve(); + + int status = 0; + const double error_tol = 0.002; + if (sim.errorNorm_ > error_tol) { + status = 1; + } + return status; +} From b55e9335ca11632c22f5d47e8ea708c149680341 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 8 Oct 2025 22:34:47 +0000 Subject: [PATCH 11/34] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/hydro/HLLD.hpp | 2 +- .../HydroBalsaraVortex/test_hydro_balsara_vortex.cpp | 4 ++-- src/problems/MHDBalsaraVortex/test_mhd_balsara_vortex.cpp | 6 +++--- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/hydro/HLLD.hpp b/src/hydro/HLLD.hpp index c5000437d9..b58f544bb0 100644 --- a/src/hydro/HLLD.hpp +++ b/src/hydro/HLLD.hpp @@ -132,7 +132,7 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState::computeReferenceSolution(amrex::Multi auto problem_main() -> int { amrex::ParmParse const hpp("setup"); - + int advection_int = 0; int num_orbits = 1; hpp.query("vortex_Mach", vortex_Mach); @@ -147,7 +147,7 @@ auto problem_main() -> int auto BCs_cc = quokka::BC(quokka::BCType::int_dir); // periodic QuokkaSimulation sim(BCs_cc); - + double stop_time = 0.0; if (is_advection_enabled) { const double advection_speed = vortex_speed; diff --git a/src/problems/MHDBalsaraVortex/test_mhd_balsara_vortex.cpp b/src/problems/MHDBalsaraVortex/test_mhd_balsara_vortex.cpp index 514cde24a7..058e639e77 100644 --- a/src/problems/MHDBalsaraVortex/test_mhd_balsara_vortex.cpp +++ b/src/problems/MHDBalsaraVortex/test_mhd_balsara_vortex.cpp @@ -130,7 +130,7 @@ template <> void QuokkaSimulation::setInitialConditionsOnGridF template <> void QuokkaSimulation::computeReferenceSolution(amrex::MultiFab &ref, amrex::GpuArray const &dx, - amrex::GpuArray const &prob_lo) + amrex::GpuArray const &prob_lo) { for (amrex::MFIter iter(ref); iter.isValid(); ++iter) { const amrex::Box &indexRange = iter.validbox(); @@ -149,7 +149,7 @@ void QuokkaSimulation::computeReferenceSolution(amrex::MultiFa auto problem_main() -> int { amrex::ParmParse const hpp("setup"); - + int advection_int = 0; int num_orbits = 1; hpp.query("vortex_Mach", vortex_Mach); @@ -170,7 +170,7 @@ auto problem_main() -> int } QuokkaSimulation sim(BCs_cc, BCs_fc); - + double stop_time = 0.0; if (is_advection_enabled) { const double advection_speed = vortex_speed; From 8b6c3c06f7b2dcf90d95cbac7241655465206426 Mon Sep 17 00:00:00 2001 From: neco kriel Date: Tue, 23 Dec 2025 10:11:26 +1100 Subject: [PATCH 12/34] remove hydro test + rename tests to be inline with new standard --- ...{balsara_vortex.in => MHDBalsaraVortex.in} | 0 src/problems/CMakeLists.txt | 1 - .../HydroBalsaraVortex/CMakeLists.txt | 8 - .../test_hydro_balsara_vortex.cpp | 181 ------------------ src/problems/MHDBalsaraVortex/CMakeLists.txt | 20 +- ...alsara_vortex.cpp => testMHDBalsaraVortex} | 2 +- 6 files changed, 15 insertions(+), 197 deletions(-) rename inputs/{balsara_vortex.in => MHDBalsaraVortex.in} (100%) delete mode 100644 src/problems/HydroBalsaraVortex/CMakeLists.txt delete mode 100644 src/problems/HydroBalsaraVortex/test_hydro_balsara_vortex.cpp rename src/problems/MHDBalsaraVortex/{test_mhd_balsara_vortex.cpp => testMHDBalsaraVortex} (99%) diff --git a/inputs/balsara_vortex.in b/inputs/MHDBalsaraVortex.in similarity index 100% rename from inputs/balsara_vortex.in rename to inputs/MHDBalsaraVortex.in diff --git a/src/problems/CMakeLists.txt b/src/problems/CMakeLists.txt index bee63dd61c..563bd62449 100644 --- a/src/problems/CMakeLists.txt +++ b/src/problems/CMakeLists.txt @@ -17,7 +17,6 @@ add_subdirectory(ParticleSF) add_subdirectory(HydroBlast2D) add_subdirectory(HydroBlast3D) add_subdirectory(HydroContact) -add_subdirectory(HydroBalsaraVortex) add_subdirectory(HydroKelvinHelmholz) add_subdirectory(HydroLeblanc) add_subdirectory(HydroRichtmeyerMeshkov) diff --git a/src/problems/HydroBalsaraVortex/CMakeLists.txt b/src/problems/HydroBalsaraVortex/CMakeLists.txt deleted file mode 100644 index 5f796d2dfc..0000000000 --- a/src/problems/HydroBalsaraVortex/CMakeLists.txt +++ /dev/null @@ -1,8 +0,0 @@ -if (AMReX_SPACEDIM EQUAL 3) - add_executable(test_hydro_balsara_vortex test_hydro_balsara_vortex.cpp ${QuokkaObjSources}) - if(AMReX_GPU_BACKEND MATCHES "CUDA") - setup_target_for_cuda_compilation(test_hydro_balsara_vortex) - endif() - - add_test(NAME HydroBalsaraVortex COMMAND test_hydro_balsara_vortex ../inputs/balsara_vortex.in ${QuokkaTestParams} WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) -endif() diff --git a/src/problems/HydroBalsaraVortex/test_hydro_balsara_vortex.cpp b/src/problems/HydroBalsaraVortex/test_hydro_balsara_vortex.cpp deleted file mode 100644 index cf30e3f749..0000000000 --- a/src/problems/HydroBalsaraVortex/test_hydro_balsara_vortex.cpp +++ /dev/null @@ -1,181 +0,0 @@ -//============================================================================== -// Copyright 2025 Neco Kriel. -// Released under the MIT license. See LICENSE file included in the GitHub repo. -//============================================================================== -/// \file test_hydro_balsara_vortex.cpp -/// \brief hydro (HLLC) Balsara vortex (constant-ρ, prescribed p). -/// - -#include -#include -#include -#include - -#include "AMReX_Array.H" -#include "AMReX_Array4.H" -#include "AMReX_Gpu.H" -#include "AMReX_ParmParse.H" -#include "AMReX_REAL.H" - -#include "QuokkaSimulation.hpp" -#include "grid.hpp" -#include "hydro/EOS.hpp" -#include "physics_info.hpp" -#include "util/BC.hpp" - -struct HydroBalsaraVortex { -}; - -template <> struct quokka::EOS_Traits { - static constexpr double gamma = 5.0 / 3.0; - static constexpr double mean_molecular_weight = C::m_u; - static constexpr double boltzmann_constant = C::k_B; -}; - -template <> struct Physics_Traits { - static constexpr bool is_hydro_enabled = true; - static constexpr int numMassScalars = 0; - static constexpr int numPassiveScalars = numMassScalars + 0; - static constexpr bool is_self_gravity_enabled = false; - static constexpr bool is_radiation_enabled = false; - static constexpr bool is_mhd_enabled = false; - static constexpr int nGroups = 1; - static constexpr UnitSystem unit_system = UnitSystem::CGS; -}; - -// background state -constexpr double gamma_gas = quokka::EOS_Traits::gamma; -constexpr double bg_density = 1.0; -constexpr double bg_pressure = 1.0; -constexpr double sound_speed = gcem::sqrt(gamma_gas * bg_pressure / bg_density); -// vortex parameters -AMREX_GPU_MANAGED double vortex_Mach = 0.01; // NOLINT -// domain extends over [-5, 5] by default -constexpr double vortex_center_x1 = 0.0; -constexpr double vortex_center_x2 = 0.0; -// drift is off by default -AMREX_GPU_MANAGED double vortex_drift_x1 = 0.0; // NOLINT -AMREX_GPU_MANAGED double vortex_drift_x2 = 0.0; // NOLINT - -AMREX_GPU_DEVICE -inline void computeVortexSolution(int i, int j, int k, amrex::Array4 const &state, amrex::GpuArray const &dx, - amrex::GpuArray const &prob_lo) -{ - const amrex::Real x1_L = prob_lo[0] + i * dx[0]; - const amrex::Real x2_L = prob_lo[1] + j * dx[1]; - const amrex::Real x1_C = x1_L + 0.5 * dx[0]; - const amrex::Real x2_C = x2_L + 0.5 * dx[1]; - - const double delta_x1_from_center = static_cast(x1_C) - vortex_center_x1; - const double delta_x2_from_center = static_cast(x2_C) - vortex_center_x2; - const double radius_sq = delta_x1_from_center * delta_x1_from_center + delta_x2_from_center * delta_x2_from_center; - const double radial_profile = std::exp(0.5 * (1.0 - radius_sq)); - const double radial_profile_sq = radial_profile * radial_profile; - - const double density = bg_density; - const double vortex_speed = vortex_Mach * sound_speed; - const double pressure = bg_pressure - 0.5 * density * vortex_speed * vortex_speed * radial_profile_sq; - - const double delta_vel_x1 = -delta_x2_from_center * vortex_speed * radial_profile; - const double delta_vel_x2 = delta_x1_from_center * vortex_speed * radial_profile; - const double vel_x1 = vortex_drift_x1 + delta_vel_x1; - const double vel_x2 = vortex_drift_x2 + delta_vel_x2; - const double vel_x3 = 0.0; - - const double mom_x1 = density * vel_x1; - const double mom_x2 = density * vel_x2; - const double mom_x3 = density * vel_x3; - - const double Eint = pressure / (gamma_gas - 1.0); - const double Ekin = 0.5 * density * (vel_x1 * vel_x1 + vel_x2 * vel_x2 + vel_x3 * vel_x3); - const double Etot = Eint + Ekin; - - state(i, j, k, HydroSystem::density_index) = density; - state(i, j, k, HydroSystem::x1Momentum_index) = mom_x1; - state(i, j, k, HydroSystem::x2Momentum_index) = mom_x2; - state(i, j, k, HydroSystem::x3Momentum_index) = mom_x3; - state(i, j, k, HydroSystem::energy_index) = Etot; - state(i, j, k, HydroSystem::internalEnergy_index) = Eint; -} - -template <> void QuokkaSimulation::setInitialConditionsOnGrid(quokka::grid const &grid_elem) -{ - const amrex::GpuArray dx = grid_elem.dx_; - const amrex::GpuArray prob_lo = grid_elem.prob_lo_; - const amrex::Array4 &state_cc = grid_elem.array_; - const amrex::Box &indexRange = grid_elem.indexRange_; - const int ncomp_cc = Physics_Indices::nvarTotal_cc; - - amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { - for (int icomp = 0; icomp < ncomp_cc; ++icomp) { - state_cc(i, j, k, icomp) = 0.0; - } - computeVortexSolution(i, j, k, state_cc, dx, prob_lo); - }); -} - -template <> -void QuokkaSimulation::computeReferenceSolution(amrex::MultiFab &ref, amrex::GpuArray const &dx, - amrex::GpuArray const &prob_lo) -{ - 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, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - for (int icomp = 0; icomp < ncomp; ++icomp) { - stateExact(i, j, k, icomp) = 0.0; - } - computeVortexSolution(i, j, k, stateExact, dx, prob_lo); - }); - } -} - -auto problem_main() -> int -{ - amrex::ParmParse const hpp("setup"); - - int advection_int = 0; - int num_orbits = 1; - hpp.query("vortex_Mach", vortex_Mach); - hpp.query("advection", advection_int); - hpp.query("num_orbits", num_orbits); - const double vortex_speed = vortex_Mach * sound_speed; - const bool is_advection_enabled = (advection_int != 0); - - auto BCs_cc = quokka::BC(quokka::BCType::int_dir); // periodic - - QuokkaSimulation sim(BCs_cc); - - double stop_time = 0.0; - if (is_advection_enabled) { - const double advection_speed = vortex_speed; - vortex_drift_x2 = vortex_drift_x1 = advection_speed / std::sqrt(2.0); - const double length_x1 = sim.geom[0].ProbLength(0); - const double length_x2 = sim.geom[0].ProbLength(1); - if (std::abs(length_x1 - length_x2) > 1e-12) { - amrex::Abort("The domain must be square for advection."); - } - const double advection_distance = std::sqrt(length_x1 * length_x1 + length_x2 * length_x2); - const double advection_duration = advection_distance / vortex_speed; - stop_time = static_cast(num_orbits) * advection_duration; - } else { - vortex_drift_x1 = vortex_drift_x2 = 0.0; - const double orbital_duration = 2 * M_PI / vortex_speed; - stop_time = static_cast(num_orbits) * orbital_duration; - } - - sim.stopTime_ = stop_time; - sim.computeReferenceSolution_ = true; - - sim.setInitialConditions(); - sim.evolve(); - - int status = 0; - const double error_tol = 0.002; - if (sim.errorNorm_ > error_tol) { - status = 1; - } - return status; -} diff --git a/src/problems/MHDBalsaraVortex/CMakeLists.txt b/src/problems/MHDBalsaraVortex/CMakeLists.txt index 6fdc38c8b6..80067f3a8a 100644 --- a/src/problems/MHDBalsaraVortex/CMakeLists.txt +++ b/src/problems/MHDBalsaraVortex/CMakeLists.txt @@ -1,8 +1,16 @@ -if (AMReX_SPACEDIM EQUAL 3) - add_executable(test_mhd_balsara_vortex test_mhd_balsara_vortex.cpp ${QuokkaObjSources}) - if(AMReX_GPU_BACKEND MATCHES "CUDA") - setup_target_for_cuda_compilation(test_mhd_balsara_vortex) - endif() +if(AMReX_SPACEDIM EQUAL 3) + set(JOB_NAME MHDBalsaraVortex) - add_test(NAME HydroBalsaraVortex COMMAND test_mhd_balsara_vortex ../inputs/balsara_vortex.in ${QuokkaTestParams} WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) + add_executable(${JOB_NAME} test${JOB_NAME}.cpp ${QuokkaObjSources}) + + if(AMReX_GPU_BACKEND MATCHES "CUDA") + setup_target_for_cuda_compilation(${JOB_NAME}) + endif() + + add_test( + NAME ${JOB_NAME} + COMMAND ${JOB_NAME} ../inputs/${JOB_NAME}.in ${QuokkaTestParams} + WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) + + set_tests_properties(MHDBalsaraVortex PROPERTIES LABELS "MHD") endif() diff --git a/src/problems/MHDBalsaraVortex/test_mhd_balsara_vortex.cpp b/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex similarity index 99% rename from src/problems/MHDBalsaraVortex/test_mhd_balsara_vortex.cpp rename to src/problems/MHDBalsaraVortex/testMHDBalsaraVortex index 058e639e77..673e18aa0c 100644 --- a/src/problems/MHDBalsaraVortex/test_mhd_balsara_vortex.cpp +++ b/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex @@ -2,7 +2,7 @@ // Copyright 2025 Neco Kriel. // Released under the MIT license. See LICENSE file included in the GitHub repo. //============================================================================== -/// \file test_hydro_balsara_vortex.cpp +/// \file testMHDBalsaraVortex.cpp /// \brief hydro (HLLC) Balsara vortex (constant-ρ, prescribed p). /// From 46fa51b9ed48053b7425e3bc7fa0d8dd1ebd598e Mon Sep 17 00:00:00 2001 From: neco kriel Date: Tue, 23 Dec 2025 13:21:03 +1100 Subject: [PATCH 13/34] draft: extend test to MHD --- .../MHDBalsaraVortex/testMHDBalsaraVortex | 173 ++++++++++++------ 1 file changed, 122 insertions(+), 51 deletions(-) diff --git a/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex b/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex index 673e18aa0c..4e204051a9 100644 --- a/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex +++ b/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex @@ -3,7 +3,7 @@ // Released under the MIT license. See LICENSE file included in the GitHub repo. //============================================================================== /// \file testMHDBalsaraVortex.cpp -/// \brief hydro (HLLC) Balsara vortex (constant-ρ, prescribed p). +/// \brief advecting MHD Balsara vortex. /// #include @@ -49,82 +49,134 @@ constexpr double bg_pressure = 1.0; constexpr double sound_speed = gcem::sqrt(gamma_gas * bg_pressure / bg_density); // vortex parameters AMREX_GPU_MANAGED double vortex_Mach = 0.01; // NOLINT +AMREX_GPU_MANAGED double vortex_b_magn = 0.01; // NOLINT // domain extends over [-5, 5] by default constexpr double vortex_center_x1 = 0.0; constexpr double vortex_center_x2 = 0.0; -// drift is off by default +// no drift by default AMREX_GPU_MANAGED double vortex_drift_x1 = 0.0; // NOLINT AMREX_GPU_MANAGED double vortex_drift_x2 = 0.0; // NOLINT +AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto computeRadiusSq(const double x1, const double x2) -> double +{ + const double delta_x1_from_center = x1 - vortex_center_x1; + const double delta_x2_from_center = x2 - vortex_center_x2; + return delta_x1_from_center * delta_x1_from_center + delta_x2_from_center * delta_x2_from_center; +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto computeRadialProfile(const double radius_sq) -> double +{ + return std::exp(0.5 * (1.0 - radius_sq)); +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto Az(const double x1, const double x2) -> double +{ + // vec(a) = (0, 0, a_z), with: + // a_z = b_magn * exp((1 - r^2)/2) + // so that: + // b_x = d(a_z)/dy = -y * b_magn * exp((1 - r^2)/2) + // b_y = -d(a_z)/dx = x * b_magn * exp((1 - r^2)/2) + const double radius_sq = computeRadiusSq(x1, x2); + const double radial_profile = computeRadialProfile(radius_sq); + return vortex_b_magn * radial_profile; +} + AMREX_GPU_DEVICE -inline void computeVortexSolution(int i, int j, int k, amrex::Array4 const &state, amrex::GpuArray const &dx, - amrex::GpuArray const &prob_lo) +void computeVortexSolution(int i, int j, int k, amrex::Array4 const &state, amrex::GpuArray const &dx, amrex::GpuArray const &prob_lo, quokka::centering cen, quokka::direction dir) { + // lower-left corner const amrex::Real x1_L = prob_lo[0] + i * dx[0]; const amrex::Real x2_L = prob_lo[1] + j * dx[1]; - const amrex::Real x1_C = x1_L + 0.5 * dx[0]; - const amrex::Real x2_C = x2_L + 0.5 * dx[1]; - - const double delta_x1_from_center = static_cast(x1_C) - vortex_center_x1; - const double delta_x2_from_center = static_cast(x2_C) - vortex_center_x2; - const double radius_sq = delta_x1_from_center * delta_x1_from_center + delta_x2_from_center * delta_x2_from_center; - const double radial_profile = std::exp(0.5 * (1.0 - radius_sq)); - const double radial_profile_sq = radial_profile * radial_profile; - - const double density = bg_density; - const double vortex_speed = vortex_Mach * sound_speed; - const double pressure = bg_pressure - 0.5 * density * vortex_speed * vortex_speed * radial_profile_sq; - - const double delta_vel_x1 = -delta_x2_from_center * vortex_speed * radial_profile; - const double delta_vel_x2 = delta_x1_from_center * vortex_speed * radial_profile; - const double vel_x1 = vortex_drift_x1 + delta_vel_x1; - const double vel_x2 = vortex_drift_x2 + delta_vel_x2; - const double vel_x3 = 0.0; - - const double mom_x1 = density * vel_x1; - const double mom_x2 = density * vel_x2; - const double mom_x3 = density * vel_x3; - - const double Eint = pressure / (gamma_gas - 1.0); - const double Ekin = 0.5 * density * (vel_x1 * vel_x1 + vel_x2 * vel_x2 + vel_x3 * vel_x3); - const double Etot = Eint + Ekin; - - state(i, j, k, HydroSystem::density_index) = density; - state(i, j, k, HydroSystem::x1Momentum_index) = mom_x1; - state(i, j, k, HydroSystem::x2Momentum_index) = mom_x2; - state(i, j, k, HydroSystem::x3Momentum_index) = mom_x3; - state(i, j, k, HydroSystem::energy_index) = Etot; - state(i, j, k, HydroSystem::internalEnergy_index) = Eint; + + if (cen == quokka::centering::cc) { + const amrex::Real x1_C = x1_L + static_cast(0.5) * dx[0]; + const amrex::Real x2_C = x2_L + static_cast(0.5) * dx[1]; + + const double radius_sq = computeRadiusSq(x1_C, x2_C); + const double radial_profile = computeRadialProfile(radius_sq); + const double radial_profile_sq = radial_profile * radial_profile; + const double delta_x1_from_center = x1_C - vortex_center_x1; + const double delta_x2_from_center = x2_C - vortex_center_x2; + + const double vortex_u_magn = vortex_Mach * sound_speed; + + const double pressure = bg_pressure + 0.5 * (vortex_b_magn * vortex_b_magn * (1.0 - radius_sq) - vortex_u_magn * vortex_u_magn) * radial_profile_sq; + const double Eint = pressure / (gamma_gas - 1.0); + + const double delta_vel_x1 = -delta_x2_from_center * vortex_u_magn * radial_profile; + const double delta_vel_x2 = delta_x1_from_center * vortex_u_magn * radial_profile; + const double vel_x1 = vortex_drift_x1 + delta_vel_x1; + const double vel_x2 = vortex_drift_x2 + delta_vel_x2; + const double vel_x3 = 0.0; + const double mom_x1 = bg_density * vel_x1; + const double mom_x2 = bg_density * vel_x2; + const double mom_x3 = bg_density * vel_x3; + const double Ekin = 0.5 * bg_density * (vel_x1 * vel_x1 + vel_x2 * vel_x2 + vel_x3 * vel_x3); + + const double Bx_cc = -delta_x2_from_center * vortex_b_magn * radial_profile; + const double By_cc = delta_x1_from_center * vortex_b_magn * radial_profile; + const double Bz_cc = 0.0; + const double Emag = 0.5 * (Bx_cc * Bx_cc + By_cc * By_cc + Bz_cc * Bz_cc); + + const double Etot = Eint + Ekin + Emag; + + state(i, j, k, HydroSystem::density_index) = bg_density; + state(i, j, k, HydroSystem::x1Momentum_index) = mom_x1; + state(i, j, k, HydroSystem::x2Momentum_index) = mom_x2; + state(i, j, k, HydroSystem::x3Momentum_index) = mom_x3; + state(i, j, k, HydroSystem::energy_index) = Etot; + state(i, j, k, HydroSystem::internalEnergy_index) = Eint; + + } else if (cen == quokka::centering::fc) { + const auto b_x1_L = (Az(x1_L, x2_L + static_cast(dx[1])) - Az(x1_L, x2_L)) / static_cast(dx[1]); + const auto b_x2_L = (Az(x1_L, x2_L) - Az(x1_L + static_cast(dx[0]), x2_L)) / static_cast(dx[0]); + const auto b_x3_L = 0.0; + + if (dir == quokka::direction::x) { + state(i, j, k, MHDSystem::bfield_index) = b_x1_L; + } else if (dir == quokka::direction::y) { + state(i, j, k, MHDSystem::bfield_index) = b_x2_L; + } else if (dir == quokka::direction::z) { + state(i, j, k, MHDSystem::bfield_index) = b_x3_L; + } + } } template <> void QuokkaSimulation::setInitialConditionsOnGrid(quokka::grid const &grid_elem) { const amrex::GpuArray dx = grid_elem.dx_; const amrex::GpuArray prob_lo = grid_elem.prob_lo_; - const amrex::Array4 &state_cc = grid_elem.array_; + const amrex::Array4 &state_cc = grid_elem.array_; const amrex::Box &indexRange = grid_elem.indexRange_; + const quokka::centering cen = grid_elem.cen_; + const quokka::direction dir = grid_elem.dir_; + const int ncomp_cc = Physics_Indices::nvarTotal_cc; amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { for (int icomp = 0; icomp < ncomp_cc; ++icomp) { - state_cc(i, j, k, icomp) = 0.0; + state_cc(i, j, k, icomp) = 0.0; // fill unused quantities with zeros } - computeVortexSolution(i, j, k, state_cc, dx, prob_lo); + computeVortexSolution(i, j, k, state_cc, dx, prob_lo, cen, dir); }); } template <> void QuokkaSimulation::setInitialConditionsOnGridFaceVars(quokka::grid const &grid_elem) { - // extract grid information - const amrex::Array4 &state_fc = grid_elem.array_; + const amrex::GpuArray dx = grid_elem.dx_; + const amrex::GpuArray prob_lo = grid_elem.prob_lo_; + const amrex::Array4 &state_fc = grid_elem.array_; const amrex::Box &indexRange = grid_elem.indexRange_; + const quokka::centering cen = grid_elem.cen_; + const quokka::direction dir = grid_elem.dir_; const int ncomp_fc = Physics_Indices::nvarPerDim_fc; - // loop over the grid and set the initial condition + amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { - for (int n = 0; n < ncomp_fc; ++n) { - state_fc(i, j, k, n) = 0; // fill all b-field quantities with zeros + for (int icomp = 0; icomp < ncomp_fc; ++icomp) { + state_fc(i, j, k, icomp) = 0.0; // fill unused quantities with zeros } + computeVortexSolution(i, j, k, state_fc, dx, prob_lo, cen, dir); }); } @@ -139,9 +191,27 @@ void QuokkaSimulation::computeReferenceSolution(amrex::MultiFa amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { for (int icomp = 0; icomp < ncomp; ++icomp) { - stateExact(i, j, k, icomp) = 0.0; + stateExact(i, j, k, icomp) = 0.0; // fill unused quantities with zeros + } + computeVortexSolution(i, j, k, stateExact, dx, prob_lo, quokka::centering::cc, quokka::direction::na); + }); + } +} + +template <> +void QuokkaSimulation::computeReferenceSolution_fc(amrex::MultiFab &ref, amrex::GpuArray const &dx, + amrex::GpuArray const &prob_lo, quokka::direction const dir) +{ + 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, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + for (int icomp = 0; icomp < ncomp; ++icomp) { + stateExact(i, j, k, icomp) = 0.0; // fill unused quantities with zeros } - computeVortexSolution(i, j, k, stateExact, dx, prob_lo); + computeVortexSolution(i, j, k, stateExact, dx, prob_lo, quokka::centering::fc, dir); }); } } @@ -153,9 +223,10 @@ auto problem_main() -> int int advection_int = 0; int num_orbits = 1; hpp.query("vortex_Mach", vortex_Mach); + hpp.query("vortex_b_magn", vortex_b_magn); hpp.query("advection", advection_int); hpp.query("num_orbits", num_orbits); - const double vortex_speed = vortex_Mach * sound_speed; + const double vortex_u_magn = vortex_Mach * sound_speed; const bool is_advection_enabled = (advection_int != 0); auto BCs_cc = quokka::BC(quokka::BCType::int_dir); @@ -173,7 +244,7 @@ auto problem_main() -> int double stop_time = 0.0; if (is_advection_enabled) { - const double advection_speed = vortex_speed; + const double advection_speed = vortex_u_magn; vortex_drift_x2 = vortex_drift_x1 = advection_speed / std::sqrt(2.0); const double length_x1 = sim.geom[0].ProbLength(0); const double length_x2 = sim.geom[0].ProbLength(1); @@ -181,11 +252,11 @@ auto problem_main() -> int amrex::Abort("The domain must be square for advection."); } const double advection_distance = std::sqrt(length_x1 * length_x1 + length_x2 * length_x2); - const double advection_duration = advection_distance / vortex_speed; + const double advection_duration = advection_distance / vortex_u_magn; stop_time = static_cast(num_orbits) * advection_duration; } else { vortex_drift_x1 = vortex_drift_x2 = 0.0; - const double orbital_duration = 2 * M_PI / vortex_speed; + const double orbital_duration = 2.0 * M_PI / vortex_u_magn; stop_time = static_cast(num_orbits) * orbital_duration; } From 80e418d7704fbb773b97fafc7f072da41349a85d Mon Sep 17 00:00:00 2001 From: neco kriel Date: Tue, 23 Dec 2025 15:09:30 +1100 Subject: [PATCH 14/34] draft: correct low Mach pressure scaling --- src/hydro/HLLD.hpp | 40 +++++++++++++++++++++++++--------------- 1 file changed, 25 insertions(+), 15 deletions(-) diff --git a/src/hydro/HLLD.hpp b/src/hydro/HLLD.hpp index 833520c3f6..ed3478aa8d 100644 --- a/src/hydro/HLLD.hpp +++ b/src/hydro/HLLD.hpp @@ -15,7 +15,8 @@ namespace quokka::Riemann { constexpr double DELTA = 1.0e-4; -// HLLD solver following Miyoshi and Kusano (2005), hereafter MK5. +// HLLD solver based on Miyoshi & Kusano (2005), hereafter MK5. +// Corrections based on Minoshima & Miyoshi (2021), hereafter MM21. template AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState const &sL, quokka::HydroState const &sR, const double gamma, const double bx, const double perp_v_jump) @@ -76,10 +77,10 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState compression - double theta = 1.0; // tp := shock anisotropy, clamped to [0, 1], with theta = tp^4 - const double denom_tp = std::max(1e-14, max_spd - std::min(perp_v_jump, 0.0)); - double tp = (max_spd - std::min(para_v_jump, 0.0)) / denom_tp; + const double denom_tp = std::max(1e-14, spd_fms_max - std::min(perp_v_jump, 0.0)); + double tp = (spd_fms_max - std::min(para_v_jump, 0.0)) / denom_tp; tp = amrex::Clamp(tp, 0.0, 1.0); - theta = SQUARE(SQUARE(tp)); - // modified middle speed S_M from MK5 eqn (38) + const double theta = SQUARE(SQUARE(tp)); + // modified middle speed S_M from MK5 eqn 38 with theta from MM21 eqn 9 const double sm_denom = (siui_R * u_R.rho - siui_L * u_L.rho); spds[2] = (siui_R * u_R.mx - siui_L * u_L.mx + theta * (ptot_L - ptot_R)) / sm_denom; // S_i - S_M (for i=L or R) const double sism_L = spds[0] - spds[2]; - const double sism_R = spds[4] - spds[2]; const double sism_inv_L = 1.0 / sism_L; const double sism_inv_R = 1.0 / sism_R; @@ -156,11 +155,22 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState 0 for very low Mach, phi -> 1 as advection becomes comparable to wave speeds) + const double spd_fss_L = std::abs(sL.u) + spd_fms_L; + const double spd_fss_R = std::abs(sR.u) + spd_fms_R; + const double spd_fss_max = std::max(spd_fss_L, spd_fss_R); + const double chi = std::min(1.0, spd_fss_max / std::max(1e-14, spd_fms_max)); + const double phi = chi * (2.0 - chi); + // MK5 star total pressures estimate (eqn 23) const double ptot_star_L = ptot_L + u_L.rho * siui_L * (spds[2] - sL.u); const double ptot_star_R = ptot_R + u_R.rho * siui_R * (spds[2] - sR.u); - const double ptot_star = 0.5 * (ptot_star_L + ptot_star_R); + const double ptot_star_avg_LR = 0.5 * (ptot_star_L + ptot_star_R); + // correct the excess (compressive) contributions in the low Mach regime (MM21) + const double ptot_avg_LR = 0.5 * (ptot_L + ptot_R); + const double ptot_star_excess = ptot_star_avg_LR - ptot_avg_LR; + const double ptot_star = ptot_avg_LR + phi * ptot_star_excess; // MK5: u_L^(star, dstar) from, eqn (39) u_star_L.mx = u_star_L.rho * spds[2]; From f3fb0f013cebb5dcf7c181d73a20f38753543870 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 23 Dec 2025 23:23:20 +0000 Subject: [PATCH 15/34] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/hydro/HLLD.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/hydro/HLLD.hpp b/src/hydro/HLLD.hpp index ec93f3bc41..2dce38406e 100644 --- a/src/hydro/HLLD.hpp +++ b/src/hydro/HLLD.hpp @@ -220,7 +220,6 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState 0.0 ? 1.0 : -1.0); @@ -234,7 +233,7 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState Date: Wed, 24 Dec 2025 10:23:51 +1100 Subject: [PATCH 16/34] minor style fixes post merge --- src/hydro/HLLD.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/hydro/HLLD.hpp b/src/hydro/HLLD.hpp index ec93f3bc41..2dce38406e 100644 --- a/src/hydro/HLLD.hpp +++ b/src/hydro/HLLD.hpp @@ -220,7 +220,6 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState 0.0 ? 1.0 : -1.0); @@ -234,7 +233,7 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState Date: Wed, 24 Dec 2025 10:48:41 +1100 Subject: [PATCH 17/34] fix: add cpp extension --- .../{testMHDBalsaraVortex => testMHDBalsaraVortex.cpp} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/problems/MHDBalsaraVortex/{testMHDBalsaraVortex => testMHDBalsaraVortex.cpp} (100%) diff --git a/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex b/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp similarity index 100% rename from src/problems/MHDBalsaraVortex/testMHDBalsaraVortex rename to src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp From a21e9b3de0cbd09c75252afa006d176e7783f5eb Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 23 Dec 2025 23:48:55 +0000 Subject: [PATCH 18/34] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../MHDBalsaraVortex/testMHDBalsaraVortex.cpp | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp b/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp index 4e204051a9..f9e7e8b3c3 100644 --- a/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp +++ b/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp @@ -48,7 +48,7 @@ constexpr double bg_density = 1.0; constexpr double bg_pressure = 1.0; constexpr double sound_speed = gcem::sqrt(gamma_gas * bg_pressure / bg_density); // vortex parameters -AMREX_GPU_MANAGED double vortex_Mach = 0.01; // NOLINT +AMREX_GPU_MANAGED double vortex_Mach = 0.01; // NOLINT AMREX_GPU_MANAGED double vortex_b_magn = 0.01; // NOLINT // domain extends over [-5, 5] by default constexpr double vortex_center_x1 = 0.0; @@ -64,10 +64,7 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto computeRadiusSq(const double x1, const return delta_x1_from_center * delta_x1_from_center + delta_x2_from_center * delta_x2_from_center; } -AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto computeRadialProfile(const double radius_sq) -> double -{ - return std::exp(0.5 * (1.0 - radius_sq)); -} +AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto computeRadialProfile(const double radius_sq) -> double { return std::exp(0.5 * (1.0 - radius_sq)); } AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto Az(const double x1, const double x2) -> double { @@ -82,7 +79,8 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto Az(const double x1, const double x2) -> } AMREX_GPU_DEVICE -void computeVortexSolution(int i, int j, int k, amrex::Array4 const &state, amrex::GpuArray const &dx, amrex::GpuArray const &prob_lo, quokka::centering cen, quokka::direction dir) +void computeVortexSolution(int i, int j, int k, amrex::Array4 const &state, amrex::GpuArray const &dx, + amrex::GpuArray const &prob_lo, quokka::centering cen, quokka::direction dir) { // lower-left corner const amrex::Real x1_L = prob_lo[0] + i * dx[0]; @@ -99,8 +97,9 @@ void computeVortexSolution(int i, int j, int k, amrex::Array4 const const double delta_x2_from_center = x2_C - vortex_center_x2; const double vortex_u_magn = vortex_Mach * sound_speed; - - const double pressure = bg_pressure + 0.5 * (vortex_b_magn * vortex_b_magn * (1.0 - radius_sq) - vortex_u_magn * vortex_u_magn) * radial_profile_sq; + + const double pressure = + bg_pressure + 0.5 * (vortex_b_magn * vortex_b_magn * (1.0 - radius_sq) - vortex_u_magn * vortex_u_magn) * radial_profile_sq; const double Eint = pressure / (gamma_gas - 1.0); const double delta_vel_x1 = -delta_x2_from_center * vortex_u_magn * radial_profile; From a54108958a904e48aa84211e7dbddb84f8947c33 Mon Sep 17 00:00:00 2001 From: neco kriel Date: Wed, 24 Dec 2025 11:08:10 +1100 Subject: [PATCH 19/34] minor update: use spaces for indentation --- src/problems/MHDBalsaraVortex/CMakeLists.txt | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/problems/MHDBalsaraVortex/CMakeLists.txt b/src/problems/MHDBalsaraVortex/CMakeLists.txt index 80067f3a8a..4b078cb703 100644 --- a/src/problems/MHDBalsaraVortex/CMakeLists.txt +++ b/src/problems/MHDBalsaraVortex/CMakeLists.txt @@ -9,8 +9,7 @@ if(AMReX_SPACEDIM EQUAL 3) add_test( NAME ${JOB_NAME} - COMMAND ${JOB_NAME} ../inputs/${JOB_NAME}.in ${QuokkaTestParams} + COMMAND ${JOB_NAME} ../inputs/${JOB_NAME}.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) - - set_tests_properties(MHDBalsaraVortex PROPERTIES LABELS "MHD") + set_tests_properties(MHDBalsaraVortex PROPERTIES LABELS "MHD") endif() From 0bf6cf980594698bf33b8b9e13aa6993a70750e5 Mon Sep 17 00:00:00 2001 From: neco kriel Date: Wed, 24 Dec 2025 11:08:51 +1100 Subject: [PATCH 20/34] add b magn param to input file --- inputs/MHDBalsaraVortex.in | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/inputs/MHDBalsaraVortex.in b/inputs/MHDBalsaraVortex.in index 46e8fd0bbe..6d17b3678c 100644 --- a/inputs/MHDBalsaraVortex.in +++ b/inputs/MHDBalsaraVortex.in @@ -13,12 +13,12 @@ amr.v = 1 # ***************************************************************** # Resolution and refinement # ***************************************************************** -amr.n_cell = 64 64 16 +amr.n_cell = 32 32 8 amr.max_level = 0 -amr.max_grid_size = 64 -amr.blocking_factor_x = 64 -amr.blocking_factor_y = 64 -amr.blocking_factor_z = 16 +amr.max_grid_size = 32 +amr.blocking_factor_x = 32 +amr.blocking_factor_y = 32 +amr.blocking_factor_z = 8 do_reflux = 0 do_subcycle = 0 @@ -34,5 +34,6 @@ hydro.reconstruction_order = 5 hydro.use_dual_energy = 0 setup.vortex_Mach = 0.1 +setup.vortex_b_magn = 0.1 setup.advection = 1 setup.num_orbits = 3 From c73fa1c800a477680ae45a59c80d69280ea64c08 Mon Sep 17 00:00:00 2001 From: neco kriel Date: Wed, 24 Dec 2025 11:09:31 +1100 Subject: [PATCH 21/34] fix: add req'ed params to turn dust off at compile time --- .../MHDBalsaraVortex/testMHDBalsaraVortex.cpp | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp b/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp index 4e204051a9..267fefe463 100644 --- a/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp +++ b/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp @@ -37,6 +37,8 @@ template <> struct Physics_Traits { static constexpr int numPassiveScalars = numMassScalars + 0; static constexpr bool is_self_gravity_enabled = false; static constexpr bool is_radiation_enabled = false; + static constexpr bool is_dust_enabled = false; + static constexpr int nDustGroups = 1; // number of dust groups static constexpr bool is_mhd_enabled = true; static constexpr int nGroups = 1; static constexpr UnitSystem unit_system = UnitSystem::CGS; @@ -261,15 +263,21 @@ auto problem_main() -> int } sim.stopTime_ = stop_time; - sim.computeReferenceSolution_ = true; sim.setInitialConditions(); sim.evolve(); - int status = 0; - const double error_tol = 0.002; - if (sim.errorNorm_ > error_tol) { - status = 1; + int status = 1; + const double error_tol = 0.005; + amrex::Real const error_norm = sim.computeErrorNorm(); + if (error_norm < error_tol) { + status = 0; + amrex::Print() << "Error norm = " << error_norm << "\n"; + amrex::Print() << "test passed\n"; + } else { + amrex::Print() << "Error norm = " << error_norm << "\n"; + amrex::Print() << "test failed\n"; } + return status; } From 8b3a0d655e7ca6127da1832afe9a6038d75685f8 Mon Sep 17 00:00:00 2001 From: neco kriel Date: Wed, 24 Dec 2025 11:10:42 +1100 Subject: [PATCH 22/34] use standard formula --- src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp b/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp index 2c85e54619..64e9397199 100644 --- a/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp +++ b/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp @@ -246,7 +246,7 @@ auto problem_main() -> int double stop_time = 0.0; if (is_advection_enabled) { const double advection_speed = vortex_u_magn; - vortex_drift_x2 = vortex_drift_x1 = advection_speed / std::sqrt(2.0); + vortex_drift_x2 = vortex_drift_x1 = advection_speed / std::numbers::sqrt2; const double length_x1 = sim.geom[0].ProbLength(0); const double length_x2 = sim.geom[0].ProbLength(1); if (std::abs(length_x1 - length_x2) > 1e-12) { From f2f350a73ecf984a98f6fc28675e5c363003a142 Mon Sep 17 00:00:00 2001 From: neco kriel Date: Fri, 26 Dec 2025 13:25:14 +1100 Subject: [PATCH 23/34] update: use consistent var name for velocity field --- .../MHDBalsaraVortex/testMHDBalsaraVortex.cpp | 30 +++++++++---------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp b/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp index 64e9397199..e46c032ad3 100644 --- a/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp +++ b/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp @@ -50,7 +50,7 @@ constexpr double bg_density = 1.0; constexpr double bg_pressure = 1.0; constexpr double sound_speed = gcem::sqrt(gamma_gas * bg_pressure / bg_density); // vortex parameters -AMREX_GPU_MANAGED double vortex_Mach = 0.01; // NOLINT +AMREX_GPU_MANAGED double vortex_Mach = 0.01; // NOLINT AMREX_GPU_MANAGED double vortex_b_magn = 0.01; // NOLINT // domain extends over [-5, 5] by default constexpr double vortex_center_x1 = 0.0; @@ -104,20 +104,20 @@ void computeVortexSolution(int i, int j, int k, amrex::Array4 const bg_pressure + 0.5 * (vortex_b_magn * vortex_b_magn * (1.0 - radius_sq) - vortex_u_magn * vortex_u_magn) * radial_profile_sq; const double Eint = pressure / (gamma_gas - 1.0); - const double delta_vel_x1 = -delta_x2_from_center * vortex_u_magn * radial_profile; - const double delta_vel_x2 = delta_x1_from_center * vortex_u_magn * radial_profile; - const double vel_x1 = vortex_drift_x1 + delta_vel_x1; - const double vel_x2 = vortex_drift_x2 + delta_vel_x2; - const double vel_x3 = 0.0; - const double mom_x1 = bg_density * vel_x1; - const double mom_x2 = bg_density * vel_x2; - const double mom_x3 = bg_density * vel_x3; - const double Ekin = 0.5 * bg_density * (vel_x1 * vel_x1 + vel_x2 * vel_x2 + vel_x3 * vel_x3); - - const double Bx_cc = -delta_x2_from_center * vortex_b_magn * radial_profile; - const double By_cc = delta_x1_from_center * vortex_b_magn * radial_profile; - const double Bz_cc = 0.0; - const double Emag = 0.5 * (Bx_cc * Bx_cc + By_cc * By_cc + Bz_cc * Bz_cc); + const double delta_u_x1 = -delta_x2_from_center * vortex_u_magn * radial_profile; + const double delta_u_x2 = delta_x1_from_center * vortex_u_magn * radial_profile; + const double u_x1 = vortex_drift_x1 + delta_u_x1; + const double u_x2 = vortex_drift_x2 + delta_u_x2; + const double u_x3 = 0.0; + const double mom_x1 = bg_density * u_x1; + const double mom_x2 = bg_density * u_x2; + const double mom_x3 = bg_density * u_x3; + const double Ekin = 0.5 * bg_density * (u_x1 * u_x1 + u_x2 * u_x2 + u_x3 * u_x3); + + const double b_x1 = -delta_x2_from_center * vortex_b_magn * radial_profile; + const double b_x2 = delta_x1_from_center * vortex_b_magn * radial_profile; + const double b_x3 = 0.0; + const double Emag = 0.5 * (b_x1 * b_x1 + b_x2 * b_x2 + b_x3 * b_x3); const double Etot = Eint + Ekin + Emag; From 5c7e81bec5fce17d372546c8f394a12d926ca701 Mon Sep 17 00:00:00 2001 From: neco kriel Date: Fri, 26 Dec 2025 13:25:55 +1100 Subject: [PATCH 24/34] update: reformat correction eval --- src/hydro/HLLD.hpp | 55 ++++++++++++++++++++++++++++++---------------- 1 file changed, 36 insertions(+), 19 deletions(-) diff --git a/src/hydro/HLLD.hpp b/src/hydro/HLLD.hpp index 2dce38406e..8db04fe5a6 100644 --- a/src/hydro/HLLD.hpp +++ b/src/hydro/HLLD.hpp @@ -41,16 +41,20 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState u_star_R{}; // frequently used term + const double vel_magn_sq_L = SQUARE(sL.u) + (SQUARE(sL.v) + SQUARE(sL.w)); + const double vel_magn_sq_R = SQUARE(sR.u) + (SQUARE(sR.v) + SQUARE(sR.w)); const double bx_sq = SQUARE(bx); + const double b_magn_sq_L = bx_sq + (SQUARE(sL.by) + SQUARE(sL.bz)); + const double b_magn_sq_R = bx_sq + (SQUARE(sR.by) + SQUARE(sR.bz)); // compute L/R states for select conserved variables // (group transverse vector components for floating-point associativity symmetry) // magnetic pressure - const double pb_L = 0.5 * (bx_sq + (SQUARE(sL.by) + SQUARE(sL.bz))); - const double pb_R = 0.5 * (bx_sq + (SQUARE(sR.by) + SQUARE(sR.bz))); + const double pb_L = 0.5 * b_magn_sq_L; + const double pb_R = 0.5 * b_magn_sq_R; // kinetic energy - const double ke_L = 0.5 * sL.rho * (SQUARE(sL.u) + (SQUARE(sL.v) + SQUARE(sL.w))); - const double ke_R = 0.5 * sR.rho * (SQUARE(sR.u) + (SQUARE(sR.v) + SQUARE(sR.w))); + const double ke_L = 0.5 * (sL.rho * vel_magn_sq_L); + const double ke_R = 0.5 * (sR.rho * vel_magn_sq_R); // set left conserved states u_L.rho = sL.rho; u_L.mx = sL.u * u_L.rho; @@ -155,22 +159,35 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState 0 for very low Mach, phi -> 1 as advection becomes comparable to wave speeds) - const double spd_fss_L = std::abs(sL.u) + spd_fms_L; - const double spd_fss_R = std::abs(sR.u) + spd_fms_R; - const double spd_fss_max = std::max(spd_fss_L, spd_fss_R); - const double chi = std::min(1.0, spd_fss_max / std::max(1e-14, spd_fms_max)); + // // total pressure (no correction) + // // MK5: eqn (41) can be calculated (more explicitly) via eqn (23) + // double ptot_star_L = ptot_L + u_L.rho * siui_L * (spds[2] - sL.u); + // double ptot_star_R = ptot_R + u_R.rho * siui_R * (spds[2] - sR.u); + // double const ptot_star = 0.5 * (ptot_star_L + ptot_star_R); + + // total pressure (w/ MM21 low-Mach correction) + // MM21 eqn 8 + const double spd_ca_sq_L = b_magn_sq_L / sL.rho; + const double spd_ca_sq_R = b_magn_sq_R / sR.rho; + const double spd_cax_sq_L = bx_sq / sL.rho; + const double spd_cax_sq_R = bx_sq / sR.rho; + // cu := modified fast magnetosonic speed (MM21 eqn 14) + const double cu_sqrt_arg_L = std::max(0.0, SQUARE(spd_ca_sq_L + vel_magn_sq_L) - 4.0 * vel_magn_sq_L * spd_cax_sq_L); + const double cu_sqrt_arg_R = std::max(0.0, SQUARE(spd_ca_sq_R + vel_magn_sq_R) - 4.0 * vel_magn_sq_R * spd_cax_sq_R); + const double spd_cu_sq_L = 0.5 * ((spd_ca_sq_L + vel_magn_sq_L) + std::sqrt(cu_sqrt_arg_L)); + const double spd_cu_sq_R = 0.5 * ((spd_ca_sq_R + vel_magn_sq_R) + std::sqrt(cu_sqrt_arg_R)); + const double spd_cu_L = std::sqrt(std::max(0.0, spd_cu_sq_L)); + const double spd_cu_R = std::sqrt(std::max(0.0, spd_cu_sq_R)); + const double spd_cu_max = std::max(spd_cu_L, spd_cu_R); + // limiter (MM21 eqn 16) + const double chi = std::min(1.0, spd_cu_max / std::max(1e-14, spd_fms_max)); const double phi = chi * (2.0 - chi); - // MK5 star total pressures estimate (eqn 23) - const double ptot_star_L = ptot_L + u_L.rho * siui_L * (spds[2] - sL.u); - const double ptot_star_R = ptot_R + u_R.rho * siui_R * (spds[2] - sR.u); - const double ptot_star_avg_LR = 0.5 * (ptot_star_L + ptot_star_R); - // correct the excess (compressive) contributions in the low Mach regime (MM21) - const double ptot_avg_LR = 0.5 * (ptot_L + ptot_R); - const double ptot_star_excess = ptot_star_avg_LR - ptot_avg_LR; - const double ptot_star = ptot_avg_LR + phi * ptot_star_excess; + // corrected total star pressure (MM21 eqn 15) + const double ptot_star_numer_term1 = (siui_R * u_R.rho) * ptot_L; + const double ptot_star_numer_term2 = (siui_L * u_L.rho) * ptot_R; + const double ptot_star_numer_term3 = phi * (u_L.rho * u_R.rho) * (siui_R * siui_L) * (sR.u - sL.u); + const double ptot_star_numer = ptot_star_numer_term1 - ptot_star_numer_term2 + ptot_star_numer_term3; + const double ptot_star = ptot_star_numer / sm_denom; // MK5: u_L^(star, dstar) from, eqn (39) u_star_L.mx = u_star_L.rho * spds[2]; From 12182ef913cf8b32f14f2b92e187e358f6ee87fb Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 26 Dec 2025 02:26:13 +0000 Subject: [PATCH 25/34] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp b/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp index e46c032ad3..7afdf8c25e 100644 --- a/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp +++ b/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp @@ -50,7 +50,7 @@ constexpr double bg_density = 1.0; constexpr double bg_pressure = 1.0; constexpr double sound_speed = gcem::sqrt(gamma_gas * bg_pressure / bg_density); // vortex parameters -AMREX_GPU_MANAGED double vortex_Mach = 0.01; // NOLINT +AMREX_GPU_MANAGED double vortex_Mach = 0.01; // NOLINT AMREX_GPU_MANAGED double vortex_b_magn = 0.01; // NOLINT // domain extends over [-5, 5] by default constexpr double vortex_center_x1 = 0.0; From 1baf5e85e9c69bc49b4b510e0f2fd4c63a1a0732 Mon Sep 17 00:00:00 2001 From: neco kriel Date: Fri, 26 Dec 2025 18:23:53 +1100 Subject: [PATCH 26/34] minor consistency updates (ty gemini) --- src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp b/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp index 7afdf8c25e..87b4cf2479 100644 --- a/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp +++ b/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp @@ -9,6 +9,7 @@ #include #include #include +#include #include "AMReX_Array.H" #include "AMReX_Array4.H" @@ -129,8 +130,8 @@ void computeVortexSolution(int i, int j, int k, amrex::Array4 const state(i, j, k, HydroSystem::internalEnergy_index) = Eint; } else if (cen == quokka::centering::fc) { - const auto b_x1_L = (Az(x1_L, x2_L + static_cast(dx[1])) - Az(x1_L, x2_L)) / static_cast(dx[1]); - const auto b_x2_L = (Az(x1_L, x2_L) - Az(x1_L + static_cast(dx[0]), x2_L)) / static_cast(dx[0]); + const auto b_x1_L = (Az(x1_L, x2_L + dx[1]) - Az(x1_L, x2_L)) / dx[1]; + const auto b_x2_L = (Az(x1_L, x2_L) - Az(x1_L + dx[0], x2_L)) / dx[0]; const auto b_x3_L = 0.0; if (dir == quokka::direction::x) { @@ -257,7 +258,7 @@ auto problem_main() -> int stop_time = static_cast(num_orbits) * advection_duration; } else { vortex_drift_x1 = vortex_drift_x2 = 0.0; - const double orbital_duration = 2.0 * M_PI / vortex_u_magn; + const double orbital_duration = 2.0 * std::numbers::pi / vortex_u_magn; stop_time = static_cast(num_orbits) * orbital_duration; } From c98d8be444fea5e67173e367bd373c9acfc47c78 Mon Sep 17 00:00:00 2001 From: neco kriel Date: Fri, 26 Dec 2025 18:35:48 +1100 Subject: [PATCH 27/34] update: add vortex radius param --- .../MHDBalsaraVortex/testMHDBalsaraVortex.cpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp b/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp index 87b4cf2479..a6d7909e3f 100644 --- a/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp +++ b/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp @@ -51,6 +51,7 @@ constexpr double bg_density = 1.0; constexpr double bg_pressure = 1.0; constexpr double sound_speed = gcem::sqrt(gamma_gas * bg_pressure / bg_density); // vortex parameters +AMREX_GPU_MANAGED double vortex_radius = 1.0; // NOLINT AMREX_GPU_MANAGED double vortex_Mach = 0.01; // NOLINT AMREX_GPU_MANAGED double vortex_b_magn = 0.01; // NOLINT // domain extends over [-5, 5] by default @@ -102,11 +103,11 @@ void computeVortexSolution(int i, int j, int k, amrex::Array4 const const double vortex_u_magn = vortex_Mach * sound_speed; const double pressure = - bg_pressure + 0.5 * (vortex_b_magn * vortex_b_magn * (1.0 - radius_sq) - vortex_u_magn * vortex_u_magn) * radial_profile_sq; + bg_pressure + 0.5 * (vortex_b_magn * vortex_b_magn * (1.0 - radius_sq) - bg_density * vortex_u_magn * vortex_u_magn) * radial_profile_sq; const double Eint = pressure / (gamma_gas - 1.0); - const double delta_u_x1 = -delta_x2_from_center * vortex_u_magn * radial_profile; - const double delta_u_x2 = delta_x1_from_center * vortex_u_magn * radial_profile; + const double delta_u_x1 = -(delta_x2_from_center / vortex_radius) * vortex_u_magn * radial_profile; + const double delta_u_x2 = (delta_x1_from_center / vortex_radius) * vortex_u_magn * radial_profile; const double u_x1 = vortex_drift_x1 + delta_u_x1; const double u_x2 = vortex_drift_x2 + delta_u_x2; const double u_x3 = 0.0; @@ -115,8 +116,8 @@ void computeVortexSolution(int i, int j, int k, amrex::Array4 const const double mom_x3 = bg_density * u_x3; const double Ekin = 0.5 * bg_density * (u_x1 * u_x1 + u_x2 * u_x2 + u_x3 * u_x3); - const double b_x1 = -delta_x2_from_center * vortex_b_magn * radial_profile; - const double b_x2 = delta_x1_from_center * vortex_b_magn * radial_profile; + const double b_x1 = -(delta_x2_from_center / vortex_radius) * vortex_b_magn * radial_profile; + const double b_x2 = (delta_x1_from_center / vortex_radius) * vortex_b_magn * radial_profile; const double b_x3 = 0.0; const double Emag = 0.5 * (b_x1 * b_x1 + b_x2 * b_x2 + b_x3 * b_x3); From 415bbcddee41dbbbbbde525c3964ba9c7f546b6b Mon Sep 17 00:00:00 2001 From: neco kriel Date: Fri, 26 Dec 2025 19:56:29 +1100 Subject: [PATCH 28/34] fix: clamp fms speeds + sgn of shock detection --- src/hydro/HLLD.hpp | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/hydro/HLLD.hpp b/src/hydro/HLLD.hpp index 8db04fe5a6..174edab491 100644 --- a/src/hydro/HLLD.hpp +++ b/src/hydro/HLLD.hpp @@ -83,8 +83,11 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState compression + const double para_v_jump = sR.u - sL.u; // negative -> compression // tp := shock anisotropy, clamped to [0, 1], with theta = tp^4 const double denom_tp = std::max(1e-14, spd_fms_max - std::min(perp_v_jump, 0.0)); double tp = (spd_fms_max - std::min(para_v_jump, 0.0)) / denom_tp; From c3d531f2aa2f951dc085e59096d1f288d86d6f8e Mon Sep 17 00:00:00 2001 From: neco kriel Date: Sun, 28 Dec 2025 12:19:29 +1100 Subject: [PATCH 29/34] update: remove redundant min/max calls --- src/hydro/HLLD.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/hydro/HLLD.hpp b/src/hydro/HLLD.hpp index 174edab491..612c1becf5 100644 --- a/src/hydro/HLLD.hpp +++ b/src/hydro/HLLD.hpp @@ -88,8 +88,8 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState Date: Sun, 28 Dec 2025 12:20:18 +1100 Subject: [PATCH 30/34] update: remove old estimate for middle state pressure --- src/hydro/HLLD.hpp | 6 ------ 1 file changed, 6 deletions(-) diff --git a/src/hydro/HLLD.hpp b/src/hydro/HLLD.hpp index 612c1becf5..4c88fae17d 100644 --- a/src/hydro/HLLD.hpp +++ b/src/hydro/HLLD.hpp @@ -161,12 +161,6 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState Date: Sun, 28 Dec 2025 12:21:29 +1100 Subject: [PATCH 31/34] update: remove balsara vortex problem from ASAN suite: it takes too long --- src/problems/MHDBalsaraVortex/CMakeLists.txt | 1 - 1 file changed, 1 deletion(-) diff --git a/src/problems/MHDBalsaraVortex/CMakeLists.txt b/src/problems/MHDBalsaraVortex/CMakeLists.txt index 4b078cb703..9641393a5f 100644 --- a/src/problems/MHDBalsaraVortex/CMakeLists.txt +++ b/src/problems/MHDBalsaraVortex/CMakeLists.txt @@ -11,5 +11,4 @@ if(AMReX_SPACEDIM EQUAL 3) NAME ${JOB_NAME} COMMAND ${JOB_NAME} ../inputs/${JOB_NAME}.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) - set_tests_properties(MHDBalsaraVortex PROPERTIES LABELS "MHD") endif() From 0ff84a161662663e6b785945fd33aca0fb989817 Mon Sep 17 00:00:00 2001 From: neco kriel Date: Sun, 28 Dec 2025 20:44:08 +1100 Subject: [PATCH 32/34] minor update: add vortex radius param around the place for consistency --- .../MHDBalsaraVortex/testMHDBalsaraVortex.cpp | 23 +++++++++++++------ 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp b/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp index a6d7909e3f..6e7f3c8309 100644 --- a/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp +++ b/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp @@ -65,21 +65,27 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto computeRadiusSq(const double x1, const { const double delta_x1_from_center = x1 - vortex_center_x1; const double delta_x2_from_center = x2 - vortex_center_x2; - return delta_x1_from_center * delta_x1_from_center + delta_x2_from_center * delta_x2_from_center; + const double vortex_radius_inv = 1.0 / vortex_radius; + const double delta_x1_scaled = delta_x1_from_center * vortex_radius_inv; + const double delta_x2_scaled = delta_x2_from_center * vortex_radius_inv; + return delta_x1_scaled * delta_x1_scaled + delta_x2_scaled * delta_x2_scaled; // (r/R)^2 } -AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto computeRadialProfile(const double radius_sq) -> double { return std::exp(0.5 * (1.0 - radius_sq)); } +AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto computeRadialProfile(const double radius_sq) -> double +{ + return std::exp(0.5 * (1.0 - radius_sq)); +} AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto Az(const double x1, const double x2) -> double { // vec(a) = (0, 0, a_z), with: - // a_z = b_magn * exp((1 - r^2)/2) - // so that: - // b_x = d(a_z)/dy = -y * b_magn * exp((1 - r^2)/2) - // b_y = -d(a_z)/dx = x * b_magn * exp((1 - r^2)/2) + // a_z = b_magn * R * exp((1 - (r/R)^2)/2) + // so: + // b_x = d(a_z)/dy = -(y/R) * b_magn * exp((1 - (r/R)^2)/2) + // b_y = -d(a_z)/dx = (x/R) * b_magn * exp((1 - (r/R)^2)/2) const double radius_sq = computeRadiusSq(x1, x2); const double radial_profile = computeRadialProfile(radius_sq); - return vortex_b_magn * radial_profile; + return vortex_b_magn * vortex_radius * radial_profile; } AMREX_GPU_DEVICE @@ -231,6 +237,9 @@ auto problem_main() -> int hpp.query("num_orbits", num_orbits); const double vortex_u_magn = vortex_Mach * sound_speed; const bool is_advection_enabled = (advection_int != 0); + if (vortex_radius <= 0.0) { + amrex::Abort("vortex_radius must be > 0."); + } auto BCs_cc = quokka::BC(quokka::BCType::int_dir); From 11e31fb47d450052a3f192b84c0ac0a5bcfdc554 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Sun, 28 Dec 2025 09:44:22 +0000 Subject: [PATCH 33/34] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp b/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp index 6e7f3c8309..e1ade14ad1 100644 --- a/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp +++ b/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp @@ -71,10 +71,7 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto computeRadiusSq(const double x1, const return delta_x1_scaled * delta_x1_scaled + delta_x2_scaled * delta_x2_scaled; // (r/R)^2 } -AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto computeRadialProfile(const double radius_sq) -> double -{ - return std::exp(0.5 * (1.0 - radius_sq)); -} +AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto computeRadialProfile(const double radius_sq) -> double { return std::exp(0.5 * (1.0 - radius_sq)); } AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto Az(const double x1, const double x2) -> double { From b961b05a1ad657e34c10a5bc69fe64958c106e9b Mon Sep 17 00:00:00 2001 From: Neco Kriel Date: Tue, 30 Dec 2025 14:26:47 +1100 Subject: [PATCH 34/34] tmp update: add old pressure estimate for comparison --- src/hydro/HLLD.hpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/hydro/HLLD.hpp b/src/hydro/HLLD.hpp index 4c88fae17d..612c1becf5 100644 --- a/src/hydro/HLLD.hpp +++ b/src/hydro/HLLD.hpp @@ -161,6 +161,12 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState