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

[Christophe Ortiz]

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.



Christophe
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20131010/1f9a7213/attachment.html>