76. Parameter Dependent Problems

We have observed that adding the constraint to the A-block may improve the convergence of saddle-point solvers. If we add the constraint with a large factor \(1/\varepsilon\), we may consider the augmented problem as penalty formulation for the constrained minimization problem, and skip the second equation and the Lagrange parameter completely. We call this a parameter dependent problem:

\[ \big(A + \frac{1}{\varepsilon} B^T C^{-1} B \big) u = f \]

Let \(V_0 = \operatorname{ker} B\) the null-space of the \(B\)-matrix.

On the null-space \(V_0\) only the \(A\)-operator is seen. On its complement \(V_0^\bot\), the second part \(\frac{1}{\varepsilon} B^T C^{-1} B\) dominates for small \(\varepsilon\). Unless \(V_0\) is trivial, the condition number degenerates with \(O(\varepsilon^{-1})\).

The goal is to develop preconditioners for the parameter dependent problems which are robust in \(\varepsilon\).

76.1. Dirichlet boundary conditions by penalty:

Adding the Dirichlet condition \(u = u_D\) on \(\Gamma_D\):

\[ \int \nabla u \nabla v + \int_{\Gamma_D} \alpha u v = \int f v + \int_{\Gamma_D} \alpha u_D v \]

with large parameter \(\alpha = \varepsilon^{-1}\), or depending on the mesh-size \(\alpha = h^{-1}\).

The null-space \(V_0\) is the finite element sub-space of \(H_{0,D}^1 = \{ v \in H^1 : v = 0 \text{ on } \Gamma_D \}\).

76.2. Penalty formulation for the flux:

Approximating the mixed formulation \(\sigma = \nabla u\), and \(\operatorname{div} \sigma = -f\) by penalty leads to:

\[ \int \sigma \tau + \frac{1}{\varepsilon} \int \operatorname{div} \sigma \operatorname{div} \tau = \frac{-1}{\varepsilon} \int f \operatorname{div} \tau \]

This is a variational formulation on the Hilbert-space \(H(\operatorname{div}) := \{ \tau \in [L_2]^d : \operatorname{div} \tau \in L_2\}\).

The null-space \(V_0\) consists of all \(\operatorname{div}\)-free functions.

76.3. Maxwell equations:

Very similar is the equation for the magnetic field, a special case of the Maxwell equations. Given a divergence-free current density \(j\), search for a divergence-free magnetic field \(B\) such that

\[ \operatorname{curl} \mu^{-1} B = j, \]

where \(\mu\) is the so called permeability, a kind of conductivity for the magnetic field. Since we search for a \(\operatorname{div}\)-free \(B\), and assuming a simply connected domain, we may introduce a vector-potential \(u\) such that \(B = \operatorname{curl} u\). Together we have the second order equation

\[ \operatorname{curl} \mu^{-1} \operatorname{curl} u = j. \]

Its weak formulation (skipping discussion of boundary conditions) is

\[ \int \mu^{-1} \operatorname{curl} u \, \operatorname{curl} v = \int j v \qquad \forall \, v \]

There is no unique solution: adding an arbitrary gradient field to a solution is also a solution. One way out is Coloumb - gauging, where we add the constraint \(\int u \nabla \psi = 0\), i.e. a weak formulation of \(\operatorname{div} u = 0\). Another, very simple method is to regularize with a small \(L_2\)-term:

\[ \int \mu^{-1} \operatorname{curl} u \, \operatorname{curl} v + \int \varepsilon u v = \int j v \]

Again, we have a parameter dependent problem, with the fictitious regularization parameter \(\varepsilon\). The canonical space is

\[ H(\operatorname{curl}) = \{ u \in [L_2]^3 : \operatorname{curl} u \in [L_2]^3 \}, \]

equipped with the norm

\[ \| u \|_{H(\operatorname{curl})}^2 = \| u \|_{L_2}^2 + \| \operatorname{curl} u \|_{L_2}^2 \]

The null-space is

