I just committed a new, accurate integration method "Driscoll&Healy" in thorn EinsteinAnalysis/Multipole. This method is described in J.R. Driscoll and D.M. Healy Jr., Computing Fourier transforms and convolutions on the 2-sphere, Advances in Applied Mathematics, 15(2):202–250, 1994. This method uses the interpolated values on the sphere to define a spectral reconstruction of the function on the sphere, and integrates this function, yielding essentially exponential convergence.
Of course, accuracy in a production simulation will still be limited by the accuracy of the interpolation onto the sphere, as well as the accuracy of the time evolution itself.
-erik