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)
