On Fri, Jun 10, 2011 at 17:56, Derek Gaston <firstname.lastname@example.org>
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().