  We present the multi-GPU realization of the StePS (Stereographically
Projected Cosmological Simulations) algorithm with MPI-OpenMP-CUDA hybrid
parallelization and nearly ideal scale-out to multiple compute nodes. Our new
zoom-in cosmological direct N-body simulation method simulates the infinite
universe with unprecedented dynamic range for a given amount of memory and, in
contrast to traditional periodic simulations, its fundamental geometry and
topology match observations. By using a spherical geometry instead of periodic
boundary conditions, and gradually decreasing the mass resolution with radius,
our code is capable of running simulations with a few gigaparsecs in diameter
and with a mass resolution of $\sim 10^{9}M_{\odot}$ in the center in four days
on three compute nodes with four GTX 1080Ti GPUs in each. The code can also be
used to run extremely fast simulations with reasonable resolution for fitting
cosmological parameters. These simulations are useful for prediction needs of
large surveys. The StePS code is publicly available for the research community.
