<font size=2 face="sans-serif">Thanks for responding Jed.</font>
<br>
<br><font size=2 face="sans-serif">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.</font>
<br>
<br><font size=2 face="sans-serif">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.</font>
<br>
<br><font size=2 face="sans-serif">For theta=1/2, x_0 is input, and </font>
<br>
<br><font size=2 face="sans-serif">f(x_1-x_0)/dt,(x_1+x_0)/2)=0 is solved
for x_1.</font>
<br>
<br><font size=2 face="sans-serif">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</font>
<br>
<br><font size=2 face="sans-serif">1+(1-0)/dt*dt = 2??</font>
<br>
<br><font size=2 face="sans-serif">Next step is </font>
<br>
<br><font size=2 face="sans-serif">1+(1-2)/dt*dt = 0</font>
<br>
<br><font size=2 face="sans-serif">and then it oscillates.</font>
<br>
<br><font size=2 face="sans-serif">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.</font>
<br>
<br><font size=2 face="sans-serif">Does that make sense?</font>
<br>
<br><font size=2 face="sans-serif">Thanks,</font>
<br>
<br><font size=2 face="sans-serif">Mike<br>
------------------------------------------------------------------------------------------------------------------------------------<br>
Mathematical Sciences, TJ Watson Research Center<br>
mhender@watson.ibm.com<br>
</font><a href=http://www.research.ibm.com/people/h/henderson/><font size=2 face="sans-serif">http://www.research.ibm.com/people/h/henderson/</font></a><font size=2 face="sans-serif"><br>
</font><a href=http://multifario.sourceforge.net/><font size=2 face="sans-serif">http://multifario.sourceforge.net/</font></a><font size=2 face="sans-serif"><br>
</font>
<br>
<br>
<br>
<br><font size=1 color=#5f5f5f face="sans-serif">From:
</font><font size=1 face="sans-serif">Jed Brown <jed@59A2.org></font>
<br><font size=1 color=#5f5f5f face="sans-serif">To:
</font><font size=1 face="sans-serif">For users of the development
version of PETSc <petsc-dev@mcs.anl.gov></font>
<br><font size=1 color=#5f5f5f face="sans-serif">Date:
</font><font size=1 face="sans-serif">12/21/2010 01:53 PM</font>
<br><font size=1 color=#5f5f5f face="sans-serif">Subject:
</font><font size=1 face="sans-serif">Re: [petsc-dev]
Fw: TSTheta for DAE's</font>
<br><font size=1 color=#5f5f5f face="sans-serif">Sent by:
</font><font size=1 face="sans-serif">petsc-dev-bounces@mcs.anl.gov</font>
<br>
<hr noshade>
<br>
<br>
<br><font size=3>On Tue, Dec 21, 2010 at 18:48, Michael E Henderson <</font><a href=mailto:mhender@us.ibm.com><font size=3 color=blue><u>mhender@us.ibm.com</u></font></a><font size=3>>
wrote:</font>
<br><font size=2 face="sans-serif">Who can I talk to about the implementation
of the TSTheta Implicit time stepper?</font><font size=3> <br>
</font><font size=2 face="sans-serif"><br>
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.</font><font size=3> </font>
<br>
<br><font size=3>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.</font>
<br>
<br><font size=3>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.</font>
<br>