[petsc-dev] Fw: TSTheta for DAE's
Michael E Henderson
mhender at us.ibm.com
Tue Dec 21 13:56:29 CST 2010
Fine, Your call.
------------------------------------------------------------------------------------------------------------------------------------
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 <jed at 59A2.org>
To: For users of the development version of PETSc
<petsc-dev at mcs.anl.gov>
Date: 12/21/2010 02:28 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 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/a9c6c2c3/attachment.html>
More information about the petsc-dev
mailing list