The only way I know is to export the matrix to petsc4py and use the interface there.
The following code example uses all of
(1) MUMPS direct solver
(2) FEniCS implementation of a simple block preconditioner
(3) PETSc fieldsplit preconditioner
to solve the Stokes demo.
from dolfin import *
from petsc4py import *
import numpy as np
# basic setup copied from the fenics stokes demo
mesh = UnitCubeMesh(8, 8, 8)
V = VectorFunctionSpace(mesh, 'CG', 2)
Q = FunctionSpace(mesh, 'CG', 1)
W = V * Q
def right(x, on_boundary): return x[0] > (1.0 - DOLFIN_EPS)
def left(x, on_boundary): return x[0] < DOLFIN_EPS
def top_bottom(x, on_boundary):
return x[1] > 1.0 - DOLFIN_EPS or x[1] < DOLFIN_EPS
noslip = Constant((0.0, 0.0, 0.0))
bc0 = DirichletBC(W.sub(0), noslip, top_bottom)
inflow = Expression(('-sin(x[1]*pi)', '0.0', '0.0'))
bc1 = DirichletBC(W.sub(0), inflow, right)
zero = Constant(0)
bc2 = DirichletBC(W.sub(1), zero, left)
bcs = [bc0, bc1, bc2]
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
f = Constant((0.0, 0.0, 0.0))
a = inner(grad(u), grad(v))*dx + div(v)*p*dx + q*div(u)*dx
L = inner(f, v)*dx
m = inner(grad(u), grad(v))*dx + p*q*dx
(A, b) = assemble_system(a, L, bcs)
(M, _) = assemble_system(m, L, bcs)
# fenics direct solve
solver = LUSolver('mumps')
solver.set_operator(A)
U = Function(W)
solver.solve(U.vector(), b)
# fenics iterative solve
solver = KrylovSolver("tfqmr", "amg")
solver.set_operators(A, M)
U1 = Function(W)
solver.solve(U1.vector(), b)
diff = np.max(np.abs(U.vector().array()-U1.vector().array()))
# export to petsc4py
A = as_backend_type(A).mat()
b = as_backend_type(b).vec()
M = as_backend_type(M).mat()
# fieldsplit solve
ksp = PETSc.KSP().create()
ksp.setType(PETSc.KSP.Type.TFQMR)
pc = ksp.getPC()
pc.setType(PETSc.PC.Type.FIELDSPLIT)
is0 = PETSc.IS().createGeneral(W.sub(0).dofmap().dofs())
is1 = PETSc.IS().createGeneral(W.sub(1).dofmap().dofs())
pc.setFieldSplitIS(('u', is0), ('p', is1))
pc.setFieldSplitType(0) # 0=additive
subksps = pc.getFieldSplitSubKSP()
subksps[0].setType("preonly")
subksps[0].getPC().setType("hypre")
subksps[1].setType("preonly")
subksps[1].getPC().setType("hypre")
ksp.setOperators(A, M)
ksp.setFromOptions()
(x, _) = A.getVecs()
ksp.solve(b, x)
print(np.max(np.abs(U1.vector().array()-x.array)))
