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
