We present ExaPF.jl, a solver for power flow on GPUs, entirely written in Julia. It implements a highly parallel Newton-Raphson solver for nonlinear equations. We exploit Julia packages for kernel and array abstractions at the modeling level, and generate efficient codes at runtime for both CPU and NVIDIA GPU using Julia's inherent metaprogramming capabilities. In the future, this infrastructure will allow us to leverage this machinery for AMD and Intel GPUs by targeting the ROCm and oneAPI frameworks. The composable design of the Julia language allows us to apply automatic differentiation that alleviates the user from providing derivatives. We also detail a GPU implementation of the iterative solver BiCGSTAB to solve efficiently the linear systems arising in the Newton-Raphson algorithm, and show how to improve its performance with a block-Jacobi preconditioner tailored to the batched matrix inversion capabilities of modern GPUs. The Newton-Raphson algorithm will eventually serve as the foundation of a reduced-space optimization method that will run entirely on the GPU.