TS grand plan

[Barry Smith]

    Jed,

     I suggest you "start from scratch" with "f(t,x,x')" and design  
what makes sense for you and ignore the current Arhs etc.
Then we can try to transition the simple ones we have now as you  
suggest.  Please look at the explicit RK in TS and see how
that may fit into your new grand plan.

    TS NEVER had a unified plan, which is a big reason why it is no  
used much and is pretty much a POS.


    Barry

On Sep 20, 2009, at 11:33 PM, Jed Brown wrote:

> I'm putting in support for DAEs written in fully implicit form,
> i.e. F(t,x,x')=0.  (I brought this up ages ago, but just got to work
> implementing it.)  TS uses a frustrating mixed design in an attempt to
> reuse code between the different implementations.  Is there clear  
> intent
> behind A, B, Alhs, Arhs?  My expectation was that we would be calling
> KSPSetOperators(ksp,A,B,...) and SNESSetJacobian(snes,A,B,...) where  
> the
> other matrices are for the user, but the usage is much less clear.
>
> TSSetUp:
>  KSPSetOperators(ksp,Arhs,B,...)
>  SNESSetJacobian(snes,Arhs,B,...)
>
> TS_BEULER:
>  A = Arhs; (aliasing)
>  KSPSetOperators(ksp,A,A,...)
>  SNESSetJacobian(snes,Arhs,B,...)  [shift applied to A=Arhs]
>
>  * -snes_mf_operator doesn't do anything matrix-free, see
>    tutorials/ex2 -nox -ts_max_steps 1 -ts_type beuler - 
> snes_mf_operator -snes_view
>    There was a hack to detect MFFD which was introduced in 2004 and
>    removed a day later (although the associated comment remains)
>  * No support for shifting a preconditioning matrix
>    (also made an appearance in 2004)
>  * No way to get a different preconditioning matrix to the KSP
>
> TS_CN:
>  KSPSetOperators(ksp,A,A,...)
>  SNESSetJacobian(snes,A,B,...)
>
>  * linear variable matrix calls ts->ops->rhsmatrix with NULL  
> preconditioning
>  matrix
>  * For linear problems with variable matrix, A is destroyed and
>  recreated every step
>  * For nonlinear problems, only B is shifted
>
>
> Is there some consistent grand plan?
>
>
> More trivially, I ripped isExplicit and Iindex out of _p_TS because  
> they
> are never used in PETSc.
>
>
> As a very basic DAE integrator, I wrote an implementation of the theta
> method as a 1-step Runge-Kutta which only uses the implicit form
> F(t,x,xdot)=0.  With the implicit formulation, the matrices Alhs, Arhs
> never appear and instead the user assembles exactly the matrix that is
> needed by SNES.  Of course, we can serve up the various explicit
> representations (xdot = f(t,x) and the linear cases) with light
> wrappers, but I feel like I'm fighting the use of A,B,Alhs,Arhs in
> interface/ts.c.
>
> It seems to me that with a little caching, we should be able to make  
> all
> the implicit integrators use only the F(t,x,xdot)=0 form and still get
> comparable efficiency to the specialized versions that we have now.
> Since TS_THETA includes the functionality of TS_BEULER and TS_CN as
> special cases, we could remove the latter if we can demonstrate that  
> the
> fully implicit formulation doesn't lose efficiency for the RHS and
> linear cases (I expect it will be more efficient when the user  
> provides
> the implicit form because we get to skip scaling and shifting).
>
> Jed
>