[petsc-users] issue with ghost points

[Gideon Simpson]

Not sure if this is a bug or not, but, for the following code:

#include <petscvec.h>
#include <petscdm.h>
#include <petscdmda.h>

typedef struct{
  PetscScalar u;
  PetscScalar v;
}b_node;


int main(int argc, char **args)
{
  char dataname[PETSC_MAX_PATH_LEN] = "testvec.bin";
  PetscInt N = 8;
  DM dm;
  Vec b;
  PetscScalar mass, eng;
  PetscViewer viewer;
  DMDALocalInfo info;
  b_node * b_array;
  Vec b_local;
  PetscInt local_index, local_width, i;
  PetscScalar u, v;

  
  PetscInitialize(&argc,&args,NULL,NULL);

  PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);
  
  DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_GHOSTED, N, 2 , 1 , NULL, &dm);
  DMCreateGlobalVector(dm ,&b);

  /* Load vector from disk */
  PetscViewerBinaryOpen(PETSC_COMM_WORLD,dataname, FILE_MODE_READ, &viewer);
  VecLoad(b, viewer);
  PetscViewerDestroy(&viewer);

  /* Inspect with viewer */
  VecView(b, PETSC_VIEWER_STDOUT_WORLD);

  DMGetLocalVector(dm,&b_local);
  DMGlobalToLocalBegin(dm,b,INSERT_VALUES,b_local);
  DMGlobalToLocalEnd(dm,b,INSERT_VALUES,b_local);

  DMDAVecGetArray(dm, b_local, &b_array);
  DMDAGetCorners (dm, &local_index, NULL, NULL,&local_width, NULL, NULL);

  /* Zero out the ghost points */
  if(info.xs == 0){
    b_array[-1].u =0.0;
    b_array[-1].v =0.0;    
  }
  if(info.xs + info.xm == info.mx){
    b_array[info.mx].u =0.0;
    b_array[info.mx].v =0.0;    
  }

  /* Manually inspect */
  for(i=local_index;i<local_index+local_width;i++){
    u = b_array[i].u;
    v = b_array[i].v;
    PetscPrintf(PETSC_COMM_WORLD," %i: u = %g, v = %g\n", i, u, v);
  }

  DMDAVecRestoreArray (dm,b_local,&b_array);
  DMRestoreLocalVector(dm,&b_local);
  
  VecDestroy(&b);
  DMDestroy(&dm) ;

  PetscFinalize();
  return 0;
  
}

I am getting the output 

Vec Object: 1 MPI processes
  type: seq
Vec Object:Vec_0x84000000_0 1 MPI processes
  type: mpi
Process [0]
1.
0.
0.707107
0.707107
6.12323e-17
1.
-0.707107
0.707107
-1.
1.22465e-16
-0.707107
-0.707107
-1.83697e-16
-1.
0.707107
-0.707107
 0: u = 0., v = 0.
 1: u = 0.707107, v = 0.707107
 2: u = 6.12323e-17, v = 1.
 3: u = -0.707107, v = 0.707107
 4: u = -1., v = 1.22465e-16
 5: u = -0.707107, v = -0.707107
 6: u = -1.83697e-16, v = -1.
 7: u = 0.707107, v = -0.707107

The input vector is an alternating sequence of cos( pi / 4 * j), sin(pi/4 * j), for j = 0,1,2…7.  What I want to point out is that the first element of this vector, as it is read in,  is 1.  But when we do the manually looping, it spits out u = 0 at index i = 0.  If i turn off the piece of code that zeros out the ghost points, the problem disappears.  The reason I am concerned is that this kind of looping through the local vector with ghost points shows up in a piece of my code, and I am now concerned about the output.  Do I have a bug here?  



-gideon

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