[petsc-users] Suggestions for code with dof >> 1 ?

[Christophe Ortiz]

On Thu, Oct 10, 2013 at 11:50 AM, Christophe Ortiz <
christophe.ortiz at ciemat.es> wrote:

>
> On Wed, Oct 9, 2013 at 10:51 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>
>>
>> On Oct 9, 2013, at 2:04 PM, Christophe Ortiz <christophe.ortiz at ciemat.es>
>> wrote:
>>
>> >
>> > On Wed, Oct 9, 2013 at 3:13 PM, Christophe Ortiz <
>> christophe.ortiz at ciemat.es> wrote:
>> >
>> > > For instance, if we take the example from petsc (MatSetValuesBlocked
>> page),
>> > > what should I do to pass the block of values 3, 4, 7 and 8 with
>> > > MatSetValuesBlocked ?
>> > >
>> > >   1  2  | 3  4
>> > >   5  6  | 7  8
>> > >     - - - | - - -
>> > >  9  10 | 11 12
>> > > 13 14 | 15 16
>> >
>> > Block row 0, block column 1.  (Count from 0.)
>> >
>> >
>> >
>> > And if with the previous case I would like to pass the values 1, 2, 3,
>> 4, 5, 6, 7 and 8 (the two first rows) what is the matrix val[ ] that I
>> should pass ? Should it be a matrix val[2][4] ?
>>
>>    It takes a one dimensional array as input; by default the values
>> should be ordered by row so you would pass in the numerical values 1, 2, 3,
>> 4,5, 6,7,8 in that order. If you wish to have the values ordered by column
>> then you would pass 1 5 2 6 3 7 4 8 after calling
>> MatSetOptions(mat,MAT_ROW_ORIENTED,PETSC_FALSE).
>>
>>    Note that the ordering you use on the numerical values is the same for
>> both MatSetValues() and MatSetValuesBlocked().
>>
>
>
> I am facing some problems with this when filling in the Jacobian for a
> problem of 2 diffusion equations (dof=2). See ex25.c.
>
> If I'm right, then for node i, there should be 3 blocks of 2x2:
>
> X  0  |  Y  0  |  X  0
> 0  X  |  0  Y  |  0  X
>
> Is it correct ?
>
> First I tried what you said, to pass a one dimensional array:
>    row = i;
>    cols[0] = i-1; cols[1] = i; cols[2] = i+1;
>
>    j = 0;
>    vallin[j] = X; vallin[j++] = 0; vallin[j++] = Y; vallin[j++] = 0;
>  vallin[j++] = X; vallin[j++] = 0;
>
>    vallin[j++] = 0; vallin[j++] = X; vallin[j++] = 0; vallin[j++] = Y;
> vallin[j++] = 0; vallin[j++] = X;
>
>   ierr = MatSetValuesBlocked(*Jpre,1,&row,3,cols,&vallin[0],INSERT_VALUES);
>
> I hit the following message:
>
> [0]PETSC ERROR: --------------------- Error Message
> ------------------------------------
> [0]PETSC ERROR: Floating point exception!
> [0]PETSC ERROR: Vec entry at local location 0 is not-a-number or infinite
> at end of function: Parameter number 3!
>
>
> Then, I tried a 2-dimensional array. I saw in some examples of
> ts/tutorials that it is also possible to pass a multi-dimensional arrays.
> After all, in memory the rows are stored contiguously, so I guess it is
> equivalent to a one-dimensional array.
>
>       j = 0; k = 0;
>       valdiff[j][k] = X; valdiff[j][k++] = 0;
>       valdiff[j][k++] = Y; valdiff[j][k++] = 0;
>       valdiff[j][k++] = X; valdiff[j][k++] = 0;
>
>       j= 1; k = 0;
>       valdiff[j][k] = 0; valdiff[j][k++] = X;
>       valdiff[j][k++] = 0; valdiff[j][k++] = Y;
>       valdiff[j][k++] = 0; valdiff[j][k++] = X;
>
>   ierr =
> MatSetValuesBlocked(*Jpre,1,&row,3,cols,&valdiff[0][0],INSERT_VALUES);
>
> I get the same error message.
> What did I do wrong ?
>
> In ex25.c they use a 3-dimensional array:
>
>  const PetscScalar vals[2][3][2] = {{{dxxL,0},{a *hx+dxx0,0},{dxxR,0}},
>
>  {{0,dxxL},{0,a*hx+dxx0},{0,dxxR}}};
>
> and &vals[0][0][0] is passed in MatSetValuesBlocked (....,....,
> &vals[0][0][0],...), and it works, though I don't understand how.
>
>
>
>
I found my mistake ! In the indices of the arrays I should have written ++j
and not j++. A good printf of each value helped me to find out...:-).

Now it works.
Christophe
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20131010/744be8da/attachment-0001.html>