GPU-based pseudo-transient finite difference solution for 3-D gravity- and shear-driven power-law viscous flow

Power-law viscous flow describes the first-order features of long-term lithosphere deformation. Due to the ellipticity of the Earth, the lithosphere is mechanically analogous to a shell, characterized by a double curvature. The mechanical characteristics of a shell are fundamentally different to the characteristics of plates, having no curvature in their undeformed state. The systematic quantification of the magnitude and the spatiotemporal distribution of strain, strain-rate and stress inside a deforming lithospheric shell is thus of major importance: stress is, for example, a key physical quantity that controls geodynamic processes such as metamorphic reactions, decompression melting, lithospheric flexure, subduction initiation or earthquakes. Calculating these stresses in a three-dimensional (3-D), geometrically and mechanically heterogeneous lithosphere requires high-resolution and high-performance computing.

 

Here, we present numerical simulations of 3-D power-law viscous flow. We employ the pseudo-transient finite difference (PTFD) method, which enables efficient simulations of high-resolution 3-D deformation processes by implementing an iterative implicit solution strategy of the governing equations. The main challenges for the PTFD method are to guarantee convergence, minimize the required iteration count and speed-up the iterations. We implemented the PTFD algorithm using the Julia language (julialang.org) to enable optimal parallel execution on multiple CPUs and GPUs using the ParallelStencil.jl module (https://github.com/omlins/ParallelStencil.jl). ParallelStencil.jl enables execution on multi-threaded CPUs and Nvidia GPUs using a single switch.

 

We present PTFD simulations of mechanically heterogeneous (weak and less dense spherical inclusion), incompressible 3-D power-law viscous flow under gravity in cartesian, cylindrical and spherical coordinates systems. The viscous flow is described by a linear combination of a linear viscous and a power-law viscous flow law, representing diffusion and dislocation creep, respectively. The iterative solution strategy builds upon pseudo-viscoelastic behavior to minimize the iteration count by exploiting the fundamental characteristics of viscoelastic wave propagation. We performed systematic numerical simulations to investigate the impact of (i) buoyancy versus shear forces and (ii) linear versus power-law viscous flow on the vertical velocity of the spherical inclusion under bulk strike-slip shearing. We report the systematic results using the controlling dimensionless numbers and compare the numerical results with analytical predictions for buoyancy-driven flow of inclusions in a power-law matrix. We also aim to unveil preliminary results for a vertically and locally loaded power-law viscous lithosphere showing the impact of different lithosphere curvatures on the resulting stress field.