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

[Matthew Knepley]

I could not understand the post. We will have to go more slowly. To begin,
I cannot understand why you cannot use TSSetMatrices()


http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-dev/docs/manualpages/TS/TSSetMatrices.html#TSSetMatrices

  Thanks,

     Matt

On Tue, Mar 17, 2009 at 9:15 PM, Jed Brown <jed at 59a2.org> wrote:

> I have a use for high-order implicit strong stability preserving
> integrators and plan to write a TS implementation based on some known
> optimal methods (cf. Gottlieb, Ketcheson, Shu 2009).  There are two
> cases of interest that I don't think TS can currently handle.
>
> Suppose we are using a Galerkin discretization so that the mass matrix
> is not the identity, but we still want to use high-order temporal
> discretization.  In particular, for high-order spatial discretizations,
> lumping the mass matrix is usually not acceptable so it would have dense
> inverse.  In these cases, MatShift is not sufficient [*].  One option is
> to allow shifts by an assembled mass matrix, but this is awkward with
> JFNK and requires significant storage when the mass matrix is not
> lumped.
>
> Perhaps a better option is for function evaluation and Jacobian assembly
> to observe the time-step (or a_{kk}*dt for Runge-Kutta methods).  This
> also works in the more general setting of semi-explicit index-1 DAE,
> presented in the form
>
>  V u_t = F(u,t)
>
> where F is the spatial discretization and V is typically block diagonal
> with mass matrices for differential variables and zero for algebraic
> variables.  To integrate such systems, the user would provide functions
> to evaluate
>
>  G(u,u_i,t,dt) = V u - F(dt*u+u_i,t)
>
> (u_i is a known inhomogeneous part) and assemble the
> Jacobian/preconditioner G'(u,t,dt).  This is sufficient for multistep
> and Diagonally Implicit Runge-Kutta (DIRK) methods.  For example, in the
> case of DIRK, stage k is computed by solving
>
>  G(u^k, u_0 + dt \sum_{l=0}^{k-1} a_{kl} u^l, t^k, dt a_{kk}) = 0
>
> for u^k.  For singly implicit Runge-Kutta (SIRK, a_{kk} constant)
> methods, lagging the preconditioner is perfectly acceptable since the
> time-step does not change frequently.  For JFNK with non-SIRK methods
> and with high-order multistep methods, the preconditioner would often
> have to be refactored to accommodate the shift, hence the inability to
> effectively lag the preconditioner is an unavoidable weakness.  Of
> course it is common that assembly is more expensive than PC setup in
> which case assembling F' and shifting by the mass matrix is appealing
> (assuming that storage for an assembled mass matrix is available).
>
> To summarize, asking the user to provide G,G' has the pros:
>  + consistent mass matrices with minimal cost
>  + DAE
>  + moving meshes
> and cons:
>  - not compatible with current F,F'
>  - shifts requiring reassembly are more expensive
>
> Asking the user to observe the time-step is at odds with the current
> interface where the user only provides F and F'.  Are the
> performance/flexibility benefits worth it?  How about a middle ground
> where the application can provide a lumped mass matrix (identity by
> default) and either functions F,F' or G,G' (where if G' is provided and
> the preconditioning matrix is different from the Jacobian, the
> preconditioner would be shifted using the lumped mass matrix)?
>
> I'd love to hear comments or suggestions.
>
>
> Jed
>
>
> [*] For my purposes, I don't actually need an assembled consistent mass
> matrix because my Jacobian is available via MatShell or MFFD (i.e. since
> shifting by the lumped mass matrix still produces a good preconditioner,
> the generalization of MatShift would be sufficient).  However, some
> people use the same matrix for the Jacobian and for preconditioning,
> hence a general scheme probably ought to be able to handle this.
>



-- 
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/20090318/b8e57f5e/attachment.html>