  Understanding the interaction of vortices with inclusions in type-II
superconductors is a major outstanding challenge both for fundamental science
and energy applications. At application-relevant scales, the long-range
interactions between a dense configuration of vortices and the dependence of
their behavior on external parameters, such as temperature and an applied
magnetic field, are all important to the net response of the superconductor.
Capturing these features, in general, precludes analytical description of
vortex dynamics and has also made numerical simulation prohibitively expensive.
Here we report on a highly optimized iterative implicit solver for the
time-dependent Ginzburg-Landau equations suitable for investigations of type-II
superconductors on massively parallel architectures. Its main purpose is to
study vortex dynamics in disordered or geometrically confined mesoscopic
systems. In this work, we present the discretization and time integration
scheme in detail for two types of boundary conditions. We describe the
necessary conditions for a stable and physically accurate integration of the
equations of motion. Using an inclusion pattern generator, we can simulate
complex pinning landscapes and the effect of geometric confinement. We show
that our algorithm, implemented on a GPU, can provide static and dynamic
solutions of the Ginzburg-Landau equations for mesoscopically large systems
over thousands of time steps in a matter of hours. Using our formulation,
studying scientifically-relevant problems is a computationally reasonable task.
