[petsc-users] VecScatter

[Barry Smith]

   You need to first convert the Vec to the "natural" ordering and then bring that vector down to one rank. Something like

DMDACreateNaturalVector(dm,&n);
DMDAGlobalToNaturalBegin/End(dm,   n);
> VecScatterCreateToZero(n, &scat, &Xseq);
>   VecScatterBegin(scat, n Xseq, INSERT_VALUES, SCATTER_FORWARD);
>   VecScatterEnd(scat, n,  Xseq, INSERT_VALUES, SCATTER_FORWARD);



> On Oct 11, 2022, at 7:45 PM, Jorti, Zakariae via petsc-users <petsc-users at mcs.anl.gov> wrote:
> 
> Hello, 
> 
> I have a code that handles a PETSc Vec on many procs and I would like to use VecScatterCreateToZero to have all elements of this vector on a single proc, call VecGetArrayRead on this proc to get the corresponding array to carry out some diagnostics.
> Unfortunately, what I noticed is that the ordering of the initial Vec is not preserved after the VecScatterCreateToZero call. Is there a way to have the same ordering for both Vecs? 
> You will find below a minimal example that demonstrates the issue.
> Many thanks,
> 
> Zakariae   
> 
> -------------------------------------------------------------
> 
> static char help[] = "Demonstrates ordering change after VecScatterCreateToZero call.\n\n";
> 
> 
> #include <petscpf.h>
> #include <petscdm.h>
> #include <petscdmda.h>
> 
> int main(int argc,char **argv)
> {
>   Vec            xy;
>   DM             da;
>   PetscErrorCode ierr;
>   PetscInt       m = 11, n = 11, dof = 2;
>   PetscMPIInt rank;
>   DMDACoor2d **coors;
> 
>   ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
>   ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,m,n,PETSC_DECIDE,PETSC_DECIDE,dof,1,0,0,&da);CHKERRQ(ierr);
>   ierr = DMSetFromOptions(da);CHKERRQ(ierr);
>   ierr = DMSetUp(da);CHKERRQ(ierr);
>   ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);CHKERRQ(ierr);
>   ierr = DMGetCoordinates(da,&xy);CHKERRQ(ierr);
> 
>   PetscInt      i, j, ixs, ixm, iys, iym;
>   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
> 
>   DMDAGetCorners(da, &ixs, &iys, 0, &ixm, &iym, 0);
>   DMDAVecGetArray(da, xy, &coors);
>   for (j = iys; j < iys + iym; j++) {
>     for (i = ixs; i < ixs + ixm; i++) {
>         PetscPrintf(PETSC_COMM_SELF, "rank=%d, %d, %d, (%g, %g)\n",rank, 
>                     i, j,coors[j][i].x,coors[j][i].y);
>     }
>   }
>   DMDAVecRestoreArray(da, xy, &coors);
> 
>   VecScatter  scat;
>   Vec         Xseq;
>   const PetscScalar *array; 
>    
>   /* create scater to zero */
>   VecScatterCreateToZero(xy, &scat, &Xseq);
>   VecScatterBegin(scat, xy, Xseq, INSERT_VALUES, SCATTER_FORWARD);
>   VecScatterEnd(scat, xy, Xseq, INSERT_VALUES, SCATTER_FORWARD);
> 
>   if (rank == 0) {
>     PetscInt sizeX;
>     VecGetSize(Xseq, &sizeX); 
>     PetscPrintf(PETSC_COMM_SELF,"The size of Xseq is %d, and the grid size is %d\n",sizeX,11*11);
>     VecGetArrayRead(Xseq, &array);
> 
>     for (j = 0; j < 11; j++) {
>       for (i = 0; i < 11; i++) {
>         PetscPrintf(PETSC_COMM_SELF, "%d,%d, (%g,%g)\n", i, j, (double)array[2*(j*11+i)], (double)array[1+2*(j*11+i)]);
>       }
>     }
>     VecRestoreArrayRead(Xseq, &array);
>   }
> 
>   /*
>      Free work space.  All PETSc objects should be destroyed when they
>      are no longer needed.
>   */
>   ierr = DMDestroy(&da);CHKERRQ(ierr);
>   ierr = PetscFinalize();
>   return ierr;
> }

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20221011/d1be13d5/attachment.html>