75. A Small Number of Constraints

There are many examples where we have only a small number of constraints. An example is a three-phase transformer, where the voltage is prescribed for each coil. This gives three scalar constraints to the field problem.

We consider now saddle-point problems

\[\begin{split} \left( \begin{array}{cc} A & B^T \\ B & \end{array} \right) \left( \begin{array}{c} u \\ p \end{array} \right) = \left( \begin{array}{c} f \\ g \end{array} \right), \end{split}\]

where \(B \in {\mathbb R}^{m \times n}\) with \(m\) small.

Let us reduce the problem to the Schur complement system

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

Now, the Schur complement

\[ S = B A^{-1} B^T \]

is a small matrix, which can be actually computed by solving \(m\) problems with the row vectors of \(B\).

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.geom2d import CSG2d, Circle, Rectangle
geo = CSG2d()
rect = Rectangle( pmin=(0,0), pmax=(10,1), bottom="free", top="free", left="left", right="right" )
geo.Add(rect)
mesh = Mesh(geo.GenerateMesh(maxh=0.3))
Draw(mesh);

We consider an elasticity problem on a rectangular domain, fixed at the left end, and loaded by a force at the right boundary:

fes = VectorH1(mesh, order=3, dirichlet="left")
u,v = fes.TnT()
a = BilinearForm(fes)
a += InnerProduct(Sym(Grad(u)), Sym(Grad(v)))*dx + div(u)*div(v)*dx
a.Assemble()

f = LinearForm(1e-3*v[1]*ds("right")).Assemble()
gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec
Draw (gfu, deformation=True);

Now, instead of the loading, we want to prescribe an average displacement at the right end:

\[ \frac{1}{|\Gamma_r|} \int_{\Gamma_r} u_x = U_x \qquad \text{and} \qquad \frac{1}{|\Gamma_r|} \int_{\Gamma_r} u_y = U_y \]

these constraints can be described by linear-forms

\[ b_x(u) = U_x \qquad \text{and} \qquad b_y(u) = U_y \]

Assemble the two constraint-vectors, and write into a MultiVector:

bx = LinearForm(v[0]*ds("right")).Assemble()
by = LinearForm(v[1]*ds("right")).Assemble()

B = MultiVector(bx.vec, 2)
B[0] = bx.vec
B[1] = by.vec

Compute the Schur complement:

S = InnerProduct(B, a.mat.Inverse(fes.FreeDofs()) * B)
print (S)

  6.6441 -0.0189391
 -0.0189391 2676.93

Solve for the Lagrange parameter. It has the physical meaning of a force:

Uxy = Vector( (1, 3) )
force = S.I * Uxy
print (force)

 0.150513
 0.00112175

and reconstruct the solution \(u\) from the computed Lagrange parameter:

fvec = B * force
gfu.vec.data = a.mat.Inverse(fes.FreeDofs()) * fvec
Draw (gfu, deformation=True);

This method requires to solve one linear system with right hand side \(f\), and \(m\) linerar systems for the Schur complement, one for each scalar constraint. We can use a preconditioner for \(A\), and solve these \(m+1\) systems by conjugate gradient iterations.

75.1. Projected preconditioner

Next we present a method to solve these problem at the cost of one unconstrained problem, plus \(m\) applications of the preconditioner.

It consists of two steps:

• Reduction to homogeneous constraints

• Construct a preconditioner on the kernel of \(B\)

Reduction to homogeneous constraints

We solve the problem in two steps: First we find an \(u_0\) such that \(B u_0 = g\). We could do so by solving

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

using the Schur complement. However, this still requires the solution of \(m\) problems for the Schur complement.

A cheaper way is to solve

\[\begin{split} \left( \begin{array}{cc} \hat{A} & B^T \\ B & \end{array} \right) \left( \begin{array}{c} u_0 \\ p \end{array} \right) = \left( \begin{array}{c} 0 \\ g \end{array} \right), \end{split}\]

where \(\hat A\) is a preconditioner to \(A\). This can be solved by the inexact Schur complement, where \(A^{-1}\) is replace by \(\hat A^{-1}\). We use the very cheap Jacobi preconditioner:

prea = a.mat.CreateSmoother(fes.FreeDofs())
Shat = InnerProduct(B, prea * B)
print (Shat)

 0.190039       0
       0 0.203514

force = Shat.I * Uxy

fvec = B * force
gfu0 = GridFunction(fes)
gfu0.vec.data = prea * fvec
Draw (gfu0, deformation=True);

Now we have to solve the correction problem \(u = u_0 + u_1\), where \(u_1\) solves the saddle-point problem with homogeneous constraint:

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

Let \(V_0 = \operatorname{ker} \{ B \}\) be the null-space of \(B\). Instead of solving the problem with constraints, we could solve the variational problem on the linear sub-space \(V_0\): find \(u_1 \in V_0\) such that

\[ v^T A u_1 = v^T (f - A u_0) \qquad v \in V_0 \]

We solve this problem on the sub-space \(V_0\) by a preconditioned conjugate gradient method, where the range of the projected preconditioner-action \(\hat A^{-1}\) spans the sub-space \(V_0\):

\[ Proj(\hat A^{-1}) : d \mapsto w \qquad \text{such that} \]

\[ w \in V_0 : \quad v^T \hat A w = v^T d \qquad \forall v \in V_0 \]

The projected problem for the preconditioned system can be solved by the saddle point problem

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

With the inexact Schur complement \(\hat S = B \hat A^{-1} B^T\) this is

\[ \hat S p = B \hat A^{-1} d \]

and then

\[ w = \hat A^{-1} (d - B^T p) = \hat A^{-1} (I - B^T S^{-1} B \hat A^{-1}) \]

Bmat = BaseMatrix(B)
projpre = prea @ (IdentityMatrix() - Bmat @ BaseMatrix(Shat.I) @ Bmat.T @ prea)

res = -a.mat * gfu0.vec
res = res.Evaluate()

from ngsolve.krylovspace import CGSolver
projinv = CGSolver(a.mat, projpre, printrates=False, maxiter=1000)

gfu1 = GridFunction(fes)
gfu1.vec.data = projinv * res
Draw (gfu1, deformation=True)

gfu1.vec.data += gfu0.vec
Draw (gfu1, deformation=True);

Some remarks:

• instead of applying the preconditioner \(\hat A\) twice, one could store the matrix \(\hat A^{-1} B\)

• the Lagrange parameter \(p\) can be computed by multiplying the equation \(A u + B^T p = f\) from left with \(B \hat A^{-1}\), and then solving with the inexact Schur complement

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