diff --git a/lib/OrdinaryDiffEqBDF/Project.toml b/lib/OrdinaryDiffEqBDF/Project.toml new file mode 100644 index 0000000000..4c0e23c6f0 --- /dev/null +++ b/lib/OrdinaryDiffEqBDF/Project.toml @@ -0,0 +1,27 @@ +name = "OrdinaryDiffEqBDF" +uuid = "6ad6398a-0878-4a85-9266-38940aa047c8" +authors = ["ParamThakkar123 "] +version = "1.0.0" + +[deps] +OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +TruncatedStacktraces = "781d530d-4396-4725-bb49-402e4bee1e77" +MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" +MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" +FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898" +RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" +ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" + +[extras] +DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[compat] +julia = "1.10" + +[targets] +test = ["DiffEqDevTools", "Random", "SafeTestsets", "Test"] \ No newline at end of file diff --git a/lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl b/lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl new file mode 100644 index 0000000000..2080050161 --- /dev/null +++ b/lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl @@ -0,0 +1,41 @@ +module OrdinaryDiffEqBDF + +import OrdinaryDiffEq: alg_order, calculate_residuals!, + initialize!, perform_step!, @unpack, unwrap_alg, + calculate_residuals, alg_extrapolates, + OrdinaryDiffEqAlgorithm, + OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, + OrdinaryDiffEqNewtonAdaptiveAlgorithm, + OrdinaryDiffEqNewtonAlgorithm, + AbstractController, DEFAULT_PRECS, + CompiledFloats, uses_uprev, + NLNewton, alg_cache, _vec, _reshape, @cache, + isfsal, full_cache, build_nlsolver, + nlsolve!, nlsolvefail, isnewton, + constvalue, _unwrap_val, + DIRK, set_new_W!, + du_alias_or_new, trivial_limiter!, + ImplicitEulerConstantCache, + ImplicitEulerCache, COEFFICIENT_MULTISTEP, + markfirststage!, UJacobianWrapper, mul!, + issplit, qsteady_min_default, qsteady_max_default, + get_current_alg_order, get_current_adaptive_order, + default_controller, stepsize_controller!, step_accept_controller!, + step_reject_controller!, post_newton_controller!, + u_modified! +using TruncatedStacktraces, MuladdMacro, MacroTools, FastBroadcast, RecursiveArrayTools +import StaticArrays: SArray, MVector, SVector, @SVector, StaticArray, MMatrix, SA +using LinearAlgebra: I +using ArrayInterface + +include("algorithms.jl") +include("alg_utils.jl") +include("bdf_utils.jl") +include("bdf_caches.jl") +include("controllers.jl") +include("bdf_perform_step.jl") + +export ABDF2, QNDF1, QBDF1, QNDF2, QBDF2, QNDF, QBDF, FBDF, + SBDF2, SBDF3, SBDF4, MEBDF2, IMEXEuler, IMEXEulerARK + +end diff --git a/lib/OrdinaryDiffEqBDF/src/alg_utils.jl b/lib/OrdinaryDiffEqBDF/src/alg_utils.jl new file mode 100644 index 0000000000..aea50d0df0 --- /dev/null +++ b/lib/OrdinaryDiffEqBDF/src/alg_utils.jl @@ -0,0 +1,26 @@ +alg_extrapolates(alg::ABDF2) = true +alg_extrapolates(alg::SBDF) = true +alg_extrapolates(alg::MEBDF2) = true + +alg_order(alg::ABDF2) = 2 +alg_order(alg::SBDF) = alg.order +alg_order(alg::QNDF1) = 1 +alg_order(alg::QNDF2) = 2 +alg_order(alg::QNDF) = 1 #dummy value +alg_order(alg::MEBDF2) = 2 +alg_order(alg::FBDF) = 1 #dummy value + +issplit(alg::SBDF) = true + +qsteady_min_default(alg::FBDF) = 9 // 10 + +qsteady_max_default(alg::QNDF) = 2 // 1 +qsteady_max_default(alg::QNDF2) = 2 // 1 +qsteady_max_default(alg::QNDF1) = 2 // 1 +qsteady_max_default(alg::FBDF) = 2 // 1 + +get_current_alg_order(alg::QNDF, cache) = cache.order +get_current_alg_order(alg::FBDF, cache) = cache.order + +get_current_adaptive_order(alg::QNDF, cache) = cache.order +get_current_adaptive_order(alg::FBDF, cache) = cache.order diff --git a/lib/OrdinaryDiffEqBDF/src/algorithms.jl b/lib/OrdinaryDiffEqBDF/src/algorithms.jl new file mode 100644 index 0000000000..0a91dce0b6 --- /dev/null +++ b/lib/OrdinaryDiffEqBDF/src/algorithms.jl @@ -0,0 +1,390 @@ +""" +E. Alberdi Celayaa, J. J. Anza Aguirrezabalab, P. Chatzipantelidisc. Implementation of +an Adaptive BDF2 Formula and Comparison with The MATLAB Ode15s. Procedia Computer Science, +29, pp 1014-1026, 2014. doi: https://doi.org/10.1016/j.procs.2014.05.091 + +ABDF2: Multistep Method +An adaptive order 2 L-stable fixed leading coefficient multistep BDF method. +""" +struct ABDF2{CS, AD, F, F2, P, FDT, ST, CJ, K, T, StepLimiter} <: + OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + κ::K + tol::T + smooth_est::Bool + extrapolant::Symbol + controller::Symbol + step_limiter!::StepLimiter +end +function ABDF2(; chunk_size = Val{0}(), autodiff = true, standardtag = Val{true}(), + concrete_jac = nothing, diff_type = Val{:forward}, + κ = nothing, tol = nothing, linsolve = nothing, precs = DEFAULT_PRECS, + nlsolve = NLNewton(), + smooth_est = true, extrapolant = :linear, + controller = :Standard, step_limiter! = trivial_limiter!) + ABDF2{ + _unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve), + typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac), + typeof(κ), typeof(tol), typeof(step_limiter!)}(linsolve, nlsolve, precs, κ, tol, + smooth_est, extrapolant, controller, step_limiter!) +end + +""" +Uri M. Ascher, Steven J. Ruuth, Brian T. R. Wetton. Implicit-Explicit Methods for Time- +Dependent Partial Differential Equations. 1995 Society for Industrial and Applied Mathematics +Journal on Numerical Analysis, 32(3), pp 797-823, 1995. doi: https://doi.org/10.1137/0732037 +""" +struct SBDF{CS, AD, F, F2, P, FDT, ST, CJ, K, T} <: + OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + κ::K + tol::T + extrapolant::Symbol + order::Int + ark::Bool +end + +function SBDF(order; chunk_size = Val{0}(), autodiff = Val{true}(), + standardtag = Val{true}(), concrete_jac = nothing, diff_type = Val{:forward}, + linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), κ = nothing, + tol = nothing, + extrapolant = :linear, ark = false) + SBDF{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve), + typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac), + typeof(κ), typeof(tol)}(linsolve, + nlsolve, + precs, + κ, + tol, + extrapolant, + order, + ark) +end + +# All keyword form needed for remake +function SBDF(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(), + concrete_jac = nothing, diff_type = Val{:forward}, + linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), κ = nothing, + tol = nothing, + extrapolant = :linear, + order, ark = false) + SBDF{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve), + typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac), + typeof(κ), typeof(tol)}(linsolve, + nlsolve, + precs, + κ, + tol, + extrapolant, + order, + ark) +end + +""" + SBDF2(;kwargs...) + +The two-step version of the IMEX multistep methods of + + - Uri M. Ascher, Steven J. Ruuth, Brian T. R. Wetton. + Implicit-Explicit Methods for Time-Dependent Partial Differential Equations. + Society for Industrial and Applied Mathematics. + Journal on Numerical Analysis, 32(3), pp 797-823, 1995. + doi: [https://doi.org/10.1137/0732037](https://doi.org/10.1137/0732037) + +See also `SBDF`. +""" +SBDF2(; kwargs...) = SBDF(2; kwargs...) + +""" + SBDF3(;kwargs...) + +The three-step version of the IMEX multistep methods of + + - Uri M. Ascher, Steven J. Ruuth, Brian T. R. Wetton. + Implicit-Explicit Methods for Time-Dependent Partial Differential Equations. + Society for Industrial and Applied Mathematics. + Journal on Numerical Analysis, 32(3), pp 797-823, 1995. + doi: [https://doi.org/10.1137/0732037](https://doi.org/10.1137/0732037) + +See also `SBDF`. +""" +SBDF3(; kwargs...) = SBDF(3; kwargs...) + +""" + SBDF4(;kwargs...) + +The four-step version of the IMEX multistep methods of + + - Uri M. Ascher, Steven J. Ruuth, Brian T. R. Wetton. + Implicit-Explicit Methods for Time-Dependent Partial Differential Equations. + Society for Industrial and Applied Mathematics. + Journal on Numerical Analysis, 32(3), pp 797-823, 1995. + doi: [https://doi.org/10.1137/0732037](https://doi.org/10.1137/0732037) + +See also `SBDF`. +""" +SBDF4(; kwargs...) = SBDF(4; kwargs...) + +""" +QNDF1: Multistep Method +An adaptive order 1 quasi-constant timestep L-stable numerical differentiation function (NDF) method. +Optional parameter kappa defaults to Shampine's accuracy-optimal -0.1850. + +See also `QNDF`. +""" +struct QNDF1{CS, AD, F, F2, P, FDT, ST, CJ, κType, StepLimiter} <: + OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + extrapolant::Symbol + kappa::κType + controller::Symbol + step_limiter!::StepLimiter +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 = -0.1850, + controller = :Standard, step_limiter! = trivial_limiter!) + QNDF1{ + _unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve), + typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac), + typeof(kappa), typeof(step_limiter!)}(linsolve, + nlsolve, + precs, + extrapolant, + kappa, + controller, + step_limiter!) +end + +""" +QNDF2: Multistep Method +An adaptive order 2 quasi-constant timestep L-stable numerical differentiation function (NDF) method. + +See also `QNDF`. +""" +struct QNDF2{CS, AD, F, F2, P, FDT, ST, CJ, κType, StepLimiter} <: + OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + extrapolant::Symbol + kappa::κType + controller::Symbol + step_limiter!::StepLimiter +end + +function QNDF2(; 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 = -1 // 9, + controller = :Standard, step_limiter! = trivial_limiter!) + QNDF2{ + _unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve), + typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac), + typeof(kappa), typeof(step_limiter!)}(linsolve, + nlsolve, + precs, + extrapolant, + kappa, + controller, + step_limiter!) +end + +""" +QNDF: Multistep Method +An adaptive order quasi-constant timestep NDF method. +Utilizes Shampine's accuracy-optimal kappa values as defaults (has a keyword argument for a tuple of kappa coefficients). + +@article{shampine1997matlab, +title={The matlab ode suite}, +author={Shampine, Lawrence F and Reichelt, Mark W}, +journal={SIAM journal on scientific computing}, +volume={18}, +number={1}, +pages={1--22}, +year={1997}, +publisher={SIAM} +} +""" +struct QNDF{MO, CS, AD, F, F2, P, FDT, ST, CJ, K, T, κType, StepLimiter} <: + OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} + max_order::Val{MO} + linsolve::F + nlsolve::F2 + precs::P + κ::K + tol::T + extrapolant::Symbol + kappa::κType + controller::Symbol + step_limiter!::StepLimiter +end + +function QNDF(; max_order::Val{MO} = Val{5}(), chunk_size = Val{0}(), + autodiff = Val{true}(), standardtag = Val{true}(), concrete_jac = nothing, + diff_type = Val{:forward}, + linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), κ = nothing, + tol = nothing, + extrapolant = :linear, kappa = promote(-0.1850, -1 // 9, -0.0823, -0.0415, 0), + 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), + _unwrap_val(concrete_jac), + typeof(κ), typeof(tol), typeof(kappa), typeof(step_limiter!)}( + max_order, linsolve, nlsolve, precs, κ, tol, + extrapolant, kappa, controller, step_limiter!) +end + +TruncatedStacktraces.@truncate_stacktrace QNDF + +""" +MEBDF2: Multistep Method +The second order Modified Extended BDF method, which has improved stability properties over the standard BDF. +Fixed timestep only. +""" +struct MEBDF2{CS, AD, F, F2, P, FDT, ST, CJ} <: + OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + extrapolant::Symbol +end +function MEBDF2(; chunk_size = Val{0}(), autodiff = true, standardtag = Val{true}(), + concrete_jac = nothing, diff_type = Val{:forward}, + linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), + extrapolant = :constant) + MEBDF2{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), + _unwrap_val(concrete_jac)}(linsolve, + nlsolve, + precs, + extrapolant) +end + +""" +FBDF: Fixed leading coefficient BDF + +An adaptive order quasi-constant timestep NDF method. +Utilizes Shampine's accuracy-optimal kappa values as defaults (has a keyword argument for a tuple of kappa coefficients). + +@article{shampine2002solving, +title={Solving 0= F (t, y (t), y′(t)) in Matlab}, +author={Shampine, Lawrence F}, +year={2002}, +publisher={Walter de Gruyter GmbH \\& Co. KG} +} +""" +struct FBDF{MO, CS, AD, F, F2, P, FDT, ST, CJ, K, T, StepLimiter} <: + OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} + max_order::Val{MO} + linsolve::F + nlsolve::F2 + precs::P + κ::K + tol::T + extrapolant::Symbol + controller::Symbol + step_limiter!::StepLimiter +end + +function FBDF(; max_order::Val{MO} = Val{5}(), chunk_size = Val{0}(), + autodiff = Val{true}(), standardtag = Val{true}(), concrete_jac = nothing, + diff_type = Val{:forward}, + linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), κ = nothing, + tol = nothing, + extrapolant = :linear, controller = :Standard, step_limiter! = trivial_limiter!) where {MO} + FBDF{MO, _unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), + _unwrap_val(concrete_jac), + typeof(κ), typeof(tol), typeof(step_limiter!)}( + max_order, linsolve, nlsolve, precs, κ, tol, extrapolant, + controller, step_limiter!) +end + +TruncatedStacktraces.@truncate_stacktrace FBDF + +""" +QBDF1: Multistep Method + +An alias of `QNDF1` with κ=0. +""" +QBDF1(; kwargs...) = QNDF1(; kappa = 0, kwargs...) + +""" +QBDF2: Multistep Method + +An alias of `QNDF2` with κ=0. +""" +QBDF2(; kwargs...) = QNDF2(; kappa = 0, kwargs...) + +""" +QBDF: Multistep Method + +An alias of `QNDF` with κ=0. +""" +QBDF(; kwargs...) = QNDF(; kappa = tuple(0 // 1, 0 // 1, 0 // 1, 0 // 1, 0 // 1), kwargs...) + +""" + IMEXEuler(;kwargs...) + +The one-step version of the IMEX multistep methods of + + - Uri M. Ascher, Steven J. Ruuth, Brian T. R. Wetton. + Implicit-Explicit Methods for Time-Dependent Partial Differential Equations. + Society for Industrial and Applied Mathematics. + Journal on Numerical Analysis, 32(3), pp 797-823, 1995. + doi: [https://doi.org/10.1137/0732037](https://doi.org/10.1137/0732037) + +When applied to a `SplitODEProblem` of the form + +``` +u'(t) = f1(u) + f2(u) +``` + +The default `IMEXEuler()` method uses an update of the form + +``` +unew = uold + dt * (f1(unew) + f2(uold)) +``` + +See also `SBDF`, `IMEXEulerARK`. +""" +IMEXEuler(; kwargs...) = SBDF(1; kwargs...) + +""" + IMEXEulerARK(;kwargs...) + +The one-step version of the IMEX multistep methods of + + - Uri M. Ascher, Steven J. Ruuth, Brian T. R. Wetton. + Implicit-Explicit Methods for Time-Dependent Partial Differential Equations. + Society for Industrial and Applied Mathematics. + Journal on Numerical Analysis, 32(3), pp 797-823, 1995. + doi: [https://doi.org/10.1137/0732037](https://doi.org/10.1137/0732037) + +When applied to a `SplitODEProblem` of the form + +``` +u'(t) = f1(u) + f2(u) +``` + +A classical additive Runge-Kutta method in the sense of +[Araújo, Murua, Sanz-Serna (1997)](https://doi.org/10.1137/S0036142995292128) +consisting of the implicit and the explicit Euler method given by + +``` +y1 = uold + dt * f1(y1) +unew = uold + dt * (f1(unew) + f2(y1)) +``` + +See also `SBDF`, `IMEXEuler`. +""" +IMEXEulerARK(; kwargs...) = SBDF(1; ark = true, kwargs...) diff --git a/src/caches/bdf_caches.jl b/lib/OrdinaryDiffEqBDF/src/bdf_caches.jl similarity index 100% rename from src/caches/bdf_caches.jl rename to lib/OrdinaryDiffEqBDF/src/bdf_caches.jl diff --git a/src/perform_step/bdf_perform_step.jl b/lib/OrdinaryDiffEqBDF/src/bdf_perform_step.jl similarity index 100% rename from src/perform_step/bdf_perform_step.jl rename to lib/OrdinaryDiffEqBDF/src/bdf_perform_step.jl diff --git a/src/bdf_utils.jl b/lib/OrdinaryDiffEqBDF/src/bdf_utils.jl similarity index 100% rename from src/bdf_utils.jl rename to lib/OrdinaryDiffEqBDF/src/bdf_utils.jl diff --git a/lib/OrdinaryDiffEqBDF/src/controllers.jl b/lib/OrdinaryDiffEqBDF/src/controllers.jl new file mode 100644 index 0000000000..4c1e8b9274 --- /dev/null +++ b/lib/OrdinaryDiffEqBDF/src/controllers.jl @@ -0,0 +1,264 @@ +struct DummyController <: AbstractController +end + +function default_controller(alg::Union{QNDF, FBDF}, args...) + DummyController() +end + +# QNBDF +stepsize_controller!(integrator, alg::QNDF) = nothing + +# this stepsize and order controller is taken from +# Implementation of an Adaptive BDF2 Formula and Comparison with the MATLAB Ode15s paper +# E. Alberdi Celaya, J. J. Anza Aguirrezabala, and P. Chatzipantelidis + +function step_accept_controller!(integrator, alg::QNDF{max_order}, q) where {max_order} + #step is accepted, reset count of consecutive failed steps + integrator.cache.consfailcnt = 0 + integrator.cache.nconsteps += 1 + if iszero(integrator.EEst) + return integrator.dt * integrator.opts.qmax + else + est = integrator.EEst + estₖ₋₁ = integrator.cache.EEst1 + estₖ₊₁ = integrator.cache.EEst2 + h = integrator.dt + k = integrator.cache.order + cache = integrator.cache + prefer_const_step = integrator.cache.nconsteps < integrator.cache.order + 2 + zₛ = 1.2 # equivalent to integrator.opts.gamma + zᵤ = 0.1 + Fᵤ = 10 + expo = 1 / (k + 1) + z = zₛ * ((est)^expo) + F = inv(z) + hₙ = h + kₙ = k + if z <= zᵤ + hₖ = Fᵤ * h + else + hₖ = F * h + end + hₖ₋₁ = 0.0 + hₖ₊₁ = 0.0 + + if k > 1 + expo = 1 / k + zₖ₋₁ = 1.3 * ((estₖ₋₁)^expo) + Fₖ₋₁ = inv(zₖ₋₁) + if zₖ₋₁ <= 0.1 + hₖ₋₁ = 10 * h + elseif 1 / 10 < zₖ₋₁ <= 1.3 + hₖ₋₁ = Fₖ₋₁ * h + end + if hₖ₋₁ > hₖ + hₙ = hₖ₋₁ + kₙ = k - 1 + else + hₙ = hₖ + kₙ = k + end + else + hₙ = hₖ + kₙ = k + end + + if k < max_order + expo = 1 / (k + 2) + zₖ₊₁ = 1.4 * ((estₖ₊₁)^expo) + Fₖ₊₁ = inv(zₖ₊₁) + + if zₖ₊₁ <= 0.1 + hₖ₊₁ = 10 * h + elseif 0.1 < zₖ₊₁ <= 1.4 + hₖ₊₁ = Fₖ₊₁ * h + end + if hₖ₊₁ > hₙ + hₙ = hₖ₊₁ + kₙ = k + 1 + end + end + cache.order = kₙ + q = integrator.dt / hₙ + end + if prefer_const_step + if q < 1.2 && q > 0.6 + return integrator.dt + end + end + if q <= integrator.opts.qsteady_max && q >= integrator.opts.qsteady_min + return integrator.dt + end + return integrator.dt / q +end + +function step_reject_controller!(integrator, ::QNDF) + bdf_step_reject_controller!(integrator, integrator.cache.EEst1) +end + +function step_reject_controller!(integrator, ::FBDF) + bdf_step_reject_controller!(integrator, integrator.cache.terkm1) +end + +function bdf_step_reject_controller!(integrator, EEst1) + k = integrator.cache.order + h = integrator.dt + integrator.cache.consfailcnt += 1 + integrator.cache.nconsteps = 0 + if integrator.cache.consfailcnt > 1 + h = h / 2 + end + zₛ = 1.2 # equivalent to integrator.opts.gamma + expo = 1 / (k + 1) + z = zₛ * ((integrator.EEst)^expo) + F = inv(z) + if z <= 10 + hₖ = F * h + else # z > 10 + hₖ = 0.1 * h + end + hₙ = hₖ + kₙ = k + if k > 1 + expo = 1 / k + zₖ₋₁ = 1.3 * (EEst1^expo) + Fₖ₋₁ = inv(zₖ₋₁) + if zₖ₋₁ <= 10 + hₖ₋₁ = Fₖ₋₁ * h + elseif zₖ₋₁ > 10 + hₖ₋₁ = 0.1 * h + end + if integrator.cache.consfailcnt > 2 || hₖ₋₁ > hₖ + hₙ = min(h, hₖ₋₁) + kₙ = k - 1 + end + end + # Restart BDf (clear history) when we failed repeatedly + if kₙ == 1 && integrator.cache.consfailcnt > 3 + u_modified!(integrator, true) + end + integrator.dt = hₙ + integrator.cache.order = kₙ +end + +function post_newton_controller!(integrator, alg::FBDF) + @unpack cache = integrator + if cache.order > 1 && cache.nlsolver.nfails >= 3 + cache.order -= 1 + end + integrator.dt = integrator.dt / integrator.opts.failfactor + integrator.cache.consfailcnt += 1 + integrator.cache.nconsteps = 0 + nothing +end + +function choose_order!(alg::FBDF, integrator, + cache::OrdinaryDiffEqMutableCache, + ::Val{max_order}) where {max_order} + @unpack t, dt, u, cache, uprev = integrator + @unpack atmp, ts_tmp, terkm2, terkm1, terk, terkp1, terk_tmp, u_history = cache + k = cache.order + # only when the order of amount of terk follows the order of step size, and achieve enough constant step size, the order could be increased. + if k < max_order && integrator.cache.nconsteps >= integrator.cache.order + 2 && + ((k == 1 && terk > terkp1) || + (k == 2 && terkm1 > terk > terkp1) || + (k > 2 && terkm2 > terkm1 > terk > terkp1)) + k += 1 + terk = terkp1 + else + while !(terkm2 > terkm1 > terk > terkp1) && k > 2 + terkp1 = terk + terk = terkm1 + terkm1 = terkm2 + fd_weights = calc_finite_difference_weights(ts_tmp, t + dt, k - 2, + Val(max_order)) + terk_tmp = @.. broadcast=false fd_weights[k - 2, 1]*u + vc = _vec(terk_tmp) + for i in 2:(k - 2) + @.. broadcast=false @views vc += fd_weights[i, k - 2] * u_history[:, i - 1] + end + @.. broadcast=false terk_tmp*=abs(dt^(k - 2)) + calculate_residuals!(atmp, _vec(terk_tmp), _vec(uprev), _vec(u), + integrator.opts.abstol, integrator.opts.reltol, + integrator.opts.internalnorm, t) + terkm2 = integrator.opts.internalnorm(atmp, t) + k -= 1 + end + end + return k, terk +end + +function choose_order!(alg::FBDF, integrator, + cache::OrdinaryDiffEqConstantCache, + ::Val{max_order}) where {max_order} + @unpack t, dt, u, cache, uprev = integrator + @unpack ts_tmp, terkm2, terkm1, terk, terkp1, u_history = cache + k = cache.order + if k < max_order && integrator.cache.nconsteps >= integrator.cache.order + 2 && + ((k == 1 && terk > terkp1) || + (k == 2 && terkm1 > terk > terkp1) || + (k > 2 && terkm2 > terkm1 > terk > terkp1)) + k += 1 + terk = terkp1 + else + while !(terkm2 > terkm1 > terk > terkp1) && k > 2 + terkp1 = terk + terk = terkm1 + terkm1 = terkm2 + fd_weights = calc_finite_difference_weights(ts_tmp, t + dt, k - 2, + Val(max_order)) + terk_tmp = @.. broadcast=false fd_weights[k - 2, 1]*u + if u isa Number + for i in 2:(k - 2) + terk_tmp += fd_weights[i, k - 2] * u_history[i - 1] + end + terk_tmp *= abs(dt^(k - 2)) + else + vc = _vec(terk_tmp) + for i in 2:(k - 2) + @.. broadcast=false @views vc += fd_weights[i, k - 2] * + u_history[:, i - 1] + end + terk_tmp = reshape(vc, size(terk_tmp)) + terk_tmp *= @.. broadcast=false abs(dt^(k - 2)) + end + atmp = calculate_residuals(_vec(terk_tmp), _vec(uprev), _vec(u), + integrator.opts.abstol, integrator.opts.reltol, + integrator.opts.internalnorm, t) + terkm2 = integrator.opts.internalnorm(atmp, t) + k -= 1 + end + end + return k, terk +end + +function stepsize_controller!(integrator, + alg::FBDF{max_order}) where { + max_order, +} + @unpack cache = integrator + cache.prev_order = cache.order + k, terk = choose_order!(alg, integrator, cache, Val(max_order)) + if k != cache.order + integrator.cache.nconsteps = 0 + cache.order = k + end + if iszero(terk) + q = inv(integrator.opts.qmax) + else + q = ((2 * terk / (k + 1))^(1 / (k + 1))) + end + integrator.qold = q + q +end + +function step_accept_controller!(integrator, alg::FBDF{max_order}, + q) where {max_order} + integrator.cache.consfailcnt = 0 + if q <= integrator.opts.qsteady_max && q >= integrator.opts.qsteady_min + q = one(q) + end + integrator.cache.nconsteps += 1 + integrator.cache.iters_from_event += 1 + return integrator.dt / q +end diff --git a/lib/OrdinaryDiffEqSDIRK/Project.toml b/lib/OrdinaryDiffEqSDIRK/Project.toml new file mode 100644 index 0000000000..7d39bed0cb --- /dev/null +++ b/lib/OrdinaryDiffEqSDIRK/Project.toml @@ -0,0 +1,26 @@ +name = "OrdinaryDiffEqSDIRK" +uuid = "2d112036-d095-4a1e-ab9a-08536f3ecdbf" +authors = ["ParamThakkar123 "] +version = "1.0.0" + +[deps] +OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +TruncatedStacktraces = "781d530d-4396-4725-bb49-402e4bee1e77" +MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" +MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" +FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898" +RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" +SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" + +[compat] +julia = "1.10" + +[extras] +DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["DiffEqDevTools", "Random", "SafeTestsets", "Test"] \ No newline at end of file diff --git a/lib/OrdinaryDiffEqSDIRK/src/OrdinaryDiffEqSDIRK.jl b/lib/OrdinaryDiffEqSDIRK/src/OrdinaryDiffEqSDIRK.jl new file mode 100644 index 0000000000..30c0bf3b7b --- /dev/null +++ b/lib/OrdinaryDiffEqSDIRK/src/OrdinaryDiffEqSDIRK.jl @@ -0,0 +1,38 @@ +module OrdinaryDiffEqSDIRK + +import OrdinaryDiffEq: alg_order, calculate_residuals!, + initialize!, perform_step!, @unpack, unwrap_alg, + calculate_residuals, alg_extrapolates, + OrdinaryDiffEqAlgorithm, + OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, + OrdinaryDiffEqNewtonAdaptiveAlgorithm, + OrdinaryDiffEqNewtonAlgorithm, + DEFAULT_PRECS, NLNewton, + OrdinaryDiffEqAdaptiveAlgorithm, CompiledFloats, uses_uprev, + alg_cache, _vec, _reshape, @cache, isfsal, full_cache, + constvalue, _unwrap_val, du_alias_or_new, _ode_interpolant, + markfirststage!, trivial_limiter!, _ode_interpolant!, + UJacobianWrapper, set_new_W!, dolinsolve, get_W, + build_nlsolver, nlsolve!, nlsolvefail, isnewton, + COEFFICIENT_MULTISTEP, mul!, isesdirk, issplit, + ssp_coefficient +using TruncatedStacktraces, MuladdMacro, MacroTools, FastBroadcast, RecursiveArrayTools +using SciMLBase: SplitFunction +using LinearAlgebra: I + +include("algorithms.jl") +include("alg_utils.jl") +include("sdirk_caches.jl") +include("kencarp_kvaerno_caches.jl") +include("sdirk_perform_step.jl") +include("kencarp_kvaerno_perform_step.jl") +include("sdirk_tableaus.jl") + +export ImplicitEuler, ImplicitMidpoint, Trapezoid, TRBDF2, SDIRK2, SDIRK22, + Kvaerno3, KenCarp3, Cash4, Hairer4, Hairer42, SSPSDIRK2, Kvaerno4, + Kvaerno5, KenCarp4, KenCarp47, KenCarp5, KenCarp58, ESDIRK54I8L2SA, SFSDIRK4, + SFSDIRK5, CFNLIRK3, SFSDIRK6, SFSDIRK7, SFSDIRK8, Kvaerno5, KenCarp4, KenCarp5, + SFSDIRK4, SFSDIRK5, CFNLIRK3, SFSDIRK6, + SFSDIRK7, SFSDIRK8, ESDIRK436L2SA2, ESDIRK437L2SA, ESDIRK547L2SA2, ESDIRK659L2SA + +end diff --git a/lib/OrdinaryDiffEqSDIRK/src/alg_utils.jl b/lib/OrdinaryDiffEqSDIRK/src/alg_utils.jl new file mode 100644 index 0000000000..a4ac5fee51 --- /dev/null +++ b/lib/OrdinaryDiffEqSDIRK/src/alg_utils.jl @@ -0,0 +1,55 @@ +alg_extrapolates(alg::ImplicitEuler) = true +alg_extrapolates(alg::Trapezoid) = true +alg_extrapolates(alg::SDIRK22) = true + +alg_order(alg::Trapezoid) = 2 +alg_order(alg::ImplicitEuler) = 1 +alg_order(alg::ImplicitMidpoint) = 2 +alg_order(alg::TRBDF2) = 2 +alg_order(alg::SSPSDIRK2) = 2 +alg_order(alg::SDIRK2) = 2 +alg_order(alg::SDIRK22) = 2 +alg_order(alg::Kvaerno3) = 3 +alg_order(alg::Kvaerno4) = 4 +alg_order(alg::Kvaerno5) = 5 +alg_order(alg::ESDIRK54I8L2SA) = 5 +alg_order(alg::ESDIRK436L2SA2) = 4 +alg_order(alg::ESDIRK437L2SA) = 4 +alg_order(alg::ESDIRK547L2SA2) = 5 +alg_order(alg::ESDIRK659L2SA) = 6 +alg_order(alg::KenCarp3) = 3 +alg_order(alg::CFNLIRK3) = 3 +alg_order(alg::KenCarp4) = 4 +alg_order(alg::KenCarp47) = 4 +alg_order(alg::KenCarp5) = 5 +alg_order(alg::KenCarp58) = 5 +alg_order(alg::Cash4) = 4 +alg_order(alg::SFSDIRK4) = 4 +alg_order(alg::SFSDIRK5) = 4 +alg_order(alg::SFSDIRK6) = 4 +alg_order(alg::SFSDIRK7) = 4 +alg_order(alg::SFSDIRK8) = 4 +alg_order(alg::Hairer4) = 4 +alg_order(alg::Hairer42) = 4 + +function isesdirk(alg::Union{KenCarp3, KenCarp4, KenCarp5, KenCarp58, + Kvaerno3, Kvaerno4, Kvaerno5, ESDIRK437L2SA, + ESDIRK54I8L2SA, ESDIRK436L2SA2, ESDIRK547L2SA2, + ESDIRK659L2SA, CFNLIRK3}) + true +end + +alg_adaptive_order(alg::Trapezoid) = 1 +alg_adaptive_order(alg::ImplicitMidpoint) = 1 +alg_adaptive_order(alg::ImplicitEuler) = 0 + +ssp_coefficient(alg::SSPSDIRK2) = 4 + +isesdirk(alg::TRBDF2) = true + +issplit(alg::KenCarp3) = true +issplit(alg::KenCarp4) = true +issplit(alg::KenCarp47) = true +issplit(alg::KenCarp5) = true +issplit(alg::KenCarp58) = true +issplit(alg::CFNLIRK3) = true diff --git a/lib/OrdinaryDiffEqSDIRK/src/algorithms.jl b/lib/OrdinaryDiffEqSDIRK/src/algorithms.jl new file mode 100644 index 0000000000..01f796f650 --- /dev/null +++ b/lib/OrdinaryDiffEqSDIRK/src/algorithms.jl @@ -0,0 +1,862 @@ +# SDIRK Methods +""" +ImplicitEuler: SDIRK Method +A 1st order implicit solver. A-B-L-stable. Adaptive timestepping through a divided differences estimate via memory. +Strong-stability preserving (SSP). +""" +struct ImplicitEuler{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <: + OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + extrapolant::Symbol + controller::Symbol + step_limiter!::StepLimiter +end + +function ImplicitEuler(; 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 = :constant, + controller = :PI, step_limiter! = trivial_limiter!) + ImplicitEuler{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), + _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, + nlsolve, precs, extrapolant, controller, step_limiter!) +end +""" +ImplicitMidpoint: SDIRK Method +A second order A-stable symplectic and symmetric implicit solver. +Good for highly stiff equations which need symplectic integration. +""" +struct ImplicitMidpoint{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <: + OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + extrapolant::Symbol + step_limiter!::StepLimiter +end + +function ImplicitMidpoint(; 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, step_limiter! = trivial_limiter!) + ImplicitMidpoint{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), + _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, + nlsolve, + precs, + extrapolant, + step_limiter!) +end + +""" +Andre Vladimirescu. 1994. The Spice Book. John Wiley & Sons, Inc., New York, +NY, USA. + +Trapezoid: SDIRK Method +A second order A-stable symmetric ESDIRK method. +"Almost symplectic" without numerical dampening. +Also known as Crank-Nicolson when applied to PDEs. Adaptive timestepping via divided +differences approximation to the second derivative terms in the local truncation error +estimate (the SPICE approximation strategy). +""" +struct Trapezoid{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <: + OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + extrapolant::Symbol + controller::Symbol + step_limiter!::StepLimiter +end + +function Trapezoid(; 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, + controller = :PI, step_limiter! = trivial_limiter!) + Trapezoid{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), + _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, + nlsolve, + precs, + extrapolant, + controller, + step_limiter!) +end + +""" +@article{hosea1996analysis, +title={Analysis and implementation of TR-BDF2}, +author={Hosea, ME and Shampine, LF}, +journal={Applied Numerical Mathematics}, +volume={20}, +number={1-2}, +pages={21--37}, +year={1996}, +publisher={Elsevier} +} + +TRBDF2: SDIRK Method +A second order A-B-L-S-stable one-step ESDIRK method. +Includes stiffness-robust error estimates for accurate adaptive timestepping, smoothed derivatives for highly stiff and oscillatory problems. +""" +struct TRBDF2{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <: + OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + smooth_est::Bool + extrapolant::Symbol + controller::Symbol + step_limiter!::StepLimiter +end + +function TRBDF2(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(), + concrete_jac = nothing, diff_type = Val{:forward}, + linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), + smooth_est = true, extrapolant = :linear, + controller = :PI, step_limiter! = trivial_limiter!) + TRBDF2{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), + _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, nlsolve, precs, + smooth_est, extrapolant, controller, step_limiter!) +end + +TruncatedStacktraces.@truncate_stacktrace TRBDF2 + +""" +@article{hindmarsh2005sundials, +title={{SUNDIALS}: Suite of nonlinear and differential/algebraic equation solvers}, +author={Hindmarsh, Alan C and Brown, Peter N and Grant, Keith E and Lee, Steven L and Serban, Radu and Shumaker, Dan E and Woodward, Carol S}, +journal={ACM Transactions on Mathematical Software (TOMS)}, +volume={31}, +number={3}, +pages={363--396}, +year={2005}, +publisher={ACM} +} + +SDIRK2: SDIRK Method +An A-B-L stable 2nd order SDIRK method +""" +struct SDIRK2{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <: + OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + smooth_est::Bool + extrapolant::Symbol + controller::Symbol + step_limiter!::StepLimiter +end + +function SDIRK2(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(), + concrete_jac = nothing, diff_type = Val{:forward}, + linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), + smooth_est = true, extrapolant = :linear, + controller = :PI, step_limiter! = trivial_limiter!) + SDIRK2{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), + _unwrap_val(concrete_jac), typeof(step_limiter!)}( + linsolve, nlsolve, precs, smooth_est, extrapolant, + controller, + step_limiter!) +end + +struct SDIRK22{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <: + OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + extrapolant::Symbol + controller::Symbol + step_limiter!::StepLimiter +end + +function SDIRK22(; + 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, + controller = :PI, step_limiter! = trivial_limiter!) + Trapezoid{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), + _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, + nlsolve, + precs, + extrapolant, + controller, + step_limiter!) +end + +struct SSPSDIRK2{CS, AD, F, F2, P, FDT, ST, CJ} <: + OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ} # Not adaptive + linsolve::F + nlsolve::F2 + precs::P + smooth_est::Bool + extrapolant::Symbol + controller::Symbol +end + +function SSPSDIRK2(; chunk_size = Val{0}(), autodiff = Val{true}(), + standardtag = Val{true}(), concrete_jac = nothing, + diff_type = Val{:forward}, + linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), + smooth_est = true, extrapolant = :constant, + controller = :PI) + SSPSDIRK2{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), + _unwrap_val(concrete_jac)}(linsolve, nlsolve, precs, smooth_est, extrapolant, + controller) +end + +""" +@article{kvaerno2004singly, +title={Singly diagonally implicit Runge--Kutta methods with an explicit first stage}, +author={Kv{\\ae}rn{\\o}, Anne}, +journal={BIT Numerical Mathematics}, +volume={44}, +number={3}, +pages={489--502}, +year={2004}, +publisher={Springer} +} + +Kvaerno3: SDIRK Method +An A-L stable stiffly-accurate 3rd order ESDIRK method +""" +struct Kvaerno3{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <: + OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + smooth_est::Bool + extrapolant::Symbol + controller::Symbol + step_limiter!::StepLimiter +end +function Kvaerno3(; chunk_size = Val{0}(), autodiff = Val{true}(), + standardtag = Val{true}(), concrete_jac = nothing, + diff_type = Val{:forward}, + linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), + smooth_est = true, extrapolant = :linear, + controller = :PI, step_limiter! = trivial_limiter!) + Kvaerno3{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), + _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, nlsolve, precs, + smooth_est, extrapolant, controller, step_limiter!) +end + +""" +@book{kennedy2001additive, +title={Additive Runge-Kutta schemes for convection-diffusion-reaction equations}, +author={Kennedy, Christopher Alan}, +year={2001}, +publisher={National Aeronautics and Space Administration, Langley Research Center} +} + +KenCarp3: SDIRK Method +An A-L stable stiffly-accurate 3rd order ESDIRK method with splitting +""" +struct KenCarp3{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <: + OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + smooth_est::Bool + extrapolant::Symbol + controller::Symbol + step_limiter!::StepLimiter +end +function KenCarp3(; chunk_size = Val{0}(), autodiff = Val{true}(), + standardtag = Val{true}(), concrete_jac = nothing, + diff_type = Val{:forward}, + linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), + smooth_est = true, extrapolant = :linear, + controller = :PI, step_limiter! = trivial_limiter!) + KenCarp3{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), + _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, nlsolve, precs, + smooth_est, extrapolant, controller, step_limiter!) +end + +struct CFNLIRK3{CS, AD, F, F2, P, FDT, ST, CJ} <: + OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + extrapolant::Symbol +end +function CFNLIRK3(; 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) + CFNLIRK3{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), + _unwrap_val(concrete_jac)}(linsolve, + nlsolve, + precs, + extrapolant) +end + +""" +@article{hindmarsh2005sundials, +title={{SUNDIALS}: Suite of nonlinear and differential/algebraic equation solvers}, +author={Hindmarsh, Alan C and Brown, Peter N and Grant, Keith E and Lee, Steven L and Serban, Radu and Shumaker, Dan E and Woodward, Carol S}, +journal={ACM Transactions on Mathematical Software (TOMS)}, +volume={31}, +number={3}, +pages={363--396}, +year={2005}, +publisher={ACM} +} + +Cash4: SDIRK Method +An A-L stable 4th order SDIRK method +""" +struct Cash4{CS, AD, F, F2, P, FDT, ST, CJ} <: + OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + smooth_est::Bool + extrapolant::Symbol + embedding::Int + controller::Symbol +end +function Cash4(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(), + concrete_jac = nothing, diff_type = Val{:forward}, + linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), + smooth_est = true, extrapolant = :linear, + controller = :PI, embedding = 3) + Cash4{ + _unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve), + typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac)}( + linsolve, + nlsolve, + precs, + smooth_est, + extrapolant, + embedding, + controller) +end + +struct SFSDIRK4{CS, AD, F, F2, P, FDT, ST, CJ} <: + OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + extrapolant::Symbol +end +function SFSDIRK4(; 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) + SFSDIRK4{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), + _unwrap_val(concrete_jac)}(linsolve, + nlsolve, + precs, + extrapolant) +end + +struct SFSDIRK5{CS, AD, F, F2, P, FDT, ST, CJ} <: + OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + extrapolant::Symbol +end + +function SFSDIRK5(; 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) + SFSDIRK5{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), + _unwrap_val(concrete_jac)}(linsolve, + nlsolve, + precs, + extrapolant) +end + +struct SFSDIRK6{CS, AD, F, F2, P, FDT, ST, CJ} <: + OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + extrapolant::Symbol +end + +function SFSDIRK6(; 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) + SFSDIRK6{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), + _unwrap_val(concrete_jac)}(linsolve, + nlsolve, + precs, + extrapolant) +end + +struct SFSDIRK7{CS, AD, F, F2, P, FDT, ST, CJ} <: + OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + extrapolant::Symbol +end + +function SFSDIRK7(; 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) + SFSDIRK7{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), + _unwrap_val(concrete_jac)}(linsolve, + nlsolve, + precs, + extrapolant) +end + +struct SFSDIRK8{CS, AD, F, F2, P, FDT, ST, CJ} <: + OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + extrapolant::Symbol +end + +function SFSDIRK8(; 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) + SFSDIRK8{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), + _unwrap_val(concrete_jac)}(linsolve, + nlsolve, + precs, + extrapolant) +end + +""" +E. Hairer, G. Wanner, Solving ordinary differential equations II, stiff and +differential-algebraic problems. Computational mathematics (2nd revised ed.), +Springer (1996) + +Hairer4: SDIRK Method +An A-L stable 4th order SDIRK method +""" +struct Hairer4{CS, AD, F, F2, P, FDT, ST, CJ} <: + OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + smooth_est::Bool + extrapolant::Symbol + controller::Symbol +end +function Hairer4(; + chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(), + concrete_jac = nothing, diff_type = Val{:forward}, + linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), + smooth_est = true, extrapolant = :linear, + controller = :PI) + Hairer4{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), + _unwrap_val(concrete_jac)}(linsolve, nlsolve, precs, smooth_est, extrapolant, + controller) +end + +""" +E. Hairer, G. Wanner, Solving ordinary differential equations II, stiff and +differential-algebraic problems. Computational mathematics (2nd revised ed.), +Springer (1996) + +Hairer42: SDIRK Method +An A-L stable 4th order SDIRK method +""" +struct Hairer42{CS, AD, F, F2, P, FDT, ST, CJ} <: + OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + smooth_est::Bool + extrapolant::Symbol + controller::Symbol +end +function Hairer42(; chunk_size = Val{0}(), autodiff = Val{true}(), + standardtag = Val{true}(), concrete_jac = nothing, + diff_type = Val{:forward}, + linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), + smooth_est = true, extrapolant = :linear, + controller = :PI) + Hairer42{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), + _unwrap_val(concrete_jac)}(linsolve, nlsolve, precs, smooth_est, extrapolant, + controller) +end + +""" +@article{kvaerno2004singly, +title={Singly diagonally implicit Runge--Kutta methods with an explicit first stage}, +author={Kv{\\ae}rn{\\o}, Anne}, +journal={BIT Numerical Mathematics}, +volume={44}, +number={3}, +pages={489--502}, +year={2004}, +publisher={Springer} +} + +Kvaerno4: SDIRK Method +An A-L stable stiffly-accurate 4th order ESDIRK method. +""" +struct Kvaerno4{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <: + OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + smooth_est::Bool + extrapolant::Symbol + controller::Symbol + step_limiter!::StepLimiter +end +function Kvaerno4(; chunk_size = Val{0}(), autodiff = Val{true}(), + standardtag = Val{true}(), concrete_jac = nothing, + diff_type = Val{:forward}, + linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), + smooth_est = true, extrapolant = :linear, + controller = :PI, step_limiter! = trivial_limiter!) + Kvaerno4{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), + _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, nlsolve, precs, + smooth_est, extrapolant, controller, step_limiter!) +end + +""" +@article{kvaerno2004singly, +title={Singly diagonally implicit Runge--Kutta methods with an explicit first stage}, +author={Kv{\\ae}rn{\\o}, Anne}, +journal={BIT Numerical Mathematics}, +volume={44}, +number={3}, +pages={489--502}, +year={2004}, +publisher={Springer} +} + +Kvaerno5: SDIRK Method +An A-L stable stiffly-accurate 5th order ESDIRK method +""" +struct Kvaerno5{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <: + OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + smooth_est::Bool + extrapolant::Symbol + controller::Symbol + step_limiter!::StepLimiter +end +function Kvaerno5(; chunk_size = Val{0}(), autodiff = Val{true}(), + standardtag = Val{true}(), concrete_jac = nothing, + diff_type = Val{:forward}, + linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), + smooth_est = true, extrapolant = :linear, + controller = :PI, step_limiter! = trivial_limiter!) + Kvaerno5{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), + _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, nlsolve, precs, + smooth_est, extrapolant, controller, step_limiter!) +end + +""" +@book{kennedy2001additive, +title={Additive Runge-Kutta schemes for convection-diffusion-reaction equations}, +author={Kennedy, Christopher Alan}, +year={2001}, +publisher={National Aeronautics and Space Administration, Langley Research Center} +} + +KenCarp4: SDIRK Method +An A-L stable stiffly-accurate 4th order ESDIRK method with splitting +""" +struct KenCarp4{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <: + OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + smooth_est::Bool + extrapolant::Symbol + controller::Symbol + step_limiter!::StepLimiter +end +function KenCarp4(; chunk_size = Val{0}(), autodiff = Val{true}(), + standardtag = Val{true}(), concrete_jac = nothing, + diff_type = Val{:forward}, + linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), + smooth_est = true, extrapolant = :linear, + controller = :PI, step_limiter! = trivial_limiter!) + KenCarp4{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), + _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, nlsolve, precs, + smooth_est, extrapolant, controller, step_limiter!) +end + +TruncatedStacktraces.@truncate_stacktrace KenCarp4 + +""" +@article{kennedy2019higher, +title={Higher-order additive Runge--Kutta schemes for ordinary differential equations}, +author={Kennedy, Christopher A and Carpenter, Mark H}, +journal={Applied Numerical Mathematics}, +volume={136}, +pages={183--205}, +year={2019}, +publisher={Elsevier} +} + +KenCarp47: SDIRK Method +An A-L stable stiffly-accurate 4th order seven-stage ESDIRK method with splitting +""" +struct KenCarp47{CS, AD, F, F2, P, FDT, ST, CJ} <: + OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + smooth_est::Bool + extrapolant::Symbol + controller::Symbol +end +function KenCarp47(; chunk_size = Val{0}(), autodiff = Val{true}(), + standardtag = Val{true}(), concrete_jac = nothing, + diff_type = Val{:forward}, + linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), + smooth_est = true, extrapolant = :linear, + controller = :PI) + KenCarp47{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), + _unwrap_val(concrete_jac)}(linsolve, nlsolve, precs, smooth_est, extrapolant, + controller) +end + +""" +@book{kennedy2001additive, +title={Additive Runge-Kutta schemes for convection-diffusion-reaction equations}, +author={Kennedy, Christopher Alan}, +year={2001}, +publisher={National Aeronautics and Space Administration, Langley Research Center} +} + +KenCarp5: SDIRK Method +An A-L stable stiffly-accurate 5th order ESDIRK method with splitting +""" +struct KenCarp5{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <: + OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + smooth_est::Bool + extrapolant::Symbol + controller::Symbol + step_limiter!::StepLimiter +end +function KenCarp5(; chunk_size = Val{0}(), autodiff = Val{true}(), + standardtag = Val{true}(), concrete_jac = nothing, + diff_type = Val{:forward}, + linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), + smooth_est = true, extrapolant = :linear, + controller = :PI, step_limiter! = trivial_limiter!) + KenCarp5{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), + _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, nlsolve, precs, + smooth_est, extrapolant, controller, step_limiter!) +end +""" +@article{kennedy2019higher, +title={Higher-order additive Runge--Kutta schemes for ordinary differential equations}, +author={Kennedy, Christopher A and Carpenter, Mark H}, +journal={Applied Numerical Mathematics}, +volume={136}, +pages={183--205}, +year={2019}, +publisher={Elsevier} +} + +KenCarp58: SDIRK Method +An A-L stable stiffly-accurate 5th order eight-stage ESDIRK method with splitting +""" +struct KenCarp58{CS, AD, F, F2, P, FDT, ST, CJ} <: + OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + smooth_est::Bool + extrapolant::Symbol + controller::Symbol +end +function KenCarp58(; chunk_size = Val{0}(), autodiff = Val{true}(), + standardtag = Val{true}(), concrete_jac = nothing, + diff_type = Val{:forward}, + linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), + smooth_est = true, extrapolant = :linear, + controller = :PI) + KenCarp58{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), + _unwrap_val(concrete_jac)}(linsolve, nlsolve, precs, smooth_est, extrapolant, + controller) +end + +# `smooth_est` is not necessary, as the embedded method is also L-stable +struct ESDIRK54I8L2SA{CS, AD, F, F2, P, FDT, ST, CJ} <: + OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + extrapolant::Symbol + controller::Symbol +end +function ESDIRK54I8L2SA(; 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, controller = :PI) + ESDIRK54I8L2SA{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), + _unwrap_val(concrete_jac)}(linsolve, nlsolve, precs, extrapolant, + controller) +end + +""" +@article{Kennedy2019DiagonallyIR, +title={Diagonally implicit Runge–Kutta methods for stiff ODEs}, +author={Christopher A. Kennedy and Mark H. Carpenter}, +journal={Applied Numerical Mathematics}, +year={2019}, +volume={146}, +pages={221-244} +} +""" +struct ESDIRK436L2SA2{CS, AD, F, F2, P, FDT, ST, CJ} <: + OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + extrapolant::Symbol + controller::Symbol +end +function ESDIRK436L2SA2(; 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, controller = :PI) + ESDIRK436L2SA2{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), + _unwrap_val(concrete_jac)}(linsolve, nlsolve, precs, extrapolant, + controller) +end + +""" +@article{Kennedy2019DiagonallyIR, +title={Diagonally implicit Runge–Kutta methods for stiff ODEs}, +author={Christopher A. Kennedy and Mark H. Carpenter}, +journal={Applied Numerical Mathematics}, +year={2019}, +volume={146}, +pages={221-244} +} +""" +struct ESDIRK437L2SA{CS, AD, F, F2, P, FDT, ST, CJ} <: + OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + extrapolant::Symbol + controller::Symbol +end +function ESDIRK437L2SA(; 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, controller = :PI) + ESDIRK437L2SA{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), + _unwrap_val(concrete_jac)}(linsolve, nlsolve, precs, extrapolant, + controller) +end + +""" +@article{Kennedy2019DiagonallyIR, +title={Diagonally implicit Runge–Kutta methods for stiff ODEs}, +author={Christopher A. Kennedy and Mark H. Carpenter}, +journal={Applied Numerical Mathematics}, +year={2019}, +volume={146}, +pages={221-244} +} +""" +struct ESDIRK547L2SA2{CS, AD, F, F2, P, FDT, ST, CJ} <: + OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + extrapolant::Symbol + controller::Symbol +end +function ESDIRK547L2SA2(; 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, controller = :PI) + ESDIRK547L2SA2{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), + _unwrap_val(concrete_jac)}(linsolve, nlsolve, precs, extrapolant, + controller) +end + +""" +@article{Kennedy2019DiagonallyIR, +title={Diagonally implicit Runge–Kutta methods for stiff ODEs}, +author={Christopher A. Kennedy and Mark H. Carpenter}, +journal={Applied Numerical Mathematics}, +year={2019}, +volume={146}, +pages={221-244} + +Currently has STABILITY ISSUES, causing it to fail the adaptive tests. +Check issue https://github.com/SciML/OrdinaryDiffEq.jl/issues/1933 for more details. +} +""" +struct ESDIRK659L2SA{CS, AD, F, F2, P, FDT, ST, CJ} <: + OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + extrapolant::Symbol + controller::Symbol +end +function ESDIRK659L2SA(; 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, controller = :PI) + ESDIRK659L2SA{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), + _unwrap_val(concrete_jac)}(linsolve, nlsolve, precs, extrapolant, + controller) +end diff --git a/src/caches/kencarp_kvaerno_caches.jl b/lib/OrdinaryDiffEqSDIRK/src/kencarp_kvaerno_caches.jl similarity index 100% rename from src/caches/kencarp_kvaerno_caches.jl rename to lib/OrdinaryDiffEqSDIRK/src/kencarp_kvaerno_caches.jl diff --git a/src/perform_step/kencarp_kvaerno_perform_step.jl b/lib/OrdinaryDiffEqSDIRK/src/kencarp_kvaerno_perform_step.jl similarity index 100% rename from src/perform_step/kencarp_kvaerno_perform_step.jl rename to lib/OrdinaryDiffEqSDIRK/src/kencarp_kvaerno_perform_step.jl diff --git a/src/caches/sdirk_caches.jl b/lib/OrdinaryDiffEqSDIRK/src/sdirk_caches.jl similarity index 100% rename from src/caches/sdirk_caches.jl rename to lib/OrdinaryDiffEqSDIRK/src/sdirk_caches.jl diff --git a/src/perform_step/sdirk_perform_step.jl b/lib/OrdinaryDiffEqSDIRK/src/sdirk_perform_step.jl similarity index 100% rename from src/perform_step/sdirk_perform_step.jl rename to lib/OrdinaryDiffEqSDIRK/src/sdirk_perform_step.jl diff --git a/src/tableaus/sdirk_tableaus.jl b/lib/OrdinaryDiffEqSDIRK/src/sdirk_tableaus.jl similarity index 100% rename from src/tableaus/sdirk_tableaus.jl rename to lib/OrdinaryDiffEqSDIRK/src/sdirk_tableaus.jl diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index b832c88f8c..b674894e81 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -142,22 +142,18 @@ include("nlsolve/type.jl") include("nlsolve/utils.jl") include("nlsolve/nlsolve.jl") include("nlsolve/functional.jl") -include("nlsolve/newton.jl") include("generic_rosenbrock.jl") include("caches/basic_caches.jl") include("caches/low_order_rk_caches.jl") include("caches/high_order_rk_caches.jl") -include("caches/sdirk_caches.jl") include("caches/firk_caches.jl") -include("caches/kencarp_kvaerno_caches.jl") include("caches/linear_caches.jl") include("caches/linear_nonlinear_caches.jl") include("caches/rosenbrock_caches.jl") include("caches/adams_bashforth_moulton_caches.jl") include("caches/nordsieck_caches.jl") -include("caches/bdf_caches.jl") include("caches/prk_caches.jl") include("caches/pdirk_caches.jl") include("caches/dae_caches.jl") @@ -166,7 +162,6 @@ include("caches/qprk_caches.jl") include("tableaus/low_order_rk_tableaus.jl") include("tableaus/high_order_rk_tableaus.jl") include("tableaus/rosenbrock_tableaus.jl") -include("tableaus/sdirk_tableaus.jl") include("tableaus/firk_tableaus.jl") include("tableaus/qprk_tableaus.jl") @@ -186,17 +181,13 @@ include("perform_step/exponential_rk_perform_step.jl") include("perform_step/explicit_rk_perform_step.jl") include("perform_step/low_order_rk_perform_step.jl") include("perform_step/high_order_rk_perform_step.jl") -include("perform_step/sdirk_perform_step.jl") -include("perform_step/kencarp_kvaerno_perform_step.jl") include("perform_step/firk_perform_step.jl") include("perform_step/rosenbrock_perform_step.jl") include("perform_step/composite_perform_step.jl") include("perform_step/adams_bashforth_moulton_perform_step.jl") include("perform_step/nordsieck_perform_step.jl") -include("perform_step/bdf_perform_step.jl") include("perform_step/prk_perform_step.jl") include("perform_step/pdirk_perform_step.jl") -include("perform_step/dae_perform_step.jl") include("perform_step/qprk_perform_step.jl") include("dense/generic_dense.jl") @@ -209,7 +200,6 @@ include("dense/high_order_rk_addsteps.jl") include("derivative_utils.jl") include("nordsieck_utils.jl") include("adams_utils.jl") -include("bdf_utils.jl") include("derivative_wrappers.jl") include("iterator_interface.jl") include("constants.jl") @@ -275,10 +265,33 @@ include("../lib/OrdinaryDiffEqVerner/src/OrdinaryDiffEqVerner.jl") using ..OrdinaryDiffEqVerner export Vern6, Vern7, Vern8, Vern9 +include("../lib/OrdinaryDiffEqSDIRK/src/OrdinaryDiffEqSDIRK.jl") +using ..OrdinaryDiffEqSDIRK +import .OrdinaryDiffEqSDIRK: ImplicitEulerConstantCache, ImplicitEulerCache +export ImplicitEuler, ImplicitMidpoint, Trapezoid, TRBDF2, SDIRK2, SDIRK22, + Kvaerno3, KenCarp3, Cash4, Hairer4, Hairer42, SSPSDIRK2, Kvaerno4, + Kvaerno5, KenCarp4, KenCarp47, KenCarp5, KenCarp58, ESDIRK54I8L2SA, SFSDIRK4, + SFSDIRK5, CFNLIRK3, SFSDIRK6, SFSDIRK7, SFSDIRK8, Kvaerno5, KenCarp4, KenCarp5, + SFSDIRK4, SFSDIRK5, CFNLIRK3, SFSDIRK6, + SFSDIRK7, SFSDIRK8, ESDIRK436L2SA2, ESDIRK437L2SA, ESDIRK547L2SA2, ESDIRK659L2SA + +include("../lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl") +using ..OrdinaryDiffEqBDF +export ABDF2, QNDF1, QBDF1, QNDF2, QBDF2, QNDF, QBDF, FBDF, + SBDF2, SBDF3, SBDF4, MEBDF2, IMEXEuler, IMEXEulerARK + include("../lib/OrdinaryDiffEqDefault/src/OrdinaryDiffEqDefault.jl") using ..OrdinaryDiffEqDefault export DefaultODEAlgorithm +using ..OrdinaryDiffEqBDF: reinitFBDF!, error_constant, estimate_terk!, + calc_Lagrange_interp!, + calc_finite_difference_weights, estimate_terk, + calc_Lagrange_interp, + bdf_step_reject_controller! +include("nlsolve/newton.jl") +include("perform_step/dae_perform_step.jl") + import PrecompileTools PrecompileTools.@compile_workload begin @@ -417,13 +430,6 @@ export FunctionMap, Euler, Heun, Ralston, Midpoint, RK4, ExplicitRK, OwrenZen3, export RadauIIA3, RadauIIA5, RadauIIA9 -export ImplicitEuler, ImplicitMidpoint, Trapezoid, TRBDF2, SDIRK2, SDIRK22, - Kvaerno3, KenCarp3, Cash4, Hairer4, Hairer42, SSPSDIRK2, Kvaerno4, - Kvaerno5, KenCarp4, KenCarp47, KenCarp5, KenCarp58, ESDIRK54I8L2SA, SFSDIRK4, - SFSDIRK5, CFNLIRK3, SFSDIRK6, SFSDIRK7, SFSDIRK8, Kvaerno5, KenCarp4, KenCarp5, - SFSDIRK4, SFSDIRK5, CFNLIRK3, SFSDIRK6, - SFSDIRK7, SFSDIRK8, ESDIRK436L2SA2, ESDIRK437L2SA, ESDIRK547L2SA2, ESDIRK659L2SA - export MagnusMidpoint, LinearExponential, MagnusLeapfrog, LieEuler, CayleyEuler, MagnusGauss4, MagnusNC6, MagnusGL6, MagnusGL8, MagnusNC8, MagnusGL4, MagnusAdapt4, RKMK2, RKMK4, LieRK4, CG2, CG3, CG4a @@ -449,16 +455,10 @@ export VCAB3, VCAB4, VCAB5, VCABM3, VCABM4, VCABM5 export VCABM -export IMEXEuler, IMEXEulerARK, CNAB2, CNLF2 +export CNAB2, CNLF2 export AN5, JVODE, JVODE_Adams, JVODE_BDF -export ABDF2, QNDF1, QBDF1, QNDF2, QBDF2, QNDF, QBDF, FBDF - -export SBDF2, SBDF3, SBDF4 - -export MEBDF2 - export Alshina2, Alshina3, Alshina6 export AutoSwitch, AutoTsit5, AutoDP5, diff --git a/src/alg_utils.jl b/src/alg_utils.jl index 3c4e3ece28..d2eab67539 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -1,5 +1,4 @@ ## SciMLBase Trait Definitions - function SciMLBase.isautodifferentiable(alg::Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm, FunctionMap}) true @@ -347,14 +346,8 @@ end alg_extrapolates(alg::Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}) = false alg_extrapolates(alg::CompositeAlgorithm) = any(alg_extrapolates.(alg.algs)) -alg_extrapolates(alg::ImplicitEuler) = true alg_extrapolates(alg::DImplicitEuler) = true alg_extrapolates(alg::DABDF2) = true -alg_extrapolates(alg::Trapezoid) = true -alg_extrapolates(alg::SDIRK22) = true -alg_extrapolates(alg::ABDF2) = true -alg_extrapolates(alg::SBDF) = true -alg_extrapolates(alg::MEBDF2) = true alg_extrapolates(alg::MagnusLeapfrog) = true function alg_order(alg::Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}) @@ -374,10 +367,6 @@ function get_current_adaptive_order(alg::OrdinaryDiffEqAdamsVarOrderVarStepAlgor cache.order end get_current_alg_order(alg::JVODE, cache) = get_current_adaptive_order(alg, cache) -get_current_alg_order(alg::QNDF, cache) = cache.order -get_current_alg_order(alg::FBDF, cache) = cache.order -get_current_adaptive_order(alg::QNDF, cache) = cache.order -get_current_adaptive_order(alg::FBDF, cache) = cache.order #alg_adaptive_order(alg::OrdinaryDiffEqAdaptiveAlgorithm) = error("Algorithm is adaptive with no order") function get_current_adaptive_order(alg::Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}, @@ -459,34 +448,6 @@ alg_order(alg::MagnusGL4) = 4 alg_order(alg::MagnusAdapt4) = 4 alg_order(alg::LinearExponential) = 1 alg_order(alg::MagnusLeapfrog) = 2 -alg_order(alg::Trapezoid) = 2 -alg_order(alg::ImplicitMidpoint) = 2 -alg_order(alg::TRBDF2) = 2 -alg_order(alg::SSPSDIRK2) = 2 -alg_order(alg::SDIRK2) = 2 -alg_order(alg::SDIRK22) = 2 -alg_order(alg::Kvaerno3) = 3 -alg_order(alg::Kvaerno4) = 4 -alg_order(alg::Kvaerno5) = 5 -alg_order(alg::ESDIRK54I8L2SA) = 5 -alg_order(alg::ESDIRK436L2SA2) = 4 -alg_order(alg::ESDIRK437L2SA) = 4 -alg_order(alg::ESDIRK547L2SA2) = 5 -alg_order(alg::ESDIRK659L2SA) = 6 -alg_order(alg::KenCarp3) = 3 -alg_order(alg::CFNLIRK3) = 3 -alg_order(alg::KenCarp4) = 4 -alg_order(alg::KenCarp47) = 4 -alg_order(alg::KenCarp5) = 5 -alg_order(alg::KenCarp58) = 5 -alg_order(alg::Cash4) = 4 -alg_order(alg::SFSDIRK4) = 4 -alg_order(alg::SFSDIRK5) = 4 -alg_order(alg::SFSDIRK6) = 4 -alg_order(alg::SFSDIRK7) = 4 -alg_order(alg::SFSDIRK8) = 4 -alg_order(alg::Hairer4) = 4 -alg_order(alg::Hairer42) = 4 alg_order(alg::PFRK87) = 8 alg_order(alg::ROS2) = 2 @@ -547,16 +508,6 @@ alg_order(alg::CNLF2) = 2 alg_order(alg::AN5) = 5 alg_order(alg::JVODE) = 1 #dummy value -alg_order(alg::ABDF2) = 2 -alg_order(alg::QNDF1) = 1 -alg_order(alg::QNDF2) = 2 - -alg_order(alg::QNDF) = 1 #dummy value -alg_order(alg::FBDF) = 1 #dummy value - -alg_order(alg::SBDF) = alg.order - -alg_order(alg::MEBDF2) = 2 alg_order(alg::PDIRK44) = 4 alg_order(alg::DImplicitEuler) = 1 @@ -582,11 +533,8 @@ alg_adaptive_order(alg::RadauIIA3) = 1 alg_adaptive_order(alg::RadauIIA5) = 3 alg_adaptive_order(alg::RadauIIA9) = 7 -alg_adaptive_order(alg::ImplicitEuler) = 0 -alg_adaptive_order(alg::Trapezoid) = 1 # this is actually incorrect and is purposefully decreased as this tends # to track the real error much better -alg_adaptive_order(alg::ImplicitMidpoint) = 1 # this is actually incorrect and is purposefully decreased as this tends # to track the real error much better @@ -620,7 +568,7 @@ function _digest_beta1_beta2(alg, cache, ::Val{QT}, _beta1, _beta2) where {QT} end # other special cases in controllers.jl -function default_controller(alg::Union{JVODE, QNDF, FBDF}, args...) +function default_controller(alg::Union{JVODE}, args...) DummyController() end @@ -646,17 +594,12 @@ gamma_default(alg::CompositeAlgorithm) = maximum(gamma_default, alg.algs) fac_default_gamma(alg) = false qsteady_min_default(alg::Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}) = 1 -qsteady_min_default(alg::FBDF) = 9 // 10 qsteady_max_default(alg::Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}) = 1 qsteady_max_default(alg::OrdinaryDiffEqAdaptiveImplicitAlgorithm) = 6 // 5 # But don't re-use Jacobian if not adaptive: too risky and cannot pull back qsteady_max_default(alg::OrdinaryDiffEqImplicitAlgorithm) = isadaptive(alg) ? 1 // 1 : 0 qsteady_max_default(alg::AN5) = 3 // 2 qsteady_max_default(alg::JVODE) = 3 // 2 -qsteady_max_default(alg::QNDF1) = 2 // 1 -qsteady_max_default(alg::QNDF2) = 2 // 1 -qsteady_max_default(alg::QNDF) = 2 // 1 -qsteady_max_default(alg::FBDF) = 2 // 1 #TODO #DiffEqBase.nlsolve_default(::QNDF, ::Val{κ}) = 1//2 @@ -671,7 +614,6 @@ ssp_coefficient(alg::Euler) = 1 # We shouldn't do this probably. #ssp_coefficient(alg::ImplicitEuler) = Inf -ssp_coefficient(alg::SSPSDIRK2) = 4 # stability regions alg_stability_size(alg::ExplicitRK) = alg.tableau.stability_size @@ -752,21 +694,9 @@ isWmethod(alg::ROS34PRw) = true isWmethod(alg::ROK4a) = true isWmethod(alg::RosenbrockW6S4OS) = true -isesdirk(alg::TRBDF2) = true -function isesdirk(alg::Union{KenCarp3, KenCarp4, KenCarp5, KenCarp58, - Kvaerno3, Kvaerno4, Kvaerno5, ESDIRK437L2SA, - ESDIRK54I8L2SA, ESDIRK436L2SA2, ESDIRK547L2SA2, - ESDIRK659L2SA, CFNLIRK3}) - true -end isesdirk(alg::Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}) = false is_mass_matrix_alg(alg::Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}) = false is_mass_matrix_alg(alg::CompositeAlgorithm) = all(is_mass_matrix_alg, alg.algs) is_mass_matrix_alg(alg::RosenbrockAlgorithm) = true is_mass_matrix_alg(alg::NewtonAlgorithm) = !isesdirk(alg) -# hack for the default alg -function is_mass_matrix_alg(alg::CompositeAlgorithm{ - <:Any, <:Tuple{Tsit5, Rosenbrock23, Rodas5P, FBDF, FBDF}}) - true -end diff --git a/src/algorithms.jl b/src/algorithms.jl index abeea83bdb..1cc85fe0a1 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -490,340 +490,6 @@ function CNLF2(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Va extrapolant) end -""" -QNDF1: Multistep Method -An adaptive order 1 quasi-constant timestep L-stable numerical differentiation function (NDF) method. -Optional parameter kappa defaults to Shampine's accuracy-optimal -0.1850. - -See also `QNDF`. -""" -struct QNDF1{CS, AD, F, F2, P, FDT, ST, CJ, κType, StepLimiter} <: - OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - extrapolant::Symbol - kappa::κType - controller::Symbol - step_limiter!::StepLimiter -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 = -0.1850, - controller = :Standard, step_limiter! = trivial_limiter!) - QNDF1{ - _unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve), - typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac), - typeof(kappa), typeof(step_limiter!)}(linsolve, - nlsolve, - precs, - extrapolant, - kappa, - controller, - step_limiter!) -end - -""" -QBDF1: Multistep Method - -An alias of `QNDF1` with κ=0. -""" -QBDF1(; kwargs...) = QNDF1(; kappa = 0, kwargs...) - -""" -QNDF2: Multistep Method -An adaptive order 2 quasi-constant timestep L-stable numerical differentiation function (NDF) method. - -See also `QNDF`. -""" -struct QNDF2{CS, AD, F, F2, P, FDT, ST, CJ, κType, StepLimiter} <: - OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - extrapolant::Symbol - kappa::κType - controller::Symbol - step_limiter!::StepLimiter -end - -function QNDF2(; 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 = -1 // 9, - controller = :Standard, step_limiter! = trivial_limiter!) - QNDF2{ - _unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve), - typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac), - typeof(kappa), typeof(step_limiter!)}(linsolve, - nlsolve, - precs, - extrapolant, - kappa, - controller, - step_limiter!) -end - -""" -QBDF2: Multistep Method - -An alias of `QNDF2` with κ=0. -""" -QBDF2(; kwargs...) = QNDF2(; kappa = 0, kwargs...) - -""" -QNDF: Multistep Method -An adaptive order quasi-constant timestep NDF method. -Utilizes Shampine's accuracy-optimal kappa values as defaults (has a keyword argument for a tuple of kappa coefficients). - -@article{shampine1997matlab, -title={The matlab ode suite}, -author={Shampine, Lawrence F and Reichelt, Mark W}, -journal={SIAM journal on scientific computing}, -volume={18}, -number={1}, -pages={1--22}, -year={1997}, -publisher={SIAM} -} -""" -struct QNDF{MO, CS, AD, F, F2, P, FDT, ST, CJ, K, T, κType, StepLimiter} <: - OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} - max_order::Val{MO} - linsolve::F - nlsolve::F2 - precs::P - κ::K - tol::T - extrapolant::Symbol - kappa::κType - controller::Symbol - step_limiter!::StepLimiter -end - -function QNDF(; max_order::Val{MO} = Val{5}(), chunk_size = Val{0}(), - autodiff = Val{true}(), standardtag = Val{true}(), concrete_jac = nothing, - diff_type = Val{:forward}, - linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), κ = nothing, - tol = nothing, - extrapolant = :linear, kappa = promote(-0.1850, -1 // 9, -0.0823, -0.0415, 0), - 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), - _unwrap_val(concrete_jac), - typeof(κ), typeof(tol), typeof(kappa), typeof(step_limiter!)}( - max_order, linsolve, nlsolve, precs, κ, tol, - extrapolant, kappa, controller, step_limiter!) -end - -TruncatedStacktraces.@truncate_stacktrace QNDF - -""" -QBDF: Multistep Method - -An alias of `QNDF` with κ=0. -""" -QBDF(; kwargs...) = QNDF(; kappa = tuple(0 // 1, 0 // 1, 0 // 1, 0 // 1, 0 // 1), kwargs...) - -""" -FBDF: Fixed leading coefficient BDF - -An adaptive order quasi-constant timestep NDF method. -Utilizes Shampine's accuracy-optimal kappa values as defaults (has a keyword argument for a tuple of kappa coefficients). - -@article{shampine2002solving, -title={Solving 0= F (t, y (t), y′(t)) in Matlab}, -author={Shampine, Lawrence F}, -year={2002}, -publisher={Walter de Gruyter GmbH \\& Co. KG} -} -""" -struct FBDF{MO, CS, AD, F, F2, P, FDT, ST, CJ, K, T, StepLimiter} <: - OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} - max_order::Val{MO} - linsolve::F - nlsolve::F2 - precs::P - κ::K - tol::T - extrapolant::Symbol - controller::Symbol - step_limiter!::StepLimiter -end - -function FBDF(; max_order::Val{MO} = Val{5}(), chunk_size = Val{0}(), - autodiff = Val{true}(), standardtag = Val{true}(), concrete_jac = nothing, - diff_type = Val{:forward}, - linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), κ = nothing, - tol = nothing, - extrapolant = :linear, controller = :Standard, step_limiter! = trivial_limiter!) where {MO} - FBDF{MO, _unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac), - typeof(κ), typeof(tol), typeof(step_limiter!)}( - max_order, linsolve, nlsolve, precs, κ, tol, extrapolant, - controller, step_limiter!) -end - -TruncatedStacktraces.@truncate_stacktrace FBDF - -""" -Uri M. Ascher, Steven J. Ruuth, Brian T. R. Wetton. Implicit-Explicit Methods for Time- -Dependent Partial Differential Equations. 1995 Society for Industrial and Applied Mathematics -Journal on Numerical Analysis, 32(3), pp 797-823, 1995. doi: https://doi.org/10.1137/0732037 -""" -struct SBDF{CS, AD, F, F2, P, FDT, ST, CJ, K, T} <: - OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - κ::K - tol::T - extrapolant::Symbol - order::Int - ark::Bool -end - -function SBDF(order; chunk_size = Val{0}(), autodiff = Val{true}(), - standardtag = Val{true}(), concrete_jac = nothing, diff_type = Val{:forward}, - linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), κ = nothing, - tol = nothing, - extrapolant = :linear, ark = false) - SBDF{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve), - typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac), - typeof(κ), typeof(tol)}(linsolve, - nlsolve, - precs, - κ, - tol, - extrapolant, - order, - ark) -end - -# All keyword form needed for remake -function SBDF(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(), - concrete_jac = nothing, diff_type = Val{:forward}, - linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), κ = nothing, - tol = nothing, - extrapolant = :linear, - order, ark = false) - SBDF{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve), - typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac), - typeof(κ), typeof(tol)}(linsolve, - nlsolve, - precs, - κ, - tol, - extrapolant, - order, - ark) -end - -""" - IMEXEuler(;kwargs...) - -The one-step version of the IMEX multistep methods of - - - Uri M. Ascher, Steven J. Ruuth, Brian T. R. Wetton. - Implicit-Explicit Methods for Time-Dependent Partial Differential Equations. - Society for Industrial and Applied Mathematics. - Journal on Numerical Analysis, 32(3), pp 797-823, 1995. - doi: [https://doi.org/10.1137/0732037](https://doi.org/10.1137/0732037) - -When applied to a `SplitODEProblem` of the form - -``` -u'(t) = f1(u) + f2(u) -``` - -The default `IMEXEuler()` method uses an update of the form - -``` -unew = uold + dt * (f1(unew) + f2(uold)) -``` - -See also `SBDF`, `IMEXEulerARK`. -""" -IMEXEuler(; kwargs...) = SBDF(1; kwargs...) - -""" - IMEXEulerARK(;kwargs...) - -The one-step version of the IMEX multistep methods of - - - Uri M. Ascher, Steven J. Ruuth, Brian T. R. Wetton. - Implicit-Explicit Methods for Time-Dependent Partial Differential Equations. - Society for Industrial and Applied Mathematics. - Journal on Numerical Analysis, 32(3), pp 797-823, 1995. - doi: [https://doi.org/10.1137/0732037](https://doi.org/10.1137/0732037) - -When applied to a `SplitODEProblem` of the form - -``` -u'(t) = f1(u) + f2(u) -``` - -A classical additive Runge-Kutta method in the sense of -[Araújo, Murua, Sanz-Serna (1997)](https://doi.org/10.1137/S0036142995292128) -consisting of the implicit and the explicit Euler method given by - -``` -y1 = uold + dt * f1(y1) -unew = uold + dt * (f1(unew) + f2(y1)) -``` - -See also `SBDF`, `IMEXEuler`. -""" -IMEXEulerARK(; kwargs...) = SBDF(1; ark = true, kwargs...) - -""" - SBDF2(;kwargs...) - -The two-step version of the IMEX multistep methods of - - - Uri M. Ascher, Steven J. Ruuth, Brian T. R. Wetton. - Implicit-Explicit Methods for Time-Dependent Partial Differential Equations. - Society for Industrial and Applied Mathematics. - Journal on Numerical Analysis, 32(3), pp 797-823, 1995. - doi: [https://doi.org/10.1137/0732037](https://doi.org/10.1137/0732037) - -See also `SBDF`. -""" -SBDF2(; kwargs...) = SBDF(2; kwargs...) - -""" - SBDF3(;kwargs...) - -The three-step version of the IMEX multistep methods of - - - Uri M. Ascher, Steven J. Ruuth, Brian T. R. Wetton. - Implicit-Explicit Methods for Time-Dependent Partial Differential Equations. - Society for Industrial and Applied Mathematics. - Journal on Numerical Analysis, 32(3), pp 797-823, 1995. - doi: [https://doi.org/10.1137/0732037](https://doi.org/10.1137/0732037) - -See also `SBDF`. -""" -SBDF3(; kwargs...) = SBDF(3; kwargs...) - -""" - SBDF4(;kwargs...) - -The four-step version of the IMEX multistep methods of - - - Uri M. Ascher, Steven J. Ruuth, Brian T. R. Wetton. - Implicit-Explicit Methods for Time-Dependent Partial Differential Equations. - Society for Industrial and Applied Mathematics. - Journal on Numerical Analysis, 32(3), pp 797-823, 1995. - doi: [https://doi.org/10.1137/0732037](https://doi.org/10.1137/0732037) - -See also `SBDF`. -""" -SBDF4(; kwargs...) = SBDF(4; kwargs...) - # Adams/BDF methods in Nordsieck forms """ AN5: Adaptive step size Adams explicit Method @@ -871,1047 +537,184 @@ for Alg in [ m::Int iop::Int end - @eval $Alg(; krylov = false, m = 30, iop = 0) = $Alg(krylov, m, iop) -end - -struct MagnusAdapt4 <: OrdinaryDiffEqAdaptiveAlgorithm end - -struct LinearExponential <: - OrdinaryDiffEqExponentialAlgorithm{1, false, Val{:forward}, Val{true}, nothing} - krylov::Symbol - m::Int - iop::Int -end -LinearExponential(; krylov = :off, m = 10, iop = 0) = LinearExponential(krylov, m, iop) - -struct CayleyEuler <: OrdinaryDiffEqAlgorithm end - -################################################################################ - -# FIRK Methods - -""" -@article{hairer1999stiff, -title={Stiff differential equations solved by Radau methods}, -author={Hairer, Ernst and Wanner, Gerhard}, -journal={Journal of Computational and Applied Mathematics}, -volume={111}, -number={1-2}, -pages={93--111}, -year={1999}, -publisher={Elsevier} -} - -RadauIIA3: Fully-Implicit Runge-Kutta Method -An A-B-L stable fully implicit Runge-Kutta method with internal tableau complex basis transform for efficiency. -""" -struct RadauIIA3{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <: - OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - precs::P - extrapolant::Symbol - κ::Tol - maxiters::Int - fast_convergence_cutoff::C1 - new_W_γdt_cutoff::C2 - controller::Symbol - step_limiter!::StepLimiter -end - -function RadauIIA3(; chunk_size = Val{0}(), autodiff = Val{true}(), - standardtag = Val{true}(), concrete_jac = nothing, - diff_type = Val{:forward}, - linsolve = nothing, precs = DEFAULT_PRECS, - extrapolant = :dense, fast_convergence_cutoff = 1 // 5, - new_W_γdt_cutoff = 1 // 5, - controller = :Predictive, κ = nothing, maxiters = 10, - step_limiter! = trivial_limiter!) - RadauIIA3{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac), - typeof(κ), typeof(fast_convergence_cutoff), - typeof(new_W_γdt_cutoff), typeof(step_limiter!)}(linsolve, - precs, - extrapolant, - κ, - maxiters, - fast_convergence_cutoff, - new_W_γdt_cutoff, - controller, - step_limiter!) -end - -TruncatedStacktraces.@truncate_stacktrace RadauIIA3 - -""" -@article{hairer1999stiff, -title={Stiff differential equations solved by Radau methods}, -author={Hairer, Ernst and Wanner, Gerhard}, -journal={Journal of Computational and Applied Mathematics}, -volume={111}, -number={1-2}, -pages={93--111}, -year={1999}, -publisher={Elsevier} -} - -RadauIIA5: Fully-Implicit Runge-Kutta Method -An A-B-L stable fully implicit Runge-Kutta method with internal tableau complex basis transform for efficiency. -""" -struct RadauIIA5{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <: - OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - precs::P - smooth_est::Bool - extrapolant::Symbol - κ::Tol - maxiters::Int - fast_convergence_cutoff::C1 - new_W_γdt_cutoff::C2 - controller::Symbol - step_limiter!::StepLimiter -end - -function RadauIIA5(; chunk_size = Val{0}(), autodiff = Val{true}(), - standardtag = Val{true}(), concrete_jac = nothing, - diff_type = Val{:forward}, - linsolve = nothing, precs = DEFAULT_PRECS, - extrapolant = :dense, fast_convergence_cutoff = 1 // 5, - new_W_γdt_cutoff = 1 // 5, - controller = :Predictive, κ = nothing, maxiters = 10, smooth_est = true, - step_limiter! = trivial_limiter!) - RadauIIA5{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac), - typeof(κ), typeof(fast_convergence_cutoff), - typeof(new_W_γdt_cutoff), typeof(step_limiter!)}(linsolve, - precs, - smooth_est, - extrapolant, - κ, - maxiters, - fast_convergence_cutoff, - new_W_γdt_cutoff, - controller, - step_limiter!) -end -TruncatedStacktraces.@truncate_stacktrace RadauIIA5 - -""" -@article{hairer1999stiff, -title={Stiff differential equations solved by Radau methods}, -author={Hairer, Ernst and Wanner, Gerhard}, -journal={Journal of Computational and Applied Mathematics}, -volume={111}, -number={1-2}, -pages={93--111}, -year={1999}, -publisher={Elsevier} -} - -RadauII97: Fully-Implicit Runge-Kutta Method -An A-B-L stable fully implicit Runge-Kutta method with internal tableau complex basis transform for efficiency. -""" -struct RadauIIA9{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <: - OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - precs::P - smooth_est::Bool - extrapolant::Symbol - κ::Tol - maxiters::Int - fast_convergence_cutoff::C1 - new_W_γdt_cutoff::C2 - controller::Symbol - step_limiter!::StepLimiter -end - -function RadauIIA9(; chunk_size = Val{0}(), autodiff = Val{true}(), - standardtag = Val{true}(), concrete_jac = nothing, - diff_type = Val{:forward}, - linsolve = nothing, precs = DEFAULT_PRECS, - extrapolant = :dense, fast_convergence_cutoff = 1 // 5, - new_W_γdt_cutoff = 1 // 5, - controller = :Predictive, κ = nothing, maxiters = 10, smooth_est = true, - step_limiter! = trivial_limiter!) - RadauIIA9{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac), - typeof(κ), typeof(fast_convergence_cutoff), - typeof(new_W_γdt_cutoff), typeof(step_limiter!)}(linsolve, - precs, - smooth_est, - extrapolant, - κ, - maxiters, - fast_convergence_cutoff, - new_W_γdt_cutoff, - controller, - step_limiter!) -end -TruncatedStacktraces.@truncate_stacktrace RadauIIA9 - -################################################################################ - -# SDIRK Methods -""" -ImplicitEuler: SDIRK Method -A 1st order implicit solver. A-B-L-stable. Adaptive timestepping through a divided differences estimate via memory. -Strong-stability preserving (SSP). -""" -struct ImplicitEuler{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <: - OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - extrapolant::Symbol - controller::Symbol - step_limiter!::StepLimiter -end - -function ImplicitEuler(; 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 = :constant, - controller = :PI, step_limiter! = trivial_limiter!) - ImplicitEuler{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, - nlsolve, precs, extrapolant, controller, step_limiter!) -end -""" -ImplicitMidpoint: SDIRK Method -A second order A-stable symplectic and symmetric implicit solver. -Good for highly stiff equations which need symplectic integration. -""" -struct ImplicitMidpoint{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <: - OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - extrapolant::Symbol - step_limiter!::StepLimiter -end - -function ImplicitMidpoint(; 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, step_limiter! = trivial_limiter!) - ImplicitMidpoint{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, - nlsolve, - precs, - extrapolant, - step_limiter!) -end - -""" -Andre Vladimirescu. 1994. The Spice Book. John Wiley & Sons, Inc., New York, -NY, USA. - -Trapezoid: SDIRK Method -A second order A-stable symmetric ESDIRK method. -"Almost symplectic" without numerical dampening. -Also known as Crank-Nicolson when applied to PDEs. Adaptive timestepping via divided -differences approximation to the second derivative terms in the local truncation error -estimate (the SPICE approximation strategy). -""" -struct Trapezoid{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <: - OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - extrapolant::Symbol - controller::Symbol - step_limiter!::StepLimiter -end - -function Trapezoid(; 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, - controller = :PI, step_limiter! = trivial_limiter!) - Trapezoid{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, - nlsolve, - precs, - extrapolant, - controller, - step_limiter!) -end - -""" -@article{hosea1996analysis, -title={Analysis and implementation of TR-BDF2}, -author={Hosea, ME and Shampine, LF}, -journal={Applied Numerical Mathematics}, -volume={20}, -number={1-2}, -pages={21--37}, -year={1996}, -publisher={Elsevier} -} - -TRBDF2: SDIRK Method -A second order A-B-L-S-stable one-step ESDIRK method. -Includes stiffness-robust error estimates for accurate adaptive timestepping, smoothed derivatives for highly stiff and oscillatory problems. -""" -struct TRBDF2{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <: - OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - smooth_est::Bool - extrapolant::Symbol - controller::Symbol - step_limiter!::StepLimiter -end - -function TRBDF2(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(), - concrete_jac = nothing, diff_type = Val{:forward}, - linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), - smooth_est = true, extrapolant = :linear, - controller = :PI, step_limiter! = trivial_limiter!) - TRBDF2{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, nlsolve, precs, - smooth_est, extrapolant, controller, step_limiter!) -end - -TruncatedStacktraces.@truncate_stacktrace TRBDF2 - -""" -@article{hindmarsh2005sundials, -title={{SUNDIALS}: Suite of nonlinear and differential/algebraic equation solvers}, -author={Hindmarsh, Alan C and Brown, Peter N and Grant, Keith E and Lee, Steven L and Serban, Radu and Shumaker, Dan E and Woodward, Carol S}, -journal={ACM Transactions on Mathematical Software (TOMS)}, -volume={31}, -number={3}, -pages={363--396}, -year={2005}, -publisher={ACM} -} - -SDIRK2: SDIRK Method -An A-B-L stable 2nd order SDIRK method -""" -struct SDIRK2{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <: - OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - smooth_est::Bool - extrapolant::Symbol - controller::Symbol - step_limiter!::StepLimiter -end - -function SDIRK2(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(), - concrete_jac = nothing, diff_type = Val{:forward}, - linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), - smooth_est = true, extrapolant = :linear, - controller = :PI, step_limiter! = trivial_limiter!) - SDIRK2{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac), typeof(step_limiter!)}( - linsolve, nlsolve, precs, smooth_est, extrapolant, - controller, - step_limiter!) -end - -struct SDIRK22{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <: - OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - extrapolant::Symbol - controller::Symbol - step_limiter!::StepLimiter -end - -function SDIRK22(; - 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, - controller = :PI, step_limiter! = trivial_limiter!) - Trapezoid{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, - nlsolve, - precs, - extrapolant, - controller, - step_limiter!) -end - -struct SSPSDIRK2{CS, AD, F, F2, P, FDT, ST, CJ} <: - OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ} # Not adaptive - linsolve::F - nlsolve::F2 - precs::P - smooth_est::Bool - extrapolant::Symbol - controller::Symbol -end - -function SSPSDIRK2(; chunk_size = Val{0}(), autodiff = Val{true}(), - standardtag = Val{true}(), concrete_jac = nothing, - diff_type = Val{:forward}, - linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), - smooth_est = true, extrapolant = :constant, - controller = :PI) - SSPSDIRK2{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac)}(linsolve, nlsolve, precs, smooth_est, extrapolant, - controller) -end - -""" -@article{kvaerno2004singly, -title={Singly diagonally implicit Runge--Kutta methods with an explicit first stage}, -author={Kv{\\ae}rn{\\o}, Anne}, -journal={BIT Numerical Mathematics}, -volume={44}, -number={3}, -pages={489--502}, -year={2004}, -publisher={Springer} -} - -Kvaerno3: SDIRK Method -An A-L stable stiffly-accurate 3rd order ESDIRK method -""" -struct Kvaerno3{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <: - OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - smooth_est::Bool - extrapolant::Symbol - controller::Symbol - step_limiter!::StepLimiter -end -function Kvaerno3(; chunk_size = Val{0}(), autodiff = Val{true}(), - standardtag = Val{true}(), concrete_jac = nothing, - diff_type = Val{:forward}, - linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), - smooth_est = true, extrapolant = :linear, - controller = :PI, step_limiter! = trivial_limiter!) - Kvaerno3{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, nlsolve, precs, - smooth_est, extrapolant, controller, step_limiter!) -end - -""" -@book{kennedy2001additive, -title={Additive Runge-Kutta schemes for convection-diffusion-reaction equations}, -author={Kennedy, Christopher Alan}, -year={2001}, -publisher={National Aeronautics and Space Administration, Langley Research Center} -} - -KenCarp3: SDIRK Method -An A-L stable stiffly-accurate 3rd order ESDIRK method with splitting -""" -struct KenCarp3{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <: - OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - smooth_est::Bool - extrapolant::Symbol - controller::Symbol - step_limiter!::StepLimiter -end -function KenCarp3(; chunk_size = Val{0}(), autodiff = Val{true}(), - standardtag = Val{true}(), concrete_jac = nothing, - diff_type = Val{:forward}, - linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), - smooth_est = true, extrapolant = :linear, - controller = :PI, step_limiter! = trivial_limiter!) - KenCarp3{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, nlsolve, precs, - smooth_est, extrapolant, controller, step_limiter!) -end - -struct CFNLIRK3{CS, AD, F, F2, P, FDT, ST, CJ} <: - OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - extrapolant::Symbol -end -function CFNLIRK3(; 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) - CFNLIRK3{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac)}(linsolve, - nlsolve, - precs, - extrapolant) -end - -""" -@article{hindmarsh2005sundials, -title={{SUNDIALS}: Suite of nonlinear and differential/algebraic equation solvers}, -author={Hindmarsh, Alan C and Brown, Peter N and Grant, Keith E and Lee, Steven L and Serban, Radu and Shumaker, Dan E and Woodward, Carol S}, -journal={ACM Transactions on Mathematical Software (TOMS)}, -volume={31}, -number={3}, -pages={363--396}, -year={2005}, -publisher={ACM} -} - -Cash4: SDIRK Method -An A-L stable 4th order SDIRK method -""" -struct Cash4{CS, AD, F, F2, P, FDT, ST, CJ} <: - OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - smooth_est::Bool - extrapolant::Symbol - embedding::Int - controller::Symbol -end -function Cash4(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(), - concrete_jac = nothing, diff_type = Val{:forward}, - linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), - smooth_est = true, extrapolant = :linear, - controller = :PI, embedding = 3) - Cash4{ - _unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve), - typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac)}( - linsolve, - nlsolve, - precs, - smooth_est, - extrapolant, - embedding, - controller) -end - -struct SFSDIRK4{CS, AD, F, F2, P, FDT, ST, CJ} <: - OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - extrapolant::Symbol -end -function SFSDIRK4(; 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) - SFSDIRK4{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac)}(linsolve, - nlsolve, - precs, - extrapolant) -end - -struct SFSDIRK5{CS, AD, F, F2, P, FDT, ST, CJ} <: - OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - extrapolant::Symbol -end - -function SFSDIRK5(; 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) - SFSDIRK5{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac)}(linsolve, - nlsolve, - precs, - extrapolant) -end - -struct SFSDIRK6{CS, AD, F, F2, P, FDT, ST, CJ} <: - OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - extrapolant::Symbol -end - -function SFSDIRK6(; 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) - SFSDIRK6{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac)}(linsolve, - nlsolve, - precs, - extrapolant) -end - -struct SFSDIRK7{CS, AD, F, F2, P, FDT, ST, CJ} <: - OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - extrapolant::Symbol -end - -function SFSDIRK7(; 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) - SFSDIRK7{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac)}(linsolve, - nlsolve, - precs, - extrapolant) -end - -struct SFSDIRK8{CS, AD, F, F2, P, FDT, ST, CJ} <: - OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - extrapolant::Symbol -end - -function SFSDIRK8(; 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) - SFSDIRK8{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac)}(linsolve, - nlsolve, - precs, - extrapolant) -end - -""" -E. Hairer, G. Wanner, Solving ordinary differential equations II, stiff and -differential-algebraic problems. Computational mathematics (2nd revised ed.), -Springer (1996) - -Hairer4: SDIRK Method -An A-L stable 4th order SDIRK method -""" -struct Hairer4{CS, AD, F, F2, P, FDT, ST, CJ} <: - OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - smooth_est::Bool - extrapolant::Symbol - controller::Symbol -end -function Hairer4(; - chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(), - concrete_jac = nothing, diff_type = Val{:forward}, - linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), - smooth_est = true, extrapolant = :linear, - controller = :PI) - Hairer4{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac)}(linsolve, nlsolve, precs, smooth_est, extrapolant, - controller) -end - -""" -E. Hairer, G. Wanner, Solving ordinary differential equations II, stiff and -differential-algebraic problems. Computational mathematics (2nd revised ed.), -Springer (1996) - -Hairer42: SDIRK Method -An A-L stable 4th order SDIRK method -""" -struct Hairer42{CS, AD, F, F2, P, FDT, ST, CJ} <: - OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - smooth_est::Bool - extrapolant::Symbol - controller::Symbol -end -function Hairer42(; chunk_size = Val{0}(), autodiff = Val{true}(), - standardtag = Val{true}(), concrete_jac = nothing, - diff_type = Val{:forward}, - linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), - smooth_est = true, extrapolant = :linear, - controller = :PI) - Hairer42{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac)}(linsolve, nlsolve, precs, smooth_est, extrapolant, - controller) -end - -""" -@article{kvaerno2004singly, -title={Singly diagonally implicit Runge--Kutta methods with an explicit first stage}, -author={Kv{\\ae}rn{\\o}, Anne}, -journal={BIT Numerical Mathematics}, -volume={44}, -number={3}, -pages={489--502}, -year={2004}, -publisher={Springer} -} - -Kvaerno4: SDIRK Method -An A-L stable stiffly-accurate 4th order ESDIRK method. -""" -struct Kvaerno4{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <: - OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - smooth_est::Bool - extrapolant::Symbol - controller::Symbol - step_limiter!::StepLimiter -end -function Kvaerno4(; chunk_size = Val{0}(), autodiff = Val{true}(), - standardtag = Val{true}(), concrete_jac = nothing, - diff_type = Val{:forward}, - linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), - smooth_est = true, extrapolant = :linear, - controller = :PI, step_limiter! = trivial_limiter!) - Kvaerno4{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, nlsolve, precs, - smooth_est, extrapolant, controller, step_limiter!) -end - -""" -@article{kvaerno2004singly, -title={Singly diagonally implicit Runge--Kutta methods with an explicit first stage}, -author={Kv{\\ae}rn{\\o}, Anne}, -journal={BIT Numerical Mathematics}, -volume={44}, -number={3}, -pages={489--502}, -year={2004}, -publisher={Springer} -} - -Kvaerno5: SDIRK Method -An A-L stable stiffly-accurate 5th order ESDIRK method -""" -struct Kvaerno5{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <: - OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - smooth_est::Bool - extrapolant::Symbol - controller::Symbol - step_limiter!::StepLimiter -end -function Kvaerno5(; chunk_size = Val{0}(), autodiff = Val{true}(), - standardtag = Val{true}(), concrete_jac = nothing, - diff_type = Val{:forward}, - linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), - smooth_est = true, extrapolant = :linear, - controller = :PI, step_limiter! = trivial_limiter!) - Kvaerno5{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, nlsolve, precs, - smooth_est, extrapolant, controller, step_limiter!) -end - -""" -@book{kennedy2001additive, -title={Additive Runge-Kutta schemes for convection-diffusion-reaction equations}, -author={Kennedy, Christopher Alan}, -year={2001}, -publisher={National Aeronautics and Space Administration, Langley Research Center} -} - -KenCarp4: SDIRK Method -An A-L stable stiffly-accurate 4th order ESDIRK method with splitting -""" -struct KenCarp4{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <: - OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - smooth_est::Bool - extrapolant::Symbol - controller::Symbol - step_limiter!::StepLimiter + @eval $Alg(; krylov = false, m = 30, iop = 0) = $Alg(krylov, m, iop) end -function KenCarp4(; chunk_size = Val{0}(), autodiff = Val{true}(), - standardtag = Val{true}(), concrete_jac = nothing, - diff_type = Val{:forward}, - linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), - smooth_est = true, extrapolant = :linear, - controller = :PI, step_limiter! = trivial_limiter!) - KenCarp4{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, nlsolve, precs, - smooth_est, extrapolant, controller, step_limiter!) + +struct MagnusAdapt4 <: OrdinaryDiffEqAdaptiveAlgorithm end + +struct LinearExponential <: + OrdinaryDiffEqExponentialAlgorithm{1, false, Val{:forward}, Val{true}, nothing} + krylov::Symbol + m::Int + iop::Int end +LinearExponential(; krylov = :off, m = 10, iop = 0) = LinearExponential(krylov, m, iop) -TruncatedStacktraces.@truncate_stacktrace KenCarp4 +struct CayleyEuler <: OrdinaryDiffEqAlgorithm end -""" -@article{kennedy2019higher, -title={Higher-order additive Runge--Kutta schemes for ordinary differential equations}, -author={Kennedy, Christopher A and Carpenter, Mark H}, -journal={Applied Numerical Mathematics}, -volume={136}, -pages={183--205}, -year={2019}, -publisher={Elsevier} -} +################################################################################ -KenCarp47: SDIRK Method -An A-L stable stiffly-accurate 4th order seven-stage ESDIRK method with splitting -""" -struct KenCarp47{CS, AD, F, F2, P, FDT, ST, CJ} <: - OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - smooth_est::Bool - extrapolant::Symbol - controller::Symbol -end -function KenCarp47(; chunk_size = Val{0}(), autodiff = Val{true}(), - standardtag = Val{true}(), concrete_jac = nothing, - diff_type = Val{:forward}, - linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), - smooth_est = true, extrapolant = :linear, - controller = :PI) - KenCarp47{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac)}(linsolve, nlsolve, precs, smooth_est, extrapolant, - controller) -end +# FIRK Methods """ -@book{kennedy2001additive, -title={Additive Runge-Kutta schemes for convection-diffusion-reaction equations}, -author={Kennedy, Christopher Alan}, -year={2001}, -publisher={National Aeronautics and Space Administration, Langley Research Center} +@article{hairer1999stiff, +title={Stiff differential equations solved by Radau methods}, +author={Hairer, Ernst and Wanner, Gerhard}, +journal={Journal of Computational and Applied Mathematics}, +volume={111}, +number={1-2}, +pages={93--111}, +year={1999}, +publisher={Elsevier} } -KenCarp5: SDIRK Method -An A-L stable stiffly-accurate 5th order ESDIRK method with splitting +RadauIIA3: Fully-Implicit Runge-Kutta Method +An A-B-L stable fully implicit Runge-Kutta method with internal tableau complex basis transform for efficiency. """ -struct KenCarp5{CS, AD, F, F2, P, FDT, ST, CJ, StepLimiter} <: +struct RadauIIA3{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} linsolve::F - nlsolve::F2 precs::P - smooth_est::Bool extrapolant::Symbol + κ::Tol + maxiters::Int + fast_convergence_cutoff::C1 + new_W_γdt_cutoff::C2 controller::Symbol step_limiter!::StepLimiter end -function KenCarp5(; chunk_size = Val{0}(), autodiff = Val{true}(), + +function RadauIIA3(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(), concrete_jac = nothing, diff_type = Val{:forward}, - linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), - smooth_est = true, extrapolant = :linear, - controller = :PI, step_limiter! = trivial_limiter!) - KenCarp5{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac), typeof(step_limiter!)}(linsolve, nlsolve, precs, - smooth_est, extrapolant, controller, step_limiter!) + linsolve = nothing, precs = DEFAULT_PRECS, + extrapolant = :dense, fast_convergence_cutoff = 1 // 5, + new_W_γdt_cutoff = 1 // 5, + controller = :Predictive, κ = nothing, maxiters = 10, + step_limiter! = trivial_limiter!) + RadauIIA3{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac), + typeof(κ), typeof(fast_convergence_cutoff), + typeof(new_W_γdt_cutoff), typeof(step_limiter!)}(linsolve, + precs, + extrapolant, + κ, + maxiters, + fast_convergence_cutoff, + new_W_γdt_cutoff, + controller, + step_limiter!) end + +TruncatedStacktraces.@truncate_stacktrace RadauIIA3 + """ -@article{kennedy2019higher, -title={Higher-order additive Runge--Kutta schemes for ordinary differential equations}, -author={Kennedy, Christopher A and Carpenter, Mark H}, -journal={Applied Numerical Mathematics}, -volume={136}, -pages={183--205}, -year={2019}, +@article{hairer1999stiff, +title={Stiff differential equations solved by Radau methods}, +author={Hairer, Ernst and Wanner, Gerhard}, +journal={Journal of Computational and Applied Mathematics}, +volume={111}, +number={1-2}, +pages={93--111}, +year={1999}, publisher={Elsevier} } -KenCarp58: SDIRK Method -An A-L stable stiffly-accurate 5th order eight-stage ESDIRK method with splitting +RadauIIA5: Fully-Implicit Runge-Kutta Method +An A-B-L stable fully implicit Runge-Kutta method with internal tableau complex basis transform for efficiency. """ -struct KenCarp58{CS, AD, F, F2, P, FDT, ST, CJ} <: +struct RadauIIA5{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} linsolve::F - nlsolve::F2 precs::P smooth_est::Bool extrapolant::Symbol + κ::Tol + maxiters::Int + fast_convergence_cutoff::C1 + new_W_γdt_cutoff::C2 controller::Symbol -end -function KenCarp58(; chunk_size = Val{0}(), autodiff = Val{true}(), - standardtag = Val{true}(), concrete_jac = nothing, - diff_type = Val{:forward}, - linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), - smooth_est = true, extrapolant = :linear, - controller = :PI) - KenCarp58{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac)}(linsolve, nlsolve, precs, smooth_est, extrapolant, - controller) + step_limiter!::StepLimiter end -# `smooth_est` is not necessary, as the embedded method is also L-stable -struct ESDIRK54I8L2SA{CS, AD, F, F2, P, FDT, ST, CJ} <: - OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - extrapolant::Symbol - controller::Symbol -end -function ESDIRK54I8L2SA(; chunk_size = Val{0}(), autodiff = Val{true}(), +function RadauIIA5(; 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, controller = :PI) - ESDIRK54I8L2SA{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac)}(linsolve, nlsolve, precs, extrapolant, - controller) + linsolve = nothing, precs = DEFAULT_PRECS, + extrapolant = :dense, fast_convergence_cutoff = 1 // 5, + new_W_γdt_cutoff = 1 // 5, + controller = :Predictive, κ = nothing, maxiters = 10, smooth_est = true, + step_limiter! = trivial_limiter!) + RadauIIA5{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac), + typeof(κ), typeof(fast_convergence_cutoff), + typeof(new_W_γdt_cutoff), typeof(step_limiter!)}(linsolve, + precs, + smooth_est, + extrapolant, + κ, + maxiters, + fast_convergence_cutoff, + new_W_γdt_cutoff, + controller, + step_limiter!) end +TruncatedStacktraces.@truncate_stacktrace RadauIIA5 """ -@article{Kennedy2019DiagonallyIR, -title={Diagonally implicit Runge–Kutta methods for stiff ODEs}, -author={Christopher A. Kennedy and Mark H. Carpenter}, -journal={Applied Numerical Mathematics}, -year={2019}, -volume={146}, -pages={221-244} +@article{hairer1999stiff, +title={Stiff differential equations solved by Radau methods}, +author={Hairer, Ernst and Wanner, Gerhard}, +journal={Journal of Computational and Applied Mathematics}, +volume={111}, +number={1-2}, +pages={93--111}, +year={1999}, +publisher={Elsevier} } -""" -struct ESDIRK436L2SA2{CS, AD, F, F2, P, FDT, ST, CJ} <: - OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - extrapolant::Symbol - controller::Symbol -end -function ESDIRK436L2SA2(; 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, controller = :PI) - ESDIRK436L2SA2{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac)}(linsolve, nlsolve, precs, extrapolant, - controller) -end +RadauII97: Fully-Implicit Runge-Kutta Method +An A-B-L stable fully implicit Runge-Kutta method with internal tableau complex basis transform for efficiency. """ -@article{Kennedy2019DiagonallyIR, -title={Diagonally implicit Runge–Kutta methods for stiff ODEs}, -author={Christopher A. Kennedy and Mark H. Carpenter}, -journal={Applied Numerical Mathematics}, -year={2019}, -volume={146}, -pages={221-244} -} -""" -struct ESDIRK437L2SA{CS, AD, F, F2, P, FDT, ST, CJ} <: +struct RadauIIA9{CS, AD, F, P, FDT, ST, CJ, Tol, C1, C2, StepLimiter} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} linsolve::F - nlsolve::F2 precs::P + smooth_est::Bool extrapolant::Symbol + κ::Tol + maxiters::Int + fast_convergence_cutoff::C1 + new_W_γdt_cutoff::C2 controller::Symbol -end -function ESDIRK437L2SA(; 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, controller = :PI) - ESDIRK437L2SA{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac)}(linsolve, nlsolve, precs, extrapolant, - controller) + step_limiter!::StepLimiter end -""" -@article{Kennedy2019DiagonallyIR, -title={Diagonally implicit Runge–Kutta methods for stiff ODEs}, -author={Christopher A. Kennedy and Mark H. Carpenter}, -journal={Applied Numerical Mathematics}, -year={2019}, -volume={146}, -pages={221-244} -} -""" -struct ESDIRK547L2SA2{CS, AD, F, F2, P, FDT, ST, CJ} <: - OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - extrapolant::Symbol - controller::Symbol -end -function ESDIRK547L2SA2(; chunk_size = Val{0}(), autodiff = Val{true}(), +function RadauIIA9(; 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, controller = :PI) - ESDIRK547L2SA2{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac)}(linsolve, nlsolve, precs, extrapolant, - controller) + linsolve = nothing, precs = DEFAULT_PRECS, + extrapolant = :dense, fast_convergence_cutoff = 1 // 5, + new_W_γdt_cutoff = 1 // 5, + controller = :Predictive, κ = nothing, maxiters = 10, smooth_est = true, + step_limiter! = trivial_limiter!) + RadauIIA9{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), + typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac), + typeof(κ), typeof(fast_convergence_cutoff), + typeof(new_W_γdt_cutoff), typeof(step_limiter!)}(linsolve, + precs, + smooth_est, + extrapolant, + κ, + maxiters, + fast_convergence_cutoff, + new_W_γdt_cutoff, + controller, + step_limiter!) end +TruncatedStacktraces.@truncate_stacktrace RadauIIA9 -""" -@article{Kennedy2019DiagonallyIR, -title={Diagonally implicit Runge–Kutta methods for stiff ODEs}, -author={Christopher A. Kennedy and Mark H. Carpenter}, -journal={Applied Numerical Mathematics}, -year={2019}, -volume={146}, -pages={221-244} - -Currently has STABILITY ISSUES, causing it to fail the adaptive tests. -Check issue https://github.com/SciML/OrdinaryDiffEq.jl/issues/1933 for more details. -} -""" -struct ESDIRK659L2SA{CS, AD, F, F2, P, FDT, ST, CJ} <: - OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - extrapolant::Symbol - controller::Symbol -end -function ESDIRK659L2SA(; 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, controller = :PI) - ESDIRK659L2SA{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac)}(linsolve, nlsolve, precs, extrapolant, - controller) -end +################################################################################ ################################################################################ @@ -2188,39 +991,6 @@ struct ETD2 <: ######################################### -""" -E. Alberdi Celayaa, J. J. Anza Aguirrezabalab, P. Chatzipantelidisc. Implementation of -an Adaptive BDF2 Formula and Comparison with The MATLAB Ode15s. Procedia Computer Science, -29, pp 1014-1026, 2014. doi: https://doi.org/10.1016/j.procs.2014.05.091 - -ABDF2: Multistep Method -An adaptive order 2 L-stable fixed leading coefficient multistep BDF method. -""" -struct ABDF2{CS, AD, F, F2, P, FDT, ST, CJ, K, T, StepLimiter} <: - OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - κ::K - tol::T - smooth_est::Bool - extrapolant::Symbol - controller::Symbol - step_limiter!::StepLimiter -end -function ABDF2(; chunk_size = Val{0}(), autodiff = true, standardtag = Val{true}(), - concrete_jac = nothing, diff_type = Val{:forward}, - κ = nothing, tol = nothing, linsolve = nothing, precs = DEFAULT_PRECS, - nlsolve = NLNewton(), - smooth_est = true, extrapolant = :linear, - controller = :Standard, step_limiter! = trivial_limiter!) - ABDF2{ - _unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve), - typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac), - typeof(κ), typeof(tol), typeof(step_limiter!)}(linsolve, nlsolve, precs, κ, tol, - smooth_est, extrapolant, controller, step_limiter!) -end - ######################################### struct CompositeAlgorithm{CS, T, F} <: OrdinaryDiffEqCompositeAlgorithm @@ -2294,29 +1064,6 @@ struct AutoSwitch{nAlg, sAlg, tolType, T} end ################################################################################ -""" -MEBDF2: Multistep Method -The second order Modified Extended BDF method, which has improved stability properties over the standard BDF. -Fixed timestep only. -""" -struct MEBDF2{CS, AD, F, F2, P, FDT, ST, CJ} <: - OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - extrapolant::Symbol -end -function MEBDF2(; chunk_size = Val{0}(), autodiff = true, standardtag = Val{true}(), - concrete_jac = nothing, diff_type = Val{:forward}, - linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), - extrapolant = :constant) - MEBDF2{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), - typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag), - _unwrap_val(concrete_jac)}(linsolve, - nlsolve, - precs, - extrapolant) -end ################################################# """ @@ -2342,12 +1089,11 @@ function PDIRK44(; chunk_size = Val{0}(), autodiff = true, standardtag = Val{tru end ### Algorithm Groups +# ABDF2 const MultistepAlgorithms = Union{ - ABDF2, AB3, AB4, AB5, ABM32, ABM43, ABM54} -const SplitAlgorithms = Union{CNAB2, CNLF2, SBDF, - KenCarp3, KenCarp4, KenCarp47, KenCarp5, KenCarp58, CFNLIRK3} +const SplitAlgorithms = Union{CNAB2, CNLF2} #= struct DBDF{CS,AD,F,F2,P,FDT,ST,CJ} <: DAEAlgorithm{CS,AD,FDT,ST,CJ} diff --git a/src/integrators/controllers.jl b/src/integrators/controllers.jl index 5db779bb62..3965cffec7 100644 --- a/src/integrators/controllers.jl +++ b/src/integrators/controllers.jl @@ -57,6 +57,9 @@ the predicted step size. struct IController <: AbstractController end +struct DummyController <: AbstractController +end + @inline function stepsize_controller!(integrator, controller::IController, alg) @unpack qmin, qmax, gamma = integrator.opts EEst = DiffEqBase.value(integrator.EEst) @@ -448,8 +451,6 @@ end # Dummy controller without any method implementations. # This is used to transfer the special controllers associated to certain # algorithms to the new controller infrastructure with -struct DummyController <: AbstractController -end # JVODE function stepsize_controller!(integrator, alg::JVODE) @@ -474,148 +475,16 @@ function step_reject_controller!(integrator, alg::JVODE) integrator.dt *= integrator.qold end -# QNBDF -stepsize_controller!(integrator, alg::QNDF) = nothing - -# this stepsize and order controller is taken from -# Implementation of an Adaptive BDF2 Formula and Comparison with the MATLAB Ode15s paper -# E. Alberdi Celaya, J. J. Anza Aguirrezabala, and P. Chatzipantelidis - -function step_accept_controller!(integrator, alg::QNDF{max_order}, q) where {max_order} - #step is accepted, reset count of consecutive failed steps - integrator.cache.consfailcnt = 0 - integrator.cache.nconsteps += 1 - if iszero(integrator.EEst) - return integrator.dt * integrator.opts.qmax - else - est = integrator.EEst - estₖ₋₁ = integrator.cache.EEst1 - estₖ₊₁ = integrator.cache.EEst2 - h = integrator.dt - k = integrator.cache.order - cache = integrator.cache - prefer_const_step = integrator.cache.nconsteps < integrator.cache.order + 2 - zₛ = 1.2 # equivalent to integrator.opts.gamma - zᵤ = 0.1 - Fᵤ = 10 - expo = 1 / (k + 1) - z = zₛ * ((est)^expo) - F = inv(z) - hₙ = h - kₙ = k - if z <= zᵤ - hₖ = Fᵤ * h - else - hₖ = F * h - end - hₖ₋₁ = 0.0 - hₖ₊₁ = 0.0 - - if k > 1 - expo = 1 / k - zₖ₋₁ = 1.3 * ((estₖ₋₁)^expo) - Fₖ₋₁ = inv(zₖ₋₁) - if zₖ₋₁ <= 0.1 - hₖ₋₁ = 10 * h - elseif 1 / 10 < zₖ₋₁ <= 1.3 - hₖ₋₁ = Fₖ₋₁ * h - end - if hₖ₋₁ > hₖ - hₙ = hₖ₋₁ - kₙ = k - 1 - else - hₙ = hₖ - kₙ = k - end - else - hₙ = hₖ - kₙ = k - end - - if k < max_order - expo = 1 / (k + 2) - zₖ₊₁ = 1.4 * ((estₖ₊₁)^expo) - Fₖ₊₁ = inv(zₖ₊₁) - - if zₖ₊₁ <= 0.1 - hₖ₊₁ = 10 * h - elseif 0.1 < zₖ₊₁ <= 1.4 - hₖ₊₁ = Fₖ₊₁ * h - end - if hₖ₊₁ > hₙ - hₙ = hₖ₊₁ - kₙ = k + 1 - end - end - cache.order = kₙ - q = integrator.dt / hₙ - end - if prefer_const_step - if q < 1.2 && q > 0.6 - return integrator.dt - end - end - if q <= integrator.opts.qsteady_max && q >= integrator.opts.qsteady_min - return integrator.dt - end - return integrator.dt / q -end - -function step_reject_controller!(integrator, ::QNDF) - bdf_step_reject_controller!(integrator, integrator.cache.EEst1) -end - -function step_reject_controller!(integrator, ::Union{FBDF, DFBDF}) +function step_reject_controller!(integrator, ::DFBDF) bdf_step_reject_controller!(integrator, integrator.cache.terkm1) end -function bdf_step_reject_controller!(integrator, EEst1) - k = integrator.cache.order - h = integrator.dt - integrator.cache.consfailcnt += 1 - integrator.cache.nconsteps = 0 - if integrator.cache.consfailcnt > 1 - h = h / 2 - end - zₛ = 1.2 # equivalent to integrator.opts.gamma - expo = 1 / (k + 1) - z = zₛ * ((integrator.EEst)^expo) - F = inv(z) - if z <= 10 - hₖ = F * h - else # z > 10 - hₖ = 0.1 * h - end - hₙ = hₖ - kₙ = k - if k > 1 - expo = 1 / k - zₖ₋₁ = 1.3 * (EEst1^expo) - Fₖ₋₁ = inv(zₖ₋₁) - if zₖ₋₁ <= 10 - hₖ₋₁ = Fₖ₋₁ * h - elseif zₖ₋₁ > 10 - hₖ₋₁ = 0.1 * h - end - if integrator.cache.consfailcnt > 2 || hₖ₋₁ > hₖ - hₙ = min(h, hₖ₋₁) - kₙ = k - 1 - end - end - # Restart BDf (clear history) when we failed repeatedly - if kₙ == 1 && integrator.cache.consfailcnt > 3 - u_modified!(integrator, true) - end - integrator.dt = hₙ - integrator.cache.order = kₙ -end - function post_newton_controller!(integrator, alg) integrator.dt = integrator.dt / integrator.opts.failfactor nothing end -function post_newton_controller!(integrator, alg::Union{FBDF, DFBDF}) +function post_newton_controller!(integrator, alg::DFBDF) @unpack cache = integrator if cache.order > 1 && cache.nlsolver.nfails >= 3 cache.order -= 1 @@ -626,7 +495,7 @@ function post_newton_controller!(integrator, alg::Union{FBDF, DFBDF}) nothing end -function choose_order!(alg::Union{FBDF, DFBDF}, integrator, +function choose_order!(alg::DFBDF, integrator, cache::OrdinaryDiffEqMutableCache, ::Val{max_order}) where {max_order} @unpack t, dt, u, cache, uprev = integrator @@ -662,7 +531,7 @@ function choose_order!(alg::Union{FBDF, DFBDF}, integrator, return k, terk end -function choose_order!(alg::Union{FBDF, DFBDF}, integrator, +function choose_order!(alg::DFBDF, integrator, cache::OrdinaryDiffEqConstantCache, ::Val{max_order}) where {max_order} @unpack t, dt, u, cache, uprev = integrator @@ -707,7 +576,7 @@ function choose_order!(alg::Union{FBDF, DFBDF}, integrator, end function stepsize_controller!(integrator, - alg::Union{FBDF{max_order}, DFBDF{max_order}}) where { + alg::DFBDF{max_order}) where { max_order, } @unpack cache = integrator @@ -726,7 +595,7 @@ function stepsize_controller!(integrator, q end -function step_accept_controller!(integrator, alg::Union{FBDF{max_order}, DFBDF{max_order}}, +function step_accept_controller!(integrator, alg::DFBDF{max_order}, q) where {max_order} integrator.cache.consfailcnt = 0 if q <= integrator.opts.qsteady_max && q >= integrator.opts.qsteady_min diff --git a/src/perform_step/rosenbrock_perform_step.jl b/src/perform_step/rosenbrock_perform_step.jl index 8cd2ffd1f3..67db01d3fc 100644 --- a/src/perform_step/rosenbrock_perform_step.jl +++ b/src/perform_step/rosenbrock_perform_step.jl @@ -1947,7 +1947,7 @@ end if integrator.opts.adaptive if (integrator.alg isa Rodas5Pe) - du = 0.2606326497975715 * k1 - 0.005158627295444251 * k2 + + @. du = 0.2606326497975715 * k1 - 0.005158627295444251 * k2 + 1.3038988631109731 * k3 + 1.235000722062074 * k4 + -0.7931985603795049 * k5 - 1.005448461135913 * k6 - 0.18044626132120234 * k7 + 0.17051519239113755 * k8 diff --git a/test/runtests.jl b/test/runtests.jl index 1dbd745f77..2bce643dd9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -71,6 +71,12 @@ function activate_symplectic_rk() Pkg.instantiate() end +function activate_sdirk() + Pkg.activate("../lib/OrdinaryDiffEqSDIRK") + Pkg.develop(PackageSpec(path = dirname(@__DIR__))) + Pkg.instantiate() +end + #Start Test Script @time begin