Skip to content

Change initialization of tolerances of fixed point iterations #95

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.

Sign up for GitHub

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

Merged
merged 4 commits into from
Jan 20, 2019
Merged

Change initialization of tolerances of fixed point iterations #95

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
81a2ec6
Change initialization of tolerances of fixed point iterations
devmotion Jan 20, 2019
d0802ff
Test default tolerances
devmotion Jan 20, 2019
61ce257
Fix typo
ChrisRackauckas Jan 20, 2019
ea2d2a3
Copy default tolerances
devmotion Jan 20, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
37 changes: 22 additions & 15 deletions src/solve.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,32 +92,39 @@ function DiffEqBase.__init(
mass_matrix = prob.f.mass_matrix)
end

# absolut tolerance for fixed-point iterations has to be of same type as elements of u
# in particular important for calculations with units
if typeof(alg.fixedpoint_abstol) <: Nothing
fixedpoint_abstol_internal = map(eltype(uType), integrator.opts.abstol)
# define absolute tolerance for fixed-point iterations
if alg.fixedpoint_abstol === nothing
abstol = integrator.opts.abstol
if typeof(abstol) <: Number
fixedpoint_abstol_internal = abstol
else
Copy link
Member Author

Choose a reason for hiding this comment

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

I'm not sure whether this specialization is required. To me it seems recursivecopy(::Number) falls back to deepcopy which both give a more complex output with e.g. @code_native recursivecopy(3) and @code_native deepcopy(3) than @code_native copy(3)? Should we just define recursivecopy(x::Number) = copy(x) in RecursiveArrayTools which by default just returns x (https://github.com/JuliaLang/julia/blob/master/base/number.jl#L88)?

Copy link
Member

Choose a reason for hiding this comment

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

Yes. There shouldn't be a need to deepcopy a number anymore. This could improve performance in some places as well.

Copy link
Member

Choose a reason for hiding this comment

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

We should make sure Measurements.jl still works after changing to copy. It should.

fixedpoint_abstol_internal = recursivecopy(abstol)
end
else
fixedpoint_abstol_internal = map(eltype(uType), alg.fixedpoint_abstol)
fixedpoint_abstol_internal = real.(alg.fixedpoint_abstol)
end

# use norm of ODE integrator if no norm for fixed-point iterations is specified
if typeof(alg.fixedpoint_norm) <: Nothing
fixedpoint_norm = integrator.opts.internalnorm
end

# derive unitless types
uEltypeNoUnits = typeof(recursive_one(integrator.u))
tTypeNoUnits = typeof(recursive_one(prob.tspan[1]))

# relative tolerance for fixed-point iterations has to be of same type as elements of u
# without units
# in particular important for calculations with units
if typeof(alg.fixedpoint_reltol) <: Nothing
fixedpoint_reltol_internal = convert.(eltype(uType), integrator.opts.reltol)
# define relative tolerance for fixed-point iterations
if alg.fixedpoint_reltol === nothing
reltol = integrator.opts.reltol
if typeof(reltol) <: Number
fixedpoint_reltol_internal = reltol
else
fixedpoint_reltol_internal = recursivecopy(reltol)
end
else
fixedpoint_reltol_internal = convert.(eltype(uType), oneunit.(integrator.u).*alg.fixedpoint_reltol)
fixedpoint_reltol_internal = real.(alg.fixedpoint_reltol)
end

# derive unitless types
uEltypeNoUnits = recursive_unitless_eltype(prob.u0)
tTypeNoUnits = typeof(one(tType))

# create separate copies u and uprev, not pointing to integrator.u or integrator.uprev,
# containers for residuals and to cache uprev with correct dimensions and types
# in particular for calculations with units residuals have to be unitless
Expand Down
32 changes: 18 additions & 14 deletions test/units.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,40 +6,44 @@ using Unitful
DiffEqProblemLibrary.DDEProblemLibrary.build_prob_dde_1delay_long_scalar_notinplace(1.0u"N", 1.0u"s")
prob_inplace = DiffEqProblemLibrary.DDEProblemLibrary.build_prob_dde_1delay_long(1.0u"N", 1.0u"s")

@testset for prob in (prob_notinplace, prob_inplace), constrained in (false, true)
@testset for prob in (prob_notinplace, prob_inplace)
@testset "correct" begin
# without units
alg1 = MethodOfSteps(Tsit5(), constrained=constrained, max_fixedpoint_iters=100,
fixedpoint_abstol=1e-10, fixedpoint_reltol=1e-4)
# default tolerances
alg1 = MethodOfSteps(Tsit5(), constrained=false, max_fixedpoint_iters=100)
sol1 = solve(prob, alg1)

# with correct units
alg2 = MethodOfSteps(Tsit5(), constrained=constrained, max_fixedpoint_iters=100,
fixedpoint_abstol=1e-7u"mN", fixedpoint_reltol=1e-4)
alg2 = MethodOfSteps(Tsit5(), constrained=false, max_fixedpoint_iters=100,
fixedpoint_abstol=1e-6u"N", fixedpoint_reltol=1e-3u"N")
sol2 = solve(prob, alg2)

@test sol1.t == sol2.t && sol1.u == sol2.u

# with correct units as vectors
if typeof(prob.u0) <: AbstractArray
alg3 = MethodOfSteps(Tsit5(), constrained=constrained, max_fixedpoint_iters=100,
fixedpoint_abstol=[1e-7u"mN"], fixedpoint_reltol=[1e-4])
alg3 = MethodOfSteps(Tsit5(), constrained=false, max_fixedpoint_iters=100,
fixedpoint_abstol=[1e-6u"N"], fixedpoint_reltol=[1e-3u"N"])
sol3 = solve(prob, alg3)

@test sol1.t == sol3.t && sol1.u == sol3.u
end
end

@testset "incorrect" begin
# with incorrect units for absolute tolerance
alg1 = MethodOfSteps(Tsit5(), constrained=constrained, max_fixedpoint_iters=100,
fixedpoint_abstol=1e-10u"s", fixedpoint_reltol=1e-4)
# without units
alg1 = MethodOfSteps(Tsit5(), constrained=false, max_fixedpoint_iters=100,
fixedpoint_abstol=1e-6, fixedpoint_reltol=1e-3)
@test_throws Unitful.DimensionError solve(prob, alg1)

# with incorrect units for relative tolerance
alg2 = MethodOfSteps(Tsit5(), constrained=constrained, max_fixedpoint_iters=100,
fixedpoint_abstol=1e-10, fixedpoint_reltol=1e-4u"N")
# with incorrect units for absolute tolerance
alg2 = MethodOfSteps(Tsit5(), constrained=false, max_fixedpoint_iters=100,
fixedpoint_abstol=1e-6u"s", fixedpoint_reltol=1e-3u"N")
@test_throws Unitful.DimensionError solve(prob, alg2)

# with incorrect units for relative tolerance
alg3 = MethodOfSteps(Tsit5(), constrained=false, max_fixedpoint_iters=100,
fixedpoint_abstol=1e-6u"N", fixedpoint_reltol=1e-3u"s")
@test_throws Unitful.DimensionError solve(prob, alg3)
end
end
end