GPU-based pseudo-transient finite difference solution for power-law viscous flow in cartesian, polar and spherical coordinates

Power-law viscous flow describes well 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.

Stress calculations in a geometrically and mechanically heterogeneous 3-D lithospheric shell require high-resolution and high-performance computing. The pseudo-transient finite difference (PTFD) method recently enabled efficient simulations of high-resolution 3-D deformation processes, implementing an iterative implicit solution strategy of the governing equations for power-law viscous flow. Main challenges for the PTFD method is to guarantee convergence, minimize the required iteration count and speed-up the iterations.

Here, we present PTFD simulations for simple mechanically heterogeneous (weak circular inclusion) incompressible 2-D power-law viscous flow in cartesian and cylindrical coordinates. The flow laws employ a pseudo-viscoelastic behavior to optimize the iterative solution by exploiting the fundamental characteristics of viscoelastic wave propagation.

The developed PTFD algorithm executes in parallel on CPUs and GPUs. The development was done in Matlab (mathworks.com), then translated into the Julia language (julialang.org), and finally made compatible for parallel GPU architectures using the ParallelStencil.jl package (https://github.com/omlins/ParallelStencil.jl). We may unveil preliminary results for 3-D spherical configurations including gravity-controlled lithospheric stress distributions around continental plateaus.