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

[Jed Brown]

On Thu 2009-03-19 08:04, Matthew Knepley wrote:
> 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.

Great.  I also prefer having the nonlinear DAE case be the bottom of the
hierarchy.  I think most specializations (and backward compatibility
with the current interface, if desired) could be achieved with thin
wrappers similar in spirit to SNESDAFormFunction.


The interface I have in mind is quite similar to IDA.  The DAE/ODE is
written in the form

  F(x_t,x,t) = 0

With an implicit system, the nonlinear equation that needs to be solved is

  G(x) = F(y + a x, x, t) = 0

(y is known and specific to the method) which has Jacobian

  J(x) = dG/dx = a dF/dx_t + dF/dx


The decision in IDA, translated to PETSc conventions, is that the user
provides F

  PetscErrorCode TSFormFunction_1(TS,Vec x_t,Vec x,PetscReal t,Vec f,void*);

and J(a,x_t,x,t)

  PetscErrorCode TSFormJacobian_1(TS,PetscReal a,Vec x_t,Vec x,PetscReal t,Mat jac,Mat jpre,MatStructure*,void*);

An alternative would be to provide G and G' as

  PetscErrorCode TSFormFunction_2(TS,PetscReal a,Vec y,Vec x,PetscReal t,Vec f,void*);
  PetscErrorCode TSFormJacobian_2(TS,PetscReal a,Vec y,Vec x,PetscReal t,Mat jac,Mat jpre,MatStructure*,void*);

in which case the SNES wrappers would not need to do the VecWAXPY
implied by interface 1.  I don't know if there is a compelling reason
not to choose interface 2.

Explicit methods can use the same function evaluation (calling with
a=0), but they need a wrapper when used with a mass matrix other than
the identity (which only makes sense if it's trivial to solve with,
e.g. block diagonal as with DG).  The wrapper would just turn

  M x_t = - F(0,x,t)

into

  x_t = G(x,t) = - M^{-1} F(0,x,t)



From reading a bit about Sundials, a useful operation for error
estimation is the "weighted RMS norm" defined as

  \sqrt{(1/N) \sum_{i=1}^N (x_i w_i)^2}

This may not make sense except for NORM_2, but perhaps it's worth

  PetscErrorCode VecWRMSNorm(Vec x,Vec w,NormType,PetscScalar*);

Thoughts?

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

Yes, a wiki would be handy.

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/20090320/4e32e80a/attachment.sig>