High Performance Computing (HPC) is a powerful tool in scientific research for experimental simulation and data analysis. However, it has not been well studied and applied for earthquake cycle simulation. In this work, we explore computational methods and HPC techniques for earthquake cycle simulations with a focus on iterative methods and preconditioning techniques. The problem that we target is solving Poisson's equation numerically, with boundary conditions derived from geophysics. Specifically, we use summation-by-parts (SBP) finite difference operators with the simultaneous-approximation-term (SAT) method. The combined SBP-SAT scheme has high-order accuracy and provable stability. The purpose of this work is to examine existing research on problems similar to ours, so we can design algorithms for large-scale 3D earthquake cycle simulations suitable for HPC platforms in the future. We begin with a mathematical introduction to our problem, and then provide an overview of iterative methods for solving the problem numerically. Finally, we explore the HPC aspects of these different iterative methods.