Highly sparse surface couplings for subdomain-wise isoviscous Stokes finite element discretizations

TL;DR

This paper introduces a novel discretization approach for the Stokes system with variable viscosity, optimizing computational efficiency and accuracy in large-scale, multi-phase flow simulations by localizing complex calculations.

Contribution

The authors propose a new method that maintains the benefits of gradient-based decoupling in isoviscous regions while correctly handling cross derivatives in variable viscosity scenarios.

Findings

Faster multigrid convergence in parallel computations.

Maintains physical accuracy with local stencil modifications.

Demonstrated effectiveness in geophysical flow simulations.

Abstract

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.