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

Wed Oct 9 07:50:56 CDT 2013

Hi,

I am building my Jacobian (IJacobian) with dof > 1 and I am considering to
use either MatSetValuesLocal (what you did in ex10.c from
advection-diffusion) or MatSetValuesBlocked.

Which one do you recommend ?

For the moment I am trying to understand how MatSetValuesBlocked works from
the various examples given in ts tutorials, but there is a mistery for
me...In these examples, the block size is not set. So when m and idxm[] (n
and idxn[]) are passed in, how does PETSc know the correspondance with
global row and column indices in the matrix ?

For instance, ex22.c

for (i=info.xs; i<info.xs+info.xm; i++) {
const PetscReal *k = user->k;
PetscScalar     v[2][2];

v[0][0] = a + k[0]; v[0][1] =  -k[1];
v[1][0] =    -k[0]; v[1][1] = a+k[1];
MatSetValuesBlocked(*Jpre,1,&i,1,&i,&v[0][0],INSERT_VALUES);
}

Here 4 terms due to chemical reactions are passed to the matrix at once.
&i is passed for indice of row and column, which is supposed to be the
indice of the block. But if the blocksize is not set, how PETSc knows to
which row/colum it really corresponds in the matrix ?

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

Christophe

On Tue, Oct 8, 2013 at 10:57 AM, Christophe Ortiz <
christophe.ortiz at ciemat.es> wrote:

>
> On Mon, Oct 7, 2013 at 11:14 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>
>>
>> On Oct 7, 2013, at 2:49 PM, Christophe Ortiz <christophe.ortiz at ciemat.es>
>> wrote:
>>
>> > Thanks. I got it.
>> >
>> > For the moment I have some general questions.
>> >
>> > - I see that in this example all terms--diffusion and reactions--are
>> put into RHSFunction. There is no IFunction. Is there a reason for that ?
>>  In which cases one should use IFunction or RHSFunction ? In ts/.../ex25.c,
>> diffusion terms are defined in IFunction and reaction terms in RHSFunction.
>> Why not here ?
>>
>>    RHS function is for non-stiff parts of the problem. IFunction for
>> stiff parts.
>>
>>    There reason they are not separated in ex10.c is because I don't know
>> for this problem what is what, so I just put them all in RHSFunction. But
>> use a fully implicit solver.
>>
>
> I understand.
>
>
>> >
>> > - If only RHSFunction is used, how is it possible to define boundary
>> conditions ? For instance, a Dirichlet boundary condition, u(i=0)=3. In
>> ts/.../ex2.c, only RHSFunction is used and the boundary conditions are
>> defined through u_t,
>>
>>    If, for example u is zero on the boundary for all time then you can
>> just set u(t=0) to 0 on the boundary and u_t = 0 on the boundary and the
>> boundary values of u will remain zero.  If you want to do something like u
>> = f(t) on the boundary you can put that into the IFunction.
>>
>
> I see. So I could put either everything in the IFunction, or everything in
> RHS but the two boundaries that cannot be defined in RHS (u=f(t)), in
> IFunction. Correct ?
>
>
>> >
>> > - In this type of problem, which TS scheme is recommended ? TSARKIMEX ?
>>
>>    Beats me.
>>
>> >
>> > - The option TSARKIMEXSetFullyImplicit is used. Why ? Does it help ?
>>
>>     I just wanted to use a fully implicit solver; since I just used
>> RHSFunction by default TSARKIMEX would treat it as explicit which I didn't
>> want.
>>
>
> It makes sense.
>
> Thanks a lot. All this is clear now. I'll probably come back to you soon
> with more specific questions about ordering and matrices.
>
> Christophe
>
>
>>
>>    Barry
>>
>> >
>> > Christophe
>> >
>> > CIEMAT
>> > Laboratorio Nacional de Fusión por Confinamiento Magnético
>> > Unidad de Materiales
>> > Edificio 2 - Planta 0 - Despacho 28m
>> > Avenida Complutense 40,
>> > 28040 Madrid, Spain
>> > Tel: +34 91496 2582
>> > Fax: +34 91346 6442
>> >
>> > --
>> > Q
>> > Por favor, piense en el medio ambiente antes de imprimir este mensaje.
>> > Please consider the environment before printing this email.
>> >
>> >
>> > On Mon, Oct 7, 2013 at 3:23 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>> >
>> >   I have put it into the next branch. So just download petsc from the
>> bitbucket site with git and then do git checkout next
>> >
>> >    The code won't work with the current PETSc release.
>> >
>> > On Oct 7, 2013, at 8:16 AM, Christophe Ortiz <
>> christophe.ortiz at ciemat.es> wrote:
>> >
>> > > Hi Barry,
>> > >
>> > > Thanks for the reply.
>> > >
>> > > I had a look at the example ex10.c. But I see it does not fully
>> corresponds to what one can see in the branch. For instance, in the example
>> that is found in the tutorials, there is no RHS defined, no Jacobian, just
>> the IFunction.
>> > >
>> > > Could you please send me the latest version of ex10.c ? Many thanks
>> in advance.
>> > >
>> > > Christophe
>> > >
>> > > CIEMAT
>> > > Laboratorio Nacional de Fusión por Confinamiento Magnético
>> > > Unidad de Materiales
>> > > Edificio 2 - Planta 0 - Despacho 28m
>> > > Avenida Complutense 40,
>> > > 28040 Madrid, Spain
>> > > Tel: +34 91496 2582
>> > > Fax: +34 91346 6442
>> > >
>> > > --
>> > > Q
>> > > Por favor, piense en el medio ambiente antes de imprimir este mensaje.
>> > > Please consider the environment before printing this email.
>> > >
>> > >
>> > > On Mon, Oct 7, 2013 at 2:00 PM, Barry Smith <bsmith at mcs.anl.gov>
>> wrote:
>> > >
>> > >   Chris,
>> > >
>> > >     Take a look in the petsc branch barry/wirth-fusion-materials
>> https://bitbucket.org/petsc/petsc/branch/barry%2Fwirth-fusion-materials
>> > > in the file
>> src/ts/examples/tutorials/advection-diffusion-reaction/ex10.c   feel free
>> to ask any questions.
>> > >
>> > >     Barry
>> > >
>> > > On Oct 7, 2013, at 3:03 AM, Christophe Ortiz <
>> christophe.ortiz at ciemat.es> wrote:
>> > >
>> > > > Hi all
>> > > >
>> > > > I need some suggestions to design a code with PETSc.
>> > > >
>> > > > I want to solve a 1D problem composed of several diffusion
>> equations and a lot of ODEs (thousands).
>> > > >
>> > > > - After discretization, the diffusion equations have terms in i-1,
>> i and i+1.
>> > > >
>> > > > - The diffusion equations have additional terms due to couping with
>> all the ODEs. These terms are non-linear.
>> > > >
>> > > > - The terms of the ODEs are local (only depend on node i) but are
>> non-linear.
>> > > >
>> > > > Any advice to design the code ?
>> > > >
>> > > > How should I distribute the terms between IFunction and RHSFunction
>> ?
>> > > >
>> > > > Any special attention to DMDA ? Should I declare just one DMDA with
>> dof >>1 ?
>> > > >
>> > > > Many thanks in advance !
>> > > > Christophe
>> > > >
>> > >
>> > >
>> >
>> >
>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20131009/3b480a8a/attachment-0001.html>