Towards locally refined 3D staggered grids for the Stokes equations

Many problems in computational geophysics can be answered by solving a discrete problem on a staggered grid - these define quantities on the elements, faces, edges, and vertices of a regular cell complex, embedded in space with an orthogonal coordinate system. A key example is the Stokes equations for creeping or highly-viscous fluid flow in Cartesian or spherical polar coordinates. These discretizations are attractive in many ways, both practically and computationally. They provide an intuitive way to pose the discrete problem, as they can be motivated as low-order finite volume methods or even as finite difference methods. They provide very regular data layout and relative easy of implementation, allowing for optimized and portable implementations in software. They couple particularly well with particle systems (as in the classical MAC scheme) which define material coefficients with arbitrary, non-grid-aligned interfaces; they provide first-order convergence with an extremely narrow stencil - in solving the discrete system, unknowns are coupled to fewer neighbors than in higher-order finite element or finite volume methods. However, especially in the context of geodynamical and seismic simulations, a severe limitation exists. Many open scientific questions relate to resolving localized features like shear bands, faults, subduction zone boundaries, and entrainment zones, and other boundary layers. With first order methods on regular grids, to resolve these fine spatial features, one must uniformly refine the entire domain (or resort to "Swiss cross" meshes, which are limited in the geometry they can resolve). These quickly becomes computationally intractable. One approach to resolving this tension is to define and use a non-uniform mesh. This is the basis for adaptive mesh refinement (AMR) methods, which are well-established in the context of finite element methods. Here, we present ongoing work in defining a stable, non-uniform, analogue to the MAC scheme for the 3D Stokes equations, attempting to retain the narrow stencil width characteristic of the uniform method while supporting arbitrary 2-to-1 refinement in a mesh described by an octree. Several researchers have made steps in this direction, so in this contribution we review those methods and discuss how they may be generalized or modified to provide a 3D scheme suitable for use in variable-viscosity Stokes flow.