TS_RK enhancements

[Kai Germaschewski]

Hi everyone,

I've been using petsc for quite a while, and for a general disclaimer, I
think there are parts which I understand quite well, but it's also quite
possible that there's other functionality I totally missed, so I may end up
asking stupid questions at some point ;) Anyway, I think it's long overdue
to see whether I can contribute some pieces, if, of course, they sound like
a good idea.

Anyway, I'll start with something simple:

The attached patch fixes two issues (IMO) with TS_RK:

- the step numbering is slightly messed up: Before the timeloop starts,
TSMonitor is called with steps=0, then after the first step, it is called
again with steps=0, since steps is only incremented later. (Arguably, the
step numbering is more messed up than this because failed steps are counted,
but TSMonitor is not called. I would actually prefer to use ts->steps to
count the number of successful steps, and leave the number of failed steps
as an internal diagnostics only (shown in TSView()).

- TS_RK doesn't honor ts->max_steps. Whereas one could argue that setting a
number of steps doesn't make too much sense for a timestepper with variable
step size, it's actually useful for someone like me, see below.

- As a side effect, this patch also takes care of the warning "Very small
steps: 0.0000" that normally happened just as the final time is reached. (In
the next-to-last step, the step size would be set to be the difference
between final time and current, time, the last step would then update the
final time to current time and set the step size to 0, triggering the
warning).

Anyway, these two little changes should be rather uncontroversial ;)

The reason why I care about max_steps is that I'm using petsc to drive an
AMR code (not just TS_RK, I alternatively use TS_CN). Now this is pretty
much a much bigger beast and I guess it's open for discussion as to whether
petsc wants to add the complexity to support this better. In an AMR code,
timestepping as such shouldn't be that different in principle, but it really
does add complexity: One of the fundamental operations in AMR is to
refine/derefine the grid, which means that the state vector changes, and
with it all the operators / nonlinear function etc. The hack I'm currently
using is to use TS to integrate the system for a number of timesteps, then I
check whether regridding is necessary, and if so, I build a new state vector
etc, throw away the old TS and create a new one. As I'm in general in favor
of avoiding bloat unless it's well justified, I can't really see a good
reason to put all this logic into TS, unless there are other users who need
it.

It woulde be nice, though, to extend the TSSetUpdate() hook to indicate to
TS that it should break out of timestepping now. Then I could put the
regridding check into the TSUpdate(), and if regridding is necessary, have
it break out of the timestepping. One potential way of doing this without
adding yet another flag would be to break out if TSUpdate() sets dt to zero,
assuming there's no legitimate reason why an application would want to do
that currently. (or one could use a negative dt as a flag, but who knows
whether there might be a valid reason for a negative dt...)

If this sounds reasonable, I'd give it a shot.

As I've just now started to look at the code, it appears to me I may be
opening up the usual can of worms. Unless I'm missing it in some ungrepable
header, I don't see ops->update being called currently in the first place,
so I'd fix that, too. In addition, I see "reform()" and "reallocate()"
methods in TS, which sound like someone might have worked on something like
the more complex AMR scenario I described above, but they're never actually
used at all AFAICS.

More of that: prestep() and poststep() are actually implemented, but don't
seem to match the description as they're called before / after each
timestep, not once in the beginning / the end. I think the bug here is the
description, the names match watch they do. Except that really these
functions should have the *dt argument to set a timestep (prestep, at
least). update(), if it were implemented, would be called before each r.h.s.
evaluation which for RK would be each substep, and changing dt at that time
is likely to lead to garbage. Does anyone feel responsible for the TS
interface, or should I have a go at proposing something I consider sensible?

--Kai


-- 
Kai Germaschewski
Assistant Professor, Dept of Physics / Space Science Center
University of New Hampshire, Durham, NH 03824
office: Morse Hall 245F
phone:  +1-603-862-2912
fax: +1-603-862-2771
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20090511/0e0887d7/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: pp-ts_rk
Type: application/octet-stream
Size: 1590 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20090511/0e0887d7/attachment.obj>