diff --git a/lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl b/lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl index 2080050161..33078d809b 100644 --- a/lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl +++ b/lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl @@ -16,13 +16,14 @@ import OrdinaryDiffEq: alg_order, calculate_residuals!, DIRK, set_new_W!, du_alias_or_new, trivial_limiter!, ImplicitEulerConstantCache, + compute_step!, 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! + u_modified!, DAEAlgorithm, _unwrap_val using TruncatedStacktraces, MuladdMacro, MacroTools, FastBroadcast, RecursiveArrayTools import StaticArrays: SArray, MVector, SVector, @SVector, StaticArray, MMatrix, SA using LinearAlgebra: I @@ -32,10 +33,13 @@ include("algorithms.jl") include("alg_utils.jl") include("bdf_utils.jl") include("bdf_caches.jl") +include("dae_caches.jl") include("controllers.jl") +include("dae_perform_step.jl") include("bdf_perform_step.jl") export ABDF2, QNDF1, QBDF1, QNDF2, QBDF2, QNDF, QBDF, FBDF, - SBDF2, SBDF3, SBDF4, MEBDF2, IMEXEuler, IMEXEulerARK + SBDF2, SBDF3, SBDF4, MEBDF2, IMEXEuler, IMEXEulerARK, + DABDF2, DImplicitEuler, DFBDF end diff --git a/lib/OrdinaryDiffEqBDF/src/alg_utils.jl b/lib/OrdinaryDiffEqBDF/src/alg_utils.jl index aea50d0df0..8400ef7ea7 100644 --- a/lib/OrdinaryDiffEqBDF/src/alg_utils.jl +++ b/lib/OrdinaryDiffEqBDF/src/alg_utils.jl @@ -24,3 +24,16 @@ 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 + +isadaptive(alg::DImplicitEuler) = true +isadaptive(alg::DABDF2) = true +isadaptive(alg::DFBDF) = true + +alg_extrapolates(alg::DImplicitEuler) = true +alg_extrapolates(alg::DABDF2) = true + +alg_order(alg::DImplicitEuler) = 1 +alg_order(alg::DABDF2) = 2 +alg_order(alg::DFBDF) = 1#dummy value + +isfsal(alg::DImplicitEuler) = false \ No newline at end of file diff --git a/lib/OrdinaryDiffEqBDF/src/algorithms.jl b/lib/OrdinaryDiffEqBDF/src/algorithms.jl index 0a91dce0b6..0ed740a489 100644 --- a/lib/OrdinaryDiffEqBDF/src/algorithms.jl +++ b/lib/OrdinaryDiffEqBDF/src/algorithms.jl @@ -388,3 +388,79 @@ unew = uold + dt * (f1(unew) + f2(y1)) See also `SBDF`, `IMEXEuler`. """ IMEXEulerARK(; kwargs...) = SBDF(1; ark = true, kwargs...) + +struct DImplicitEuler{CS, AD, F, F2, P, FDT, ST, CJ} <: DAEAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + extrapolant::Symbol + controller::Symbol +end +function DImplicitEuler(; + 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, + controller = :Standard) + DImplicitEuler{_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 + +struct DABDF2{CS, AD, F, F2, P, FDT, ST, CJ} <: DAEAlgorithm{CS, AD, FDT, ST, CJ} + linsolve::F + nlsolve::F2 + precs::P + extrapolant::Symbol + controller::Symbol +end +function DABDF2(; 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 = :Standard) + DABDF2{_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 + +#= +struct DBDF{CS,AD,F,F2,P,FDT,ST,CJ} <: DAEAlgorithm{CS,AD,FDT,ST,CJ} + linsolve::F + nlsolve::F2 + precs::P + extrapolant::Symbol +end + +DBDF(;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) = + DBDF{_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) +=# + +struct DFBDF{MO, CS, AD, F, F2, P, FDT, ST, CJ, K, T} <: DAEAlgorithm{CS, AD, FDT, ST, CJ} + max_order::Val{MO} + linsolve::F + nlsolve::F2 + precs::P + κ::K + tol::T + extrapolant::Symbol + controller::Symbol +end +function DFBDF(; 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) where {MO} + DFBDF{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)}(max_order, linsolve, nlsolve, precs, κ, tol, extrapolant, + controller) +end + +TruncatedStacktraces.@truncate_stacktrace DFBDF \ No newline at end of file diff --git a/lib/OrdinaryDiffEqBDF/src/controllers.jl b/lib/OrdinaryDiffEqBDF/src/controllers.jl index 4c1e8b9274..23b9aab655 100644 --- a/lib/OrdinaryDiffEqBDF/src/controllers.jl +++ b/lib/OrdinaryDiffEqBDF/src/controllers.jl @@ -262,3 +262,135 @@ function step_accept_controller!(integrator, alg::FBDF{max_order}, integrator.cache.iters_from_event += 1 return integrator.dt / q end + + +function step_reject_controller!(integrator, ::DFBDF) + bdf_step_reject_controller!(integrator, integrator.cache.terkm1) +end + +function post_newton_controller!(integrator, alg) + integrator.dt = integrator.dt / integrator.opts.failfactor + nothing +end + +function post_newton_controller!(integrator, alg::DFBDF) + @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::DFBDF, 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::DFBDF, 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::DFBDF{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::DFBDF{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 \ No newline at end of file diff --git a/src/caches/dae_caches.jl b/lib/OrdinaryDiffEqBDF/src/dae_caches.jl similarity index 99% rename from src/caches/dae_caches.jl rename to lib/OrdinaryDiffEqBDF/src/dae_caches.jl index cc53b10a61..c432408952 100644 --- a/src/caches/dae_caches.jl +++ b/lib/OrdinaryDiffEqBDF/src/dae_caches.jl @@ -253,4 +253,4 @@ function alg_cache(alg::DFBDF{MO}, du, u, res_prototype, rate_prototype, uEltype u_corrector, u₀, bdf_coeffs, Val(5), nconsteps, consfailcnt, tmp, atmp, terkm2, terkm1, terk, terkp1, terk_tmp, terkp1_tmp, r, weights, equi_ts, iters_from_event) -end +end \ No newline at end of file diff --git a/src/perform_step/dae_perform_step.jl b/lib/OrdinaryDiffEqBDF/src/dae_perform_step.jl similarity index 99% rename from src/perform_step/dae_perform_step.jl rename to lib/OrdinaryDiffEqBDF/src/dae_perform_step.jl index a89ddc0001..4b4eab7e71 100644 --- a/src/perform_step/dae_perform_step.jl +++ b/lib/OrdinaryDiffEqBDF/src/dae_perform_step.jl @@ -457,4 +457,4 @@ function perform_step!(integrator, cache::DFBDFCache{max_order}, @.. broadcast=false integrator.fsallast=integrator.du = (nlsolver.α * z + nlsolver.tmp) * inv(nlsolver.γ * dt) #TODO Lorenz plot seems not smooth -end +end \ No newline at end of file diff --git a/test/interface/dae_ad_tests.jl b/lib/OrdinaryDiffEqBDF/test/dae_ad_tests.jl similarity index 99% rename from test/interface/dae_ad_tests.jl rename to lib/OrdinaryDiffEqBDF/test/dae_ad_tests.jl index 1040152adb..398aa310be 100644 --- a/test/interface/dae_ad_tests.jl +++ b/lib/OrdinaryDiffEqBDF/test/dae_ad_tests.jl @@ -57,4 +57,4 @@ sol1 = solve(prob, DFBDF(), dt = 1e-5, abstol = 1e-8, reltol = 1e-8) sum(sol) end @test ForwardDiff.gradient(f, [0.04, 3e7, 1e4])≈[0, 0, 0] atol=1e-8 -end +end \ No newline at end of file diff --git a/test/algconvergence/dae_convergence_tests.jl b/lib/OrdinaryDiffEqBDF/test/dae_convergence_tests.jl similarity index 99% rename from test/algconvergence/dae_convergence_tests.jl rename to lib/OrdinaryDiffEqBDF/test/dae_convergence_tests.jl index 859f194d68..e3ef2dd474 100644 --- a/test/algconvergence/dae_convergence_tests.jl +++ b/lib/OrdinaryDiffEqBDF/test/dae_convergence_tests.jl @@ -53,4 +53,4 @@ prob_dae_linear_oop = DAEProblem( @test sim24.𝒪est[:final]≈2 atol=testTol @test_nowarn solve(prob, DFBDF()) -end +end \ No newline at end of file diff --git a/test/integrators/dae_event.jl b/lib/OrdinaryDiffEqBDF/test/dae_event.jl similarity index 95% rename from test/integrators/dae_event.jl rename to lib/OrdinaryDiffEqBDF/test/dae_event.jl index 7c32b8d46e..c20b600670 100644 --- a/test/integrators/dae_event.jl +++ b/lib/OrdinaryDiffEqBDF/test/dae_event.jl @@ -27,4 +27,4 @@ sol = solve(prob, DFBDF(), callback = cb, tstops = [50.0], abstol = 1e-12, relto @test sol.t[end] == 100.0 @test sol[end][1]≈0.686300529575259 atol=1e-7 @test sol[end][2]≈2.0797982209353813e-6 atol=1e-7 -@test sol[end][3]≈1.31369739062652 atol=1e-7 +@test sol[end][3]≈1.31369739062652 atol=1e-7 \ No newline at end of file diff --git a/test/integrators/dae_initialization_tests.jl b/lib/OrdinaryDiffEqBDF/test/dae_initialization_tests.jl similarity index 99% rename from test/integrators/dae_initialization_tests.jl rename to lib/OrdinaryDiffEqBDF/test/dae_initialization_tests.jl index b3cdd2a858..5cf37c27f4 100644 --- a/test/integrators/dae_initialization_tests.jl +++ b/lib/OrdinaryDiffEqBDF/test/dae_initialization_tests.jl @@ -265,4 +265,4 @@ f = ODEFunction(f2, mass_matrix = Diagonal([1.0, 1.0, 0.0])) prob = ODEProblem(f, ones(3), (0.0, 1.0)) integrator = init(prob, Rodas5P(), initializealg = ShampineCollocationInit(1.0, BrokenNLSolve())) -@test all(isequal(reinterpret(Float64, 0xDEADBEEFDEADBEEF)), integrator.u) +@test all(isequal(reinterpret(Float64, 0xDEADBEEFDEADBEEF)), integrator.u) \ No newline at end of file diff --git a/test/interface/dae_initialize_integration.jl b/lib/OrdinaryDiffEqBDF/test/dae_initialize_integration.jl similarity index 98% rename from test/interface/dae_initialize_integration.jl rename to lib/OrdinaryDiffEqBDF/test/dae_initialize_integration.jl index 2d09cb1a41..33913925e2 100644 --- a/test/interface/dae_initialize_integration.jl +++ b/lib/OrdinaryDiffEqBDF/test/dae_initialize_integration.jl @@ -75,4 +75,4 @@ sol = solve(prob, Rodas5P(), dt = 1e-10) @test SciMLBase.successful_retcode(sol) @test sol[1] == [1.0] @test sol[2] ≈ [0.9999999998] -@test sol[end] ≈ [-1.0] +@test sol[end] ≈ [-1.0] \ No newline at end of file diff --git a/test/regression/hard_dae.jl b/lib/OrdinaryDiffEqBDF/test/hard_dae.jl similarity index 99% rename from test/regression/hard_dae.jl rename to lib/OrdinaryDiffEqBDF/test/hard_dae.jl index 302de87d2a..f47b226ce1 100644 --- a/test/regression/hard_dae.jl +++ b/lib/OrdinaryDiffEqBDF/test/hard_dae.jl @@ -272,4 +272,4 @@ for prob in [prob1, prob2], alg in [simple_implicit_euler, alg_switch] @test sol(0, idxs = 1) == 5 @test sol(2, idxs = 1) == 0 @test sol(4, idxs = 1) > 10 -end +end \ No newline at end of file diff --git a/lib/OrdinaryDiffEqBDF/test/runtest.jl b/lib/OrdinaryDiffEqBDF/test/runtest.jl new file mode 100644 index 0000000000..eaddde84ae --- /dev/null +++ b/lib/OrdinaryDiffEqBDF/test/runtest.jl @@ -0,0 +1 @@ +using SafeTestsets \ No newline at end of file diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index 6c66031e7b..57c9ee3aa1 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -20,6 +20,8 @@ using LinearAlgebra import StaticArrayInterface +using PrecompileTools + import InteractiveUtils using LinearSolve, SimpleNonlinearSolve @@ -155,7 +157,6 @@ include("caches/adams_bashforth_moulton_caches.jl") include("caches/nordsieck_caches.jl") include("caches/prk_caches.jl") include("caches/pdirk_caches.jl") -include("caches/dae_caches.jl") include("caches/qprk_caches.jl") include("tableaus/low_order_rk_tableaus.jl") @@ -165,7 +166,6 @@ include("tableaus/qprk_tableaus.jl") include("integrators/type.jl") include("integrators/controllers.jl") -include("integrators/integrator_utils.jl") include("integrators/integrator_interface.jl") include("cache_utils.jl") @@ -275,7 +275,8 @@ export ImplicitEuler, ImplicitMidpoint, Trapezoid, TRBDF2, SDIRK2, SDIRK22, include("../lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl") using ..OrdinaryDiffEqBDF export ABDF2, QNDF1, QBDF1, QNDF2, QBDF2, QNDF, QBDF, FBDF, - SBDF2, SBDF3, SBDF4, MEBDF2, IMEXEuler, IMEXEulerARK + SBDF2, SBDF3, SBDF4, MEBDF2, IMEXEuler, IMEXEulerARK, + DImplicitEuler, DABDF2, DFBDF include("../lib/OrdinaryDiffEqDefault/src/OrdinaryDiffEqDefault.jl") using ..OrdinaryDiffEqDefault @@ -290,11 +291,9 @@ using ..OrdinaryDiffEqBDF: reinitFBDF!, error_constant, estimate_terk!, 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 +using ..OrdinaryDiffEqBDF: post_newton_controller! +include("integrators/integrator_utils.jl") PrecompileTools.@compile_workload begin function lorenz(du, u, p, t) @@ -464,7 +463,7 @@ export Alshina2, Alshina3, Alshina6 export AutoSwitch, AutoTsit5, AutoDP5, AutoVern6, AutoVern7, AutoVern8, AutoVern9 -export KuttaPRK2p5, PDIRK44, DImplicitEuler, DABDF2, DFBDF +export KuttaPRK2p5, PDIRK44 export ShampineCollocationInit, BrownFullBasicInit, NoInit diff --git a/src/alg_utils.jl b/src/alg_utils.jl index 4120c7652b..80d21e2b8b 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -41,7 +41,6 @@ isfsal(alg::Rodas4P) = false isfsal(alg::Rodas4P2) = false # Pseudo Non-FSAL isfsal(alg::PDIRK44) = false -isfsal(alg::DImplicitEuler) = false isfsal(alg::RKO65) = false isfsal(alg::FRK65) = true #isfsal(alg::RKM) = false @@ -155,9 +154,6 @@ ismultistep(alg::ETD2) = true isadaptive(alg::Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}) = false isadaptive(alg::OrdinaryDiffEqAdaptiveAlgorithm) = true isadaptive(alg::OrdinaryDiffEqCompositeAlgorithm) = all(isadaptive.(alg.algs)) -isadaptive(alg::DImplicitEuler) = true -isadaptive(alg::DABDF2) = true -isadaptive(alg::DFBDF) = true anyadaptive(alg::Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}) = isadaptive(alg) anyadaptive(alg::OrdinaryDiffEqCompositeAlgorithm) = any(isadaptive, alg.algs) @@ -347,8 +343,6 @@ end alg_extrapolates(alg::Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}) = false alg_extrapolates(alg::CompositeAlgorithm) = any(alg_extrapolates.(alg.algs)) -alg_extrapolates(alg::DImplicitEuler) = true -alg_extrapolates(alg::DABDF2) = true alg_extrapolates(alg::MagnusLeapfrog) = true function alg_order(alg::Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}) @@ -507,10 +501,6 @@ alg_order(alg::JVODE) = 1 #dummy value alg_order(alg::PDIRK44) = 4 -alg_order(alg::DImplicitEuler) = 1 -alg_order(alg::DABDF2) = 2 -alg_order(alg::DFBDF) = 1#dummy value - alg_order(alg::Alshina2) = 2 alg_order(alg::Alshina3) = 3 alg_order(alg::Alshina6) = 6 diff --git a/src/algorithms.jl b/src/algorithms.jl index 238926aa4c..52985f65fb 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -931,80 +931,4 @@ end const MultistepAlgorithms = Union{ AB3, AB4, AB5, ABM32, ABM43, ABM54} -const SplitAlgorithms = Union{CNAB2, CNLF2} - -#= -struct DBDF{CS,AD,F,F2,P,FDT,ST,CJ} <: DAEAlgorithm{CS,AD,FDT,ST,CJ} - linsolve::F - nlsolve::F2 - precs::P - extrapolant::Symbol -end - -DBDF(;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) = - DBDF{_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) -=# - -struct DImplicitEuler{CS, AD, F, F2, P, FDT, ST, CJ} <: DAEAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - extrapolant::Symbol - controller::Symbol -end -function DImplicitEuler(; - 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, - controller = :Standard) - DImplicitEuler{_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 - -struct DABDF2{CS, AD, F, F2, P, FDT, ST, CJ} <: DAEAlgorithm{CS, AD, FDT, ST, CJ} - linsolve::F - nlsolve::F2 - precs::P - extrapolant::Symbol - controller::Symbol -end -function DABDF2(; 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 = :Standard) - DABDF2{_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 - -struct DFBDF{MO, CS, AD, F, F2, P, FDT, ST, CJ, K, T} <: DAEAlgorithm{CS, AD, FDT, ST, CJ} - max_order::Val{MO} - linsolve::F - nlsolve::F2 - precs::P - κ::K - tol::T - extrapolant::Symbol - controller::Symbol -end -function DFBDF(; 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) where {MO} - DFBDF{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)}(max_order, linsolve, nlsolve, precs, κ, tol, extrapolant, - controller) -end - -TruncatedStacktraces.@truncate_stacktrace DFBDF +const SplitAlgorithms = Union{CNAB2, CNLF2} \ No newline at end of file diff --git a/src/integrators/controllers.jl b/src/integrators/controllers.jl index 10af114a40..87b6d71fc1 100644 --- a/src/integrators/controllers.jl +++ b/src/integrators/controllers.jl @@ -422,135 +422,4 @@ end function step_reject_controller!(integrator, alg::JVODE) integrator.dt *= integrator.qold -end - -function step_reject_controller!(integrator, ::DFBDF) - bdf_step_reject_controller!(integrator, integrator.cache.terkm1) -end - -function post_newton_controller!(integrator, alg) - integrator.dt = integrator.dt / integrator.opts.failfactor - nothing -end - -function post_newton_controller!(integrator, alg::DFBDF) - @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::DFBDF, 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::DFBDF, 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::DFBDF{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::DFBDF{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 +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 3327354672..337ac4fa19 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -77,8 +77,14 @@ function activate_sdirk() Pkg.instantiate() end -function activate_firk() - Pkg.activate("../lib/OrdinaryDiffEqFIRK") +function activate_dae() + Pkg.activate("../lib/OrdinaryDiffEqDAE") + Pkg.develop(PackageSpec(path = dirname(@__DIR__))) + Pkg.instantiate() +end + +function activate_dae() + Pkg.activate("../lib/OrdinaryDiffEqDAE") Pkg.develop(PackageSpec(path = dirname(@__DIR__))) Pkg.instantiate() end @@ -147,9 +153,9 @@ end if !is_APPVEYOR && (GROUP == "All" || GROUP == "InterfaceV" || GROUP == "Interface") @time @safetestset "Interpolation Derivative Error Tests" include("interface/interpolation_derivative_error_tests.jl") @time @safetestset "AD Tests" include("interface/ad_tests.jl") - @time @safetestset "DAE AD Tests" include("interface/dae_ad_tests.jl") + @time @safetestset "DAE AD Tests" include("../lib/OrdinaryDiffEqBDF/test/dae_ad_tests.jl") @time @safetestset "Newton Tests" include("interface/newton_tests.jl") - @time @safetestset "DAE Initialize Integration" include("interface/dae_initialize_integration.jl") + @time @safetestset "DAE Initialize Integration" include("../lib/OrdinaryDiffEqBDF/test/dae_initialize_integration.jl") end if !is_APPVEYOR && @@ -173,8 +179,8 @@ end @time @safetestset "Reverse Directioned Event Tests" include("integrators/rev_events_tests.jl") @time @safetestset "Differentiation Direction Tests" include("integrators/diffdir_tests.jl") @time @safetestset "Resize Tests" include("integrators/resize_tests.jl") - @time @safetestset "DAE Initialization Tests" include("integrators/dae_initialization_tests.jl") - @time @safetestset "DAE Event Tests" include("integrators/dae_event.jl") + @time @safetestset "DAE Initialization Tests" include("../lib/OrdinaryDiffEqBDF/test/dae_initialization_tests.jl") + @time @safetestset "DAE Event Tests" include("../lib/OrdinaryDiffEqBDF/test/dae_event.jl") @time @safetestset "Cache Tests" include("integrators/ode_cache_tests.jl") @time @safetestset "Add Steps Tests" include("integrators/ode_add_steps_tests.jl") @time @safetestset "IMEX Split Function Tests" include("integrators/split_ode_tests.jl") @@ -185,7 +191,7 @@ end @time @safetestset "Special Interp Tests" include("regression/special_interps.jl") @time @safetestset "Inplace Tests" include("regression/ode_inplace_tests.jl") @time @safetestset "Adaptive Tests" include("regression/ode_adaptive_tests.jl") - @time @safetestset "Hard DAE Tests" include("regression/hard_dae.jl") + @time @safetestset "Hard DAE Tests" include("../lib/OrdinaryDiffEqBDF/test/hard_dae.jl") end if !is_APPVEYOR && (GROUP == "All" || GROUP == "Regression_II" || GROUP == "Regression") @@ -199,7 +205,7 @@ end if !is_APPVEYOR && GROUP == "AlgConvergence_I" @time @safetestset "Partitioned Methods Tests" include("algconvergence/partitioned_methods_tests.jl") @time @safetestset "Convergence Tests" include("algconvergence/ode_convergence_tests.jl") - @time @safetestset "DAE Convergence Tests" include("algconvergence/dae_convergence_tests.jl") + @time @safetestset "DAE Convergence Tests" include("../lib/OrdinaryDiffEqBDF/test/dae_convergence_tests.jl") @time @safetestset "Non-autonomous Convergence Tests" include("algconvergence/non-autonomous_convergence_tests.jl") @time @safetestset "Adams Variable Coefficients Tests" include("algconvergence/adams_tests.jl") @time @safetestset "Nordsieck Tests" include("algconvergence/nordsieck_tests.jl")