On Fri, Jun 10, 2011 at 17:56, Derek Gaston <[email protected]> wrote:
> This happens during the call to update() in the Petsc callback to compute
> the residual (and other callbacks). It's not from user code at all. Petsc
> is actually giving us a vector that isn't closed.
This comes from an impedance mismatch between the interfaces. Libmesh uses
the term "closed" to refer to the combination of PETSc's "assembled" (values
in the global vector have been communicated) and "ghost values are current".
The PETSc interfaces were designed with the idea that you would call
VecGhostUpdateBegin/End when you wanted ghost values to be current,
otherwise they could be stale indefinitely. The ghost values are not used
for any linear algebra operations (norms, dot products, matix-vector
multiplies, etc). They are semantically independent of the mathematical
object "vector", and could be stored somewhere else entirely without
affecting semantics.
An unassembled vector, on the other hand, is useless in terms of linear
algebra. If you perform local modification, then you have to assemble before
the solver can do anything useful with it.
The problem here (at least from my perspective) is that Libmesh provides
operator() as a non-collective operation that is supposed to be defined
without the user doing anything. If it was collective on the first call,
then you could do the VecGhostUpdate() at that time, but that would be a
horrible API because then you would get deadlock if one process had an empty
subdomain or was otherwise degenerate so that the updates were called out of
order.
PETSc makes no guarantee about what is stored in the ghost values when it
comes to you from the solver. The solvers don't even know about the
existence of ghost points and it would be inefficient if they did. The
vector you get is assembled, but the ghost values are not current. According
to closed(), you have to assemble and then communicate ghost values.
Communicating ghost values is a simpler operation than assembling a vector
because assembly requires all-to-all communication to find out who needs to
talk to who (because the API allows VecSetValues() to set any value on any
process from any other process). Modifying a local vector (e.g. what you get
with VecGhostGetLocalForm()) and calling
VecGhostUpdateBegin/End(...,SCATTER_REVERSE) is more efficient than
VecAssemblyBegin/End().
Note that a vector only loses its "assembled" status through VecSetValues().
In the vast majority of finite element computations, you only need access to
neighbors so you could just use VecGhostGetLocalForm(), set the entries
directly, and then VecGhostUpdateBegin/End().
------------------------------------------------------------------------------
EditLive Enterprise is the world's most technically advanced content
authoring tool. Experience the power of Track Changes, Inline Image
Editing and ensure content is compliant with Accessibility Checking.
http://p.sf.net/sfu/ephox-dev2dev
_______________________________________________
Libmesh-devel mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-devel
