[WIP] Relaxation Runge-Kutta

[Theozeud]

Start of the implementation of relaxation Runge-Kutta methods (see #1029)

Instead of writing things with a callback as suggested in the issue, I prefer to propose, I think, a more suitable solution by decoupling perform_step! into several steps one of which will allow the user to do modifications on the computations of the method. I think the problem with the callback is when it occurs : in the loopfooter function where things that have already modified have been registered as dtnew, t, tprev and others things related to tstops, which could be changed by the relaxation methods.

For the moment, to test the new structure, I copy the Tsit5 method with another name Tsit5_for_relaxation. I will certainly change some parts of the code and its structure. I am not sure yet if the idea is good, but I think it is a good starting point.

To do

From the relaxation method point of view :

• [x] Implement a copy of Tsit5 with the new structure explained above,
• [x] Test if the code run without the relaxation,
• [x] Implement a structure to perform a relaxation and look for (a) good package(s) and method(s) to do the needed optimization step,
• [x] Implement FSAL-R,
• [x] See if it is possible to implement easily R-FSAL,
• [x] Look for suitable problems for the relaxation method,
• [x] Make the relaxation method running,
• [x] Make the relaxation method working,
• [x] Make benchmark.

From an interface point of view :

• [ ] See which variables are really needed to be added to the structure integrator,
• [x] See if the splitting in three steps of the perform_step! function is a good choice or if it is better to do only two,
• [x] Investigate all things needed to meticulously take care of time as matching tstops ect,
• [x] Find out what need to be changed in DiffEqBase.jl if we want to add a solver option (because it looks like something needs to be changed)

[ChrisRackauckas]

Alerting @ranocha who may be interested in keeping on top of this.

[Theozeud]

Two comments on the current implementation :

• The new structure seems to work at least without any modification of the calculation by a user but there is a small bug with interpolation that I do not understand for the moment

• I divided perform_step ! into three steps:

  • computation!
  • modify_step !

  • finalize_step !

    The first is to do the mathematical side of the method, the second is to allow the user to make modifications to the calculation and the solver to track it, and the last is to do all the work related to the DiffEq solver. But it is not easy to divide the tasks between computation ! and finalize_step ! In particular, we need to reload some objects to do this, so I am not sure it is really necessary, except perhaps to recalculate fsal or things like that if u changes.

[JoshuaLampert]

Regarding the step

Implement a structure to perform a relaxation and look for a good package/method to do the needed optimization step

in your TODO list above, I have an implementation of the relaxation step (there it is implemented as a callback) in https://github.com/JoshuaLampert/DispersiveShallowWater.jl/blob/a7eb1307169d1278f0c4e1afeae55b3e9a8cea64/src/callbacks_step/relaxation.jl, which might or might not help you.

[Theozeud]

Interesting development. If you want to achieve best performance with step size control, you can have a look at our recent preprint https://arxiv.org/abs/2311.14050

Feel free to ask me if you have any specific questions while working on this

Regarding the step

Implement a structure to perform a relaxation and look for a good package/method to do the needed optimization step

in your TODO list above, I have an implementation of the relaxation step (there it is implemented as a callback) in https://github.com/JoshuaLampert/DispersiveShallowWater.jl/blob/a7eb1307169d1278f0c4e1afeae55b3e9a8cea64/src/callbacks_step/relaxation.jl, which might or might not help you.

Ok, thanks to both of you, I will take a look!

[Theozeud]

Update on my work :

• [x] The new structure with relaxation runs and works partially. I need to fix things related to tstops matching.
• [x] I am reading the article mentionned by Hendrik Ranocha https://arxiv.org/abs/2311.14050 about better step size control.

[Theozeud]

@ChrisRackauckas, @ranocha I have a small question if someone can help me. I compare my implementation Tsit5_for_relaxtion with the existing Tsit5 on a simple problem without relaxation such that results must be the same. I encounter a small problem of interpolation which could be due to my code, but I wonder if it could be due to the fact I do not define _ode_interpolant for Tsit5_for_relaxation.

[oscardssmith]

I wonder if it could be due to the fact I do not define _ode_interpolant for Tsit5_for_relaxation.

that is the reason.

[oscardssmith]

The fact that this requires ~1000 lines to add relaxation to a single algorithm seems like it should be improved. It seems like we should be able to duplicate less of the steps between the regular solvers and the relaxation ones. Currently, it seems like the chance of us making improvements to the regular and those changes not getting reflected in the relaxation solver is pretty high.

[Theozeud]

The fact that this requires ~1000 lines to add relaxation to a single algorithm seems like it should be improved. It seems like we should be able to duplicate less of the steps between the regular solvers and the relaxation ones. Currently, it seems like the chance of us making improvements to the regular and those changes not getting reflected in the relaxation solver is pretty high.

Thanks for your review. In fact, it does not require as many lines at all. I have replicated the Tsit5 solver as Tsit5_for_relaxation so as not to "break" the initial solver and make comparisons. But what I had in mind would be to modify the perform_step! function for all the methods for which we can do relaxations, into a general perform_step! function for all the methods. Instead, each method would have a computation! function and other functions such as the calculation of error estimates. I am aware that this is a big change to the general solver, but I can't see any other way of including relaxation before the loopfooter! function. (So we wouldn't be able to use all the variants of relaxation, including the R-FSAL method). The other solution would be to add a condtion like "if relaxation do relaxation end" inside each method, with a lot of copy and paste, but this does not seem viable for me in the long term, for example if there are other different methods than relaxations. I therefore propose this structure, which of course may not be desirable for reasons I am unaware of due to lack of experience, so I wanted to discuss it with you before going any further.

In the 1000 lines, there are several test files that I have made so that you can test them and see that the code works with relaxation, as well as a jupiter notebook that I should remove that shows that the code is not slowed down by the new perform_step! structure. For the tests, I have taken the examples I found in the https://arxiv.org/abs/2311.14050 article.