On Sun, Aug 9, 2020 at 2:37 PM Nidish <[email protected]> wrote:
>
> On 8/9/20 1:45 AM, Barry Smith wrote:
>
> On Aug 9, 2020, at 1:14 AM, Nidish <[email protected]> <[email protected]> wrote:
>
> Hello,
>
> My question is related to solving a generic nonlinear algebraic system in
> parallel. I'm using the word generic to point to the fact that the system may
> or may not be from a finite element/mesh type object - it's just a residual
> function and a Jacobian function that I can provide.
>
> I followed the serial example in $PETSC_DIR/src/snes/tutorials/ex1.c pretty
> much verbatim (albeit with a slightly different objective function whose size
> can be controlled), and have been trying to extend it to a parallel
> application (the given example is strictly serial).
>
> When I run my code in serial everything works well and I get the solution I
> expect. However, even if I use just two cores, the code starts returning
> garbage. I had a feeling it has something to do with index ordering, but I
> want to get this code working without using a DM object. I'm attaching my
> code with this email.
>
> TL;DR Version: What modifications must be made in
> "$PETSC_DIR/src/snes/tutorials/ex1.c" to make it work in Parallel? Is it
> possible to make minimal modifications to make it run in parallel too?
>
> This is wrong:
>
> VecGetSize(x, &N);
> VecGetOwnershipRange(x, &is, &ie);
> /* VecGetArrayRead(x, &xx); */
> VecGetArray(x, &xx);
> for (int i=is; i<ie; ++i)
>
> You appear to be trying to mimic a DM example but without a DM.
>
> VecGetArray() always returns an array with indices from 0 to
> VecGetLocalSize()
>
> Thus you need to make your loops over local indices, not global indices.
>
> In addition you are likely to need "ghost" values of the input vectors. If
> you are not using DM then you need to have code that determines what ghost
> values are needed and create a VecScatter to manage communicating those ghost
> values. Look for examples that use VecScatterBegin in the formfunction
> routine for how this may be done.
>
> Barry
>
> I wasn't accounting for the ghost values! My compiler wasn't giving
> segfaults when I was accessing the arrays returned by VecGetArray beyond
> the allocation size. (the loop indices was me trying to be smart, I was
> accessing xx by xx[i -is] and so on).
>
> I've changed my implementation with VecCreateGhost - the code's slightly
> longer but it works, thank you for this!
>
>
> I also had another related question about the storage of the arrays from
> VecGetArray() calls. Do these functions allocate new locations in memory
> (and can be expected to be slow)? Even for the ghosted vectors I found
> distinct addresses for the ghost values. I'm assuming the addresses are
> distinct to accommodate for cluster programming, but should a user be
> judicious in calls to VecGetArray() ?
>
No. For the vectors you are using, this is a no-copy operation.
Thanks,
Matt
> Thank you,
> Nidish
>
> Thank You,
> Nidish
> <sness.c>
>
> --
> Nidish
>
--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
