Deep convection occurs in the outer one-third of the solar interior and transports energy generated by nuclear reactions to the surface. It leads to a characteristic pattern of time-averaged differential rotation, with the poles rotating approximately 20% slower than the equator. A particularly notable feature of the solar differential rotation is that it departs significantly from the Taylor-Proudman state of rotation constant on cylinders aligned with the rotation axis. Although this observation contrasts with results from early numerical simulations, such simulations provide the best hope of understanding the observations. Many studies have adopted the DNS (Direct Numerical Simulation) approach and justified the artificially large viscosities and thermal diffusivities used as modelling transport by unresolved eddies. LES (Large Eddy Simulation) techniques, which use a suitable turbulence closure model, offer a superior alternative but face the problem of choosing an appropriate turbulence closure; this can be difficult in the face of complicating factors such as stratification and rotation. An alternative approach is to shift responsibility for truncating the turbulent cascade to the numerical scheme itself. Since this approach abandons the rigorous notions of the LES approach, we refer to it as a VLES (Very Large Eddy Simulation). This paper compares results of DNS simulations carried out with a spherical harmonic code, and preliminary results obtained using a VLES-type code. Both make the anelastic approximation.