© The Authors 2026

1 Introduction

Radiation plays a central role in the dynamics of many astrophysical systems, often dominating or significantly altering their evolution. Whether generated in situ or irradiating the system from the outside, radiation can often be the dominant force at work, even counteracting gravity in the most extreme cases. This is the case in accreting black holes (Hirose et al. 2009; Jiang et al. 2019), stellar winds (Kudritzki & Puls 2000) and disk winds (Proga et al. 1998), galactic disks (Thompson et al. 2005), and star-forming regions (Murray et al. 2010), to name only a few examples. Yet, even when it is not dynamically dominant, radiation remains an indispensable component of astrophysical fluid dynamics because it regulates the thermodynamics of the flow, thereby determining the relative importance of thermal, magnetic, and gravitational energy. Consequently, radiation must be evolved self-consistently with the fluid dynamics, what is known as radiation magnetohydrodynamics (RMHD).

Coupling radiative transfer to fluid dynamics, however, requires substantial approximations, as solving the full frequency- and angle-dependent radiative transfer equation is often computationally prohibitive, though it is feasible in some cases using Monte Carlo radiative transfer (Miller et al. 2019, 2020; Dexter et al. 2021). A common simplification is the gray approximation, in which opacities are averaged over frequency. With the advent of exascale computing, simulations that employ multiple frequency groups have become feasible, although such multigroup radiative transfer calculations remain relatively rare (Jiang et al. 2025; Melon Fuksman et al. 2025; Roth et al. 2025). Another widely used simplification is to take angular moments of the radiative transfer equation; this yields the moment equations of radiative transfer, which require a closure relation. The most common closure, flux-limited diffusion (FLD), assumes that the radiative flux is proportional to the gradient of the radiative energy density; while computationally efficient and widely used (Hirose et al. 2009; Flock et al. 2017), this approach loses directional information and is inadequate in many astrophysical regimes. A higher-order closure is the M1 approximation, in which the radiation pressure tensor is expressed as a nonlinear function of the radiative energy and flux, ensuring correct behavior in both the optically thick (Eddington) and optically thin (free-streaming) limits (Minerbo 1978; Levermore 1984; Ripoll et al. 2001; González et al. 2007; Rosdahl et al. 2013; Sądowski et al. 2013; Skinner & Ostriker 2013). The M1 closure has been successfully applied in a variety of astrophysical contexts, including accretion disks (Sądowski et al. 2013; McKinney et al. 2014; Melon Fuksman et al. 2021; Liska et al. 2022), galaxy evolution (Rosdahl et al. 2013; Rosdahl & Teyssier 2015), and core collapse supernovae and neutron star mergers, where the radiative fluid is replaced by a neutrino fluid (Mezzacappa et al. 2020; Foucart 2023). While there are higher-order approaches that keep more information about the directionality of the fluid, such as the variable Eddington tensor with a closure made using a Monte Carlo scheme (Foucart 2018) or a short-characteristic scheme (Jiang et al. 2012), we chose here to adopt the M1 approximation for its robustness, performance, and ease of implementation.

In this paper, we present the radiation module we implemented in the magnetohydrodynamic (MHD) public code IDE-FIX (Lesur et al. 2023). The main advantage of IDEFIX is that it uses the Kokkos library to ensure performance portability across a large variety of architectures, including accelerated architectures, in anticipation of the exascale era. We detail the method that we developed for the radiative transfer module in Sect. 2, the tests of the method in Sect. 3, and the performance of the radiative version of IDEFIX in Sect. 4, before concluding in Sect. 5.

2 Radiative transfer module of IDEFIX

In this section we present the radiative module that we implemented in the IDEFIX code solving for the radiative transfer moment equations in the M1 approximation. The module closely follows the implementation from Melon Fuksman & Mignone (2019) and Melon Fuksman et al. (2021) except for a few notable differences that we highlight below.

2.1 Equations

The full system for the RMHD equations is $${\partial }_t\rho +\nabla \cdot [\rho v]=0,$$(1) $${\partial }_t(\rho v)+\nabla \cdot [\rho v\otimes v-B\otimes B+\mathbb{I}]=-\rho \nabla \psi +G,$$(2) $${\partial }_{t}B+\nabla \times [-v\times B]=0,$$(3) $$\begin{array}{r}\hfill {\partial }_{t}E+\nabla \cdot [(E+)v-B(B\cdot v)]\\ \hfill =-\rho v\cdot \nabla \psi +c{G}^0-\nabla \cdot F_{\text{irr}},\end{array}$$(4) $${\partial }_t{E}_r+\nabla \cdot (\widehat{c}{F}_r)=-\widehat{c}{G}^0,$$(5) $${\partial }_t{F}_r+\nabla \cdot (\widehat{c}{P}_r)=-\widehat{c}G,$$(6)

where ρ is the density, v is the velocity, B is the magnetic field, 𝒫 ≡ P + B°/2 is the total pressure with P the thermal pressure, I is the identity tensor, ψ is the gravitational potential, F_irr is the irradiation flux, and E is the total gas energy density: $$E=e+\rho v^2/2+B^2/2,$$(7)

with e the internal energy density.

In the last two equations, c is the speed of light, ĉ is the reduced speed of light (see Sect. 2.7), E_r is the radiative energy, F_r the radiative energy flux and ℙ_r the radiative energy tensor defined, respectively, from the direction and frequency-dependent specific intensity, I_v(t, x, n) as $$E_r=\frac{1}{c}{\displaystyle {\int }_0^{\infty }d\nu {{\displaystyle \oint }}^d\Omega I_{\nu }(t,x,n)}$$(8) $$F_r^i=\frac{1}{c}{\displaystyle {\int }_0^{\infty }d\nu {{\displaystyle \oint }}^d\Omega I_{\nu }(t,x,n)n^i}$$(9) $$P_r^{ij}=\frac{1}{c}{\displaystyle {\int }_0^{\infty }d\nu {{\displaystyle \oint }}^d\Omega I_{\nu }(t,x,n)n^i{n}^j.}$$(10)

Note that Eq. (9) contains a factor, 1/c, that is usually not present in the definition of the radiative flux (see Mihalas & Mihalas 2013). The expression of Eqs. (5) and (6) takes into account this extra factor.

To close the set of the radiation moment equations, we used the closure from Levermore (1984) relating the radiation pressure to the radiation energy and the radiation flux: $$P_r^{ij}=D^{ij}{E}_r,$$(11)

with $$D^{ij}=\frac{1-\xi }{2}{\delta }^{ij}+\frac{3\xi -1}{2}{n}^i{n}^j,$$(12)

and $$\xi =\frac{3+4f^2}{5+2\sqrt{4-3f^2}},$$(13)

where δ^ij is the Kronecker delta symbol, n = F_r/∥F_r∥ is a unit vector pointing in the direction of the radiative flux, and f is the so-called reduced flux defined as $$f\equiv \frac{\parallel F_{\text{r}}\parallel }{E_r}\le 1.$$(14)

The condition that f ≤ 1 must be ensured at all times to keep a physical solution. We discuss this further in Sect. 2.5.

