[petsc-users] use of VecGetArray() with dof>1 (block size > 1)

[Dave May]

Hi Ed,


On 28 May 2016 at 00:18, Ed Bueler <elbueler at alaska.edu> wrote:

> Dave --
>
> You are right that I had a cut-paste-edit error in my VecGetArray()
> invocation in the email.   Sorry about that.  I should of cut and pasted
> from my functioning, and DMDA-free, code.
>
> In this context I meant
>
> Vec v;
> ...  // create or load the Vec
> Node *u;
> VecGetArray(v,(PetscScalar**)&u);
>
> This *is* correct, and it works fine in the right context, without memory
> leaks.  I am *not* using a DMDA in this case.  At all.
>

Ah okay.


>
> My original quote from the PETSc User Manual should be read *as is*,
> however.  It *does* refer to DMDAVecGetArray().  And you are right that
> *it* has an error: should be [j][i] not [i][j].
>

Oh crap, that indexing error is in the manual (multiple times)!
Yikes.


> And DMDAVecGetArray() does a cast like the above, but internally.
>
> Again, my original question was about distinquishing/explaining the
> different types of the returned pointers from DMDAVecGetArray() and
> VecGetArray().
>

Okay.

VecGetArray() knows nothing about any DM. It only provides access to the
underlying entries in the vector which are PetscScalar's. So the last arg
is naturally PetscScalar**.

DMDAVecGet{Array,ArrayDOF}() exploits the DMDA ijk structure.
DMDAVecGet{Array,ArrayDOF}() exist solely for the convenience of the user
who wants to write an FD code, and wishes to express the discrete operators
with multi-dimensional arrays, e.g. which can be indexed like x[j+1][i-1].

DMDAVecGetArray() maps entries from a Vec to a users data type, which can
be indexed with the dimension of the DMDA (ijk). Since the dimension of the
DMDA can be 1,2,3 and the blocksize (e.g. number of members in your user
struct) is defined by the user - the last arg must be void*.

DMDAVecGetArrayDOF() maps entries from a Vec to multi-dimensional array
indexed by the dimension of the DMDA AND the number of fields  (defined via
the block size). Since the dimension of the DMDA can be 1,2 or 3, again,
the last arg must be void* unless a separate API is introduced for 1D, 2D
and 3D.

Why do both exist? Well, Jed provided one reasons - the block size may be a
runtime choice and thus the definition of the struct cannot be changed at
runtime. Another reason could be the user just doesn't think it is useful
to attach names (i.e. though members in a struct) to their DOFs - hence
they want DMDAVecGetArrayDOF(). This could arise if you used the DMDA to
describe a DG discretization. The DOFs would then just represent
coefficients associated with your basis function. Maybe the user just
prefers to write out loops.

I see DMDAVec{Get,Restore}XXX as tools to help the user.
The user can pick the API they prefer.

I use the DMDA, but I always use plain old VecGetArray() and obtain the
result in a variable declared as PetscScalar*. I don't bother with mapping
the entries into a struct. I see no advantage to this in terms of code
clarity. I prefer not to use DMDAVecGetArray() and DMDAVecGetArrayDOF() as
these methods allocate additional memory.

In my opinion, the line you refer to from the manual regarding
multi-component PDEs should only applied in the context of usage with the
DMDA. Others may disagree.

I hope I finally helped answered your question.


Cheers,
  Dave




