WIP: implement fifth order Radau IIA #612

YingboMa · 2019-01-20T04:38:37Z

julia> using OrdinaryDiffEq

julia> function lorenz(u,p,t)
        du = similar(u)
        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]
        du
       end
lorenz (generic function with 1 method)

julia> solve(prob, Vern9(), dt=0.01, adaptive=false)[end]
3-element Array{Float64,1}:
 -5.857685382424004
 -5.831082486429051
 23.932132987022584

julia> solve(prob, RadauIIA5(), dt=0.01, adaptive=false)[end]
3-element Array{Float64,1}:
 -5.857676300291801
 -5.831062113159235
 23.932140720000515

coveralls · 2019-01-20T06:07:44Z

Coverage decreased (-9.1%) to 53.271% when pulling 8947511 on myb/radau into ef6b965 on master.

codecov · 2019-01-20T08:14:42Z

Codecov Report

Merging #612 into master will increase coverage by 0.17%.
The diff coverage is 82.35%.

@@            Coverage Diff             @@
##           master     #612      +/-   ##
==========================================
+ Coverage   68.45%   68.62%   +0.17%     
==========================================
  Files          83       86       +3     
  Lines       24721    24883     +162     
==========================================
+ Hits        16923    17077     +154     
- Misses       7798     7806       +8

Impacted Files	Coverage Δ
src/OrdinaryDiffEq.jl	`100% <ø> (ø)`	⬆️
src/alg_utils.jl	`30.63% <0%> (-0.21%)`	⬇️
src/algorithms.jl	`98.05% <100%> (+0.01%)`	⬆️
src/caches/firk_caches.jl	`100% <100%> (ø)`
src/integrators/integrator_utils.jl	`78.07% <54.54%> (+6.59%)`	⬆️
src/perform_step/firk_perform_step.jl	`81.73% <81.73%> (ø)`
src/tableaus/firk_tableaus.jl	`96.87% <96.87%> (ø)`
... and 1 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 ef6b965...8947511. Read the comment docs.

YingboMa · 2019-01-20T12:03:42Z

The adaptivity of the current implementation is very weird.

julia> using OrdinaryDiffEq, DiffEqDevTools, Test, LinearAlgebra

julia> using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()

julia> import DiffEqProblemLibrary.ODEProblemLibrary: van

julia> solve(remake(vanstiff, p=1e7), RadauIIA5(), abstol=1e-5, reltol=1e-7) |> length
668

julia> solve(remake(vanstiff, p=1e7), RadauIIA5()) |> length
┌ Warning: Interrupted. Larger maxiters is needed.
└ @ DiffEqBase ~/.julia/packages/DiffEqBase/8usQ9/src/integrator_interface.jl:150
95

julia> solve(remake(vanstiff, p=1e7), Kvaerno5()) |> length
110

julia> solve(remake(vanstiff, p=1e7), Kvaerno5(), abstol=1e-5, reltol=1e-7) |> length
401

src/alg_utils.jl

src/perform_step/firk_perform_step.jl

ChrisRackauckas · 2019-01-20T13:20:49Z

src/perform_step/firk_perform_step.jl

+
+  if adaptive
+    e1dt, e2dt, e3dt = e1/dt, e2/dt, e3/dt
+    tmp = @. e1dt*z1 + e2dt*z2 + e3dt*z3


is this not the err tilde term?

The err tilde term is integrator.fsalfirst + tmp.

https://github.com/luchr/ODEInterface.jl/blob/0bd134a5a358c4bc13e0fb6a90e27e4ee79e0115/src/dc_lapack.f#L1244-L1254

ChrisRackauckas · 2019-01-20T13:39:14Z

src/perform_step/firk_perform_step.jl

+    e1dt, e2dt, e3dt = e1/dt, e2/dt, e3/dt
+    tmp = @. e1dt*z1 + e2dt*z2 + e3dt*z3
+    mass_matrix != I && (tmp = mass_matrix*tmp)
+    err = @. integrator.fsalfirst + tmp


shouldn't this not be err tilde?

ChrisRackauckas · 2019-01-21T02:47:15Z

src/perform_step/firk_perform_step.jl

+  γdt, αdt, βdt = γ/dt, α/dt, β/dt
+  J = calc_J(integrator, cache, is_compos)
+  if u isa Number
+    LU1 = γdt*mass_matrix - J


this term should be mass_matrix/(gamma*dt), but it's not

wait, yes it is.

In all cases you should have gamma*dt. The transformation is to divide both sides by gdt = gamma*dt. The reason is because J is O(n^2) but I is O(n) (without mass matrices of course), so this decreases the overall calculations. This is exactly the same as the invW_t transformation on the Rosenbrock and SDIRK methods.

