Numerical simulations were done to find low energy states in frustrated large square Josephson Junction arrays in a perpendicular magnetic field using simulated annealing on the coupled RSJ model. These simulations were made possible by a new algorithm suitable for parallel gpu computing and reduced complexity. Free boundary conditions were used so that values of the frustration factor f that are incommensurate with the array size are permitted. The resulting energy as a function of f is continuous with logarithmic discontinuities in the derivative dE/df at rational frustration factors f=p/q with small q, substantiating the mathematical proof that this curve is continuous and further showing that the staircase state hypothesis is incorrect. The solution shows qualitative similarities with the lowest energy branch of the Hofstadter butterfly, which is a closely related problem. Furthermore, it is found that at the edge of an array there are either extra vortices or missing vortices depending the frustration factor, and the width of this region is independent of the array size.