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

Adding Alshina2, Alshina3, Alshin6 from Alshina, Zaks, and Kalitkin (2008) #1898

Merged
merged 21 commits into from
Mar 31, 2023

Conversation

ayushinav
Copy link
Contributor

Adding the 2nd order from Alshina, Zaks, and Kalitkin (2008), Optimal first- to sixth-order accurate Runge-Kutta schemes.
This is not a new algo per se (trying to get familiar with the code structure), but if it works fine, I'll add the other orders soon.

@ayushinav ayushinav closed this Mar 11, 2023
@ayushinav ayushinav reopened this Mar 11, 2023
@ayushinav ayushinav marked this pull request as draft March 11, 2023 00:39
Copy link
Member

@ErikQQY ErikQQY left a comment

Choose a reason for hiding this comment

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

Here are some suggestions:

  1. Please make sure the new scheme passes the convergence tests to ensure the implementation is correct, for example, you can put this new scheme in a self-contained example to easily test the implementation. Also, you can learn from existing schemes or merged PRs to get familiar with the structure.
  2. Add convergence plots so that it would be convenient for committers to be sure this new scheme "behaves" well, see devdocs about the convergence plots.
  3. Since all of the SciML projects follow the SciML style, it would be nice if you can format your commit before pushing, see the CONTRIBUTING.md to know how.

@ayushinav ayushinav marked this pull request as ready for review March 14, 2023 07:39
@ayushinav
Copy link
Contributor Author

Alshina2
Alshina3
Alshina4
Alshina6

Thanks for the suggestions, @ErikQQY. Updated the codes and added the other methods.

@ayushinav ayushinav requested a review from ErikQQY March 14, 2023 19:48
step_limiter!::StepLimiter
thread::Thread
end
export Alshina2
Copy link
Member

Choose a reason for hiding this comment

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

Since we already export this algorithm in the OrdinaryDiffEq.jl file, this export is redundant.

step_limiter!::StepLimiter
thread::Thread
end
export Alshina3
Copy link
Member

Choose a reason for hiding this comment

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

Same here.

step_limiter!::StepLimiter
thread::Thread
end
export Alshina4
Copy link
Member

Choose a reason for hiding this comment

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

Here

step_limiter!::StepLimiter
thread::Thread
end
export Alshina6
Copy link
Member

Choose a reason for hiding this comment

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

And here


sim6 = test_convergence(dts1, prob, Alshina6())
@test sim6.𝒪est[:l2]≈4 atol=testTol

Copy link
Member

Choose a reason for hiding this comment

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

The convergence tests here look weird. What we are doing here is testing whether the order of the convergence estimate is equal to the order of this method. If it doesn't equal, you can adjust dts so that it doesn't saturate, if it still doesn't converge, then there are bugs in your implementation, you need to carefully compare your implementation with this paper, probably some coefficients or the initialize! and perform_step! are wrong.

Copy link
Contributor Author

@ayushinav ayushinav Mar 16, 2023

Choose a reason for hiding this comment

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

Right, my implementation had bugs. I was so confused by FSAL being used inspite of the algorithm not being FSAL'd. It's clear to me now.

The convergence tests make so much sense now. Also, am I correct in understanding that integrator.stats.nf stores the number of function calls?

@ayushinav ayushinav requested a review from ErikQQY March 16, 2023 01:22
Comment on lines 1938 to 1976
struct Alshina4ConstantCache{T, T2} <: OrdinaryDiffEqConstantCache
a21::T
a31::T
a32::T
a41::T
a42::T
a43::T


b1::T
b2::T
b3::T
b4::T

c2::T2
c3::T2
c4::T2

end

function Alshina4ConstantCache(T::Type{<:Float64}, T2::Type{<:Float64})
a21 = convert(T, 0.5)
a31 = convert(T, 0.0)
a32 = convert(T, 0.5)
a41 = convert(T, 0.0)
a42 = convert(T, 0.0)
a43 = convert(T, 1.0)

b1 = convert(T, 0.16666666666666666)
b2 = convert(T, 0.3333333333333333)
b3 = convert(T, 0.3333333333333333)
b4 = convert(T, 0.16666666666666666)

c2 = convert(T2, 0.5)
c3 = convert(T2, 0.5)
c4 = convert(T2, 1.0)

Alshina4ConstantCache(a21, a31, a32, a41, a42, a43, b1, b2, b3, b4, c2, c3, c4)
end
Copy link
Member

Choose a reason for hiding this comment

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

This is just RK4?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Do I need to resolve anything here? This is indeed just RK4.

Copy link
Member

Choose a reason for hiding this comment

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

Delete it because it's not useful.

Comment on lines 1898 to 1906
function Alshina2ConstantCache(T::Type{<:Float64}, T2::Type{<:Float64})
a21 = convert(T, 0.6666666666666666)

b1 = convert(T, 0.25)
b2 = convert(T, 0.75)

c2 = convert(T2, 0.6666666666666666)

Alshina2ConstantCache(a21, b1, b2, c2)
Copy link
Member

Choose a reason for hiding this comment

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

the non float64 dispatches are missing


