This article considers the iterative solution of a finite element discretisation of the magma dynamics equations. In simplified form, the magma dynamics equations share some features of the Stokes equations. We therefore formulate, analyse and numerically test a Elman, Silvester and Wathen-type block preconditioner for magma dynamics. We prove analytically and demonstrate numerically the optimality of the preconditioner. The presented analysis highlights the dependence of the preconditioner on parameters in the magma dynamics equations that can affect convergence of iterative linear solvers. The analysis is verified through a range of two- and three-dimensional numerical examples on unstructured grids, from simple illustrative problems through to large problems on subduction zone-like geometries. The computer code to reproduce all numerical examples is freely available as supporting material.