photoelectric heating

src/QuokkaSimulation.hpp

@@ -1617,10 +1617,10 @@ void QuokkaSimulation<problem_t>::subcycleRadiationAtLevel(int lev, amrex::Real
 		// new radiation state is stored in state_new_cc_
 		// new hydro state is stored in state_new_cc_ (always the case during radiation update)
 
-		amrex::Gpu::Buffer<int> iteration_failure_counter(
-		    {0, 0, 0}); // failure counter for: matter-radiation coupling, dust temperature, outer iteration
-		amrex::Gpu::Buffer<int> iteration_counter(
-		    {0, 0, 0}); // iteration counter for: radiation update, Newton-Raphson iterations, max Newton-Raphson iterations
+		// failure counter for: matter-radiation coupling, dust temperature, outer iteration
+		amrex::Gpu::Buffer<int> iteration_failure_counter({0, 0, 0});
+		// iteration counter for: radiation update, Newton-Raphson iterations, max Newton-Raphson iterations, decoupled gas-dust update
+		amrex::Gpu::Buffer<int> iteration_counter({0, 0, 0, 0});
 		int *p_iteration_failure_counter = iteration_failure_counter.data();
 		int *p_iteration_counter = iteration_counter.data();
 
@@ -1659,13 +1659,15 @@ void QuokkaSimulation<problem_t>::subcycleRadiationAtLevel(int lev, amrex::Real
 
 		if (print_rad_counter_) {
 			auto h_iteration_counter = iteration_counter.copyToHost();
-			long global_solver_count = h_iteration_counter[0];  // number of Newton-Raphson solvings
-			long global_iteration_sum = h_iteration_counter[1]; // sum of Newton-Raphson iterations
-			int global_iteration_max = h_iteration_counter[2];  // max number of Newton-Raphson iterations
+			long global_solver_count = h_iteration_counter[0];	      // number of Newton-Raphson solvings
+			long global_iteration_sum = h_iteration_counter[1];	      // sum of Newton-Raphson iterations
+			int global_iteration_max = h_iteration_counter[2];	      // max number of Newton-Raphson iterations
+			long global_decoupled_iteration_sum = h_iteration_counter[3]; // sum of decoupled gas-dust Newton-Raphson iterations
 
 			amrex::ParallelDescriptor::ReduceLongSum(global_solver_count);
 			amrex::ParallelDescriptor::ReduceLongSum(global_iteration_sum);
 			amrex::ParallelDescriptor::ReduceIntMax(global_iteration_max);
+			amrex::ParallelDescriptor::ReduceLongSum(global_decoupled_iteration_sum);
 
 			if (amrex::ParallelDescriptor::IOProcessor()) {
 				const auto n_cells = CountCells(lev);
@@ -1675,9 +1677,15 @@ void QuokkaSimulation<problem_t>::subcycleRadiationAtLevel(int lev, amrex::Real
 					    static_cast<double>(global_iteration_sum) / static_cast<double>(global_solver_count);
 					const double global_solving_mean =
 					    static_cast<double>(global_solver_count) / static_cast<double>(n_cells) / 2.0; // 2 stages
-					amrex::Print() << "average number of Newton-Raphson solvings per IMEX stage is " << global_solving_mean
+					const double global_decoupled_iteration_mean =
+					    static_cast<double>(global_decoupled_iteration_sum) / static_cast<double>(global_solver_count);
+					amrex::Print() << "The average number of Newton-Raphson solvings per IMEX stage is " << global_solving_mean
 						       << ", (mean, max) number of Newton-Raphson iterations are " << global_iteration_mean << ", "
-						       << global_iteration_max << "\n";
+						       << global_iteration_max << ".\n";
+					if constexpr (RadSystem_Traits<problem_t>::enable_dust_gas_thermal_coupling_model) {
+						amrex::Print() << "The fraction of gas-dust interactions that are decoupled is "
+							       << global_decoupled_iteration_mean << "\n";
+					}
 				}
 			}
 		}

src/problems/CMakeLists.txt

@@ -36,6 +36,7 @@ add_subdirectory(RadTube)
 add_subdirectory(RadDust)
 add_subdirectory(RadDustMG)
 add_subdirectory(RadMarshakDust)
+add_subdirectory(RadMarshakDustPE)
 
 add_subdirectory(RadhydroShell)
 add_subdirectory(RadhydroShock)

src/problems/RadDustMG/test_rad_dust_MG.hpp

@@ -14,7 +14,7 @@
 // internal headers
 
 #include "math/interpolate.hpp"
-#include "radiation/radiation_system.hpp"
+#include "radiation/radiation_dust_system.hpp"
 
 // function definitions
 

src/problems/RadMarshakDust/test_radiation_marshak_dust.cpp

@@ -13,7 +13,7 @@
 #include "util/fextract.hpp"
 #include "util/valarray.hpp"
 
-struct StreamingProblem {
+struct MarshakProblem {
 };
 
 AMREX_GPU_MANAGED double kappa1 = NAN; // dust opacity at IR
@@ -42,13 +42,13 @@ constexpr int n_group_ = 2;
 static constexpr amrex::GpuArray<double, n_group_ + 1> radBoundaries_{1e-10, 100, 1e4};
 static constexpr OpacityModel opacity_model_ = OpacityModel::piecewise_constant_opacity;
 
-template <> struct quokka::EOS_Traits<StreamingProblem> {
+template <> struct quokka::EOS_Traits<MarshakProblem> {
 	static constexpr double mean_molecular_weight = mu;
 	static constexpr double boltzmann_constant = 1.0;
 	static constexpr double gamma = 5. / 3.;
 };
 
-template <> struct Physics_Traits<StreamingProblem> {
+template <> struct Physics_Traits<MarshakProblem> {
 	// cell-centred
 	static constexpr bool is_hydro_enabled = false;
 	static constexpr int numMassScalars = 0;		     // number of mass scalars
@@ -59,7 +59,7 @@ template <> struct Physics_Traits<StreamingProblem> {
 	static constexpr int nGroups = n_group_; // number of radiation groups
 };
 
-template <> struct RadSystem_Traits<StreamingProblem> {
+template <> struct RadSystem_Traits<MarshakProblem> {
 	static constexpr double c_light = c;
 	static constexpr double c_hat = chat;
 	static constexpr double radiation_constant = a_rad;
@@ -71,20 +71,25 @@ template <> struct RadSystem_Traits<StreamingProblem> {
 	static constexpr OpacityModel opacity_model = opacity_model_;
 };
 
-template <> AMREX_GPU_HOST_DEVICE auto RadSystem<StreamingProblem>::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real
+template <> struct ISM_Traits<MarshakProblem> {
+	static constexpr double gas_dust_coupling_threshold = 1.0e-5;
+	static constexpr bool enable_photoelectric_heating = false;
+};
+
+template <> AMREX_GPU_HOST_DEVICE auto RadSystem<MarshakProblem>::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real
 {
 	return kappa1;
 }
 
-template <> AMREX_GPU_HOST_DEVICE auto RadSystem<StreamingProblem>::ComputeFluxMeanOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real
+template <> AMREX_GPU_HOST_DEVICE auto RadSystem<MarshakProblem>::ComputeFluxMeanOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real
 {
 	return kappa1;
 }
 
 template <>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto
-RadSystem<StreamingProblem>::DefineOpacityExponentsAndLowerValues(amrex::GpuArray<double, nGroups_ + 1> /*rad_boundaries*/, const double /*rho*/,
-								  const double /*Tgas*/) -> amrex::GpuArray<amrex::GpuArray<double, nGroups_ + 1>, 2>
+RadSystem<MarshakProblem>::DefineOpacityExponentsAndLowerValues(amrex::GpuArray<double, nGroups_ + 1> /*rad_boundaries*/, const double /*rho*/,
+								const double /*Tgas*/) -> amrex::GpuArray<amrex::GpuArray<double, nGroups_ + 1>, 2>
 {
 	amrex::GpuArray<amrex::GpuArray<double, nGroups_ + 1>, 2> exponents_and_values{};
 	for (int i = 0; i < nGroups_ + 1; ++i) {
@@ -98,36 +103,36 @@ RadSystem<StreamingProblem>::DefineOpacityExponentsAndLowerValues(amrex::GpuArra
 	return exponents_and_values;
 }
 
-template <> void QuokkaSimulation<StreamingProblem>::setInitialConditionsOnGrid(quokka::grid const &grid_elem)
+template <> void QuokkaSimulation<MarshakProblem>::setInitialConditionsOnGrid(quokka::grid const &grid_elem)
 {
 	const amrex::Box &indexRange = grid_elem.indexRange_;
 	const amrex::Array4<double> &state_cc = grid_elem.array_;
 
 	const auto Egas0 = initial_T * CV;
-	const auto Erads = RadSystem<StreamingProblem>::ComputeThermalRadiationMultiGroup(initial_Trad, radBoundaries_);
+	const auto Erads = RadSystem<MarshakProblem>::ComputeThermalRadiationMultiGroup(initial_Trad, radBoundaries_);
 
 	// loop over the grid and set the initial condition
 	amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
-		for (int g = 0; g < Physics_Traits<StreamingProblem>::nGroups; ++g) {
-			state_cc(i, j, k, RadSystem<StreamingProblem>::radEnergy_index + Physics_NumVars::numRadVars * g) = Erads[g];
-			state_cc(i, j, k, RadSystem<StreamingProblem>::x1RadFlux_index + Physics_NumVars::numRadVars * g) = 0;
-			state_cc(i, j, k, RadSystem<StreamingProblem>::x2RadFlux_index + Physics_NumVars::numRadVars * g) = 0;
-			state_cc(i, j, k, RadSystem<StreamingProblem>::x3RadFlux_index + Physics_NumVars::numRadVars * g) = 0;
+		for (int g = 0; g < Physics_Traits<MarshakProblem>::nGroups; ++g) {
+			state_cc(i, j, k, RadSystem<MarshakProblem>::radEnergy_index + Physics_NumVars::numRadVars * g) = Erads[g];
+			state_cc(i, j, k, RadSystem<MarshakProblem>::x1RadFlux_index + Physics_NumVars::numRadVars * g) = 0;
+			state_cc(i, j, k, RadSystem<MarshakProblem>::x2RadFlux_index + Physics_NumVars::numRadVars * g) = 0;
+			state_cc(i, j, k, RadSystem<MarshakProblem>::x3RadFlux_index + Physics_NumVars::numRadVars * g) = 0;
 		}
-		state_cc(i, j, k, RadSystem<StreamingProblem>::gasEnergy_index) = Egas0;
-		state_cc(i, j, k, RadSystem<StreamingProblem>::gasDensity_index) = rho0;
-		state_cc(i, j, k, RadSystem<StreamingProblem>::gasInternalEnergy_index) = Egas0;
-		state_cc(i, j, k, RadSystem<StreamingProblem>::x1GasMomentum_index) = 0.;
-		state_cc(i, j, k, RadSystem<StreamingProblem>::x2GasMomentum_index) = 0.;
-		state_cc(i, j, k, RadSystem<StreamingProblem>::x3GasMomentum_index) = 0.;
+		state_cc(i, j, k, RadSystem<MarshakProblem>::gasEnergy_index) = Egas0;
+		state_cc(i, j, k, RadSystem<MarshakProblem>::gasDensity_index) = rho0;
+		state_cc(i, j, k, RadSystem<MarshakProblem>::gasInternalEnergy_index) = Egas0;
+		state_cc(i, j, k, RadSystem<MarshakProblem>::x1GasMomentum_index) = 0.;
+		state_cc(i, j, k, RadSystem<MarshakProblem>::x2GasMomentum_index) = 0.;
+		state_cc(i, j, k, RadSystem<MarshakProblem>::x3GasMomentum_index) = 0.;
 	});
 }
 
 template <>
 AMREX_GPU_DEVICE AMREX_FORCE_INLINE void
-AMRSimulation<StreamingProblem>::setCustomBoundaryConditions(const amrex::IntVect &iv, amrex::Array4<amrex::Real> const &consVar, int /*dcomp*/,
-							     int /*numcomp*/, amrex::GeometryData const &geom, const amrex::Real /*time*/,
-							     const amrex::BCRec * /*bcr*/, int /*bcomp*/, int /*orig_comp*/)
+AMRSimulation<MarshakProblem>::setCustomBoundaryConditions(const amrex::IntVect &iv, amrex::Array4<amrex::Real> const &consVar, int /*dcomp*/, int /*numcomp*/,
+							   amrex::GeometryData const &geom, const amrex::Real /*time*/, const amrex::BCRec * /*bcr*/,
+							   int /*bcomp*/, int /*orig_comp*/)
 {
 #if (AMREX_SPACEDIM == 1)
 	auto i = iv.toArray()[0];
@@ -145,7 +150,7 @@ AMRSimulation<StreamingProblem>::setCustomBoundaryConditions(const amrex::IntVec
 	amrex::Box const &box = geom.Domain();
 	amrex::GpuArray<int, 3> lo = box.loVect3d();
 
-	// const auto Erads = RadSystem<StreamingProblem>::ComputeThermalRadiation(T_rad_L, radBoundaries_);
+	// const auto Erads = RadSystem<MarshakProblem>::ComputeThermalRadiation(T_rad_L, radBoundaries_);
 	quokka::valarray<double, 2> const Erads = {erad_floor, EradL};
 	const double c_light = c;
 	const auto Frads = Erads * c_light;
@@ -154,22 +159,22 @@ AMRSimulation<StreamingProblem>::setCustomBoundaryConditions(const amrex::IntVec
 		// streaming inflow boundary
 		// multigroup radiation
 		// x1 left side boundary (Marshak)
-		for (int g = 0; g < Physics_Traits<StreamingProblem>::nGroups; ++g) {
-			consVar(i, j, k, RadSystem<StreamingProblem>::radEnergy_index + Physics_NumVars::numRadVars * g) = Erads[g];
-			consVar(i, j, k, RadSystem<StreamingProblem>::x1RadFlux_index + Physics_NumVars::numRadVars * g) = Frads[g];
-			consVar(i, j, k, RadSystem<StreamingProblem>::x2RadFlux_index + Physics_NumVars::numRadVars * g) = 0;
-			consVar(i, j, k, RadSystem<StreamingProblem>::x3RadFlux_index + Physics_NumVars::numRadVars * g) = 0;
+		for (int g = 0; g < Physics_Traits<MarshakProblem>::nGroups; ++g) {
+			consVar(i, j, k, RadSystem<MarshakProblem>::radEnergy_index + Physics_NumVars::numRadVars * g) = Erads[g];
+			consVar(i, j, k, RadSystem<MarshakProblem>::x1RadFlux_index + Physics_NumVars::numRadVars * g) = Frads[g];
+			consVar(i, j, k, RadSystem<MarshakProblem>::x2RadFlux_index + Physics_NumVars::numRadVars * g) = 0;
+			consVar(i, j, k, RadSystem<MarshakProblem>::x3RadFlux_index + Physics_NumVars::numRadVars * g) = 0;
 		}
 	}
 
 	// gas boundary conditions are the same everywhere
 	const double Egas = initial_T * CV;
-	consVar(i, j, k, RadSystem<StreamingProblem>::gasEnergy_index) = Egas;
-	consVar(i, j, k, RadSystem<StreamingProblem>::gasDensity_index) = rho0;
-	consVar(i, j, k, RadSystem<StreamingProblem>::gasInternalEnergy_index) = Egas;
-	consVar(i, j, k, RadSystem<StreamingProblem>::x1GasMomentum_index) = 0.;
-	consVar(i, j, k, RadSystem<StreamingProblem>::x2GasMomentum_index) = 0.;
-	consVar(i, j, k, RadSystem<StreamingProblem>::x3GasMomentum_index) = 0.;
+	consVar(i, j, k, RadSystem<MarshakProblem>::gasEnergy_index) = Egas;
+	consVar(i, j, k, RadSystem<MarshakProblem>::gasDensity_index) = rho0;
+	consVar(i, j, k, RadSystem<MarshakProblem>::gasInternalEnergy_index) = Egas;
+	consVar(i, j, k, RadSystem<MarshakProblem>::x1GasMomentum_index) = 0.;
+	consVar(i, j, k, RadSystem<MarshakProblem>::x2GasMomentum_index) = 0.;
+	consVar(i, j, k, RadSystem<MarshakProblem>::x3GasMomentum_index) = 0.;
 }
 
 auto problem_main() -> int
@@ -187,7 +192,7 @@ auto problem_main() -> int
 	pp.query("kappa2", kappa2);
 
 	// Boundary conditions
-	constexpr int nvars = RadSystem<StreamingProblem>::nvar_;
+	constexpr int nvars = RadSystem<MarshakProblem>::nvar_;
 	amrex::Vector<amrex::BCRec> BCs_cc(nvars);
 	for (int n = 0; n < nvars; ++n) {
 		BCs_cc[n].setLo(0, amrex::BCType::ext_dir);  // Dirichlet x1
@@ -199,7 +204,7 @@ auto problem_main() -> int
 	}
 
 	// Problem initialization
-	QuokkaSimulation<StreamingProblem> sim(BCs_cc);
+	QuokkaSimulation<MarshakProblem> sim(BCs_cc);
 
 	sim.radiationReconstructionOrder_ = 3; // PPM
 	// sim.stopTime_ = tmax; // set with runtime parameters
@@ -230,14 +235,14 @@ auto problem_main() -> int
 	for (int i = 0; i < nx; ++i) {
 		amrex::Real const x = position[i];
 		xs.at(i) = x;
-		erad1.at(i) = values.at(RadSystem<StreamingProblem>::radEnergy_index + Physics_NumVars::numRadVars * 0)[i];
+		erad1.at(i) = values.at(RadSystem<MarshakProblem>::radEnergy_index + Physics_NumVars::numRadVars * 0)[i];
 		erad.at(i) = erad1.at(i);
 		if (n_group_ > 1) {
-			erad2.at(i) = values.at(RadSystem<StreamingProblem>::radEnergy_index + Physics_NumVars::numRadVars * 1)[i];
+			erad2.at(i) = values.at(RadSystem<MarshakProblem>::radEnergy_index + Physics_NumVars::numRadVars * 1)[i];
 			erad.at(i) += erad2.at(i);
 		}
-		const double e_gas = values.at(RadSystem<StreamingProblem>::gasInternalEnergy_index)[i];
-		T.at(i) = quokka::EOS<StreamingProblem>::ComputeTgasFromEint(rho0, e_gas);
+		const double e_gas = values.at(RadSystem<MarshakProblem>::gasInternalEnergy_index)[i];
+		T.at(i) = quokka::EOS<MarshakProblem>::ComputeTgasFromEint(rho0, e_gas);
 		T_exact.at(i) = T_end_exact;
 
 		erad2_exact.at(i) = x < sim.tNew_[0] ? EradL * std::exp(-x * rho0 * kappa2) : erad_floor;

src/problems/RadMarshakDust/test_radiation_marshak_dust.hpp

@@ -17,7 +17,7 @@
 
 // internal headers
 
-#include "radiation/radiation_system.hpp"
+#include "radiation/radiation_dust_system.hpp"
 
 // function definitions
 

src/problems/RadMarshakDustPE/CMakeLists.txt

@@ -0,0 +1,8 @@
+add_executable(test_radiation_marshak_dust_PE test_radiation_marshak_dust_and_PE.cpp ../../util/fextract.cpp ${QuokkaObjSources})
+
+if(AMReX_GPU_BACKEND MATCHES "CUDA")
+    setup_target_for_cuda_compilation(test_radiation_marshak_dust_PE)
+endif(AMReX_GPU_BACKEND MATCHES "CUDA")
+
+add_test(NAME RadiationMarshakDustPE-coupled COMMAND test_radiation_marshak_dust_PE RadMarshakDustPEcoupled.in ${QuokkaTestParams} WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests)
+add_test(NAME RadiationMarshakDustPE-decoupled COMMAND test_radiation_marshak_dust_PE RadMarshakDustPEdecoupled.in ${QuokkaTestParams} WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests)