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

[Christophe Ortiz]

On Thu, Oct 10, 2013 at 2:33 PM, Jed Brown <jedbrown at mcs.anl.gov> wrote:

> Christophe Ortiz <christophe.ortiz at ciemat.es> writes:
> > 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;
>
> The expression j++ in vallin[j++] is post-increment, so the two
> statements above both write to vallin[0] and vallin[11] is undefined
> (whatever garbage was in the array before).
>

Thanks.
I just found out few minutes ago :-). Hehe. The point is that I come from
15 years of Fortran...I have to refresh the C that I learnt at Faculty and
that is deeply buried in my memory !



>
> > 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;
>
> Same problem here.
>
> >       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.
>
> That grouping in the 3D array is the same ordering, but (perhaps) makes
> it easier to see which expression refers to which entry of the matrix.
>

Yes, probably.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20131010/f36be5a4/attachment.html>