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>
