Skip to content

Commit

Permalink
Merge pull request SciML#2286 from ParamThakkar123/FIRK
Browse files Browse the repository at this point in the history
Added FIRK solvers
  • Loading branch information
ChrisRackauckas authored Jul 26, 2024
2 parents acce9cc + 66d3a86 commit 4fd5afb
Show file tree
Hide file tree
Showing 19 changed files with 308 additions and 236 deletions.
2 changes: 1 addition & 1 deletion lib/OrdinaryDiffEqExtrapolation/test/runtests.jl
Original file line number Diff line number Diff line change
@@ -1,3 +1,3 @@
using SafeTestsets

@time @safetestset "Extrapolation Tests" include("ode_extrapolation_tests.jl")
@time @safetestset "Extrapolation Tests" include("ode_extrapolation_tests.jl")
24 changes: 24 additions & 0 deletions lib/OrdinaryDiffEqFIRK/Project.toml
Original file line number Diff line number Diff line change
@@ -0,0 +1,24 @@
name = "OrdinaryDiffEqFIRK"
uuid = "5960d6e9-dd7a-4743-88e7-cf307b64f125"
authors = ["ParamThakkar123 <[email protected]>"]
version = "1.0.0"

[deps]
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
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"]
35 changes: 35 additions & 0 deletions lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,35 @@
module OrdinaryDiffEqFIRK

import OrdinaryDiffEq: alg_order, calculate_residuals!,
initialize!, perform_step!, @unpack, unwrap_alg,
calculate_residuals,
OrdinaryDiffEqAlgorithm, OrdinaryDiffEqNewtonAdaptiveAlgorithm,
OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache,
OrdinaryDiffEqAdaptiveAlgorithm, CompiledFloats, uses_uprev,
alg_cache, _vec, _reshape, @cache, isfsal, full_cache,
constvalue, _unwrap_val, du_alias_or_new,
explicit_rk_docstring, trivial_limiter!,
_ode_interpolant!, _ode_addsteps!, AbstractController,
NLStatus, qmax_default, alg_adaptive_order, DEFAULT_PRECS,
UJacobianWrapper, build_J_W, build_jac_config, UDerivativeWrapper,
Convergence, calc_J!, dolinsolve, FastConvergence, calc_J,
stepsize_controller!, islinearfunction, step_accept_controller!, step_reject_controller!,
PredictiveController, alg_can_repeat_jac, NewtonAlgorithm, fac_default_gamma,
get_new_W_γdt_cutoff, get_current_adaptive_order, VerySlowConvergence,
Divergence, isfirk
using MuladdMacro, DiffEqBase, RecursiveArrayTools
using SciMLOperators: AbstractSciMLOperator
using LinearAlgebra: I, UniformScaling, mul!, lu
import LinearSolve
import FastBroadcast: @..
include("algorithms.jl")
include("alg_utils.jl")
include("controllers.jl")
include("firk_caches.jl")
include("firk_tableaus.jl")
include("firk_perform_step.jl")
include("integrator_interface.jl")

export RadauIIA3, RadauIIA5, RadauIIA7

end
13 changes: 13 additions & 0 deletions lib/OrdinaryDiffEqFIRK/src/alg_utils.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,13 @@
qmax_default(alg::Union{RadauIIA3, RadauIIA5, RadauIIA7}) = 8

alg_order(alg::RadauIIA3) = 3
alg_order(alg::RadauIIA5) = 5
alg_order(alg::RadauIIA7) = 7

isfirk(alg::RadauIIA3) = true
isfirk(alg::RadauIIA5) = true
isfirk(alg::RadauIIA7) = true

alg_adaptive_order(alg::RadauIIA3) = 1
alg_adaptive_order(alg::RadauIIA5) = 3
alg_adaptive_order(alg::RadauIIA7) = 5
155 changes: 155 additions & 0 deletions lib/OrdinaryDiffEqFIRK/src/algorithms.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,155 @@
# 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

