On Tue, Sep 22, 2009 at 12:02 PM, Lisandro Dalcin <span dir="ltr"><<a href="mailto:dalcinl@gmail.com">dalcinl@gmail.com</a>></span> wrote:<br><div class="gmail_quote"><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
I commented about this in the past, but let's do it once more...<br>
<br>
Just for reference, in petsc4py (accesible from core PETSc as TSPYTHON<br>
type when configured --with-petsc4py) more or less uses Jed's model:<br>
the user routines are in charge of computing F(...)=0 and the<br>
Jacobian. Moreover, this code can also manage some basic time-step<br>
control.<br>
<br>
BTW, I proposed something similar to this in the past (look for<br>
subject "flexible TS implementation for user-defined timestepping").<br>
Matt (strongly?) resisted the idea, and then I did not continued<br>
advocating for it<br><div><div class="h5"></div></div></blockquote><div><br>I think my objections were based upon my conviction that methods that need<br>specific operators or restrictions of operators should have a high level way<br>
of expressing this which incorporates the discretization. Since this seems far<br>away (unfortunately), I have stopped opposing the "good enough" stuff.<br><br>  Matt<br> </div><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
<div><div class="h5">
On Mon, Sep 21, 2009 at 12:01 PM, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>> wrote:<br>
><br>
>   Jed,<br>
><br>
>    I suggest you "start from scratch" with "f(t,x,x')" and design what makes<br>
> sense for you and ignore the current Arhs etc.<br>
> Then we can try to transition the simple ones we have now as you suggest.<br>
>  Please look at the explicit RK in TS and see how<br>
> that may fit into your new grand plan.<br>
><br>
>   TS NEVER had a unified plan, which is a big reason why it is no used much<br>
> and is pretty much a POS.<br>
><br>
><br>
>   Barry<br>
><br>
> On Sep 20, 2009, at 11:33 PM, Jed Brown wrote:<br>
><br>
>> I'm putting in support for DAEs written in fully implicit form,<br>
>> i.e. F(t,x,x')=0.  (I brought this up ages ago, but just got to work<br>
>> implementing it.)  TS uses a frustrating mixed design in an attempt to<br>
>> reuse code between the different implementations.  Is there clear intent<br>
>> behind A, B, Alhs, Arhs?  My expectation was that we would be calling<br>
>> KSPSetOperators(ksp,A,B,...) and SNESSetJacobian(snes,A,B,...) where the<br>
>> other matrices are for the user, but the usage is much less clear.<br>
>><br>
>> TSSetUp:<br>
>>  KSPSetOperators(ksp,Arhs,B,...)<br>
>>  SNESSetJacobian(snes,Arhs,B,...)<br>
>><br>
>> TS_BEULER:<br>
>>  A = Arhs; (aliasing)<br>
>>  KSPSetOperators(ksp,A,A,...)<br>
>>  SNESSetJacobian(snes,Arhs,B,...)  [shift applied to A=Arhs]<br>
>><br>
>>  * -snes_mf_operator doesn't do anything matrix-free, see<br>
>>   tutorials/ex2 -nox -ts_max_steps 1 -ts_type beuler -snes_mf_operator<br>
>> -snes_view<br>
>>   There was a hack to detect MFFD which was introduced in 2004 and<br>
>>   removed a day later (although the associated comment remains)<br>
>>  * No support for shifting a preconditioning matrix<br>
>>   (also made an appearance in 2004)<br>
>>  * No way to get a different preconditioning matrix to the KSP<br>
>><br>
>> TS_CN:<br>
>>  KSPSetOperators(ksp,A,A,...)<br>
>>  SNESSetJacobian(snes,A,B,...)<br>
>><br>
>>  * linear variable matrix calls ts->ops->rhsmatrix with NULL<br>
>> preconditioning<br>
>>  matrix<br>
>>  * For linear problems with variable matrix, A is destroyed and<br>
>>  recreated every step<br>
>>  * For nonlinear problems, only B is shifted<br>
>><br>
>><br>
>> Is there some consistent grand plan?<br>
>><br>
>><br>
>> More trivially, I ripped isExplicit and Iindex out of _p_TS because they<br>
>> are never used in PETSc.<br>
>><br>
>><br>
>> As a very basic DAE integrator, I wrote an implementation of the theta<br>
>> method as a 1-step Runge-Kutta which only uses the implicit form<br>
>> F(t,x,xdot)=0.  With the implicit formulation, the matrices Alhs, Arhs<br>
>> never appear and instead the user assembles exactly the matrix that is<br>
>> needed by SNES.  Of course, we can serve up the various explicit<br>
>> representations (xdot = f(t,x) and the linear cases) with light<br>
>> wrappers, but I feel like I'm fighting the use of A,B,Alhs,Arhs in<br>
>> interface/ts.c.<br>
>><br>
>> It seems to me that with a little caching, we should be able to make all<br>
>> the implicit integrators use only the F(t,x,xdot)=0 form and still get<br>
>> comparable efficiency to the specialized versions that we have now.<br>
>> Since TS_THETA includes the functionality of TS_BEULER and TS_CN as<br>
>> special cases, we could remove the latter if we can demonstrate that the<br>
>> fully implicit formulation doesn't lose efficiency for the RHS and<br>
>> linear cases (I expect it will be more efficient when the user provides<br>
>> the implicit form because we get to skip scaling and shifting).<br>
>><br>
>> Jed<br>
>><br>
><br>
><br>
<br>
<br>
<br>
</div></div><font color="#888888">--<br>
Lisandro Dalcín<br>
---------------<br>
Centro Internacional de Métodos Computacionales en Ingeniería (CIMEC)<br>
Instituto de Desarrollo Tecnológico para la Industria Química (INTEC)<br>
Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET)<br>
PTLC - Güemes 3450, (3000) Santa Fe, Argentina<br>
Tel/Fax: +54-(0)342-451.1594<br>
</font></blockquote></div><br><br clear="all"><br>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener<br>