SUPG, residual stabilization and variational multiscale methods

Fix an integer $N\ge 1$, put $h=1\mathrm{/}N$, and let $x_j=jh$ for $j=0,\dots ,N$. The cells are $K_j=[x_{j-1},x_j]$. For $g_0,g_1\in \mathbb{R}$, let $S_h(g_0,g_1)$ be the finite-dimensional affine space of continuous functions on $[0,1]$ that are affine on each cell and have endpoint values $g_0,g_1$. Write $V_h=S_h(0,0)$ for the homogeneous test and trial vector space; explicitly, $S_h(g_0,g_1)=\mathrm{\ell }+V_h$ with $\mathrm{\ell }(x)=(1-x)g_0+x{g}_1$. “P1” means polynomial degree at most one on each interval.

Let $ϵ>0$, $b\in \mathbb{R}\setminus \{0\}$, and $f\in C([0,1])$. The stationary boundary-value problem is to find $u\in C^2([0,1])$ such that $-ϵu^{\mathrm{\prime }\mathrm{\prime }}+b{u}^{\mathrm{\prime }}=f$ on $(0,1)$, with prescribed endpoint values $u(0)=g_0$ and $u(1)=g_1$. Its dimensionless cell Péclet number is $\mathrm{P}\mathrm{e}=\mathrm{\mid }b\mathrm{\mid }h\mathrm{/}(2ϵ)>0$. The stability theorem on this page uses the homogeneous branch $g_0=g_1=0$; the general endpoint notation also permits the boundary-layer example.

Fix $d\ge 2$ and a bounded connected polyhedral open set $\mathrm{\Omega }\subset {\mathbb{R}}^d$. Let $\{{\mathcal{T}}_h\}$ be a shape-regular family of finite matching affine simplicial meshes in the sense of Matching affine simplicial meshes and shape regularity: one finite constant $\sigma$ satisfies $h_K\mathrm{/}{\rho }_K\le \sigma$ for every cell of every mesh. Put $h={\max }_{K\in {\mathcal{T}}_h}{h}_K$. Let $V_h$ be the space of continuous functions on $\bar{\mathrm{\Omega }}$ that are affine on every cell and vanish on $\mathrm{\partial }\mathrm{\Omega }$.

Let $ϵ>0$, $b\in C^1(\bar{\mathrm{\Omega }};{\mathbb{R}}^d)$, $c\in C(\bar{\mathrm{\Omega }})$, and $f\in C(\bar{\mathrm{\Omega }})$. Define $B:=\mathrm{\parallel }b{\mathrm{\parallel }}_{L^{\mathrm{\infty }}(\mathrm{\Omega })}$ and $C:=\mathrm{\parallel }c{\mathrm{\parallel }}_{L^{\mathrm{\infty }}(\mathrm{\Omega })}$. Assume a fixed $\omega >0$ satisfies $c-\frac{1}{2}\mathrm{\nabla }\cdot b\ge \omega$ on $\mathrm{\Omega }$. The boundary condition for this branch is $u=0$ on the whole boundary, including the inflow part $\{x\in \mathrm{\partial }\mathrm{\Omega }:b(x)\cdot n(x)<0\}$. A classical exact solution is a function $u\in C^2(\bar{\mathrm{\Omega }})$ satisfying $-ϵ\mathrm{\Delta }u+b\cdot \mathrm{\nabla }u+cu=f$ in $\mathrm{\Omega }$ and this boundary condition.

