Use IMEX PD-ARS scheme for radiation update

src/RadMarshakAsymptotic/CMakeLists.txt

@@ -4,4 +4,4 @@ if(AMReX_GPU_BACKEND MATCHES "CUDA")
     setup_target_for_cuda_compilation(test_radiation_marshak_asymptotic)
 endif(AMReX_GPU_BACKEND MATCHES "CUDA")
 
-#add_test(NAME MarshakWaveAsymptoticDiffusion COMMAND test_radiation_marshak_asymptotic MarshakAsymptotic.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests)
+add_test(NAME MarshakWaveAsymptoticDiffusion COMMAND test_radiation_marshak_asymptotic MarshakAsymptotic.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests)

src/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp

@@ -19,12 +19,15 @@ struct SuOlsonProblemCgs {
 constexpr double kappa = 300.0;		      // cm^-1 (opacity)
 constexpr double rho0 = 2.0879373766122384;   // g cm^-3 (matter density)
 constexpr double T_hohlraum = 1.1604448449e7; // K (1 keV)
-constexpr double T_initial = 300.;	      // K
+// constexpr double T_initial = 300.;	      // K
+constexpr double T_initial = T_hohlraum * 0.001;
 
 // constexpr double kelvin_to_eV = 8.617385e-5;
 constexpr double a_rad = radiation_constant_cgs_;
 constexpr double c_v = (C::k_B / C::m_u) / (5. / 3. - 1.);
 
+constexpr double Erad_floor_ = a_rad * T_initial * T_initial * T_initial * T_initial;
+
 template <> struct quokka::EOS_Traits<SuOlsonProblemCgs> {
 	static constexpr double mean_molecular_weight = C::m_u;
 	static constexpr double boltzmann_constant = C::k_B;
@@ -35,8 +38,8 @@ template <> struct RadSystem_Traits<SuOlsonProblemCgs> {
 	static constexpr double c_light = c_light_cgs_;
 	static constexpr double c_hat = c_light_cgs_;
 	static constexpr double radiation_constant = radiation_constant_cgs_;
-	static constexpr double Erad_floor = 0.;
-	static constexpr bool compute_v_over_c_terms = true;
+	static constexpr double Erad_floor = Erad_floor_;
+	static constexpr bool compute_v_over_c_terms = false;
 };
 
 template <> struct Physics_Traits<SuOlsonProblemCgs> {
@@ -130,7 +133,10 @@ AMRSimulation<SuOlsonProblemCgs>::setCustomBoundaryConditions(const amrex::IntVe
 		//		      (1. / 12.) * (c * E_1 + 2.0 * F_1);
 
 		// use value at interface to solve for F_rad in the ghost zones
+		// const double F_bdry = 0.5 * c * E_inc - 0.5 * (c * E_0 + 2.0 * F_0);
 		const double F_bdry = 0.5 * c * E_inc - 0.5 * (c * E_0 + 2.0 * F_0);
+		// F_bdry = std::max(F_bdry, 0.0);
+		// AMREX_ASSERT(F_bdry >= 0.0);
 
 		AMREX_ASSERT(std::abs(F_bdry / (c * E_inc)) < 1.0);
 
@@ -328,14 +334,16 @@ auto problem_main() -> int
 #ifdef HAVE_PYTHON
 	// plot results
 	matplotlibcpp::clf();
-	std::map<std::string, std::string> Tgas_args;
+	std::unordered_map<std::string, std::string> Tgas_args;
 	std::map<std::string, std::string> Tgas_exact_args;
-	Tgas_args["label"] = "gas temperature";
-	Tgas_args["marker"] = ".";
+	Tgas_args["label"] = "gas temperature (numerical)";
+	Tgas_args["marker"] = "o";
+	Tgas_args["color"] = "C1";
 	Tgas_exact_args["label"] = "gas temperature (exact)";
-	Tgas_exact_args["marker"] = "x";
-	matplotlibcpp::plot(xs, Tgas_keV, Tgas_args);
+	Tgas_exact_args["color"] = "C0";
+	// Tgas_exact_args["marker"] = "x";
 	matplotlibcpp::plot(xs_exact, Tmat_exact, Tgas_exact_args);
+	matplotlibcpp::scatter(xs, Tgas_keV, 20.0, Tgas_args);
 
 	matplotlibcpp::ylim(0.0, 1.0);	// keV
 	matplotlibcpp::xlim(0.0, 0.55); // cm

src/RadhydroSimulation.hpp

@@ -247,10 +247,15 @@ template <typename problem_t> class RadhydroSimulation : public AMRSimulation<pr
 	auto computeNumberOfRadiationSubsteps(int lev, amrex::Real dt_lev_hydro) -> int;
 	void advanceRadiationSubstepAtLevel(int lev, amrex::Real time, amrex::Real dt_radiation, int iter_count, int nsubsteps,
 					    amrex::YAFluxRegister *fr_as_crse, amrex::YAFluxRegister *fr_as_fine);
+	void advanceRadiationForwardEuler(int lev, amrex::Real time, amrex::Real dt_radiation, int iter_count, int nsubsteps, amrex::YAFluxRegister *fr_as_crse,
+					  amrex::YAFluxRegister *fr_as_fine);
+	void advanceRadiationMidpointRK2(int lev, amrex::Real time, amrex::Real dt_radiation, int iter_count, int nsubsteps, amrex::YAFluxRegister *fr_as_crse,
+					 amrex::YAFluxRegister *fr_as_fine);
+
 	void subcycleRadiationAtLevel(int lev, amrex::Real time, amrex::Real dt_lev_hydro, amrex::YAFluxRegister *fr_as_crse,
 				      amrex::YAFluxRegister *fr_as_fine);
 
-	void operatorSplitSourceTerms(amrex::Array4<amrex::Real> const &stateNew, const amrex::Box &indexRange, amrex::Real time, double dt,
+	void operatorSplitSourceTerms(amrex::Array4<amrex::Real> const &stateNew, const amrex::Box &indexRange, amrex::Real time, double dt, int stage,
 				      amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const &dx, amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const &prob_lo,
 				      amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const &prob_hi);
 
@@ -1503,20 +1508,42 @@ void RadhydroSimulation<problem_t>::subcycleRadiationAtLevel(int lev, amrex::Rea
 			swapRadiationState(state_old_cc_[lev], state_new_cc_[lev]);
 		}
 
-		// advance hyperbolic radiation subsystem starting from state_old_cc_ to state_new_cc_
-		advanceRadiationSubstepAtLevel(lev, time_subcycle, dt_radiation, i, nsubSteps, fr_as_crse, fr_as_fine);
+		// We use the IMEX PD-ARS scheme to evolve the radiation subsystem and radiation-matter coupling.
+
+		// Stage 1: advance hyperbolic radiation subsystem using Forward Euler method, starting from state_old_cc_ to state_new_cc_
+		advanceRadiationForwardEuler(lev, time_subcycle, dt_radiation, i, nsubSteps, fr_as_crse, fr_as_fine);
+
+		// new radiation state is stored in state_new_cc_
+		// new hydro state is stored in state_new_cc_ (always the case during radiation update)
+
+		if constexpr (IMEX_a22 > 0.0) {
+			// matter-radiation exchange source terms of stage 1
+			for (amrex::MFIter iter(state_new_cc_[lev]); iter.isValid(); ++iter) {
+				const amrex::Box &indexRange = iter.validbox();
+				auto const &stateNew = state_new_cc_[lev].array(iter);
+				auto const &prob_lo = geom[lev].ProbLoArray();
+				auto const &prob_hi = geom[lev].ProbHiArray();
+				// update state_new_cc_[lev] in place (updates both radiation and hydro vars)
+				// Note that only a fraction (IMEX_a32) of the matter-radiation exchange source terms are added to hydro. This ensures that the
+				// hydro properties get to t + IMEX_a32 dt in terms of matter-radiation exchange.
+				operatorSplitSourceTerms(stateNew, indexRange, time_subcycle, dt_radiation, 1, dx, prob_lo, prob_hi);
+			}
+		}
+
+		// Stage 2: advance hyperbolic radiation subsystem using midpoint RK2 method, starting from state_old_cc_ to state_new_cc_
+		advanceRadiationMidpointRK2(lev, time_subcycle, dt_radiation, i, nsubSteps, fr_as_crse, fr_as_fine);
 
 		// new radiation state is stored in state_new_cc_
 		// new hydro state is stored in state_new_cc_ (always the case during radiation update)
 
-		// matter-radiation exchange source terms
+		// Add the matter-radiation exchange source terms to the radiation subsystem and evolve by (1 - IMEX_a32) * dt
 		for (amrex::MFIter iter(state_new_cc_[lev]); iter.isValid(); ++iter) {
 			const amrex::Box &indexRange = iter.validbox();
 			auto const &stateNew = state_new_cc_[lev].array(iter);
 			auto const &prob_lo = geom[lev].ProbLoArray();
 			auto const &prob_hi = geom[lev].ProbHiArray();
 			// update state_new_cc_[lev] in place (updates both radiation and hydro vars)
-			operatorSplitSourceTerms(stateNew, indexRange, time_subcycle, dt_radiation, dx, prob_lo, prob_hi);
+			operatorSplitSourceTerms(stateNew, indexRange, time_subcycle, dt_radiation, 2, dx, prob_lo, prob_hi);
 		}
 
 		// new hydro+radiation state is stored in state_new_cc_
@@ -1596,9 +1623,86 @@ void RadhydroSimulation<problem_t>::advanceRadiationSubstepAtLevel(int lev, amre
 	}
 }
 
+template <typename problem_t>
+void RadhydroSimulation<problem_t>::advanceRadiationForwardEuler(int lev, amrex::Real time, amrex::Real dt_radiation, int const /*iter_count*/,
+								 int const /*nsubsteps*/, amrex::YAFluxRegister *fr_as_crse, amrex::YAFluxRegister *fr_as_fine)
+{
+	// get cell sizes
+	auto const &dx = geom[lev].CellSizeArray();
+
+	// update ghost zones [old timestep]
+	fillBoundaryConditions(state_old_cc_[lev], state_old_cc_[lev], lev, time, quokka::centering::cc, quokka::direction::na, PreInterpState,
+			       PostInterpState);
+
+	// advance all grids on local processor (Stage 1 of integrator)
+	for (amrex::MFIter iter(state_new_cc_[lev]); iter.isValid(); ++iter) {
+		const amrex::Box &indexRange = iter.validbox();
+		auto const &stateOld = state_old_cc_[lev].const_array(iter);
+		auto const &stateNew = state_new_cc_[lev].array(iter);
+		auto [fluxArrays, fluxDiffusiveArrays] = computeRadiationFluxes(stateOld, indexRange, ncompHyperbolic_, dx);
+
+		// Stage 1 of RK2-SSP
+		RadSystem<problem_t>::PredictStep(
+		    stateOld, stateNew, {AMREX_D_DECL(fluxArrays[0].array(), fluxArrays[1].array(), fluxArrays[2].array())},
+		    {AMREX_D_DECL(fluxDiffusiveArrays[0].const_array(), fluxDiffusiveArrays[1].const_array(), fluxDiffusiveArrays[2].const_array())},
+		    dt_radiation, dx, indexRange, ncompHyperbolic_);
+
+		if (do_reflux) {
+			// increment flux registers
+			// WARNING: as written, diffusive flux correction is not compatible with reflux!!
+			auto expandedFluxes = expandFluxArrays(fluxArrays, nstartHyperbolic_, state_new_cc_[lev].nComp());
+			incrementFluxRegisters(iter, fr_as_crse, fr_as_fine, expandedFluxes, lev, 0.5 * dt_radiation);
+		}
+	}
+
+	// update ghost zones [intermediate stage stored in state_new_cc_]
+	fillBoundaryConditions(state_new_cc_[lev], state_new_cc_[lev], lev, (time + dt_radiation), quokka::centering::cc, quokka::direction::na, PreInterpState,
+			       PostInterpState);
+}
+
+template <typename problem_t>
+void RadhydroSimulation<problem_t>::advanceRadiationMidpointRK2(int lev, amrex::Real time, amrex::Real dt_radiation, int const /*iter_count*/,
+								int const /*nsubsteps*/, amrex::YAFluxRegister *fr_as_crse, amrex::YAFluxRegister *fr_as_fine)
+{
+	auto const &dx = geom[lev].CellSizeArray();
+
+	// update ghost zones [intermediate stage stored in state_new_cc_]
+	fillBoundaryConditions(state_new_cc_[lev], state_new_cc_[lev], lev, (time + dt_radiation), quokka::centering::cc, quokka::direction::na, PreInterpState,
+			       PostInterpState);
+
+	// advance all grids on local processor (Stage 2 of integrator)
+	for (amrex::MFIter iter(state_new_cc_[lev]); iter.isValid(); ++iter) {
+		const amrex::Box &indexRange = iter.validbox();
+		auto const &stateOld = state_old_cc_[lev].const_array(iter);
+		auto const &stateInter = state_new_cc_[lev].const_array(iter);
+		auto const &stateNew = state_new_cc_[lev].array(iter);
+		auto [fluxArraysOld, fluxDiffusiveArraysOld] = computeRadiationFluxes(stateOld, indexRange, ncompHyperbolic_, dx);
+		auto [fluxArrays, fluxDiffusiveArrays] = computeRadiationFluxes(stateInter, indexRange, ncompHyperbolic_, dx);
+
+		// Stage 2 of RK2-SSP
+		RadSystem<problem_t>::AddFluxesRK2(
+		    stateNew, stateOld, stateInter, {AMREX_D_DECL(fluxArraysOld[0].array(), fluxArraysOld[1].array(), fluxArraysOld[2].array())},
+		    {AMREX_D_DECL(fluxArrays[0].array(), fluxArrays[1].array(), fluxArrays[2].array())},
+		    {AMREX_D_DECL(fluxDiffusiveArraysOld[0].const_array(), fluxDiffusiveArraysOld[1].const_array(), fluxDiffusiveArraysOld[2].const_array())},
+		    {AMREX_D_DECL(fluxDiffusiveArrays[0].const_array(), fluxDiffusiveArrays[1].const_array(), fluxDiffusiveArrays[2].const_array())},
+		    dt_radiation, dx, indexRange, ncompHyperbolic_);
+
+		if (do_reflux) {
+			// increment flux registers
+			// WARNING: as written, diffusive flux correction is not compatible with reflux!!
+			auto expandedFluxes = expandFluxArrays(fluxArrays, nstartHyperbolic_, state_new_cc_[lev].nComp());
+			incrementFluxRegisters(iter, fr_as_crse, fr_as_fine, expandedFluxes, lev, 0.5 * dt_radiation);
+		}
+	}
+
+	// update ghost zones [intermediate stage stored in state_new_cc_]
+	fillBoundaryConditions(state_new_cc_[lev], state_new_cc_[lev], lev, (time + dt_radiation), quokka::centering::cc, quokka::direction::na, PreInterpState,
+			       PostInterpState);
+}
+
 template <typename problem_t>
 void RadhydroSimulation<problem_t>::operatorSplitSourceTerms(amrex::Array4<amrex::Real> const &stateNew, const amrex::Box &indexRange, const amrex::Real time,
-							     const double dt, amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const &dx,
+							     const double dt, const int stage, amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const &dx,
 							     amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const &prob_lo,
 							     amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const &prob_hi)
 {
@@ -1614,7 +1718,7 @@ void RadhydroSimulation<problem_t>::operatorSplitSourceTerms(amrex::Array4<amrex
 	RadSystem<problem_t>::SetRadEnergySource(radEnergySource.array(), indexRange, dx, prob_lo, prob_hi, time + dt);
 
 	// cell-centered source terms
-	RadSystem<problem_t>::AddSourceTerms(stateNew, radEnergySource.const_array(), advectionFluxes.const_array(), indexRange, dt);
+	RadSystem<problem_t>::AddSourceTerms(stateNew, radEnergySource.const_array(), advectionFluxes.const_array(), indexRange, dt, stage);
 }
 
 template <typename problem_t>