Skip to content

Conversation

@utkarsh530
Copy link
Member

@utkarsh530 utkarsh530 commented Feb 14, 2020

Hi, I have currently implemented RK6(5) based on the paper (It does not account for zero dissipation yet, some basic changes needs to be done) here which is the basis of paper mentioned in #867. However, since this is my first attempt, I may have made some errors and sim tests are failing. Can anyone help where to possibly search for errors. Thanks.

@utkarsh530 utkarsh530 changed the title Zero Order Dissipation ERK6(5) Zero Dissipation ERK6(5) Feb 14, 2020
@coveralls
Copy link

coveralls commented Feb 14, 2020

Coverage Status

Coverage increased (+0.07%) to 86.797% when pulling c15a279 on utkarsh530:FRK65 into 46b4763 on JuliaDiffEq:master.

@ChrisRackauckas
Copy link
Member

What does the convergence plot look like?

@utkarsh530
Copy link
Member Author

1-D out of place
plot
2-D in place
plot1

@codecov
Copy link

codecov bot commented Feb 15, 2020

Codecov Report

Merging #1033 into master will decrease coverage by 41.11%.
The diff coverage is 0%.

Impacted file tree graph

@@             Coverage Diff             @@
##           master    #1033       +/-   ##
===========================================
- Coverage   76.74%   35.62%   -41.12%     
===========================================
  Files          95       95               
  Lines       30621    30044      -577     
===========================================
- Hits        23500    10704    -12796     
- Misses       7121    19340    +12219
Impacted Files Coverage Δ
src/OrdinaryDiffEq.jl 100% <ø> (ø) ⬆️
src/caches/low_order_rk_caches.jl 34.67% <0%> (-61.07%) ⬇️
src/algorithms.jl 44.25% <0%> (-39.63%) ⬇️
src/alg_utils.jl 7.39% <0%> (-20.1%) ⬇️
src/perform_step/low_order_rk_perform_step.jl 67.93% <0%> (-30.63%) ⬇️
src/caches/rkn_caches.jl 0% <0%> (-100%) ⬇️
src/iterator_interface.jl 0% <0%> (-100%) ⬇️
src/caches/symplectic_caches.jl 0% <0%> (-100%) ⬇️
src/caches/prk_caches.jl 0% <0%> (-100%) ⬇️
src/tableaus/rkc_tableaus.jl 0% <0%> (-100%) ⬇️
... and 66 more

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update b762a2d...c15a279. Read the comment docs.

@ChrisRackauckas
Copy link
Member

This isn't even converging. Do the b_i sum to 1? To the rows of a_i sum to c_i?

@unpack t, dt, uprev, u, f, p = integrator
@unpack α21, α31, α41, α51, α61, α71, α81, α91, α32, α43, α53, α63, α73, α83, α54, α64, α74, α84, α94, α65, α75, α85, α95, α76, α86, α96, α87, α97, α98, β1, β4, β5, β6, β7, β8, β1tilde, β4tilde, β5tilde, β6tilde, β7tilde, β8tilde, β9tilde, c2, c3, c4, c5, c6, c7, c8, c9 = cache
alg = unwrap_alg(integrator, true)
omega = alg.omega
Copy link
Member

Choose a reason for hiding this comment

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

looks like omega isn't used?

Copy link
Member Author

Choose a reason for hiding this comment

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

Actually, I am testing for convergence of paper I mentioned in my earlier comment, which had constant variables. The paper for zero dissipation requires only to change b4 b5 b6 to a function of omega which I will implement after this just to be sure that original coefficients are leading to convergence or not.

@kanav99
Copy link
Contributor

kanav99 commented Feb 15, 2020

Hey! The b_i coefficients were bad. Check the tests now with the suggestions

@utkarsh530
Copy link
Member Author

utkarsh530 commented Feb 15, 2020

Yep, I corrected another coefficient a_i & btilde_i. The 2-D in place plot looks like this:

plot3

Order estimate:
1-D out of place
Expression: ≈(sim3.�est[:l∞], 6, atol=testTol) Evaluated: 2.7735194214179764 ≈ 6 (atol=0.2)
2-D in place:
Expression: ≈(sim3.�est[:l∞], 6, atol=testTol) Evaluated: 2.865412492134244 ≈ 6 (atol=0.2)
Coefficients are for Periodic IVPs and hence result does not follow?

@ChrisRackauckas
Copy link
Member

Make the dts a little bit larger: it's saturating at floating point epsilon.

@utkarsh530
Copy link
Member Author

Yep, that worked, thanks! I will add code for zero dissipation coefficients.

@utkarsh530 utkarsh530 force-pushed the FRK65 branch 3 times, most recently from 945f1ae to 17c83f3 Compare February 21, 2020 06:40
@utkarsh530
Copy link
Member Author

I have added zero-dissipation fitted coefficients.
Question: should we set omega = 0 if not specified? (It coincides with RK6(5))

@ChrisRackauckas
Copy link
Member

Question: should we set omega = 0 if not specified? (It coincides with RK6(5))

sure

@utkarsh530
Copy link
Member Author

utkarsh530 commented Feb 21, 2020

Made required changes to struct FRK65. Please review.

lazy::Bool
Vern9(;lazy=true) = new(lazy)
end
struct FRK65 <: OrdinaryDiffEqAdaptiveAlgorithm
Copy link
Member

Choose a reason for hiding this comment

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

Parameterize this

lazy::Bool
Vern9(;lazy=true) = new(lazy)
end
struct FRK65 <: OrdinaryDiffEqAdaptiveAlgorithm
Copy link
Member

Choose a reason for hiding this comment

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

Suggested change
struct FRK65 <: OrdinaryDiffEqAdaptiveAlgorithm
struct FRK65{T} <: OrdinaryDiffEqAdaptiveAlgorithm

Vern9(;lazy=true) = new(lazy)
end
struct FRK65 <: OrdinaryDiffEqAdaptiveAlgorithm
omega::Float64
Copy link
Member

Choose a reason for hiding this comment

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

Suggested change
omega::Float64
omega::T

@utkarsh530
Copy link
Member Author

Made the required changes.

sim3 = test_convergence(dts4,prob,FRK65(0)) #standard value for non zero-dissipation tests
@test sim3.𝒪est[:l∞] 6 atol=testTol
sim3 = test_convergence(dts4,prob,FRK65())
@test sim3.𝒪est[:l∞] 6 atol=testTol1
Copy link
Member

Choose a reason for hiding this comment

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

instead, make this 6.5 and keep the same testtol. Looks like super convergence. What does plot(sim3) look like?

Copy link
Member Author

@utkarsh530 utkarsh530 Feb 21, 2020

Choose a reason for hiding this comment

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

Here's the plot:
plot
sim3.�est[:l∞]=6.547286289246398
plot1
sim3.�est[:l∞]=5.862512678209683
A test is failing if setting it to 6.5 and same testtol.

Copy link
Member

Choose a reason for hiding this comment

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

huh, that's weird. Instead of defining at testTol1, just put the larger value right in the code there and I'll accept.

Copy link
Member Author

Choose a reason for hiding this comment

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

Done. By the way, what is super convergence?

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.

4 participants