Clean up code for single-group radiation

src/RadhydroSimulation.hpp

@@ -129,6 +129,8 @@ template <typename problem_t> class RadhydroSimulation : public AMRSimulation<pr
 	amrex::Real errorNorm_ = NAN;
 	amrex::Real pressureFloor_ = 0.;
 
+	bool use_wavespeed_correction_ = false;
+
 	int lowLevelDebuggingOutput_ = 0;	// 0 == do nothing; 1 == output intermediate multifabs used in hydro each timestep (ONLY USE FOR DEBUGGING)
 	int integratorOrder_ = 2;		// 1 == forward Euler; 2 == RK2-SSP (default)
 	int reconstructionOrder_ = 3;		// 1 == donor cell; 2 == PLM; 3 == PPM (default)
@@ -1882,8 +1884,7 @@ void RadhydroSimulation<problem_t>::fluxFunction(amrex::Array4<const amrex::Real
 	// interface-centered kernel
 	amrex::Box const &x1FluxRange = amrex::surroundingNodes(indexRange, dir);
 	RadSystem<problem_t>::template ComputeFluxes<DIR>(x1Flux.array(), x1FluxDiffusive.array(), x1LeftState.array(), x1RightState.array(), x1FluxRange,
-							  consState,
-							  dx); // watch out for argument order!!
+							  consState, dx, use_wavespeed_correction_); // watch out for argument order!!
 }
 
 #endif // RADIATION_SIMULATION_HPP_

src/problems/CMakeLists.txt

@@ -23,6 +23,7 @@ add_subdirectory(RadForce)
 add_subdirectory(RadMarshak)
 add_subdirectory(RadMarshakCGS)
 add_subdirectory(RadMarshakAsymptotic)
+add_subdirectory(RadMarshakAsymptoticCorr)
 add_subdirectory(RadMarshakVaytet)
 add_subdirectory(RadMatterCoupling)
 add_subdirectory(RadMatterCouplingRSLA)

src/problems/RadBeam/test_radiation_beam.cpp

@@ -51,20 +51,14 @@ template <> struct Physics_Traits<BeamProblem> {
 	static constexpr int nGroups = 1; // number of radiation groups
 };
 
-template <>
-AMREX_GPU_HOST_DEVICE auto RadSystem<BeamProblem>::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> quokka::valarray<double, nGroups_>
+template <> AMREX_GPU_HOST_DEVICE auto RadSystem<BeamProblem>::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real
 {
-	quokka::valarray<double, nGroups_> kappaPVec{};
-	for (int g = 0; g < nGroups_; ++g) {
-		kappaPVec[g] = kappa0;
-	}
-	return kappaPVec;
+	return kappa0;
 }
 
-template <>
-AMREX_GPU_HOST_DEVICE auto RadSystem<BeamProblem>::ComputeFluxMeanOpacity(const double /*rho*/, const double /*Tgas*/) -> quokka::valarray<double, nGroups_>
+template <> AMREX_GPU_HOST_DEVICE auto RadSystem<BeamProblem>::ComputeFluxMeanOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real
 {
-	return ComputePlanckOpacity(0., 0.);
+	return kappa0;
 }
 
 template <>

src/problems/RadForce/test_radiation_force.cpp

@@ -13,6 +13,7 @@
 #include "AMReX.H"
 #include "AMReX_BC_TYPES.H"
 #include "AMReX_BLassert.H"
+#include "AMReX_Box.H"
 #include "AMReX_REAL.H"
 
 #include "RadhydroSimulation.hpp"
@@ -30,7 +31,7 @@
 struct TubeProblem {
 };
 
-constexpr double kappa0 = 5.0;	     // cm^2 g^-1
+constexpr amrex::Real kappa0 = 5.0;  // cm^2 g^-1
 constexpr double mu = 2.33 * C::m_u; // g
 constexpr double gamma_gas = 1.0;    // isothermal gas EOS
 constexpr double a0 = 0.2e5;	     // cm s^-1
@@ -72,26 +73,11 @@ template <> struct RadSystem_Traits<TubeProblem> {
 	static constexpr int beta_order = 1;
 };
 