\[ V_0 = \{ v : \operatorname{curl} v = 0 \} = \nabla H^1 \]

76.4. Penalty formulation for the Stokes equation:

We approximate the div-free condition by penalty:

Find \(u \in [H_{0,D}]^d\) such that

\[ \int \nu \nabla u : \nabla v + \frac{1}{\varepsilon} \int \operatorname{div} u \operatorname{div} v = \int f v \]

It looks similar as the penalty formulation for the flux, but now the A-form is a \(H^1\)-inner product instead of an \(L_2\)-inner product. Now, conforming low order methods lead to bad approximations (called locking effect). Robust discretizations are constructed and analyzed by by mixed formulations. If the pressure space consists of discontinuous finite elements, and the \(C\) matrix is chosen as \(L_2\)-inner product, it can be cheaply inverted. We end up with the discretization of the parameter dependent problem. However, the detour over the mixed mixed formulation leads to a projection operator \(R_h\):

Find \(u_h \in V_h \subset [H_{0,D}]^d\) such that

\[ \int \nu \nabla u_h : \nabla v_h + \frac{1}{\varepsilon} \int R_h \operatorname{div} u_h \operatorname{div} v_h = \int f v_h \]

A stable example is to choose \(V_h\) a \(P^2\) finite element space. Then \(\operatorname{div} u_h\) is a \(P^1\) function. Choosing a \(P^0\) pressure leads to a projection operator \(R_h\) which is the \(L_2\)-orthogonal projection from \(P^1\) onto \(P^0\).

76.5. Robust two-level methods for parameter dependent problems

Let

\[ A_h^\varepsilon = A + \frac{1}{\varepsilon} B^T C^{-1} B \]

be the discretization matrix of the parameter dependent problem on the finite element space \(V_h\). We are going to analyze two-level preconditioners

\[ C_{2L}^{-1} = D_h^{-1} + P \, ( A_H^\varepsilon ) ^{-1} P^T, \]

where \(A_H^\varepsilon\) is the discretization matrix on the coarse space \(V_H\), \(P : V_H \rightarrow V_h\) is the prolongation matrix, and \(D_h\) is a local block-Jacobi preconditioner on \(V_h\).

76.5.1. Robust smoothers

By the ASM - lemma, we can represent the norm generated by the smoother \(D_h\) as

\[ \| u \|_{D_h}^2 = \inf_{u = \sum u_i} \sum \| u_i \|_{A + \frac{1}{\varepsilon} B^T C^{-1} B}^2, \]

where we have the sub-space splitting

\[ V_h = \sum V_i. \]

If \(u\) is a function in the finite-element null-space \(V_0\), then its norm \(\| u \|_{A + \frac{1}{\varepsilon} B^T C^{-1} B}^2\) is independent of the small parameter \(\varepsilon\). Only if we can split \(u\) into local null-space functions \(u_i\), the \(D_h\) norm is also independent of \(\varepsilon\). Thus, we want a splitting compatible with the null-space:

\[ V_0 = \sum V_i \cap V_0 \]

Example: Maxwell equations

For a compatible discretization with Nedelec finite elements, the discrete null-space consists of discrete gradients:

\[ \{ v_h \in V_h : \operatorname{curl} v_h = 0 \} = \{ \nabla \psi_h : \psi_h \in W_h \subset H^1 \} \]

Let \(u_h \in V_0\), and \(\psi_h\) such that \(\nabla \psi_h = u_h\). Now, decompose \(\psi_h = \sum \psi_i\). Then, \(\sum \nabla \psi_i\) is a decomposition of \(u_h\) into null-space functions \(u_i = \nabla \psi_i\). The block-Jacobi must consist of blocks large enough such that each component \(\nabla \psi_i\) is contained in one block. If \(\psi_i\) is a vertex-hat basis function, then its gradient is contained in the vertex patch. These block-smoothers are called AFW-blocks, due to Arnold, Falk, and Winther: Multigrid in H(div) and H(curl), Numerische Mathematik, 2000.

