On Fri, Oct 25, 2013 at 10:04 PM, Jed Brown <[email protected]> wrote:
> Christophe Ortiz <[email protected]> writes: > > > 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 > > > > ierr = TSSetEquationType(ts,TS_EQ_DAE_IMPLICIT_INDEX1); > > Thanks for tracking this down. Indeed, that example enforces boundary > conditions using an algebraic expression: > > if (i == 0) { > f[i].u = hx * (x[i].u - user->uleft); > f[i].v = hx * (x[i].v - user->vleft); > > For this reason, we can only obtain the stage derivative by > differencing, rather than by evaluating the function. > > We need to audit the examples and use TSSetEquationType everywhere that > it's needed, which is probably in several of the examples. > Hope that it will help. But just to understand what I'm doing, what does TSSetEquationType do internally ? 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 _INDEXHI ? Christophe > > Note that Dirichlet boundary conditions can also be chosen by setting an > appropriate initial condition and using the equation du/dt = 0. > > > Well...the results is totally different, and in my opinion much more > > consistent. The value at the boundary is now the expected one: > > > > 0 10.000000 <----- > > 1 9.888647 > > 2 9.777382 > > 3 9.666293 > > 4 9.555467 > > > > ... > > ... > > 112 5.011530 > > 113 5.010665 > > 114 5.009880 > > 115 5.009169 > > ... > > ... > > 252 9.666293 > > 253 9.777382 > > 254 9.888647 > > 255 10.000000 <---- > > > > > > 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. > > > > What is the effect of the > TSSetEquationType(ts,TS_EQ_DAE_IMPLICIT_INDEX1); ? > > > > Christophe > > > > > > On Fri, Oct 25, 2013 at 11:47 AM, Christophe Ortiz < > > [email protected]> wrote: > > > >> Hi guys, > >> > >> 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. > >> > >> First, to have a better mesh I set 256 points (instead of the 11) in > >> > >> ierr = > >> DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,256,2,2,NULL,&da); > >> > >> Then, I did something simple. I removed reaction terms by commenting > >> > >> ierr = TSSetRHSFunction(ts,NULL,FormRHSFunction,&user);CHKERRQ(ierr); > >> > >> So we are left with 2 diffusion equations: > >> > >> u_t -alpha u_xx = 0 > >> v_t - alpha v_xx = 0 > >> > >> Now, do something simple: > >> 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. > >> > >> Now, for boundary conditions (Dirichlet here), set vleft and vright to > 10, > >> for instance. > >> > >> 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. > >> 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: > >> > >> 0 5.000000 > >> 1 8.951153 > >> 2 8.414078 > >> 3 7.964193 > >> 4 7.581196 > >> ... > >> ... > >> 77 5.000235 > >> 78 5.000207 > >> 79 5.000182 > >> 80 5.000161 > >> ... > >> ... > >> 252 7.964193 > >> 253 8.414078 > >> 254 8.951153 > >> 255 5.000000 > >> > >> 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. > >> Some output: > >> > >> TSAdapt 'basic': step 324 accepted t=9.73302 + 9.496e-02 > >> wlte=0.798 family='arkimex' scheme=0:'1bee' dt=9.569e-02 > >> 5.000000 > >> 10.000000 > >> 10.000000 > >> 10.000000 > >> 10.000000 > >> > >> TSAdapt 'basic': step 325 accepted t=9.82798 + 9.569e-02 > >> wlte=0.798 family='arkimex' scheme=0:'1bee' dt=9.642e-02 > >> 5.000000 > >> 10.000000 > >> 10.000000 > >> 10.000000 > >> 10.000000 > >> > >> 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. > >> > >> > >> I use petsc 3.4.1and configured it with > >> --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. > >> > >> The options I used to run the example are > >> -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. > >> > >> > >> Christophe > >> > >> > >> > >> > >> > >> > >> >