The source terms in the radiation equations are defined as $${\tilde{G}}^0={\kappa }_{\text{p}}\rho \left({\tilde{E}}_r-a_{\text{R}}{T}^4\right)$$(15) $$\tilde{G}=\chi \rho {\tilde{F}}_r,$$(16)

where the tildes denotes quantities in the fluid rest frame. κ_P is the Planck absorption opacity, χ ≡ κ_R + σ where κ_R and σ are the Rosseland absorption and scattering opacities, respectively. It is customary to Lorentz-transform the fluid rest frame source terms in the coordinate frame to find G° and G. However, since our main application is for protoplanetary disks, where we typically have β ≡ ∥v∥/c < 10⁻³, we neglected the Lorentz transform and set $G^0\equiv {\tilde{G}}^0\text{\hspace{0.17em}and\hspace{0.17em}}G\equiv \tilde{G}.$.

Note that in the nonrelativistic regime of protoplanetary disks, most relativistic corrections arising from the Lorentz transform (typically of order 𝒪(β) or 𝒪(β°)) can be safely neglected. However, a correction term of order 𝒪(τβ) may still be relevant. This term becomes significant in the dynamic diffusion limit, where the timescales for photon diffusion and advection are comparable. In this regime, photons diffusing through an optically thick medium are advected with the flow.

In this work, we neglected this 𝒪(τβ) correction, as its magnitude is generally comparable to the numerical error of our scheme for τβ ≈ 1 (Skinner & Ostriker 2013). Future studies focusing on accretion disks around compact objects will incorporate all relativistic corrections to address their potential impact in more extreme environments.

2.2 Units

The MHD version of IDEFIX is dimensionless. However, radiation naturally introduces units into the system of equation. To keep it conceptually simple for a user, we ask the user to set units for the density, velocity and length. These units are then used to explicitly re-dimension the radiation energy, the radiation flux and the temperature in the source terms. The hyperbolic part of the radiation equation remains dimensionless.

2.3 Implementation details

IDEFIX leverages abstract C++ class templates to enhance modularity, improve performance, and minimize code duplication. The framework employs a single, generic class template to represent all fluid types, with the template parameter defining the specific physics to be solved: whether for classical hydrodynamics, MHD, dust fluids (zero pressure), or radiative fluids. This design eliminates redundancy in critical components such as the Riemann solver, reconstruction scheme, and boundary conditions. Additionally, the use of class templates enables seamless extensibility to multigroup radiative transfer. This is achieved by instantiating multiple radiative fluid objects within a vector, requiring no structural changes to the existing implementation.

2.4 Integration scheme

We use a split operator scheme to integrate the system of equations from (1) to (6). Namely, we perform a (magneto)-hydrodynamic step, then we do a radiation transport step and finally add the radiation source terms. This general method is followed for each sub-step of the time integration loop of IDEFIX. For RK2, this is equivalent to the IMEX1 method that is used in Melon Fuksman & Mignone (2019). We followed Melon Fuksman et al. (2021) by using an explicit integration for the radiation transport (the hyperbolic part of the radiation equations) and an implicit integration for the radiation source terms solving radiation-matter coupling.

For a Euler time integration scheme, this is what a step evolving the radiation field looks like:

1. Set the boundary conditions

2. Convert the primitive variables to conservative variables

3. Reconstruct the primitive variables on the faces

4. Enforce f ≤ 1 condition on the faces

5. Compute wave speed at cell interface

6. Compute Riemann flux at cell interface

7. Evolve cell-centered conservative variables from flux divergence in each direction

8. Evolve cell-centered conservative variables from local source terms

9. Convert conservative to primitive variables

10. Update time step.

We note that in radiation hydrodynamics the sets of primitive variables and conservative variables for the radiation field are identical, namely (E_r, F_r). Hence, the primitive to conservative (or vice-versa) transformations are trivial.

2.5 Hyperbolic solver

Radiation transport focuses on the left-hand side of Eqs. (5) and (6), which corresponds to the hyperbolic component of the system. Analogous to the MHD equations, we solved this hyperbolic part using a high-order finite-volume Godunov method. In each computational cell, the primitive variables are first reconstructed at the cell faces using either a slope-limited linear scheme or higher-order methods so that $$E_{r,f}=ℛ(E_{r,c})$$(17) $$F_{r,f}^i=ℛ(F_{r,f}^i),$$(18)

where E_r,f, E_r,c, $F_{r,f}^i$ and $F_{r,c}^i$ are the radiative energy and the radiative flux in the direction i on the face and cells, respectively, and where ℛ() denotes the standard reconstruction method of IDEFIX such as a piecewise linear reconstruction using the van Leer slope limiter (PLM), compact third-order reconstruction (LimO3; Čada & Torrilhon 2009) or piecewise parabolic reconstruction using an extremum-preserving limiter (PPM; Colella & Woodward 1984, Colella & Sekora 2008). Using these reconstructed variables, a Riemann problem is solved at each interface to compute the inter-cell fluxes of radiation energy and radiation flux in all spatial directions.

In our current implementation, the user choice of reconstruction order and of slope limiter are applied to both the radiation field and hydrodynamic fluid. However, we find that radiation requires a special reconstruction scheme. Indeed, the physical constraint f ≤ 1 can be violated on the (reconstructed) faces even if it is met at the cell centers. In a more general way, the values of f computed on the cell centers and on the cell faces are likely to differ. We find that this discrepancy between cells and faces can lead to unwanted behaviors in the free-streaming limit (see Sects. 3.1 and 3.2). To avoid this issue, Melon Fuksman & Mignone (2019) switch to a flat reconstruction scheme so that E_r,f = E_r,c and $F_{r,f}^i=F_{r,c}^i$ whenever $\parallel ℛ(F_{r,c}^i)\parallel >ℛ(E_{r,c})$. In the remainder of the paper, we call this reconstruction the “PLUTO reconstruction scheme.”

However, we find that the PLUTO reconstruction scheme tends to render the scheme more diffusive and lead to numerical artifacts (see Sects. 3.1 and 3.2). In order to partially keep the benefit of high-order reconstruction scheme, we propose here an “ f -preserving reconstruction scheme.” The idea is that, after the reconstruction of E_r and $F_r^i$ on the faces, we check if $\parallel ℛ(F_{r,c}^i)\parallel >ℛ(E_{r,c})$ and if $\parallel ℛ(F_{r,c}^i)\parallel /ℛ(E_{r,c})<\parallel F_{r,c}\parallel /E_{r,c}$. If one of these conditions is met, we recompute the fluxes on the faces in each direction with the constraint that the value of f must be equal on the faces and at the center of the cell so that $$E_{r,f}=ℛ(E_{r,c}),$$(19) $$F_{r,f}^i=ℛ(F_{r,c}^i)\text{if\hspace{0.17em}}\frac{\parallel ℛ(F_{r,c}^i)\parallel }{ℛ(E_{r,c})}\le 1\text{\hspace{0.17em}or\hspace{0.17em}}\frac{\parallel ℛ(F_{r,c}^i)\parallel }{ℛ(E_{r,c})}>\frac{\parallel F_{r,c}\parallel }{E_{r,c}},$$(20) $$F_{r,f}^i=F_{r,c}^i\frac{ℛ(E_{r,c})}{E_{r,c}}\text{otherwise.}$$(21)

