  We present radiation hydrodynamics simulations of the collapse of massive
pre-stellar cores. We treat frequency dependent radiative feedback from stellar
evolution and accretion luminosity at a numerical resolution down to 1.27 AU.
In the 2D approximation of axially symmetric simulations, it is possible for
the first time to simulate the whole accretion phase of several 10^5 yr for the
forming massive star and to perform a comprehensive scan of the parameter
space. Our simulation series show evidently the necessity to incorporate the
dust sublimation front to preserve the high shielding property of massive
accretion disks. Our disk accretion models show a persistent high anisotropy of
the corresponding thermal radiation field, yielding to the growth of the
highest-mass stars ever formed in multi-dimensional radiation hydrodynamics
simulations. Non-axially symmetric effects are not necessary to sustain
accretion. The radiation pressure launches a stable bipolar outflow, which
grows in angle with time as presumed from observations. For an initial mass of
the pre-stellar host core of 60, 120, 240, and 480 Msol the masses of the final
stars formed in our simulations add up to 28.2, 56.5, 92.6, and at least 137.2
Msol respectively.
