Some Examples of ASM preconditioners

78. Some Examples of ASM preconditioners#

78.1. Diagonal preconditioner for $L_2$-norm#

Let $\Omega \subset {\mathbb R}^d$, and $M$ be the finite element discretization matrix of

\[ M(u,v) = \int_\Omega u v \, dx \]

on a conforming finite element space $V_h$ of order $p$. Then there hold the spectral estimates

\[ c_1 \, u^T C_{Jac} u \leq u^T M u \leq c_2 \, u^T C_{Jac} u. \]

with constants independent of the mesh-size $h$, but dependent on $p$.

Proof: The sub-space decomposition is

\[ V_h = \sum_{i=1}^N \operatorname{span} \{ \varphi_i \} \]

With the basis function $\varphi_i$ we have

\[ u = \sum_{i=1}^N u_i \varphi_i, \]

and the norms are

\[ \| u \|_M^2 = \big\| \sum_i u_i \varphi_i \big\|_{L_2(\Omega)}^2 \]

and

\[ \| u \|_{C_{Jac}}^2 = \sum_i \| u_i \varphi_i \|_{L_2(\Omega)}^2. \]

We use $\Omega = \cup \, T$ to split the integral into integrals on elements, and use that elements are mapped from the reference element $T = \Phi_T (\widehat T)$. We assume that $J_T = \operatorname{det} \Phi_T^\prime \approx h^d$ is element-wise constant (otherwise we have to replace it by the element-wise min/max). For any function $v \in L_2(\Omega)$ there holds by means of the transformation rule of integrals:

\[ \| v \|_{L_2(\Omega)}^2 = \sum_T \| v \|_{L_2(T)}^2 = \sum_T J_T \| v \circ \Phi_T \|_{L_2(\widehat T)}^2 \]

Thus we have

\[\begin{eqnarray*} \| u \|_M^2 & = & \sum_T J_T \, \big\| \sum_\alpha u_{T,\alpha} \widehat \varphi_\alpha \big\|_{L_2(\widehat T)}^2, \\ \| u \|_{C_{Jac}}^2 & = & \sum_T J_T \, \sum_\alpha \| u_{T,\alpha} \widehat \varphi_\alpha \|_{L_2(\widehat T)}^2, \end{eqnarray*}\]

where the basis function $\varphi_{i|T} \circ \Phi_T = \widehat \varphi_\alpha$, and $u_{T,\alpha}$ the corresponding coefficient.

Since all norms in ${\mathbb R}^n$ (with a fix n !!!) are equivalent we get

\[ \big\| \sum_\alpha u_{T,\alpha} \widehat \varphi_\alpha \big\|_{L_2(\widehat T)}^2 \approx \sum_\alpha \| u_{T,\alpha} \widehat \varphi_\alpha \|_{L_2(\widehat T)}^2. \]

This holds since both sides are squared norms for $(u_{T,\alpha}) \in {\mathbb R}^{N_T}$

from ngsolve import *
from ngsolve.la import EigenValues_Preconditioner
mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))
fes = H1(mesh, order=3)
u,v = fes.TnT()
a = BilinearForm(u*v*dx).Assemble()
c = a.mat.CreateSmoother()

lam = list(EigenValues_Preconditioner(a.mat, c))
print ("lammin, lammax=", lam[0], lam[-1])
print ("kappa=", lam[-1]/lam[0])

lammin, lammax= 0.020831301863592395 4.762472029308668
kappa= 228.6209503608706

Exercise: Do some experiments with the mesh-size $h$ and the order $p$.

78.2. Diagonal preconditioner for the $H^1$ norm#

Now let $V_h \subset H_0^1(\Omega)$, and the bilinear-form

\[ A(u,v) = \int_\Omega \nabla u \nabla v \, dx. \]

Then there holds

\[ c_1 h^2 \, u^T C_{Jac} u \leq u^T A u \leq c_2 \, u^T C_{Jac} u, \]

which implies $\kappa(C_{Jac}^{-1} A) \preceq h^{-2}$.

Proof: We prove that

\[ \| u \|_M^2 \preceq \| u \|_A^2 \preceq h^{-2} \, \| u \|_M^2. \]

Then there follows immediately via the ASM-lemma

\[ \| u \|_{C_{Jac}^M}^2 \preceq \| u \|_{C_{Jac}^A}^2 \preceq h^{-2} \, \| u \|_{C_{Jac}^M}^2 \]

and