The second condition in Eq. (21) is introduced to prevent what we refer to as “beam widening” in the free-streaming regime. Without this constraint, the reconstruction scheme may artificially introduce a transverse pressure-like component into an otherwise free-streaming beam, thereby causing an unphysical broadening of the beam (see Sect. 3.2 for further details).

We implemented three Riemann solvers for radiation, Lax-Friedrichs-Rusanov, HLL (González et al. 2007), and HLLC Melon Fuksman & Mignone (2019). In practice, for our applications we only use HLL but all three solvers have been tested on all test cases. We do not describe here the details of the solvers and refer the reader to Melon Fuksman & Mignone (2019) for details.

2.6 Signal speed

We computed the wave velocities as in Skinner & Ostriker (2013). The M1 system of equations has four eigenvalues, two of which are degenerate; they are $${\lambda }_{1,3}=\frac{f\cos (\theta )\pm \zeta (f,\theta )}{\sqrt{4-3f^2}}$$(22) $${\lambda }_2=\frac{3\xi (f)-1}{2f}\cos (\theta ),$$(23)

where $$\zeta (f,\theta )=[\frac{2}{3}(4-3f^2-\sqrt{4-3f^2})$$(24) $${+2{\cos }^2(\theta )\left(2-f^2-\sqrt{4-3f^2}\right)]}^{1/2}$$(25)

and θ is the angle between F_r and the direction in which we are evaluating the Riemann flux.

In the free-streaming-limit (f → 1), we recover that all eigenvalues λ_1,2,3 → cos(θ) so that there are equal to ±1 along the direction of the radiative flux and zero in the direction perpendicular to it. In the diffusion limit ( f → 0), we recover ${\lambda }_{1,3}\to \mp 1/\sqrt{3}$ and λ₂ → 0 in all directions, in the Eddington limit. In our scheme, only λ_1,3 are computed. However, λ₂ is useful for understanding the structure of the Riemann problem tests that we present in Sect. 3.1. It is the propagation velocity of the equivalent of a contact discontinuity in hydrodynamics (meaning it is canceled out in the discontinuity comoving frame).

In the very optically thick limit, the diffusion timescale can be much larger than λ_1,3 leading to excessive numerical diffusion. To circumvent this issue, we adopted the method used in Sądowski et al. (2013), Rosdahl & Teyssier (2015), and Melon Fuksman & Mignone (2019) and corrected the signal speed as follows: $${\lambda }_L=\max \left({\lambda }_L,-\frac{4\widehat{c}}{3{\tau }_{\text{cell}}}\right)$$(26) $${\lambda }_R=\min \left({\lambda }_R,\frac{4\widehat{c}}{3{\tau }_{\text{cell}}}\right),$$(27)

where λ_L and λ_R are respectively the maximum and minimum signal velocities across the left and right faces, and τ_cell = ρχ∆x is the optical depth within a cell of maximum size ∆x.

2.7 Reduced speed of light

The use of an explicit integration scheme has the drawback of drastically reducing the timestep when the characteristic MHD wave speeds are much slower than the speed of light. To alleviate this issue, we used the reduced speed of light approximation (Gnedin & Abel 2001), where ĉ is an input parameter taken to be constant across the domain. A direct consequence of the reduced speed of light approximation is the violation of the usual form of the energy conservation (Melon Fuksman et al. 2021). Without gravity or dissipative source terms, the rescaled quantities that are effectively conserved are $$E_{\text{tot}}=E+(c/\widehat{c})E_r,$$(28) $$m_{\text{tot}}=\rho v+(1/\widehat{c})F_r.$$(29)

The reduced speed of light approximation does guarantee the convergence to the correct steady-state solution but can lead to artifacts in the temporal evolution leading to steady-state. To stay as close as possible to the real solution, one needs to keep the timescale ordering, i.e., the reduced speed of light should still exceed all other characteristic speeds of the system and the diffusion time should be shorter than the dynamical time (see Sect. 3.5 for more details).

2.8 Radiation-matter coupling terms

Radiation matter-coupling deals with the right hand-side of Eqs. (5) and (6). As stated above, it is customary to Lorentz transform the source terms from the fluid-frame to the laboratory-frame so that the source terms usually involve both hydrodynamic and radiation-related variables. However, since our main application is for protoplanetary disks, we decided to neglect the Lorentz transform as β ≪ 1. This means that the source terms only depend on the density and temperature for the hydrodynamic variables and not on the velocity anymore.

To solve Eqs. (5) and (6), we used an implicit method. Using Eqs. (15) and (16), we can discretize and rewrite Eqs. (5) and (6) (omitting the hyperbolic part of the equation) as $$E_r^{n+1}(1+\text{Δ}t/t_{\text{abs,red}})-\text{Δ}t/t_{\text{abs,red}}{a}_R{\left(T^{n+1}\right)}^4=E_r^n,$$(30) $$F_r^{n+1}(1+\text{Δ}t/t_{\text{scatt,red}})=F_r^n,$$(31)

where n denotes the integration step, t_abs,red ≡ 1/(ĉκ_Pρⁿ) is the characteristic timescale for radiation absorption associated with the reduced speed of light and t_scatt,red ≡ 1/(ĉχρⁿ) is the typical timescale for radiation scattering associated with the reduced speed of light.

Equation (31) is trivially solved, but Eq. (30) requires an equation on Tⁿ⁺¹ to close the system. To this end, we use the conservation of internal energy, $$C_v{\partial }_{t}T=c{\kappa }_{\text{P}}\rho \left(E_r-a_{\text{R}}{T}^4\right),$$(32)

which we discretize as $$C_v{T}^{n+1}+\text{Δ}t/t_{\text{abs}}{a}_R{\left(T^{n+1}\right)}^4-\text{Δ}t/t_{\text{abs}}{E}_r^{n+1}=C_v{T}^n,$$(33)

where t_abs ≡ 1/(cκ_Pρⁿ) and C_v is the calorific capacity.

We then followed Hayes & Norman (2003) and Commerçon et al. (2011) and assumed that T does not change significantly over one time step so that we can write (Tⁿ⁺¹)° = 4Tⁿ⁺¹(Tⁿ)° − 3(Tⁿ)°. We can then rewrite Eqs. (30) and (32) in the following way: $$\left[\begin{array}{cc}1+{\gamma }_{\text{abs,red}}& -4{\gamma }_{\text{abs,red}}{a}_R{T}^{n3}\\ -{\gamma }_{\text{abs}}& C_v+4{\gamma }_{\text{abs}}{a}_R{T}^{n3}\end{array}\right]\left[\begin{array}{c}{E}_r^{n+1}\\ T^{n+1}\end{array}\right]=\left[\begin{array}{c}{E}_r^n-3{\gamma }_{\text{abs,red}}{a}_R{T}^{n4}\\ C_v{T}^n+3{\gamma }_{\text{abs}}{a}_R{T}^{n4}\end{array}\right],$$(34)

where γ_abs,red ≡ ∆t/t_abs,red and γ_abs ≡ ∆t/t_abs. To find $E_r^{n+1}$, we inverted this 2° matrix for each cell, which we could do analytically. Finally, to complete the implicit step, we computed the new gas energy and momentum in each direction by using the rescaled conservation of energy, Eq. (28). While the implicit part of our scheme is currently written for single group radiative transfer, our implementation can trivially be extended to multi-group by solving for a (m + 1)° matrix instead of a 2°, where m is the number of frequency groups (see Sect. 2.8).

