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

[Michael E Henderson]

Fine, Your call. 

------------------------------------------------------------------------------------------------------------------------------------
Mathematical Sciences, TJ Watson Research Center
mhender at watson.ibm.com
http://www.research.ibm.com/people/h/henderson/
http://multifario.sourceforge.net/




From:   Jed Brown <jed at 59A2.org>
To:     For users of the development version of PETSc 
<petsc-dev at mcs.anl.gov>
Date:   12/21/2010 02:28 PM
Subject:        Re: [petsc-dev] Fw: TSTheta for DAE's
Sent by:        petsc-dev-bounces at mcs.anl.gov



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