[petsc-users] issue with ghost points

[Barry Smith]

  You are using the info struct without ever having set values into it  You are, at a minimum, missing a call to DMDAGetLocalInfo(). Valgrind would immediately have found the use of uninitialized values. 

   Note that DMDAGetCorners() provides some of the information that DMDAGetLocalInfo() provides but in a slightly different format.


  Barry



> On Feb 27, 2017, at 9:18 PM, Gideon Simpson <gideon.simpson at gmail.com> wrote:
> 
> 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
>