\[ \| u \|_A^2 \preceq h^{-2} \, \| u \|_M^2 \preceq h^{-2} \, \| u \|_{C_{Jac}^M}^2 \preceq h^{-2} \, \| u \|_{C_{Jac}^A}^2. \]

The claim $\| u \|_M^2 \preceq \| u \|_A^2$ is exactly the Friedrichs inequality. The other way around is an inverse inequality (we call an inequality on a finite element space inverse inequality if the constant deteriorates for $h \rightarrow 0$). Since

\[ \| u \|_A^2 = \sum_T \| \nabla u_{|T} \|_{L_2(T)}^2 \qquad \text{and} \qquad \| u \|_M^2 = \sum_T \| u_{|T} \|_{L_2(T)}^2 \]

it is enough to show

\[ \| \nabla u \|_{L_2(T)}^2 \preceq h^{-2} \| u \|_{L_2(T)}^2. \]

Again we use $T = \Phi(\widehat T)$ and write $\widehat u = u \circ \Phi$. From the chain-rule we get $\nabla u = { \Phi^\prime }^{-1} \nabla \widehat u$, and together with the transformation rule we get

\[ \| \nabla u \|_{L_2(T)}^2 \approx h^{d-2} \, \| \nabla \widehat u \|_{L_2(\hat T)}^2 \]

On the reference element, $\| \nabla \widehat u \|_{L_2(\hat T)}$ is a semi-norm on a fixed, finite dimensional space, and thus bounded by any norm and so also for the $L_2$-norm. Together we have proven that

\[ \| \nabla u \|_{L_2(T)}^2 \approx h^{d-2} \| \nabla \widehat u \|_{L_2(\widehat T)}^2 \preceq h^{d-2} \, \| \widehat u \|_{L_2(\widehat T)}^2 \preceq h^{-2} \, \| u \|_{L_2(T)}^2 \]

mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))
fes = H1(mesh, order=1, dirichlet=".*")
u,v = fes.TnT()
a = BilinearForm(grad(u)*grad(v)*dx).Assemble()
c = a.mat.CreateSmoother(fes.FreeDofs())

lam = list(EigenValues_Preconditioner(a.mat, c))
print ("lammin, lammax=", lam[0], lam[-1])
print ("kappa=", lam[-1]/lam[0])

lammin, lammax= 0.04897732520653299 1.4967019608058028
kappa= 30.559079216644538

Exercise:

is the Dirichlet boundary important ?
add some term $\int_\Gamma u v \, ds$ to the bilinear-form

Numerical experiments indicate that the estimate $\kappa \{C_{Jac}^{-1}A \} = O(h^{-2})$ is asymptotically sharp. To verify it, take a smooth function such that $\| \nabla u \|_{\Omega} \approx 1$. Then $\| u \|_{C}^2 = \sum u_i^2 \| \nabla \varphi_i \|^2$. There holds $u_i \approx 1$ and $\| \nabla \varphi_i \|^2 \approx h^{d-2}$, and there are $h^{-d}$ summands, which verifies that the estimate is sharp.

78.3. $H^1$-norm with small $L_2$-term#

We now consider $V_h \subset H^1$, and the bilinear-form

\[ A(u,v) = \int_\Omega \nabla u \nabla v \, dx + \varepsilon \int_\Omega u v \, dx. \]

In the limit $\varepsilon = 0$ the bilinear-form has a non-trivial kernel, namely the space of constant functions. Then the matrix is only semi-definite. For small $\varepsilon$ we expect a bad condition number.

mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))
fes = H1(mesh, order=1)
u,v = fes.TnT()
eps = 0.01
a = BilinearForm(grad(u)*grad(v)*dx + eps*u*v*dx ).Assemble()
c = a.mat.CreateSmoother(fes.FreeDofs())

lam = list(EigenValues_Preconditioner(a.mat, c))
print ("lammin, lam2, lammax=", lam[0], lam[1], lam[-1])
print ("kappa=", lam[-1]/lam[0])

lammin, lam2, lammax= 2.4456256298405192e-05 0.024384741600227022 1.6917198310276222
kappa= 69173.29497965478

The largest eigenvalue is good, but the smallest deteriorates with $\varepsilon$. This is conform with the analysis of the $H^1$-case, since now $$ \| u \|_{L_2}^2 \preceq \varepsilon^{-1} \, \| u \|_A^2 $$

We can fix the problem by adding the null-space function to the sub-space decomposition:

\[ V_h = \operatorname{span} \{ 1 \} + \sum \operatorname{span} \{\varphi_i\} \]

