PETScViewer and FEniCS: dump/load vectors and matrices

Question

Dear FEniCS community,
I need to dump and load PETSc vectors and matrices in my FEniCS application. Do you think that it is possible/meaningful to do so using PETScViewer?

I have tried the following (for a vector), in serial. PETSc error is enclosed at the end of the code.

import sys, petsc4py
petsc4py.init(sys.argv)
from petsc4py import PETSc
from dolfin import *

mesh = UnitSquareMesh(5,5)
V = FunctionSpace(mesh, "CG", 2)

v = Function(V)
v.vector()[:] = 1

viewer = PETSc.Viewer().createBinary("petsc4py-viewer.dat", "w")
viewer(as_backend_type(v.vector()).vec())

viewer = PETSc.Viewer().createBinary("petsc4py-viewer.dat", "r")
data = PETSc.Vec().load(viewer)
vec = PETScVector(data)

print vec.array() # Fine, vector of [1, 1, .... 1]

# a not so smart way to double the vector, which however
# is representative of the (more complex) logic of my final program
output = vec.copy()
#output.apply("") # or
#vec.apply("")    # do not fix the issue
output.add_local(vec.array()) # Error (below)
# Error:   Unable to successfully call PETSc function 'VecSetValuesLocal'.
# Reason:  PETSc error code is: 73.
# Where:   This error was encountered inside /scratch/fballarin/src/fenics/dolfin/dolfin/la/PETScVector.cpp.
# Process: 0

Moreover, would this work with different number of processors when reading and writing? If not, how should I take into account the different ordering?

Thanks,
Francesco

Comments

A direct call to petsc4py on the variable data reveals a more helpful error message:

[0] VecSetValuesLocal() line 1055 in petsc/vec/vec/interface/rvector.c
[0] Object is in wrong state
[0] Local to global never set with VecSetLocalToGlobalMapping()

Do you have any suggestion on how to set the local to global mapping? How would that change for a different numbers of processors when reading/writing?