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 <[email protected]> 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>
