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

304

views

0

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/

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))
pc.setFieldSplitType(0) # 0=additive
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

Please login to add an answer/comment or follow this question.

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)