Recrystallization models and simulations have been the subject of much attention in materials community in the past decades due to this process having a significant effect on many technologically important materials characteristics. Statistical analysis performed close to the steady state requires large-scale simulations, which are often prohibitively expensive from computational point of view. Graphical Processing Unit (GPU)-based realizations provide a viable approach to addressing this challenge, yet they remain relatively under-explored in this context. In the present manuscript, we develop a fully-parallelizable matrix-free GPU-based algorithm for implementing a two-dimensional vertex model of recrystallization based on the stored energy formalism. Nucleation is assumed to take place at triple junctions and obeys a Metropolis-type criterion. We include a complete mathematical analysis of the nucleation model deriving conditions under which nucleation is successful. Stability analysis of the dynamics of a triple junction under the presence of bulk energy is provided. On the computational side, we propose a novel polling system for handling topological transitions to ensure robust GPU implementation. Single grain tests are performed for benchmarking purposes. Finally, a set of numerical experiments for large scale systems is presented to explore the effect of initial distributions of stored energy on several statistical characteristics.