Fix radiation transport error in Y and Z directions. (#633)

### Description
Fix bug #632 . This is a dumb mistake I made on using `AMREX_D_TERM` in
`AddFluxesRK2`, causing a wrong calculation of the HLL flux in 2D/3D
cases. The mistake was introduced in the "Use IMEX PD-ARS ..." commit
(78c5970) 6 months ago. Because of this
mistake, all 2D/3D simulations since then with non-zero y/z-component
radiation fluxes are wrong.

To avoid mistakes like this in the future, I added a new test of
radiation streaming in the Y direction. I plan to add a better test
(say, the RadShock test) in both Y and Z direction, with the following
grid configuration: [4, 128, 4] and [4, 4, 128]. Or, alternatively, we
should have a non-trivial 3D test with an exact solution.

### Related issues
Fixed #632 

### Checklist
_Before this pull request can be reviewed, all of these tasks should be
completed. Denote completed tasks with an `x` inside the square brackets
`[ ]` in the Markdown source below:_
- [x] I have added a description (see above).
- [x] I have added a link to any related issues see (see above).
- [x] I have read the [Contributing
Guide](https://github.com/quokka-astro/quokka/blob/development/CONTRIBUTING.md).
- [x] I have added tests for any new physics that this PR adds to the
code.
- [x] I have tested this PR on my local computer and all tests pass.
- [x] I have manually triggered the GPU tests with the magic comment
`/azp run`.
- [x] I have requested a reviewer for this PR.

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

src/CMakeLists.txt

@@ -164,6 +164,7 @@ add_subdirectory(RadMatterCouplingRSLA)
 add_subdirectory(RadPulse)
 add_subdirectory(RadShadow)
 add_subdirectory(RadStreaming)
+add_subdirectory(RadStreamingY)
 add_subdirectory(RadSuOlson)
 add_subdirectory(RadTophat)
 add_subdirectory(RadTube)

src/RadStreamingY/CMakeLists.txt

@@ -0,0 +1,9 @@
+if (AMReX_SPACEDIM GREATER_EQUAL 2)
+	add_executable(test_radiation_streaming_y test_radiation_streaming_y.cpp ../fextract.cpp ${QuokkaObjSources})
+
+	if(AMReX_GPU_BACKEND MATCHES "CUDA")
+			setup_target_for_cuda_compilation(test_radiation_streaming_y)
+	endif(AMReX_GPU_BACKEND MATCHES "CUDA")
+
+	add_test(NAME RadiationStreamingY COMMAND test_radiation_streaming_y RadStreamingY.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests)
+endif()

src/RadStreamingY/test_radiation_streaming_y.cpp

@@ -0,0 +1,270 @@
+//==============================================================================
+// 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_streaming.cpp
+/// \brief Defines a test problem for radiation in the free-streaming regime.
+///
+
+#include "test_radiation_streaming_y.hpp"
+#include "AMReX.H"
+#include "AMReX_BLassert.H"
+#include "RadhydroSimulation.hpp"
+#include "fextract.hpp"
+#include "valarray.hpp"
+
+struct StreamingProblem {
+};
+
+constexpr int direction = 1;
+
+constexpr double initial_Erad = 1.0e-5;
+constexpr double initial_Egas = 1.0e-5;
+constexpr double c = 1.0;	   // speed of light
+constexpr double chat = c;	   // reduced speed of light
+constexpr double kappa0 = 1.0e-10; // opacity
+constexpr double rho = 1.0;
+
+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 = 1.0;
+	static constexpr double Erad_floor = initial_Erad;
+	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_>
+{
+	quokka::valarray<double, nGroups_> kappaPVec{};
+	for (int g = 0; g < nGroups_; ++g) {
+		kappaPVec[g] = kappa0;
+	}
+	return kappaPVec;
+}
+
+template <>
+AMREX_GPU_HOST_DEVICE auto RadSystem<StreamingProblem>::ComputeFluxMeanOpacity(const double /*rho*/, const double /*Tgas*/)
+    -> quokka::valarray<double, nGroups_>
+{
+	return ComputePlanckOpacity(0.0, 0.0);
+}
+
+template <> void RadhydroSimulation<StreamingProblem>::setInitialConditionsOnGrid(quokka::grid 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;
+
+	// loop over the grid and set the initial condition
+	amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
+		state_cc(i, j, k, RadSystem<StreamingProblem>::radEnergy_index) = Erad0;
+		state_cc(i, j, k, RadSystem<StreamingProblem>::x1RadFlux_index) = 0;
+		state_cc(i, j, k, RadSystem<StreamingProblem>::x2RadFlux_index) = 0;
+		state_cc(i, j, k, RadSystem<StreamingProblem>::x3RadFlux_index) = 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();
+
+	if constexpr (direction == 0) {
+		if (i < lo[0]) {
+			// streaming inflow boundary
+			const double Erad = 1.0;
+			const double Frad = c * Erad;
+
+			// x1 left side boundary (Marshak)
+			consVar(i, j, k, RadSystem<StreamingProblem>::radEnergy_index) = Erad;
+			consVar(i, j, k, RadSystem<StreamingProblem>::x1RadFlux_index) = Frad;
+			consVar(i, j, k, RadSystem<StreamingProblem>::x2RadFlux_index) = 0.;
+			consVar(i, j, k, RadSystem<StreamingProblem>::x3RadFlux_index) = 0.;
+		} else if (i >= hi[0]) {
+			// right-side boundary -- constant
+			const double Erad = initial_Erad;
+			consVar(i, j, k, RadSystem<StreamingProblem>::radEnergy_index) = Erad;
+			consVar(i, j, k, RadSystem<StreamingProblem>::x1RadFlux_index) = 0;
+			consVar(i, j, k, RadSystem<StreamingProblem>::x2RadFlux_index) = 0;
+			consVar(i, j, k, RadSystem<StreamingProblem>::x3RadFlux_index) = 0;
+		}
+	} else {
+		if (j < lo[1]) {
+			// streaming inflow boundary
+			const double Erad = 1.0;
+			const double Frad = c * Erad;
+
+			// x1 left side boundary (Marshak)
+			consVar(i, j, k, RadSystem<StreamingProblem>::radEnergy_index) = Erad;
+			consVar(i, j, k, RadSystem<StreamingProblem>::x1RadFlux_index) = 0.;
+			consVar(i, j, k, RadSystem<StreamingProblem>::x2RadFlux_index) = Frad;
+			consVar(i, j, k, RadSystem<StreamingProblem>::x3RadFlux_index) = 0.;
+		} else if (j >= hi[1]) {
+			// right-side boundary -- constant
+			const double Erad = initial_Erad;
+			consVar(i, j, k, RadSystem<StreamingProblem>::radEnergy_index) = Erad;
+			consVar(i, j, k, RadSystem<StreamingProblem>::x1RadFlux_index) = 0;
+			consVar(i, j, k, RadSystem<StreamingProblem>::x2RadFlux_index) = 0;
+			consVar(i, j, k, RadSystem<StreamingProblem>::x3RadFlux_index) = 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;
+	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) {
+		// assert at compile time
+		if constexpr (direction == 0) {
+			BCs_cc[n].setLo(0, amrex::BCType::ext_dir);  // Dirichlet x1
+			BCs_cc[n].setHi(0, amrex::BCType::foextrap); // extrapolate x1
+			BCs_cc[n].setLo(1, amrex::BCType::int_dir);  // periodic
+			BCs_cc[n].setHi(1, amrex::BCType::int_dir);
+		} else {
+			BCs_cc[n].setLo(0, amrex::BCType::int_dir); // periodic
+			BCs_cc[n].setHi(0, amrex::BCType::int_dir);
+			BCs_cc[n].setLo(1, amrex::BCType::ext_dir);  // Dirichlet x1
+			BCs_cc[n].setHi(1, amrex::BCType::foextrap); // extrapolate x1
+		}
+	}
+
+	// Problem initialization
+	RadhydroSimulation<StreamingProblem> sim(BCs_cc);
+
+	// read tmax from inputs file
+	amrex::ParmParse pp; // NOLINT
+	pp.query("max_time", tmax);
+
+	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), direction, 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> xs(nx);
+	for (int i = 0; i < nx; ++i) {
+		amrex::Real const x = position[i];
+		xs.at(i) = x;
+		erad_exact.at(i) = (x <= chat * tmax) ? 1.0 : 0.0;
+		erad.at(i) = values.at(RadSystem<StreamingProblem>::radEnergy_index)[i];
+	}
+
+	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.05;
+	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.1, 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]));
+	if constexpr (direction == 0) {
+		matplotlibcpp::save("./radiation_streaming_x.pdf");
+	} else {
+		matplotlibcpp::save("./radiation_streaming_y.pdf");
+	}
+#endif // HAVE_PYTHON
+
+	// Cleanup and exit
+	amrex::Print() << "Finished." << std::endl;
+	return status;
+}