Let $d\ge 1$ and $\mathrm{\Omega }\subset {\mathbb{R}}^d$ have a finite cell partition $\mathcal{T}$. Let $V_h$ be a finite-dimensional vector space of continuous vector-valued functions, each cellwise polynomial, and let $Q_h$ be a finite-dimensional scalar space of continuous cellwise polynomials with zero mean. Require all velocity components, their cellwise first derivatives and cellwise Laplacians, and all pressures and their cellwise first derivatives to belong to $L^2(\mathrm{\Omega })$. Boundary conditions are stated with each application; the Stokes theorem below uses zero velocity boundary values, and the periodic example below uses periodic identification. The pair is equal order when velocity components and pressure use the same polynomial degree on each cell. On a cell $K$, derivatives in residuals are classical derivatives restricted to the cell interior. For well-defined cell integrals, take $\sigma \in L^{\mathrm{\infty }}(\mathrm{\Omega })$, $\beta \in L^{\mathrm{\infty }}(\mathrm{\Omega };{\mathbb{R}}^d)$, $f\in L^2(\mathrm{\Omega };{\mathbb{R}}^d)$, and parameters ${\tau }_K\ge 0$, $\gamma \ge 0$.

For viscosity $\nu >0$, reaction coefficient $\sigma \ge 0$, advective field $\beta$, force $f$, and trial pair $(u_h,p_h)\in V_h\times Q_h$, define the strong momentum residual $R_{M,K}(u_h,p_h):=\sigma u_h+\beta \cdot \mathrm{\nabla }{u}_h-\nu {\mathrm{\Delta }}_K{u}_h+\mathrm{\nabla }{p}_h-f\mathrm{.}$ For test pairs $(v_h,q_h)$, the full residual PSPG mixed form is $B_h((u_h,p_h),(v_h,q_h)):=\nu (\mathrm{\nabla }{u}_h,\mathrm{\nabla }{v}_h{)}_{\mathrm{\Omega }}+(\beta \cdot \mathrm{\nabla }{u}_h+\sigma u_h,v_h{)}_{\mathrm{\Omega }}-(p_h,\mathrm{\nabla }\cdot v_h{)}_{\mathrm{\Omega }}+(q_h,\mathrm{\nabla }\cdot u_h{)}_{\mathrm{\Omega }}+{\sum }_K{\tau }_K(\mathrm{\nabla }{q}_h,\sigma u_h+\beta \cdot \mathrm{\nabla }{u}_h-\nu {\mathrm{\Delta }}_K{u}_h+\mathrm{\nabla }{p}_h{)}_K,$ with load $(f,v_h{)}_{\mathrm{\Omega }}+{\sum }_K{\tau }_K(\mathrm{\nabla }{q}_h,f{)}_K$. The grad-div addend in velocity tests $v_h$ is $\gamma {\sum }_K(\mathrm{\nabla }\cdot u_h,\mathrm{\nabla }\cdot v_h{)}_K$ with $\gamma \ge 0$. This definition alone asserts no inf-sup condition, pressure bound, energy estimate, or convergence result.

For piecewise affine velocities, ${\mathrm{\Delta }}_K{u}_h=0$ on each cell. Therefore the standard P1 PSPG residual omits the cellwise viscous Laplacian of the exact velocity; this fact does not justify a pressure-accuracy or convergence claim. Such conclusions require separate consistency and stabilization analysis.

Let $X_c$ and $X_f$ be finite-dimensional real vector spaces and set $X=X_c\oplus X_f$; elements are written $(u_c,u_f)$, with $u_c$ the resolved component and $u_f$ the unresolved component. For a real vector space $Y$, let $Y^{\ast }$ be its space of real linear functionals. Let $A_{cc}:X_c\to X_c^{\ast }$, $A_{cf}:X_f\to X_c^{\ast }$, $A_{fc}:X_c\to X_f^{\ast }$, and $A_{ff}:X_f\to X_f^{\ast }$ be linear maps, with $A_{ff}$ invertible, and let $F_c\in X_c^{\ast }$ and $F_f\in X_f^{\ast }$. The linear equation on $X$ has block form $\left(\begin{array}{cc}{A}_{cc}& A_{cf}\\ A_{fc}& A_{ff}\end{array}\right)\left(\begin{array}{c}{u}_c\\ u_f\end{array}\right)=\left(\begin{array}{c}{F}_c\\ F_f\end{array}\right)\mathrm{.}$ The fine residual for a resolved vector $v_c$ is $r_f(v_c):=F_f-A_{fc}{v}_c$. Exact elimination gives $u_f=A_{ff}^{-1}{r}_f(u_c)$ and the Schur equation $(A_{cc}-A_{cf}{A}_{ff}^{-1}{A}_{fc})u_c=F_c-A_{cf}{A}_{ff}^{-1}{F}_f\mathrm{.}$