-template <>
-AMREX_GPU_HOST_DEVICE auto RadSystem<TubeProblem>::ComputePlanckOpacity(const double /*rho*/,
-									const double /*Tgas*/) -> quokka::valarray<double, Physics_Traits<TubeProblem>::nGroups>
-{
-	quokka::valarray<double, Physics_Traits<TubeProblem>::nGroups> kappaPVec{};
-	for (int g = 0; g < nGroups_; ++g) {
-		kappaPVec[g] = 0.; // no heating/cooling
-	}
-	return kappaPVec;
-}
+template <> AMREX_GPU_HOST_DEVICE auto RadSystem<TubeProblem>::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real { return 0.; }
 
-template <>
-AMREX_GPU_HOST_DEVICE auto
-RadSystem<TubeProblem>::ComputeFluxMeanOpacity(const double /*rho*/, const double /*Tgas*/) -> quokka::valarray<double, Physics_Traits<TubeProblem>::nGroups>
+template <> AMREX_GPU_HOST_DEVICE auto RadSystem<TubeProblem>::ComputeFluxMeanOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real
 {
-	quokka::valarray<double, Physics_Traits<TubeProblem>::nGroups> kappaFVec{};
-	for (int g = 0; g < nGroups_; ++g) {
-		kappaFVec[g] = kappa0;
-	}
-	return kappaFVec;
+	return kappa0;
 }
 
 template <> void RadhydroSimulation<TubeProblem>::setInitialConditionsOnGrid(quokka::grid grid_elem)

src/problems/RadMarshak/test_radiation_marshak.cpp

@@ -55,20 +55,14 @@ template <> struct Physics_Traits<SuOlsonProblem> {
 	static constexpr int nGroups = 1; // number of radiation groups
 };
 
-template <>
-AMREX_GPU_HOST_DEVICE auto RadSystem<SuOlsonProblem>::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> quokka::valarray<double, nGroups_>
+template <> AMREX_GPU_HOST_DEVICE auto RadSystem<SuOlsonProblem>::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real
 {
-	quokka::valarray<double, nGroups_> kappaPVec{};
-	for (int i = 0; i < nGroups_; ++i) {
-		kappaPVec[i] = kappa;
-	}
-	return kappaPVec;
+	return kappa;
 }
 
-template <>
-AMREX_GPU_HOST_DEVICE auto RadSystem<SuOlsonProblem>::ComputeFluxMeanOpacity(const double /*rho*/, const double /*Tgas*/) -> quokka::valarray<double, nGroups_>
+template <> AMREX_GPU_HOST_DEVICE auto RadSystem<SuOlsonProblem>::ComputeFluxMeanOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real
 {
-	return ComputePlanckOpacity(0.0, 0.0);
+	return kappa;
 }
 
 static constexpr int nmscalars_ = Physics_Traits<SuOlsonProblem>::numMassScalars;

src/problems/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp

@@ -52,21 +52,14 @@ template <> struct Physics_Traits<SuOlsonProblemCgs> {
 	static constexpr int nGroups = 1; // number of radiation groups
 };
 
-template <>
-AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto RadSystem<SuOlsonProblemCgs>::ComputePlanckOpacity(const double rho,
-												 const double Tgas) -> quokka::valarray<double, nGroups_>
+template <> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto RadSystem<SuOlsonProblemCgs>::ComputePlanckOpacity(const double rho, const double Tgas) -> amrex::Real
 {
 	auto sigma = kappa * std::pow(Tgas / T_hohlraum, -3); // cm^-1
-	quokka::valarray<double, nGroups_> kappaPVec{};
-	for (int i = 0; i < nGroups_; ++i) {
-		kappaPVec[i] = sigma / rho; // cm^2 g^-1
-	}
-	return kappaPVec;
+	return sigma / rho;				      // cm^2 g^-1
 }
 
 template <>
-AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto RadSystem<SuOlsonProblemCgs>::ComputeFluxMeanOpacity(const double rho,
-												   const double Tgas) -> quokka::valarray<double, nGroups_>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto RadSystem<SuOlsonProblemCgs>::ComputeFluxMeanOpacity(const double rho, const double Tgas) -> amrex::Real
 {
 	return ComputePlanckOpacity(rho, Tgas);
 }
@@ -244,6 +237,9 @@ auto problem_main() -> int
 	sim.maxTimesteps_ = max_timesteps;
 	sim.plotfileInterval_ = -1;
 
+	amrex::ParmParse pp("marshak");
+	pp.query("use_wavespeed_correction", sim.use_wavespeed_correction_);
+
 	// initialize
 	sim.setInitialConditions();
 

src/problems/RadMarshakAsymptoticCorr/CMakeLists.txt

