[petsc-dev] Strange behaviour from TSARKIMEX

[Zhang, Hong]

> On Nov 15, 2016, at 10:56 AM, Emil Constantinescu <emconsta at mcs.anl.gov> wrote:
> 
> 
> 
> 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.

Like +1. 

Given the equation M*udot=f(u), would F(u,dot,t)=M*udot, G(t,u)=f(u) be a valid partition for IMEX? Obviously not. This is like pushing everything to the explicit part, leaving nothing in the implicit part. Essentially this partition makes PETSc solve the implicit formulation M*udot=f(u) with an explicit time integration method, which does not make sense. Also it is not a good idea for the solver to try to invert the mass matrix M because M could be singular.

This is a good example for an invalid partitioning. Maybe we should put it in the manual.

>> 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
>>