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

[Matthew Knepley]

Okay, I understand now. You are correct that this formulation should be
available.
In fact, I think it should be possible to make this the bottom of the
hierarchy, with
the linear case a specialization, and the identity case a specialization of
that (with
a MF application). This is a good chance to redesign TS.

Next question. How should we do it? We can continue the discussion on email,
but
that makes it hard to refer back to pieces of code, or suggested interfaces.
I think
for a collaborative job, we need to augment the email discussion with a Wiki
of some
type (I was very impressed with Tim Gowers' polymath experiment).

Satish, do we have a Wiki, or could we setup one in short order?

  Thanks,

     Matt

On Thu, Mar 19, 2009 at 5:28 AM, Jed Brown <jed at 59a2.org> wrote:

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



-- 
What most experimenters take for granted before they begin their experiments
is infinitely more interesting than any results to which their experiments
lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20090319/05188041/attachment.html>