@@ -0,0 +1,7 @@
+add_executable(test_radiation_marshak_asymptotic_corr ../RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp ../../util/fextract.cpp ../../math/interpolate.cpp ${QuokkaObjSources})
+
+if(AMReX_GPU_BACKEND MATCHES "CUDA")
+    setup_target_for_cuda_compilation(test_radiation_marshak_asymptotic_corr)
+endif(AMReX_GPU_BACKEND MATCHES "CUDA")
+
+add_test(NAME MarshakWaveAsymptoticDiffusion COMMAND test_radiation_marshak_asymptotic_corr MarshakAsymptoticCorr.in ${QuokkaTestParams} WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests)

src/problems/RadMarshakCGS/test_radiation_marshak_cgs.cpp

@@ -57,21 +57,14 @@ template <> struct Physics_Traits<SuOlsonProblemCgs> {
 	static constexpr int nGroups = 1; // number of radiation groups
 };
 
-template <>
-AMREX_GPU_HOST_DEVICE auto RadSystem<SuOlsonProblemCgs>::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> quokka::valarray<double, nGroups_>
+template <> AMREX_GPU_HOST_DEVICE auto RadSystem<SuOlsonProblemCgs>::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real
 {
-	quokka::valarray<double, nGroups_> kappaPVec{};
-	for (int i = 0; i < nGroups_; ++i) {
-		kappaPVec[i] = kappa;
-	}
-	return kappaPVec;
+	return kappa;
 }
 
-template <>
-AMREX_GPU_HOST_DEVICE auto RadSystem<SuOlsonProblemCgs>::ComputeFluxMeanOpacity(const double /*rho*/,
-										const double /*Tgas*/) -> quokka::valarray<double, nGroups_>
+template <> AMREX_GPU_HOST_DEVICE auto RadSystem<SuOlsonProblemCgs>::ComputeFluxMeanOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real
 {
-	return ComputePlanckOpacity(0.0, 0.0);
+	return kappa;
 }
 
 static constexpr int nmscalars_ = Physics_Traits<SuOlsonProblemCgs>::numMassScalars;

src/problems/RadMatterCoupling/test_radiation_matter_coupling.cpp

@@ -55,20 +55,14 @@ template <> struct Physics_Traits<CouplingProblem> {
 	static constexpr int nGroups = 1; // number of radiation groups
 };
 
-template <>
-AMREX_GPU_HOST_DEVICE auto RadSystem<CouplingProblem>::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> quokka::valarray<double, nGroups_>
+template <> AMREX_GPU_HOST_DEVICE auto RadSystem<CouplingProblem>::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real
 {
-	quokka::valarray<double, nGroups_> kappaPVec{};
-	for (int i = 0; i < nGroups_; ++i) {
-		kappaPVec[i] = 1.0;
-	}
-	return kappaPVec;
+	return 1.0;
 }
 
-template <>
-AMREX_GPU_HOST_DEVICE auto RadSystem<CouplingProblem>::ComputeFluxMeanOpacity(const double /*rho*/, const double /*Tgas*/) -> quokka::valarray<double, nGroups_>
+template <> AMREX_GPU_HOST_DEVICE auto RadSystem<CouplingProblem>::ComputeFluxMeanOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real
 {
-	return ComputePlanckOpacity(0.0, 0.0);
+	return 1.0;
 }
 
 static constexpr int nmscalars_ = Physics_Traits<CouplingProblem>::numMassScalars;
@@ -194,7 +188,7 @@ auto problem_main() -> int
 		std::vector<double> t_exact(nmax);
 		std::vector<double> Tgas_exact(nmax);
 		const double initial_Tgas = quokka::EOS<CouplingProblem>::ComputeTgasFromEint(rho0, Egas0);
