### Using fieldsplit on a petsc matnest object via cbc.block

304
views
0
5 months ago by
I'm trying to use petsc fieldsplit with a block matrix assembled via cbc.block. MWE below is a modified "Stokes_mixed.py" from cbc.block demos following the approach here:
https://github.com/MiroK/cbc.block-eig/blob/master/stokes.py

and the fieldsplit approach here:
https://fenicsproject.org/qa/5287/using-the-petsc-pcfieldsplit-in-fenics/

from dolfin import *
from block import *
from block.algebraic.petsc import *
from block.iterative import *

mesh = UnitSquareMesh(32, 32)

P2 = VectorElement("Lagrange", triangle, 2)
P1 = FiniteElement("Lagrange", triangle, 1)
TH = MixedElement([P2, P1])

W = FunctionSpace(mesh, TH)

f = Constant((0., 0.))

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

a = inner(grad(u), grad(v)) * dx \
+ p * div(v) * dx \
+ q * div(u) * dx

b = inner(grad(u), grad(v)) * dx + inner(u, v) * dx \
+ p * q * dx

L = inner(f, v) * dx

bcs = [DirichletBC(W.sub(0), (0., 0.), "on_boundary&&(x[1]<1-DOLFIN_EPS)"),
DirichletBC(W.sub(0), (1., 0.), "on_boundary&&(x[1]>1-DOLFIN_EPS)")]

# assemble as block matrices
A, y = block_assemble(a, L, bcs)

_, x = block_assemble(a, L, bcs) #doing this to get a usable solution vector

from petsc4py import PETSc

[[A00, A01], [A10, A11]] = A
[Y0,Y1] = y
[X0,X1] = x

A00a = A00
A00, A01, A10, A11 = [as_backend_type(mat).mat() for mat in (A00, A01, A10, A11)]
Y0,Y1 = [as_backend_type(b).vec() for b in (Y0,Y1)]
X0,X1 = [as_backend_type(b).vec() for b in (X0,X1)]

comm = mpi_comm_world()
AA = PETSc.Mat().createNest([[A00, A01], [A10, A11]], comm=comm)
bb = PETSc.Vec().createNest([Y0,Y1], comm=comm)
xx = PETSc.Vec().createNest([X0,X1], comm=comm)

ksp = PETSc.KSP().create()
ksp.setType(PETSc.KSP.Type.GMRES)

pc = ksp.getPC()
if False:
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))
subksps = pc.getFieldSplitSubKSP()
subksps[0].setType("preonly")
subksps[0].getPC().setType("lu")
subksps[1].setType("preonly")
subksps[1].getPC().setType("lu")

ksp.view()

ksp.setOperators(AA)
ksp.setFromOptions()
ksp.solve(bb,xx)

print 'done'
​

Note the code runs when the if statement is 'False', but encounters a petsc error:
petsc4py.PETSc.Error: error code 75
( PETSC_ERR_ARG_INCOMP 75 /* two arguments are incompatible */)
when enabled.

Community: FEniCS Project
As far as I know PETSc fieldsplit is intended for dealing with block systems which are assembled monolithically. So either assemble everything together (without cbc.block) and use fieldsplit, or assemble separately as matnest and use relevant functionality (not fieldsplit).
written 5 months ago by Jan Blechta
Really? I thought my problem might have been from the way I was identifying my 'index set' from the dof maps instead of functions like 'MatNestFindIS()' that I found reference too googling around...

As an interim solution I'd be happy to be able to assemble the system matrix by blocks monolithically so I can put that through a normal PETSc KSP solve if someone can point me to an example of that? (although I think MatNest would be more memory efficient)
written 5 months ago by Mike Welland
There's the matnest dev branch Stokes demo which may be of some help here.
written 5 months ago by Nate
Thanks. Doing some digging, it looks like this is for a new MatNest interface in FEniCS that relies on new MatNestGetISs functionality in petsc4py... maybe that's the solution that I'd have to use developmental tools. I'm hoping to run it on an old build (1.6.0) so it would be nice if I didn't have to recompile everything.
written 5 months ago by Mike Welland