Lets try:

from ngsolve import *
from netgen.csg import *
from ngsolve.webgui import Draw

geo = CSGeometry()
magnet = Cylinder(Pnt(-1,0,0), Pnt(1,0,0), 0.3) \
    * OrthoBrick(Pnt(-1,-1,-1), Pnt(1,1,1))
box = OrthoBrick(Pnt(-3, -3, -3), Pnt(3,3,3))
geo.Add (magnet.mat("magnet"))
geo.Add ((box-magnet).mat("air"))
mesh = Mesh(geo.GenerateMesh(maxh=1))
# mesh.ngmesh.Refine()

clipping = { "pnt" : (0,0,0), "function" : False}
Draw (mesh, clipping=clipping);

fes = HCurl(mesh, order=0)   # one dof per edge
u,v = fes.TnT()
eps = 1e-3

a = BilinearForm(curl(u)*curl(v)*dx + eps*u*v*dx).Assemble()
M = mesh.MaterialCF( { "magnet" : (1,0,0)} )
f = LinearForm(M*curl(v)*dx("magnet")).Assemble()

gfu = GridFunction(fes)

pre = Preconditioner(a, "local", block=True)
pre.Update()
from ngsolve.krylovspace import CGSolver

inv = CGSolver(a.mat, pre.mat, printrates=True, maxiter=200)
gfu.vec.data = inv * f.vec

