$\texttt{GPUmonty}$: A GPU-accelerated relativistic Monte Carlo radiative transfer code

TL;DR

By significantly reducing computational costs, GPUmonty enables the extensive parameter space surveys and faster spectra modeling required to interpret horizon-scale observations of supermassive black holes.

Abstract

We introduce $\texttt{GPUmonty}$, a CUDA/C-based Monte Carlo radiative transfer code accelerated using graphics processing units (GPUs). $\texttt{GPUmonty}$ derives from the CPU-based code $\texttt{grmonty}$ and offloads the most computationally expensive stages of the calculation -- superphoton generation, sampling, tracking, and scattering -- to the GPU. Whereas $\texttt{grmonty}$ handles photons sequentially, $\texttt{GPUmonty}$ processes large numbers of superphotons concurrently, leveraging the single-instruction, multiple-thread (SIMT) execution model of modern GPUs. Benchmarks demonstrate a speedup of about $12\times$ relative to the original CPU implementation on a single GPU, with runtime limited primarily by register pressure rather than compute or memory bandwidth saturation. We validate the implementation through analytic tests for a optically thin synchrotron sphere, as well as comparisons with $\texttt{igrmonty}$ for scattering synchrotron sphere and GRMHD simulation data. Relative errors remain below a percent level and convergence is consistent with the expected $N_{\rm s}^{-1/2}$ Monte Carlo scaling. By significantly reducing computational costs, GPUmonty enables the extensive parameter space surveys and faster spectra modeling required to interpret horizon-scale observations of supermassive black holes. $\texttt{GPUmonty}$ is publicly available under the GNU General Public License.

Figures (6)

• Figure 1: Panel (a): Optically thin synchrotron sphere spectrum considering $N_{\rm s} = 10^8$. The red line represents the synthetic spectrum generated with GPUmonty and the purple markers are the analytical results from the angle-integrated emissivity calculated with equation \ref{['eq:Lnu_analytic_thin_test']}. Panel (b): Residuals computed with equation \ref{['eq:relative_error']}.
• Figure 2: Normalized integrated error for different values of $N_{\rm s}$, computed with equation \ref{['eq:conv_parameter']}. The blue dashed line represents the convergence $\epsilon_{\rm err} \propto N_{\rm s}^{-1/2}$ as expected in Monte Carlo simulations.
• Figure 3: Panel (a): Scattering test for an uniformly spherical cloud of gas. The GPUmonty spectrum is shown with the red solid line and igrmonty's with the blue dashed line, both with $N_{\rm s} = 10^8$.Panel (b): The frequency-dependent residuals, computed as in Figure \ref{['fig:thin_synch_test']}, for different values of $N_{\rm s}$. For clarity, the curves have been smoothed using a moving average over the $10$ nearest frequency bins.
• Figure 4: Normalized integrated error (equation \ref{['eq:conv_parameter']}) for different $N_{\rm s}$ sizes. The green dashed line represents $\epsilon_{\rm err} \propto N_{\rm s}^{-1/2}$.
• Figure 5: Panel (a): Comparison of spectra computed with GPUmonty (red solid line) and igrmonty (blue dashed line) for a 3D SANE accretion flow GRMHD simulation around a Kerr black hole. Both spectra use the same parameters and $N_{\rm s} = 10^6$. Panel (b): Residuals computed as in Figure \ref{['fig:thin_synch_test']} with igrmonty as baseline.