diff --git a/Project.toml b/Project.toml index 9cc17d5a6c..9194b91f37 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "OrdinaryDiffEq" uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" authors = ["Chris Rackauckas ", "Yingbo Ma "] -version = "6.87.0" +version = "6.88.1" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" @@ -98,19 +98,33 @@ MacroTools = "0.5" MuladdMacro = "0.2.1" NonlinearSolve = "3" OrdinaryDiffEqAdamsBashforthMoulton = "1" +OrdinaryDiffEqBDF = "1" OrdinaryDiffEqCore = "1" +OrdinaryDiffEqDefault = "1" OrdinaryDiffEqDifferentiation = "1" +OrdinaryDiffEqExplicitRK = "1" +OrdinaryDiffEqExponentialRK = "1" OrdinaryDiffEqExtrapolation = "1" +OrdinaryDiffEqFIRK = "1" OrdinaryDiffEqFeagin = "1" OrdinaryDiffEqFunctionMap = "1" +OrdinaryDiffEqHighOrderRK = "1" +OrdinaryDiffEqIMEXMultistep = "1" +OrdinaryDiffEqLinear = "1" OrdinaryDiffEqLowOrderRK = "1" OrdinaryDiffEqLowStorageRK = "1" -OrdinaryDiffEqHighOrderRK = "1" OrdinaryDiffEqNonlinearSolve = "1" +OrdinaryDiffEqNordsieck = "1" +OrdinaryDiffEqPDIRK = "1" OrdinaryDiffEqPRK = "1" OrdinaryDiffEqQPRK = "1" OrdinaryDiffEqRKN = "1" +OrdinaryDiffEqRosenbrock = "1" +OrdinaryDiffEqSDIRK = "1" +OrdinaryDiffEqStabilizedIRK = "1" +OrdinaryDiffEqSSPRK = "1" OrdinaryDiffEqStabilizedRK = "1" +OrdinaryDiffEqSymplecticRK = "1" OrdinaryDiffEqTsit5 = "1" OrdinaryDiffEqVerner = "1" Polyester = "0.7" diff --git a/lib/OrdinaryDiffEqAdamsBashforthMoulton/src/adams_bashforth_moulton_caches.jl b/lib/OrdinaryDiffEqAdamsBashforthMoulton/src/adams_bashforth_moulton_caches.jl index fda72a8cd8..d417f9cdad 100644 --- a/lib/OrdinaryDiffEqAdamsBashforthMoulton/src/adams_bashforth_moulton_caches.jl +++ b/lib/OrdinaryDiffEqAdamsBashforthMoulton/src/adams_bashforth_moulton_caches.jl @@ -1,7 +1,9 @@ abstract type ABMMutableCache <: OrdinaryDiffEqMutableCache end abstract type ABMVariableCoefficientMutableCache <: OrdinaryDiffEqMutableCache end -get_fsalfirstlast(cache::ABMMutableCache,u) = (cache.fsalfirst, cache.k) -get_fsalfirstlast(cache::ABMVariableCoefficientMutableCache,u) = (cache.fsalfirst, cache.k4) +get_fsalfirstlast(cache::ABMMutableCache, u) = (cache.fsalfirst, cache.k) +function get_fsalfirstlast(cache::ABMVariableCoefficientMutableCache, u) + (cache.fsalfirst, cache.k4) +end @cache mutable struct AB3Cache{uType, rateType} <: ABMMutableCache u::uType uprev::uType diff --git a/lib/OrdinaryDiffEqBDF/src/algorithms.jl b/lib/OrdinaryDiffEqBDF/src/algorithms.jl index 1a4c960c7b..4f25e268a7 100644 --- a/lib/OrdinaryDiffEqBDF/src/algorithms.jl +++ b/lib/OrdinaryDiffEqBDF/src/algorithms.jl @@ -150,7 +150,7 @@ end function QNDF1(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(), concrete_jac = nothing, diff_type = Val{:forward}, linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), - extrapolant = :linear, kappa = -37//200, + extrapolant = :linear, kappa = -37 // 200, controller = :Standard, step_limiter! = trivial_limiter!) QNDF1{ _unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve), @@ -233,7 +233,8 @@ function QNDF(; max_order::Val{MO} = Val{5}(), chunk_size = Val{0}(), diff_type = Val{:forward}, linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), κ = nothing, tol = nothing, - extrapolant = :linear, kappa = (-37//200, -1//9, -823//10000, -83//2000, 0//1), + extrapolant = :linear, kappa = ( + -37 // 200, -1 // 9, -823 // 10000, -83 // 2000, 0 // 1), controller = :Standard, step_limiter! = trivial_limiter!) where {MO} QNDF{MO, _unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), diff --git a/lib/OrdinaryDiffEqBDF/src/bdf_caches.jl b/lib/OrdinaryDiffEqBDF/src/bdf_caches.jl index 9d60023c13..00c6875cef 100644 --- a/lib/OrdinaryDiffEqBDF/src/bdf_caches.jl +++ b/lib/OrdinaryDiffEqBDF/src/bdf_caches.jl @@ -1,5 +1,7 @@ abstract type BDFMutableCache <: OrdinaryDiffEqMutableCache end -get_fsalfirstlast(cache::BDFMutableCache,u) = (cache.fsalfirst, du_alias_or_new(cache.nlsolver, cache.fsalfirst)) +function get_fsalfirstlast(cache::BDFMutableCache, u) + (cache.fsalfirst, du_alias_or_new(cache.nlsolver, cache.fsalfirst)) +end @cache mutable struct ABDF2ConstantCache{N, dtType, rate_prototype} <: OrdinaryDiffEqConstantCache diff --git a/lib/OrdinaryDiffEqBDF/src/controllers.jl b/lib/OrdinaryDiffEqBDF/src/controllers.jl index c001e6a586..fc13dbd91a 100644 --- a/lib/OrdinaryDiffEqBDF/src/controllers.jl +++ b/lib/OrdinaryDiffEqBDF/src/controllers.jl @@ -206,7 +206,7 @@ function choose_order!(alg::FBDF, integrator, Val(max_order)) local terk_tmp if u isa Number - terk_tmp = fd_weights[k - 2, 1]*u + terk_tmp = fd_weights[k - 2, 1] * u for i in 2:(k - 2) terk_tmp += fd_weights[i, k - 2] * u_history[i - 1] end @@ -215,7 +215,7 @@ function choose_order!(alg::FBDF, integrator, # we need terk_tmp to be mutable. # so it can be updated terk_tmp = similar(u) - @.. terk_tmp = fd_weights[k - 2, 1]*_vec(u) + @.. terk_tmp = fd_weights[k - 2, 1] * _vec(u) for i in 2:(k - 2) @.. @views terk_tmp += fd_weights[i, k - 2] * u_history[:, i - 1] end diff --git a/lib/OrdinaryDiffEqBDF/src/dae_caches.jl b/lib/OrdinaryDiffEqBDF/src/dae_caches.jl index 00bfb4ca5a..708a3d119f 100644 --- a/lib/OrdinaryDiffEqBDF/src/dae_caches.jl +++ b/lib/OrdinaryDiffEqBDF/src/dae_caches.jl @@ -1,5 +1,7 @@ abstract type DAEBDFMutableCache <: OrdinaryDiffEqMutableCache end -get_fsalfirstlast(cache::DAEBDFMutableCache,u) = (cache.fsalfirst, du_alias_or_new(cache.nlsolver, cache.fsalfirst)) +function get_fsalfirstlast(cache::DAEBDFMutableCache, u) + (cache.fsalfirst, du_alias_or_new(cache.nlsolver, cache.fsalfirst)) +end @cache mutable struct DImplicitEulerCache{uType, rateType, uNoUnitsType, N} <: DAEBDFMutableCache @@ -13,7 +15,7 @@ get_fsalfirstlast(cache::DAEBDFMutableCache,u) = (cache.fsalfirst, du_alias_or_n end # Not FSAL -get_fsalfirstlast(cache::DImplicitEulerCache,u) = (u,u) +get_fsalfirstlast(cache::DImplicitEulerCache, u) = (u, u) mutable struct DImplicitEulerConstantCache{N} <: OrdinaryDiffEqConstantCache nlsolver::N diff --git a/lib/OrdinaryDiffEqBDF/src/dae_perform_step.jl b/lib/OrdinaryDiffEqBDF/src/dae_perform_step.jl index fd9b007963..d372a89a5e 100644 --- a/lib/OrdinaryDiffEqBDF/src/dae_perform_step.jl +++ b/lib/OrdinaryDiffEqBDF/src/dae_perform_step.jl @@ -160,8 +160,6 @@ end end function initialize!(integrator, cache::DABDF2Cache) - - integrator.kshortsize = 2 @unpack k₁, k₂ = cache.eulercache resize!(integrator.k, integrator.kshortsize) diff --git a/lib/OrdinaryDiffEqCore/Project.toml b/lib/OrdinaryDiffEqCore/Project.toml index f5228b185c..9a950e8230 100644 --- a/lib/OrdinaryDiffEqCore/Project.toml +++ b/lib/OrdinaryDiffEqCore/Project.toml @@ -1,7 +1,7 @@ name = "OrdinaryDiffEqCore" uuid = "bbf590c4-e513-4bbe-9b18-05decba2e5d8" authors = ["ParamThakkar123 "] -version = "1.2.0" +version = "1.3.0" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" diff --git a/lib/OrdinaryDiffEqCore/ext/OrdinaryDiffEqCoreEnzymeCoreExt.jl b/lib/OrdinaryDiffEqCore/ext/OrdinaryDiffEqCoreEnzymeCoreExt.jl index c43dd98c3e..b8bc336b57 100644 --- a/lib/OrdinaryDiffEqCore/ext/OrdinaryDiffEqCoreEnzymeCoreExt.jl +++ b/lib/OrdinaryDiffEqCore/ext/OrdinaryDiffEqCoreEnzymeCoreExt.jl @@ -1,12 +1,20 @@ module OrdinaryDiffEqCoreEnzymeCoreExt import OrdinaryDiffEqCore, EnzymeCore -Enzyme.EnzymeCore.EnzymeRules.inactive(::typeof(OrdinaryDiffEqCore.increment_nf!), args...) = true -Enzyme.EnzymeCore.EnzymeRules.inactive(::typeof(OrdinaryDiffEqCore.fixed_t_for_floatingpoint_error!), args...) = true -Enzyme.EnzymeCore.EnzymeRules.inactive(::typeof(OrdinaryDiffEqCore.increment_accept!), args...) = true -Enzyme.EnzymeCore.EnzymeRules.inactive(::typeof(OrdinaryDiffEqCore.increment_reject!), args...) = true -Enzyme.EnzymeCore.EnzymeRules.inactive(::typeof(OrdinaryDiffEqCore.increment_nf_perform_step!), args...) = true -Enzyme.EnzymeCore.EnzymeRules.inactive(::typeof(OrdinaryDiffEqCore.check_error!), args...) = true -Enzyme.EnzymeCore.EnzymeRules.inactive(::typeof(OrdinaryDiffEqCore.log_step!), args...) = true +EnzymeCore.EnzymeRules.inactive(::typeof(OrdinaryDiffEqCore.increment_nf!), args...) = true +function EnzymeCore.EnzymeRules.inactive( + ::typeof(OrdinaryDiffEqCore.fixed_t_for_floatingpoint_error!), args...) + true +end +function EnzymeCore.EnzymeRules.inactive( + ::typeof(OrdinaryDiffEqCore.increment_accept!), args...) + true +end +function EnzymeCore.EnzymeRules.inactive( + ::typeof(OrdinaryDiffEqCore.increment_reject!), args...) + true +end +EnzymeCore.EnzymeRules.inactive(::typeof(OrdinaryDiffEqCore.check_error!), args...) = true +EnzymeCore.EnzymeRules.inactive(::typeof(OrdinaryDiffEqCore.log_step!), args...) = true -end \ No newline at end of file +end diff --git a/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl b/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl index 6cbb4ec20a..bc36e5c7a0 100644 --- a/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl +++ b/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl @@ -41,7 +41,7 @@ using SimpleUnPack, RecursiveArrayTools, DataStructures, ArrayInterface import TruncatedStacktraces import StaticArraysCore: SArray, MVector, SVector, StaticArray, MMatrix, - StaticMatrix + StaticMatrix # Integrator Interface import DiffEqBase: resize!, deleteat!, addat!, full_cache, user_cache, u_cache, du_cache, @@ -54,7 +54,6 @@ import DiffEqBase: resize!, deleteat!, addat!, full_cache, user_cache, u_cache, isautodifferentiable, get_tstops, get_tstops_array, get_tstops_max - using DiffEqBase: check_error!, @def, _vec, _reshape using FastBroadcast: @.., True, False diff --git a/lib/OrdinaryDiffEqCore/src/alg_utils.jl b/lib/OrdinaryDiffEqCore/src/alg_utils.jl index b8d1fdd9ba..955540f2d6 100644 --- a/lib/OrdinaryDiffEqCore/src/alg_utils.jl +++ b/lib/OrdinaryDiffEqCore/src/alg_utils.jl @@ -345,6 +345,18 @@ qsteady_max_default(alg::OrdinaryDiffEqImplicitAlgorithm) = isadaptive(alg) ? 1 #DiffEqBase.nlsolve_default(::QNDF, ::Val{κ}) = 1//2 # SSP coefficients + +""" + ssp_coefficient(alg) +Return the SSP coefficient of the ODE algorithm `alg`. If one time step of size +`dt` with `alg` can be written as a convex combination of explicit Euler steps +with step sizes `cᵢ * dt`, the SSP coefficient is the minimal value of `1/cᵢ`. +# Examples +```julia-repl +julia> ssp_coefficient(SSPRK104()) +6 +``` +""" ssp_coefficient(alg) = error("$alg is not a strong stability preserving method.") # We shouldn't do this probably. diff --git a/lib/OrdinaryDiffEqCore/src/caches/basic_caches.jl b/lib/OrdinaryDiffEqCore/src/caches/basic_caches.jl index f8aef04cda..088fcc4b76 100644 --- a/lib/OrdinaryDiffEqCore/src/caches/basic_caches.jl +++ b/lib/OrdinaryDiffEqCore/src/caches/basic_caches.jl @@ -5,7 +5,7 @@ struct ODEEmptyCache <: OrdinaryDiffEqConstantCache end struct ODEChunkCache{CS} <: OrdinaryDiffEqConstantCache end # Don't worry about the potential alloc on a constant cache -get_fsalfirstlast(cache::OrdinaryDiffEqConstantCache,u) = (zero(u), zero(u)) +get_fsalfirstlast(cache::OrdinaryDiffEqConstantCache, u) = (zero(u), zero(u)) mutable struct CompositeCache{T, F} <: OrdinaryDiffEqCache caches::T @@ -13,7 +13,7 @@ mutable struct CompositeCache{T, F} <: OrdinaryDiffEqCache current::Int end -get_fsalfirstlast(cache::CompositeCache,u) = get_fsalfirstlast(cache.caches[1],u) +get_fsalfirstlast(cache::CompositeCache, u) = get_fsalfirstlast(cache.caches[1], u) mutable struct DefaultCache{T1, T2, T3, T4, T5, T6, A, F, uType} <: OrdinaryDiffEqCache args::A @@ -28,12 +28,13 @@ mutable struct DefaultCache{T1, T2, T3, T4, T5, T6, A, F, uType} <: OrdinaryDiff cache6::T6 function DefaultCache{T1, T2, T3, T4, T5, T6, F, uType}( args, choice_function, current, u) where {T1, T2, T3, T4, T5, T6, F, uType} - new{T1, T2, T3, T4, T5, T6, typeof(args), F, uType}(args, choice_function, current, u) + new{T1, T2, T3, T4, T5, T6, typeof(args), F, uType}( + args, choice_function, current, u) end end -function get_fsalfirstlast(cache::DefaultCache,u) - (cache.u,cache.u) # will be overwritten by the cache choice +function get_fsalfirstlast(cache::DefaultCache, u) + (cache.u, cache.u) # will be overwritten by the cache choice end function alg_cache(alg::CompositeAlgorithm, u, rate_prototype, ::Type{uEltypeNoUnits}, diff --git a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl index 3024f40827..fc56ed8d50 100644 --- a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl +++ b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl @@ -255,8 +255,8 @@ function _loopfooter!(integrator) end if integrator.opts.progress && integrator.iter % integrator.opts.progress_steps == 0 log_step!(integrator.opts.progress_name, integrator.opts.progress_id, - integrator.opts.progress_message, integrator.dt, integrator.u, - integrator.p, integrator.t, integrator.sol.prob.tspan) + integrator.opts.progress_message, integrator.dt, integrator.u, + integrator.p, integrator.t, integrator.sol.prob.tspan) end # Take value because if t is dual then maxeig can be dual @@ -280,18 +280,18 @@ end function log_step!(progress_name, progress_id, progress_message, dt, u, p, t, tspan) t1, t2 = tspan - @logmsg(LogLevel(-1),progress_name, - _id=progress_id, - message=progress_message(dt, u, p, t), - progress=(t - t1) / (t2 - t1)) + @logmsg(LogLevel(-1), progress_name, + _id=progress_id, + message=progress_message(dt, u, p, t), + progress=(t - t1) / (t2 - t1)) end function fixed_t_for_floatingpoint_error!(integrator, ttmp) if has_tstop(integrator) tstop = integrator.tdir * first_tstop(integrator) if abs(ttmp - tstop) < - 100eps(float(max(integrator.t, tstop) / oneunit(integrator.t))) * - oneunit(integrator.t) + 100eps(float(max(integrator.t, tstop) / oneunit(integrator.t))) * + oneunit(integrator.t) tstop else ttmp diff --git a/lib/OrdinaryDiffEqCore/src/interp_func.jl b/lib/OrdinaryDiffEqCore/src/interp_func.jl index f0c29d4967..a2bcd6b96c 100644 --- a/lib/OrdinaryDiffEqCore/src/interp_func.jl +++ b/lib/OrdinaryDiffEqCore/src/interp_func.jl @@ -62,7 +62,6 @@ end # strip interpolation of function information function SciMLBase.strip_interpolation(id::InterpolationData) - cache = strip_cache(id.cache) InterpolationData(nothing, id.timeseries, @@ -78,7 +77,7 @@ end function strip_cache(cache) if hasfield(typeof(cache), :jac_config) || hasfield(typeof(cache), :grad_config) fieldnums = length(fieldnames(typeof(cache))) - noth_list = fill(nothing,fieldnums) + noth_list = fill(nothing, fieldnums) cache_type_name = Base.typename(typeof(cache)).wrapper cache_type_name(noth_list...) else diff --git a/lib/OrdinaryDiffEqCore/src/perform_step/composite_perform_step.jl b/lib/OrdinaryDiffEqCore/src/perform_step/composite_perform_step.jl index c3b277a7bc..91d6d2a781 100644 --- a/lib/OrdinaryDiffEqCore/src/perform_step/composite_perform_step.jl +++ b/lib/OrdinaryDiffEqCore/src/perform_step/composite_perform_step.jl @@ -32,40 +32,40 @@ function initialize!(integrator, cache::DefaultCache) init_ith_default_cache(cache, algs, cache.current) u = integrator.u if cache.current == 1 - fsalfirst, fsallast = get_fsalfirstlast(cache.cache1,u) + fsalfirst, fsallast = get_fsalfirstlast(cache.cache1, u) integrator.fsalfirst = fsalfirst integrator.fsallast = fsallast initialize!(integrator, cache.cache1) elseif cache.current == 2 - fsalfirst, fsallast = get_fsalfirstlast(cache.cache2,u) + fsalfirst, fsallast = get_fsalfirstlast(cache.cache2, u) integrator.fsalfirst = fsalfirst integrator.fsallast = fsallast initialize!(integrator, cache.cache2) # the controller was initialized by default for algs[1] reset_alg_dependent_opts!(integrator.opts.controller, algs[1], algs[2]) elseif cache.current == 3 - fsalfirst, fsallast = get_fsalfirstlast(cache.cache3,u) + fsalfirst, fsallast = get_fsalfirstlast(cache.cache3, u) integrator.fsalfirst = fsalfirst integrator.fsallast = fsallast initialize!(integrator, cache.cache3) # the controller was initialized by default for algs[1] reset_alg_dependent_opts!(integrator.opts.controller, algs[1], algs[3]) elseif cache.current == 4 - fsalfirst, fsallast = get_fsalfirstlast(cache.cache4,u) + fsalfirst, fsallast = get_fsalfirstlast(cache.cache4, u) integrator.fsalfirst = fsalfirst integrator.fsallast = fsallast initialize!(integrator, cache.cache4) # the controller was initialized by default for algs[1] reset_alg_dependent_opts!(integrator.opts.controller, algs[1], algs[4]) elseif cache.current == 5 - fsalfirst, fsallast = get_fsalfirstlast(cache.cache5,u) + fsalfirst, fsallast = get_fsalfirstlast(cache.cache5, u) integrator.fsalfirst = fsalfirst integrator.fsallast = fsallast initialize!(integrator, cache.cache5) # the controller was initialized by default for algs[1] reset_alg_dependent_opts!(integrator.opts.controller, algs[1], algs[5]) elseif cache.current == 6 - fsalfirst, fsallast = get_fsalfirstlast(cache.cache6,u) + fsalfirst, fsallast = get_fsalfirstlast(cache.cache6, u) integrator.fsalfirst = fsalfirst integrator.fsallast = fsallast initialize!(integrator, cache.cache6) @@ -79,12 +79,12 @@ function initialize!(integrator, cache::CompositeCache) cache.current = cache.choice_function(integrator) u = integrator.u if cache.current == 1 - fsalfirst, fsallast = get_fsalfirstlast(cache.caches[1],u) + fsalfirst, fsallast = get_fsalfirstlast(cache.caches[1], u) integrator.fsalfirst = fsalfirst integrator.fsallast = fsallast initialize!(integrator, @inbounds(cache.caches[1])) elseif cache.current == 2 - fsalfirst, fsallast = get_fsalfirstlast(cache.caches[2],u) + fsalfirst, fsallast = get_fsalfirstlast(cache.caches[2], u) integrator.fsalfirst = fsalfirst integrator.fsallast = fsallast initialize!(integrator, @inbounds(cache.caches[2])) @@ -92,7 +92,7 @@ function initialize!(integrator, cache::CompositeCache) reset_alg_dependent_opts!(integrator.opts.controller, integrator.alg.algs[1], integrator.alg.algs[2]) else - fsalfirst, fsallast = get_fsalfirstlast(cache.caches[cache.current],u) + fsalfirst, fsallast = get_fsalfirstlast(cache.caches[cache.current], u) integrator.fsalfirst = fsalfirst integrator.fsallast = fsallast initialize!(integrator, @inbounds(cache.caches[cache.current])) @@ -106,12 +106,12 @@ function initialize!(integrator, cache::CompositeCache{Tuple{T1, T2}, F}) where cache.current = cache.choice_function(integrator) u = integrator.u if cache.current == 1 - fsalfirst, fsallast = get_fsalfirstlast(cache.caches[1],u) + fsalfirst, fsallast = get_fsalfirstlast(cache.caches[1], u) integrator.fsalfirst = fsalfirst integrator.fsallast = fsallast initialize!(integrator, @inbounds(cache.caches[1])) elseif cache.current == 2 - fsalfirst, fsallast = get_fsalfirstlast(cache.caches[2],u) + fsalfirst, fsallast = get_fsalfirstlast(cache.caches[2], u) integrator.fsalfirst = fsalfirst integrator.fsallast = fsallast initialize!(integrator, @inbounds(cache.caches[2])) @@ -172,12 +172,12 @@ function choose_algorithm!(integrator, @inbounds if new_current != old_current cache.current = new_current if new_current == 1 - fsalfirst, fsallast = get_fsalfirstlast(cache.caches[1],u) + fsalfirst, fsallast = get_fsalfirstlast(cache.caches[1], u) integrator.fsalfirst = fsalfirst integrator.fsallast = fsallast initialize!(integrator, @inbounds(cache.caches[1])) elseif new_current == 2 - fsalfirst, fsallast = get_fsalfirstlast(cache.caches[2],u) + fsalfirst, fsallast = get_fsalfirstlast(cache.caches[2], u) integrator.fsalfirst = fsalfirst integrator.fsallast = fsallast initialize!(integrator, @inbounds(cache.caches[2])) @@ -205,37 +205,37 @@ function choose_algorithm!(integrator, cache::DefaultCache) cache.current = new_current init_ith_default_cache(cache, algs, new_current) if new_current == 1 - fsalfirst, fsallast = get_fsalfirstlast(cache.cache1,u) + fsalfirst, fsallast = get_fsalfirstlast(cache.cache1, u) integrator.fsalfirst = fsalfirst integrator.fsallast = fsallast initialize!(integrator, @inbounds(cache.cache1)) new_cache = cache.cache1 elseif new_current == 2 - fsalfirst, fsallast = get_fsalfirstlast(cache.cache2,u) + fsalfirst, fsallast = get_fsalfirstlast(cache.cache2, u) integrator.fsalfirst = fsalfirst integrator.fsallast = fsallast initialize!(integrator, @inbounds(cache.cache2)) new_cache = cache.cache2 elseif new_current == 3 - fsalfirst, fsallast = get_fsalfirstlast(cache.cache3,u) + fsalfirst, fsallast = get_fsalfirstlast(cache.cache3, u) integrator.fsalfirst = fsalfirst integrator.fsallast = fsallast initialize!(integrator, @inbounds(cache.cache3)) new_cache = cache.cache3 elseif new_current == 4 - fsalfirst, fsallast = get_fsalfirstlast(cache.cache4,u) + fsalfirst, fsallast = get_fsalfirstlast(cache.cache4, u) integrator.fsalfirst = fsalfirst integrator.fsallast = fsallast initialize!(integrator, @inbounds(cache.cache4)) new_cache = cache.cache4 elseif new_current == 5 - fsalfirst, fsallast = get_fsalfirstlast(cache.cache5,u) + fsalfirst, fsallast = get_fsalfirstlast(cache.cache5, u) integrator.fsalfirst = fsalfirst integrator.fsallast = fsallast initialize!(integrator, @inbounds(cache.cache5)) new_cache = cache.cache5 elseif new_current == 6 - fsalfirst, fsallast = get_fsalfirstlast(cache.cache6,u) + fsalfirst, fsallast = get_fsalfirstlast(cache.cache6, u) integrator.fsalfirst = fsalfirst integrator.fsallast = fsallast initialize!(integrator, @inbounds(cache.cache6)) diff --git a/lib/OrdinaryDiffEqCore/src/precompilation_setup.jl b/lib/OrdinaryDiffEqCore/src/precompilation_setup.jl index 655db1737c..bfb61ff69a 100644 --- a/lib/OrdinaryDiffEqCore/src/precompilation_setup.jl +++ b/lib/OrdinaryDiffEqCore/src/precompilation_setup.jl @@ -12,16 +12,16 @@ PrecompileTools.@compile_workload begin ODEProblem(lorenz, [1.0; 0.0; 0.0], (0.0, 1.0)) ODEProblem(lorenz, [1.0; 0.0; 0.0], (0.0, 1.0), Float64[]) ODEProblem{true, SciMLBase.AutoSpecialize}(lorenz, [1.0; 0.0; 0.0], - (0.0, 1.0)) + (0.0, 1.0)) ODEProblem{true, SciMLBase.AutoSpecialize}(lorenz, [1.0; 0.0; 0.0], - (0.0, 1.0), Float64[]) + (0.0, 1.0), Float64[]) ODEProblem{true, SciMLBase.FunctionWrapperSpecialize}(lorenz, [1.0; 0.0; 0.0], - (0.0, 1.0)) + (0.0, 1.0)) ODEProblem{true, SciMLBase.FunctionWrapperSpecialize}(lorenz, [1.0; 0.0; 0.0], - (0.0, 1.0), Float64[]) + (0.0, 1.0), Float64[]) ODEProblem{true, SciMLBase.NoSpecialize}(lorenz, [1.0; 0.0; 0.0], (0.0, 1.0)) ODEProblem{true, SciMLBase.NoSpecialize}(lorenz, [1.0; 0.0; 0.0], (0.0, 1.0), - Float64[]) + Float64[]) lorenz([1.0; 0.0; 0.0], [1.0; 0.0; 0.0], DiffEqBase.NullParameters(), 0.0) lorenz([1.0; 0.0; 0.0], [1.0; 0.0; 0.0], Float64[], 0.0) diff --git a/lib/OrdinaryDiffEqCore/src/solve.jl b/lib/OrdinaryDiffEqCore/src/solve.jl index cef15ee62b..563c7654db 100644 --- a/lib/OrdinaryDiffEqCore/src/solve.jl +++ b/lib/OrdinaryDiffEqCore/src/solve.jl @@ -469,7 +469,7 @@ function DiffEqBase.__init( reinitiailize = true saveiter = 0 # Starts at 0 so first save is at 1 saveiter_dense = 0 - faslfirst, fsallast = get_fsalfirstlast(cache,rate_prototype) + faslfirst, fsallast = get_fsalfirstlast(cache, rate_prototype) integrator = ODEIntegrator{typeof(_alg), isinplace(prob), uType, typeof(du), tType, typeof(p), diff --git a/lib/OrdinaryDiffEqDifferentiation/src/alg_utils.jl b/lib/OrdinaryDiffEqDifferentiation/src/alg_utils.jl index 50e605e409..cfc14219c5 100644 --- a/lib/OrdinaryDiffEqDifferentiation/src/alg_utils.jl +++ b/lib/OrdinaryDiffEqDifferentiation/src/alg_utils.jl @@ -76,4 +76,4 @@ end @generated function pick_static_chunksize(::Val{chunksize}) where {chunksize} x = ForwardDiff.pickchunksize(chunksize) :(Val{$x}()) -end +end \ No newline at end of file diff --git a/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl b/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl index d16df97fa7..0914b355b0 100644 --- a/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl +++ b/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl @@ -12,7 +12,7 @@ struct StaticWOperator{isinv, T, F} <: AbstractSciMLOperator{T} # doing to how StaticArrays and StaticArraysCore are split up StaticArrays.LU(LowerTriangular(W), UpperTriangular(W), SVector{n}(1:n)) else - lu(W, check=false) + lu(W, check = false) end # when constructing W for the first time for the type # inv(W) can be singular @@ -932,28 +932,28 @@ function LinearSolve.init_cacheval( end for alg in [LinearSolve.AppleAccelerateLUFactorization, - LinearSolve.BunchKaufmanFactorization, - LinearSolve.CHOLMODFactorization, - LinearSolve.CholeskyFactorization, - LinearSolve.CudaOffloadFactorization, - LinearSolve.DiagonalFactorization, - LinearSolve.FastLUFactorization, - LinearSolve.FastQRFactorization, - LinearSolve.GenericFactorization, - LinearSolve.GenericLUFactorization, - LinearSolve.KLUFactorization, - LinearSolve.LDLtFactorization, - LinearSolve.LUFactorization, - LinearSolve.MKLLUFactorization, - LinearSolve.MetalLUFactorization, - LinearSolve.NormalBunchKaufmanFactorization, - LinearSolve.NormalCholeskyFactorization, - LinearSolve.QRFactorization, - LinearSolve.RFLUFactorization, - LinearSolve.SVDFactorization, - LinearSolve.SimpleLUFactorization, - LinearSolve.SparspakFactorization, - LinearSolve.UMFPACKFactorization] + LinearSolve.BunchKaufmanFactorization, + LinearSolve.CHOLMODFactorization, + LinearSolve.CholeskyFactorization, + LinearSolve.CudaOffloadFactorization, + LinearSolve.DiagonalFactorization, + LinearSolve.FastLUFactorization, + LinearSolve.FastQRFactorization, + LinearSolve.GenericFactorization, + LinearSolve.GenericLUFactorization, + LinearSolve.KLUFactorization, + LinearSolve.LDLtFactorization, + LinearSolve.LUFactorization, + LinearSolve.MKLLUFactorization, + LinearSolve.MetalLUFactorization, + LinearSolve.NormalBunchKaufmanFactorization, + LinearSolve.NormalCholeskyFactorization, + LinearSolve.QRFactorization, + LinearSolve.RFLUFactorization, + LinearSolve.SVDFactorization, + LinearSolve.SimpleLUFactorization, + LinearSolve.SparspakFactorization, + LinearSolve.UMFPACKFactorization] @eval function LinearSolve.init_cacheval(alg::$alg, A::WOperator, b, u, Pl, Pr, maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) @@ -997,4 +997,4 @@ function resize_J_W!(cache, integrator, i) end nothing -end +end \ No newline at end of file diff --git a/lib/OrdinaryDiffEqExplicitRK/src/explicit_rk_caches.jl b/lib/OrdinaryDiffEqExplicitRK/src/explicit_rk_caches.jl index 853e216b06..943bbb4d64 100644 --- a/lib/OrdinaryDiffEqExplicitRK/src/explicit_rk_caches.jl +++ b/lib/OrdinaryDiffEqExplicitRK/src/explicit_rk_caches.jl @@ -11,7 +11,7 @@ tab::TabType end -get_fsalfirstlast(cache::ExplicitRKCache,u) = (cache.kk[1], cache.fsallast) +get_fsalfirstlast(cache::ExplicitRKCache, u) = (cache.kk[1], cache.fsallast) function alg_cache(alg::ExplicitRK, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, diff --git a/lib/OrdinaryDiffEqExponentialRK/src/exponential_rk_caches.jl b/lib/OrdinaryDiffEqExponentialRK/src/exponential_rk_caches.jl index 65b662f825..6a2217e3e0 100644 --- a/lib/OrdinaryDiffEqExponentialRK/src/exponential_rk_caches.jl +++ b/lib/OrdinaryDiffEqExponentialRK/src/exponential_rk_caches.jl @@ -2,7 +2,7 @@ # Classical ExpRK method caches abstract type ExpRKCache <: OrdinaryDiffEqMutableCache end abstract type ExpRKConstantCache <: OrdinaryDiffEqConstantCache end -get_fsalfirstlast(cache::ExpRKCache,u) = (zero(cache.rtmp), zero(cache.rtmp)) +get_fsalfirstlast(cache::ExpRKCache, u) = (zero(cache.rtmp), zero(cache.rtmp)) # Precomputation of exponential-like operators """ @@ -872,7 +872,7 @@ struct ETD2ConstantCache{expType} <: OrdinaryDiffEqConstantCache B0::expType # -ϕ2(hA) end -get_fsalfirstlast(cache::ETD2ConstantCache,u) = (ETD2Fsal(u), ETD2Fsal(u)) +get_fsalfirstlast(cache::ETD2ConstantCache, u) = (ETD2Fsal(u), ETD2Fsal(u)) function alg_cache(alg::ETD2, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, @@ -894,7 +894,7 @@ end B1::expType # ϕ1(hA) + ϕ2(hA) B0::expType # -ϕ2(hA) end -get_fsalfirstlast(cache::ETD2Cache,u) = (ETD2Fsal(cache.rtmp1), ETD2Fsal(cache.rtmp1)) +get_fsalfirstlast(cache::ETD2Cache, u) = (ETD2Fsal(cache.rtmp1), ETD2Fsal(cache.rtmp1)) function alg_cache(alg::ETD2, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, diff --git a/lib/OrdinaryDiffEqExtrapolation/src/extrapolation_caches.jl b/lib/OrdinaryDiffEqExtrapolation/src/extrapolation_caches.jl index 7ede4e5b1a..30c8fc16c9 100644 --- a/lib/OrdinaryDiffEqExtrapolation/src/extrapolation_caches.jl +++ b/lib/OrdinaryDiffEqExtrapolation/src/extrapolation_caches.jl @@ -1,5 +1,5 @@ abstract type ExtrapolationMutableCache <: OrdinaryDiffEqMutableCache end -get_fsalfirstlast(cache::ExtrapolationMutableCache,u) = (cache.fsalfirst, cache.k) +get_fsalfirstlast(cache::ExtrapolationMutableCache, u) = (cache.fsalfirst, cache.k) @cache mutable struct AitkenNevilleCache{ uType, @@ -128,7 +128,7 @@ end diff1::Array{uType, 1} diff2::Array{uType, 1} end -get_fsalfirstlast(cache::ImplicitEulerExtrapolationCache,u) = (zero(u), zero(u)) +get_fsalfirstlast(cache::ImplicitEulerExtrapolationCache, u) = (zero(u), zero(u)) @cache mutable struct ImplicitEulerExtrapolationConstantCache{QType, dtType, arrayType, TF, UF, sequenceType} <: diff --git a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl index 5d7ccba8c4..e0ed1c2eb3 100644 --- a/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl +++ b/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl @@ -1,5 +1,5 @@ abstract type FIRKMutableCache <: OrdinaryDiffEqMutableCache end -get_fsalfirstlast(cache::FIRKMutableCache,u) = (cache.fsalfirst, cache.k) +get_fsalfirstlast(cache::FIRKMutableCache, u) = (cache.fsalfirst, cache.k) mutable struct RadauIIA3ConstantCache{F, Tab, Tol, Dt, U, JType} <: OrdinaryDiffEqConstantCache diff --git a/lib/OrdinaryDiffEqFeagin/src/feagin_caches.jl b/lib/OrdinaryDiffEqFeagin/src/feagin_caches.jl index 927d4778b2..e350b35725 100644 --- a/lib/OrdinaryDiffEqFeagin/src/feagin_caches.jl +++ b/lib/OrdinaryDiffEqFeagin/src/feagin_caches.jl @@ -1,5 +1,5 @@ abstract type FeaginCache <: OrdinaryDiffEqMutableCache end -get_fsalfirstlast(cache::FeaginCache,u) = (cache.fsalfirst, cache.k) +get_fsalfirstlast(cache::FeaginCache, u) = (cache.fsalfirst, cache.k) @cache struct Feagin10Cache{uType, uNoUnitsType, rateType, TabType, StepLimiter} <: FeaginCache diff --git a/lib/OrdinaryDiffEqFeagin/src/feagin_rk_perform_step.jl b/lib/OrdinaryDiffEqFeagin/src/feagin_rk_perform_step.jl index 296f5ece4f..ddba6505b9 100644 --- a/lib/OrdinaryDiffEqFeagin/src/feagin_rk_perform_step.jl +++ b/lib/OrdinaryDiffEqFeagin/src/feagin_rk_perform_step.jl @@ -86,7 +86,6 @@ end end function initialize!(integrator, cache::Feagin10Cache) - integrator.kshortsize = 2 resize!(integrator.k, integrator.kshortsize) integrator.k[1] = integrator.fsalfirst @@ -417,7 +416,6 @@ end end function initialize!(integrator, cache::Feagin12Cache) - integrator.kshortsize = 2 resize!(integrator.k, integrator.kshortsize) integrator.k[1] = integrator.fsalfirst @@ -907,7 +905,6 @@ end end function initialize!(integrator, cache::Feagin14Cache) - integrator.kshortsize = 2 resize!(integrator.k, integrator.kshortsize) integrator.k[1] = integrator.fsalfirst diff --git a/lib/OrdinaryDiffEqFunctionMap/src/functionmap_caches.jl b/lib/OrdinaryDiffEqFunctionMap/src/functionmap_caches.jl index 401a55adeb..38d7afe1d8 100644 --- a/lib/OrdinaryDiffEqFunctionMap/src/functionmap_caches.jl +++ b/lib/OrdinaryDiffEqFunctionMap/src/functionmap_caches.jl @@ -3,7 +3,7 @@ uprev::uType tmp::rateType end -get_fsalfirstlast(cache::FunctionMapCache,u) = (cache.u, cache.uprev) +get_fsalfirstlast(cache::FunctionMapCache, u) = (cache.u, cache.uprev) function alg_cache(alg::FunctionMap, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, diff --git a/lib/OrdinaryDiffEqHighOrderRK/src/high_order_rk_caches.jl b/lib/OrdinaryDiffEqHighOrderRK/src/high_order_rk_caches.jl index 2f8a0f4244..3aed829fb2 100644 --- a/lib/OrdinaryDiffEqHighOrderRK/src/high_order_rk_caches.jl +++ b/lib/OrdinaryDiffEqHighOrderRK/src/high_order_rk_caches.jl @@ -1,5 +1,5 @@ abstract type HighOrderRKMutableCache <: OrdinaryDiffEqMutableCache end -get_fsalfirstlast(cache::HighOrderRKMutableCache,u) = (cache.fsalfirst, cache.k) +get_fsalfirstlast(cache::HighOrderRKMutableCache, u) = (cache.fsalfirst, cache.k) @cache struct TanYam7Cache{uType, rateType, uNoUnitsType, TabType, StageLimiter, StepLimiter, Thread} <: HighOrderRKMutableCache @@ -92,7 +92,7 @@ end step_limiter!::StepLimiter thread::Thread end -get_fsalfirstlast(cache::DP8Cache,u) = (cache.k1, cache.k13) +get_fsalfirstlast(cache::DP8Cache, u) = (cache.k1, cache.k13) function alg_cache(alg::DP8, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, diff --git a/lib/OrdinaryDiffEqIMEXMultistep/src/imex_multistep_caches.jl b/lib/OrdinaryDiffEqIMEXMultistep/src/imex_multistep_caches.jl index 20d33b9d1b..d30425975a 100644 --- a/lib/OrdinaryDiffEqIMEXMultistep/src/imex_multistep_caches.jl +++ b/lib/OrdinaryDiffEqIMEXMultistep/src/imex_multistep_caches.jl @@ -1,6 +1,8 @@ # IMEX Multistep methods abstract type IMEXMutableCache <: OrdinaryDiffEqMutableCache end -get_fsalfirstlast(cache::IMEXMutableCache,u) = (cache.fsalfirst, du_alias_or_new(cache.nlsolver, cache.fsalfirst)) +function get_fsalfirstlast(cache::IMEXMutableCache, u) + (cache.fsalfirst, du_alias_or_new(cache.nlsolver, cache.fsalfirst)) +end # CNAB2 diff --git a/lib/OrdinaryDiffEqLinear/src/linear_caches.jl b/lib/OrdinaryDiffEqLinear/src/linear_caches.jl index b18605f1fa..47a1d1bde1 100644 --- a/lib/OrdinaryDiffEqLinear/src/linear_caches.jl +++ b/lib/OrdinaryDiffEqLinear/src/linear_caches.jl @@ -1,5 +1,5 @@ abstract type LinearMutableCache <: OrdinaryDiffEqMutableCache end -get_fsalfirstlast(cache::LinearMutableCache,u) = (cache.fsalfirst, cache.k) +get_fsalfirstlast(cache::LinearMutableCache, u) = (cache.fsalfirst, cache.k) @cache struct MagnusMidpointCache{uType, rateType, WType, expType} <: LinearMutableCache @@ -565,7 +565,7 @@ end exp_cache::expType end -get_fsalfirstlast(cache::LinearExponentialCache,u) = (zero(u), zero(u)) +get_fsalfirstlast(cache::LinearExponentialCache, u) = (zero(u), zero(u)) function _phiv_timestep_caches(u_prototype, maxiter::Int, p::Int) n = length(u_prototype) diff --git a/lib/OrdinaryDiffEqLowOrderRK/Project.toml b/lib/OrdinaryDiffEqLowOrderRK/Project.toml index 737d7e9e1f..75b4acea91 100644 --- a/lib/OrdinaryDiffEqLowOrderRK/Project.toml +++ b/lib/OrdinaryDiffEqLowOrderRK/Project.toml @@ -1,7 +1,7 @@ name = "OrdinaryDiffEqLowOrderRK" uuid = "1344f307-1e59-4825-a18e-ace9aa3fa4c6" authors = ["ParamThakkar123 "] -version = "1.1.0" +version = "1.2.0" [deps] DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" diff --git a/lib/OrdinaryDiffEqLowOrderRK/src/fixed_timestep_perform_step.jl b/lib/OrdinaryDiffEqLowOrderRK/src/fixed_timestep_perform_step.jl index c4ab6f12a6..cfccbcb25a 100644 --- a/lib/OrdinaryDiffEqLowOrderRK/src/fixed_timestep_perform_step.jl +++ b/lib/OrdinaryDiffEqLowOrderRK/src/fixed_timestep_perform_step.jl @@ -21,7 +21,7 @@ function perform_step!(integrator, cache::EulerConstantCache, repeat_step = fals integrator.u = u end -get_fsalfirstlast(cache::EulerCache,u) = (cache.fsalfirst, cache.k) +get_fsalfirstlast(cache::EulerCache, u) = (cache.fsalfirst, cache.k) function initialize!(integrator, cache::EulerCache) integrator.kshortsize = 2 @unpack k, fsalfirst = cache @@ -97,7 +97,7 @@ end integrator.u = u end -get_fsalfirstlast(cache::Union{HeunCache, RalstonCache},u) = (cache.fsalfirst, cache.k) +get_fsalfirstlast(cache::Union{HeunCache, RalstonCache}, u) = (cache.fsalfirst, cache.k) function initialize!(integrator, cache::Union{HeunCache, RalstonCache}) integrator.kshortsize = 2 @unpack k, fsalfirst = cache @@ -189,7 +189,7 @@ end integrator.u = u end -get_fsalfirstlast(cache::MidpointCache,u) = (cache.fsalfirst, cache.k) +get_fsalfirstlast(cache::MidpointCache, u) = (cache.fsalfirst, cache.k) function initialize!(integrator, cache::MidpointCache) @unpack k, fsalfirst = cache integrator.fsalfirst = fsalfirst @@ -288,7 +288,7 @@ end integrator.u = u end -get_fsalfirstlast(cache::RK4Cache,u) = (cache.fsalfirst, cache.k) +get_fsalfirstlast(cache::RK4Cache, u) = (cache.fsalfirst, cache.k) function initialize!(integrator, cache::RK4Cache) @unpack tmp, fsalfirst, k₂, k₃, k₄, k = cache integrator.fsalfirst = fsalfirst @@ -420,7 +420,7 @@ end integrator.u = u end -get_fsalfirstlast(cache::Anas5Cache,u) = (cache.k1, cache.k7) +get_fsalfirstlast(cache::Anas5Cache, u) = (cache.k1, cache.k7) function initialize!(integrator, cache::Anas5Cache) integrator.kshortsize = 7 resize!(integrator.k, integrator.kshortsize) diff --git a/lib/OrdinaryDiffEqLowOrderRK/src/low_order_rk_perform_step.jl b/lib/OrdinaryDiffEqLowOrderRK/src/low_order_rk_perform_step.jl index b7f8bdd7f0..9d9928f678 100644 --- a/lib/OrdinaryDiffEqLowOrderRK/src/low_order_rk_perform_step.jl +++ b/lib/OrdinaryDiffEqLowOrderRK/src/low_order_rk_perform_step.jl @@ -33,7 +33,7 @@ end integrator.u = u end -get_fsalfirstlast(cache::BS3Cache,u) = (cache.fsalfirst, cache.k4) +get_fsalfirstlast(cache::BS3Cache, u) = (cache.fsalfirst, cache.k4) function initialize!(integrator, cache::BS3Cache) integrator.kshortsize = 2 resize!(integrator.k, integrator.kshortsize) @@ -113,7 +113,7 @@ end integrator.u = u end -get_fsalfirstlast(cache::OwrenZen3Cache,u) = (cache.k1, cache.k4) +get_fsalfirstlast(cache::OwrenZen3Cache, u) = (cache.k1, cache.k4) function initialize!(integrator, cache::OwrenZen3Cache) integrator.kshortsize = 4 resize!(integrator.k, integrator.kshortsize) @@ -195,7 +195,7 @@ end integrator.u = u end -get_fsalfirstlast(cache::OwrenZen4Cache,u) = (cache.k1, cache.k6) +get_fsalfirstlast(cache::OwrenZen4Cache, u) = (cache.k1, cache.k6) function initialize!(integrator, cache::OwrenZen4Cache) integrator.kshortsize = 6 resize!(integrator.k, integrator.kshortsize) @@ -297,7 +297,7 @@ end integrator.u = u end -get_fsalfirstlast(cache::OwrenZen5Cache,u) = (cache.k1, cache.k8) +get_fsalfirstlast(cache::OwrenZen5Cache, u) = (cache.k1, cache.k8) function initialize!(integrator, cache::OwrenZen5Cache) integrator.kshortsize = 8 resize!(integrator.k, integrator.kshortsize) @@ -454,7 +454,7 @@ end end end -get_fsalfirstlast(cache::BS5Cache,u) = (cache.k1, cache.k8) +get_fsalfirstlast(cache::BS5Cache, u) = (cache.k1, cache.k8) function initialize!(integrator, cache::BS5Cache) alg = unwrap_alg(integrator, false) alg.lazy ? (integrator.kshortsize = 8) : (integrator.kshortsize = 11) @@ -619,7 +619,7 @@ end integrator.u = u end -get_fsalfirstlast(cache::DP5Cache,u) = (cache.k1, cache.k7) +get_fsalfirstlast(cache::DP5Cache, u) = (cache.k1, cache.k7) function initialize!(integrator, cache::DP5Cache) integrator.kshortsize = 4 resize!(integrator.k, integrator.kshortsize) @@ -736,7 +736,7 @@ end integrator.u = u end -get_fsalfirstlast(cache::RKO65Cache,u) = (cache.k1, cache.k6) +get_fsalfirstlast(cache::RKO65Cache, u) = (cache.k1, cache.k6) function initialize!(integrator, cache::RKO65Cache) @unpack k, fsalfirst = cache integrator.kshortsize = 6 @@ -863,7 +863,7 @@ end integrator.u = u end -get_fsalfirstlast(cache::FRK65Cache,u) = (cache.k1, cache.k9) +get_fsalfirstlast(cache::FRK65Cache, u) = (cache.k1, cache.k9) function initialize!(integrator, cache::FRK65Cache) integrator.kshortsize = 9 @@ -991,7 +991,7 @@ end integrator.u = u end -get_fsalfirstlast(cache::RKMCache,u) = (cache.k1, zero(cache.k1)) +get_fsalfirstlast(cache::RKMCache, u) = (cache.k1, zero(cache.k1)) function initialize!(integrator, cache::RKMCache) @unpack k, fsalfirst = cache integrator.kshortsize = 6 @@ -1081,7 +1081,7 @@ function perform_step!(integrator, cache::PSRK4p7q6ConstantCache, repeat_step = integrator.u = u end -get_fsalfirstlast(cache::PSRK4p7q6Cache,u) = (cache.k1, cache.k6) +get_fsalfirstlast(cache::PSRK4p7q6Cache, u) = (cache.k1, cache.k6) function initialize!(integrator, cache::PSRK4p7q6Cache) @unpack uprev, f, p, t = integrator @@ -1169,7 +1169,7 @@ function perform_step!(integrator, cache::PSRK3p6q5ConstantCache, repeat_step = integrator.u = u end -get_fsalfirstlast(cache::PSRK3p6q5Cache,u) = (cache.k1, cache.k5) +get_fsalfirstlast(cache::PSRK3p6q5Cache, u) = (cache.k1, cache.k5) function initialize!(integrator, cache::PSRK3p6q5Cache) @unpack uprev, f, p, t = integrator @@ -1246,7 +1246,7 @@ function perform_step!(integrator, cache::PSRK3p5q4ConstantCache, repeat_step = integrator.u = u end -get_fsalfirstlast(cache::PSRK3p5q4Cache,u) = (cache.k1, cache.k4) +get_fsalfirstlast(cache::PSRK3p5q4Cache, u) = (cache.k1, cache.k4) function initialize!(integrator, cache::PSRK3p5q4Cache) @unpack uprev, f, p, t = integrator @@ -1335,7 +1335,7 @@ function perform_step!(integrator, cache::MSRK5ConstantCache, repeat_step = fals integrator.u = u end -get_fsalfirstlast(cache::MSRK5Cache,u) = (cache.k1, cache.k9) +get_fsalfirstlast(cache::MSRK5Cache, u) = (cache.k1, cache.k9) function initialize!(integrator, cache::MSRK5Cache) @unpack uprev, f, p, t = integrator @@ -1452,7 +1452,7 @@ function perform_step!(integrator, cache::MSRK6ConstantCache, repeat_step = fals integrator.u = u end -get_fsalfirstlast(cache::MSRK6Cache,u) = (cache.k1, cache.k9) +get_fsalfirstlast(cache::MSRK6Cache, u) = (cache.k1, cache.k9) function initialize!(integrator, cache::MSRK6Cache) @unpack uprev, f, p, t = integrator @@ -1572,7 +1572,7 @@ function perform_step!(integrator, cache::Stepanov5ConstantCache, repeat_step = integrator.u = u end -get_fsalfirstlast(cache::Stepanov5Cache,u) = (cache.k1, cache.k7) +get_fsalfirstlast(cache::Stepanov5Cache, u) = (cache.k1, cache.k7) function initialize!(integrator, cache::Stepanov5Cache) @unpack uprev, f, p, t = integrator @@ -1693,7 +1693,7 @@ function perform_step!(integrator, cache::SIR54ConstantCache, repeat_step = fals integrator.u = u end -get_fsalfirstlast(cache::SIR54Cache,u) = (cache.k1, cache.k8) +get_fsalfirstlast(cache::SIR54Cache, u) = (cache.k1, cache.k8) function initialize!(integrator, cache::SIR54Cache) @unpack uprev, f, p, t = integrator @@ -1801,7 +1801,7 @@ function perform_step!(integrator, cache::Alshina2ConstantCache, repeat_step = f integrator.u = u end -get_fsalfirstlast(cache::Alshina2Cache,u) = (cache.k1, cache.k2) +get_fsalfirstlast(cache::Alshina2Cache, u) = (cache.k1, cache.k2) function initialize!(integrator, cache::Alshina2Cache) @unpack uprev, f, p, t = integrator @@ -1882,7 +1882,7 @@ function perform_step!(integrator, cache::Alshina3ConstantCache, repeat_step = f integrator.u = u end -get_fsalfirstlast(cache::Alshina3Cache,u) = (cache.k1, cache.k3) +get_fsalfirstlast(cache::Alshina3Cache, u) = (cache.k1, cache.k3) function initialize!(integrator, cache::Alshina3Cache) @unpack uprev, f, p, t = integrator @@ -1978,7 +1978,7 @@ function perform_step!(integrator, cache::Alshina6ConstantCache, repeat_step = f integrator.u = u end -get_fsalfirstlast(cache::Alshina6Cache,u) = (cache.k1, cache.k7) +get_fsalfirstlast(cache::Alshina6Cache, u) = (cache.k1, cache.k7) function initialize!(integrator, cache::Alshina6Cache) @unpack uprev, f, p, t = integrator diff --git a/lib/OrdinaryDiffEqLowOrderRK/src/split_perform_step.jl b/lib/OrdinaryDiffEqLowOrderRK/src/split_perform_step.jl index 997610988f..aeb8511f0c 100644 --- a/lib/OrdinaryDiffEqLowOrderRK/src/split_perform_step.jl +++ b/lib/OrdinaryDiffEqLowOrderRK/src/split_perform_step.jl @@ -24,7 +24,7 @@ end integrator.u = u end -get_fsalfirstlast(cache::SplitEulerCache,u) = (cache.fsalfirst, cache.k) +get_fsalfirstlast(cache::SplitEulerCache, u) = (cache.fsalfirst, cache.k) function initialize!(integrator, cache::SplitEulerCache) integrator.kshortsize = 2 @unpack k, fsalfirst = cache diff --git a/lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl b/lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl index 2e7f9fb8e2..867d79097d 100644 --- a/lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl +++ b/lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl @@ -103,6 +103,6 @@ export ORK256, CarpenterKennedy2N54, SHLDDRK64, HSLDDRK64, DGLDDRK73_C, DGLDDRK8 ParsaniKetchesonDeconinck3S53, ParsaniKetchesonDeconinck3S173, ParsaniKetchesonDeconinck3S94, ParsaniKetchesonDeconinck3S184, ParsaniKetchesonDeconinck3S105, ParsaniKetchesonDeconinck3S205, - RDPK3Sp35, RDPK3SpFSAL35, RDPK3Sp49, RDPK3SpFSAL49, RDPK3Sp510, RDPK3SpFSAL510, - KYK2014DGSSPRK_3S2, RK46NL + RDPK3Sp35, RDPK3SpFSAL35, RDPK3Sp49, RDPK3SpFSAL49, RDPK3Sp510, RDPK3SpFSAL510, + RK46NL, SHLDDRK_2N, SHLDDRK52 end diff --git a/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl b/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl index 787fa9ba00..1dc1a9e81a 100644 --- a/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl +++ b/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl @@ -1,4 +1,3 @@ -alg_order(alg::KYK2014DGSSPRK_3S2) = 2 alg_order(alg::ORK256) = 2 alg_order(alg::CarpenterKennedy2N54) = 4 alg_order(alg::NDBLSRK124) = 4 @@ -41,6 +40,8 @@ alg_order(alg::CKLLSRK54_3N_4R) = 4 alg_order(alg::CKLLSRK75_4M_5R) = 5 alg_order(alg::CKLLSRK54_3M_4R) = 4 alg_order(alg::CKLLSRK65_4M_4R) = 5 +alg_order(alg::SHLDDRK_2N) = 4 +alg_order(alg::SHLDDRK52) = 2 isfsal(alg::ORK256) = false isfsal(alg::CarpenterKennedy2N54) = false @@ -83,22 +84,6 @@ uses_uprev(alg::CKLLSRK65_4M_4R, adaptive::Bool) = adaptive uses_uprev(alg::CKLLSRK85_4FM_4R, adaptive::Bool) = adaptive uses_uprev(alg::CKLLSRK75_4M_5R, adaptive::Bool) = adaptive -""" - ssp_coefficient(alg) - -Return the SSP coefficient of the ODE algorithm `alg`. If one time step of size -`dt` with `alg` can be written as a convex combination of explicit Euler steps -with step sizes `cᵢ * dt`, the SSP coefficient is the minimal value of `1/cᵢ`. - -# Examples - -```julia-repl -julia> ssp_coefficient(SSPRK104()) -6 -``` -""" -ssp_coefficient(alg::KYK2014DGSSPRK_3S2) = 0.8417 - function default_controller(alg::RDPK3Sp35, cache, qoldinit, args...) QT = typeof(qoldinit) return PIDController(map(Base.Fix1(convert, QT), (0.64, -0.31, 0.04))...) diff --git a/lib/OrdinaryDiffEqLowStorageRK/src/algorithms.jl b/lib/OrdinaryDiffEqLowStorageRK/src/algorithms.jl index 58b3b3808f..d0f69f7711 100644 --- a/lib/OrdinaryDiffEqLowStorageRK/src/algorithms.jl +++ b/lib/OrdinaryDiffEqLowStorageRK/src/algorithms.jl @@ -51,6 +51,32 @@ function DGLDDRK73_C(stage_limiter!, step_limiter! = trivial_limiter!; williamson_condition) end +@doc explicit_rk_docstring("TBD", "SHLDDRK_2N") +Base.@kwdef struct SHLDDRK_2N{StageLimiter, StepLimiter, Thread} <: OrdinaryDiffEqAlgorithm + stage_limiter!::StageLimiter = trivial_limiter! + step_limiter!::StepLimiter = trivial_limiter! + thread::Thread = False() +end +# for backwards compatibility +function SHLDDRK_2N(stage_limiter!, step_limiter! = trivial_limiter!) + SHLDDRK_2N(stage_limiter!, + step_limiter!, + False()) +end + +@doc explicit_rk_docstring("TBD", "SHLDDRK52") +Base.@kwdef struct SHLDDRK52{StageLimiter, StepLimiter, Thread} <: OrdinaryDiffEqAlgorithm + stage_limiter!::StageLimiter = trivial_limiter! + step_limiter!::StepLimiter = trivial_limiter! + thread::Thread = False() +end +# for backwards compatibility +function SHLDDRK52(stage_limiter!, step_limiter! = trivial_limiter!) + SHLDDRK52(stage_limiter!, + step_limiter!, + False()) +end + @doc explicit_rk_docstring( "A fourth-order, five-stage explicit low-storage method of Carpenter and Kennedy (free 3rd order Hermite interpolant). Fixed timestep only. Designed for @@ -888,19 +914,4 @@ function NDBLSRK134(stage_limiter!, step_limiter! = trivial_limiter!; williamson_condition) end -#SSP Optimized Runge-Kutta Methods - -@doc explicit_rk_docstring("TBD", - "KYK2014DGSSPRK_3S2") -Base.@kwdef struct KYK2014DGSSPRK_3S2{StageLimiter, StepLimiter, Thread} <: - OrdinaryDiffEqAlgorithm - stage_limiter!::StageLimiter = trivial_limiter! - step_limiter!::StepLimiter = trivial_limiter! - thread::Thread = False() -end -# for backwards compatibility -function KYK2014DGSSPRK_3S2(stage_limiter!, step_limiter! = trivial_limiter!) - KYK2014DGSSPRK_3S2(stage_limiter!, - step_limiter!, - False()) -end +#SSP Optimized Runge-Kutta Methods \ No newline at end of file diff --git a/lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_caches.jl b/lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_caches.jl index c2f4ca4623..6c813472d3 100644 --- a/lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_caches.jl +++ b/lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_caches.jl @@ -1,5 +1,5 @@ abstract type LowStorageRKMutableCache <: OrdinaryDiffEqMutableCache end -get_fsalfirstlast(cache::LowStorageRKMutableCache,u) = (cache.fsalfirst, cache.k) +get_fsalfirstlast(cache::LowStorageRKMutableCache, u) = (cache.fsalfirst, cache.k) # 2N low storage methods introduced by Williamson @cache struct LowStorageRK2NCache{uType, rateType, TabType, StageLimiter, StepLimiter, @@ -119,83 +119,6 @@ struct RK46NLConstantCache{T, T2} <: OrdinaryDiffEqConstantCache end end -@cache struct KYK2014DGSSPRK_3S2_Cache{uType, rateType, TabType, StageLimiter, StepLimiter, - Thread} <: - LowStorageRKMutableCache - u::uType - uprev::uType - k::rateType - fsalfirst::rateType - tab::TabType - #temporary values for Shu-Osher - u_1::uType - u_2::uType - kk_1::rateType - kk_2::rateType - stage_limiter!::StageLimiter - step_limiter!::StepLimiter - thread::Thread -end - -struct KYK2014DGSSPRK_3S2_ConstantCache{T, T2} <: OrdinaryDiffEqConstantCache - #These are not α and β for RK but for Shu-Osher - #see top of page 317 in - #Optimal Strong-Stability-Preserving Runge–Kutta Time Discretizations for - #Discontinuous Garlekin Methods, Kubatko, Yaeger, Ketcheson 2014 - α_10::T - α_20::T - α_21::T - α_30::T - α_32::T - β_10::T - β_21::T - β_30::T - β_32::T - #Shu-Osher is normally stated for autonomous systems, the times - #are calculated by hand for this scheme - c_1::T - c_2::T - - function KYK2014DGSSPRK_3S2_ConstantCache(T, T2) - α_10 = T(1.0) - α_20 = T(0.087353119859156) - α_21 = T(0.912646880140844) - α_30 = T(0.344956917166841) - α_32 = T(0.655043082833159) - β_10 = T(0.528005024856522) - β_21 = T(0.481882138633993) - β_30 = T(0.022826837460491) - β_32 = T(0.345866039233415) - c_1 = β_10 - c_2 = α_21 * β_10 + β_21 # ==0.96376427726 - new{T, T2}(α_10, α_20, α_21, α_30, α_32, β_10, β_21, β_30, β_32, c_1, c_2) - end -end - -function alg_cache(alg::KYK2014DGSSPRK_3S2, u, rate_prototype, ::Type{uEltypeNoUnits}, - ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, - dt, reltol, p, calck, - ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} - u_1 = zero(u) - u_2 = zero(u) - kk_1 = zero(rate_prototype) - kk_2 = zero(rate_prototype) - k = zero(rate_prototype) - fsalfirst = zero(rate_prototype) - tab = KYK2014DGSSPRK_3S2_ConstantCache(constvalue(uBottomEltypeNoUnits), - constvalue(tTypeNoUnits)) - KYK2014DGSSPRK_3S2_Cache(u, uprev, k, fsalfirst, tab, u_1, u_2, kk_1, kk_2, - alg.stage_limiter!, alg.step_limiter!, alg.thread) -end - -function alg_cache(alg::KYK2014DGSSPRK_3S2, u, rate_prototype, ::Type{uEltypeNoUnits}, - ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, - dt, reltol, p, calck, - ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} - KYK2014DGSSPRK_3S2_ConstantCache(constvalue(uBottomEltypeNoUnits), - constvalue(tTypeNoUnits)) -end - function alg_cache(alg::RK46NL, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck, @@ -244,6 +167,176 @@ function CarpenterKennedy2N54ConstantCache(T, T2) LowStorageRK2NConstantCache{4, T, T2}(A2end, B1, B2end, c2end) end +@cache mutable struct SHLDDRK_2NCache{uType, rateType, TabType, StageLimiter, StepLimiter, + Thread} <: LowStorageRKMutableCache + u::uType + uprev::uType + k::rateType + tmp::uType + fsalfirst::rateType + tab::TabType + step::Int + stage_limiter!::StageLimiter + step_limiter!::StepLimiter + thread::Thread +end + +mutable struct SHLDDRK_2NConstantCache{T1, T2} <: OrdinaryDiffEqConstantCache + α21::T1 + α31::T1 + α41::T1 + α51::T1 + β11::T1 + β21::T1 + β31::T1 + β41::T1 + β51::T1 + c21::T2 + c31::T2 + c41::T2 + c51::T2 + + α22::T1 + α32::T1 + α42::T1 + α52::T1 + α62::T1 + β12::T1 + β22::T1 + β32::T1 + β42::T1 + β52::T1 + β62::T1 + c22::T2 + c32::T2 + c42::T2 + c52::T2 + c62::T2 + + step::Int +end + +function SHLDDRK_2NConstantCache(T1, T2) + α21 = T1(-0.6051226) + α31 = T1(-2.0437564) + α41 = T1(-0.7406999) + α51 = T1(-4.4231765) + β11 = T1(0.2687454) + β21 = T1(0.8014706) + β31 = T1(0.5051570) + β41 = T1(0.5623568) + β51 = T1(0.0590065) + c21 = T2(0.2687454) + c31 = T2(0.5852280) + c41 = T2(0.6827066) + c51 = T2(1.1646854) + + α22 = T1(-0.4412737) + α32 = T1(-1.0739820) + α42 = T1(-1.7063570) + α52 = T1(-2.7979293) + α62 = T1(-4.0913537) + β12 = T1(0.1158488) + β22 = T1(0.3728769) + β32 = T1(0.7379536) + β42 = T1(0.5798110) + β52 = T1(1.0312849) + β62 = T1(0.15) + c22 = T2(0.1158485) + c32 = T2(0.3241850) + c42 = T2(0.6193208) + c52 = T2(0.8034472) + c62 = T2(0.9184166) + SHLDDRK_2NConstantCache( + α21, α31, α41, α51, β11, β21, β31, β41, β51, c21, c31, c41, c51, + α22, α32, α42, α52, α62, β12, β22, β32, β42, β52, β62, c22, c32, + c42, c52, c62, 1) +end + +function alg_cache(alg::SHLDDRK_2N, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, + dt, reltol, p, calck, + ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + SHLDDRK_2NConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits)) +end + +function alg_cache(alg::SHLDDRK_2N, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, + dt, reltol, p, calck, + ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + tmp = zero(u) + k = zero(rate_prototype) + fsalfirst = zero(rate_prototype) + tab = SHLDDRK_2NConstantCache(constvalue(uBottomEltypeNoUnits), + constvalue(tTypeNoUnits)) + SHLDDRK_2NCache(u, uprev, k, tmp, fsalfirst, tab, 1, alg.stage_limiter!, + alg.step_limiter!, alg.thread) +end + +@cache struct SHLDDRK52Cache{uType, rateType, TabType, StageLimiter, StepLimiter, Thread} <: LowStorageRKMutableCache + u::uType + uprev::uType + k::rateType + tmp::uType + fsalfirst::rateType + tab::TabType + stage_limiter!::StageLimiter + step_limiter!::StepLimiter + thread::Thread +end + +struct SHLDDRK52ConstantCache{T1, T2} <: OrdinaryDiffEqConstantCache + α2::T1 + α3::T1 + α4::T1 + α5::T1 + β1::T1 + β2::T1 + β3::T1 + β4::T1 + β5::T1 + c2::T2 + c3::T2 + c4::T2 + c5::T2 +end + +function SHLDDRK52ConstantCache(T1, T2) + α2 = T1(-0.6913065) + α3 = T1(-2.655155) + α4 = T1(-0.8147688) + α5 = T1(-0.6686587) + β1 = T1(0.1) + β2 = T1(0.75) + β3 = T1(0.7) + β4 = T1(0.479313) + β5 = T1(0.310392) + c2 = T2(0.1) + c3 = T2(0.3315201) + c4 = T2(0.4577796) + c5 = T2(0.8666528) + SHLDDRK52ConstantCache(α2, α3, α4, α5, β1, β2, β3, β4, β5, c2, c3, c4, c5) +end + +function alg_cache(alg::SHLDDRK52, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, + dt, reltol, p, calck, + ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + SHLDDRK52ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits)) +end + +function alg_cache(alg::SHLDDRK52, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, + dt, reltol, p, calck, + ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + tmp = zero(u) + k = zero(rate_prototype) + fsalfirst = zero(rate_prototype) + tab = SHLDDRK52ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits)) + SHLDDRK52Cache(u, uprev, k, tmp, fsalfirst, tab, alg.stage_limiter!, alg.step_limiter!, + alg.thread) +end + function alg_cache(alg::CarpenterKennedy2N54, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck, diff --git a/lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_perform_step.jl b/lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_perform_step.jl index f76fcb4eb6..8805ca5b2f 100644 --- a/lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_perform_step.jl +++ b/lib/OrdinaryDiffEqLowStorageRK/src/low_storage_rk_perform_step.jl @@ -33,7 +33,7 @@ end integrator.u = u end -get_fsalfirstlast(cache::LowStorageRK2NCache,u) = (cache.k, cache.k) +get_fsalfirstlast(cache::LowStorageRK2NCache, u) = (cache.k, cache.k) function initialize!(integrator, cache::LowStorageRK2NCache) @unpack k, tmp, williamson_condition = cache @@ -932,34 +932,46 @@ end integrator.u = u end -function initialize!(integrator, cache::KYK2014DGSSPRK_3S2_ConstantCache) - integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal - OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) +function initialize!(integrator, cache::SHLDDRK52ConstantCache) 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) - return nothing + integrator.k[1] = integrator.fsalfirst + integrator.k[2] = integrator.fsallast end -@muladd function perform_step!(integrator, cache::KYK2014DGSSPRK_3S2_ConstantCache, +@muladd function perform_step!(integrator, cache::SHLDDRK52ConstantCache, repeat_step = false) @unpack t, dt, uprev, u, f, p = integrator - @unpack α_10, α_20, α_21, α_30, α_32, β_10, β_21, β_30, β_32, c_1, c_2 = cache - u_1 = α_10 * uprev + dt * β_10 * integrator.fsalfirst - u_2 = (α_20 * uprev + - α_21 * u_1 + dt * β_21 * f(u_1, p, t + c_1 * dt)) - integrator.u = (α_30 * uprev + dt * β_30 * integrator.fsalfirst + - α_32 * u_2 + dt * β_32 * f(u_2, p, t + c_2 * dt)) + @unpack α2, α3, α4, α5, β1, β2, β3, β4, β5, c2, c3, c4, c5 = cache + + # u1 + tmp = dt * integrator.fsalfirst + u = uprev + β1 * tmp + # u2 + tmp = α2 * tmp + dt * f(u, p, t + c2 * dt) + u = u + β2 * tmp + # u3 + tmp = α3 * tmp + dt * f(u, p, t + c3 * dt) + u = u + β3 * tmp + # u4 + tmp = α4 * tmp + dt * f(u, p, t + c4 * dt) + u = u + β4 * tmp + # u5 = u + tmp = α5 * tmp + dt * f(u, p, t + c5 * dt) + u = u + β5 * tmp + + integrator.fsallast = f(u, p, t + dt) # For interpolation, then FSAL'd integrator.k[1] = integrator.fsalfirst - integrator.k[2] = f(integrator.u, p, t + dt) # For interpolation, then FSAL'd - OrdinaryDiffEqCore.increment_nf!(integrator.stats, 3) - integrator.fsallast = integrator.k[2] - return nothing + integrator.k[2] = integrator.fsallast + integrator.u = u end -function initialize!(integrator, cache::KYK2014DGSSPRK_3S2_Cache) +function initialize!(integrator, cache::SHLDDRK52Cache) @unpack k, fsalfirst = cache integrator.kshortsize = 2 @@ -968,28 +980,188 @@ function initialize!(integrator, cache::KYK2014DGSSPRK_3S2_Cache) integrator.k[2] = integrator.fsallast integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # FSAL for interpolation OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) - return nothing end -@muladd function perform_step!(integrator, cache::KYK2014DGSSPRK_3S2_Cache, - repeat_step = false) +@muladd function perform_step!(integrator, cache::SHLDDRK52Cache, repeat_step = false) @unpack t, dt, uprev, u, f, p = integrator - @unpack k, fsalfirst, u_1, u_2, kk_1, kk_2, stage_limiter!, step_limiter!, thread = cache - @unpack α_10, α_20, α_21, α_30, α_32, β_10, β_21, β_30, β_32, c_1, c_2 = cache.tab - - @.. broadcast=false thread=thread u_1=α_10 * uprev + dt * β_10 * integrator.fsalfirst - stage_limiter!(u_1, integrator, p, t + c_1 * dt) - f(kk_1, u_1, p, t + c_1 * dt) - @.. broadcast=false thread=thread u_2=(α_20 * uprev + - α_21 * u_1 + dt * β_21 * kk_1) - stage_limiter!(u_2, integrator, p, t + c_2 * dt) - f(kk_2, u_2, p, t + c_2 * dt) - @.. broadcast=false thread=thread u=(α_30 * uprev + - dt * β_30 * integrator.fsalfirst + - α_32 * u_2 + dt * β_32 * kk_2) + @unpack k, fsalfirst, tmp, stage_limiter!, step_limiter!, thread = cache + @unpack α2, α3, α4, α5, β1, β2, β3, β4, β5, c2, c3, c4, c5 = cache.tab + + # u1 + @.. thread=thread tmp=dt * fsalfirst + @.. thread=thread u=uprev + β1 * tmp + stage_limiter!(u, integrator, p, t + c2 * dt) + # u2 + f(k, u, p, t + c2 * dt) + @.. thread=thread tmp=α2 * tmp + dt * k + @.. thread=thread u=u + β2 * tmp + stage_limiter!(u, integrator, p, t + c3 * dt) + # u3 + f(k, u, p, t + c3 * dt) + @.. thread=thread tmp=α3 * tmp + dt * k + @.. thread=thread u=u + β3 * tmp + stage_limiter!(u, integrator, p, t + c4 * dt) + # u4 + f(k, u, p, t + c4 * dt) + @.. thread=thread tmp=α4 * tmp + dt * k + @.. thread=thread u=u + β4 * tmp + stage_limiter!(u, integrator, p, t + c5 * dt) + # u5 = u + f(k, u, p, t + c5 * dt) + @.. thread=thread tmp=α5 * tmp + dt * k + @.. thread=thread u=u + β5 * tmp stage_limiter!(u, integrator, p, t + dt) step_limiter!(u, integrator, p, t + dt) - f(integrator.k[2], u, p, t + dt) # For interpolation, then FSAL'd - OrdinaryDiffEqCore.increment_nf!(integrator.stats, 3) - return nothing + + f(k, u, p, t + dt) end + +function initialize!(integrator, cache::SHLDDRK_2NConstantCache) + 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 +end + +@muladd function perform_step!(integrator, cache::SHLDDRK_2NConstantCache, + repeat_step = false) + @unpack t, dt, uprev, u, f, p = integrator + @unpack α21, α31, α41, α51, β11, β21, β31, β41, β51, c21, c31, c41, c51, α22, α32, α42, α52, α62, β12, β22, β32, β42, β52, β62, c22, c32, c42, c52, c62 = cache + + if integrator.u_modified + cache.step = 1 + end + # cnt = cache.step + + if cache.step % 2 == 1 + cache.step += 1 + # u1 + tmp = dt * integrator.fsalfirst + u = uprev + β11 * tmp + # u2 + tmp = α21 * tmp + dt * f(u, p, t + c21 * dt) + u = u + β21 * tmp + # u3 + tmp = α31 * tmp + dt * f(u, p, t + c31 * dt) + u = u + β31 * tmp + # u4 + tmp = α41 * tmp + dt * f(u, p, t + c41 * dt) + u = u + β41 * tmp + # u5 = u + tmp = α51 * tmp + dt * f(u, p, t + c51 * dt) + u = u + β51 * tmp + + else + cache.step += 1 + # u1 + tmp = dt * integrator.fsalfirst + u = uprev + β12 * tmp + # u2 + tmp = α22 * tmp + dt * f(u, p, t + c22 * dt) + u = u + β22 * tmp + # u3 + tmp = α32 * tmp + dt * f(u, p, t + c32 * dt) + u = u + β32 * tmp + # u4 + tmp = α42 * tmp + dt * f(u, p, t + c42 * dt) + u = u + β42 * tmp + # u5 = u + tmp = α52 * tmp + dt * f(u, p, t + c52 * dt) + u = u + β52 * tmp + tmp = α62 * tmp + dt * f(u, p, t + c62 * dt) + u = u + β62 * tmp + end + + integrator.fsallast = f(u, p, t + dt) # For interpolation, then FSAL'd + integrator.k[1] = integrator.fsalfirst + integrator.k[2] = integrator.fsallast + integrator.u = u +end + +function initialize!(integrator, cache::SHLDDRK_2NCache) + @unpack k, fsalfirst = cache + + 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) # FSAL for interpolation +end + +@muladd function perform_step!(integrator, cache::SHLDDRK_2NCache, repeat_step = false) + @unpack t, dt, uprev, u, f, p = integrator + @unpack k, fsalfirst, tmp, stage_limiter!, step_limiter!, thread = cache + @unpack α21, α31, α41, α51, β11, β21, β31, β41, β51, c21, c31, c41, c51, α22, α32, α42, α52, α62, β12, β22, β32, β42, β52, β62, c22, c32, c42, c52, c62 = cache.tab + + if integrator.u_modified + cache.step = 1 + end + + if cache.step % 2 == 1 + # u1 + @.. thread=thread tmp=dt * fsalfirst + @.. thread=thread u=uprev + β11 * tmp + stage_limiter!(u, integrator, p, t + c21 * dt) + # u2 + f(k, u, p, t + c21 * dt) + @.. thread=thread tmp=α21 * tmp + dt * k + @.. thread=thread u=u + β21 * tmp + stage_limiter!(u, integrator, p, t + c31 * dt) + # u3 + f(k, u, p, t + c31 * dt) + @.. thread=thread tmp=α31 * tmp + dt * k + @.. thread=thread u=u + β31 * tmp + stage_limiter!(u, integrator, p, t + c41 * dt) + # u4 + f(k, u, p, t + c41 * dt) + @.. thread=thread tmp=α41 * tmp + dt * k + @.. thread=thread u=u + β41 * tmp + stage_limiter!(u, integrator, p, t + c51 * dt) + # u5 = u + f(k, u, p, t + c51 * dt) + @.. thread=thread tmp=α51 * tmp + dt * k + @.. thread=thread u=u + β51 * tmp + stage_limiter!(u, integrator, p, t + dt) + step_limiter!(u, integrator, p, t + dt) + + f(k, u, p, t + dt) + else + # u1 + @.. thread=thread tmp=dt * fsalfirst + @.. thread=thread u=uprev + β12 * tmp + stage_limiter!(u, integrator, p, t + c22 * dt) + # u2 + f(k, u, p, t + c22 * dt) + @.. thread=thread tmp=α22 * tmp + dt * k + @.. thread=thread u=u + β22 * tmp + stage_limiter!(u, integrator, p, t + c32 * dt) + # u3 + f(k, u, p, t + c32 * dt) + @.. thread=thread tmp=α32 * tmp + dt * k + @.. thread=thread u=u + β32 * tmp + stage_limiter!(u, integrator, p, t + c42 * dt) + # u4 + f(k, u, p, t + c42 * dt) + @.. thread=thread tmp=α42 * tmp + dt * k + @.. thread=thread u=u + β42 * tmp + stage_limiter!(u, integrator, p, t + c52 * dt) + # u5 = u + f(k, u, p, t + c52 * dt) + @.. thread=thread tmp=α52 * tmp + dt * k + @.. thread=thread u=u + β52 * tmp + stage_limiter!(u, integrator, p, t + c62 * dt) + # u6 = u + f(k, u, p, t + c62 * dt) + @.. thread=thread tmp=α62 * tmp + dt * k + @.. thread=thread u=u + β62 * tmp + stage_limiter!(u, integrator, p, t + dt) + step_limiter!(u, integrator, p, t + dt) + + f(k, u, p, t + dt) + end +end \ No newline at end of file diff --git a/lib/OrdinaryDiffEqLowStorageRK/test/ode_low_storage_rk_tests.jl b/lib/OrdinaryDiffEqLowStorageRK/test/ode_low_storage_rk_tests.jl index 109aedb3e8..2ddc6f4401 100644 --- a/lib/OrdinaryDiffEqLowStorageRK/test/ode_low_storage_rk_tests.jl +++ b/lib/OrdinaryDiffEqLowStorageRK/test/ode_low_storage_rk_tests.jl @@ -41,6 +41,40 @@ test_problems_nonlinear = [prob_ode_nonlinear, prob_ode_nonlinear_inplace] u0_large = rand(10^6) prob_ode_large = ODEProblem((du, u, p, t) -> du .= u, u0_large, (0.0, 1.0)) +println("SHLDDRK_2N") +dts_SHLDDRK_2N = (1 / 2) .^ (0:3) +alg = SHLDDRK_2N() +for prob in test_problems_only_time + sim = test_convergence(dts_SHLDDRK_2N, prob, alg) + @test sim.𝒪est[:final]≈4 atol=0.46 +end +for prob in test_problems_linear + sim = test_convergence(dts_SHLDDRK_2N, prob, alg) + @test sim.𝒪est[:final]≈4 atol=0.46 +end +for prob in test_problems_nonlinear + sim = test_convergence(dts_SHLDDRK_2N, prob, alg) + @test sim.𝒪est[:final]≈4 atol=1 + # due to unusual saturation towards high dts(0.5 and onwards) and + # saturation towards low dts due to less precision in the provided values of weights , tolerance is kept so high +end + +println("SHLDDRK52") +dts = 1 .// 2 .^ (8:-1:4) +alg = SHLDDRK52() +for prob in test_problems_only_time + sim = test_convergence(dts, prob, alg) + @test_broken sim.𝒪est[:final]≈OrdinaryDiffEqLowStorageRK.alg_order(alg) atol=testTol +end +for prob in test_problems_linear + sim = test_convergence(dts, prob, alg) + @test_broken sim.𝒪est[:final]≈OrdinaryDiffEqLowStorageRK.alg_order(alg) atol=testTol +end +for prob in test_problems_nonlinear + sim = test_convergence(dts, prob, alg) + @test_broken sim.𝒪est[:final]≈OrdinaryDiffEqLowStorageRK.alg_order(alg) atol=testTol +end + @testset "ORK256" begin alg = ORK256() alg2 = ORK256(; williamson_condition = false) diff --git a/lib/OrdinaryDiffEqNordsieck/src/nordsieck_caches.jl b/lib/OrdinaryDiffEqNordsieck/src/nordsieck_caches.jl index fb175b5e38..a9ecb45215 100644 --- a/lib/OrdinaryDiffEqNordsieck/src/nordsieck_caches.jl +++ b/lib/OrdinaryDiffEqNordsieck/src/nordsieck_caches.jl @@ -265,4 +265,6 @@ function alg_cache(alg::JVODE, u, rate_prototype, ::Type{uEltypeNoUnits}, dts, Δ, atmp, tsit5cache, 2, 1, 1, 2, η, η, η, η, η) end -get_fsalfirstlast(cache::Union{JVODECache, AN5Cache},u) = get_fsalfirstlast(cache.tsit5cache,u) \ No newline at end of file +function get_fsalfirstlast(cache::Union{JVODECache, AN5Cache}, u) + get_fsalfirstlast(cache.tsit5cache, u) +end diff --git a/lib/OrdinaryDiffEqPDIRK/src/pdirk_caches.jl b/lib/OrdinaryDiffEqPDIRK/src/pdirk_caches.jl index 0b1d86d995..2cbbbf828f 100644 --- a/lib/OrdinaryDiffEqPDIRK/src/pdirk_caches.jl +++ b/lib/OrdinaryDiffEqPDIRK/src/pdirk_caches.jl @@ -8,7 +8,7 @@ end # Non-FSAL -get_fsalfirstlast(cache::PDIRK44Cache,u) = (cache.u, cache.uprev) +get_fsalfirstlast(cache::PDIRK44Cache, u) = (cache.u, cache.uprev) struct PDIRK44ConstantCache{N, TabType} <: OrdinaryDiffEqConstantCache nlsolver::N diff --git a/lib/OrdinaryDiffEqPRK/src/prk_caches.jl b/lib/OrdinaryDiffEqPRK/src/prk_caches.jl index 17c5fa8ba5..059fffcbaa 100644 --- a/lib/OrdinaryDiffEqPRK/src/prk_caches.jl +++ b/lib/OrdinaryDiffEqPRK/src/prk_caches.jl @@ -11,7 +11,7 @@ fsalfirst::rateType tab::TabType end -get_fsalfirstlast(cache::KuttaPRK2p5Cache,u) = (cache.fsalfirst, cache.k) +get_fsalfirstlast(cache::KuttaPRK2p5Cache, u) = (cache.fsalfirst, cache.k) struct KuttaPRK2p5ConstantCache{T, T2} <: OrdinaryDiffEqConstantCache α21::T diff --git a/lib/OrdinaryDiffEqQPRK/src/qprk_caches.jl b/lib/OrdinaryDiffEqQPRK/src/qprk_caches.jl index 667d45823b..8ceb2bc451 100644 --- a/lib/OrdinaryDiffEqQPRK/src/qprk_caches.jl +++ b/lib/OrdinaryDiffEqQPRK/src/qprk_caches.jl @@ -30,7 +30,7 @@ struct QPRK98ConstantCache <: OrdinaryDiffEqConstantCache end thread::Thread end -get_fsalfirstlast(cache::QPRK98Cache,u) = (cache.fsalfirst, cache.k) +get_fsalfirstlast(cache::QPRK98Cache, u) = (cache.fsalfirst, cache.k) function alg_cache(alg::QPRK98, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, diff --git a/lib/OrdinaryDiffEqRKN/src/rkn_caches.jl b/lib/OrdinaryDiffEqRKN/src/rkn_caches.jl index 649052ce8c..e259a1e6ed 100644 --- a/lib/OrdinaryDiffEqRKN/src/rkn_caches.jl +++ b/lib/OrdinaryDiffEqRKN/src/rkn_caches.jl @@ -1,5 +1,5 @@ abstract type NystromMutableCache <: OrdinaryDiffEqMutableCache end -get_fsalfirstlast(cache::NystromMutableCache,u) = (cache.fsalfirst, cache.k) +get_fsalfirstlast(cache::NystromMutableCache, u) = (cache.fsalfirst, cache.k) @cache struct Nystrom4Cache{uType, rateType, reducedRateType} <: NystromMutableCache u::uType diff --git a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl index 4a266e77f3..0a70a7fdde 100644 --- a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl +++ b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl @@ -2,7 +2,7 @@ abstract type RosenbrockMutableCache <: OrdinaryDiffEqMutableCache end abstract type RosenbrockConstantCache <: OrdinaryDiffEqConstantCache end # Fake values since non-FSAL -get_fsalfirstlast(cache::RosenbrockMutableCache,u) = (zero(u), zero(u)) +get_fsalfirstlast(cache::RosenbrockMutableCache, u) = (zero(u), zero(u)) ################################################################################ @@ -1132,10 +1132,13 @@ function alg_cache( constvalue(tTypeNoUnits)), J, W, linsolve) end - -get_fsalfirstlast(cache::Union{Rosenbrock23Cache,Rosenbrock32Cache, Rosenbrock33Cache, -Rosenbrock34Cache, -Rosenbrock4Cache},u) = (cache.fsalfirst, cache.fsallast) +function get_fsalfirstlast( + cache::Union{Rosenbrock23Cache, Rosenbrock32Cache, Rosenbrock33Cache, + Rosenbrock34Cache, + Rosenbrock4Cache}, + u) + (cache.fsalfirst, cache.fsallast) +end ################################################################################ diff --git a/lib/OrdinaryDiffEqSDIRK/src/sdirk_caches.jl b/lib/OrdinaryDiffEqSDIRK/src/sdirk_caches.jl index 21c4f878ee..75a6453fcf 100644 --- a/lib/OrdinaryDiffEqSDIRK/src/sdirk_caches.jl +++ b/lib/OrdinaryDiffEqSDIRK/src/sdirk_caches.jl @@ -1,6 +1,8 @@ abstract type SDIRKMutableCache <: OrdinaryDiffEqMutableCache end abstract type SDIRKConstantCache <: OrdinaryDiffEqConstantCache end -get_fsalfirstlast(cache::SDIRKMutableCache,u) = (cache.fsalfirst, du_alias_or_new(cache.nlsolver, cache.fsalfirst)) +function get_fsalfirstlast(cache::SDIRKMutableCache, u) + (cache.fsalfirst, du_alias_or_new(cache.nlsolver, cache.fsalfirst)) +end @cache mutable struct ImplicitEulerCache{ uType, rateType, uNoUnitsType, N, AV, StepLimiter} <: diff --git a/lib/OrdinaryDiffEqSSPRK/Project.toml b/lib/OrdinaryDiffEqSSPRK/Project.toml index 4875988c14..faccc02874 100644 --- a/lib/OrdinaryDiffEqSSPRK/Project.toml +++ b/lib/OrdinaryDiffEqSSPRK/Project.toml @@ -1,7 +1,7 @@ name = "OrdinaryDiffEqSSPRK" uuid = "669c94d9-1f4b-4b64-b377-1aa079aa2388" authors = ["ParamThakkar123 "] -version = "1.1.0" +version = "1.2.0" [deps] DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" diff --git a/lib/OrdinaryDiffEqSSPRK/src/OrdinaryDiffEqSSPRK.jl b/lib/OrdinaryDiffEqSSPRK/src/OrdinaryDiffEqSSPRK.jl index fde64f8ef7..527bae3c4f 100644 --- a/lib/OrdinaryDiffEqSSPRK/src/OrdinaryDiffEqSSPRK.jl +++ b/lib/OrdinaryDiffEqSSPRK/src/OrdinaryDiffEqSSPRK.jl @@ -46,7 +46,6 @@ PrecompileTools.@compile_workload begin ] low_storage_nonadaptive = [ - ] if Preferences.@load_preference("PrecompileLowStorage", false) @@ -99,6 +98,6 @@ end export SSPRK53_2N2, SSPRK22, SSPRK53, SSPRK63, SSPRK83, SSPRK43, SSPRK432, SSPRKMSVS32, SSPRK54, SSPRK53_2N1, SSPRK104, SSPRK932, SSPRKMSVS43, SSPRK73, SSPRK53_H, - SSPRK33, SHLDDRK_2N, KYKSSPRK42, SHLDDRK52 + SSPRK33, KYKSSPRK42, KYK2014DGSSPRK_3S2 end diff --git a/lib/OrdinaryDiffEqSSPRK/src/alg_utils.jl b/lib/OrdinaryDiffEqSSPRK/src/alg_utils.jl index f68b6e1b18..81102dd71e 100644 --- a/lib/OrdinaryDiffEqSSPRK/src/alg_utils.jl +++ b/lib/OrdinaryDiffEqSSPRK/src/alg_utils.jl @@ -13,6 +13,7 @@ isfsal(alg::SSPRK932) = false isfsal(alg::SSPRK54) = false isfsal(alg::SSPRK104) = false +alg_order(alg::KYK2014DGSSPRK_3S2) = 2 alg_order(alg::KYKSSPRK42) = 2 alg_order(alg::SSPRKMSVS32) = 2 alg_order(alg::SSPRK33) = 3 @@ -30,8 +31,6 @@ 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 @@ -50,3 +49,4 @@ ssp_coefficient(alg::SSPRK22) = 1 ssp_coefficient(alg::SSPRKMSVS32) = 0.5 ssp_coefficient(alg::SSPRKMSVS43) = 0.33 ssp_coefficient(alg::KYKSSPRK42) = 2.459 +ssp_coefficient(alg::KYK2014DGSSPRK_3S2) = 0.8417 diff --git a/lib/OrdinaryDiffEqSSPRK/src/algorithms.jl b/lib/OrdinaryDiffEqSSPRK/src/algorithms.jl index 68c97a9b10..efe6a1ef75 100644 --- a/lib/OrdinaryDiffEqSSPRK/src/algorithms.jl +++ b/lib/OrdinaryDiffEqSSPRK/src/algorithms.jl @@ -320,19 +320,6 @@ function SSPRK33(stage_limiter!, step_limiter! = trivial_limiter!) step_limiter!, False()) end -@doc explicit_rk_docstring("TBD", "SHLDDRK_2N") -Base.@kwdef struct SHLDDRK_2N{StageLimiter, StepLimiter, Thread} <: OrdinaryDiffEqAlgorithm - stage_limiter!::StageLimiter = trivial_limiter! - step_limiter!::StepLimiter = trivial_limiter! - thread::Thread = False() -end -# for backwards compatibility -function SHLDDRK_2N(stage_limiter!, step_limiter! = trivial_limiter!) - SHLDDRK_2N(stage_limiter!, - step_limiter!, - False()) -end - @doc explicit_rk_docstring("TBD", "KYKSSPRK42") Base.@kwdef struct KYKSSPRK42{StageLimiter, StepLimiter, Thread} <: OrdinaryDiffEqAlgorithm @@ -347,15 +334,17 @@ function KYKSSPRK42(stage_limiter!, step_limiter! = trivial_limiter!) False()) end -@doc explicit_rk_docstring("TBD", "SHLDDRK52") -Base.@kwdef struct SHLDDRK52{StageLimiter, StepLimiter, Thread} <: OrdinaryDiffEqAlgorithm +@doc explicit_rk_docstring("TBD", + "KYK2014DGSSPRK_3S2") +Base.@kwdef struct KYK2014DGSSPRK_3S2{StageLimiter, StepLimiter, Thread} <: + OrdinaryDiffEqAlgorithm stage_limiter!::StageLimiter = trivial_limiter! step_limiter!::StepLimiter = trivial_limiter! thread::Thread = False() end # for backwards compatibility -function SHLDDRK52(stage_limiter!, step_limiter! = trivial_limiter!) - SHLDDRK52(stage_limiter!, +function KYK2014DGSSPRK_3S2(stage_limiter!, step_limiter! = trivial_limiter!) + KYK2014DGSSPRK_3S2(stage_limiter!, step_limiter!, False()) -end +end \ No newline at end of file diff --git a/lib/OrdinaryDiffEqSSPRK/src/ssprk_caches.jl b/lib/OrdinaryDiffEqSSPRK/src/ssprk_caches.jl index cdd41f568b..443166b9fd 100644 --- a/lib/OrdinaryDiffEqSSPRK/src/ssprk_caches.jl +++ b/lib/OrdinaryDiffEqSSPRK/src/ssprk_caches.jl @@ -217,176 +217,6 @@ function alg_cache(alg::SSPRK53, u, rate_prototype, ::Type{uEltypeNoUnits}, SSPRK53ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits)) end -@cache struct SHLDDRK52Cache{uType, rateType, TabType, StageLimiter, StepLimiter, Thread} <: SSPRKMutableCache - u::uType - uprev::uType - k::rateType - tmp::uType - fsalfirst::rateType - tab::TabType - stage_limiter!::StageLimiter - step_limiter!::StepLimiter - thread::Thread -end - -struct SHLDDRK52ConstantCache{T1, T2} <: SSPRKConstantCache - α2::T1 - α3::T1 - α4::T1 - α5::T1 - β1::T1 - β2::T1 - β3::T1 - β4::T1 - β5::T1 - c2::T2 - c3::T2 - c4::T2 - c5::T2 -end - -function SHLDDRK52ConstantCache(T1, T2) - α2 = T1(-0.6913065) - α3 = T1(-2.655155) - α4 = T1(-0.8147688) - α5 = T1(-0.6686587) - β1 = T1(0.1) - β2 = T1(0.75) - β3 = T1(0.7) - β4 = T1(0.479313) - β5 = T1(0.310392) - c2 = T2(0.1) - c3 = T2(0.3315201) - c4 = T2(0.4577796) - c5 = T2(0.8666528) - SHLDDRK52ConstantCache(α2, α3, α4, α5, β1, β2, β3, β4, β5, c2, c3, c4, c5) -end - -function alg_cache(alg::SHLDDRK52, u, rate_prototype, ::Type{uEltypeNoUnits}, - ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, - dt, reltol, p, calck, - ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} - SHLDDRK52ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits)) -end - -function alg_cache(alg::SHLDDRK52, u, rate_prototype, ::Type{uEltypeNoUnits}, - ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, - dt, reltol, p, calck, - ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} - tmp = zero(u) - k = zero(rate_prototype) - fsalfirst = zero(rate_prototype) - tab = SHLDDRK52ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits)) - SHLDDRK52Cache(u, uprev, k, tmp, fsalfirst, tab, alg.stage_limiter!, alg.step_limiter!, - alg.thread) -end - -@cache mutable struct SHLDDRK_2NCache{uType, rateType, TabType, StageLimiter, StepLimiter, - Thread} <: SSPRKMutableCache - u::uType - uprev::uType - k::rateType - tmp::uType - fsalfirst::rateType - tab::TabType - step::Int - stage_limiter!::StageLimiter - step_limiter!::StepLimiter - thread::Thread -end - -mutable struct SHLDDRK_2NConstantCache{T1, T2} <: SSPRKConstantCache - α21::T1 - α31::T1 - α41::T1 - α51::T1 - β11::T1 - β21::T1 - β31::T1 - β41::T1 - β51::T1 - c21::T2 - c31::T2 - c41::T2 - c51::T2 - - α22::T1 - α32::T1 - α42::T1 - α52::T1 - α62::T1 - β12::T1 - β22::T1 - β32::T1 - β42::T1 - β52::T1 - β62::T1 - c22::T2 - c32::T2 - c42::T2 - c52::T2 - c62::T2 - - step::Int -end - -function SHLDDRK_2NConstantCache(T1, T2) - α21 = T1(-0.6051226) - α31 = T1(-2.0437564) - α41 = T1(-0.7406999) - α51 = T1(-4.4231765) - β11 = T1(0.2687454) - β21 = T1(0.8014706) - β31 = T1(0.5051570) - β41 = T1(0.5623568) - β51 = T1(0.0590065) - c21 = T2(0.2687454) - c31 = T2(0.5852280) - c41 = T2(0.6827066) - c51 = T2(1.1646854) - - α22 = T1(-0.4412737) - α32 = T1(-1.0739820) - α42 = T1(-1.7063570) - α52 = T1(-2.7979293) - α62 = T1(-4.0913537) - β12 = T1(0.1158488) - β22 = T1(0.3728769) - β32 = T1(0.7379536) - β42 = T1(0.5798110) - β52 = T1(1.0312849) - β62 = T1(0.15) - c22 = T2(0.1158485) - c32 = T2(0.3241850) - c42 = T2(0.6193208) - c52 = T2(0.8034472) - c62 = T2(0.9184166) - SHLDDRK_2NConstantCache( - α21, α31, α41, α51, β11, β21, β31, β41, β51, c21, c31, c41, c51, - α22, α32, α42, α52, α62, β12, β22, β32, β42, β52, β62, c22, c32, - c42, c52, c62, 1) -end - -function alg_cache(alg::SHLDDRK_2N, u, rate_prototype, ::Type{uEltypeNoUnits}, - ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, - dt, reltol, p, calck, - ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} - SHLDDRK_2NConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits)) -end - -function alg_cache(alg::SHLDDRK_2N, u, rate_prototype, ::Type{uEltypeNoUnits}, - ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, - dt, reltol, p, calck, - ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} - tmp = zero(u) - k = zero(rate_prototype) - fsalfirst = zero(rate_prototype) - tab = SHLDDRK_2NConstantCache(constvalue(uBottomEltypeNoUnits), - constvalue(tTypeNoUnits)) - SHLDDRK_2NCache(u, uprev, k, tmp, fsalfirst, tab, 1, alg.stage_limiter!, - alg.step_limiter!, alg.thread) -end - @cache struct SSPRK53_2N1Cache{ uType, rateType, @@ -1255,3 +1085,80 @@ function alg_cache(alg::SSPRK104, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} SSPRK104ConstantCache() end + +@cache struct KYK2014DGSSPRK_3S2_Cache{uType, rateType, TabType, StageLimiter, StepLimiter, + Thread} <: + SSPRKMutableCache + u::uType + uprev::uType + k::rateType + fsalfirst::rateType + tab::TabType + #temporary values for Shu-Osher + u_1::uType + u_2::uType + kk_1::rateType + kk_2::rateType + stage_limiter!::StageLimiter + step_limiter!::StepLimiter + thread::Thread +end + +struct KYK2014DGSSPRK_3S2_ConstantCache{T, T2} <: OrdinaryDiffEqConstantCache + #These are not α and β for RK but for Shu-Osher + #see top of page 317 in + #Optimal Strong-Stability-Preserving Runge–Kutta Time Discretizations for + #Discontinuous Garlekin Methods, Kubatko, Yaeger, Ketcheson 2014 + α_10::T + α_20::T + α_21::T + α_30::T + α_32::T + β_10::T + β_21::T + β_30::T + β_32::T + #Shu-Osher is normally stated for autonomous systems, the times + #are calculated by hand for this scheme + c_1::T + c_2::T + + function KYK2014DGSSPRK_3S2_ConstantCache(T, T2) + α_10 = T(1.0) + α_20 = T(0.087353119859156) + α_21 = T(0.912646880140844) + α_30 = T(0.344956917166841) + α_32 = T(0.655043082833159) + β_10 = T(0.528005024856522) + β_21 = T(0.481882138633993) + β_30 = T(0.022826837460491) + β_32 = T(0.345866039233415) + c_1 = β_10 + c_2 = α_21 * β_10 + β_21 # ==0.96376427726 + new{T, T2}(α_10, α_20, α_21, α_30, α_32, β_10, β_21, β_30, β_32, c_1, c_2) + end +end + +function alg_cache(alg::KYK2014DGSSPRK_3S2, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, + dt, reltol, p, calck, + ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + u_1 = zero(u) + u_2 = zero(u) + kk_1 = zero(rate_prototype) + kk_2 = zero(rate_prototype) + k = zero(rate_prototype) + fsalfirst = zero(rate_prototype) + tab = KYK2014DGSSPRK_3S2_ConstantCache(constvalue(uBottomEltypeNoUnits), + constvalue(tTypeNoUnits)) + KYK2014DGSSPRK_3S2_Cache(u, uprev, k, fsalfirst, tab, u_1, u_2, kk_1, kk_2, + alg.stage_limiter!, alg.step_limiter!, alg.thread) +end + +function alg_cache(alg::KYK2014DGSSPRK_3S2, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, + dt, reltol, p, calck, + ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + KYK2014DGSSPRK_3S2_ConstantCache(constvalue(uBottomEltypeNoUnits), + constvalue(tTypeNoUnits)) +end \ No newline at end of file diff --git a/lib/OrdinaryDiffEqSSPRK/src/ssprk_perform_step.jl b/lib/OrdinaryDiffEqSSPRK/src/ssprk_perform_step.jl index cbdf80f982..2d5ed7ca7d 100644 --- a/lib/OrdinaryDiffEqSSPRK/src/ssprk_perform_step.jl +++ b/lib/OrdinaryDiffEqSSPRK/src/ssprk_perform_step.jl @@ -116,239 +116,6 @@ end integrator.u = u end -function initialize!(integrator, cache::SHLDDRK52ConstantCache) - 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 -end - -@muladd function perform_step!(integrator, cache::SHLDDRK52ConstantCache, - repeat_step = false) - @unpack t, dt, uprev, u, f, p = integrator - @unpack α2, α3, α4, α5, β1, β2, β3, β4, β5, c2, c3, c4, c5 = cache - - # u1 - tmp = dt * integrator.fsalfirst - u = uprev + β1 * tmp - # u2 - tmp = α2 * tmp + dt * f(u, p, t + c2 * dt) - u = u + β2 * tmp - # u3 - tmp = α3 * tmp + dt * f(u, p, t + c3 * dt) - u = u + β3 * tmp - # u4 - tmp = α4 * tmp + dt * f(u, p, t + c4 * dt) - u = u + β4 * tmp - # u5 = u - tmp = α5 * tmp + dt * f(u, p, t + c5 * dt) - u = u + β5 * tmp - - integrator.fsallast = f(u, p, t + dt) # For interpolation, then FSAL'd - integrator.k[1] = integrator.fsalfirst - integrator.k[2] = integrator.fsallast - integrator.u = u -end - -function initialize!(integrator, cache::SHLDDRK52Cache) - @unpack k, fsalfirst = cache - - 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) # FSAL for interpolation -end - -@muladd function perform_step!(integrator, cache::SHLDDRK52Cache, repeat_step = false) - @unpack t, dt, uprev, u, f, p = integrator - @unpack k, fsalfirst, tmp, stage_limiter!, step_limiter!, thread = cache - @unpack α2, α3, α4, α5, β1, β2, β3, β4, β5, c2, c3, c4, c5 = cache.tab - - # u1 - @.. thread=thread tmp=dt * fsalfirst - @.. thread=thread u=uprev + β1 * tmp - stage_limiter!(u, integrator, p, t + c2 * dt) - # u2 - f(k, u, p, t + c2 * dt) - @.. thread=thread tmp=α2 * tmp + dt * k - @.. thread=thread u=u + β2 * tmp - stage_limiter!(u, integrator, p, t + c3 * dt) - # u3 - f(k, u, p, t + c3 * dt) - @.. thread=thread tmp=α3 * tmp + dt * k - @.. thread=thread u=u + β3 * tmp - stage_limiter!(u, integrator, p, t + c4 * dt) - # u4 - f(k, u, p, t + c4 * dt) - @.. thread=thread tmp=α4 * tmp + dt * k - @.. thread=thread u=u + β4 * tmp - stage_limiter!(u, integrator, p, t + c5 * dt) - # u5 = u - f(k, u, p, t + c5 * dt) - @.. thread=thread tmp=α5 * tmp + dt * k - @.. thread=thread u=u + β5 * tmp - stage_limiter!(u, integrator, p, t + dt) - step_limiter!(u, integrator, p, t + dt) - - f(k, u, p, t + dt) -end - -function initialize!(integrator, cache::SHLDDRK_2NConstantCache) - 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 -end - -@muladd function perform_step!(integrator, cache::SHLDDRK_2NConstantCache, - repeat_step = false) - @unpack t, dt, uprev, u, f, p = integrator - @unpack α21, α31, α41, α51, β11, β21, β31, β41, β51, c21, c31, c41, c51, α22, α32, α42, α52, α62, β12, β22, β32, β42, β52, β62, c22, c32, c42, c52, c62 = cache - - if integrator.u_modified - cache.step = 1 - end - # cnt = cache.step - - if cache.step % 2 == 1 - cache.step += 1 - # u1 - tmp = dt * integrator.fsalfirst - u = uprev + β11 * tmp - # u2 - tmp = α21 * tmp + dt * f(u, p, t + c21 * dt) - u = u + β21 * tmp - # u3 - tmp = α31 * tmp + dt * f(u, p, t + c31 * dt) - u = u + β31 * tmp - # u4 - tmp = α41 * tmp + dt * f(u, p, t + c41 * dt) - u = u + β41 * tmp - # u5 = u - tmp = α51 * tmp + dt * f(u, p, t + c51 * dt) - u = u + β51 * tmp - - else - cache.step += 1 - # u1 - tmp = dt * integrator.fsalfirst - u = uprev + β12 * tmp - # u2 - tmp = α22 * tmp + dt * f(u, p, t + c22 * dt) - u = u + β22 * tmp - # u3 - tmp = α32 * tmp + dt * f(u, p, t + c32 * dt) - u = u + β32 * tmp - # u4 - tmp = α42 * tmp + dt * f(u, p, t + c42 * dt) - u = u + β42 * tmp - # u5 = u - tmp = α52 * tmp + dt * f(u, p, t + c52 * dt) - u = u + β52 * tmp - tmp = α62 * tmp + dt * f(u, p, t + c62 * dt) - u = u + β62 * tmp - end - - integrator.fsallast = f(u, p, t + dt) # For interpolation, then FSAL'd - integrator.k[1] = integrator.fsalfirst - integrator.k[2] = integrator.fsallast - integrator.u = u -end - -function initialize!(integrator, cache::SHLDDRK_2NCache) - @unpack k, fsalfirst = cache - - 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) # FSAL for interpolation -end - -@muladd function perform_step!(integrator, cache::SHLDDRK_2NCache, repeat_step = false) - @unpack t, dt, uprev, u, f, p = integrator - @unpack k, fsalfirst, tmp, stage_limiter!, step_limiter!, thread = cache - @unpack α21, α31, α41, α51, β11, β21, β31, β41, β51, c21, c31, c41, c51, α22, α32, α42, α52, α62, β12, β22, β32, β42, β52, β62, c22, c32, c42, c52, c62 = cache.tab - - if integrator.u_modified - cache.step = 1 - end - - if cache.step % 2 == 1 - # u1 - @.. thread=thread tmp=dt * fsalfirst - @.. thread=thread u=uprev + β11 * tmp - stage_limiter!(u, integrator, p, t + c21 * dt) - # u2 - f(k, u, p, t + c21 * dt) - @.. thread=thread tmp=α21 * tmp + dt * k - @.. thread=thread u=u + β21 * tmp - stage_limiter!(u, integrator, p, t + c31 * dt) - # u3 - f(k, u, p, t + c31 * dt) - @.. thread=thread tmp=α31 * tmp + dt * k - @.. thread=thread u=u + β31 * tmp - stage_limiter!(u, integrator, p, t + c41 * dt) - # u4 - f(k, u, p, t + c41 * dt) - @.. thread=thread tmp=α41 * tmp + dt * k - @.. thread=thread u=u + β41 * tmp - stage_limiter!(u, integrator, p, t + c51 * dt) - # u5 = u - f(k, u, p, t + c51 * dt) - @.. thread=thread tmp=α51 * tmp + dt * k - @.. thread=thread u=u + β51 * tmp - stage_limiter!(u, integrator, p, t + dt) - step_limiter!(u, integrator, p, t + dt) - - f(k, u, p, t + dt) - else - # u1 - @.. thread=thread tmp=dt * fsalfirst - @.. thread=thread u=uprev + β12 * tmp - stage_limiter!(u, integrator, p, t + c22 * dt) - # u2 - f(k, u, p, t + c22 * dt) - @.. thread=thread tmp=α22 * tmp + dt * k - @.. thread=thread u=u + β22 * tmp - stage_limiter!(u, integrator, p, t + c32 * dt) - # u3 - f(k, u, p, t + c32 * dt) - @.. thread=thread tmp=α32 * tmp + dt * k - @.. thread=thread u=u + β32 * tmp - stage_limiter!(u, integrator, p, t + c42 * dt) - # u4 - f(k, u, p, t + c42 * dt) - @.. thread=thread tmp=α42 * tmp + dt * k - @.. thread=thread u=u + β42 * tmp - stage_limiter!(u, integrator, p, t + c52 * dt) - # u5 = u - f(k, u, p, t + c52 * dt) - @.. thread=thread tmp=α52 * tmp + dt * k - @.. thread=thread u=u + β52 * tmp - stage_limiter!(u, integrator, p, t + c62 * dt) - # u6 = u - f(k, u, p, t + c62 * dt) - @.. thread=thread tmp=α62 * tmp + dt * k - @.. thread=thread u=u + β62 * tmp - stage_limiter!(u, integrator, p, t + dt) - step_limiter!(u, integrator, p, t + dt) - - f(k, u, p, t + dt) - end -end - function initialize!(integrator, cache::SSPRK33ConstantCache) integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) @@ -1688,3 +1455,65 @@ end step_limiter!(u, integrator, p, t + dt) OrdinaryDiffEqCore.increment_nf!(integrator.stats, 10) end + +function initialize!(integrator, cache::KYK2014DGSSPRK_3S2_ConstantCache) + integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal + OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) + integrator.kshortsize = 2 + integrator.k = typeof(integrator.k)(undef, integrator.kshortsize) + + # Avoid undefined entries if k is an array of arrays + integrator.fsallast = zero(integrator.fsalfirst) + return nothing +end + +@muladd function perform_step!(integrator, cache::KYK2014DGSSPRK_3S2_ConstantCache, + repeat_step = false) + @unpack t, dt, uprev, u, f, p = integrator + @unpack α_10, α_20, α_21, α_30, α_32, β_10, β_21, β_30, β_32, c_1, c_2 = cache + u_1 = α_10 * uprev + dt * β_10 * integrator.fsalfirst + u_2 = (α_20 * uprev + + α_21 * u_1 + dt * β_21 * f(u_1, p, t + c_1 * dt)) + integrator.u = (α_30 * uprev + dt * β_30 * integrator.fsalfirst + + α_32 * u_2 + dt * β_32 * f(u_2, p, t + c_2 * dt)) + integrator.k[1] = integrator.fsalfirst + integrator.k[2] = f(integrator.u, p, t + dt) # For interpolation, then FSAL'd + OrdinaryDiffEqCore.increment_nf!(integrator.stats, 3) + integrator.fsallast = integrator.k[2] + return nothing +end + +function initialize!(integrator, cache::KYK2014DGSSPRK_3S2_Cache) + @unpack k, fsalfirst = cache + + 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) # FSAL for interpolation + OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) + return nothing +end + +@muladd function perform_step!(integrator, cache::KYK2014DGSSPRK_3S2_Cache, + repeat_step = false) + @unpack t, dt, uprev, u, f, p = integrator + @unpack k, fsalfirst, u_1, u_2, kk_1, kk_2, stage_limiter!, step_limiter!, thread = cache + @unpack α_10, α_20, α_21, α_30, α_32, β_10, β_21, β_30, β_32, c_1, c_2 = cache.tab + + @.. broadcast=false thread=thread u_1=α_10 * uprev + dt * β_10 * integrator.fsalfirst + stage_limiter!(u_1, integrator, p, t + c_1 * dt) + f(kk_1, u_1, p, t + c_1 * dt) + @.. broadcast=false thread=thread u_2=(α_20 * uprev + + α_21 * u_1 + dt * β_21 * kk_1) + stage_limiter!(u_2, integrator, p, t + c_2 * dt) + f(kk_2, u_2, p, t + c_2 * dt) + @.. broadcast=false thread=thread u=(α_30 * uprev + + dt * β_30 * integrator.fsalfirst + + α_32 * u_2 + dt * β_32 * kk_2) + stage_limiter!(u, integrator, p, t + dt) + step_limiter!(u, integrator, p, t + dt) + f(integrator.k[2], u, p, t + dt) # For interpolation, then FSAL'd + OrdinaryDiffEqCore.increment_nf!(integrator.stats, 3) + return nothing +end \ No newline at end of file diff --git a/lib/OrdinaryDiffEqSSPRK/test/ode_ssprk_tests.jl b/lib/OrdinaryDiffEqSSPRK/test/ode_ssprk_tests.jl index 31facc3d6a..ec0c73f11b 100644 --- a/lib/OrdinaryDiffEqSSPRK/test/ode_ssprk_tests.jl +++ b/lib/OrdinaryDiffEqSSPRK/test/ode_ssprk_tests.jl @@ -102,39 +102,6 @@ sol = solve(test_problem_ssp_long, alg, dt = OrdinaryDiffEqSSPRK.ssp_coefficient dense = false) @test all(sol.u .>= 0) -println("SHLDDRK52") -alg = SHLDDRK52() -for prob in test_problems_only_time - sim = test_convergence(dts, prob, alg) - @test_broken sim.𝒪est[:final]≈OrdinaryDiffEqSSPRK.alg_order(alg) atol=testTol -end -for prob in test_problems_linear - sim = test_convergence(dts, prob, alg) - @test_broken sim.𝒪est[:final]≈OrdinaryDiffEqSSPRK.alg_order(alg) atol=testTol -end -for prob in test_problems_nonlinear - sim = test_convergence(dts, prob, alg) - @test_broken sim.𝒪est[:final]≈OrdinaryDiffEqSSPRK.alg_order(alg) atol=testTol -end - -println("SHLDDRK_2N") -dts_SHLDDRK_2N = (1 / 2) .^ (0:3) -alg = SHLDDRK_2N() -for prob in test_problems_only_time - sim = test_convergence(dts_SHLDDRK_2N, prob, alg) - @test sim.𝒪est[:final]≈4 atol=0.46 -end -for prob in test_problems_linear - sim = test_convergence(dts_SHLDDRK_2N, prob, alg) - @test sim.𝒪est[:final]≈4 atol=0.46 -end -for prob in test_problems_nonlinear - sim = test_convergence(dts_SHLDDRK_2N, prob, alg) - @test sim.𝒪est[:final]≈4 atol=1 - # due to unusual saturation towards high dts(0.5 and onwards) and - # saturation towards low dts due to less precision in the provided values of weights , tolerance is kept so high -end - println("SSPRK33") alg = SSPRK33() for prob in test_problems_only_time @@ -507,7 +474,7 @@ integ = init(prob_ode_large, alg, dt = 1.e-2, save_start = false, save_end = fal @test Base.summarysize(integ) ÷ Base.summarysize(u0_large) <= 5 println("KYK2014DGSSPRK_3S2") -alg = OrdinaryDiffEqLowStorageRK.KYK2014DGSSPRK_3S2() +alg = KYK2014DGSSPRK_3S2() for prob in test_problems_only_time sim = test_convergence(dts, prob, alg) @test abs(sim.𝒪est[:final] - OrdinaryDiffEqSSPRK.alg_order(alg)) < testTol diff --git a/lib/OrdinaryDiffEqStabilizedIRK/src/irkc_caches.jl b/lib/OrdinaryDiffEqStabilizedIRK/src/irkc_caches.jl index 80f92326f9..60239258c9 100644 --- a/lib/OrdinaryDiffEqStabilizedIRK/src/irkc_caches.jl +++ b/lib/OrdinaryDiffEqStabilizedIRK/src/irkc_caches.jl @@ -23,7 +23,9 @@ end constantcache::C end -get_fsalfirstlast(cache::IRKCCache,u) = (cache.fsalfirst, du_alias_or_new(cache.nlsolver, cache.fsalfirst)) +function get_fsalfirstlast(cache::IRKCCache, u) + (cache.fsalfirst, du_alias_or_new(cache.nlsolver, cache.fsalfirst)) +end function alg_cache(alg::IRKC, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, diff --git a/lib/OrdinaryDiffEqStabilizedRK/src/rkc_caches.jl b/lib/OrdinaryDiffEqStabilizedRK/src/rkc_caches.jl index cbdaa5cc2e..e598d1243c 100644 --- a/lib/OrdinaryDiffEqStabilizedRK/src/rkc_caches.jl +++ b/lib/OrdinaryDiffEqStabilizedRK/src/rkc_caches.jl @@ -1,5 +1,5 @@ abstract type StabilizedRKMutableCache <: OrdinaryDiffEqMutableCache end -get_fsalfirstlast(cache::StabilizedRKMutableCache,u) = (cache.fsalfirst, cache.k) +get_fsalfirstlast(cache::StabilizedRKMutableCache, u) = (cache.fsalfirst, cache.k) mutable struct ROCK2ConstantCache{T, T2, zType} <: OrdinaryDiffEqConstantCache ms::SVector{46, Int} diff --git a/lib/OrdinaryDiffEqSymplecticRK/src/symplectic_caches.jl b/lib/OrdinaryDiffEqSymplecticRK/src/symplectic_caches.jl index 6112a89792..dca14b6d28 100644 --- a/lib/OrdinaryDiffEqSymplecticRK/src/symplectic_caches.jl +++ b/lib/OrdinaryDiffEqSymplecticRK/src/symplectic_caches.jl @@ -421,4 +421,7 @@ function alg_cache(alg::SofSpa10, u, rate_prototype, ::Type{uEltypeNoUnits}, SofSpa10ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits)) end -get_fsalfirstlast(cache::Union{HamiltonMutableCache, VelocityVerletCache, SymplecticEulerCache},u) = (cache.fsalfirst, cache.k) +function get_fsalfirstlast( + cache::Union{HamiltonMutableCache, VelocityVerletCache, SymplecticEulerCache}, u) + (cache.fsalfirst, cache.k) +end diff --git a/lib/OrdinaryDiffEqSymplecticRK/src/symplectic_perform_step.jl b/lib/OrdinaryDiffEqSymplecticRK/src/symplectic_perform_step.jl index 06173ff5f2..2e581be283 100644 --- a/lib/OrdinaryDiffEqSymplecticRK/src/symplectic_perform_step.jl +++ b/lib/OrdinaryDiffEqSymplecticRK/src/symplectic_perform_step.jl @@ -126,7 +126,6 @@ end function initialize!(integrator, cache::C) where {C <: Union{HamiltonMutableCache, VelocityVerletCache}} - integrator.kshortsize = 2 resize!(integrator.k, integrator.kshortsize) integrator.k[1] = integrator.fsalfirst diff --git a/lib/OrdinaryDiffEqTsit5/src/tsit_caches.jl b/lib/OrdinaryDiffEqTsit5/src/tsit_caches.jl index b8bc6b611b..040f9a9869 100644 --- a/lib/OrdinaryDiffEqTsit5/src/tsit_caches.jl +++ b/lib/OrdinaryDiffEqTsit5/src/tsit_caches.jl @@ -36,7 +36,7 @@ function alg_cache(alg::Tsit5, u, rate_prototype, ::Type{uEltypeNoUnits}, alg.stage_limiter!, alg.step_limiter!, alg.thread) end -get_fsalfirstlast(cache::Tsit5Cache,u) = (cache.k1, cache.k7) +get_fsalfirstlast(cache::Tsit5Cache, u) = (cache.k1, cache.k7) function alg_cache(alg::Tsit5, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, diff --git a/lib/OrdinaryDiffEqVerner/src/verner_caches.jl b/lib/OrdinaryDiffEqVerner/src/verner_caches.jl index a118fc9f24..19289eb2c2 100644 --- a/lib/OrdinaryDiffEqVerner/src/verner_caches.jl +++ b/lib/OrdinaryDiffEqVerner/src/verner_caches.jl @@ -23,7 +23,7 @@ lazy::Bool end -get_fsalfirstlast(cache::Vern6Cache,u) = (cache.k1, cache.k9) +get_fsalfirstlast(cache::Vern6Cache, u) = (cache.k1, cache.k9) function alg_cache(alg::Vern6, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, @@ -87,7 +87,7 @@ end end # fake values since non-FSAL method -get_fsalfirstlast(cache::Vern7Cache,u) = (cache.k1, cache.k2) +get_fsalfirstlast(cache::Vern7Cache, u) = (cache.k1, cache.k2) function alg_cache(alg::Vern7, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, @@ -153,7 +153,7 @@ end end # fake values since non-FSAL method -get_fsalfirstlast(cache::Vern8Cache,u) = (cache.k1, cache.k2) +get_fsalfirstlast(cache::Vern8Cache, u) = (cache.k1, cache.k2) function alg_cache(alg::Vern8, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, @@ -227,7 +227,7 @@ end end # fake values since non-FSAL method -get_fsalfirstlast(cache::Vern9Cache,u) = (cache.k1, cache.k2) +get_fsalfirstlast(cache::Vern9Cache, u) = (cache.k1, cache.k2) function alg_cache(alg::Vern9, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index dfafe9b9d3..503484ac33 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -58,7 +58,7 @@ import OrdinaryDiffEqCore: trivial_limiter!, CompositeAlgorithm, alg_order, _change_t_via_interpolation!, ODEIntegrator, _ode_interpolant!, current_interpolant, resize_nlsolver!, _ode_interpolant, handle_tstop!, _postamble!, update_uprev!, resize_J_W!, - DAEAlgorithm, get_fsalfirstlast, strip_cache, strip_interpolation + DAEAlgorithm, get_fsalfirstlast, strip_cache export CompositeAlgorithm, ShampineCollocationInit, BrownFullBasicInit, NoInit AutoSwitch @@ -101,12 +101,12 @@ export ORK256, CarpenterKennedy2N54, SHLDDRK64, HSLDDRK64, DGLDDRK73_C, DGLDDRK8 ParsaniKetchesonDeconinck3S94, ParsaniKetchesonDeconinck3S184, ParsaniKetchesonDeconinck3S105, ParsaniKetchesonDeconinck3S205, RDPK3Sp35, RDPK3SpFSAL35, RDPK3Sp49, RDPK3SpFSAL49, RDPK3Sp510, RDPK3SpFSAL510, - KYK2014DGSSPRK_3S2, RK46NL + RK46NL, SHLDDRK_2N, SHLDDRK52 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 + SSPRK33, KYKSSPRK42, KYK2014DGSSPRK_3S2 using OrdinaryDiffEqFeagin export Feagin10, Feagin12, Feagin14 diff --git a/test/interface/default_solver_tests.jl b/test/interface/default_solver_tests.jl index 4cdef37e80..3c56f9f7ea 100644 --- a/test/interface/default_solver_tests.jl +++ b/test/interface/default_solver_tests.jl @@ -112,9 +112,9 @@ SA_ode_problem = ODEProblem((u, p, t) -> zero(u), SA[0], 2) @test solve(SA_ode_problem; callback = cb).retcode == ReturnCode.Success # test Complex numbers -H(s) = (1-s) * complex([0 1; 1 0]) + s * complex([1 0; 0 -1]) +H(s) = (1 - s) * complex([0 1; 1 0]) + s * complex([1 0; 0 -1]) schrod_eq(state, time, s) = -im * time * H(s) * state -prob_complex = ODEProblem(schrod_eq, complex([1, -1]/sqrt(2)), (0,1), 100) +prob_complex = ODEProblem(schrod_eq, complex([1, -1] / sqrt(2)), (0, 1), 100) complex_sol = solve(prob_complex) @test complex_sol.retcode == ReturnCode.Success diff --git a/test/interface/linear_solver_test.jl b/test/interface/linear_solver_test.jl index 12a449d603..a615deb515 100644 --- a/test/interface/linear_solver_test.jl +++ b/test/interface/linear_solver_test.jl @@ -161,7 +161,7 @@ end using OrdinaryDiffEq, StaticArrays, LinearSolve, ParameterizedFunctions hires = @ode_def Hires begin - dy1 = -1.71f0 * y1 + 0.43f0 * y2 + 8.32f0 * y3 + 0.0007f0 + 1f-18*t + dy1 = -1.71f0 * y1 + 0.43f0 * y2 + 8.32f0 * y3 + 0.0007f0 + 1.0f-18 * t dy2 = 1.71f0 * y1 - 8.75f0 * y2 dy3 = -10.03f0 * y3 + 0.43f0 * y4 + 0.035f0 * y5 dy4 = 8.32f0 * y2 + 1.71f0 * y3 - 1.12f0 * y4 @@ -178,11 +178,11 @@ u0[8] = 0.0057 probiip = ODEProblem{true}(hires, u0, (0.0, 10.0)) proboop = ODEProblem{false}(hires, u0, (0.0, 10.0)) probstatic = ODEProblem{false}(hires, SVector{8}(u0), (0.0, 10.0)) -probiipf32 = ODEProblem{true}(hires, Float32.(u0), (0f0, 10f0)) -proboopf32 = ODEProblem{false}(hires, Float32.(u0), (0f0, 10f0)) -probstaticf32 = ODEProblem{false}(hires, SVector{8}(Float32.(u0)), (0f0, 10f0)) +probiipf32 = ODEProblem{true}(hires, Float32.(u0), (0.0f0, 10.0f0)) +proboopf32 = ODEProblem{false}(hires, Float32.(u0), (0.0f0, 10.0f0)) +probstaticf32 = ODEProblem{false}(hires, SVector{8}(Float32.(u0)), (0.0f0, 10.0f0)) probs = (; probiip, proboop, probstatic) -probsf32 = (;probiipf32, proboopf32, probstaticf32) +probsf32 = (; probiipf32, proboopf32, probstaticf32) qndf = QNDF() krylov_qndf = QNDF(linsolve = KrylovJL_GMRES()) fbdf = FBDF() diff --git a/test/interface/noindex_tests.jl b/test/interface/noindex_tests.jl index f013c59fd9..3aca1a635d 100644 --- a/test/interface/noindex_tests.jl +++ b/test/interface/noindex_tests.jl @@ -56,17 +56,18 @@ for alg in algs @test_nowarn sol(similar(prob.u0), 0.1) end - struct CustomArray{T, N} x::Array{T, N} end Base.size(x::CustomArray) = size(x.x) Base.axes(x::CustomArray) = axes(x.x) Base.ndims(x::CustomArray) = ndims(x.x) -Base.ndims(::Type{<:CustomArray{T,N}}) where {T,N} = N +Base.ndims(::Type{<:CustomArray{T, N}}) where {T, N} = N Base.zero(x::CustomArray) = CustomArray(zero(x.x)) -Base.zero(::Type{<:CustomArray{T,N}}) where {T,N} = CustomArray(zero(Array{T,N})) -Base.similar(x::CustomArray, dims::Union{Integer, AbstractUnitRange}...) = CustomArray(similar(x.x, dims...)) +Base.zero(::Type{<:CustomArray{T, N}}) where {T, N} = CustomArray(zero(Array{T, N})) +function Base.similar(x::CustomArray, dims::Union{Integer, AbstractUnitRange}...) + CustomArray(similar(x.x, dims...)) +end Base.copyto!(x::CustomArray, y::CustomArray) = CustomArray(copyto!(x.x, y.x)) Base.copy(x::CustomArray) = CustomArray(copy(x.x)) Base.length(x::CustomArray) = length(x.x) @@ -82,21 +83,28 @@ Base.any(f::Function, x::CustomArray; kwargs...) = any(f, x.x; kwargs...) Base.all(f::Function, x::CustomArray; kwargs...) = all(f, x.x; kwargs...) Base.similar(x::CustomArray, t) = CustomArray(similar(x.x, t)) Base.:(==)(x::CustomArray, y::CustomArray) = x.x == y.x -Base.:(*)(x::Number, y::CustomArray) = CustomArray(x*y.x) -Base.:(/)(x::CustomArray, y::Number) = CustomArray(x.x/y) +Base.:(*)(x::Number, y::CustomArray) = CustomArray(x * y.x) +Base.:(/)(x::CustomArray, y::Number) = CustomArray(x.x / y) LinearAlgebra.norm(x::CustomArray) = norm(x.x) struct CustomStyle{N} <: Broadcast.BroadcastStyle where {N} end -CustomStyle(::Val{N}) where N = CustomStyle{N}() -CustomStyle{M}(::Val{N}) where {N,M} = NoIndexStyle{N}() -Base.BroadcastStyle(::Type{<:CustomArray{T,N}}) where {T,N} = CustomStyle{N}() -Broadcast.BroadcastStyle(::CustomStyle{N}, ::Broadcast.DefaultArrayStyle{0}) where {N} = CustomStyle{N}() -Base.similar(bc::Base.Broadcast.Broadcasted{CustomStyle{N}}, ::Type{ElType}) where {N, ElType} = CustomArray(similar(Array{ElType, N}, axes(bc))) +CustomStyle(::Val{N}) where {N} = CustomStyle{N}() +CustomStyle{M}(::Val{N}) where {N, M} = NoIndexStyle{N}() +Base.BroadcastStyle(::Type{<:CustomArray{T, N}}) where {T, N} = CustomStyle{N}() +function Broadcast.BroadcastStyle( + ::CustomStyle{N}, ::Broadcast.DefaultArrayStyle{0}) where {N} + CustomStyle{N}() +end +function Base.similar( + bc::Base.Broadcast.Broadcasted{CustomStyle{N}}, ::Type{ElType}) where {N, ElType} + CustomArray(similar(Array{ElType, N}, axes(bc))) +end Base.Broadcast._broadcast_getindex(x::CustomArray, i) = x.x[i] Base.Broadcast.extrude(x::CustomArray) = x Base.Broadcast.broadcastable(x::CustomArray) = x -@inline function Base.copyto!(dest::CustomArray, bc::Base.Broadcast.Broadcasted{<:CustomStyle}) +@inline function Base.copyto!( + dest::CustomArray, bc::Base.Broadcast.Broadcasted{<:CustomStyle}) axes(dest) == axes(bc) || throwdm(axes(dest), axes(bc)) bc′ = Base.Broadcast.preprocess(dest, bc) dest′ = dest.x @@ -127,7 +135,7 @@ RecursiveArrayTools.recursivefill!(x::CustomArray, a) = fill!(x, a) Base.show_vector(io::IO, x::CustomArray) = Base.show_vector(io, x.x) -Base.show(io::IO, x::CustomArray) = (print(io, "CustomArray");show(io, x.x)) +Base.show(io::IO, x::CustomArray) = (print(io, "CustomArray"); show(io, x.x)) function Base.show(io::IO, ::MIME"text/plain", x::CustomArray) println(io, Base.summary(x), ":") Base.print_array(io, x.x) @@ -139,4 +147,4 @@ for alg in algs sol = @test_nowarn solve(prob, alg) @test_nowarn sol(0.1) @test_nowarn sol(similar(prob.u0), 0.1) -end \ No newline at end of file +end diff --git a/test/interface/ode_strip_test.jl b/test/interface/ode_strip_test.jl index ed5a448051..e476c024ee 100644 --- a/test/interface/ode_strip_test.jl +++ b/test/interface/ode_strip_test.jl @@ -14,4 +14,4 @@ prob = ODEProblem(lorenz!, u0, tspan) sol = solve(prob, Rosenbrock23()) @test isnothing(SciMLBase.strip_interpolation(sol.interp).f) -@test isnothing(SciMLBase.strip_interpolation(sol.interp).cache.jac_config) \ No newline at end of file +@test isnothing(SciMLBase.strip_interpolation(sol.interp).cache.jac_config)