[petsc-dev] 1D heat equation decay and implicit TS results

[Chetan Jhurani]

> -----Original Message-----
> From: Mark F. Adams
> Sent: Monday, October 01, 2012 2:21 PM
> To: For users of the development version of PETSc
> Subject: Re: [petsc-dev] 1D heat equation decay and implicit TS results
> 
> It looks like you are hitting the floating point limits someplace.
>
> The first thing to notice is that just because you ask for an rtol of 1.e-200 does not mean that you
> are going to get it.  Use -ksp_monitor to see what you are really getting.

I did some more investigation.  Also saw similar behavior when coded
this up in MATLAB for backward Euler.  It is an x-dependent
floating point error accumulation issue (not a tolerance issue or
a typical bug).  What was puzzling was the lack of chaotic behavior
that is usually present when the relative error reaches machine epsilon.
Instead we had a gradual change in slope.

Read the details and a fix if this sounds interesting.

The floating point error accumulation is x-dependent because the typical
linear solvers have a preferred DOF direction (left to right or right to
left in 1D).

What happens is that as iterations proceed, the original solution gradually
loses anti-symmetry (it had at t=0: sin(pi x) in [-1,1]).  One lobe starts
becoming smaller.  The kink happens when both lobes "merge" and the solution
becomes a curve with a single extremum.  This changes the character of
the convergence curve.  The same happens when SVD is used to solve
the implicit system. 

This also qualitatively matches with the facts that the new slope is approx
two times smaller than the "correct" slope and that the undesired shape
(accumulated over many time-steps) is in the near null-space of the matrix
being inverted.

If after each time-step one makes the solution explicitly anti-symmetric,
by forming u = 0.5*(u - u(end:-1:1)), then one can even reach denormals
(< 1E-308 ish) before the FEM solution starts deviating from the exact
solution.  Of course, this is a hack that is only suitable for some
specific initial functions like sin(pi x).  A more general way is to shift
u to maintain a zero average, which works here too.

Thanks for being an audience.

Chetan

> On Oct 1, 2012, at 5:06 PM, Chetan Jhurani <chetan.jhurani at gmail.com> wrote:
> 
> > Hi,
> >
> > I'm using petsc-dev to solve the 1D heat equation
> >
> > 	u_t = u_xx in [-1,1]
> >
> > with zero boundary conditions and u_initial = sin(pi x).
> > The exact solution decays to 0.  The exact
> > log10(norm_inf(u)) is a linear function of t.
> >
> > I use the standard 1D hat basis functions on a uniform
> > mesh to discretize.  The mass matrix remains tri-diagonal.
> >
> > When using TS with TSTHETA (or TSBEULER) and exact linear
> > solves (PCLU), the FEM log10(norm_inf(u)) follows the exact
> > curve for around 4500 time steps and then gradually,
> > over a few hundred time-steps, changes slope.
> >
> > See attached figure.  x axis is step number. The code is
> > attached too.
> >
> > The results were obtained using
> > -ksp_type preonly -pc_type lu -ksp_atol 1e-200 -ksp_rtol 1e-200
> >
> > The 1e-200 tolerances were added just to be way below the
> > norm scales in the problem above.
> >
> > Questions:
> >
> > Is it a bug in my code or is there another good explanation?
> >
> > Do I need to turn on some tolerance related option to see
> > the exact convergence behavior for more time steps?
> >
> > Chetan
> >
> > <ts_test.cpp><log_u_inf.png>