A residual closure selects a linear map $T:X_f^{\ast }\to X_f$ to approximate the exact inverse $A_{ff}^{-1}:X_f^{\ast }\to X_f$. It defines ${\tilde{A}}_{cc}:=A_{cc}-A_{cf}T{A}_{fc}$, ${\tilde{F}}_c:=F_c-A_{cf}T{F}_f$, and the modeled fine component ${\tilde{u}}_f:=T{r}_f({\tilde{u}}_c)$, where ${\tilde{u}}_c$ solves ${\tilde{A}}_{cc}{\tilde{u}}_c={\tilde{F}}_c$. This is a finite-dimensional algebraic closure; an application to a PDE requires its own residual, spaces, invertibility, and approximation-defect analysis.

Use the uniform mesh, coefficients, and spaces of One-dimensional stationary convection–diffusion model and uniform P1 space. On every cell set $\tau =\frac{h}{2\mathrm{\mid }b\mathrm{\mid }}(\coth (\mathrm{P}\mathrm{e})-\frac{1}{\mathrm{P}\mathrm{e}}),\coth s:=\frac{\cosh s}{\sinh s}\mathrm{.}$ For $w\in S_h(g_0,g_1)+C^2([0,1])$ and $v_h\in V_h$, define $a_h(w,v_h):=ϵ{\int }_0^1{w}^{\mathrm{\prime }}{v}_h^{\mathrm{\prime }}dx+b{\int }_0^1{w}^{\mathrm{\prime }}{v}_{h}dx+{\sum }_{K\in {\mathcal{T}}_h}\tau {\int }_K(-ϵw^{\mathrm{\prime }\mathrm{\prime }}+b{w}^{\mathrm{\prime }})(b{v}_h^{\mathrm{\prime }})dx,$ $F_h(v_h):={\int }_0^{1}f{v}_{h}dx+{\sum }_{K\in {\mathcal{T}}_h}\tau {\int }_{K}f(b{v}_h^{\mathrm{\prime }})dx\mathrm{.}$ Second derivatives are taken classically on each cell. The strong equation residual is $R_K(w):=-ϵw^{\mathrm{\prime }\mathrm{\prime }}+b{w}^{\mathrm{\prime }}-f$. For P1 trial functions, $w_h^{\mathrm{\prime }\mathrm{\prime }}=0$ on each cell, so the cell term is $\tau {\int }_K(b{w}_h^{\mathrm{\prime }})(b{v}_h^{\mathrm{\prime }})$. On the smooth consistency extension it retains $-ϵw^{\mathrm{\prime }\mathrm{\prime }}$; thus the full residual pairing is $\tau {\int }_K(R_K(w)+f)(b{v}_h^{\mathrm{\prime }})$, while the load contains the matching $f$ term. The SUPG equation is $a_h(w_h,v_h)=F_h(v_h)$ and adds the same residual-weighted streamline test perturbation to the differential equation and its load.

Use the model and notation of Stationary multidimensional convection–diffusion–reaction P1 model. Define $s_b(K)=h_K\mathrm{/}\mathrm{\parallel }b{\mathrm{\parallel }}_{L^{\mathrm{\infty }}(K)}$ when that norm is positive and $s_b(K)=+\mathrm{\infty }$ otherwise; define $s_c(K)=\omega \mathrm{/}\mathrm{\parallel }c{\mathrm{\parallel }}_{L^{\mathrm{\infty }}(K)}^2$ when that norm is positive and $s_c(K)=+\mathrm{\infty }$ otherwise. For each cell set ${\delta }_K:=\theta \min \{\frac{h_K^2}{ϵ},s_b(K),s_c(K)\},0<\theta \le \frac{1}{4}\mathrm{.}$ The minimum is finite and positive because its first entry is finite and positive. This scaling gives ${\delta }_K\mathrm{\parallel }c{\mathrm{\parallel }}_{L^{\mathrm{\infty }}(K)}^2\le \omega \mathrm{/}4$.

