diff --git a/lib/OrdinaryDiffEqFIRK/Project.toml b/lib/OrdinaryDiffEqFIRK/Project.toml index 308a831ced..ef8f74b22f 100644 --- a/lib/OrdinaryDiffEqFIRK/Project.toml +++ b/lib/OrdinaryDiffEqFIRK/Project.toml @@ -6,12 +6,15 @@ version = "1.1.0" [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" SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961" diff --git a/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl b/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl index f7039b2cfc..a37ff9383b 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 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..f8ed7b9b63 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) = 9 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 d95430f5a5..4429fb78b6 100644 --- a/lib/OrdinaryDiffEqFIRK/src/algorithms.jl +++ b/lib/OrdinaryDiffEqFIRK/src/algorithms.jl @@ -153,3 +153,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 516775e085..e17f486570 100644 --- a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl +++ b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl @@ -470,3 +470,164 @@ 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::AbstractVector{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 + tab = adaptiveRadauTableau(uToltype, constvalue(tTypeNoUnits), num_stages) + + cont = Vector{typeof(u)}(undef, num_stages - 1) + for i in 1: (num_stages - 1) + 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, uNoUnitsType, rateType, JType, W1Type, W2Type, + UF, JC, F1, F2, Tab, Tol, Dt, rTol, aTol, StepLimiter} <: + OrdinaryDiffEqMutableCache + u::uType + uprev::uType + z::AbstractVector{uType} + w::AbstractVector{uType} + dw1::uType + ubuff::uType + dw2::AbstractVector{cuType} + cubuff::AbstractVector{cuType} + cont::AbstractVector{uType} + du1::rateType + fsalfirst::rateType + k::AbstractVector{rateType} + fw::AbstractVector{rateType} + J::JType + W1::W1Type #real + W2::AbstractVector{W2Type} #complex + uf::UF + tab::Tab + κ::Tol + ηold::Tol + iter::Int + tmp::AbstractVector{uType} + atmp::uNoUnitsType + jac_config::JC + linsolve1::F1 #real + linsolve2::AbstractVector{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) + alg.num_stages = num_stages + tab = adaptiveRadauTableau(uToltype, constvalue(tTypeNoUnits), num_stages) + + κ = 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:s + z[i] = w[i] = zero(u) + end + + dw1 = zero(u) + ubuff = zero(u) + dw2 = Vector{typeof(u)}(undef, floor(Int, num_stages/2)) + for i in 1 : floor(Int, num_stages/2) + dw2[i] = similar(u, Complex{eltype(u)}) + recursivefill!(dw[i], false) + end + cubuff = Vector{typeof(u)}(undef, floor(Int, num_stages/2)) + for i in 1 :floor(Int, num_stages/2) + cubuff[i] = similar(u, Complex{eltype(u)}) + recursivefill!(cubuff[i], false) + end + + cont = Vector{typeof(u)}(undef, num_stages - 1) + for i in 1: (num_stages - 1) + cont[i] = zero(u) + end + + fsalfirst = zero(rate_prototype) + fw = Vector{typeof(rate_prototype)}(undef, num_stages) + k = Vector{typeof(rate_prototype)}(undef, num_stages) + for i in 1: num_stages + k[i] = fw[i] = zero(rate_prototype) + end + + 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 RadauIIA5.") + end + W2 = vector{typeof(Complex{W1})}(undef, floor(Int, num_stages/2)) + for i in 1 : floor(Int, num_stages/2) + W2[i] = similar(J, Complex{eltype(W1)}) + recursivefill!(w2[i], false) + end + + du1 = zero(rate_prototype) + + tmp = Vector{typeof(u)}(undef, binomial(num_stages,2)) + for i in 1 : binomial(num_stages,2) + tmp[i] = zero(u) + end + + atmp = similar(u, uEltypeNoUnits) + recursivefill!(atmp, false) + + jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, 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 = Vector{typeof(linsolve1)}(undef, floor(Int, num_stages/2)) + for i in 1 : floor(int, num_stages/2) + linprob = LinearProblem(W2[i], _vec(cubuff[i]); u0 = _vec(dw2[i])) + linsolve2 = init(linprob, alg.linsolve, alias_A = true, alias_b = true, + assumptions = LinearSolve.OperatorAssumptions(true)) + end + + rtol = reltol isa Number ? reltol : zero(reltol) + atol = reltol isa Number ? reltol : zero(reltol) + + AdaptiveRadauCache(u, uprev, + z, w, dw1, ubuff, dw2, cubuff, cont, + du1, fsalfirst, 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 bb0009cdd7..1f0b418970 100644 --- a/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl +++ b/lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl @@ -26,32 +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) - 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) - integrator.k[1] = integrator.fsalfirst - integrator.k[2] = integrator.fsallast - nothing -end - -function initialize!(integrator, cache::RadauIIA9ConstantCache) +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 @@ -74,7 +49,7 @@ function initialize!(integrator, cache::RadauIIA3Cache) nothing end -function initialize!(integrator, cache::RadauIIA5Cache) +function initialize!(integrator, cache::Union{RadauIIA5Cache, RadauIIA9Cache, AdaptiveRadauCache}) integrator.kshortsize = 2 resize!(integrator.k, integrator.kshortsize) integrator.k[1] = integrator.fsalfirst @@ -94,25 +69,6 @@ function initialize!(integrator, cache::RadauIIA5Cache) nothing end -function initialize!(integrator, cache::RadauIIA9Cache) - 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^(5 / 8) / 10 - cache.atol = cache.rtol * (abstol / reltol) - else - @.. broadcast=false cache.rtol=reltol^(5 / 8) / 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 @@ -942,7 +898,7 @@ end z3 = @.. broadcast=false T31*w1+T32*w2+T33*w3+T34*w4+T35*w5 z4 = @.. broadcast=false T41*w1+T42*w2+T43*w3+T44*w4+T45*w5 z5 = @.. broadcast=false T51*w1+w2+w4 #= T52=1, T53=0, T54=1, T55=0 =# - + @show z1 # check stopping criterion iter > 1 && (η = θ / (1 - θ)) if η * ndw < κ && (iter > 1 || iszero(ndw) || !iszero(integrator.success_iter)) @@ -1018,7 +974,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, 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,26 +1028,26 @@ end c2′ = c2 * c5′ c3′ = c3 * c5′ c4′ = c4 * c5′ - z1 = @.. broadcast=false c1′*(cont1 + + @.. broadcast=false z1 = c1′*(cont1 + (c1′ - c3m1) * (cont2 + (c1′ - c2m1) * (cont3 + (c1′ - c1m1) * cont4))) - z2 = @.. broadcast=false c2′*(cont1 + + @.. broadcast=false z2 = c2′*(cont1 + (c2′ - c3m1) * (cont2 + (c2′ - c2m1) * (cont3 + (c2′ - c1m1) * cont4))) - z3 = @.. broadcast=false c3′*(cont1 + + @.. broadcast=false z3 = c3′*(cont1 + (c3′ - c3m1) * (cont2 + (c3′ - c2m1) * (cont3 + (c3′ - c1m1) * cont4))) - z4 = @.. broadcast=false c4′*(cont1 + + @.. broadcast=false z4 = c4′*(cont1 + (c4′ - c3m1) * (cont2 + (c4′ - c2m1) * (cont3 + (c4′ - c1m1) * cont4))) - z5 = @.. broadcast=false c5′*(cont1 + + @.. broadcast=false z5 = c5′*(cont1 + (c5′ - c3m1) * (cont2 + (c5′ - c2m1) * (cont3 + (c5′ - c1m1) * cont4))) - 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 + @.. broadcast=false w1 = TI11*z1+TI12*z2+TI13*z3+TI14*z4+TI15*z5 + @.. broadcast=false w2 = TI21*z1+TI22*z2+TI23*z3+TI24*z4+TI25*z5 + @.. broadcast=false w3 = TI31*z1+TI32*z2+TI33*z3+TI34*z4+TI35*z5 + @.. broadcast=false w4 = TI41*z1+TI42*z2+TI43*z3+TI44*z4+TI45*z5 + @.. broadcast=false w5 = TI51*z1+TI52*z2+TI53*z3+TI54*z4+TI55*z5 end # Newton iteration @@ -1333,3 +1289,443 @@ 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^(5 / 8)/10 + atol = @.. broadcast=false rtol*(abstol / reltol) + + γdt, αdt, βdt = γ / dt, α ./ dt, β ./ dt + + J = calc_J(integrator, cache) + LU = Vector{Any}(undef, Int((num_stages + 1) / 2)) + if u isa Number + LU[1] = -γdt * mass_matrix + J + for i in 2 : Int((num_stages + 1) / 2) + LU[i] = -(αdt[i - 1] + βdt[i - 1] * im) * mass_matrix + J + end + else + LU[1] = lu(-γdt * mass_matrix + J) + for i in 2 : Int((num_stages + 1) / 2) + LU[i] = lu(-(αdt[i - 1] + βdt[i - 1] * im) * mass_matrix + J) + end + end + integrator.stats.nw += 1 + z = w = Vector{BigFloat}(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] = w[i] = map(zero, u) + end + for i in 1 : (num_stages - 1) + cache.cont[i] = map(zero, u) + end + else + c_prime = Vector{eltype(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 - 1] * (c_prime[i] - c[1] + 1) + cont[num_stages - 2] + j = num_stages - 3 + while j > 0 + z[i] = @.. z[i] * (c_prime[i] - c[num_stages- j - 1] + 1) + cont[j] + j = j - 1 + end + z[i] = @.. z[i] * c_prime[i] + end + w = @.. TI * z + 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{eltype(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 + Mw = Vector{eltype(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{eltype(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[Int(i/2)] * Mw[i] + βdt[Int(i/2)] * Mw[i + 1] + rhs[i + 1] = @.. fw[i + 1] - βdt[Int(i/2)] * Mw[i] - αdt[Int(i/2)] * Mw[i + 1] + i += 2 + end + + dw = Vector{eltype(u)}(undef, num_stages) + dw[1] = LU[1] \ rhs[1] + for i in 2 : Int((num_stages + 1) / 2) + tmp = LU[i] \ (@.. rhs[2 * i - 2] + rhs[2 * i - 1] * im) + dw[2 * i - 2] = real(tmp) + dw[2 * i - 1] = imag(tmp) + end + integrator.stats.nsolve += Int((num_stages + 1) / 2) + + # compute norm of residuals + iter > 1 && (ndwprev = ndw) + atmp = Vector{eltype(u)}(undef, num_stages) + for i in 1 : num_stages + atmp[i] = calculate_residuals(dw[i], uprev, u, atol, rtol, internalnorm, t) + end + + ndw = 0 + for i in 1 : num_stages + ndw = ndw + internalnorm(atmp[i], 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 = @.. w - dw + + # transform `w` to `z` + z = @.. T * w + @show z[1] + + # 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 integrator.fsalfirst+tmp + alg.smooth_est && (utilde = LU[1] \ utilde; integrator.stats.nsolve += 1) + 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 f0+tmp + alg.smooth_est && (utilde = LU[1] \ utilde; integrator.stats.nsolve += 1) + 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{eltype(u)}(undef, num_stages - 1, num_stages - 1) + for i in 1 : (num_stages - 1) + for j in i : (num_stages - 1) + if i == 1 + derivatives[i, j] = @.. (z[i] - z[i + 1]) / (c[i] - c[i + 1]) #first derivatives + else + derivatives[i, j] = @.. (derivatives[i - 1, j - 1] - derivatives[i - 1, j]) / (c[j - i + 1] - c[j + 1]) #all others + end + end + end + for i in 1 : (num_stages - 1) + cache.cont[i] = derivatives[i, num_stages - 1] + 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, z, w = cache + @unpack dw1, ubuff, dw2, cubuff1, cubuff2 = cache + @unpack 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] = -(α[i]dt + β[i]dt * 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) + uzero = zero(eltype(u)) + for i in 1 : num_stages + @.. z[i] = w[i] = uzero + end + for i in 1 : (num_stages-1) + @.. cache.cont[i] = uzero + end + else + c_prime = Vector{eltype(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 - 1] * (c_prime[i] - c[1] + 1) + cont[num_stages - 2] + j = num_stages - 3 + while j > 0 + z[i] = @.. z[i] * (c_prime[i] - c[num_stages- j - 1] + 1) + cont[j] + j = j - 1 + end + z[i] = @.. z[i] * c_prime[i] + end + w = @.. TI * z + 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 + k[1] = fsallast + for i in 1 : num_stages + @.. tmp = uprev + z[i] + f(k[i], tmp, p, t + c[i] * dt) + end + integrator.stats.nf += num_stages + + @.. fw = TI * k + 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 + linres = Vector{BigFloat}(undef, num_stages) + if needfactor + linres[1] = dolinsolve(integrator, linsolve1; A = W1, b = _vec(ubuff), + linu = _vec(dw1)) + else + linres[1] = dolinsolve(integrator, linsolve1; A = nothing, b = _vec(ubuff), + linu = _vec(dw1)) + end + + cache.linsolve1 = linres1.cache + + for i in 1 : (num_stages - 1)/2 + @.. broadcast=false cubuff[i]=complex( + fw2 - αdt[i] * Mw[2 * i] + βdt[i] * Mw[2 * i + 1], fw3 - βdt[i] * Mw[2 * i] - αdt[i] * Mw[2 * i + 1]) + linsolve2[i] = cache.linsolve2[i] + if needfactor + linres[i + 1] = dolinsolve(integrator, linsolve2[i]; A = W2[i], b = _vec(cubuff[i]), + linu = _vec(dw2[i])) + else + linres[i + 1] = dolinsolve(integrator, linsolve2[i]; A = nothing, b = _vec(cubuff[i]), + linu = _vec(dw2[i])) + end + cache.linsolve2[i] = linres[i + 1].cache + end + + integrator.stats.nsolve += (num_stages + 1) / 2 + dw[1] = dw1 + i = 2 + while i <= num_stages + dw[i] = z[i] + dw[i + 1] = z[i + 1] + @.. dw[i] = real(dw2[i - 1]) + @.. dw[i + 1] = imag(dw2[i - 1]) + i += 2 + end + + # compute norm of residuals + iter > 1 && (ndwprev = ndw) + ndws = Vector{BigFloat}(undef, num_stages) + for i in 1:num_stages + calculate_residuals!(atmp, dw[i], uprev, u, atol, rtol, internalnorm, t) + ndws[i] = internalnorm(atmp, t) + end + + ndw = 0 + for i in 1 : num_stages + ndw += ndws[i] + 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 = w - dw + + # transform `w` to `z` + @.. z = T * w + # 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[s] + + step_limiter!(u, integrator, p, t + dt) + #= + if adaptive + utilde = w2 + edt = e./dt + @.. tmp= dot(edt, z) + 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 = Matrix{eltype(u)}(undef, num_stages - 1, num_stages - 1) + for i in 1 : (num_stages - 1) + for j in i : (num_stages - 1) + if i == 1 + @.. derivatives[i, j] = (z[i] - z[i + 1]) / (c[i] - c[i + 1]) #first derivatives + else + @.. derivatives[i, j] = (derivatives[i - 1, j - 1] - derivatives[i - 1, j]) / (c[j - i + 1] - c[j + 1]) #all others + end + end + end + for i in 1 : (num_stages - 1) + cache.cont[i] = derivatives[i, num_stages - 1] + 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..7860d6bfb7 100644 --- a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl +++ b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl @@ -258,3 +258,88 @@ function RadauIIA9Tableau(T, T2) γ, α1, β1, α2, β2, e1, e2, e3, e4, e5) end + +struct adaptiveRadauTableau{T, T2, Int} + T:: AbstractMatrix{T} + TI::AbstractMatrix{T} + γ::T + α::AbstractVector{T} + β::AbstractVector{T} + c::AbstractVector{T} + #e::AbstractVector{T} + num_stages::Int +end + +using Polynomials, GenericLinearAlgebra, LinearAlgebra, LinearSolve, GenericSchur + +function adaptiveRadauTableau(T, 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 + p = Polynomial{BigFloat}([tmp; tmp2]) + for i in 1:(num_stages - 1) + p = derivative(p) + end + c = roots(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(b[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) + #b_hat = Vector{BigFloat}(undef, num_stages) + #embedded = bseries(a, b_hat, c, num_stages - 2) + + #e = b_hat - b + adaptiveRadauTableau{Any, T2, Int}(T, TI, γ, α, β, c, num_stages) +end \ No newline at end of file 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..36289bdf7b 100644 --- a/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl +++ b/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl @@ -8,7 +8,12 @@ for prob in [prob_ode_linear, prob_ode_2Dlinear] @test sim21.𝒪est[:final]≈5 atol=testTol end -sim21 = test_convergence(1 ./ 2 .^ (2.5:-1:0.5), prob_ode_linear, RadauIIA9()) +sol = solve(prob_ode_linear, AdaptiveRadau(), adaptive = false, dt = 1e-2) +sol = solve(prob_ode_linear, RadauIIA9(), adaptive = false, dt = 1e-2) +sol = solve(prob_ode_linear, RadauIIA5(), adaptive = false, dt = 1e-2) + + +sim21 = test_convergence(1 ./ 10 .^ (4.5:-1:2.5), prob_ode_linear, AdaptiveRadau()) @test sim21.𝒪est[:final]≈8 atol=testTol sim21 = test_convergence(1 ./ 2 .^ (2.5:-1:0.5), prob_ode_2Dlinear, RadauIIA9())