Okay, I understand now. You are correct that this formulation should be available.<br>In fact, I think it should be possible to make this the bottom of the hierarchy, with<br>the linear case a specialization, and the identity case a specialization of that (with<br>
a MF application). This is a good chance to redesign TS.<br><br>Next question. How should we do it? We can continue the discussion on email, but<br>that makes it hard to refer back to pieces of code, or suggested interfaces. I think<br>
for a collaborative job, we need to augment the email discussion with a Wiki of some<br>type (I was very impressed with Tim Gowers' polymath experiment).<br><br>Satish, do we have a Wiki, or could we setup one in short order?<br>
<br>  Thanks,<br><br>     Matt<br><br><div class="gmail_quote">On Thu, Mar 19, 2009 at 5:28 AM, Jed Brown <span dir="ltr"><<a href="mailto:jed@59a2.org">jed@59a2.org</a>></span> wrote:<br><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
<div class="im">On Wed 2009-03-18 18:34, Matthew Knepley wrote:<br>
> Not a whole lot. Needs to go slower, since I think there are a bunch<br>
> of embedded assumptions, and I worry that the most simple, general<br>
> interface cannot be seen until they are clarified. I want it spelled<br>
> out very very explicitly. So, to begin, you are unhappy with<br>
><br>
>   A_lhs(t) u_t = F(u,t)<br>
><br>
> where A_lhs is a linear operator. Is this true?<br>
<br>
</div>Yes, I agree that this would be the natural generalization of the<br>
existing design for linear problems.  I think it is not a good option<br>
for several reasons, most notably that it will often double the required<br>
storage and/or double the required computation per time-step when used<br>
with matrix-free methods.  Also see the end of this message for the case<br>
of general index-1 DAE where it is not possible to separate A_lhs from F.<br>
<div class="im"><br>
> For now, I will assume that A_lhs(t) is actually a nonlinear operator O_lhs<br>
> (u,u_t,t). Then<br>
> we have<br>
><br>
>   O_lhs(u,u_t,t) - F(u,t) = 0<br>
><br>
> So if we do Backward-Euler,<br>
><br>
>   O_lhs(u^n+1,t^n+1) - dt F(u^n+1,t^n+1) - O_lhs(u^n,t^n) = G = 0<br>
><br>
> so for Newton<br>
><br>
>   (O'_lhs - dt F') du = G<br>
><br>
> So are you asking for a nonlinear operator O?<br>
<br>
</div>Yes, so "shifted" function evaluation is<br>
<br>
  O = O_lhs - dt F<br>
<br>
and the "shifted" Jacobian is<br>
<br>
  O' = O'_lhs - dt F' .<br>
<br>
My proposal is that there are significant space- and time-savings if the<br>
user can supply O and O' for an arbitrary "shift" dt, instead of O_lhs,<br>
O'_lhs, F, F'.  For most problems, computing O and O' is essentially the<br>
same amount of work as F and F', and computing O_lhs, O'_lhs may require<br>
a similar amount of work (in my case) and often requires a similar<br>
amount of storage (the preconditioning matrix for O'_lhs often has the<br>
same nonzero pattern as F').  The purpose of O_lhs, O'_lhs are only to<br>
compute shifts and are never used alone (indeed cannot be used alone for<br>
DAE).<br>
<br>
Use of O, O' generalizes to high-order multistep and multistage implicit<br>
methods.<br>
<div class="im"><br>
> After I understand the answers to these, I can start to understand the<br>
> comments about shifts, etc. that talk about implementing this new<br>
> model.<br>
><br>
> I just want to clarify the movement away from the original TS model.<br>
<br>
</div>Yes, of course.<br>
<br>
My motivation comes from integrating DAE using a moving-mesh scheme<br>
where the spatial discretization does not assemble the real Jacobian,<br>
only a much sparser (factor 10--100) preconditioning matrix.  Evaluating<br>
the action of the Jacobian is less expensive than function evaluation<br>
only because transcendental functions don't need to be reevaluated, and<br>
much of the cost comes from evaluating tensor-product operations at<br>
quadrature points.  Keeping O_lhs and F separate means that application<br>
of O' (this is what is actually needed by the integrator) needs to be<br>
applied by applying O'_lhs and F' separately, essentially doubling the<br>
cost (this is also the case if O'_lhs is available via MFFD).  If O'_lhs<br>
and F' are assembled, it doubles the required storage and may double the<br>
assembly time.<br>
<br>
I think the advantages of keeping O_lhs separate is only realized when<br>
it is constant and F' is assembled, or when O_lhs is diagonal.<br>
<br>
CVODE (Sundials) requires that O_lhs is the identity (there is no way to<br>
change it) while IDA (the DAE solver) essentially asks the user for the<br>
shifted versions.  Indeed, for general index-1 DAE G(u_t,u,t)=0, it is<br>
not possible to separate O_lhs and F.  There is probably no advantage to<br>
restricting ourselves to the semi-explicit case.  Maybe it's clearer in<br>
the context of general index-1 DAE, equations of the form<br>
<br>
  F(u_t, u, t) = 0 .<br>
<br>
Solution using multistep (e.g. BDF) and multistage (Implicit<br>
Runge-Kutta) methods requires solving algebraic systems of the form<br>
<br>
  G(u) = F(known + alpha u, u, t)<br>
<br>
The required Jacobian is<br>
<br>
  G'(u)  = alpha dF/du_t + dF/du .<br>
<br>
These are analogous to the "shifted" operators I've been talking about,<br>
but apply in more general circumstances.<br>
<font color="#888888"><br>
Jed<br>
</font></blockquote></div><br><br clear="all"><br>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener<br>