Merge pull request #399 from sipah00/ndf

1st order Numerical Differentiation Method (QNDF1)

src/OrdinaryDiffEq.jl

@@ -119,6 +119,7 @@ module OrdinaryDiffEq
   include("nlsolve_utils.jl")
   include("nordsieck_utils.jl")
   include("adams_utils.jl")
+  include("bdf_utils.jl")
   include("exponential_utils.jl")
   include("derivative_wrappers.jl")
   include("iterator_interface.jl")
@@ -184,7 +185,7 @@ module OrdinaryDiffEq
 
   export AN5, JVODE, JVODE_Adams
 
-  export ABDF2
+  export ABDF2, QNDF1, QBDF1
 
   export AutoSwitch, AutoTsit5, AutoDP5,
          AutoVern6, AutoVern7, AutoVern8, AutoVern9

src/alg_utils.jl

@@ -237,6 +237,7 @@ alg_order(alg::AN5) = 5
 alg_order(alg::JVODE) = 1  #dummy value
 
 alg_order(alg::ABDF2) = 2
+alg_order(alg::QNDF1) = 1
 
 alg_maximum_order(alg) = alg_order(alg)
 alg_maximum_order(alg::CompositeAlgorithm) = maximum(alg_order(x) for x in alg.algs)
@@ -284,6 +285,7 @@ qsteady_max_default(alg::OrdinaryDiffEqAdaptiveImplicitAlgorithm) = 6//5
 qsteady_max_default(alg::OrdinaryDiffEqImplicitAlgorithm) = 1//1
 qsteady_max_default(alg::AN5) = 3//2
 qsteady_max_default(alg::JVODE) = 3//2
+qsteady_max_default(alg::QNDF1) = 2//1
 
 FunctionMap_scale_by_time{scale_by_time}(alg::FunctionMap{scale_by_time}) = scale_by_time
 

src/algorithms.jl

@@ -229,6 +229,29 @@ Base.@pure CNLF2(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                       linsolve,diff_type,κ,tol,extrapolant,min_newton_iter,
                       max_newton_iter,new_jac_conv_bound)
 
+struct QNDF1{CS,AD,F,FDT,K,T,T2,κType,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}
+  linsolve::F
+  diff_type::FDT
+  κ::K
+  tol::T
+  extrapolant::Symbol
+  min_newton_iter::Int
+  max_newton_iter::Int
+  new_jac_conv_bound::T2
+  kappa::κType
+end
+Base.@pure QNDF1(;chunk_size=0,autodiff=true,diff_type=Val{:central},
+                      linsolve=DEFAULT_LINSOLVE,κ=nothing,tol=nothing,
+                      extrapolant=:linear,min_newton_iter=1,
+                      max_newton_iter=7,new_jac_conv_bound = 1e-3, kappa = -0.1850,
+                      controller = :Predictive) =
+                      QNDF1{chunk_size,autodiff,typeof(linsolve),typeof(diff_type),
+                      typeof(κ),typeof(tol),typeof(new_jac_conv_bound),typeof(kappa),controller}(
+                      linsolve,diff_type,κ,tol,extrapolant,min_newton_iter,
+                      max_newton_iter,new_jac_conv_bound,kappa)
+
+Base.@pure QBDF1(;kwargs...) = QNDF1(;kappa=0,kwargs...)
+
 # Adams/BDF methods in Nordsieck forms
 struct AN5   <: OrdinaryDiffEqAdaptiveAlgorithm end
 struct JVODE{bType,aType} <: OrdinaryDiffEqAdamsVarOrderVarStepAlgorithm

src/bdf_utils.jl

