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

[Christophe Ortiz]

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);

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/20131025/2881b446/attachment.html>