Alright, thank you.
Nidish
On 8/9/20 1:43 PM, Matthew Knepley wrote:
On Sun, Aug 9, 2020 at 2:37 PM Nidish <[email protected]
<mailto:[email protected]>> wrote:
On 8/9/20 1:45 AM, Barry Smith wrote:
On Aug 9, 2020, at 1:14 AM, Nidish<[email protected]> <mailto:[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/>
--
Nidish
