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 <[email protected]> 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 >