@@ -0,0 +1,46 @@
+# bdf_utils
+
+function U!(k, cache)
+  @unpack U = cache
+  for j = 1:k
+    for r = 1:k
+      if typeof(cache) <: OrdinaryDiffEqMutableCache
+        @. U[j,r] = inv(factorial(j)) * prod([m-r for m in 0:(j-1)])
+      else
+        U[j,r] = inv(factorial(j)) * prod([m-r for m in 0:(j-1)]) 
+      end
+    end
+  end  
+end
+
+function R!(k, ρ, cache)
+  @unpack R = cache
+  for j = 1:k
+    for r = 1:k
+      if typeof(cache) <: OrdinaryDiffEqMutableCache
+        @. R[j,r] = inv(factorial(j)) * prod([m-r*ρ for m in 0:(j-1)])
+      else
+        R[j,r] = inv(factorial(j)) * prod([m-r*ρ for m in 0:(j-1)])
+      end
+    end
+  end
+end
+
+function prev_u!(uprev, k, t, dt, cache)
+  @unpack D, uprev2 = cache
+  if typeof(cache) <: OrdinaryDiffEqMutableCache
+    @. uprev2 = uprev + (D[1,1]) * -1
+  else
+    uprev2 = uprev + (D[1,1]) * -1
+    return uprev2
+  end
+end
+
+function D2!(u, uprev, k, cache)
+  @unpack D, D2 = cache
+  if typeof(cache) <: OrdinaryDiffEqMutableCache
+    @. D2[1,1] = (u - uprev) - D[1,1] 
+  else
+    D2[1,1] = (u - uprev) - D[1,1] 
+  end
+end

src/caches/bdf_caches.jl