For a trial function $w$ that is the sum of a cellwise smooth function and a member of $V_h$, and $v_h\in V_h$, define the element residual $R_K(w):=-ϵ{\mathrm{\Delta }}_{K}w+b\cdot \mathrm{\nabla }w+cw$, and $a_h(w,v_h):=ϵ(\mathrm{\nabla }w,\mathrm{\nabla }{v}_h{)}_{\mathrm{\Omega }}+(b\cdot \mathrm{\nabla }w+cw,v_h{)}_{\mathrm{\Omega }}+{\sum }_{K\in {\mathcal{T}}_h}{\delta }_K(R_K(w),b\cdot \mathrm{\nabla }{v}_h{)}_K,$ $F_h(v_h):=(f,v_h{)}_{\mathrm{\Omega }}+{\sum }_{K\in {\mathcal{T}}_h}{\delta }_K(f,b\cdot \mathrm{\nabla }{v}_h{)}_K\mathrm{.}$ Here $(g_1,g_2{)}_D:={\int }_D{g}_1{g}_{2}dx$ and ${\mathrm{\Delta }}_K$ is the classical Laplacian restricted to the interior of $K$. For $w_h\in V_h$, ${\mathrm{\Delta }}_K{w}_h=0$ because $w_h$ is affine on $K$. No stability claim is part of this definition.

Let $0=x_0<\cdots <x_N=1$ be a finite interval partition with maximum cell length $h$, and let $u\in C^2([0,1])$. Let $I_{h}u$ be the continuous piecewise-linear interpolant with the same nodal values as $u$. Then $\mathrm{\parallel }u-I_{h}u{\mathrm{\parallel }}_{L^2(0,1)}\le h^2\mathrm{\parallel }{u}^{\mathrm{\prime }\mathrm{\prime }}{\mathrm{\parallel }}_{L^{\mathrm{\infty }}(0,1)}\mathrm{.}$

Use the fixed dimension $d\ge 2$ and the finite shape-regular affine simplicial mesh family from Stationary multidimensional convection–diffusion–reaction P1 model and Matching affine simplicial meshes and shape regularity. Let $u\in C^2(\bar{\mathrm{\Omega }})$, and on each simplex $K$ let $I_{K}u={\sum }_{i=0}^{d}u(a_i){\lambda }_i$ be its nodal affine interpolant, where $a_i$ are the vertices and ${\lambda }_i$ the barycentric coordinates. Set $I_{h}u{\mathrm{\mid }}_K=I_{K}u$ and $M_2:=\mathrm{\mid }\mathrm{\Omega }{\mathrm{\mid }}^{1\mathrm{/}2}\mathrm{\parallel }{D}^{2}u{\mathrm{\parallel }}_{L^{\mathrm{\infty }}(\mathrm{\Omega })}$, using the Euclidean matrix norm for $D^{2}u$. There is a constant $C_I$ depending only on $d$ and the shape-regularity bound such that $\mathrm{\parallel }u-I_{h}u{\mathrm{\parallel }}_{L^2(\mathrm{\Omega })}\le C_I{h}^2{M}_2,\mathrm{\parallel }\mathrm{\nabla }(u-I_{h}u){\mathrm{\parallel }}_{L^2(\mathrm{\Omega })}\le C_{I}h{M}_2\mathrm{.}$ If $u=0$ on $\mathrm{\partial }\mathrm{\Omega }$, then $I_{h}u=0$ on every boundary facet.

