Structure of Saddle-point Problems

92. Structure of Saddle-point Problems#

92.1. Constrained minimization problems#

Let $A \in {\mathbf R}^{n\times n}$ be spd, and $f \in {\mathbf R}^n$. Then the minimization problem

\[ \min_u \frac{1}{2} \, u^T A u - f^T u \]

is equivalent to the linear system

\[ A u = f. \]

Since $A$ is spd, the function to be minimized is convex.

Now let $B \in {\mathbf R}^{m\times n}$ of full rank, with $m \leq n$ and, and $g \in {\mathbf R}^m$. We consider the constrained minimization problem

\[ \min_{u \atop B u = g} \frac{1}{2} \, u^T A u - f^T u. \]

The minimizers are found among the critical points of the Lagrange function

\[ L(u, p) = \frac{1}{2} \, u^T A u - f^T u + p^T (B u - g), \]

where $p \in {\mathbb R}^m$ is the Lagrange parameter. Differentiation with respect to $u$ and $p$, and zeroing leads to the Karush-Kuhn-Tucker (KKT) equations

\[\begin{split} \begin{array}{ccccl} A u & + & B ^T p & = & f, \\ B u & & & = & g. \end{array} \end{split}\]

The solution is not a minimizer, but a saddle-point of the Lagrange function $L(.,.)$.

Theorem: The Hessian matrix of $L(.,.)$,

\[\begin{split} \left( \begin{array}{cc} A & B \\ B & \end{array} \right) \in {\mathbf R}^{(m+n) \times (m+n)} \end{split}\]

has $n$ positive and $m$ negative eigenvalues.

Proof: We define the sub-spaces

\[\begin{eqnarray*} V_1 & := & \{ (u, 0) : u \in {\mathbf R}^n \}, \\ V_2 & := & \{ (-A^{-1} B^T p, p) : p \in {\mathbf R}^m \}, \end{eqnarray*}\]

and evaluate for $u \neq 0$

\[\begin{split} \left( \begin{array}{c} u \\ 0 \end{array} \right)^T \left( \begin{array}{cc} A & B \\ B & \end{array} \right) \left( \begin{array}{c} u \\ 0 \end{array} \right) = u^T A u > 0 \end{split}\]

and for $p \neq 0$

\[\begin{split} \left( \begin{array}{c} -A^{-1}B^T p \\ p \end{array} \right)^T \left( \begin{array}{cc} A & B \\ B & \end{array} \right) \left( \begin{array}{c} -A^{-1}B^T p \\ p \end{array} \right) = - p^T B A^{-1} B^T p < 0. \end{split}\]

The matrix $S := B A^{-1} B^T$ is called the (negative) Schur complement. For an spd matrix $A$ and for the rank-$m$ matrix $B$ it is positive definit.

Since we found an $n$-dimensional sub-space where the quadratic form is positive, and an $m$-dimensional sub-space where it is negative, the $\min-\max$-theorem proves that we have $n$ positive and $m$ negative eigenvalues.

92.2. Examples#

Stokes Equation

The vector-field $u \in [H^1_{0,D}]^2$ is the velocity of an incompressible flow. It minimizes the $H^1$-norm, under the incompressibility constraint $\operatorname{div} u = 0$:

\[ \min_{u \atop \operatorname{div} u} \frac{1}{2} \int | \nabla u |^2 - \int f u \]

The variational formulation of the KKT system is

\[\begin{split} \begin{array}{ccccll} \int \nabla u : \nabla v & + & \int \operatorname{div} v \, p & = & \int f v & \qquad \forall \, v \\ \int \operatorname{div} u \, q & & & = & 0 & \qquad \forall \, q \end{array} \end{split}\]

Discretization leads to the matrix equation

\[\begin{split} \begin{array}{ccccl} A u & + & B ^T p & = & f, \\ B u & & & = & 0, \end{array} \end{split}\]

where $A$ is the discretization of the $-\Delta$-operator for both components of the velocity, and $B$ is the discretization of the $\operatorname{div}$-operator. It’s transposed matrix $B^T$ is the discretization of the negative $\nabla$-operator. The physical meaning of the Lagrange parameters $p$ is pressure.

Dirichlet boundary conditions

We can interpret a Dirichlet problem as a constrained minimization problem

\[ \min_{u \in H^1 \atop u = u_D \text{ on } \Gamma_D} \frac{1}{2} \int |\nabla u |^2 - \int f u \]

