64. Some Examples of ASM preconditioners

64.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\).

64.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.

64.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\).