We develop a fully discrete approximation scheme to fast simulate centered Gaussian random fields
(GRFs) over compact Riemannian manifolds. These GRFs can be either anisotropic or isotropic ones.
The scheme is proposed by the truncation of a series expansion depending solely on the decay of
the angular power spectrum, which is independent of the chosen space and time discretization. The
$p$th moment convergence for the scheme is proved and the error bound with its associated rate of
convergence is determined. Three simulation case studies are conducted by
applying the Laplacian operator to find the required surface harmonics for a torus endowed with
a Riemanian metric and by applying the ellipsoidal Beltrami and Laplacian operators to find the
required surface ellipsoidal harmonics respectively for anisotropic and isotropic GRFs over
ellipsoids endowed with two different Riemannian metrics.