@@ -102,3 +102,144 @@ function alg_cache(alg::ABDF2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUni
   ABDF2Cache(u,uprev,uprev2,du1,fsalfirst,fsalfirstprev,k,z,zₙ₋₁,dz,b,tmp,atmp,J,
               W,uf,jac_config,linsolve,ηold,κ,tol,10000,eulercache,dtₙ₋₁)
 end
+
+# QNDF1
+
+mutable struct QNDF1ConstantCache{F,uToltype,coefType,coefType1,dtType,uType,tType} <: OrdinaryDiffEqConstantCache
+  uf::F
+  ηold::uToltype
+  κ::uToltype
+  tol::uToltype
+  newton_iters::Int
+  eulercache::ImplicitEulerConstantCache
+  D::coefType
+  temp_D::coefType
+  D2::coefType1
+  R::coefType
+  U::coefType
+  uprev2::uType
+  uprev3::uType
+  tprev2::tType
+  dtₙ₋₁::dtType
+end
+
+mutable struct QNDF1Cache{uType,rateType,coefType,uNoUnitsType,J,UF,JC,uToltype,F,dtType} <: OrdinaryDiffEqMutableCache
+  uprev2::uType
+  du1::rateType
+  fsalfirst::rateType
+  fsalfirstprev::rateType
+  k::rateType
+  z::uType
+  zₙ₋₁::uType
+  dz::uType
+  b::uType
+  D::coefType
+  temp_D::coefType
+  temp_u::uType
+  D2::coefType
+  R::coefType
+  U::coefType
+  tmp::uType
+  atmp::uNoUnitsType
+  utilde::uType
+  J::J
+  W::J
+  uf::UF
+  jac_config::JC
+  linsolve::F
+  ηold::uToltype
+  κ::uToltype
+  tol::uToltype
+  newton_iters::Int
+  eulercache::ImplicitEulerCache
+  dtₙ₋₁::dtType
+end
+
+u_cache(c::QNDF1Cache)    = ()
+du_cache(c::QNDF1Cache)   = ()
+
+function alg_cache(alg::QNDF1,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
+  uToltype = real(uBottomEltypeNoUnits)
+  uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
+  ηold = one(uToltype)
+  uprev2 = u
+  uprev3 = u
+  tprev2 = t
+  dtₙ₋₁ = t
+
+  D = zeros(typeof(t), 1, 1)
+  temp_D = zeros(typeof(t), 1, 1)
+  D2 = zeros(typeof(t), 1, 1)
+  R = zeros(typeof(t), 1, 1)
+  U = zeros(typeof(t), 1, 1)
+
+  if alg.κ != nothing
+    κ = uToltype(alg.κ)
+  else
+    κ = uToltype(1//100)
+  end
+  if alg.tol != nothing
+    tol = uToltype(alg.tol)
+  else
+    tol = uToltype(min(0.03,first(reltol)^(0.5)))
+  end
+
+  eulercache = ImplicitEulerConstantCache(uf,ηold,κ,tol,100000)
+
+  QNDF1ConstantCache(uf,ηold,κ,tol,10000,eulercache,D,temp_D,D2,R,U,uprev2,uprev3,tprev2,dtₙ₋₁)
+end
+
+function alg_cache(alg::QNDF1,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
+  du1 = zeros(rate_prototype)
+  J = zeros(uEltypeNoUnits,length(u),length(u))
+  W = similar(J)
+  zprev = similar(u,indices(u))
+  zₙ₋₁ = similar(u,indices(u)); z = similar(u,indices(u))
+  dz = similar(u,indices(u))
+  fsalfirst = zeros(rate_prototype)
+  fsalfirstprev = zeros(rate_prototype)
+  k = zeros(rate_prototype)
+
+  D = Vector{typeof(u)}(1)
+  R = Vector{typeof(u)}(1)
+  U = Vector{typeof(u)}(1)
+  D2 = Vector{typeof(u)}(1)
+  temp_D = Vector{typeof(u)}(1)
+
+  D[1,1] = similar(u)
+  R[1,1] = similar(u)
+  U[1,1] = similar(u)
+  D2[1,1] = similar(u)
+  temp_D[1,1] = similar(u)
+  temp_u = similar(u)
+
+  tmp = similar(u); b = similar(u,indices(u))
+  atmp = similar(u,uEltypeNoUnits,indices(u))
+
+  uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
+  linsolve = alg.linsolve(Val{:init},uf,u)
+  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,dz)
+  uToltype = real(uBottomEltypeNoUnits)
+  utilde = similar(u,indices(u))
+
+  uprev2 = similar(u)
+
+  if alg.κ != nothing
+    κ = uToltype(alg.κ)
+  else
+    κ = uToltype(1//100)
+  end
+  if alg.tol != nothing
+    tol = uToltype(alg.tol)
+  else
+    tol = uToltype(min(0.03,first(reltol)^(0.5)))
+  end
+
+  ηold = one(uToltype)
+
+  eulercache = ImplicitEulerCache(u,uprev,uprev2,du1,fsalfirst,k,z,dz,b,tmp,atmp,J,W,uf,jac_config,linsolve,ηold,κ,tol,10000)
+
+  dtₙ₋₁ = one(dt)
+  QNDF1Cache(uprev2,du1,fsalfirst,fsalfirstprev,k,z,zₙ₋₁,dz,b,D,temp_D,temp_u,D2,R,U,tmp,atmp,utilde,J,
+              W,uf,jac_config,linsolve,ηold,κ,tol,10000,eulercache,dtₙ₋₁)
+end

src/perform_step/bdf_perform_step.jl

@@ -143,3 +143,139 @@ end
   end
   return
 end
+
+# QNDF1
+
+function initialize!(integrator, cache::QNDF1ConstantCache)
+  integrator.kshortsize = 2
+  integrator.k = typeof(integrator.k)(integrator.kshortsize)
+  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
+
+  # Avoid undefined entries if k is an array of arrays
+  integrator.fsallast = zero(integrator.fsalfirst)
+  integrator.k[1] = integrator.fsalfirst
+  integrator.k[2] = integrator.fsallast
+end
+
+function perform_step!(integrator,cache::QNDF1ConstantCache,repeat_step=false)
+  @unpack t,dt,uprev,u,f,p = integrator
+  @unpack uprev2,D,temp_D,D2,R,U,dtₙ₋₁,uf,κ,tol = cache
+  cnt = integrator.iter
+  if cnt == 1
+    perform_step!(integrator, cache.eulercache, repeat_step)
+    cache.uprev2 = integrator.uprev
+    cache.dtₙ₋₁ = dt
+    return
+  end
+  κ = integrator.alg.kappa
+  γₖ = 1//1
+  k = 1
+  ρ = dt/dtₙ₋₁
+  D[1,1] = uprev - uprev2 # backward diff
+  temp_D .= D
+  if ρ != 1
+    U!(k,cache)
+    R!(k,ρ,cache)
+    D[1,1] = D[1,1] * (R[1,1] * U[1,1])
+    uprev2 = prev_u!(uprev,k,t,dt,cache)
+  end
+  c1 = (1-2κ)/(1-κ)
+  c2 = (κ)/(1-κ)
+  tmp = c1*uprev + c2*uprev2
+  # precalculations
+  γ = inv(1-κ)
+  γdt = γ*dt
+  W = calc_W!(integrator, cache, γdt, repeat_step)
+
+  # initial guess
+  z = dt*integrator.fsalfirst
+
+  z, η, iter, fail_convergence = diffeq_nlsolve!(integrator, cache, W, z, tmp, γ, 1, Val{:newton})
+  fail_convergence && return
+  u = tmp + γ*z
+
+  if integrator.opts.adaptive
+    D2!(u,uprev,k,cache)
+    utilde = (κ*γₖ + inv(k+1)) * D2[1,1]
+    atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm)
+    integrator.EEst = integrator.opts.internalnorm(atmp)
+    if integrator.EEst > one(integrator.EEst)
+      cache.D .= temp_D
+      return
+    end
+  end
+  cache.dtₙ₋₁ = dt
+  cache.uprev2 = uprev
+  cache.ηold = η
+  cache.newton_iters = iter
+  integrator.fsallast = f(u, p, t+dt)
+  integrator.k[1] = integrator.fsalfirst
+  integrator.k[2] = integrator.fsallast
+  integrator.u = u
+end
+
+function initialize!(integrator, cache::QNDF1Cache)
+  integrator.kshortsize = 2
+  integrator.fsalfirst = cache.fsalfirst
+  integrator.fsallast = cache.k
+  resize!(integrator.k, integrator.kshortsize)
+  integrator.k[1] = integrator.fsalfirst
+  integrator.k[2] = integrator.fsallast
+  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # For the interpolation, needs k at the updated point
+end
+
+function perform_step!(integrator,cache::QNDF1Cache,repeat_step=false)
+  @unpack t,dt,uprev,u,f,p = integrator
+  @unpack uprev2,D,temp_D,temp_u,D2,R,U,tmp,z,W,dtₙ₋₁,uf,κ,tol,utilde,atmp = cache
+  cnt = integrator.iter
+  if cnt == 1
+    perform_step!(integrator, cache.eulercache, repeat_step)
+    cache.uprev2 .= integrator.uprev
+    cache.dtₙ₋₁ = dt
+    return
+  end
+  κ = integrator.alg.kappa
+  γₖ = 1//1
+  k = 1
+  ρ = dt/dtₙ₋₁
+  @. D[1,1] = uprev - uprev2 # backward diff
+  temp_D .= D
+  temp_u .= uprev2
+  if ρ != 1
+    U!(k,cache)
+    R!(k,ρ,cache)
+    @. D[1,1] = D[1,1] * (R[1,1] * U[1,1])
+    prev_u!(uprev,k,t,dt,cache)
+  end
+  c1 = (1-2κ)/(1-κ)
+  c2 = (κ)/(1-κ)
+  @. tmp = c1*uprev + c2*uprev2
+  # precalculations
+  γ = inv(1-κ)
+  γdt = γ*dt
+  new_W = calc_W!(integrator, cache, γdt, repeat_step)
+
+  # initial guess
+  @. z = dt*integrator.fsalfirst
+
+  z, η, iter, fail_convergence = diffeq_nlsolve!(integrator, cache, W, z, tmp, γ, 1, Val{:newton}, new_W)
+  fail_convergence && return
+  @. u = tmp + γ*z
+
+  if integrator.opts.adaptive
+    D2!(u,uprev,k,cache)
+    @. utilde = (κ*γₖ + inv(k+1)) * D2[1,1]
+    calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm)
+    integrator.EEst = integrator.opts.internalnorm(atmp)
+    if integrator.EEst > one(integrator.EEst)
+      cache.D .= temp_D
+      cache.uprev2 .= temp_u
+      return
+    end
+  end
+  cache.dtₙ₋₁ = dt
+  cache.uprev2 .= uprev
+  cache.ηold = η
+  cache.newton_iters = iter
+  f(integrator.fsallast, u, p, t+dt)
+end