The KKT system is to find $u \in H^1(\Omega)$ and $\lambda$ on the boundary such that

\[\begin{split} \begin{array}{ccccll} \int_\Omega \nabla u \nabla v & + & \int_{\Gamma_D} v \lambda & = & \int f v & \qquad \forall \, v \\ \int_{\Gamma_D} u \, \mu & & & = & \int_{\Gamma_D} u_D \, \mu & \qquad \forall \, \mu \end{array} \end{split}\]

By integration by parts in the first equation, we obtain

\[ -\int_\Omega \Delta u v + \int_{\partial \Omega} \frac{\partial u}{\partial n} v + \int_{\Gamma_D} \lambda v = \int_\Omega f v. \]

We see that the meaning of the Lagrange parameter $\lambda$ is the negative normal derivative $- \tfrac{\partial u}{\partial n}$

92.3. Schur complement iteration#

We solve the block-system

\[\begin{split} \begin{array}{ccccl} A u & + & B ^T p & = & f, \\ B u & & & = & g. \end{array} \end{split}\]

Since $A$ is spd and thus regular, we can solve for $u$ in the first block-row

\[ u = A^{-1} (f - B^T p). \]

Inserting this $u$ into the second block-row gives

\[ B A^{-1} (f - B^T p) = g, \]

or equivalently

\[ B A^{-1} B^T p = B A^{-1} f - g. \]

The Schur complement matrix $S = B A^{-1} B^T$ is positive definite. The Uzawa algorithm (and there are many variants of it) is now to use a Richardson iteration, or conjugate gradients, to solve the Schur complement equation

\[ S p = \tilde g \]

with $\tilde g = B A^{-1} f - g$. After computing the Lagrange parameter $p$, the first component $u$ is found from $u = A^{-1} (f - B^T p)$.

Of course, the drawback of this method is that it requires the inverse $A^{-1}$, which is usually not cheaply available. If the Uzawa iteration does not take many iteration steps, one can also use an (inner) iterative solver to apply $A^{-1}$ in every step.

To design preconditioners for the Uzawa iteration, it is important to understand the norm induced by the Schur-complement $S$:

\[ \| p \|_S^2 = \| B^T p \|_{A^{-1}}^2 = \left( \sup_u \frac{ u^T B^T p}{ \| u \|_{A}} \right)^2 \]

If $A$ and $B$ come from finite-element discretization of bilinear-forms $a(.,.)$ and $b(.,.)$, then the last term can be recast as

\[ \sup_{u \in V_h} \frac{b(u,p)^2}{a(u,u)} \]

If we can prove continuity of $b(.,.)$

\[ b(u,p) \leq \overline \beta \, \| u \|_A \| p \|_Q, \]

and the famous LBB-condition

\[ \sup_u \frac{b(u,p)}{\| u \|_A} \geq \underline \beta \, \| p \|_Q, \]

we immediately get

\[ \underline \beta^2 \, \| p \|_Q^2 \leq \| p \|_S^2 \leq \overline \beta^2 \, \| p \|_Q^2 \]

If we have a preconditioner for the norm-matrix of the $Q$-norm, we have a preconditioner for the Schur complement $S$.

from ngsolve import *
from ngsolve.webgui import Draw
mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))

fesu = H1(mesh, order=1)
feslam = H1(mesh, order=1, definedon=mesh.Boundaries(".*"))
u,v = fesu.TnT()
lam,mu = feslam.TnT()

a = BilinearForm(grad(u)*grad(v)*dx + u*v*dx).Assemble()
b = BilinearForm(u*mu*ds).Assemble()
f = LinearForm(10*v*dx).Assemble()
g = LinearForm((x+y)*mu*ds).Assemble()

ainv = a.mat.Inverse()
gtilde = (b.mat @ ainv * f.vec - g.vec).Evaluate()

S = b.mat @ ainv @ b.mat.T
preS = Projector(feslam.FreeDofs(), True)

from ngsolve.krylovspace import CGSolver
Sinv = CGSolver(S, preS, printrates=True, maxiter=200)

gflam = GridFunction(feslam)
gfu = GridFunction(fesu)

gflam.vec.data = Sinv * gtilde 
gfu.vec.data = ainv * (f.vec - b.mat.T * gflam.vec)

Draw (gfu);