"""
@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

"""
@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}
}
RadauIIA7: Fully-Implicit Runge-Kutta Method
An A-B-L stable fully implicit Runge-Kutta method with internal tableau complex basis transform for efficiency.
"""
struct RadauIIA7{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 RadauIIA7(; 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!)
RadauIIA7{_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
51 changes: 51 additions & 0 deletions lib/OrdinaryDiffEqFIRK/src/controllers.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,51 @@
@inline function stepsize_controller!(integrator, controller::PredictiveController, alg)
@unpack qmin, qmax, gamma = integrator.opts
EEst = DiffEqBase.value(integrator.EEst)

if iszero(EEst)
q = inv(qmax)
else
if fac_default_gamma(alg)
fac = gamma
else
if isfirk(alg)
@unpack iter = integrator.cache
@unpack maxiters = alg
else
@unpack iter, maxiters = integrator.cache.nlsolver
end
fac = min(gamma, (1 + 2 * maxiters) * gamma / (iter + 2 * maxiters))
end
expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
qtmp = DiffEqBase.fastpow(EEst, expo) / fac
@fastmath q = DiffEqBase.value(max(inv(qmax), min(inv(qmin), qtmp)))
integrator.qold = q
end
q
end

function step_accept_controller!(integrator, controller::PredictiveController, alg, q)
@unpack qmin, qmax, gamma, qsteady_min, qsteady_max = integrator.opts
EEst = DiffEqBase.value(integrator.EEst)

if integrator.success_iter > 0
expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
qgus = (integrator.dtacc / integrator.dt) *
DiffEqBase.fastpow((EEst^2) / integrator.erracc, expo)
qgus = max(inv(qmax), min(inv(qmin), qgus / gamma))
qacc = max(q, qgus)
else
qacc = q
end
if qsteady_min <= qacc <= qsteady_max
qacc = one(qacc)
end
integrator.dtacc = integrator.dt
integrator.erracc = max(1e-2, EEst)
return integrator.dt / qacc
end

function step_reject_controller!(integrator, controller::PredictiveController, alg)
@unpack dt, success_iter, qold = integrator
integrator.dt = success_iter == 0 ? 0.1 * dt : dt / qold
end
Original file line number Diff line number Diff line change
Expand Up @@ -200,7 +200,6 @@ mutable struct RadauIIA5Cache{uType, cuType, uNoUnitsType, rateType, JType, W1Ty
status::NLStatus
step_limiter!::StepLimiter
end
TruncatedStacktraces.@truncate_stacktrace RadauIIA5Cache 1

function alg_cache(alg::RadauIIA5, u, rate_prototype, ::Type{uEltypeNoUnits},
::Type{uBottomEltypeNoUnits},
Expand Down Expand Up @@ -370,7 +369,6 @@ mutable struct RadauIIA7Cache{uType, cuType, uNoUnitsType, rateType, JType, W1Ty
status::NLStatus
step_limiter!::StepLimiter
end
TruncatedStacktraces.@truncate_stacktrace RadauIIA7Cache 1

function alg_cache(alg::RadauIIA7, u, rate_prototype, ::Type{uEltypeNoUnits},
::Type{uBottomEltypeNoUnits},
Expand Down
File renamed without changes.
File renamed without changes.
4 changes: 4 additions & 0 deletions lib/OrdinaryDiffEqFIRK/src/integrator_interface.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
@inline function DiffEqBase.get_tmp_cache(integrator, alg::Union{RadauIIA3, RadauIIA5, RadauIIA7},
cache::OrdinaryDiffEqMutableCache)
(cache.tmp, cache.atmp)
end
Original file line number Diff line number Diff line change
Expand Up @@ -57,4 +57,4 @@ for iip in (true, false)
@test sol.stats.njacs < sol.stats.nw # W reuse
end
@test length(sol) < 5000 # the error estimate is not very good
end
end
3 changes: 3 additions & 0 deletions lib/OrdinaryDiffEqFIRK/test/runtests.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,3 @@
using SafeTestsets

@time @safetestset "FIRK Tests" include("ode_firk_tests.jl")
3 changes: 3 additions & 0 deletions lib/OrdinaryDiffEqFeagin/test/runtests.jl
Original file line number Diff line number Diff line change
@@ -1 +1,4 @@

using SafeTestsets

@time @safetestset "Feagin Tests" include("ode_feagin_tests.jl")
10 changes: 5 additions & 5 deletions src/OrdinaryDiffEq.jl
Original file line number Diff line number Diff line change
Expand Up @@ -148,7 +148,6 @@ 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/firk_caches.jl")
include("caches/linear_caches.jl")
include("caches/linear_nonlinear_caches.jl")
include("caches/rosenbrock_caches.jl")
Expand All @@ -162,7 +161,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/firk_tableaus.jl")
include("tableaus/qprk_tableaus.jl")

include("integrators/type.jl")
Expand All @@ -181,7 +179,6 @@ 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/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")
Expand Down Expand Up @@ -284,11 +281,16 @@ include("../lib/OrdinaryDiffEqDefault/src/OrdinaryDiffEqDefault.jl")
using ..OrdinaryDiffEqDefault
export DefaultODEAlgorithm

include("../lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl")
using ..OrdinaryDiffEqFIRK
export RadauIIA3, RadauIIA5, RadauIIA7

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")

Expand Down Expand Up @@ -428,8 +430,6 @@ export FunctionMap, Euler, Heun, Ralston, Midpoint, RK4, ExplicitRK, OwrenZen3,
FRK65, PFRK87,
RKM, MSRK5, MSRK6, Stepanov5, SIR54, QPRK98, PSRK4p7q6, PSRK3p6q5, PSRK3p5q4

export RadauIIA3, RadauIIA5, RadauIIA7

export MagnusMidpoint, LinearExponential, MagnusLeapfrog, LieEuler, CayleyEuler,
MagnusGauss4, MagnusNC6, MagnusGL6, MagnusGL8, MagnusNC8, MagnusGL4,
MagnusAdapt4, RKMK2, RKMK4, LieRK4, CG2, CG3, CG4a
Expand Down
Loading

0 comments on commit 4fd5afb

Please sign in to comment.