TS grand plan

[Jed Brown]

Lisandro Dalcin wrote:
> Just for reference, in petsc4py (accesible from core PETSc as TSPYTHON
> type when configured --with-petsc4py) more or less uses Jed's model:
> the user routines are in charge of computing F(...)=0 and the
> Jacobian. Moreover, this code can also manage some basic time-step
> control.

Neat.  What do you think of using my IFunction/IJacobian interfaces?

> BTW, I proposed something similar to this in the past (look for
> subject "flexible TS implementation for user-defined timestepping").

Interesting thread.  So the established TS interface is to write
x'=f(t,x).  Matt notes that this is screwy for Galerkin schemes because
(naively) you would have to write f(t,x) = M^{-1}g(t,x) which usually
has a dense Jacobian and is really only suitable for explicit schemes
(and inverting the mass matrix is part of the reason why FEM is
considered to be expensive for such problems).  A hack for FEM is to let
the user provide M, but then you need to build a*M+F' which might now
have a different sparsity pattern than either M,F' and you would need to
handle extra terms if M is changing (as in ALE).

I put in an alternative with the implicit formulation f(t,x,x')=0.  The
user provides (a preconditioner for) the required Jacobian which is
g'(x) where g(x)=f(t,x,x0+a*x) for some t,a,x0 determined by the
integrator.  All general linear schemes (which includes fully- and
diagonally-implicit Runge-Kutta, as well as multistep and, with slight
abuse, Rosenbrock) can get exactly what they need with this interface.
It is immediately applicable to index-1 DAE; higher index systems may
require more but it's very common to perform index reduction anyway.
Note that (discontinuous) Galerkin time-discretization can, after choice
of quadrature, be recast as a GL scheme (usually RK).

So I think we have no lack of methods that can be represented
algebraically as long as the user is happy to use an implicit method
(perhaps with their "semi-implicit" method in a preconditioner).  Time
splitting schemes are much harder to abstract, and the details of the
splitting scheme usually effects the boundary conditions and perhaps the
equations to be solved (otherwise they are more likely to produce stable
and reasonable looking garbage including very wrong steady states and
nonphysical cycles).


Jed

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 261 bytes
Desc: OpenPGP digital signature
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20090922/8ea37cfd/attachment.sig>