Radiation Hydrodynamics

Warning: Beta feature

The radiation module has been verified against the test problems in (Wibking & Krumholz, 2022), (He et al., 2024), and (He et al., 2024), but is still marked beta for science-use maturity. Please cite the relevant methods paper and record the exact commit hash used.

Quokka solves the equations of radiation hydrodynamics (RHD) with a two-moment (M1) method in the mixed-frame formulation, accurate to first order in $v/c$. Radiation may be grey (a single frequency-integrated group) or multigroup, and the solver is asymptotic-preserving: it recovers the correct diffusion solution even when the photon mean free path is far smaller than a cell. This page summarises what is solved, how, and how to set it up. The methods are presented in full in three papers, cited throughout below: the Godunov radiation solver in (Wibking & Krumholz, 2022), the grey time-integration scheme in (He et al., 2024), and the multigroup formulation in (He et al., 2024). The stage-by-stage structure of the time integrator is documented separately in the Radiation Integrator page, and the full conservation-law system including MHD, dust, and gravity in Equations.

Equations solved

The RHD system

The radiation quantities are the lab-frame group-integrated moments of the specific intensity,

$E_g = \int_{\nu_{g-}}^{\nu_{g+}} \frac{1}{c} \oint I(\boldsymbol{n}, \nu) \, d\Omega \, d\nu \, , \quad F_g^i = \int_{\nu_{g-}}^{\nu_{g+}} \oint n^i I(\boldsymbol{n}, \nu) \, d\Omega \, d\nu \, , \quad P_g^{ij} = \int_{\nu_{g-}}^{\nu_{g+}} \frac{1}{c} \oint n^i n^j I(\boldsymbol{n}, \nu) \, d\Omega \, d\nu \, ,$

for $g = 1 \ldots N_g$, where $\nu_{g-}$ and $\nu_{g+}$ are the lower and upper frequency edges of group $g$. Together with the gas variables they obey the conservative system

$\frac{\partial}{\partial t} \left[\begin{array}{c} \rho \\ \rho \boldsymbol{v} \\ E_{\rm gas} \\ E_g \\ c^{-2} \boldsymbol{F}_g \end{array}\right] + \nabla \cdot \left[\begin{array}{c} \rho \boldsymbol{v} \\ \rho \boldsymbol{v} \otimes \boldsymbol{v} + p \\ (E_{\rm gas} + p) \boldsymbol{v} \\ \boldsymbol{F}_g \\ \mathsf{P}_g \end{array}\right] = \left[\begin{array}{c} 0 \\ \sum_g \boldsymbol{G}_g \\ c \sum_g G^0_g \\ - c\, G^0_g \\ - \boldsymbol{G}_g \end{array}\right] ,$

where $(c G^0_g, \boldsymbol{G}_g)$ is the radiation four-force of group $g$. Writing the equations in the lab frame makes them manifestly conservative, so total energy and momentum are conserved to machine precision.

The system is closed with the (Levermore, 1984) M1 closure, which expresses $\mathsf{P}_g$ in terms of $E_g$ and $\boldsymbol{F}_g$ through the reduced flux $f = |\boldsymbol{F}_g| / (c E_g)$. In multigroup runs the closure is applied group by group.

Matter-radiation coupling

The four-force is written in the mixed-frame form: the opacities and emissivities are evaluated in the comoving frame, where they are isotropic and simple, while all radiation moments stay in the lab frame. Assuming the gas is in local thermodynamic equilibrium and neglecting scattering, the frequency-integrated (grey) four-force to order $v/c$ is, following (Mihalas & Mihalas, 1984) and (Krumholz et al., 2007),

$\begin{aligned} - c G^0 &= 4 \pi \chi_{0P} B - c \chi_{0E} E + c^{-1} (2 \chi_{0E} - \chi_{0F}) v^i F^i \, , \\ - G^i &= - c^{-1} \chi_{0F} F^i + 4 \pi c^{-2} v^i \chi_{0P} B + c^{-1} \chi_{0F} v^j P^{ji} + c^{-1} (\chi_{0F} - \chi_{0E}) v^i E \, , \end{aligned}$