-		const auto kappa = RadSystem<CouplingProblem>::ComputePlanckOpacity(rho0, initial_Tgas)[0];
+		const auto kappa = RadSystem<CouplingProblem>::ComputePlanckOpacity(rho0, initial_Tgas);
 
 		for (int n = 0; n < nmax; ++n) {
 			const double time_t = sim.userData_.t_vec_.at(n);

src/problems/RadMatterCouplingRSLA/test_radiation_matter_coupling_rsla.cpp

@@ -57,20 +57,14 @@ template <> struct Physics_Traits<CouplingProblem> {
 	static constexpr int nGroups = 1; // number of radiation groups
 };
 
-template <>
-AMREX_GPU_HOST_DEVICE auto RadSystem<CouplingProblem>::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> quokka::valarray<double, nGroups_>
+template <> AMREX_GPU_HOST_DEVICE auto RadSystem<CouplingProblem>::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real
 {
-	quokka::valarray<double, nGroups_> kappaPVec{};
-	for (int i = 0; i < nGroups_; ++i) {
-		kappaPVec[i] = 1.0;
-	}
-	return kappaPVec;
+	return 1.0;
 }
 
-template <>
-AMREX_GPU_HOST_DEVICE auto RadSystem<CouplingProblem>::ComputeFluxMeanOpacity(const double /*rho*/, const double /*Tgas*/) -> quokka::valarray<double, nGroups_>
+template <> AMREX_GPU_HOST_DEVICE auto RadSystem<CouplingProblem>::ComputeFluxMeanOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real
 {
-	return ComputePlanckOpacity(0.0, 0.0);
+	return 1.0;
 }
 
 static constexpr int nmscalars_ = Physics_Traits<CouplingProblem>::numMassScalars;
@@ -198,7 +192,7 @@ auto problem_main() -> int
 		std::vector<double> Tgas_rsla_exact(nmax);
 
 		const double initial_Tgas = quokka::EOS<CouplingProblem>::ComputeTgasFromEint(rho0, Egas0);
-		const auto kappa = RadSystem<CouplingProblem>::ComputePlanckOpacity(rho0, initial_Tgas)[0];
+		const auto kappa = RadSystem<CouplingProblem>::ComputePlanckOpacity(rho0, initial_Tgas);
 
 		for (int n = 0; n < nmax; ++n) {
 			const double time_t = sim.userData_.t_vec_.at(n);

src/problems/RadPulse/test_radiation_pulse.cpp

@@ -65,18 +65,13 @@ auto compute_exact_Trad(const double x, const double t) -> double
 	return 0.5 * normfac * std::exp(-(x * x) / (4.0 * width_sq));
 }
 
-template <> AMREX_GPU_HOST_DEVICE auto RadSystem<PulseProblem>::ComputePlanckOpacity(const double rho, const double Tgas) -> quokka::valarray<double, nGroups_>
+template <> AMREX_GPU_HOST_DEVICE auto RadSystem<PulseProblem>::ComputePlanckOpacity(const double rho, const double Tgas) -> amrex::Real
 {
-	quokka::valarray<double, nGroups_> kappaPVec{};
 	auto kappa = (kappa0 / rho) * std::max(std::pow(Tgas / T0, 3), 1.0);
-	for (int i = 0; i < nGroups_; ++i) {
-		kappaPVec[i] = kappa;
-	}
-	return kappaPVec;
+	return kappa;
 }
 
-template <>
-AMREX_GPU_HOST_DEVICE auto RadSystem<PulseProblem>::ComputeFluxMeanOpacity(const double rho, const double Tgas) -> quokka::valarray<double, nGroups_>
+template <> AMREX_GPU_HOST_DEVICE auto RadSystem<PulseProblem>::ComputeFluxMeanOpacity(const double rho, const double Tgas) -> amrex::Real
 {
 	return ComputePlanckOpacity(rho, Tgas);
 }

src/problems/RadShadow/test_radiation_shadow.cpp

@@ -52,18 +52,14 @@ template <> struct Physics_Traits<ShadowProblem> {
 	static constexpr int nGroups = 1;
 };
 
-template <>
-AMREX_GPU_HOST_DEVICE auto RadSystem<ShadowProblem>::ComputePlanckOpacity(const double rho, const double /*Tgas*/) -> quokka::valarray<double, nGroups_>
+template <> AMREX_GPU_HOST_DEVICE auto RadSystem<ShadowProblem>::ComputePlanckOpacity(const double rho, const double /*Tgas*/) -> amrex::Real
 {
-	quokka::valarray<double, nGroups_> kappaPVec{};
 	const amrex::Real sigma = sigma0 * std::pow(rho / rho_bg, 2);
 	const amrex::Real kappa = sigma / rho; // specific opacity [cm^2 g^-1]
-	kappaPVec.fillin(kappa);
-	return kappaPVec;
+	return kappa;
 }
 
-template <>
-AMREX_GPU_HOST_DEVICE auto RadSystem<ShadowProblem>::ComputeFluxMeanOpacity(const double rho, const double Tgas) -> quokka::valarray<double, nGroups_>
+template <> AMREX_GPU_HOST_DEVICE auto RadSystem<ShadowProblem>::ComputeFluxMeanOpacity(const double rho, const double Tgas) -> amrex::Real
 {
 	return ComputePlanckOpacity(rho, Tgas);
 }

src/problems/RadStreaming/test_radiation_streaming.cpp

@@ -48,21 +48,14 @@ template <> struct RadSystem_Traits<StreamingProblem> {
 	static constexpr int beta_order = 0;
 };
 
-template <>
-AMREX_GPU_HOST_DEVICE auto RadSystem<StreamingProblem>::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> quokka::valarray<double, nGroups_>
+template <> AMREX_GPU_HOST_DEVICE auto RadSystem<StreamingProblem>::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real
 {
-	quokka::valarray<double, nGroups_> kappaPVec{};
-	for (int g = 0; g < nGroups_; ++g) {
-		kappaPVec[g] = kappa0;
-	}
-	return kappaPVec;
+	return kappa0;
 }
 
-template <>
-AMREX_GPU_HOST_DEVICE auto RadSystem<StreamingProblem>::ComputeFluxMeanOpacity(const double /*rho*/,
-									       const double /*Tgas*/) -> quokka::valarray<double, nGroups_>
+template <> AMREX_GPU_HOST_DEVICE auto RadSystem<StreamingProblem>::ComputeFluxMeanOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real
 {
-	return ComputePlanckOpacity(0.0, 0.0);
+	return kappa0;
 }
 
 template <> void RadhydroSimulation<StreamingProblem>::setInitialConditionsOnGrid(quokka::grid grid_elem)

src/problems/RadStreamingY/test_radiation_streaming_y.cpp

@@ -51,21 +51,14 @@ template <> struct RadSystem_Traits<StreamingProblem> {
 	static constexpr int beta_order = 0;
 };
 
-template <>
-AMREX_GPU_HOST_DEVICE auto RadSystem<StreamingProblem>::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> quokka::valarray<double, nGroups_>
+template <> AMREX_GPU_HOST_DEVICE auto RadSystem<StreamingProblem>::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real
 {
-	quokka::valarray<double, nGroups_> kappaPVec{};
-	for (int g = 0; g < nGroups_; ++g) {
-		kappaPVec[g] = kappa0;
-	}
-	return kappaPVec;
+	return kappa0;
 }
 
-template <>
-AMREX_GPU_HOST_DEVICE auto RadSystem<StreamingProblem>::ComputeFluxMeanOpacity(const double /*rho*/,
-									       const double /*Tgas*/) -> quokka::valarray<double, nGroups_>
+template <> AMREX_GPU_HOST_DEVICE auto RadSystem<StreamingProblem>::ComputeFluxMeanOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real
 {
-	return ComputePlanckOpacity(0.0, 0.0);
+	return kappa0;
 }
 
 template <> void RadhydroSimulation<StreamingProblem>::setInitialConditionsOnGrid(quokka::grid grid_elem)

src/problems/RadSuOlson/test_radiation_SuOlson.cpp

@@ -70,20 +70,14 @@ template <> struct Physics_Traits<MarshakProblem> {
 	static constexpr int nGroups = 1; // number of radiation groups
 };
 
-template <>
-AMREX_GPU_HOST_DEVICE auto RadSystem<MarshakProblem>::ComputePlanckOpacity(const double rho, const double /*Tgas*/) -> quokka::valarray<double, nGroups_>
+template <> AMREX_GPU_HOST_DEVICE auto RadSystem<MarshakProblem>::ComputePlanckOpacity(const double rho, const double /*Tgas*/) -> amrex::Real
 {
-	quokka::valarray<double, nGroups_> kappaVec{};
-	for (int g = 0; g < nGroups_; ++g) {
-		kappaVec[g] = kappa / rho;
-	}
-	return kappaVec;
+	return kappa / rho;
 }
 
-template <>
-AMREX_GPU_HOST_DEVICE auto RadSystem<MarshakProblem>::ComputeFluxMeanOpacity(const double rho, const double Tgas) -> quokka::valarray<double, nGroups_>
+template <> AMREX_GPU_HOST_DEVICE auto RadSystem<MarshakProblem>::ComputeFluxMeanOpacity(const double rho, const double Tgas) -> amrex::Real
 {
-	return ComputePlanckOpacity(rho, Tgas);
+	return kappa / rho;
 }
 
 static constexpr int nmscalars_ = Physics_Traits<MarshakProblem>::numMassScalars;