>
> Ed
>
>
> On Fri, May 27, 2016 at 2:47 PM, Dave May <dave.mayhem23 at gmail.com> wrote:
>
>>
>>
>> On 27 May 2016 at 21:24, Ed Bueler <elbueler at alaska.edu> wrote:
>>
>>> Dave --
>>>
>>> Perhaps you should re-read my questions.
>>>
>>
>> Actually - maybe we got our wires crossed from the beginning.
>> I'm going back to the original email as I missed something.
>>
>>
>>>> """
>>>> The recommended approach for multi-component PDEs is to declare a
>>>> struct representing the fields defined
>>>> at each node of the grid, e.g.
>>>>
>>>> typedef struct {
>>>> PetscScalar u,v,omega,temperature;
>>>> } Node;
>>>>
>>>> and write residual evaluation using
>>>>
>>>> Node **f,**u;
>>>> DMDAVecGetArray(DM da,Vec local,&u);
>>>> DMDAVecGetArray(DM da,Vec global,&f);
>>>> ...
>>>>
>>>>
>>
>> 1)  The third argument to DMDAVec{Get,Restore}Array() is of type "void
>>>> *".  It makes the above convenient.  But the third argument of the
>>>> unstructured version Vec{Get,Restore}Array() is of type "PetscScalar **",
>>>> which means that in an unstructured case, with the same Node struct, I
>>>> would write
>>>> "VecGetArray(DM da,Vec local,(PetscScalar **)&u);"
>>>> to get the same functionality.  Why is it this way?  More specifically,
>>>> why not have the argument to VecGetArray() be of type "void *"?
>>>>
>>>
>>>
>> Is the quoted text
>>   "VecGetArray(DM da,Vec local,(PetscScalar **)&u);"
>> really what you meant?
>>
>> Sorry I didn't spot this on the first read, but probably you meant
>> something else as VecGetArray() only takes two args (Vec,PetscScalar**).
>>
>> This code
>>
>>   Node **u;
>>   VecGetArray(Vec local,(PetscScalar**)&u);
>>
>> would not be correct, neither would
>>
>>   Node ***u;
>>   VecGetArray(Vec local,(PetscScalar**)&u);
>> if the DMDA was defined in 3d.
>>
>>
>>
>>>
>>> I would say the reason why the last arg to VecGetArray() is not void* is
>>> because it is intended to give you direct access to the pointer associated
>>> with the entries within the vector - these are also PetscScalar's
>>>
>>>
>>>>
>>>> 2) Given that the "recommended approach"
>>>>
>>>
>>> I don't believe it is ever recommended anywhere to do the following:
>>>   VecGetArray(DM da,Vec local,(PetscScalar**)&u)
>>>
>>> Trying to trick the compile with such a cast is just begging for memory
>>> corruption to occur.
>>>
>>> above works just fine, why do DMDAVec{Get,Restore}ArrayDOF() exist?
>>>>  (I.e. is there something I am missing about C indexing?)
>>>>
>>>
>>> As an additional point, DMDAVec{Get,Restore}ArrayDOF() return
>>>
>>> void *array
>>>
>>> so that the same API will work for 1D, 2D and 3D DMDA's which would
>>> require PetscScalar **data, PetscScalar ***data, PetscScalar ****data
>>> respectively.
>>>
>>>
>>> Cheers,
>>>   Dave
>>>
>>>
>>>> 3) There are parts of the PETSc API that refer to "dof" and parts that
>>>> refer to "block size".  Is this a systematic distinction with an underlying
>>>> reason?  It seems "block size" is more generic, but also it seems that it
>>>> could replace "dof" everywhere.
>>>>
>>>> Thanks for addressing silly questions.
>>>>
>>>> Ed
>>>>
>>>>
>>>> --
>>>> Ed Bueler
>>>> Dept of Math and Stat and Geophysical Institute
>>>> University of Alaska Fairbanks
>>>> Fairbanks, AK 99775-6660
>>>> 301C Chapman and 410D Elvey
>>>> 907 474-7693 and 907 474-7199  (fax 907 474-5394)
>>>>
>>>
>>>
>>
>>
>> --
>> Ed Bueler
>> Dept of Math and Stat and Geophysical Institute
>> University of Alaska Fairbanks
>> Fairbanks, AK 99775-6660
>> 301C Chapman and 410D Elvey
>> 907 474-7693 and 907 474-7199  (fax 907 474-5394)
>>
>>
>
>
> --
> Ed Bueler
> Dept of Math and Stat and Geophysical Institute
> University of Alaska Fairbanks
> Fairbanks, AK 99775-6660
> 301C Chapman and 410D Elvey
> 907 474-7693 and 907 474-7199  (fax 907 474-5394)
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160528/fffefa09/attachment.html>