Numerical methods

Accurate and efficient numerical methods are of paramount importance in theoretical geophysics and fluid dynamics in general. Several reasons favor spectral methods over local grid methods for the problems I am interested in. The dynamo problem is a stability problem, so that numerical diffusivities need to be kept at bay as much as possible. Numerical diffusivity in any reasonable spectral discretization is virtually zero. In the dynamo problem, one also has to enforce a non-local boundary condition, since the magnetic field in the computational domain needs to be matched to an external potential field in the vacuum which, if possible, is not computed numerically and only appears implicitly in the boundary conditions. Matching to the vacuum field is especially convenient at a spherical interface. After decomposition into spherical harmonics, the non-local boundary condition decouples into separate boundary conditions for every spectral coefficient.

The method I used in most problems decomposes solenoidal fields into poloidal and toroidal fields, which in turn are decomposed into Chebychev polynomials in radius and into spherical harmonics. This decomposition is appropriate for spherical and also ellipsoidal domains. The time step treats diffusion implicitly and all other terms explictly. Different combinations of implicit Euler and Crank-Nicolson steps for the implicit part and Adams-Bashforth schemes of second or third order for the explicit part have been useful in different circumstances.

The code can quite easily be changed so that the radial discretization is made with a finite difference method. This simplifies the code and accelerates it a bit. Most importantly, communication problems are vastly reduced on parallel machines. However, accuracy is lost at an equal number of grid points. The figure shows the relative error on the computation of a poloidal magnetic decay mode as a function of the number Nr of gridpoints used. The circles are for centered finite differences and the error decreases as 1/Nr2. The squares are for a discretization with Chebychev polynomials. In the latter case, the error decreases much more rapidly with Nr (in theory faster than any power of Nr).

Full details of the methods are given in the references.