[petsc-users] accessing global vector in a DMMA

[Zhang]

Dear All,
Now I am confused with the way to access a global vector defined from DMDA.

Here is the codes. When I switch on _DEBUG1_ the results get wrong.  While if it's off,

the results fine. I just wonder why I gave to use the local form of ul,vl,wl to access the
values such as ul[k][j][i+1], and etc.

Thank you first for any suggestion.

Zhenyu

ierr = DMCreateGlobalVector(ctxu->grid, &ctxu->x );CHKERRQ(ierr);
  ierr = VecDuplicate(ctxu->x,&ctxu->b);CHKERRQ(ierr);
ierr = DMCreateGlobalVector(ctxv->grid, &ctxv->x );CHKERRQ(ierr);
    ierr = VecDuplicate(ctxv->x,&ctxv->b);CHKERRQ(ierr);
ierr = DMCreateGlobalVector(ctxw->grid, &ctxw->x );CHKERRQ(ierr);
    ierr = VecDuplicate(ctxw->x,&ctxw->b);CHKERRQ(ierr);
...
  VecCopy(ctxu->x,ctxu->b);
  VecCopy(ctxv->x,ctxv->b);
  VecCopy(ctxw->x,ctxw->b);

  DMDAVecGetArray( ctxu->grid, ctxu->b, &ustar );
  DMDAVecGetArray( ctxv->grid, ctxv->b, &vstar );
  DMDAVecGetArray( ctxw->grid, ctxw->b, &wstar );

#if defined(_DEBUG1_)
  DMDAVecGetArray( ctxu->grid, ctxu->x, &u );
  DMDAVecGetArray( ctxv->grid, ctxv->x, &v );
  DMDAVecGetArray( ctxw->grid, ctxw->x, &w );
#else
  DMGetLocalVector(ctxu->grid,&ctxu->local);
  DMGetLocalVector(ctxv->grid,&ctxv->local);
  DMGetLocalVector(ctxw->grid,&ctxw->local);
  DMGlobalToLocalBegin(ctxu->grid,ctxu->x,INSERT_VALUES,ctxu->local);
  DMGlobalToLocalEnd(ctxu->grid,ctxu->x,INSERT_VALUES,ctxu->local);
  DMGlobalToLocalBegin(ctxv->grid,ctxv->x,INSERT_VALUES,ctxv->local);
  DMGlobalToLocalEnd(ctxv->grid,ctxv->x,INSERT_VALUES,ctxv->local);
  DMGlobalToLocalBegin(ctxw->grid,ctxw->x,INSERT_VALUES,ctxw->local);
  DMGlobalToLocalEnd(ctxw->grid,ctxw->x,INSERT_VALUES,ctxw->local);
  DMDAVecGetArray( ctxu->grid, ctxu->local, &ul );
  DMDAVecGetArray( ctxv->grid, ctxv->local, &vl );
  DMDAVecGetArray( ctxw->grid, ctxw->local, &wl );
#endif

//----------------------------------------------------------------
  // U
  DMDAGetCorners( ctxu->grid, &is, &js, &ks, &in, &jn, &kn );
  ie = is + in - 1;
  je = js + jn - 1;
  ke = ks + kn - 1;

  is=max(is,1);
  js=max(js,1);
  ks=max(ks,1);
  ie=min(ie,ctxu->l-2);
  je=min(je,ctxu->m-2);
  ke=min(ke,ctxu->n-2);


  for (k=ks; k<=ke; k++) {
    for (j=js; j<=je; j++) {
      for (i=is; i<=ie; i++) {

#if defined(_DEBUG1_)
        ustar[k][j][i] +=
                               - dtdx*(0.25*((u[k][j][i]+u[k][j][i+1])*(u[k][j][i]+u[k][j][i+1]))
                                       - 0.25*((u[k][j][i]+u[k][j][i-1])*(u[k][j][i]+u[k][j][i-1])))
                               - dtdy*(0.25*(u [k][j][i]+u [k][j+1][i])*(v [k][j][i]+v [k][j][i+1])
                                       - 0.25*(u [k][j][i]+u [k][j-1][i])*(v [k][j-1][i]+v [k][j-1][i+1]))
                               - dtdz*(0.25*(u [k][j][i]+u [k+1][j][i])*(w [k][j][i]+w [k][j][i+1])
                                       - 0.25*(u [k][j][i]+u [k-1][j][i])*(w [k-1][j][i]+w [k-1][j][i+1]))
                               + dtdxx*(u [k][j][i-1]-2*u [k][j][i]+u [k][j][i+1])
                               + dtdyy*(u [k][j-1][i]-2*u [k][j][i]+u [k][j+1][i])
                               + dtdzz*(u [k-1][j][i]-2*u [k][j][i]+u [k+1][j][i]);
#else
        ustar[k][j][i] +=
                               - dtdx*(0.25*((ul[k][j][i]+ul[k][j][i+1])*(ul[k][j][i]+ul[k][j][i+1]))
                                       - 0.25*((ul[k][j][i]+ul[k][j][i-1])*(ul[k][j][i]+ul[k][j][i-1])))
                               - dtdy*(0.25*(ul[k][j][i]+ul[k][j+1][i])*(vl[k][j][i]+vl[k][j][i+1])
                                       - 0.25*(ul[k][j][i]+ul[k][j-1][i])*(vl[k][j-1][i]+vl[k][j-1][i+1]))
                               - dtdz*(0.25*(ul[k][j][i]+ul[k+1][j][i])*(wl[k][j][i]+wl[k][j][i+1])
                                       - 0.25*(ul[k][j][i]+ul[k-1][j][i])*(wl[k-1][j][i]+wl[k-1][j][i+1]))
                               + dtdxx*(ul[k][j][i-1]-2*ul[k][j][i]+ul[k][j][i+1])
                               + dtdyy*(ul[k][j-1][i]-2*ul[k][j][i]+ul[k][j+1][i])
                               + dtdzz*(ul[k-1][j][i]-2*ul[k][j][i]+ul[k+1][j][i]);
        #endif
      }
    }
  }

.....

#if defined(_DEBUG1_)
DMDAVecRestoreArray( ctxu->grid, ctxu->x, &u );
DMDAVecRestoreArray( ctxv->grid, ctxv->x, &v );
DMDAVecRestoreArray( ctxw->grid, ctxw->x, &w );
#else
DMDAVecRestoreArray( ctxu->grid, ctxu->local, &ul );
DMDAVecRestoreArray( ctxv->grid, ctxv->local, &vl );
DMDAVecRestoreArray( ctxw->grid, ctxw->local, &wl );
DMRestoreLocalVector(ctxu->grid,&ctxu->local);
DMRestoreLocalVector(ctxv->grid,&ctxv->local);
DMRestoreLocalVector(ctxw->grid,&ctxw->local);
#endif

DMDAVecRestoreArray( ctxu->grid, ctxu->b, &ustar );
DMDAVecRestoreArray( ctxv->grid, ctxv->b, &vstar );
DMDAVecRestoreArray( ctxw->grid, ctxw->b, &wstar );


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