We preferred to use this implicit solver rather than the fixed point method used in Melon Fuksman et al. (2021) because we find that it is more robust (less prone to overshooting or undershooting the solution) and performs better since it is not an iterative method. For completeness, we also implemented the fixed-point method described in Melon Fuksman & Mignone (2019) but refer the reader to the aforementioned paper for details as we did not use this method in the tests presented here.

2.9 Irradiation flux

To treat irradiation by an external object, we allow the user to define an irradiation flux, the divergence of which acts as a source term in the energy conservation equation, Eq. (4). This irradiation flux is not part of the radiation scheme; it is defined through a user function that updates instantaneously the irradiation flux throughout the entire domain at each time step. We find that the best way to implement this term is to incorporate it to the implicit scheme. To do so, we rewrite Eq. (28) to take into account the energy deposited by the irradiation flux so that $$E_{\text{tot}}=E+(c/\widehat{c})E_r-\nabla \cdot F_{\text{irr}}\times \text{Δ}t,$$(35)

where ∆t is the timestep. We then also rewrite the bottom term of the right hand side of Eq. (34), S₁, to $$S_1=C_v{T}^n+3{\gamma }_{\text{abs}}{a}_R{T}^{n4}-\nabla \cdot F_{\text{irr}}\times \text{Δ}t,$$(36)

solve for this new system to find $E_r^{n+1}$, and finally use Eq. (35) to find Eⁿ⁺¹. We find that adding the irradiation source terms in the implicit step provides a much better accuracy compared to a split-scheme, where irradiation would be added afterward, and prevents an overshooting of the solution in the test case (see Sect. 3.7).

3 Radiation test problems

In the first two tests, the optically thin “shock tubes” and the free-streaming beam, the choice of ĉ does not impact the results as we test only the hyperbolic part of Eqs. (5) and (6). In the third, fourth, and seventh tests, the shadow tests, the pulse tests, and the stellar irradiation tests, the hydrodynamic fluid is static, so the Courant–Friedrichs–Lewy condition is only determined by the radiation subsystem. In such cases, the choice of ĉ is also without consequences. Hence, the only tests where the choice of ĉ matters are the fifth and sixth tests, the radiative shocks and vertical diffusion tests.

3.1 Optically thin “shock tube” tests

To test the implementation of the hyperbolic solver, we performed two tests of 1D radiative shock tubes in Cartesian coordinates following Melon Fuksman & Mignone (2019). The two tests follow the evolution of two constant states L and R, located respectively to the left and right of x = 0. We neglected any interaction with matter by setting all opacities to zero. The simulation box extends from x ∈ [−20, 20].

The first test is a discontinuity along the x-direction in the tangential radiative flux. Specifically, we set as E_r,L = E_r,R = 1, F_r,x,L = F_r,x,R = 0 and F_r,y,L = 1/2 and F_r,y,R = 0. The radiation field has f = 1/2 so that it is neither in the free-streaming limit nor in the diffusion limit. We show in Fig. 1 the solution obtained at t = 20 for the HLL solver with two different reconstruction schemes, the f -preserving PLM (blue line) and f -preserving LimO3 (green line), at a resolution of 2°. We also show the result of a simulation with HLL and a f -preserving PLM reconstruction with a resolution of 2¹⁷ for reference (black line). We see in the top three panels the development of a three-wave pattern, with (1) a left-facing shock on the left, (2) a contact-like discontinuity in the middle, and (3) a rightward expansion wave on the right. The two bottom panels show that, at the highest resolution, the fields β_x ≡ (3ξ − 1)F_r,xE_r/2||F_r|| and Π ≡ (1 − ξ)E_r/2, which respectively are analogs of the hydrodynamic velocity and pressure, are constant across the contact-like discontinuity. In this test, we see that in contrast to PLM, LimO3 tends to overshoot and undershoot the solution at the discontinuity as can be seen from the right panels that are zoomed on x ≈ 11.2. We find that this behavior is quite general to all tests and can lead to negative radiative energy in situations where E_r is very small on one side of the discontinuity. In such cases, we needed to manually enforce the positivity of the radiative energy limiting the usefulness of higher-order reconstruction schemes.

The second test is that of a shock in the x direction. Specifically, we set E_r,L = 1/10, E_r,R = 1, F_r,x,L = 1/10 F_r,x,R = 0, F_r,y,L = 0, and F_r,y,R = 1. The radiative fluid has f = 1 on both sides so that we test the ability of the solver to handle the free-streaming limit. Again, we show in Fig. 2 the solution at t = 20 using the HLL solver but with different reconstruction schemes, PLM and LimO3 (both f -preserving), at a resolution of 2°. We also show the result of a simulation with a PLM f -preserving reconstruction with a resolution of 2¹⁷ and that of a simulation with a LimO3 with the PLUTO reconstruction scheme for reference (see Sect. 2.5). Again, we see in the top three panels the development of a three-wave pattern, with (1) a left-facing shock on the left, (2) a contact-like discontinuity in the middle and (3) a right-going shock on the right. This test illustrates the need for a careful treatment of the reconstruction on the faces in the free-streaming limit. Indeed, we see in the two bottom panels that our solution with the reconstruction scheme of PLUTO produces spurious oscillations on the left side of the contact-like discontinuity when using LimO3. We checked that this behavior is also observed using the PLUTO code. Our reconstruction scheme does not produce such spurious oscillations; however, it does produce a small peak in Π on the left side of the contact-like discontinuity with both PLM and LimO3. Despite this slight loss of accuracy with PLM compared to the reconstruction scheme of PLUTO, we chose to adopt our reconstruction scheme given its improvement in the solution with LimO3.

In Fig. 3, we show the result of a resolution study for both Riemann tests using second-order PLM and third-order LimO3 reconstructions schemes spanning resolutions going from 2° to 2°4. We plot the L1-norm error of E_r computed compared to the solution of reference with a resolution of 2¹⁷. We use a Courant number of 0.4 and a RK2 time integrator. We see that both tests scale as N^−1/2, which is expected for shock solutions (LeVeque 2002).

Fig. 1 First test of an optically thin shock tube. The blue and green lines show the results of simulations with 2° radial cells using a PLM reconstruction and LimO3 reconstruction, respectively. The black line shows the reference solution using PLM reconstruction with 217 radial cells. Left panels: full solution. Right panels: zoomed-in view of the left-facing shock at x ≈ 11.2.

Fig. 2 Second test of an optically thin shock tube. The blue and green lines show the results of simulations with 2° radial cells using a PLM reconstruction and a LimO3 reconstruction, respectively, with our reconstruction of f on the faces. The red line shows the result of a simulation with 2° radial cells using a PLM reconstruction with PLUTO’s reconstruction of f on the faces. The black line shows the reference solution using PLM reconstruction with 217 radial cells.

3.2 Free streaming beam