CG iteration 1, residual = 5.70374341015521     
CG iteration 2, residual = 0.3225362849898911     
CG iteration 3, residual = 0.08339010831277839     
CG iteration 4, residual = 0.01117390734935611     
CG iteration 5, residual = 0.0047927265813528675     
CG iteration 6, residual = 0.004039693809759069     
CG iteration 7, residual = 0.003477548920994166     
CG iteration 8, residual = 0.004676900022939997     
CG iteration 9, residual = 0.0028443097964922547     
CG iteration 10, residual = 0.0010427358439432045     
CG iteration 11, residual = 0.0015208669764192647     
CG iteration 12, residual = 0.0010007715432088899     
CG iteration 13, residual = 0.0008083236800528286     
CG iteration 14, residual = 0.00038646050624300133     
CG iteration 15, residual = 0.00023345852649526882     
CG iteration 16, residual = 0.00019121655931019407     
CG iteration 17, residual = 0.0002507057465019636     
CG iteration 18, residual = 7.906795223588332e-05     
CG iteration 19, residual = 6.20599823174278e-05     
CG iteration 20, residual = 5.296784849054979e-05     
CG iteration 21, residual = 3.7249857424010985e-05     
CG iteration 22, residual = 2.233427661186067e-05     
CG iteration 23, residual = 2.5759653496097197e-05     
CG iteration 24, residual = 9.522437201035332e-06     
CG iteration 25, residual = 7.208353059797299e-06     
CG iteration 26, residual = 3.650418161152721e-06     
CG iteration 27, residual = 6.150077558234375e-06     
CG iteration 28, residual = 4.63861041586411e-06     
CG iteration 29, residual = 9.558709031922435e-07     
CG iteration 30, residual = 7.351143807631565e-07     
CG iteration 31, residual = 4.395662608619207e-06     
CG iteration 32, residual = 1.728559224122473e-07     
CG iteration 33, residual = 1.9409072280009905e-07     
CG iteration 34, residual = 5.965516822381169e-08     
CG iteration 35, residual = 3.763709263634211e-08     
CG iteration 36, residual = 1.98928958686514e-08     
CG iteration 37, residual = 2.6624750020974526e-08     
CG iteration 38, residual = 7.27342755478678e-08     
CG iteration 39, residual = 3.348927688529326e-09     
CG iteration 40, residual = 2.943042808256458e-09     
CG iteration 41, residual = 1.6193400996968746e-09     
CG iteration 42, residual = 1.9205838369725923e-09     
CG iteration 43, residual = 3.6576873786252613e-09     
CG iteration 44, residual = 5.15316554116041e-10     
CG iteration 45, residual = 4.938952741759179e-10     
CG iteration 46, residual = 8.122531987166569e-10     
CG iteration 47, residual = 3.306527104434466e-10     
CG iteration 48, residual = 2.4383600045391037e-10     
CG iteration 49, residual = 1.9216564883158577e-11     
CG iteration 50, residual = 4.198500568785393e-11     
CG iteration 51, residual = 8.843758072882099e-12     
CG iteration 52, residual = 1.1584849958238028e-11     
CG iteration 53, residual = 9.66978184735845e-13     

If we experiment with the mesh-size, we observe that the iteration number doubles, if we refine the mesh by a factor of 4, which numerically indicate that

\[ \kappa \{ S \} = O(h^{-1}) \]

We will come back to this, later.

92.4. Block-preconditioning#

Now we want to solve directly the saddle-point system, without forming the Schur complement explicitly. We want to use preconditioners for both fields separately.

First, we show that the spectrum of

\[\begin{split} \left( \begin{array} {cc} A & 0 \\ 0 & S \end{array} \right)^{-1} \left( \begin{array} {cc} A & B^T \\ B & 0 \end{array} \right) \end{split}\]

consists just of a few, small eigenvalues, and thus this matrix is well conditioned.

Theorem The spectrum $\sigma$ of the generalized eigenvalue problem

\[\begin{split} \left( \begin{array} {cc} A & B^T \\ B & 0 \end{array} \right) \left( \begin{array}{cc} u \\ p \end{array} \right) = \lambda \left( \begin{array} {cc} A & 0 \\ 0 & S \end{array} \right) \left( \begin{array}{cc} u \\ p \end{array} \right) \end{split}\]

satisfies $\sigma \subset \{ 1, (1+\sqrt{5})/2, (1-\sqrt{5})/2 \}$.

Proof: The eigenvalue problem reads as

