PETSc fieldsplit: PETScPreconditioner_set_fieldsplit expected 4 arguments, got 3

Question

From updating an old code:

I used to use PETSc's fieldsplit for solving Stokes system à la

from dolfin import *

mesh = UnitSquareMesh(20, 20, 'crossed')

W = VectorElement('Lagrange', mesh.ufl_cell(), 2)
P = FiniteElement('Lagrange', mesh.ufl_cell(), 1)
WP = FunctionSpace(mesh, W*P)

u_bcs = DirichletBC(WP.sub(0), Constant((0.0, 0.0)), 'on_boundary')
p_bcs = DirichletBC(WP.sub(1), Constant(0.0), 'on_boundary')
bcs = [u_bcs, p_bcs]

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

# Create Stokes system
mu = 1.0
f = Constant((1.0, 1.0))
a = mu * inner(grad(u), grad(v))*dx \
        - p * div(v) * dx \
        - q * div(u) * dx
L = dot(f, v)*dx
A, b = assemble_system(a, L, bcs)

# solver
W = WP.sub(0)
P = WP.sub(1)
u_dofs = W.dofmap().dofs()
p_dofs = P.dofmap().dofs()
prec = PETScPreconditioner()
prec.set_fieldsplit([u_dofs, p_dofs], ['u', 'p'])

PETScOptions.set('pc_type', 'fieldsplit')
PETScOptions.set('pc_fieldsplit_type', 'additive')
PETScOptions.set('fieldsplit_u_pc_type', 'lu')
PETScOptions.set('fieldsplit_p_pc_type', 'jacobi')

solver = PETScKrylovSolver('gmres', prec)
solver.set_operator(A)

solver.parameters['monitor_convergence'] = verbose
solver.parameters['report'] = verbose
solver.parameters['absolute_tolerance'] = 0.0
solver.parameters['relative_tolerance'] = tol
solver.parameters['maximum_iterations'] = 500

# Solve
up = Function(WP)
solver.solve(up.vector(), b)

Now, the above gives the error

TypeError: PETScPreconditioner_set_fieldsplit expected 4 arguments, got 3

Is anyone using fieldsplit and knows what's going on?

Answer 1

I don't know PETScPreconditioner, but you can do it with PETSc parameters:

# solver
PETScOptions.set('ksp_view')
PETScOptions.set('ksp_monitor_true_residual')
PETScOptions.set('pc_type', 'fieldsplit')
PETScOptions.set('pc_fieldsplit_type', 'additive')
PETScOptions.set('pc_fieldsplit_detect_saddle_point')
PETScOptions.set('fieldsplit_0_ksp_type', 'preonly')
PETScOptions.set('fieldsplit_0_pc_type', 'lu')
PETScOptions.set('fieldsplit_1_ksp_type', 'preonly')
PETScOptions.set('fieldsplit_1_pc_type', 'jacobi')

solver = PETScKrylovSolver('gmres')
solver.set_operator(A)
solver.set_from_options()

solver.parameters['absolute_tolerance'] = 0.0
solver.parameters['relative_tolerance'] = 1.e-8
solver.parameters['maximum_iterations'] = 500

# Solve
up = Function(WP)
solver.solve(up.vector(), b)

pc_fieldsplit_detect_saddle_point does what the name says (looks for zeros on the diagonal of A), so you don't need to specify the dofs.
Note that you forgot solver.set_from_options(), needed to pass the PETScOptions to the solver.

Comments

Interesting. So pc_fieldsplit_detect_saddle_point magically splits up into u and p?

---

Yes. It finds the rows where the diagonal entries are zero, these are the p dofs (fieldsplit_1), the rest is u (fieldsplit_0).

---

If you want to set the fieldsplit manually, you can do it via the petsc4py interface. Maybe there is an easier way with the dolfin, I don't know. I use a Stokes fieldsplit like this:

from petsc4py import PETSc
# ...

ksp = PETSc.KSP().create()
PETScOptions.set('ksp_view')
PETScOptions.set('ksp_monitor_true_residual')
PETScOptions.set('pc_type', 'fieldsplit')
PETScOptions.set('pc_fieldsplit_type', 'additive')
# PETScOptions.set('pc_fieldsplit_detect_saddle_point')
PETScOptions.set('fieldsplit_0_ksp_type', 'preonly')
PETScOptions.set('fieldsplit_0_pc_type', 'lu')
PETScOptions.set('fieldsplit_1_ksp_type', 'preonly')
PETScOptions.set('fieldsplit_1_pc_type', 'jacobi')
ksp.setFromOptions()

ksp.setOperators(A)

pc = ksp.getPC()
is0 = PETSc.IS().createGeneral(W.sub(0).dofmap().dofs())
is1 = PETSc.IS().createGeneral(W.sub(1).dofmap().dofs())
fields = [('0', is0), ('1', is1)]
pc.setFieldSplitIS(*fields)

ksp.setUp()
ksp.solve(b, x)

A, b, x are PETScMatrix/Vectors. I find this useful, for example, if you want to use Schur preconditioners. I don't think this is possible yet with the standard dolfin interface.