Read/Write PETSc Matrices in legacy dolfin

I/O
1

Hi all,

Is there a way to read and write matrices in parallel, using for instance PETSc.Viewer() similar to what is given here. What I am trying currently:

from petsc4py import PETSc
from dolfin import *
comm = MPI.comm_world
comm_rank = comm.rank

mesh = UnitSquareMesh(comm, 10, 10)
V = FunctionSpace(mesh, "CG", 1)
du = TrialFunction(V)
v = TestFunction(V)
u = Function(V)
f = Constant(1.0)
bilinear_form = inner(grad(du), grad(v))*dx
linear_form = v*f*dx
A_ = PETScMatrix(comm)
b_ = PETScVector(comm)
A_, b_ = assemble_system(bilinear_form, linear_form, A_tensor=A_, b_tensor=b_)
solver = PETScLUSolver(comm, "mumps")
solver.solve(A_, u.vector(), b_)

# writing to file (seemingly works; not sure if correct)
viewer = PETSc.Viewer().createBinary("Amat_poisson.dat", "w", comm=comm)
viewer(A_.mat())
viewer = PETSc.Viewer().createBinary("bmat_poisson.dat", "w", comm=comm)
viewer(b_.vec())

# now loading back into the petsc matrix (doesn't work)
A_new = PETScMatrix(comm)
b_new = PETScVector(comm)
viewer = PETSc.Viewer().createBinary("Amat_poisson.dat", "r", comm=comm)
_A = PETSc.Mat().load(viewer)
_A.copy(result=A_new.mat())
viewer = PETSc.Viewer().createBinary("bmat_poisson.dat", "r", comm=comm)
_b = PETSc.Vec().load(viewer)
_b.copy(b_new.vec())

This fails with the Traceback

    _A.copy(A_new.mat())
  File "PETSc/Mat.pyx", line 686, in petsc4py.PETSc.Mat.copy
petsc4py.PETSc.Error: error code 60
2

The following works in serial but the assertion check fails in parallel

from petsc4py import PETSc
from dolfin import *
set_log_level(30)
comm = MPI.comm_world
comm_rank = comm.rank

mesh = UnitSquareMesh(comm, 10, 10)
V = FunctionSpace(mesh, "CG", 1)
du = TrialFunction(V)
v = TestFunction(V)
u = Function(V)
f = Constant(1.0)
bilinear_form = inner(grad(du), grad(v))*dx
linear_form = v*f*dx
A_ = PETScMatrix(comm)
b_ = PETScVector(comm)
A_, b_ = assemble_system(bilinear_form, linear_form, A_tensor=A_, b_tensor=b_)
solver = PETScLUSolver(comm, "mumps")
solver.solve(A_, u.vector(), b_)

viewer = PETSc.Viewer().createBinary("Amat_poisson.dat", "w", comm=comm)
viewer(A_.mat())
viewer = PETSc.Viewer().createBinary("bmat_poisson.dat", "w", comm=comm)
viewer(b_.vec())

# read values back in PETSc matrices
viewer = PETSc.Viewer().createBinary("Amat_poisson.dat", "r", comm=comm)
_A = PETSc.Mat().load(viewer)
A_new = PETScMatrix(_A)
viewer = PETSc.Viewer().createBinary("bmat_poisson.dat", "r", comm=comm)
_b = PETSc.Vec().load(viewer)
b_new = PETScVector(_b)

assert A_.mat().equal(A_new.mat()) and b_.vec().equal(b_new.vec())

with

error code '60' (Nonconforming object sizes)
1 Like