[petsc-dev] Strange behaviour from TSARKIMEX

[Emil Constantinescu]

On 11/15/16 7:35 AM, Stefano Zampini wrote:
> Emil,
>
> I have modifed the example code to also include the convection matrix in the lhs.
>
> Plain arkimex with -ts_arkimex_fully_implicit false and M*udot - K*u in the lhs works well.

ts_arkimex_fully_implicit true makes a new F(u,udot,t) <= 
F(u,udot,t)-G(u,t) and operates only with it, so that moves you in a 
different formulation.

> Trying to summarize the thread, it seems that we may have problems using arkimex with default parameters when F(u,udot,t) will just depend on udot and not on u.

This integrator is developed for udot = g(u) + f(u) as a tightly coupled 
time stepper. Other forms listed in the table are supported by doing 
some acrobatics: e.g., M udot - f(u) = 0, udot = M^(-1) g(u),... In 
other words we either use just the implicit integrator for F(u,udot,t)=0 
or assume a certain form for F(u,udot,t), such as M udot, but there is 
no direct support for a mass matrix; however, there are indirect ways of 
doing that. The point is that we cannot partition the equation randomly 
and expect it to work.

> If this is the case, is the user supposed to write G(t,u) as F_udot ^(-1) g(t,u)? So why not provide an API to allow it, instead to have the user modify the G function?

One element that could be added is direct support for a mass matrix; 
however, that would be somewhat difficult to manage (pretty much all 
solvers will need support and only a couple of methods would use it. 
That can be done on the user side. The reason is that F is broke up 
depending on it's form. There are additional (non)linear solves 
depending on what -tsarkimex_type method is used and the form of F. Note 
that the TS solver has been extended from udot = g(u) + f(u) to handling 
many other situations such as DAEs.

> Also, wouldn’t be better to have standard SDIRK the default? i.e. -ts_arkimex_fully_implicit false?

ARKIMEX is meant for IMEX methods and -ts_arkimex_fully_implicit is 
false by default. SDIRKs are nice by not great for stiff problems. So 
the explicit part is chosen to work with a good ESDIRK method. If just 
the implicit part is desired, then the user can use 
-ts_arkimex_fully_implicit true. Actually Table 12 should be updated as 
well to reflect that.

In conclusion, I agree that the documentation can benefit from 
clarifications. Currently the easiest way for the user to go about it is 
to look at the table entries and identify which partition works best in 
terms of efficiency and effectiveness - which depends in part on what's 
listed in Table 12.

Emil

> Stefano
>
> On Nov 14, 2016, at 11:11 PM, Emil Constantinescu <emconsta at mcs.anl.gov> wrote:
>
>>
>>
>> On 11/14/16 1:17 PM, Stefano Zampini wrote:
>>> I do have pasted the output for ts_monitor for the first three steps.
>>
>> I meant -ts_adapt_monitor, I wanted to see what TS solvers are used: The first step is probably an Euler method, that's why you have 10 SNES solves in place of 7 for the later ones.
>>
>>> The entry you are referring to in Table 11 actually breaks the TS
>>> assumption that the ODE can be written F(u,udot,t) = G(u,t) (I was
>>> thinking it was a typo)
>>
>> It works for certain forms as listed in the table. I'm not aware of any RK method that handles this case with G(u,t) treated explicitly.
>>
>>>> In your case, you would either use the setup corresponding to "stiff
>>>> ODE w/ mass matrix" or "nonstiff  ODE  w/  mass  matrix”.
>>>
>>> I’m not interested in solving this specific advection problem with
>>> arkimex. I came across this problem when writing a general interface to
>>> TS from a FEM library.
>>
>> Normally, I do that by solving udot=M^-1 f(u) + M^-1 g(u)
>>
>>>> It is difficult to create a decision tree for every way of writing the
>>>> splittings. Do you find Table 11 not clear about what splitting to
>>>> use? I would welcome any kind of feedback for improving it.
>>>>
>>>
>>> If the entry in Table 11 is correct, please add a comment on the manual
>>> stating that ARKIMEX does not fully support the F(u,udot,t) = G(u,t)
>>> interface.
>>
>> I see how the first statement in the caption can be confusing; The last statement in the caption states:
>> "General cases such as F(t;y;y')=G(t;y) are not amenable to IMEX Runge-Kutta, but can be solved by using fully implicit methods." A line in the table may be more clear. I'll go for that.
>>
>> Emil
>