[petsc-users] VecScatter

[Jorti, Zakariae]

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/96458d51/attachment-0001.html>