### Assigning PETSc Vector to the given function space

120
views
-1
4 months ago by

// I am trying to solve a P2P1 mixed formulation of the Stokes flow.
V_velocity   = dolfin.VectorFunctionSpace(mesh, "CG", 2)
V_pressure   = dolfin.      FunctionSpace(mesh, "CG", 1)
V = V_velocity * V_pressure
u_out = dolfin.Function(V)
// I have my Petsc compiled outside of fenics. I run MPI version of the PETSc solve
// outside fenics and now want to import the solution of the linear solve for post
// processing.
petsc_viewer = PETSc.Viewer()
petsc_viewer.createBinary(name='/tmp/solution.dat', mode=0)
solution_vec = PETSc.Vec().create()
// How do I now assign the solution_vec to the internal Generic vector of the u_out
u_out.vector() = solution_vec // this is not working.

​
Community: FEniCS Project
1

u_out.vector()[:] = solution_vec​

?

written 4 months ago by Hernán Mella

It gives me this error. Any other suggestion?

File "simpleflow.py", line 271, in <module>
sys.exit(main(sys.argv))
File "simpleflow.py", line 234, in main
u_out.vector()[:] = solution_vec
File "/usr/lib/python2.7/dist-packages/dolfin/cpp/la.py", line 1405, in __setslice__
raise IndexError("can only set full slices v[:]")
IndexError: can only set full slices v[:]​
written 4 months ago by Hardik Kabaria

The following works for me in 2017.2.  (Sorry if this appears twice; I may have run into some sort of weird glitch with the forum.)

from dolfin import *
from petsc4py import PETSc

mesh = UnitSquareMesh(10,10)

# I am trying to solve a P2P1 mixed formulation of the Stokes flow.
V_velocity   = VectorElement("CG", mesh.ufl_cell(), 2)
V_pressure   = FiniteElement("CG", mesh.ufl_cell(), 1)
VE = V_velocity * V_pressure
V = FunctionSpace(mesh,VE)

u_out = Function(V)

# I have my Petsc compiled outside of fenics. I run MPI version of the
# PETSc solve outside fenics and now want to import the solution of the linear
# solve for post processing.

# placeholder to create data file
viewer = PETSc.Viewer().createBinary("./solution.dat",'w')
viewer(as_backend_type(u_out.vector()).vec())

petsc_viewer = PETSc.Viewer()
petsc_viewer.createBinary('./solution.dat', 'r')
solution_vec = PETSc.Vec().create()
# How do I now assign the solution_vec to the internal Generic vector
# of the u_out

#u_out.vector() = solution_vec # this is not working.

# Runs without error:
u_out.vector().zero()
u_out.vector().axpy(1.0, PETScVector(solution_vec))

# Another alternative:
u_out.vector().set_local(PETScVector(solution_vec).get_local())
as_backend_type(u_out.vector()).vec().ghostUpdate()
as_backend_type(u_out.vector()).vec().assemble()

# Also runs without error, but may have issues.  See thread where axpy
# solution is suggested:
#
# https://fenicsproject.org/qa/10130/create-function-vector-containing-coefficient-finite-element/
#
# The API docs say this constructor is for internal library use only.
# Also, if I try to call set_local() after this, I get an
# "object in wrong state" PETSc error.
u_out = Function(V,PETScVector(solution_vec))


written 4 months ago by David Kamensky
What is the difference between

V_velocity = VectorElement("CG", mesh.ufl_cell(), 2)
V_velocity = VectorElement("CG", 2)
written 4 months ago by Hardik Kabaria
The second one gives me an error using version 2017.2.  I'm guessing you're using an older version of FEniCS, based on the construction of the mixed space (which is no longer valid in 2017.2).  I think this is just an API change between different versions.  I don't think this affects the vector manipulations (unless the API for that has changed too...).
written 4 months ago by David Kamensky