PETSc MatSetBlockSize For Product Function Spaces

When using a product space such as,

order = 1
Vv = VectorFunctionSpace(mesh, 'Lagrange', order)
Vs = FunctionSpace(mesh, 'Lagrange', order)
V = Vv*Vs

the matrix assemble routine does not appear to set the matrix block size appropriately. See the following exchange on the FEniCS Q&A: http://fenicsproject.org/qa/6535/setting-the-number-of-unknowns-for-hypre

I will summarize the key points here. In short,

mesh = UnitSquareMesh(4,4)
V = VectorFunctionSpace(mesh, "CG", 1)
u,v = TrialFunction(V), TestFunction(V)
A = assemble(dot(u,v)*dx)
Aback = as_backend_type(A).mat()
print Aback.block_size

returns the correct block_size, 2. But if we assemble a system using product spaces, such as:

mesh = UnitSquareMesh(4, 4)

order = 1
Vv = VectorFunctionSpace(mesh, 'Lagrange', order)
Vs = FunctionSpace(mesh, 'Lagrange', order)
V = Vv*Vs

def LSop(u, p):
    L1 = u - grad(p)
    L2 = div(u)
    L3 = curl(u)
    return [L1, L2, L3]

(u, p) = TrialFunctions(V)
(v, q) = TestFunctions(V)

Lu = LSop(u, p)
Lv = LSop(v, q)

a = dot(Lu[0], Lv[0])*dx + Lu[1]*Lv[1]*dx + dot(Lu[2], Lv[2])*dx
A = assemble(a)

Aback = as_backend_type(A).mat()
print Aback.block_size

The block_size parameter returns 1, and significantly hinders the performance of the Hypre preconditioner. If, instead, I use a "flat" VectorFunctionSpace that contains all the unknowns and decompose them manually, the correct block_size is again recovered:

mesh = UnitSquareMesh(4, 4)

order = 1

Vsp = VectorFunctionSpace(mesh, 'Lagrange', order, 3)

UU = TrialFunctions(Vsp)
VV = TestFunctions(Vsp)

u = as_vector((UU[0], UU[1]))
p = UU[2]
v = as_vector((VV[0], VV[1]))
q = VV[2]

def LSop(u, p):
    L1 = u - grad(p)
    L2 = div(u)
    L3 = curl(u)
    return [L1, L2, L3]

(u, p) = TrialFunctions(V)
(v, q) = TestFunctions(V)

Lu = LSop(u, p)
Lv = LSop(v, q)

a = dot(Lu[0], Lv[0])*dx + Lu[1]*Lv[1]*dx + dot(Lu[2], Lv[2])*dx
A = assemble(a)

Aback = as_backend_type(A).mat()
print Aback.block_size

which returns the expected size of 3, and allows for Hypre preconditioner to perform as expected. We are okay using this second formulation, but would much prefer to be able to use the more convenient product space formulation to organize our systems.

This may be difficult (or impossible) if the product space contains spaces with different orders or element types, as the unknowns can't be ordered node-wise due to a different number of DOFs per triangular element in each unknown. I propose that if using a product space where each of the base spaces is of the same order and element type, an appropriate call should be made to MatSetBlockSize(...) in order to set block_size such that Hypre will be automatically called with HYPRE_BoomerAMGSetNumFunctions(<num_unknowns>) from PETSc.

Comments (2)