  We present high-resolution, three-dimensional hydrodynamic simulations of the
fueling of supermassive black holes in elliptical galaxies from a turbulent
medium on galactic scales, taking M87* as a typical case. The simulations use a
new GPU-accelerated version of the Athena++ AMR code, and span more than 6
orders of magnitude in radius, reaching scales similar to the black hole
horizon. The key physical ingredients are radiative cooling and a
phenomenological heating model. We find that the accretion flow takes the form
of multiphase gas at radii less than about a kpc. The cold gas accretion
includes two dynamically distinct stages: the typical disk stage in which the
cold gas resides in a rotationally supported disk and relatively rare chaotic
stages ($\lesssim 10\%$ of the time) in which the cold gas inflows via chaotic
streams. Though cold gas accretion dominates the time-averaged accretion rate
at intermediate radii, accretion at the smallest radii is dominated by hot
virialized gas at most times. The accretion rate scales with radius as
$\dot{M}\propto r^{1/2}$ when hot gas dominates and we obtain
$\dot{M}\simeq10^\mathrm{-4}-10^\mathrm{-3}\,M_\odot\,\mathrm{yr^{-1}}$ near
the event horizon, similar to what is inferred from EHT observations. The
orientation of the cold gas disk can differ significantly on different spatial
scales. We propose a subgrid model for accretion in lower-resolution
simulations in which the hot gas accretion rate is suppressed relative to the
Bondi rate by $\sim (r_\mathrm{g}/r_{\rm Bondi})^{1/2}$. Our results can also
provide more realistic initial conditions for simulations of black hole
accretion at the event horizon scale.
