improve convergence rate of the radiation Newton-Raphson solver

src/QuokkaSimulation.hpp

@@ -131,6 +131,7 @@ template <typename problem_t> class QuokkaSimulation : public AMRSimulation<prob
 	amrex::Real pressureFloor_ = 0.;
 
 	bool use_wavespeed_correction_ = false;
+	bool radiation_enable_counter_ = 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)
@@ -251,8 +252,8 @@ template <typename problem_t> class QuokkaSimulation : public AMRSimulation<prob
 
 	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, int *num_failed_coupling, int *num_failed_dust,
-				      int *p_num_failed_outer_ite);
+				      amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const &prob_hi, int *p_iteration_counter, int *num_failed_coupling,
+				      int *num_failed_dust, int *p_num_failed_outer_ite);
 
 	auto computeRadiationFluxes(amrex::Array4<const amrex::Real> const &consVar, const amrex::Box &indexRange, int nvars,
 				    amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> dx)
@@ -390,6 +391,7 @@ template <typename problem_t> void QuokkaSimulation<problem_t>::readParmParse()
 		rpp.query("reconstruction_order", radiationReconstructionOrder_);
 		rpp.query("cfl", radiationCflNumber_);
 		rpp.query("dust_gas_interaction_coeff", dustGasInteractionCoeff_);
+		rpp.query("enable_counter", radiation_enable_counter_);
 	}
 }
 
@@ -1619,9 +1621,11 @@ void QuokkaSimulation<problem_t>::subcycleRadiationAtLevel(int lev, amrex::Real
 		amrex::Gpu::Buffer<int> num_failed_coupling({0});
 		amrex::Gpu::Buffer<int> num_failed_dust({0});
 		amrex::Gpu::Buffer<int> num_failed_outer({0});
+		amrex::Gpu::Buffer<int> iteration_counter({0, 0, 0});
 		int *p_num_failed_coupling = num_failed_coupling.data();
 		int *p_num_failed_dust = num_failed_dust.data();
 		int *p_num_failed_outer = num_failed_outer.data();
+		int *p_iteration_counter = iteration_counter.data();
 
 		if constexpr (IMEX_a22 > 0.0) {
 			// matter-radiation exchange source terms of stage 1
@@ -1634,8 +1638,8 @@ void QuokkaSimulation<problem_t>::subcycleRadiationAtLevel(int lev, amrex::Real
 				// 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, p_num_failed_coupling,
-							 p_num_failed_dust, p_num_failed_outer);
+				operatorSplitSourceTerms(stateNew, indexRange, time_subcycle, dt_radiation, 1, dx, prob_lo, prob_hi, p_iteration_counter,
+							 p_num_failed_coupling, p_num_failed_dust, p_num_failed_outer);
 			}
 		}
 
@@ -1652,8 +1656,16 @@ void QuokkaSimulation<problem_t>::subcycleRadiationAtLevel(int lev, amrex::Real
 			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, 2, dx, prob_lo, prob_hi, p_num_failed_coupling,
-						 p_num_failed_dust, p_num_failed_outer);
+			operatorSplitSourceTerms(stateNew, indexRange, time_subcycle, dt_radiation, 2, dx, prob_lo, prob_hi, p_iteration_counter,
+						 p_num_failed_coupling, p_num_failed_dust, p_num_failed_outer);
+		}
+
+		if (radiation_enable_counter_) {
+			// Iteration counter
+			auto h_iteration_counter = iteration_counter.copyToHost();
+			const double iteration_mean = static_cast<double>(h_iteration_counter[1]) / static_cast<double>(h_iteration_counter[0]);
+			amrex::Print() << "time_subcycle = " << time_subcycle << ", (mean, max) number of Newton-Raphson iterations are " << iteration_mean
+				       << ", " << h_iteration_counter[2] << "\n";
 		}
 
 		const int nf_coupling = *(num_failed_coupling.copyToHost());
