Merge pull request #42 from ranocha/hr/controllers

adapt to new controller infrastructure

Author: ChrisRackauckas

src/StochasticDelayDiffEq.jl

@@ -6,6 +6,7 @@ module StochasticDelayDiffEq
 # using DifferentialEquations
 using Reexport
 @reexport using StochasticDiffEq
+import StochasticDiffEq: stepsize_controller!, accept_step_controller, step_accept_controller!, step_reject_controller!, PIController
 # import StochasticDiffEq: calc_J, calc_J!, calc_tderivative!, calc_tderivative
 using LinearAlgebra, StaticArrays
 using UnPack, DataStructures
@@ -33,5 +34,6 @@ include("integrators/utils.jl")
 include("functionwrapper.jl")
 include("utils.jl")
 include("solve.jl")
+include("stepsize_controllers.jl")
 
 end # module

src/integrators/interface.jl

@@ -11,7 +11,7 @@
         if integrator.isout
           integrator.dtnew = integrator.dt*integrator.opts.qmin
         elseif !integrator.force_stepfail
-          integrator.dtnew = integrator.dt/min(inv(integrator.opts.qmin),integrator.q11/integrator.opts.gamma)
+          step_reject_controller!(integrator,getalg(integrator.alg))
         end
         StochasticDiffEq.choose_algorithm!(integrator,integrator.cache)
         StochasticDiffEq.fix_dtnew_at_bounds!(integrator)
@@ -37,13 +37,11 @@
     integrator.last_stepfail = true
     integrator.accept_step = false
   elseif integrator.opts.adaptive
-    @fastmath integrator.q11 = integrator.EEst^integrator.opts.beta1
-    @fastmath integrator.q = integrator.q11/(integrator.qold^integrator.opts.beta2)
-    @fastmath integrator.q = max(inv(integrator.opts.qmax),min(inv(integrator.opts.qmin),integrator.q/integrator.opts.gamma))
-    @fastmath integrator.dtnew = integrator.dt/integrator.q
+    stepsize_controller!(integrator,getalg(integrator.alg))
     integrator.isout = integrator.opts.isoutofdomain(integrator.u,integrator.p,ttmp)
-    integrator.accept_step = (!integrator.isout && integrator.EEst <= 1.0) || (integrator.opts.force_dtmin && integrator.dt <= integrator.opts.dtmin)
+    integrator.accept_step = (!integrator.isout && accept_step_controller(integrator, integrator.opts.controller)) || (integrator.opts.force_dtmin && integrator.dt <= integrator.opts.dtmin)
     if integrator.accept_step # Accepted
+      step_accept_controller!(integrator,getalg(integrator.alg))
       integrator.last_stepfail = false
       integrator.tprev = integrator.t
       if typeof(integrator.t)<:AbstractFloat && !isempty(integrator.opts.tstops)

src/solve.jl

@@ -32,10 +32,13 @@ function DiffEqBase.__init(prob::AbstractSDDEProblem,# TODO DiffEqBasee.Abstract
     reltol = nothing,
     qmax = StochasticDiffEq.qmax_default(getalg(alg)),
     qmin = StochasticDiffEq.qmin_default(getalg(alg)),
+    qsteady_min = StochasticDiffEq.qsteady_min_default(alg),
+    qsteady_max = StochasticDiffEq.qsteady_max_default(alg),
     qoldinit = 1 // 10^4, fullnormalize = true,
+    controller = nothing,
     failfactor = 2,
-    beta2 = StochasticDiffEq.beta2_default(getalg(alg)),
-    beta1 = StochasticDiffEq.beta1_default(getalg(alg), beta2),
+    beta2 = nothing,
+    beta1 = nothing,
     delta = StochasticDiffEq.delta_default(getalg(alg)),
     maxiters = adaptive ? 1000000 : typemax(Int),
     dtmax = eltype(prob.tspan)((prob.tspan[end] - prob.tspan[1])),
@@ -334,22 +337,41 @@ function DiffEqBase.__init(prob::AbstractSDDEProblem,# TODO DiffEqBasee.Abstract
     save_end_user = save_end
     save_end = save_end === nothing ? save_everystep || isempty(saveat) || saveat isa Number || prob.tspan[2] in saveat : save_end
 
+    # Setting up the step size controller
+    if (beta1 !== nothing || beta2 !== nothing) && controller !== nothing
+      throw(ArgumentError(
+        "Setting both the legacy PID parameters `beta1, beta2 = $((beta1, beta2))` and the `controller = $controller` is not allowed."))
+    end
+
+    if (beta1 !== nothing || beta2 !== nothing)
+      message = "Providing the legacy PID parameters `beta1, beta2` is deprecated. Use the keyword argument `controller` instead."
+      Base.depwarn(message, :init)
+      Base.depwarn(message, :solve)
+    end
+
+    if controller === nothing
+      controller = StochasticDiffEq.default_controller(getalg(alg), cache, qoldinit, beta1, beta2)
+    end
+
     opts = StochasticDiffEq.SDEOptions(maxiters, save_everystep,
                       adaptive, abstol_internal,
-                      reltol_internal, QT(gamma),
+                      reltol_internal,
+                      QT(gamma),
                       QT(qmax), QT(qmin),
+                      QT(qsteady_max),QT(qsteady_min),
+                      QT(qoldinit),
                       QT(failfactor),
-                      tType(dtmax), tType(dtmin), internalnorm, save_idxs,
+                      tType(dtmax), tType(dtmin),
+                      controller,
+                      internalnorm, save_idxs,
                       tstops_internal, saveat_internal,
                       d_discontinuities_internal,
                       tstops, saveat, d_discontinuities,
                       userdata,
                       progress, progress_steps,
                       progress_name, progress_message,
                       timeseries_errors, dense_errors,
-                      QT(beta1), QT(beta2),
                       convert.(uBottomEltypeNoUnits, delta),
-                      QT(qoldinit),
                       dense, save_on, save_start, save_end, save_end_user, save_noise,
                       callbacks_internal, isoutofdomain, unstable_check,
                       verbose, calck, force_dtmin,

src/stepsize_controllers.jl

@@ -0,0 +1,18 @@
+
+function stepsize_controller!(integrator::SDDEIntegrator, controller::PIController, alg)
+  integrator.q11 = DiffEqBase.value(DiffEqBase.fastpow(integrator.EEst,controller.beta1))
+  integrator.q = DiffEqBase.value(integrator.q11/DiffEqBase.fastpow(integrator.qold,controller.beta2))
+  @fastmath integrator.q = DiffEqBase.value(max(inv(integrator.opts.qmax),min(inv(integrator.opts.qmin),integrator.q/integrator.opts.gamma)))
+end
+
+@inline function step_accept_controller!(integrator::SDDEIntegrator, alg)
+  step_accept_controller!(integrator, integrator.opts.controller, alg)
+end
+
+function step_accept_controller!(integrator::SDDEIntegrator, controller::PIController, alg)
+  integrator.dtnew = DiffEqBase.value(integrator.dt/integrator.q) * oneunit(integrator.dt)
+end
+
+function step_reject_controller!(integrator::SDDEIntegrator, controller::PIController, alg)
+  integrator.dtnew = integrator.dt/min(inv(integrator.opts.qmin),integrator.q11/integrator.opts.gamma)
+end