Choosing between adaptive BDF and IRK integrators: adaptive meshing effect?

[kosa...@ualberta.ca]

Dear deal.ii team,

I am currently at the stage of choosing a transient solver for our deal.ii-based software OpenFCST. The solver will be applied to solving nonlinear transport and reaction problems. After days of research, I found two options:

1. Variable order, variable step BDF schemes (as seen in LSODA, VODE; they are combined with Adams methods there);
2. Variable order, variable step implicit Runge-Kutta schemes.

The performance of both in general stiff problems is more or less similar with BDF being superior when used with the method of lines. I really like the idea of variable order and step implementation of the multistep methods (see, for example, this work of L. Petzold). They have become a "must-have" for automated ODE solvers and are present in Matlab and Python, for example (RK methods are present there as well, to be fair). Their high flexibility allows to adapt to the stiffness of the problem that may change during time integration. This is important in the problems that describe coupled transport and reactions, such as the ones in our electrochemical models.

There is, however, one important drawback of the multistep methods (including BDF of order 2<k<6 that I am looking at): they require solution transfer from more than one previous steps compared to just one in one-step methods (Runge-Kutta, theta scheme). Therefore, if the mesh is continuously refined/coarsened during the simulation, this will lead to increase in computational expenses due to necessity of interpolation of several solutions onto the new grid.

I was 100 % sure in which method I was going to choose when I suddenly realized that shortcoming of the BDF method. How computationally complex is the implementation of the solution transfer in deal.ii? Will it really cause significant performance decrease in problems with, for example, rapid transients, during which the mesh would be continuously refined/coarsened and a higher order BDF scheme would be automatically chosen (the higher the order, the more solutions from previous time layers the method requires)? If I am not mistaken, IRK would also suffer from this case since, when simplified Newton iterations are used, it will require recomputation of m Jacobians (m - number of stages in the method) based on the transferred old solution. So, my question likely sums up to what will be more expensive (assuming particular implementations in deal.ii), to transfer k old solutions or to recompute k Jacobians? I do, of course, realize that complexity of computing a Jacobian depends on the particular problem. We can assume up to five arithmetic operations inside three for loops (for quadrature points and (i,j) components) for each cell.

I hope this discussion is not too broad for this message board. I will appreciate any help as I am currently unsure in how to proceed.

Regards,
Aslan Kosakian

Graduate Student
Energy Systems Design Laboratory
Department of Mechanical Engineering
University of Alberta
Edmonton, AB, Canada
T6G2G8

[Wolfgang Bangerth]

Aslan,

> There is, however, one important drawback of the multistep methods (including
> BDF of order 2<k<6 that I am looking at): they require solution transfer from
> more than one previous steps compared to just one in one-step methods
> (Runge-Kutta, theta scheme). Therefore, if the mesh is continuously
> refined/coarsened during the simulation, this will lead to increase in
> computational expenses due to necessity of interpolation of several solutions
> onto the new grid.

The way you would implement this is that every time you change the mesh, you
transfer all solution vectors from the previous to the new mesh. The old mesh,
and all solution vectors defined on it, then cease to exist. Everything only
exists on the new mesh. For all practical purposes, at any given time of the
time stepping method, you only ever work on one mesh and it looks like there
has never been another mesh. The only time when you deal with two meshes is
during the brief period where you use SolutionTransfer.

Now, transferring solution vectors is not cheap, but it's also not expensive
compared to, for example, assembling and solving a linear system. That's
particularly true because you don't typically adapt the mesh very frequently.
For most problems, I would assume that adapting the mesh every 10 or 20 time
steps is perfectly fine. For example, in our "big" code, ASPECT, the default
is every 15 time steps, and I don't know that anyone would have ever changed
that default in their applications. In any case, if you amortize the cost of
solution transfer over 10 or 20 time steps, then it's pretty much negligible.
I would say that's true regardless of whether you transfer one, or ten
solution vectors.

Best
Wolfgang

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[kosa...@ualberta.ca]

Dear Dr. Bangerth,

Thank you very much for your quick reply. In theta scheme that I currently use, I do limit refinement frequency to each n-th layer, normally n=10.

I can see that there is an implementation of IRK method in deal.ii. If you do not mind, may I ask what was the reason that that particular method was chosen in your software? I am aware that IRK and BDF schemes have overall similar efficiency, so that might have been just some personal preference of the developer(s). I am wondering how deal.ii solves the issue when the dimension of the problem is large and then it is impossible to store the matrix (I-tau*AxJ)^(-1) of the simplified Newton iteration in IRK (where x is the Kronecker product; from Hairer & Wanner). I believe you use some modification that leads to the use of simply (I-tau*J)^(-1) as it is seen from the name id_minus_tau_J_inverse of one of the variables in the code. Is it recomputed for each Newton iteration?

Regards,
Aslan Kosakian

Graduate Student
Energy Systems Design Laboratory
Department of Mechanical Engineering
University of Alberta
Edmonton, AB, Canada
T6G2G8

[Wolfgang Bangerth]

Aslan,

> I can see that there is an implementation of IRK method in deal.ii. If you do
> not mind, may I ask what was the reason that that particular method was chosen
> in your software? I am aware that IRK and BDF schemes have overall similar
> efficiency, so that might have been just some personal preference of the
> developer(s).

That would be a question for Bruno and Damien. I don't know.

> I am wondering how deal.ii solves the issue when the dimension
> of the problem is large and then it is impossible to store the matrix
> (I-tau*AxJ)^(-1) of the simplified Newton iteration in IRK (where x is the
> Kronecker product; from Hairer & Wanner). I believe you use some modification
> that leads to the use of simply (I-tau*J)^(-1) as it is seen from the name

> /id_minus_tau_J_inverse/ of one of the variables in the code. Is it recomputed
> for each Newton iteration?

You never store this matrix. What you need to provide to the integrators is a
function that multiplies a vector with this matrix. In practice, that always
means that you solve a linear system with (I - tau J) and the given vector as
the right hand side.

Take a look at the tutorial program that demonstrates the use of these time
integrators to see an example of how this is implemented in practice.

Best
W.

[kosa...@ualberta.ca]

Thank you very much. I CC Bruno Turcksin and Damien Lebrun-Grandie in hope of receiving a reply from them.

Regards,
Aslan

[Bruno Turcksin]

Aslan,

2016-11-17 19:31 GMT-05:00 <kosa...@ualberta.ca>:
> On Thursday, November 17, 2016 at 5:22:00 PM UTC-7, Wolfgang Bangerth wrote:
>> > I can see that there is an implementation of IRK method in deal.ii. If
>> > you do
>> > not mind, may I ask what was the reason that that particular method was
>> > chosen
>> > in your software? I am aware that IRK and BDF schemes have overall
>> > similar
>> > efficiency, so that might have been just some personal preference of the
>> > developer(s).
>>
>> That would be a question for Bruno and Damien. I don't know.

The reason is that all Runge-Kutta methods (including forward and
backward Euler) can be expressed using a Butcher tableau. Once you
have the code working for one explicit RK method and you have a Newton
solver, you can easily add any Runge-Kutta method by adding a new
Butcher tableau. Since this was a very small amount of work, we added
IRK schemes. There is no numerical reasons why IRK was chosen over
BDF.

Best,

Bruno

[kosa...@ualberta.ca]

Good evening,

Thank you very much for the reply. Indeed, the use of the Butcher tableau for constructing RK methods is very convenient. This discussion was very helpful, so I would like to thank again Dr. Bangerth and Dr. Turcksin for their quick responses.

Regards,
Aslan