CG iteration 1, residual = 0.4635327629266928     
CG iteration 2, residual = 0.42789858702521705     
CG iteration 3, residual = 0.3232297657475674     
CG iteration 4, residual = 0.25464403586309087     
CG iteration 5, residual = 0.18615733922752714     
CG iteration 6, residual = 0.14800181052391703     
CG iteration 7, residual = 0.11949520952077154     
CG iteration 8, residual = 0.09282291331667633     
CG iteration 9, residual = 0.07679038814089593     
CG iteration 10, residual = 0.06389531595609149     
CG iteration 11, residual = 0.05056242146491291     
CG iteration 12, residual = 0.04060539621517081     
CG iteration 13, residual = 0.03430579733377514     
CG iteration 14, residual = 0.030065910434030483     
CG iteration 15, residual = 0.025774846052082172     
CG iteration 16, residual = 0.01983965339095311     
CG iteration 17, residual = 0.013867502737470573     
CG iteration 18, residual = 0.009567470112709246     
CG iteration 19, residual = 0.006601850082293865     
CG iteration 20, residual = 0.004463903340199571     
CG iteration 21, residual = 0.002938211266571423     
CG iteration 22, residual = 0.002000033787596275     
CG iteration 23, residual = 0.0013982543263219135     
CG iteration 24, residual = 0.0009922984038274098     
CG iteration 25, residual = 0.0007815351507657772     
CG iteration 26, residual = 0.000683642014896617     
CG iteration 27, residual = 0.0006042081173435654     
CG iteration 28, residual = 0.0004933831969703863     
CG iteration 29, residual = 0.0003658814417513913     
CG iteration 30, residual = 0.00027461820027847346     
CG iteration 31, residual = 0.00021606301815433576     
CG iteration 32, residual = 0.0001852656638868365     
CG iteration 33, residual = 0.00016757613019838458     
CG iteration 34, residual = 0.00015014932580892883     
CG iteration 35, residual = 0.00012635246751698076     
CG iteration 36, residual = 0.00010095216118821378     
CG iteration 37, residual = 7.677150772435917e-05     
CG iteration 38, residual = 5.705147728680346e-05     
CG iteration 39, residual = 4.400991083489937e-05     
CG iteration 40, residual = 3.6417376943526984e-05     
CG iteration 41, residual = 3.212176823477045e-05     
CG iteration 42, residual = 2.7186175188828105e-05     
CG iteration 43, residual = 2.323767518440678e-05     
CG iteration 44, residual = 2.069312434365388e-05     
CG iteration 45, residual = 1.911627307675492e-05     
CG iteration 46, residual = 1.8718180412834895e-05     
CG iteration 47, residual = 1.842446272464718e-05     
CG iteration 48, residual = 1.7710075176906277e-05     
CG iteration 49, residual = 1.649936644773816e-05     
CG iteration 50, residual = 1.5042919233686429e-05     
CG iteration 51, residual = 1.4060059508048765e-05     
CG iteration 52, residual = 1.3294831066099182e-05     
CG iteration 53, residual = 1.2250262778494888e-05     
CG iteration 54, residual = 1.0757704058947594e-05     
CG iteration 55, residual = 8.740608113448077e-06     
CG iteration 56, residual = 6.838580733733294e-06     
CG iteration 57, residual = 5.401844236565597e-06     
CG iteration 58, residual = 4.329876676381426e-06     
CG iteration 59, residual = 3.392817659702821e-06     
CG iteration 60, residual = 2.5330312963314107e-06     
CG iteration 61, residual = 1.8791455960178543e-06     
CG iteration 62, residual = 1.512762255199913e-06     
CG iteration 63, residual = 1.400764857435027e-06     
CG iteration 64, residual = 1.4201528772553445e-06     
CG iteration 65, residual = 1.484543077553714e-06     
CG iteration 66, residual = 1.506645750744541e-06     
CG iteration 67, residual = 1.5221995110452216e-06     
CG iteration 68, residual = 1.4801364239010767e-06     
CG iteration 69, residual = 1.4110760181683581e-06     
CG iteration 70, residual = 1.2533077908298905e-06     
CG iteration 71, residual = 1.0275036694747827e-06     
CG iteration 72, residual = 8.43762385593524e-07     
CG iteration 73, residual = 7.11756578349743e-07     
CG iteration 74, residual = 6.085538139091289e-07     
CG iteration 75, residual = 5.140103064211771e-07     
CG iteration 76, residual = 4.020347788031209e-07     
CG iteration 77, residual = 3.092098918926114e-07     
CG iteration 78, residual = 2.4553367053311626e-07     
CG iteration 79, residual = 2.009547208619791e-07     
CG iteration 80, residual = 1.655442290755333e-07     
CG iteration 81, residual = 1.3974281500533848e-07     
CG iteration 82, residual = 1.20185210357764e-07     
CG iteration 83, residual = 1.015471598348598e-07     
CG iteration 84, residual = 8.26051191101564e-08     
CG iteration 85, residual = 6.635817941556107e-08     
CG iteration 86, residual = 5.120263469848235e-08     
CG iteration 87, residual = 4.142946942564037e-08     
CG iteration 88, residual = 3.4856977656683026e-08     
CG iteration 89, residual = 2.9171532625864503e-08     
CG iteration 90, residual = 2.4620885627225223e-08     
CG iteration 91, residual = 2.1143891383548138e-08     
CG iteration 92, residual = 1.8993136308760616e-08     
CG iteration 93, residual = 1.7017039819867693e-08     
CG iteration 94, residual = 1.4333454712960677e-08     
CG iteration 95, residual = 1.1557028423569548e-08     
CG iteration 96, residual = 8.95521894798872e-09     
CG iteration 97, residual = 6.947450640234518e-09     
CG iteration 98, residual = 5.438531508912747e-09     
CG iteration 99, residual = 4.346737344250734e-09     
CG iteration 100, residual = 3.6087661839121814e-09     
CG iteration 101, residual = 3.1829807261609314e-09     
CG iteration 102, residual = 2.7843641797264262e-09     
CG iteration 103, residual = 2.349635933929682e-09     
CG iteration 104, residual = 2.0666980034078006e-09     
CG iteration 105, residual = 1.8682450158839147e-09     
CG iteration 106, residual = 1.6702589191564442e-09     
CG iteration 107, residual = 1.5608872971611618e-09     
CG iteration 108, residual = 1.4141783443506158e-09     
CG iteration 109, residual = 1.1801844070219922e-09     
CG iteration 110, residual = 9.680627093063401e-10     
CG iteration 111, residual = 8.362346931360255e-10     
CG iteration 112, residual = 7.562143501253274e-10     
CG iteration 113, residual = 6.519927604494129e-10     
CG iteration 114, residual = 5.603811278315021e-10     
CG iteration 115, residual = 5.069583126200502e-10     
CG iteration 116, residual = 4.926136120458067e-10     
CG iteration 117, residual = 5.022852060313368e-10     
CG iteration 118, residual = 4.873362458040147e-10     
CG iteration 119, residual = 4.179038940459572e-10     
CG iteration 120, residual = 3.3427726344761957e-10     
CG iteration 121, residual = 2.6056347314808085e-10     
CG iteration 122, residual = 2.0391573266438795e-10     
CG iteration 123, residual = 1.6060002408897057e-10     
CG iteration 124, residual = 1.290854582112953e-10     
CG iteration 125, residual = 1.0736797747190897e-10     
CG iteration 126, residual = 9.326212478264851e-11     
CG iteration 127, residual = 8.434050375655718e-11     
CG iteration 128, residual = 7.529219410275793e-11     
CG iteration 129, residual = 6.445632525789044e-11     
CG iteration 130, residual = 5.515875985552747e-11     
CG iteration 131, residual = 4.642173216119092e-11     
CG iteration 132, residual = 3.969780154913643e-11     
CG iteration 133, residual = 3.510174891172396e-11     
CG iteration 134, residual = 3.133217721853538e-11     
CG iteration 135, residual = 2.7576109259346328e-11     
CG iteration 136, residual = 2.3543762860934322e-11     
CG iteration 137, residual = 1.8993384818526183e-11     
CG iteration 138, residual = 1.4438616067029353e-11     
CG iteration 139, residual = 1.1284839854964043e-11     
CG iteration 140, residual = 9.820233339776046e-12     
CG iteration 141, residual = 9.007363468914282e-12     
CG iteration 142, residual = 8.633510002398147e-12     
CG iteration 143, residual = 8.440755426206628e-12     
CG iteration 144, residual = 7.829576433200149e-12     
CG iteration 145, residual = 6.897243155937746e-12     
CG iteration 146, residual = 5.9525636053155235e-12     
CG iteration 147, residual = 5.264312876504194e-12     
CG iteration 148, residual = 4.8402817281136264e-12     
CG iteration 149, residual = 4.492327674466352e-12     
CG iteration 150, residual = 4.08848867205207e-12     
CG iteration 151, residual = 3.5505614653778326e-12     
CG iteration 152, residual = 2.8329546473066384e-12     
CG iteration 153, residual = 2.175796464246952e-12     
CG iteration 154, residual = 1.663515736510115e-12     
CG iteration 155, residual = 1.2634537302230126e-12     
CG iteration 156, residual = 9.361744949075666e-13     
CG iteration 157, residual = 6.973280254961185e-13     
CG iteration 158, residual = 5.230426840079289e-13     
CG iteration 159, residual = 4.1134450354834574e-13     

