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

[Jed Brown]

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.

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 --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 835 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20131025/e46af289/attachment.pgp>