To test the multidimensional hyperbolic transport of our radiative scheme, we implement an oblique free streaming radiation beam test (Richling et al. 2001; González et al. 2007). We set all opacities to zero so that there is no interaction between radiation and matter. The domain consists of a square Cartesian grid of size L = 5 in code units. We initialize the background radiation field as E_r = 10° and F_r, x = F_r, y = 0 in code units. The beam is initialized between x ∈ [0.5, 0.6] and y ∈ [0.3, 0.44] with E_r = 10¹² and $F_{r,x}=F_{r,y}=E_r/\sqrt{2}$. The boundary conditions are outflow everywhere. Note that we do not inject the beam from the boundary condition as in Melon Fuksman & Mignone (2019) to avoid the implementation of a special boundary condition for oblique injection. Finally, we use a resolution of 300 × 300, a Courant factor of 0.4, the HLL scheme and different reconstruction schemes.

We show in the first panel of Fig. 4 a color plot of this test using the PLM f -preserving reconstruction. As expected, the beam is propagating in a straight line and reaches the outer boundary. As it is propagating the beam widens due to numerical diffusion. To quantify this diffusion, we plot, in the bottom panel of Fig. 4, horizontal slices of the beam at x = 1 and x = 4.35 that are centered on the maximum value of the profile. We see that the beam maximum decreased by 8% and the beam’s full width at half maximum increased by only 7% during its propagation. This is much lower than what is reported in Melon Fuksman & Mignone (2019) where increase by as much as 40% although we do not use exactly the same metric.

In Fig. 5, we also illustrate the importance of the choice of reconstruction scheme with a high-resolution test of the free streaming beam with a resolution of 1000 × 1000. All panels use the LimO3 flux-reconstruction but enforce the physical condition f ≤ 1 in different ways. The first panel shows a reconstruction scheme where we simply enforce that f ≤ 1 on cell-faces. More precisely, we set $$E_{r,f}=ℛ(E_{r,c}),$$(37) $$F_{r,f}^i=ℛ(F_{r,c}^i)\text{if\hspace{0.17em}}\frac{\parallel ℛ(F_{r,c}^i)\parallel }{ℛ(E_{r,c})}\le 1,$$(38) $$F_{r,f}^i=ℛ(F_{r,c}^i)\frac{ℛ(E_{r,c})}{\parallel ℛ(F_{r,c}^i)\parallel }\text{otherwise.}$$(39)

We see that with this simple reconstruction scheme the beam “explodes.” We believe that this is due to the accumulation of small deviation between f on the center and faces of the cells that introduce deviation from the free-streaming limit as discussed in Sect. 2.5. The middle panels show the results of the free streaming beam test using the “PLUTO reconstruction” scheme and our f -preserving reconstruction scheme. We see that both of these methods fix the beam explosion and give similar results although the PLUTO reconstruction produces spurious features around the beam that are absent with our reconstruction.

Fig. 3 L1 norm error as a function of the number of cells for the first optically thin Riemann problem (top) and the second optically thin Riemann problem (bottom) computed from a reference solution with 217 radial cells. Blue and green points show simulations with PLM and LimO3 reconstructions, respectively.

Fig. 4 Top: color map of the radiation energy density for the free streaming beam with a resolution of 300 × 300 with a LimO3 f -preserving reconstruction scheme. Bottom: vertical cut of the radiation energy density at x = 1 and x = 4.35.

3.3 Shadow

The M1 method offers a significant advantage over the FLD approach by retaining directional information about the radiation field, thereby enabling an accurate representation of shadowing effects. To validate this capability, we conducted the benchmark test proposed in Hayes & Norman (2003), González et al. (2007), and Melon Fuksman & Mignone (2019), which consists of illuminating an opaque obstacle from one side and examining the resulting shadow formed in the downstream region. The setup consists of a rectangular box of dimension L × l = 1 cm × 0.3 cm with a resolution of 140 × 40. The radiation beam is injected from the left boundary along the x-direction with a radiation temperature of T_r = 1740 K and a flux along the x-direction F_{r, x} = E_r corresponding to the free-streaming limit. The right hand side horizontal boundary and the top vertical boundary are outflowing while the bottom vertical boundary, at y = 0, is reflective. The opaque surface is a spheroid of density ρ₁ = 10° g cm⁻³ surrounded by a background with a density of ρ₀ = 1 g cm⁻³. To allow for a smooth transition between ρ₀ and ρ₁, we set $$\rho ={\rho }_0+\frac{{\rho }_1-{\rho }_0}{1+e^{\text{Δ}}},$$(40)

where $$\text{Δ}=10\left[{\left(\frac{x}{x_0}\right)}^2+{\left(\frac{y}{y_0}\right)}^2-1\right],$$(41)

with (x₀, y₀) = (0.1, 0.06). The entire domain is set in thermal equilibrium with T = T_r = 290 K. Additionally, the radiative fluxes and gas velocities are initially set to zero. Finally, we set the absorption Rosseland opacity using Kramers’ law so that κ = 0.1(ρ/ρ₀)(T/T₀)^−3.5 cm° g⁻¹. With this choice of opacity, we ensure that the spheroid is optically thick to radiation while the rest of the domain is optically thin.

We show in Fig. 6 two color maps of the radiation energy after 10 light crossing time, by which time the beam has long passed the opaque spheroid, reached the rightmost boundary and settled to a steady state. In the top panel, we follow Skinner & Ostriker (2013) and remove the source terms responsible for emission of radiative energy by the gas so that radiation can only be absorbed. In practice, this means that the radiation-matter interaction reduces to $$E_r^{n+1}=\frac{E_r^n}{1+{\gamma }_{\text{abs,red}}}$$(42)

for the radiative energy instead of Eq. (34). The top panel of Fig. 6 confirms the expected formation of a sharp shadow behind the spheroid, with the transition between the shadowed and illuminated regions confined to a single cell. In the middle panel, we present the case where both emission and absorption terms are included so that we are solving for Eq. (34). Here, the shadow appears more diffuse, and an overshoot in radiation energy is observed at the transition between the shadow and the beam. Additionally, the radiation energy within the shadow is higher when emission terms are accounted for, as the gas located behind the spheroid contributes to the local light emission.

We attribute the observed over-density in radiation energy to the interaction between the incident beam and the emitted light originating from behind the shadowed region. Overall, our results are quantitatively consistent with those reported in González et al. (2007), Skinner & Ostriker (2013), and Melon Fuksman & Mignone (2019).

Fig. 5 Color maps of the radiation energy density for the free streaming beam test run using different reconstruction scheme. Top: reconstruction scheme where we only impose f ≤ 1 on the faces after reconstruction. Middle: PLUTO radiative reconstruction scheme. Bottom: our f -preserving reconstruction scheme.

Fig. 6 Color map of the radiation energy density without emission source terms (top) and with emission source terms (middle) for the shadow test run. Bottom: vertical slice of the radiation energy density at x = 0.5 for the run with emission (yellow line and dots) and without emission (blue line and dots).

3.4 Pulses

We next tested the ability of our scheme to handle transport of radiation in the optically thin limit (but with matter-light interactions on) in Cartesian and spherical coordinates by performing a pulse propagation test. We initialized a spherically symmetric radiation energy distribution as E_r = 4πB(T_r), where B(T) is the blackbody function and $$T_r=T_0\left(1+100e^{-r^2/w^2}\right),$$(43)

