[petsc-users] petsc4py DA with dof>1 confusion

[Lisandro Dalcin]

On 4 July 2011 15:17, Lisandro Dalcin <dalcinl at gmail.com> wrote:
> On 4 July 2011 14:25, Matthew Emmett <memmett at unc.edu> wrote:
>> Hi everyone,
>>
>> I am trying to use petsc4py in a fairly straight forward finite volume
>> shallow-water code (the PDE has three unknowns). I have a structured
>> grid and am using a DA to communicate ghost cells.  When I use dof=1,
>> everything works as I would expect.  However, when I try to use dof=3,
>> I am confused by the results.  I have probably mis-understood
>> something -- any suggestions would be appreciated.
>>
>> For example (see attached script), first I create the DA ('dof' and
>> 'width' are defined):
>>
>>  da = PETSc.DA().create(dim=2, sizes=[8, 8], dof=dof,
>>                         boundary_type='periodic',
>>                         stencil_width=width,
>>                         stencil_type='box')
>>
>> Next, I create global and local vectors:
>>
>>  gvec = da.createGlobalVec()
>>  lvec = da.createLocalVector()
>>
>> Then, put something into the local vectors, and scatter:
>>
>>  with da.getVecArray(lvec) as va:
>>    if dof == 1:
>>      va.array[width:-width,width:-width] = MPI.COMM_WORLD.rank+1
>>    else:
>>      va.array[width:-width,width:-width,0] = MPI.COMM_WORLD.rank+1
>>  da.localToGlobal(lvec, gvec)
>>
>> Finally, get the local vector again and print:
>>
>>  da.globalToLocal(gvec, lvec)
>>  with da.getVecArray(lvec) as va:
>>    if dof == 1:
>>      print va.array[:,:]
>>    else:
>>      print va.array[:,:,0]
>>
>> With 'dof=1' I get the following on the second processor:
>>
>> [[ 1.  2.  2.  2.  2.  1.]
>>  [ 1.  2.  2.  2.  2.  1.]
>>  [ 1.  2.  2.  2.  2.  1.]
>>  [ 1.  2.  2.  2.  2.  1.]
>>  [ 1.  2.  2.  2.  2.  1.]
>>  [ 1.  2.  2.  2.  2.  1.]
>>  [ 1.  2.  2.  2.  2.  1.]
>>  [ 1.  2.  2.  2.  2.  1.]
>>  [ 1.  2.  2.  2.  2.  1.]
>>  [ 1.  2.  2.  2.  2.  1.]]
>>
>> With 'dof=3' I get the following on the second processor:
>>
>> [[ 0.  0.  0.  0.  0.  0.]
>>  [ 0.  0.  0.  0.  2.  0.]
>>  [ 0.  0.  0.  0.  2.  0.]
>>  [ 0.  0.  0.  2.  2.  0.]
>>  [ 0.  0.  0.  2.  2.  0.]
>>  [ 0.  0.  0.  2.  2.  0.]
>>  [ 0.  0.  0.  2.  2.  0.]
>>  [ 0.  0.  0.  2.  2.  2.]
>>  [ 0.  0.  0.  2.  2.  2.]
>>  [ 0.  0.  0.  0.  0.  2.]]
>>
>> I was expecting to see the same as the 'dof=1' case, since I am
>> setting (and printing) va.array[:,:,0].
>>
>>
>> Again, any suggestions would be appreciated.  Thanks,
>> Matt
>>
>> PS, thanks to the developers of PETSc and petsc4py for all your hard
>> work.  I had the pleasure of meeting Jed and Matt at KAUST recently.
>>
>
> Mmm, it seems that this functionality is broken for the case of
> boundary_type='periodic'. I'll take a look.
>

No, sorry... It has nothing to do with periodicity... It is actually a
C/Fortran ordering issue I need to fix...

In Fortran 90, it seems you index a DA Vec array as A[dof,x,y,z]... ,
However, I think that for Python we should follow a more C-ish
indexing A[x,y,z,dof]. Or we could do it like PETSc in C, that is
A[z,y,x,dof] (wich is the transpose of the Fortran way) but it is
counter-intuitive to C (and likely Python) programmers ...

What do others think about this?




-- 
Lisandro Dalcin
---------------
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
3000 Santa Fe, Argentina
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169