diff --git a/Project.toml b/Project.toml index 1255e671bb..56d9074a34 100644 --- a/Project.toml +++ b/Project.toml @@ -83,7 +83,7 @@ SimpleNonlinearSolve = "1" SimpleUnPack = "1" SparseArrays = "1.9" SparseDiffTools = "2.3" -Static = "1" +Static = "0.8, 1" StaticArrayInterface = "1.2" StaticArrays = "1.0" TruncatedStacktraces = "1.2" diff --git a/lib/OrdinaryDiffEqLowStorageRK/test/runtests.jl b/lib/OrdinaryDiffEqLowStorageRK/test/runtests.jl index f823d8debf..7225e9390c 100644 --- a/lib/OrdinaryDiffEqLowStorageRK/test/runtests.jl +++ b/lib/OrdinaryDiffEqLowStorageRK/test/runtests.jl @@ -1,3 +1,3 @@ using SafeTestsets -@time @safetestset "Extrapolation Tests" include("ode_low_storage_rk_tests.jl") +@time @safetestset "Low Storage RK Tests" include("ode_low_storage_rk_tests.jl") diff --git a/lib/OrdinaryDiffEqSSPRK/Project.toml b/lib/OrdinaryDiffEqSSPRK/Project.toml new file mode 100644 index 0000000000..692af1906f --- /dev/null +++ b/lib/OrdinaryDiffEqSSPRK/Project.toml @@ -0,0 +1,23 @@ +name = "OrdinaryDiffEqSSPRK" +uuid = "669c94d9-1f4b-4b64-b377-1aa079aa2388" +authors = ["ParamThakkar123 "] +version = "1.0.0" + +[deps] +FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898" +MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" +OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" +RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" + +[compat] +julia = "1.10" + +[extras] +DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["DiffEqDevTools", "Random", "SafeTestsets", "Test"] \ No newline at end of file diff --git a/lib/OrdinaryDiffEqSSPRK/src/OrdinaryDiffEqSSPRK.jl b/lib/OrdinaryDiffEqSSPRK/src/OrdinaryDiffEqSSPRK.jl new file mode 100644 index 0000000000..658d21baef --- /dev/null +++ b/lib/OrdinaryDiffEqSSPRK/src/OrdinaryDiffEqSSPRK.jl @@ -0,0 +1,35 @@ +module OrdinaryDiffEqSSPRK + +import OrdinaryDiffEq: alg_order, calculate_residuals!, + initialize!, perform_step!, @unpack, unwrap_alg, + calculate_residuals, ssp_coefficient, + OrdinaryDiffEqAlgorithm, + OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, + OrdinaryDiffEqNewtonAdaptiveAlgorithm, + OrdinaryDiffEqRosenbrockAdaptiveAlgorithm, + OrdinaryDiffEqAdaptiveAlgorithm, uses_uprev, + alg_cache, _vec, _reshape, @cache, isfsal, full_cache, + constvalue, _unwrap_val, du_alias_or_new, + explicit_rk_docstring, trivial_limiter!, + _ode_interpolant, _ode_interpolant!, + _ode_addsteps! +using DiffEqBase, FastBroadcast, Polyester, MuladdMacro, RecursiveArrayTools +using DiffEqBase: @def +using Static: False + +import StaticArrays: SArray, MVector, SVector, @SVector, StaticArray, MMatrix, SA + +include("algorithms.jl") +include("alg_utils.jl") +include("ssprk_caches.jl") +include("interp_func.jl") +include("ssprk_perform_step.jl") +include("interpolants.jl") +include("addsteps.jl") +include("functions.jl") + +export SSPRK53_2N2, SSPRK22, SSPRK53, SSPRK63, SSPRK83, SSPRK43, SSPRK432, SSPRKMSVS32, + SSPRK54, SSPRK53_2N1, SSPRK104, SSPRK932, SSPRKMSVS43, SSPRK73, SSPRK53_H, + SSPRK33, SHLDDRK_2N, KYKSSPRK42, SHLDDRK52 + +end diff --git a/lib/OrdinaryDiffEqSSPRK/src/addsteps.jl b/lib/OrdinaryDiffEqSSPRK/src/addsteps.jl new file mode 100644 index 0000000000..5f31180d9e --- /dev/null +++ b/lib/OrdinaryDiffEqSSPRK/src/addsteps.jl @@ -0,0 +1,21 @@ +function _ode_addsteps!(k, t, uprev, u, dt, f, p, + cache::Union{SSPRK22ConstantCache, SSPRK33ConstantCache, + SSPRK43ConstantCache, SSPRK432ConstantCache}, + always_calc_begin = false, allow_calc_end = true, + force_calc_end = false) + if length(k) < 1 || always_calc_begin + copyat_or_push!(k, 1, f(uprev, p, t)) + end + nothing +end + +function _ode_addsteps!(k, t, uprev, u, dt, f, p, + cache::Union{SSPRK22Cache, SSPRK33Cache, SSPRK43Cache, + SSPRK432Cache}, always_calc_begin = false, + allow_calc_end = true, force_calc_end = false) + if length(k) < 1 || always_calc_begin + f(cache.k, uprev, p, t) + copyat_or_push!(k, 1, cache.k) + end + nothing +end \ No newline at end of file diff --git a/lib/OrdinaryDiffEqSSPRK/src/alg_utils.jl b/lib/OrdinaryDiffEqSSPRK/src/alg_utils.jl new file mode 100644 index 0000000000..e162e3a63d --- /dev/null +++ b/lib/OrdinaryDiffEqSSPRK/src/alg_utils.jl @@ -0,0 +1,52 @@ +isfsal(alg::SSPRK53_2N1) = false +isfsal(alg::SSPRK53_2N2) = false +isfsal(alg::SSPRK22) = false +isfsal(alg::SSPRK33) = false +isfsal(alg::SSPRK53) = false +isfsal(alg::SSPRK53_H) = false +isfsal(alg::SSPRK63) = false +isfsal(alg::SSPRK73) = false +isfsal(alg::SSPRK83) = false +isfsal(alg::SSPRK43) = false +isfsal(alg::SSPRK432) = false +isfsal(alg::SSPRK932) = false +isfsal(alg::SSPRK54) = false +isfsal(alg::SSPRK104) = false + +alg_order(alg::KYKSSPRK42) = 2 +alg_order(alg::SSPRKMSVS32) = 2 +alg_order(alg::SSPRK33) = 3 +alg_order(alg::SSPRK53_2N1) = 3 +alg_order(alg::SSPRK53_2N2) = 3 +alg_order(alg::SSPRK22) = 2 +alg_order(alg::SSPRK53) = 3 +alg_order(alg::SSPRK53_H) = 3 +alg_order(alg::SSPRK63) = 3 +alg_order(alg::SSPRK73) = 3 +alg_order(alg::SSPRK83) = 3 +alg_order(alg::SSPRK43) = 3 +alg_order(alg::SSPRK432) = 3 +alg_order(alg::SSPRKMSVS43) = 3 +alg_order(alg::SSPRK932) = 3 +alg_order(alg::SSPRK54) = 4 +alg_order(alg::SSPRK104) = 4 +alg_order(alg::SHLDDRK_2N) = 4 +alg_order(alg::SHLDDRK52) = 2 + +ssp_coefficient(alg::SSPRK53_2N1) = 2.18 +ssp_coefficient(alg::SSPRK53_2N2) = 2.148 +ssp_coefficient(alg::SSPRK53) = 2.65 +ssp_coefficient(alg::SSPRK53_H) = 2.65 +ssp_coefficient(alg::SSPRK63) = 3.518 +ssp_coefficient(alg::SSPRK73) = 4.2879 +ssp_coefficient(alg::SSPRK83) = 5.107 +ssp_coefficient(alg::SSPRK43) = 2 +ssp_coefficient(alg::SSPRK432) = 2 +ssp_coefficient(alg::SSPRK932) = 6 +ssp_coefficient(alg::SSPRK54) = 1.508 +ssp_coefficient(alg::SSPRK104) = 6 +ssp_coefficient(alg::SSPRK33) = 1 +ssp_coefficient(alg::SSPRK22) = 1 +ssp_coefficient(alg::SSPRKMSVS32) = 0.5 +ssp_coefficient(alg::SSPRKMSVS43) = 0.33 +ssp_coefficient(alg::KYKSSPRK42) = 2.459 \ No newline at end of file diff --git a/src/algorithms/explicit_rk_pde.jl b/lib/OrdinaryDiffEqSSPRK/src/algorithms.jl similarity index 99% rename from src/algorithms/explicit_rk_pde.jl rename to lib/OrdinaryDiffEqSSPRK/src/algorithms.jl index 68c97a9b10..579f6b29fe 100644 --- a/src/algorithms/explicit_rk_pde.jl +++ b/lib/OrdinaryDiffEqSSPRK/src/algorithms.jl @@ -1,3 +1,7 @@ +using OrdinaryDiffEq: explicit_rk_docstring +using Static: False +@inline trivial_limiter!(u, integrator, p, t) = nothing + @doc explicit_rk_docstring( "A third-order, five-stage explicit strong stability preserving (SSP) low-storage method. Fixed timestep only.", diff --git a/lib/OrdinaryDiffEqSSPRK/src/functions.jl b/lib/OrdinaryDiffEqSSPRK/src/functions.jl new file mode 100644 index 0000000000..c459cf8720 --- /dev/null +++ b/lib/OrdinaryDiffEqSSPRK/src/functions.jl @@ -0,0 +1,52 @@ +@eval @inline function DiffEqBase.get_tmp_cache(integrator, alg::SSPRK22, + cache::OrdinaryDiffEqConstantCache) + nothing +end + +@eval @inline function DiffEqBase.get_tmp_cache(integrator, alg::SSPRK33, + cache::OrdinaryDiffEqConstantCache) + nothing +end + +@eval @inline function DiffEqBase.get_tmp_cache(integrator, alg::SSPRK53_2N1, + cache::OrdinaryDiffEqConstantCache) + nothing +end + +@eval @inline function DiffEqBase.get_tmp_cache(integrator, alg::SSPRK53_2N2, + cache::OrdinaryDiffEqConstantCache) + nothing +end + +@eval @inline function DiffEqBase.get_tmp_cache(integrator, alg::SSPRK432, + cache::OrdinaryDiffEqConstantCache) + nothing +end + +@eval @inline function DiffEqBase.get_tmp_cache(integrator, alg::SSPRK932, + cache::OrdinaryDiffEqConstantCache) + nothing +end + +@eval @inline function DiffEqBase.get_tmp_cache(integrator, alg::OrdinaryDiffEqNewtonAdaptiveAlgorithm, + cache::OrdinaryDiffEqConstantCache) + nothing +end + +@eval @inline function DiffEqBase.get_tmp_cache(integrator, alg::OrdinaryDiffEqRosenbrockAdaptiveAlgorithm, + cache::OrdinaryDiffEqConstantCache) + nothing +end + +@eval @inline function DiffEqBase.get_tmp_cache(integrator, alg::OrdinaryDiffEqAlgorithm, + cache::OrdinaryDiffEqConstantCache) + nothing +end + +@inline function DiffEqBase.get_tmp_cache(integrator, + alg::Union{SSPRK22, SSPRK33, SSPRK53_2N1, + SSPRK53_2N2, SSPRK43, SSPRK432, + SSPRK932}, + cache::OrdinaryDiffEqMutableCache) +(cache.k,) +end \ No newline at end of file diff --git a/lib/OrdinaryDiffEqSSPRK/src/interp_func.jl b/lib/OrdinaryDiffEqSSPRK/src/interp_func.jl new file mode 100644 index 0000000000..062871ca4f --- /dev/null +++ b/lib/OrdinaryDiffEqSSPRK/src/interp_func.jl @@ -0,0 +1,10 @@ +function DiffEqBase.interp_summary(::Type{cacheType}, + dense::Bool) where { + cacheType <: + Union{SSPRK22, SSPRK22ConstantCache, + SSPRK33, SSPRK33ConstantCache, + SSPRK43, SSPRK43ConstantCache, + SSPRK432, SSPRK432ConstantCache +}} +dense ? "2nd order \"free\" SSP interpolation" : "1st order linear" +end \ No newline at end of file diff --git a/lib/OrdinaryDiffEqSSPRK/src/interpolants.jl b/lib/OrdinaryDiffEqSSPRK/src/interpolants.jl new file mode 100644 index 0000000000..a5258bb49c --- /dev/null +++ b/lib/OrdinaryDiffEqSSPRK/src/interpolants.jl @@ -0,0 +1,135 @@ +const NEGZERO = Float16(-0.0f0) + +@def ssprkpre0 begin + c00 = @evalpoly(Θ, 1, NEGZERO, -1) + c10 = Θ^2 + b10dt = Θ * @evalpoly(Θ, 1, -1) * dt +end + +@def ssprkpre1 begin + b10diff = @evalpoly(Θ, 1, -2) + c10diffinvdt = 2Θ * inv(dt) # = -c00diff * inv(dt) +end + +@def ssprkpre2 begin + invdt = inv(dt) + b10diff2invdt = -2 * invdt + c10diff2invdt2 = 2 * invdt^2 # = -c00diff2 * inv(dt)^2 +end + +@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, + cache::Union{SSPRK22ConstantCache, SSPRK22Cache, + SSPRK33ConstantCache, SSPRK33Cache, + SSPRK43ConstantCache, SSPRK43Cache, + SSPRK432ConstantCache, SSPRK432Cache}, + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) +@ssprkpre0 +@inbounds @.. broadcast=false y₀*c00+y₁*c10+k[1]*b10dt +end + +@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, + cache::Union{SSPRK22ConstantCache, SSPRK22Cache, + SSPRK33ConstantCache, SSPRK33Cache, + SSPRK43ConstantCache, SSPRK43Cache, + SSPRK432ConstantCache, SSPRK432Cache}, idxs, + T::Type{Val{0}}, differential_vars::Nothing) +@ssprkpre0 +@views @.. broadcast=false y₀[idxs]*c00+y₁[idxs]*c10+k[1][idxs]*b10dt +end + +@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, + cache::Union{SSPRK22ConstantCache, SSPRK22Cache, + SSPRK33ConstantCache, SSPRK33Cache, + SSPRK43ConstantCache, SSPRK43Cache, + SSPRK432ConstantCache, SSPRK432Cache}, + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) +@ssprkpre0 +@inbounds @.. broadcast=false out=y₀ * c00 + y₁ * c10 + k[1] * b10dt +out +end + +@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, + cache::Union{SSPRK22ConstantCache, SSPRK22Cache, + SSPRK33ConstantCache, SSPRK33Cache, + SSPRK43ConstantCache, SSPRK43Cache, + SSPRK432ConstantCache, SSPRK432Cache}, idxs, + T::Type{Val{0}}, differential_vars::Nothing) +@ssprkpre0 +@views @.. broadcast=false out=y₀[idxs] * c00 + y₁[idxs] * c10 + k[1][idxs] * b10dt +out +end + +@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, + cache::Union{SSPRK22ConstantCache, SSPRK22Cache, + SSPRK33ConstantCache, SSPRK33Cache, + SSPRK43ConstantCache, SSPRK43Cache, + SSPRK432ConstantCache, SSPRK432Cache}, + idxs::Nothing, T::Type{Val{1}}, differential_vars::Nothing) +@ssprkpre1 +@inbounds @.. broadcast=false (y₁ - y₀) * c10diffinvdt+k[1] * b10diff +end + +@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, + cache::Union{SSPRK22ConstantCache, SSPRK22Cache, + SSPRK33ConstantCache, SSPRK33Cache, + SSPRK43ConstantCache, SSPRK43Cache, + SSPRK432ConstantCache, SSPRK432Cache}, idxs, + T::Type{Val{1}}, differential_vars::Nothing) +@ssprkpre1 +@views @.. broadcast=false (y₁[idxs] - y₀[idxs]) * c10diffinvdt+k[1][idxs] * b10diff +end + +@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, + cache::Union{SSPRK22ConstantCache, SSPRK22Cache, + SSPRK33ConstantCache, SSPRK33Cache, + SSPRK43ConstantCache, SSPRK43Cache, + SSPRK432ConstantCache, SSPRK432Cache}, + idxs::Nothing, T::Type{Val{1}}, differential_vars::Nothing) +@ssprkpre1 +@inbounds @.. broadcast=false out=(y₁ - y₀) * c10diffinvdt + k[1] * b10diff +out +end + +@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, + cache::Union{SSPRK22ConstantCache, SSPRK22Cache, + SSPRK33ConstantCache, SSPRK33Cache, + SSPRK43ConstantCache, SSPRK43Cache, + SSPRK432ConstantCache, SSPRK432Cache}, idxs, + T::Type{Val{1}}, differential_vars::Nothing) + @ssprkpre1 + @views @.. broadcast=false out=(y₁[idxs] - y₀[idxs]) * c10diffinvdt + + k[1][idxs] * b10diff + out +end + +@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, + cache::Union{SSPRK22ConstantCache, SSPRK22Cache, + SSPRK33ConstantCache, SSPRK33Cache, + SSPRK43ConstantCache, SSPRK43Cache, + SSPRK432ConstantCache, SSPRK432Cache}, + idxs::Nothing, T::Type{Val{2}}, differential_vars::Nothing) + @ssprkpre2 + @inbounds @.. broadcast=false (y₁ - y₀) * c10diff2invdt2+k[1] * b10diff2invdt +end + +@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, + cache::Union{SSPRK22ConstantCache, SSPRK22Cache, + SSPRK33ConstantCache, SSPRK33Cache, + SSPRK43ConstantCache, SSPRK43Cache, + SSPRK432ConstantCache, SSPRK432Cache}, idxs, + T::Type{Val{2}}, differential_vars::Nothing) +@ssprkpre2 +@views @.. broadcast=false (y₁[idxs] - y₀[idxs]) * + c10diff2invdt2+k[1][idxs] * b10diff2invdt +end + +@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, + cache::Union{SSPRK22ConstantCache, SSPRK22Cache, + SSPRK33ConstantCache, SSPRK33Cache, + SSPRK43ConstantCache, SSPRK43Cache, + SSPRK432ConstantCache, SSPRK432Cache}, + idxs::Nothing, T::Type{Val{2}}, differential_vars::Nothing) +@ssprkpre2 +@inbounds @.. broadcast=false out=(y₁ - y₀) * c10diff2invdt2 + k[1] * b10diff2invdt +out +end diff --git a/src/caches/ssprk_caches.jl b/lib/OrdinaryDiffEqSSPRK/src/ssprk_caches.jl similarity index 100% rename from src/caches/ssprk_caches.jl rename to lib/OrdinaryDiffEqSSPRK/src/ssprk_caches.jl diff --git a/src/perform_step/ssprk_perform_step.jl b/lib/OrdinaryDiffEqSSPRK/src/ssprk_perform_step.jl similarity index 99% rename from src/perform_step/ssprk_perform_step.jl rename to lib/OrdinaryDiffEqSSPRK/src/ssprk_perform_step.jl index 4f39c2e83c..c332ae335a 100644 --- a/src/perform_step/ssprk_perform_step.jl +++ b/lib/OrdinaryDiffEqSSPRK/src/ssprk_perform_step.jl @@ -1704,4 +1704,4 @@ end stage_limiter!(u, integrator, p, t + dt) step_limiter!(u, integrator, p, t + dt) integrator.stats.nf += 10 -end +end \ No newline at end of file diff --git a/test/algconvergence/ode_ssprk_tests.jl b/lib/OrdinaryDiffEqSSPRK/test/ode_ssprk_tests.jl similarity index 100% rename from test/algconvergence/ode_ssprk_tests.jl rename to lib/OrdinaryDiffEqSSPRK/test/ode_ssprk_tests.jl diff --git a/lib/OrdinaryDiffEqSSPRK/test/runtests.jl b/lib/OrdinaryDiffEqSSPRK/test/runtests.jl new file mode 100644 index 0000000000..939f4eae4e --- /dev/null +++ b/lib/OrdinaryDiffEqSSPRK/test/runtests.jl @@ -0,0 +1,3 @@ +using SafeTestsets + +@time @safetestset "SSPRK Tests" include("ode_ssprk_tests.jl") \ No newline at end of file diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index 258af062ad..b511653080 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -133,7 +133,6 @@ include("misc_utils.jl") include("algorithms.jl") include("algorithms/explicit_rk.jl") -include("algorithms/explicit_rk_pde.jl") include("alg_utils.jl") @@ -149,7 +148,6 @@ include("composite_algs.jl") include("caches/basic_caches.jl") include("caches/low_order_rk_caches.jl") include("caches/high_order_rk_caches.jl") -include("caches/ssprk_caches.jl") include("caches/feagin_caches.jl") include("caches/verner_caches.jl") include("caches/sdirk_caches.jl") @@ -198,7 +196,6 @@ include("perform_step/low_order_rk_perform_step.jl") include("perform_step/high_order_rk_perform_step.jl") include("perform_step/verner_rk_perform_step.jl") include("perform_step/feagin_rk_perform_step.jl") -include("perform_step/ssprk_perform_step.jl") include("perform_step/sdirk_perform_step.jl") include("perform_step/kencarp_kvaerno_perform_step.jl") include("perform_step/firk_perform_step.jl") @@ -262,6 +259,12 @@ export ORK256, CarpenterKennedy2N54, SHLDDRK64, HSLDDRK64, DGLDDRK73_C, DGLDDRK8 RDPK3Sp35, RDPK3SpFSAL35, RDPK3Sp49, RDPK3SpFSAL49, RDPK3Sp510, RDPK3SpFSAL510, KYK2014DGSSPRK_3S2, RK46NL +include("../lib/OrdinaryDiffEqSSPRK/src/OrdinaryDiffEqSSPRK.jl") +using ..OrdinaryDiffEqSSPRK +export SSPRK53_2N2, SSPRK22, SSPRK53, SSPRK63, SSPRK83, SSPRK43, SSPRK432, SSPRKMSVS32, + SSPRK54, SSPRK53_2N1, SSPRK104, SSPRK932, SSPRKMSVS43, SSPRK73, SSPRK53_H, + SSPRK33, SHLDDRK_2N, KYKSSPRK42, SHLDDRK52 + import PrecompileTools PrecompileTools.@compile_workload begin @@ -398,10 +401,6 @@ export FunctionMap, Euler, Heun, Ralston, Midpoint, RK4, ExplicitRK, OwrenZen3, Vern9, Feagin10, Feagin12, Feagin14, CompositeAlgorithm, Anas5, RKO65, FRK65, PFRK87, RKM, MSRK5, MSRK6, Stepanov5, SIR54, QPRK98, PSRK4p7q6, PSRK3p6q5, PSRK3p5q4 -export SSPRK22, SSPRK33, KYKSSPRK42, SSPRK53, SSPRK53_2N1, SSPRK53_2N2, SSPRK53_H, SSPRK63, - SSPRK73, SSPRK83, SSPRK43, SSPRK432, - SSPRKMSVS32, SSPRKMSVS43, SSPRK932, SSPRK54, SSPRK104 - export RadauIIA3, RadauIIA5 export ImplicitEuler, ImplicitMidpoint, Trapezoid, TRBDF2, SDIRK2, SDIRK22, diff --git a/src/alg_utils.jl b/src/alg_utils.jl index 1a084d3b9b..55fd14cc4b 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -33,20 +33,6 @@ isfsal(alg::FunctionMap) = false isfsal(alg::Rodas3P) = false isfsal(alg::Rodas23W) = false isfsal(alg::Rodas5) = false -isfsal(alg::SSPRK53_2N1) = false -isfsal(alg::SSPRK53_2N2) = false -isfsal(alg::SSPRK22) = false -isfsal(alg::SSPRK33) = false -isfsal(alg::SSPRK53) = false -isfsal(alg::SSPRK53_H) = false -isfsal(alg::SSPRK63) = false -isfsal(alg::SSPRK73) = false -isfsal(alg::SSPRK83) = false -isfsal(alg::SSPRK43) = false -isfsal(alg::SSPRK432) = false -isfsal(alg::SSPRK932) = false -isfsal(alg::SSPRK54) = false -isfsal(alg::SSPRK104) = false isfsal(alg::Rodas5P) = false isfsal(alg::Rodas5Pr) = false isfsal(alg::Rodas5Pe) = false @@ -414,9 +400,6 @@ alg_order(alg::Ralston) = 2 alg_order(alg::LawsonEuler) = 1 alg_order(alg::NorsettEuler) = 1 alg_order(alg::LieEuler) = 1 -alg_order(alg::KYKSSPRK42) = 2 -alg_order(alg::SSPRKMSVS32) = 2 -alg_order(alg::SSPRK33) = 3 alg_order(alg::CayleyEuler) = 2 alg_order(alg::ETDRK2) = 2 alg_order(alg::ETDRK3) = 3 @@ -433,7 +416,6 @@ alg_order(alg::SplitEuler) = 1 alg_order(alg::ETD2) = 2 alg_order(alg::Exprb32) = 3 alg_order(alg::Exprb43) = 4 -alg_order(alg::SHLDDRK_2N) = 4 alg_order(alg::Anas5) = 5 alg_order(alg::KuttaPRK2p5) = 5 alg_order(alg::RKO65) = 5 @@ -550,20 +532,6 @@ alg_order(alg::SFSDIRK8) = 4 alg_order(alg::Hairer4) = 4 alg_order(alg::Hairer42) = 4 alg_order(alg::Feagin10) = 10 -alg_order(alg::SSPRK53_2N1) = 3 -alg_order(alg::SSPRK53_2N2) = 3 -alg_order(alg::SSPRK22) = 2 -alg_order(alg::SSPRK53) = 3 -alg_order(alg::SSPRK53_H) = 3 -alg_order(alg::SSPRK63) = 3 -alg_order(alg::SSPRK73) = 3 -alg_order(alg::SSPRK83) = 3 -alg_order(alg::SSPRK43) = 3 -alg_order(alg::SSPRK432) = 3 -alg_order(alg::SSPRKMSVS43) = 3 -alg_order(alg::SSPRK932) = 3 -alg_order(alg::SSPRK54) = 4 -alg_order(alg::SSPRK104) = 4 alg_order(alg::Feagin12) = 12 alg_order(alg::Feagin14) = 14 alg_order(alg::PFRK87) = 8 @@ -602,7 +570,6 @@ alg_order(alg::Rodas5) = 5 alg_order(alg::Rodas5P) = 5 alg_order(alg::Rodas5Pr) = 5 alg_order(alg::Rodas5Pe) = 5 -alg_order(alg::SHLDDRK52) = 2 alg_order(alg::AB3) = 3 alg_order(alg::AB4) = 4 @@ -748,23 +715,6 @@ end # SSP coefficients ssp_coefficient(alg) = error("$alg is not a strong stability preserving method.") ssp_coefficient(alg::Euler) = 1 -ssp_coefficient(alg::SSPRK53_2N1) = 2.18 -ssp_coefficient(alg::SSPRK53_2N2) = 2.148 -ssp_coefficient(alg::SSPRK53) = 2.65 -ssp_coefficient(alg::SSPRK53_H) = 2.65 -ssp_coefficient(alg::SSPRK63) = 3.518 -ssp_coefficient(alg::SSPRK73) = 4.2879 -ssp_coefficient(alg::SSPRK83) = 5.107 -ssp_coefficient(alg::SSPRK43) = 2 -ssp_coefficient(alg::SSPRK432) = 2 -ssp_coefficient(alg::SSPRK932) = 6 -ssp_coefficient(alg::SSPRK54) = 1.508 -ssp_coefficient(alg::SSPRK104) = 6 -ssp_coefficient(alg::SSPRK33) = 1 -ssp_coefficient(alg::SSPRK22) = 1 -ssp_coefficient(alg::SSPRKMSVS32) = 0.5 -ssp_coefficient(alg::SSPRKMSVS43) = 0.33 -ssp_coefficient(alg::KYKSSPRK42) = 2.459 # We shouldn't do this probably. #ssp_coefficient(alg::ImplicitEuler) = Inf diff --git a/src/dense/interpolants.jl b/src/dense/interpolants.jl index 9a658997bd..7356b8865d 100644 --- a/src/dense/interpolants.jl +++ b/src/dense/interpolants.jl @@ -1,9 +1,5 @@ RK_WITH_SPECIAL_INTERPOLATIONS = Union{FunctionMapConstantCache, FunctionMapCache, DP5ConstantCache, DP5Cache, - SSPRK22ConstantCache, SSPRK22Cache, - SSPRK33ConstantCache, SSPRK33Cache, - SSPRK43ConstantCache, SSPRK43Cache, - SSPRK432ConstantCache, SSPRK432Cache, Tsit5ConstantCache, Tsit5Cache, OwrenZen3ConstantCache, OwrenZen3Cache, OwrenZen4ConstantCache, OwrenZen4Cache, @@ -255,158 +251,6 @@ end out end -const NEGZERO = Float16(-0.0f0) -""" -Second order strong stability preserving (SSP) interpolant. - -Ketcheson, Lóczi, Jangabylova, Kusmanov: Dense output for SSP RK methods (2017). -""" -@def ssprkpre0 begin - c00 = @evalpoly(Θ, 1, NEGZERO, -1) - c10 = Θ^2 - b10dt = Θ * @evalpoly(Θ, 1, -1) * dt -end - -@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, - cache::Union{SSPRK22ConstantCache, SSPRK22Cache, - SSPRK33ConstantCache, SSPRK33Cache, - SSPRK43ConstantCache, SSPRK43Cache, - SSPRK432ConstantCache, SSPRK432Cache}, - idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) - @ssprkpre0 - @inbounds @.. broadcast=false y₀*c00+y₁*c10+k[1]*b10dt -end - -@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, - cache::Union{SSPRK22ConstantCache, SSPRK22Cache, - SSPRK33ConstantCache, SSPRK33Cache, - SSPRK43ConstantCache, SSPRK43Cache, - SSPRK432ConstantCache, SSPRK432Cache}, idxs, - T::Type{Val{0}}, differential_vars::Nothing) - @ssprkpre0 - @views @.. broadcast=false y₀[idxs]*c00+y₁[idxs]*c10+k[1][idxs]*b10dt -end - -@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, - cache::Union{SSPRK22ConstantCache, SSPRK22Cache, - SSPRK33ConstantCache, SSPRK33Cache, - SSPRK43ConstantCache, SSPRK43Cache, - SSPRK432ConstantCache, SSPRK432Cache}, - idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) - @ssprkpre0 - @inbounds @.. broadcast=false out=y₀ * c00 + y₁ * c10 + k[1] * b10dt - out -end - -@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, - cache::Union{SSPRK22ConstantCache, SSPRK22Cache, - SSPRK33ConstantCache, SSPRK33Cache, - SSPRK43ConstantCache, SSPRK43Cache, - SSPRK432ConstantCache, SSPRK432Cache}, idxs, - T::Type{Val{0}}, differential_vars::Nothing) - @ssprkpre0 - @views @.. broadcast=false out=y₀[idxs] * c00 + y₁[idxs] * c10 + k[1][idxs] * b10dt - out -end - -@def ssprkpre1 begin - b10diff = @evalpoly(Θ, 1, -2) - c10diffinvdt = 2Θ * inv(dt) # = -c00diff * inv(dt) -end - -@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, - cache::Union{SSPRK22ConstantCache, SSPRK22Cache, - SSPRK33ConstantCache, SSPRK33Cache, - SSPRK43ConstantCache, SSPRK43Cache, - SSPRK432ConstantCache, SSPRK432Cache}, - idxs::Nothing, T::Type{Val{1}}, differential_vars::Nothing) - @ssprkpre1 - @inbounds @.. broadcast=false (y₁ - y₀) * c10diffinvdt+k[1] * b10diff -end - -@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, - cache::Union{SSPRK22ConstantCache, SSPRK22Cache, - SSPRK33ConstantCache, SSPRK33Cache, - SSPRK43ConstantCache, SSPRK43Cache, - SSPRK432ConstantCache, SSPRK432Cache}, idxs, - T::Type{Val{1}}, differential_vars::Nothing) - @ssprkpre1 - @views @.. broadcast=false (y₁[idxs] - y₀[idxs]) * c10diffinvdt+k[1][idxs] * b10diff -end - -@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, - cache::Union{SSPRK22ConstantCache, SSPRK22Cache, - SSPRK33ConstantCache, SSPRK33Cache, - SSPRK43ConstantCache, SSPRK43Cache, - SSPRK432ConstantCache, SSPRK432Cache}, - idxs::Nothing, T::Type{Val{1}}, differential_vars::Nothing) - @ssprkpre1 - @inbounds @.. broadcast=false out=(y₁ - y₀) * c10diffinvdt + k[1] * b10diff - out -end - -@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, - cache::Union{SSPRK22ConstantCache, SSPRK22Cache, - SSPRK33ConstantCache, SSPRK33Cache, - SSPRK43ConstantCache, SSPRK43Cache, - SSPRK432ConstantCache, SSPRK432Cache}, idxs, - T::Type{Val{1}}, differential_vars::Nothing) - @ssprkpre1 - @views @.. broadcast=false out=(y₁[idxs] - y₀[idxs]) * c10diffinvdt + - k[1][idxs] * b10diff - out -end - -@def ssprkpre2 begin - invdt = inv(dt) - b10diff2invdt = -2 * invdt - c10diff2invdt2 = 2 * invdt^2 # = -c00diff2 * inv(dt)^2 -end - -@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, - cache::Union{SSPRK22ConstantCache, SSPRK22Cache, - SSPRK33ConstantCache, SSPRK33Cache, - SSPRK43ConstantCache, SSPRK43Cache, - SSPRK432ConstantCache, SSPRK432Cache}, - idxs::Nothing, T::Type{Val{2}}, differential_vars::Nothing) - @ssprkpre2 - @inbounds @.. broadcast=false (y₁ - y₀) * c10diff2invdt2+k[1] * b10diff2invdt -end - -@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, - cache::Union{SSPRK22ConstantCache, SSPRK22Cache, - SSPRK33ConstantCache, SSPRK33Cache, - SSPRK43ConstantCache, SSPRK43Cache, - SSPRK432ConstantCache, SSPRK432Cache}, idxs, - T::Type{Val{2}}, differential_vars::Nothing) - @ssprkpre2 - @views @.. broadcast=false (y₁[idxs] - y₀[idxs]) * - c10diff2invdt2+k[1][idxs] * b10diff2invdt -end - -@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, - cache::Union{SSPRK22ConstantCache, SSPRK22Cache, - SSPRK33ConstantCache, SSPRK33Cache, - SSPRK43ConstantCache, SSPRK43Cache, - SSPRK432ConstantCache, SSPRK432Cache}, - idxs::Nothing, T::Type{Val{2}}, differential_vars::Nothing) - @ssprkpre2 - @inbounds @.. broadcast=false out=(y₁ - y₀) * c10diff2invdt2 + k[1] * b10diff2invdt - out -end - -@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, - cache::Union{SSPRK22ConstantCache, SSPRK22Cache, - SSPRK33ConstantCache, SSPRK33Cache, - SSPRK43ConstantCache, SSPRK43Cache, - SSPRK432ConstantCache, SSPRK432Cache}, idxs, - T::Type{Val{2}}, differential_vars::Nothing) - @ssprkpre2 - @views @.. broadcast=false out=(y₁[idxs] - y₀[idxs]) * c10diff2invdt2 + - k[1][idxs] * b10diff2invdt - out -end - """ Runge–Kutta pairs of order 5(4) satisfying only the first column simplifying assumption diff --git a/src/dense/low_order_rk_addsteps.jl b/src/dense/low_order_rk_addsteps.jl index 7e93300f25..9e609808b1 100644 --- a/src/dense/low_order_rk_addsteps.jl +++ b/src/dense/low_order_rk_addsteps.jl @@ -10,28 +10,6 @@ function _ode_addsteps!(k, t, uprev, u, dt, f, p, cache::FunctionMapConstantCach nothing end -function _ode_addsteps!(k, t, uprev, u, dt, f, p, - cache::Union{SSPRK22ConstantCache, SSPRK33ConstantCache, - SSPRK43ConstantCache, SSPRK432ConstantCache}, - always_calc_begin = false, allow_calc_end = true, - force_calc_end = false) - if length(k) < 1 || always_calc_begin - copyat_or_push!(k, 1, f(uprev, p, t)) - end - nothing -end - -function _ode_addsteps!(k, t, uprev, u, dt, f, p, - cache::Union{SSPRK22Cache, SSPRK33Cache, SSPRK43Cache, - SSPRK432Cache}, always_calc_begin = false, - allow_calc_end = true, force_calc_end = false) - if length(k) < 1 || always_calc_begin - f(cache.k, uprev, p, t) - copyat_or_push!(k, 1, cache.k) - end - nothing -end - @muladd function _ode_addsteps!(k, t, uprev, u, dt, f, p, cache::OwrenZen4Cache, always_calc_begin = false, allow_calc_end = true, force_calc_end = false) diff --git a/src/integrators/integrator_interface.jl b/src/integrators/integrator_interface.jl index 140abddefb..2763656877 100644 --- a/src/integrators/integrator_interface.jl +++ b/src/integrators/integrator_interface.jl @@ -1,6 +1,7 @@ # We want to make sure that the first argument of change_t_via_interpolation! # is specialized, yet, it needs to take both ODEIntegrator and DDEIntegrator. # Hence, we need to have two separate functions. +using OrdinaryDiffEq function _change_t_via_interpolation!(integrator, t, modify_save_endpoint::Type{Val{T}}) where {T} # Can get rid of an allocation here with a function @@ -111,15 +112,12 @@ end get_tmp_cache(integrator::ODEIntegrator, integrator.alg, integrator.cache) end # avoid method ambiguity -for typ in (OrdinaryDiffEqAlgorithm, Union{RadauIIA3, RadauIIA5}, - OrdinaryDiffEqNewtonAdaptiveAlgorithm, - OrdinaryDiffEqRosenbrockAdaptiveAlgorithm, - Union{SSPRK22, SSPRK33, SSPRK53_2N1, SSPRK53_2N2, SSPRK43, SSPRK432, SSPRK932}) - @eval @inline function DiffEqBase.get_tmp_cache(integrator, alg::$typ, - cache::OrdinaryDiffEqConstantCache) - nothing - end -end +# for typ in (Union{RadauIIA3, RadauIIA5}) +# @eval @inline function DiffEqBase.get_tmp_cache(integrator, alg::$typ, +# cache::OrdinaryDiffEqConstantCache) +# nothing +# end +# end # the ordering of the cache arrays is important!!! @inline function DiffEqBase.get_tmp_cache(integrator, alg::OrdinaryDiffEqAlgorithm, @@ -144,13 +142,7 @@ end cache::OrdinaryDiffEqMutableCache) (cache.tmp, cache.linsolve_tmp) end -@inline function DiffEqBase.get_tmp_cache(integrator, - alg::Union{SSPRK22, SSPRK33, SSPRK53_2N1, - SSPRK53_2N2, SSPRK43, SSPRK432, - SSPRK932}, - cache::OrdinaryDiffEqMutableCache) - (cache.k,) -end + @inline function DiffEqBase.get_tmp_cache(integrator, alg::OrdinaryDiffEqAdaptiveExponentialAlgorithm, cache::OrdinaryDiffEqMutableCache) diff --git a/src/interp_func.jl b/src/interp_func.jl index 1f6614e46f..14c95b7aba 100644 --- a/src/interp_func.jl +++ b/src/interp_func.jl @@ -71,16 +71,7 @@ function DiffEqBase.interp_summary(::Type{cacheType}, dense ? "specialized 4rd order \"free\" stiffness-aware interpolation" : "1st order linear" end -function DiffEqBase.interp_summary(::Type{cacheType}, - dense::Bool) where { - cacheType <: - Union{SSPRK22, SSPRK22ConstantCache, - SSPRK33, SSPRK33ConstantCache, - SSPRK43, SSPRK43ConstantCache, - SSPRK432, SSPRK432ConstantCache -}} - dense ? "2nd order \"free\" SSP interpolation" : "1st order linear" -end + function DiffEqBase.interp_summary(::Type{cacheType}, dense::Bool) where { cacheType <: Union{OwrenZen3Cache, diff --git a/test/runtests.jl b/test/runtests.jl index d557bc2f1a..7e71b2c45a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -35,12 +35,24 @@ function activate_stabilized_rk() Pkg.instantiate() end +function activate_stabilized_irk() + Pkg.activate("../lib/OrdinaryDiffEqStabilizedIRK") + Pkg.develop(PackageSpec(path = dirname(@__DIR__))) + Pkg.instantiate() +end + function activate_low_storage_rk() Pkg.activate("../lib/OrdinaryDiffEqLowStorageRK") Pkg.develop(PackageSpec(path = dirname(@__DIR__))) Pkg.instantiate() end +function activate_ssprk() + Pkg.activate("../lib/OrdinaryDiffEqSSPRK") + Pkg.develop(PackageSpec(path = dirname(@__DIR__))) + Pkg.instantiate() +end + #Start Test Script @time begin @@ -164,8 +176,6 @@ end end if !is_APPVEYOR && GROUP == "AlgConvergence_II" - @time @safetestset "SSPRK Tests" include("algconvergence/ode_ssprk_tests.jl") - @time @safetestset "Low Storage RK Tests" include("../lib/OrdinaryDiffEqLowStorageRK/test/ode_low_storage_rk_tests.jl") @time @safetestset "OwrenZen Tests" include("algconvergence/owrenzen_tests.jl") @time @safetestset "Runge-Kutta-Chebyshev Tests" include("../lib/OrdinaryDiffEqStabilizedRK/test/rkc_tests.jl") end @@ -194,6 +204,15 @@ end @time @safetestset "StabilizedIRK Tests" include("../lib/OrdinaryDiffEqStabilizedIRK/test/runtests.jl") end + if !is_APPVEYOR && GROUP == "SSPRK" + @time @safetestset "SSPRK Tests" include("../lib/OrdinaryDiffEqSSPRK/test/ode_ssprk_tests.jl") + end + + if !is_APPVEYOR && GROUP == "Low Storage RK" + @time @safetestset "Low Storage RK Tests" include("../lib/OrdinaryDiffEqLowStorageRK/test/ode_low_storage_rk_tests.jl") + end + + if !is_APPVEYOR && GROUP == "Downstream" activate_downstream_env() @time @safetestset "DelayDiffEq Tests" include("downstream/delaydiffeq.jl")