  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 (up to the end of the
accretion disk epoch) for the forming massive star and to perform a broad 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. While confirming the upper mass limit of
spherically symmetric accretion, our disk accretion models show a persistent
high anisotropy of the corresponding thermal radiation field. This yields to
the growth of the highest-mass stars ever formed in multi-dimensional radiation
hydrodynamics simulations, far beyond the upper mass limit of spherical
accretion. 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 Msun the masses of the final
stars formed in our simulations add up to 28.2, 56.5, 92.6, and at least 137.2
Msun respectively.