src/RadStreamingY/test_radiation_streaming_y.hpp

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

src/radiation_system.hpp

@@ -508,8 +508,8 @@ void RadSystem<problem_t>::AddFluxesRK2(array_t &U_new, arrayconst_t &U0, arrayc
 			const double FzU_1 = (dt / dz) * (x3Flux(i, j, k, n) - x3Flux(i, j, k + 1, n));
 #endif
 			// save results in cons_new
-			cons_new[n] = (1.0 - IMEX_a32) * U_0 + IMEX_a32 * U_1 + ((0.5 - IMEX_a32) * AMREX_D_TERM(FxU_0, +FyU_0, +FzU_0)) +
-				      (0.5 * AMREX_D_TERM(FxU_1, +FyU_1, +FzU_1));
+			cons_new[n] = (1.0 - IMEX_a32) * U_0 + IMEX_a32 * U_1 + ((0.5 - IMEX_a32) * (AMREX_D_TERM(FxU_0, +FyU_0, +FzU_0))) +
+				      (0.5 * (AMREX_D_TERM(FxU_1, +FyU_1, +FzU_1)));
 		}
 
 		if (!isStateValid(cons_new)) {

tests/RadStreamingY.in

@@ -0,0 +1,24 @@
+max_time = 0.2
+
+# *****************************************************************
+# Problem size and geometry
+# *****************************************************************
+geometry.prob_lo     =  0.0  0.0  0.0 
+geometry.prob_hi     =  1.0  1.0  1.0
+geometry.is_periodic =  1    0    1
+
+# *****************************************************************
+# VERBOSITY
+# *****************************************************************
+amr.v              = 0       # verbosity in Amr
+
+# *****************************************************************
+# Resolution and refinement
+# *****************************************************************
+amr.n_cell          = 4 100 4
+amr.max_level       = 0     # number of levels = max_level + 1
+amr.blocking_factor = 4     # grid size must be divisible by this
+
+do_reflux = 0
+do_subcycle = 0
+suppress_output = 0