<div dir="ltr">Hi,<div><br></div><div>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.</div>
<div><br></div><div>Which one do you recommend ?</div><div><br></div><div>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 ?</div>
<div><br></div><div>For instance, ex22.c</div><div><br></div><div><div>for (i=info.xs; i<info.xs+info.xm; i++) {</div><div>const PetscReal *k = user->k;</div><div>PetscScalar v[2][2];</div><div><br></div><div>v[0][0] = a + k[0]; v[0][1] = -k[1];</div>
<div>v[1][0] = -k[0]; v[1][1] = a+k[1];</div><div>MatSetValuesBlocked(*Jpre,1,&i,1,&i,&v[0][0],INSERT_VALUES);</div><div>}</div></div><div><br></div><div>Here 4 terms due to chemical reactions are passed to the matrix at once.</div>
<div>&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 ?</div><div><br>
</div><div><br></div><div>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 ?</div><div><br></div><div><div> 1 2 | 3 4</div>
<div> 5 6 | 7 8</div><div> - - - | - - -</div><div> 9 10 | 11 12</div><div>13 14 | 15 16</div></div><div><br></div><div><br></div><div>Christophe</div><div><br></div><div><br></div><div class="gmail_extra"><br><div class="gmail_quote">
On Tue, Oct 8, 2013 at 10:57 AM, Christophe Ortiz <span dir="ltr"><<a href="mailto:christophe.ortiz@ciemat.es" target="_blank">christophe.ortiz@ciemat.es</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
<div dir="ltr"><br><div class="gmail_extra"><div class="gmail_quote"><div class="im">On Mon, Oct 7, 2013 at 11:14 PM, Barry Smith <span dir="ltr"><<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div><br>
On Oct 7, 2013, at 2:49 PM, Christophe Ortiz <<a href="mailto:christophe.ortiz@ciemat.es" target="_blank">christophe.ortiz@ciemat.es</a>> wrote:<br>
<br>
> Thanks. I got it.<br>
><br>
> For the moment I have some general questions.<br>
><br>
> - 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 ?<br>
<br>
</div> RHS function is for non-stiff parts of the problem. IFunction for stiff parts.<br>
<br>
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.<br></blockquote><div><br></div></div><div>
I understand.</div><div class="im">
<div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
<div>><br>
> - 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,<br>
<br>
</div> 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.<br>
</blockquote><div><br></div></div><div>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 ?</div><div class="im">
<div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
<div>><br>
> - In this type of problem, which TS scheme is recommended ? TSARKIMEX ?<br>
<br>
</div> Beats me.<br>
<div><br>
><br>
> - The option TSARKIMEXSetFullyImplicit is used. Why ? Does it help ?<br>
<br>
</div> 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.<br></blockquote><div><br></div></div><div>It makes sense.</div>
<div>
<br></div><div>Thanks a lot. All this is clear now. I'll probably come back to you soon with more specific questions about ordering and matrices.</div><span class=""><font color="#888888"><div><br></div><div>Christophe</div>
</font></span><div><div class="h5"><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
<span><font color="#888888"><br>
Barry<br>
</font></span><div><div><br>
><br>
> Christophe<br>
><br>
> CIEMAT<br>
> Laboratorio Nacional de Fusión por Confinamiento Magnético<br>
> Unidad de Materiales<br>
> Edificio 2 - Planta 0 - Despacho 28m<br>
> Avenida Complutense 40,<br>
> 28040 Madrid, Spain<br>
> Tel: +34 91496 2582<br>
> Fax: +34 91346 6442<br>
><br>
> --<br>
> Q<br>
> Por favor, piense en el medio ambiente antes de imprimir este mensaje.<br>
> Please consider the environment before printing this email.<br>
><br>
><br>
> On Mon, Oct 7, 2013 at 3:23 PM, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>> wrote:<br>
><br>
> I have put it into the next branch. So just download petsc from the bitbucket site with git and then do git checkout next<br>
><br>
> The code won't work with the current PETSc release.<br>
><br>
> On Oct 7, 2013, at 8:16 AM, Christophe Ortiz <<a href="mailto:christophe.ortiz@ciemat.es" target="_blank">christophe.ortiz@ciemat.es</a>> wrote:<br>
><br>
> > Hi Barry,<br>
> ><br>
> > Thanks for the reply.<br>
> ><br>
> > 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.<br>
> ><br>
> > Could you please send me the latest version of ex10.c ? Many thanks in advance.<br>
> ><br>
> > Christophe<br>
> ><br>
> > CIEMAT<br>
> > Laboratorio Nacional de Fusión por Confinamiento Magnético<br>
> > Unidad de Materiales<br>
> > Edificio 2 - Planta 0 - Despacho 28m<br>
> > Avenida Complutense 40,<br>
> > 28040 Madrid, Spain<br>
> > Tel: +34 91496 2582<br>
> > Fax: +34 91346 6442<br>
> ><br>
> > --<br>
> > Q<br>
> > Por favor, piense en el medio ambiente antes de imprimir este mensaje.<br>
> > Please consider the environment before printing this email.<br>
> ><br>
> ><br>
> > On Mon, Oct 7, 2013 at 2:00 PM, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>> wrote:<br>
> ><br>
> > Chris,<br>
> ><br>
> > Take a look in the petsc branch barry/wirth-fusion-materials <a href="https://bitbucket.org/petsc/petsc/branch/barry%2Fwirth-fusion-materials" target="_blank">https://bitbucket.org/petsc/petsc/branch/barry%2Fwirth-fusion-materials</a><br>
> > in the file src/ts/examples/tutorials/advection-diffusion-reaction/ex10.c feel free to ask any questions.<br>
> ><br>
> > Barry<br>
> ><br>
> > On Oct 7, 2013, at 3:03 AM, Christophe Ortiz <<a href="mailto:christophe.ortiz@ciemat.es" target="_blank">christophe.ortiz@ciemat.es</a>> wrote:<br>
> ><br>
> > > Hi all<br>
> > ><br>
> > > I need some suggestions to design a code with PETSc.<br>
> > ><br>
> > > I want to solve a 1D problem composed of several diffusion equations and a lot of ODEs (thousands).<br>
> > ><br>
> > > - After discretization, the diffusion equations have terms in i-1, i and i+1.<br>
> > ><br>
> > > - The diffusion equations have additional terms due to couping with all the ODEs. These terms are non-linear.<br>
> > ><br>
> > > - The terms of the ODEs are local (only depend on node i) but are non-linear.<br>
> > ><br>
> > > Any advice to design the code ?<br>
> > ><br>
> > > How should I distribute the terms between IFunction and RHSFunction ?<br>
> > ><br>
> > > Any special attention to DMDA ? Should I declare just one DMDA with dof >>1 ?<br>
> > ><br>
> > > Many thanks in advance !<br>
> > > Christophe<br>
> > ><br>
> ><br>
> ><br>
><br>
><br>
<br>
</div></div></blockquote></div></div></div><br></div></div>
</blockquote></div><br></div></div>