Let $d\ge 2$, let $\mathrm{\Omega }\subset {\mathbb{R}}^d$ be a bounded connected polyhedral open set, and let a finite conforming nondegenerate simplicial mesh have connected cell adjacency. Let $V_h$ be continuous vector-valued P1 functions vanishing on the boundary, and let $Q_h$ be continuous scalar P1 functions with zero mean. Let $f\in L^2(\mathrm{\Omega };{\mathbb{R}}^d)$, $\nu >0$, ${\tau }_K>0$, and $\gamma \ge 0$. For the stationary Stokes specialization of Residual PSPG and grad-div terms for equal-order flow spaces, define $B_h((u_h,p_h),(v_h,q_h)):=\nu (\mathrm{\nabla }{u}_h,\mathrm{\nabla }{v}_h{)}_{\mathrm{\Omega }}-(p_h,\mathrm{\nabla }\cdot v_h{)}_{\mathrm{\Omega }}+(q_h,\mathrm{\nabla }\cdot u_h{)}_{\mathrm{\Omega }}+\gamma (\mathrm{\nabla }\cdot u_h,\mathrm{\nabla }\cdot v_h{)}_{\mathrm{\Omega }}+{\sum }_K{\tau }_K(\mathrm{\nabla }{q}_h,-\nu {\mathrm{\Delta }}_K{u}_h+\mathrm{\nabla }{p}_h{)}_K\mathrm{.}$ For the residual load, set $F_h(v_h,q_h):=(f,v_h{)}_{\mathrm{\Omega }}+{\sum }_K{\tau }_K(\mathrm{\nabla }{q}_h,f{)}_K$. Define $\mathrm{\parallel }(v_h,q_h){\mathrm{\parallel }}_h^2:=\nu \mathrm{\parallel }\mathrm{\nabla }{v}_h{\mathrm{\parallel }}_2^2+\gamma \mathrm{\parallel }\mathrm{\nabla }\cdot v_h{\mathrm{\parallel }}_2^2+{\sum }_K{\tau }_K\mathrm{\parallel }\mathrm{\nabla }{q}_h{\mathrm{\parallel }}_{0,K}^2\mathrm{.}$ Then this is a norm, $B_h(z_h,z_h)=\mathrm{\parallel }{z}_h{\mathrm{\parallel }}_h^2$, and the discrete problem $B_h(z_h,\cdot )=F_h(\cdot )$ has a unique solution satisfying $\mathrm{\parallel }{z}_h{\mathrm{\parallel }}_h\le \mathrm{\parallel }{F}_h{\mathrm{\parallel }}_{h,\ast }$, where $\mathrm{\parallel }{F}_h{\mathrm{\parallel }}_{h,\ast }:={\sup }_{y_h\mathrm{\ne }0}\mathrm{\mid }{F}_h(y_h)\mathrm{\mid }\mathrm{/}\mathrm{\parallel }{y}_h{\mathrm{\parallel }}_h$. The estimate controls pressure in the mesh-dependent gradient stabilization norm; it does not assert a mesh-uniform $L^2$ pressure bound or convergence.

Use the block operators and residual closure of Finite-dimensional resolved–unresolved VMS split and residual closure. Equip $X_c,X_f$ and their duals with norms, and use the corresponding induced operator norms. Let $G:=A_{ff}^{-1}$, assume the exact Schur system has a solution $u_c$, and put $r_f:=F_f-A_{fc}{u}_c$. Suppose the closure map $T$ has $\tilde{S}:=A_{cc}-A_{cf}T{A}_{fc}$ invertible, and let ${\tilde{u}}_c$ solve the approximate Schur system. Then $\mathrm{\parallel }{\tilde{u}}_c-u_c\mathrm{\parallel }\le \mathrm{\parallel }{\tilde{S}}^{-1}\mathrm{\parallel }\mathrm{\parallel }{A}_{cf}\mathrm{\parallel }\mathrm{\parallel }T-G\mathrm{\parallel }\mathrm{\parallel }{r}_f\mathrm{\parallel },$ and for the modeled fine component ${\tilde{u}}_f=T(F_f-A_{fc}{\tilde{u}}_c)$, $\mathrm{\parallel }{\tilde{u}}_f-u_f\mathrm{\parallel }\le \mathrm{\parallel }T-G\mathrm{\parallel }\mathrm{\parallel }{r}_f\mathrm{\parallel }+\mathrm{\parallel }T\mathrm{\parallel }\mathrm{\parallel }{A}_{fc}\mathrm{\parallel }\mathrm{\parallel }{\tilde{u}}_c-u_c\mathrm{\parallel }\mathrm{.}$ The exact fine scale is $u_f=G{r}_f$. These bounds concern finite-dimensional algebraic blocks; they do not certify a universal PDE subgrid model.

