As a counterpoint to recent numerical methods for crystal surface evolution, which agree well with microscopic dynamics but suffer from significant stiffness that prevents simulation on fine spatial grids, we develop a new numerical method based on the macroscopic partial differential equation, leveraging its formal structure as the gradient flow of the total variation energy, with respect to a weighted $H^{-1}$ norm. This gradient flow structure relates to several metric space gradient flows of recent interest, including 2-Wasserstein flows and their generalizations to nonlinear mobilities. We develop a novel semi-implicit time discretization of the gradient flow, inspired by the classical minimizing movements scheme (known as the JKO scheme in the 2-Wasserstein case). We then use a primal dual hybrid gradient (PDHG) method to compute each element of the semi-implicit scheme. In one dimension, we prove convergence of the PDHG method to the semi-implicit scheme, under general integrability assumptions on the mobility and its reciprocal. Finally, by taking finite difference approximations of our PDHG method, we arrive at a fully discrete numerical algorithm, with iterations that converge at a rate independent of the spatial discretization: in particular, the convergence properties do not deteriorate as we refine our spatial grid. We close with several numerical examples illustrating the properties of our method, including facet formation at local maxima, pinning at local minima, and convergence as the spatial and temporal discretizations are refined.

24 pages, 6 figures, comments welcome!