diff --git a/lib/OrdinaryDiffEqFIRK/Project.toml b/lib/OrdinaryDiffEqFIRK/Project.toml index e11da0544c..2ffd2e2710 100644 --- a/lib/OrdinaryDiffEqFIRK/Project.toml +++ b/lib/OrdinaryDiffEqFIRK/Project.toml @@ -6,15 +6,20 @@ version = "1.1.1" [deps] DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898" +GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a" +GenericSchur = "c145ed77-6b09-5dd9-b285-bf645a82121e" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8" OrdinaryDiffEqDifferentiation = "4302a76b-040a-498a-8c04-15b101fed76b" OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8" +Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" +RootedTrees = "47965b36-3f3e-11e9-0dcf-4570dfd42a8c" SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961" +Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" [compat] DiffEqBase = "6.152.2" diff --git a/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl b/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl index af3e4a7dd2..40305282d7 100644 --- a/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl +++ b/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl @@ -18,6 +18,7 @@ import OrdinaryDiffEqCore: alg_order, calculate_residuals!, get_current_adaptive_order, get_fsalfirstlast, isfirk, generic_solver_docstring using MuladdMacro, DiffEqBase, RecursiveArrayTools +using Polynomials, GenericLinearAlgebra, GenericSchur using SciMLOperators: AbstractSciMLOperator using LinearAlgebra: I, UniformScaling, mul!, lu import LinearSolve @@ -42,6 +43,6 @@ include("firk_tableaus.jl") include("firk_perform_step.jl") include("integrator_interface.jl") -export RadauIIA3, RadauIIA5, RadauIIA9 +export RadauIIA3, RadauIIA5, RadauIIA9, AdaptiveRadau end diff --git a/lib/OrdinaryDiffEqFIRK/src/alg_utils.jl b/lib/OrdinaryDiffEqFIRK/src/alg_utils.jl index 692a86bb88..428393e0bf 100644 --- a/lib/OrdinaryDiffEqFIRK/src/alg_utils.jl +++ b/lib/OrdinaryDiffEqFIRK/src/alg_utils.jl @@ -3,10 +3,12 @@ qmax_default(alg::Union{RadauIIA3, RadauIIA5, RadauIIA9}) = 8 alg_order(alg::RadauIIA3) = 3 alg_order(alg::RadauIIA5) = 5 alg_order(alg::RadauIIA9) = 9 +alg_order(alg::AdaptiveRadau) = 5 isfirk(alg::RadauIIA3) = true isfirk(alg::RadauIIA5) = true isfirk(alg::RadauIIA9) = true +isfirk(alg::AdaptiveRadau) = true alg_adaptive_order(alg::RadauIIA3) = 1 alg_adaptive_order(alg::RadauIIA5) = 3 diff --git a/lib/OrdinaryDiffEqFIRK/src/algorithms.jl b/lib/OrdinaryDiffEqFIRK/src/algorithms.jl index 844f50ea96..25165422c8 100644 --- a/lib/OrdinaryDiffEqFIRK/src/algorithms.jl +++ b/lib/OrdinaryDiffEqFIRK/src/algorithms.jl @@ -150,3 +150,42 @@ function RadauIIA9(; chunk_size = Val{0}(), autodiff = Val{true}(), controller, step_limiter!) end + +struct AdaptiveRadau{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <: + OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + precs::P + smooth_est::Bool + extrapolant::Symbol + κ::Tol + maxiters::Int + fast_convergence_cutoff::C1 + new_W_γdt_cutoff::C2 + controller::Symbol + step_limiter!::StepLimiter + num_stages::Int +end + +function AdaptiveRadau(; chunk_size = Val{0}(), autodiff = Val{true}(), + standardtag = Val{true}(), concrete_jac = nothing, + diff_type = Val{:forward}, num_stages = 3, + linsolve = nothing, precs = DEFAULT_PRECS, + extrapolant = :dense, fast_convergence_cutoff = 1 // 5, + new_W_γdt_cutoff = 1 // 5, + controller = :Predictive, κ = nothing, maxiters = 10, smooth_est = true, + step_limiter! = trivial_limiter!) + AdaptiveRadau{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac), + typeof(κ), typeof(fast_convergence_cutoff), + typeof(new_W_γdt_cutoff), typeof(step_limiter!)}(linsolve, + precs, + smooth_est, + extrapolant, + κ, + maxiters, + fast_convergence_cutoff, + new_W_γdt_cutoff, + controller, + step_limiter!, num_stages) +end + diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl index 47fdad5374..af6da2a282 100644 --- a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl +++ b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl @@ -287,6 +287,7 @@ mutable struct RadauIIA9ConstantCache{F, Tab, Tol, Dt, U, JType} <: cont2::U cont3::U cont4::U + cont5::U dtprev::Dt W_γdt::Dt status::NLStatus @@ -304,7 +305,7 @@ function alg_cache(alg::RadauIIA9, u, rate_prototype, ::Type{uEltypeNoUnits}, κ = alg.κ !== nothing ? convert(uToltype, alg.κ) : convert(uToltype, 1 // 100) J = false .* _vec(rate_prototype) .* _vec(rate_prototype)' - RadauIIA9ConstantCache(uf, tab, κ, one(uToltype), 10000, u, u, u, u, dt, dt, + RadauIIA9ConstantCache(uf, tab, κ, one(uToltype), 10000, u, u, u, u, u, dt, dt, Convergence, J) end @@ -333,6 +334,7 @@ mutable struct RadauIIA9Cache{uType, cuType, uNoUnitsType, rateType, JType, W1Ty cont2::uType cont3::uType cont4::uType + cont5::uType du1::rateType fsalfirst::rateType k::rateType @@ -407,6 +409,7 @@ function alg_cache(alg::RadauIIA9, u, rate_prototype, ::Type{uEltypeNoUnits}, cont2 = zero(u) cont3 = zero(u) cont4 = zero(u) + cont5 = zero(u) fsalfirst = zero(rate_prototype) k = zero(rate_prototype) @@ -462,7 +465,7 @@ function alg_cache(alg::RadauIIA9, u, rate_prototype, ::Type{uEltypeNoUnits}, RadauIIA9Cache(u, uprev, z1, z2, z3, z4, z5, w1, w2, w3, w4, w5, - dw1, ubuff, dw23, dw45, cubuff1, cubuff2, cont1, cont2, cont3, cont4, + dw1, ubuff, dw23, dw45, cubuff1, cubuff2, cont1, cont2, cont3, cont4, cont5, du1, fsalfirst, k, k2, k3, k4, k5, fw1, fw2, fw3, fw4, fw5, J, W1, W2, W3, uf, tab, κ, one(uToltype), 10000, @@ -470,3 +473,185 @@ function alg_cache(alg::RadauIIA9, u, rate_prototype, ::Type{uEltypeNoUnits}, linsolve1, linsolve2, linsolve3, rtol, atol, dt, dt, Convergence, alg.step_limiter!) end + +mutable struct AdaptiveRadauConstantCache{F, Tab, Tol, Dt, U, JType} <: + OrdinaryDiffEqConstantCache + uf::F + tab::Tab + κ::Tol + ηold::Tol + iter::Int + cont::Vector{U} + dtprev::Dt + W_γdt::Dt + status::NLStatus + J::JType +end + +function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, + ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck, + ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + uf = UDerivativeWrapper(f, t, p) + uToltype = constvalue(uBottomEltypeNoUnits) + num_stages = alg.num_stages + + if (num_stages == 3) + tab = BigRadauIIA5Tableau(uToltype, constvalue(tTypeNoUnits)) + elseif (num_stages == 5) + tab = BigRadauIIA9Tableau(uToltype, constvalue(tTypeNoUnits)) + elseif (num_stages == 7) + tab = BigRadauIIA13Tableau(uToltype, constvalue(tTypeNoUnits)) + elseif iseven(num_stages) || num_stages <3 + error("num_stages must be odd and 3 or greater") + else + tab = adaptiveRadauTableau(uToltype, constvalue(tTypeNoUnits), num_stages) + end + + cont = Vector{typeof(u)}(undef, num_stages) + for i in 1: num_stages + cont[i] = zero(u) + end + + κ = alg.κ !== nothing ? convert(uToltype, alg.κ) : convert(uToltype, 1 // 100) + J = false .* _vec(rate_prototype) .* _vec(rate_prototype)' + + AdaptiveRadauConstantCache(uf, tab, κ, one(uToltype), 10000, cont, dt, dt, + Convergence, J) +end + +mutable struct AdaptiveRadauCache{uType, cuType, tType, uNoUnitsType, rateType, JType, W1Type, W2Type, + UF, JC, F1, F2, Tab, Tol, Dt, rTol, aTol, StepLimiter} <: + FIRKMutableCache + u::uType + uprev::uType + z::Vector{uType} + w::Vector{uType} + c_prime::Vector{tType} + dw1::uType + ubuff::uType + dw2::Vector{cuType} + cubuff::Vector{cuType} + dw::Vector{uType} + cont::Vector{uType} + derivatives:: Matrix{uType} + du1::rateType + fsalfirst::rateType + ks::Vector{rateType} + k::rateType + fw::Vector{rateType} + J::JType + W1::W1Type #real + W2::Vector{W2Type} #complex + uf::UF + tab::Tab + κ::Tol + ηold::Tol + iter::Int + tmp::uType + atmp::uNoUnitsType + jac_config::JC + linsolve1::F1 #real + linsolve2::Vector{F2} #complex + rtol::rTol + atol::aTol + dtprev::Dt + W_γdt::Dt + status::NLStatus + step_limiter!::StepLimiter +end + +function alg_cache(alg::AdaptiveRadau, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, + ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck, + ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + uf = UJacobianWrapper(f, t, p) + uToltype = constvalue(uBottomEltypeNoUnits) + num_stages = alg.num_stages + + if (num_stages == 3) + tab = BigRadauIIA5Tableau(uToltype, constvalue(tTypeNoUnits)) + elseif (num_stages == 5) + tab = BigRadauIIA9Tableau(uToltype, constvalue(tTypeNoUnits)) + elseif (num_stages == 7) + tab = BigRadauIIA13Tableau(uToltype, constvalue(tTypeNoUnits)) + elseif iseven(num_stages) || num_stages < 3 + error("num_stages must be odd and 3 or greater") + else + tab = adaptiveRadauTableau(uToltype, constvalue(tTypeNoUnits), num_stages) + end + + κ = alg.κ !== nothing ? convert(uToltype, alg.κ) : convert(uToltype, 1 // 100) + + z = Vector{typeof(u)}(undef, num_stages) + w = Vector{typeof(u)}(undef, num_stages) + for i in 1 : num_stages + z[i] = w[i] = zero(u) + end + + c_prime = Vector{typeof(t)}(undef, num_stages) #time stepping + + dw1 = zero(u) + ubuff = zero(u) + dw2 = [similar(u, Complex{eltype(u)}) for _ in 1 : (num_stages - 1) ÷ 2] + recursivefill!.(dw2, false) + cubuff = [similar(u, Complex{eltype(u)}) for _ in 1 : (num_stages - 1) ÷ 2] + recursivefill!.(cubuff, false) + dw = Vector{typeof(u)}(undef, num_stages - 1) + + cont = Vector{typeof(u)}(undef, num_stages) + for i in 1 : num_stages + cont[i] = zero(u) + end + + derivatives = Matrix{typeof(u)}(undef, num_stages, num_stages) + for i in 1 : num_stages, j in 1 : num_stages + derivatives[i, j] = zero(u) + end + + fsalfirst = zero(rate_prototype) + fw = Vector{typeof(rate_prototype)}(undef, num_stages) + ks = Vector{typeof(rate_prototype)}(undef, num_stages) + for i in 1: num_stages + ks[i] = fw[i] = zero(rate_prototype) + end + k = ks[1] + + J, W1 = build_J_W(alg, u, uprev, p, t, dt, f, uEltypeNoUnits, Val(true)) + if J isa AbstractSciMLOperator + error("Non-concrete Jacobian not yet supported by AdaptiveRadau.") + end + + W2 = [similar(J, Complex{eltype(W1)}) for _ in 1 : (num_stages - 1) ÷ 2] + recursivefill!.(W2, false) + + du1 = zero(rate_prototype) + + tmp = zero(u) + + atmp = similar(u, uEltypeNoUnits) + recursivefill!(atmp, false) + + jac_config = build_jac_config(alg, f, uf, du1, uprev, u, zero(u), dw1) + + linprob = LinearProblem(W1, _vec(ubuff); u0 = _vec(dw1)) + linsolve1 = init(linprob, alg.linsolve, alias_A = true, alias_b = true, + assumptions = LinearSolve.OperatorAssumptions(true)) + + linsolve2 = [ + init(LinearProblem(W2[i], _vec(cubuff[i]); u0 = _vec(dw2[i])), alg.linsolve, alias_A = true, alias_b = true, + assumptions = LinearSolve.OperatorAssumptions(true)) for i in 1 : (num_stages - 1) ÷ 2] + + rtol = reltol isa Number ? reltol : zero(reltol) + atol = reltol isa Number ? reltol : zero(reltol) + + AdaptiveRadauCache(u, uprev, + z, w, c_prime, dw1, ubuff, dw2, cubuff, dw, cont, derivatives, + du1, fsalfirst, ks, k, fw, + J, W1, W2, + uf, tab, κ, one(uToltype), 10000, tmp, + atmp, jac_config, + linsolve1, linsolve2, rtol, atol, dt, dt, + Convergence, alg.step_limiter!) +end + diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl b/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl index 6789b307f6..044257afad 100644 --- a/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl +++ b/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl @@ -26,19 +26,7 @@ function do_newW(integrator, nlsolver, new_jac, W_dt)::Bool # for FIRK return !smallstepchange end -function initialize!(integrator, cache::RadauIIA3ConstantCache) - integrator.kshortsize = 2 - integrator.k = typeof(integrator.k)(undef, integrator.kshortsize) - integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal - - # Avoid undefined entries if k is an array of arrays - integrator.fsallast = zero(integrator.fsalfirst) - integrator.k[1] = integrator.fsalfirst - integrator.k[2] = integrator.fsallast - nothing -end - -function initialize!(integrator, cache::RadauIIA5ConstantCache) +function initialize!(integrator, cache::Union{RadauIIA3ConstantCache, RadauIIA5ConstantCache, RadauIIA9ConstantCache,AdaptiveRadauConstantCache}) integrator.kshortsize = 2 integrator.k = typeof(integrator.k)(undef, integrator.kshortsize) integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal @@ -51,30 +39,37 @@ function initialize!(integrator, cache::RadauIIA5ConstantCache) nothing end -function initialize!(integrator, cache::RadauIIA9ConstantCache) +function initialize!(integrator, cache::RadauIIA3Cache) integrator.kshortsize = 2 - integrator.k = typeof(integrator.k)(undef, integrator.kshortsize) - integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal - OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) - - # Avoid undefined entries if k is an array of arrays - integrator.fsallast = zero(integrator.fsalfirst) + resize!(integrator.k, integrator.kshortsize) integrator.k[1] = integrator.fsalfirst integrator.k[2] = integrator.fsallast + integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) + OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) nothing end -function initialize!(integrator, cache::RadauIIA3Cache) +function initialize!(integrator, cache::RadauIIA5Cache) integrator.kshortsize = 2 resize!(integrator.k, integrator.kshortsize) integrator.k[1] = integrator.fsalfirst integrator.k[2] = integrator.fsallast integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) + if integrator.opts.adaptive + @unpack abstol, reltol = integrator.opts + if reltol isa Number + cache.rtol = reltol^(2 / 3) / 10 + cache.atol = cache.rtol * (abstol / reltol) + else + @.. broadcast=false cache.rtol=reltol^(2 / 3) / 10 + @.. broadcast=false cache.atol=cache.rtol * (abstol / reltol) + end + end nothing end -function initialize!(integrator, cache::RadauIIA5Cache) +function initialize!(integrator, cache::RadauIIA9Cache) integrator.kshortsize = 2 resize!(integrator.k, integrator.kshortsize) integrator.k[1] = integrator.fsalfirst @@ -84,17 +79,18 @@ function initialize!(integrator, cache::RadauIIA5Cache) if integrator.opts.adaptive @unpack abstol, reltol = integrator.opts if reltol isa Number - cache.rtol = reltol^(2 / 3) / 10 + cache.rtol = reltol^(3 / 5) / 10 cache.atol = cache.rtol * (abstol / reltol) else - @.. broadcast=false cache.rtol=reltol^(2 / 3) / 10 + @.. broadcast=false cache.rtol=reltol^(3 / 5) / 10 @.. broadcast=false cache.atol=cache.rtol * (abstol / reltol) end end nothing end -function initialize!(integrator, cache::RadauIIA9Cache) +function initialize!(integrator, cache::AdaptiveRadauCache) + @unpack num_stages = cache.tab integrator.kshortsize = 2 resize!(integrator.k, integrator.kshortsize) integrator.k[1] = integrator.fsalfirst @@ -104,16 +100,17 @@ function initialize!(integrator, cache::RadauIIA9Cache) if integrator.opts.adaptive @unpack abstol, reltol = integrator.opts if reltol isa Number - cache.rtol = reltol^(5 / 8) / 10 + cache.rtol = reltol^((num_stages + 1) / (2 * num_stages)) / 10 cache.atol = cache.rtol * (abstol / reltol) else - @.. broadcast=false cache.rtol=reltol^(5 / 8) / 10 + @.. broadcast=false cache.rtol=reltol^((num_stages + 1) / (2 * num_stages)) / 10 @.. broadcast=false cache.atol=cache.rtol * (abstol / reltol) end end nothing end + @muladd function perform_step!(integrator, cache::RadauIIA3ConstantCache) @unpack t, dt, uprev, u, f, p = integrator @unpack T11, T12, T21, T22, TI11, TI12, TI21, TI22 = cache.tab @@ -125,7 +122,7 @@ end mass_matrix = integrator.f.mass_matrix # precalculations - rtol = @. reltol^(2 / 3) / 10 + rtol = @. reltol^(3 / 4) / 10 atol = @. rtol * (abstol / reltol) αdt, βdt = α / dt, β / dt J = calc_J(integrator, cache) @@ -776,6 +773,7 @@ end @.. broadcast=false cache.cont3=cache.cont2 - (tmp - z1 / c1) / c2 end end + f(fsallast, u, p, t + dt) OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) return @@ -787,14 +785,14 @@ end @unpack T11, T12, T13, T14, T15, T21, T22, T23, T24, T25, T31, T32, T33, T34, T35, T41, T42, T43, T44, T45, T51 = cache.tab #= T52 = 1, T53 = 0, T54 = 1, T55 = 0=# @unpack TI11, TI12, TI13, TI14, TI15, TI21, TI22, TI23, TI24, TI25, TI31, TI32, TI33, TI34, TI35, TI41, TI42, TI43, TI44, TI45, TI51, TI52, TI53, TI54, TI55 = cache.tab @unpack c1, c2, c3, c4, γ, α1, β1, α2, β2, e1, e2, e3, e4, e5 = cache.tab - @unpack κ, cont1, cont2, cont3, cont4 = cache + @unpack κ, cont1, cont2, cont3, cont4, cont5 = cache @unpack internalnorm, abstol, reltol, adaptive = integrator.opts alg = unwrap_alg(integrator, true) @unpack maxiters = alg mass_matrix = integrator.f.mass_matrix # precalculations rtol pow is (num stages + 1)/(2*num stages) - rtol = @.. broadcast=false reltol^(5 / 8)/10 + rtol = @.. broadcast=false reltol^(3 / 5)/10 atol = @.. broadcast=false rtol*(abstol / reltol) c1m1 = c1 - 1 c2m1 = c2 - 1 @@ -832,27 +830,32 @@ end cache.cont2 = map(zero, u) cache.cont3 = map(zero, u) cache.cont4 = map(zero, u) + cache.cont5 = map(zero, u) else c5′ = dt / cache.dtprev c1′ = c1 * c5′ c2′ = c2 * c5′ c3′ = c3 * c5′ c4′ = c4 * c5′ - z1 = @.. broadcast=false c1′*(cont1 + - (c1′ - c3m1) * (cont2 + - (c1′ - c2m1) * (cont3 + (c1′ - c1m1) * cont4))) - z2 = @.. broadcast=false c2′*(cont1 + - (c2′ - c3m1) * (cont2 + - (c2′ - c2m1) * (cont3 + (c2′ - c1m1) * cont4))) - z3 = @.. broadcast=false c3′*(cont1 + - (c3′ - c3m1) * (cont2 + - (c3′ - c2m1) * (cont3 + (c3′ - c1m1) * cont4))) - z4 = @.. broadcast=false c4′*(cont1 + - (c4′ - c3m1) * (cont2 + - (c4′ - c2m1) * (cont3 + (c4′ - c1m1) * cont4))) - z5 = @.. broadcast=false c5′*(cont1 + - (c5′ - c3m1) * (cont2 + - (c5′ - c2m1) * (cont3 + (c5′ - c1m1) * cont4))) + z1 = @.. c1′ * (cont1 + + (c1′-c4m1) * (cont2 + + (c1′ - c3m1) * (cont3 + + (c1′ - c2m1) * (cont4 + (c1′ - c1m1) * cont5)))) + z2 = @.. c2′ * (cont1 + + (c2′-c4m1) * (cont2 + + (c2′ - c3m1) * (cont3 + + (c2′ - c2m1) * (cont4 + (c2′ - c1m1) * cont5)))) + z3 = @.. c3′ * (cont1 + + (c3′-c4m1) * (cont2 + + (c3′ - c3m1) * (cont3 + + (c3′ - c2m1) * (cont4 + (c3′ - c1m1) * cont5)))) + z4 = @.. c4′ * (cont1 + + (c4′-c4m1) * (cont2 + + (c4′ - c3m1) * (cont3 + + (c4′ - c2m1) * (cont4 + (c4′ - c1m1) * cont5)))) + 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 @@ -1003,6 +1006,11 @@ end 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 @@ -1019,12 +1027,12 @@ end @unpack T11, T12, T13, T14, T15, T21, T22, T23, T24, T25, T31, T32, T33, T34, T35, T41, T42, T43, T44, T45, T51 = cache.tab #= T52 = 1, T53 = 0, T54 = 1, T55 = 0=# @unpack TI11, TI12, TI13, TI14, TI15, TI21, TI22, TI23, TI24, TI25, TI31, TI32, TI33, TI34, TI35, TI41, TI42, TI43, TI44, TI45, TI51, TI52, TI53, TI54, TI55 = cache.tab @unpack c1, c2, c3, c4, γ, α1, β1, α2, β2, e1, e2, e3, e4, e5 = cache.tab - @unpack κ, cont1, cont2, cont3, cont4 = cache + @unpack κ, cont1, cont2, cont3, cont4, cont5 = cache @unpack z1, z2, z3, z4, z5, w1, w2, w3, w4, w5 = cache @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, rtol, atol, step_limiter! = cache + @unpack tmp, tmp2, tmp3, tmp4, tmp5, tmp6, 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 @@ -1072,27 +1080,32 @@ end @.. broadcast=false cache.cont2=uzero @.. broadcast=false cache.cont3=uzero @.. broadcast=false cache.cont4=uzero + @.. broadcast=false cache.cont5=uzero else c5′ = dt / cache.dtprev c1′ = c1 * c5′ c2′ = c2 * c5′ c3′ = c3 * c5′ c4′ = c4 * c5′ - z1 = @.. broadcast=false c1′*(cont1 + - (c1′ - c3m1) * (cont2 + - (c1′ - c2m1) * (cont3 + (c1′ - c1m1) * cont4))) - z2 = @.. broadcast=false c2′*(cont1 + - (c2′ - c3m1) * (cont2 + - (c2′ - c2m1) * (cont3 + (c2′ - c1m1) * cont4))) - z3 = @.. broadcast=false c3′*(cont1 + - (c3′ - c3m1) * (cont2 + - (c3′ - c2m1) * (cont3 + (c3′ - c1m1) * cont4))) - z4 = @.. broadcast=false c4′*(cont1 + - (c4′ - c3m1) * (cont2 + - (c4′ - c2m1) * (cont3 + (c4′ - c1m1) * cont4))) - z5 = @.. broadcast=false c5′*(cont1 + - (c5′ - c3m1) * (cont2 + - (c5′ - c2m1) * (cont3 + (c5′ - c1m1) * cont4))) + z1 = @.. c1′ * (cont1 + + (c1′-c4m1) * (cont2 + + (c1′ - c3m1) * (cont3 + + (c1′ - c2m1) * (cont4 + (c1′ - c1m1) * cont5)))) + z2 = @.. c2′ * (cont1 + + (c2′-c4m1) * (cont2 + + (c2′ - c3m1) * (cont3 + + (c2′ - c2m1) * (cont4 + (c2′ - c1m1) * cont5)))) + z3 = @.. c3′ * (cont1 + + (c3′-c4m1) * (cont2 + + (c3′ - c3m1) * (cont3 + + (c3′ - c2m1) * (cont4 + (c3′ - c1m1) * cont5)))) + z4 = @.. c4′ * (cont1 + + (c4′-c4m1) * (cont2 + + (c4′ - c3m1) * (cont3 + + (c4′ - c2m1) * (cont4 + (c4′ - c1m1) * cont5)))) + 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 @@ -1169,11 +1182,9 @@ end linsolve1 = cache.linsolve1 if needfactor - linres1 = dolinsolve(integrator, linsolve1; A = W1, b = _vec(ubuff), - linu = _vec(dw1)) + linres1 = dolinsolve(integrator, linsolve1; A = W1, b = _vec(ubuff), linu = _vec(dw1)) else - linres1 = dolinsolve(integrator, linsolve1; A = nothing, b = _vec(ubuff), - linu = _vec(dw1)) + linres1 = dolinsolve(integrator, linsolve1; A = nothing, b = _vec(ubuff), linu = _vec(dw1)) end cache.linsolve1 = linres1.cache @@ -1184,11 +1195,9 @@ end linsolve2 = cache.linsolve2 if needfactor - linres2 = dolinsolve(integrator, linsolve2; A = W2, b = _vec(cubuff1), - linu = _vec(dw23)) + linres2 = dolinsolve(integrator, linsolve2; A = W2, b = _vec(cubuff1), linu = _vec(dw23)) else - linres2 = dolinsolve(integrator, linsolve2; A = nothing, b = _vec(cubuff1), - linu = _vec(dw23)) + linres2 = dolinsolve(integrator, linsolve2; A = nothing, b = _vec(cubuff1), linu = _vec(dw23)) end cache.linsolve2 = linres2.cache @@ -1199,15 +1208,12 @@ end linsolve3 = cache.linsolve3 if needfactor - linres3 = dolinsolve(integrator, linsolve3; A = W3, b = _vec(cubuff2), - linu = _vec(dw45)) + linres3 = dolinsolve(integrator, linsolve3; A = W3, b = _vec(cubuff2), linu = _vec(dw45)) else - linres3 = dolinsolve(integrator, linsolve3; A = nothing, b = _vec(cubuff2), - linu = _vec(dw45)) + linres3 = dolinsolve(integrator, linsolve3; A = nothing, b = _vec(cubuff2), linu = _vec(dw45)) end cache.linsolve3 = linres3.cache - integrator.stats.nsolve += 3 dw2 = z2 dw3 = z3 @@ -1322,16 +1328,21 @@ end if integrator.EEst <= oneunit(integrator.EEst) cache.dtprev = dt if alg.extrapolant != :constant - @.. cache.cont1 = (z4 - z5) / c4m1 - @.. tmp = (z3 - z4) / c3mc4 - @.. cache.cont2 = (tmp - cache.cont1) / c3m1 - @.. tmp2 = (z2 - z3) / c2mc3 - @.. tmp3 = (tmp2 - tmp) / c2mc4 - @.. cache.cont3 = (tmp3 - cache.cont2) / c2m1 - @.. tmp4 = (z1 - z2) / c1mc2 - @.. tmp5 = (tmp4 - tmp2) / c1mc3 - @.. tmp6 = (tmp5 - tmp3) / c1mc4 - @.. cache.cont4 = (tmp6 - cache.cont3) / c1m1 + 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] end end @@ -1339,3 +1350,495 @@ end OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) return end + +@muladd function perform_step!(integrator, cache::AdaptiveRadauConstantCache, + repeat_step = false) + @unpack t, dt, uprev, u, f, p = integrator + @unpack T, TI, γ, α, β, c, e, num_stages = cache.tab + @unpack κ, cont = cache + @unpack internalnorm, abstol, reltol, adaptive = integrator.opts + alg = unwrap_alg(integrator, true) + @unpack maxiters = alg + mass_matrix = integrator.f.mass_matrix + + # precalculations rtol pow is (num stages + 1)/(2*num stages) + rtol = @.. broadcast=false reltol^((num_stages + 1) / (num_stages * 2))/10 + atol = @.. broadcast=false rtol*(abstol / reltol) + + γdt, αdt, βdt = γ / dt, α ./ dt, β ./ dt + + J = calc_J(integrator, cache) + if u isa Number + LU1 = -γdt * mass_matrix + J + tmp = -(αdt[1] + βdt[1] * im) * mass_matrix + J + else + LU1 = lu(-γdt * mass_matrix + J) + tmp = lu(-(αdt[1] + βdt[1] * im) * mass_matrix + J) + end + LU2 = Vector{typeof(tmp)}(undef, (num_stages - 1) ÷ 2) + LU2[1] = tmp + if u isa Number + for i in 2 : (num_stages - 1) ÷ 2 + LU2[i] = -(αdt[i] + βdt[i] * im) * mass_matrix + J + end + else + for i in 2 : (num_stages - 1) ÷ 2 + LU2[i] = lu(-(αdt[i] + βdt[i] * im) * mass_matrix + J) + end + end + + integrator.stats.nw += 1 + z = Vector{typeof(u)}(undef, num_stages) + w = Vector{typeof(u)}(undef, num_stages) + if integrator.iter == 1 || integrator.u_modified || alg.extrapolant == :constant + cache.dtprev = one(cache.dtprev) + for i in 1 : num_stages + z[i] = @.. map(zero, u) + w[i] = @.. map(zero, u) + cache.cont[i] = @.. map(zero, u) + end + else + c_prime = Vector{typeof(u)}(undef, num_stages) #time stepping + c_prime[num_stages] = @.. dt / cache.dtprev + for i in 1 : num_stages - 1 + c_prime[i] = @.. c[i] * c_prime[num_stages] + end + for i in 1 : num_stages # collocation polynomial + z[i] = @.. cont[num_stages] * (c_prime[i] - c[1] + 1) + cont[num_stages - 1] + j = num_stages - 2 + while j > 0 + z[i] = @.. z[i] * (c_prime[i] - c[num_stages - j] + 1) + cont[j] + j = j - 1 + end + z[i] = @.. z[i] * c_prime[i] + end + #w = TI*z + for i in 1:num_stages + w[i] = @.. zero(u) + for j in 1:num_stages + w[i] = @.. w[i] + TI[i,j] * z[j] + end + end + end + + # Newton iteration + local ndw + η = max(cache.ηold, eps(eltype(integrator.opts.reltol)))^(0.8) + fail_convergence = true + iter = 0 + while iter < maxiters + iter += 1 + integrator.stats.nnonliniter += 1 + + # evaluate function + ff = Vector{typeof(u)}(undef, num_stages) + for i in 1 : num_stages + ff[i] = f(uprev + z[i], p, t + c[i] * dt) + end + integrator.stats.nf += num_stages + + #fw = TI * ff + fw = Vector{typeof(u)}(undef, num_stages) + for i in 1 : num_stages + fw[i] = @.. zero(u) + for j in 1:num_stages + fw[i] = @.. fw[i] + TI[i,j] * ff[j] + end + end + + Mw = Vector{typeof(u)}(undef, num_stages) + if mass_matrix isa UniformScaling # `UniformScaling` doesn't play nicely with broadcast + for i in 1 : num_stages + Mw[i] = @.. mass_matrix.λ * w[i] #scaling by eigenvalue + end + else + Mw = mass_matrix * w #standard multiplication + end + + rhs = Vector{typeof(u)}(undef, num_stages) + rhs[1] = @.. fw[1] - γdt * Mw[1] + i = 2 + while i <= num_stages #block by block multiplication + rhs[i] = @.. fw[i] - αdt[i ÷ 2] * Mw[i] + βdt[i ÷ 2] * Mw[i + 1] + rhs[i + 1] = @.. fw[i + 1] - βdt[i ÷ 2] * Mw[i] - αdt[i ÷ 2] * Mw[i + 1] + i += 2 + end + + dw = Vector{typeof(u)}(undef, num_stages) + dw[1] = _reshape(LU1 \ _vec(rhs[1]), axes(u)) + for i in 2 :(num_stages + 1) ÷ 2 + tmp = _reshape(LU2[i - 1] \ _vec(@.. rhs[2 * i - 2] + rhs[2 * i - 1] * im), axes(u)) + dw[2 * i - 2] = @.. real(tmp) + dw[2 * i - 1] = @.. imag(tmp) + end + integrator.stats.nsolve +=(num_stages + 1) ÷ 2 + + # compute norm of residuals + iter > 1 && (ndwprev = ndw) + ndw = 0.0 + for i in 1 : num_stages + ndw += internalnorm(calculate_residuals(dw[i], uprev, u, atol, rtol, internalnorm, t), t) + end + + # check divergence (not in initial step) + + if iter > 1 + θ = ndw / ndwprev + (diverge = θ > 1) && (cache.status = Divergence) + (veryslowconvergence = ndw * θ^(maxiters - iter) > κ * (1 - θ)) && + (cache.status = VerySlowConvergence) + if diverge || veryslowconvergence + break + end + end + + for i in 1 : num_stages + w[i] = @.. w[i] - dw[i] + end + + # transform `w` to `z` + #z = T * w + for i in 1:num_stages - 1 + z[i] = @.. zero(u) + for j in 1:num_stages + z[i] = @.. z[i] + T[i,j] * w[j] + end + end + z[num_stages] = @.. T[num_stages, 1] * w[1] + i = 2 + while i < num_stages + z[num_stages] = @.. z[num_stages] + w[i] + i += 2 + end + + + # check stopping criterion + iter > 1 && (η = θ / (1 - θ)) + if η * ndw < κ && (iter > 1 || iszero(ndw) || !iszero(integrator.success_iter)) + # Newton method converges + cache.status = η < alg.fast_convergence_cutoff ? FastConvergence : + Convergence + fail_convergence = false + break + end + end + + if fail_convergence + integrator.force_stepfail = true + integrator.stats.nnonlinconvfail += 1 + return + end + cache.ηold = η + cache.iter = iter + + u = @.. uprev + z[num_stages] + + if adaptive + edt = e ./ dt + tmp = dot(edt, z) + mass_matrix != I && (tmp = mass_matrix * tmp) + utilde = @.. broadcast=false 1 / γ * dt * integrator.fsalfirst+tmp + if alg.smooth_est + utilde = _reshape(LU1 \ _vec(utilde), axes(u)) + integrator.stats.nsolve += 1 + end + atmp = calculate_residuals(utilde, uprev, u, atol, rtol, internalnorm, t) + integrator.EEst = internalnorm(atmp, t) + + if !(integrator.EEst < oneunit(integrator.EEst)) && integrator.iter == 1 || + integrator.u_modified + f0 = f(uprev .+ utilde, p, t) + integrator.stats.nf += 1 + utilde = @.. broadcast=false 1 / γ * dt * f0 + tmp + if alg.smooth_est + utilde = _reshape(LU1 \ _vec(utilde), axes(u)) + integrator.stats.nsolve += 1 + end + atmp = calculate_residuals(utilde, uprev, u, atol, rtol, internalnorm, t) + integrator.EEst = internalnorm(atmp, t) + end + end + + if integrator.EEst <= oneunit(integrator.EEst) + cache.dtprev = dt + if alg.extrapolant != :constant + derivatives = Matrix{typeof(u)}(undef, num_stages, num_stages) + derivatives[1, 1] = @.. z[1] / c[1] + for j in 2 : num_stages + derivatives[1, j] = @.. (z[j - 1] - z[j]) / (c[j - 1] - c[j]) #first derivatives + end + for i in 2 : num_stages + derivatives[i, i] = @.. (derivatives[i - 1, i] - derivatives[i - 1, i - 1]) / c[i] + for j in i+1 : num_stages + derivatives[i, j] = @.. (derivatives[i - 1, j - 1] - derivatives[i - 1, j]) / (c[j - i] - c[j]) #all others + end + end + for i in 1 : num_stages + cache.cont[i] = @.. derivatives[i, num_stages] + end + end + end + + integrator.fsallast = f(u, p, t + dt) + integrator.stats.nf += 1 + integrator.k[1] = integrator.fsalfirst + integrator.k[2] = integrator.fsallast + integrator.u = u + return +end + +@muladd function perform_step!(integrator, cache::AdaptiveRadauCache, repeat_step = false) + @unpack t, dt, uprev, u, f, p, fsallast, fsalfirst = integrator + @unpack T, TI, γ, α, β, c, e, num_stages = cache.tab + @unpack κ, cont, derivatives, z, w, c_prime = 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 + @unpack internalnorm, abstol, reltol, adaptive = integrator.opts + alg = unwrap_alg(integrator, true) + @unpack maxiters = alg + mass_matrix = integrator.f.mass_matrix + + # precalculations + γdt, αdt, βdt = γ / dt, α ./ dt, β ./ dt + (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)) + @inbounds for II in CartesianIndices(J) + W1[II] = -γdt * mass_matrix[Tuple(II)...] + J[II] + for i in 1 :(num_stages - 1) ÷ 2 + W2[i][II] = -(αdt[i] + βdt[i] * im) * mass_matrix[Tuple(II)...] + J[II] + end + end + integrator.stats.nw += 1 + end + + # TODO better initial guess + if integrator.iter == 1 || integrator.u_modified || alg.extrapolant == :constant + cache.dtprev = one(cache.dtprev) + for i in 1 : num_stages + @.. z[i] = map(zero, u) + @.. w[i] = map(zero, u) + @.. cache.cont[i] = map(zero, u) + end + else + c_prime[num_stages] = dt / cache.dtprev + for i in 1 : num_stages - 1 + c_prime[i] = c[i] * c_prime[num_stages] + end + for i in 1 : num_stages # collocation polynomial + @.. z[i] = cont[num_stages] * (c_prime[i] - c[1] + 1) + cont[num_stages - 1] + j = num_stages - 2 + while j > 0 + @.. z[i] *= (c_prime[i] - c[num_stages - j] + 1) + cont[j] + j = j - 1 + end + @.. z[i] *= c_prime[i] + end + #mul!(w, TI, z) + for i in 1:num_stages + @.. w[i] = zero(u) + for j in 1:num_stages + @.. w[i] += TI[i,j] * z[j] + end + end + end + + # Newton iteration + local ndw + η = max(cache.ηold, eps(eltype(integrator.opts.reltol)))^(0.8) + fail_convergence = true + iter = 0 + while iter < maxiters + iter += 1 + integrator.stats.nnonliniter += 1 + + # evaluate function + ks[1] = fsallast + for i in 1 : num_stages + @.. tmp = uprev + z[i] + f(ks[i], tmp, p, t + c[i] * dt) + end + integrator.stats.nf += num_stages + + #mul!(fw, TI, ks) + for i in 1:num_stages + fw[i] = @.. zero(u) + for j in 1:num_stages + fw[i] = @.. fw[i] + TI[i,j] * ks[j] + end + end + + if mass_matrix === I + Mw = w + elseif mass_matrix isa UniformScaling + for i in 1 : num_stages + mul!(z[i], mass_matrix.λ, w[i]) + end + Mw = z + else + for i in 1 : num_stages + mul!(z[i], mass_matrix, w[i]) + end + Mw = z + 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)) + else + linres = dolinsolve(integrator, linsolve1; A = nothing, b = _vec(ubuff), linu = _vec(dw1)) + end + + cache.linsolve1 = linres.cache + + linres2 = Vector{Any}(undef,(num_stages - 1) ÷ 2) + + 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]) + if needfactor + linres2[i] = dolinsolve(integrator, linsolve2[i]; A = W2[i], b = _vec(cubuff[i]), linu = _vec(dw2[i])) + else + linres2[i] = dolinsolve(integrator, linsolve2[i]; A = nothing, b = _vec(cubuff[i]), linu = _vec(dw2[i])) + end + cache.linsolve2[i] = linres2[i].cache + end + + integrator.stats.nsolve += (num_stages + 1) / 2 + + for i in 1 : (num_stages - 1) ÷ 2 + dw[2 * i - 1] = z[2 * i - 1] + dw[2 * i] = z[2 * i] + dw[2 * i - 1] = real(dw2[i]) + dw[2 * i] = imag(dw2[i]) + end + + # compute norm of residuals + iter > 1 && (ndwprev = ndw) + calculate_residuals!(atmp, dw1, uprev, u, atol, rtol, internalnorm, t) + ndw = internalnorm(atmp, t) + for i in 2 : num_stages + calculate_residuals!(atmp, dw[i - 1], uprev, u, atol, rtol, internalnorm, t) + ndw += internalnorm(atmp, t) + end + + # check divergence (not in initial step) + + if iter > 1 + θ = ndw / ndwprev + (diverge = θ > 1) && (cache.status = Divergence) + (veryslowconvergence = ndw * θ^(maxiters - iter) > κ * (1 - θ)) && + (cache.status = VerySlowConvergence) + if diverge || veryslowconvergence + break + end + end + + @.. w[1] = w[1] - dw1 + for i in 2 : num_stages + @.. w[i] = w[i] - dw[i - 1] + end + + # transform `w` to `z` + #mul!(z, T, w) + for i in 1:num_stages - 1 + z[i] = @.. zero(u) + for j in 1:num_stages + z[i] = @.. z[i] + T[i,j] * w[j] + end + end + z[num_stages] = @.. T[num_stages, 1] * w[1] + i = 2 + while i < num_stages + z[num_stages] = @.. z[num_stages] + w[i] + i += 2 + end + + # check stopping criterion + iter > 1 && (η = θ / (1 - θ)) + if η * ndw < κ && (iter > 1 || iszero(ndw) || !iszero(integrator.success_iter)) + # Newton method converges + cache.status = η < alg.fast_convergence_cutoff ? FastConvergence : + Convergence + fail_convergence = false + break + end + end + if fail_convergence + integrator.force_stepfail = true + integrator.stats.nnonlinconvfail += 1 + return + end + + cache.ηold = η + cache.iter = iter + + @.. broadcast=false u=uprev + z[num_stages] + + step_limiter!(u, integrator, p, t + dt) + + if adaptive + utilde = w2 + edt = e./dt + @.. tmp = dot(edt, z) + 1 / γ * dt * fsalfirst + mass_matrix != I && (mul!(w1, mass_matrix, tmp); copyto!(tmp, w1)) + @.. ubuff=integrator.fsalfirst + tmp + + if alg.smooth_est + linres1 = dolinsolve(integrator, linres1.cache; b = _vec(ubuff), + linu = _vec(utilde)) + cache.linsolve1 = linres1.cache + integrator.stats.nsolve += 1 + end + + # RadauIIA5 needs a transformed rtol and atol see + # https://github.com/luchr/ODEInterface.jl/blob/0bd134a5a358c4bc13e0fb6a90e27e4ee79e0115/src/radau5.f#L399-L421 + calculate_residuals!(atmp, utilde, uprev, u, atol, rtol, internalnorm, t) + integrator.EEst = internalnorm(atmp, t) + + if !(integrator.EEst < oneunit(integrator.EEst)) && integrator.iter == 1 || + integrator.u_modified + @.. broadcast=false utilde=uprev + utilde + f(fsallast, utilde, p, t) + integrator.stats.nf += 1 + @.. broadcast=false ubuff=fsallast + tmp + + if alg.smooth_est + linres1 = dolinsolve(integrator, linres1.cache; b = _vec(ubuff), + linu = _vec(utilde)) + cache.linsolve1 = linres1.cache + integrator.stats.nsolve += 1 + end + + calculate_residuals!(atmp, utilde, uprev, u, atol, rtol, internalnorm, t) + integrator.EEst = internalnorm(atmp, t) + end + end + + if integrator.EEst <= oneunit(integrator.EEst) + cache.dtprev = dt + if alg.extrapolant != :constant + @.. derivatives[1, 1] = z[1] / c[1] + for j in 2 : num_stages + @.. derivatives[1, j] = (z[j - 1] - z[j]) / (c[j - 1] - c[j]) #first derivatives + end + for i in 2 : num_stages + @.. derivatives[i, i] = (derivatives[i - 1, i] - derivatives[i - 1, i - 1]) / c[i] + for j in i+1 : num_stages + @.. derivatives[i, j] = (derivatives[i - 1, j - 1] - derivatives[i - 1, j]) / (c[j - i] - c[j]) #all others + end + end + for i in 1 : num_stages + cache.cont[i] = derivatives[i, num_stages] + end + end + end + + f(fsallast, u, p, t + dt) + integrator.stats.nf += 1 + return +end diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl index b2dae1c5fa..750a9a0b34 100644 --- a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl +++ b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl @@ -43,8 +43,8 @@ struct RadauIIA5Tableau{T, T2} T22::T T23::T T31::T - #T32::T = 1 - #T33::T = 0 + #T32::T + #T33::T TI11::T TI12::T TI13::T @@ -56,7 +56,7 @@ struct RadauIIA5Tableau{T, T2} TI33::T c1::T2 c2::T2 - #c3::T2 = 1 + #c3::T2 γ::T α::T β::T @@ -111,6 +111,141 @@ function RadauIIA5Tableau(T, T2) e1, e2, e3) end +struct RadauIIATableau{T1, T2} + T::Matrix{T1} + TI::Matrix{T1} + c::Vector{T2} + γ::T1 + α::Vector{T1} + β::Vector{T1} + e::Vector{T1} + num_stages::Int +end + +function BigRadauIIA5Tableau(T1, T2) + γ = convert(T1, big"3.63783425274449573220841851357777579794593608687391153215117488565841871456727143375130115708511223004183651123208497057248238260532214672028700625775335843") + α = Vector{T1}(undef, 1) + β = Vector{T1}(undef, 1) + α[1] = big"2.68108287362775213389579074321111210102703195656304423392441255717079064271636428312434942145744388497908174438395751471375880869733892663985649687112332242" + β[1] = big"3.05043019924741056942637762478756790444070419917947659226291744751211727051786694870515117615266028855554735929171362769761399150862332538376382934625577549" + + c = Vector{T2}(undef, 3) + c[1] = big"0.155051025721682190180271592529410860803405251934332987156730743274903962254268497346014056689535976518140539877338581087514113454016224265837421604876272084" + c[2] = big"0.644948974278317809819728407470589139196594748065667012843269256725096037745731502653985943310464023481859460122661418912485886545983775734162578395123729143" + c[3] = big"1" + + e = Vector{T1}(undef, 3) + e[1] = big"-0.428298294115368098113417591057340987284723986878723769598436582629558514867595819404541110575367601847354683220540647741034052880983125451920841851066713815" + e[2] = big"0.245039074384916438547779982042524963814308131639809054920681599475247366387530452311400517264433163448821319862243292031101333835037542080594323438762615605" + e[3] = big"-0.0916296098652259003910911359018870983677439465358993542342383692664256286871842771687905889462502859518430848371143253189423645209875225522045944109790661043" + + TI = Matrix{T1}(undef, 3, 3) + TI[1, 1] = big"4.32557989006315535102435095295614882731995158490590784287320458848019483341979047442263696495019938973156007686663488090615420049217658854859024016717169837" + TI[1, 2] = big"0.339199251815809869542824974053410987511771566126056902312311333553438988409693737874718833892037643701271502187763370262948704203562215007824701228014200056" + TI[1, 3] = big"0.541770539935874871186523033492089631898841317849243944095021379289933921771713116368931784890546144473788347538203807242114936998948954098533375649163016612" + TI[2, 1] = big"-4.17871859155190472734646265851205623000038388214686525896709481539843195209360778128456932548583273459040707932166364293012713818843609182148794380267482041" + TI[2, 2] = big"-0.327682820761062387082533272429616234245791838308340887801415258608836530255609335712523838667242449344879454518796849992049787172023800373390124427898159896" + TI[2, 3] = big"0.476623554500550451960069084091012497939942928625055897109833707684876604712862299049343675491204859381277636585708398915065951363736337328178192801074535132" + TI[3, 1] = big"-0.502872634945786875951247343139544292859248429570937886791036339034110181540695221500843782634464164585836226038438397328726973424362168221527501738985822875" + TI[3, 2] = big"2.57192694985560542918678535360167505469448742842178326395573566888176471664393761903447163100353067504020263109067033226021288356347565113471227052083596358" + TI[3, 3] = big"-0.596039204828224924968821911099302403289857517521591823052174732952989090998130905722763344484798508456930766594977798579939415052669401095404149917833710127" + + T = Matrix{T1}(undef, 3, 3) + T[1, 1] = big"0.091232394870892942791548135249436196118684699372210280712184363514099824021240149574725365814781580305065489937969163922775110463056339192206701819661425186" + T[1, 2] = big"-0.141255295020954208427990383807797309409263248498594798844289981408804297900674604638610419147468875667691398225003133444988034605081071965848437945842767211" + T[1 ,3] = big"-0.0300291941051474244918611170890538666683842974606300802563717702200388818691214144173874588956764952224874407424115249418136547481236684478531215095064078994" + T[2, 1] = big"0.241717932707107018957474779310148232884879540532595279746187345714229132659465207414913313803429072060469564350914390845001169448350326344874859416624577348" + T[2, 2] = big"0.204129352293799931995990810298338174086540402523315938937516234649384944528706774788799548853122282827246947911905379230680096946800308693162079538975632443" + T[2, 3] = big"0.382942112757261937795438233599873210357792575012007744255205163027042915338009760005422153613194350161760232119048691964499888989151661861236831969497483828" + T[3, 1] = big"0.966048182615092936190567080794590794996748754810883844283183333914131408744555961195911605614405476210484499875001737558078500322423463946527349731087504518" + T[3, 2] = big"1.0" + T[3, 3] = big"0.0" + RadauIIATableau{T1, T2}(T, TI, + c, γ, α, β, e, 3) +end + +function BigRadauIIA9Tableau(T1, T2) + γ = convert(T1, big"6.28670475172927664517315334186940904959068186655567041187229167532923622489525703260842273089261139845280626287956099768662193453067483410165932355981736786") + α = Vector{T1}(undef, 2) + β = Vector{T1}(undef, 2) + α[1] = big"3.65569432546357225824320796009543385435699888857815445045567025741630720509235614026228963385258117304229337679733945535812317372403535763551850772878775217" + α[2] = big"5.70095329867178941917021536896986162084766017814401034360818390491907468246001534343349900070111312773130349176288004579856585901062722531365183049130382405" + β[1] = big"6.5437368993600772940210715093936863183637851728134458820202187133882261290012752452972782843700946890488789462524897903624959996932392239962196563965573345" + β[2] = big"3.21026560030854988842501065297211721232153653493981008029923647488964744732168461657389754087826565709085773529539707072244537983491480773006949966789260925" + + c = Vector{T2}(undef, 5) + c[1] = big"0.0571041961145176821931211925541156212350779455987501643278082929309346782020731645861138168198427368635148018903413155731609901559772929443100370500757072557" + c[2] = big"0.276843013638123827680045997685625141110889169695030468349442048831121339683708036772541528564051130879197377136636984534220758899839905855114024309075271826" + c[3] = big"0.583590432368916820056697668662917248693432639896771640176293841831747501961831012005632277467456299345321045569611992496682381919275766424103024358378365496" + c[4] = big"0.8602401356562194478479129188751197667383780225872255049242335941839742579301655644134901549264276106897445531811874851737136468026848125542506920602484255" + c[5] = big"1.0" + + e = Vector{T1}(undef, 5) + e[1] = big"-0.239909571163200476817707991076962793618683493830916562279975225042872448414819259070978815977101189851237591634144816820425592764432710089981892023441549743" + e[2] = big"0.125293484229223300606887443525929336197638450194973243323835481816625682669684634271160851996902445722310139152840852195000603355301430153561918341655137084" + e[3] = big"-0.0688288849083782089370741375422279772873399871312158026536514369967825732836040693366396751988019623495452730460577336089848105495733304937016477629621990433" + e[4] = big"0.0372433600301293198267284480468961585078575499935477902539140092558239369583143189611737274891328175087350382605569395937945872776839987254373869550943049195" + e[5] = big"-0.012863950751139890646895902730137465239579845437088427263654957605488133042638147254426913683751171160691603649073170415735165315443538347036196755548109703" + + TI = Matrix{T1}(undef, 5, 5) + TI[1, 1] = big"30.0415677215444016277146611632467970747634862837368422955138470463852339244593400023985957753164599415374157317627305099177616927640413043608408838747985125" + TI[1, 2] = big"13.8651078562714131651762946846279728486098595017962436746405940971751244384714668104145151259298432908422191238542910724677205181071665482818120092330632702" + TI[1 ,3] = big"3.48000277479518556182840016971955819123081637245954095062693470191383865922357339844125383481645392882289968250993872221445874555610460465838129969397069557" + TI[1, 4] = big"-1.03200879782526342277108071214631493513824682491749273908106331923801396656058254294323988505859654767877050109789490714699847664805679842903430004696170252" + TI[1, 5] = big"0.804303045073989917475330383606196086089578671788707543063308602519859970319818304759856653218877415405946945572102875643297890954688508528143272905631829894" + TI[2, 1] = big"5.34418643783491159889531030409736033885455686563071401172022718575590068536629704134603404624953791012861634674294690788961703408019660066685859393456498931" + TI[2, 2] = big"4.59361556775916100445407449817656238428260055301676371438973411021009514435572975394999086474831271997070798032181411537895658457000537727156665947774751386" + TI[2, 3] = big"-3.03636032345942429864615756872018980250277648141683630832856906288036929718223473102394179699607901856890769270810252103326382063852039607285826867723587514" + TI[2, 4] = big"1.05066019023145886385983615715299311307615150447133905233370933194949591737765763708886464382722316727972166443876395823044171403663404254906698768838255919" + TI[2, 5] = big"-0.272778611864296270538614649997366804891835224042737605275699398413256470423268908248569612750117948720141667949532252500428432062582365619208502333677907158" + TI[3, 1] = big"3.74805980743980486005103450189256983678052751095791526209741655305580351377124372457009580386663275146166007984852101733055495783906881063060757645038080343" + TI[3, 2] = big"-3.98496573634388466725226385805351110838575115293851360514636734529255361185420464416807882769853298186283398369873418552760618971047757002216338511286260041" + TI[3, 3] = big"-1.04441564160801879294224732309562532189841624726401645191058551173485917137499204844819781779667611903670073971659834929382224472890100209497741235960707456" + TI[3, 4] = big"1.18409856813794848723102038838340482030291345603197522521517834943166421242518751666675199211369552058487095283489346390066317584532997854692445653563909898" + TI[3, 5] = big"-0.449917770156780368898811918314095435942113881883174152777026977062686286863549565130412864190301081537983106397709991028107600781961279985605930655683680139" + TI[4, 1] = big"-33.0418802135190000080614469426109507742858088371383868670878639187564531424382858814386742148456699143328462132296293097447566408853495288807407929988004676" + TI[4, 2] = big"-17.3769534790635670194549806058987105852733409102703844354448800193942184746909147697382687117638715195698950138089979798321855885541817752366521518811413713" + TI[4, 3] = big"-0.172129063254005561151528806427751383749451500597823574207174433146207178559871803504021077429693091164540897873472803934375603405253541639437370184767553293" + TI[4, 4] = big"-0.0991697779825426425881662214017368584726354746776989845479783944003623924121748016326495070834800297497011104846871751430208559227945252758721362340763610828" + TI[4, 5] = big"0.531228115838306667184911422606024795426589562580669892779793097035561488973256023529352389498509937781553683467106048413485632583844632286562240161995145055" + TI[5, 1] = big"-8.61144397987529197770008251257034851950485933115010902789613925540488896812417081206983938638600226846804467531843522104806738090683710882069500386691775154" + TI[5, 2] = big"9.69999140952880823133589405342003266497120753048627084327055311528684684237122654108691149692242002085965723391934376924400492239317026460192827344970015484" + TI[5, 3] = big"1.91472863969687428485137560339172471528025297511003983469957355306260543484472462223194401768126877615795915146192537091374017807611943419264038682143890747" + TI[5, 4] = big"2.41869200608494002642656343408298350771199306961305597858229870375990977712805399625496435641846363295393762353024017195444763964531237381728801981679934304" + TI[5, 5] = big"-1.0474634879353374186944329992117360176590042540536055452919974336199826846201614544718272622833822842591012529895091659029452542118642301415759073410771819" + + T = Matrix{T1}(undef, 5, 5) + T[1, 1] = big"0.0125175862205010458901356760368001462557655123420858705973577952199246108029451084239310924615007306721702298573083400752464277227557045438770401832498107968" + T[1, 2] = big"-0.0102420478179088270700863300668590125015813934827825923708366359399562125950804289592272678367034071306578383319296130180550178248531589487456925441921649293" + T[1 ,3] = big"0.0476738772902957238631839478592069782970238490568258436986723993118380988311441474394156362952631834786373081794857384127209450988829840886524135970873769918" + T[1, 4] = big"-0.0114785152552295147079415554121555049385506204591245712490409384029671974157542450636658532835395855844059342442518520033304129991000509527123870917346017759" + T[1, 5] = big"-0.0140198588928754102810778942934959307831026572823203692568448424056201483917805257790275956734469193171917730378117501915144713896813544630288006687542182225" + T[2, 1] = big"0.00149167015189538242900444775236282223594625052328927847572623038484966999313257893341818287477809424303168766872838075463220122499449382436194198620498144296" + T[2, 2] = big"0.050172864517371058162991380262646513853120568882725793734131676894272706020317186004736779675826101816279321643304301437029912742375638648226701787880031719" + T[2, 3] = big"-0.0943318191816114369806569003363724471884924328367212069321438749304281980331334016578193750445513659941246363262225907407726099492713722343006925656625258579" + T[2, 4] = big"-0.00766883074918016288515687679203608074116106558796378201472238095295554979920808799930579174190884587422912077296093093698836937450535804218413704866981728518" + T[2, 5] = big"0.024708578426518526812525205377780382655366504554979744093019395818934704623702078004474076773426928900579988063099593288435684744957695210778788200213260272" + T[3, 1] = big"0.072981876388087148622657299703669587832652508881663282287850495621401398441897288250625556038835308015912409648841893161563884759791665776933761278383553608" + T[3, 2] = big"-0.230539534043417946721421862180000422679228296568599014834226319726930529322581417981617275287468418138394077987361681288909676234537699721082090802790143303" + T[3, 3] = big"0.102703045380125899792210456947141185148813233939327773583525878521508211077874610560448598369259541346968946573971195783374996178436435357335759255990489434" + T[3, 4] = big"0.0193984639988289509112232896408330872285824216708905773930244363652651247181543158008567311548336143384128605013911312875018664026371225431993252265128272262" + T[3, 5] = big"0.0818003537037511708363908122287572533071340646031113975848869261019231448226334426630664318901554550460201409321555775999869184033436795623062614812355590017" + T[4, 1] = big"0.380091440003568104126439184355215575526619121262253024859378518379910007234696730891540745160675744992320824590679292148769326540463161583672773762554445506" + T[4, 2] = big"0.377893902248861249543862293745933995234687511602719536459666284734445918178134851270924212812363352965391508894581698067329905034837778770261095647458874628" + T[4, 3] = big"0.466744130332494359289559582964906703283968612669234331018678042733321473730897217606173184300477207393539851157929838664168404778962779344509707214938022808" + T[4, 4] = big"0.40760117128019906662166237021895987274626181127101561893104166874567447589187790736078997321464949349935802836110699884016973990503134772720646054039223561" + T[4, 5] = big"0.199682427886802525936540566022390695167018315867216115995143539347975271751460199398235415129329119718414206048034051939441434136353381864781262773401023899" + T[5, 1] = big"0.921978973681210488488254647415676321266345412943047462855852351388222898143904205962703147998267738964059170225806964893009202287585991334322032058414768529" + T[5, 2] = big"1.0" + T[5, 3] = big"0.0" + T[5, 4] = big"1.0" + T[5, 5] = big"0.0" + + RadauIIATableau{T1, T2}(T, TI, + c, γ, α, β, e, 5) +end + + struct RadauIIA9Tableau{T, T2} T11::T T12::T @@ -258,3 +393,241 @@ function RadauIIA9Tableau(T, T2) γ, α1, β1, α2, β2, e1, e2, e3, e4, e5) end + +function BigRadauIIA13Tableau(T1, T2) + γ = convert(T1, big"8.93683278840521633730209691330107970355008194433956657198414191417624969654351559268800871286734194720118970058657997472527299153742511021973612156231867783") + α = Vector{T1}(undef, 3) + β = Vector{T1}(undef, 3) + α[1] = big"4.37869356150680600252334919268856129165763746518197948235657247177701087073069907016715245914093899486193202405685779803686971216417800783050995450529391908" + α[2] = big"7.14105521918764010577498142571556804318193862372238812855726792587872300446315860222917039505087745633962330233504078264632719519730762016919715839787116038" + α[3] = big"8.51183482510294572305062092494533081338538293892584910309408864525614127653438453125967278937451257519784982331481143195416659686980181689042482631568989031" + β[1] = big"10.1696932837950116273183544188477298930096536824510223588525334625762336174947183926243705927725260475934351162622185429326813205432867247703480391692806137" + β[2] = big"6.62304592263927597062055811591186110468148199066707542227575094761515104946479159063603447729283770429494038962408904312215452856333028405675512985803584472" + β[3] = big"3.2810136243250588300359425270393915846791621918405321383787427650552081712406957205287551182809705166989352673500472974040971593568323836675590314648604458" + + c = Vector{T2}(undef, 7) + c[1] = big"0.0293164271597848919720502769131649103737303925637149277869106839449360382416657787486309483651843695097273923248526200112627747993405898353736305552306269904" + c[2] = big"0.148078599668484291849976852495979212230248774808594461412594641801598386090878321806369397661747576057906341132861865305306667654594593138746653233717241913" + c[3] = big"0.336984690281154299097052972080775705197568750028473347122562968073691350512784060852409141173654482529393236826516171319486086447256539582972346127980810124" + c[4] = big"0.558671518771550132081393341805521940074368288965407825555747226117350122897421078323820052012282581935200398463518265914564420109615277886000739200777932339" + c[5] = big"0.769233862030054500916883360115645451837142143322295416166948169636548130573953285685200211542774367652885154701431860087378103033801830280742146083476036669" + c[6] = big"0.926945671319741114851873965819682011056172419542283252724467079656645202452528243814339480013587391545656707320049986592771178724621938506933715568048004783" + c[7] = big"1.0" + + e = Vector{T1}(undef, 7) + e[1] = big"-0.171003707892600662399316289094713451418682228802554671094543075316220821099884263713032871578607576486535539488632407200766379971279557791263375813632287147" + e[2] = big"0.0934967172358652400317534533028674569308657324394316331629203486361371292312231403973668315582870547753526899857449840409175610009437530537219068836721686211" + e[3] = big"-0.0538908303114758775848180855518003793385120454908028879947132475960997563222416509683480350817114056356343433378213334826358824351525243402758389616051172681" + e[4] = big"0.03036786965048439581219923157368590250090409822952169396779792168990510618756404452728392288892998298088147691907128776816545685599760715439221674418662785" + e[5] = big"-0.0169792974425458224044481617230998766694942329759644144506911809530904808476835739189122151558601434810772520378036293579816345384682687602578758514350075723" + e[6] = big"0.00942688256820236884916415666439281573527695349338346787075550606528435808025071461023926432308932314282041885090975780812273515256740094297311541275151861292" + e[7] = big"-0.00331409873565629283448601269346047459594635696619041493081994712789731442072563377354737487903843138987115421498455722938021358621090485566426506726181005806" + + TI = Matrix{T1}(undef, 7, 7) + TI[1, 1] = big"258.131926319982229276108947425184471333411128774462923076434633414645220927977539758484670571338176678808837829326061674950321562391576244286310404028770676" + TI[1, 2] = big"189.073763081398508951976143411165126555759459745371576264125287430947886413126866952443113984840310549596923934762141954737541643761162558070450614795561734" + TI[1, 3] = big"49.0873148179301311944474703372633419330229683717897887664283914712555334645741343066714059043135343948204451450061803442374878045458955826422757210762412997" + TI[1, 4] = big"4.11064746966142841811238518636124668078589358089581133578005291508858571621836624121708112101643343488669794287298806656198949715476379639435093560435010553" + TI[1, 5] = big"4.05344788931556330417512803837862541661144275947069236866476426664242632965376171604053865483440478823853326237912519148507906655855071507442222711969825069" + TI[1, 6] = big"-3.11275536660734607655357698925636361735741304308245452106573904595716690770542970584435712650159533448326091358879097717388530116398450168049097806992817596" + TI[1, 7] = big"1.64677491355844465016894934800942442334612077828885771793164268655566366462165061862443368822544695623147966149765223644798045399342853834086413561960176148" + TI[2, 1] = big"-3.00739016945129213173149353792169083141834116044470099212013728771587881480191343754504173052952073006187734389002396348355357273701343509199048972794392147" + TI[2, 2] = big"-11.0158660787657713291120393664792067595453921824881213620299497076376976067619617086470844707815815293102862568459526162951253770377715406520772358338647188" + TI[2, 3] = big"1.48779945613165628148618248664965038886474377325027865838645297753993182317594482435706956176392903188004580583104018591540474622009639200188521283880201225" + TI[2, 4] = big"2.13038815955928245943197208332824475219642634294808813866153957342980992047877237670079423767538654092424134276380826377135080667266661637001176204430488753" + TI[2, 5] = big"-1.81614108681756562482220455159496741723359999245934818387747079566312917815672128128449281415737713177900591942282975861961228230314168417307836619006791605" + TI[2, 6] = big"1.13432558789516110008277908420532415765361628740656810686297793967986689714948610119162966211301325316623863222505219543867472186257492829970663316956377323" + TI[2, 7] = big"-0.414699045943303531993049422295928526684402022493736427543557958358387925728160703636844863663828153394608981043415378230601486738224597324364079320598162815" + TI[3, 1] = big"-8.44196318832108468175691559413731210343158392484322786670758421404507417209484447031645790366021837365786640573614305718894911853549168061902141351516580451" + TI[3, 2] = big"-0.650525274057515002816904045893485631294530894981669254094573985727348985809697093879080285963063573837365484483755274668080611163704039179328960851461387071" + TI[3, 3] = big"6.94067073036987647880408175445008301222030789462375109942012235845495260572570799226646472429196555932436186979400567616504159564738984233922289782922787445" + TI[3, 4] = big"-3.20504752559789843156502799159713971965747774043426947358779973217345866996463287674334224123932879873323284636947452187683408110992957222808611161423213549" + TI[3, 5] = big"1.07128094354647858978279562700457911254627057919002861801894953308482120936700881726232902304000322718645130593907512149815870969208873216470962770569998532" + TI[3, 6] = big"-0.354850749121622187972972761073874956531274189535504546398851680169235702590362534883357256681588685608802983372517893712333972644320006895019178184808028042" + TI[3, 7] = big"0.0919854913278655415440864884207305663999562250023079120516746551750254082665966708567906888946992351083964961208132558221142585217674963218388224937302473142" + TI[4, 1] = big"74.6783322350226997715286176267232500441551583987525066913719852490109364599462546293112601362342028584101507709386240000804692470037564789980905370400509214" + TI[4, 2] = big"87.4085889799008164020396362924136436577534600993283836959398121813667403209890699914314446222016952621954817633686823685774595935180374571416781238038364186" + TI[4, 3] = big"4.02415873737999787701407840793921059156554118449220356776918749072220128918152906578385457943212213189933447495921754693186811343717296680238755923076427455" + TI[4, 4] = big"-3.7148063151583641866387382381081795406061842159003055897302686185198568522128509989890869602984467843559169959313018612449354703104270603001605170037725663" + TI[4, 5] = big"-3.43009398598231735074090769130593476067104938465255451803266927011738721835297930406017172365070584279715308905584391225176154776278518922912169890517961929" + TI[4, 6] = big"2.69660480976531237885262500230842013033719691844775548640355919138284680959979836353143310081338215041119022648809147361433752919265159399610746756470853959" + TI[4, 7] = big"-0.938692743607546193356785681771531136814109179879957291315724533839534255667763099330792864148293396694586387338161584706252944483821135344465739888811338788" + TI[5, 1] = big"58.3565288519065772423731088606544342599129168115273649928818622008651860145833895668543250775742899696760389837877193028417145182338484929599333810581515993" + TI[5, 2] = big"-10.0687739578001809632495544545749228539542767485211306078205622876595603032162891608453826862136355989387474454697691529766293644115682409173741730758425432" + TI[5, 3] = big"-30.3663888425666712081087189214021522992426235463582449811325590575576319489955157279473313224901192335775884848736150180108985558310423628914140477437063457" + TI[5, 4] = big"-1.02002086518486598502718784312141857841892430616701325398305811243769008274372077411348691412296276168896198187688441456921700292037247387330560786140723416" + TI[5, 5] = big"-0.112417500378424962126670249921897816128157398591725875330925039631874967429838848482089690872916638698820411392685501889126627650123714184027159547685248056" + TI[5, 6] = big"1.89064083100037762279966919417932484200269828564004442737723486475878958135985745266991261770924069476112679285337233931312540904735632744873728510014970829" + TI[5, 7] = big"-0.971648639383148228217233127548943147296423534674266405843322723719694664032217172325052282800290275002731997713145411340983758516166807609661717915219518127" + TI[6, 1] = big"-299.18624802825209667863642523944728107942141534516550178278869311293354511449399684666660494133688445719285752471650937062695632169114367079856135650539072" + TI[6, 2] = big"-243.040745368744791181900565230083092669143049316165122405971394775932180012728275256467636352341415340547177922968547123544546515287229215470481168446631934" + TI[6, 3] = big"-48.7771040780378692121909344887388032694629956594617430615510915251995189158287187599892740037773277403958100797917560590738598108409472582147091119440886778" + TI[6, 4] = big"-2.03867190574193440528015205293433905622043272233073734690244789947707827347049413187234402189062846366658963666461334786306660732097114011309282331323116958" + TI[6, 5] = big"1.67356023986108494426829042309213202110891938292923077616474877079402040904687073610625868939896244842053999572446723558562427506280564629528151134946587118" + TI[6, 6] = big"-1.0873740320571061644555969255032311107358443063278089996181949045168433801494845898897631535619158410753032807069032950523487601457868753453652745002841107" + TI[6, 7] = big"0.901938249296099373842715514839004052963355800714627971724094542443991299921284427589690820402982448873149676210397055957126153220340909284180014056386791594" + TI[7, 1] = big"-93.076502897435305911571945263737383854569504715670989865831914555937966339933932282945955570244055882294556430466422133231853008314991630740535709028417842" + TI[7, 2] = big"23.8816310562811442770319002318043863376962876994405756649585750650966186536576789769674007990310112890015051984278059899811178135726914390958188405071290871" + TI[7, 3] = big"39.2788807308138438271015646136760366834412493325456249795727722130258444051594274416196392795817449902122139076648927894476044063388859377757097127385794539" + TI[7, 4] = big"14.3889156854910800698761307424979534708984169042483973564042387223013868069040933228077604321320066763752720714195604903398768371784013771964086553618150626" + TI[7, 5] = big"-3.51043839939936122108708432480845734972162782563284715495715984978907792386567906732993553255070093796782368160341757151292477304975079070782335737053297468" + TI[7, 6] = big"4.86328488556618070121491058699734313503568312572977577331134555924656926935558698308076704662503608259898740028814153544991114426972747448736702277116049277" + TI[7, 7] = big"-2.24648272959123991640046924839711232278867381637608763335081676684616443569602032178385937243819174902544136208243971053224668691848283004752869023074006745" + + T = Matrix{T1}(undef, 7, 7) + T[1, 1] = big"0.00215375462731052642282751906550204337272018200721827917615061640312650856312529840445028048591986867096756005142895325420603307041594804305862850861253757163" + T[1, 2] = big"0.021567551351320773386914226953811992365459277376204369162736830595700124529879508417849062386878143122032508776691627063229415272329484156789207145821702462" + T[1, 3] = big"0.00878356792514414440732555660043326940873333657406338685620618347939710728032290406426688328221296324998146697730909767495361893387567339044816921837538988154" + T[1, 4] = big"-0.00405516145233102389819844704090310382485225922827010954643577855973533421255114497764957587851178840064428149215351434824919490696577563849929483184955933965" + T[1, 5] = big"0.00442723275326828547967807873499027629097834766201549949492135358632150336069311115075327876323707841703727317338755331613570950287342825020738596326021052902" + T[1, 6] = big"-0.00123864618795287405637686870391105285581324510790128485733529975336279476721707053186563729417080236061385260749762448518679294700311105630290083016823761156" + T[1, 7] = big"-0.00276061748054385249954800379096675592021481213358861974911688001011761550911589157738523818859000828996335817774948428177282421412491830529445501318154035024" + T[2, 1] = big"-0.00160002507788042852683067347985080829550105638728462477214069614397009338180775134535418790113854904464693278677067195562013777079470430165035085043732753352" + T[2, 2] = big"-0.0381316481344115466944201512445271892551007922443248010648630183723114657457789198582213862424187595732944781586531399310738197517976083499508550510483478779" + T[2, 3] = big"-0.0215255605940068755238494349163503963236812065771639056145559371805737876208350036328339608215271680572576146954552666030277743869132676140541472724370558091" + T[2, 4] = big"0.00841556827655958923717700333156546206587781542530241328710392714333753219743181540077241302321588065650704924760060316717877095134935044662592211744890794666" + T[2, 5] = big"-0.00403194957022454949230429372587008587329606687054571010486662485715979240183165499902791387008699068626978608835015342675934092134962673636484308565473356683" + T[2, 6] = big"-6.6666353393963381817604789740257628821376819567901071737415235834331307484818353061850936507762955342131861918219584166678095273744210157164382779907235669e-05" + T[2, 7] = big"0.00318547482516620984874835878222687621122035448401205459368674257818574765593899794870819769668503869906022860261901897250913569265553156976061140932045107432" + T[3, 1] = big"0.00405910730194768309165024146216588597640781263680870767202041411242133338742562561902630276038676420444232405079851555753917806998064489819308813790494788924" + T[3, 2] = big"0.0573965089393817153975680203880753938458832782600090443030839643350468249623833638779578474891654213594195393636829414422184571666256857425091138479371917574" + T[3, 3] = big"0.0588505292084267910561208969865829735901655409220388105109199298038946675765714122525765330769443473927581930134049676200572930797370286476504623214740871248" + T[3, 4] = big"-0.00856043106160343206017727185390754992573940897343949944649743606465705403614377469754987858631901604547097801042861815249197647886051332362774581709381720893" + T[3, 5] = big"-0.00692321266502390892414068519049460069371592099748070119636478595631451405094203293036429762819458535062492059219566837532157551782305886338773933077463475632" + T[3, 6] = big"-0.00235218098294333834053519532555529491776729377182703234025085030409255592197086839142988525473684138901264206886166295186155491132922909402254443843846019141" + T[3, 7] = big"0.00041690777252975626914088803059940941342549922756308931704215701350026719541939053570614368159222367707113801117750298289694571643601584878405615892432648487" + T[4, 1] = big"0.0157504880793768442034586734054915501004520506405808322686493022779655453114657621318660532381583918124125360276320121127974912393389579826125529804830864399" + T[4, 2] = big"-0.0382146935969683504846411337659300127514788882892071252172987515109399372135899067290947441850340146027892665775682097051548343529370733593281856326317259999" + T[4, 3] = big"-0.165736811272943851241241116255535218556011122333381899790277357803281567727036568454939356458468926429537927937619042817050400333625919290585510785057955509" + T[4, 4] = big"-0.0373712423023844574190702119163246888117181457309185176497005310822879226235861373253125139016964433591381638592353617347369492240160809914228784174846477722" + T[4, 5] = big"0.00823900729850771940449868235563938395546999707236910359131464615707125576979409087864780171789078059526539789661318173387826643385244974406562622466790754233" + T[4, 6] = big"0.00311507115234617525272547086289315208054441921705361129575617631104650731644437585122142710666234276633544335552925569262424677362146587776195531866754755781" + T[4, 7] = big"0.025116604913438821928363823471446698278976101918753236732238210724710282378748917637317846485853317873304329580245705683618093593158791190832004186288367408" + T[5, 1] = big"0.112977661024220807608615842313106352633973778091080400075534257952348289641328709240673869677499013004285003126194992176632265223545565047727637631580337111" + T[5, 2] = big"-0.249174212465263686330825594009221950347570740813751325091913985975498424569678307894304962660904874986611526140914403971840496728150916599999921976188547708" + T[5, 3] = big"0.273563305798662321213236935135336593478278696397012151365678540099566245199777083242808233574654642014215983653810819494932091426330017240672955510133726276" + T[5, 4] = big"0.00536676137918177009427930181087914853701809128264121101773394730339300080525157052081366996826642003169044168721911822166683675089051631342776752635189343996" + T[5, 5] = big"0.193211116101262014431211225620266980060733605289133050251158448403922545905872373640500736693735926480983370235582910255756813799388364741420161359961401418" + T[5, 6] = big"0.101717732481715146808078931323995112561027763392448195424858681165964478003318758266672250034474900552688318026734856778296896546916272032434282368222825518" + T[5, 7] = big"0.0950450203560462282103892144485647895183175432965514336285840628832838918715022627077373617151475963061484489345238022187829573892306346658797861719620799413" + T[6, 1] = big"0.458381043183931501028085939964292092908293295595258886425372669820276128937720150467378912424378376379185138190017965370589550781979145790869568608776861466" + T[6, 2] = big"0.5315846490836284292050500994300107341125728347976407285397462896004659632807779347307732180848765709277026749725126234633983063167374333425454720010026876" + T[6, 3] = big"0.486322836617572894056685295353340203321316764127126557475136642083389075853199222650975554544550110757249234979120491845825690852575400863926535437662617201" + T[6, 4] = big"0.526574226458449262914091192639271913456008564881594253716678163127743947224108435833618497118891017505982561930788522171455486058320589875335702474378251931" + T[6, 5] = big"0.275534394989625814192875938762525038291639319966986287664787801569471609648366101593885546008609962622035890891754680149203464179471952105174480329668882489" + T[6, 6] = big"0.521751945274765285294609453181807034209434470364856664246194441011327338299794536726049398636575212016960129143954076748520870645966241492966592488607495009" + T[6, 7] = big"0.128071944635543894414114939510913357662538610722706228789484435811417614332529416514635125851744500940930818246509599119254761178392202724896572159336577251" + T[7, 1] = big"0.881391578353818376313498879127399181693003124999819194603124949551827789004545406999549226388170693806014968936224161749923163222614460424501073405017519348" + T[7, 2] = big"1.0" + T[7, 3] = big"0.0" + T[7, 4] = big"1.0" + T[7, 5] = big"0.0" + T[7, 6] = big"1.0" + T[7, 7] = big"0.0" + + RadauIIATableau{T1, T2}(T, TI, + c, γ, α, β, e, 7) +end + +using Polynomials, LinearAlgebra, GenericSchur, RootedTrees, Symbolics +using Symbolics: variables, variable, unwrap + +function adaptiveRadauTableau(T1, T2, num_stages::Int) + tmp = Vector{BigFloat}(undef, num_stages - 1) + for i in 1:(num_stages - 1) + tmp[i] = 0 + end + tmp2 = Vector{BigFloat}(undef, num_stages + 1) + for i in 1:(num_stages + 1) + tmp2[i]=(-1)^(num_stages + 1 - i) * binomial(num_stages , num_stages + 1 - i) + end + radau_p = Polynomial{BigFloat}([tmp; tmp2]) + for i in 1:(num_stages - 1) + radau_p = derivative(radau_p) + end + c = real(roots(radau_p)) + c[num_stages] = 1 + c_powers = Matrix{BigFloat}(undef, num_stages, num_stages) + for i in 1 : num_stages + for j in 1 : num_stages + c_powers[i,j] = c[i]^(j - 1) + end + end + inverse_c_powers = inv(c_powers) + c_q = Matrix{BigFloat}(undef, num_stages, num_stages) + for i in 1 : num_stages + for j in 1 : num_stages + c_q[i,j] = c[i]^(j) / j + end + end + a = c_q * inverse_c_powers + a_inverse = inv(a) + b = Vector{BigFloat}(undef, num_stages) + for i in 1 : num_stages + b[i] = a[num_stages, i] + end + vals = eigvals(a_inverse) + γ = real(vals[num_stages]) + α = Vector{BigFloat}(undef, floor(Int, num_stages/2)) + β = Vector{BigFloat}(undef, floor(Int, num_stages/2)) + index = 1 + i = 1 + while i <= (num_stages - 1) + α[index] = real(vals[i]) + β[index] = imag(vals[i + 1]) + index = index + 1 + i = i + 2 + end + eigvec = eigvecs(a) + vecs = Vector{Vector{BigFloat}}(undef, num_stages) + i = 1 + index = 2 + while i < num_stages + vecs[index] = real(eigvec[:, i] ./ eigvec[num_stages, i]) + vecs[index + 1] = -imag(eigvec[:, i] ./ eigvec[num_stages, i]) + index += 2 + i += 2 + end + vecs[1] = real(eigvec[:, num_stages]) + tmp3 = vcat(vecs) + T = Matrix{BigFloat}(undef, num_stages, num_stages) + for j in 1 : num_stages + for i in 1 : num_stages + T[i, j] = tmp3[j][i] + end + end + TI = inv(T) + + if (num_stages == 9) + e = Vector{BigFloat}(undef, 9) + e[1] = big"-0.133101731359431287515066981129913748644705107621439651956220186897253838380345034218235538734967567153163030284540660584311040323114847240173627907922903296" + e[2] = big"0.0754476228408557299650196603226967248368445025181771896522057250989188754588885465998346476425502117889420021664297319179240040109156780754680742172762707621" + e[3] = big"-0.0458369394236156144604575482137179697005739995740615341890112217655441769701945378217626766299683076189687755618065050383493055018324395934911567207485032988" + e[4] = big"0.0271430329153098694457979735602502142083095152399102869109830450899844979409229538982100527256348792152825816553434603418662939944133319974874915933773657075" + e[5] = big"-0.0156126300301219212217568535995825232086423550686814635293876744035364259647929167763641353639085929285192729729570945658304937255929114458885296622493040224" + e[6] = big"0.00890598154557403928205152521539967562877335780940124672915181111908317890891659158654221736499522823959933517986673010006749138291836676520080172845444352328" + e[7] = big"-0.00514824122639241252178399021479378841872099572255461304439292434131750195489022869965968028106854978547414579491205935930595041763060069987112580994637398395" + e[8] = big"0.00296533914055503317169967748114188676589522458557982039693426239853498956125735811263087631479968309978854200615027412311940897061471388689986239742919640848" + e[9] = big"-0.0010634368308888065260482548541946175520274736959410047497431569257848032902381738362547705844630238841535652230832162703806430112125115777122361837311714267" + else + p = num_stages + eb = variables(:b, 1:num_stages + 1) + @variables y + zz = zeros(size(a, 1) + 1) + zz2 = zeros(size(a, 1)) + eA = [zz' + zz2 a] + ec = [0; c] + constraints = map(Iterators.flatten(RootedTreeIterator(i) for i in 1:2*p-3)) do t + residual_order_condition(t, RungeKuttaMethod(eA, eb, ec)) + end + AA, bb, islinear = Symbolics.linear_expansion(Symbolics.substitute.(constraints, (eb[1]=>1/γ,)), eb[2:end]) + AA = Float64.(map(unwrap, AA)) + idxs = qr(AA', ColumnNorm()).p[1:num_stages] + @assert rank(AA[idxs, :]) == num_stages + @assert islinear + b_hat = Symbolics.expand.((AA \ -bb)) + e = [Symbolics.symbolic_to_float(b_hat[i] - b[i]) for i in 1 : num_stages] + end + RadauIIATableau{T1, T2}(T, TI, c, γ, α, β, e, num_stages) +end diff --git a/lib/OrdinaryDiffEqFIRK/src/integrator_interface.jl b/lib/OrdinaryDiffEqFIRK/src/integrator_interface.jl index 54178b9923..dc4bbb1e88 100644 --- a/lib/OrdinaryDiffEqFIRK/src/integrator_interface.jl +++ b/lib/OrdinaryDiffEqFIRK/src/integrator_interface.jl @@ -1,5 +1,4 @@ -@inline function DiffEqBase.get_tmp_cache( - integrator, alg::Union{RadauIIA3, RadauIIA5, RadauIIA9}, +@inline function DiffEqBase.get_tmp_cache(integrator, alg::Union{RadauIIA3, RadauIIA5, RadauIIA9, AdaptiveRadau}, cache::OrdinaryDiffEqMutableCache) (cache.tmp, cache.atmp) end diff --git a/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl b/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl index b2da07a091..2eb857827d 100644 --- a/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl +++ b/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl @@ -9,10 +9,19 @@ for prob in [prob_ode_linear, prob_ode_2Dlinear] end sim21 = test_convergence(1 ./ 2 .^ (2.5:-1:0.5), prob_ode_linear, RadauIIA9()) -@test sim21.𝒪est[:final]≈8 atol=testTol +@test sim21.𝒪est[:final]≈8.5 atol=testTol sim21 = test_convergence(1 ./ 2 .^ (2.5:-1:0.5), prob_ode_2Dlinear, RadauIIA9()) -@test sim21.𝒪est[:final]≈8 atol=testTol +@test sim21.𝒪est[:final]≈8.5 atol=testTol + +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, 9], prob in [prob_ode_linear_big, prob_ode_2Dlinear_big] + dts = 1 ./ 2 .^ (4.25:-1:0.25) + sim21 = test_convergence(dts, prob, AdaptiveRadau(num_stages = i)) + @test sim21.𝒪est[:final]≈ (2 * i - 1) atol=testTol +end # test adaptivity for iip in (true, false)