[petsc-dev] [petsc-maint #42811] TSPSEUDO and snes ex27

[Barry Smith]

On Mar 8, 2010, at 3:34 AM, Jed Brown wrote:

> I was just looking at TSPSEUDO to put in DAE support (which basically
> amounts to using TSComputeIFunction/TSComputeIJacobian instead of
> rolling this logic in the implementation with TSScaleShiftMatrices)  
> and
> notice that it is not doing pseudotransient continuation as  
> described in
> the literature.  For example, see Eq 2.3 of
>
>  Todd S. Coffey and C. T. Kelley and David E. Keyes, Pseudotransient
>  Continuation and Differential-Algebraic Equations, 2003.
>
> To explain, I'll write the following in DAE form where the transient
> system is
>
>  F(X,Xdot) = 0
>
> and pseudotransient continuation should be solving
>
>  X += inv(G') F(X,0)
>
> where
>
>  G(Y) = F(Y,(Y-X)/dt)
>
> and dt is chosen with switched evolution relaxation (SER) which  
> relates
> residual norms at each step.  This is linearly implicit Euler (i.e. a
> Rosenbrock method) except that the residual is the steady state  
> residual
> instead of the transient residual. TSPSEUDO does nonlinearly implicit
> Euler with the transient residual and SER.
>
> I then looked at snes ex27.c to see how it was implemented for the
> paper, and it turns out that they are doing the same thing  
> (nonlinearly
> implicit Euler with transient residual).  The runex27 and runex27_p
> targets in the makefile date from this checkin
>
>  changeset:   1534:85d178697bd4
>  parent:      1532:4f36db04e928
>  user:        kaushik at manu.mcs.anl.gov
>  date:        Sat Nov 02 16:29:16 2002 -0600
>  files:       src/snes/examples/tutorials/ex27.c src/snes/examples/ 
> tutorials/makefile
>  description:
>  bk-changeset-1.866.4.1
>  kaushik at manu.mcs.anl.gov|ChangeSet|20021102222916|48153
>  ChangeSet
>    1.866.4.1 02/11/02 16:29:16 kaushik at manu.mcs.anl.gov +2 -0
>    runex27 and runex27_p in the makefile correpond to David Keyes's  
> runtime options (used
>    in Fig 3.1 of his paper at http://www.icase.edu/~keyes/papers/ptc2.pdf) 
> .
>
>    src/snes/examples/tutorials/makefile
>      1.16 02/11/02 16:29:13 kaushik at manu.mcs.anl.gov +4 -0
>      Added runex27 (DAE) and runex27_p (parabolic) with David Keyes's
>      runtime options
>
>    src/snes/examples/tutorials/ex27.c
>      1.3 02/11/02 16:29:13 kaushik at manu.mcs.anl.gov +11 -6
>      Put the command line option to choose between parabolic and DAE.
>
> and running these targets today reproduces Figure 3.1 of the paper,  
> but
> is definitely not doing the algorithm stated in the paper.  I'm a bit
> unsure how to resolve this.
>
> As for making TSPSEUDO do what the papers assert, one alternative  
> would
> be to use KSP directly, but this requires reproducing some code to
> support matrix free (and with coloring).  More seriously, it makes the
> scheme far less robust because a useful mode (mentioned in some of the
> papers) is to move all the way to the steady state problem if dt  
> exceeds
> some threshold, in which case we want a real Newton with line search
> once that occurs.
>
> If an analytic Jacobian is available, then we can just run Newton for
> one step with no line search (because it's sometimes correct for
> residuals to increase after a step of pseudotransient continuation)
> until the threshold is exceeded at which point we reactivate the line
> search.  But with an MFFD Jacobian, we have to be able to distinguish
> between when residual evaluation is being called for Jacobian purposes
> (transient term should be included) and when it is being called for
> steady-state residual puproses (transient term should not be  
> included).

    The matrix-free finite differencing in SNES uses the MATMFFD. You  
could overwrite the use of the SNES function to use the time-dependent  
one by calling
MatMFFDSetFunction() "after the fact" (that is to overwrite the  
default one from SNES. Similarly in the case of  
SNESDefaultComputeJacobianColoring() the function  
MatFDColoringSetFunction() could be called to use time-dependent  
function instead of the original. Maybe.

The drawback to this is, of course, "Jacobian" specific code in the  
pseudo transient solver.

    BTW: I think the right thing to do is to bag the pseudo transient  
solver support in TS and "somehow" put the support into SNES.  If we  
had a model that let us compose SNES solvers, like we have for PCs  
this would be straightforward. So the first step is to fix up the SNES  
model to allow composing them "somehow".

   Moving this discussion to petsc-dev


    Barry




>
> In any case, unlike supporting DAE, these changes are not going to  
> make
> it into the release.
>
> Jed
>