TS grand plan

[Lisandro Dalcin]

On Tue, Sep 22, 2009 at 3:32 PM, Jed Brown <jed at 59a2.org> wrote:
> 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?
>

I need some time to review and think about it (just noticed your
commits in petsc-dev!), but the general idea you discussed in previous
posts was OK for me.

>> 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).
>

Indeed. I was "forced" to invent my TSPYTHON implementation just
because for the problems I have to solve at home (full implicit,
nonlinear, theta method), I had no other facilities in core PETSc. In
short, at every time step, I need to solve F(t,x,x')=0, so your
proposal is what I was always wanted to have in PETSc.

> 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).
>




-- 
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