where $\chi_{0P}$, $\chi_{0E}$, and $\chi_{0F}$ are the comoving-frame Planck-, energy-, and flux-mean absorption coefficients and $B$ is the Planck function at the gas temperature. Quokka works with mass opacities $\kappa = \chi / \rho$ in $\mathrm{cm^2\,g^{-1}}$, which is what a problem generator supplies. The leading terms are the familiar emission, absorption, and radiation force; the terms in $v/c$ carry the work done by the radiation force on the gas and the frame-transformation (“frame-dragging”) effects that become order unity in the dynamic diffusion regime. Terms of order $v^2/c^2$ and higher are not documented here; the single-group solver can optionally include them (see beta_order below).

The multigroup four-force

For multigroup, the expressions above are integrated over each group, from $\nu_{g-}$ to $\nu_{g+}$. Quokka solves the group-integrated four-force derived in (He et al., 2024),

$\begin{aligned} - c G_g^0 &= \underbrace{4 \pi \chi_{0B,g} B_g}_{\text{emission}} - \underbrace{c \, \chi_{0E,g} E_g}_{\text{absorption}} + \underbrace{c^{-1} (1 + \alpha_{\chi_0, g}) \chi_{0F,g} \, v^i F_g^i}_{\text{work on the gas}} \, , \\[4pt] - G_g^i &= \underbrace{- c^{-1} \chi_{0F,g} F_g^i}_{\text{radiation force}} + \underbrace{\frac{4 \pi}{c^{2}} v^i \chi_{0B,g} B_g}_{\text{momentum of emission}} - \underbrace{\frac{4 \pi}{3 c^{2}} v^i \, \Delta_g (\nu \chi_0 B_{\nu})}_{\text{group coupling}} + \underbrace{c^{-1} (1 + \alpha_{\chi_0, g}) \chi_{0E,g} \, v^j P_g^{ji}}_{\text{frame dragging}} \, , \end{aligned}$

where $B_g$, $E_g$, $\boldsymbol{F}_g$, and $\mathsf{P}_g$ are the Planck function and the three radiation moments integrated over group $g$; $\chi_{0B,g}$, $\chi_{0E,g}$, and $\chi_{0F,g}$ are the comoving-frame absorption coefficients averaged across the group weighted by $B_\nu$, $E_\nu$, and $F_\nu$ respectively; $\alpha_{\chi_0,g}$ is the power-law index of the opacity across the group; and

$\Delta_g(Q) \equiv Q(\nu_{g+}) - Q(\nu_{g-})$

is the difference of a frequency-dependent quantity between the upper and lower edges of the group. See Multigroup opacity models below for how the three mean opacities and $\alpha_{\chi_0,g}$ are evaluated, and Multigroup opacities for how to supply them from a problem generator.

Term by term:

Emission. The rate at which the gas radiates thermally into group $g$. Only the part of the blackbody spectrum that falls inside the group contributes, so $B_g$ depends on the gas temperature and on where the group edges sit relative to the spectral peak.
Absorption. The rate at which group-$g$ radiation is absorbed by the gas. Emission and absorption balance when the group is in radiative equilibrium with the matter, which is what forces $E_g \to 4 \pi B_g / c$ at high optical depth.
Work on the gas. The $O(v/c)$ energy exchange that accompanies the radiation force: as the gas is pushed by the flux, the radiation does work on it. The factor $(1 + \alpha_{\chi_0,g})$ arises from the $\nu \, \partial \chi_0 / \partial \nu$ term in the lab-frame opacity — a moving observer sees Doppler-shifted frequencies, and if the opacity varies across the group, that shift changes how strongly the group is absorbed. It reduces to unity for an opacity that is constant across the group.
Radiation force. The momentum the gas absorbs from the group’s flux. This is the leading term of the momentum exchange and the one responsible for radiation pressure on matter.
Momentum of emission. Thermal emission from moving matter is beamed forward by the Doppler effect, so it carries net momentum even though it is isotropic in the comoving frame. This term is the corresponding recoil on the gas.
Group coupling. This term has no counterpart in the grey equations. The same Doppler shift that gives the thermal emission net momentum also spreads that momentum over frequency differently from the way it spreads the energy; this term is that difference, and it is what shifts momentum between groups. Because it telescopes, $\sum_g \Delta_g (\nu \chi_0 B_\nu)$ collapses to the value of $\nu \chi_0 B_\nu$ at the two ends of the whole frequency grid, so it contributes nothing to the total momentum exchange and the grey expressions are recovered exactly — it only redistributes momentum among groups. This does place a requirement on radBoundaries: the grid must be wide enough that $\nu \chi_0 B_\nu$ is negligible at both ends, otherwise the cancellation is incomplete and the total momentum exchange through $-G_g^i$ is in error.
Frame dragging. The $O(v/c)$ transformation of the group’s radiation pressure between the lab and comoving frames. Together with the work term it becomes order unity in the dynamic diffusion regime, which is why both must be retained.

