  The study of fractional order differential operators is receiving renewed
attention in many scientific fields. In order to accommodate researchers doing
work in these areas, there is a need for highly scalable numerical methods for
solving partial differential equations that involve fractional order operators
on complex geometries. These operators have desirable special properties that
also change the computational considerations in such a way that undermines
traditional methods and makes certain other approaches more appealing. We have
developed a scalable framework for solving fractional diffusion equations using
one such method, specifically the method of eigenfunction expansion. In this
paper, we will discuss the specific parallelization strategies used to
efficiently compute the full set of eigenvalues and eigenvectors for a
discretized Laplace eigenvalue problem and apply them to construct approximate
solutions to our fractional order model problems. Additionally, we demonstrate
the performance of the method on the Frontera computing cluster and the
accuracy of the method on simple geometries using known exact solutions.
