A scalable, parallel matrix-free Stokes solver for geodynamic applications

Here I describe a numerical method suitable for studying non-linear, large deformation processes in crustal and lithopspheric dynamics. The method utilizes a hybrid spatial discretisation which consists of mixed finite elements for the Stokes flow problem, coupled to a Lagrangian marker based discretisation to represent the material properties (viscosity and density). This approach is akin to the classical Marker-And-Cell (MAC) scheme of Harlow and the subsequently developed Material Point Method (MPM) of Sulsky and co-workers. The geometric flexibility and ease of modelling large deformation processes afforded by such mesh-particle methods has been exploited by the lithospheric dynamics community over the last 20 years. The strength of the Stokes preconditioner fundamentally controls the scientific throughput achievable and represents the largest bottleneck in the development of our understanding of geodynamic processes. The possibility to develop a 'cheap' and efficient preconditioning methodology which is suitable for the mixed Q2-P1 element is explored here. I describe a flexible strategy, which aims to address the Stokes preconditioning issue using an upper block triangular preconditioner, together with a geometric multi-grid preconditioner for the viscous block. The key to the approach is to utilize algorithms and data-structures that exploit current multi-core hardware and avoid the need for excessive global reductions. In order to develop a scalable method, special consideration is given to; the definition of the coarse grid operator, the smoother and the coarse grid solver. The performance characteristics of this hybrid matrix-free / partially assembled multi-level preconditioning strategy is examined. The robustness of the preconditioner with respect to the viscosity contrast and the topology of the viscosity field, together with the parallel scalability is demonstrated.