From b150b8874d90a6804c98e0232db2e3a958589edf Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Tue, 23 Jul 2024 02:03:40 +0530 Subject: [PATCH 01/11] Added DAE solvers --- lib/OrdinaryDiffEqDAE/Project.toml | 20 +++ .../src/OrdinaryDiffEqDAE.jl | 16 +++ lib/OrdinaryDiffEqDAE/src/alg_utils.jl | 12 ++ lib/OrdinaryDiffEqDAE/src/algorithms.jl | 75 ++++++++++ lib/OrdinaryDiffEqDAE/src/controllers.jl | 130 +++++++++++++++++ .../OrdinaryDiffEqDAE/src}/dae_caches.jl | 2 +- .../src}/dae_perform_step.jl | 2 +- src/OrdinaryDiffEq.jl | 6 +- src/alg_utils.jl | 10 -- src/algorithms.jl | 78 +--------- src/integrators/controllers.jl | 133 +----------------- 11 files changed, 261 insertions(+), 223 deletions(-) create mode 100644 lib/OrdinaryDiffEqDAE/Project.toml create mode 100644 lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl create mode 100644 lib/OrdinaryDiffEqDAE/src/alg_utils.jl create mode 100644 lib/OrdinaryDiffEqDAE/src/algorithms.jl create mode 100644 lib/OrdinaryDiffEqDAE/src/controllers.jl rename {src/caches => lib/OrdinaryDiffEqDAE/src}/dae_caches.jl (99%) rename {src/perform_step => lib/OrdinaryDiffEqDAE/src}/dae_perform_step.jl (99%) diff --git a/lib/OrdinaryDiffEqDAE/Project.toml b/lib/OrdinaryDiffEqDAE/Project.toml new file mode 100644 index 0000000000..ee3e8c74e7 --- /dev/null +++ b/lib/OrdinaryDiffEqDAE/Project.toml @@ -0,0 +1,20 @@ +name = "OrdinaryDiffEqDAE" +uuid = "7d2d7016-1c68-401d-834f-6be2c6dc7d88" +authors = ["ParamThakkar123 "] +version = "1.0.0" + +[deps] +OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +TruncatedStacktraces = "781d530d-4396-4725-bb49-402e4bee1e77" + +[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" + +[compat] +julia = "1.10" + +[targets] +test = ["DiffEqDevTools", "Random", "SafeTestsets", "Test"] \ No newline at end of file diff --git a/lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl b/lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl new file mode 100644 index 0000000000..0be59f3cbf --- /dev/null +++ b/lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl @@ -0,0 +1,16 @@ +module OrdinaryDiffEqDAE + +import OrdinaryDiffEq: _unwrap_val, NLNewton, DAEAlgorithm, + DEFAULT_PRECS, OrdinaryDiffEqConstantCache, + OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache +using TruncatedStacktraces + +include("algorithms.jl") +include("alg_utils.jl") +include("controllers.jl") +include("dae_caches.jl") +include("dae_perform_step.jl") + +export DABDF2, DImplicitEuler, DFBDF + +end # module OrdinaryDiffEqDAE diff --git a/lib/OrdinaryDiffEqDAE/src/alg_utils.jl b/lib/OrdinaryDiffEqDAE/src/alg_utils.jl new file mode 100644 index 0000000000..ff0d60daac --- /dev/null +++ b/lib/OrdinaryDiffEqDAE/src/alg_utils.jl @@ -0,0 +1,12 @@ +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/OrdinaryDiffEqDAE/src/algorithms.jl b/lib/OrdinaryDiffEqDAE/src/algorithms.jl new file mode 100644 index 0000000000..f8584348a5 --- /dev/null +++ b/lib/OrdinaryDiffEqDAE/src/algorithms.jl @@ -0,0 +1,75 @@ +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/OrdinaryDiffEqDAE/src/controllers.jl b/lib/OrdinaryDiffEqDAE/src/controllers.jl new file mode 100644 index 0000000000..4ae1a93a28 --- /dev/null +++ b/lib/OrdinaryDiffEqDAE/src/controllers.jl @@ -0,0 +1,130 @@ +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/OrdinaryDiffEqDAE/src/dae_caches.jl similarity index 99% rename from src/caches/dae_caches.jl rename to lib/OrdinaryDiffEqDAE/src/dae_caches.jl index cc53b10a61..c432408952 100644 --- a/src/caches/dae_caches.jl +++ b/lib/OrdinaryDiffEqDAE/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/OrdinaryDiffEqDAE/src/dae_perform_step.jl similarity index 99% rename from src/perform_step/dae_perform_step.jl rename to lib/OrdinaryDiffEqDAE/src/dae_perform_step.jl index a89ddc0001..4b4eab7e71 100644 --- a/src/perform_step/dae_perform_step.jl +++ b/lib/OrdinaryDiffEqDAE/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/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index 6214e87808..4ba78958c1 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -156,7 +156,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") @@ -289,8 +288,11 @@ using ..OrdinaryDiffEqBDF: reinitFBDF!, error_constant, estimate_terk!, calc_finite_difference_weights, estimate_terk, calc_Lagrange_interp, bdf_step_reject_controller! + +include("../lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl") +using ..OrdinaryDiffEqDAE +export DImplicitEuler, DABDF2, DFBDF include("nlsolve/newton.jl") -include("perform_step/dae_perform_step.jl") import PrecompileTools diff --git a/src/alg_utils.jl b/src/alg_utils.jl index 77bf1df0f4..34e9722c1d 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 @@ -153,9 +152,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) @@ -346,8 +342,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}) @@ -509,10 +503,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 1b921906bb..1ec6d9e17b 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -1093,80 +1093,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 69169566c3..e4cab76b04 100644 --- a/src/integrators/controllers.jl +++ b/src/integrators/controllers.jl @@ -473,135 +473,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 From 09f9af8ce03445da96734237d70a1403a94086a3 Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Tue, 23 Jul 2024 02:17:41 +0530 Subject: [PATCH 02/11] @unpack --- lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl | 3 ++- .../OrdinaryDiffEqDAE/test}/dae_ad_tests.jl | 2 +- .../test}/dae_convergence_tests.jl | 2 +- .../OrdinaryDiffEqDAE/test}/dae_event.jl | 2 +- .../test}/dae_initialization_tests.jl | 2 +- .../test}/dae_initialize_integration.jl | 2 +- .../OrdinaryDiffEqDAE/test}/hard_dae.jl | 2 +- lib/OrdinaryDiffEqDAE/test/runtest.jl | 0 src/OrdinaryDiffEq.jl | 2 +- test/runtests.jl | 18 ++++++++++++------ 10 files changed, 21 insertions(+), 14 deletions(-) rename {test/interface => lib/OrdinaryDiffEqDAE/test}/dae_ad_tests.jl (99%) rename {test/algconvergence => lib/OrdinaryDiffEqDAE/test}/dae_convergence_tests.jl (99%) rename {test/integrators => lib/OrdinaryDiffEqDAE/test}/dae_event.jl (95%) rename {test/integrators => lib/OrdinaryDiffEqDAE/test}/dae_initialization_tests.jl (99%) rename {test/interface => lib/OrdinaryDiffEqDAE/test}/dae_initialize_integration.jl (98%) rename {test/regression => lib/OrdinaryDiffEqDAE/test}/hard_dae.jl (99%) create mode 100644 lib/OrdinaryDiffEqDAE/test/runtest.jl diff --git a/lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl b/lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl index 0be59f3cbf..ed77fee311 100644 --- a/lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl +++ b/lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl @@ -2,7 +2,8 @@ module OrdinaryDiffEqDAE import OrdinaryDiffEq: _unwrap_val, NLNewton, DAEAlgorithm, DEFAULT_PRECS, OrdinaryDiffEqConstantCache, - OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache + OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, + @unpack using TruncatedStacktraces include("algorithms.jl") diff --git a/test/interface/dae_ad_tests.jl b/lib/OrdinaryDiffEqDAE/test/dae_ad_tests.jl similarity index 99% rename from test/interface/dae_ad_tests.jl rename to lib/OrdinaryDiffEqDAE/test/dae_ad_tests.jl index 1040152adb..398aa310be 100644 --- a/test/interface/dae_ad_tests.jl +++ b/lib/OrdinaryDiffEqDAE/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/OrdinaryDiffEqDAE/test/dae_convergence_tests.jl similarity index 99% rename from test/algconvergence/dae_convergence_tests.jl rename to lib/OrdinaryDiffEqDAE/test/dae_convergence_tests.jl index 859f194d68..e3ef2dd474 100644 --- a/test/algconvergence/dae_convergence_tests.jl +++ b/lib/OrdinaryDiffEqDAE/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/OrdinaryDiffEqDAE/test/dae_event.jl similarity index 95% rename from test/integrators/dae_event.jl rename to lib/OrdinaryDiffEqDAE/test/dae_event.jl index 7c32b8d46e..c20b600670 100644 --- a/test/integrators/dae_event.jl +++ b/lib/OrdinaryDiffEqDAE/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/OrdinaryDiffEqDAE/test/dae_initialization_tests.jl similarity index 99% rename from test/integrators/dae_initialization_tests.jl rename to lib/OrdinaryDiffEqDAE/test/dae_initialization_tests.jl index b3cdd2a858..5cf37c27f4 100644 --- a/test/integrators/dae_initialization_tests.jl +++ b/lib/OrdinaryDiffEqDAE/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/OrdinaryDiffEqDAE/test/dae_initialize_integration.jl similarity index 98% rename from test/interface/dae_initialize_integration.jl rename to lib/OrdinaryDiffEqDAE/test/dae_initialize_integration.jl index 2d09cb1a41..33913925e2 100644 --- a/test/interface/dae_initialize_integration.jl +++ b/lib/OrdinaryDiffEqDAE/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/OrdinaryDiffEqDAE/test/hard_dae.jl similarity index 99% rename from test/regression/hard_dae.jl rename to lib/OrdinaryDiffEqDAE/test/hard_dae.jl index 302de87d2a..f47b226ce1 100644 --- a/test/regression/hard_dae.jl +++ b/lib/OrdinaryDiffEqDAE/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/OrdinaryDiffEqDAE/test/runtest.jl b/lib/OrdinaryDiffEqDAE/test/runtest.jl new file mode 100644 index 0000000000..e69de29bb2 diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index 4ba78958c1..cc134875bd 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -466,7 +466,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/test/runtests.jl b/test/runtests.jl index 2bce643dd9..5779113d6b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -77,6 +77,12 @@ function activate_sdirk() Pkg.instantiate() end +function activate_dae() + Pkg.activate("../lib/OrdinaryDiffEqDAE") + Pkg.develop(PackageSpec(path = dirname(@__DIR__))) + Pkg.instantiate() +end + #Start Test Script @time begin @@ -141,9 +147,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/OrdinaryDiffEqDAE/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/OrdinaryDiffEqDAE/test/dae_initialize_integration.jl") end if !is_APPVEYOR && @@ -167,8 +173,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/OrdinaryDiffEqDAE/test/dae_initialization_tests.jl") + @time @safetestset "DAE Event Tests" include("../lib/OrdinaryDiffEqDAE/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") @@ -179,7 +185,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/OrdinaryDiffEqDAE/test/hard_dae.jl") end if !is_APPVEYOR && (GROUP == "All" || GROUP == "Regression_II" || GROUP == "Regression") @@ -193,7 +199,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/OrdinaryDiffEqDAE/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") From dbb4e5fa28cc685df0100047ef0e55fb518c3de5 Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Tue, 23 Jul 2024 10:19:15 +0530 Subject: [PATCH 03/11] @.. --- lib/OrdinaryDiffEqDAE/Project.toml | 1 + lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl | 1 + 2 files changed, 2 insertions(+) diff --git a/lib/OrdinaryDiffEqDAE/Project.toml b/lib/OrdinaryDiffEqDAE/Project.toml index ee3e8c74e7..b15e3bbe39 100644 --- a/lib/OrdinaryDiffEqDAE/Project.toml +++ b/lib/OrdinaryDiffEqDAE/Project.toml @@ -6,6 +6,7 @@ version = "1.0.0" [deps] OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" TruncatedStacktraces = "781d530d-4396-4725-bb49-402e4bee1e77" +FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898" [extras] DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d" diff --git a/lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl b/lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl index ed77fee311..82a1d7d496 100644 --- a/lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl +++ b/lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl @@ -5,6 +5,7 @@ import OrdinaryDiffEq: _unwrap_val, NLNewton, DAEAlgorithm, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, @unpack using TruncatedStacktraces +import FastBroadCast: @.. include("algorithms.jl") include("alg_utils.jl") From 70e3240518eab50f367d01e08c0a6177ea0324f1 Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Tue, 23 Jul 2024 10:41:13 +0530 Subject: [PATCH 04/11] @.. --- lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl | 2 +- lib/OrdinaryDiffEqDAE/test/runtest.jl | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl b/lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl index 82a1d7d496..9643b52c2d 100644 --- a/lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl +++ b/lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl @@ -5,7 +5,7 @@ import OrdinaryDiffEq: _unwrap_val, NLNewton, DAEAlgorithm, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, @unpack using TruncatedStacktraces -import FastBroadCast: @.. +import FastBroadcast: @.. include("algorithms.jl") include("alg_utils.jl") diff --git a/lib/OrdinaryDiffEqDAE/test/runtest.jl b/lib/OrdinaryDiffEqDAE/test/runtest.jl index e69de29bb2..a639da4027 100644 --- a/lib/OrdinaryDiffEqDAE/test/runtest.jl +++ b/lib/OrdinaryDiffEqDAE/test/runtest.jl @@ -0,0 +1,2 @@ +using SafeTestsets + From c40127efaf84f06746bea3cdfa227129cbe5e90f Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Tue, 23 Jul 2024 10:43:51 +0530 Subject: [PATCH 05/11] @cache --- lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl b/lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl index 9643b52c2d..374e6b772b 100644 --- a/lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl +++ b/lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl @@ -3,7 +3,7 @@ module OrdinaryDiffEqDAE import OrdinaryDiffEq: _unwrap_val, NLNewton, DAEAlgorithm, DEFAULT_PRECS, OrdinaryDiffEqConstantCache, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, - @unpack + @unpack, @cache using TruncatedStacktraces import FastBroadcast: @.. From d9769abf17aedbfaa9b666163ee9911034c6b6b1 Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Tue, 23 Jul 2024 10:46:02 +0530 Subject: [PATCH 06/11] MuladdMacro --- lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl b/lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl index 374e6b772b..0657ec30b9 100644 --- a/lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl +++ b/lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl @@ -4,7 +4,7 @@ import OrdinaryDiffEq: _unwrap_val, NLNewton, DAEAlgorithm, DEFAULT_PRECS, OrdinaryDiffEqConstantCache, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, @unpack, @cache -using TruncatedStacktraces +using TruncatedStacktraces, MuladdMacro import FastBroadcast: @.. include("algorithms.jl") From 44e3a2a6e07e645624b502bcf0810be15ca039f5 Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Tue, 23 Jul 2024 10:56:06 +0530 Subject: [PATCH 07/11] controllers --- lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl b/lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl index 0657ec30b9..9017917bb8 100644 --- a/lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl +++ b/lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl @@ -3,7 +3,9 @@ module OrdinaryDiffEqDAE import OrdinaryDiffEq: _unwrap_val, NLNewton, DAEAlgorithm, DEFAULT_PRECS, OrdinaryDiffEqConstantCache, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, - @unpack, @cache + @unpack, @cache, default_controller, stepsize_controller!, + step_accept_controller!, + step_reject_controller!, post_newton_controller! using TruncatedStacktraces, MuladdMacro import FastBroadcast: @.. From 79a15979908064270e55c7bf1822d1df2d6ffb89 Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Tue, 23 Jul 2024 11:13:21 +0530 Subject: [PATCH 08/11] post_newton_controller! --- src/OrdinaryDiffEq.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index cc134875bd..18d1dcc5dc 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -166,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") @@ -292,6 +291,10 @@ using ..OrdinaryDiffEqBDF: reinitFBDF!, error_constant, estimate_terk!, include("../lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl") using ..OrdinaryDiffEqDAE export DImplicitEuler, DABDF2, DFBDF + + +using ..OrdinaryDiffEqDAE: post_newton_controller! +include("integrators/integrator_utils.jl") include("nlsolve/newton.jl") import PrecompileTools From 24d3c61a4b8897ea0cc2206a9aa0353838ebca57 Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Fri, 26 Jul 2024 22:20:47 +0530 Subject: [PATCH 09/11] DAE added to BDF --- .../src/OrdinaryDiffEqBDF.jl | 7 +- lib/OrdinaryDiffEqBDF/src/alg_utils.jl | 13 ++ lib/OrdinaryDiffEqBDF/src/algorithms.jl | 76 ++++++++++ lib/OrdinaryDiffEqBDF/src/controllers.jl | 132 ++++++++++++++++++ .../src/dae_caches.jl | 0 .../src/dae_perform_step.jl | 0 .../test/dae_ad_tests.jl | 0 .../test/dae_convergence_tests.jl | 0 .../test/dae_event.jl | 0 .../test/dae_initialization_tests.jl | 0 .../test/dae_initialize_integration.jl | 0 .../test/hard_dae.jl | 0 .../test/runtest.jl | 0 lib/OrdinaryDiffEqDAE/Project.toml | 21 --- .../src/OrdinaryDiffEqDAE.jl | 20 --- lib/OrdinaryDiffEqDAE/src/alg_utils.jl | 12 -- lib/OrdinaryDiffEqDAE/src/algorithms.jl | 75 ---------- lib/OrdinaryDiffEqDAE/src/controllers.jl | 130 ----------------- src/OrdinaryDiffEq.jl | 10 +- test/runtests.jl | 12 +- 20 files changed, 235 insertions(+), 273 deletions(-) rename lib/{OrdinaryDiffEqDAE => OrdinaryDiffEqBDF}/src/dae_caches.jl (100%) rename lib/{OrdinaryDiffEqDAE => OrdinaryDiffEqBDF}/src/dae_perform_step.jl (100%) rename lib/{OrdinaryDiffEqDAE => OrdinaryDiffEqBDF}/test/dae_ad_tests.jl (100%) rename lib/{OrdinaryDiffEqDAE => OrdinaryDiffEqBDF}/test/dae_convergence_tests.jl (100%) rename lib/{OrdinaryDiffEqDAE => OrdinaryDiffEqBDF}/test/dae_event.jl (100%) rename lib/{OrdinaryDiffEqDAE => OrdinaryDiffEqBDF}/test/dae_initialization_tests.jl (100%) rename lib/{OrdinaryDiffEqDAE => OrdinaryDiffEqBDF}/test/dae_initialize_integration.jl (100%) rename lib/{OrdinaryDiffEqDAE => OrdinaryDiffEqBDF}/test/hard_dae.jl (100%) rename lib/{OrdinaryDiffEqDAE => OrdinaryDiffEqBDF}/test/runtest.jl (100%) delete mode 100644 lib/OrdinaryDiffEqDAE/Project.toml delete mode 100644 lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl delete mode 100644 lib/OrdinaryDiffEqDAE/src/alg_utils.jl delete mode 100644 lib/OrdinaryDiffEqDAE/src/algorithms.jl delete mode 100644 lib/OrdinaryDiffEqDAE/src/controllers.jl diff --git a/lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl b/lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl index 2080050161..740ced5261 100644 --- a/lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl +++ b/lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl @@ -22,7 +22,7 @@ import OrdinaryDiffEq: alg_order, calculate_residuals!, 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 +32,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/lib/OrdinaryDiffEqDAE/src/dae_caches.jl b/lib/OrdinaryDiffEqBDF/src/dae_caches.jl similarity index 100% rename from lib/OrdinaryDiffEqDAE/src/dae_caches.jl rename to lib/OrdinaryDiffEqBDF/src/dae_caches.jl diff --git a/lib/OrdinaryDiffEqDAE/src/dae_perform_step.jl b/lib/OrdinaryDiffEqBDF/src/dae_perform_step.jl similarity index 100% rename from lib/OrdinaryDiffEqDAE/src/dae_perform_step.jl rename to lib/OrdinaryDiffEqBDF/src/dae_perform_step.jl diff --git a/lib/OrdinaryDiffEqDAE/test/dae_ad_tests.jl b/lib/OrdinaryDiffEqBDF/test/dae_ad_tests.jl similarity index 100% rename from lib/OrdinaryDiffEqDAE/test/dae_ad_tests.jl rename to lib/OrdinaryDiffEqBDF/test/dae_ad_tests.jl diff --git a/lib/OrdinaryDiffEqDAE/test/dae_convergence_tests.jl b/lib/OrdinaryDiffEqBDF/test/dae_convergence_tests.jl similarity index 100% rename from lib/OrdinaryDiffEqDAE/test/dae_convergence_tests.jl rename to lib/OrdinaryDiffEqBDF/test/dae_convergence_tests.jl diff --git a/lib/OrdinaryDiffEqDAE/test/dae_event.jl b/lib/OrdinaryDiffEqBDF/test/dae_event.jl similarity index 100% rename from lib/OrdinaryDiffEqDAE/test/dae_event.jl rename to lib/OrdinaryDiffEqBDF/test/dae_event.jl diff --git a/lib/OrdinaryDiffEqDAE/test/dae_initialization_tests.jl b/lib/OrdinaryDiffEqBDF/test/dae_initialization_tests.jl similarity index 100% rename from lib/OrdinaryDiffEqDAE/test/dae_initialization_tests.jl rename to lib/OrdinaryDiffEqBDF/test/dae_initialization_tests.jl diff --git a/lib/OrdinaryDiffEqDAE/test/dae_initialize_integration.jl b/lib/OrdinaryDiffEqBDF/test/dae_initialize_integration.jl similarity index 100% rename from lib/OrdinaryDiffEqDAE/test/dae_initialize_integration.jl rename to lib/OrdinaryDiffEqBDF/test/dae_initialize_integration.jl diff --git a/lib/OrdinaryDiffEqDAE/test/hard_dae.jl b/lib/OrdinaryDiffEqBDF/test/hard_dae.jl similarity index 100% rename from lib/OrdinaryDiffEqDAE/test/hard_dae.jl rename to lib/OrdinaryDiffEqBDF/test/hard_dae.jl diff --git a/lib/OrdinaryDiffEqDAE/test/runtest.jl b/lib/OrdinaryDiffEqBDF/test/runtest.jl similarity index 100% rename from lib/OrdinaryDiffEqDAE/test/runtest.jl rename to lib/OrdinaryDiffEqBDF/test/runtest.jl diff --git a/lib/OrdinaryDiffEqDAE/Project.toml b/lib/OrdinaryDiffEqDAE/Project.toml deleted file mode 100644 index b15e3bbe39..0000000000 --- a/lib/OrdinaryDiffEqDAE/Project.toml +++ /dev/null @@ -1,21 +0,0 @@ -name = "OrdinaryDiffEqDAE" -uuid = "7d2d7016-1c68-401d-834f-6be2c6dc7d88" -authors = ["ParamThakkar123 "] -version = "1.0.0" - -[deps] -OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" -TruncatedStacktraces = "781d530d-4396-4725-bb49-402e4bee1e77" -FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898" - -[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" - -[compat] -julia = "1.10" - -[targets] -test = ["DiffEqDevTools", "Random", "SafeTestsets", "Test"] \ No newline at end of file diff --git a/lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl b/lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl deleted file mode 100644 index 9017917bb8..0000000000 --- a/lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl +++ /dev/null @@ -1,20 +0,0 @@ -module OrdinaryDiffEqDAE - -import OrdinaryDiffEq: _unwrap_val, NLNewton, DAEAlgorithm, - DEFAULT_PRECS, OrdinaryDiffEqConstantCache, - OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, - @unpack, @cache, default_controller, stepsize_controller!, - step_accept_controller!, - step_reject_controller!, post_newton_controller! -using TruncatedStacktraces, MuladdMacro -import FastBroadcast: @.. - -include("algorithms.jl") -include("alg_utils.jl") -include("controllers.jl") -include("dae_caches.jl") -include("dae_perform_step.jl") - -export DABDF2, DImplicitEuler, DFBDF - -end # module OrdinaryDiffEqDAE diff --git a/lib/OrdinaryDiffEqDAE/src/alg_utils.jl b/lib/OrdinaryDiffEqDAE/src/alg_utils.jl deleted file mode 100644 index ff0d60daac..0000000000 --- a/lib/OrdinaryDiffEqDAE/src/alg_utils.jl +++ /dev/null @@ -1,12 +0,0 @@ -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/OrdinaryDiffEqDAE/src/algorithms.jl b/lib/OrdinaryDiffEqDAE/src/algorithms.jl deleted file mode 100644 index f8584348a5..0000000000 --- a/lib/OrdinaryDiffEqDAE/src/algorithms.jl +++ /dev/null @@ -1,75 +0,0 @@ -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/OrdinaryDiffEqDAE/src/controllers.jl b/lib/OrdinaryDiffEqDAE/src/controllers.jl deleted file mode 100644 index 4ae1a93a28..0000000000 --- a/lib/OrdinaryDiffEqDAE/src/controllers.jl +++ /dev/null @@ -1,130 +0,0 @@ -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/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index 69c6504a62..68cc4991e1 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -273,7 +273,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 @@ -289,12 +290,7 @@ using ..OrdinaryDiffEqBDF: reinitFBDF!, error_constant, estimate_terk!, calc_Lagrange_interp, bdf_step_reject_controller! -include("../lib/OrdinaryDiffEqDAE/src/OrdinaryDiffEqDAE.jl") -using ..OrdinaryDiffEqDAE -export DImplicitEuler, DABDF2, DFBDF - - -using ..OrdinaryDiffEqDAE: post_newton_controller! +using ..OrdinaryDiffEqBDF: post_newton_controller! include("integrators/integrator_utils.jl") include("nlsolve/newton.jl") diff --git a/test/runtests.jl b/test/runtests.jl index fd9dbd5378..d3f88cdec2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -153,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("../lib/OrdinaryDiffEqDAE/test/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("../lib/OrdinaryDiffEqDAE/test/dae_initialize_integration.jl") + @time @safetestset "DAE Initialize Integration" include("../lib/OrdinaryDiffEqBDF/test/dae_initialize_integration.jl") end if !is_APPVEYOR && @@ -179,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("../lib/OrdinaryDiffEqDAE/test/dae_initialization_tests.jl") - @time @safetestset "DAE Event Tests" include("../lib/OrdinaryDiffEqDAE/test/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") @@ -191,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("../lib/OrdinaryDiffEqDAE/test/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") @@ -205,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("../lib/OrdinaryDiffEqDAE/test/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") From cc67c34b1dabf5769a41c7fcd979f4288f0dba00 Mon Sep 17 00:00:00 2001 From: ParamThakkar123 Date: Fri, 26 Jul 2024 22:31:11 +0530 Subject: [PATCH 10/11] added tests --- lib/OrdinaryDiffEqBDF/test/runtest.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/lib/OrdinaryDiffEqBDF/test/runtest.jl b/lib/OrdinaryDiffEqBDF/test/runtest.jl index a639da4027..aeacfd70d2 100644 --- a/lib/OrdinaryDiffEqBDF/test/runtest.jl +++ b/lib/OrdinaryDiffEqBDF/test/runtest.jl @@ -1,2 +1,8 @@ using SafeTestsets +@time @safetestset "DAE AD Tests" include("dae_ad_tests.jl") +@time @safetestset "DAE Initialize Integration" include("dae_initialize_integration.jl") +@time @safetestset "DAE Initialization Tests" include("dae_initialization_tests.jl") +@time @safetestset "DAE Event Tests" include("dae_event.jl") +@time @safetestset "Hard DAE Tests" include("hard_dae.jl") +@time @safetestset "DAE Convergence Tests" include("dae_convergence_tests.jl") \ No newline at end of file From 1f834d6fb40fbbde3d4eb410403ef96f0b14b806 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Fri, 26 Jul 2024 23:33:54 -0400 Subject: [PATCH 11/11] Update lib/OrdinaryDiffEqBDF/src/algorithms.jl --- lib/OrdinaryDiffEqBDF/src/algorithms.jl | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/lib/OrdinaryDiffEqBDF/src/algorithms.jl b/lib/OrdinaryDiffEqBDF/src/algorithms.jl index 0ed740a489..71900e4925 100644 --- a/lib/OrdinaryDiffEqBDF/src/algorithms.jl +++ b/lib/OrdinaryDiffEqBDF/src/algorithms.jl @@ -426,19 +426,6 @@ function DABDF2(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = V 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}