Replace the basic `kexpm` routines with a variable time step version

[Simkern]

In order to speed up the exponential integrator, we need to improve the fixed-tau Krylov exponential integrator.

Design requirements:

• Re-use (as far as possible) of the computed Krylov-subspace when changing tau
• tau should be intent(inout) such that subsequent steps can avoid the (possibly expensive) first iteration requiring many Krylov vectors
• In order to allow for time step increases when possible, the input tau for the subsequent time step should be larger than the output of the last time step.
• The final routine should, if possible, conform to the existing interface.

[Simkern]

Notes on commit e315fc48d9123b8ffdde37b3e71100c50eba9f83 :

• The current implementation consistently passes the dp test but also consistently fails the sp test. I don't know why. For diagostics, switch the verbosity flag on in the corresponding test in test/TestExpmLib.fypp.
• The interface now has the additional dt (inout) parameter.
• The routine includes an internal parameter scale determining the time step adaptation factor each step. The current value is scale = 1.5 which seems to work well.
• kdim is fixed.

Reconsidering the problem, I think there can be substantial gains by also considering variable kdim. In typical use cases for LightKrylov, the by far most expensive part of the time stepping is the call to A%matvec. We don't know how expensive such a call is since it depends on the specific implementation of the user so we can't directly compare with the other parts of the algorithm. Nevertheless, the total number of matvecs should be a pretty effective proxy for optimization.

Choose dt,kdim such that dt/tau * kdim is minimized.

The difficulty compared with only varying dt is that we need some estimation of the error as a function of dt and kdim to adapt the subpsace size it since we don't want to do it willy-nilly. There is a reference by Jitse Niese and Will Wright (A Krylov subspace algorithm for evaluating the ϕ-functions appearing in exponential integrators) who give both theoretical and empirical suggestions for this estimation.

As we do not have access to A, we need to change a few details. We base our choice on an indicator given by

with the estimate of the error eps (using the last column of the dense exponential of the Hessenberg matrix) and the integration tolerance tol.

Error estimation dt

Assuming that the time-stepping method as order q, the error is approximately proportional to dt^(q+1). If no priors are available, a heuristic is q = kdim/4. If we have to reject a time step dt, an estimate for the heuristic order q is

based on the old time step dt_old and the error estimates for the old and current steps.

With this estimate, we predict the best new timestep:

Error estimation kdim

I heuristic estimate of the error as a function of the Krylov subspace size is proportional to k^kdim. If we have no reference, the authors suggest k = 2 or, given a previous estimate with an previous kdim_old

to choose the next kdim as

Balancing dt and kdim

In practice, we change only the time step or the Krylov subspace size, the choice being based on the expected cost to complete the calculation, in our case as a function of the number of calls to A%matvec.

We choose the smaller side:

Initialization

The paper proposes an initialization for the time step based on the norm of A. Instead, starting from either the input value of dt or the simple dt = tau, we can use successive halving of the time step at a fixed kdim (10 seems reasonable) until the tolerance is achieved with that subspace size. Once this reasonable initial time step is found, the algorithm above can be used ensuring that the changes are not too drastic and the error estimates are applicable.