Under the piecewise constant opacity model, $\alpha_{\chi_0,g} = 0$ and the three mean opacities collapse to a single value $\chi_{0,g}$, leaving

$\begin{aligned} - c G_g^0 &= \chi_{0,g} \left( 4 \pi B_g - c E_g + c^{-1} v^i F_g^i \right) , \\[4pt] - G_g^i &= \chi_{0,g} \left[ - c^{-1} F_g^i + \frac{4 \pi}{c^{2}} v^i \left( B_g - \frac{1}{3} \Delta_g (\nu B_{\nu}) \right) + c^{-1} v^j P_g^{ji} \right] . \end{aligned}$

Summing either form over all groups recovers the grey four-force of Matter-radiation coupling above.

Reduced speed of light

To relax the radiation timestep, the radiation subsystem may be solved with a reduced speed of light $\hat{c} < c$ (the RSLA), set through c_hat_over_c. This scales the transport term by $\hat{c}/c$ and leaves the equations exact when $\hat{c} = c$ (c_hat_over_c = 1). $\hat{c}$ must remain much larger than every hydrodynamic speed in the problem. Energy and momentum are conserved to machine precision only for $\hat{c} = c$.

Numerical method

The scheme is developed and tested in (He et al., 2024) for grey radiation and extended to multigroup in (He et al., 2024); the underlying Godunov radiation solver is that of (Wibking & Krumholz, 2022). A full timestep is operator-split into two parts.

Hydrodynamic transport, advanced explicitly with the PPM + RK2-SSP Godunov scheme described in (Wibking & Krumholz, 2022).
Radiation transport and matter-radiation coupling, advanced with an implicit-explicit (IMEX) scheme, subcycled with respect to the hydro step at radiation.cfl (see Runtime parameters).

Within the radiation step, the transport terms $\nabla \cdot \boldsymbol{F}_g$ and $\nabla \cdot \mathsf{P}_g$ are treated explicitly and the four-force terms implicitly. This split is what keeps the method cheap on GPUs: the implicit part contains no spatial derivatives, so every cell is solved independently and the radiation update needs no more communication than a pure hydro update.

IMEX PD-ARS

The two parts are combined with the asymptotic-preserving IMEX PD-ARS integrator, which advances the state $\boldsymbol{U}$ over $\Delta t$ in two stages,

$\begin{aligned} \boldsymbol{U}^{(n+1/2)} &= \boldsymbol{U}^{(n)} + \Delta t \, \mathsf{T}(\boldsymbol{U}^{(n)}) + \Delta t \, \mathsf{S}(\boldsymbol{U}^{(n+1/2)}) \, , \\ \boldsymbol{U}^{(n+1)} &= \boldsymbol{U}^{(n)} + \frac{\Delta t}{2} \left[ \mathsf{T}(\boldsymbol{U}^{(n)}) + \mathsf{T}(\boldsymbol{U}^{(n+1/2)}) \right] + \frac{\Delta t}{2} \left[ \mathsf{S}(\boldsymbol{U}^{(n+1/2)}) + \mathsf{S}(\boldsymbol{U}^{(n+1)}) \right] , \end{aligned}$

