Petsc4py.IS matrices

General
1

Hello,

I am interested in using petsc4py.IS matrices as they are needed for KSPFETIDP. The matrix must be subassembled locally with the Neumann problems and PETSc takes care of the rest. However, I have not been able to find proper documentation on how to do this.

Moreover, should you be aware of a tutorial on the KSPFETIDP different than the C++ one, it would be greatly appreciated.

Thanks a lot.

2

To follow-up on this topic, I have the following MWE:

from petsc4py import PETSc
import dolfinx, ufl
from mpi4py import MPI
from ufl import inner, grad, dx
import numpy as np

comm = MPI.COMM_WORLD
mesh = dolfinx.mesh.create_rectangle(comm, [[0, 0], [1, 1]], [10, 10], dolfinx.mesh.CellType.quadrilateral)

V = dolfinx.fem.FunctionSpace(mesh, ('CG', 1))
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)

f = 1.0
a = inner(grad(u), grad(v))
L = inner(f, v)

zero_dofs = dolfinx.fem.locate_dofs_geometrical(V, lambda x: np.logical_or(np.logical_or(np.isclose(x[0], 0), np.isclose(x[0], 1)), np.logical_or(np.isclose(x[1], 0), np.isclose(x[1], 1))))
bc = dolfinx.fem.dirichletbc(dolfinx.fem.Constant(mesh, 0.0), zero_dofs, V)

A = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(a * dx), bcs=[bc])
A.assemble()
b = dolfinx.fem.petsc.assemble_vector(dolfinx.fem.form(L * dx))

dolfinx.fem.petsc.apply_lifting(b, [dolfinx.fem.form(a * dx)], bcs=[[bc]])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
b.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
dolfinx.fem.petsc.set_bc(b, [bc])

imap = V.dofmap.index_map
lgmap = PETSc.LGMap().create(imap.local_to_global(np.arange(imap.size_local).tolist()).tolist(), comm=comm)
mat = PETSc.Mat().createIS((imap.size_global, imap.size_global), lgmap, lgmap, comm=comm)
print(comm.rank, lgmap.size, mat.getLocalSize())

ownership_ = A.getOwnershipRange()
local_vals = A.getValues(np.arange(ownership_[0], ownership_[1], dtype=np.int32), np.arange(ownership_[0], ownership_[1], dtype=np.int32))
A_local = PETSc.Mat().createAIJ(size=(imap.size_local, imap.size_local), comm=MPI.COMM_SELF)
A_local.setUp()

A_local.setValues(np.arange(imap.size_local, dtype=np.int32), np.arange(imap.size_local, dtype=np.int32), local_vals)
mat.setISLocalMat(A_local)
mat.assemble()

ksp = PETSc.KSP().create(comm)
ksp.setOperators(mat)
ksp.setType('fetidp')
vec = b.copy()
t = dolfinx.common.Timer()
ksp.solve(b, vec)

print(f'Solve took {t.elapsed()[0]}')
print(f'Matrix size {mat.size[0]}')

vec.ghostUpdate(
addv=PETSc.InsertMode.INSERT,
mode=PETSc.ScatterMode.FORWARD,
)

uh = dolfinx.fem.Function(V)
uh.vector.setArray(vec.array)
uh.name = 'uh_FOM'

file = dolfinx.io.XDMFFile(comm, 'example.xdmf', 'w')
file.write_mesh(mesh)
file.write_function(uh)

Something weird happens when initializing the mat object, indeed the output of the above is:

0 59 (61, 61)
1 62 (60, 60)

and then throws an error in the ksp.solve because there are (of course) inconsistent arrays. It’s as if there is a switch between DOFs assigned. It seems to me that, however, one cannot pass a list of LGMap(s) to the createIS function, so how should one go about it?

Thank you.