Separate single- and multi-group radiation modules (#731)

### Description

The single-group radiation scheme has complex high-order terms, while
the multi-group scheme has complexy opacity functions and Jacobian
matrix. Previously, these two schemes are mixed together. It's a mess.
In this PR, I splitted `AddSourceTerms` into
`AddSourceTermsSingleGroup` and `AddSourceTermsMultiGroup`. In
`AddSourceTermsSingleGroup`, whether or not the dust model is enabled,
the Jacobian is a 2x2 matrix and can be inverted with a simple
analytical expression. Second or higher order terms are supported. In
`AddSourceTermsMultiGroup`, complexy multigroup opacity models are
included but only 0th and 1st order terms are supported. The Jacobian is
a (Ng + 1) by (Ng + 1) sparse matrix.

### Related issues

None.

### 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)](https://github.com/quokka-astro/quokka/blob/development/CONTRIBUTING.md).
- [ ] 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/QuokkaSimulation.hpp

@@ -1858,8 +1858,13 @@ 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_iteration_counter,
-					     p_num_failed_coupling, p_num_failed_dust, p_num_failed_outer_ite);
+	if constexpr (Physics_Traits<problem_t>::nGroups <= 1) {
+		RadSystem<problem_t>::AddSourceTermsSingleGroup(stateNew, radEnergySource.const_array(), indexRange, dt, stage, dustGasInteractionCoeff_,
+								p_iteration_counter, p_num_failed_coupling, p_num_failed_dust, p_num_failed_outer_ite);
+	} else {
+		RadSystem<problem_t>::AddSourceTermsMultiGroup(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/RadDust/test_rad_dust.cpp

@@ -74,23 +74,16 @@ template <> AMREX_GPU_HOST_DEVICE auto RadSystem<DustProblem>::ComputeFluxMeanOp
 	return ComputePlanckOpacity(rho, Tgas);
 }
 
-template <>
-AMREX_GPU_HOST_DEVICE auto RadSystem<DustProblem>::ComputeThermalRadiation(amrex::Real temperature, amrex::GpuArray<double, nGroups_ + 1> const &boundaries)
-    -> quokka::valarray<amrex::Real, nGroups_>
+template <> AMREX_GPU_HOST_DEVICE auto RadSystem<DustProblem>::ComputeThermalRadiationSingleGroup(amrex::Real temperature) -> amrex::Real
 {
-	auto radEnergyFractions = ComputePlanckEnergyFractions(boundaries, temperature);
-	const double power = radiation_constant_ * temperature;
-	auto Erad_g = power * radEnergyFractions;
-	return Erad_g;
+	// We assume the thermal emission is proportional to T_d in order to linearize the problem so that we can derive an analytical solution.
+	return radiation_constant_ * temperature;
 }
 
-template <>
-AMREX_GPU_HOST_DEVICE auto RadSystem<DustProblem>::ComputeThermalRadiationTempDerivative(
-    amrex::Real temperature, amrex::GpuArray<double, nGroups_ + 1> const &boundaries) -> quokka::valarray<amrex::Real, nGroups_>
+template <> AMREX_GPU_HOST_DEVICE auto RadSystem<DustProblem>::ComputeThermalRadiationTempDerivativeSingleGroup(amrex::Real /*temperature*/) -> amrex::Real
 {
-	auto radEnergyFractions = ComputePlanckEnergyFractions(boundaries, temperature);
-	const double d_power_dt = radiation_constant_;
-	return d_power_dt * radEnergyFractions;
+	// Same. We assume B(T_d) = a_rad * T_d.
+	return radiation_constant_;
 }
 
 template <> void QuokkaSimulation<DustProblem>::setInitialConditionsOnGrid(quokka::grid const &grid_elem)

src/problems/RadDustMG/test_rad_dust_MG.cpp

@@ -82,8 +82,9 @@ RadSystem<DustProblem>::DefineOpacityExponentsAndLowerValues(amrex::GpuArray<dou
 }
 
 template <>
-AMREX_GPU_HOST_DEVICE auto RadSystem<DustProblem>::ComputeThermalRadiation(amrex::Real temperature, amrex::GpuArray<double, nGroups_ + 1> const &boundaries)
-    -> quokka::valarray<amrex::Real, nGroups_>
+AMREX_GPU_HOST_DEVICE auto
+RadSystem<DustProblem>::ComputeThermalRadiationMultiGroup(amrex::Real temperature,
+							  amrex::GpuArray<double, nGroups_ + 1> const &boundaries) -> quokka::valarray<amrex::Real, nGroups_>
 {
 	auto radEnergyFractions = ComputePlanckEnergyFractions(boundaries, temperature);
 	const double power = radiation_constant_ * temperature;
@@ -92,7 +93,7 @@ AMREX_GPU_HOST_DEVICE auto RadSystem<DustProblem>::ComputeThermalRadiation(amrex
 }
 
 template <>
-AMREX_GPU_HOST_DEVICE auto RadSystem<DustProblem>::ComputeThermalRadiationTempDerivative(
+AMREX_GPU_HOST_DEVICE auto RadSystem<DustProblem>::ComputeThermalRadiationTempDerivativeMultiGroup(
     amrex::Real temperature, amrex::GpuArray<double, nGroups_ + 1> const &boundaries) -> quokka::valarray<amrex::Real, nGroups_>
 {
 	auto radEnergyFractions = ComputePlanckEnergyFractions(boundaries, temperature);

src/problems/RadMarshakVaytet/test_radiation_marshak_Vaytet.cpp

@@ -200,7 +200,7 @@ AMRSimulation<SuOlsonProblemCgs>::setCustomBoundaryConditions(const amrex::IntVe
 			T_H = T_R;
 		}
 
-		auto Erad_g = RadSystem<SuOlsonProblemCgs>::ComputeThermalRadiation(T_H, radBoundaries_g);
+		auto Erad_g = RadSystem<SuOlsonProblemCgs>::ComputeThermalRadiationMultiGroup(T_H, radBoundaries_g);
 		const double Egas = quokka::EOS<SuOlsonProblemCgs>::ComputeEintFromTgas(rho0, T_initial);
 
 		for (int g = 0; g < Physics_Traits<SuOlsonProblemCgs>::nGroups; ++g) {
@@ -231,7 +231,7 @@ template <> void QuokkaSimulation<SuOlsonProblemCgs>::setInitialConditionsOnGrid
 	amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
 		const double Egas = quokka::EOS<SuOlsonProblemCgs>::ComputeEintFromTgas(rho0, T_initial);
 		// const double Erad = a_rad * std::pow(T_initial, 4);
-		auto Erad_g = RadSystem<SuOlsonProblemCgs>::ComputeThermalRadiation(T_initial, radBoundaries_g);
+		auto Erad_g = RadSystem<SuOlsonProblemCgs>::ComputeThermalRadiationMultiGroup(T_initial, radBoundaries_g);
 
 		for (int g = 0; g < Physics_Traits<SuOlsonProblemCgs>::nGroups; ++g) {
 			state_cc(i, j, k, RadSystem<SuOlsonProblemCgs>::radEnergy_index + Physics_NumVars::numRadVars * g) = Erad_g[g];

src/problems/RadhydroPulseMGconst/test_radhydro_pulse_MG_const_kappa.cpp

@@ -220,7 +220,7 @@ template <> void QuokkaSimulation<MGproblem>::setInitialConditionsOnGrid(quokka:
 		const double rho = compute_exact_rho(x - x0);
 		const double Egas = quokka::EOS<MGproblem>::ComputeEintFromTgas(rho, Trad);
 
-		auto Erad_g = RadSystem<MGproblem>::ComputeThermalRadiation(Trad, radBoundaries_g);
+		auto Erad_g = RadSystem<MGproblem>::ComputeThermalRadiationMultiGroup(Trad, radBoundaries_g);
 
 		for (int g = 0; g < Physics_Traits<MGproblem>::nGroups; ++g) {
 			state_cc(i, j, k, RadSystem<MGproblem>::radEnergy_index + Physics_NumVars::numRadVars * g) = Erad_g[g];

src/problems/RadhydroPulseMGint/test_radhydro_pulse_MG_int.cpp

@@ -238,7 +238,7 @@ template <> void QuokkaSimulation<MGProblem>::setInitialConditionsOnGrid(quokka:
 		const double Egas = quokka::EOS<MGProblem>::ComputeEintFromTgas(rho, Trad);
 		const double v0 = v0_adv;
 
-		auto Erad_g = RadSystem<MGProblem>::ComputeThermalRadiation(Trad, radBoundaries_g);
+		auto Erad_g = RadSystem<MGProblem>::ComputeThermalRadiationMultiGroup(Trad, radBoundaries_g);
 		auto Frad_g = RadSystem<MGProblem>::ComputeFluxInDiffusionLimit(radBoundaries_g, Trad, v0);
 
 		for (int g = 0; g < Physics_Traits<MGProblem>::nGroups; ++g) {

src/problems/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp

@@ -126,7 +126,7 @@ AMRSimulation<ShockProblem>::setCustomBoundaryConditions(const amrex::IntVect &i
 		const double px_L = rho0 * v0;
 		const double Egas_L = Egas0;
 
-		auto Erad_g = RadSystem<ShockProblem>::ComputeThermalRadiation(T0, radBoundaries_g);
+		auto Erad_g = RadSystem<ShockProblem>::ComputeThermalRadiationMultiGroup(T0, radBoundaries_g);
 
 		// x1 left side boundary -- constant
 		consVar(i, j, k, RadSystem<ShockProblem>::gasDensity_index) = rho0;
@@ -146,7 +146,7 @@ AMRSimulation<ShockProblem>::setCustomBoundaryConditions(const amrex::IntVect &i
 		const double px_R = rho1 * v1;
 		const double Egas_R = Egas1;
 
-		auto Erad_g = RadSystem<ShockProblem>::ComputeThermalRadiation(T1, radBoundaries_g);
+		auto Erad_g = RadSystem<ShockProblem>::ComputeThermalRadiationMultiGroup(T1, radBoundaries_g);
 
 		// x1 right-side boundary -- constant
 		consVar(i, j, k, RadSystem<ShockProblem>::gasDensity_index) = rho1;
@@ -199,7 +199,7 @@ template <> void QuokkaSimulation<ShockProblem>::setInitialConditionsOnGrid(quok
 			temp = T1;
 		}
 
-		auto Erad_g = RadSystem<ShockProblem>::ComputeThermalRadiation(temp, radBoundaries_g);
+		auto Erad_g = RadSystem<ShockProblem>::ComputeThermalRadiationMultiGroup(temp, radBoundaries_g);
 
 		state_cc(i, j, k, RadSystem<ShockProblem>::gasDensity_index) = density;
 		state_cc(i, j, k, RadSystem<ShockProblem>::gasInternalEnergy_index) = 0.;