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