Draw (curl(gfu), mesh, clipping=clipping, \
      vectors={"grid_size" : 40 }, draw_surf=False);

from ngsolve.la import EigenValues_Preconditioner
lam = EigenValues_Preconditioner(a.mat, pre.mat)
print ("lammin, lammax=", lam[0], lam[-1])
print ("kappa=", lam[-1]/lam[0])

lammin, lammax= 0.02952541052556933 5.916054099105543
kappa= 200.37161190298016

We observe that the condition number of the block-Jacobi preconditioner

• is robust in \(\varepsilon\)

• depends on the mesh-size

We observe the the condition number of the Jacobi preconditioner

• depends on \(\varepsilon\)

• depends on the mesh-size

76.5.2. Robust coarse-grid correction

If the penalty term involves a projection, then the coarse-grid null-space \(V_{H,0}\) is not a sub-space of the fine-grid null-space \(V_{h,0}\). An example is Stokes with penalty:

\[ \| u_H \|_{A_H^\varepsilon}^2 = \| u_H \|_A^2 + \frac{1}{\varepsilon} \| R_H \operatorname{div} u_H \|_{L_2}^2 \]

The operator-norm of the prolongation is

\[ \| P \|_{(V_H, A_H^\varepsilon) \rightarrow (V_h, A_h^\varepsilon)} = \sup_{v_H} \frac{ \| P v_H \|_{A_h^\varepsilon}}{\| v_H \|_{A_H^\varepsilon}} \]

