Summary: We present a simulation code which can solve a broad range of partial differential equations in a full sphere. The code expands tensorial variables in a spectral series of spin-weighted spherical harmonics in the angular directions and a scaled Jacobi polynomial basis in the radial direction, as described in Vasil et al. (2018; hereafter, Part-I). Nonlinear terms are calculated by transforming from the coefficients of the spectral series to the value of each quantity on the physical grid, where it is easy to calculate products and perform other local operations. The expansion makes it straightforward to solve equations in tensor form (i.e., without decomposition into scalars). We propose and study several unit tests which demonstrate the code can accurately solve linear problems, implement boundary conditions, and transform between spectral and physical space. We then run a series of benchmark problems proposed in [P. Marti et al., Geophys. J. Int. 197, No. 1, 119–134 (2014; doi:10.1093/gji/ggt518)], implementing the hydrodynamic and magnetohydrodynamic equations. We are able to calculate more accurate solutions than reported in [Marti et al., loc. cit.] by running at higher spatial resolution and using a higher-order timestepping scheme. We find the rotating convection and convective dynamo benchmark problems depend sensitively on details of timestepping and data analysis. We also demonstrate that in low resolution simulations of the dynamo problem, small changes in a numerical scheme can lead to large changes in the solution. To aid future comparison to these benchmarks, we include the source code used to generate the data, as well as the data and analysis scripts used to generate the figures.
For Part I, see [G. M. Vasil et al., J. Comput. Phys.: X 3, Article ID 100013, 31 p. (2019; Zbl 07785508)].