### Schur complement for the Stokes problem

390
views
0
9 months ago by
Dear FEniCS users and developers,

I would like to implement a Schur complement solver for the Stokes problem. I tried the approach below that is based on the BlockMatrix class. However, it does not works out. I get PETSc error code 73 in line 114.

I also think that the size of the sub-block are not correct. Their sizes are all the same which should not be the case, since the velocity and pressure degrees of freedom are not identical.

Does anyone have idea how to implement a Schur complement? Are there any methods in FEniCS or the linear backends (PETSc or Trilinos) that support the Schur complement.

Best wishes,
Sebastian

# -*- coding: utf-8 -*-
import dolfin as dlfn
#======================================================================
# setting discretization parameters
p_deg = 2 # polynomial degree
q_deg = p_deg + 2 # quadrature degree
#======================================================================
# setting linear solver parameters
solver_method = "bicgstab"
abs_tol = 1e-12 # absolute tolerance for krylov solver
rel_tol = 1e-9 # relative tolerance for krylov solver
max_iter = 100 # maxinum number of outer krylov iterations
#======================================================================
# mesh
mesh = dlfn.UnitSquareMesh(10, 10)
space_dim = mesh.geometry().dim()
n_cells = mesh.num_cells()
#======================================================================
# element formulation
c = mesh.ufl_cell()
elemV = dlfn.VectorElement("CG", c, p_deg)
elemP = dlfn.FiniteElement("CG", c, p_deg - 1)
Wh = dlfn.FunctionSpace(mesh, elemV * elemP)
n_dofs = Wh.dim()
#======================================================================
# boundary subdomains
left_bndry = dlfn.CompiledSubDomain("near(x[0], 0.0) && on_boundary")
right_bndry = dlfn.CompiledSubDomain("near(x[0], l_x) && on_boundary", l_x = 1.0)
top_bndry = dlfn.CompiledSubDomain("near(x[1], l_y) && on_boundary", l_y = 1.0)
bot_bndry = dlfn.CompiledSubDomain("near(x[1], 0.0) && on_boundary")
#======================================================================
# dirichlet boundary conditions
Vh, _ = Wh.split()
null_vector = dlfn.Constant((0., ) *  space_dim)
dirichlet_bcs = dict(v = [], p = [])
vel_bcs = dirichlet_bcs["v"]
vel_bcs.append(dlfn.DirichletBC(Vh, null_vector, left_bndry))
vel_bcs.append(dlfn.DirichletBC(Vh, null_vector, right_bndry))
vel_bcs.append(dlfn.DirichletBC(Vh, null_vector, bot_bndry))
vel_bcs.append(dlfn.DirichletBC(Vh.sub(0), dlfn.Constant(1.), top_bndry))
#======================================================================
# creating test and trial functions
(v, p) = dlfn.TrialFunctions(Wh)
(w, q) = dlfn.TestFunctions(Wh)
#======================================================================
# surface and volume element
dA = dlfn.Measure("ds", domain = mesh)
dV = dlfn.Measure("dx", domain = mesh)
#======================================================================
# bilinear form
b = - p * dlfn.div(w) * dV
c = - dlfn.div(v) * q * dV
#======================================================================
# linear form
f = dlfn.inner(dlfn.Constant((0.0, ) * space_dim), w) * dV
#======================================================================
# assembling block matrix and right-hand side components
A = dlfn.assemble(a)
B = dlfn.assemble(b)
C = dlfn.assemble(c)
F = dlfn.assemble(f)
for bc in dirichlet_bcs["v"]:
bc.apply(A)
bc.apply(B)
bc.apply(C)
bc.apply(F)
# assigning block matrix components
M = dlfn.BlockMatrix(2,2)
M[0,0] = A
M[0,1] = B
M[1,0] = C
# assigning block right-hand side components
y = dlfn.BlockVector(2)
G = dlfn.Vector()
C.init_vector(G, 0)
y[0] = F
y[1] = G
# creating solution vector
x = dlfn.BlockVector(2)
x[0] = A.init_vector(dlfn.Vector(), 1)
x[1] = B.init_vector(dlfn.Vector(), 1)
#======================================================================
class SchurComplement(dlfn.GenericLinearOperator):
"""
Schur complement of a block system.
"""
def __init__(self, system_matrix, Ainv):
assert isinstance(system_matrix, dlfn.BlockMatrix)
assert system_matrix.size(0) == 2
assert system_matrix.size(1) == 2
self._sysmatrix = system_matrix
assert hasattr(Ainv, "solve")
self._Ainv = Ainv
def mult(self, x):
"""
Simple version for D = 0.
"""
tmp01 = dlfn.Vector()
tmp02 = dlfn.Vector()
y = dlfn.Vector()
self._system_matrix[0,1].mult(x, tmp01)
self._Ainv.solve(tmp02, tmp01)
self._sysmatrix[1,0].mult(tmp02, y)
return y
# creating schur complement
invA = dlfn.LUSolver(M[0,0])
S = SchurComplement(M, invA)
# creating right hand side of schur complement
tmp, schur_rhs = (dlfn.Vector(), ) * 2
M[1,0].init_vector(schur_rhs, 1)
invA.solve(tmp, y[0])
M[1,0].mult(tmp, schur_rhs)
#======================================================================
slvr = dlfn.KrylovSolver("cg")
slvr.solve(S, x[1], schur_rhs)​
Community: FEniCS Project
Have a look at cbc.block (https://bitbucket.org/fenics-apps/cbc.block). This library probably offers the functionality you are looking for.
written 9 months ago by Jakob Maljaars
Dear Jakob,

thanks for your reply. I installed cbc.block and adapt the example stokes.py to my problem, see code below. The method apply defined in block_bc throws an error saying that the dimension do not match. From the documentation and the examples I cannot tell what is exactly wrong.

Two questions: Is cbc.block actively developed and fully tested? Will it be integrated into FEniCS at some time?

# -*- coding: utf-8 -*-
import dolfin as dlfn
import block
from block.iterative import MinRes
from block.algebraic.petsc import ML
#======================================================================
# mesh
mesh = dlfn.UnitSquareMesh(10, 10)
space_dim = mesh.geometry().dim()
n_cells = mesh.num_cells()
#======================================================================
# element formulation
p_deg = 2 # polynomial degree
Vh = dlfn.VectorFunctionSpace(mesh, "CG", p_deg)
Ph = dlfn.FunctionSpace(mesh, "CG", p_deg - 1)
#======================================================================
# boundary subdomains
left_bndry = dlfn.CompiledSubDomain("near(x[0], 0.0) && on_boundary")
right_bndry = dlfn.CompiledSubDomain("near(x[0], l_x) && on_boundary", l_x = 1.0)
top_bndry = dlfn.CompiledSubDomain("near(x[1], l_y) && on_boundary", l_y = 1.0)
bot_bndry = dlfn.CompiledSubDomain("near(x[1], 0.0) && on_boundary")
#======================================================================
# dirichlet boundary conditions
null_vector = dlfn.Constant((0., ) *  space_dim)
vel_bcs = []
vel_bcs.append(dlfn.DirichletBC(Vh, null_vector, left_bndry))
vel_bcs.append(dlfn.DirichletBC(Vh, null_vector, right_bndry))
vel_bcs.append(dlfn.DirichletBC(Vh, null_vector, bot_bndry))
vel_bcs.append(dlfn.DirichletBC(Vh.sub(0), dlfn.Constant(1.), top_bndry))
#======================================================================
# creating test and trial functions
v, w = dlfn.TrialFunction(Vh), dlfn.TestFunction(Vh)
p, q = dlfn.TrialFunction(Ph), dlfn.TestFunction(Ph)
#======================================================================
# surface and volume element
dA = dlfn.Measure("ds", domain = mesh)
dV = dlfn.Measure("dx", domain = mesh)
#======================================================================
# bilinear form
b = - p * dlfn.div(w) * dV
c = - dlfn.div(v) * q * dV
#======================================================================
# linear form
f = dlfn.inner(dlfn.Constant((0.0, ) * space_dim), w) * dV
#======================================================================
# assembling block matrix and right-hand side components
AA = block.block_assemble([[a, b], [c, 0]])
bb = block.block_assemble([f, 0])
bcs = block.block_bc([vel_bcs, []], True)
bcs.apply(AA)
bcs.apply(bb)
#======================================================================
# preconditioner
m = p * q * dV
M = dlfn.assemble(m)
A = AA[0,0]
BB = block.block_mat([[ML(A),  0], [0, ML(M)]])
#======================================================================
# solver
AAinv = MinRes(AA, precond=BB, tolerance=1e-8, show=1)
v, p = AAinv * bb​
written 9 months ago by sebastian_g
> Will it be integrated into FEniCS at some time?

It is very improbable that cbc.block will be integrated into vanilla FEniCS. Instead FEniCS is planned to get its own support for block problems.
written 9 months ago by Jan Blechta
You can approach the problem by using PETSc functionality more directly, e.g. using petsc4py. Have a look into FENaPack (https://fenapack.readthedocs.io/en/latest/) implementing a Schur complement appoximation for Navier-Stokes (PCD preconditioner).
written 9 months ago by Jan Blechta
Thank you. I had to update my FEniCS version from 2016.2 to 2017.1. The demo files in the repository now work on my machine. I somehow had to re-install cmake as well. I will try to use the functionality provided for the Stokes.
written 9 months ago by sebastian_g