To be robust in the parameter, the prolongation operator \(P : V_H \rightarrow V_h\) must map the coarse-grid null-space into the fine-grid null-space:

\[ P V_{H,0} \subset V_{h,0} \]

If the null-spaces are nested (as for Maxwell equations), the prolongation is canonical embedding.

For the ASM - analysis we need the existence of stable interpolation operators from fine to coarse: \(I_H : V_h \rightarrow V_H\). It must be compatible with null-spaces, i.e. \(I_H V_{h,0} \subset V_{H,0}\).

76.5.3. Two-level analysis for Maxwell equations

Let \(V_h\) be a Nedelec finite element sub-space of \(H(\operatorname{curl})\), and consider the bilinear-form

\[ A^\varepsilon(u,v) = \int \operatorname{curl} u \operatorname{curl} v + \varepsilon u v \]

We prove that a two-level method with an AFW block-Jacobi smoother provides an optimal preconditioner.

Let \(W_H \subset W_h \subset H^1\) be compatible finite element spaces such that \(\nabla W_h = V_{h,0} = \{ v \in V_h : \operatorname{curl} v_h = 0 \}\), and \(\nabla W_H = W_{H,0}\).

We need a few technical tools:

Theorem: (Regular Decomposition): There holds

\[ H(\operatorname{curl}) = [H^1]^3 + \nabla H^1 \]

For every function \(u\) from \(H(\operatorname{curl})\) there exists a decomposition

\[ u = z + \nabla \phi \]

with \(z \in [H^1]^3\) and \(\phi \in H^1\), and which is stable in the sense

\[\begin{eqnarray*} \| \phi \|_{H^1} & \preceq & \| u \|_{H(\operatorname{curl})} \\ \| z \|_{H^1} & \preceq & \| \operatorname{curl} u \|_{L_2} \\ \end{eqnarray*}\]

On the other hand, every function from \([H^1]^3\), and every gradient from \(H^1\) is in \(H(\operatorname{curl})\).

Proof: Pasciak, J. E.; Zhao, J.: Overlapping Schwarz methods in H(curl) on polyhedral domains. J. Numer. Math. 10 (2002)

Theorem: (Commuting quasi-interpolation operators) There exist projections \(\pi_{V_h} : H(\operatorname{curl}) \rightarrow V_h \) and \(\pi_{W_h} : H^1 \rightarrow W_h\) which commute with the gradient:

\[ \pi_{V_h} \nabla = \nabla \pi_{W_h} \]

They are stable in norms and semi-norms, and provide optimal approximation error estimates.

Proof: Schöberl, A multilevel decomposition result in H(curl). in Multigrid, Multilevel and Multiscale Methods, EMG 2005

We are ready to prove that the condition-number of the two-level preconditioner is robust in the mesh-size \(h\), as well as the regularization parameter \(\varepsilon\). We assume that \(H = ch\). We transfer results from the \(H^1\)-case to the \(H(\operatorname{curl})\) case.

In the \(H^1\)-case, we have shown that the coarse grid interpolation error satisfies

\[ \| u_h - \pi_{W_H} u_h \|_{L_2} \preceq h \| \nabla u_h \| \]

The corresponding lemma for the \(H(\operatorname{curl})\) case is:

Lemma: For the coarse-grid interpolation error there exists a decomposition