This sub-space splitting leads to the preconditioner

\[ C_{ASM}^{-1} = E_0 \, (E_0^T A E_0)^{-1} E_0^T + C_{Jac}^{-1}, \]

where the column vector of $E_0 \in {\mathbf R}^{N \times 1}$ is the coefficient vector of the constant function $u=1$.

We obtain the vector $E_0$ by setting some GridFunction to the function 1.

mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))
fes = H1(mesh, order=1)
u,v = fes.TnT()
eps = 1e-2
a = BilinearForm(grad(u)*grad(v)*dx + eps*u*v*dx ).Assemble()
cjac = a.mat.CreateSmoother(fes.FreeDofs())

gfconst = GridFunction(fes)
gfconst.Set(1)
e0 = BaseMatrix(gfconst.vec)   # N*1 matrix
a0 = InnerProduct(a.mat*gfconst.vec, gfconst.vec)
print ("a0=", a0)
c = cjac + 1/a0 * e0 @ e0.T

lam = list(EigenValues_Preconditioner(a.mat, c))
print ("lammin, lam2, lammax=", lam[0], lam[1], lam[-1])
print ("kappa=", lam[-1]/lam[0])

a0= 0.010000000000006035
lammin, lam2, lammax= 0.024547976218031113 0.042529997182085555 1.6917198260952784
kappa= 68.9148388881307

The rank-1 correction fixed the problem !

The analysis uses the ASM representation

\[ \| u \|_{C_{ASM}}^2 = \inf_{u = u_0 + \sum u_i } \| u_0 \|_{A}^2 + \sum_i \| u_i \|_A^2, \]

where $u_0$ is a constant function, and $u_i \in \operatorname{span} \{ \varphi_i \}$. We perform the decomposition in two steps, first we decompose $u = u_0 + u_f$, and then we decompose $u_f = \sum u_i$ into the one-dimensional spaces of the basis functions. Thus

\[ \| u \|_{C_{ASM}}^2 = \inf_{u = u_0 + u_f } \left\{ \| u_0 \|_{A}^2 + \inf_{u_f = \sum u_i } \sum_i \| u_i \|_A^2 \right\} \]

The second term is exactly the norm induced by the Jacobi-preconditioner of $u_f$:

\[ \| u \|_{C_{ASM}}^2 = \inf_{u = u_0 + u_f } \| u_0 \|_{A}^2 + \| u_f \| _{C_{Jac}}^2 \]

The estimate $A \preceq C_{ASM}$ follows from the triangle inequality, and the result of the Jacobi preconditioner as follows: Let $u_0$ be an arbitrary constant function. Then

\[ \| u \|_A^2 \leq 2 \big( \| u_0 \|_A^2 + \| u - u_0 \|_A^2 \big) \preceq \| u_0 \|_A^2 + \| u - u_0 \|_{C_{Jac}}^2. \]

Since $u_0$ is an arbitrarily chosen function, it also holds for the infimum.

The more difficult part is to show the existence of such a decomposition. Let $u \in V_h$ be given. Define the mean value

\[ u_0 = \frac{1}{|\Omega|} \int_\Omega u \, dx \]

There holds $\| u_0 \|_{L_2} \leq \| u \|_{L_2}$, and thus also

\[ \| u_0 \|_A^2 \leq \| u \|_A^2. \]

Setting $u_f = u - u_0$, we also get from the triangle inequality

\[ \| u_f \|_A^2 \preceq \| u \|_A^2. \]

The $A$-norm has only weak control of the mean value. But now we have that $u_f$ has mean-value zero. We apply the Poincaré inequality to obtain also control of the $L_2$-norm of $u_f$:

\[ \| u_f \|_{L_2}^2 + \| u_f \|_A^2 \preceq \| u \|_A^2 \]

We proceed as in the previous case:

\[ \| u_f \|_{C_{Jac}^A}^2 \preceq h^{-2} \| u_f \|_{C_{Jac}^M}^2 \preceq h^{-2} \, \| u_f \|_M^2 \preceq h^{-2} \, \| u \|_A^2 \]

Thus, all terms of the ASM-splitting are controlled by the $A$-norm of $u$.

Some Examples of ASM preconditioners

Contents

78. Some Examples of ASM preconditioners#

78.1. Diagonal preconditioner for \(L_2\)-norm#

78.2. Diagonal preconditioner for the \(H^1\) norm#

78.3. \(H^1\)-norm with small \(L_2\)-term#