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)
-
-
reporter Thanks for the quick reply and for pointing us towards
DofMapBuilder
, we will have a look and see if we can sort it out. I'm happy to hear that you are planning to tackle this beast, best of luck on the project! - Log in to comment
Nice to hear that someone is happy about the dof reordering I added, rather than complaining about!
The code for detecting dof block structure is rather simple and conservative. @chris_richardson and I have a new project starting next month which will likely include improved detection of dof blocking. In the meantime, if you want to improve the detection of dof blocks, take a look at the code in
DofMapBuilder
.