We perform general relativistic radiation magnetohydrodynamics (MHD) simulations of super-Eddington accretion flows around a neutron star with a dipole magnetic field for modeling the galactic ultra-luminous X-ray source (ULX) exhibiting X-ray pulsations, Swift J0243.6+6124. Our simulations show the accretion columns near the magnetic poles, the accretion disk outside the magnetosphere, and the outflows from the disk. It is revealed that the effectively optically thick outflows, consistent with the observed thermal emission at \(\sim10^7\) K, are generated if the mass accretion rate is much higher than the Eddington rate \(\dot{M}_{\rm Edd}\) and the magnetospheric radius is smaller than the spherization radius. In order to explain the blackbody radius (\(\sim 100-500\) km) without contradicting the reported spin period (\(9.8~{\rm s}\)) and spin-up rate (\(\dot{P}=-2.22\times10^{-8}~{\rm s~s^{-1}}\)), the mass accretion rate of \((200-1200)\dot{M}_{\rm Edd}\) is required. Since the thermal emission was detected in two observations with \(\dot{P}\) of \(-2.22\times10^{-8}~{\rm s~s^{-1}}\) and \(-1.75\times10^{-8}~{\rm s~s^{-1}}\) but not in another with \(\dot{P}=-6.8 \times10^{-9}~{\rm s~s^{-1}}\), the surface magnetic field strength of the neutron star in Swift J0243.6+6124 is estimated to be between \(3\times10^{11}~{\rm G}\) and \(4\times10^{12}~{\rm G}\). From this restricted range of magnetic field strength, the accretion rate would be \((200-500)\dot{M}_{\rm Edd}\) when the thermal emission appears and \((60-100)\dot{M}_{\rm Edd}\) when it is not detected. Our results support the hypothesis that the super-Eddington phase in the 2017-2018 giant outburst of Swift J0243.6+6124 is powered by highly super-Eddington accretion flows onto a magnetized neutron star.