@@ -1820,8 +1832,8 @@ template <typename problem_t>
 void QuokkaSimulation<problem_t>::operatorSplitSourceTerms(amrex::Array4<amrex::Real> const &stateNew, const amrex::Box &indexRange, const amrex::Real time,
 							   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, int *p_num_failed_coupling,
-							   int *p_num_failed_dust, int *p_num_failed_outer_ite)
+							   amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const &prob_hi, int *p_iteration_counter,
+							   int *p_num_failed_coupling, int *p_num_failed_dust, int *p_num_failed_outer_ite)
 {
 	amrex::FArrayBox radEnergySource(indexRange, Physics_Traits<problem_t>::nGroups,
 					 amrex::The_Async_Arena()); // cell-centered scalar
@@ -1832,8 +1844,8 @@ void QuokkaSimulation<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(), indexRange, dt, stage, dustGasInteractionCoeff_, p_num_failed_coupling,
-					     p_num_failed_dust, p_num_failed_outer_ite);
+	RadSystem<problem_t>::AddSourceTerms(stateNew, radEnergySource.const_array(), indexRange, dt, stage, dustGasInteractionCoeff_, p_iteration_counter,
+					     p_num_failed_coupling, p_num_failed_dust, p_num_failed_outer_ite);
 }
 
 template <typename problem_t>

src/problems/CMakeLists.txt

@@ -35,6 +35,7 @@ add_subdirectory(RadTophat)
 add_subdirectory(RadTube)
 add_subdirectory(RadDust)
 add_subdirectory(RadDustMG)
+add_subdirectory(RadMarshakHot)
 
 add_subdirectory(RadhydroShell)
 add_subdirectory(RadhydroShock)

src/problems/RadMarshakHot/CMakeLists.txt

@@ -0,0 +1,7 @@
+add_executable(test_radiation_marshak_hot test_radiation_marshak_hot.cpp ../../util/fextract.cpp ${QuokkaObjSources})
+
+if(AMReX_GPU_BACKEND MATCHES "CUDA")
+    setup_target_for_cuda_compilation(test_radiation_marshak_hot)
+endif(AMReX_GPU_BACKEND MATCHES "CUDA")
+
+add_test(NAME RadiationMarshakHot COMMAND test_radiation_marshak_hot RadMarshakHot.in ${QuokkaTestParams} WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests)

src/problems/RadMarshakHot/test_radiation_marshak_hot.cpp

