[petsc-users] issue with ghost points

[Gideon Simpson]

Thanks, as always. 

Interesting that -Wall (with clang) didn’t catch that.  

-gideon

> On Feb 27, 2017, at 10:28 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> 
> 
>  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
>> 
> 

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