Fix resampled cooling internal energy with MHD (#1886)

## Summary
- Pass face-centered MHD state into `ResampledCooling::computeCooling`.
- Recover gas internal energy in the resampled cooling module as total
gas energy minus kinetic and magnetic energy.
- Convert `ResampledCoolingTest` to exercise the MHD path in 3D with a
uniform magnetic field initialized at plasma beta ~ 1.

Fixes #1885.

## Validation
- `scripts/bash/quokka config -d 3d-debug`
- `cmake --build build/3d-debug --target ResampledCoolingTest`
- `ctest --test-dir build/3d-debug -R ResampledCoolingTest
--output-on-failure`

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

Author: BenWibking

src/QuokkaSimulation.hpp

@@ -1006,8 +1006,8 @@ auto QuokkaSimulation<problem_t>::addStrangSplitSourcesWithBuiltin(amrex::MultiF
 				coolingTableType_ = "resampled";
 			}
 			if (coolingTableType_ == "resampled") {
-				cool_success =
-				    quokka::ResampledCooling::computeCooling<problem_t>(state, dt, resampledTables_, tempFloor_, const_heating_rate_per_H);
+				cool_success = quokka::ResampledCooling::computeCooling<problem_t>(state, state_fc, dt, resampledTables_, tempFloor_,
+												   const_heating_rate_per_H);
 			} else {
 				amrex::Abort("Invalid cooling table type! Only 'resampled' is supported.");
 			}

src/cooling/ResampledCooling.hpp

@@ -19,7 +19,6 @@
 #include "math/FastMath.hpp"
 #include "math/ODEIntegrate.hpp"
 #include "math/root_finding.hpp"
-#include "radiation/radiation_system.hpp"
 #include "util/DataTable.hpp"
 #include <format>
 
@@ -195,8 +194,8 @@ struct ResampledCoolingFunctor {
 
 // const_heating_rate_per_H: unit erg/s/H
 template <typename problem_t>
-auto computeCooling(amrex::MultiFab &mf, const Real dt_in, resampled_tables &resampledTables, const Real temp_floor, const Real const_heating_rate_per_H)
-    -> bool
+auto computeCooling(amrex::MultiFab &mf, std::array<amrex::MultiFab, AMREX_SPACEDIM> const &mf_fc, const Real dt_in, resampled_tables &resampledTables,
+		    const Real temp_floor, const Real const_heating_rate_per_H) -> bool
 {
 	const BL_PROFILE("quokka::ResampledCooling::computeCooling()");
 
@@ -219,21 +218,26 @@ auto computeCooling(amrex::MultiFab &mf, const Real dt_in, resampled_tables &res
 		auto const &state = mf.array(iter);
 		auto const &nsubsteps = nsubstepsMF.array(iter);
 
+		std::array<amrex::Array4<const amrex::Real>, AMREX_SPACEDIM> state_fc{};
+		if constexpr (Physics_Traits<problem_t>::is_mhd_enabled) {
+			for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
+				state_fc[dir] = mf_fc[dir].const_array(iter);
+			}
+		}
+
 		amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+			// cooling function
 			const Real rho = state(i, j, k, HydroSystem<problem_t>::density_index);
-			const Real x1Mom = state(i, j, k, HydroSystem<problem_t>::x1Momentum_index);
-			const Real x2Mom = state(i, j, k, HydroSystem<problem_t>::x2Momentum_index);
-			const Real x3Mom = state(i, j, k, HydroSystem<problem_t>::x3Momentum_index);
-			const Real Egas = state(i, j, k, HydroSystem<problem_t>::energy_index);
+			const Real nH = rho * tables.cloudy_H_mass_fraction / C::m_p;			    // unit: cm^-3
+			const ResampledCoolingFunctor user_rhs(rho, tables, const_heating_rate_per_H * nH); // unit: erg/cm^3/s
 
-			// number density
-			const Real nH = rho * tables.cloudy_H_mass_fraction / C::m_p; // unit: cm^-3
+			// state vector
+			const Real Eint = HydroSystem<problem_t>::ComputeInternalEnergy(state, i, j, k, &state_fc);
+			quokka::valarray<Real, 1> y = {Eint};
 
-			const Real Eint = RadSystem<problem_t>::ComputeEintFromEgas(rho, x1Mom, x2Mom, x3Mom, Egas);
+			// integration tolerance
 			const Real Eint_floor = (temp_floor > 0.0) ? ComputeEgasFromTgas(rho, temp_floor, tables) : 0.0;
 			const Real abstol_floor = amrex::max(Eint_floor, std::numeric_limits<Real>::min());
-			const ResampledCoolingFunctor user_rhs(rho, tables, const_heating_rate_per_H * nH); // unit: erg/cm^3/s
-			quokka::valarray<Real, 1> y = {Eint};
 			quokka::valarray<Real, 1> const abstol = {reltol_floor * abstol_floor};
 
 			// do integration with RK2 (Heun's method)

src/problems/ResampledCoolingTest/testResampledCoolingTest.cpp

@@ -15,9 +15,13 @@
 #include "cooling/ResampledCooling.hpp"
 #include "math/interpolate.hpp"
 #include "util/BC.hpp"
+#include <algorithm>
+#include <array>
+#include <cmath>
 #include <format>
 #include <fstream>
 #include <iomanip>
+#include <limits>
 #include <sstream>
 #include <vector>
 
@@ -97,14 +101,25 @@ template <> struct Physics_Traits<ResampledCoolingTest> {
 	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 = false;
+	static constexpr bool is_mhd_enabled = (AMREX_SPACEDIM == 3);
 	static constexpr int nGroups = 1; // number of radiation groups
 	static constexpr UnitSystem unit_system = UnitSystem::CGS;
 };
 
 // Initial conditions: hot gas that will cool down
 constexpr double T_initial = 1.0e7;	// K
 constexpr double rho_initial = 1.0e-26; // g cm^-3 (constant density for isochoric)
+constexpr double mu_initial = 0.6 * quokka::EOS_Traits<ResampledCoolingTest>::mean_molecular_weight;
+constexpr double pressure_initial = rho_initial * C::k_B * T_initial / mu_initial;
+constexpr double Bx_initial = 5.264491941623788e-06; // chosen so plasma beta = P_gas / (B^2 / 2) ~= 1
+constexpr double magnetic_energy_initial = 0.5 * Bx_initial * Bx_initial;
+constexpr double active_magnetic_energy_initial = Physics_Traits<ResampledCoolingTest>::is_mhd_enabled ? magnetic_energy_initial : 0.0;
+
+auto computeInitialInternalEnergy() -> double
+{
+	constexpr double gamma = quokka::EOS_Traits<ResampledCoolingTest>::gamma;
+	return rho_initial * C::k_B * T_initial / ((gamma - 1.0) * mu_initial);
+}
 
 template <> void QuokkaSimulation<ResampledCoolingTest>::setInitialConditionsOnGrid(quokka::grid const &grid_elem)
 {
@@ -114,25 +129,40 @@ template <> void QuokkaSimulation<ResampledCoolingTest>::setInitialConditionsOnG
 	// Compute initial internal energy from temperature
 	const double k_B = C::k_B;
 	const double gamma = quokka::EOS_Traits<ResampledCoolingTest>::gamma;
-	const double mu = 0.6 * quokka::EOS_Traits<ResampledCoolingTest>::mean_molecular_weight;
 
 	// For ideal gas: P = (gamma - 1) * rho * e_int
 	// and P = rho * k_B * T / (mu * m_u)
 	// Therefore: e_int = k_B * T / ((gamma - 1) * mu * m_u)
-	const double e_int_initial = k_B * T_initial / ((gamma - 1.0) * mu);
+	const double e_int_initial = k_B * T_initial / ((gamma - 1.0) * mu_initial);
 	const double Eint_initial = rho_initial * e_int_initial;
+	const double Egas_initial = Eint_initial + active_magnetic_energy_initial;
 
 	// loop over the grid and set the initial condition
 	amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
 		state_cc(i, j, k, HydroSystem<ResampledCoolingTest>::density_index) = rho_initial;
 		state_cc(i, j, k, HydroSystem<ResampledCoolingTest>::x1Momentum_index) = 0.;
 		state_cc(i, j, k, HydroSystem<ResampledCoolingTest>::x2Momentum_index) = 0.;
 		state_cc(i, j, k, HydroSystem<ResampledCoolingTest>::x3Momentum_index) = 0.;
-		state_cc(i, j, k, HydroSystem<ResampledCoolingTest>::energy_index) = Eint_initial;
+		state_cc(i, j, k, HydroSystem<ResampledCoolingTest>::energy_index) = Egas_initial;
 		state_cc(i, j, k, HydroSystem<ResampledCoolingTest>::internalEnergy_index) = Eint_initial;
 	});
 }
 
+template <> void QuokkaSimulation<ResampledCoolingTest>::setInitialConditionsOnGridFaceVars(quokka::grid const &grid_elem)
+{
+	const amrex::Array4<double> &state_fc = grid_elem.array_;
+	const amrex::Box &indexRange = grid_elem.indexRange_;
+	const quokka::direction dir = grid_elem.dir_;
+
+	amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
+		if (dir == quokka::direction::x) {
+			state_fc(i, j, k, Physics_Indices<ResampledCoolingTest>::mhdFirstIndex) = Bx_initial;
+		} else {
+			state_fc(i, j, k, Physics_Indices<ResampledCoolingTest>::mhdFirstIndex) = 0.0;
+		}
+	});
+}
+
 template <> void QuokkaSimulation<ResampledCoolingTest>::computeAfterTimestep()
 {
 	// Extract solution at center of domain
@@ -143,8 +173,8 @@ template <> void QuokkaSimulation<ResampledCoolingTest>::computeAfterTimestep()
 
 		const amrex::Real Etot = values.at(HydroSystem<ResampledCoolingTest>::energy_index)[0];
 		const amrex::Real rho = values.at(HydroSystem<ResampledCoolingTest>::density_index)[0];
-		// For isochoric cooling with no kinetic energy, Eint = Etot
-		const amrex::Real Eint = Etot;
+		// For isochoric MHD cooling with no kinetic energy and a uniform magnetic field, subtract the constant magnetic energy.
+		const amrex::Real Eint = Etot - active_magnetic_energy_initial;
 
 		// Get temperature from tables
 		amrex::Real T = NAN;