@@ -0,0 +1,266 @@
+//==============================================================================
+// 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 test_radiation_marshak_hot.cpp
+/// \brief Defines a test Marshak wave problem with extremely high specific heat capacity in radiation.
+///
+
+#include "test_radiation_marshak_hot.hpp"
+#include "AMReX.H"
+#include "QuokkaSimulation.hpp"
+#include "util/fextract.hpp"
+#include "util/valarray.hpp"
+
+struct StreamingProblem {
+};
+
+constexpr double initial_Erad = 1.0e-5;
+constexpr double erad_floor = 1.0e-15;
+constexpr double initial_Egas = 1.0e-5;
+constexpr double c = 1.0;	 // speed of light
+constexpr double chat = 1.0;	 // reduced speed of light
+constexpr double kappa0 = 1.0e4; // opacity
+constexpr double rho = 1.0;
+constexpr double a_rad = 1.0e15;
+constexpr double EradL = a_rad;
+
+template <> struct quokka::EOS_Traits<StreamingProblem> {
+	static constexpr double mean_molecular_weight = 1.0;
+	static constexpr double boltzmann_constant = 1.0;
+	static constexpr double gamma = 5. / 3.;
+};
+
+template <> struct Physics_Traits<StreamingProblem> {
+	// cell-centred
+	static constexpr bool is_hydro_enabled = false;
+	static constexpr int numMassScalars = 0;		     // number of mass scalars
+	static constexpr int numPassiveScalars = numMassScalars + 0; // number of passive scalars
+	static constexpr bool is_radiation_enabled = true;
+	// face-centred
+	static constexpr bool is_mhd_enabled = false;
+	static constexpr int nGroups = 1; // number of radiation groups
+};
+
+template <> struct RadSystem_Traits<StreamingProblem> {
+	static constexpr double c_light = c;
+	static constexpr double c_hat = chat;
+	static constexpr double radiation_constant = a_rad;
+	static constexpr double Erad_floor = erad_floor;
+	static constexpr int beta_order = 0;
+	static constexpr bool enable_dust_gas_thermal_coupling_model = false;
+};
+
+template <> AMREX_GPU_HOST_DEVICE auto RadSystem<StreamingProblem>::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real
+{
+	return kappa0;
+}
+
+template <> AMREX_GPU_HOST_DEVICE auto RadSystem<StreamingProblem>::ComputeFluxMeanOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real
+{
+	return kappa0;
+}
+
+template <> void QuokkaSimulation<StreamingProblem>::setInitialConditionsOnGrid(quokka::grid const &grid_elem)
+{
+	const amrex::Box &indexRange = grid_elem.indexRange_;
+	const amrex::Array4<double> &state_cc = grid_elem.array_;
+
+	const auto Erad0 = initial_Erad;
+	const auto Egas0 = initial_Egas;
+
+	// calculate radEnergyFractions
+	quokka::valarray<amrex::Real, Physics_Traits<StreamingProblem>::nGroups> radEnergyFractions{};
+	for (int g = 0; g < Physics_Traits<StreamingProblem>::nGroups; ++g) {
+		radEnergyFractions[g] = 1.0 / Physics_Traits<StreamingProblem>::nGroups;
+	}
+
+	// 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) = Erad0 * radEnergyFractions[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;
+		}
+		state_cc(i, j, k, RadSystem<StreamingProblem>::gasEnergy_index) = Egas0;
+		state_cc(i, j, k, RadSystem<StreamingProblem>::gasDensity_index) = rho;
+		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.;
+	});
+}
+
+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*/)
+{
+#if (AMREX_SPACEDIM == 1)
+	auto i = iv.toArray()[0];
+	int j = 0;
+	int k = 0;
+#endif
+#if (AMREX_SPACEDIM == 2)
+	auto [i, j] = iv.toArray();
+	int k = 0;
+#endif
+#if (AMREX_SPACEDIM == 3)
+	auto [i, j, k] = iv.toArray();
+#endif
+
+	amrex::Box const &box = geom.Domain();
+	amrex::GpuArray<int, 3> lo = box.loVect3d();
+	amrex::GpuArray<int, 3> hi = box.hiVect3d();
+
+	// calculate radEnergyFractions
+	quokka::valarray<amrex::Real, Physics_Traits<StreamingProblem>::nGroups> radEnergyFractions{};
+	for (int g = 0; g < Physics_Traits<StreamingProblem>::nGroups; ++g) {
+		radEnergyFractions[g] = 1.0 / Physics_Traits<StreamingProblem>::nGroups;
+	}
+
+	if (i < lo[0]) {
+		// streaming inflow boundary
+		const double Erad = EradL;
+		const double Frad = c * Erad;
+
+		// 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) = Erad * radEnergyFractions[g];
+			consVar(i, j, k, RadSystem<StreamingProblem>::x1RadFlux_index + Physics_NumVars::numRadVars * g) = Frad * radEnergyFractions[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;
+		}
+	} else if (i >= hi[0]) {
+		// right-side boundary -- constant
+		const double Erad = initial_Erad;
+		for (int g = 0; g < Physics_Traits<StreamingProblem>::nGroups; ++g) {
+			auto const Erad_g = Erad * radEnergyFractions[g];
+			consVar(i, j, k, RadSystem<StreamingProblem>::radEnergy_index + Physics_NumVars::numRadVars * g) = Erad_g;
+			consVar(i, j, k, RadSystem<StreamingProblem>::x1RadFlux_index + Physics_NumVars::numRadVars * g) = 0;
+			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;
+		}
+	}
+
+	// gas boundary conditions are the same everywhere
+	const double Egas = initial_Egas;
+	consVar(i, j, k, RadSystem<StreamingProblem>::gasEnergy_index) = Egas;
+	consVar(i, j, k, RadSystem<StreamingProblem>::gasDensity_index) = rho;
+	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.;
+}
+
+auto problem_main() -> int
+{
+	// Problem parameters
+	// const int nx = 1000;
+	// const double Lx = 1.0;
+	const double CFL_number = 0.8;
+	const double dt_max = 1e-2;
+	const double tmax = 1.0;
+	const int max_timesteps = 5000;
+
+	// Boundary conditions
+	constexpr int nvars = RadSystem<StreamingProblem>::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
+		BCs_cc[n].setHi(0, amrex::BCType::foextrap); // extrapolate x1
+		for (int i = 1; i < AMREX_SPACEDIM; ++i) {
+			BCs_cc[n].setLo(i, amrex::BCType::int_dir); // periodic
+			BCs_cc[n].setHi(i, amrex::BCType::int_dir);
+		}
+	}
+
+	// Problem initialization
+	QuokkaSimulation<StreamingProblem> sim(BCs_cc);
+
+	sim.radiationReconstructionOrder_ = 3; // PPM
+	sim.stopTime_ = tmax;
+	sim.radiationCflNumber_ = CFL_number;
+	sim.maxDt_ = dt_max;
+	sim.maxTimesteps_ = max_timesteps;
+	sim.plotfileInterval_ = -1;
+
+	// initialize
+	sim.setInitialConditions();
+
+	// evolve
+	sim.evolve();
+
+	// read output variables
+	auto [position, values] = fextract(sim.state_new_cc_[0], sim.Geom(0), 0, 0.0);
+	const int nx = static_cast<int>(position.size());
+
+	// compute error norm
+	std::vector<double> erad(nx);
+	std::vector<double> erad_exact(nx);
+	std::vector<double> T(nx);
+	std::vector<double> xs(nx);
+	for (int i = 0; i < nx; ++i) {
+		amrex::Real const x = position[i];
+		xs.at(i) = x;
+		double erad_sim = 0.0;
+		for (int g = 0; g < Physics_Traits<StreamingProblem>::nGroups; ++g) {
+			erad_sim += values.at(RadSystem<StreamingProblem>::radEnergy_index + Physics_NumVars::numRadVars * g)[i];
+		}
+		erad.at(i) = erad_sim;
+		const double e_gas = values.at(RadSystem<StreamingProblem>::gasInternalEnergy_index)[i];
+		T.at(i) = quokka::EOS<StreamingProblem>::ComputeTgasFromEint(rho, e_gas);
+
+		// compute exact solution
+		const double tau = kappa0 * rho * x;
+		erad_exact.at(i) = (x <= chat * tmax) ? EradL * std::exp(-tau) : 0.0;
+	}
+
+	double err_norm = 0.;
+	double sol_norm = 0.;
+	for (int i = 0; i < nx; ++i) {
+		err_norm += std::abs(erad[i] - erad_exact[i]);
+		sol_norm += std::abs(erad_exact[i]);
+	}
+
+	const double rel_err_norm = err_norm / sol_norm;
+	const double rel_err_tol = 0.02;
+	int status = 1;
+	if (rel_err_norm < rel_err_tol) {
+		status = 0;
+	}
+	amrex::Print() << "Relative L1 norm = " << rel_err_norm << std::endl;
+
+#ifdef HAVE_PYTHON
+	// Plot results
+	matplotlibcpp::clf();
+	// matplotlibcpp::ylim(0.0, 1.1);
+
+	std::map<std::string, std::string> erad_args;
+	std::map<std::string, std::string> erad_exact_args;
+	erad_args["label"] = "numerical solution";
+	erad_exact_args["label"] = "exact solution";
+	erad_exact_args["linestyle"] = "--";
+	matplotlibcpp::plot(xs, erad, erad_args);
+	matplotlibcpp::plot(xs, erad_exact, erad_exact_args);
+
+	matplotlibcpp::legend();
+	matplotlibcpp::title(fmt::format("t = {:.4f}", sim.tNew_[0]));
+	matplotlibcpp::save("./radiation_marshak_hot_Erad.pdf");
+
+	// plot temperature
+	matplotlibcpp::clf();
+	matplotlibcpp::plot(xs, T);
+	matplotlibcpp::title("Temperature");
+	matplotlibcpp::save("./radiation_marshak_hot_temperature.pdf");
+#endif // HAVE_PYTHON
+
+	// Cleanup and exit
+	amrex::Print() << "Finished." << std::endl;
+	// return status;
+	return 0;
+}

src/problems/RadMarshakHot/test_radiation_marshak_hot.hpp

@@ -0,0 +1,24 @@
+#ifndef TEST_RADIATION_MARSHAK_HOT_HPP_ // NOLINT
+#define TEST_RADIATION_MARSHAK_HOT_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 test_radiation_marshak_hot.hpp
+/// \brief Defines a test problem for radiation in the free-streaming regime.
+///
+
+// external headers
+#ifdef HAVE_PYTHON
+#include "util/matplotlibcpp.h"
+#endif
+#include <fmt/format.h>
+
+// internal headers
+
+#include "radiation/radiation_system.hpp"
+
+// function definitions
+
+#endif // TEST_RADIATION_MARSHAK_HOT_HPP_