with $\mathsf{T}$ the transport terms and $\mathsf{S}$ the source terms. Transport and source terms enter symmetrically at both stages, which is what allows the near-exact cancellation between them that the diffusion limit requires. The scheme is second-order accurate and reduces to RK2-SSP in the streaming limit. (He et al., 2024) give the Butcher tableaux, the formal asymptotic analysis of the diffusion limit, and the numerical tests; the stage-by-stage mapping onto the code is in Radiation Integrator.

One practical consequence, demonstrated in (He et al., 2024): because the scheme is asymptotic-preserving by construction, the ad hoc correction to the HLL wavespeeds that earlier explicit radiation schemes needed in order to recover the diffusion limit is no longer required. Quokka uses the uncorrected HLL fluxes with PPM reconstruction by default. The correction survives only as the diagnostic flag use_wavespeed_correction_, used by RadMarshakAsymptotic to demonstrate what it does.

The implicit solve

Each implicit stage solves, cell by cell, a system of $4 + 4 N_g$ equations for the gas energy, the gas momentum, and the energy and flux of every group. Following (Howell & Greenough, 2003) and (Wibking & Krumholz, 2022), and as generalised to multigroup in (He et al., 2024), it is split into two nested iterations:

an inner Newton-Raphson iteration over the $1 + N_g$ energy variables (gas energy and the group exchange terms $R_g$), with $\boldsymbol{v}$ and $\boldsymbol{F}_g$ frozen;
an outer iteration that updates $\boldsymbol{F}_g$ and the gas momentum analytically, then returns to the inner solve if the velocity-dependent terms have changed.

The inner Jacobian is sparse — groups couple to the gas but not directly to each other — so (He et al., 2024) invert it by Gauss-Jordan elimination in $O(N_g)$ operations rather than $O(N_g^3)$. Outside the dynamic diffusion limit the outer loop almost always converges in one pass. The gas energy is recovered from the converged exchange terms rather than solved for independently, so energy is conserved only to the accuracy set by radiation.iteration_tolerance. Convergence tolerances are set by radiation.iteration_tolerance and radiation.iteration_tolerance_rel; the choice of per-group unknown and the round-off floor on the residual are discussed in Radiation Integrator.

Multigroup opacity models

The group-integrated four-force given in The multigroup four-force involves integrals of the opacity multiplied by a radiation quantity over each group, so a model for the frequency dependence of $\chi_0$ within a group is needed. Quokka offers two, both introduced in (He et al., 2024) and selected with the opacity_model trait. Which model to use gives the recommendation, and Multigroup opacities shows how to supply either from a problem generator.

Piecewise constant (PC)

The opacity is taken to be constant across each group, $\chi_0(\nu) = \chi_{0,g}$. The Planck-, energy-, and flux-mean opacities of a group are then all equal to that one value, and the $\partial \chi_0 / \partial \nu$ terms vanish. This is the simplest model and the one used by most earlier multigroup codes. It is accurate when the frequency grid is fine enough to resolve the variation of the opacity, and inaccurate when it is not. It is the case for which the four-force reduces to the compact form at the end of The multigroup four-force.

Piecewise power law (PPL)

The opacity is instead assumed to be a power law across each group,

$\chi_{0}(\nu) = \chi_{0,g-} \left( \frac{\nu}{\nu_{g-}} \right)^{\alpha_{\chi_0, g}} \, , \qquad \nu_{g-} \le \nu \le \nu_{g+} \, ,$

specified by its value at the lower group edge and a power-law index. This costs very little extra and is much more accurate than PC at low frequency resolution, which is the common case in practice. Power laws are also the natural description of many real opacities — infrared dust opacity and free-free opacity are both close to power laws in frequency.

Evaluating the group means also requires the shape of the radiation spectrum within each group. Assuming both $E_\nu$ and $B_\nu$ are power laws of index $\alpha_{Q,g}$ within the group gives the closed form

$\chi_{0Q,g} = \chi_{0,g-} \left[ \frac{ r_g^{\alpha_{Q,g} + 1} - 1}{\alpha_{Q,g} + 1} \right]^{-1} \left[ \frac{r_g^{\alpha_{\chi_0,g} + \alpha_{Q,g} + 1} - 1}{\alpha_{\chi_0,g} + \alpha_{Q,g} + 1} \right] , \qquad r_g \equiv \frac{\nu_{g+}}{\nu_{g-}} \, ,$

