From 393a506549358b54a44bb6dd2a23ce71faf82f8e Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 17 Dec 2023 12:41:23 -0500 Subject: [PATCH 01/53] Update StaticArray tests for eltype requirements --- test/interface/static_array_tests.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/interface/static_array_tests.jl b/test/interface/static_array_tests.jl index cae40bc4f4..3622051ed7 100644 --- a/test/interface/static_array_tests.jl +++ b/test/interface/static_array_tests.jl @@ -2,7 +2,7 @@ using StaticArrays, Test using OrdinaryDiffEq using RecursiveArrayTools -u0 = [fill(2, MVector{2, Float64}), ones(MVector{2, Float64})] +u0 = VectorOfAray([fill(2, MVector{2, Float64}), ones(MVector{2, Float64})]) g(u, p, t) = SA[u[1] + u[2], u[1]] f = (du, u, p, t) -> begin for i in 1:2 @@ -17,7 +17,7 @@ sol = solve(ode, Tsit5()) sol = solve(ode, Vern9()) @test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0) -u0 = [fill(2, SVector{2, Float64}), ones(SVector{2, Float64})] +u0 = VectorOfAray([fill(2, SVector{2, Float64}), ones(SVector{2, Float64})]) ode = ODEProblem(f, u0, (0.0, 1.0)) sol = solve(ode, Euler(), dt = 1e-2) @test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0) @@ -28,14 +28,14 @@ sol = solve(ode, SSPRK22(), dt = 1e-2) sol = solve(ode, ROCK4()) @test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0) -u0 = ones(MVector{2, Float64}) +u0 = VectorOfAray(ones(MVector{2, Float64})) ode = ODEProblem(g, u0, (0.0, 1.0)) sol = solve(ode, Euler(), dt = 1e-2) @test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0) sol = solve(ode, Tsit5(), dt = 1e-2) @test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0) -u0 = ones(SVector{2, Float64}) +u0 = VectorOfAray(ones(SVector{2, Float64})) f = (u, p, t) -> u ode = ODEProblem(f, u0, (0.0, 1.0)) sol = solve(ode, Euler(), dt = 1e-2) From f35f470be08b690beeb34271012b97a362c89705 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 17 Dec 2023 13:12:52 -0500 Subject: [PATCH 02/53] typo --- test/interface/static_array_tests.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/interface/static_array_tests.jl b/test/interface/static_array_tests.jl index 3622051ed7..5ad6c07064 100644 --- a/test/interface/static_array_tests.jl +++ b/test/interface/static_array_tests.jl @@ -2,7 +2,7 @@ using StaticArrays, Test using OrdinaryDiffEq using RecursiveArrayTools -u0 = VectorOfAray([fill(2, MVector{2, Float64}), ones(MVector{2, Float64})]) +u0 = VectorOfArray([fill(2, MVector{2, Float64}), ones(MVector{2, Float64})]) g(u, p, t) = SA[u[1] + u[2], u[1]] f = (du, u, p, t) -> begin for i in 1:2 @@ -17,7 +17,7 @@ sol = solve(ode, Tsit5()) sol = solve(ode, Vern9()) @test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0) -u0 = VectorOfAray([fill(2, SVector{2, Float64}), ones(SVector{2, Float64})]) +u0 = VectorOfArray([fill(2, SVector{2, Float64}), ones(SVector{2, Float64})]) ode = ODEProblem(f, u0, (0.0, 1.0)) sol = solve(ode, Euler(), dt = 1e-2) @test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0) @@ -28,14 +28,14 @@ sol = solve(ode, SSPRK22(), dt = 1e-2) sol = solve(ode, ROCK4()) @test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0) -u0 = VectorOfAray(ones(MVector{2, Float64})) +u0 = VectorOfArray(ones(MVector{2, Float64})) ode = ODEProblem(g, u0, (0.0, 1.0)) sol = solve(ode, Euler(), dt = 1e-2) @test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0) sol = solve(ode, Tsit5(), dt = 1e-2) @test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0) -u0 = VectorOfAray(ones(SVector{2, Float64})) +u0 = VectorOfArray(ones(SVector{2, Float64})) f = (u, p, t) -> u ode = ODEProblem(f, u0, (0.0, 1.0)) sol = solve(ode, Euler(), dt = 1e-2) From 1ed5ccce54b54d77d0f7e17e68f0467d70d78cc5 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 18 Dec 2023 08:20:02 -0500 Subject: [PATCH 03/53] Update test/interface/static_array_tests.jl --- test/interface/static_array_tests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/interface/static_array_tests.jl b/test/interface/static_array_tests.jl index 5ad6c07064..c9e2216486 100644 --- a/test/interface/static_array_tests.jl +++ b/test/interface/static_array_tests.jl @@ -28,7 +28,7 @@ sol = solve(ode, SSPRK22(), dt = 1e-2) sol = solve(ode, ROCK4()) @test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0) -u0 = VectorOfArray(ones(MVector{2, Float64})) +u0 = ones(MVector{2, Float64}) ode = ODEProblem(g, u0, (0.0, 1.0)) sol = solve(ode, Euler(), dt = 1e-2) @test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0) From 8296bc3d2538f4eb85941c13928883965f430624 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 18 Dec 2023 08:20:16 -0500 Subject: [PATCH 04/53] Update test/interface/static_array_tests.jl --- test/interface/static_array_tests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/interface/static_array_tests.jl b/test/interface/static_array_tests.jl index c9e2216486..ff2bcbb3cf 100644 --- a/test/interface/static_array_tests.jl +++ b/test/interface/static_array_tests.jl @@ -35,7 +35,7 @@ sol = solve(ode, Euler(), dt = 1e-2) sol = solve(ode, Tsit5(), dt = 1e-2) @test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0) -u0 = VectorOfArray(ones(SVector{2, Float64})) +u0 = ones(SVector{2, Float64}) f = (u, p, t) -> u ode = ODEProblem(f, u0, (0.0, 1.0)) sol = solve(ode, Euler(), dt = 1e-2) From 6e20ceeff26f4c7ead4a84a3d0c595e0f132d59f Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 24 Dec 2023 14:42:22 -0500 Subject: [PATCH 05/53] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 5b5eb4085c..3f310a7194 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "OrdinaryDiffEq" uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" authors = ["Chris Rackauckas ", "Yingbo Ma "] -version = "6.63.0" +version = "6.64.0" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From 092be52cf4d6547d19ad3577419111fc18621262 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 26 Dec 2023 16:01:17 -0500 Subject: [PATCH 06/53] Test master --- README.md | 1 + 1 file changed, 1 insertion(+) diff --git a/README.md b/README.md index b4e9bf8ad4..7b07bdbe0c 100644 --- a/README.md +++ b/README.md @@ -15,6 +15,7 @@ ordinary differential equation solvers and utilities. While completely independe and usable on its own, users interested in using this functionality should check out [DifferentialEquations.jl](https://github.com/SciML/DifferentialEquations.jl). + ## Installation Assuming that you already have Julia correctly installed, it suffices to import From 0306da0ef4b9036f43e154cb2255ba75820e917a Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 26 Dec 2023 16:34:20 -0500 Subject: [PATCH 07/53] don't differentiate the convergence check in the nonlinear solve --- src/nlsolve/nlsolve.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/nlsolve/nlsolve.jl b/src/nlsolve/nlsolve.jl index b522239ae3..140fbcec7d 100644 --- a/src/nlsolve/nlsolve.jl +++ b/src/nlsolve/nlsolve.jl @@ -89,7 +89,7 @@ function nlsolve!(nlsolver::AbstractNLSolver, integrator::DiffEqBase.DEIntegrato apply_step!(nlsolver, integrator) # check for convergence - iter > 1 && (η = θ / (1 - θ)) + iter > 1 && (η = DiffEqBase.value(θ / (1 - θ))) if (iter == 1 && ndz < 1e-5) || (iter > 1 && (η >= zero(η) && η * ndz < κ)) nlsolver.status = Convergence nlsolver.nfails = 0 From 4749c38351933e2c5d085ff7f6a4f6d2506b1034 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 26 Dec 2023 17:03:11 -0500 Subject: [PATCH 08/53] `@test_broken` on inference is changed in v1.1.0 --- test/regression/inference.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/test/regression/inference.jl b/test/regression/inference.jl index dc28b6bc33..453c4037c8 100644 --- a/test/regression/inference.jl +++ b/test/regression/inference.jl @@ -7,13 +7,13 @@ using Test inferred = [BS3(), Tsit5(), RK4(), Vern6()] for alg in inferred - @test_broken init(prob, alg) - @test_broken init(prob2D, alg) + @inferred init(prob, alg) + @inferred init(prob2D, alg) end - notinferred = [SDIRK2(), TRBDF2(), KenCarp4(), Rosenbrock23(), Rodas4()] - for alg in notinferred - @test_broken @inferred init(prob, alg) - @test_broken @inferred init(prob2D, alg) - end + #notinferred = [SDIRK2(), TRBDF2(), KenCarp4(), Rosenbrock23(), Rodas4()] + #for alg in notinferred + # @test_broken @inferred init(prob, alg).t[1] == 0.0 + # @test_broken @inferred init(prob2D, alg).t[1] == 0.0 + #end end From cecce4008e393f2e060af58a8ad1be51f0f440f9 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 26 Dec 2023 17:03:19 -0500 Subject: [PATCH 09/53] simplify --- test/regression/inference.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/regression/inference.jl b/test/regression/inference.jl index 453c4037c8..9e62015f54 100644 --- a/test/regression/inference.jl +++ b/test/regression/inference.jl @@ -13,7 +13,7 @@ using Test #notinferred = [SDIRK2(), TRBDF2(), KenCarp4(), Rosenbrock23(), Rodas4()] #for alg in notinferred - # @test_broken @inferred init(prob, alg).t[1] == 0.0 - # @test_broken @inferred init(prob2D, alg).t[1] == 0.0 + # @test_broken @inferred init(prob, alg) + # @test_broken @inferred init(prob2D, alg) #end end From e91367e7a2f402d0a3d0c4cf69475018bdeb03ac Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 26 Dec 2023 18:47:03 -0500 Subject: [PATCH 10/53] Workaround v1.10 function redefinition issue https://github.com/JuliaLang/julia/issues/52635 --- src/derivative_wrappers.jl | 4 ++-- test/interface/autodiff_error_tests.jl | 3 ++- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/derivative_wrappers.jl b/src/derivative_wrappers.jl index 75f3484407..0123e2e578 100644 --- a/src/derivative_wrappers.jl +++ b/src/derivative_wrappers.jl @@ -6,7 +6,7 @@ const FIRST_AUTODIFF_TGRAD_MESSAGE = """ 1. Turn off automatic differentiation (e.g. Rosenbrock23() becomes Rosenbrock23(autodiff=false)). More details can be found at https://docs.sciml.ai/DiffEqDocs/stable/features/performance_overloads/ - 2. Improving the compatibility of `f` with ForwardDiff.jl automatic + 2. Improving the compatibility of `f` with ForwardDiff.jl automatic differentiation (using tools like PreallocationTools.jl). More details can be found at https://docs.sciml.ai/DiffEqDocs/stable/basics/faq/#Autodifferentiation-and-Dual-Numbers 3. Defining analytical Jacobians and time gradients. More details can be @@ -48,7 +48,7 @@ const FIRST_AUTODIFF_JAC_MESSAGE = """ 1. Turn off automatic differentiation (e.g. Rosenbrock23() becomes Rosenbrock23(autodiff=false)). More details can befound at https://docs.sciml.ai/DiffEqDocs/stable/features/performance_overloads/ - 2. Improving the compatibility of `f` with ForwardDiff.jl automatic + 2. Improving the compatibility of `f` with ForwardDiff.jl automatic differentiation (using tools like PreallocationTools.jl). More details can be found at https://docs.sciml.ai/DiffEqDocs/stable/basics/faq/#Autodifferentiation-and-Dual-Numbers 3. Defining analytical Jacobians. More details can be diff --git a/test/interface/autodiff_error_tests.jl b/test/interface/autodiff_error_tests.jl index c87589a428..e374166b54 100644 --- a/test/interface/autodiff_error_tests.jl +++ b/test/interface/autodiff_error_tests.jl @@ -34,12 +34,13 @@ tspan = (0.0, 1.0) prob = ODEProblem(lorenz!, u0, tspan) @test_throws OrdinaryDiffEq.FirstAutodiffJacError solve(prob, Rosenbrock23()) -function lorenz!(du, u, p, t) +function lorenz2!(du, u, p, t) du[1] = 10.0(u[2] - u[1]) a[1] = t du[2] = u[1] * (28.0 - u[3]) - u[2] du[3] = u[1] * u[2] - (8 / 3) * u[3] end +prob = ODEProblem(lorenz2!, u0, tspan) @test_throws OrdinaryDiffEq.FirstAutodiffTgradError solve(prob, Rosenbrock23()) ## Test that nothing is using duals when autodiff=false From 44064304dcb43ae1c108f45ce00e7283d8a1829a Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 26 Dec 2023 19:44:35 -0500 Subject: [PATCH 11/53] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 3f310a7194..ca3a7c2b4e 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "OrdinaryDiffEq" uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" authors = ["Chris Rackauckas ", "Yingbo Ma "] -version = "6.64.0" +version = "6.65.0" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From 44a9d078270ae68a438ad754e90f09d58ab98b47 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 27 Dec 2023 02:16:37 -0500 Subject: [PATCH 12/53] Fix interpolation output types for dynamical ODEs Fixes https://github.com/SciML/OrdinaryDiffEq.jl/issues/2086 --- src/dense/generic_dense.jl | 42 ++++++++++++++++---- test/interface/interpolation_output_types.jl | 22 ++++++++++ test/runtests.jl | 1 + 3 files changed, 58 insertions(+), 7 deletions(-) create mode 100644 test/interface/interpolation_output_types.jl diff --git a/src/dense/generic_dense.jl b/src/dense/generic_dense.jl index 22fc4860ba..68a490b8f7 100644 --- a/src/dense/generic_dense.jl +++ b/src/dense/generic_dense.jl @@ -684,6 +684,18 @@ Hairer Norsett Wanner Solving Ordinary Differential Euations I - Nonstiff Proble Herimte Interpolation, chosen if no other dispatch for ode_interpolant """ +@muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{false}}, idxs::Nothing, + T::Type{Val{0}}, differential_vars) # Default interpolant is Hermite + #@.. broadcast=false (1-Θ)*y₀+Θ*y₁+Θ*(Θ-1)*((1-2Θ)*(y₁-y₀)+(Θ-1)*dt*k[1] + Θ*dt*k[2]) + if all(differential_vars) + @inbounds (1 - Θ) * y₀ + Θ * y₁ + + (Θ * (Θ - 1) * ((1 - 2Θ) * (y₁ - y₀) + (Θ - 1) * dt * k[1] + Θ * dt * k[2])) + else + @inbounds (1 - Θ) * y₀ + Θ * y₁ + + differential_vars .* (Θ * (Θ - 1) * ((1 - 2Θ) * (y₁ - y₀) + (Θ - 1) * dt * k[1] + Θ * dt * k[2])) + end +end + @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{false}}, idxs::Nothing, T::Type{Val{0}}, differential_vars) # Default interpolant is Hermite #@.. broadcast=false (1-Θ)*y₀+Θ*y₁+Θ*(Θ-1)*((1-2Θ)*(y₁-y₀)+(Θ-1)*dt*k[1] + Θ*dt*k[2]) @@ -755,10 +767,17 @@ Herimte Interpolation, chosen if no other dispatch for ode_interpolant @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{false}}, idxs::Nothing, T::Type{Val{1}}, differential_vars) # Default interpolant is Hermite #@.. broadcast=false k[1] + Θ*(-4*dt*k[1] - 2*dt*k[2] - 6*y₀ + Θ*(3*dt*k[1] + 3*dt*k[2] + 6*y₀ - 6*y₁) + 6*y₁)/dt - @inbounds (.!differential_vars).*(y₁ - y₀)/dt + differential_vars .*( - k[1] + - Θ * (-4 * dt * k[1] - 2 * dt * k[2] - 6 * y₀ + - Θ * (3 * dt * k[1] + 3 * dt * k[2] + 6 * y₀ - 6 * y₁) + 6 * y₁) / dt) + if all(differential_vars) + @inbounds (y₁ - y₀)/dt + ( + k[1] + + Θ * (-4 * dt * k[1] - 2 * dt * k[2] - 6 * y₀ + + Θ * (3 * dt * k[1] + 3 * dt * k[2] + 6 * y₀ - 6 * y₁) + 6 * y₁) / dt) + else + @inbounds (.!differential_vars).*(y₁ - y₀)/dt + differential_vars .*( + k[1] + + Θ * (-4 * dt * k[1] - 2 * dt * k[2] - 6 * y₀ + + Θ * (3 * dt * k[1] + 3 * dt * k[2] + 6 * y₀ - 6 * y₁) + 6 * y₁) / dt) + end end @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{true}}, idxs::Nothing, @@ -826,8 +845,13 @@ Herimte Interpolation, chosen if no other dispatch for ode_interpolant @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{false}}, idxs::Nothing, T::Type{Val{2}}, differential_vars) # Default interpolant is Hermite #@.. broadcast=false (-4*dt*k[1] - 2*dt*k[2] - 6*y₀ + Θ*(6*dt*k[1] + 6*dt*k[2] + 12*y₀ - 12*y₁) + 6*y₁)/(dt*dt) - @inbounds differential_vars .* (-4 * dt * k[1] - 2 * dt * k[2] - 6 * y₀ + - Θ * (6 * dt * k[1] + 6 * dt * k[2] + 12 * y₀ - 12 * y₁) + 6 * y₁) / (dt * dt) + if all(differential_vars) + @inbounds (-4 * dt * k[1] - 2 * dt * k[2] - 6 * y₀ + + Θ * (6 * dt * k[1] + 6 * dt * k[2] + 12 * y₀ - 12 * y₁) + 6 * y₁) / (dt * dt) + else + @inbounds differential_vars .* (-4 * dt * k[1] - 2 * dt * k[2] - 6 * y₀ + + Θ * (6 * dt * k[1] + 6 * dt * k[2] + 12 * y₀ - 12 * y₁) + 6 * y₁) / (dt * dt) + end end @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{true}}, idxs::Nothing, @@ -887,7 +911,11 @@ Herimte Interpolation, chosen if no other dispatch for ode_interpolant @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{false}}, idxs::Nothing, T::Type{Val{3}}, differential_vars) # Default interpolant is Hermite #@.. broadcast=false (6*dt*k[1] + 6*dt*k[2] + 12*y₀ - 12*y₁)/(dt*dt*dt) - @inbounds differential_vars .* (6 * dt * k[1] + 6 * dt * k[2] + 12 * y₀ - 12 * y₁) / (dt * dt * dt) + if all(differential_vars) + @inbounds (6 * dt * k[1] + 6 * dt * k[2] + 12 * y₀ - 12 * y₁) / (dt * dt * dt) + else + @inbounds differential_vars .* (6 * dt * k[1] + 6 * dt * k[2] + 12 * y₀ - 12 * y₁) / (dt * dt * dt) + end end @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{true}}, idxs::Nothing, diff --git a/test/interface/interpolation_output_types.jl b/test/interface/interpolation_output_types.jl new file mode 100644 index 0000000000..ea6e9cac1c --- /dev/null +++ b/test/interface/interpolation_output_types.jl @@ -0,0 +1,22 @@ +using OrdinaryDiffEq, Test + +# in terms of the voltage across all three elements +rlc1!(v′,v,(R,L,C),t) = -(v′/R + v/L)/C +identity_f(v,u,p,t) = v # needed to form second order dynamical ODE + +setup_rlc(R,L,C;v_init=0.0,v′_init=0.0,tspan=(0.0,50.0)) = + DynamicalODEProblem{false}(rlc1!,identity_f,v′_init,v_init,tspan,(R,L,C)) + +# simulate voltage impulse +R,L,C = 10, 0.3, 2 + +prob = setup_rlc(R,L,C,v_init=2.0) + +res1 = solve(prob,Vern8(),dt=1/10,saveat=1/10) # success +res3 = solve(prob,CalvoSanz4(),dt=1/10,saveat=1/10) # fail + +sol = solve(prob,CalvoSanz4(),dt=1/10) +@test sol(0.32) isa OrdinaryDiffEq.ArrayPartition +@test sol(0.32, Val{1}) isa OrdinaryDiffEq.ArrayPartition +@test sol(0.32, Val{2}) isa OrdinaryDiffEq.ArrayPartition +@test sol(0.32, Val{3}) isa OrdinaryDiffEq.ArrayPartition diff --git a/test/runtests.jl b/test/runtests.jl index 395c139a84..a3b27ccfbf 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -44,6 +44,7 @@ end @time @safetestset "Complex Tests" include("interface/complex_tests.jl") @time @safetestset "Ndim Complex Tests" include("interface/ode_ndim_complex_tests.jl") @time @safetestset "Number Type Tests" include("interface/ode_numbertype_tests.jl") + @time @safetestset "Interpolation Output Type Tests" include("interface/interpolation_output_types.jl") @time @safetestset "Stiffness Detection Tests" include("interface/stiffness_detection_test.jl") @time @safetestset "Composite Interpolation Tests" include("interface/composite_interpolation.jl") @time @safetestset "Export tests" include("interface/export_tests.jl") From 9cefb66fc14e68c6418d4f4b6ac228bfdb430ccb Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 27 Dec 2023 02:35:28 -0500 Subject: [PATCH 13/53] Fix duplicate function --- src/dense/generic_dense.jl | 7 ------- test/interface/interpolation_output_types.jl | 4 ++-- 2 files changed, 2 insertions(+), 9 deletions(-) diff --git a/src/dense/generic_dense.jl b/src/dense/generic_dense.jl index 68a490b8f7..1047cea661 100644 --- a/src/dense/generic_dense.jl +++ b/src/dense/generic_dense.jl @@ -696,13 +696,6 @@ Herimte Interpolation, chosen if no other dispatch for ode_interpolant end end -@muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{false}}, idxs::Nothing, - T::Type{Val{0}}, differential_vars) # Default interpolant is Hermite - #@.. broadcast=false (1-Θ)*y₀+Θ*y₁+Θ*(Θ-1)*((1-2Θ)*(y₁-y₀)+(Θ-1)*dt*k[1] + Θ*dt*k[2]) - @inbounds (1 - Θ) * y₀ + Θ * y₁ + - differential_vars .* (Θ * (Θ - 1) * ((1 - 2Θ) * (y₁ - y₀) + (Θ - 1) * dt * k[1] + Θ * dt * k[2])) -end - @muladd function hermite_interpolant(Θ, dt, y₀, y₁, k, ::Type{Val{true}}, idxs::Nothing, T::Type{Val{0}}, differential_vars) # Default interpolant is Hermite #@.. broadcast=false (1-Θ)*y₀+Θ*y₁+Θ*(Θ-1)*((1-2Θ)*(y₁-y₀)+(Θ-1)*dt*k[1] + Θ*dt*k[2]) diff --git a/test/interface/interpolation_output_types.jl b/test/interface/interpolation_output_types.jl index ea6e9cac1c..66f0599baf 100644 --- a/test/interface/interpolation_output_types.jl +++ b/test/interface/interpolation_output_types.jl @@ -12,8 +12,8 @@ R,L,C = 10, 0.3, 2 prob = setup_rlc(R,L,C,v_init=2.0) -res1 = solve(prob,Vern8(),dt=1/10,saveat=1/10) # success -res3 = solve(prob,CalvoSanz4(),dt=1/10,saveat=1/10) # fail +res1 = solve(prob,Vern8(),dt=1/10,saveat=1/10) +res3 = solve(prob,CalvoSanz4(),dt=1/10,saveat=1/10) sol = solve(prob,CalvoSanz4(),dt=1/10) @test sol(0.32) isa OrdinaryDiffEq.ArrayPartition From dbb3e02dab28d85e90cb77695b479ff0d08b8050 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 27 Dec 2023 03:07:57 -0500 Subject: [PATCH 14/53] Update src/dense/generic_dense.jl --- src/dense/generic_dense.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dense/generic_dense.jl b/src/dense/generic_dense.jl index 1047cea661..a943579eb2 100644 --- a/src/dense/generic_dense.jl +++ b/src/dense/generic_dense.jl @@ -761,7 +761,7 @@ Herimte Interpolation, chosen if no other dispatch for ode_interpolant T::Type{Val{1}}, differential_vars) # Default interpolant is Hermite #@.. broadcast=false k[1] + Θ*(-4*dt*k[1] - 2*dt*k[2] - 6*y₀ + Θ*(3*dt*k[1] + 3*dt*k[2] + 6*y₀ - 6*y₁) + 6*y₁)/dt if all(differential_vars) - @inbounds (y₁ - y₀)/dt + ( + @inbounds ( k[1] + Θ * (-4 * dt * k[1] - 2 * dt * k[2] - 6 * y₀ + Θ * (3 * dt * k[1] + 3 * dt * k[2] + 6 * y₀ - 6 * y₁) + 6 * y₁) / dt) From 40dfbd10f7181b716d3d1c9bd79a448e747454a9 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 28 Dec 2023 08:10:45 -0500 Subject: [PATCH 15/53] Fix save_end overriding behavior Fixes https://github.com/SciML/OrdinaryDiffEq.jl/issues/1842 --- src/integrators/integrator_utils.jl | 5 +++++ src/solve.jl | 13 ++++++++++--- test/interface/ode_saveat_tests.jl | 25 +++++++++++++++++++++++++ 3 files changed, 40 insertions(+), 3 deletions(-) diff --git a/src/integrators/integrator_utils.jl b/src/integrators/integrator_utils.jl index 4bc032d1a6..da1ef70034 100644 --- a/src/integrators/integrator_utils.jl +++ b/src/integrators/integrator_utils.jl @@ -79,6 +79,10 @@ function _savevalues!(integrator, force_save, reduce_size)::Tuple{Bool, Bool} integrator.cache.current) end else # ==t, just save + if curt == integrator.sol.prob.tspan[2] && !integrator.opts.save_end + integrator.saveiter -= 1 + continue + end savedexactly = true copyat_or_push!(integrator.sol.t, integrator.saveiter, integrator.t) if integrator.opts.save_idxs === nothing @@ -145,6 +149,7 @@ postamble!(integrator::ODEIntegrator) = _postamble!(integrator) function _postamble!(integrator) DiffEqBase.finalize!(integrator.opts.callback, integrator.u, integrator.t, integrator) solution_endpoint_match_cur_integrator!(integrator) + save resize!(integrator.sol.t, integrator.saveiter) resize!(integrator.sol.u, integrator.saveiter) if !(integrator.sol isa DAESolution) diff --git a/src/solve.jl b/src/solve.jl index acb6ffee54..0661c72434 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -293,9 +293,16 @@ function DiffEqBase.__init(prob::Union{DiffEqBase.AbstractODEProblem, sizehint!(ts, 50) sizehint!(ks, 50) elseif !isempty(saveat_internal) - sizehint!(timeseries, length(saveat_internal) + 1) - sizehint!(ts, length(saveat_internal) + 1) - sizehint!(ks, length(saveat_internal) + 1) + savelength = length(saveat_internal) + 1 + if save_start == false + savelength -= 1 + end + if save_end == false && prob.tspan[2] in saveat_internal.valtree + savelength -= 1 + end + sizehint!(timeseries, savelength) + sizehint!(ts, savelength) + sizehint!(ks, savelength) else sizehint!(timeseries, 2) sizehint!(ts, 2) diff --git a/test/interface/ode_saveat_tests.jl b/test/interface/ode_saveat_tests.jl index 858fada876..52df1a9276 100644 --- a/test/interface/ode_saveat_tests.jl +++ b/test/interface/ode_saveat_tests.jl @@ -187,3 +187,28 @@ prob = ODEProblem(SIR!, [0.99, 0.01, 0.0], (t_obs[1], t_obs[end]), [0.20, 0.15]) sol = solve(prob, DP5(), reltol = 1e-6, abstol = 1e-6, saveat = t_obs) @test maximum(sol) <= 1 @test minimum(sol) >= 0 + +@testset "Proper save_start and save_end behavior" begin + function f2(du, u, p, t) + du[1] = -cos(u[1]) * u[1] + end + prob = ODEProblem(f2, [10], (0.0, 0.4)) + + @test solve(prob, Tsit5(); saveat = 0:.1:.4).t == [0.0; 0.1; 0.2; 0.3; 0.4] + @test solve(prob, Tsit5(); saveat = 0:.1:.4, save_start = true, save_end = true).t == [0.0; 0.1; 0.2; 0.3; 0.4] + @test solve(prob, Tsit5(); saveat = 0:.1:.4, save_start = false, save_end = false).t == [0.1; 0.2; 0.3] + + ts = solve(prob, Tsit5()).t + @test 0.0 in ts + @test 0.4 in ts + ts = solve(prob, Tsit5(); save_start = true, save_end = true).t + @test 0.0 in ts + @test 0.4 in ts + ts = solve(prob, Tsit5(); save_start = false, save_end = false).t + @test 0.0 ∉ ts + @test 0.4 ∉ ts + + @test solve(prob, Tsit5(); saveat = [.2]).t == [0.2] + @test solve(prob, Tsit5(); saveat = [.2], save_start = true, save_end = true).t == [0.0; 0.2; 0.4] + @test solve(prob, Tsit5(); saveat = [.2], save_start = false, save_end = false).t == [0.2] +end \ No newline at end of file From 8643fdbb670f7aa51929432de7f66063dacbebe3 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 28 Dec 2023 08:14:07 -0500 Subject: [PATCH 16/53] Update src/integrators/integrator_utils.jl --- src/integrators/integrator_utils.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/integrators/integrator_utils.jl b/src/integrators/integrator_utils.jl index da1ef70034..8edec505a3 100644 --- a/src/integrators/integrator_utils.jl +++ b/src/integrators/integrator_utils.jl @@ -149,7 +149,6 @@ postamble!(integrator::ODEIntegrator) = _postamble!(integrator) function _postamble!(integrator) DiffEqBase.finalize!(integrator.opts.callback, integrator.u, integrator.t, integrator) solution_endpoint_match_cur_integrator!(integrator) - save resize!(integrator.sol.t, integrator.saveiter) resize!(integrator.sol.u, integrator.saveiter) if !(integrator.sol isa DAESolution) From f3751847255d069eb5fdcbccabb95150760f0788 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 28 Dec 2023 08:14:33 -0500 Subject: [PATCH 17/53] fix save forcing --- src/integrators/integrator_utils.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/integrators/integrator_utils.jl b/src/integrators/integrator_utils.jl index da1ef70034..eab7044f36 100644 --- a/src/integrators/integrator_utils.jl +++ b/src/integrators/integrator_utils.jl @@ -111,7 +111,9 @@ function _savevalues!(integrator, force_save, reduce_size)::Tuple{Bool, Bool} end end if force_save || (integrator.opts.save_everystep && - (isempty(integrator.sol.t) || (integrator.t !== integrator.sol.t[end]))) + (isempty(integrator.sol.t) || (integrator.t !== integrator.sol.t[end]) && + (integrator.opts.save_end || integrator.t !== integrator.sol.prob.tspan[2]) + )) integrator.saveiter += 1 saved, savedexactly = true, true if integrator.opts.save_idxs === nothing From d054a0e3db98e437cedbc8cca8f6a12751c77643 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 28 Dec 2023 09:47:48 -0500 Subject: [PATCH 18/53] Update ode_saveat_tests.jl --- test/interface/ode_saveat_tests.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/interface/ode_saveat_tests.jl b/test/interface/ode_saveat_tests.jl index 52df1a9276..5e1a6de97e 100644 --- a/test/interface/ode_saveat_tests.jl +++ b/test/interface/ode_saveat_tests.jl @@ -160,7 +160,7 @@ integ = init(ODEProblem((u, p, t) -> u, 0.0, (0.0, 1.0)), Tsit5(), saveat = _sav save_end = false) add_tstop!(integ, 2.0) solve!(integ) -@test integ.sol.t == _saveat +@test integ.sol.t == _saveat[1:end-1] # Catch save for maxiters ode = ODEProblem((u, p, t) -> u, 1.0, (0.0, 1.0)) @@ -211,4 +211,4 @@ sol = solve(prob, DP5(), reltol = 1e-6, abstol = 1e-6, saveat = t_obs) @test solve(prob, Tsit5(); saveat = [.2]).t == [0.2] @test solve(prob, Tsit5(); saveat = [.2], save_start = true, save_end = true).t == [0.0; 0.2; 0.4] @test solve(prob, Tsit5(); saveat = [.2], save_start = false, save_end = false).t == [0.2] -end \ No newline at end of file +end From f5e0454a489611dfd3e58aa1872315903defa80c Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Fri, 29 Dec 2023 01:21:59 -0500 Subject: [PATCH 19/53] Fix McAte5 time dependence Fixes https://github.com/SciML/OrdinaryDiffEq.jl/issues/2066 --- src/perform_step/symplectic_perform_step.jl | 4 ++-- test/algconvergence/symplectic_tests.jl | 26 +++++++++++++++++++++ 2 files changed, 28 insertions(+), 2 deletions(-) diff --git a/src/perform_step/symplectic_perform_step.jl b/src/perform_step/symplectic_perform_step.jl index 44a16bb110..3e4f6f849f 100644 --- a/src/perform_step/symplectic_perform_step.jl +++ b/src/perform_step/symplectic_perform_step.jl @@ -568,7 +568,7 @@ end kdu = f.f1(du, u, p, tnew) du = du + dt * a5 * kdu - tnew = tnew + t + a5 * dt + tnew = tnew + a5 * dt ku = f.f2(du, u, p, tnew) u = u + dt * b6 * ku @@ -625,7 +625,7 @@ end f.f1(kdu, du, u, p, tnew) @.. broadcast=false du=du + dt * a5 * kdu - tnew = tnew + t + a5 * dt + tnew = tnew + a5 * dt f.f2(ku, du, u, p, tnew) @.. broadcast=false u=u + dt * b6 * ku diff --git a/test/algconvergence/symplectic_tests.jl b/test/algconvergence/symplectic_tests.jl index 3d496ee4c8..c3cf9145b9 100644 --- a/test/algconvergence/symplectic_tests.jl +++ b/test/algconvergence/symplectic_tests.jl @@ -67,3 +67,29 @@ end @test_throws ArgumentError solve(prob, alg(); dt = dt) end end + +function motionfuncDirect1(dv,v,u,p,t) + # 1:Electron, 2: Be + ω_1,ω_2,γ,m_1,m_2,η,ω_d=p + dv[1]=-ω_1^2*u[1]*(1+η*cos(ω_d*t))-γ*u[2]/m_1 + dv[2]=-ω_2^2*u[2]-γ*u[1]/m_2 +end + +function motionfuncDirect1(v,u,p,t) + # 1:Electron, 2: Be + ω_1,ω_2,γ,m_1,m_2,η,ω_d=p + [-ω_1^2*u[1]*(1+η*cos(ω_d*t))-γ*u[2]/m_1,-ω_2^2*u[2]-γ*u[1]/m_2] +end + +param=[90386.15717208837, 3938.9288690708827, 8560.718748264337, 0.000544617021484666, 8.947079933513658, 0.7596480420227258, 78778.57738141765] +u0_direct=zeros(2) # mm, mm +v0_direct = [0.0, 135.83668926684385] +tspan=(0.0, 1.321179076090661) +prob_direct=SecondOrderODEProblem(motionfuncDirect1,v0_direct,u0_direct,tspan,param) +dt=2e-8 +ref=solve(prob_direct,DPRKN12(),abstol=1e-12,reltol=1e-12,maxiters=1e7,saveat=0.01) + +@testset "symplectic time-dependent $alg-$iip-$pa" for (alg, x, d) in ALGOS + sol=solve(prob_direct,alg,dt=dt,saveat=0.01) + @test maximum(ref-sol) < 1e-3 +end From 4d5c0a753d8f2e5b325bb901b605dea6298fd609 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Fri, 29 Dec 2023 19:29:43 -0500 Subject: [PATCH 20/53] Set dtmin to zero by default MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This has been a long time coming. Ever since the issue https://github.com/SciML/Sundials.jl/pull/423 I've been like "could you actually do that?", since Sundials does seem to default to dtmin=0, so I've been testing it around. We started down this path because the Hairer codes use it, and there's a justification that I gave as to why it might make more sense for Runge-Kutta codes than it does for BDF codes (https://github.com/SciML/Sundials.jl/pull/423). But... I've come to the realization that it is a good idea to have dtmin=0 for a few reasons. It's hard to pinpoint the codes since I've seen a few thousand models since September, but a few things have stood out to me. 1. We've had a lot of false negatives with dtmin. https://github.com/SciML/OrdinaryDiffEq.jl/pull/1830, fixes https://github.com/SciML/DifferentialEquations.jl/issues/924, fixes https://github.com/SciML/OrdinaryDiffEq.jl/issues/1616, fixes https://github.com/SciML/OrdinaryDiffEq.jl/issues/1879. There's just a lot of these cases that can pop up. It's hard to individually root them all out. There was a whole bunch of work that @oscardssmith did last year https://github.com/SciML/OrdinaryDiffEq.jl/pull/1762 in order to try and root out the false positives, and I think there are a lot less, but it led to a separate concern I will mention next. This basically means that any dtmin>0 has a downside at some point getting false negative exits. It's good to exit early in some cases, but in many cases it's "hyperactive". 2. Setting dtmin to a reasonable default to not get false positives ends up being too conservative for it to be useful in a lot of scenarios. To see what I mean, here's some stiff ODE solves: ```julia function rober!(du, u, p, t) y₁, y₂, y₃ = u k₁, k₂, k₃ = p du[1] = -k₁ * y₁ + k₃ * y₂ * y₃ du[2] = k₁ * y₁ - k₂ * y₂^2 - k₃ * y₂ * y₃ du[3] = k₂ * y₂^2 nothing end prob = ODEProblem(rober!, [1.0, 0.0, 0.0], (0.0, 1e5), [0.04, 3e7, 1e4]) sol = solve(prob, Tsit5()) const k1=.35e0 const k2=.266e2 const k3=.123e5 const k4=.86e-3 const k5=.82e-3 const k6=.15e5 const k7=.13e-3 const k8=.24e5 const k9=.165e5 const k10=.9e4 const k11=.22e-1 const k12=.12e5 const k13=.188e1 const k14=.163e5 const k15=.48e7 const k16=.35e-3 const k17=.175e-1 const k18=.1e9 const k19=.444e12 const k20=.124e4 const k21=.21e1 const k22=.578e1 const k23=.474e-1 const k24=.178e4 const k25=.312e1 function f2(dy,y,p,t) r1 = k1 *y[1] r2 = k2 *y[2]*y[4] r3 = k3 *y[5]*y[2] r4 = k4 *y[7] r5 = k5 *y[7] r6 = k6 *y[7]*y[6] r7 = k7 *y[9] r8 = k8 *y[9]*y[6] r9 = k9 *y[11]*y[2] r10 = k10*y[11]*y[1] r11 = k11*y[13] r12 = k12*y[10]*y[2] r13 = k13*y[14] r14 = k14*y[1]*y[6] r15 = k15*y[3] r16 = k16*y[4] r17 = k17*y[4] r18 = k18*y[16] r19 = k19*y[16] r20 = k20*y[17]*y[6] r21 = k21*y[19] r22 = k22*y[19] r23 = k23*y[1]*y[4] r24 = k24*y[19]*y[1] r25 = k25*y[20] dy[1] = -r1-r10-r14-r23-r24+ r2+r3+r9+r11+r12+r22+r25 dy[2] = -r2-r3-r9-r12+r1+r21 dy[3] = -r15+r1+r17+r19+r22 dy[4] = -r2-r16-r17-r23+r15 dy[5] = -r3+r4+r4+r6+r7+r13+r20 dy[6] = -r6-r8-r14-r20+r3+r18+r18 dy[7] = -r4-r5-r6+r13 dy[8] = r4+r5+r6+r7 dy[9] = -r7-r8 dy[10] = -r12+r7+r9 dy[11] = -r9-r10+r8+r11 dy[12] = r9 dy[13] = -r11+r10 dy[14] = -r13+r12 dy[15] = r14 dy[16] = -r18-r19+r16 dy[17] = -r20 dy[18] = r20 dy[19] = -r21-r22-r24+r23+r25 dy[20] = -r25+r24 end function fjac(J,y,p,t) J .= 0.0 J[1,1] = -k1-k10*y[11]-k14*y[6]-k23*y[4]-k24*y[19] J[1,11] = -k10*y[1]+k9*y[2] J[1,6] = -k14*y[1] J[1,4] = -k23*y[1]+k2*y[2] J[1,19] = -k24*y[1]+k22 J[1,2] = k2*y[4]+k9*y[11]+k3*y[5]+k12*y[10] J[1,13] = k11 J[1,20] = k25 J[1,5] = k3*y[2] J[1,10] = k12*y[2] J[2,4] = -k2*y[2] J[2,5] = -k3*y[2] J[2,11] = -k9*y[2] J[2,10] = -k12*y[2] J[2,19] = k21 J[2,1] = k1 J[2,2] = -k2*y[4]-k3*y[5]-k9*y[11]-k12*y[10] J[3,1] = k1 J[3,4] = k17 J[3,16] = k19 J[3,19] = k22 J[3,3] = -k15 J[4,4] = -k2*y[2]-k16-k17-k23*y[1] J[4,2] = -k2*y[4] J[4,1] = -k23*y[4] J[4,3] = k15 J[5,5] = -k3*y[2] J[5,2] = -k3*y[5] J[5,7] = 2k4+k6*y[6] J[5,6] = k6*y[7]+k20*y[17] J[5,9] = k7 J[5,14] = k13 J[5,17] = k20*y[6] J[6,6] = -k6*y[7]-k8*y[9]-k14*y[1]-k20*y[17] J[6,7] = -k6*y[6] J[6,9] = -k8*y[6] J[6,1] = -k14*y[6] J[6,17] = -k20*y[6] J[6,2] = k3*y[5] J[6,5] = k3*y[2] J[6,16] = 2k18 J[7,7] = -k4-k5-k6*y[6] J[7,6] = -k6*y[7] J[7,14] = k13 J[8,7] = k4+k5+k6*y[6] J[8,6] = k6*y[7] J[8,9] = k7 J[9,9] = -k7-k8*y[6] J[9,6] = -k8*y[9] J[10,10] = -k12*y[2] J[10,2] = -k12*y[10]+k9*y[11] J[10,9] = k7 J[10,11] = k9*y[2] J[11,11] = -k9*y[2]-k10*y[1] J[11,2] = -k9*y[11] J[11,1] = -k10*y[11] J[11,9] = k8*y[6] J[11,6] = k8*y[9] J[11,13] = k11 J[12,11] = k9*y[2] J[12,2] = k9*y[11] J[13,13] = -k11 J[13,11] = k10*y[1] J[13,1] = k10*y[11] J[14,14] = -k13 J[14,10] = k12*y[2] J[14,2] = k12*y[10] J[15,1] = k14*y[6] J[15,6] = k14*y[1] J[16,16] = -k18-k19 J[16,4] = k16 J[17,17] = -k20*y[6] J[17,6] = -k20*y[17] J[18,17] = k20*y[6] J[18,6] = k20*y[17] J[19,19] = -k21-k22-k24*y[1] J[19,1] = -k24*y[19]+k23*y[4] J[19,4] = k23*y[1] J[19,20] = k25 J[20,20] = -k25 J[20,1] = k24*y[19] J[20,19] = k24*y[1] return end u0 = zeros(20) u0[2] = 0.2 u0[4] = 0.04 u0[7] = 0.1 u0[8] = 0.3 u0[9] = 0.01 u0[17] = 0.007 prob = ODEProblem(ODEFunction{true, SciMLBase.FullSpecialize}(f2, jac=fjac),u0,(0.0,60.0)) sol = solve(prob,Tsit5()) const N = 32 xyd_brusselator = range(0,stop=1,length=N) brusselator_f(x, y, t) = (((x-0.3)^2 + (y-0.6)^2) <= 0.1^2) * (t >= 1.1) * 5. limit(a, N) = a == N+1 ? 1 : a == 0 ? N : a function brusselator_2d_loop(du, u, p, t) A, B, alpha, dx = p alpha = alpha/dx^2 @inbounds for I in CartesianIndices((N, N)) i, j = Tuple(I) x, y = xyd_brusselator[I[1]], xyd_brusselator[I[2]] ip1, im1, jp1, jm1 = limit(i+1, N), limit(i-1, N), limit(j+1, N), limit(j-1, N) du[i,j,1] = alpha*(u[im1,j,1] + u[ip1,j,1] + u[i,jp1,1] + u[i,jm1,1] - 4u[i,j,1]) + B + u[i,j,1]^2*u[i,j,2] - (A + 1)*u[i,j,1] + brusselator_f(x, y, t) du[i,j,2] = alpha*(u[im1,j,2] + u[ip1,j,2] + u[i,jp1,2] + u[i,jm1,2] - 4u[i,j,2]) + A*u[i,j,1] - u[i,j,1]^2*u[i,j,2] end end p = (3.4, 1., 10., step(xyd_brusselator)) input = rand(N,N,2) output = similar(input) f = ODEFunction{true, SciMLBase.FullSpecialize}(brusselator_2d_loop) function init_brusselator_2d(xyd) N = length(xyd) u = zeros(N, N, 2) for I in CartesianIndices((N, N)) x = xyd[I[1]] y = xyd[I[2]] u[I,1] = 22*(y*(1-y))^(3/2) u[I,2] = 27*(x*(1-x))^(3/2) end u end u0 = init_brusselator_2d(xyd_brusselator) prob = ODEProblem(f,u0,(0.,11.5),p); sol = solve(prob,Tsit5()) ``` The first two stiff ODEs exist with maxiters. The last one actually completes, but used to `dtmin`, so it was a false negative and now just a slow way to compute it but it works. So these are the cases where in theory `dtmin` should cause earlier exits, non-stiff ODE solver on a stiff ODE, and it's not even the thing causing exits here because it's small enough to not be "the" factor. 3. One of the arguments that I have seen for dtmin was that it can reduce the possibility of numerical error. But numerical error with dt being super small isn't compounding, the issue is that you get `x + dt*y = x`. If you look at what happens with very small steps on exponential growth (exponential growth is an ODE which has exponential error growth and is thus a hard not easy test for ODE solver error behavior), then you can see what I mean: ```julia using OrdinaryDiffEq, Plots f = (u,p,t) -> -p.a*u prob = ODEProblem(f, [10.0], (0.0, 0.1), (a=-100.0,)) sol1 = solve(prob, Vern7(), dtmin=0); sol2 = solve(prob, Vern7(), tstops=[sol1.t[2] + eps()], dtmin=0) sol3 = solve(prob, Vern7(), tstops=[sol1.t[2] + 0.01eps()], dtmin=0) sol4 = solve(prob, Vern7(), tstops=[sol1.t[2] + 0.001eps()], dtmin=0) plot(sol1) plot!(sol2) plot!(sol3) plot!(sol4) t = 0.0:1e-5:0.1 @show maximum((sol1(t) - sol2(t)) ./ sol1(t)) # 8.678693835607108e-7 @show maximum((sol1(t) - sol3(t)) ./ sol1(t)) # 2.680284044683277e-13 @show maximum((sol1(t) - sol4(t)) ./ sol1(t)) # 0.0 ``` Small dt's forced by tstops are safe. I've been playing around in many of these examples that have come up, it's safe. The worst that happens is you take a 1e-18 sized step and don't update `u`, but that's because the `u` update is below floating point resolution. You resolve the tstop just fine! So I don't think dtmin is enforcing a safety behavior. And finally... 4. It can do good things to have dtmin=0, and dt=eps() # to prevent infinite loops. - abs(dt) < DiffEqBase.prob2dtmin(prob) && + abs(dt) < dtmin && throw(ArgumentError("Supplied dt is smaller than dtmin")) steps = ceil(Int, internalnorm((tspan[2] - tspan[1]) / dt, tspan[1])) end @@ -358,8 +358,6 @@ function DiffEqBase.__init(prob::Union{DiffEqBase.AbstractODEProblem, controller = default_controller(_alg, cache, qoldinit, beta1, beta2) end - dtmin === nothing && (dtmin = DiffEqBase.prob2dtmin(prob)) - save_end_user = save_end save_end = save_end === nothing ? save_everystep || isempty(saveat) || saveat isa Number || From f6a42e26a8c45686da3acc56df4c63ff5a937d2e Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 30 Dec 2023 05:04:31 -0500 Subject: [PATCH 21/53] test fixes --- test/algconvergence/symplectic_tests.jl | 2 +- test/integrators/check_error.jl | 2 +- test/interface/ode_initdt_tests.jl | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/test/algconvergence/symplectic_tests.jl b/test/algconvergence/symplectic_tests.jl index c3cf9145b9..57de7c0d0e 100644 --- a/test/algconvergence/symplectic_tests.jl +++ b/test/algconvergence/symplectic_tests.jl @@ -89,7 +89,7 @@ prob_direct=SecondOrderODEProblem(motionfuncDirect1,v0_direct,u0_direct,tspan,pa dt=2e-8 ref=solve(prob_direct,DPRKN12(),abstol=1e-12,reltol=1e-12,maxiters=1e7,saveat=0.01) -@testset "symplectic time-dependent $alg-$iip-$pa" for (alg, x, d) in ALGOS +@testset "symplectic time-dependent $alg" for (alg, x, d) in ALGOS sol=solve(prob_direct,alg,dt=dt,saveat=0.01) @test maximum(ref-sol) < 1e-3 end diff --git a/test/integrators/check_error.jl b/test/integrators/check_error.jl index 9287e35dad..50c114e41c 100644 --- a/test/integrators/check_error.jl +++ b/test/integrators/check_error.jl @@ -5,7 +5,7 @@ u0 = 0.0 # explosion time is 1.0 tspan = (0.0, 10.0) prob = ODEProblem(f_ec, u0, tspan) options = [:reltol => 1e-8, :abstol => 1e-8, :verbose => false] -desired_code = ReturnCode.DtLessThanMin +desired_code = ReturnCode.MaxIters # Test that sol.retcode is set to the correct value by various ways to # invoke integrator. diff --git a/test/interface/ode_initdt_tests.jl b/test/interface/ode_initdt_tests.jl index 04f8bd8f04..9e75bbf252 100644 --- a/test/interface/ode_initdt_tests.jl +++ b/test/interface/ode_initdt_tests.jl @@ -28,7 +28,7 @@ prob = remake(prob, u0 = u0, tspan = tspan) tspan = T.((2000, 2100)) prob = remake(prob, tspan = tspan) # set maxiters to prevent infinite loop on test failure -@test_throws ArgumentError solve(prob, Euler(); dt = T(0.0001), maxiters = 10) +@test solve(prob, Euler(); dt = T(0.0001), maxiters = 10).retcode == SciMLBase.ReturnCode.MaxIters function rober(du, u, p, t) y₁, y₂, y₃ = u From 11f33f8e63f259bac631392c370218bae50bd929 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 30 Dec 2023 05:10:46 -0500 Subject: [PATCH 22/53] Fix FunctionMap idxs handling in interpolation Fixes https://github.com/SciML/SciMLBase.jl/issues/575 --- src/dense/interpolants.jl | 16 +++++++++++++-- test/interface/interpolation_output_types.jl | 21 +++++++++++++++++++- 2 files changed, 34 insertions(+), 3 deletions(-) diff --git a/src/dense/interpolants.jl b/src/dense/interpolants.jl index 060423ab81..164e178b24 100644 --- a/src/dense/interpolants.jl +++ b/src/dense/interpolants.jl @@ -30,16 +30,28 @@ end @muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{FunctionMapConstantCache, FunctionMapCache}, - idxs, T::Type{Val{0}}, differential_vars::Nothing) + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) y₀ end -@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, +@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Union{FunctionMapConstantCache, FunctionMapCache}, idxs, T::Type{Val{0}}, differential_vars::Nothing) + y₀[idxs] +end + +@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, + cache::Union{FunctionMapConstantCache, FunctionMapCache}, + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) recursivecopy!(out, y₀) end +@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, + cache::Union{FunctionMapConstantCache, FunctionMapCache}, + idxs, T::Type{Val{0}}, differential_vars::Nothing) + @views out[idxs] .= y₀[idxs] +end + """ Hairer Norsett Wanner Solving Ordinary Differential Euations I - Nonstiff Problems Page 192 """ diff --git a/test/interface/interpolation_output_types.jl b/test/interface/interpolation_output_types.jl index 66f0599baf..f16a0fcdee 100644 --- a/test/interface/interpolation_output_types.jl +++ b/test/interface/interpolation_output_types.jl @@ -1,4 +1,4 @@ -using OrdinaryDiffEq, Test +using OrdinaryDiffEq, RecursiveArrayTools, Test # in terms of the voltage across all three elements rlc1!(v′,v,(R,L,C),t) = -(v′/R + v/L)/C @@ -20,3 +20,22 @@ sol = solve(prob,CalvoSanz4(),dt=1/10) @test sol(0.32, Val{1}) isa OrdinaryDiffEq.ArrayPartition @test sol(0.32, Val{2}) isa OrdinaryDiffEq.ArrayPartition @test sol(0.32, Val{3}) isa OrdinaryDiffEq.ArrayPartition + +function f(du,u,p,t) + du .= u + nothing +end +dprob = DiscreteProblem(f, [1,2,3], (0,100)) +sol = solve(dprob, FunctionMap()) +@test sol(0:0.1:100;idxs=1) isa DiffEqArray +@test length(sol(0:0.1:100;idxs=1)) == length(0:0.1:100) +@test length(sol(0:0.1:100;idxs=1).u[1]) == 1 +sol(0:0.1:100;idxs=[1,2]) + +@test sol(0:0.1:100;idxs=[1,2]) isa DiffEqArray +@test length(sol(0:0.1:100;idxs=[1,2])) == length(0:0.1:100) +@test length(sol(0:0.1:100;idxs=[1,2]).u[1]) == 2 + +@test sol(0:0.1:100) isa DiffEqArray +@test length(sol(0:0.1:100)) == length(0:0.1:100) +@test length(sol(0:0.1:100).u[1]) == 3 From 42ee4c8fc717ca344967a24e386398707cf7adcd Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 30 Dec 2023 05:36:59 -0500 Subject: [PATCH 23/53] typo --- test/algconvergence/symplectic_tests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/algconvergence/symplectic_tests.jl b/test/algconvergence/symplectic_tests.jl index 57de7c0d0e..193f5b46bb 100644 --- a/test/algconvergence/symplectic_tests.jl +++ b/test/algconvergence/symplectic_tests.jl @@ -90,6 +90,6 @@ dt=2e-8 ref=solve(prob_direct,DPRKN12(),abstol=1e-12,reltol=1e-12,maxiters=1e7,saveat=0.01) @testset "symplectic time-dependent $alg" for (alg, x, d) in ALGOS - sol=solve(prob_direct,alg,dt=dt,saveat=0.01) + sol=solve(prob_direct,alg(),dt=dt,saveat=0.01) @test maximum(ref-sol) < 1e-3 end From 0e4888dc99813363b2668aac7c6f6ef4cb46760d Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 30 Dec 2023 06:30:17 -0500 Subject: [PATCH 24/53] tweak some tests --- test/algconvergence/symplectic_tests.jl | 2 +- test/downstream/delaydiffeq.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/test/algconvergence/symplectic_tests.jl b/test/algconvergence/symplectic_tests.jl index 193f5b46bb..cae9004bb5 100644 --- a/test/algconvergence/symplectic_tests.jl +++ b/test/algconvergence/symplectic_tests.jl @@ -91,5 +91,5 @@ ref=solve(prob_direct,DPRKN12(),abstol=1e-12,reltol=1e-12,maxiters=1e7,saveat=0. @testset "symplectic time-dependent $alg" for (alg, x, d) in ALGOS sol=solve(prob_direct,alg(),dt=dt,saveat=0.01) - @test maximum(ref-sol) < 1e-3 + @test maximum(ref[4,:]-sol[4,:]) < 1e-3 end diff --git a/test/downstream/delaydiffeq.jl b/test/downstream/delaydiffeq.jl index 413fe86740..1c6d845c85 100644 --- a/test/downstream/delaydiffeq.jl +++ b/test/downstream/delaydiffeq.jl @@ -24,7 +24,7 @@ using Test @test sol.errors[:l∞] < error sol_scalar = solve(prob_scalar, ddealg) - @test sol.t ≈ sol_scalar.t + @test sol.t ≈ sol_scalar.t atol=1e-6 @test sol[1, :] ≈ sol_scalar.u end end From 8fdf4bfd2885adde5787a6ac43dcf09b4cea6a3b Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 30 Dec 2023 07:11:34 -0500 Subject: [PATCH 25/53] bump tolerance a bit and release --- Project.toml | 2 +- test/algconvergence/symplectic_tests.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index ca3a7c2b4e..45e81e6835 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "OrdinaryDiffEq" uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" authors = ["Chris Rackauckas ", "Yingbo Ma "] -version = "6.65.0" +version = "6.66.0" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" diff --git a/test/algconvergence/symplectic_tests.jl b/test/algconvergence/symplectic_tests.jl index cae9004bb5..894dd4e3dc 100644 --- a/test/algconvergence/symplectic_tests.jl +++ b/test/algconvergence/symplectic_tests.jl @@ -91,5 +91,5 @@ ref=solve(prob_direct,DPRKN12(),abstol=1e-12,reltol=1e-12,maxiters=1e7,saveat=0. @testset "symplectic time-dependent $alg" for (alg, x, d) in ALGOS sol=solve(prob_direct,alg(),dt=dt,saveat=0.01) - @test maximum(ref[4,:]-sol[4,:]) < 1e-3 + @test maximum(ref[4,:]-sol[4,:]) < 3e-3 end From 97a7cb307a9e8e21b67a99b6a31cc567b1e234a1 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 31 Dec 2023 11:16:17 -0500 Subject: [PATCH 26/53] Update NonlinearSolve ForwardDiff support for DAE Initialization diff NonlinearSolve.jl no longer differentiates the solver when used with ForwardDiff and instead solves and applies the chain rule at the solved point. This of course is a good thing, it's much more efficient. However, the DAE initialization algorithms were setup assuming that the element type of the nonlinear solve matches that of the ODE solver, which is violated by this change. This PR thus changes the caches to be non-dual'd to match what NonlinearSolve is expecting, and improves the testing of AD on initialization and solving of DAEs --- Project.toml | 226 ++++++++++++++++----------------- src/initialize_dae.jl | 100 +++++++++++---- test/interface/ad_tests.jl | 18 --- test/interface/dae_ad_tests.jl | 56 ++++++++ test/runtests.jl | 3 +- 5 files changed, 246 insertions(+), 157 deletions(-) create mode 100644 test/interface/dae_ad_tests.jl diff --git a/Project.toml b/Project.toml index 45e81e6835..7c73ec81b2 100644 --- a/Project.toml +++ b/Project.toml @@ -1,113 +1,113 @@ -name = "OrdinaryDiffEq" -uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" -authors = ["Chris Rackauckas ", "Yingbo Ma "] -version = "6.66.0" - -[deps] -ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" -Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" -ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" -DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" -DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" -DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18" -FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898" -FastClosures = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a" -FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" -FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" -ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" -FunctionWrappersWrappers = "77dc65aa-8811-40c2-897b-53d922fa7daf" -IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173" -InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" -LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" -LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" -Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" -LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890" -MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" -MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" -NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" -Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" -PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46" -PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" -Preferences = "21216c6a-2e73-6563-6e65-726566657250" -RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" -Reexport = "189a3867-3050-52da-a836-e630ba90ab69" -SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" -SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961" -SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7" -SimpleUnPack = "ce78b400-467f-4804-87d8-8f486da07d0a" -SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" -SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804" -StaticArrayInterface = "0d7ed370-da01-4f52-bd93-41d350b8b718" -StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" -TruncatedStacktraces = "781d530d-4396-4725-bb49-402e4bee1e77" - -[compat] -ADTypes = "0.2" -Adapt = "3.0, 4" -ArrayInterface = "7" -DataStructures = "0.18" -DiffEqBase = "6.128.2" -DocStringExtensions = "0.9" -ExponentialUtilities = "1.22" -FastBroadcast = "0.2" -FastClosures = "0.3" -FillArrays = "1.9" -FiniteDiff = "2" -ForwardDiff = "0.10.3" -FunctionWrappersWrappers = "0.1" -IfElse = "0.1" -InteractiveUtils = "1.9" -LineSearches = "7" -LinearAlgebra = "1.9" -LinearSolve = "2.1.10" -Logging = "1.9" -LoopVectorization = "0.12" -MacroTools = "0.5" -MuladdMacro = "0.2.1" -NLsolve = "4" -NonlinearSolve = "3" -Polyester = "0.7" -PreallocationTools = "0.4" -PrecompileTools = "1" -Preferences = "1.3" -RecursiveArrayTools = "2.36, 3" -Reexport = "1.0" -SciMLBase = "2" -SciMLOperators = "0.3" -SimpleNonlinearSolve = "1" -SimpleUnPack = "1" -SparseArrays = "1.9" -SparseDiffTools = "2.3" -StaticArrayInterface = "1.2" -StaticArrays = "1.0" -TruncatedStacktraces = "1.2" -julia = "1.9" - -[extras] -AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c" -Calculus = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9" -ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" -DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" -DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d" -ElasticArrays = "fdbdab4c-e67f-52f5-8c3f-e7b388dad3d4" -IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895" -InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" -ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" -NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" -ODEProblemLibrary = "fdc4e326-1af4-4b90-96e7-779fcce2daa5" -Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" -PoissonRandom = "e409e4f3-bfea-5376-8464-e040bb5c01ab" -Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" -Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" -ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" -SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" -SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" -Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" -Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" - -[targets] -test = ["Calculus", "ComponentArrays", "Symbolics", "AlgebraicMultigrid", "IncompleteLU", "DiffEqCallbacks", "DiffEqDevTools", "ODEProblemLibrary", "ElasticArrays", "InteractiveUtils", "PoissonRandom", "Printf", "Random", "ReverseDiff", "SafeTestsets", "SparseArrays", "Statistics", "Test", "Unitful", "ModelingToolkit", "Pkg", "NLsolve"] +name = "OrdinaryDiffEq" +uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +authors = ["Chris Rackauckas ", "Yingbo Ma "] +version = "6.66.0" + +[deps] +ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" +Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" +ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" +DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" +DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" +DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18" +FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898" +FastClosures = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a" +FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" +FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +FunctionWrappersWrappers = "77dc65aa-8811-40c2-897b-53d922fa7daf" +IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173" +InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" +Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" +LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890" +MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" +MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" +NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" +Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" +PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46" +PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" +Preferences = "21216c6a-2e73-6563-6e65-726566657250" +RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" +Reexport = "189a3867-3050-52da-a836-e630ba90ab69" +SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" +SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961" +SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7" +SimpleUnPack = "ce78b400-467f-4804-87d8-8f486da07d0a" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804" +StaticArrayInterface = "0d7ed370-da01-4f52-bd93-41d350b8b718" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +TruncatedStacktraces = "781d530d-4396-4725-bb49-402e4bee1e77" + +[compat] +ADTypes = "0.2" +Adapt = "3.0, 4" +ArrayInterface = "7" +DataStructures = "0.18" +DiffEqBase = "6.128.2" +DocStringExtensions = "0.9" +ExponentialUtilities = "1.22" +FastBroadcast = "0.2" +FastClosures = "0.3" +FillArrays = "1.9" +FiniteDiff = "2" +ForwardDiff = "0.10.3" +FunctionWrappersWrappers = "0.1" +IfElse = "0.1" +InteractiveUtils = "1.9" +LineSearches = "7" +LinearAlgebra = "1.9" +LinearSolve = "2.1.10" +Logging = "1.9" +LoopVectorization = "0.12" +MacroTools = "0.5" +MuladdMacro = "0.2.1" +NLsolve = "4" +NonlinearSolve = "3.3" +Polyester = "0.7" +PreallocationTools = "0.4" +PrecompileTools = "1" +Preferences = "1.3" +RecursiveArrayTools = "2.36, 3" +Reexport = "1.0" +SciMLBase = "2" +SciMLOperators = "0.3" +SimpleNonlinearSolve = "1" +SimpleUnPack = "1" +SparseArrays = "1.9" +SparseDiffTools = "2.3" +StaticArrayInterface = "1.2" +StaticArrays = "1.0" +TruncatedStacktraces = "1.2" +julia = "1.9" + +[extras] +AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c" +Calculus = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9" +ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" +DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" +DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d" +ElasticArrays = "fdbdab4c-e67f-52f5-8c3f-e7b388dad3d4" +IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895" +InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" +NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" +ODEProblemLibrary = "fdc4e326-1af4-4b90-96e7-779fcce2daa5" +Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +PoissonRandom = "e409e4f3-bfea-5376-8464-e040bb5c01ab" +Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" +SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" + +[targets] +test = ["Calculus", "ComponentArrays", "Symbolics", "AlgebraicMultigrid", "IncompleteLU", "DiffEqCallbacks", "DiffEqDevTools", "ODEProblemLibrary", "ElasticArrays", "InteractiveUtils", "PoissonRandom", "Printf", "Random", "ReverseDiff", "SafeTestsets", "SparseArrays", "Statistics", "Test", "Unitful", "ModelingToolkit", "Pkg", "NLsolve"] diff --git a/src/initialize_dae.jl b/src/initialize_dae.jl index bb71431648..72a1de57f5 100644 --- a/src/initialize_dae.jl +++ b/src/initialize_dae.jl @@ -141,6 +141,19 @@ function _initialize_dae!(integrator, prob::ODEProblem, alg::ShampineCollocation failed = nlsolvefail(nlsolver) @.. broadcast=false integrator.u=integrator.uprev + z else + + # _u0 should be non-dual since NonlinearSolve does not differentiate the solver + # These non-dual values are thus used to make the caches + #_du = DiffEqBase.value.(du) + _u0 = DiffEqBase.value.(u0) + + # If not doing auto-diff of the solver, save an allocation + if typeof(u0) === typeof(_u0) + tmp = get_tmp_cache(integrator)[1] + else + tmp = copy(_u0) + end + isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff if isAD chunk = ForwardDiff.pickchunksize(length(tmp)) @@ -150,10 +163,11 @@ function _initialize_dae!(integrator, prob::ODEProblem, alg::ShampineCollocation end nlequation! = @closure (out, u, p) -> begin + T = one(Base.promote_type(eltype(u),eltype(p))) update_coefficients!(M, u, p, t) # f(u,p,t) + M * (u0 - u)/dt - tmp = isAD ? PreallocationTools.get_tmp(_tmp, u) : _tmp - @. tmp = (u0 - u) / dt + tmp = isAD ? PreallocationTools.get_tmp(_tmp, T) : _tmp + @. tmp = (_u0 - u) / dt mul!(_vec(out), M, _vec(tmp)) f(tmp, u, p, t) out .+= tmp @@ -277,7 +291,6 @@ function _initialize_dae!(integrator, prob::DAEProblem, u0 = integrator.u dtmax = integrator.opts.dtmax - tmp = get_tmp_cache(integrator)[1] resid = get_tmp_cache(integrator)[2] dt = t != 0 ? min(t / 1000, dtmax / 10) : dtmax / 10 # Haven't implemented norm reduction @@ -285,6 +298,18 @@ function _initialize_dae!(integrator, prob::DAEProblem, f(resid, integrator.du, u0, p, t) integrator.opts.internalnorm(resid, t) <= integrator.opts.abstol && return + # _du and _u should be non-dual since NonlinearSolve does not differentiate the solver + # These non-dual values are thus used to make the caches + #_du = DiffEqBase.value.(du) + _u0 = DiffEqBase.value.(u0) + + # If not doing auto-diff of the solver, save an allocation + if typeof(u0) === typeof(_u0) + tmp = get_tmp_cache(integrator)[1] + else + tmp = copy(_u0) + end + isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff if isAD chunk = ForwardDiff.pickchunksize(length(tmp)) @@ -294,9 +319,10 @@ function _initialize_dae!(integrator, prob::DAEProblem, end nlequation! = @closure (out, u, p) -> begin - tmp = isAD ? PreallocationTools.get_tmp(_tmp, u) : _tmp + T = one(Base.promote_type(eltype(u),eltype(p))) + tmp = isAD ? PreallocationTools.get_tmp(_tmp, T) : _tmp #M * (u-u0)/dt - f(u,p,t) - @. tmp = (u - u0) / dt + @. tmp = (u - _u0) / dt f(out, tmp, u, p, t) nothing end @@ -414,9 +440,17 @@ function _initialize_dae!(integrator, prob::ODEProblem, integrator.opts.internalnorm(tmp, t) <= alg.abstol && return alg_u = @view u[algebraic_vars] + # These non-dual values are thus used to make the caches + _u = DiffEqBase.value.(u) + + # If auto-diff of the solver, should be non-dual since NonlinearSolve does not differentiate the solver + if typeof(u) !== typeof(_u) + tmp = DiffEqBase.value.(tmp) + end + isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff if isAD - chunk = ForwardDiff.pickchunksize(count(algebraic_vars)) + chunk = ForwardDiff.pickchunksize(max(count(algebraic_vars),length(p))) _tmp = PreallocationTools.dualcache(tmp, chunk) _du_tmp = PreallocationTools.dualcache(similar(tmp), chunk) else @@ -424,9 +458,10 @@ function _initialize_dae!(integrator, prob::ODEProblem, end nlequation! = @closure (out, x, p) -> begin - uu = isAD ? PreallocationTools.get_tmp(_tmp, x) : _tmp - du_tmp = isAD ? PreallocationTools.get_tmp(_du_tmp, x) : _du_tmp - copyto!(uu, integrator.u) + T = one(Base.promote_type(eltype(x),eltype(p))) + uu = isAD ? PreallocationTools.get_tmp(_tmp, T) : _tmp + du_tmp = isAD ? PreallocationTools.get_tmp(_du_tmp, T) : _du_tmp + copyto!(uu, _u) alg_uu = @view uu[algebraic_vars] alg_uu .= x f(du_tmp, uu, p, t) @@ -531,31 +566,46 @@ function _initialize_dae!(integrator, prob::DAEProblem, u = integrator.u du = integrator.du - tmp = get_tmp_cache(integrator)[1] - du_tmp = get_tmp_cache(integrator)[2] - f(tmp, du, u, p, t) + # _du and _u should be non-dual since NonlinearSolve does not differentiate the solver + # These non-dual values are thus used to make the caches + _du = DiffEqBase.value.(du) + _u = DiffEqBase.value.(u) - if integrator.opts.internalnorm(tmp, t) <= alg.abstol + # If not doing auto-diff of the solver, save an allocation + if typeof(u) === typeof(_u) + tmp = get_tmp_cache(integrator)[1] + du_tmp = get_tmp_cache(integrator)[2] + else + tmp = copy(_u) + du_tmp = copy(_du) + end + + # Can be the same as tmp + normtmp = get_tmp_cache(integrator)[1] + f(normtmp, du, u, p, t) + + if integrator.opts.internalnorm(normtmp, t) <= alg.abstol return elseif differential_vars === nothing error("differential_vars must be set for DAE initialization to occur. Either set consistent initial conditions, differential_vars, or use a different initialization algorithm.") end - isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff - if isAD - chunk = ForwardDiff.pickchunksize(length(tmp)) - _tmp = PreallocationTools.dualcache(tmp, chunk) - _du_tmp = PreallocationTools.dualcache(du_tmp, chunk) - else - _tmp, _du_tmp = tmp, similar(tmp) - end + isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff + if isAD + chunk = ForwardDiff.pickchunksize(length(tmp)) + _tmp = PreallocationTools.dualcache(tmp, chunk) + _du_tmp = PreallocationTools.dualcache(du_tmp, chunk) + else + _tmp, _du_tmp = tmp, _du_tmp + end nlequation! = @closure (out, x, p) -> begin - du_tmp = isAD ? PreallocationTools.get_tmp(_du_tmp, x) : _du_tmp - uu = isAD ? PreallocationTools.get_tmp(_tmp, x) : _tmp + T = one(Base.promote_type(eltype(x),eltype(p))) + du_tmp = isAD ? PreallocationTools.get_tmp(_du_tmp, T) : _du_tmp + uu = isAD ? PreallocationTools.get_tmp(_tmp, T) : _tmp - @. du_tmp = ifelse(differential_vars, x, du) - @. uu = ifelse(differential_vars, u, x) + @. du_tmp = ifelse(differential_vars, x, _du) + @. uu = ifelse(differential_vars, _u, x) f(out, du_tmp, uu, p, t) end diff --git a/test/interface/ad_tests.jl b/test/interface/ad_tests.jl index a8cd6ec1e6..1a85c97ad8 100644 --- a/test/interface/ad_tests.jl +++ b/test/interface/ad_tests.jl @@ -271,24 +271,6 @@ end @test ForwardDiff.derivative(g, 0.0)≈_u0 / 2 * exp(-0.5) rtol=rtol end -function f(out, du, u, p, t) - out[1] = -p[1] * u[1] + 1e4 * u[2] * u[3] - du[1] - out[2] = +p[1] * u[1] - 3e7 * u[2]^2 - 1e4 * u[2] * u[3] - du[2] - out[3] = u[1] + u[2] + u[3] - 1.0 -end - -p = [0.5] -u₀ = [1.0, 0, 0] -du₀ = [-0.04, 0.04, 0.0] -tspan = (0.0, 100000.0) -differential_vars = [true, true, false] -prob = DAEProblem(f, du₀, u₀, tspan, p, differential_vars = differential_vars) - -function f(p) - sum(solve(remake(prob, p = p), DABDF2(), saveat = 0.1, abstol = 1e-6, reltol = 1e-6)) -end -@test ForwardDiff.gradient(f, [0.5])[1]≈0 atol=1e-2 - # https://github.com/SciML/DifferentialEquations.jl/issues/903 #ode function diff --git a/test/interface/dae_ad_tests.jl b/test/interface/dae_ad_tests.jl new file mode 100644 index 0000000000..0b7cb7d817 --- /dev/null +++ b/test/interface/dae_ad_tests.jl @@ -0,0 +1,56 @@ +using OrdinaryDiffEq, ForwardDiff, Test + +function rober(du, u, p, t) + y₁, y₂, y₃ = u + k₁, k₂, k₃ = p + du[1] = -k₁ * y₁ + k₃ * y₂ * y₃ + du[2] = k₁ * y₁ - k₃ * y₂ * y₃ - k₂ * y₂^2 + du[3] = y₁ + y₂ + y₃ - 1 + nothing +end +function rober(u, p, t) + y₁, y₂, y₃ = u + k₁, k₂, k₃ = p + [-k₁ * y₁ + k₃ * y₂ * y₃, + k₁ * y₁ - k₃ * y₂ * y₃ - k₂ * y₂^2, + y₁ + y₂ + y₃ - 1] +end +M = [1.0 0 0 + 0 1.0 0 + 0 0 0] +roberf = ODEFunction(rober, mass_matrix = M) +roberf_oop = ODEFunction{false}(rober, mass_matrix = M) +prob_mm = ODEProblem(roberf, [1.0, 0.0, 0.2], (0.0, 1e5), (0.04, 3e7, 1e4)) +prob_mm_oop = ODEProblem(roberf_oop, [1.0, 0.0, 0.2], (0.0, 1e5), (0.04, 3e7, 1e4)) +sol = solve(prob_mm, Rodas5P(), reltol = 1e-8, abstol = 1e-8) +sol = solve(prob_mm_oop, Rodas5P(), reltol = 1e-8, abstol = 1e-8) + +function f(out, du, u, p, t) + out[1] = -p[1] * u[1] + p[3] * u[2] * u[3] - du[1] + out[2] = +p[1] * u[1] - p[2] * u[2]^2 - p[3] * u[2] * u[3] - du[2] + out[3] = u[1] + u[2] + u[3] - 1.0 +end +function f(du, u, p, t) + [-p[1] * u[1] + p[3] * u[2] * u[3] - du[1], + +p[1] * u[1] - p[2] * u[2]^2 - p[3] * u[2] * u[3] - du[2], + u[1] + u[2] + u[3] - 1.0] +end +p = [0.04, 3e7, 1e4] +u₀ = [1.0, 0, 0] +du₀ = [-0.04, 0.04, 0.0] +tspan = (0.0, 100000.0) +differential_vars = [true, true, false] +prob = DAEProblem(f, du₀, u₀, tspan, p, differential_vars = differential_vars) +prob_oop = DAEProblem{false}(f, du₀, u₀, tspan, p, differential_vars = differential_vars) +sol1 = solve(prob, DFBDF(), dt=1e-5, abstol = 1e-8, reltol = 1e-8) + +# These tests flex differentiation of the solver and through the initialization +# To only test the solver part and isolate potential issues, set the initialization to consistent +@testset "Inplace: $(isinplace(_prob)), DAEProblem: $(_prob isa DAEProblem), BrownBasic: $(initalg isa BrownFullBasicInit)" for _prob in [prob, prob_mm, prob_oop, prob_mm_oop], initalg in [BrownFullBasicInit(), ShampineCollocationInit()] + alg = _prob isa DAEProblem ? DFBDF() : Rodas5P() + function f(p) + sol = solve(remake(_prob, p = p), alg, abstol = 1e-14, reltol = 1e-14, initializealg = initalg) + sum(sol) + end + @test ForwardDiff.gradient(f, [0.04, 3e7, 1e4]) ≈ [0,0,0] atol=1e-8 +end diff --git a/test/runtests.jl b/test/runtests.jl index a3b27ccfbf..183e3a4b6f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -84,8 +84,9 @@ end end if !is_APPVEYOR && (GROUP == "All" || GROUP == "InterfaceV" || GROUP == "Interface") - @time @safetestset "AD Tests" include("interface/interpolation_derivative_error_tests.jl") + @time @safetestset "Interpolation Derivative Error Tests" include("interface/interpolation_derivative_error_tests.jl") @time @safetestset "AD Tests" include("interface/ad_tests.jl") + @time @safetestset "DAE AD Tests" include("interface/dae_ad_tests.jl") @time @safetestset "Newton Tests" include("interface/newton_tests.jl") @time @safetestset "DAE Initialize Integration" include("interface/dae_initialize_integration.jl") end From 41808f8ff74d650a19867f1767a2edd9eb238f65 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 31 Dec 2023 12:54:48 -0500 Subject: [PATCH 27/53] bump PAT --- .buildkite/pipeline.yml | 1 - .github/workflows/CI.yml | 1 - Project.toml | 4 ++-- 3 files changed, 2 insertions(+), 4 deletions(-) diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index cb1566aa93..8feaee42f1 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -20,7 +20,6 @@ steps: matrix: setup: version: - - "1.9" - "1" group: - "Regression_I" diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 1524481770..ab02143c1e 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -25,7 +25,6 @@ jobs: - ODEInterfaceRegression - Multithreading version: - - '1.9' - '1' steps: - uses: actions/checkout@v4 diff --git a/Project.toml b/Project.toml index 7c73ec81b2..a83af9c7d7 100644 --- a/Project.toml +++ b/Project.toml @@ -69,7 +69,7 @@ MuladdMacro = "0.2.1" NLsolve = "4" NonlinearSolve = "3.3" Polyester = "0.7" -PreallocationTools = "0.4" +PreallocationTools = "0.4.14" PrecompileTools = "1" Preferences = "1.3" RecursiveArrayTools = "2.36, 3" @@ -83,7 +83,7 @@ SparseDiffTools = "2.3" StaticArrayInterface = "1.2" StaticArrays = "1.0" TruncatedStacktraces = "1.2" -julia = "1.9" +julia = "1.10" [extras] AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c" From dbd807925a56bb7571878ea48ef1a0079cc05688 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 31 Dec 2023 13:25:53 -0500 Subject: [PATCH 28/53] fix non-diff case --- src/initialize_dae.jl | 7 +- test/algconvergence/symplectic_tests.jl | 194 ++++++++++++------------ 2 files changed, 105 insertions(+), 96 deletions(-) diff --git a/src/initialize_dae.jl b/src/initialize_dae.jl index 72a1de57f5..f29e1c187a 100644 --- a/src/initialize_dae.jl +++ b/src/initialize_dae.jl @@ -450,7 +450,12 @@ function _initialize_dae!(integrator, prob::ODEProblem, isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff if isAD - chunk = ForwardDiff.pickchunksize(max(count(algebraic_vars),length(p))) + csize = if p isa SciMLBase.NullParameters || typeof(_u) === typeof(u) + count(algebraic_vars) + else + max(count(algebraic_vars),length(p)) + end + chunk = ForwardDiff.pickchunksize(csize) _tmp = PreallocationTools.dualcache(tmp, chunk) _du_tmp = PreallocationTools.dualcache(similar(tmp), chunk) else diff --git a/test/algconvergence/symplectic_tests.jl b/test/algconvergence/symplectic_tests.jl index 894dd4e3dc..8bae2f9de7 100644 --- a/test/algconvergence/symplectic_tests.jl +++ b/test/algconvergence/symplectic_tests.jl @@ -1,95 +1,99 @@ - -using Test, LinearAlgebra -using OrdinaryDiffEq, DiffEqBase - -# algorithm, dq(p) != p, convergence order -const ALGOS = ((SymplecticEuler, true, 1), - (VelocityVerlet, false, 2), - (VerletLeapfrog, true, 2), - (PseudoVerletLeapfrog, true, 2), - (McAte2, true, 2), - (Ruth3, true, 3), - (McAte3, true, 3), - (CandyRoz4, true, 4), - (McAte4, true, 4), - (CalvoSanz4, true, 4), - (McAte42, true, 1), # known to be broken - (McAte5, true, 5), - (Yoshida6, true, 6), - (KahanLi6, true, 6), - (McAte8, true, 8), - (KahanLi8, true, 8), - (SofSpa10, true, 10)) - -function dp(p, q, pa, t) - 0q .+ pa[2] -end - -function dq(p, q, pa, t) - p .* pa[1] -end - -dp(res, p, q, pa, t) = (res .= dp(p, q, pa, t)) -dq(res, p, q, pa, t) = (res .= dq(p, q, pa, t)) - -dynode(iip, dp, dq) = DynamicalODEFunction{iip}(dp, dq) - -# [0:1] used in dp, dq; [3:4] start values for p0, q0 -const PARAMS = ((1.0, 0.1, 1.0, 0.0), (0.1, 1.0, 1.0, -1.0)) -const IIPS = (true, false) -const TSPAN = (0.0, 1.0) - -solution(t, w) = (w[2] * t + w[3], (w[2] / 2 * t + w[3]) * w[1] * t + w[4]) -apa(iip::Bool, x) = iip ? vcat.(x) : x -errorbound(dt, d, x) = 100 * abs(dt)^d + 1000 * eps(norm(x)) -function printerrors(text, calc, solution, pa, t1) - print(text, ": ") - print(norm(calc[1] - solution(t1, pa)[1]), " ") - print(norm(calc[2] - solution(t1, pa)[2])) - println() -end - -@testset "symplectic $alg-$iip-$pa" for (alg, x, d) in ALGOS, iip in IIPS, pa in PARAMS - dt = 0.01 - tspan = TSPAN - t0, t1 = tspan - dynfun = dynode(iip, dp, dq) - p0, q0 = apa(iip, solution(t0, pa)) - prob = DynamicalODEProblem(dynfun, p0, q0, tspan, pa) - - if x || pa[1] == 1 - sol = solve(prob, alg(); dt = dt) - calc = sol(t1) - # printerrors("$alg-$iip-$pa", calc, solution, pa, t1) - @test calc[1]≈solution(t1, pa)[1] rtol=errorbound(dt, d, calc[1]) - @test calc[2]≈solution(t1, pa)[2] rtol=errorbound(dt, d, calc[2]) - else - @test_throws ArgumentError solve(prob, alg(); dt = dt) - end -end - -function motionfuncDirect1(dv,v,u,p,t) - # 1:Electron, 2: Be - ω_1,ω_2,γ,m_1,m_2,η,ω_d=p - dv[1]=-ω_1^2*u[1]*(1+η*cos(ω_d*t))-γ*u[2]/m_1 - dv[2]=-ω_2^2*u[2]-γ*u[1]/m_2 -end - -function motionfuncDirect1(v,u,p,t) - # 1:Electron, 2: Be - ω_1,ω_2,γ,m_1,m_2,η,ω_d=p - [-ω_1^2*u[1]*(1+η*cos(ω_d*t))-γ*u[2]/m_1,-ω_2^2*u[2]-γ*u[1]/m_2] -end - -param=[90386.15717208837, 3938.9288690708827, 8560.718748264337, 0.000544617021484666, 8.947079933513658, 0.7596480420227258, 78778.57738141765] -u0_direct=zeros(2) # mm, mm -v0_direct = [0.0, 135.83668926684385] -tspan=(0.0, 1.321179076090661) -prob_direct=SecondOrderODEProblem(motionfuncDirect1,v0_direct,u0_direct,tspan,param) -dt=2e-8 -ref=solve(prob_direct,DPRKN12(),abstol=1e-12,reltol=1e-12,maxiters=1e7,saveat=0.01) - -@testset "symplectic time-dependent $alg" for (alg, x, d) in ALGOS - sol=solve(prob_direct,alg(),dt=dt,saveat=0.01) - @test maximum(ref[4,:]-sol[4,:]) < 3e-3 -end + +using Test, LinearAlgebra +using OrdinaryDiffEq, DiffEqBase + +# algorithm, dq(p) != p, convergence order +const ALGOS = ((SymplecticEuler, true, 1), + (VelocityVerlet, false, 2), + (VerletLeapfrog, true, 2), + (PseudoVerletLeapfrog, true, 2), + (McAte2, true, 2), + (Ruth3, true, 3), + (McAte3, true, 3), + (CandyRoz4, true, 4), + (McAte4, true, 4), + (CalvoSanz4, true, 4), + (McAte42, true, 1), # known to be broken + (McAte5, true, 5), + (Yoshida6, true, 6), + (KahanLi6, true, 6), + (McAte8, true, 8), + (KahanLi8, true, 8), + (SofSpa10, true, 10)) + +function dp(p, q, pa, t) + 0q .+ pa[2] +end + +function dq(p, q, pa, t) + p .* pa[1] +end + +dp(res, p, q, pa, t) = (res .= dp(p, q, pa, t)) +dq(res, p, q, pa, t) = (res .= dq(p, q, pa, t)) + +dynode(iip, dp, dq) = DynamicalODEFunction{iip}(dp, dq) + +# [0:1] used in dp, dq; [3:4] start values for p0, q0 +const PARAMS = ((1.0, 0.1, 1.0, 0.0), (0.1, 1.0, 1.0, -1.0)) +const IIPS = (true, false) +const TSPAN = (0.0, 1.0) + +solution(t, w) = (w[2] * t + w[3], (w[2] / 2 * t + w[3]) * w[1] * t + w[4]) +apa(iip::Bool, x) = iip ? vcat.(x) : x +errorbound(dt, d, x) = 100 * abs(dt)^d + 1000 * eps(norm(x)) +function printerrors(text, calc, solution, pa, t1) + print(text, ": ") + print(norm(calc[1] - solution(t1, pa)[1]), " ") + print(norm(calc[2] - solution(t1, pa)[2])) + println() +end + +@testset "symplectic $alg-$iip-$pa" for (alg, x, d) in ALGOS, iip in IIPS, pa in PARAMS + dt = 0.01 + tspan = TSPAN + t0, t1 = tspan + dynfun = dynode(iip, dp, dq) + p0, q0 = apa(iip, solution(t0, pa)) + prob = DynamicalODEProblem(dynfun, p0, q0, tspan, pa) + + if x || pa[1] == 1 + sol = solve(prob, alg(); dt = dt) + calc = sol(t1) + # printerrors("$alg-$iip-$pa", calc, solution, pa, t1) + @test calc[1]≈solution(t1, pa)[1] rtol=errorbound(dt, d, calc[1]) + @test calc[2]≈solution(t1, pa)[2] rtol=errorbound(dt, d, calc[2]) + else + @test_throws ArgumentError solve(prob, alg(); dt = dt) + end +end + +function motionfuncDirect1(dv,v,u,p,t) + # 1:Electron, 2: Be + ω_1,ω_2,γ,m_1,m_2,η,ω_d=p + dv[1]=-ω_1^2*u[1]*(1+η*cos(ω_d*t))-γ*u[2]/m_1 + dv[2]=-ω_2^2*u[2]-γ*u[1]/m_2 +end + +function motionfuncDirect1(v,u,p,t) + # 1:Electron, 2: Be + ω_1,ω_2,γ,m_1,m_2,η,ω_d=p + [-ω_1^2*u[1]*(1+η*cos(ω_d*t))-γ*u[2]/m_1,-ω_2^2*u[2]-γ*u[1]/m_2] +end + +param=[90386.15717208837, 3938.9288690708827, 8560.718748264337, 0.000544617021484666, 8.947079933513658, 0.7596480420227258, 78778.57738141765] +u0_direct=zeros(2) # mm, mm +v0_direct = [0.0, 135.83668926684385] +tspan=(0.0, 1.321179076090661) +prob_direct=SecondOrderODEProblem(motionfuncDirect1,v0_direct,u0_direct,tspan,param) +dt=2e-8 +ref=solve(prob_direct,DPRKN12(),abstol=1e-12,reltol=1e-12,maxiters=1e7,saveat=0.01) + +@testset "symplectic time-dependent $alg" for (alg, x, d) in ALGOS + sol=solve(prob_direct,alg(),dt=dt,saveat=0.01) + if alg <: Yoshida + @test maximum(ref[4,:]-sol[4,:]) < 9e-3 + else + @test maximum(ref[4,:]-sol[4,:]) < 3e-3 + end +end From bbf205d207c4920b2fe88841f02a1ba80597017d Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 31 Dec 2023 13:39:32 -0500 Subject: [PATCH 29/53] just dispatch on the eltype for AD initialize --- src/initialize_dae.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/initialize_dae.jl b/src/initialize_dae.jl index f29e1c187a..6155fb32fe 100644 --- a/src/initialize_dae.jl +++ b/src/initialize_dae.jl @@ -463,7 +463,7 @@ function _initialize_dae!(integrator, prob::ODEProblem, end nlequation! = @closure (out, x, p) -> begin - T = one(Base.promote_type(eltype(x),eltype(p))) + T = Base.promote_type(eltype(x),eltype(p)) uu = isAD ? PreallocationTools.get_tmp(_tmp, T) : _tmp du_tmp = isAD ? PreallocationTools.get_tmp(_du_tmp, T) : _du_tmp copyto!(uu, _u) @@ -605,7 +605,7 @@ function _initialize_dae!(integrator, prob::DAEProblem, end nlequation! = @closure (out, x, p) -> begin - T = one(Base.promote_type(eltype(x),eltype(p))) + T = Base.promote_type(eltype(x),eltype(p)) du_tmp = isAD ? PreallocationTools.get_tmp(_du_tmp, T) : _du_tmp uu = isAD ? PreallocationTools.get_tmp(_tmp, T) : _tmp From 243892a621987aee6a588fd58f060ba2bbd725f6 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 31 Dec 2023 15:41:47 -0500 Subject: [PATCH 30/53] another PreallocationTools bump --- .github/workflows/Downstream.yml | 2 +- Project.toml | 2 +- test/algconvergence/symplectic_tests.jl | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/Downstream.yml b/.github/workflows/Downstream.yml index 61beadf612..ae25ad90da 100644 --- a/.github/workflows/Downstream.yml +++ b/.github/workflows/Downstream.yml @@ -14,7 +14,7 @@ jobs: strategy: fail-fast: false matrix: - julia-version: [1.9, 1] + julia-version: [1] os: [ubuntu-latest] package: - {user: SciML, repo: DelayDiffEq.jl, group: Interface} diff --git a/Project.toml b/Project.toml index a83af9c7d7..04914aa60a 100644 --- a/Project.toml +++ b/Project.toml @@ -69,7 +69,7 @@ MuladdMacro = "0.2.1" NLsolve = "4" NonlinearSolve = "3.3" Polyester = "0.7" -PreallocationTools = "0.4.14" +PreallocationTools = "0.4.15" PrecompileTools = "1" Preferences = "1.3" RecursiveArrayTools = "2.36, 3" diff --git a/test/algconvergence/symplectic_tests.jl b/test/algconvergence/symplectic_tests.jl index 8bae2f9de7..c0442ca591 100644 --- a/test/algconvergence/symplectic_tests.jl +++ b/test/algconvergence/symplectic_tests.jl @@ -91,7 +91,7 @@ ref=solve(prob_direct,DPRKN12(),abstol=1e-12,reltol=1e-12,maxiters=1e7,saveat=0. @testset "symplectic time-dependent $alg" for (alg, x, d) in ALGOS sol=solve(prob_direct,alg(),dt=dt,saveat=0.01) - if alg <: Yoshida + if alg <: Yoshida6 @test maximum(ref[4,:]-sol[4,:]) < 9e-3 else @test maximum(ref[4,:]-sol[4,:]) < 3e-3 From 4b97107691945e919b4bc0c2f3053b04c7445f84 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 31 Dec 2023 17:26:50 -0500 Subject: [PATCH 31/53] Fix eltype computation when parameters is null --- src/initialize_dae.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/initialize_dae.jl b/src/initialize_dae.jl index 6155fb32fe..2d1dbe72a2 100644 --- a/src/initialize_dae.jl +++ b/src/initialize_dae.jl @@ -463,7 +463,7 @@ function _initialize_dae!(integrator, prob::ODEProblem, end nlequation! = @closure (out, x, p) -> begin - T = Base.promote_type(eltype(x),eltype(p)) + T = p isa SciMLBase.NullParameters ? eltype(x) : Base.promote_type(eltype(x),eltype(p)) uu = isAD ? PreallocationTools.get_tmp(_tmp, T) : _tmp du_tmp = isAD ? PreallocationTools.get_tmp(_du_tmp, T) : _du_tmp copyto!(uu, _u) @@ -605,7 +605,7 @@ function _initialize_dae!(integrator, prob::DAEProblem, end nlequation! = @closure (out, x, p) -> begin - T = Base.promote_type(eltype(x),eltype(p)) + T = p isa SciMLBase.NullParameters ? eltype(x) : Base.promote_type(eltype(x),eltype(p)) du_tmp = isAD ? PreallocationTools.get_tmp(_du_tmp, T) : _du_tmp uu = isAD ? PreallocationTools.get_tmp(_tmp, T) : _tmp From 205ae9f7c9322825677547d16630a8f22c42f64a Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 31 Dec 2023 17:55:48 -0500 Subject: [PATCH 32/53] fix other eltype computation --- src/initialize_dae.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/initialize_dae.jl b/src/initialize_dae.jl index 2d1dbe72a2..7c1fbc4da4 100644 --- a/src/initialize_dae.jl +++ b/src/initialize_dae.jl @@ -163,7 +163,7 @@ function _initialize_dae!(integrator, prob::ODEProblem, alg::ShampineCollocation end nlequation! = @closure (out, u, p) -> begin - T = one(Base.promote_type(eltype(u),eltype(p))) + T = p isa SciMLBase.NullParameters ? eltype(x) : Base.promote_type(eltype(x),eltype(p)) update_coefficients!(M, u, p, t) # f(u,p,t) + M * (u0 - u)/dt tmp = isAD ? PreallocationTools.get_tmp(_tmp, T) : _tmp @@ -319,7 +319,7 @@ function _initialize_dae!(integrator, prob::DAEProblem, end nlequation! = @closure (out, u, p) -> begin - T = one(Base.promote_type(eltype(u),eltype(p))) + T = p isa SciMLBase.NullParameters ? eltype(x) : Base.promote_type(eltype(x),eltype(p)) tmp = isAD ? PreallocationTools.get_tmp(_tmp, T) : _tmp #M * (u-u0)/dt - f(u,p,t) @. tmp = (u - _u0) / dt From e6cb2d02ba4c4880ae0d350b469f5ab8e0ab5f08 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 31 Dec 2023 18:18:26 -0500 Subject: [PATCH 33/53] fix typo --- src/initialize_dae.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/initialize_dae.jl b/src/initialize_dae.jl index 7c1fbc4da4..30ef15bdb7 100644 --- a/src/initialize_dae.jl +++ b/src/initialize_dae.jl @@ -163,7 +163,7 @@ function _initialize_dae!(integrator, prob::ODEProblem, alg::ShampineCollocation end nlequation! = @closure (out, u, p) -> begin - T = p isa SciMLBase.NullParameters ? eltype(x) : Base.promote_type(eltype(x),eltype(p)) + T = p isa SciMLBase.NullParameters ? eltype(u) : Base.promote_type(eltype(u),eltype(p)) update_coefficients!(M, u, p, t) # f(u,p,t) + M * (u0 - u)/dt tmp = isAD ? PreallocationTools.get_tmp(_tmp, T) : _tmp @@ -319,7 +319,7 @@ function _initialize_dae!(integrator, prob::DAEProblem, end nlequation! = @closure (out, u, p) -> begin - T = p isa SciMLBase.NullParameters ? eltype(x) : Base.promote_type(eltype(x),eltype(p)) + T = p isa SciMLBase.NullParameters ? eltype(u) : Base.promote_type(eltype(u),eltype(p)) tmp = isAD ? PreallocationTools.get_tmp(_tmp, T) : _tmp #M * (u-u0)/dt - f(u,p,t) @. tmp = (u - _u0) / dt From f2ca03e6a8e9ebcbd8b595a9462a8ef59f95c470 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 1 Jan 2024 01:10:21 -0500 Subject: [PATCH 34/53] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 04914aa60a..328d2dfdd2 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "OrdinaryDiffEq" uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" authors = ["Chris Rackauckas ", "Yingbo Ma "] -version = "6.66.0" +version = "6.67.0" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From ff4e252511c8bb920b4ca3fe587b70e620d22a3a Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 1 Jan 2024 02:26:11 -0500 Subject: [PATCH 35/53] Rodas5(P) 2nd and third derivatives Fixes https://github.com/SciML/OrdinaryDiffEq.jl/issues/2101 Second derivative looks fine, third derivative looks weird. ```julia using OrdinaryDiffEq, ODEProblemLibrary prob = prob_ode_linear sol1 = solve(prob, Rodas5P(), dt = 1 // 2^(8), dense = true) sol2 = solve(prob, Vern9(), dt = 1 // 2^(8), dense = true) using Plots t = 0.0:0.01:1.0 plot(t, Array(sol1(t, Val{1}))) plot!(t, Array(sol2(t, Val{1}))) plot(t, Array(sol1(t, Val{2}))) plot!(t, Array(sol2(t, Val{2}))) plot(t, Array(sol1(t, Val{3}))) plot!(t, Array(sol2(t, Val{3}))) ``` --- src/dense/rosenbrock_interpolants.jl | 68 +++++++++++++++++++++++++++- test/regression/ode_dense_tests.jl | 4 +- 2 files changed, 69 insertions(+), 3 deletions(-) diff --git a/src/dense/rosenbrock_interpolants.jl b/src/dense/rosenbrock_interpolants.jl index 7717f34d6c..866dccd172 100644 --- a/src/dense/rosenbrock_interpolants.jl +++ b/src/dense/rosenbrock_interpolants.jl @@ -311,4 +311,70 @@ end y₀[idxs] + y₁[idxs]) / dt out end -#- + +# Second Derivative +@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Rosenbrock5ConstantCache, + idxs::Nothing, T::Type{Val{2}}, differential_vars) + @inbounds (-2 * k[1] + 2 * k[2] + Θ * (-6 * k[2] + 6 * k[3] - 12 * Θ * k[3])) / dt^2 +end + +@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Rosenbrock5Cache, idxs::Nothing, + T::Type{Val{2}}, differential_vars) + @inbounds @.. broadcast=false (-2 * k[1] + 2 * k[2] + + Θ * (-6 * k[2] + 6 * k[3] - 12 * Θ * k[3]))/dt^2 +end + +@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, + cache::Union{Rosenbrock5ConstantCache, Rosenbrock5Cache}, + idxs, T::Type{Val{2}}, differential_vars) + @.. broadcast=false (-2 * k[1][idxs] + 2 * k[2][idxs] + + Θ * (-6 * k[2][idxs] + 6 * k[3][idxs] - 12 * Θ * k[3][idxs]))/dt^2 +end + +@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, + cache::Union{Rosenbrock5ConstantCache, Rosenbrock5Cache}, + idxs::Nothing, T::Type{Val{2}}, differential_vars) + @.. broadcast=false out= (-2 * k[1] + 2 * k[2] + + Θ * (-6 * k[2] + 6 * k[3] - 12 * Θ * k[3])) / dt^2 + out +end + +@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, + cache::Union{Rosenbrock5ConstantCache, Rosenbrock5Cache}, + idxs, T::Type{Val{2}}, differential_vars) + @views @.. broadcast=false out= (-2 * k[1][idxs] + 2 * k[2][idxs] + + Θ * + (-6 * k[2][idxs] + 6 * k[3][idxs] - 12 * Θ * k[3][idxs])) / dt^2 + out +end + +# Third Derivative +@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Rosenbrock5ConstantCache, + idxs::Nothing, T::Type{Val{3}}, differential_vars) + @inbounds (-6 * k[2] + 6 * k[3] - 24 * Θ * k[3]) / dt^3 +end + +@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, cache::Rosenbrock5Cache, idxs::Nothing, + T::Type{Val{3}}, differential_vars) + @inbounds @.. broadcast=false (-6 * k[2] + 6 * k[3] - 24 * Θ * k[3])/dt^3 +end + +@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, + cache::Union{Rosenbrock5ConstantCache, Rosenbrock5Cache}, + idxs, T::Type{Val{3}}, differential_vars) + @.. broadcast=false (-6 * k[2][idxs] + 6 * k[3][idxs] - 24 * Θ * k[3][idxs])/dt^3 +end + +@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, + cache::Union{Rosenbrock5ConstantCache, Rosenbrock5Cache}, + idxs::Nothing, T::Type{Val{3}}, differential_vars) + @.. broadcast=false out= (-6 * k[2] + 6 * k[3] - 24 * Θ * k[3]) / dt^3 + out +end + +@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, + cache::Union{Rosenbrock5ConstantCache, Rosenbrock5Cache}, + idxs, T::Type{Val{3}}, differential_vars) + @views @.. broadcast=false out= (-6 * k[2][idxs] + 6 * k[3][idxs] - 24 * Θ * k[3][idxs]) / dt^3 + out +end diff --git a/test/regression/ode_dense_tests.jl b/test/regression/ode_dense_tests.jl index dc60c17c0b..5eb1441d6b 100644 --- a/test/regression/ode_dense_tests.jl +++ b/test/regression/ode_dense_tests.jl @@ -537,10 +537,10 @@ regression_test(Rodas4P(), 4e-5, 6e-5, test_diff1 = true, nth_der = 1, dertol = regression_test(Rodas4P2(), 2e-5, 3e-5, test_diff1 = true, nth_der = 1, dertol = 1e-13) # Rodas5 -regression_test(Rodas5(), 2e-6, 3e-6, test_diff1 = true, nth_der = 1, dertol = 1e-13) +regression_test(Rodas5(), 2e-6, 3e-6, test_diff1 = true, nth_der = 3, dertol = 5e-1) # Rodas5P -regression_test(Rodas5P(), 2e-5, 3e-5, test_diff1 = true, nth_der = 1, dertol = 1e-13) +regression_test(Rodas5P(), 2e-5, 3e-5, test_diff1 = true, nth_der = 3, dertol = 5e-1) # ExplicitRK regression_test(ExplicitRK(), 7e-5, 2e-4) From a9a13e0d5a6292ce4f2c5826d90722ee2f43874e Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Wed, 3 Jan 2024 13:10:30 +0530 Subject: [PATCH 36/53] refactor: remove getproperty method for ODEIntegrator This is already handled in SciMLBase --- src/integrators/type.jl | 9 --------- 1 file changed, 9 deletions(-) diff --git a/src/integrators/type.jl b/src/integrators/type.jl index 3229ff80df..1300ac956a 100644 --- a/src/integrators/type.jl +++ b/src/integrators/type.jl @@ -186,12 +186,3 @@ end # When this is changed, DelayDiffEq.jl must be changed as well! TruncatedStacktraces.@truncate_stacktrace ODEIntegrator 2 1 3 4 - -function Base.getproperty(integ::ODEIntegrator, s::Symbol) - if s === :destats - @warn "destats has been deprecated for stats" - getfield(integ, :stats) - else - getfield(integ, s) - end -end From ff72443e91ac7bcc1bb702ce0562f9ea2515cec7 Mon Sep 17 00:00:00 2001 From: Chris Elrod Date: Wed, 3 Jan 2024 08:12:26 -0500 Subject: [PATCH 37/53] No need to depend on LV --- Project.toml | 224 +++++++++++++++++++++--------------------- src/OrdinaryDiffEq.jl | 2 - 2 files changed, 111 insertions(+), 115 deletions(-) diff --git a/Project.toml b/Project.toml index 328d2dfdd2..1ea6225034 100644 --- a/Project.toml +++ b/Project.toml @@ -1,113 +1,111 @@ -name = "OrdinaryDiffEq" -uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" -authors = ["Chris Rackauckas ", "Yingbo Ma "] -version = "6.67.0" - -[deps] -ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" -Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" -ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" -DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" -DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" -DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18" -FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898" -FastClosures = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a" -FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" -FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" -ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" -FunctionWrappersWrappers = "77dc65aa-8811-40c2-897b-53d922fa7daf" -IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173" -InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" -LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" -LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" -Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" -LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890" -MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" -MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" -NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" -Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" -PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46" -PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" -Preferences = "21216c6a-2e73-6563-6e65-726566657250" -RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" -Reexport = "189a3867-3050-52da-a836-e630ba90ab69" -SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" -SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961" -SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7" -SimpleUnPack = "ce78b400-467f-4804-87d8-8f486da07d0a" -SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" -SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804" -StaticArrayInterface = "0d7ed370-da01-4f52-bd93-41d350b8b718" -StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" -TruncatedStacktraces = "781d530d-4396-4725-bb49-402e4bee1e77" - -[compat] -ADTypes = "0.2" -Adapt = "3.0, 4" -ArrayInterface = "7" -DataStructures = "0.18" -DiffEqBase = "6.128.2" -DocStringExtensions = "0.9" -ExponentialUtilities = "1.22" -FastBroadcast = "0.2" -FastClosures = "0.3" -FillArrays = "1.9" -FiniteDiff = "2" -ForwardDiff = "0.10.3" -FunctionWrappersWrappers = "0.1" -IfElse = "0.1" -InteractiveUtils = "1.9" -LineSearches = "7" -LinearAlgebra = "1.9" -LinearSolve = "2.1.10" -Logging = "1.9" -LoopVectorization = "0.12" -MacroTools = "0.5" -MuladdMacro = "0.2.1" -NLsolve = "4" -NonlinearSolve = "3.3" -Polyester = "0.7" -PreallocationTools = "0.4.15" -PrecompileTools = "1" -Preferences = "1.3" -RecursiveArrayTools = "2.36, 3" -Reexport = "1.0" -SciMLBase = "2" -SciMLOperators = "0.3" -SimpleNonlinearSolve = "1" -SimpleUnPack = "1" -SparseArrays = "1.9" -SparseDiffTools = "2.3" -StaticArrayInterface = "1.2" -StaticArrays = "1.0" -TruncatedStacktraces = "1.2" -julia = "1.10" - -[extras] -AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c" -Calculus = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9" -ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" -DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" -DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d" -ElasticArrays = "fdbdab4c-e67f-52f5-8c3f-e7b388dad3d4" -IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895" -InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" -ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" -NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" -ODEProblemLibrary = "fdc4e326-1af4-4b90-96e7-779fcce2daa5" -Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" -PoissonRandom = "e409e4f3-bfea-5376-8464-e040bb5c01ab" -Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" -Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" -ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" -SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" -SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" -Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" -Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" - -[targets] -test = ["Calculus", "ComponentArrays", "Symbolics", "AlgebraicMultigrid", "IncompleteLU", "DiffEqCallbacks", "DiffEqDevTools", "ODEProblemLibrary", "ElasticArrays", "InteractiveUtils", "PoissonRandom", "Printf", "Random", "ReverseDiff", "SafeTestsets", "SparseArrays", "Statistics", "Test", "Unitful", "ModelingToolkit", "Pkg", "NLsolve"] +name = "OrdinaryDiffEq" +uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +authors = ["Chris Rackauckas ", "Yingbo Ma "] +version = "6.67.0" + +[deps] +ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" +Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" +ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" +DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" +DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" +DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18" +FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898" +FastClosures = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a" +FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" +FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +FunctionWrappersWrappers = "77dc65aa-8811-40c2-897b-53d922fa7daf" +IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173" +InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" +Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" +MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" +MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" +NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" +Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" +PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46" +PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" +Preferences = "21216c6a-2e73-6563-6e65-726566657250" +RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" +Reexport = "189a3867-3050-52da-a836-e630ba90ab69" +SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" +SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961" +SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7" +SimpleUnPack = "ce78b400-467f-4804-87d8-8f486da07d0a" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804" +StaticArrayInterface = "0d7ed370-da01-4f52-bd93-41d350b8b718" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +TruncatedStacktraces = "781d530d-4396-4725-bb49-402e4bee1e77" + +[compat] +ADTypes = "0.2" +Adapt = "3.0, 4" +ArrayInterface = "7" +DataStructures = "0.18" +DiffEqBase = "6.128.2" +DocStringExtensions = "0.9" +ExponentialUtilities = "1.22" +FastBroadcast = "0.2" +FastClosures = "0.3" +FillArrays = "1.9" +FiniteDiff = "2" +ForwardDiff = "0.10.3" +FunctionWrappersWrappers = "0.1" +IfElse = "0.1" +InteractiveUtils = "1.9" +LineSearches = "7" +LinearAlgebra = "1.9" +LinearSolve = "2.1.10" +Logging = "1.9" +MacroTools = "0.5" +MuladdMacro = "0.2.1" +NLsolve = "4" +NonlinearSolve = "3.3" +Polyester = "0.7" +PreallocationTools = "0.4.15" +PrecompileTools = "1" +Preferences = "1.3" +RecursiveArrayTools = "2.36, 3" +Reexport = "1.0" +SciMLBase = "2" +SciMLOperators = "0.3" +SimpleNonlinearSolve = "1" +SimpleUnPack = "1" +SparseArrays = "1.9" +SparseDiffTools = "2.3" +StaticArrayInterface = "1.2" +StaticArrays = "1.0" +TruncatedStacktraces = "1.2" +julia = "1.10" + +[extras] +AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c" +Calculus = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9" +ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" +DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" +DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d" +ElasticArrays = "fdbdab4c-e67f-52f5-8c3f-e7b388dad3d4" +IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895" +InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" +NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" +ODEProblemLibrary = "fdc4e326-1af4-4b90-96e7-779fcce2daa5" +Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +PoissonRandom = "e409e4f3-bfea-5376-8464-e040bb5c01ab" +Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" +SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" + +[targets] +test = ["Calculus", "ComponentArrays", "Symbolics", "AlgebraicMultigrid", "IncompleteLU", "DiffEqCallbacks", "DiffEqDevTools", "ODEProblemLibrary", "ElasticArrays", "InteractiveUtils", "PoissonRandom", "Printf", "Random", "ReverseDiff", "SafeTestsets", "SparseArrays", "Statistics", "Test", "Unitful", "ModelingToolkit", "Pkg", "NLsolve"] diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index 2cea945101..91af679eb7 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -18,8 +18,6 @@ using MuladdMacro, SparseArrays, FastClosures using LinearAlgebra -using LoopVectorization - import StaticArrayInterface import InteractiveUtils From 1bb77f54b62eaffaedcba9478fada4f71c7dae93 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Tue, 2 Jan 2024 18:10:05 -0500 Subject: [PATCH 38/53] make dae initialization work for non eltypable p --- src/initialize_dae.jl | 24 +++++++++++++++----- test/integrators/dae_initialization_tests.jl | 9 ++++++++ 2 files changed, 27 insertions(+), 6 deletions(-) diff --git a/src/initialize_dae.jl b/src/initialize_dae.jl index 30ef15bdb7..6f9d920ea3 100644 --- a/src/initialize_dae.jl +++ b/src/initialize_dae.jl @@ -163,7 +163,11 @@ function _initialize_dae!(integrator, prob::ODEProblem, alg::ShampineCollocation end nlequation! = @closure (out, u, p) -> begin - T = p isa SciMLBase.NullParameters ? eltype(u) : Base.promote_type(eltype(u),eltype(p)) + if p isa AbstractArray{<:Dual} + T = Base.promote_type(eltype(u), eltype(p)) + else + T = eltype(u) + end update_coefficients!(M, u, p, t) # f(u,p,t) + M * (u0 - u)/dt tmp = isAD ? PreallocationTools.get_tmp(_tmp, T) : _tmp @@ -319,7 +323,11 @@ function _initialize_dae!(integrator, prob::DAEProblem, end nlequation! = @closure (out, u, p) -> begin - T = p isa SciMLBase.NullParameters ? eltype(u) : Base.promote_type(eltype(u),eltype(p)) + if p isa AbstractArray{<:Dual} + T = Base.promote_type(eltype(u), eltype(p)) + else + T = eltype(u) + end tmp = isAD ? PreallocationTools.get_tmp(_tmp, T) : _tmp #M * (u-u0)/dt - f(u,p,t) @. tmp = (u - _u0) / dt @@ -450,10 +458,10 @@ function _initialize_dae!(integrator, prob::ODEProblem, isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff if isAD - csize = if p isa SciMLBase.NullParameters || typeof(_u) === typeof(u) - count(algebraic_vars) - else + csize = if p isa SciMLBase.NullParameters || typeof(_u) !== typeof(u) max(count(algebraic_vars),length(p)) + else + count(algebraic_vars) end chunk = ForwardDiff.pickchunksize(csize) _tmp = PreallocationTools.dualcache(tmp, chunk) @@ -463,7 +471,11 @@ function _initialize_dae!(integrator, prob::ODEProblem, end nlequation! = @closure (out, x, p) -> begin - T = p isa SciMLBase.NullParameters ? eltype(x) : Base.promote_type(eltype(x),eltype(p)) + if p isa AbstractArray{<:Dual} + T = Base.promote_type(eltype(x), eltype(p)) + else + T = eltype(x) + end uu = isAD ? PreallocationTools.get_tmp(_tmp, T) : _tmp du_tmp = isAD ? PreallocationTools.get_tmp(_du_tmp, T) : _du_tmp copyto!(uu, _u) diff --git a/test/integrators/dae_initialization_tests.jl b/test/integrators/dae_initialization_tests.jl index 77e13227bd..a621510f1a 100644 --- a/test/integrators/dae_initialization_tests.jl +++ b/test/integrators/dae_initialization_tests.jl @@ -197,6 +197,14 @@ integrator = init(prob, DABDF2()) @test integrator.du[1]≈1.0 atol=1e-9 @test integrator.du[2]≈1.0 atol=1e-9 +# test initialization with parameters without eltype/length +struct A +end +probp = DAEProblem(f, du₀, u₀, A(), tspan, differential_vars = differential_vars) +for initializealg in (ShampineCollocationInit(), BrownFullBasicInit()) + @test isapprox(init(probp, DABDF2(); initializealg).u, init(prob, DABDF2(); initializealg).u) + + # to test that we get the right NL solve we need a broken solver. struct BrokenNLSolve <: SciMLBase.AbstractNonlinearAlgorithm BrokenNLSolve(; kwargs...) = new() @@ -216,3 +224,4 @@ prob = ODEProblem(f, ones(3), (0.0, 1.0)) integrator = init(prob, Rodas5P(), initializealg = ShampineCollocationInit(1.0, BrokenNLSolve())) @test all(isequal(reinterpret(Float64, 0xDEADBEEFDEADBEEF)), integrator.u) + From 363bbcb2005629bf39d3df379fdfcf50acfc90e0 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Wed, 3 Jan 2024 11:20:55 -0500 Subject: [PATCH 39/53] typo --- test/integrators/dae_initialization_tests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/integrators/dae_initialization_tests.jl b/test/integrators/dae_initialization_tests.jl index a621510f1a..77a34e4097 100644 --- a/test/integrators/dae_initialization_tests.jl +++ b/test/integrators/dae_initialization_tests.jl @@ -203,7 +203,7 @@ end probp = DAEProblem(f, du₀, u₀, A(), tspan, differential_vars = differential_vars) for initializealg in (ShampineCollocationInit(), BrownFullBasicInit()) @test isapprox(init(probp, DABDF2(); initializealg).u, init(prob, DABDF2(); initializealg).u) - +end # to test that we get the right NL solve we need a broken solver. struct BrokenNLSolve <: SciMLBase.AbstractNonlinearAlgorithm From 97a0606b26947c9dc4833d143baae01e1ca175e1 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Wed, 3 Jan 2024 17:14:56 -0500 Subject: [PATCH 40/53] fixes --- src/initialize_dae.jl | 22 ++++++---- test/integrators/dae_initialization_tests.jl | 44 +++++++++++++++++--- 2 files changed, 53 insertions(+), 13 deletions(-) diff --git a/src/initialize_dae.jl b/src/initialize_dae.jl index 6f9d920ea3..d0256df80a 100644 --- a/src/initialize_dae.jl +++ b/src/initialize_dae.jl @@ -163,7 +163,7 @@ function _initialize_dae!(integrator, prob::ODEProblem, alg::ShampineCollocation end nlequation! = @closure (out, u, p) -> begin - if p isa AbstractArray{<:Dual} + if DiffEqBase.anyeltypedual(p) <: Dual T = Base.promote_type(eltype(u), eltype(p)) else T = eltype(u) @@ -323,7 +323,7 @@ function _initialize_dae!(integrator, prob::DAEProblem, end nlequation! = @closure (out, u, p) -> begin - if p isa AbstractArray{<:Dual} + if DiffEqBase.anyeltypedual(p) <: Dual T = Base.promote_type(eltype(u), eltype(p)) else T = eltype(u) @@ -458,10 +458,12 @@ function _initialize_dae!(integrator, prob::ODEProblem, isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff if isAD - csize = if p isa SciMLBase.NullParameters || typeof(_u) !== typeof(u) - max(count(algebraic_vars),length(p)) - else - count(algebraic_vars) + csize = count(algebraic_vars) + if !(p isa SciMLBase.NullParameters) && typeof(_u) !== typeof(u) + try + csize = max(csize,length(p)) + catch + end end chunk = ForwardDiff.pickchunksize(csize) _tmp = PreallocationTools.dualcache(tmp, chunk) @@ -471,7 +473,7 @@ function _initialize_dae!(integrator, prob::ODEProblem, end nlequation! = @closure (out, x, p) -> begin - if p isa AbstractArray{<:Dual} + if DiffEqBase.anyeltypedual(p) <: Dual T = Base.promote_type(eltype(x), eltype(p)) else T = eltype(x) @@ -617,7 +619,11 @@ function _initialize_dae!(integrator, prob::DAEProblem, end nlequation! = @closure (out, x, p) -> begin - T = p isa SciMLBase.NullParameters ? eltype(x) : Base.promote_type(eltype(x),eltype(p)) + if DiffEqBase.anyeltypedual(p) <: Dual + T = Base.promote_type(eltype(x), eltype(p)) + else + T = eltype(x) + end du_tmp = isAD ? PreallocationTools.get_tmp(_du_tmp, T) : _du_tmp uu = isAD ? PreallocationTools.get_tmp(_tmp, T) : _tmp diff --git a/test/integrators/dae_initialization_tests.jl b/test/integrators/dae_initialization_tests.jl index 77a34e4097..ae81100e0f 100644 --- a/test/integrators/dae_initialization_tests.jl +++ b/test/integrators/dae_initialization_tests.jl @@ -63,6 +63,37 @@ for alg in [Rodas5(autodiff = false), Trapezoid()] @test sum(sol[1]) ≈ 1 end +function rober_no_p(du, u, p, t) + y₁, y₂, y₃ = u + ;(k₁, k₂, k₃) = (0.04, 3e7, 1e4) + du[1] = -k₁ * y₁ + k₃ * y₂ * y₃ + du[2] = k₁ * y₁ - k₃ * y₂ * y₃ - k₂ * y₂^2 + du[3] = y₁ + y₂ + y₃ - 1 + nothing +end + +function rober_oop_no_p(du, u, p, t) + y₁, y₂, y₃ = u + ;(k₁, k₂, k₃) = (0.04, 3e7, 1e4) + du1 = -k₁ * y₁ + k₃ * y₂ * y₃ + du2 = k₁ * y₁ - k₃ * y₂ * y₃ - k₂ * y₂^2 + du3 = y₁ + y₂ + y₃ - 1 + [du1, du2, du3] +end + +# test oop and iip ODE initialization with parameters without eltype/length +struct UnusedParam +end +for f in (ODEFunction(rober_no_p, mass_matrix = M), ODEFunction(rober_oop_no_p, mass_matrix = M)) + local prob, probp + prob = ODEProblem(f, [1.0, 0.0, 1.0], (0.0, 1e5)) + probp = ODEProblem(f, [1.0, 0.0, 1.0], (0.0, 1e5), UnusedParam) + for initializealg in (ShampineCollocationInit(), BrownFullBasicInit()) + isapprox(init(prob, Rodas5(), abstol=1e-10; initializealg).u, + init(prob, Rodas5(), abstol=1e-10; initializealg).u) + end +end + ## DAEProblem f = function (du, u, p, t) @@ -170,6 +201,12 @@ integrator = init(prob, DABDF2(); initializealg = ShampineCollocationInit()) @test_broken integrator.du[2]≈-1.0 atol=1e-9 @test_broken integrator.u[3]≈2.0 atol=1e-9 +# test iip dae initialization with parameters without eltype/length +probp = DAEProblem(f, du₀, u₀, tspan, UnusedParam(), differential_vars = differential_vars) +for initializealg in (ShampineCollocationInit(), BrownFullBasicInit()) + @test isapprox(init(probp, DABDF2(); initializealg).u, init(prob, DABDF2(); initializealg).u) +end + f = function (du, u, p, t) du - u end @@ -196,11 +233,8 @@ integrator = init(prob, DABDF2()) @test integrator.du[1]≈1.0 atol=1e-9 @test integrator.du[2]≈1.0 atol=1e-9 - -# test initialization with parameters without eltype/length -struct A -end -probp = DAEProblem(f, du₀, u₀, A(), tspan, differential_vars = differential_vars) +# test oop DAE initialization with parameters without eltype/length +probp = DAEProblem(f, du₀, u₀, tspan, UnusedParam(), differential_vars = differential_vars) for initializealg in (ShampineCollocationInit(), BrownFullBasicInit()) @test isapprox(init(probp, DABDF2(); initializealg).u, init(prob, DABDF2(); initializealg).u) end From ed3f222379798a6e81486fe800301d6cd8f2d32d Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 4 Jan 2024 00:44:10 -0500 Subject: [PATCH 41/53] Use the computed ForwardDiff type for preallocation choices in init --- src/initialize_dae.jl | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/src/initialize_dae.jl b/src/initialize_dae.jl index d0256df80a..e8574c0a03 100644 --- a/src/initialize_dae.jl +++ b/src/initialize_dae.jl @@ -163,8 +163,9 @@ function _initialize_dae!(integrator, prob::ODEProblem, alg::ShampineCollocation end nlequation! = @closure (out, u, p) -> begin - if DiffEqBase.anyeltypedual(p) <: Dual - T = Base.promote_type(eltype(u), eltype(p)) + TP = DiffEqBase.anyeltypedual(p) + if TP <: Dual + T = Base.promote_type(eltype(u), TP) else T = eltype(u) end @@ -323,8 +324,9 @@ function _initialize_dae!(integrator, prob::DAEProblem, end nlequation! = @closure (out, u, p) -> begin - if DiffEqBase.anyeltypedual(p) <: Dual - T = Base.promote_type(eltype(u), eltype(p)) + TP = DiffEqBase.anyeltypedual(p) + if TP <: Dual + T = Base.promote_type(eltype(u), TP) else T = eltype(u) end @@ -473,8 +475,9 @@ function _initialize_dae!(integrator, prob::ODEProblem, end nlequation! = @closure (out, x, p) -> begin - if DiffEqBase.anyeltypedual(p) <: Dual - T = Base.promote_type(eltype(x), eltype(p)) + TP = DiffEqBase.anyeltypedual(p) + if TP <: Dual + T = Base.promote_type(eltype(x), TP) else T = eltype(x) end @@ -619,8 +622,9 @@ function _initialize_dae!(integrator, prob::DAEProblem, end nlequation! = @closure (out, x, p) -> begin - if DiffEqBase.anyeltypedual(p) <: Dual - T = Base.promote_type(eltype(x), eltype(p)) + TP = DiffEqBase.anyeltypedual(p) + if TP <: Dual + T = Base.promote_type(eltype(x), TP) else T = eltype(x) end From ef267c8e31ef587a156f2ff5ce8a726dacdbc02d Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 4 Jan 2024 01:15:57 -0500 Subject: [PATCH 42/53] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 1ea6225034..edcff5d0bd 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "OrdinaryDiffEq" uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" authors = ["Chris Rackauckas ", "Yingbo Ma "] -version = "6.67.0" +version = "6.68.0" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From cc31776421a585eec95d95bd57f2f1893044e193 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 4 Jan 2024 01:16:29 -0500 Subject: [PATCH 43/53] Test CUDA 5 --- test/gpu/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/gpu/Project.toml b/test/gpu/Project.toml index 5a49c5b551..39a8540125 100644 --- a/test/gpu/Project.toml +++ b/test/gpu/Project.toml @@ -3,4 +3,4 @@ CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" [compat] -CUDA = "4.2" +CUDA = "5" From e07c3778f4ff38b799bf63f373669381f1b9dc15 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Fri, 5 Jan 2024 19:13:19 -0500 Subject: [PATCH 44/53] Fix DAE initialization AD when autodiff=false Simply needs to still create the dual cache --- src/initialize_dae.jl | 11 +++++++---- test/interface/dae_ad_tests.jl | 6 +++--- 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/src/initialize_dae.jl b/src/initialize_dae.jl index e8574c0a03..34945a60a4 100644 --- a/src/initialize_dae.jl +++ b/src/initialize_dae.jl @@ -154,7 +154,7 @@ function _initialize_dae!(integrator, prob::ODEProblem, alg::ShampineCollocation tmp = copy(_u0) end - isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff + isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff || typeof(u0) !== typeof(_u0) if isAD chunk = ForwardDiff.pickchunksize(length(tmp)) _tmp = PreallocationTools.dualcache(tmp, chunk) @@ -315,7 +315,7 @@ function _initialize_dae!(integrator, prob::DAEProblem, tmp = copy(_u0) end - isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff + isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff || typeof(u0) !== typeof(_u0) if isAD chunk = ForwardDiff.pickchunksize(length(tmp)) _tmp = PreallocationTools.dualcache(tmp, chunk) @@ -458,7 +458,7 @@ function _initialize_dae!(integrator, prob::ODEProblem, tmp = DiffEqBase.value.(tmp) end - isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff + isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff || typeof(u) !== typeof(_u) if isAD csize = count(algebraic_vars) if !(p isa SciMLBase.NullParameters) && typeof(_u) !== typeof(u) @@ -483,6 +483,9 @@ function _initialize_dae!(integrator, prob::ODEProblem, end uu = isAD ? PreallocationTools.get_tmp(_tmp, T) : _tmp du_tmp = isAD ? PreallocationTools.get_tmp(_du_tmp, T) : _du_tmp + @show isAD + @show T, TP + @show typeof(x), typeof(p), typeof(du_tmp) copyto!(uu, _u) alg_uu = @view uu[algebraic_vars] alg_uu .= x @@ -612,7 +615,7 @@ function _initialize_dae!(integrator, prob::DAEProblem, error("differential_vars must be set for DAE initialization to occur. Either set consistent initial conditions, differential_vars, or use a different initialization algorithm.") end - isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff + isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff || typeof(u) !== typeof(_u) if isAD chunk = ForwardDiff.pickchunksize(length(tmp)) _tmp = PreallocationTools.dualcache(tmp, chunk) diff --git a/test/interface/dae_ad_tests.jl b/test/interface/dae_ad_tests.jl index 0b7cb7d817..d8976af291 100644 --- a/test/interface/dae_ad_tests.jl +++ b/test/interface/dae_ad_tests.jl @@ -1,4 +1,4 @@ -using OrdinaryDiffEq, ForwardDiff, Test +using OrdinaryDiffEq, LinearAlgebra, ForwardDiff, Test function rober(du, u, p, t) y₁, y₂, y₃ = u @@ -46,8 +46,8 @@ sol1 = solve(prob, DFBDF(), dt=1e-5, abstol = 1e-8, reltol = 1e-8) # These tests flex differentiation of the solver and through the initialization # To only test the solver part and isolate potential issues, set the initialization to consistent -@testset "Inplace: $(isinplace(_prob)), DAEProblem: $(_prob isa DAEProblem), BrownBasic: $(initalg isa BrownFullBasicInit)" for _prob in [prob, prob_mm, prob_oop, prob_mm_oop], initalg in [BrownFullBasicInit(), ShampineCollocationInit()] - alg = _prob isa DAEProblem ? DFBDF() : Rodas5P() +@testset "Inplace: $(isinplace(_prob)), DAEProblem: $(_prob isa DAEProblem), BrownBasic: $(initalg isa BrownFullBasicInit), Autodiff: $autodiff" for _prob in [prob, prob_mm, prob_oop, prob_mm_oop], initalg in [BrownFullBasicInit(), ShampineCollocationInit()], autodiff in [true, false] + alg = _prob isa DAEProblem ? DFBDF(;autodiff) : Rodas5P(;autodiff) function f(p) sol = solve(remake(_prob, p = p), alg, abstol = 1e-14, reltol = 1e-14, initializealg = initalg) sum(sol) From ab89c3e566633bbd750badddcdd678263feff3ea Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Fri, 5 Jan 2024 19:49:21 -0500 Subject: [PATCH 45/53] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index edcff5d0bd..1633b2187c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "OrdinaryDiffEq" uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" authors = ["Chris Rackauckas ", "Yingbo Ma "] -version = "6.68.0" +version = "6.68.1" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From 50af2be647beeef9b803fc06ae81eb686828d0c5 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 6 Jan 2024 11:06:05 -0500 Subject: [PATCH 46/53] hotfix show --- src/initialize_dae.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/initialize_dae.jl b/src/initialize_dae.jl index 34945a60a4..9d7537168b 100644 --- a/src/initialize_dae.jl +++ b/src/initialize_dae.jl @@ -483,9 +483,6 @@ function _initialize_dae!(integrator, prob::ODEProblem, end uu = isAD ? PreallocationTools.get_tmp(_tmp, T) : _tmp du_tmp = isAD ? PreallocationTools.get_tmp(_du_tmp, T) : _du_tmp - @show isAD - @show T, TP - @show typeof(x), typeof(p), typeof(du_tmp) copyto!(uu, _u) alg_uu = @view uu[algebraic_vars] alg_uu .= x From a62c4d11ffa9b5b1979d85338d98796cf2ca4e79 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 6 Jan 2024 11:06:35 -0500 Subject: [PATCH 47/53] patch --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 1633b2187c..c71c063238 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "OrdinaryDiffEq" uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" authors = ["Chris Rackauckas ", "Yingbo Ma "] -version = "6.68.1" +version = "6.68.2" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From b4bdf24333f948d6d47e95983f35eed6a6bd1168 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 7 Jan 2024 13:55:56 -0500 Subject: [PATCH 48/53] Single type for sensitivity interpolation and composites --- Project.toml | 2 +- src/dense/generic_dense.jl | 4 ++++ src/interp_func.jl | 40 ++++++++++++-------------------------- src/solve.jl | 19 ++++++------------ 4 files changed, 23 insertions(+), 42 deletions(-) diff --git a/Project.toml b/Project.toml index c71c063238..9e0c1dd9be 100644 --- a/Project.toml +++ b/Project.toml @@ -72,7 +72,7 @@ PrecompileTools = "1" Preferences = "1.3" RecursiveArrayTools = "2.36, 3" Reexport = "1.0" -SciMLBase = "2" +SciMLBase = "2.17" SciMLOperators = "0.3" SimpleNonlinearSolve = "1" SimpleUnPack = "1" diff --git a/src/dense/generic_dense.jl b/src/dense/generic_dense.jl index a943579eb2..3683c25bed 100644 --- a/src/dense/generic_dense.jl +++ b/src/dense/generic_dense.jl @@ -360,6 +360,7 @@ function ode_interpolation(tvals, id::I, idxs, deriv::D, p, i₋ = max(1, _searchsortedlast(ts, t, i₋, tdir > 0)) i₊ = i₋ < lastindex(ts) ? i₋ + 1 : i₋ end + id.sensitivitymode && error(SENSITIVITY_INTERP_MESSAGE) i₋₊ref[] = (i₋, i₊) dt = ts[i₊] - ts[i₋] Θ = iszero(dt) ? oneunit(t) / oneunit(dt) : (t - ts[i₋]) / dt @@ -405,6 +406,7 @@ function ode_interpolation!(vals, tvals, id::I, idxs, deriv::D, p, i₋ = max(1, _searchsortedlast(ts, t, i₋, tdir > 0)) i₊ = i₋ < lastindex(ts) ? i₋ + 1 : i₋ end + id.sensitivitymode && error(SENSITIVITY_INTERP_MESSAGE) dt = ts[i₊] - ts[i₋] Θ = iszero(dt) ? oneunit(t) / oneunit(dt) : (t - ts[i₋]) / dt @@ -477,6 +479,7 @@ function ode_interpolation(tval::Number, id::I, idxs, deriv::D, p, i₋ = max(1, _searchsortedlast(ts, tval, 1, tdir > 0)) i₊ = i₋ < lastindex(ts) ? i₋ + 1 : i₋ end + id.sensitivitymode && error(SENSITIVITY_INTERP_MESSAGE) @inbounds begin dt = ts[i₊] - ts[i₋] @@ -525,6 +528,7 @@ function ode_interpolation!(out, tval::Number, id::I, idxs, deriv::D, p, i₋ = max(1, _searchsortedlast(ts, tval, 1, tdir > 0)) i₊ = i₋ < lastindex(ts) ? i₋ + 1 : i₋ end + id.sensitivitymode && error(SENSITIVITY_INTERP_MESSAGE) @inbounds begin dt = ts[i₊] - ts[i₋] diff --git a/src/interp_func.jl b/src/interp_func.jl index d6be19dddb..ea302b7bc1 100644 --- a/src/interp_func.jl +++ b/src/interp_func.jl @@ -1,27 +1,25 @@ abstract type OrdinaryDiffEqInterpolation{cacheType} <: DiffEqBase.AbstractDiffEqInterpolation end -struct InterpolationData{F, uType, tType, kType, cacheType, DV} <: +struct InterpolationData{F, uType, tType, kType, algType <: Union{Nothing, Vector{Int}}, cacheType, DV} <: OrdinaryDiffEqInterpolation{cacheType} f::F timeseries::uType ts::tType ks::kType + alg_choice::algType dense::Bool cache::cacheType differential_vars::DV + sensitivitymode::Bool end -struct CompositeInterpolationData{F, uType, tType, kType, cacheType, DV} <: - OrdinaryDiffEqInterpolation{cacheType} - f::F - timeseries::uType - ts::tType - ks::kType - alg_choice::Vector{Int} - dense::Bool - cache::cacheType - differential_vars::DV +@static if isdefined(SciMLBase, :enable_interpolation_sensitivitymode) + function SciMLBase.enable_interpolation_sensitivitymode(interp::InterpolationData) + InterpolationData(interp.f,interp.timeseries,interp.ts,interp.ks, + interp.alg_choice, interp.dense, interp.cache, + interp.differential_vars, true) + end end function DiffEqBase.interp_summary(interp::OrdinaryDiffEqInterpolation{ @@ -169,31 +167,17 @@ end function (interp::InterpolationData)(tvals, idxs, deriv, p, continuity::Symbol = :left) ode_interpolation(tvals, interp, idxs, deriv, p, continuity) end -function (interp::CompositeInterpolationData)(tvals, idxs, deriv, p, - continuity::Symbol = :left) - ode_interpolation(tvals, interp, idxs, deriv, p, continuity) -end function (interp::InterpolationData)(val, tvals, idxs, deriv, p, continuity::Symbol = :left) ode_interpolation!(val, tvals, interp, idxs, deriv, p, continuity) end -function (interp::CompositeInterpolationData)(val, tvals, idxs, deriv, p, - continuity::Symbol = :left) - ode_interpolation!(val, tvals, interp, idxs, deriv, p, continuity) -end function InterpolationData(id::InterpolationData, f) InterpolationData(f, id.timeseries, - id.ts, - id.ks, - id.dense, - id.cache) -end - -function CompositeInterpolationData(id::CompositeInterpolationData, f) - CompositeInterpolationData(f, id.timeseries, id.ts, id.ks, id.alg_choice, id.dense, - id.cache) + id.cache, + id.differential_vars, + id.sensitivitymode) end diff --git a/src/solve.jl b/src/solve.jl index c0a01ad126..b04e3b665f 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -410,19 +410,12 @@ function DiffEqBase.__init(prob::Union{DiffEqBase.AbstractODEProblem, stats = SciMLBase.DEStats(0) differential_vars = prob isa DAEProblem ? prob.differential_vars : get_differential_vars(f, u) - - if _alg isa OrdinaryDiffEqCompositeAlgorithm - id = CompositeInterpolationData(f, timeseries, ts, ks, alg_choice, dense, cache, differential_vars) - sol = DiffEqBase.build_solution(prob, _alg, ts, timeseries, - dense = dense, k = ks, interp = id, - alg_choice = alg_choice, - calculate_error = false, stats = stats) - else - id = InterpolationData(f, timeseries, ts, ks, dense, cache, differential_vars) - sol = DiffEqBase.build_solution(prob, _alg, ts, timeseries, - dense = dense, k = ks, interp = id, - calculate_error = false, stats = stats) - end + + id = InterpolationData(f, timeseries, ts, ks, alg_choice, dense, cache, differential_vars) + sol = DiffEqBase.build_solution(prob, _alg, ts, timeseries, + dense = dense, k = ks, interp = id, + alg_choice = alg_choice, + calculate_error = false, stats = stats) if recompile_flag == true FType = typeof(f) From 0891c277955682ab3a8201659a2468b4186235d2 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 7 Jan 2024 14:10:17 -0500 Subject: [PATCH 49/53] fix constructor --- src/solve.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solve.jl b/src/solve.jl index b04e3b665f..1962b02327 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -411,7 +411,7 @@ function DiffEqBase.__init(prob::Union{DiffEqBase.AbstractODEProblem, stats = SciMLBase.DEStats(0) differential_vars = prob isa DAEProblem ? prob.differential_vars : get_differential_vars(f, u) - id = InterpolationData(f, timeseries, ts, ks, alg_choice, dense, cache, differential_vars) + id = InterpolationData(f, timeseries, ts, ks, alg_choice, dense, cache, differential_vars, false) sol = DiffEqBase.build_solution(prob, _alg, ts, timeseries, dense = dense, k = ks, interp = id, alg_choice = alg_choice, From 24d5738fbabe4f4ab25913de52ee79b5ce2c811e Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 7 Jan 2024 14:34:58 -0500 Subject: [PATCH 50/53] fix type definition --- src/interp_func.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/interp_func.jl b/src/interp_func.jl index ea302b7bc1..07613e2f26 100644 --- a/src/interp_func.jl +++ b/src/interp_func.jl @@ -1,7 +1,7 @@ abstract type OrdinaryDiffEqInterpolation{cacheType} <: DiffEqBase.AbstractDiffEqInterpolation end -struct InterpolationData{F, uType, tType, kType, algType <: Union{Nothing, Vector{Int}}, cacheType, DV} <: +struct InterpolationData{F, uType, tType, kType, algType <: Union{Tuple{}, Vector{Int}}, cacheType, DV} <: OrdinaryDiffEqInterpolation{cacheType} f::F timeseries::uType From 03d7857a9cd320ce21f63ed84f8c5c873b990efa Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 7 Jan 2024 14:54:40 -0500 Subject: [PATCH 51/53] No alg_choice is nothing --- src/interp_func.jl | 2 +- src/solve.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/interp_func.jl b/src/interp_func.jl index 07613e2f26..ea302b7bc1 100644 --- a/src/interp_func.jl +++ b/src/interp_func.jl @@ -1,7 +1,7 @@ abstract type OrdinaryDiffEqInterpolation{cacheType} <: DiffEqBase.AbstractDiffEqInterpolation end -struct InterpolationData{F, uType, tType, kType, algType <: Union{Tuple{}, Vector{Int}}, cacheType, DV} <: +struct InterpolationData{F, uType, tType, kType, algType <: Union{Nothing, Vector{Int}}, cacheType, DV} <: OrdinaryDiffEqInterpolation{cacheType} f::F timeseries::uType diff --git a/src/solve.jl b/src/solve.jl index 1962b02327..6b8e2ddd34 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -272,7 +272,7 @@ function DiffEqBase.__init(prob::Union{DiffEqBase.AbstractODEProblem, ts = ts_init === () ? tType[] : convert(Vector{tType}, ts_init) ks = ks_init === () ? ksEltype[] : convert(Vector{ksEltype}, ks_init) - alg_choice = _alg isa CompositeAlgorithm ? Int[] : () + alg_choice = _alg isa CompositeAlgorithm ? Int[] : nothing if (!adaptive || !isadaptive(_alg)) && save_everystep && tspan[2] - tspan[1] != Inf if dt == 0 From ffcad28475f24a8082df8a46d78c5f10ac302109 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 7 Jan 2024 15:14:48 -0500 Subject: [PATCH 52/53] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 9e0c1dd9be..029cfde0cf 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "OrdinaryDiffEq" uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" authors = ["Chris Rackauckas ", "Yingbo Ma "] -version = "6.68.2" +version = "6.69.0" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From eab9e653d4ed56c2d5b52dfdec1a2b4c7e805341 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Wed, 17 Jan 2024 18:35:01 +0200 Subject: [PATCH 53/53] use ForwardDiff.value when computing internal tolerances --- src/solve.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/solve.jl b/src/solve.jl index 6b8e2ddd34..c999c84f66 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -195,10 +195,10 @@ function DiffEqBase.__init(prob::Union{DiffEqBase.AbstractODEProblem, reltol_internal = false elseif reltol === nothing if uBottomEltypeNoUnits == uBottomEltype - reltol_internal = real(convert(uBottomEltype, - oneunit(uBottomEltype) * 1 // 10^3)) + reltol_internal = ForwardDiff.value(real(convert(uBottomEltype, + oneunit(uBottomEltype) * 1 // 10^3))) else - reltol_internal = real.(oneunit.(u) .* 1 // 10^3) + reltol_internal = ForwardDiff.value.(real.(oneunit.(u) .* 1 // 10^3)) end else reltol_internal = real.(reltol)