Function -> PETScVector -> Vec -> PETScVector -> Function is losing information

[Patrick Farrell]

Consider the following diff to the petsc4py demo:

diff --git demo/undocumented/petsc4py/python/demo_petsc4py.py demo/undocumented/petsc4py/python/demo_petsc4py.py
index 4449f7c..6dd87a9 100644
--- demo/undocumented/petsc4py/python/demo_petsc4py.py
+++ demo/undocumented/petsc4py/python/demo_petsc4py.py
@@ -82,6 +82,12 @@ A_petsc = as_backend_type(A).mat()
 b_petsc = as_backend_type(b).vec()
 x_petsc = as_backend_type(u.vector()).vec()

+output = File("/tmp/output.pvd")
+output << u
+
+u_again = Function(V, PETScVector(x_petsc))
+output << u_again
+
 # Create solver, apply preconditioner and solve system
 ksp = PETSc.KSP().create()
 ksp.setOperators(A_petsc)

We try to recreate the Function from its Vec, and output it for visualisation. This runs fine in parallel, but with 2 processors:

Traceback (most recent call last):
  File "demo_petsc4py.py", line 89, in <module>
    output << u_again
*** -------------------------------------------------------------------------
*** Error:   Unable to successfully call PETSc function 'VecGetValues'.
*** Reason:  PETSc error code is: 63.
*** Where:   This error was encountered inside /home/pef/src/dolfin/git/dolfin/la/PETScVector.cpp.
*** Process: unknown
*** 
*** DOLFIN version: 1.3.0+
*** Git changeset:  a846a5677af85623582b381a737ca4a2d2256859
*** -------------------------------------------------------------------------

It seems that the chain Function -> PETScVector -> Vec -> PETScVector -> Function is losing information somehow.

[Jan Blechta]

Can't replicate with your branch pefarrell/petscmatrix-remove-use-gpu merged into master.

[Jan Blechta]

Sorry, here PETScMatrix-fix plays no role. Even though, can't replicate.

[Prof Garth Wells]

@blechta Do you have petsc4py installed? I've seen some issues in the Python interface that only show up when petsc4py is installed.

[Jan Blechta]

@garth-wells Yes.

[Jan Blechta]

It seems to me that problem is that PETScVector::x_ghosted, ghost_global_to_local are not initialized when constructing from Vec. (Maybe, I don't reproduce the error because partitioning differes on my machine and x_ghosted is not needed when gathering values for VTK output.)

[Patrick Farrell]

Garth: Can you replicate?

Jan: I'm on petsc-dev, perhaps that's a difference?

[Jan Blechta]

@pefarrell I think I've suggested a correct explanation in my last post - uninitialized x_ghosted, ghost_global_to_local...

[Joachim Berdal Haga]

See also issue #270, which is basically the same problem for Trilinos.

[Prof Garth Wells]

The upcoming switch to local dof will fix this as ghost_global_to_local will be redundant.

[Prof Garth Wells]

Will be taken care of with switch to local dof indices in 1.5.

[Jan Blechta]

petsc4py demo now plots and outputs wrong ghost values. Adding

as_backend_type(u.vector()).update_ghost_values()

before plot command fixes the problem.

[Jan Blechta]

Or

u.vector().apply('insert')

helps which looks more compact.

[Prof Garth Wells]

Fixed with switch to local indices.

[Jan Blechta]

Ok, I will open new issue.

[Jan Blechta]

See #425.