From f445f25d000b4dc2f3a4f2ee766029a9a3a21585 Mon Sep 17 00:00:00 2001 From: CompatHelper Julia Date: Sun, 3 Nov 2024 20:15:29 +0000 Subject: [PATCH 01/12] CompatHelper: bump compat for SimpleNonlinearSolve to 2, (keep existing compat) --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 69e498bc76..2dcb48233d 100644 --- a/Project.toml +++ b/Project.toml @@ -136,7 +136,7 @@ Reexport = "1.0" SciMLBase = "2.53" SciMLOperators = "0.3" SciMLStructures = "1" -SimpleNonlinearSolve = "1" +SimpleNonlinearSolve = "1, 2" SimpleUnPack = "1" SparseDiffTools = "2" Static = "0.8, 1" From 21a5a66e24f8dd6e8b0cdc2fe89617fa4794b2fd Mon Sep 17 00:00:00 2001 From: CompatHelper Julia Date: Sun, 3 Nov 2024 22:15:12 +0000 Subject: [PATCH 02/12] CompatHelper: bump compat for NonlinearSolve to 4, (keep existing compat) --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 69e498bc76..59e4c810a7 100644 --- a/Project.toml +++ b/Project.toml @@ -96,7 +96,7 @@ LinearSolve = "2" Logging = "1.9" MacroTools = "0.5" MuladdMacro = "0.2.1" -NonlinearSolve = "3" +NonlinearSolve = "3, 4" OrdinaryDiffEqAdamsBashforthMoulton = "1" OrdinaryDiffEqBDF = "1" OrdinaryDiffEqCore = "1" From e50bf9a4d2baf1111dc5d7dd2b3be4d9fba46d1a Mon Sep 17 00:00:00 2001 From: Shreyas Ekanathan Date: Wed, 13 Nov 2024 17:19:20 -0500 Subject: [PATCH 03/12] rename things and fix broadcast --- lib/OrdinaryDiffEqFIRK/src/algorithms.jl | 4 ++-- lib/OrdinaryDiffEqFIRK/src/controllers.jl | 6 +++--- lib/OrdinaryDiffEqFIRK/src/firk_caches.jl | 19 ++++++++++++------- .../src/firk_perform_step.jl | 13 +++++++++---- 4 files changed, 26 insertions(+), 16 deletions(-) diff --git a/lib/OrdinaryDiffEqFIRK/src/algorithms.jl b/lib/OrdinaryDiffEqFIRK/src/algorithms.jl index 75ccae073c..ba45841ac6 100644 --- a/lib/OrdinaryDiffEqFIRK/src/algorithms.jl +++ b/lib/OrdinaryDiffEqFIRK/src/algorithms.jl @@ -169,7 +169,7 @@ end function AdaptiveRadau(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(), concrete_jac = nothing, - diff_type = Val{:forward}, min_stages = 3, max_stages = 7, + diff_type = Val{:forward}, min_order = 5, max_order = 13, linsolve = nothing, precs = DEFAULT_PRECS, extrapolant = :dense, fast_convergence_cutoff = 1 // 5, new_W_γdt_cutoff = 1 // 5, @@ -187,6 +187,6 @@ function AdaptiveRadau(; chunk_size = Val{0}(), autodiff = Val{true}(), fast_convergence_cutoff, new_W_γdt_cutoff, controller, - step_limiter!, min_stages, max_stages) + step_limiter!, min_order, max_order) end diff --git a/lib/OrdinaryDiffEqFIRK/src/controllers.jl b/lib/OrdinaryDiffEqFIRK/src/controllers.jl index 6de40914d0..95816c8b32 100644 --- a/lib/OrdinaryDiffEqFIRK/src/controllers.jl +++ b/lib/OrdinaryDiffEqFIRK/src/controllers.jl @@ -23,11 +23,11 @@ function step_accept_controller!(integrator, controller::PredictiveController, a hist_iter = hist_iter * 0.8 + iter * 0.2 cache.hist_iter = hist_iter if (step > 10) - if (hist_iter < 2.6 && num_stages < alg.max_stages) + if (hist_iter < 2.6 && num_stages < (alg.max_order + 1) ÷ 2) cache.num_stages += 2 cache.step = 1 cache.hist_iter = iter - elseif ((hist_iter > 8 || cache.status == VerySlowConvergence || cache.status == Divergence) && num_stages > alg.min_stages) + elseif ((hist_iter > 8 || cache.status == VerySlowConvergence || cache.status == Divergence) && num_stages > (alg.min_order + 1) ÷ 2) cache.num_stages -= 2 cache.step = 1 cache.hist_iter = iter @@ -45,7 +45,7 @@ function step_reject_controller!(integrator, controller::PredictiveController, a hist_iter = hist_iter * 0.8 + iter * 0.2 cache.hist_iter = hist_iter if (step > 10) - if ((hist_iter > 8 || cache.status == VerySlowConvergence || cache.status == Divergence) && num_stages > alg.min_stages) + if ((hist_iter > 8 || cache.status == VerySlowConvergence || cache.status == Divergence) && num_stages > (alg.min_order + 1) ÷ 2) cache.num_stages -= 2 cache.step = 1 cache.hist_iter = iter diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl index 8dbcd690c7..e1dbe08ad8 100644 --- a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl +++ b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl @@ -497,12 +497,12 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits} ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} uf = UDerivativeWrapper(f, t, p) uToltype = constvalue(uBottomEltypeNoUnits) - num_stages = alg.min_stages - max = alg.max_stages + max = (alg.max_order + 1) ÷ 2 + num_stages = (alg.min_order + 1) ÷ 2 tabs = [BigRadauIIA5Tableau(uToltype, constvalue(tTypeNoUnits)), BigRadauIIA9Tableau(uToltype, constvalue(tTypeNoUnits)), BigRadauIIA13Tableau(uToltype, constvalue(tTypeNoUnits))] i = 9 - while i <= alg.max_stages + while i <= max push!(tabs, adaptiveRadauTableau(uToltype, constvalue(tTypeNoUnits), i)) i += 2 end @@ -525,6 +525,8 @@ mutable struct AdaptiveRadauCache{uType, cuType, tType, uNoUnitsType, rateType, z::Vector{uType} w::Vector{uType} c_prime::Vector{tType} + αdt::Vector{tType} + βdt::Vector{tType} dw1::uType ubuff::uType dw2::Vector{cuType} @@ -568,8 +570,8 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits} uf = UJacobianWrapper(f, t, p) uToltype = constvalue(uBottomEltypeNoUnits) - max = alg.max_stages - num_stages = alg.min_stages + max = (alg.max_order + 1) ÷ 2 + num_stages = (alg.min_order + 1) ÷ 2 tabs = [BigRadauIIA5Tableau(uToltype, constvalue(tTypeNoUnits)), BigRadauIIA9Tableau(uToltype, constvalue(tTypeNoUnits)), BigRadauIIA13Tableau(uToltype, constvalue(tTypeNoUnits))] i = 9 @@ -583,9 +585,12 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits} z = Vector{typeof(u)}(undef, max) w = Vector{typeof(u)}(undef, max) for i in 1 : max - z[i] = w[i] = zero(u) + z[i] = zero(u) + w[i] = zero(u) end + αdt = [zero(t) for i in 1:max] + βdt = [zero(t) for i in 1:max] c_prime = Vector{typeof(t)}(undef, max) #time stepping for i in 1 : max c_prime[i] = zero(t) @@ -641,7 +646,7 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits} atol = reltol isa Number ? reltol : zero(reltol) AdaptiveRadauCache(u, uprev, - z, w, c_prime, dw1, ubuff, dw2, cubuff, dw, cont, derivatives, + z, w, c_prime, αdt, βdt, dw1, ubuff, dw2, cubuff, dw, cont, derivatives, du1, fsalfirst, ks, k, fw, J, W1, W2, uf, tabs, κ, one(uToltype), 10000, tmp, diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl b/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl index 5e42aa5b94..874d968cb6 100644 --- a/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl +++ b/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl @@ -1598,7 +1598,7 @@ end @unpack num_stages, tabs = cache tab = tabs[(num_stages - 1) ÷ 2] @unpack T, TI, γ, α, β, c, e = tab - @unpack κ, cont, derivatives, z, w, c_prime = cache + @unpack κ, cont, derivatives, z, w, c_prime, αdt, βdt= cache @unpack dw1, ubuff, dw2, cubuff, dw = cache @unpack ks, k, fw, J, W1, W2 = cache @unpack tmp, atmp, jac_config, linsolve1, linsolve2, rtol, atol, step_limiter! = cache @@ -1608,7 +1608,12 @@ end mass_matrix = integrator.f.mass_matrix # precalculations - γdt, αdt, βdt = γ / dt, α ./ dt, β ./ dt + γdt = γ / dt + for i in 1 : (num_stages - 1) ÷ 2 + αdt[i] = α[i]/dt + βdt[i] = β[i]/dt + end + (new_jac = do_newJ(integrator, alg, cache, repeat_step)) && (calc_J!(J, integrator, cache); cache.W_γdt = dt) if (new_W = do_newW(integrator, alg, new_jac, cache.W_γdt)) @@ -1750,12 +1755,12 @@ end # transform `w` to `z` #mul!(z, T, w) for i in 1:num_stages - 1 - z[i] = zero(u) + @.. z[i] = zero(u) for j in 1:num_stages @.. z[i] += T[i,j] * w[j] end end - z[num_stages] = T[num_stages, 1] * w[1] + @.. z[num_stages] = T[num_stages, 1] * w[1] i = 2 while i < num_stages @.. z[num_stages] += w[i] From 937641f9cbf10234152a12300f82c1338ecbe1c1 Mon Sep 17 00:00:00 2001 From: Shreyas Ekanathan Date: Mon, 18 Nov 2024 21:31:51 -0500 Subject: [PATCH 04/12] speed ups --- lib/OrdinaryDiffEqFIRK/src/algorithms.jl | 4 +- lib/OrdinaryDiffEqFIRK/src/controllers.jl | 9 ++- lib/OrdinaryDiffEqFIRK/src/firk_caches.jl | 39 ++++++++-- .../src/firk_perform_step.jl | 75 +++++++++---------- 4 files changed, 75 insertions(+), 52 deletions(-) diff --git a/lib/OrdinaryDiffEqFIRK/src/algorithms.jl b/lib/OrdinaryDiffEqFIRK/src/algorithms.jl index ba45841ac6..e98389ec7c 100644 --- a/lib/OrdinaryDiffEqFIRK/src/algorithms.jl +++ b/lib/OrdinaryDiffEqFIRK/src/algorithms.jl @@ -163,8 +163,8 @@ struct AdaptiveRadau{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <: new_W_γdt_cutoff::C2 controller::Symbol step_limiter!::StepLimiter - min_stages::Int - max_stages::Int + min_order::Int + max_order::Int end function AdaptiveRadau(; chunk_size = Val{0}(), autodiff = Val{true}(), diff --git a/lib/OrdinaryDiffEqFIRK/src/controllers.jl b/lib/OrdinaryDiffEqFIRK/src/controllers.jl index 95816c8b32..3849d8114c 100644 --- a/lib/OrdinaryDiffEqFIRK/src/controllers.jl +++ b/lib/OrdinaryDiffEqFIRK/src/controllers.jl @@ -22,12 +22,14 @@ function step_accept_controller!(integrator, controller::PredictiveController, a cache.step = step + 1 hist_iter = hist_iter * 0.8 + iter * 0.2 cache.hist_iter = hist_iter + max_stages = (alg.max_order - 1) ÷ 4 * 2 + 1 + min_stages = (alg.min_order - 1) ÷ 4 * 2 + 1 if (step > 10) - if (hist_iter < 2.6 && num_stages < (alg.max_order + 1) ÷ 2) + if (hist_iter < 2.6 && num_stages <= max_stages) cache.num_stages += 2 cache.step = 1 cache.hist_iter = iter - elseif ((hist_iter > 8 || cache.status == VerySlowConvergence || cache.status == Divergence) && num_stages > (alg.min_order + 1) ÷ 2) + elseif ((hist_iter > 8 || cache.status == VerySlowConvergence || cache.status == Divergence) && num_stages >= min_stages) cache.num_stages -= 2 cache.step = 1 cache.hist_iter = iter @@ -44,8 +46,9 @@ function step_reject_controller!(integrator, controller::PredictiveController, a cache.step = step + 1 hist_iter = hist_iter * 0.8 + iter * 0.2 cache.hist_iter = hist_iter + min_stages = (alg.min_order - 1) ÷ 4 * 2 + 1 if (step > 10) - if ((hist_iter > 8 || cache.status == VerySlowConvergence || cache.status == Divergence) && num_stages > (alg.min_order + 1) ÷ 2) + if ((hist_iter > 8 || cache.status == VerySlowConvergence || cache.status == Divergence) && num_stages >= min_stages) cache.num_stages -= 2 cache.step = 1 cache.hist_iter = iter diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl index e1dbe08ad8..b971b3cc4f 100644 --- a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl +++ b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl @@ -362,6 +362,10 @@ mutable struct RadauIIA9Cache{uType, cuType, uNoUnitsType, rateType, JType, W1Ty tmp4::uType tmp5::uType tmp6::uType + tmp7::uType + tmp8::uType + tmp9::uType + tmp10::uType atmp::uNoUnitsType jac_config::JC linsolve1::F1 @@ -440,6 +444,10 @@ function alg_cache(alg::RadauIIA9, u, rate_prototype, ::Type{uEltypeNoUnits}, tmp4 = zero(u) tmp5 = zero(u) tmp6 = zero(u) + tmp7 = zero(u) + tmp8 = zero(u) + tmp9 = zero(u) + tmp10 = zero(u) atmp = similar(u, uEltypeNoUnits) recursivefill!(atmp, false) jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, dw1) @@ -469,7 +477,7 @@ function alg_cache(alg::RadauIIA9, u, rate_prototype, ::Type{uEltypeNoUnits}, du1, fsalfirst, k, k2, k3, k4, k5, fw1, fw2, fw3, fw4, fw5, J, W1, W2, W3, uf, tab, κ, one(uToltype), 10000, - tmp, tmp2, tmp3, tmp4, tmp5, tmp6, atmp, jac_config, + tmp, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9, tmp10, atmp, jac_config, linsolve1, linsolve2, linsolve3, rtol, atol, dt, dt, Convergence, alg.step_limiter!) end @@ -497,17 +505,26 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits} ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} uf = UDerivativeWrapper(f, t, p) uToltype = constvalue(uBottomEltypeNoUnits) - max = (alg.max_order + 1) ÷ 2 - num_stages = (alg.min_order + 1) ÷ 2 + + max_order = alg.max_order + min_order = alg.min_order + max = (max_order - 1) ÷ 4 * 2 + 1 + min = (min_order - 1) ÷ 4 * 2 + 1 + if (alg.min_order < 5) + error("min_order choice $min_order below 5 is not compatible with the algorithm") + elseif (max < min) + error("max_order $max_order is below min_order $min_order") + end + num_stages = min + tabs = [BigRadauIIA5Tableau(uToltype, constvalue(tTypeNoUnits)), BigRadauIIA9Tableau(uToltype, constvalue(tTypeNoUnits)), BigRadauIIA13Tableau(uToltype, constvalue(tTypeNoUnits))] - i = 9 while i <= max push!(tabs, adaptiveRadauTableau(uToltype, constvalue(tTypeNoUnits), i)) i += 2 end cont = Vector{typeof(u)}(undef, max) - for i in 1: max + for i in 1:max cont[i] = zero(u) end @@ -570,8 +587,16 @@ function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits} uf = UJacobianWrapper(f, t, p) uToltype = constvalue(uBottomEltypeNoUnits) - max = (alg.max_order + 1) ÷ 2 - num_stages = (alg.min_order + 1) ÷ 2 + max_order = alg.max_order + min_order = alg.min_order + max = (max_order - 1) ÷ 4 * 2 + 1 + min = (min_order - 1) ÷ 4 * 2 + 1 + if (alg.min_order < 5) + error("min_order choice $min_order below 5 is not compatible with the algorithm") + elseif (max < min) + error("max_order $max_order is below min_order $min_order") + end + num_stages = min tabs = [BigRadauIIA5Tableau(uToltype, constvalue(tTypeNoUnits)), BigRadauIIA9Tableau(uToltype, constvalue(tTypeNoUnits)), BigRadauIIA13Tableau(uToltype, constvalue(tTypeNoUnits))] i = 9 diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl b/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl index 874d968cb6..2058c4fb15 100644 --- a/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl +++ b/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl @@ -1032,7 +1032,7 @@ end @unpack dw1, ubuff, dw23, dw45, cubuff1, cubuff2 = cache @unpack k, k2, k3, k4, k5, fw1, fw2, fw3, fw4, fw5 = cache @unpack J, W1, W2, W3 = cache - @unpack tmp, tmp2, tmp3, tmp4, tmp5, tmp6, atmp, jac_config, linsolve1, linsolve2, linsolve3, rtol, atol, step_limiter! = cache + @unpack tmp, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9, tmp10, atmp, jac_config, linsolve1, linsolve2, linsolve3, rtol, atol, step_limiter! = cache @unpack internalnorm, abstol, reltol, adaptive = integrator.opts alg = unwrap_alg(integrator, true) @unpack maxiters = alg @@ -1087,30 +1087,30 @@ end c2′ = c2 * c5′ c3′ = c3 * c5′ c4′ = c4 * c5′ - z1 = @.. c1′ * (cont1 + + @.. z1 = c1′ * (cont1 + (c1′-c4m1) * (cont2 + (c1′ - c3m1) * (cont3 + (c1′ - c2m1) * (cont4 + (c1′ - c1m1) * cont5)))) - z2 = @.. c2′ * (cont1 + + @.. z2 = c2′ * (cont1 + (c2′-c4m1) * (cont2 + (c2′ - c3m1) * (cont3 + (c2′ - c2m1) * (cont4 + (c2′ - c1m1) * cont5)))) - z3 = @.. c3′ * (cont1 + + @.. z3 = c3′ * (cont1 + (c3′-c4m1) * (cont2 + (c3′ - c3m1) * (cont3 + (c3′ - c2m1) * (cont4 + (c3′ - c1m1) * cont5)))) - z4 = @.. c4′ * (cont1 + + @.. z4 = c4′ * (cont1 + (c4′-c4m1) * (cont2 + (c4′ - c3m1) * (cont3 + (c4′ - c2m1) * (cont4 + (c4′ - c1m1) * cont5)))) - z5 = @.. c5′ * (cont1 + + @.. z5 = c5′ * (cont1 + (c5′-c4m1) * (cont2 + (c5′ - c3m1) * (cont3 + (c5′ - c2m1) * (cont4 + (c5′ - c1m1) * cont5)))) - w1 = @.. broadcast=false TI11*z1+TI12*z2+TI13*z3+TI14*z4+TI15*z5 - w2 = @.. broadcast=false TI21*z1+TI22*z2+TI23*z3+TI24*z4+TI25*z5 - w3 = @.. broadcast=false TI31*z1+TI32*z2+TI33*z3+TI34*z4+TI35*z5 - w4 = @.. broadcast=false TI41*z1+TI42*z2+TI43*z3+TI44*z4+TI45*z5 - w5 = @.. broadcast=false TI51*z1+TI52*z2+TI53*z3+TI54*z4+TI55*z5 + @.. w1 = TI11*z1+TI12*z2+TI13*z3+TI14*z4+TI15*z5 + @.. w2 = TI21*z1+TI22*z2+TI23*z3+TI24*z4+TI25*z5 + @.. w3 = TI31*z1+TI32*z2+TI33*z3+TI34*z4+TI35*z5 + @.. w4 = TI41*z1+TI42*z2+TI43*z3+TI44*z4+TI45*z5 + @.. w5 = TI51*z1+TI52*z2+TI53*z3+TI54*z4+TI55*z5 end # Newton iteration @@ -1328,21 +1328,21 @@ end if integrator.EEst <= oneunit(integrator.EEst) cache.dtprev = dt if alg.extrapolant != :constant - cache.cont1 = @.. (z4 - z5) / c4m1 # first derivative on [c4, 1] - tmp1 = @.. (z3 - z4) / c3mc4 # first derivative on [c3, c4] - cache.cont2 = @.. (tmp1 - cache.cont1) / c3m1 # second derivative on [c3, 1] - tmp2 = @.. (z2 - z3) / c2mc3 # first derivative on [c2, c3] - tmp3 = @.. (tmp2 - tmp1) / c2mc4 # second derivative on [c2, c4] - cache.cont3 = @.. (tmp3 - cache.cont2) / c2m1 # third derivative on [c2, 1] - tmp4 = @.. (z1 - z2) / c1mc2 # first derivative on [c1, c2] - tmp5 = @.. (tmp4 - tmp2) / c1mc3 # second derivative on [c1, c3] - tmp6 = @.. (tmp5 - tmp3) / c1mc4 # third derivative on [c1, c4] - cache.cont4 = @.. (tmp6 - cache.cont3) / c1m1 #fourth derivative on [c1, 1] - tmp7 = @.. z1 / c1 #first derivative on [0, c1] - tmp8 = @.. (tmp4 - tmp7) / c2 #second derivative on [0, c2] - tmp9 = @.. (tmp5 - tmp8) / c3 #third derivative on [0, c3] - tmp10 = @.. (tmp6 - tmp9) / c4 #fourth derivative on [0,c4] - cache.cont5 = @.. cache.cont4 - tmp10 #fifth derivative on [0,1] + @.. cache.cont1 = (z4 - z5) / c4m1 # first derivative on [c4, 1] + @.. tmp = (z3 - z4) / c3mc4 # first derivative on [c3, c4] + @.. cache.cont2 = (tmp - cache.cont1) / c3m1 # second derivative on [c3, 1] + @.. tmp2 = (z2 - z3) / c2mc3 # first derivative on [c2, c3] + @.. tmp3 = (tmp2 - tmp) / c2mc4 # second derivative on [c2, c4] + @.. cache.cont3 = (tmp3 - cache.cont2) / c2m1 # third derivative on [c2, 1] + @.. tmp4 = (z1 - z2) / c1mc2 # first derivative on [c1, c2] + @.. tmp5 = (tmp4 - tmp2) / c1mc3 # second derivative on [c1, c3] + @.. tmp6 = (tmp5 - tmp3) / c1mc4 # third derivative on [c1, c4] + @.. cache.cont4 = (tmp6 - cache.cont3) / c1m1 #fourth derivative on [c1, 1] + @.. tmp7 = z1 / c1 #first derivative on [0, c1] + @.. tmp8 = (tmp4 - tmp7) / c2 #second derivative on [0, c2] + @.. tmp9 = (tmp5 - tmp8) / c3 #third derivative on [0, c3] + @.. tmp10 = (tmp6 - tmp9) / c4 #fourth derivative on [0,c4] + @.. cache.cont5 = cache.cont4 - tmp10 #fifth derivative on [0,1] end end @@ -1437,7 +1437,7 @@ end for i in 1 : num_stages z[i] = f(uprev + z[i], p, t + c[i] * dt) end - OrdinaryDiffEqCore.increment_nf!(integrator.stats, 5) + OrdinaryDiffEqCore.increment_nf!(integrator.stats, num_stages) #fw = TI * ff fw = Vector{typeof(u)}(undef, num_stages) @@ -1619,7 +1619,7 @@ end if (new_W = do_newW(integrator, alg, new_jac, cache.W_γdt)) @inbounds for II in CartesianIndices(J) W1[II] = -γdt * mass_matrix[Tuple(II)...] + J[II] - for i in 1 :(num_stages - 1) ÷ 2 + for i in 1 : (num_stages - 1) ÷ 2 W2[i][II] = -(αdt[i] + βdt[i] * im) * mass_matrix[Tuple(II)...] + J[II] end end @@ -1673,7 +1673,7 @@ end @.. tmp = uprev + z[i] f(ks[i], tmp, p, t + c[i] * dt) end - OrdinaryDiffEqCore.increment_nf!(integrator.stats, 5) + OrdinaryDiffEqCore.increment_nf!(integrator.stats, num_stages) #mul!(fw, TI, ks) for i in 1:num_stages @@ -1700,15 +1700,12 @@ end @.. ubuff = fw[1] - γdt * Mw[1] needfactor = iter == 1 && new_W - linsolve1 = cache.linsolve1 if needfactor - linres = dolinsolve(integrator, linsolve1; A = W1, b = _vec(ubuff), linu = _vec(dw1)) + cache.linsolve1 = dolinsolve(integrator, linsolve1; A = W1, b = _vec(ubuff), linu = _vec(dw1)).cache else - linres = dolinsolve(integrator, linsolve1; A = nothing, b = _vec(ubuff), linu = _vec(dw1)) + cache.linsolve1 = dolinsolve(integrator, linsolve1; A = nothing, b = _vec(ubuff), linu = _vec(dw1)).cache end - cache.linsolve1 = linres.cache - for i in 1 :(num_stages - 1) ÷ 2 @.. cubuff[i]=complex( fw[2 * i] - αdt[i] * Mw[2 * i] + βdt[i] * Mw[2 * i + 1], fw[2 * i + 1] - βdt[i] * Mw[2 * i] - αdt[i] * Mw[2 * i + 1]) @@ -1801,9 +1798,8 @@ end @.. broadcast=false ubuff=integrator.fsalfirst + tmp if alg.smooth_est - linres = dolinsolve(integrator, linres.cache; b = _vec(ubuff), - linu = _vec(utilde)) - cache.linsolve1 = linres.cache + cache.linsolve1 = dolinsolve(integrator, linsolve1; b = _vec(ubuff), + linu = _vec(utilde)).cache integrator.stats.nsolve += 1 end @@ -1821,9 +1817,8 @@ end @.. broadcast=false ubuff=fsallast + tmp if alg.smooth_est - linres = dolinsolve(integrator, linres.cache; b = _vec(ubuff), - linu = _vec(utilde)) - cache.linsolve1 = linres.cache + cache.linsolve1 = dolinsolve(integrator, linsolve1; b = _vec(ubuff), + linu = _vec(utilde)).cache integrator.stats.nsolve += 1 end From b6b86de50a9b8eec97f61ee109d1127c34463976 Mon Sep 17 00:00:00 2001 From: Shreyas Ekanathan Date: Tue, 19 Nov 2024 07:42:41 -0500 Subject: [PATCH 05/12] fix tests --- lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl | 6 +++--- lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl b/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl index 72e0be40f4..cbbddf9467 100644 --- a/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl +++ b/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl @@ -17,10 +17,10 @@ sim21 = test_convergence(1 ./ 2 .^ (2.5:-1:0.5), prob_ode_2Dlinear, RadauIIA9()) prob_ode_linear_big = remake(prob_ode_linear, u0 = big.(prob_ode_linear.u0), tspan = big.(prob_ode_linear.tspan)) prob_ode_2Dlinear_big = remake(prob_ode_2Dlinear, u0 = big.(prob_ode_2Dlinear.u0), tspan = big.(prob_ode_2Dlinear.tspan)) -for i in [3, 5, 7], prob in [prob_ode_linear_big, prob_ode_2Dlinear_big] +for i in [5, 9, 13], prob in [prob_ode_linear_big, prob_ode_2Dlinear_big] dts = 1 ./ 2 .^ (4.25:-1:0.25) - sim21 = test_convergence(dts, prob, AdaptiveRadau(min_stages = i, max_stages = i)) - @test sim21.𝒪est[:final]≈ (2 * i - 1) atol=testTol + sim21 = test_convergence(dts, prob, AdaptiveRadau(min_order = i, max_order = i)) + @test sim21.𝒪est[:final]≈ i atol=testTol end # test adaptivity diff --git a/lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl b/lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl index d3c19ebbb5..9b2c5b6ca8 100644 --- a/lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl +++ b/lib/OrdinaryDiffEqFIRKGenerator/test/ode_firk_tests.jl @@ -6,8 +6,8 @@ testTol = 0.5 prob_ode_linear_big = remake(prob_ode_linear, u0 = big.(prob_ode_linear.u0), tspan = big.(prob_ode_linear.tspan)) prob_ode_2Dlinear_big = remake(prob_ode_2Dlinear, u0 = big.(prob_ode_2Dlinear.u0), tspan = big.(prob_ode_2Dlinear.tspan)) -for i in [9], prob in [prob_ode_linear_big, prob_ode_2Dlinear_big] +for i in [17, 21], prob in [prob_ode_linear_big, prob_ode_2Dlinear_big] dts = 1 ./ 2 .^ (4.25:-1:0.25) - sim21 = test_convergence(dts, prob, AdaptiveRadau(min_stages = i, max_stages = i)) - @test sim21.𝒪est[:final]≈ (2 * i - 1) atol=testTol + sim21 = test_convergence(dts, prob, AdaptiveRadau(min_order = i, max_order = i)) + @test sim21.𝒪est[:final]≈ i atol=testTol end From f86d0ea385aefad0839a10d522c46aced9218a51 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 11 Nov 2024 23:09:46 +0530 Subject: [PATCH 06/12] refactor: generalize `_initialize_dae!` to use SciMLBase implementations --- .../src/OrdinaryDiffEqCore.jl | 2 +- lib/OrdinaryDiffEqCore/src/initialize_dae.jl | 123 +++--------------- 2 files changed, 18 insertions(+), 107 deletions(-) diff --git a/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl b/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl index a7cbd95167..ac2e671b21 100644 --- a/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl +++ b/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl @@ -60,7 +60,7 @@ using DiffEqBase: check_error!, @def, _vec, _reshape using FastBroadcast: @.., True, False -using SciMLBase: NoInit, CheckInit, _unwrap_val +using SciMLBase: NoInit, CheckInit, OverrideInit, AbstractDEProblem, _unwrap_val import SciMLBase: alg_order diff --git a/lib/OrdinaryDiffEqCore/src/initialize_dae.jl b/lib/OrdinaryDiffEqCore/src/initialize_dae.jl index ebefa3d91c..8ad99e7f1e 100644 --- a/lib/OrdinaryDiffEqCore/src/initialize_dae.jl +++ b/lib/OrdinaryDiffEqCore/src/initialize_dae.jl @@ -20,16 +20,6 @@ function BrownFullBasicInit(; abstol = 1e-10, nlsolve = nothing) end BrownFullBasicInit(abstol) = BrownFullBasicInit(; abstol = abstol, nlsolve = nothing) -struct OverrideInit{T, F} <: DiffEqBase.DAEInitializationAlgorithm - abstol::T - nlsolve::F -end - -function OverrideInit(; abstol = 1e-10, nlsolve = nothing) - OverrideInit(abstol, nlsolve) -end -OverrideInit(abstol) = OverrideInit(; abstol = abstol, nlsolve = nothing) - ## Notes #= @@ -143,19 +133,15 @@ end ## NoInit -function _initialize_dae!(integrator, prob::Union{ODEProblem, DAEProblem}, +function _initialize_dae!(integrator, prob::AbstractDEProblem, alg::NoInit, x::Union{Val{true}, Val{false}}) end ## OverrideInit -function _initialize_dae!(integrator, prob::Union{ODEProblem, DAEProblem}, +function _initialize_dae!(integrator, prob::AbstractDEProblem, alg::OverrideInit, isinplace::Union{Val{true}, Val{false}}) - initializeprob = prob.f.initializeprob - - if SciMLBase.has_update_initializeprob!(prob.f) - prob.f.update_initializeprob!(initializeprob, prob) - end + initializeprob = prob.f.initialization_data.initializeprob # If it doesn't have autodiff, assume it comes from symbolic system like ModelingToolkit # Since then it's the case of not a DAE but has initializeprob @@ -168,105 +154,30 @@ function _initialize_dae!(integrator, prob::Union{ODEProblem, DAEProblem}, true end - alg = default_nlsolve(alg.nlsolve, isinplace, initializeprob.u0, initializeprob, isAD) - nlsol = solve(initializeprob, alg, abstol = integrator.opts.abstol, reltol = integrator.opts.reltol) + nlsolve_alg = default_nlsolve(alg.nlsolve, isinplace, initializeprob.u0, initializeprob, isAD) + + u0, p, success = SciMLBase.get_initial_values(prob, prob.f, integrator, alg, isinplace; nlsolve_alg, abstol = integrator.opts.abstol, reltol = integrator.opts.reltol) + if isinplace === Val{true}() - integrator.u .= prob.f.initializeprobmap(nlsol) + integrator.u .= u0 elseif isinplace === Val{false}() - integrator.u = prob.f.initializeprobmap(nlsol) + integrator.u = u0 else error("Unreachable reached. Report this error.") end - if SciMLBase.has_initializeprobpmap(prob.f) - integrator.p = prob.f.initializeprobpmap(prob, nlsol) - sol = integrator.sol - @reset sol.prob.p = integrator.p - integrator.sol = sol - end + integrator.p = p + sol = integrator.sol + @reset sol.prob.p = integrator.p + integrator.sol = sol - if nlsol.retcode != ReturnCode.Success + if !success integrator.sol = SciMLBase.solution_new_retcode(integrator.sol, ReturnCode.InitialFailure) end end ## CheckInit -struct CheckInitFailureError <: Exception - normresid::Any - abstol::Any -end - -function Base.showerror(io::IO, e::CheckInitFailureError) - print(io, - "CheckInit specified but initialization not satisifed. normresid = $(e.normresid) > abstol = $(e.abstol)") -end - -function _initialize_dae!(integrator, prob::ODEProblem, alg::CheckInit, - isinplace::Val{true}) - @unpack p, t, f = integrator - M = integrator.f.mass_matrix - tmp = first(get_tmp_cache(integrator)) - u0 = integrator.u - - algebraic_vars = [all(iszero, x) for x in eachcol(M)] - algebraic_eqs = [all(iszero, x) for x in eachrow(M)] - (iszero(algebraic_vars) || iszero(algebraic_eqs)) && return - update_coefficients!(M, u0, p, t) - f(tmp, u0, p, t) - tmp .= ArrayInterface.restructure(tmp, algebraic_eqs .* _vec(tmp)) - - normresid = integrator.opts.internalnorm(tmp, t) - if normresid > integrator.opts.abstol - throw(CheckInitFailureError(normresid, integrator.opts.abstol)) - end -end - -function _initialize_dae!(integrator, prob::ODEProblem, alg::CheckInit, - isinplace::Val{false}) - @unpack p, t, f = integrator - u0 = integrator.u - M = integrator.f.mass_matrix - - algebraic_vars = [all(iszero, x) for x in eachcol(M)] - algebraic_eqs = [all(iszero, x) for x in eachrow(M)] - (iszero(algebraic_vars) || iszero(algebraic_eqs)) && return - update_coefficients!(M, u0, p, t) - du = f(u0, p, t) - resid = _vec(du)[algebraic_eqs] - - normresid = integrator.opts.internalnorm(resid, t) - if normresid > integrator.opts.abstol - throw(CheckInitFailureError(normresid, integrator.opts.abstol)) - end -end - -function _initialize_dae!(integrator, prob::DAEProblem, - alg::CheckInit, isinplace::Val{true}) - @unpack p, t, f = integrator - u0 = integrator.u - resid = get_tmp_cache(integrator)[2] - - f(resid, integrator.du, u0, p, t) - normresid = integrator.opts.internalnorm(resid, t) - if normresid > integrator.opts.abstol - throw(CheckInitFailureError(normresid, integrator.opts.abstol)) - end -end - -function _initialize_dae!(integrator, prob::DAEProblem, - alg::CheckInit, isinplace::Val{false}) - @unpack p, t, f = integrator - u0 = integrator.u - - nlequation_oop = u -> begin - f((u - u0) / dt, u, p, t) - end - - nlequation = (u, _) -> nlequation_oop(u) - - resid = f(integrator.du, u0, p, t) - normresid = integrator.opts.internalnorm(resid, t) - if normresid > integrator.opts.abstol - throw(CheckInitFailureError(normresid, integrator.opts.abstol)) - end +function _initialize_dae!(integrator, prob::AbstractDEProblem, alg::CheckInit, + isinplace::Union{Val{true}, Val{false}}) + SciMLBase.get_initial_values(prob, integrator, prob.f, alg, isinplace; abstol = integrator.opts.abstol) end From 8f9a1b248d298fe7b5357374eaec9e22ede9599e Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Tue, 19 Nov 2024 11:23:15 +0530 Subject: [PATCH 07/12] build: bump SciMLBase compat --- lib/OrdinaryDiffEqCore/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/OrdinaryDiffEqCore/Project.toml b/lib/OrdinaryDiffEqCore/Project.toml index ea0eed4466..1c1476d5e1 100644 --- a/lib/OrdinaryDiffEqCore/Project.toml +++ b/lib/OrdinaryDiffEqCore/Project.toml @@ -70,7 +70,7 @@ Random = "<0.0.1, 1" RecursiveArrayTools = "2.36, 3" Reexport = "1.0" SafeTestsets = "0.1.0" -SciMLBase = "2.60" +SciMLBase = "2.62" SciMLOperators = "0.3" SciMLStructures = "1" SimpleUnPack = "1" From 349bfcdd5f7d7ded70a7504d84b5e78aa8d8c611 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Wed, 20 Nov 2024 13:13:10 +0530 Subject: [PATCH 08/12] test: fix `CheckInit` tests --- test/interface/checkinit_tests.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/interface/checkinit_tests.jl b/test/interface/checkinit_tests.jl index 8320449151..e71b412e2c 100644 --- a/test/interface/checkinit_tests.jl +++ b/test/interface/checkinit_tests.jl @@ -24,9 +24,9 @@ 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)) -@test_throws OrdinaryDiffEqCore.CheckInitFailureError solve( +@test_throws SciMLBase.CheckInitFailureError solve( prob_mm, Rodas5P(), reltol = 1e-8, abstol = 1e-8, initializealg = SciMLBase.CheckInit()) -@test_throws OrdinaryDiffEqCore.CheckInitFailureError solve( +@test_throws SciMLBase.CheckInitFailureError solve( prob_mm_oop, Rodas5P(), reltol = 1e-8, abstol = 1e-8, initializealg = SciMLBase.CheckInit()) @@ -49,7 +49,7 @@ tspan = (0.0, 100000.0) differential_vars = [true, true, false] prob = DAEProblem(f, du₀, u₀, tspan, differential_vars = differential_vars) prob_oop = DAEProblem(f_oop, du₀, u₀, tspan, differential_vars = differential_vars) -@test_throws OrdinaryDiffEqCore.CheckInitFailureError solve( +@test_throws SciMLBase.CheckInitFailureError solve( prob, DFBDF(), reltol = 1e-8, abstol = 1e-8, initializealg = SciMLBase.CheckInit()) -@test_throws OrdinaryDiffEqCore.CheckInitFailureError solve( +@test_throws SciMLBase.CheckInitFailureError solve( prob_oop, DFBDF(), reltol = 1e-8, abstol = 1e-8, initializealg = SciMLBase.CheckInit()) From c7b28ff3d56941539e1b7324e93c741595faef1f Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 20 Nov 2024 06:10:39 -0500 Subject: [PATCH 09/12] Update CI.yml --- .github/workflows/CI.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 7d11cf320e..1c05dbf029 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -95,4 +95,4 @@ jobs: with: token: ${{ secrets.CODECOV_TOKEN }} file: lcov.info - fail_ci_if_error: true + fail_ci_if_error: false From e01d9c29f16237cb330e8766ea81733b64784984 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 20 Nov 2024 09:36:17 -0500 Subject: [PATCH 10/12] Update Project.toml --- lib/OrdinaryDiffEqCore/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/OrdinaryDiffEqCore/Project.toml b/lib/OrdinaryDiffEqCore/Project.toml index 1c1476d5e1..f6d4f5e9b2 100644 --- a/lib/OrdinaryDiffEqCore/Project.toml +++ b/lib/OrdinaryDiffEqCore/Project.toml @@ -1,7 +1,7 @@ name = "OrdinaryDiffEqCore" uuid = "bbf590c4-e513-4bbe-9b18-05decba2e5d8" authors = ["ParamThakkar123 "] -version = "1.11.0" +version = "1.12.0" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From 7e4c00544cd779417600996b50a8beeba3c16b55 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 20 Nov 2024 09:49:07 -0500 Subject: [PATCH 11/12] Update Project.toml --- lib/OrdinaryDiffEqFIRKGenerator/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/OrdinaryDiffEqFIRKGenerator/Project.toml b/lib/OrdinaryDiffEqFIRKGenerator/Project.toml index 3e4156ddea..d212271061 100644 --- a/lib/OrdinaryDiffEqFIRKGenerator/Project.toml +++ b/lib/OrdinaryDiffEqFIRKGenerator/Project.toml @@ -1,7 +1,7 @@ name = "OrdinaryDiffEqFIRKGenerator" uuid = "677d4f02-548a-44fa-8eaf-26579094acaf" authors = ["ParamThakkar123 "] -version = "1.0.0" +version = "1.1.0" [deps] GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a" From 79ff531b8e081defb85b25756a2359155528f146 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 20 Nov 2024 09:50:21 -0500 Subject: [PATCH 12/12] Update Project.toml --- lib/OrdinaryDiffEqFIRK/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/OrdinaryDiffEqFIRK/Project.toml b/lib/OrdinaryDiffEqFIRK/Project.toml index 0525c0059e..9e8d3373ba 100644 --- a/lib/OrdinaryDiffEqFIRK/Project.toml +++ b/lib/OrdinaryDiffEqFIRK/Project.toml @@ -1,7 +1,7 @@ name = "OrdinaryDiffEqFIRK" uuid = "5960d6e9-dd7a-4743-88e7-cf307b64f125" authors = ["ParamThakkar123 "] -version = "1.4.0" +version = "1.5.0" [deps] DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"