diff --git a/inputs/MHDBalsaraVortex.in b/inputs/MHDBalsaraVortex.in
new file mode 100644
index 0000000000..6d17b3678c
--- /dev/null
+++ b/inputs/MHDBalsaraVortex.in
@@ -0,0 +1,39 @@
+# *****************************************************************
+# Problem size and geometry
+# *****************************************************************
+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
+# *****************************************************************
+amr.v = 1
+
+# *****************************************************************
+# Resolution and refinement
+# *****************************************************************
+amr.n_cell = 32 32 8
+amr.max_level = 0
+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
+
+plotfile_interval = 10
+plotfile_prefix = balsara_vortex_plt
+checkpoint_interval = -1
+cfl = 0.3
+max_timesteps = 10000
+
+hydro.rk_integrator_order = 2
+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
diff --git a/src/hydro/HLLD.hpp b/src/hydro/HLLD.hpp
index 7ea1782f3e..612c1becf5 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 <typename problem_t, int N_scalars, int N_mscalars, int fluxdim>
 AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState<N_scalars, N_mscalars> const &sL, quokka::HydroState<N_scalars, N_mscalars> const &sR,
 					      const double gamma, const double bx, const double perp_v_jump)
@@ -40,16 +41,20 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState<N_scalars, N_ms
 	ConsHydro1D<N_scalars> u_star_R{};
 
 	// frequently used term
-	double const bx_sq = SQUARE(bx);
+	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
-	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 * b_magn_sq_L;
+	const double pb_R = 0.5 * b_magn_sq_R;
 	// 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 * 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;
@@ -76,18 +81,21 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState<N_scalars, N_ms
 
 	//--- Step 2. Compute L & R wave speeds according to MK5, eqn. (67)
 
-	const double fspd_L = FastMagnetoSonicSpeed(gamma, sL, bx);
-	const double fspd_R = FastMagnetoSonicSpeed(gamma, sR, bx);
-	spds[0] = std::min(sL.u - fspd_L, sR.u - fspd_R);
-	spds[4] = std::max(sL.u + fspd_L, sR.u + fspd_R);
-	const double fspd_m = -std::min(0.0, spds[0]);
-	const double fspd_p = std::max(0.0, spds[4]);
+	const double spd_fms_L = FastMagnetoSonicSpeed(gamma, sL, bx);
+	const double spd_fms_R = FastMagnetoSonicSpeed(gamma, sR, bx);
+	const double spd_fms_max = std::max(spd_fms_L, spd_fms_R);
+	const double u_min = std::min(sL.u, sR.u);
+	const double u_max = std::max(sL.u, sR.u);
+	spds[0] = std::min(0.0, u_min - spd_fms_max);
+	spds[4] = std::max(0.0, u_max + spd_fms_max);
+	const double fspd_m = -spds[0];
+	const double fspd_p = spds[4];
 
 	//--- Step 3. Compute L/R fluxes
 
 	// total pressure
-	double ptot_L = sL.P + pb_L;
-	double ptot_R = sR.P + pb_R;
+	const double ptot_L = sL.P + pb_L;
+	const double ptot_R = sR.P + pb_R;
 	// fluxes on the left side of the interface
 	f_L.rho = u_L.mx;
 	f_L.mx = u_L.mx * sL.u + ptot_L - bx_sq;
@@ -116,25 +124,23 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState<N_scalars, N_ms
 	//--- Step 4. Compute middle and Alfven wave speeds
 
 	// MK5: S_i - u_i (for i=L or R)
-	double siui_L = spds[0] - sL.u;
-	double siui_R = spds[4] - sR.u;
+	const double siui_L = spds[0] - sL.u;
+	const double siui_R = spds[4] - sR.u;
 	// carbuncle detector
-	const double max_spd = std::max(fspd_L, fspd_R);
-	const double para_v_jump = sL.u - sR.u; // negative -> compression
-	double theta = 1.0;
+	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, 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)
-	double sism_L = spds[0] - spds[2];
-	double sism_R = spds[4] - spds[2];
-	double sism_inv_L = 1.0 / sism_L;
-	double sism_inv_R = 1.0 / sism_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;
 	// MK5: rho_i from eqn (43)
 	u_star_L.rho = u_L.rho * siui_L * sism_inv_L;
 	u_star_R.rho = u_R.rho * siui_R * sism_inv_R;
