LOD Methodology for Heterogeneous Stokes

• Localized Orthogonal Decomposition is a computational framework that constructs locally supported, problem-adapted basis functions for simulating PDEs with heterogeneous coefficients.
• It employs a coarse-fine scale decomposition with localized corrector functions to capture small-scale variations, ensuring accurate multiscale approximation.
• For heterogeneous Stokes problems, the LOD approach achieves optimal convergence rates and robust error control independent of coefficient regularity.

A Localized Orthogonal Decomposition (LOD) methodology is a computational framework for multiscale numerical approximation of partial differential equations (PDEs) with rough, heterogeneous coefficients. By constructing problem-adapted, locally supported basis functions on a coarse mesh, LOD enables accurate simulation of phenomena characterized by small-scale variability without the need to resolve all fine-scale features globally. For heterogeneous Stokes problems, the approach achieves optimal convergence rates in velocity and pressure approximations with errors independent of coefficient regularity, based on localized enrichment of coarse velocity spaces and a piecewise constant pressure treatment.

1. Variational Principle for Heterogeneous Stokes Problems

Consider a Lipschitz domain $\Omega \subset \mathbb{R}^n$ ($n=2,3$), with heterogeneous viscosity $\nu$ ($0 < \nu_{\min} \leq \nu \leq \nu_{\max}$) and an optional drag coefficient $\sigma$ ($0 \leq \sigma \leq \sigma_{\max}$). The strong form of the stationary Stokes equations is

$$\begin{cases} -\operatorname{div}( \nu \nabla u ) + \sigma u + \nabla p = f,\ \operatorname{div} u = 0,\ u|_{\partial\Omega} = 0,\ \int_\Omega p = 0, \end{cases}$$

where $u: \Omega \rightarrow \mathbb{R}^n$ is the velocity, $p: \Omega \rightarrow \mathbb{R}$ is the pressure, and $f \in L^2(\Omega)^n$ is the force.

The weak (variational) formulation: Find $n=2,3$0 with $n=2,3$1 and $n=2,3$2 such that

$n=2,3$3

where

$n=2,3$4

Well-posedness follows from the coercivity of $n=2,3$5 and the inf–sup condition on $n=2,3$6.

2. Coarse and Fine Spaces, Interpolation, and the Fine-Scale Decomposition

The LOD construction employs a quasi-uniform, simplicial coarse mesh $n=2,3$7 (diameter $n=2,3$8). The coarse velocity space is

$n=2,3$9

where $\nu$0 denotes the space of piecewise constant functions on $\nu$1. The discrete pressure space is

$\nu$2

A linear interpolation operator $\nu$3 is defined to preserve all interior face averages, ensuring that $\nu$4 is compatible with the divergence structure. The fine-scale supplement is

$\nu$5

In practical computations over unresolved oscillations, one uses a sufficiently refined fine mesh $\nu$6 and constructs $\nu$7 as a Crouzeix–Raviart discrete space.

3. Corrector Construction and Localized Enrichment

Given a coarse function $\nu$8, the LOD defines the global (ideal) corrector $\nu$9 via

$0 < \nu_{\min} \leq \nu \leq \nu_{\max}$0

This is the unique $0 < \nu_{\min} \leq \nu \leq \nu_{\max}$1-orthogonal projection of $0 < \nu_{\min} \leq \nu \leq \nu_{\max}$2 onto $0 < \nu_{\min} \leq \nu \leq \nu_{\max}$3. The enriched, problem-adapted trial function is $0 < \nu_{\min} \leq \nu \leq \nu_{\max}$4. However, as $0 < \nu_{\min} \leq \nu \leq \nu_{\max}$5 is globally supported, direct computation is infeasible.

To localize, decompose $0 < \nu_{\min} \leq \nu \leq \nu_{\max}$6 into a sum of local basis functions (associated to interior faces $0 < \nu_{\min} \leq \nu \leq \nu_{\max}$7 and coordinate directions), and solve for the correctors of these:

• For coarse element $0 < \nu_{\min} \leq \nu \leq \nu_{\max}$8, define the $0 < \nu_{\min} \leq \nu \leq \nu_{\max}$9-layer patch $\sigma$0 via

$\sigma$1

• The localized corrector $\sigma$2 solves the same variational equation as $\sigma$3 but with test functions $\sigma$4 supported in $\sigma$5.
• The localized corrector $\sigma$6 is then defined on $\sigma$7 by $\sigma$8.

Exponential decay of the correctors is established: For each local basis $\sigma$9,

$0 \leq \sigma \leq \sigma_{\max}$0

enabling sharp control of localization error by adjusting the oversampling parameter $0 \leq \sigma \leq \sigma_{\max}$1.

4. Multiscale Space, Discrete Problem, and Divergence-Free Constraint

The multiscale velocity space is constructed as

$0 \leq \sigma \leq \sigma_{\max}$2

where $0 \leq \sigma \leq \sigma_{\max}$3 is the set of interior faces and $0 \leq \sigma \leq \sigma_{\max}$4 is a local face-based basis function in direction $0 \leq \sigma \leq \sigma_{\max}$5. The full discrete LOD approximation seeks $0 \leq \sigma \leq \sigma_{\max}$6 satisfying

$0 \leq \sigma \leq \sigma_{\max}$7

This structure guarantees $0 \leq \sigma \leq \sigma_{\max}$8-free velocities in the discrete sense due to face-average conformity.

5. Convergence Analysis and Localization Parameter Selection

