Code benchmarking: why should we spend time on it?

Benchmarking of a numerical code means comparing the numerical solution obtained with solving the system of linear equations with (i) analytical solutions (ii) results of physical (analogue) experiments (iii) numerical results from other (well-established) codes and (iv) general physical considerations. Benchmarking of newly created numerical tools is sometimes very tedious, but an absolutely necessary stage of code development as its purpose is to test the code's robustness in a broad range of situations relevant to geodynamic modelling applications. For instance, if you plan to model with your code shear heating processes in deforming rocks – make sure that your code provides the correct temperature changes related to mechanical energy dissipation; if you model subduction – make sure that your code handles correctly large viscosity contrasts and has no notable numerical diffusion of the temperature field and composition; if you intend to model self-gravitating planetary bodies – make sure that your code computes correct gravity field etc. We should not be lazy and limit ourselves to one or two common benchmarks, such as the Rayleigh–Taylor instability and convection with constant viscosity hoping that everything else will work automatically. No, it will not! Therefore, test your code on a broad range of challenging cases (several of which are discussed below), to explore its limitations and be creative in inventing and calibrating new numerical benchmarks.

Multigrid – what is it?

The use of direct methods to solve the system of equations places a strong limitation on the maximum possible number of nodal points within a numerical grid due to limitations in computer memory and computational speed. This limitation is particularly critical in 3D where, to reach the same resolution as in 2D (tens and hundreds of grid points in each direction), the amount of linear equations to be solved increases by at least two orders of magnitude and the number of computational operations required to solve the same equations grows by at least three orders of magnitude. Try to increase resolution in your visco-elasto-plastic code (Exercise 13.2) by two to three times in each direction and you will see how much slower it will compute …

Therefore, for very high resolution in 2D and for a moderate resolution in 3D, we are forced to use iterative methods which do not have such strong memory limitations. However, many iterative methods also have the problem that they require an increasing number of iterations with an increasing number of grid points (Fig. 14.1).

Grids

As we already discussed, the numerical solution of partial differential equations (PDEs) requires the definition of a grid of nodal points within the numerical model. The choice of this grid depends strongly on the type of equations to be solved. Discretisation schemes for these equations will also change with the changing types of numerical grids. The following types of numerical grid exist in numerical geodynamic modelling:

• Depending on the dimension of the problem, the numerical grid can be one, two and three dimensional (1D, 2D, 3D) (Fig. 7.1).

• Depending on the shape of the basic elements, the grid can be rectangular and triangular (Fig. 7.2).

• Depending on the distribution of nodal points, the grid can be regular and non-regular (irregular) (Fig. 7.3).

• Depending on the distribution of different variables within the grid, it can be non-staggered or staggered (Fig. 7.4).

The simplest grids are non-staggered. All variables are defined at the same nodal points (Fig. 7.5). When using finite differences (FD) with such a grid, all equations are formulated at the same nodal points.

Introduction

Program 1: Visualisation_is_important.m (Exercise Introduction.2) – visualisation of ‘sin’ and ‘cos’ functions with ‘plot’, ‘pcolor’, ‘contour’ and ‘surf’.

Chapter 1

Program 2: Divergence.m (Exercise 1.2) – computation and visualisation of velocity, divergence of velocity and time derivatives of density with ‘pcolor’ and ‘quiver’.

Chapter 2

Program 3: Periclase_EOS.m (Exercise 2.2) – computation and visualisation of density, thermal expansion and compressibility for periclase (MgO) using external Gibbs free energy function G_periclase.m.

Program 4: Density_map.m (Exercise 2.3) – loading from data files (m895_ro, morn_ro) and visualising density maps for pyrolite (m895_ro) and MORB (morn_ro) and density difference between pyrolite and MORB.

Chapter 3

Program 5: Poisson1D.m (Exercise 3.1) – solution of 1D Poisson equation with finite differences on a regular grid using direct solver ‘\’.

Program 6: Poisson2D_direct.m (Exercise 3.2) – solution of 2D Poisson equation with finite differences on a regular grid using direct solver ‘\’.

Program 7: Poisson2D_Gauss_Seidel.m (Exercise 3.3) – solution of 2D Poisson equation with finite differences on a regular grid using Gauss–Seidel iteration.

Program 8: Poisson2D_Jacobi.m (Exercise 3.4) – solution of 2D Poisson equation with finite differences on a regular grid using Jacobi iteration.

Chapter 4

Program 9: Strain_rate.m (Exercise 4.2) – computation and visualisation of velocity field, strain rate, deviatoric strain rate, and second strain rate invariant.

What is this book?

This book is a practical, hands-on introduction to numerical geodynamic modelling for inexperienced people, i.e. for young students and newcomers from other fields. It does not require much background in mathematics or physics and is therefore written with a maximum amount of simple technical details. If you are inexperienced – this book is for you!

What this book is not

This book is not a treatise or a compendium of knowledge for experienced researchers. It does not contain large overviews of existing numerical techniques and only simple approaches are explained. If you are experienced in numerical methods, read Chapter 17 first and then decide if you wish to read about the technical details presented in previous chapters.

Get started

Already decided?! Then let us get started! In recent decades numerical modelling has become an essential approach in geosciences in general and in geodynamics in particular. This is a very natural process (‘instinctive evolution’) since direct human observation scales are extremely limited in both time and space (depth) and rapid progress in computer technology offers every day new and exceptional possibilities to explore sophisticated mathematical models and this is true in every discipline, and even industrial applications.