where r is the spherical radius, T₀ = 1, 05 × 10⁻⁷ K, and w = 5 in code units. The gas is initially at a temperature of T₀ so that it is in thermal equilibrium with the radiation away from the pulse. We set F_r = 0, ρ = 1 in code units, κ = 0, Γ = 5/3, C_a = 0.4. Our units are unit_v = c, unit_ρ = 1.67 × 10⁻²⁴ g cm⁻³ and unit_L = 1.496 × 10¹³ cm.

To ensure an optically thin domain, we set a constant scattering opacity of σ = 3.9 × 10° cm° g⁻¹. The top panel of Fig. 7 shows the energy density of the optically thin pulse at t = 35 in code units when propagating on 3D on a Cartesian grid of 200×200×200 with (x, y, z) ∈ [−50, 50]×[−50, 50]×[−50, 50]. We see that the pulse is spherically expanding as expected. For this resolution, we do not see, with the naked eye, artifacts along each axis due to the Cartesian grid as can be seen in Melon Fuksman & Mignone (2019). In the bottom panel, we plot two 1D profiles of the radiation energy density for the 3D Cartesian simulation, one along the x-axis at y = 0 and one along the x = y axis. We also plot the radiation energy density from a 1D spherical simulation for comparison. We see that the spherical simulation decays as 1/r° as expected. The profile of the pulse along the x = y axis in the Cartesian run is almost indistinguishable from the spherical profile. However, along the x-axis of the Cartesian simulation the pulse has a larger amplitude than the spherical case. This is an artifact due to propagation along the grid. For consistency, we checked that the total energy of the pulse is the same in the two simulations and that in each case it is conserved to machine precision.

3.5 Shocks

Radiative shocks differ from hydrodynamical shocks as radiation is able to redistribute the thermal energy deposited in the post-shock region, which is heated to a temperature of T₂, back into the pre-shock region. In an optically thick medium, the pre-shock region will rise to a temperature T₋ ≤ T₂ before reaching the front. As a result, the temperature on the shocked side of the front is T₊ > T₂ (see Fig. 8) and relax to T₂ away from the front. There are two types of radiative shocks depending on the strength of the shock: (1) subcritical shock, where the post-shock temperature T₂ is not high enough for radiation to dramatically affect the pre-shock region so that T₋ < T₂, (2) supercritical shock, where T₂ is so large that the radiation escaping from the post-shock region is able to heat the precursor to T₋ = T₂.

We perform two radiative shocks simulations, a subcritical one and a supercritical one, to test the dependence of our results with the reduced speed of light. Indeed, since radiative shocks are dynamic problems involving light-matter interaction, they provide a good illustration of how our choice of reduced speed of light can affect the final solution. As noted in Sect. 2.7, the use of the reduced speed of light approximation introduces an artificial loss of energy such that $$\text{Δ}(E+E_r)=(1-\widehat{c}/c)\text{Δ}{E}_r,$$(44)

which is particularly important in problems where the radiative energy changes suddenly such as in the case of radiative shocks.

To initialize the simulation, we reproduced the setup of Ensman (1994), Hayes & Norman (2003), González et al. (2007), and Melon Fuksman et al. (2021). We initialized a 1D Cartesian grid extending from x ∈ [0, 7 × 10¹⁰ cm]. We initialized a constant density of ρ = 7.78 × 10⁻¹⁰ g cm⁻³ and a constant pressure and radiation field using a temperature of T = 10 K, µ = 1 and Γ = 7/5. We used a constant absorption opacity of κ = 0.39 and a null scattering opacity so that the domain is optically thick to absorption opacity with ρκL ≈ 20. A rightward moving shock was created by setting v_x < 0 and a reflective boundary condition on the left boundary. Our choice of v_x determines the shock regime. We used v_x = −6 and −20 km s⁻¹ for the subcritical and supercritical regime, respectively. Finally, we used the HLL Riemann solver with the PLM f -preserving reconstruction.

We see in Fig. 8 that in both cases for ĉ/c = 1, we retrieve a very similar solution as González et al. (2007) and Melon Fuksman et al. (2021), with T₂ = 834 K, T₋ = 358 K, and T₊ = 1080 K for the subcritical shock and T₂ = 4262 K and T₊ = 5314 K for the supercritical shock. We also see that our choice of ĉ/c critically affect the solution, especially for the supercritical shock. For the subcritical shock, we find T₂ = 800K for ĉ/c = 10⁻⁴ quite close to the solution with ĉ/c = 1. However, for the supercritical shock, we find T₂ = 2720 K for ĉ/c = 10⁻⁴ much lower than for ĉ/c = 1. Interestingly, we see that a naive prediction of the maximum ĉ we can use for the subcritical shock would lead to ĉ ≫ (|v_x| + c_s) ≈ 9 × 10⁻⁵ while our departure from the ĉ solution is already very large at ĉ = 10⁻³. However, as noted by Melon Fuksman et al. (2021), this simple estimate does not take into account the timescale on which energy is injected into the system. We propose that a better estimate of ĉ is given by $\max (E_r/\left|\rho v_x^3\right|)$, which gives the typical velocity at which radiation can escape in order to compensate for kinetic energy injection. We find that at the shock location of our simulation with ĉ/c = 1, this estimate gives ĉ ≫ 1230 and ĉ ≫ 1 for the subcritical and supercritical shocks. This is consistent with the fact that the subcritical shock solution departs from the reference solution only for ĉ/c = 10⁻⁴ while the supercritical shock departs from the reference solution for ĉ/c = 10 already. This highlights the importance of examining all relevant timescales of the simulations and of performing convergence tests with different reduced speeds of light.

Fig. 7 Top: color map of the radiation energy density for the optically thin 3D Cartesian pulse test. Bottom: slices along the x axis and x = y axis for the 3D Cartesian pulse test (solid and dash-dotted lines, respectively) at different times compared to the 1D spherical case (dashed lines). The dotted black line shows the expected decrease in the energy density as 1/r°.

Fig. 8 Gas temperature (solid lines) and radiation temperature (dashed lines) as a function of x for the subcritical radiation shock (top) and supercritical shock (bottom) for different choices of reduced speed of light ĉ/c = 1, 10−2, 10−3, and 10−4 as blue, salmon, yellow, and brown lines, respectively.

3.6 Vertical diffusion

We perform a test of radiative energy diffusion along the vertical extent of a disk under constant injection of viscous energy (see Melon Fuksman et al. 2021). With this test, we check the ability of our code to perform in the diffusion limit. We also tested the diffusion of an optically thick pulse (as in Melon Fuksman & Mignone 2019) but do not present it here as the two tests are testing similar features of the code.

We model a 1D vertical slice of a disk in hydrostatic equilibrium as $$\rho (x)={\rho }_0\exp \left(-x^2/2H^2\right)+{\rho }_{\text{min}},$$(45)

where ρ₀ = 10⁻¹⁰ g cm⁻³ is the midplane density, ρ_min = 10⁻¹⁰ρ₀ is the density floor, H = 0.05R is the pressure scale height and x ∈ [−1, 1] is the height of the disk. The pressure is initialized using a constant temperature T₀ = 1000 K. We used a mean molecular weight µ = 2.35 and an adiabatic index Γ = 1.41. We used a resolution of 201 points in x, an HLL solver with PLM f -preserving reconstruction, and a RK2 time-integrator. Finally, we set local thermodynamic equilibrium with T = T₀ at the boundaries for E_r as well as zero-gradient boundary conditions on F_r.