A priori error estimates for the LOD Stokes method (under minimal regularity and no scale separation) are: $0 \leq \sigma \leq \sigma_{\max}$9 where $$\begin{cases} -\operatorname{div}( \nu \nabla u ) + \sigma u + \nabla p = f,\ \operatorname{div} u = 0,\ u|_{\partial\Omega} = 0,\ \int_\Omega p = 0, \end{cases}$$0 is a locally post-processed pressure using the fine-scale pressures from corrector solves. By choosing $$\begin{cases} -\operatorname{div}( \nu \nabla u ) + \sigma u + \nabla p = f,\ \operatorname{div} u = 0,\ u|_{\partial\Omega} = 0,\ \int_\Omega p = 0, \end{cases}$$1, one ensures that the exponential terms become higher order in $$\begin{cases} -\operatorname{div}( \nu \nabla u ) + \sigma u + \nabla p = f,\ \operatorname{div} u = 0,\ u|_{\partial\Omega} = 0,\ \int_\Omega p = 0, \end{cases}$$2; optimal convergence rates ($$\begin{cases} -\operatorname{div}( \nu \nabla u ) + \sigma u + \nabla p = f,\ \operatorname{div} u = 0,\ u|_{\partial\Omega} = 0,\ \int_\Omega p = 0, \end{cases}$$3 in $$\begin{cases} -\operatorname{div}( \nu \nabla u ) + \sigma u + \nabla p = f,\ \operatorname{div} u = 0,\ u|_{\partial\Omega} = 0,\ \int_\Omega p = 0, \end{cases}$$4 for velocity, $$\begin{cases} -\operatorname{div}( \nu \nabla u ) + \sigma u + \nabla p = f,\ \operatorname{div} u = 0,\ u|_{\partial\Omega} = 0,\ \int_\Omega p = 0, \end{cases}$$5 in $$\begin{cases} -\operatorname{div}( \nu \nabla u ) + \sigma u + \nabla p = f,\ \operatorname{div} u = 0,\ u|_{\partial\Omega} = 0,\ \int_\Omega p = 0, \end{cases}$$6 for velocity, $$\begin{cases} -\operatorname{div}( \nu \nabla u ) + \sigma u + \nabla p = f,\ \operatorname{div} u = 0,\ u|_{\partial\Omega} = 0,\ \int_\Omega p = 0, \end{cases}$$7 in $$\begin{cases} -\operatorname{div}( \nu \nabla u ) + \sigma u + \nabla p = f,\ \operatorname{div} u = 0,\ u|_{\partial\Omega} = 0,\ \int_\Omega p = 0, \end{cases}$$8 for pressure) are then achieved.

6. Implementation and Computational Aspects

• Fine mesh $$\begin{cases} -\operatorname{div}( \nu \nabla u ) + \sigma u + \nabla p = f,\ \operatorname{div} u = 0,\ u|_{\partial\Omega} = 0,\ \int_\Omega p = 0, \end{cases}$$9 must resolve the smallest scale of coefficient variation ($u: \Omega \rightarrow \mathbb{R}^n$0).
• Commonly: Crouzeix–Raviart elements for velocity and $u: \Omega \rightarrow \mathbb{R}^n$1 for pressure on $u: \Omega \rightarrow \mathbb{R}^n$2.
• For each face $u: \Omega \rightarrow \mathbb{R}^n$3 and direction $u: \Omega \rightarrow \mathbb{R}^n$4, assemble and solve the local saddle-point problem for $u: \Omega \rightarrow \mathbb{R}^n$5 on the patch $u: \Omega \rightarrow \mathbb{R}^n$6 with homogeneous boundary conditions.
• All patch problems are independent (embarrassingly parallel).
• Global stiffness matrix assembly uses the basis $u: \Omega \rightarrow \mathbb{R}^n$7 (size $u: \Omega \rightarrow \mathbb{R}^n$8).
• Complexity per patch: $u: \Omega \rightarrow \mathbb{R}^n$9 fine unknowns; for $p: \Omega \rightarrow \mathbb{R}$0 this scales quasi-linearly in $p: \Omega \rightarrow \mathbb{R}$1.
• The global (coarse) solve in the multiscale space is relatively low-dimensional and well-conditioned due to face-based basis structure.

Typical choices and recommendations:

• For the interpolation $p: \Omega \rightarrow \mathbb{R}$2, admissible options include face-preserving, Clément, or Scott–Zhang operators.
• The choice of fine finite element (e.g., CR, $p: \Omega \rightarrow \mathbb{R}$3–$p: \Omega \rightarrow \mathbb{R}$4) for the local Stokes saddle-point problem can be tailored as needed.
• Energy-orthogonal splitting and exponential corrector decay are invariant across these implementation choices.

7. Numerical Evidence and Observed Properties

Numerical experiments confirm that the LOD for Stokes problems achieves the following:

• Robust approximation independent of the heterogeneity structure; no dependence on coefficient smoothness or scale separation.
• Observed rates in practice match theoretical predictions: $p: \Omega \rightarrow \mathbb{R}$5 in $p: \Omega \rightarrow \mathbb{R}$6-norm for velocity, $p: \Omega \rightarrow \mathbb{R}$7 in $p: \Omega \rightarrow \mathbb{R}$8-norm, and $p: \Omega \rightarrow \mathbb{R}$9 for (post-processed) pressure.
• Localization error is sharply controlled: moderate oversampling, $f \in L^2(\Omega)^n$0, suffices even in the presence of highly oscillatory and high-contrast coefficients.
• Computational workload is parallelizable on element faces, with trivial coordination overhead.

---

The Localized Orthogonal Decomposition methodology, as applied to the Stokes problem, establishes a framework for multiscale velocity-pressure approximation with optimal rates and robust error control, independent of coefficient regularity and suitable for high-performance, parallel computation (Hauck et al., 2024). This approach generalizes to related vector-valued, nonscalar, and saddle-point PDEs with strong theoretical and practical guarantees.