[petsc-users] Correct way to access a sequential version of a DMDA Vec?

[Smith, Barry F.]

   This won't work as written for two reasons

1) the VecScatterCreateToAll() will just concatenate the values from each process in a long array on each process, thus the resulting values will be "scrambled" and it won't be practical to access the values (because the parallel layout of DMDA vectors is scrambled from the logical global 2d array)

2) As you note, DMDAVecGetArrayRead() cannot handle a large sequential vector not from a DM.

   Fortunately it is doable. 

1) You first need to convert the parallel vector from PETSc ordering (by process) to natural ordering (matching the underlying parallel grid) use 

DMDACreateNaturalVector
DMDAGlobalToNaturalBegin
DMDAGlobalToNaturalEnd

then 

> VecScatterCreateToAll(natural &scatter, &xs);
> VecScatterBegin(scatter, natural xs, INSERT_VALUES, SCATTER_FORWARD);
> VecScatterEnd(scatter, natural, xs, INSERT_VALUES, SCATTER_FORWARD);

Note you don't need to convert to local at all.

2)  Use VecGetArray2d() to access the full array on each process. You simply need to pass in the global sizes of the 2d array to get is "shape" right.

If any of the steps seem to cause grief I would debug with a tiny DMDA on 2 processes fill the initial vector up with easy to recognize values like 1 2 3 ...  and print the vectors as you move through each stage 

  Barry


> On Oct 27, 2019, at 6:25 PM, Ellen M. Price via petsc-users <petsc-users at mcs.anl.gov> wrote:
> 
> Hi PETSc users!
> 
> I'm trying to access a DMDA Vec's values after doing a scatter to all
> processors. I've run into some trouble, however, because the call to
> DMDAVecGetArray doesn't seem to care that I'm using a sequential vector.
> The calling sequence that "works" (runs without error *until* I try to
> use the values) is:
> 
> // da is the 2-dimensional DMDA
> // x is a global, input vector
> // xs is the sequential vector
> // values is an array of structs suitable for this DMDA
> 
> DMGetLocalVector(da, &xloc);
> DMGlobalToLocalBegin(da, x, INSERT_VALUES, xloc);
> DMGlobalToLocalEnd(da, x, INSERT_VALUES, xloc);
> VecScatterCreateToAll(xloc, &scatter, &xs);
> VecScatterBegin(scatter, xloc, xs, INSERT_VALUES, SCATTER_FORWARD);
> VecScatterEnd(scatter, xloc, xs, INSERT_VALUES, SCATTER_FORWARD);
> DMDAVecGetArrayRead(da, xs, &values);
> 
> From a skim of the PETSc source code, it looks like the start parameters
> that get handed off to VecGetArray2d only depend on the DMDA and do not
> take into account that the vector might already be sequential. (I might
> be wrong about this, but that's what it looks like at first glance.) Is
> there a way around this that is more elegant than just modifying the
> output pointer?
> 
> Feel free to point me to an applicable message board post if there is
> one, but Google hasn't brought up anything useful yet.
> 
> Thanks,
> Ellen Price