\[ u_h - \pi_{V_H} u_h = r + \nabla \psi, \]

where

\[ \| r \|_{L_2}^2 \preceq h^2 \| \operatorname{curl} u \|^2 \]

and

\[ \varepsilon \| \psi \|_{L_2}^2 \preceq h^2 \| u \|_{A^\varepsilon}^2 \]

Proof: Let \(u_h \in V_h \subset H(\operatorname{curl})\). There exists a regular decomposition

\[ u_h = z + \nabla \Phi, \]

where

\[ \| z \|_{H^1}^2 \preceq \| \operatorname{curl} u \|_{L_2}^2 \]

and

\[ \varepsilon \| \Phi \|_{H^1}^2 \preceq \| u \|_{A^\varepsilon}^2 \]

By the commutativity of the interpolation operator we obtain

\[\begin{eqnarray*} u_h - \pi_{V_H} u_h & = & z - \pi_{V_H} z + \nabla \Phi - \pi_{V_H} \nabla \Phi \\ & = & z - \pi_{V_H} z + \nabla ( \Phi - \pi_{W_H} \Phi ) \end{eqnarray*}\]

By setting

\[ r = z - \pi_{V_H} z \qquad \text{and} \qquad \psi = \Phi - \pi_{W_H} \Phi \]

and using the approximation properties of the interpolation operators we obtain the splitting.

Lemma: For the coarse-grid interpolation error there exists a discrete decomposition

\[ u_h - \pi_{V_H} u_h = r_h + \nabla \psi_h. \]

Proof: Apply the commuting projection operator \(\pi_{V_h}\) to the previous lemma to obtain

(76.1)\[\begin{eqnarray} u_h - \pi_{V_H} u_h & = & \pi_{V_h} ( u_h - \pi_{V_H} u_h ) \\ & = & \pi_{V_h} r + \pi_{V_h} \nabla \psi \\ & = & \pi_{V_h} r + \nabla \pi_{W_h} \psi \end{eqnarray}\]

The discrete decomposition is now \(r_h = \pi_{V_h} r\) and \(\psi_h = \pi_{W_h} \psi\).

In the \(H^1\)-case, we can apply the ASM lemma to analyze the Jacobi preconditioner \(D\) as follows:

\[ \| u \|_D^2 = \inf_{u = \sum u_i} \sum \| u_i \|_{H^1}^2 \preceq \inf_{u = \sum u_i} h^{-2} \sum \| u_i \|_{L_2}^2 \preceq h^{-2} \| u \|_{L_2}^2 \]

Lemma: For the AFW block-Jacobi preconditioner \(D\) there holds

\[ \| u_h \|_D^2 \preceq \inf_{u_h = r_h + \nabla \psi_h} h^{-2} \| r_h \|_{L_2}^2 + \varepsilon h^{-2} \| \psi_h \|_{L_2}^2 \]

Proof: Let \(u_h = r_h + \nabla \psi_h\) be an arbitrary decomposition. Then the triangle inequality implies

\[ \| u_h \|_D \leq \| r_h \|_D + \| \nabla \psi_h \|_D \]

Now we apply the ASM lemma for both terms:

\[\begin{eqnarray*} \| r_h \|_D^2 & = & \inf_{r = \sum r_i} \sum \| r_i \|_{A^\varepsilon}^2 \\ & = & \inf_{r = \sum r_i} \sum \| \operatorname{curl} r_i \|_{L_2}^2 + \varepsilon \| r_i \|_{L_2}^2 \\ & \preceq & \inf_{r = \sum r_i} h^{-2} \| r_i \|_{L_2}^2 \\ & \approx & h^{-2} \| r \|_{L_2}^2 \end{eqnarray*}\]

For the decomposition of \(\nabla \psi\) we use that we can decompose \(\psi = \sum \psi_i\) with \(\psi_i \in W_i\), where \(\nabla W_i \subset V_i\):

