  Black holes are believed to be one of the key ingredients of galaxy formation
models, but it has been notoriously challenging to simulate them due to the
very complex physics and large dynamical range of spatial scales involved. Here
we address a significant shortcoming of a Bondi-Hoyle-like prescription
commonly invoked to estimate black hole accretion in cosmological hydrodynamic
simulations of galaxy formation, namely that the Bondi-Hoyle radius is
frequently unresolved. We describe and implement a novel super-Lagrangian
refinement scheme to increase, adaptively and 'on the fly', the mass and
spatial resolution in targeted regions around the accreting black holes at
limited computational cost. While our refinement scheme is generically
applicable and flexible, for the purpose of this paper we select the smallest
resolvable scales to match black holes' instantaneous Bondi radii, thus
effectively resolving Bondi-Hoyle-like accretion in full galaxy formation
simulations. This permits us to not only estimate gas properties close to the
Bondi radius much more accurately, but also allows us to improve black hole
accretion and feedback implementations. We thus devise a more generic feedback
model where accretion and feedback depend on the geometry of the local gas
distribution and where mass, energy and momentum loading are followed
simultaneously. We present a series of tests of our refinement and feedback
methods and apply them to models of isolated disc galaxies. Our simulations
demonstrate that resolving gas properties in the vicinity of black holes is
necessary to follow black hole accretion and feedback with a higher level of
realism and that doing so allows us to incorporate important physical processes
so far neglected in cosmological simulations.
