Immersed boundary methods (IBMs) make the simulation of flow around complex, moving, and deforming bodies possible on the relatively simple Cartesian grids. The sharp-interface IBM even provides high-order accuracy and conservation that are close to ones with body-fitted grid methods. Although the sharp-interface IBM is highly scalable for the parallel computation and applicable to the large-scale problems for realistic, engineering configurations, it still requires huge amount of computational resources and time. In this work, therefore, we improve the performance of our in-house, sharp-interface IBM flow solver, ViCar3D by using graphical processing units (GPUs). The implementations are done by using the OpenACC, CUDA, and CUDA-aware MPI to port the current solver to multi-GPU architectures. Verification and scalability studies are performed for various benchmark cases including a direct numerical simulation (DNS) of flow past a finite span rectangular wing. We observe an approximately 20 times speedup (node-to-node comparison) relative to the CPU-based computation. The GPU powered solver is then applied to various complex moving body problems in the application to aero/hydro-dynamics and bio-inspired engineering. The GPU powered solver can simulate complex 3D flows with up to 200 million mesh points on a single node equipped with four GPUs, and strong as well as weak scaling tests demonstrate maximum scaling efficiencies of 92% and 93%, respectively, on multi-GPU systems relative to the single GPU system.