[petsc-users] Bug in boundary conditions (ex25.c) ?

[Christophe Ortiz]

On Fri, Oct 25, 2013 at 10:04 PM, Jed Brown <jedbrown at mcs.anl.gov> wrote:

> Christophe Ortiz <christophe.ortiz at ciemat.es> 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 <
> > christophe.ortiz at ciemat.es> 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
> >>
> >>
> >>
> >>
> >>
> >>
> >>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20131028/1cc8e3a1/attachment.html>