Black hole accretion disks serve as the central engines driving high-energy phenomena such as X-ray binaries, active galactic nuclei, and gamma-ray bursts, and give rise to complex physical processes such as turbulence induced by magnetorotational instability, intense heating, and the formation of relativistic jets. Traditional one-dimensional theoretical models and general relativistic magnetohydrodynamic (GR-MHD) simulations have successfully captured many aspects of these phenomena, but they fall short in reproducing the effects of radiative cooling and radiation pressure, which are critical in high accretion rate environments. To overcome these limitations, general relativistic radiation magnetohydrodynamics (GR-RMHD) simulations have been developed to incorporate the interplay between radiation and plasma dynamics. However, the inclusion of radiative transfer significantly increases the computational complexity, requiring advanced computational strategies for high-fidelity multidimensional simulations. In this study, we present a GR-RMHD code implemented using CUDA Fortran and MPI that achieves up to 39-fold speedup over conventional CPU implementations, with nearly linear scalability as the number of GPUs increases. This acceleration enables detailed analysis of the complex dynamics of black hole accretion disks and jet formation, providing new insights into the energy production mechanisms in high-energy astrophysical systems. We discuss the evolution of theoretical models and numerical simulations from GR-MHD to GR-RMHD, describe the algorithms and GPU implementation of our approach, and provide quantitative performance evaluations.
KOBAYASHI et al. (Thu,) studied this question.