<font size=2 face="sans-serif">Fine, Your call. </font>
<br>
<br><font size=2 face="sans-serif">------------------------------------------------------------------------------------------------------------------------------------<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 02:28 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 20:07, 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">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>
<br>
<br><font size=3>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.</font>
<br><font size=3> </font>
<br><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><font size=3>
<br>
</font><font size=2 face="sans-serif"><br>
For theta=1/2, x_0 is input, and </font><font size=3><br>
</font><font size=2 face="sans-serif"><br>
f(x_1-x_0)/dt,(x_1+x_0)/2)=0 is solved for x_1.</font><font size=3> <br>
</font><font size=2 face="sans-serif"><br>
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><font size=3> <br>
</font><font size=2 face="sans-serif"><br>
1+(1-0)/dt*dt = 2??</font><font size=3> <br>
</font><font size=2 face="sans-serif"><br>
Next step is </font><font size=3><br>
</font><font size=2 face="sans-serif"><br>
1+(1-2)/dt*dt = 0</font><font size=3> <br>
</font><font size=2 face="sans-serif"><br>
and then it oscillates.</font><font size=3> <br>
</font><font size=2 face="sans-serif"><br>
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=3>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.</font>
<br>
<br><font size=3>There are several reasons I do not want to place the monitor
at the quadrature point.</font>
<br>
<br><font size=3>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.</font>
<br>
<br><font size=3>2. The solution at quadrature points is generally lower
accuracy and not symplectic even if the propagated solution is.</font>
<br>
<br><font size=3>3. It is inconsistent with other methods where solutions
at quadrature points may not even be available (some Runge-Kutta and IMEX
methods).</font>
<br>
<br><font size=3>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)?</font>
<br>
<br><font size=3>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.</font>
<br>
<br>
<br><font size=3>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.</font>
<br>