with the bracketed factors replaced by $\ln r_g$ when the corresponding exponent is $-1$. Two ways of choosing $\alpha_{Q,g}$ are implemented:

Fixed slope (PPL_opacity_fixed_slope_spectrum): assume $\nu Q_\nu$ is constant across each group, i.e. $\alpha_{Q,g} = -1$. This is not arbitrary — weighted by the spectrum itself, the mean power-law index of any spectrum that integrates to a finite energy is exactly $-1$. It costs almost nothing and performs well even with a handful of groups. This is the PPL variant to use.
Full spectrum (PPL_opacity_full_spectrum): fit $\alpha_{E,g}$ and $\alpha_{B,g}$ to the actual radiation and Planck spectra on the fly, refitting during the first few Newton iterations. This is significantly more expensive and, in the tests of (He et al., 2024), never more accurate than the fixed-slope method at any number of groups. It is kept for testing and is not recommended for production at any frequency resolution.

Note that $\alpha_{Q,g}$ matters only when $\alpha_{\chi_0,g} \ne 0$: the spectrum shape within a group is relevant only if the opacity varies across that group. Setting all exponents to zero reduces PPL to PC.

Which model to use

The choice is governed by frequency resolution, and the deciding comparison is between PC and PPL fixed slope — not between the two PPL variants. (He et al., 2024) measure the convergence of all three models on a Marshak wave with a continuously varying opacity $\chi_0 \propto \nu^{-2}$ at 4, 8, and 16 groups, and recommend:

Frequency resolution	Recommended model	Why
$N_g \lesssim 10$	`PPL_opacity_fixed_slope_spectrum`	Best accuracy. With few groups the opacity varies strongly across each bin, which is exactly what PC cannot represent.
$N_g \gtrsim 10$	`piecewise_constant_opacity`	Bins are narrow enough that a constant opacity is a good approximation, so PPL buys no accuracy, and PC is slightly cheaper because no slopes are evaluated.

The sharper form of the criterion is stated in terms of bin width rather than group count: use PC whenever the bin width in logarithmic frequency is smaller than about $W/10$, where $W$ is the logarithmic width of the whole frequency range, and PPL fixed slope otherwise.

Two things are worth knowing about the ends of this range. At high resolution PC is not merely adequate but marginally the most accurate of the three, because the flux-mean opacity relation of Flux-mean opacity introduces a small error that the PPL group means do not cancel. At low resolution the advantage of PPL is large: in the advecting radiation pulse test of (He et al., 2024), PPL fixed slope with 4 groups has errors two to three times smaller than PC with 4 groups — comparable to PC with 8 — and 4 groups with PPL fixed slope is already enough to keep every quantity accurate to better than 10 per cent. The full-spectrum variant is not recommended in either regime.

Flux-mean opacity

The flux-mean opacity $\chi_{0F,g}$ that appears in the four-force is not supplied by the user. It is computed internally from $\chi_{0E,g}$, $\chi_{0B,g}$, and the gas temperature so that $\boldsymbol{G}_g \to 0$ in an optically thick moving medium. As shown in (He et al., 2024), enforcing this relation is what guarantees that the multigroup scheme reaches the correct diffusion limit; in the PC case it reduces to $\chi_{0F,g} = \chi_{0,g}$, as expected.

Setting up a problem

Compile-time traits

Enable radiation in Physics_Traits and set the number of groups:

template <> struct Physics_Traits<MyProblem> : DefaultPhysicsTraits {
	static constexpr bool is_hydro_enabled = true;
	static constexpr bool is_radiation_enabled = true;
	static constexpr int nGroups = 4; // 1 (default) means grey radiation
};

Then specialise RadSystem_Traits. A specialisation does not inherit defaults, so c_hat_over_c, Erad_floor, and beta_order must always be set; omitting one is a compile error.