\[\begin{split} \left( \begin{array} {ccc} A u & + & B^T p \\ B u & & \end{array} \right) = \lambda \left( \begin{array} {c} A u \\ S p \end{array} \right), \end{split}\]

with $(u,p) \neq 0$. The first row reads as

\[ A u + B^T p = \lambda A u \]

or $$ B^T p = (\lambda-1) A u. $$ This implies $$ B A^{-1} B^T p = (\lambda-1) B u $$ or $$ S p = (\lambda-1) B u $$

Inserting the second equation $Bu = \lambda S p$, we obtain

\[ S p = (\lambda-1) \lambda S p \]

If $p$ is non-zero, the eigenvalue $\lambda$ must satisfy the quadratic equation

\[ 1 = (\lambda - 1) \lambda, \]

i.e. $\lambda = \frac{1}{2} (1 \pm \sqrt{5})$. If $p = 0$, but $u \neq 0$, the equation $B^T p = (\lambda-1) A u$ can be satisfied only if $\lambda = 1$.

Since

\[\begin{split} \left(\begin{array}{cc} A&B^T \\ B&0\end{array}\right) \end{split}\]

is symmetric, and

\[\begin{split} \left(\begin{array}{cc} A&0 \\ 0&S\end{array}\right) \end{split}\]

provides a norm, the fraction $$ \frac{\max | \sigma |}{\min | \sigma |} = \frac{1+\sqrt{5}}{| 1-\sqrt{5} |} = 2.618... $$ bounds the norm of the preconditioned matrix. It also bounds the number of iterations of a Krylov-space method. However, since eigenvalues are positive and negative, we cannot use conjugate gradients. One can apply MinRes, or GMRes.

from ngsolve import *
from ngsolve.webgui import Draw

mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))

fesu = H1(mesh, order=1, autoupdate=True)
feslam = H1(mesh, order=1, definedon=mesh.Boundaries(".*"), autoupdate=True)
u,v = fesu.TnT()
lam,mu = feslam.TnT()

a = BilinearForm(grad(u)*grad(v)*dx + u*v*dx)
b = BilinearForm(u*mu*ds)
pre = Preconditioner(a, type="multigrid")
a.Assemble()
b.Assemble()

for i in range(2):
    mesh.Refine()
    a.Assemble()
    b.Assemble()

f = LinearForm(10*v*dx).Assemble()
g = LinearForm((x+y)*mu*ds).Assemble()

K = BlockMatrix( [[a.mat, b.mat.T], [b.mat, None]] )
F = BlockVector( [f.vec, g.vec])
U = K.CreateRowVector()
prelam = Projector(feslam.FreeDofs(), True)
PRE = BlockMatrix( [[pre.mat, None], [None, prelam]] )

from ngsolve.krylovspace import MinRes

U.data = MinRes(K, pre=PRE, rhs=F, maxsteps=400)

