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

[Michael E Henderson]

Thanks for responding Jed.

I'm suggesting that the monitor be passed th-X after the SNES solver 
solves the implicit equations, with whatever time value that corresponds 
to.

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

Does that make sense?

Thanks,

Mike
------------------------------------------------------------------------------------------------------------------------------------
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 01:53 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 18:48, Michael E Henderson <mhender at us.ibm.com> 
wrote:
Who can I talk to about the implementation of the TSTheta Implicit time 
stepper? 

The implementation computes ts->vec_sol which solves the implicit system. 
But the monitor returns ts->vec_sol+dt*th->Xdot, which does not 
necessarily satisfy the equations. At least that's how I read the code. I 
think that the call to VecAXPY which changes ts->vec_sol needs to be moved 
to after the call to TSMonitor in src/ts/impls/implicit/theta/theta.c 
Purely a matter of what is passed to the monitor routine. 

The state propagated by the method is the one that the monitor sees. 
 Calling the monitor at the quadrature point would be inconsistent, in my 
opinion.

The Theta method is not L-stable, it is not a robust method for stiff 
systems (but it is efficient and happens to work fairly frequently).  The 
standard alternative is BDF-2 which we don't currently have an 
implementation of, but would be fairly easy to add.  The downside of BDF-2 
is that it's less efficient and less accurate.  Note that theta=1 is 
implicit Euler which has every stability property you could want, but is 
only first order accurate.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20101221/29a209a7/attachment.html>