Trait	Type	Default	Meaning
`c_hat_over_c`	`double`	required	Reduced speed of light $\hat{c}/c$. Use `1.0` unless the radiation timestep is prohibitive.
`Erad_floor`	`double`	required	Floor on the total radiation energy density, in $\mathrm{erg\,cm^{-3}}$; each group gets `Erad_floor / nGroups`.
`energy_unit`	`double`	`C::ev2erg`	Unit in which `radBoundaries` is expressed. Use `C::hplanck` to give group boundaries in Hz, `C::ev2erg` to give them in eV.
`radBoundaries`	`GpuArray<double, nGroups+1>`	`{0., inf}` (grey only)	Group boundaries, monotonically increasing, in units of `energy_unit`. For `nGroups > 1`, give finite, positive end points (not `0` or `inf`) that span the full range where $\nu B_\nu$ matters, so that $\nu \chi_0 B_\nu$ is negligible at both.
`beta_order`	`int`	required	Highest order of $v/c$ retained in the four-force. `0` drops all velocity terms; `1` is the documented scheme. The single-group solver also accepts `2` and `3`.
`opacity_model`	`OpacityModel`	`single_group`	`single_group` for grey radiation; when `nGroups > 1`, one of the multigroup models — see Which model to use.

Grey opacities

For nGroups = 1, define the mean opacities $\chi_{0P}$, $\chi_{0F}$, and $\chi_{0E}$ of Matter-radiation coupling as functions of density and gas temperature. Each returns a mass opacity in $\mathrm{cm^2\,g^{-1}}$. Only ComputePlanckOpacity is mandatory; the other two default to it.

template <> struct RadSystem_Traits<MyProblem> {
	static constexpr double c_hat_over_c = 1.0;
	static constexpr double Erad_floor = 1.0e-20;
	static constexpr int beta_order = 1;
};

template <>
AMREX_GPU_HOST_DEVICE auto RadSystem<MyProblem>::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real
{
	return kappa0;
}

template <>
AMREX_GPU_HOST_DEVICE auto RadSystem<MyProblem>::ComputeFluxMeanOpacity(const double rho, const double Tgas) -> amrex::Real
{
	return ComputePlanckOpacity(rho, Tgas);
}

ComputeEnergyMeanOpacity may be specialised in the same way.

Multigroup opacities

For nGroups > 1 the grey hooks are not used. Instead specialise a single function, DefineOpacityExponentsAndLowerValues, which returns two arrays of length nGroups + 1: the power-law exponents $\alpha_{\chi_0,g}$ in element [0], and the opacity at the lower edge of each group, $\kappa_{0,g-}$, in element [1]. These are the two quantities that define the power law in Piecewise power law (PPL); under the PC model [1][g] is used directly as the constant opacity of group g, and the exponents are used only to extrapolate the opacity to the upper group edge in the group-coupling term, so they should be set consistently (0 for a constant opacity). The three group means $\chi_{0B,g}$, $\chi_{0E,g}$, and $\chi_{0F,g}$ that enter the four-force are derived from them internally.

A piecewise constant opacity, uniform across all groups:

template <> struct RadSystem_Traits<MyProblem> {
	static constexpr double c_hat_over_c = 1.0;
	static constexpr double Erad_floor = 1.0e-20;
	static constexpr double energy_unit = C::hplanck;      // radBoundaries given in Hz
	static constexpr amrex::GpuArray<double, n_groups_ + 1> radBoundaries = rad_boundaries_;
	static constexpr int beta_order = 1;
	static constexpr OpacityModel opacity_model = OpacityModel::piecewise_constant_opacity;
};

template <>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto
RadSystem<MyProblem>::DefineOpacityExponentsAndLowerValues(amrex::GpuArray<double, nGroups_ + 1> /*rad_boundaries*/, const double /*rho*/,
							   const double /*Tgas*/) -> amrex::GpuArray<amrex::GpuArray<double, nGroups_ + 1>, 2>
{
	amrex::GpuArray<amrex::GpuArray<double, nGroups_ + 1>, 2> exponents_and_values{};
	for (int g = 0; g < nGroups_ + 1; ++g) {
		exponents_and_values[0][g] = 0.0;    // constant opacity within each group
		exponents_and_values[1][g] = kappa0; // cm^2 g^-1
	}
	return exponents_and_values;
}

