  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.
