<div dir="ltr"><div class="gmail_extra"><br></div><div class="gmail_extra">After playing around, I found something interesting, an option that I have never seen in any example: TSSetEquationType. I did a TSGetEquationType and I got TS_EQ_UNSPECIFIED. Then I thought, ok, let's try to specify the equation type and I added the following line to the code</div>
<div class="gmail_extra"><br></div><div class="gmail_extra"> ierr = TSSetEquationType(ts,TS_EQ_DAE_IMPLICIT_INDEX1);</div><div class="gmail_extra"><br></div><div class="gmail_extra">Well...the results is totally different, and in my opinion much more consistent. The value at the boundary is now the expected one:</div>
<div class="gmail_extra"><br></div><div class="gmail_extra"><div class="gmail_extra">0 10.000000 <-----</div><div class="gmail_extra">1 9.888647 </div><div class="gmail_extra">2 9.777382 </div><div class="gmail_extra">
3 9.666293 </div><div class="gmail_extra">4 9.555467 </div><div class="gmail_extra"><br></div><div class="gmail_extra">...</div><div class="gmail_extra">...</div><div class="gmail_extra"><div class="gmail_extra">112 5.011530 </div>
<div class="gmail_extra">113 5.010665 </div><div class="gmail_extra">114 5.009880 </div><div class="gmail_extra">115 5.009169 </div><div>...</div><div>...</div><div>252 9.666293 </div><div>253 9.777382 </div>
<div>254 9.888647 </div><div>255 10.000000 <----</div><div><br></div></div><div><br></div></div><div class="gmail_extra">Moreover, I see much more effect of diffusion. Without this option, diffusion is always very weak, whatever the diffusion coefficient one puts, which always surprised me.</div>
<div class="gmail_extra"><br></div><div class="gmail_extra">What is the effect of the TSSetEquationType(ts,TS_EQ_DAE_IMPLICIT_INDEX1); ?</div><div class="gmail_extra"><br></div><div class="gmail_extra">Christophe<br>
<br><br><div class="gmail_quote">On Fri, Oct 25, 2013 at 11:47 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">Hi guys,
<div><br></div><div>I was having some troubles with my code (dof >> 1) and I decided to go back to simple things. Therefore, I took ex25.c, which is a nice case for me: diffusion of 2 species and 1 reaction.</div><div>
<br></div><div><div>First, to have a better mesh I set 256 points (instead of the 11) in </div><div><br></div><div> ierr = DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,256,2,2,NULL,&da);</div></div><div><br></div>
<div>Then, I did something simple. I removed reaction terms by commenting </div><div><br></div><div>ierr = TSSetRHSFunction(ts,NULL,FormRHSFunction,&user);CHKERRQ(ierr);<br></div><div><br></div><div>So we are left with 2 diffusion equations:</div>
<div><br></div><div>u_t -alpha u_xx = 0</div><div>v_t - alpha v_xx = 0</div><div><br></div><div>Now, do something simple:</div><div>In FormInitialSolution, for variable v, change it and simply set it to 5 (for instance), ie, v will have uniform values from i=0 to i=mx-1.</div>
<div><br></div><div>Now, for boundary conditions (Dirichlet here), set vleft and vright to 10, for instance.</div><div><br></div><div>Since values at the boundaries (10.0) are higher than values of initial solution (5.0), what you expect is some diffusion towards the center from the boundaries. Since we use Dirichlet boundary conditions here, you expect the values of v at i=0 and i=mx-1 to be 10.0. It can't be something else.</div>
<div>Well, that's not the case. I see that the value of v at i=0 and at i=mx-1 is always the initial value (t=0). Here is a sample of the final values after diffusion:</div><div><br></div><div><div>0 5.000000 </div>
<div>1 8.951153 </div><div>2 8.414078 </div><div>3 7.964193 </div><div>4 7.581196 </div></div><div>...</div><div>...</div><div><div>77 5.000235 </div><div>78 5.000207 </div><div>79 5.000182 </div><div>
80 5.000161 </div></div><div>...</div><div>...</div><div><div>252 7.964193 </div><div>253 8.414078 </div><div>254 8.951153 </div><div>255 5.000000</div></div><div><br></div><div>Then I checked the value of v inside IFunction at each timestep. I found out that each time IFunction is called, the first value of v that is used to evaluate F is the initial value (5.0), instead of of the value set by Dirichlet (10). For next evaluations within the same timestep, it is 10.0, the value imposed by Dirichlet boundary conditions. </div>
<div>Some output:</div><div><br></div><div><div> TSAdapt 'basic': step 324 accepted t=9.73302 + 9.496e-02 wlte=0.798 family='arkimex' scheme=0:'1bee' dt=9.569e-02 </div><div>5.000000</div>
<div>
10.000000</div><div>10.000000</div><div>10.000000</div><div>10.000000</div></div><div><br></div><div><div> TSAdapt 'basic': step 325 accepted t=9.82798 + 9.569e-02 wlte=0.798 family='arkimex' scheme=0:'1bee' dt=9.642e-02 </div>
<div>5.000000</div><div>10.000000</div><div>10.000000</div><div>10.000000</div><div>10.000000</div></div><div><br></div><div>Is it a bug or did I miss something ? I would expect that to evaluate F within IFunction, it would always use the value set by Dirichlet.</div>
<div><br></div><div><br></div><div>I use petsc 3.4.1and configured it with</div><div> --with-mpi=0 --with-fc=ifort --with-cc=icc --with-debugging=1 --with-scalar-type=real --with-precision=double --with-blas-lapack-dir=/opt/intel/mkl.</div>
<div><br></div><div>The options I used to run the example are</div><div>-ts_adapt_monitor -ts_adapt_basic_clip 0.1,1.1 -draw_pause -2 -ts_arkimex_type 1bee -ts_max_snes_failures -1 -snes_type newtonls -snes_linesearch_type bt.</div>
<span class=""><font color="#888888">
<div><br></div><div><br></div><div>Christophe</div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div></font></span></div>
</blockquote></div><br></div></div>