[petsc-dev] Fw: TSTheta for DAE's

[Jed Brown]

On Tue, Dec 21, 2010 at 20:07, Michael E Henderson <mhender at us.ibm.com>wrote:

> If my initial condition does not satisfy the algebraic constraints the
> difference between it and th->X, which is used to estimate the "time"
> derivative is not correct, and so the solution after the VecAXPY which adds
> dt*th->Xdot to it also does not satisfy the implicit eqs. This continues,
> and we see oscillations in the algebraic variables.


Whether or not we place the monitor where these oscillations are visible,
they _are_ being propagated.  Theta=0.5 is not a good method for DAE, just
because the oscillations vanish at the quadrature point does not mean that
they are not present and contaminating the solution.


> On the otherhand, th->X after SNES does satisfy the algebraic equations.
> I'm not concerned with the traditional type of instability so much, just
> that the algebraic constraints be satisfied.
>
> For theta=1/2, x_0 is input, and
>
> f(x_1-x_0)/dt,(x_1+x_0)/2)=0 is solved for x_1.
>
> Then monitor is passed x_1+(x_1-x_0)/dt *dt.    If my algebraic constraint
> is x[0]-1=0, and x_0[0]=0., SNES should return x_1[0]=1., so monitor sees
>
> 1+(1-0)/dt*dt = 2??
>
> Next step is
>
> 1+(1-2)/dt*dt = 0
>
> and then it oscillates.
>
> You could say (probably correctly) that the initial condition should
> satisfy the algebraic constraints, but in my case this is proving hard to do
> (I don't explicitly know which variables appear without a time derivative),
> and even a small error leads to visible oscillations.


The correct way to handle inconsistent initial conditions is for the library
to provide explicit support for solving for them.  This is not implemented
yet, though a single short step of implicit Euler would be enough in simple
cases.  But the problem is deeper than just initial conditions.
 Oscillations will arise any time the algebraic solution changes rapidly.
 The right thing to do is to use an L-stable integrator.

There are several reasons I do not want to place the monitor at the
quadrature point.

1. It would hide oscillations that are real so the user could think they
have a decent method when in fact the accuracy is garbage.

2. The solution at quadrature points is generally lower accuracy and not
symplectic even if the propagated solution is.

3. It is inconsistent with other methods where solutions at quadrature
points may not even be available (some Runge-Kutta and IMEX methods).

4. The solution at quadrature points normally becomes available long before
it is decided whether to accept or reject a step.  Is it okay to call the
monitor anyway, or do we have to save the solution in case the monitor needs
it (turning a many-stage low-memory method into a full-memory method)?

5. Quadrature points may lie outside the interval or not increase
monotonically in time, so the monitor would be called out of order.  The
distance between points is usually irregular even if the step size is
regular.


I do not see a similarly strong case for placing the monitor at quadrature
points.  Now you might ask for a different sort of monitor that is called at
quadrature points (immediately after each SNESSolve, more-or-less), but in
light of the points above, I don't see a way to define it in a useful way.
 Maybe a special one only for TSTheta?  We could do that, but I'm still
skeptical that it's what you want.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20101221/778dfb2c/attachment.html>