<div dir="ltr"><br><div class="gmail_extra"><div class="gmail_quote">On Fri, Oct 25, 2013 at 10:04 PM, Jed Brown <span dir="ltr"><<a href="mailto:jedbrown@mcs.anl.gov" target="_blank">jedbrown@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>Christophe Ortiz <<a href="mailto:christophe.ortiz@ciemat.es" target="_blank">christophe.ortiz@ciemat.es</a>> writes:<br>
<br>
> After playing around, I found something interesting, an option that I have<br>
> never seen in any example: TSSetEquationType. I did a TSGetEquationType and<br>
> I got TS_EQ_UNSPECIFIED. Then I thought, ok, let's try to specify the<br>
> equation type and I added the following line to the code<br>
><br>
> ierr = TSSetEquationType(ts,TS_EQ_DAE_IMPLICIT_INDEX1);<br>
<br>
</div>Thanks for tracking this down. Indeed, that example enforces boundary<br>
conditions using an algebraic expression:<br>
<br>
if (i == 0) {<br>
f[i].u = hx * (x[i].u - user->uleft);<br>
f[i].v = hx * (x[i].v - user->vleft);<br>
<br>
For this reason, we can only obtain the stage derivative by<br>
differencing, rather than by evaluating the function.<br>
<br>
We need to audit the examples and use TSSetEquationType everywhere that<br>
it's needed, which is probably in several of the examples.<br></blockquote><div><br></div><div>Hope that it will help.</div><div><br></div><div>But just to understand what I'm doing, what does TSSetEquationType do internally ? </div>
<div><br></div>
<div>For the type of problem I'm trying to solve, is TS_EQ_DAE_IMPLICIT_INDEX1 ok ? What is the difference with TS_EQ_DAE_IMPLICIT_INDEX2,3 or _<span style="color:rgb(0,0,0)">INDEXHI ?</span></div><div><span style="color:rgb(0,0,0)"><br>
</span></div><div><span style="color:rgb(0,0,0)">Christophe</span></div><div><br></div><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">
<br>
Note that Dirichlet boundary conditions can also be chosen by setting an<br>
appropriate initial condition and using the equation du/dt = 0.<br>
<div><div><br>
> Well...the results is totally different, and in my opinion much more<br>
> consistent. The value at the boundary is now the expected one:<br>
><br>
> 0 10.000000 <-----<br>
> 1 9.888647<br>
> 2 9.777382<br>
> 3 9.666293<br>
> 4 9.555467<br>
><br>
> ...<br>
> ...<br>
> 112 5.011530<br>
> 113 5.010665<br>
> 114 5.009880<br>
> 115 5.009169<br>
> ...<br>
> ...<br>
> 252 9.666293<br>
> 253 9.777382<br>
> 254 9.888647<br>
> 255 10.000000 <----<br>
><br>
><br>
> Moreover, I see much more effect of diffusion. Without this option,<br>
> diffusion is always very weak, whatever the diffusion coefficient one puts,<br>
> which always surprised me.<br>
><br>
> What is the effect of the TSSetEquationType(ts,TS_EQ_DAE_IMPLICIT_INDEX1); ?<br>
><br>
> Christophe<br>
><br>
><br>
> On Fri, Oct 25, 2013 at 11:47 AM, Christophe Ortiz <<br>
> <a href="mailto:christophe.ortiz@ciemat.es" target="_blank">christophe.ortiz@ciemat.es</a>> wrote:<br>
><br>
>> Hi guys,<br>
>><br>
>> I was having some troubles with my code (dof >> 1) and I decided to go<br>
>> back to simple things. Therefore, I took ex25.c, which is a nice case for<br>
>> me: diffusion of 2 species and 1 reaction.<br>
>><br>
>> First, to have a better mesh I set 256 points (instead of the 11) in<br>
>><br>
>> ierr =<br>
>> DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,256,2,2,NULL,&da);<br>
>><br>
>> Then, I did something simple. I removed reaction terms by commenting<br>
>><br>
>> ierr = TSSetRHSFunction(ts,NULL,FormRHSFunction,&user);CHKERRQ(ierr);<br>
>><br>
>> So we are left with 2 diffusion equations:<br>
>><br>
>> u_t -alpha u_xx = 0<br>
>> v_t - alpha v_xx = 0<br>
>><br>
>> Now, do something simple:<br>
>> In FormInitialSolution, for variable v, change it and simply set it to 5<br>
>> (for instance), ie, v will have uniform values from i=0 to i=mx-1.<br>
>><br>
>> Now, for boundary conditions (Dirichlet here), set vleft and vright to 10,<br>
>> for instance.<br>
>><br>
>> Since values at the boundaries (10.0) are higher than values of initial<br>
>> solution (5.0), what you expect is some diffusion towards the center from<br>
>> the boundaries. Since we use Dirichlet boundary conditions here, you expect<br>
>> the values of v at i=0 and i=mx-1 to be 10.0. It can't be something else.<br>
>> Well, that's not the case. I see that the value of v at i=0 and at i=mx-1<br>
>> is always the initial value (t=0). Here is a sample of the final values<br>
>> after diffusion:<br>
>><br>
>> 0 5.000000<br>
>> 1 8.951153<br>
>> 2 8.414078<br>
>> 3 7.964193<br>
>> 4 7.581196<br>
>> ...<br>
>> ...<br>
>> 77 5.000235<br>
>> 78 5.000207<br>
>> 79 5.000182<br>
>> 80 5.000161<br>
>> ...<br>
>> ...<br>
>> 252 7.964193<br>
>> 253 8.414078<br>
>> 254 8.951153<br>
>> 255 5.000000<br>
>><br>
>> Then I checked the value of v inside IFunction at each timestep. I found<br>
>> out that each time IFunction is called, the first value of v that is used<br>
>> to evaluate F is the initial value (5.0), instead of of the value set by<br>
>> Dirichlet (10). For next evaluations within the same timestep, it is 10.0,<br>
>> the value imposed by Dirichlet boundary conditions.<br>
>> Some output:<br>
>><br>
>> TSAdapt 'basic': step 324 accepted t=9.73302 + 9.496e-02<br>
>> wlte=0.798 family='arkimex' scheme=0:'1bee' dt=9.569e-02<br>
>> 5.000000<br>
>> 10.000000<br>
>> 10.000000<br>
>> 10.000000<br>
>> 10.000000<br>
>><br>
>> TSAdapt 'basic': step 325 accepted t=9.82798 + 9.569e-02<br>
>> wlte=0.798 family='arkimex' scheme=0:'1bee' dt=9.642e-02<br>
>> 5.000000<br>
>> 10.000000<br>
>> 10.000000<br>
>> 10.000000<br>
>> 10.000000<br>
>><br>
>> Is it a bug or did I miss something ? I would expect that to evaluate F<br>
>> within IFunction, it would always use the value set by Dirichlet.<br>
>><br>
>><br>
>> I use petsc 3.4.1and configured it with<br>
>> --with-mpi=0 --with-fc=ifort --with-cc=icc --with-debugging=1<br>
>> --with-scalar-type=real --with-precision=double<br>
>> --with-blas-lapack-dir=/opt/intel/mkl.<br>
>><br>
>> The options I used to run the example are<br>
>> -ts_adapt_monitor -ts_adapt_basic_clip 0.1,1.1 -draw_pause -2<br>
>> -ts_arkimex_type 1bee -ts_max_snes_failures -1 -snes_type newtonls<br>
>> -snes_linesearch_type bt.<br>
>><br>
>><br>
>> Christophe<br>
>><br>
>><br>
>><br>
>><br>
>><br>
>><br>
>><br>
</div></div></blockquote></div><br></div></div>