Advances in 3D geodynamic modeling
Anton Popov
This projects aims at bringing the 3D lithospheric deformation modeling to a qualitatively different level. Since this a very challenging task we have to take a number of correct decision steps regarding the discretization and solution algorithms. Our new code LaMEM (Lithosphere and Mantle Evolution Model) is based on the following solid building blocks:
* Massively-parallel data-distributed implementation model based on PETSc library
* Light, stable and accurate staggered-grid finite difference spatial discretization
* Pedictor-corector time discretization with Runge-Kutta 4-th order
* Staircase-type internal free surface boundary condition without artificial viscosity contrast
* Geodynamically relevant visco-elasto-plastic rheology
* Local nonlinear solver based on FZERO algorithm
* Elastic stress rotation algorithm based on the time integration of the vorticity pseudo-vector
* Coupled velocity-pressure-temperature Newton-Raphson nonlinear solver
* Coupled velocity-pressure geometric multigrid preconditioner with Galerkin coarsening
Staggered grid finite difference, being inherently Eulerian and rather complicated discretization method (in terms of spatial arrangement of the unknowns) poses a number of challenges.
There is no natural treatment for free surface boundary condition. The solution based on the artificial quasi-viscous air phase introduces significant viscosity contrasts and spoils the convergence of the iterative solvers. In LaMEM we are currently implementing an approximate stair-case type of the free surface boundary condition which excludes the empty cells and restores the solver convergence.
Because of the mutual dependence of the stress and strain-rate tensor components, and their different spatial locations in the grid, there is no straightforward way of implementing the nonlinear rheology. In LaMEM we have developed and implemented an efficient interpolation scheme for the second invariant of the strain-rate tensor, that solves this problem.
Scalable efficient linear solvers are the key components of the successful nonlinear problem solution. In LaMEM we have a range of PETSc-based preconditioning techniques that either employ a block factorization of the velocity-pressure matrix, or treat it as a monolithic piece. In particular we have implemented the custom restriction-interpolation operators for the coupled Galerkin multigrid preconditioner. We have found that this type of algorithm is very robust with respect to the high grid resolutions and realistic viscosity variations. The coupled Galerking geometric multigrid implemented with the custom restriction-interpolation operators currently enables LaMEM to run efficiently with the grid sizes up to 1024*1024*512 cells on 16384 cores of the Blue Gene/Q JUQUEEN machine at the Jülich supercomputer center (Germany).
We also figured out the importance of the consistent linearization of the nonlinear residuals, efficient local iterations and accurate stress rotation algorithm. We are currently integrating all the building blocks into a single well-organized user-friendly package.