Purpose: To enable accelerated Bloch simulations by enhancing the open-source multi-purpose MRI simulation tool JEMRIS with graphic processing units (GPU) parallelization. Methods: A GPU-compatible version of JEMRIS was built by shifting the computationally expensive parallelizable processes to the GPU to benefit from heterogeneous computing and by adding asynchronous communication and mixed precision support. With key classes reimplemented in CUDA C++, the developed GPU-JEMRIS framework was tested on simulations of common MRI artifacts in numerical phantoms. The accuracy and performance of the GPU-parallelized JEMRIS simulator were benchmarked against the CPU-parallelized JEMRIS and GPU-enabled KomaMRI.jl simulators. Additionally, an example of liver fat quantification errors due to respiratory motion artifacts was simulated in a multi-echo gradient echo (MEGRE) acquisition. Results: The GPU-accelerated JEMRIS achieved speed-up factors 3–12 and 7–65 using double and single precision numerical integrators, respectively, when compared to the parallelized CPU implementation in the investigated numerical phantom scenarios. While double precision GPU simulations negligibly differ (<0.1% NRMSE) from double precision CPU simulations, the single precision simulations still present small errors of up to 1% k-space signal NRMSE. The developed a GPU extension enabled computationally demanding motion simulations with a multi-species abdominal phantom and a MEGRE sequence, showing significant and spatially varying fat fraction bias in the presence of motion. Conclusion: By solving the Bloch equations in parallel on device, accelerated Bloch simulations can be performed on any GPU-equipped device with CUDA support using the developed GPU-JEMRIS. This would enable further insights into more realistic large spin pool MR simulations such as experiments with large multi-dimensional phantoms, multiple chemical species and dynamic effects.