<div class="gmail_quote">On Tue, Dec 21, 2010 at 20:07, Michael E Henderson <span dir="ltr"><<a href="mailto:mhender@us.ibm.com">mhender@us.ibm.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
<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.</font></blockquote><div><br></div><div>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.</div>
<div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;"><font size="2" face="sans-serif"> 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></blockquote></div><br><div>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.</div>
<div><br></div><div>There are several reasons I do not want to place the monitor at the quadrature point.</div><div><br></div><div>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.</div>
<div><br></div><div>2. The solution at quadrature points is generally lower accuracy and not symplectic even if the propagated solution is.</div><div><br></div><div>3. It is inconsistent with other methods where solutions at quadrature points may not even be available (some Runge-Kutta and IMEX methods).</div>
<div><br></div><div>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)?</div>
<div><br></div><div>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.</div>
<div><br></div><div><br></div><div>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.</div>