LinearSolver iteration 1, residual = 10.006045236381954     
LinearSolver iteration 2, residual = 2.6923699593826087     
LinearSolver iteration 3, residual = 1.5965294423392078     
LinearSolver iteration 4, residual = 0.1792227859371181     
LinearSolver iteration 5, residual = 0.165531846610526     
LinearSolver iteration 6, residual = 0.1288692895408873     
LinearSolver iteration 7, residual = 0.04192901526139578     
LinearSolver iteration 8, residual = 0.0409233253352542     
LinearSolver iteration 9, residual = 0.03804945087632947     
LinearSolver iteration 10, residual = 0.007618029409365915     
LinearSolver iteration 11, residual = 0.005474471005680294     
LinearSolver iteration 12, residual = 0.005471921389164381     
LinearSolver iteration 13, residual = 0.004824883063804512     
LinearSolver iteration 14, residual = 0.0019283965718214849     
LinearSolver iteration 15, residual = 0.0018171867877555712     
LinearSolver iteration 16, residual = 0.0018122634640930398     
LinearSolver iteration 17, residual = 0.0013064909694558882     
LinearSolver iteration 18, residual = 0.0011031206144912565     
LinearSolver iteration 19, residual = 0.0011026808151518698     
LinearSolver iteration 20, residual = 0.0008766395383136351     
LinearSolver iteration 21, residual = 0.0007898291412045709     
LinearSolver iteration 22, residual = 0.0007895774460520821     
LinearSolver iteration 23, residual = 0.0006500874349254983     
LinearSolver iteration 24, residual = 0.0005612430886061916     
LinearSolver iteration 25, residual = 0.0005610940596745528     
LinearSolver iteration 26, residual = 0.0004940141768467384     
LinearSolver iteration 27, residual = 0.00047399237453467195     
LinearSolver iteration 28, residual = 0.0004724799260642711     
LinearSolver iteration 29, residual = 0.0003364149752795049     
LinearSolver iteration 30, residual = 0.0003002716303972033     
LinearSolver iteration 31, residual = 0.00029653895051867045     
LinearSolver iteration 32, residual = 0.00028592842807910197     
LinearSolver iteration 33, residual = 0.0002859281127930502     
LinearSolver iteration 34, residual = 0.0002573518232922922     
LinearSolver iteration 35, residual = 0.0001663262957824955     
LinearSolver iteration 36, residual = 0.00016535339144486697     
LinearSolver iteration 37, residual = 0.0001649625882264077     
LinearSolver iteration 38, residual = 0.00012967853147466594     
LinearSolver iteration 39, residual = 0.00010578988303286085     
LinearSolver iteration 40, residual = 0.00010540545205318694     
LinearSolver iteration 41, residual = 0.00010522759094698594     
LinearSolver iteration 42, residual = 8.449756676256749e-05     
LinearSolver iteration 43, residual = 7.218279182715383e-05     
LinearSolver iteration 44, residual = 7.212874971234854e-05     
LinearSolver iteration 45, residual = 7.164508232782922e-05     
LinearSolver iteration 46, residual = 5.951829264189257e-05     
LinearSolver iteration 47, residual = 5.812467620630036e-05     
LinearSolver iteration 48, residual = 5.8099204887225615e-05     
LinearSolver iteration 49, residual = 5.2564680785120206e-05     
LinearSolver iteration 50, residual = 4.827573870951947e-05     
LinearSolver iteration 51, residual = 4.8263555139361436e-05     
LinearSolver iteration 52, residual = 4.738533691525023e-05     
LinearSolver iteration 53, residual = 4.364952832136796e-05     
LinearSolver iteration 54, residual = 4.3606994302290244e-05     
LinearSolver iteration 55, residual = 4.3353504172289444e-05     
LinearSolver iteration 56, residual = 4.2991659399057094e-05     
LinearSolver iteration 57, residual = 4.2275383542573695e-05     
LinearSolver iteration 58, residual = 3.4519934307075106e-05     
LinearSolver iteration 59, residual = 3.425085125135383e-05     
LinearSolver iteration 60, residual = 3.422906040051716e-05     
LinearSolver iteration 61, residual = 3.1073849695662394e-05     
LinearSolver iteration 62, residual = 2.970720481660876e-05     
LinearSolver iteration 63, residual = 2.970716286972366e-05     
LinearSolver iteration 64, residual = 2.865277986967327e-05     
LinearSolver iteration 65, residual = 2.6888064585467143e-05     
LinearSolver iteration 66, residual = 2.688392177842838e-05     
LinearSolver iteration 67, residual = 2.6125103340677654e-05     
LinearSolver iteration 68, residual = 2.3221242888964587e-05     
LinearSolver iteration 69, residual = 2.3205611741745814e-05     
LinearSolver iteration 70, residual = 2.2490258574118317e-05     
LinearSolver iteration 71, residual = 2.210309025415418e-05     
LinearSolver iteration 72, residual = 2.2102972968135286e-05     
LinearSolver iteration 73, residual = 2.1263318781088735e-05     
LinearSolver iteration 74, residual = 1.7979853058349303e-05     
LinearSolver iteration 75, residual = 1.7930695236106596e-05     
LinearSolver iteration 76, residual = 1.7920411676795458e-05     
LinearSolver iteration 77, residual = 1.7153156518584036e-05     
LinearSolver iteration 78, residual = 1.6630502840545385e-05     
LinearSolver iteration 79, residual = 1.6629786499658186e-05     
LinearSolver iteration 80, residual = 1.639852908928904e-05     
LinearSolver iteration 81, residual = 1.6148025231338725e-05     
LinearSolver iteration 82, residual = 1.6026599527235916e-05     
LinearSolver iteration 83, residual = 1.575553706216768e-05     
LinearSolver iteration 84, residual = 1.5735726883430885e-05     
LinearSolver iteration 85, residual = 1.3687309393695147e-05     
LinearSolver iteration 86, residual = 1.3177274508807713e-05     
LinearSolver iteration 87, residual = 1.3177015434573429e-05     
LinearSolver iteration 88, residual = 1.2992148726230931e-05     
LinearSolver iteration 89, residual = 1.1934563299011639e-05     
LinearSolver iteration 90, residual = 1.1929186016552033e-05     
LinearSolver iteration 91, residual = 1.0992054989320825e-05     
LinearSolver iteration 92, residual = 1.0460447875170974e-05     
LinearSolver iteration 93, residual = 1.0459882321381914e-05     
LinearSolver iteration 94, residual = 1.0247339264672588e-05     
LinearSolver iteration 95, residual = 8.892348053345243e-06     
LinearSolver iteration 96, residual = 8.866382161878328e-06     
LinearSolver iteration 97, residual = 8.854617309016814e-06     
LinearSolver iteration 98, residual = 7.795460155435322e-06     
LinearSolver iteration 99, residual = 7.342648617939514e-06     
LinearSolver iteration 100, residual = 7.340513058100612e-06     
LinearSolver iteration 101, residual = 7.286468677777898e-06     
LinearSolver iteration 102, residual = 6.182855069827605e-06     
LinearSolver iteration 103, residual = 6.124788015211296e-06     
LinearSolver iteration 104, residual = 6.123974093789975e-06     
LinearSolver iteration 105, residual = 5.868808550906637e-06     
LinearSolver iteration 106, residual = 5.767181139583967e-06     
LinearSolver iteration 107, residual = 5.3281460872830935e-06     
LinearSolver iteration 108, residual = 5.026212578904594e-06     
LinearSolver iteration 109, residual = 5.02579998159349e-06     
LinearSolver iteration 110, residual = 4.922902758822715e-06     
LinearSolver iteration 111, residual = 4.031479316625231e-06     
LinearSolver iteration 112, residual = 4.000867552588764e-06     
LinearSolver iteration 113, residual = 4.000859648129893e-06     
LinearSolver iteration 114, residual = 3.91225720577679e-06     
LinearSolver iteration 115, residual = 3.259728634039026e-06     
LinearSolver iteration 116, residual = 3.247815353847782e-06     
LinearSolver iteration 117, residual = 3.2231676483122137e-06     
LinearSolver iteration 118, residual = 2.6579440575624015e-06     
LinearSolver iteration 119, residual = 2.602391512135414e-06     
LinearSolver iteration 120, residual = 2.602373769854909e-06     
LinearSolver iteration 121, residual = 2.5595543911213636e-06     
LinearSolver iteration 122, residual = 2.189030868636511e-06     
LinearSolver iteration 123, residual = 2.181684323444141e-06     
LinearSolver iteration 124, residual = 2.143191617394749e-06     
LinearSolver iteration 125, residual = 1.8159294730406653e-06     
LinearSolver iteration 126, residual = 1.8104776968147678e-06     
LinearSolver iteration 127, residual = 1.8079437136171905e-06     
LinearSolver iteration 128, residual = 1.7011307312035322e-06     
LinearSolver iteration 129, residual = 1.6949912529183096e-06     
LinearSolver iteration 130, residual = 1.6687746383460705e-06     
LinearSolver iteration 131, residual = 1.6644067604214153e-06     
LinearSolver iteration 132, residual = 1.622609931392713e-06     
LinearSolver iteration 133, residual = 1.562511031078445e-06     
LinearSolver iteration 134, residual = 1.562482826797393e-06     
LinearSolver iteration 135, residual = 1.5339463436135662e-06     
LinearSolver iteration 136, residual = 1.5079174967231562e-06     
LinearSolver iteration 137, residual = 1.5076700637248655e-06     
LinearSolver iteration 138, residual = 1.382422284986777e-06     
LinearSolver iteration 139, residual = 1.253845664137486e-06     
LinearSolver iteration 140, residual = 1.2535335656750144e-06     
LinearSolver iteration 141, residual = 1.2345160251500424e-06     
LinearSolver iteration 142, residual = 1.016189917332574e-06     
LinearSolver iteration 143, residual = 1.0105933525959124e-06     
LinearSolver iteration 144, residual = 1.0095249823151575e-06     
LinearSolver iteration 145, residual = 8.889765385466757e-07     

gfu = GridFunction(fesu)
gfu.vec.data = U[0]
Draw (gfu);