Merge pull request #462 from sipah00/sbdf2

SBDF2, SBDF3 and SBDF4.
SciML · Aug 2, 2018 · fe202c8 · fe202c8
2 parents 1716b77 + b3bb1b0
commit fe202c8
Show file tree

Hide file tree

Showing 6 changed files with 851 additions and 6 deletions.
diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl
@@ -66,6 +66,7 @@ module OrdinaryDiffEq
   include("caches/adams_bashforth_moulton_caches.jl")
   include("caches/nordsieck_caches.jl")
   include("caches/bdf_caches.jl")
+  include("caches/sbdf_caches.jl")
   include("caches/rkc_caches.jl")
   include("caches/euler_imex_caches.jl")
 
@@ -107,6 +108,7 @@ module OrdinaryDiffEq
   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/sbdf_perform_step.jl")
   include("perform_step/rkc_perform_step.jl")
   include("perform_step/euler_imex_perform_step.jl")
 
@@ -194,6 +196,8 @@ module OrdinaryDiffEq
 
   export ABDF2, QNDF1, QBDF1, QNDF2, QBDF2, QNDF, QBDF
 
+  export SBDF2, SBDF3, SBDF4
+
   export AutoSwitch, AutoTsit5, AutoDP5,
          AutoVern6, AutoVern7, AutoVern8, AutoVern9
 end # module
diff --git a/src/alg_utils.jl b/src/alg_utils.jl
@@ -251,6 +251,10 @@ alg_order(alg::QNDF2) = 2
 
 alg_order(alg::QNDF) = 1 #dummy value
 
+alg_order(alg::SBDF2) = 2
+alg_order(alg::SBDF3) = 3
+alg_order(alg::SBDF4) = 4
+
 alg_order(alg::ROCK2) = 2
 
 alg_maximum_order(alg) = alg_order(alg)

diff --git a/src/algorithms.jl b/src/algorithms.jl
@@ -326,6 +326,66 @@ Base.@pure QNDF(;chunk_size=0,autodiff=true,diff_type=Val{:central},
 
 Base.@pure QBDF(;kwargs...) = QNDF(;kappa=tuple(0,0,0,0,0),kwargs...)
 
+struct SBDF2{CS,AD,F,F2,FDT,K,T,T2} <: OrdinaryDiffEqNewtonAlgorithm{CS,AD}
+  linsolve::F
+  nonlinsolve::F2
+  diff_type::FDT
+  κ::K
+  tol::T
+  extrapolant::Symbol
+  min_newton_iter::Int
+  max_newton_iter::Int
+  new_jac_conv_bound::T2
+end
+SBDF2(;chunk_size=0,autodiff=true,diff_type=Val{:central},
+                      linsolve=DEFAULT_LINSOLVE,nonlinsolve=Val{:newton},κ=nothing,tol=nothing,
+                      extrapolant=:linear,min_newton_iter=1,
+                      max_newton_iter=7,new_jac_conv_bound = 1e-3) =
+                      SBDF2{chunk_size,autodiff,typeof(linsolve),typeof(nonlinsolve),typeof(diff_type),
+                      typeof(κ),typeof(tol),typeof(new_jac_conv_bound)}(
+                      linsolve,nonlinsolve,diff_type,κ,tol,extrapolant,min_newton_iter,
+                      max_newton_iter,new_jac_conv_bound)
+
+struct SBDF3{CS,AD,F,F2,FDT,K,T,T2} <: OrdinaryDiffEqNewtonAlgorithm{CS,AD}
+  linsolve::F
+  nonlinsolve::F2
+  diff_type::FDT
+  κ::K
+  tol::T
+  extrapolant::Symbol
+  min_newton_iter::Int
+  max_newton_iter::Int
+  new_jac_conv_bound::T2
+end
+SBDF3(;chunk_size=0,autodiff=true,diff_type=Val{:central},
+                      linsolve=DEFAULT_LINSOLVE,nonlinsolve=Val{:newton},κ=nothing,tol=nothing,
+                      extrapolant=:linear,min_newton_iter=1,
+                      max_newton_iter=7,new_jac_conv_bound = 1e-3) =
+                      SBDF3{chunk_size,autodiff,typeof(linsolve),typeof(nonlinsolve),typeof(diff_type),
+                      typeof(κ),typeof(tol),typeof(new_jac_conv_bound)}(
+                      linsolve,nonlinsolve,diff_type,κ,tol,extrapolant,min_newton_iter,
+                      max_newton_iter,new_jac_conv_bound)
+
+struct SBDF4{CS,AD,F,F2,FDT,K,T,T2} <: OrdinaryDiffEqNewtonAlgorithm{CS,AD}
+  linsolve::F
+  nonlinsolve::F2
+  diff_type::FDT
+  κ::K
+  tol::T
+  extrapolant::Symbol
+  min_newton_iter::Int
+  max_newton_iter::Int
+  new_jac_conv_bound::T2
+end
+SBDF4(;chunk_size=0,autodiff=true,diff_type=Val{:central},
+                      linsolve=DEFAULT_LINSOLVE,nonlinsolve=Val{:newton},κ=nothing,tol=nothing,
+                      extrapolant=:linear,min_newton_iter=1,
+                      max_newton_iter=7,new_jac_conv_bound = 1e-3) =
+                      SBDF4{chunk_size,autodiff,typeof(linsolve),typeof(nonlinsolve),typeof(diff_type),
+                      typeof(κ),typeof(tol),typeof(new_jac_conv_bound)}(
+                      linsolve,nonlinsolve,diff_type,κ,tol,extrapolant,min_newton_iter,
+                      max_newton_iter,new_jac_conv_bound)
+
 # Adams/BDF methods in Nordsieck forms
 struct AN5   <: OrdinaryDiffEqAdaptiveAlgorithm end
 struct JVODE{bType,aType} <: OrdinaryDiffEqAdamsVarOrderVarStepAlgorithm