TS: high-order with Galerkin spatial discretization, semi-explicit index-1 DAE

[Jed Brown]

On Wed 2009-03-18 18:34, Matthew Knepley wrote:
> Not a whole lot. Needs to go slower, since I think there are a bunch
> of embedded assumptions, and I worry that the most simple, general
> interface cannot be seen until they are clarified. I want it spelled
> out very very explicitly. So, to begin, you are unhappy with
> 
>   A_lhs(t) u_t = F(u,t)
> 
> where A_lhs is a linear operator. Is this true?

Yes, I agree that this would be the natural generalization of the
existing design for linear problems.  I think it is not a good option
for several reasons, most notably that it will often double the required
storage and/or double the required computation per time-step when used
with matrix-free methods.  Also see the end of this message for the case
of general index-1 DAE where it is not possible to separate A_lhs from F.

> For now, I will assume that A_lhs(t) is actually a nonlinear operator O_lhs
> (u,u_t,t). Then
> we have
> 
>   O_lhs(u,u_t,t) - F(u,t) = 0
> 
> So if we do Backward-Euler,
> 
>   O_lhs(u^n+1,t^n+1) - dt F(u^n+1,t^n+1) - O_lhs(u^n,t^n) = G = 0
> 
> so for Newton
> 
>   (O'_lhs - dt F') du = G
> 
> So are you asking for a nonlinear operator O?

Yes, so "shifted" function evaluation is

  O = O_lhs - dt F

and the "shifted" Jacobian is

  O' = O'_lhs - dt F' .

My proposal is that there are significant space- and time-savings if the
user can supply O and O' for an arbitrary "shift" dt, instead of O_lhs,
O'_lhs, F, F'.  For most problems, computing O and O' is essentially the
same amount of work as F and F', and computing O_lhs, O'_lhs may require
a similar amount of work (in my case) and often requires a similar
amount of storage (the preconditioning matrix for O'_lhs often has the
same nonzero pattern as F').  The purpose of O_lhs, O'_lhs are only to
compute shifts and are never used alone (indeed cannot be used alone for
DAE).

Use of O, O' generalizes to high-order multistep and multistage implicit
methods.

> After I understand the answers to these, I can start to understand the
> comments about shifts, etc. that talk about implementing this new
> model.
> 
> I just want to clarify the movement away from the original TS model.

Yes, of course.

My motivation comes from integrating DAE using a moving-mesh scheme
where the spatial discretization does not assemble the real Jacobian,
only a much sparser (factor 10--100) preconditioning matrix.  Evaluating
the action of the Jacobian is less expensive than function evaluation
only because transcendental functions don't need to be reevaluated, and
much of the cost comes from evaluating tensor-product operations at
quadrature points.  Keeping O_lhs and F separate means that application
of O' (this is what is actually needed by the integrator) needs to be
applied by applying O'_lhs and F' separately, essentially doubling the
cost (this is also the case if O'_lhs is available via MFFD).  If O'_lhs
and F' are assembled, it doubles the required storage and may double the
assembly time.

I think the advantages of keeping O_lhs separate is only realized when
it is constant and F' is assembled, or when O_lhs is diagonal.

CVODE (Sundials) requires that O_lhs is the identity (there is no way to
change it) while IDA (the DAE solver) essentially asks the user for the
shifted versions.  Indeed, for general index-1 DAE G(u_t,u,t)=0, it is
not possible to separate O_lhs and F.  There is probably no advantage to
restricting ourselves to the semi-explicit case.  Maybe it's clearer in
the context of general index-1 DAE, equations of the form

  F(u_t, u, t) = 0 .

Solution using multistep (e.g. BDF) and multistage (Implicit
Runge-Kutta) methods requires solving algebraic systems of the form

  G(u) = F(known + alpha u, u, t)

The required Jacobian is

  G'(u)  = alpha dF/du_t + dF/du .

These are analogous to the "shifted" operators I've been talking about,
but apply in more general circumstances.

Jed
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 197 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20090319/9c197403/attachment.sig>