Isotropic constant thermal conduction

inputs/ThermalConduction.toml

@@ -0,0 +1,33 @@
+# *****************************************************************
+# Problem size and geometry
+# *****************************************************************
+geometry.prob_lo = [-1.5428e18, -1.5428e18, -1.5428e18]
+geometry.prob_hi = [1.5428e18, 1.5428e18, 1.5428e18]
+quokka.bc = ["foextrap", "foextrap", "foextrap"]
+
+# *****************************************************************
+# VERBOSITY
+# *****************************************************************
+amr.v = 1      # verbosity in Amr
+
+# *****************************************************************
+# Resolution and refinement
+# *****************************************************************
+amr.n_cell = [128, 128, 128]
+amr.max_level = 0  # number of levels = max_level + 1
+amr.blocking_factor = 32  # grid size must be divisible by this
+
+cfl = 0.3
+
+cooling.enabled = 0
+cooling.read_tables_even_if_disabled = 1
+cooling.cooling_table_type = "resampled"
+cooling.hdf5_data_file = "../extern/cooling/CloudyData_UVB=HM2012_resampled.h5"
+
+
+conduction.enabled = 1
+conduction.conductivity_prefactor = 0.6113916490935586e12
+conduction.flux_limiter_phi = 1
+conduction.saturation_factor = 1.2e20
+conduction.conduction_cfl = 0.1
+conduction.eos_flag = 1

scripts/hpc_profiles/gadi_hopper.profile

@@ -1,4 +1,5 @@
-module load gcc/12.2.0
+module purge
+module load gcc/14.2.0
 module load cuda/12.9.0
 
 module unload openmpi

src/QuokkaSimulation.hpp

@@ -68,6 +68,7 @@ namespace filesystem = experimental::filesystem;
 
 #include "SimulationData.hpp"
 #include "chemistry/Chemistry.hpp"
+#include "conduction/ElectronConduction.hpp"
 #include "cooling/ResampledCooling.hpp"
 #include "dust/DustDrag.hpp"
 #include "dust/dust_system.hpp"
@@ -141,6 +142,10 @@ template <typename problem_t> class QuokkaSimulation : public AMRSimulation<prob
 	using AMRSimulation<problem_t>::sfh_interval_;
 	using AMRSimulation<problem_t>::sfh_time_interval_;
 
+	using AMRSimulation<problem_t>::enableElectronConduction_;
+	using AMRSimulation<problem_t>::electronConductionKappa0_;
+	using AMRSimulation<problem_t>::conductionCFL;
+
 #if AMREX_SPACEDIM == 3
 	using AMRSimulation<problem_t>::luminosityTables_;
 #endif // AMREX_SPACEDIM == 3
@@ -167,6 +172,10 @@ template <typename problem_t> class QuokkaSimulation : public AMRSimulation<prob
 	std::string coolingTableType_;
 	std::string coolingTableFilename_;
 
+	amrex::Real electronConductionFluxLimiterPhi_ = 1.0;
+	amrex::Real electronConductionSaturationFactor_ = 5.0;
+	int eosFlagForElectronConduction_ = 1; // 1 == use quokka::EOS; 0 == use resampled cooling
+
 	std::map<std::string, std::string> turbParams_;
 
 	static constexpr int nvarTotal_cc_ = Physics_Indices<problem_t>::nvarTotal_cc;
@@ -616,6 +625,17 @@ template <typename problem_t> void QuokkaSimulation<problem_t>::readParmParse()
 		}
 	}
 
+	// set electron thermal conduction runtime parameters
+	{
+		amrex::ParmParse const hpp("conduction");
+		hpp.query("enabled", enableElectronConduction_);
+		hpp.query("conductivity_prefactor", electronConductionKappa0_);
+		hpp.query("conduction_cfl", conductionCFL);
+		hpp.query("flux_limiter_phi", electronConductionFluxLimiterPhi_);
+		hpp.query("saturation_factor", electronConductionSaturationFactor_);
+		hpp.query("eos_flag", eosFlagForElectronConduction_);
+	}
+
 	// set turbulence runtime parameters
 	{
 		amrex::ParmParse const hpp("turbulence");
@@ -1032,6 +1052,19 @@ auto QuokkaSimulation<problem_t>::addStrangSplitSourcesWithBuiltin(amrex::MultiF
 		}
 	};
 