We set all velocity components to zero and set the mass and momentum fluxes to zero so as to keep the density and velocity structure fixed. Hence, the energy is evolved only under the influence of radiative cooling and viscous heating defined as $$S_E=\frac{9}{4}\alpha {\Omega }_K{c}_s^2\rho ,$$(46)

where α = 10⁻³, Ω_K is the Keplerian angular velocity at 5 au, and c_s is the isothermal sound speed.

We plot in Fig. 9, the semi-analytical solution to this problem, which we computed following the procedure of Melon Fuksman et al. (2021), as well as solutions for ĉ/c = 1, 10⁻² and 10⁻⁴. We see that our solution is very close to the analytical solution for all ĉ/c. However, we find that find cell-centered variables always exhibit a flattening of the radiative flux near the origin that can be clearly seen in the right panels of Fig. 9. This flattening of the radiative flux is actually an artifact due to the reconstruction of cell-centered variables on the cell faces when solving the Riemann problem. Indeed, we plot the radiative flux from the Riemann problem and find that it does not show a flattening. We also plot in the bottom panel of Fig. 9 the relative L1 norm error on the radiative flux as a function of time. We see that the final error on our solution is identical for all of our values of reduced speed of light and is around 10⁻³. The only difference between simulations using different reduced speed of light is the time they take to get to the final solution. We note that in this test, although we do inject energy all the time as in the radiative shock test, the characteristic time of energy injection is min(E_r/S_E) ≈ 5 × 10° s, which is ten times slower than the typical light crossing timescale for ĉ/c = 10⁻⁴.

Fig. 9 Top left and center left: radiative energy and radiative flux as a function of x for the vertical diffusion test. Empty circles show the cell-centered values. Crosses show the values of radiative flux from the Riemann solver and radiative energy reconstructed from the Riemann fluxes of energy and flux. The black line shows the analytical solution. Top right and center right: difference between our simulation and the analytical solution. Bottom: L1 relative norm for three different reduced speed of light ĉ/c = 10−2, 10−3, and 10−4 as blue, red, and yellow lines, respectively.

3.7 Stellar irradiation

In this last test, we checked the ability of our code to handle irradiation of a passive disk by a central star (Pascucci et al. 2004). We defined our disk density as $$\rho (R,z)={\rho }_0\left(\frac{500\text{\hspace{0.17em}AU}}{R}\right)\exp \left(-\frac{\pi }{4}{\left(\frac{z}{h(r)}\right)}^2\right),$$(47)

where (R, z) = r(cos θ, sin θ) with r ∈ [1, 1000] AU and θ ∈ [0, π] and where h(R) = 125 AU × (R/500 AU)^1.125. We used a spherical grid with 240 cells spaced logarithmically in r and 100 cells spaced uniformly in θ. We used outflowing boundary conditions in all directions and used an HLL solver with PLM reconstruction and a RK2 time integrator. We used a reduced speed of light of ĉ/c = 10⁻⁴ and ran the simulation for 10 years of physical time until the temperature reached equilibrium out to R ≈ 300 AU. We set all velocity components to zero and we cancel the fluxes of density, momentum, and energy so that the disk temperature is only changed by the irradiation source term and the absorption source term.

We computed the opacities for absorption and irradiation from the opacities provided by RADMC3D. We used the “astronomical silicate” opacity of Draine (2003), with a grain size of 0.12 µm and a density of 3.6 g cm⁻³. We set ρ₀ = 6.66 × 10⁻¹⁷ g cm⁻³ so that the optical depth of a ray at 550 nm going through the midplane experiences an optical depth of τ = 100. From these opacities as a function of wavelength, we pre-computed 1D Planck and Rosseland opacity tables as a function of temperature that we use to compute the opacity by interpolation in each cell. We neglect scattering by setting σ = 0.

The irradiation flux is given by $$F_{\text{irr}}(r,\theta )=\pi {\left(\frac{R_s}{r}\right)}^2{\displaystyle {\int }_{{\nu }_{\text{min}}}^{{\nu }_{\text{max}}}d\nu B_{\nu }(T_s)e^{-{\tau }_{\text{irr}}(r,\theta ,\nu )}{r}_{,}}$$(48)

where R_s = R_⊙ is the radius of the central star,ν_min = 1.5× 10¹¹ Hz and ν_max = 1.5 × 10¹⁵ Hz, B_ν is the Planck function and T_s = 5800 K is the temperature of the central star. We define the optical depth seen by the irradiation flux as $${\tau }_{\text{irr}}(r,\theta ,\nu )={\displaystyle {\int }_{r_{\text{in}}}^r\kappa (\nu )\rho (r^{\prime },\theta )d{r}^{\prime }.}$$(49)

In practice, we pre-compute tables of the integral on frequency in Eq. (48) as a function of τ_irr from which we interpolate the flux in each cell.

We show in Fig. 10 a color map of the temperature as well as a cut in the midplane as a function of radius and a cut at a radius of 2 AU as a function of θ. We retrieve the typical structure of an irradiated disk, with hot upper layers and a cold disk midplane where stellar irradiation gets absorbed as it propagates radially outward. We compare our solution with a multifrequency multidimensional Monte Carlo simulation done with RADMC3D for the same grid and the same opacity tables using 10° photons. We see that the agreement is quite good between our simulation and the Monte Carlo simulation. We typically overestimate the temperature near the pole, because of beam crossing at the axis, and we underestimate the temperature in the mid-plane, but these deviations between the two methods do not exceed 5%. Nonetheless, we find that at higher optical depth the disagreement between the two method gets worse. We plot in Fig. 11 cuts of the temperature at 2 AU for increasing dust to gas ratios (so increasing optical depth). We see that as we increase the optical depth, the midplane temperature flattens and eventually gets a “Mexican hat” shape. This problem was raised in Melon Fuksman et al. (2025) where it was attributed to beam crossing at the midplane. In Melon Fuksman et al. (2025), the authors proposed a modified M1 method, called the half moment method, where the equations are integrated over hemispheres instead of the full solid angle when computing the moment of the radiative transfer equations. However, implementing this method is outside the scope of this paper.

Fig. 10 Top: color map of the gas temperature as a function of r and z/r for the stellar irradiation test. Middle and bottom: comparison of the temperature for the Monte Carlo and M1 solution in the midplane as a function of r (middle) and at r = 2 AU as a function of θ (bottom).

4 Performance

To test the performance of our radiative version of IDEFIX, we perform the pulse test described in Sect. 3.4 on a 3D Cartesian grid. We let the simulation run for a thousand cycles without performing any outputs. We use the HLL solver for radiation and the hydrodynamics with the PLM f -preserving scheme and RK2 time-integration. When we increase the number of nodes, we increase proportionally the size of the domain as well as the number of points in the simulation. In this way, we maintain a constant time step in the simulation and so a fixed number of cycles.

