The In-Core Fuel Management or Loading Pattern Optimization (LPO) is the problem of finding an optimal configuration of nuclear fuel assemblies in a reactor core. Over the years Optimization Metaheuristics based on populations have been used successfully for solving the LPO. However, for such methods, thousands of evaluations of candidate loading patterns with reactor physics codes are required. In addition, more precise calculations are also desirable, which in principle may conflict with faster evaluations. For coping with both criteria simultaneously, a code based on the Galerkin Finite Element Method (GFEM) for solving the neutron diffusion equation designed for Graphics Processing Unit (GPU) is proposed in the present work. The resulting program GFEM_GPU has been implemented in Python language with a library for GPU-accelerated computing. A comparison between GFEM_GPU and the CPU version of GFEM regarding three nuclear reactor benchmarks (IAEA-2D, BIBLIS-2D, and ZION-2D) is provided and the accuracy of the results is compared to the Nodal Expansion Method (NEM). The implementation of the GFEM solve-phase with GFEM_GPU provides iterative solutions up to 50% faster than the CPU version depending on the discretization of the problem. This approach presents several advantages, such as a high-level implementation, faster prototyping, and application to LPO problems, while keeping the high accuracy provided by fine mesh methods.