Magnetohydrodynamic (MHD) simulations in the domain of spherical shell are a crucial and challenging subject in many fields such as geophysics and solar-terrestrial physics, due to the complication of the MHD equations and the specificity of the domain. Besides, due to the real-time requirement, accelerating the heavy computation is proposed in many practical problems, of which the space weather simulation and forecast from the Sun to Earth is a typical case. Considering these factors, we first develop a new, spatially second-order accurate finite volume (FV) solver for three-dimensional (3D) MHD simulations with the multiple time steps strategy, which is based on the six-component grid for spherical shell domain. Then to speed up the simulation, we implement the solver on multiple GPUs with optimizations of CUDA and establish an effective multi-GPU FV solver on the spherical shell domain. A MHD manufactured solution is used to validate the solvers’ spatial accuracy, and to measure their performances. Results show that both solvers have nice scalability, and speedup ratios of 27.7x to 30.06x are obtained on GPUs. Then we utilize them to study the ambient solar wind for Carrington rotation (CR) 2060. The multi-GPU FV solver can not only obtain speedup ratios of about 29.0x, but capture main features of the solar corona and the mapped in-situ solar wind measurements. © 2018