@@ -145,21 +151,45 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState<N_scalars, N_ms
 		u_star_R.scalar[n] = u_R.scalar[n] * siui_R * sism_inv_R;
 	}
 
-	double u_star_rho_inv_L = 1.0 / u_star_L.rho;
-	double u_star_rho_inv_R = 1.0 / u_star_R.rho;
-	double rho_sqrt_L = std::sqrt(u_star_L.rho);
-	double rho_sqrt_R = std::sqrt(u_star_R.rho);
+	const double u_star_rho_inv_L = 1.0 / u_star_L.rho;
+	const double u_star_rho_inv_R = 1.0 / u_star_R.rho;
+	const double rho_sqrt_L = std::sqrt(u_star_L.rho);
+	const double rho_sqrt_R = std::sqrt(u_star_R.rho);
 	// MK5: eqn (51)
 	spds[1] = spds[2] - std::abs(bx) / rho_sqrt_L;
 	spds[3] = spds[2] + std::abs(bx) / rho_sqrt_R;
 
 	//--- Step 5. Compute intermediate states
 
-	// compute total pressure
-	// 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 (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);
+	// 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];
@@ -205,7 +235,7 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState<N_scalars, N_ms
 	}
 	// vec(v_R^star) dot vec(b_R^star)
 	// group transverse momenta-components for floating-point associativity
-	double vb_star_R = (u_star_R.mx * bx + (u_star_R.my * u_star_R.by + u_star_R.mz * u_star_R.bz)) * u_star_rho_inv_R;
+	const double vb_star_R = (u_star_R.mx * bx + (u_star_R.my * u_star_R.by + u_star_R.mz * u_star_R.bz)) * u_star_rho_inv_R;
 	// MK5: eqn (48)
 	u_star_R.E = (siui_R * u_R.E - ptot_R * sR.u + ptot_star * spds[2] + bx * (sR.u * bx + (sR.v * u_R.by + sR.w * u_R.bz) - vb_star_R)) * sism_inv_R;
 