ChrisRackauckas · 2019-01-21T02:58:12Z

src/integrators/integrator_utils.jl

    qtmp = (integrator.EEst^expo)/fac
    @fastmath q = max(inv(integrator.opts.qmax),min(inv(integrator.opts.qmin),qtmp))
+    integrator.qold = q
    if q <= integrator.opts.qsteady_max && q >= integrator.opts.qsteady_min


This is supposed to be handling

But it looks slightly incorrect. It's really only used with a few recent algorithms:

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/alg_utils.jl#L312-L321

This might be the main difference between us and CVODE. Check with LSODE stuff and see if they are using qold here like Hairer.

ChrisRackauckas · 2019-01-21T03:00:07Z

src/perform_step/firk_perform_step.jl

+  u = @. uprev + z3
+
+  if adaptive
+    e1dt, e2dt, e3dt = e1/dt, e2/dt, e3/dt


These should all be divided by γdt

e1, e2, and e3 are already divided by γ0 in the tableau

Add basic oop RadauII5

8287fa7

YingboMa added 7 commits January 20, 2019 01:18

Add convergence test

66492d7

Differentiate (γ′, α′, β′) and (γ, α, β)

7ab70b8

Correct comments in RadauIIA5Tableau

41a62ba

Add precalculations and handle both scalar and array u

a9d2a2a

Add Newton iteration

73a9906

Add LTE calculation

457b4ee

Add adaptivity test for RadauIIA5

67e747f

YingboMa added 3 commits January 20, 2019 06:51

Add smooth_est option

ed2c06f

add nl_iters field in cache

668d6aa

Fix controller dispatch

26b2204

ChrisRackauckas reviewed Jan 20, 2019

View reviewed changes

src/alg_utils.jl Show resolved Hide resolved

ChrisRackauckas reviewed Jan 20, 2019

View reviewed changes

src/perform_step/firk_perform_step.jl Outdated Show resolved Hide resolved

ChrisRackauckas reviewed Jan 20, 2019

View reviewed changes

src/perform_step/firk_perform_step.jl Outdated Show resolved Hide resolved

ChrisRackauckas reviewed Jan 20, 2019

View reviewed changes

src/perform_step/firk_perform_step.jl Show resolved Hide resolved

ChrisRackauckas reviewed Jan 20, 2019

View reviewed changes

YingboMa force-pushed the myb/radau branch from bddf1fd to 7521c23 Compare January 20, 2019 13:48

Set RadauIIA5's alg_adaptive_order to 3

04daad1

YingboMa force-pushed the myb/radau branch from 7521c23 to 04daad1 Compare January 20, 2019 14:31

YingboMa added 3 commits January 20, 2019 10:46

update the adaptivity test

c58c1b6

Add better initial guess for RadauIIA5

7f4541e

Fix predictive controller

e17fc18

ChrisRackauckas reviewed Jan 21, 2019

View reviewed changes

Move "round to 1" to step accept controllers

70a8e08

YingboMa added 3 commits January 20, 2019 23:09

Rename RHS variables in RadauIIA's Newton iteration

9b8c8ee

Only take the last nl_iters

8dfdfa0

Merge branch 'master' into myb/radau

8947511

ChrisRackauckas merged commit 09f48d9 into master Jan 21, 2019

ChrisRackauckas deleted the myb/radau branch January 21, 2019 12:04

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

WIP: implement fifth order Radau IIA #612

WIP: implement fifth order Radau IIA #612

YingboMa commented Jan 20, 2019

coveralls commented Jan 20, 2019 •

edited

codecov bot commented Jan 20, 2019 •

edited

YingboMa commented Jan 20, 2019

ChrisRackauckas Jan 20, 2019

YingboMa Jan 20, 2019

YingboMa Jan 20, 2019

ChrisRackauckas Jan 20, 2019

ChrisRackauckas Jan 21, 2019

ChrisRackauckas Jan 21, 2019

ChrisRackauckas Jan 21, 2019

ChrisRackauckas Jan 21, 2019

ChrisRackauckas Jan 21, 2019

ChrisRackauckas Jan 21, 2019

ChrisRackauckas Jan 21, 2019

YingboMa Jan 21, 2019

WIP: implement fifth order Radau IIA #612

WIP: implement fifth order Radau IIA #612

Conversation

YingboMa commented Jan 20, 2019

coveralls commented Jan 20, 2019 • edited

codecov bot commented Jan 20, 2019 • edited

Codecov Report

YingboMa commented Jan 20, 2019

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

coveralls commented Jan 20, 2019 •

edited

codecov bot commented Jan 20, 2019 •

edited