\[\begin{eqnarray*} \varepsilon \| \nabla \psi_h \|_D^2 & = & \inf_{\nabla \psi_h = \sum v_i} \sum \| v_i \|_{A^\varepsilon}^2 \\ & \leq & \inf_{\psi_h = \sum \psi_i} \sum \| \nabla \psi_i \|_{A^\varepsilon}^2 \\ & = & \inf_{\psi_h = \sum \psi_i} \sum \varepsilon \| \nabla \psi_i \|_{L_2}^2 \\ & \preceq & \inf_{\psi_h = \sum \psi_i} \sum \varepsilon h^{-2} \| \psi_i \|_{L_2}^2 \\ & \approx & \varepsilon h^{-2} \| \psi_h \|_{L_2}^2 \end{eqnarray*}\]

Since \(\| u \|_D\) is bounded by any such decomposition of \(u_h\), it holds also for the infimum, and the lemma is proven.

Theorem The two-level preconditioner with AFW block-smoothers provide a optimal condition number \(\kappa = O(1)\).

Proof: We have to verify that

\[ \inf_{u_h = u_H + \sum u_i} \left\{ \| u_H \|_{A^\varepsilon}^2 + \sum \| u_i \|_{A^\varepsilon}^2 \right\} \preceq \| u_h \|_{A^\varepsilon}^2 \]

We choose \(u_H := \pi_{V_H} u_h\). Continuity in energy-norm of the projector bounds the first term. The interpolation rest provides a stable splitting

\[ u_f = u_h - u_H = r_h + \nabla \psi_h, \]

where the components \(r_h\) and \(\psi_h\) are bounded in \(L_2\)-norms. This \(u_f\) can be further split into local functions as \(u_f = \sum u_i\).

from ngsolve import *
from netgen.csg import *
from ngsolve.webgui import Draw

geo = CSGeometry()
magnet = Cylinder(Pnt(-1,0,0), Pnt(1,0,0), 0.3) \
    * OrthoBrick(Pnt(-1,-1,-1), Pnt(1,1,1))
box = OrthoBrick(Pnt(-3, -3, -3), Pnt(3,3,3))
geo.Add (magnet.mat("magnet"))
geo.Add ((box-magnet).mat("air"))
mesh = Mesh(geo.GenerateMesh(maxh=1))

fes = HCurl(mesh, order=0, autoupdate=True)
print ("ndof =", fes.ndof)

u,v = fes.TnT()
eps = 1e-3

a = BilinearForm(curl(u)*curl(v)*dx + eps*u*v*dx)
pre = Preconditioner(a, "multigrid", smoother="block")

M = mesh.MaterialCF( { "magnet" : (1,0,0)} )
f = LinearForm(M*curl(v)*dx("magnet"))

with TaskManager():
    a.Assemble()
    f.Assemble()
    for l in range(1):  # try also 2
        mesh.Refine()
        print ("ndof =", fes.ndof)
        a.Assemble()
        f.Assemble()

gfu = GridFunction(fes)

ndof = 8800

ndof = 77372

from ngsolve.krylovspace import CGSolver

inv = CGSolver(a.mat, pre.mat, printrates=True, maxiter=200)
with TaskManager():
    gfu.vec.data = inv * f.vec

CG iteration 1, residual = 0.6907017058419193     

CG iteration 2, residual = 0.016142253725977114     
CG iteration 3, residual = 0.00046060801777727964     
CG iteration 4, residual = 1.0543115072670755e-05     

CG iteration 5, residual = 3.8277013571751647e-07     

CG iteration 6, residual = 3.371148838672236e-08     
CG iteration 7, residual = 2.477239892603522e-09     
CG iteration 8, residual = 1.0905496474263337e-10     

CG iteration 9, residual = 1.1153690820937868e-11     

CG iteration 10, residual = 5.439845270680144e-13     

clipping = { "pnt" : (0,0,0), "function" : False}
Draw (curl(gfu), mesh, order=1, clipping=clipping, \
      vectors={"grid_size" : 40 }, draw_surf=False);