+	if (enableElectronConduction_ == 1) {
+		if (max_level > 0) {
+			amrex::Abort("Electron conduction not implemented for > 0 levels.");
+		}
+		fillBoundaryConditions(state, state, lev, time, quokka::centering::cc, quokka::direction::na, PreInterpState, PostInterpState);
+		const quokka::conduction::ElectronConductionParams conduction_params{.conductivity_prefactor = electronConductionKappa0_,
+										     .flux_limiter_phi = electronConductionFluxLimiterPhi_,
+										     .saturation_factor = electronConductionSaturationFactor_,
+										     .min_temperature = tempFloor_,
+										     .eos_flag = eosFlagForElectronConduction_};
+		quokka::conduction::ElectronConduction<problem_t>::ComputeExplicit(state, state_fc, geom[lev], dt, conduction_params, resampledTables_);
+	}
+
 	auto const applyUserSources = [&]() {
 		// compute user-specified sources
 		addStrangSplitSources(state, lev, time, dt);

src/conduction/ElectronConduction.hpp

@@ -0,0 +1,214 @@
+#ifndef ELECTRON_CONDUCTION_HPP_ // NOLINT
+#define ELECTRON_CONDUCTION_HPP_
+
+//==============================================================================
+// TwoMomentRad - a radiation transport library for patch-based AMR codes
+// Copyright 2020 Benjamin Wibking.
+// Released under the MIT license. See LICENSE file included in the GitHub repo.
+//==============================================================================
+/// \file ElectronConduction.hpp
+/// \brief Explicit flux-limited electron thermal conduction update.
+
+#include <array>
+#include <cmath>
+#include <limits>
+
+#include "AMReX_Array4.H"
+#include "AMReX_Geometry.H"
+#include "AMReX_GpuQualifiers.H"
+#include "AMReX_MultiFab.H"
+#include "AMReX_REAL.H"
+#include "AMReX_SPACE.H"
+#include "AMReX_Vector.H"
+#include "cooling/ResampledCooling.hpp"
+#include "hydro/hydro_system.hpp"
+
+namespace quokka::conduction
+{
+
+struct ElectronConductionParams {
+	amrex::Real conductivity_prefactor = 3.e34; // units of erg cm^-1 s^-1 K^-1
+	amrex::Real flux_limiter_phi = 0.1;
+	amrex::Real saturation_factor = 5.0; // refer to equation 8 of Cowee & McKee 1977
+	amrex::Real min_temperature = 0.0;   // default value will be overwritten by tempFloor_ during initialization
+	int eos_flag = 1;		     // 1 == use quokka::EOS; 0 == use resampled cooling
+};
+
+template <typename problem_t> class ElectronConduction
+{
+      public:
+	static void ComputeExplicit(amrex::MultiFab &state, std::array<amrex::MultiFab, AMREX_SPACEDIM> const &state_fc, amrex::Geometry const &geom,
+				    amrex::Real dt, ElectronConductionParams const &params, const quokka::ResampledCooling::resampled_tables &tables)
+	{
+		static_assert(Physics_Traits<problem_t>::is_hydro_enabled, "Electron conduction requires hydro to be enabled.");
+
+		if ((dt <= 0.0) || (params.conductivity_prefactor <= 0.0)) {
+			return;
+		}
+
+		if constexpr (HydroSystem<problem_t>::is_eos_isothermal()) {
+			amrex::ignore_unused(state_fc, geom, params);
+			return;
+		}
+
+		AMREX_ALWAYS_ASSERT_WITH_MESSAGE(state.nGrow() >= 1, "Electron conduction requires at least 1 ghost cell.");
+
+		const auto dx = geom.CellSizeArray();
+		const amrex::Real flux_limiter_phi = params.flux_limiter_phi;
+		const amrex::Real saturation_factor = params.saturation_factor;
+		const amrex::Real t_min = params.min_temperature;
+		const amrex::Real small = std::numeric_limits<amrex::Real>::min();
+
+		amrex::MultiFab temperature(state.boxArray(), state.DistributionMap(), 1, state.nGrow());
+		temperature.setVal(0.0);
+		amrex::MultiFab conductivity(state.boxArray(), state.DistributionMap(), 1, state.nGrow());
+		conductivity.setVal(0.0);
+		amrex::MultiFab saturated_flux(state.boxArray(), state.DistributionMap(), 1, state.nGrow());
+		saturated_flux.setVal(0.0);
+
+		auto const &state_x0 = state.const_arrays();
+		auto const &state_fc_x0 = state_fc[0].const_arrays();
+#if AMREX_SPACEDIM >= 2
+		auto const &state_fc_x1 = state_fc[1].const_arrays();
+#endif
+#if AMREX_SPACEDIM == 3
+		auto const &state_fc_x2 = state_fc[2].const_arrays();
+#endif
+		auto temperature_arr = temperature.arrays();
+		auto conductivity_arr = conductivity.arrays();
+		auto saturated_flux_arr = saturated_flux.arrays();
+		amrex::IntVect ng = amrex::IntVect(AMREX_D_DECL(2, 2, 2));
+		std::optional<decltype(tables.const_tables())> tables_dev;
+		if (params.eos_flag == 0) {
+			tables_dev = tables.const_tables();
+		} else if (params.eos_flag != 0 && params.eos_flag != 1) {
+			amrex::Abort("Invalid eos_flag value in ElectronConduction. Must be 0 (resampled cooling) or 1 (quokka::EOS).");
+		}
+
+		amrex::ParallelFor(state, ng, [=] AMREX_GPU_DEVICE(int bx, int i, int j, int k) noexcept {
+			std::array<amrex::Array4<const amrex::Real>, AMREX_SPACEDIM> local_state_fc{};
+			if (Physics_Traits<problem_t>::is_mhd_enabled) {
+				local_state_fc[0] = state_fc_x0[bx];
+#if AMREX_SPACEDIM >= 2
+				local_state_fc[1] = state_fc_x1[bx];
+#endif
+#if AMREX_SPACEDIM == 3
+				local_state_fc[2] = state_fc_x2[bx];
+#endif
+			}
+
+			auto const &cons = state_x0[bx];
+			const amrex::Real rho = cons(i, j, k, HydroSystem<problem_t>::density_index);
+			const amrex::Real Eint = HydroSystem<problem_t>::ComputeInternalEnergy(cons, i, j, k, &local_state_fc);
+			amrex::Real Tgas = NAN;
+			amrex::Real cs = NAN;
+			if (params.eos_flag == 0) {
+				Tgas = quokka::ResampledCooling::ComputeTgasFromEgas(rho, Eint, *tables_dev);
+				cs = quokka::ResampledCooling::ComputeSoundSpeedFromRhoEint(rho, Eint, *tables_dev);
+			} else if (params.eos_flag == 1) {
+				Tgas = quokka::EOS<problem_t>::ComputeTgasFromEint(rho, Eint);
+				amrex::Real Pgas = quokka::EOS<problem_t>::ComputePressure(rho, Eint);
+				cs = quokka::EOS<problem_t>::ComputeSoundSpeed(rho, Pgas);
+			}
+
+			const amrex::Real Tuse = amrex::max(Tgas, t_min);
+			const amrex::Real kappa = params.conductivity_prefactor;
+			const amrex::Real qsat = amrex::max(saturation_factor * flux_limiter_phi * rho * cs * cs * cs, small);
+
+			temperature_arr[bx](i, j, k) = Tuse;
+			conductivity_arr[bx](i, j, k) = kappa;
+			saturated_flux_arr[bx](i, j, k) = qsat;
+		});
+
+		std::array<amrex::MultiFab, AMREX_SPACEDIM> heat_flux;
+		for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+			amrex::BoxArray const ba_face = amrex::convert(state.boxArray(), amrex::IntVect::TheDimensionVector(idim));
+			heat_flux[idim].define(ba_face, state.DistributionMap(), 1, 0);
+			heat_flux[idim].setVal(0.0);
+		}
+
+		auto const &temp = temperature.const_arrays();
+		auto const &kappa = conductivity.const_arrays();
+		auto const &qsat = saturated_flux.const_arrays();
+		auto flux_x = heat_flux[0].arrays();
+		amrex::ParallelFor(heat_flux[0], [=] AMREX_GPU_DEVICE(int bx, int i, int j, int k) noexcept {
+			const amrex::Real gradT = (temp[bx](i, j, k) - temp[bx](i - 1, j, k)) / dx[0];
+			const amrex::Real kappa_face = 0.5 * (kappa[bx](i, j, k) + kappa[bx](i - 1, j, k));
+			const amrex::Real q_classical = -kappa_face * gradT;
+			const amrex::Real q_sat_face = 0.5 * (qsat[bx](i, j, k) + qsat[bx](i - 1, j, k));
+			const amrex::Real limiter = 1.0 + std::abs(q_classical) / amrex::max(q_sat_face, small);
+			flux_x[bx](i, j, k) = q_classical / limiter;
+		});
+
+#if AMREX_SPACEDIM >= 2
+		auto flux_y = heat_flux[1].arrays();
+		amrex::ParallelFor(heat_flux[1], [=] AMREX_GPU_DEVICE(int bx, int i, int j, int k) noexcept {
+			const amrex::Real gradT = (temp[bx](i, j, k) - temp[bx](i, j - 1, k)) / dx[1];
+			const amrex::Real kappa_face = 0.5 * (kappa[bx](i, j, k) + kappa[bx](i, j - 1, k));
+			const amrex::Real q_classical = -kappa_face * gradT;
+			const amrex::Real q_sat_face = 0.5 * (qsat[bx](i, j, k) + qsat[bx](i, j - 1, k));
+			const amrex::Real limiter = 1.0 + std::abs(q_classical) / amrex::max(q_sat_face, small);
+			flux_y[bx](i, j, k) = q_classical / limiter;
+		});
+#endif
+
+#if AMREX_SPACEDIM == 3
+		auto flux_z = heat_flux[2].arrays();
+		amrex::ParallelFor(heat_flux[2], [=] AMREX_GPU_DEVICE(int bx, int i, int j, int k) noexcept {
+			const amrex::Real gradT = (temp[bx](i, j, k) - temp[bx](i, j, k - 1)) / dx[2];
+			const amrex::Real kappa_face = 0.5 * (kappa[bx](i, j, k) + kappa[bx](i, j, k - 1));
+			const amrex::Real q_classical = -kappa_face * gradT;
+			const amrex::Real q_sat_face = 0.5 * (qsat[bx](i, j, k) + qsat[bx](i, j, k - 1));
+			const amrex::Real limiter = 1.0 + std::abs(q_classical) / amrex::max(q_sat_face, small);
+			flux_z[bx](i, j, k) = q_classical / limiter;
+		});
+#endif
+
+		auto state_out = state.arrays();
+		auto const &flux_x_const = heat_flux[0].const_arrays();
+#if AMREX_SPACEDIM >= 2
+		auto const &flux_y_const = heat_flux[1].const_arrays();
+#endif
+#if AMREX_SPACEDIM == 3
+		auto const &flux_z_const = heat_flux[2].const_arrays();
+#endif
+
+		amrex::ParallelFor(state, [=] AMREX_GPU_DEVICE(int bx, int i, int j, int k) noexcept {
+			std::array<amrex::Array4<const amrex::Real>, AMREX_SPACEDIM> local_state_fc{};
+			if (Physics_Traits<problem_t>::is_mhd_enabled) {
+				local_state_fc[0] = state_fc_x0[bx];
+#if AMREX_SPACEDIM >= 2
+				local_state_fc[1] = state_fc_x1[bx];
+#endif
+#if AMREX_SPACEDIM == 3
+				local_state_fc[2] = state_fc_x2[bx];
+#endif
+			}
+
+			const amrex::Real rho = state_out[bx](i, j, k, HydroSystem<problem_t>::density_index);
+			const amrex::Real px = state_out[bx](i, j, k, HydroSystem<problem_t>::x1Momentum_index);
+			const amrex::Real py = state_out[bx](i, j, k, HydroSystem<problem_t>::x2Momentum_index);
+			const amrex::Real pz = state_out[bx](i, j, k, HydroSystem<problem_t>::x3Momentum_index);
+
+			const amrex::Real Ekin = 0.5 * (px * px + py * py + pz * pz) / rho;
+			const amrex::Real Eint_old = HydroSystem<problem_t>::ComputeInternalEnergy(state_out[bx], i, j, k, &local_state_fc);
+			const amrex::Real Emag = HydroSystem<problem_t>::ComputeMagneticEnergy(i, j, k, &local_state_fc);
+			amrex::Real div_flux = (flux_x_const[bx](i + 1, j, k) - flux_x_const[bx](i, j, k)) / dx[0];
+#if AMREX_SPACEDIM >= 2
+			div_flux += (flux_y_const[bx](i, j + 1, k) - flux_y_const[bx](i, j, k)) / dx[1];
+#endif
+#if AMREX_SPACEDIM == 3
+			div_flux += (flux_z_const[bx](i, j, k + 1) - flux_z_const[bx](i, j, k)) / dx[2];
+#endif
+
+			amrex::Real Eint_new = Eint_old - dt * div_flux;
+
+			state_out[bx](i, j, k, HydroSystem<problem_t>::energy_index) = Eint_new + Ekin + Emag;
+			state_out[bx](i, j, k, HydroSystem<problem_t>::internalEnergy_index) = Eint_new;
+		});
+	}
+};
+
+} // namespace quokka::conduction
+
+#endif // ELECTRON_CONDUCTION_HPP_

src/problems/ThermalConduction/CMakeLists.txt

@@ -0,0 +1 @@
+quokka_add_problem(JOB_NAME ThermalConduction PRIORITY 100)