diff --git a/lib/OrdinaryDiffEqCore/Project.toml b/lib/OrdinaryDiffEqCore/Project.toml index 5f29c546a8..aa198fcdbc 100644 --- a/lib/OrdinaryDiffEqCore/Project.toml +++ b/lib/OrdinaryDiffEqCore/Project.toml @@ -14,6 +14,7 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56" FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898" FastClosures = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a" +FastPower = "a4df4552-cc26-4903-aec0-212e50a0e84b" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" FunctionWrappersWrappers = "77dc65aa-8811-40c2-897b-53d922fa7daf" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" @@ -53,6 +54,7 @@ DocStringExtensions = "0.9" EnumX = "1" FastBroadcast = "0.2, 0.3" FastClosures = "0.3" +FastPower = "1" FillArrays = "1.9" FunctionWrappersWrappers = "0.1" InteractiveUtils = "1.9" diff --git a/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl b/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl index 9276a467f9..6376270182 100644 --- a/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl +++ b/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl @@ -19,6 +19,8 @@ using PrecompileTools import FillArrays: Trues, Falses +import FastPower + # Interfaces import DiffEqBase: solve!, step!, initialize!, isadaptive diff --git a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl index a13c3daf88..1f47f61ebc 100644 --- a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl +++ b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl @@ -68,7 +68,7 @@ end q = inv(qmax) else expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1) - qtmp = DiffEqBase.fastpow(EEst, expo) / gamma + qtmp = FastPower.fastpower(EEst, expo) / gamma @fastmath q = DiffEqBase.value(max(inv(qmax), min(inv(qmin), qtmp))) # TODO: Shouldn't this be in `step_accept_controller!` as for the PI controller? integrator.qold = DiffEqBase.value(integrator.dt) / q @@ -141,8 +141,8 @@ end if iszero(EEst) q = inv(qmax) else - q11 = DiffEqBase.fastpow(EEst, float(beta1)) - q = q11 / DiffEqBase.fastpow(qold, float(beta2)) + q11 = FastPower.fastpower(EEst, float(beta1)) + q = q11 / FastPower.fastpower(qold, float(beta2)) integrator.q11 = q11 @fastmath q = max(inv(qmax), min(inv(qmin), q / gamma)) end diff --git a/lib/OrdinaryDiffEqExtrapolation/Project.toml b/lib/OrdinaryDiffEqExtrapolation/Project.toml index 8889d5355e..5e6abc5c66 100644 --- a/lib/OrdinaryDiffEqExtrapolation/Project.toml +++ b/lib/OrdinaryDiffEqExtrapolation/Project.toml @@ -6,6 +6,7 @@ version = "1.1.0" [deps] DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898" +FastPower = "a4df4552-cc26-4903-aec0-212e50a0e84b" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8" @@ -18,6 +19,7 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69" DiffEqBase = "6.152.2" DiffEqDevTools = "2.44.4" FastBroadcast = "0.3.5" +FastPower = "1" LinearSolve = "2.32.0" MuladdMacro = "0.2.4" OrdinaryDiffEqCore = "1.1" diff --git a/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl b/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl index 29bbc6ca8b..3482022699 100644 --- a/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl +++ b/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl @@ -20,6 +20,7 @@ import OrdinaryDiffEqCore: alg_order, alg_maximum_order, get_current_adaptive_or differentiation_rk_docstring using DiffEqBase, FastBroadcast, Polyester, MuladdMacro, RecursiveArrayTools, LinearSolve import OrdinaryDiffEqCore +import FastPower import OrdinaryDiffEqDifferentiation: TimeDerivativeWrapper, UDerivativeWrapper, calc_J, WOperator, TimeGradientWrapper, UJacobianWrapper, build_grad_config, diff --git a/lib/OrdinaryDiffEqExtrapolation/src/controllers.jl b/lib/OrdinaryDiffEqExtrapolation/src/controllers.jl index 84a280f565..a6ae7ea315 100644 --- a/lib/OrdinaryDiffEqExtrapolation/src/controllers.jl +++ b/lib/OrdinaryDiffEqExtrapolation/src/controllers.jl @@ -31,10 +31,10 @@ function stepsize_controller_internal!(integrator, else # Update gamma and beta1 controller.beta1 = typeof(controller.beta1)(1 // (2integrator.cache.n_curr + 1)) - integrator.opts.gamma = DiffEqBase.fastpow(typeof(integrator.opts.gamma)(1 // 4), + integrator.opts.gamma = FastPower.fastpower(typeof(integrator.opts.gamma)(1 // 4), controller.beta1) # Compute new stepsize scaling - qtmp = DiffEqBase.fastpow(integrator.EEst, controller.beta1) / integrator.opts.gamma + qtmp = FastPower.fastpower(integrator.EEst, controller.beta1) / integrator.opts.gamma @fastmath q = max(inv(integrator.opts.qmax), min(inv(integrator.opts.qmin), qtmp)) end integrator.cache.Q[integrator.cache.n_curr - alg.min_order + 1] = q @@ -57,10 +57,10 @@ function stepsize_predictor!(integrator, s_new = stage_number[n_new - alg.min_order + 1] # Update gamma and beta1 controller.beta1 = typeof(controller.beta1)(1 // (2integrator.cache.n_curr + 1)) - integrator.opts.gamma = DiffEqBase.fastpow(typeof(integrator.opts.gamma)(1 // 4), + integrator.opts.gamma = FastPower.fastpower(typeof(integrator.opts.gamma)(1 // 4), controller.beta1) # Compute new stepsize scaling - qtmp = EEst * DiffEqBase.fastpow(DiffEqBase.fastpow(tol, (1.0 - s_curr / s_new)), + qtmp = EEst * FastPower.fastpower(FastPower.fastpower(tol, (1.0 - s_curr / s_new)), controller.beta1) / integrator.opts.gamma @fastmath q = max(inv(integrator.opts.qmax), min(inv(integrator.opts.qmin), qtmp)) end @@ -158,12 +158,12 @@ function stepsize_controller_internal!(integrator, controller.beta1 = typeof(controller.beta1)(1 // (integrator.cache.n_curr - 1)) end - integrator.opts.gamma = DiffEqBase.fastpow( + integrator.opts.gamma = FastPower.fastpower( typeof(integrator.opts.gamma)(65 // 100), controller.beta1) # Compute new stepsize scaling - qtmp = DiffEqBase.fastpow(integrator.EEst, controller.beta1) / + qtmp = FastPower.fastpower(integrator.EEst, controller.beta1) / (integrator.opts.gamma) @fastmath q = max(inv(integrator.opts.qmax), min(inv(integrator.opts.qmin), qtmp)) @@ -175,12 +175,12 @@ function stepsize_controller_internal!(integrator, else # Update gamma and beta1 controller.beta1 = typeof(controller.beta1)(1 // (2integrator.cache.n_curr + 1)) - integrator.opts.gamma = DiffEqBase.fastpow( + integrator.opts.gamma = FastPower.fastpower( typeof(integrator.opts.gamma)(65 // 100), controller.beta1) # Compute new stepsize scaling - qtmp = DiffEqBase.fastpow(integrator.EEst, controller.beta1) / + qtmp = FastPower.fastpower(integrator.EEst, controller.beta1) / integrator.opts.gamma @fastmath q = max(inv(integrator.opts.qmax), min(inv(integrator.opts.qmin), qtmp)) diff --git a/lib/OrdinaryDiffEqFIRK/Project.toml b/lib/OrdinaryDiffEqFIRK/Project.toml index 2ffd2e2710..a3b3fb2a4c 100644 --- a/lib/OrdinaryDiffEqFIRK/Project.toml +++ b/lib/OrdinaryDiffEqFIRK/Project.toml @@ -6,6 +6,7 @@ version = "1.1.1" [deps] DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898" +FastPower = "a4df4552-cc26-4903-aec0-212e50a0e84b" GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a" GenericSchur = "c145ed77-6b09-5dd9-b285-bf645a82121e" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -25,6 +26,7 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" DiffEqBase = "6.152.2" DiffEqDevTools = "2.44.4" FastBroadcast = "0.3.5" +FastPower = "1" LinearAlgebra = "<0.0.1, 1" LinearSolve = "2.32.0" MuladdMacro = "0.2.4" diff --git a/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl b/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl index 40305282d7..753f094704 100644 --- a/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl +++ b/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl @@ -24,6 +24,7 @@ using LinearAlgebra: I, UniformScaling, mul!, lu import LinearSolve import FastBroadcast: @.. import OrdinaryDiffEqCore +import FastPower using OrdinaryDiffEqDifferentiation: UJacobianWrapper, build_J_W, build_jac_config, UDerivativeWrapper, calc_J!, dolinsolve, calc_J, diff --git a/lib/OrdinaryDiffEqFIRK/src/controllers.jl b/lib/OrdinaryDiffEqFIRK/src/controllers.jl index df6cf153ab..2887bacd4f 100644 --- a/lib/OrdinaryDiffEqFIRK/src/controllers.jl +++ b/lib/OrdinaryDiffEqFIRK/src/controllers.jl @@ -17,7 +17,7 @@ fac = min(gamma, (1 + 2 * maxiters) * gamma / (iter + 2 * maxiters)) end expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1) - qtmp = DiffEqBase.fastpow(EEst, expo) / fac + qtmp = FastPower.fastpower(EEst, expo) / fac @fastmath q = DiffEqBase.value(max(inv(qmax), min(inv(qmin), qtmp))) integrator.qold = q end @@ -31,7 +31,7 @@ function step_accept_controller!(integrator, controller::PredictiveController, a if integrator.success_iter > 0 expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1) qgus = (integrator.dtacc / integrator.dt) * - DiffEqBase.fastpow((EEst^2) / integrator.erracc, expo) + FastPower.fastpower((EEst^2) / integrator.erracc, expo) qgus = max(inv(qmax), min(inv(qmin), qgus / gamma)) qacc = max(q, qgus) else diff --git a/test/interface/composite_algorithm_test.jl b/test/interface/composite_algorithm_test.jl index be3a496fcf..01bedd5953 100644 --- a/test/interface/composite_algorithm_test.jl +++ b/test/interface/composite_algorithm_test.jl @@ -47,11 +47,11 @@ v = @inferred OrdinaryDiffEqCore.ode_extrapolant( alg_mixed_r = CompositeAlgorithm((ABM54(), Tsit5()), reverse_choice) alg_mixed2 = CompositeAlgorithm((Tsit5(), ABM54()), reverse_choice) - @test_throws ErrorException solve(prob_ode_linear, alg_mixed) + @test_throws ArgumentError solve(prob_ode_linear, alg_mixed) sol2 = solve(prob_ode_linear, Tsit5()) - sol3 = solve(prob_ode_linear, alg_mixed; dt = 0.05) - sol4 = solve(prob_ode_linear, alg_mixed_r; dt = 0.05) - sol5 = solve(prob_ode_linear, alg_mixed2; dt = 0.05) + sol3 = solve(prob_ode_linear, alg_mixed; dt = 0.05, adaptive=false) + sol4 = solve(prob_ode_linear, alg_mixed_r; dt = 0.05, adaptive=false) + sol5 = solve(prob_ode_linear, alg_mixed2; dt = 0.05, adaptive=false) @test sol3.t == sol4.t && sol3.u == sol4.u @test sol3(0.8)≈sol2(0.8) atol=1e-4 @test sol5(0.8)≈sol2(0.8) atol=1e-4