function Alshina3ConstantCache(T::Type{<:Float64}, T2::Type{<:Float64})
a21 = convert(T, 0.5)
a31 = convert(T, 0.0)
Copy link
Member

Choose a reason for hiding this comment

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

remove the zeros from the calculation

Comment on lines 1968 to 1970
b3 = convert(T, 0.3333333333333333)
b4 = convert(T, 0.16666666666666666)

Copy link
Member

Choose a reason for hiding this comment

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

Embedded schemes are missing?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Ahh, I guess I did not include them because the paper had only for RK2 and RK3. My bad, will include them shortly.

@ayushinav ayushinav reopened this Mar 18, 2023
@ayushinav
Copy link
Contributor Author

Updated

@ayushinav
Copy link
Contributor Author

ayushinav commented Mar 21, 2023

Also, I guess there's a small bug in here, utilde is not brought in. I get the utilde not defined error when I run the in place lorenz fn using SIR54().

@ayushinav
Copy link
Contributor Author

ayushinav commented Mar 21, 2023

Added the regression tests as I understood. Please let me know if anything needs to be changed.

prob_linear= ODEProblemLibrary.prob_ode_linear

function lorenz!(du,u,p,t)
    du[1] = 10.0(u[2]-u[1])
    du[2] = u[1]*(28.0-u[3]) - u[2]
    du[3] = u[1]*u[2] - (8/3)*u[3]
end
u0 = [1.0;0.0;0.0]
tspan = (0.0,1.0)
prob_lorenz = ODEProblem(lorenz!,u0,tspan)

sol_linear2 = solve(prob_linear, Alshina2())
@show length(sol_linear2.u)

sol_lorenz2 = solve(prob_lorenz,Alshina2())
@show length(sol_lorenz2.t)

sol_linear3 = solve(prob_linear, Alshina3())
@show length(sol_linear3.u)

sol_lorenz3 = solve(prob_lorenz,Alshina3())
@show length(sol_lorenz3.t)

gives

length(sol_linear2.u) = 2040
length(sol_lorenz2.t) = 64110
length(sol_linear3.u) = 1286
length(sol_lorenz3.t) = 40524

@ayushinav ayushinav requested review from ChrisRackauckas and ErikQQY and removed request for ErikQQY and ChrisRackauckas March 22, 2023 02:41
@ChrisRackauckas
Copy link
Member

Test failure looks real.

@ayushinav
Copy link
Contributor Author

ayushinav commented Mar 22, 2023

I was not clear what the failures meant. Made the changes. Hope it looks right now.

@ayushinav ayushinav changed the title 2nd order of Alshina, Zaks, and Kalitkin (2008) Adding Alshina2, Alshina3, Alshin6 from Alshina, Zaks, and Kalitkin (2008) Mar 23, 2023
Comment on lines 131 to 138
function lorenz!(du,u,p,t)
du[1] = 10.0(u[2]-u[1])
du[2] = u[1]*(28.0-u[3]) - u[2]
du[3] = u[1]*u[2] - (8/3)*u[3]
end
u0 = [1.0;0.0;0.0]
tspan = (0.0,1.0)
prob_lorenz = ODEProblem(lorenz!,u0,tspan)
Copy link
Member

Choose a reason for hiding this comment

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

already defined above

Copy link
Contributor Author

Choose a reason for hiding this comment

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

True. The same function is defined on line 33. While I was trying to understand the tests, I realized that ESDIRK tests use the lorenz defined on line 59, which does not update the du and will always give 8 as the length of the soln, irrespective of the adaptive method used. Though ESDIRK pass the tests on right lorenz too, I did it to explicitly show what's happening.


sol_linear = solve(prob_linear, Alshina2())
sol_lorenz = solve(prob_lorenz, Alshina2())
@test length(sol_linear.u) < length(sol_lorenz.u)
Copy link
Member

Choose a reason for hiding this comment

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

I'm not sure this relative test makes sense.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I am not entirely surely sure about what to put in. Could you please point to what exactly I might be putting in? Thanks

Copy link
Member

Choose a reason for hiding this comment

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

Make it a regression test like the ones on linear above.

Copy link
Contributor Author

@ayushinav ayushinav Mar 26, 2023

Choose a reason for hiding this comment

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

Would this be fine then?

for prob in [prob_ode_2Dlinear, prob_ode_linear]
    sol = solve(prob, Alshina2())
    val= maximum(abs.(sol.u[end] - sol.u_analytic[end]))
    @test val<1e-6
    
    sol = solve(prob, Alshina3())
    val= maximum(abs.(sol.u[end] - sol.u_analytic[end]))
    @test val<1e-6
end

@ChrisRackauckas
Copy link
Member

Rebase onto master.

@ayushinav
Copy link
Contributor Author

Rebased again. Please let me know if I have to make any changes in here.

@ChrisRackauckas ChrisRackauckas merged commit c77b13d into SciML:master Mar 31, 2023
@ErikQQY
Copy link
Member

ErikQQY commented Apr 2, 2023

After adding these algorithms, it would be nice if you can record these new methods in solvers docs: DiffEqDocs.jl @ayushinav

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.

None yet

3 participants