[petsc-users] accessing global vector in a DMMA

Wed Aug 21 22:43:47 CDT 2013

On Aug 21, 2013, at 10:04 PM, Zhang <zyzhang at nuaa.edu.cn> wrote:

> Thank your help, Barry,
> 
> If I assign the sum of local values (, such as ul[i+1]-2*ul[i]+ul[i-1],) to a global value ustar[i], 
> 1. should I use the local value for ustar as well ?

   Generally with finite differences there is no reason to use the local ghosted vector for the result since each process calculates its own values and doesn't calculate any values that to into the ghosted points.

> 2. Once I finished the work of assigning values for ustar and restored them, should I use VecAssemblyBegin/End to them later?

   No, if you use any of the vecgetarray functions you do not use VecAssemblyBegin/End() (though it is harmless to call the VecAssemblyBegin/end because they end up doing nothing.

> Because I find such a usage in the examples. I am not sure it necessary or not.
> 3. AS a test if I choose the incorrect way, say, ustar[i] comes straightly from u[i+/-1]..., why I did not get any error, such as 
> memory access out of range? Or even this is dangerous, such errors will not be shown in a serial run?

   No errors in serial runs because there are no ghost points between processes. In parallel you will get wrong answers but may or may not get crashes due to memory access out of range. We highly recommend making some parallel runs with valgrind because it finds almost all memory access bugs even ones that don't crash the code. http://www.mcs.anl.gov/petsc/petsc-as/documentation/faq.html#valgrind

  Barry

> 
> Cheers,
> 
> Zhenyu
> 
> 
>> -----原始邮件-----
>> 发件人: "Barry Smith" <bsmith at mcs.anl.gov>
>> 发送时间: 2013-08-22 11:02:13 (星期四)
>> 收件人: Zhang <zyzhang at nuaa.edu.cn>
>> 抄送: petsc-users at mcs.anl.gov
>> 主题: Re: [petsc-users] accessing global vector in a DMMA
>> 
>> 
>>   In PETSc the "local vector" or "local representation" refers to vectors WITH ghost points. The global vector or global representation refers to vectors WITHOUT ghost points.  Hence to access locations like [i+1] which extend over to the next process you need to use the local ghosted representation. 
>> 
>>   Barry
>> 
>> One could argue that it is not the clearest names to use. Just remember local == ghosted and global == not ghosted.
>> 
>> 
>> On Aug 21, 2013, at 8:25 PM, Zhang <zyzhang at nuaa.edu.cn> wrote:
>> 
>>> 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);
>>> & nbsp; 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]))
>>>                             &nbs p;         - 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])
>>>          &nb sp;                            - 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 );
>>> 
>>> 
>> 
> 