Let $s>0$ and define $m(s):=\coth s-1\mathrm{/}s$ and $c_{\ast }:=\min \{(3\sinh 1{)}^{-1},\coth 1-1\}>0$. Then $c_{\ast }\min \{s,1\}\le m(s)<1,m(s)\le s\mathrm{/}2.$ Consequently, for the parameter in One-dimensional residual SUPG form with the exact cell parameter, $0<\frac{c_{\ast }}{2}\min \{\frac{h^2}{2ϵ},\frac{h}{\mathrm{\mid }b\mathrm{\mid }}\}\le \tau <\frac{h}{2\mathrm{\mid }b\mathrm{\mid }},\tau \le \frac{h^2}{8ϵ}\mathrm{.}$

Use the model of Stationary multidimensional convection–diffusion–reaction P1 model and form/parameter of Multidimensional P1 SUPG residual form and parameter scaling. Assume Countable Choice, and let $u\in C^2(\bar{\mathrm{\Omega }})$ solve the stated strong problem. For every continuous, zero-boundary function $w$ whose restriction to each cell has square-integrable values and square-integrable weak first derivatives, define $\mathrm{\parallel }w{\mathrm{\parallel }}_{\mathrm{S}\mathrm{D},h}^2:=ϵ\mathrm{\parallel }\mathrm{\nabla }w{\mathrm{\parallel }}_{L^2(\mathrm{\Omega })}^2+\omega \mathrm{\parallel }w{\mathrm{\parallel }}_{L^2(\mathrm{\Omega })}^2+{\sum }_K{\delta }_K\mathrm{\parallel }b\cdot \mathrm{\nabla }w{\mathrm{\parallel }}_{L^2(K)}^2\mathrm{.}$ $M_2:=\mathrm{\mid }\mathrm{\Omega }{\mathrm{\mid }}^{1\mathrm{/}2}\mathrm{\parallel }{D}^{2}u{\mathrm{\parallel }}_{L^{\mathrm{\infty }}(\mathrm{\Omega })},B:=\mathrm{\parallel }b{\mathrm{\parallel }}_{L^{\mathrm{\infty }}(\mathrm{\Omega })},C:=\mathrm{\parallel }c{\mathrm{\parallel }}_{L^{\mathrm{\infty }}(\mathrm{\Omega })},{\delta }_{\ast }:={\max }_K{\delta }_K\mathrm{.}$ Then the discrete form is coercive: $a_h(v_h,v_h)\ge \frac{3}{4}\mathrm{\parallel }{v}_h{\mathrm{\parallel }}_{\mathrm{S}\mathrm{D},h}^2(v_h\in V_h),$ the discrete problem has a unique solution $u_h$, and the method is consistent, $a_h(u,v_h)=F_h(v_h)$. Moreover, with $C_{\ast }$ depending only on dimension, shape regularity, and the interpolation constant, $\mathrm{\parallel }u-u_h{\mathrm{\parallel }}_{\mathrm{S}\mathrm{D},h}\le C_{\ast }{M}_2[(\sqrt{ϵ}+\frac{B}{\sqrt{\omega }})h+(\sqrt{\omega }+\frac{C}{\sqrt{\omega }})h^2+B\sqrt{{\delta }_{\ast }}h+\sqrt{{\delta }_{\ast }}(\sqrt{d}ϵ+Bh+C{h}^2)]\mathrm{.}$ This is a P1 estimate for the stated strong regularity and parameter branch; it does not claim stability for unrestricted multidimensional parameters or weaker regularity.