[petsc-dev] Vector operations on ghosted vectors

[Jed Brown]

On Fri, Mar 1, 2013 at 10:37 AM, John Mousel <elafint.john at gmail.com> wrote:

> Let me be even more specific. I'm implementing a velocity correction
> scheme in rotational form (reference given below). In the first step,
> forgetting even a sub-iterative loop, a projecting pressure is solved for
> from
>
> L(p(n+1)) = -div(beta*u(n) + gamma*u(n-1) + exci*u(n-2) + H(n+1) - H(n) +
> f(n+1) - f(n))
>
> Here, L and div are Laplacian and divergence operators, and H(n) and f(n)
> are the non-linear and body force at time level n. In my mind, to form part
> of the RHS, I can use VecAXPY on all the terms in the div(), and then
> perform the divergence on that. To avoid communication, I can just use the
> VecAXPY on the ghosted portion as well, which since I've copied the entire
> ghosted vector during the VecCopy(u(n),u(n-1)) operation, should make sense.
> Is that not a wise thing to do?
>

It is common to be able to assemble the RHS without updating ghost points
in the local form, but that's not the case here. What you are doing is
reasonable, but it couples the time integration method to the physics more
tightly than I tend to like. In any case, you need not run the VecAXPY on
the global form at all. Since you used the local form in the solve for
u(n), it will still be valid on the next iteration. The result of the
VecAXPY is a pure work vector so it can be a local vector (rather than the
local form of some global vector).

I predict that this is all premature optimization because the solve does
much more communication and that cutting corners here will only serve to
make your code more brittle without providing meaningful performance gains.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20130302/a0f6c143/attachment.html>