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; }