diff --git a/src/problems/MHDBalsaraVortex/CMakeLists.txt b/src/problems/MHDBalsaraVortex/CMakeLists.txt
new file mode 100644
index 0000000000..9641393a5f
--- /dev/null
+++ b/src/problems/MHDBalsaraVortex/CMakeLists.txt
@@ -0,0 +1,14 @@
+if(AMReX_SPACEDIM EQUAL 3)
+  set(JOB_NAME MHDBalsaraVortex)
+
+  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
+    WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests)
+endif()
diff --git a/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp b/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp
new file mode 100644
index 0000000000..e1ade14ad1
--- /dev/null
+++ b/src/problems/MHDBalsaraVortex/testMHDBalsaraVortex.cpp
@@ -0,0 +1,290 @@
+//==============================================================================
+// Copyright 2025 Neco Kriel.
+// Released under the MIT license. See LICENSE file included in the GitHub repo.
+//==============================================================================
+/// \file testMHDBalsaraVortex.cpp
+/// \brief advecting MHD Balsara vortex.
+///
+
+#include <cassert>
+#include <cmath>
+#include <gcem.hpp>
+#include <numbers>
+
+#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<MHDBalsaraVortex> {
+	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<MHDBalsaraVortex> {
+	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_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;
+};
+
+// background state
+constexpr double gamma_gas = quokka::EOS_Traits<MHDBalsaraVortex>::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_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
+constexpr double vortex_center_x1 = 0.0;
+constexpr double vortex_center_x2 = 0.0;
+// 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;
+	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 Az(const double x1, const double x2) -> double
+{
+	// vec(a) = (0, 0, a_z), with:
+	//   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 * vortex_radius * radial_profile;
+}
+
+AMREX_GPU_DEVICE
+void computeVortexSolution(int i, int j, int k, amrex::Array4<amrex::Real> const &state, amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const &dx,
+			   amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> 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];
+
+	if (cen == quokka::centering::cc) {
+		const amrex::Real x1_C = x1_L + static_cast<amrex::Real>(0.5) * dx[0];
+		const amrex::Real x2_C = x2_L + static_cast<amrex::Real>(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) - 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_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;
+		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_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);
+
+		const double Etot = Eint + Ekin + Emag;
+
+		state(i, j, k, HydroSystem<MHDBalsaraVortex>::density_index) = bg_density;
+		state(i, j, k, HydroSystem<MHDBalsaraVortex>::x1Momentum_index) = mom_x1;
+		state(i, j, k, HydroSystem<MHDBalsaraVortex>::x2Momentum_index) = mom_x2;
+		state(i, j, k, HydroSystem<MHDBalsaraVortex>::x3Momentum_index) = mom_x3;
+		state(i, j, k, HydroSystem<MHDBalsaraVortex>::energy_index) = Etot;
+		state(i, j, k, HydroSystem<MHDBalsaraVortex>::internalEnergy_index) = Eint;
+
+	} else if (cen == quokka::centering::fc) {
+		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) {
+			state(i, j, k, MHDSystem<MHDBalsaraVortex>::bfield_index) = b_x1_L;
+		} else if (dir == quokka::direction::y) {
+			state(i, j, k, MHDSystem<MHDBalsaraVortex>::bfield_index) = b_x2_L;
+		} else if (dir == quokka::direction::z) {
+			state(i, j, k, MHDSystem<MHDBalsaraVortex>::bfield_index) = b_x3_L;
+		}
+	}
+}
+
+template <> void QuokkaSimulation<MHDBalsaraVortex>::setInitialConditionsOnGrid(quokka::grid const &grid_elem)
+{
+	const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> dx = grid_elem.dx_;
+	const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> prob_lo = grid_elem.prob_lo_;
+	const amrex::Array4<amrex::Real> &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<MHDBalsaraVortex>::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; // fill unused quantities with zeros
+		}
+		computeVortexSolution(i, j, k, state_cc, dx, prob_lo, cen, dir);
+	});
+}
+
+template <> void QuokkaSimulation<MHDBalsaraVortex>::setInitialConditionsOnGridFaceVars(quokka::grid const &grid_elem)
+{
+	const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> dx = grid_elem.dx_;
+	const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> prob_lo = grid_elem.prob_lo_;
+	const amrex::Array4<amrex::Real> &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<MHDBalsaraVortex>::nvarPerDim_fc;
+
+	amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
+		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);
+	});
+}
+
+template <>
+void QuokkaSimulation<MHDBalsaraVortex>::computeReferenceSolution(amrex::MultiFab &ref, amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const &dx,
+								  amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> 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; // fill unused quantities with zeros
+			}
+			computeVortexSolution(i, j, k, stateExact, dx, prob_lo, quokka::centering::cc, quokka::direction::na);
+		});
+	}
+}
+
+template <>
+void QuokkaSimulation<MHDBalsaraVortex>::computeReferenceSolution_fc(amrex::MultiFab &ref, amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const &dx,
+								     amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> 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, quokka::centering::fc, dir);
+		});
+	}
+}
+
+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("vortex_b_magn", vortex_b_magn);
+	hpp.query("advection", advection_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<MHDBalsaraVortex>(quokka::BCType::int_dir);
+
+	const int nvars_fc = Physics_Indices<MHDBalsaraVortex>::nvarTotal_fc;
+	amrex::Vector<amrex::BCRec> 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<MHDBalsaraVortex> sim(BCs_cc, BCs_fc);
+
+	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::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) {
+			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_u_magn;
+		stop_time = static_cast<double>(num_orbits) * advection_duration;
+	} else {
+		vortex_drift_x1 = vortex_drift_x2 = 0.0;
+		const double orbital_duration = 2.0 * std::numbers::pi / vortex_u_magn;
+		stop_time = static_cast<double>(num_orbits) * orbital_duration;
+	}
+
+	sim.stopTime_ = stop_time;
+
+	sim.setInitialConditions();
+	sim.evolve();
+
+	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;
+}