Vertex patch smoothing for augmented Lagrangian formulations of the Oseen problem

In this tutorial, we will see how to use an augmented Lagrangian formulation to precondition the Oseen problem, i.e.

\[ \begin{align}\begin{aligned}\text{Given } \vec{\beta} \in \mathbb{R}^3 \text{ find } (\vec{u}, p) \in [H^1_{0}(\Omega)]^d \times L^2(\Omega) \text{ s.t. }\\\begin{split}\begin{cases} \nu (\nabla \vec{u}, \nabla \vec{v})_{L^2(\Omega)} + (\nabla \cdot \vec{v}, p)_{L^2(\Omega)} - (\nabla \vec{u} \vec{\beta}, \vec{v})_{L^2(\Omega)} + \gamma (\text{div}(\vec{u}), \text{div}(\vec{v}))_{L^2(\Omega)} = (\vec{f}, \vec{v})_{L^2(\Omega)} \quad \forall v \in H^1_{0}(\Omega) \\ (\nabla \cdot \vec{u}, q)_{L^2(\Omega)} = 0 \quad \forall q \in L^2(\Omega) \end{cases}\end{split}\end{aligned}\end{align} \]

Let us begin defining the parameters of the problem.

from ngsolve import *
from ngsPETSc import *
nu = 1e-4; gamma = 1e8
nu = Parameter(nu)
gamma = Parameter(gamma)
b = CoefficientFunction((4*(2*y-1)*(1-x)*x, -4*(2*x-1)*(1-y)*y))
uin = CoefficientFunction((1,0))

In particular, we will consider a high-order Hood-Taylor discretisation of the problem. Such a discretisation can easily be constructed using NGSolve as follows:

from netgen.occ import *
import netgen.gui
import netgen.meshing as ngm
from ngsolve.meshes import MakeStructured2DMesh
from mpi4py.MPI import COMM_WORLD

if COMM_WORLD.rank == 0:
   ngmesh = MakeStructured2DMesh(False, 40,40).ngmesh
   mesh = Mesh(ngmesh.Distribute(COMM_WORLD))
else:
   mesh = Mesh(ngm.Mesh.Receive(COMM_WORLD))
V = H1(mesh, order=4, hoprolongation=True, dirichlet=".*")
V = V**mesh.dim
Q = L2(mesh, order=3)
u,v = V.TnT(); p,q = Q.TnT()
a = BilinearForm(nu*InnerProduct(Grad(u),Grad(v))*dx)
a += InnerProduct(Grad(u)*b,v)*dx
a += gamma*div(u)*div(v)*dx
b = BilinearForm(div(u)*q*dx)
c = BilinearForm(-1e-10*p*q*dx)
f = LinearForm(V)
g = LinearForm(Q)

In Augmented Lagrangian preconditioning for P2-P0 discretisation of the Stokes problem we have seen that augmenting the Stokes linear system with a \(\gamma(div(\vec{u}),div(\vec{v}))_{L^2(\Omega)}\) we can obtain better converge result for a field-split preconditioner where we use the pressure mass matrix instead of the Schur complement.

mG = BilinearForm((1/nu+gamma)*p*q*dx)

We can now assemble the matrices and the right-hand side of the problem.

from ngsPETSc.pc import *
a.Assemble()
aCoarsePre = PETScPreconditioner(a.mat, V.FreeDofs(),
                                 solverParameters={"pc_type":"lu"})
for l in range(2):
   if l != 0: mesh.Refine()
V.Update(); Q.Update()
dofs = BitArray(V.ndof+Q.ndof); dofs[:] = True
a.Assemble(); b.Assemble(); mG.Assemble();
f.Assemble(); g.Assemble(); c.Assemble();
prol = V.Prolongation().Operator(1)
mGpre = PETScPreconditioner(mG.mat, Q.FreeDofs(),
                                    solverParameters={"pc_type":"lu"})
K = BlockMatrix( [ [a.mat, b.mat.T], [b.mat, c.mat] ] )

As discussed in Augmented Lagrangian preconditioning for P2-P0 discretisation of the Stokes problem, the hard part remains the construction of \((A+\gamma B^TM^{-1}B)^{-1}\) to precondition the (1,1) block. We will use a two-level additive Schwarz preconditioner made of an exact coarse correction and a vertex patch smoother, similar to what we presented in Vertex patch smoothing for p-multigrid preconditioners for the Poisson problem. Notice that while the smoother is very similar to the one used in Vertex patch smoothing for p-multigrid preconditioners for the Poisson problem, for the coarse correction we are here using h-multigrid and not p-multigrid.

def VertexStarPatchBlocks(mesh, fes):
   blocks = []
   freedofs = fes.FreeDofs()
   for v in mesh.vertices:
      vdofs = set(d for d in fes.GetDofNrs(v) if freedofs[d])
      for ed in mesh[v].edges:
         vdofs |= set(d for d in fes.GetDofNrs(ed) if freedofs[d])
      blocks.append(vdofs)
   return blocks

blocks = VertexStarPatchBlocks(mesh, V)
dofs = BitArray(V.ndof); dofs[:] = True
smoother = ASMPreconditioner(a.mat, dofs, blocks=blocks,
                             solverParameters={"pc_type": "asm",
                                               "sub_ksp_type": "preonly",
                                               "sub_pc_type": "lu"})
prol = V.Prolongation().Operator(1)
two_lv = prol@ aCoarsePre @ prol.T + smoother
C = BlockMatrix( [ [two_lv, None], [None, mGpre] ] )
print("-----------|Additive h-Multigird + Vertex star smoothing|-----------")
gfu = GridFunction(V, name='AdditiveVel'); gfp = GridFunction(Q, name='AdditivePres')
gfu.vec.data[:] = 0; gfp.vec.data[:] = 0
gfu.Set(uin, definedon=mesh.Boundaries("top"))
rhs = BlockVector( [f.vec, g.vec] )
sol = BlockVector( [gfu.vec, gfp.vec] )
rhs -= K * sol
dofs = BitArray(V.ndof+Q.ndof); dofs[:] = True
solver = KrylovSolver(K,dofs, p=C,
                      solverParameters={"ksp_type": "gmres",
                                        "ksp_max_it":100,
                                        "ksp_monitor_true_residual": None,
                                        "ksp_rtol": 1e-6,
                                        "pc_type": "mat"
                                        })
solver.solve(rhs, sol)
gfu0 = GridFunction(V, name="PETSc0"); gfp0 = GridFunction(Q)
gfu0.vec.data[:]= 0
gfu0.Set(uin, definedon=mesh.Boundaries("top"))
sol0 = BlockVector( [gfu0.vec, gfp0.vec] )
sol += sol0
gfu.vec.data = sol[0]
Draw(gfu)
vtk = VTKOutput(ma=mesh, coefs=[gfu],
             names = ["velocity"],
             filename="output/Oseen_{}".format(nu.Get()),
             subdivision=0)
vtk.Do()
Preconditioner performance for different values of the Reynolds number, for a fixed penalty parameter \(\gamma=10^8\)

Raynolds number

1e-2

1e-3

1e-4

Augmented Lagrangian preconditioner

2 (1.57e-5)

3 (7.44e-6)

5 (6.59e-6)