From 3d09fa24078c6b58d2c6c162d050196603138a18 Mon Sep 17 00:00:00 2001 From: Shreyas Ekanathan Date: Tue, 20 Aug 2024 20:53:53 -0400 Subject: [PATCH] lots of edits --- lib/OrdinaryDiffEqFIRK/Project.toml | 3 + .../src/OrdinaryDiffEqFIRK.jl | 1 + lib/OrdinaryDiffEqFIRK/src/alg_utils.jl | 1 + lib/OrdinaryDiffEqFIRK/src/algorithms.jl | 2 +- lib/OrdinaryDiffEqFIRK/src/firk_caches.jl | 202 +++++++++--------- .../src/firk_perform_step.jl | 121 ++++------- lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl | 12 +- lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl | 7 +- 8 files changed, 157 insertions(+), 192 deletions(-) diff --git a/lib/OrdinaryDiffEqFIRK/Project.toml b/lib/OrdinaryDiffEqFIRK/Project.toml index 63aa2f368c..02c7ace7b9 100644 --- a/lib/OrdinaryDiffEqFIRK/Project.toml +++ b/lib/OrdinaryDiffEqFIRK/Project.toml @@ -15,6 +15,9 @@ MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8" +GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a" +GenericSchur = "c145ed77-6b09-5dd9-b285-bf645a82121e" +Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" [compat] julia = "1.10" diff --git a/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl b/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl index 5d071e8777..e26805bc6e 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, isfirk using MuladdMacro, DiffEqBase, RecursiveArrayTools +using Polynomials, GenericLinearAlgebra, GenericSchur using SciMLOperators: AbstractSciMLOperator using LinearAlgebra: I, UniformScaling, mul!, lu import LinearSolve diff --git a/lib/OrdinaryDiffEqFIRK/src/alg_utils.jl b/lib/OrdinaryDiffEqFIRK/src/alg_utils.jl index 5953ecffc1..7b0e8cdc1a 100644 --- a/lib/OrdinaryDiffEqFIRK/src/alg_utils.jl +++ b/lib/OrdinaryDiffEqFIRK/src/alg_utils.jl @@ -3,6 +3,7 @@ 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 diff --git a/lib/OrdinaryDiffEqFIRK/src/algorithms.jl b/lib/OrdinaryDiffEqFIRK/src/algorithms.jl index ba4c51c0d0..4429fb78b6 100644 --- a/lib/OrdinaryDiffEqFIRK/src/algorithms.jl +++ b/lib/OrdinaryDiffEqFIRK/src/algorithms.jl @@ -171,7 +171,7 @@ end function AdaptiveRadau(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(), concrete_jac = nothing, - diff_type = Val{:forward}, num_stages = 5, + 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, diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl index 40f6631765..8b531cda7d 100644 --- a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl +++ b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl @@ -468,41 +468,42 @@ function alg_cache(alg::RadauIIA9, u, rate_prototype, ::Type{uEltypeNoUnits}, Convergence, alg.step_limiter!) end -mutable struct adaptiveRadauConstantCache{S, F, Tab, Tol, Dt, U, JType} <: +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 + 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) -tab = adaptiveRadau(uToltype, constvalue(tTypeNoUnits), alg.num_stages) - -cont = Vector{typeof(u)}(undef, num_stages - 1) -for i in 1: (num_stages - 1) - cont[i] = zero(u) -end + ::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)' + κ = 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) + AdaptiveRadauConstantCache(uf, tab, κ, one(uToltype), 10000, cont, dt, dt, + Convergence, J) end -mutable struct adaptiveRadauCache{uType, cuType, uNoUnitsType, rateType, JType, W1Type, W2Type, +mutable struct AdaptiveRadauCache{uType, cuType, uNoUnitsType, rateType, JType, W1Type, W2Type, UF, JC, F1, F2, Tab, Tol, Dt, rTol, aTol, StepLimiter} <: OrdinaryDiffEqMutableCache u::uType @@ -540,89 +541,90 @@ mutable struct adaptiveRadauCache{uType, cuType, uNoUnitsType, rateType, 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{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} -uf = UJacobianWrapper(f, t, p) -uToltype = constvalue(uBottomEltypeNoUnits) -tab = RadauIIA9Tableau(uToltype, constvalue(tTypeNoUnits)) - -κ = 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 + ::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) -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 + κ = alg.κ !== nothing ? convert(uToltype, alg.κ) : convert(uToltype, 1 // 100) -cont = Vector{typeof(u)}(undef, num_stages - 1) -for i in 1: (num_stages - 1) - cont[i] = zero(u) -end + 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 -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 + 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 -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 + cont = Vector{typeof(u)}(undef, num_stages - 1) + for i in 1: (num_stages - 1) + cont[i] = zero(u) + end -du1 = zero(rate_prototype) + 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 -tmp = Vector{typeof(u)}(undef, binomial(num_stages,2)) -for i in 1 : binomial(num_stages,2) - tmp[i] = zero(u) -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 -atmp = similar(u, uEltypeNoUnits) -recursivefill!(atmp, false) + du1 = zero(rate_prototype) -jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, dw1) + tmp = Vector{typeof(u)}(undef, binomial(num_stages,2)) + for i in 1 : binomial(num_stages,2) + tmp[i] = zero(u) + end -linprob = LinearProblem(W1, _vec(ubuff); u0 = _vec(dw1)) -linsolve1 = init(linprob, alg.linsolve, alias_A = true, alias_b = true, - assumptions = LinearSolve.OperatorAssumptions(true)) + atmp = similar(u, uEltypeNoUnits) + recursivefill!(atmp, false) -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 + 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) -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!) + 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 51313203a6..1f52f50aca 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 - integrator.stats.nf += 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 @@ -76,7 +51,7 @@ function initialize!(integrator, cache::RadauIIA3Cache) nothing end -function initialize!(integrator, cache::RadauIIA5Cache) +function initialize!(integrator, cache::Union{RadauIIA5Cache, RadauIIA9Cache, AdaptiveRadauCache}) integrator.kshortsize = 2 integrator.fsalfirst = cache.fsalfirst integrator.fsallast = cache.k @@ -98,27 +73,6 @@ function initialize!(integrator, cache::RadauIIA5Cache) nothing end -function initialize!(integrator, cache::RadauIIA9Cache) - integrator.kshortsize = 2 - integrator.fsalfirst = cache.fsalfirst - integrator.fsallast = cache.k - 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) - integrator.stats.nf += 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 @@ -948,7 +902,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)) @@ -1340,10 +1294,10 @@ end return end -@muladd function perform_step!(integrator, cache::adaptiveRadauConstantCache, +@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 T, TI, γ, α, β, c, #=e,=# num_stages = cache.tab @unpack κ, cont = cache @unpack internalnorm, abstol, reltol, adaptive = integrator.opts alg = unwrap_alg(integrator, true) @@ -1357,16 +1311,16 @@ end γdt, αdt, βdt = γ / dt, α ./ dt, β ./ dt J = calc_J(integrator, cache) - LU = Vector{Any}(undef, (num_stages + 1) / 2) + LU = Vector{Any}(undef, Int((num_stages + 1) / 2)) if u isa Number LU[1] = -γdt * mass_matrix + J - for i in 2 : (num_stages + 1) / 2 - LU[i] = -(α[i - 1]dt + β[i - 1]dt * im) * 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 : (num_stages + 1) / 2 - LU[i] = lu(-(α[i - 1]dt + β[i - 1]dt * im) * 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 @@ -1377,21 +1331,22 @@ end z[i] = w[i] = map(zero, u) end for i in 1 : (num_stages - 1) - cont[i] = map(zero, u) + cache.cont[i] = map(zero, u) end else - c' = Vector{eltype(u)}(undef, num_stages) #time stepping - c'[num_stages] = dt / cache.dtprev + 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'[i] = c[i] * c'[num_stages] + 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'[i] - c[1] + 1) + cont[num_stages - 1] - j = s - 2 + 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'[i] - c[num_stages-j] + 1) + cont[j] + z[i] = @.. z[i] * (c_prime[i] - c[num_stages- j - 1] + 1) + cont[j] + j = j - 1 end - z[i] = @.. z[i] * c'[i] + z[i] = @.. z[i] * c_prime[i] end w = @.. TI * z end @@ -1410,7 +1365,7 @@ end for i in 1 : num_stages ff[i] = f(uprev + z[i], p, t + c[i] * dt) end - integrator.stats.nf += 5 + integrator.stats.nf += num_stages fw = @.. TI * ff Mw = Vector{eltype(u)}(undef, num_stages) @@ -1426,19 +1381,19 @@ end 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] + 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] = LU1 \ rhs[1] - for i in 2 : (num_stages + 1) / 2 + 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.nlsolve += (num_stages + 1) / 2 + integrator.stats.nsolve += Int((num_stages + 1) / 2) # compute norm of residuals iter > 1 && (ndwprev = ndw) @@ -1467,6 +1422,7 @@ end # transform `w` to `z` z = @.. T * w + @show z[1] # check stopping criterion iter > 1 && (η = θ / (1 - θ)) @@ -1536,9 +1492,9 @@ end return end -@muladd function perform_step!(integrator, cache::adaptiveRadauCache, repeat_step = false) +@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 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 @@ -1573,21 +1529,22 @@ end for i in 1 : (num_stages-1) @.. cache.cont[i] = uzero end - else - c' = Vector{eltype(u)}(undef, num_stages) #time stepping - c'[num_stages] = dt / cache.dtprev - for i in 1 : (num_stages-1) - c'[i] = c[i] * c'[num_stages] + 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'[i] - c[1] + 1) + cont[num_stages - 1] - j = num_stages - 2 + 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'[i] - c[num_stages - j] + 1) + cont[j] + z[i] = @.. z[i] * (c_prime[i] - c[num_stages- j - 1] + 1) + cont[j] + j = j - 1 end - @.. z[i] = z[i] * c'[i] + z[i] = @.. z[i] * c_prime[i] end - @.. w = TI * z + w = @.. TI * z end # Newton iteration diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl index 3e36fd9363..7860d6bfb7 100644 --- a/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl +++ b/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl @@ -266,11 +266,11 @@ struct adaptiveRadauTableau{T, T2, Int} α::AbstractVector{T} β::AbstractVector{T} c::AbstractVector{T} - e::AbstractVector{T} - S::Int + #e::AbstractVector{T} + num_stages::Int end -using Polynomials, GenericLinearAlgebra, LinearAlgebra, LinearSolve, GenericSchur, BSeries +using Polynomials, GenericLinearAlgebra, LinearAlgebra, LinearSolve, GenericSchur function adaptiveRadauTableau(T, T2, num_stages::Int) tmp = Vector{BigFloat}(undef, num_stages - 1) @@ -341,7 +341,5 @@ function adaptiveRadauTableau(T, T2, num_stages::Int) #embedded = bseries(a, b_hat, c, num_stages - 2) #e = b_hat - b - #adaptiveRadautableau(T, TI, γ, α, β, c, e, s) -end - -adaptiveRadauTableau(0, 0, 3) + adaptiveRadauTableau{Any, T2, Int}(T, TI, γ, α, β, c, num_stages) +end \ No newline at end of file diff --git a/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl b/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl index 0c960a61df..36289bdf7b 100644 --- a/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl +++ b/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl @@ -8,9 +8,12 @@ for prob in [prob_ode_linear, prob_ode_2Dlinear] @test sim21.𝒪est[:final]≈5 atol=testTol end -sol = solve(prob_ode_linear, AdaptiveRadau()) +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 ./ 2 .^ (2.5:-1:0.5), prob_ode_linear, RadauIIA9()) + +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())