We measure the performance of IDEFIX on the French supercomputer AdAstra on two types of GPUs, AMD Mi250X and AMD Mi300. We plot in the top panels of Fig. 12 the number of cell updates per second per node as a function of the number of nodes. We see that for the largest sub-domains of 256°, the performance goes as high as 7.6 × 10° cell updates per second per node on Mi250X and 1.2 × 10° cell updates per second per node on Mi300. This compares to 1.24 × 10° cell updates per second per node on Mi250X for a test run on MHD (see Lesur et al. 2023). This shows that the radiative transfer scheme is roughly 1.6 times more expensive than the MHD scheme, because of the complexity of the operations needed to compute the closure in the former case. Interestingly, we see that although the performance on the Mi300 GPUs is better than on the Mi250X GPUs for sub-domains of 256°, it is the contrary for sub-domains of 64° and 32°. This behavior arises because IDEFIX employs two MPI processes per GPU on the Mi250 (one process per graphic compute die), effectively increasing the domain size per GPU. As a result, the Mi250 can better hide memory access latency compared to the Mi300 for small domain sizes.

We also plot the weak scaling of our simulations in the lower panels of Fig. 12. We see that on both Mi250x and Mi300 our radiative version of IDEFIX performs with almost perfect scaling, staying higher than 90% even when running on 64 nodes. This is consistent with the test presented on the MHD version of IDEFIX and is quite natural given how close is the implementation of the radiative module and the hydrodynamic module.

Fig. 11 Temperature at r = 2 AU as a function of θ for four different dust-to-gas ratios of 10−2, 10−1, 1, and 10, shown as brown, yellow, salmon, and blue lines, respectively. Dashed lines show the results from the Monte Carlo simulations, and solid lines show the results from the M1 simulations.

Fig. 12 Top: performance-in-cell update per second per cell on AdAstra’s Mi250x and Mi300 as a function of nodes. Bottom: weak scaling efficiency.

5 Conclusion

In this paper, we have presented the new radiative transfer module of IDEFIX. We used the M1 approximation, effectively treating radiation as a fluid, which is evolved using a high-order finite-volume Godunov method. We solved for the hyperbolic part of the equations using an explicit method with radiation Riemann solvers. To compute the reconstructed inter-cell fluxes, we propose an improved reconstruction scheme that provides more accurate results in the free-streaming limit. Because the scheme is nonrelativistic, we used a reduced speed of light to reduce the difference between the hydrodynamic and radiative timescales. The source terms are treated implicitly by inverting a matrix of (n + 1)°, where n is the number of frequency groups. For now, the scheme uses only one frequency group, but we intend to use several in future work. Radiation can be used on Cartesian, cylindrical and spherical grids that are uniform or stretched and in multiple dimensions. The scheme is also compatible with the MHD module.

The radiative transfer version of IDEFIX performs well, with as many as 7.6 × 10° cell updates per second per node on an AMD Mi250x on AdAstra. It keeps a weak scaling efficiency of 90% up to 64 nodes. As such, radiation only costs 1.6 times an MHD simulation.

Acknowledgements

The authors acknowledge support from the European Research Council (ERC) under the European Union Horizon 2020 research and innovation program (Grant agreement No. 815559 (MHDiscs)). This work was supported by the MHD@Exascale project (reference 22-EXOR-0015) of PEPR Origins (PI: Morbidelli). This project was provided with computing HPC and storage resources by GENCI at CINES thanks to the grant 2025-A0180402231 on the supercomputer Adastra Mi250 and Mi300 partitions.

References

All Figures

	Fig. 1 First test of an optically thin shock tube. The blue and green lines show the results of simulations with 2° radial cells using a PLM reconstruction and LimO3 reconstruction, respectively. The black line shows the reference solution using PLM reconstruction with 217 radial cells. Left panels: full solution. Right panels: zoomed-in view of the left-facing shock at x ≈ 11.2.
In the text

Fig. 2 Second test of an optically thin shock tube. The blue and green lines show the results of simulations with 2° radial cells using a PLM reconstruction and a LimO3 reconstruction, respectively, with our reconstruction of f on the faces. The red line shows the result of a simulation with 2° radial cells using a PLM reconstruction with PLUTO’s reconstruction of f on the faces. The black line shows the reference solution using PLM reconstruction with 217 radial cells.

In the text

	Fig. 3 L1 norm error as a function of the number of cells for the first optically thin Riemann problem (top) and the second optically thin Riemann problem (bottom) computed from a reference solution with 217 radial cells. Blue and green points show simulations with PLM and LimO3 reconstructions, respectively.
In the text

	Fig. 4 Top: color map of the radiation energy density for the free streaming beam with a resolution of 300 × 300 with a LimO3 f -preserving reconstruction scheme. Bottom: vertical cut of the radiation energy density at x = 1 and x = 4.35.
In the text

	Fig. 5 Color maps of the radiation energy density for the free streaming beam test run using different reconstruction scheme. Top: reconstruction scheme where we only impose f ≤ 1 on the faces after reconstruction. Middle: PLUTO radiative reconstruction scheme. Bottom: our f -preserving reconstruction scheme.
In the text

	Fig. 6 Color map of the radiation energy density without emission source terms (top) and with emission source terms (middle) for the shadow test run. Bottom: vertical slice of the radiation energy density at x = 0.5 for the run with emission (yellow line and dots) and without emission (blue line and dots).
In the text

	Fig. 7 Top: color map of the radiation energy density for the optically thin 3D Cartesian pulse test. Bottom: slices along the x axis and x = y axis for the 3D Cartesian pulse test (solid and dash-dotted lines, respectively) at different times compared to the 1D spherical case (dashed lines). The dotted black line shows the expected decrease in the energy density as 1/r°.
In the text

	Fig. 8 Gas temperature (solid lines) and radiation temperature (dashed lines) as a function of x for the subcritical radiation shock (top) and supercritical shock (bottom) for different choices of reduced speed of light ĉ/c = 1, 10−2, 10−3, and 10−4 as blue, salmon, yellow, and brown lines, respectively.
In the text

Fig. 9 Top left and center left: radiative energy and radiative flux as a function of x for the vertical diffusion test. Empty circles show the cell-centered values. Crosses show the values of radiative flux from the Riemann solver and radiative energy reconstructed from the Riemann fluxes of energy and flux. The black line shows the analytical solution. Top right and center right: difference between our simulation and the analytical solution. Bottom: L1 relative norm for three different reduced speed of light ĉ/c = 10−2, 10−3, and 10−4 as blue, red, and yellow lines, respectively.

In the text

	Fig. 10 Top: color map of the gas temperature as a function of r and z/r for the stellar irradiation test. Middle and bottom: comparison of the temperature for the Monte Carlo and M1 solution in the midplane as a function of r (middle) and at r = 2 AU as a function of θ (bottom).
In the text

	Fig. 11 Temperature at r = 2 AU as a function of θ for four different dust-to-gas ratios of 10−2, 10−1, 1, and 10, shown as brown, yellow, salmon, and blue lines, respectively. Dashed lines show the results from the Monte Carlo simulations, and solid lines show the results from the M1 simulations.
In the text

	Fig. 12 Top: performance-in-cell update per second per cell on AdAstra’s Mi250x and Mi300 as a function of nodes. Bottom: weak scaling efficiency.
In the text