Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[WIP] Relaxation Runge-Kutta #2227

Closed
wants to merge 46 commits into from
Closed

Conversation

Theozeud
Copy link
Contributor

@Theozeud Theozeud commented May 27, 2024

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 :

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

From an interface point of view :

  • See which variables are really needed to be added to the structure integrator,
  • 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,
  • Investigate all things needed to meticulously take care of time as matching tstops ect,
  • 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
Copy link
Member

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

@Theozeud
Copy link
Contributor Author

Theozeud commented May 30, 2024

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.

Copy link
Member

@ranocha ranocha left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

@JoshuaLampert
Copy link
Contributor

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
Copy link
Contributor Author

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
Copy link
Contributor Author

Theozeud commented Jun 11, 2024

Update on my work :

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

@Theozeud
Copy link
Contributor Author

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

image

@oscardssmith
Copy link
Contributor

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

that is the reason.

@@ -77,6 +77,8 @@ import DiffEqBase: resize!, deleteat!, addat!, full_cache, user_cache, u_cache,
add_saveat!, set_reltol!,
set_abstol!, postamble!, last_step_failed,
isautodifferentiable

export change_u!, change_dt!, apriori_bounds_dt
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

do these need to be exported?

thread::Thread = False()
end

# for backwards compatibility
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

backwards compatibility with what?

@@ -448,6 +448,44 @@ if isdefined(Base, :Experimental) && isdefined(Base.Experimental, :silence!)
Base.Experimental.silence!(Tsit5Cache)
end

@cache struct Tsit5Cache_for_relaxation{uType, rateType, uNoUnitsType, StageLimiter, StepLimiter,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

could this just use the Tsit5Cache?

integrator.dt_changed = dt
end

function change_u!(integrator::ODEIntegrator, u)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

should this set integrator.umodified?

end

# Calculate f(uₙ₊₁) if the aglortihm has the FSAL property
if isfsal(integrator.alg)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

won't this always be true?

@@ -0,0 +1,255 @@
{
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

shouldn't the benchmark be in SciMLBenchmarks rather than OrdinaryDiffEq?

@oscardssmith
Copy link
Contributor

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
Copy link
Contributor Author

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.

@oscardssmith
Copy link
Contributor

Closing in favor of #2283

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants