Relativistic radiative gas dynamics models are widely implemented in simulations of radiation-driven astrophysical processes. We investigate a gamma-ray burst (GRB) mathematical model to interpret observed radiation spectra, considering a three-dimensional relativistic radiative gas dynamics model with quasi-stationary radiative transfer equations and multiple scattering kernels. The numerical scheme combines physical process splitting with unstructured tetrahedral meshes, using a Godunov-type HLLC method for the hyperbolic gas dynamics and the short characteristics method for radiation field approximation. The parallel algorithm is designed for GPU-equipped multiprocessor clusters, with collision integral calculations distributed across GPU nodes. We analyze the algorithm’s convergence and scalability, simulating both relativistic jet evolution in interstellar media during GRBs and photon spectrum redistribution via Compton scattering on relativistic steady-state flows.