  The Stokes system with constant viscosity can be cast into different
formulations by exploiting the incompressibility constraint. For instance the
strain in the weak formulation can be replaced by the gradient to decouple the
velocity components in the different coordinate directions. Thus the
discretization of the simplified problem leads to fewer nonzero entries in the
stiffness matrix. This is of particular interest in large scale simulations
where a reduced memory bandwidth requirement can help to significantly
accelerate the computations. In the case of a piecewise constant viscosity, as
it typically arises in multi-phase flows, or when the boundary conditions
involve traction, the situation is more complex, and one has to treat the cross
derivatives in the original Stokes system with care. A naive application of the
standard vectorial Laplacian results in a physically incorrect solution, while
formulations based on the strain increase the computational effort everywhere,
even when the inconsistencies arise only from an incorrect treatment in a small
fraction of the computational domain. Here we propose a new approach that is
consistent with the strain-based formulation and preserves the decoupling
advantages of the gradient-based formulation in isoviscous subdomains. The
modification is equivalent to locally changing the discretization stencils,
hence the more expensive discretization is restricted to a lower dimensional
interface, making the additional computational cost asymptotically negligible.
We demonstrate the consistency and convergence properties of the method and
show that in a massively parallel setup, the multigrid solution of the
resulting discrete systems is faster than for the classical strain-based
formulation. Moreover, we give an application example which is inspired by
geophysical research.
