TS_RK enhancements

Thu May 14 17:00:29 CDT 2009

Not sure if this is relevant to this discussion, but you could take a
view to my custom TS implementation for petsc4py. The idea is to being
able to do time-stepping regarding the time-step problem as a
nonlinear problem (think of the Jacobian being the sum of the
mass-matrix and the other linear/nonlinear contributions).

http://petsc.cs.iit.edu/petsc4py/petsc4py-dev/file/tip/src/source/src/ts/impls/python/python.c

In particular, see this:

    41 typedef struct _TSPyOps *TSPyOps;
    42
    43 struct _TSPyOps {
    44   PetscErrorCode (*presolve)  (TS);
    45   PetscErrorCode (*postsolve) (TS);
    46   PetscErrorCode (*prestep)   (TS,PetscReal,Vec);
    47   PetscErrorCode (*poststep)  (TS,PetscReal,Vec);
    48
    49   PetscErrorCode (*start)     (TS,PetscReal,Vec);
    50   PetscErrorCode (*step)      (TS,PetscReal,Vec);
    51   PetscErrorCode (*verify)    (TS,PetscReal,Vec,PetscTruth*,PetscReal*);
    52
    53   PetscErrorCode (*monitor)(TS,PetscInt,PetscReal,Vec);
    54
    55 };

The presolve/postsolve are ran only once for the entire timestepping,
prestep/poststep for each timestep, and the start/step/verify are
intended for implementing strategies for controlling the size of the
time-step.

On Thu, May 14, 2009 at 4:36 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>
>   Kai,
>
>    Thanks for the patch. I have applied it and pushed to petsc-dev
>>
>> 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?
>
>   There is no one who is coordinating the TS interface at the moment. So you
> are welcome to propose changes/fixes/improvements and are in as good as a
> position as any of us.
>
>   I see there are four methods: prestep, poststep, update, postupdate; as
> you noticed the naming and usage of them is a bit inconsistent. Even the
> TSStep() function has a confusing name since it actually takes (potentially)
> many steps.  I suggest starting the following possible changes.
>
>   change TSStep to TSSolve(),
>   move the prestep() and poststep() methods to be the ones used AT the
> beginning and end of each timestep (not called in TSSolve())
>   eliminate update() and postupdate() confusing names anyways
>   determine if there is a need for a pre and post called in TSSolve() could
> call them presolve and postsolve if you put them in.
>   do we really want prestep and poststep deciding to kick out of the
> iteration or do we want yet another method (corresponding to the converged()
> method in KSP and SNES) that is used to exit early from the solver?
>   I think using a flag to exit early from the solver is better than setting
> dt to zero.
>   having the steps and ptime output arguments from TSSolve() I think is
> wrong. They should be kept in the TS object and then calls like
> TSGetSolveTime() and TSGetNumberSteps()
>   to access them.
>
>  So some of the reorganization would be to make the use of TS match more
> closely SNES and KSP.
>
>   Barry
>
>
>
>
>
>
>
>>
>>
>> --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
>>
>> <pp-ts_rk>
>
>

-- 
Lisandro Dalcín
---------------
Centro Internacional de Métodos Computacionales en Ingeniería (CIMEC)
Instituto de Desarrollo Tecnológico para la Industria Química (INTEC)
Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET)
PTLC - Güemes 3450, (3000) Santa Fe, Argentina
Tel/Fax: +54-(0)342-451.1594