A piecewise power-law opacity $\kappa \propto \nu^{-2}$ differs only in the choice of model and the two lines that fill the arrays:

	static constexpr OpacityModel opacity_model = OpacityModel::PPL_opacity_fixed_slope_spectrum;

	for (int g = 0; g < nGroups_ + 1; ++g) {
		exponents_and_values[0][g] = -2.0;
		exponents_and_values[1][g] = kappa0 * std::pow(rad_boundaries[g] / nu_ref, -2.0);
	}

The function is called on the device with the current cell density and gas temperature, so the opacity may depend on both. If the opacity is tabulated rather than analytic, a convenient choice of exponent is the secant slope across the group,

$\alpha_{\chi_0,g} = \frac{\ln \left[ \chi_0(\nu_{g+}) / \chi_0(\nu_{g-}) \right]}{\ln (\nu_{g+} / \nu_{g-})} \, .$

Radiation sources

A problem can inject radiation directly — for example from stellar sources — by specialising RadSystem<problem_t>::AddRadSource, which is called before each implicit solve. See the Radiation Integrator page for the buffer conventions.

Runtime parameters

Parameter	Type	Default	Description
`radiation.cfl`	Float	`0.3`	CFL number for the radiation substeps, based on $\hat{c}$. Independent of the hydro CFL number.
`radiation.reconstruction_order`	Integer	`3`	Spatial reconstruction for the radiation variables: 1 (donor cell), 2 (PLM), 3 (PPM), 5 (extremum-preserving PPM).
`radiation.iteration_tolerance`	Float	`1e-11`	Relative tolerance on the Newton-Raphson residuals of the implicit solve.
`radiation.iteration_tolerance_rel`	Float	`-1.0` (disabled)	Optional tolerance on the relative change between consecutive Newton iterations.
`radiation.print_iteration_counts`	Boolean (0/1)	`0`	Print the number of Newton iterations per step. Useful when diagnosing a stiff or non-converging problem.
`radiation.dust_gas_interaction_coeff`	Float	`2.5e-34`	Coefficient of the dust-gas thermal coupling term, used when `ISM_Traits::enable_dust_gas_thermal_coupling_model` is enabled.

The number of radiation substeps per hydro step is computed automatically from radiation.cfl and $\hat{c}$; it is not set directly. This holds only for hydro+radiation runs with an adaptive timestep: with constant_dt set, or without hydro, there is a single radiation step per hydro step, and radiation.cfl and $\hat{c}$ do not limit it.

Test problems

The following test problems exercise the solver across the streaming, static diffusion, and dynamic diffusion regimes:

Radiative shock test — non-equilibrium radiating shock.
Matter-radiation temperature equilibrium test — the implicit coupling in isolation.
Advecting radiation pulse test — static and dynamic diffusion, single-group and multigroup.
Uniform advecting radiation in diffusive limit — the $v/c$ terms in the dynamic diffusion limit.
1D H II region and dust reprocessing test — multigroup radiation with dust.

References

Because the three Quokka methods papers share authors and year, the short citations above do not distinguish them on sight. They are, in the order a reader should approach them:

(Wibking & Krumholz, 2022) — the original Quokka paper: the Godunov radiation solver, PPM reconstruction and HLL fluxes for the radiation moments, and the M1 closure. Start here.
(He et al., 2024) — An asymptotically correct implicit-explicit time integration scheme for finite volume radiation-hydrodynamics. The IMEX PD-ARS scheme of Numerical method, its asymptotic analysis in the static and dynamic diffusion limits, and the removal of the wavespeed correction.
(He et al., 2024) — A novel numerical method for mixed-frame multigroup radiation-hydrodynamics with GPU acceleration implemented in the QUOKKA code. Everything multigroup: the group-integrated four-force of The multigroup four-force, the PC and PPL opacity models, and the sparse Newton solve.

The mixed-frame formulation itself follows (Mihalas & Mihalas, 1984) and (Krumholz et al., 2007), and the inner Newton-Raphson iteration follows (Howell & Greenough, 2003). If you use the radiation module, please cite the papers that apply to your work — see Citation.

Keyboard shortcuts

Quokka