[petsc-users] 2nd order time-stepping in PETSc

[Jed Brown]

On Wed, 26 May 2010 21:48:56 -0300, Lisandro Dalcin <dalcinl at gmail.com> wrote:
> On 26 May 2010 18:58, Matthew Knepley <knepley at gmail.com> wrote:
> > It sounds like you could formulate it to solve for timestep n+1, and
> > then the RHS would depend on timestep n (which PETSc would give
> > you), and timestep n-1.  You could save t^{n-1} yourself with a
> > custom Monitor and get it in your RHS function through the user
> > context.     Matt
> >
> 
> What about using TSSet{Pre|Post}Step() ?

The problem is that this sort of ad-hoc dependence on the last step
causes the application to assume things about the integrator in use.
You can of course reduce the second-order problem to a first-order
system, but at the expense of increasing the length of the Krylov
vectors and ruining potential symmetry.  This is currently the only
general method supported by TS, but I don't like it.  An alternative is
to have the method perform differencing in auxiliary variables, as I
proposed a few weeks ago in a different context:

  http://lists.mcs.anl.gov/pipermail/petsc-dev/2010-April/002671.html

Some integrators, like RADAU, recognize special orderings and a flag for
this, but I find this especially awkward in parallel and we don't live
in an F77 world, so I would prefer to write this capability into the
API.

Specifically in this case, we have

  x_tt + f(x) = 0

which with time step \Delta t ~ 1/a and mass matrix M, yields a Jacobian

  J(x) = a^2 M + f'(x)

(I am ignoring the details of the time discretization which involves
precise definition of  a  and the right hand side.)  We can transform
this as

  x_t - y    = 0
  y_t + f(x) = 0

but we now have a Jacobian

  J(x,y) = [aM, -M; f'(x), aM]

which is somewhat unfortunate.  The proposal in the message referenced
above is not sufficient to transform this into the system you want, but
consider

  /* Set by the user before TSSolve(), called by the integrator at each
   * stage (maybe lazily). */
  typedef PetscErrorCode (*TSAuxiliaryFunction)(TS,Vec X,Vec X_t,Vec Y,void *ctx);
  PetscErrorCode TSSetAuxiliaryFunction(TS,TSAuxiliaryFunction,void*);

  /* Called by the user in TSIFunction(), Y is provided directly be the
   * user's function, Y_t is determined by the integrator using some
   * differencing of the Y values at prior steps/stages.  (This is a
   * virtual function to be defined by supported implementations.) */
  PetscErrorCode TSGetAuxiliaryVector(TS,Vec *Y,Vec *Y_t);

The material energy in the referenced message would by Y = Y(X), for
Hamiltonian systems you might have Y = Y(X_t) = X_t.  I don't think it
is possible to avoid leaving the user responsible for differentiating
through their auxiliary function when computing a Jacobian, but this is
easy for the easy cases (like the wave equation) and when it is actually
difficult, I think (naively) it is true complexity and not API
awkwardness.  What do you think?

Jed