From 2ad1812716d36fb9347f62aa4048356e51a6f00c Mon Sep 17 00:00:00 2001 From: Gaurav Arya Date: Sun, 26 Mar 2023 18:05:54 -0400 Subject: [PATCH 01/54] Support SparseDiffTools v2 --- Project.toml | 2 +- src/OrdinaryDiffEq.jl | 2 ++ src/derivative_utils.jl | 11 ++--------- src/integrators/integrator_interface.jl | 13 ++++--------- src/misc_utils.jl | 2 +- 5 files changed, 10 insertions(+), 20 deletions(-) diff --git a/Project.toml b/Project.toml index b330b06b3e..62947fffe1 100644 --- a/Project.toml +++ b/Project.toml @@ -69,7 +69,7 @@ SciMLNLSolve = "0.1" SimpleNonlinearSolve = "0.1.4" SimpleUnPack = "1" SnoopPrecompile = "1" -SparseDiffTools = "1.26.2" +SparseDiffTools = "1.26.2,2" StaticArrayInterface = "1.2" StaticArrays = "0.11, 0.12, 1.0" TruncatedStacktraces = "1.2" diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index 5074891f58..87af471e9e 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -38,6 +38,8 @@ using DiffEqBase: TimeGradientWrapper, UJacobianWrapper, TimeDerivativeWrapper, using DiffEqBase: DEIntegrator +import SciMLBase: update_coefficients! + import RecursiveArrayTools: chain, recursivecopy! using SimpleUnPack, ForwardDiff, RecursiveArrayTools, diff --git a/src/derivative_utils.jl b/src/derivative_utils.jl index daedbf0d11..e81c67ddb8 100644 --- a/src/derivative_utils.jl +++ b/src/derivative_utils.jl @@ -294,19 +294,13 @@ SciMLBase.isinplace(::WOperator{IIP}, i) where {IIP} = IIP Base.eltype(W::WOperator) = eltype(W.J) set_gamma!(W::WOperator, gamma) = (W.gamma = gamma; W) -function DiffEqBase.update_coefficients!(W::WOperator, u, p, t) +function update_coefficients!(W::WOperator, u, p, t) update_coefficients!(W.J, u, p, t) update_coefficients!(W.mass_matrix, u, p, t) W.jacvec !== nothing && update_coefficients!(W.jacvec, u, p, t) W end -function DiffEqBase.update_coefficients!(J::SparseDiffTools.JacVec, u, p, t) - copyto!(J.x, u) - J.f.t = t - J.f.p = p -end - function Base.convert(::Type{AbstractMatrix}, W::WOperator{IIP}) where {IIP} if !IIP # Allocating @@ -678,8 +672,7 @@ function calc_W!(W, integrator, nlsolver::Union{Nothing, AbstractNLSolver}, cach isnewton(nlsolver) || DiffEqBase.update_coefficients!(W, uprev, p, t) # we will call `update_coefficients!` in NLNewton W.transform = W_transform set_gamma!(W, dtgamma) - if W.J !== nothing && !(W.J isa SparseDiffTools.JacVec) && - !(W.J isa AbstractSciMLOperator) + if W.J !== nothing && !(W.J isa AbstractSciMLOperator) islin, isode = islinearfunction(integrator) islin ? (J = isode ? f.f : f.f1.f) : (new_jac && (calc_J!(W.J, integrator, lcache, next_step))) diff --git a/src/integrators/integrator_interface.jl b/src/integrators/integrator_interface.jl index d452119d5d..92cb6df745 100644 --- a/src/integrators/integrator_interface.jl +++ b/src/integrators/integrator_interface.jl @@ -234,20 +234,15 @@ function resize_J_W!(cache, integrator, i) nf = nlsolve_f(f, integrator.alg) islin = f isa Union{ODEFunction, SplitFunction} && islinear(nf.f) if !islin - if isa(cache.J, AbstractSciMLOperator) + if cache.J isa AbstractSciMLOperator resize!(cache.J, i) elseif f.jac_prototype !== nothing J = similar(f.jac_prototype, i, i) J = DiffEqArrayOperator(J; update_func = f.jac) - elseif cache.J isa SparseDiffTools.JacVec - resize!(cache.J.cache1, i) - resize!(cache.J.cache2, i) - resize!(cache.J.x, i) end - if cache.W.jacvec !== nothing - resize!(cache.W.jacvec.cache1, i) - resize!(cache.W.jacvec.cache2, i) - resize!(cache.W.jacvec.x, i) + if cache.W.jacvec isa AbstractSciMLOperator + # TODO: resize! will need to be implemented upstream to handle previously handled cases + resize!(cache.W.jacvec, i) end cache.W = WOperator{DiffEqBase.isinplace(integrator.sol.prob)}(f.mass_matrix, integrator.dt, diff --git a/src/misc_utils.jl b/src/misc_utils.jl index c83d20d1b8..7045e86269 100644 --- a/src/misc_utils.jl +++ b/src/misc_utils.jl @@ -108,7 +108,7 @@ function dolinsolve(integrator, linsolve; A = nothing, linu = nothing, b = nothi # TODO: this ignores the add of the `f` count for add_steps! if integrator isa SciMLBase.DEIntegrator && _alg.linsolve !== nothing && !LinearSolve.needs_concrete_A(_alg.linsolve) && - linsolve.A isa WOperator && linsolve.A.J isa SparseDiffTools.JacVec + linsolve.A isa WOperator && linsolve.A.J isa AbstractSciMLOperator if alg_autodiff(_alg) integrator.stats.nf += linres.iters else From ea9ef6599d066d7ba3c4c80db904ea9dd90e3eb8 Mon Sep 17 00:00:00 2001 From: Gaurav Arya Date: Sun, 26 Mar 2023 18:08:07 -0400 Subject: [PATCH 02/54] Tweak version format --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 62947fffe1..1dde0dbd7c 100644 --- a/Project.toml +++ b/Project.toml @@ -69,7 +69,7 @@ SciMLNLSolve = "0.1" SimpleNonlinearSolve = "0.1.4" SimpleUnPack = "1" SnoopPrecompile = "1" -SparseDiffTools = "1.26.2,2" +SparseDiffTools = "1.26.2, 2" StaticArrayInterface = "1.2" StaticArrays = "0.11, 0.12, 1.0" TruncatedStacktraces = "1.2" From 2ef6999829fe2cb374f378fd58a5c54d830d8bd5 Mon Sep 17 00:00:00 2001 From: Gaurav Arya Date: Sun, 26 Mar 2023 18:10:19 -0400 Subject: [PATCH 03/54] Drop support for SparseDiffTools v1 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 1dde0dbd7c..7b3ce829ab 100644 --- a/Project.toml +++ b/Project.toml @@ -69,7 +69,7 @@ SciMLNLSolve = "0.1" SimpleNonlinearSolve = "0.1.4" SimpleUnPack = "1" SnoopPrecompile = "1" -SparseDiffTools = "1.26.2, 2" +SparseDiffTools = "2" StaticArrayInterface = "1.2" StaticArrays = "0.11, 0.12, 1.0" TruncatedStacktraces = "1.2" From 9c8c24bfc348e3200abfb276c3b6d20144bddabc Mon Sep 17 00:00:00 2001 From: Gaurav Arya Date: Tue, 28 Mar 2023 00:13:49 -0400 Subject: [PATCH 04/54] Update autodiff selection --- Project.toml | 1 + src/OrdinaryDiffEq.jl | 1 + src/alg_utils.jl | 26 ++++++++++++++++++++------ 3 files changed, 22 insertions(+), 6 deletions(-) diff --git a/Project.toml b/Project.toml index 7b3ce829ab..e184af8d8f 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["Chris Rackauckas ", "Yingbo Ma Date: Tue, 28 Mar 2023 00:16:06 -0400 Subject: [PATCH 05/54] Fix JacVec calls --- src/derivative_utils.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/derivative_utils.jl b/src/derivative_utils.jl index e81c67ddb8..cde69376e8 100644 --- a/src/derivative_utils.jl +++ b/src/derivative_utils.jl @@ -847,7 +847,7 @@ function build_J_W(alg, u, uprev, p, t, dt, f::F, ::Type{uEltypeNoUnits}, _f = islin ? (isode ? f.f : f.f1.f) : f jacvec = SparseDiffTools.JacVec(UJacobianWrapper(_f, t, p), copy(u), - OrdinaryDiffEqTag(), autodiff = alg_autodiff(alg)) + p, t; autodiff = alg_autodiff(alg)) J = jacvec W = WOperator{IIP}(f.mass_matrix, dt, J, u, jacvec) @@ -863,7 +863,7 @@ function build_J_W(alg, u, uprev, p, t, dt, f::F, ::Type{uEltypeNoUnits}, deepcopy(f.jac_prototype) end jacvec = SparseDiffTools.JacVec(UJacobianWrapper(_f, t, p), copy(u), - OrdinaryDiffEqTag(), autodiff = alg_autodiff(alg)) + p, t; autodiff = alg_autodiff(alg)) W = WOperator{IIP}(f.mass_matrix, dt, J, u, jacvec) elseif islin || (!IIP && DiffEqBase.has_jac(f)) From 646684583a364f5fa0de05396be9a82ed472372b Mon Sep 17 00:00:00 2001 From: Gaurav Arya Date: Mon, 10 Apr 2023 15:01:14 -0400 Subject: [PATCH 06/54] Update usage of alg_autodiff --- src/OrdinaryDiffEq.jl | 1 + src/derivative_utils.jl | 7 +++++ src/derivative_wrappers.jl | 37 +++++++++++++++++-------- src/misc_utils.jl | 6 ++-- src/perform_step/linear_perform_step.jl | 1 + 5 files changed, 38 insertions(+), 14 deletions(-) diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index e09f4f73d1..526042de71 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -31,6 +31,7 @@ import DiffEqBase: solve!, step!, initialize!, isadaptive import DiffEqBase: ODE_DEFAULT_NORM, ODE_DEFAULT_ISOUTOFDOMAIN, ODE_DEFAULT_PROG_MESSAGE, ODE_DEFAULT_UNSTABLE_CHECK +# TODO: adjust all uses of the below two using DiffEqBase: DiffEqArrayOperator, DEFAULT_UPDATE_FUNC using DiffEqBase: TimeGradientWrapper, UJacobianWrapper, TimeDerivativeWrapper, diff --git a/src/derivative_utils.jl b/src/derivative_utils.jl index cde69376e8..4ca71a3372 100644 --- a/src/derivative_utils.jl +++ b/src/derivative_utils.jl @@ -301,6 +301,12 @@ function update_coefficients!(W::WOperator, u, p, t) W end +function DiffEqBase.update_coefficients!(J::FunctionOperator{UJacobianWrapper}, u, p, t) + copyto!(J.x, u) + J.f.t = t + J.f.p = p +end + function Base.convert(::Type{AbstractMatrix}, W::WOperator{IIP}) where {IIP} if !IIP # Allocating @@ -445,6 +451,7 @@ function do_newJW(integrator, alg, nlsolver, repeat_step)::NTuple{2, Bool} isfreshJ = isJcurrent(nlsolver, integrator) && !integrator.u_modified iszero(nlsolver.fast_convergence_cutoff) && return isfs && !isfreshJ, isfs mm = integrator.f.mass_matrix + # TODO: adjust is_varying_mm = mm isa DiffEqArrayOperator && mm.update_func !== SciMLBase.DEFAULT_UPDATE_FUNC if isfreshJ diff --git a/src/derivative_wrappers.jl b/src/derivative_wrappers.jl index 8e4cb0bd2a..209e58f263 100644 --- a/src/derivative_wrappers.jl +++ b/src/derivative_wrappers.jl @@ -80,7 +80,7 @@ function derivative!(df::AbstractArray{<:Number}, f, integrator, grad_config) alg = unwrap_alg(integrator, true) tmp = length(x) # We calculate derivtive for all elements in gradient - if alg_autodiff(alg) + if alg_autodiff(alg) isa AutoForwardDiff T = if standardtag(alg) typeof(ForwardDiff.Tag(OrdinaryDiffEqTag(), eltype(df))) else @@ -102,7 +102,7 @@ function derivative!(df::AbstractArray{<:Number}, f, df .= first.(ForwardDiff.partials.(grad_config)) integrator.stats.nf += 1 - else + elseif alg_autodiff(alg) isa AutoFiniteDiff FiniteDiff.finite_difference_gradient!(df, f, x, grad_config, dir = diffdir(integrator)) fdtype = alg_difftype(alg) @@ -113,6 +113,8 @@ function derivative!(df::AbstractArray{<:Number}, f, end end integrator.stats.nf += tmp + else + error("$alg_autodiff not yet supported in derivative! function") end nothing end @@ -122,7 +124,7 @@ function derivative(f, x::Union{Number, AbstractArray{<:Number}}, local d tmp = length(x) # We calculate derivtive for all elements in gradient alg = unwrap_alg(integrator, true) - if alg_autodiff(alg) + if alg_autodiff(alg) isa AutoForwardDiff integrator.stats.nf += 1 if integrator.iter == 1 try @@ -133,7 +135,7 @@ function derivative(f, x::Union{Number, AbstractArray{<:Number}}, else d = ForwardDiff.derivative(f, x) end - else + elseif alg_autodiff(alg) isa AutoFiniteDiff d = FiniteDiff.finite_difference_derivative(f, x, alg_difftype(alg), dir = diffdir(integrator)) if alg_difftype(alg) === Val{:central} || alg_difftype(alg) === Val{:forward} @@ -141,6 +143,8 @@ function derivative(f, x::Union{Number, AbstractArray{<:Number}}, end integrator.stats.nf += tmp d + else + error("$alg_autodiff not yet supported in derivative function") end end @@ -186,7 +190,7 @@ end function jacobian(f, x, integrator) alg = unwrap_alg(integrator, true) local tmp - if alg_autodiff(alg) + if alg_autodiff(alg) isa AutoForwardDiff if integrator.iter == 1 try J, tmp = jacobian_autodiff(f, x, integrator.f, alg) @@ -196,12 +200,14 @@ function jacobian(f, x, integrator) else J, tmp = jacobian_autodiff(f, x, integrator.f, alg) end - else + elseif alg_autodiff(alg) isa AutoFiniteDiff jac_prototype = integrator.f.jac_prototype sparsity, colorvec = sparsity_colorvec(integrator.f, x) dir = diffdir(integrator) J, tmp = jacobian_finitediff(f, x, alg_difftype(alg), dir, colorvec, sparsity, jac_prototype) + else + bleh end integrator.stats.nf += tmp J @@ -222,7 +228,7 @@ function jacobian!(J::AbstractMatrix{<:Number}, f, x::AbstractArray{<:Number}, fx::AbstractArray{<:Number}, integrator::DiffEqBase.DEIntegrator, jac_config) alg = unwrap_alg(integrator, true) - if alg_autodiff(alg) + if alg_autodiff(alg) isa AutoForwardDiff if integrator.iter == 1 try forwarddiff_color_jacobian!(J, f, x, jac_config) @@ -233,7 +239,7 @@ function jacobian!(J::AbstractMatrix{<:Number}, f, x::AbstractArray{<:Number}, forwarddiff_color_jacobian!(J, f, x, jac_config) end integrator.stats.nf += 1 - else + elseif alg_autodiff(alg) isa AutoFiniteDiff isforward = alg_difftype(alg) === Val{:forward} if isforward forwardcache = get_tmp_cache(integrator, alg, unwrap_cache(integrator, true))[2] @@ -245,6 +251,8 @@ function jacobian!(J::AbstractMatrix{<:Number}, f, x::AbstractArray{<:Number}, tmp = jacobian_finitediff!(J, f, x, jac_config, integrator) end integrator.stats.nf += tmp + else + error("$alg_autodiff not yet supported in jacobian! function") end nothing end @@ -272,7 +280,8 @@ function build_jac_config(alg, f::F1, uf::F2, du1, uprev, u, tmp, du2, end sparsity, colorvec = sparsity_colorvec(f, u) - if alg_autodiff(alg) + # TODO: more generc, do we need this? + if alg_autodiff(alg) isa AutoForwardDiff _chunksize = get_chunksize(alg) === Val(0) ? nothing : get_chunksize(alg) # SparseDiffEq uses different convection... T = if standardtag(alg) @@ -282,7 +291,7 @@ function build_jac_config(alg, f::F1, uf::F2, du1, uprev, u, tmp, du2, end jac_config = ForwardColorJacCache(uf, uprev, _chunksize; colorvec = colorvec, sparsity = sparsity, tag = T) - else + elseif alg_autodiff(alg) isa AutoFiniteDiff if alg_difftype(alg) !== Val{:complex} jac_config = FiniteDiff.JacobianCache(tmp, du1, du2, alg_difftype(alg), colorvec = colorvec, @@ -294,6 +303,8 @@ function build_jac_config(alg, f::F1, uf::F2, du1, uprev, u, tmp, du2, colorvec = colorvec, sparsity = sparsity) end + else + error("$alg_autodiff not yet supported in build_jac_config function") end else jac_config = nothing @@ -343,7 +354,7 @@ end function build_grad_config(alg, f::F1, tf::F2, du1, t) where {F1, F2} if !DiffEqBase.has_tgrad(f) - if alg_autodiff(alg) + if alg_autodiff(alg) isa AutoForwardDiff T = if standardtag(alg) typeof(ForwardDiff.Tag(OrdinaryDiffEqTag(), eltype(du1))) else @@ -362,8 +373,10 @@ function build_grad_config(alg, f::F1, tf::F2, du1, t) where {F1, F2} (ForwardDiff.Partials((one(eltype(du1)),)),)) .* false) end - else + elseif alg_autodiff(alg) isa AutoFiniteDiff grad_config = FiniteDiff.GradientCache(du1, t, alg_difftype(alg)) + else + error("$alg_autodiff not yet supported in build_grad_config function") end else grad_config = nothing diff --git a/src/misc_utils.jl b/src/misc_utils.jl index 7045e86269..3b7f692f24 100644 --- a/src/misc_utils.jl +++ b/src/misc_utils.jl @@ -109,10 +109,12 @@ function dolinsolve(integrator, linsolve; A = nothing, linu = nothing, b = nothi if integrator isa SciMLBase.DEIntegrator && _alg.linsolve !== nothing && !LinearSolve.needs_concrete_A(_alg.linsolve) && linsolve.A isa WOperator && linsolve.A.J isa AbstractSciMLOperator - if alg_autodiff(_alg) + if alg_autodiff(_alg) isa AutoForwardDiff integrator.stats.nf += linres.iters - else + elseif alg_autodiff(_alg) isa AutoFiniteDiff integrator.stats.nf += 2 * linres.iters + else + error("$alg_autodiff not yet supported in dolinsolve function") end end diff --git a/src/perform_step/linear_perform_step.jl b/src/perform_step/linear_perform_step.jl index 168c5304cb..44d138a1f7 100644 --- a/src/perform_step/linear_perform_step.jl +++ b/src/perform_step/linear_perform_step.jl @@ -778,6 +778,7 @@ function perform_step!(integrator, cache::CayleyEulerConstantCache, repeat_step A = f.f end + # TODO: this is not in place, think about this L = update_coefficients(A, uprev, p, t) V = cay(L * dt) u = V * uprev * transpose(V) From 6025e3780e93484d70f02ccbb3da92c909bbf2f9 Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Mon, 8 May 2023 11:58:57 -0400 Subject: [PATCH 07/54] added scimloperators v0.2.2 --- Project.toml | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index c71196b7e3..ad2adbeb61 100644 --- a/Project.toml +++ b/Project.toml @@ -32,6 +32,7 @@ RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SciMLNLSolve = "e9a6253c-8580-4d32-9898-8661bb511710" +SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961" SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7" SimpleUnPack = "ce78b400-467f-4804-87d8-8f486da07d0a" SnoopPrecompile = "66db9d55-30c0-4569-8b51-7e840670fc0c" @@ -66,6 +67,7 @@ Preferences = "1.3" RecursiveArrayTools = "2.36" Reexport = "0.2, 1.0" SciMLBase = "1.90" +SciMLOperators = "0.2.2" SciMLNLSolve = "0.1" SimpleNonlinearSolve = "0.1.4" SimpleUnPack = "1" @@ -100,4 +102,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [targets] -test = ["Calculus", "ComponentArrays", "Symbolics", "AlgebraicMultigrid", "IncompleteLU", "DiffEqCallbacks", "DiffEqDevTools", "ODEProblemLibrary", "ElasticArrays", "InteractiveUtils", "PoissonRandom", "Printf", "Random", "ReverseDiff", "SafeTestsets", "SparseArrays", "Statistics", "Test", "Unitful", "ModelingToolkit", "Pkg"] \ No newline at end of file +test = ["Calculus", "ComponentArrays", "Symbolics", "AlgebraicMultigrid", "IncompleteLU", "DiffEqCallbacks", "DiffEqDevTools", "ODEProblemLibrary", "ElasticArrays", "InteractiveUtils", "PoissonRandom", "Printf", "Random", "ReverseDiff", "SafeTestsets", "SparseArrays", "Statistics", "Test", "Unitful", "ModelingToolkit", "Pkg"] From c2e13f33ae3c43a5c3288f7293fce9f573e87248 Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Mon, 8 May 2023 11:59:21 -0400 Subject: [PATCH 08/54] diffeqarrayop -> matrixop --- src/OrdinaryDiffEq.jl | 8 +++----- src/alg_utils.jl | 2 +- src/derivative_utils.jl | 15 +++++++-------- src/integrators/integrator_interface.jl | 2 +- 4 files changed, 12 insertions(+), 15 deletions(-) diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index ac0910c3b8..0c02264387 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -31,16 +31,14 @@ import DiffEqBase: solve!, step!, initialize!, isadaptive import DiffEqBase: ODE_DEFAULT_NORM, ODE_DEFAULT_ISOUTOFDOMAIN, ODE_DEFAULT_PROG_MESSAGE, ODE_DEFAULT_UNSTABLE_CHECK -# TODO: adjust all uses of the below two -using DiffEqBase: DiffEqArrayOperator, DEFAULT_UPDATE_FUNC +using SciMLOperators +using SciMLOperators: AbstractSciMLOperator, DEFAULT_UPDATE_FUNC using DiffEqBase: TimeGradientWrapper, UJacobianWrapper, TimeDerivativeWrapper, UDerivativeWrapper using DiffEqBase: DEIntegrator -import SciMLBase: update_coefficients! - import RecursiveArrayTools: chain, recursivecopy! using SimpleUnPack, ForwardDiff, RecursiveArrayTools, @@ -77,7 +75,7 @@ using FastBroadcast: @.., True, False using IfElse -using SciMLBase: NoInit, _unwrap_val, AbstractSciMLOperator +using SciMLBase: NoInit, _unwrap_val import DiffEqBase: calculate_residuals, calculate_residuals!, unwrap_cache, @tight_loop_macros, diff --git a/src/alg_utils.jl b/src/alg_utils.jl index 13ce463f2e..36f54f9f30 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -250,7 +250,7 @@ function DiffEqBase.prepare_alg(alg::Union{ !(typeof(prob.f.jac_prototype) <: AbstractSciMLOperator))) linsolve = LinearSolve.defaultalg(prob.f.jac_prototype, u0) else - # If mm is a sparse matrix and A is a DiffEqArrayOperator, then let linear + # If mm is a sparse matrix and A is a MatrixOperator, then let linear # solver choose things later linsolve = nothing end diff --git a/src/derivative_utils.jl b/src/derivative_utils.jl index 00b077add2..b63c69b28d 100644 --- a/src/derivative_utils.jl +++ b/src/derivative_utils.jl @@ -210,9 +210,9 @@ mutable struct WOperator{IIP, T, end _func_cache = nothing else - AJ = J isa DiffEqArrayOperator ? convert(AbstractMatrix, J) : J + AJ = J isa MatrixOperator ? convert(AbstractMatrix, J) : J if AJ isa AbstractMatrix - mm = mass_matrix isa DiffEqArrayOperator ? + mm = mass_matrix isa MatrixOperator ? convert(AbstractMatrix, mass_matrix) : mass_matrix if AJ isa AbstractSparseMatrix @@ -223,7 +223,7 @@ mutable struct WOperator{IIP, T, # # Constant operators never refactorize so always use the correct values there # as well - if gamma == 0 && !(J isa DiffEqArrayOperator && SciMLBase.isconstant(J)) + if gamma == 0 && !(J isa MatrixOperator && SciMLBase.isconstant(J)) # Workaround https://github.com/JuliaSparse/SparseArrays.jl/issues/190 # Hopefully `rand()` does not match any value in the array (prob ~ 0, with a check) # Then `one` is required since gamma is zero @@ -285,7 +285,7 @@ function WOperator{IIP}(f, u, gamma; transform = false) where {IIP} J = deepcopy(f.jac_prototype) if J isa AbstractMatrix @assert DiffEqBase.has_jac(f) "f needs to have an associated jacobian" - J = DiffEqArrayOperator(J; update_func = f.jac) + J = MatrixOperator(J; update_func! = f.jac) end return WOperator{IIP}(mass_matrix, gamma, J, u; transform = transform) end @@ -452,8 +452,7 @@ function do_newJW(integrator, alg, nlsolver, repeat_step)::NTuple{2, Bool} iszero(nlsolver.fast_convergence_cutoff) && return isfs && !isfreshJ, isfs mm = integrator.f.mass_matrix # TODO: adjust - is_varying_mm = mm isa DiffEqArrayOperator && - mm.update_func !== SciMLBase.DEFAULT_UPDATE_FUNC + is_varying_mm = isconstant(mm) if isfreshJ jbad = false smallstepchange = true @@ -740,7 +739,7 @@ end else if !isa(J, AbstractSciMLOperator) && (!isnewton(nlsolver) || nlsolver.cache.W.J isa AbstractSciMLOperator) - J = DiffEqArrayOperator(J) + J = MatrixOperator(J) end W = WOperator{false}(mass_matrix, dtgamma, J, uprev, cache.W.jacvec; transform = W_transform) @@ -876,7 +875,7 @@ function build_J_W(alg, u, uprev, p, t, dt, f::F, ::Type{uEltypeNoUnits}, elseif islin || (!IIP && DiffEqBase.has_jac(f)) J = islin ? (isode ? f.f : f.f1.f) : f.jac(uprev, p, t) # unwrap the Jacobian accordingly if !isa(J, AbstractSciMLOperator) - J = DiffEqArrayOperator(J) + J = MatrixOperator(J) end W = WOperator{IIP}(f.mass_matrix, dt, J, u) else diff --git a/src/integrators/integrator_interface.jl b/src/integrators/integrator_interface.jl index 92cb6df745..133698de45 100644 --- a/src/integrators/integrator_interface.jl +++ b/src/integrators/integrator_interface.jl @@ -238,7 +238,7 @@ function resize_J_W!(cache, integrator, i) resize!(cache.J, i) elseif f.jac_prototype !== nothing J = similar(f.jac_prototype, i, i) - J = DiffEqArrayOperator(J; update_func = f.jac) + J = MatrixOperator(J; update_func! = f.jac) end if cache.W.jacvec isa AbstractSciMLOperator # TODO: resize! will need to be implemented upstream to handle previously handled cases From f6ce742241b6301e8dc04c4d4120c94ce84d9fb4 Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Mon, 8 May 2023 12:26:07 -0400 Subject: [PATCH 09/54] comments and fixes --- src/OrdinaryDiffEq.jl | 3 ++- src/derivative_utils.jl | 21 ++++++++++----------- src/integrators/integrator_interface.jl | 1 - src/perform_step/linear_perform_step.jl | 2 +- 4 files changed, 13 insertions(+), 14 deletions(-) diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index 0c02264387..1e0a3d0fd1 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -75,6 +75,7 @@ using FastBroadcast: @.., True, False using IfElse +using SciMLBase using SciMLBase: NoInit, _unwrap_val import DiffEqBase: calculate_residuals, calculate_residuals!, unwrap_cache, @@ -87,7 +88,7 @@ else struct OrdinaryDiffEqTag end end -import SparseDiffTools +using SparseDiffTools import SparseDiffTools: matrix_colors, forwarddiff_color_jacobian!, forwarddiff_color_jacobian, ForwardColorJacCache, default_chunk_size, getsize diff --git a/src/derivative_utils.jl b/src/derivative_utils.jl index b63c69b28d..e0be5bad20 100644 --- a/src/derivative_utils.jl +++ b/src/derivative_utils.jl @@ -1,6 +1,6 @@ const ROSENBROCK_INV_CUTOFF = 7 # https://github.com/SciML/OrdinaryDiffEq.jl/pull/1539 -struct StaticWOperator{isinv, T} +struct StaticWOperator{isinv, T} <: AbstractSciMLOperator{T} W::T function StaticWOperator(W::T, callinv = true) where {T} isinv = size(W, 1) <= ROSENBROCK_INV_CUTOFF @@ -294,14 +294,14 @@ SciMLBase.isinplace(::WOperator{IIP}, i) where {IIP} = IIP Base.eltype(W::WOperator) = eltype(W.J) set_gamma!(W::WOperator, gamma) = (W.gamma = gamma; W) -function update_coefficients!(W::WOperator, u, p, t) +function SciMLOperators.update_coefficients!(W::WOperator, u, p, t) update_coefficients!(W.J, u, p, t) update_coefficients!(W.mass_matrix, u, p, t) - W.jacvec !== nothing && update_coefficients!(W.jacvec, u, p, t) + !isnothing(W.jacvec) && update_coefficients!(W.jacvec, u, p, t) W end -function DiffEqBase.update_coefficients!(J::FunctionOperator{UJacobianWrapper}, u, p, t) +function SciMLOperators.update_coefficients!(J::FunctionOperator{UJacobianWrapper}, u, p, t) copyto!(J.x, u) J.f.t = t J.f.p = p @@ -451,7 +451,6 @@ function do_newJW(integrator, alg, nlsolver, repeat_step)::NTuple{2, Bool} isfreshJ = isJcurrent(nlsolver, integrator) && !integrator.u_modified iszero(nlsolver.fast_convergence_cutoff) && return isfs && !isfreshJ, isfs mm = integrator.f.mass_matrix - # TODO: adjust is_varying_mm = isconstant(mm) if isfreshJ jbad = false @@ -675,7 +674,7 @@ function calc_W!(W, integrator, nlsolver::Union{Nothing, AbstractNLSolver}, cach # calculate W if W isa WOperator - isnewton(nlsolver) || DiffEqBase.update_coefficients!(W, uprev, p, t) # we will call `update_coefficients!` in NLNewton + isnewton(nlsolver) || update_coefficients!(W, uprev, p, t) # we will call `update_coefficients!` in NLNewton W.transform = W_transform set_gamma!(W, dtgamma) if W.J !== nothing && !(W.J isa AbstractSciMLOperator) @@ -766,7 +765,7 @@ end end end (W isa WOperator && unwrap_alg(integrator, true) isa NewtonAlgorithm) && - (W = DiffEqBase.update_coefficients!(W, uprev, p, t)) # we will call `update_coefficients!` in NLNewton + (W = update_coefficients!(W, uprev, p, t)) # we will call `update_coefficients!` in NLNewton is_compos && (integrator.eigen_est = isarray ? constvalue(opnorm(J, Inf)) : integrator.opts.internalnorm(J, t)) return W @@ -852,8 +851,8 @@ function build_J_W(alg, u, uprev, p, t, dt, f::F, ::Type{uEltypeNoUnits}, # be overridden with concrete_jac. _f = islin ? (isode ? f.f : f.f1.f) : f - jacvec = SparseDiffTools.JacVec(UJacobianWrapper(_f, t, p), copy(u), - p, t; autodiff = alg_autodiff(alg)) + jacvec = JacVec(UJacobianWrapper(_f, t, p), copy(u), + p, t; autodiff = alg_autodiff(alg)) J = jacvec W = WOperator{IIP}(f.mass_matrix, dt, J, u, jacvec) @@ -868,8 +867,8 @@ function build_J_W(alg, u, uprev, p, t, dt, f::F, ::Type{uEltypeNoUnits}, else deepcopy(f.jac_prototype) end - jacvec = SparseDiffTools.JacVec(UJacobianWrapper(_f, t, p), copy(u), - p, t; autodiff = alg_autodiff(alg)) + jacvec = JacVec(UJacobianWrapper(_f, t, p), copy(u), + p, t; autodiff = alg_autodiff(alg)) W = WOperator{IIP}(f.mass_matrix, dt, J, u, jacvec) elseif islin || (!IIP && DiffEqBase.has_jac(f)) diff --git a/src/integrators/integrator_interface.jl b/src/integrators/integrator_interface.jl index 133698de45..3cde50a262 100644 --- a/src/integrators/integrator_interface.jl +++ b/src/integrators/integrator_interface.jl @@ -241,7 +241,6 @@ function resize_J_W!(cache, integrator, i) J = MatrixOperator(J; update_func! = f.jac) end if cache.W.jacvec isa AbstractSciMLOperator - # TODO: resize! will need to be implemented upstream to handle previously handled cases resize!(cache.W.jacvec, i) end cache.W = WOperator{DiffEqBase.isinplace(integrator.sol.prob)}(f.mass_matrix, diff --git a/src/perform_step/linear_perform_step.jl b/src/perform_step/linear_perform_step.jl index 804c9595b8..7cc9577f4b 100644 --- a/src/perform_step/linear_perform_step.jl +++ b/src/perform_step/linear_perform_step.jl @@ -819,7 +819,7 @@ function perform_step!(integrator, cache::CayleyEulerConstantCache, repeat_step A = f.f end - # TODO: this is not in place, think about this + # OOP update_coefficients() L = update_coefficients(A, uprev, p, t) V = cay(L * dt) u = V * uprev * transpose(V) From f125b59d501919fc66d8bfadaca242a12b88850b Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Mon, 8 May 2023 12:39:21 -0400 Subject: [PATCH 10/54] diffeqarrayop -> matrixop in tests --- test/algconvergence/linear_method_tests.jl | 46 +++++++++---------- .../linear_nonlinear_convergence_tests.jl | 10 ++-- .../linear_nonlinear_krylov_tests.jl | 6 +-- test/integrators/ode_event_tests.jl | 2 +- test/interface/ad_tests.jl | 2 +- .../interface/linear_solver_split_ode_test.jl | 4 +- test/interface/linear_solver_test.jl | 4 +- test/interface/mass_matrix_tests.jl | 8 ++-- test/interface/nonfulldiagonal_sparse.jl | 6 +-- test/interface/utility_tests.jl | 7 +-- 10 files changed, 46 insertions(+), 49 deletions(-) diff --git a/test/algconvergence/linear_method_tests.jl b/test/algconvergence/linear_method_tests.jl index 7e3f44bfa3..237c85be62 100644 --- a/test/algconvergence/linear_method_tests.jl +++ b/test/algconvergence/linear_method_tests.jl @@ -2,7 +2,7 @@ using OrdinaryDiffEq, Test, DiffEqDevTools using LinearAlgebra, Random # Linear exponential solvers -A = DiffEqArrayOperator([2.0 -1.0; -1.0 2.0]) +A = MatrixOperator([2.0 -1.0; -1.0 2.0]) u0 = ones(2) prob = ODEProblem(A, u0, (0.0, 1.0)) solve(prob, LinearExponential(krylov = :off)) @@ -19,13 +19,13 @@ sol_analytic = exp(1.0 * Matrix(A)) * u0 @test isapprox(sol4, sol_analytic, rtol = 1e-8) # u' = A(t)u solvers -function update_func(A, u, p, t) +function update_func!(A, u, p, t) A[1, 1] = 0 A[2, 1] = sin(u[1]) A[1, 2] = -1 A[2, 2] = 0 end -A = DiffEqArrayOperator(ones(2, 2), update_func = update_func) +A = MatrixOperator(ones(2, 2), update_func! = update_func!) prob = ODEProblem(A, ones(2), (10, 50.0)) sol1 = solve(prob, OrdinaryDiffEq.Vern9(), dt = 1 / 4) sol2 = solve(prob, OrdinaryDiffEq.RKMK2(), dt = 1 / 4) @@ -34,7 +34,7 @@ test_setup = Dict(:alg => Vern9(), :reltol => 1e-14, :abstol => 1e-14) sim = analyticless_test_convergence(dts, prob, RKMK2(), test_setup) @test sim.𝒪est[:l2]≈2 atol=0.2 -A = DiffEqArrayOperator(ones(2, 2), update_func = update_func) +A = MatrixOperator(ones(2, 2), update_func! = update_func!) prob = ODEProblem(A, ones(2), (0, 30.0)) sol1 = solve(prob, OrdinaryDiffEq.Vern9(), dt = 1 / 4) sol2 = solve(prob, OrdinaryDiffEq.RKMK4(), dt = 1 / 4) @@ -43,7 +43,7 @@ test_setup = Dict(:alg => Vern9(), :reltol => 1e-14, :abstol => 1e-14) sim = analyticless_test_convergence(dts, prob, RKMK4(), test_setup) @test sim.𝒪est[:l2]≈4 atol=0.22 -A = DiffEqArrayOperator(ones(2, 2), update_func = update_func) +A = MatrixOperator(ones(2, 2), update_func! = update_func!) prob = ODEProblem(A, ones(2), (0, 30.0)) sol1 = solve(prob, OrdinaryDiffEq.Vern9(), dt = 1 / 4) sol2 = solve(prob, OrdinaryDiffEq.LieRK4(), dt = 1 / 4) @@ -52,7 +52,7 @@ test_setup = Dict(:alg => Vern9(), :reltol => 1e-14, :abstol => 1e-14) sim = analyticless_test_convergence(dts, prob, LieRK4(), test_setup) @test sim.𝒪est[:l2]≈5 atol=0.2 -A = DiffEqArrayOperator(ones(2, 2), update_func = update_func) +A = MatrixOperator(ones(2, 2), update_func! = update_func!) prob = ODEProblem(A, ones(2), (0, 30.0)) sol1 = solve(prob, OrdinaryDiffEq.Vern9(), dt = 1 / 4) sol2 = solve(prob, OrdinaryDiffEq.CG2(), dt = 1 / 4) @@ -61,7 +61,7 @@ test_setup = Dict(:alg => Vern9(), :reltol => 1e-14, :abstol => 1e-14) sim = analyticless_test_convergence(dts, prob, CG2(), test_setup) @test sim.𝒪est[:l2]≈2 atol=0.2 -A = DiffEqArrayOperator(ones(2, 2), update_func = update_func) +A = MatrixOperator(ones(2, 2), update_func! = update_func!) prob = ODEProblem(A, ones(2), (0, 20.0)) sol1 = solve(prob, OrdinaryDiffEq.Vern6(), dt = 1 / 8) sol2 = solve(prob, OrdinaryDiffEq.CG3(), dt = 1 / 8) @@ -70,7 +70,7 @@ test_setup = Dict(:alg => Vern6(), :reltol => 1e-14, :abstol => 1e-14) sim = analyticless_test_convergence(dts, prob, CG3(), test_setup) @test sim.𝒪est[:l2]≈3 atol=0.2 -A = DiffEqArrayOperator(ones(2, 2), update_func = update_func) +A = MatrixOperator(ones(2, 2), update_func! = update_func!) prob = ODEProblem(A, ones(2), (0, 30.0)) sol1 = solve(prob, Vern9(), dt = 1 / 4) sol2 = solve(prob, CG4a(), dt = 1 / 4) @@ -79,26 +79,26 @@ test_setup = Dict(:alg => Vern9(), :reltol => 1e-14, :abstol => 1e-14) sim = analyticless_test_convergence(dts, prob, CG4a(), test_setup) @test sim.𝒪est[:l2]≈4 atol=0.2 -function update_func(A, u, p, t) +function update_func!(A, u, p, t) A[1, 1] = 0 A[2, 1] = 1 A[1, 2] = -2 * (1 - cos(u[2]) - u[2] * sin(u[2])) A[2, 2] = 0 end -A = DiffEqArrayOperator(ones(2, 2), update_func = update_func) +A = MatrixOperator(ones(2, 2), update_func! = update_func!) prob = ODEProblem(A, ones(2), (30, 150.0)) dts = 1 ./ 2 .^ (7:-1:1) test_setup = Dict(:alg => Tsit5(), :reltol => 1e-14, :abstol => 1e-14) sim = analyticless_test_convergence(dts, prob, MagnusAdapt4(), test_setup) @test sim.𝒪est[:l2]≈4 atol=0.2 -function update_func(A, u, p, t) +function update_func!(A, u, p, t) A[1, 1] = cos(t) A[2, 1] = sin(t) A[1, 2] = -sin(t) A[2, 2] = cos(t) end -A = DiffEqArrayOperator(ones(2, 2), update_func = update_func) +A = MatrixOperator(ones(2, 2), update_func! = update_func!) prob = ODEProblem(A, ones(2), (0.0, 5.0)) dts = 1 ./ 2 .^ (10:-1:1) sol = solve(prob, OrdinaryDiffEq.MagnusMidpoint(), dt = 1 / 4) @@ -110,7 +110,7 @@ sim = analyticless_test_convergence(dts, prob, MagnusMidpoint(), test_setup) sim = analyticless_test_convergence(dts, prob, MagnusMidpoint(krylov = true), test_setup) @test sim.𝒪est[:l2]≈2 atol=0.2 -A = DiffEqArrayOperator(ones(2, 2), update_func = update_func) +A = MatrixOperator(ones(2, 2), update_func! = update_func!) prob = ODEProblem(A, ones(2), (0.0, 5.0)) dts = 1 ./ 2 .^ (10:-1:1) sol = solve(prob, OrdinaryDiffEq.MagnusLeapfrog(), dt = 1 / 4) @@ -122,7 +122,7 @@ sim = analyticless_test_convergence(dts, prob, MagnusLeapfrog(), test_setup) sim = analyticless_test_convergence(dts, prob, MagnusLeapfrog(krylov = true), test_setup) @test sim.𝒪est[:l2]≈2 atol=0.2 -A = DiffEqArrayOperator(ones(2, 2), update_func = update_func) +A = MatrixOperator(ones(2, 2), update_func! = update_func!) prob = ODEProblem(A, ones(2), (0.5, 5.0)) dts = 1 ./ 2 .^ (10:-1:1) sol = solve(prob, OrdinaryDiffEq.LieEuler(), dt = 1 / 4) @@ -134,7 +134,7 @@ sim = analyticless_test_convergence(dts, prob, LieEuler(), test_setup) sim = analyticless_test_convergence(dts, prob, LieEuler(krylov = true), test_setup) @test sim.𝒪est[:l2]≈1 atol=0.2 -A = DiffEqArrayOperator(ones(2, 2), update_func = update_func) +A = MatrixOperator(ones(2, 2), update_func! = update_func!) prob = ODEProblem(A, ones(2), (1.0, 6.0)) dts = 1 ./ 2 .^ (10:-1:1) sol = solve(prob, MagnusGauss4(), dt = 1 / 4) @@ -146,7 +146,7 @@ sim = analyticless_test_convergence(dts, prob, MagnusGauss4(), test_setup) sim = analyticless_test_convergence(dts, prob, MagnusGauss4(krylov = true), test_setup) @test sim.𝒪est[:l2]≈4 atol=0.2 -A = DiffEqArrayOperator(ones(2, 2), update_func = update_func) +A = MatrixOperator(ones(2, 2), update_func! = update_func!) prob = ODEProblem(A, ones(2), (1.0, 6.0)) dts = 1 ./ 2 .^ (4:-1:1) test_setup = Dict(:alg => Vern9(), :reltol => 1e-14, :abstol => 1e-14) @@ -155,7 +155,7 @@ sim = analyticless_test_convergence(dts, prob, MagnusNC6(), test_setup) sim = analyticless_test_convergence(dts, prob, MagnusNC6(krylov = true), test_setup) @test sim.𝒪est[:l2]≈6 atol=0.2 -A = DiffEqArrayOperator(ones(2, 2), update_func = update_func) +A = MatrixOperator(ones(2, 2), update_func! = update_func!) prob = ODEProblem(A, ones(2), (1.0, 6.0)) sol = solve(prob, MagnusGL6(), dt = 1 / 10) dts = 1 ./ 2 .^ (4:-1:1) @@ -165,7 +165,7 @@ sim = analyticless_test_convergence(dts, prob, MagnusGL6(), test_setup) sim = analyticless_test_convergence(dts, prob, MagnusGL6(krylov = true), test_setup) @test sim.𝒪est[:l2]≈6 atol=0.3 -A = DiffEqArrayOperator(ones(2, 2), update_func = update_func) +A = MatrixOperator(ones(2, 2), update_func! = update_func!) prob = ODEProblem(A, ones(2), (0.0, 100.0)) dts = 1.775 .^ (5:-1:0) test_setup = Dict(:alg => Vern9(), :reltol => 1e-14, :abstol => 1e-14) @@ -174,7 +174,7 @@ sim = analyticless_test_convergence(dts, prob, MagnusGL8(), test_setup) sim = analyticless_test_convergence(dts, prob, MagnusGL8(krylov = true), test_setup) @test sim.𝒪est[:l2]≈8 atol=0.2 -A = DiffEqArrayOperator(ones(2, 2), update_func = update_func) +A = MatrixOperator(ones(2, 2), update_func! = update_func!) prob = ODEProblem(A, ones(2), (0.0, 100.0)) dts = 1.773 .^ (5:-1:0) test_setup = Dict(:alg => Vern9(), :reltol => 1e-14, :abstol => 1e-14) @@ -183,7 +183,7 @@ sim = analyticless_test_convergence(dts, prob, MagnusNC8(), test_setup) sim = analyticless_test_convergence(dts, prob, MagnusNC8(krylov = true), test_setup) @test sim.𝒪est[:l2]≈8 atol=0.2 -A = DiffEqArrayOperator(ones(2, 2), update_func = update_func) +A = MatrixOperator(ones(2, 2), update_func! = update_func!) prob = ODEProblem(A, ones(2), (1.0, 6.0)) dts = 1 ./ 2 .^ (7:-1:1) test_setup = Dict(:alg => Vern6(), :reltol => 1e-14, :abstol => 1e-14) @@ -192,7 +192,7 @@ sim = analyticless_test_convergence(dts, prob, MagnusGL4(), test_setup) sim = analyticless_test_convergence(dts, prob, MagnusGL4(krylov = true), test_setup) @test sim.𝒪est[:l2]≈4 atol=0.2 -A = DiffEqArrayOperator(ones(2, 2), update_func = update_func) +A = MatrixOperator(ones(2, 2), update_func! = update_func!) prob = ODEProblem(A, ones(2), (0, 20.0)) sol1 = solve(prob, OrdinaryDiffEq.Vern6(), dt = 1 / 8) sol2 = solve(prob, OrdinaryDiffEq.CG3(), dt = 1 / 8) @@ -222,13 +222,13 @@ function B(y::AbstractMatrix) return b end -function update_func(A, u, p, t) +function update_func!(A, u, p, t) A .= B(u) return nothing end η = diagm([1.0, 2, 3, 4, 5]) -A = DiffEqArrayOperator(Matrix{eltype(η)}(I(size(η, 1))), update_func = update_func) +A = MatrixOperator(Matrix{eltype(η)}(I(size(η, 1))), update_func! = update_func!) dts = 1 ./ 2 .^ (10:-1:2) tspan = (0.0, 20.0) diff --git a/test/algconvergence/linear_nonlinear_convergence_tests.jl b/test/algconvergence/linear_nonlinear_convergence_tests.jl index 06d361fe64..e6f3864020 100644 --- a/test/algconvergence/linear_nonlinear_convergence_tests.jl +++ b/test/algconvergence/linear_nonlinear_convergence_tests.jl @@ -37,7 +37,7 @@ end μ = 1.01 u0 = rand(2) A = [2.0 -1.0; -1.0 2.0] - linnonlin_f1 = DiffEqArrayOperator(A) + linnonlin_f1 = MatrixOperator(A) linnonlin_f2 = (du, u, p, t) -> du .= μ .* u linnonlin_fun_iip = SplitFunction(linnonlin_f1, linnonlin_f2; analytic = (u0, p, t) -> exp((A + μ * I) * t) * u0) @@ -102,8 +102,8 @@ end # Setup nonlinear problem A = [-2.0 1.0; 1.0 -2.0] f = (du, u, p, t) -> (mul!(du, A, u); du .-= u .^ 3) - jac_update = (J, u, p, t) -> (copyto!(J, A); J[1, 1] -= 3u[1]^2; J[2, 2] -= 3u[2]^2) - jac_prototype = DiffEqArrayOperator(zeros(2, 2); update_func = jac_update) + jac_update! = (J, u, p, t) -> (copyto!(J, A); J[1, 1] -= 3u[1]^2; J[2, 2] -= 3u[2]^2) + jac_prototype = MatrixOperator(zeros(2, 2); update_func! = jac_update!) fun = ODEFunction(f; jac_prototype = jac_prototype) Random.seed!(0) u0 = rand(2) @@ -181,8 +181,8 @@ end # Setup nonlinear problem A = [-2.0 1.0; 1.0 -2.0] f = (du, u, p, t) -> (mul!(du, A, u); du .-= u .^ 3) - jac_update = (J, u, p, t) -> (copyto!(J, A); J[1, 1] -= 3u[1]^2; J[2, 2] -= 3u[2]^2) - jac_prototype = DiffEqArrayOperator(zeros(2, 2); update_func = jac_update) + jac_update! = (J, u, p, t) -> (copyto!(J, A); J[1, 1] -= 3u[1]^2; J[2, 2] -= 3u[2]^2) + jac_prototype = MatrixOperator(zeros(2, 2); update_func! = jac_update!) fun = ODEFunction(f; jac_prototype = jac_prototype) Random.seed!(0) u0 = rand(2) diff --git a/test/algconvergence/linear_nonlinear_krylov_tests.jl b/test/algconvergence/linear_nonlinear_krylov_tests.jl index 144cbef607..0d46c76149 100644 --- a/test/algconvergence/linear_nonlinear_krylov_tests.jl +++ b/test/algconvergence/linear_nonlinear_krylov_tests.jl @@ -8,15 +8,15 @@ let N = 20 _f = (u, p, t) -> A * u - u .^ 3 _f_ip = (du, u, p, t) -> (mul!(du, A, u); du .-= u .^ 3) _jac = (u, p, t) -> A - 3 * diagm(0 => u .^ 2) - _jac_ip = (J, u, p, t) -> begin + _jac_ip! = (J, u, p, t) -> begin copyto!(J, A) @inbounds for i in 1:N J[i, i] -= 3 * u[i]^2 end end # f = ODEFunction(_f; jac=_jac) - # f_ip = ODEFunction(_f_ip; jac=_jac_ip, jac_prototype=zeros(N,N)) - jac_prototype = DiffEqArrayOperator(zeros(N, N); update_func = _jac_ip) + # f_ip = ODEFunction(_f_ip; jac=_jac_ip!, jac_prototype=zeros(N,N)) + jac_prototype = MatrixOperator(zeros(N, N); update_func! = _jac_ip!) f = ODEFunction(_f; jac_prototype = jac_prototype) f_ip = ODEFunction(_f_ip; jac_prototype = jac_prototype) prob = ODEProblem(f, u0, (0.0, 1.0)) diff --git a/test/integrators/ode_event_tests.jl b/test/integrators/ode_event_tests.jl index 891c0510ce..1b5e410347 100644 --- a/test/integrators/ode_event_tests.jl +++ b/test/integrators/ode_event_tests.jl @@ -422,7 +422,7 @@ step!(integrator, 1e-5, true) cb = SavingCallback(save_func, saved_values, saveat = t_l) u0 = normalize(rand(ComplexF64, 100)) - A = DiffEqArrayOperator(-1im * A) + A = MatrixOperator(-1im * A) prob = ODEProblem(A, u0, (0, 1.0)) solve(prob, LinearExponential(), dt = t_l[2] - t_l[1], callback = cb) @test length(saved_values.saveval) == length(t_l) diff --git a/test/interface/ad_tests.jl b/test/interface/ad_tests.jl index f6d9fe05b0..2af42e88f3 100644 --- a/test/interface/ad_tests.jl +++ b/test/interface/ad_tests.jl @@ -329,7 +329,7 @@ fn(x, solver) = objfun(x, prob, data, solver, reltol, abstol) # test AD with LinearExponential() function f(x) - K = DiffEqArrayOperator(x) + K = MatrixOperator(x) u0 = eltype(x).([1.0, 0.0]) prob = ODEProblem(K, u0, (0.0, 10.0)) sol = solve(prob, LinearExponential(), tstops = [0.0, 10.0])[2, :] diff --git a/test/interface/linear_solver_split_ode_test.jl b/test/interface/linear_solver_split_ode_test.jl index b4b023aa7f..a3f2ddadbd 100644 --- a/test/interface/linear_solver_split_ode_test.jl +++ b/test/interface/linear_solver_split_ode_test.jl @@ -12,8 +12,8 @@ tspan = (0.0, 1.0) M1 = 2ones(n) |> Diagonal #|> Array M2 = 2ones(n) |> Diagonal #|> Array -f1 = M1 |> DiffEqArrayOperator -f2 = M2 |> DiffEqArrayOperator +f1 = M1 |> MatrixOperator +f2 = M2 |> MatrixOperator prob = SplitODEProblem(f1, f2, u0, tspan) for algname in (:SBDF2, diff --git a/test/interface/linear_solver_test.jl b/test/interface/linear_solver_test.jl index f76809670a..374845cfee 100644 --- a/test/interface/linear_solver_test.jl +++ b/test/interface/linear_solver_test.jl @@ -5,8 +5,8 @@ using LinearAlgebra, Random N = 30 AA = sprand(MersenneTwister(12), N, N, 0.5) mm = sprand(MersenneTwister(123), N, N, 0.5) -A = DiffEqArrayOperator(AA) -M = DiffEqArrayOperator(mm'mm) +A = MatrixOperator(AA) +M = MatrixOperator(mm'mm) u0 = ones(N) prob = ODEProblem(ODEFunction(A; mass_matrix = M), u0, (0.0, 1.0)) diff --git a/test/interface/mass_matrix_tests.jl b/test/interface/mass_matrix_tests.jl index 3715910487..83f8e6f82f 100644 --- a/test/interface/mass_matrix_tests.jl +++ b/test/interface/mass_matrix_tests.jl @@ -35,7 +35,7 @@ function _norm_dsol(alg, prob, prob2, dt = 1 / 10) norm(sol .- sol2) end -function update_func1(A, u, p, t) +function update_func1!(A, u, p, t) A[1, 1] = cos(t) A[2, 1] = sin(t) * u[1] A[3, 1] = t^2 @@ -47,7 +47,7 @@ function update_func1(A, u, p, t) A[3, 3] = t * cos(t) + 1 end -function update_func2(A, u, p, t) +function update_func2!(A, u, p, t) A[1, 1] = cos(t) A[2, 1] = sin(t) A[3, 1] = t^2 @@ -61,8 +61,8 @@ end almost_I = Matrix{Float64}(1.01I, 3, 3) mm_A = Float64[-2 1 4; 4 -2 1; 2 1 3] -dependent_M1 = DiffEqArrayOperator(ones(3, 3), update_func = update_func1) -dependent_M2 = DiffEqArrayOperator(ones(3, 3), update_func = update_func2) +dependent_M1 = MatrixOperator(ones(3, 3), update_func! = update_func1!) +dependent_M2 = MatrixOperator(ones(3, 3), update_func! = update_func2!) @testset "Mass Matrix Accuracy Tests" for mm in (almost_I, mm_A) # test each method for exactness diff --git a/test/interface/nonfulldiagonal_sparse.jl b/test/interface/nonfulldiagonal_sparse.jl index 6f15ef050e..27018a334b 100644 --- a/test/interface/nonfulldiagonal_sparse.jl +++ b/test/interface/nonfulldiagonal_sparse.jl @@ -19,7 +19,7 @@ function enclosethetimedifferential(parameters::NamedTuple)::Function lower[1] = -1.0 upper[end] = 1.0 M = hcat(lower, sparse(diagm(-1 => du, 0 => diag, 1 => du2)), upper) - DiffEqArrayOperator(1 / dx * M) + MatrixOperator(1 / dx * M) end function second_deriv(N) @@ -32,7 +32,7 @@ function enclosethetimedifferential(parameters::NamedTuple)::Function lower[1] = 1.0 upper[end] = 1.0 M = hcat(lower, sparse(diagm(-1 => du, 0 => diag, 1 => du2)), upper) - DiffEqArrayOperator(1 / dx^2 * M) + MatrixOperator(1 / dx^2 * M) end function extender(N) @@ -45,7 +45,7 @@ function enclosethetimedifferential(parameters::NamedTuple)::Function M = vcat(transpose(lower), sparse(diagm(diag)), transpose(upper)) - DiffEqArrayOperator(1 / dx^2 * M) + MatrixOperator(1 / dx^2 * M) end bc_handler = extender(N) diff --git a/test/interface/utility_tests.jl b/test/interface/utility_tests.jl index 934892203d..44fc551812 100644 --- a/test/interface/utility_tests.jl +++ b/test/interface/utility_tests.jl @@ -24,7 +24,7 @@ using OrdinaryDiffEq, LinearAlgebra, SparseArrays, Random, Test, LinearSolve # In-place fun = ODEFunction((du, u, p, t) -> mul!(du, A, u); mass_matrix = mm, - jac_prototype = DiffEqArrayOperator(A)) + jac_prototype = MatrixOperator(A)) integrator = init(ODEProblem(fun, u0, tspan), ImplicitEuler(); adaptive = false, dt = dt) calc_W!(integrator.cache.nlsolver.cache.W, integrator, integrator.cache.nlsolver, @@ -56,10 +56,7 @@ end fun2 = ODEFunction(_f; mass_matrix = mm, jac = (u, p, t) -> t * A) fun1_ip = ODEFunction(_f_ip; mass_matrix = mm) fun2_ip = ODEFunction(_f_ip; mass_matrix = mm, - jac_prototype = DiffEqArrayOperator(similar(A); - update_func = (J, u, p, t) -> (J .= t .* - A; - J))) + jac_prototype = MatrixOperator(similar(A); update_func! = (J, u, p, t) -> J .= t .* A)) for Alg in [ImplicitEuler, Rosenbrock23, Rodas5] println(Alg) From 060e4befb02b6bca4d57ffa3e2192455e9e1e695 Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Mon, 8 May 2023 12:44:13 -0400 Subject: [PATCH 11/54] diffeqlinearop -> scimlop in newton.jl --- src/nlsolve/newton.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/nlsolve/newton.jl b/src/nlsolve/newton.jl index 49d3f23b52..f9389f1f0a 100644 --- a/src/nlsolve/newton.jl +++ b/src/nlsolve/newton.jl @@ -169,7 +169,7 @@ end b, ustep = _compute_rhs!(nlsolver, integrator, f, z) # update W - if W isa DiffEqBase.AbstractDiffEqLinearOperator + if W isa AbstractSciMLOperator update_coefficients!(W, ustep, p, tstep) end From 67af2b2973a436a38914673aaadd2e108e3407df Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Mon, 8 May 2023 13:11:18 -0400 Subject: [PATCH 12/54] rosenbroock accepts ADTypes --- src/alg_utils.jl | 17 +++++++++-------- src/caches/rosenbrock_caches.jl | 12 ++++++------ 2 files changed, 15 insertions(+), 14 deletions(-) diff --git a/src/alg_utils.jl b/src/alg_utils.jl index eaacc43345..97b2b495f0 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -320,14 +320,15 @@ end function _alg_autodiff(alg::OrdinaryDiffEqAlgorithm) error("This algorithm does not have an autodifferentiation option defined.") end -_alg_autodiff(alg::OrdinaryDiffEqAdaptiveImplicitAlgorithm{CS, AD}) where {CS, AD} = Val{AD}() -_alg_autodiff(alg::DAEAlgorithm{CS, AD}) where {CS, AD} = Val{AD}() -_alg_autodiff(alg::OrdinaryDiffEqImplicitAlgorithm{CS, AD}) where {CS, AD} = Val{AD}() -function _alg_autodiff(alg::Union{OrdinaryDiffEqExponentialAlgorithm{CS, AD}, - OrdinaryDiffEqAdaptiveExponentialAlgorithm{CS, AD}}) where { - CS, - AD - } +_alg_autodiff(::OrdinaryDiffEqAdaptiveImplicitAlgorithm{CS, AD}) where {CS, AD} = Val{AD}() +_alg_autodiff(::DAEAlgorithm{CS, AD}) where {CS, AD} = Val{AD}() +_alg_autodiff(::OrdinaryDiffEqImplicitAlgorithm{CS, AD}) where {CS, AD} = Val{AD}() +function _alg_autodiff(::Union{OrdinaryDiffEqExponentialAlgorithm{CS, AD}, + OrdinaryDiffEqAdaptiveExponentialAlgorithm{CS, AD} + } + ) where { + CS, AD, + } Val{AD}() end diff --git a/src/caches/rosenbrock_caches.jl b/src/caches/rosenbrock_caches.jl index 1cc59af055..ba1f6da0d8 100644 --- a/src/caches/rosenbrock_caches.jl +++ b/src/caches/rosenbrock_caches.jl @@ -152,7 +152,7 @@ function alg_cache(alg::Rosenbrock32, u, rate_prototype, ::Type{uEltypeNoUnits}, grad_config, reltol, alg, algebraic_vars) end -struct Rosenbrock23ConstantCache{T, TF, UF, JType, WType, F} <: OrdinaryDiffEqConstantCache +struct Rosenbrock23ConstantCache{T, TF, UF, JType, WType, F, AD} <: OrdinaryDiffEqConstantCache c₃₂::T d::T tf::TF @@ -160,7 +160,7 @@ struct Rosenbrock23ConstantCache{T, TF, UF, JType, WType, F} <: OrdinaryDiffEqCo J::JType W::WType linsolve::F - autodiff::Bool + autodiff::AD end function Rosenbrock23ConstantCache(::Type{T}, tf, uf, J, W, linsolve, autodiff) where {T} @@ -181,7 +181,7 @@ function alg_cache(alg::Rosenbrock23, u, rate_prototype, ::Type{uEltypeNoUnits}, alg_autodiff(alg)) end -struct Rosenbrock32ConstantCache{T, TF, UF, JType, WType, F} <: OrdinaryDiffEqConstantCache +struct Rosenbrock32ConstantCache{T, TF, UF, JType, WType, F, AD} <: OrdinaryDiffEqConstantCache c₃₂::T d::T tf::TF @@ -189,7 +189,7 @@ struct Rosenbrock32ConstantCache{T, TF, UF, JType, WType, F} <: OrdinaryDiffEqCo J::JType W::WType linsolve::F - autodiff::Bool + autodiff::AD end function Rosenbrock32ConstantCache(::Type{T}, tf, uf, J, W, linsolve, autodiff) where {T} @@ -416,14 +416,14 @@ jac_cache(c::Rosenbrock4Cache) = (c.J, c.W) ### Rodas methods -struct Rodas4ConstantCache{TF, UF, Tab, JType, WType, F} <: OrdinaryDiffEqConstantCache +struct Rodas4ConstantCache{TF, UF, Tab, JType, WType, F, AD} <: OrdinaryDiffEqConstantCache tf::TF uf::UF tab::Tab J::JType W::WType linsolve::F - autodiff::Bool + autodiff::AD end @cache mutable struct Rodas4Cache{uType, rateType, uNoUnitsType, JType, WType, TabType, From 2e47fe5564647d2c0f10d944e601ae9a3ad75b18 Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Mon, 8 May 2023 13:16:13 -0400 Subject: [PATCH 13/54] autodiff fix in newton.jl --- src/nlsolve/utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/nlsolve/utils.jl b/src/nlsolve/utils.jl index f940c261d5..97a1d8d717 100644 --- a/src/nlsolve/utils.jl +++ b/src/nlsolve/utils.jl @@ -74,7 +74,7 @@ mutable struct DAEResidualJacobianWrapper{AD, F, pType, duType, uType, alphaType t::tType function DAEResidualJacobianWrapper(alg, f, p, α, invγdt, tmp, uprev, t) isautodiff = alg_autodiff(alg) - if isautodiff + if isautodiff isa AutoForwardDiff tmp_du = PreallocationTools.dualcache(uprev) tmp_u = PreallocationTools.dualcache(uprev) else From a17a77e435fd1198ca56278bc6cbc12757c78d0f Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Mon, 8 May 2023 13:25:40 -0400 Subject: [PATCH 14/54] fix ad error in nlsolve/utisl.jl --- src/nlsolve/utils.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/nlsolve/utils.jl b/src/nlsolve/utils.jl index 97a1d8d717..4cb0444dfb 100644 --- a/src/nlsolve/utils.jl +++ b/src/nlsolve/utils.jl @@ -61,7 +61,7 @@ function du_alias_or_new(nlsolver::AbstractNLSolver, rate_prototype) end end -mutable struct DAEResidualJacobianWrapper{AD, F, pType, duType, uType, alphaType, gammaType, +mutable struct DAEResidualJacobianWrapper{isAD, F, pType, duType, uType, alphaType, gammaType, tmpType, uprevType, tType} <: Function f::F p::pType @@ -73,8 +73,8 @@ mutable struct DAEResidualJacobianWrapper{AD, F, pType, duType, uType, alphaType uprev::uprevType t::tType function DAEResidualJacobianWrapper(alg, f, p, α, invγdt, tmp, uprev, t) - isautodiff = alg_autodiff(alg) - if isautodiff isa AutoForwardDiff + isautodiff = alg_autodiff(alg) isa AutoForwardDiff + if isautodiff tmp_du = PreallocationTools.dualcache(uprev) tmp_u = PreallocationTools.dualcache(uprev) else @@ -87,7 +87,7 @@ mutable struct DAEResidualJacobianWrapper{AD, F, pType, duType, uType, alphaType end end -is_autodiff(m::DAEResidualJacobianWrapper{AD}) where {AD} = AD +is_autodiff(m::DAEResidualJacobianWrapper{isAD}) where {isAD} = isAD function (m::DAEResidualJacobianWrapper)(out, x) if is_autodiff(m) From 4a3a4deaab5983804de19d035eb1f8f9f11b34fb Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Mon, 8 May 2023 14:13:57 -0400 Subject: [PATCH 15/54] bump scimlops to 0.2.3 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 652cda5382..6a925fbf75 100644 --- a/Project.toml +++ b/Project.toml @@ -69,7 +69,7 @@ Preferences = "1.3" RecursiveArrayTools = "2.36" Reexport = "0.2, 1.0" SciMLBase = "1.90" -SciMLOperators = "0.2.2" +SciMLOperators = "0.2.3" SciMLNLSolve = "0.1" SimpleNonlinearSolve = "0.1.4" SimpleUnPack = "1" From f50d52cb6d918784c4a2c21b7a9db4736863773b Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Mon, 8 May 2023 14:14:12 -0400 Subject: [PATCH 16/54] AD bug in initialize_dae --- src/initialize_dae.jl | 44 +++++++++++++++++++++---------------------- 1 file changed, 22 insertions(+), 22 deletions(-) diff --git a/src/initialize_dae.jl b/src/initialize_dae.jl index 95dea03533..23c25105cc 100644 --- a/src/initialize_dae.jl +++ b/src/initialize_dae.jl @@ -142,8 +142,8 @@ function _initialize_dae!(integrator, prob::ODEProblem, alg::ShampineCollocation failed = nlsolvefail(nlsolver) @.. broadcast=false integrator.u=integrator.uprev + z else - isad = alg_autodiff(integrator.alg) - if isad + isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff + if isAD chunk = ForwardDiff.pickchunksize(length(tmp)) _tmp = PreallocationTools.dualcache(tmp, chunk) else @@ -153,7 +153,7 @@ function _initialize_dae!(integrator, prob::ODEProblem, alg::ShampineCollocation nlequation! = @closure (out, u, p) -> begin update_coefficients!(M, u, p, t) #M * (u-u0)/dt - f(u,p,t) - tmp = isad ? PreallocationTools.get_tmp(_tmp, u) : _tmp + tmp = isAD ? PreallocationTools.get_tmp(_tmp, u) : _tmp @. tmp = (u - u0) / dt mul!(_vec(out), M, _vec(tmp)) f(tmp, u, p, t) @@ -161,7 +161,7 @@ function _initialize_dae!(integrator, prob::ODEProblem, alg::ShampineCollocation nothing end - nlsolve = default_nlsolve(alg.nlsolve, isinplace, u0, isad) + nlsolve = default_nlsolve(alg.nlsolve, isinplace, u0, isAD) nlfunc = NonlinearFunction(nlequation!; jac_prototype = f.jac_prototype) @@ -265,8 +265,8 @@ function _initialize_dae!(integrator, prob::DAEProblem, f(resid, integrator.du, u0, p, t) integrator.opts.internalnorm(resid, t) <= integrator.opts.abstol && return - isad = alg_autodiff(integrator.alg) - if isad + isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff + if isAD chunk = ForwardDiff.pickchunksize(length(tmp)) _tmp = PreallocationTools.dualcache(tmp, chunk) else @@ -274,14 +274,14 @@ function _initialize_dae!(integrator, prob::DAEProblem, end nlequation! = @closure (out, u, p) -> begin - tmp = isad ? PreallocationTools.get_tmp(_tmp, u) : _tmp + tmp = isAD ? PreallocationTools.get_tmp(_tmp, u) : _tmp #M * (u-u0)/dt - f(u,p,t) @. tmp = (u - u0) / dt f(out, tmp, u, p, t) nothing end - nlsolve = default_nlsolve(alg.nlsolve, isinplace, u0, isad) + nlsolve = default_nlsolve(alg.nlsolve, isinplace, u0, isAD) nlfunc = NonlinearFunction(nlequation!; jac_prototype = f.jac_prototype) nlprob = NonlinearProblem(nlfunc, u0, p) @@ -372,8 +372,8 @@ function _initialize_dae!(integrator, prob::ODEProblem, integrator.opts.internalnorm(tmp, t) <= alg.abstol && return alg_u = @view u[algebraic_vars] - isad = alg_autodiff(integrator.alg) - if isad + isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff + if isAD chunk = ForwardDiff.pickchunksize(count(algebraic_vars)) _tmp = PreallocationTools.dualcache(tmp, chunk) _du_tmp = PreallocationTools.dualcache(similar(tmp), chunk) @@ -382,8 +382,8 @@ function _initialize_dae!(integrator, prob::ODEProblem, end nlequation! = @closure (out, x, p) -> begin - uu = isad ? PreallocationTools.get_tmp(_tmp, x) : _tmp - du_tmp = isad ? PreallocationTools.get_tmp(_du_tmp, x) : _du_tmp + uu = isAD ? PreallocationTools.get_tmp(_tmp, x) : _tmp + du_tmp = isAD ? PreallocationTools.get_tmp(_du_tmp, x) : _du_tmp copyto!(uu, integrator.u) alg_uu = @view uu[algebraic_vars] alg_uu .= x @@ -394,7 +394,7 @@ function _initialize_dae!(integrator, prob::ODEProblem, J = algebraic_jacobian(f.jac_prototype, algebraic_eqs, algebraic_vars) - nlsolve = default_nlsolve(alg.nlsolve, isinplace, u, isad) + nlsolve = default_nlsolve(alg.nlsolve, isinplace, u, isAD) nlfunc = NonlinearFunction(nlequation!; jac_prototype = J) nlprob = NonlinearProblem(nlfunc, alg_u, p) @@ -429,8 +429,8 @@ function _initialize_dae!(integrator, prob::ODEProblem, integrator.opts.internalnorm(resid, t) <= alg.abstol && return - isad = alg_autodiff(integrator.alg) - if isad + isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff + if isAD chunk = ForwardDiff.pickchunksize(count(algebraic_vars)) _tmp = PreallocationTools.dualcache(similar(u0), chunk) else @@ -445,7 +445,7 @@ function _initialize_dae!(integrator, prob::ODEProblem, end nlequation = @closure (x, _) -> begin - uu = isad ? PreallocationTools.get_tmp(_tmp, x) : _tmp + uu = isAD ? PreallocationTools.get_tmp(_tmp, x) : _tmp copyto!(uu, integrator.u) alg_u = @view uu[algebraic_vars] alg_u .= x @@ -455,7 +455,7 @@ function _initialize_dae!(integrator, prob::ODEProblem, J = algebraic_jacobian(f.jac_prototype, algebraic_eqs, algebraic_vars) - nlsolve = default_nlsolve(alg.nlsolve, isinplace, u0, isad) + nlsolve = default_nlsolve(alg.nlsolve, isinplace, u0, isAD) nlfunc = NonlinearFunction(nlequation; jac_prototype = J) nlprob = NonlinearProblem(nlfunc, u0[algebraic_vars]) @@ -499,8 +499,8 @@ function _initialize_dae!(integrator, prob::DAEProblem, error("differential_vars must be set for DAE initialization to occur. Either set consistent initial conditions, differential_vars, or use a different initialization algorithm.") end - isad = alg_autodiff(integrator.alg) - if isad + isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff + if isAD chunk = ForwardDiff.pickchunksize(length(tmp)) _tmp = PreallocationTools.dualcache(tmp, chunk) _du_tmp = PreallocationTools.dualcache(du_tmp, chunk) @@ -509,8 +509,8 @@ function _initialize_dae!(integrator, prob::DAEProblem, end nlequation! = @closure (out, x, p) -> begin - du_tmp = isad ? PreallocationTools.get_tmp(_du_tmp, x) : _du_tmp - uu = isad ? PreallocationTools.get_tmp(_tmp, x) : _tmp + du_tmp = isAD ? PreallocationTools.get_tmp(_du_tmp, x) : _du_tmp + uu = isAD ? PreallocationTools.get_tmp(_tmp, x) : _tmp @. du_tmp = ifelse(differential_vars, x, du) @. uu = ifelse(differential_vars, u, x) @@ -521,7 +521,7 @@ function _initialize_dae!(integrator, prob::DAEProblem, if alg.nlsolve !== nothing nlsolve = alg.nlsolve else - nlsolve = NewtonRaphson(autodiff = isad) + nlsolve = NewtonRaphson(autodiff = isAD) end nlfunc = NonlinearFunction(nlequation!; jac_prototype = f.jac_prototype) From c86a0515499d64d376ba1dc910f0268bb144f9c4 Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Tue, 9 May 2023 07:05:11 -0400 Subject: [PATCH 17/54] using/import fixes in ODE.jl --- src/OrdinaryDiffEq.jl | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index 003fb2f852..deb5e7fa4e 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -33,8 +33,9 @@ import DiffEqBase: solve!, step!, initialize!, isadaptive import DiffEqBase: ODE_DEFAULT_NORM, ODE_DEFAULT_ISOUTOFDOMAIN, ODE_DEFAULT_PROG_MESSAGE, ODE_DEFAULT_UNSTABLE_CHECK -using SciMLOperators -using SciMLOperators: AbstractSciMLOperator, DEFAULT_UPDATE_FUNC +import SciMLOperators: SciMLOperators, AbstractSciMLOperator, + MatrixOperator, FunctionOperator, + update_coefficients, update_coefficients!, DEFAULT_UPDATE_FUNC using DiffEqBase: TimeGradientWrapper, UJacobianWrapper, TimeDerivativeWrapper, UDerivativeWrapper @@ -77,7 +78,6 @@ using FastBroadcast: @.., True, False using IfElse -using SciMLBase using SciMLBase: NoInit, _unwrap_val import DiffEqBase: calculate_residuals, calculate_residuals!, unwrap_cache, @@ -90,11 +90,12 @@ else struct OrdinaryDiffEqTag end end -using SparseDiffTools import SparseDiffTools: matrix_colors, forwarddiff_color_jacobian!, forwarddiff_color_jacobian, ForwardColorJacCache, - default_chunk_size, getsize -using ADTypes + default_chunk_size, getsize, JacVec + +import ADTypes: AbstractADType, AutoFiniteDiff, AutoForwardDiff, AutoReverseDiff, + AutoTracker, AutoZygote, AutoEnzyme import Polyester using MacroTools, Adapt From 15c92ac6e8953d0d8f4e300995302e9eea2a5dc8 Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Tue, 9 May 2023 07:08:17 -0400 Subject: [PATCH 18/54] import/using fixes in ODE.jl --- src/OrdinaryDiffEq.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index deb5e7fa4e..cbdcd5c459 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -90,7 +90,7 @@ else struct OrdinaryDiffEqTag end end -import SparseDiffTools: matrix_colors, forwarddiff_color_jacobian!, +import SparseDiffTools: SparseDiffTools, matrix_colors, forwarddiff_color_jacobian!, forwarddiff_color_jacobian, ForwardColorJacCache, default_chunk_size, getsize, JacVec From ea3e9b8b1300a81b345cb1a32b49e9783e6d3a89 Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Tue, 9 May 2023 08:34:46 -0400 Subject: [PATCH 19/54] fix autodiff error in stiff_addsteps.jl --- src/dense/stiff_addsteps.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/dense/stiff_addsteps.jl b/src/dense/stiff_addsteps.jl index e7c497baa8..5e3f68a348 100644 --- a/src/dense/stiff_addsteps.jl +++ b/src/dense/stiff_addsteps.jl @@ -7,7 +7,7 @@ function _ode_addsteps!(k, t, uprev, u, dt, f, p, @unpack tf, uf, d = cache γ = dt * d tf.u = uprev - if cache.autodiff + if cache.autodiff isa AutoForwardDiff dT = ForwardDiff.derivative(tf, t) else dT = FiniteDiff.finite_difference_derivative(tf, t, dir = sign(dt)) @@ -180,7 +180,7 @@ function _ode_addsteps!(k, t, uprev, u, dt, f, p, cache::Rodas4ConstantCache, # Time derivative tf.u = uprev - if cache.autodiff + if cache.autodiff isa AutoForwardDiff dT = ForwardDiff.derivative(tf, t) else dT = FiniteDiff.finite_difference_derivative(tf, t, dir = sign(dt)) @@ -566,7 +566,7 @@ function _ode_addsteps!(k, t, uprev, u, dt, f, p, cache::Rosenbrock5ConstantCach # Time derivative tf.u = uprev - # if cache.autodiff + # if cache.autodiff isa AutoForwardDiff # dT = ForwardDiff.derivative(tf, t) # else dT = FiniteDiff.finite_difference_derivative(tf, t, dir = sign(dt)) From 98b01e8fb8f86e91dc6939555b4173be9d5569d4 Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Tue, 9 May 2023 09:15:23 -0400 Subject: [PATCH 20/54] tests unexpectedly passing --- test/regression/iipvsoop_tests.jl | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/test/regression/iipvsoop_tests.jl b/test/regression/iipvsoop_tests.jl index dafc2750b6..6860063f5b 100644 --- a/test/regression/iipvsoop_tests.jl +++ b/test/regression/iipvsoop_tests.jl @@ -109,8 +109,13 @@ end sol_ip = solve(prob_ip, alg, dt = 0.0125) sol_scalar = solve(prob_scalar, alg, dt = 0.0125) - @test_broken sol_ip(ts, idxs = 1) ≈ sol_scalar(ts) - @test_broken sol_ip.t ≈ sol_scalar.t && sol_ip[1, :] ≈ sol_scalar.u + if alg isa Kvaerno4 + @test_broken sol_ip(ts, idxs = 1) ≈ sol_scalar(ts) + @test_broken sol_ip.t ≈ sol_scalar.t && sol_ip[1, :] ≈ sol_scalar.u + else + @test sol_ip(ts, idxs = 1) ≈ sol_scalar(ts) + @test sol_ip.t ≈ sol_scalar.t && sol_ip[1, :] ≈ sol_scalar.u + end end working_rosenbrock_algs = [Rosenbrock23(), ROS3P(), Rodas3(), From caa7cd3c707631a24e40d2f9f2c05f351391fb4b Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Tue, 9 May 2023 10:04:52 -0400 Subject: [PATCH 21/54] TODO - fix CayleyEuler cache to use SciMLOps --- src/perform_step/linear_perform_step.jl | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/src/perform_step/linear_perform_step.jl b/src/perform_step/linear_perform_step.jl index 7cc9577f4b..b781cb723e 100644 --- a/src/perform_step/linear_perform_step.jl +++ b/src/perform_step/linear_perform_step.jl @@ -794,7 +794,9 @@ function perform_step!(integrator, cache::LinearExponentialCache, repeat_step = # integrator.k is automatically set due to aliasing end +# TODO - CayleyEuler needs to be fixed for SciMLOps cay!(tmp, A) = mul!(tmp, inv(I - 1 / 2 * A), (I + 1 / 2 * A)) +cay!(tmp, A::AbstractSciMLOperator) = @error "cannot multiply two SciMLOperators with mul!" cay(A) = inv(I - 1 / 2 * A) * (I + 1 / 2 * A) function initialize!(integrator, cache::CayleyEulerConstantCache) @@ -821,6 +823,12 @@ function perform_step!(integrator, cache::CayleyEulerConstantCache, repeat_step # OOP update_coefficients() L = update_coefficients(A, uprev, p, t) + + # TODO - find general scimloperator solution + if L isa MatrixOperator + L = L.A + end + V = cay(L * dt) u = V * uprev * transpose(V) @@ -856,6 +864,11 @@ function perform_step!(integrator, cache::CayleyEulerCache, repeat_step = false) update_coefficients!(L, uprev, p, t) + # TODO - InvertibleMatrixOperator, general solution for CaylerEuler + if L isa MatrixOperator + L = L.A + end + cay!(V, L * dt) mul!(tmp, uprev, transpose(V)) mul!(u, V, tmp) From 70166a621a9a7e8d26f58dae52b9a6065a7c0847 Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Tue, 9 May 2023 10:17:01 -0400 Subject: [PATCH 22/54] fix CayleyEuler --- src/perform_step/linear_perform_step.jl | 15 ++++----------- 1 file changed, 4 insertions(+), 11 deletions(-) diff --git a/src/perform_step/linear_perform_step.jl b/src/perform_step/linear_perform_step.jl index b781cb723e..4619d664a5 100644 --- a/src/perform_step/linear_perform_step.jl +++ b/src/perform_step/linear_perform_step.jl @@ -794,7 +794,6 @@ function perform_step!(integrator, cache::LinearExponentialCache, repeat_step = # integrator.k is automatically set due to aliasing end -# TODO - CayleyEuler needs to be fixed for SciMLOps cay!(tmp, A) = mul!(tmp, inv(I - 1 / 2 * A), (I + 1 / 2 * A)) cay!(tmp, A::AbstractSciMLOperator) = @error "cannot multiply two SciMLOperators with mul!" cay(A) = inv(I - 1 / 2 * A) * (I + 1 / 2 * A) @@ -823,14 +822,10 @@ function perform_step!(integrator, cache::CayleyEulerConstantCache, repeat_step # OOP update_coefficients() L = update_coefficients(A, uprev, p, t) - - # TODO - find general scimloperator solution - if L isa MatrixOperator - L = L.A - end + L = convert(AbstractMatrix, L) V = cay(L * dt) - u = V * uprev * transpose(V) + u = (V * uprev) * transpose(V) # Update integrator state integrator.fsallast = f(u, p, t + dt) @@ -864,10 +859,8 @@ function perform_step!(integrator, cache::CayleyEulerCache, repeat_step = false) update_coefficients!(L, uprev, p, t) - # TODO - InvertibleMatrixOperator, general solution for CaylerEuler - if L isa MatrixOperator - L = L.A - end + # TODO CayleyEuler cache. even if we store and cache cay(L * dt), we dont know how to apply inv(AddedOperator) + L = convert(AbstractMatrix, L) cay!(V, L * dt) mul!(tmp, uprev, transpose(V)) From 67e66f0b34d4a6b25992425cbcc1fa998f0463bc Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Tue, 9 May 2023 16:19:18 -0400 Subject: [PATCH 23/54] isconstant fix --- src/OrdinaryDiffEq.jl | 3 ++- src/derivative_utils.jl | 6 +++--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index cbdcd5c459..9c280c2ef0 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -35,7 +35,8 @@ import DiffEqBase: ODE_DEFAULT_NORM, ODE_DEFAULT_ISOUTOFDOMAIN, ODE_DEFAULT_PROG import SciMLOperators: SciMLOperators, AbstractSciMLOperator, MatrixOperator, FunctionOperator, - update_coefficients, update_coefficients!, DEFAULT_UPDATE_FUNC + update_coefficients, update_coefficients!, DEFAULT_UPDATE_FUNC, + isconstant using DiffEqBase: TimeGradientWrapper, UJacobianWrapper, TimeDerivativeWrapper, UDerivativeWrapper diff --git a/src/derivative_utils.jl b/src/derivative_utils.jl index e0be5bad20..4d24901e15 100644 --- a/src/derivative_utils.jl +++ b/src/derivative_utils.jl @@ -202,7 +202,7 @@ mutable struct WOperator{IIP, T, transform = false) where {IIP} # TODO: there is definitely a missing interface. # Tentative interface: `has_concrete` and `concertize(A)` - if J isa Union{Number, DiffEqScalar} + if J isa Union{Number, ScalarOperator} if transform _concrete_form = -mass_matrix / gamma + convert(Number, J) else @@ -223,7 +223,7 @@ mutable struct WOperator{IIP, T, # # Constant operators never refactorize so always use the correct values there # as well - if gamma == 0 && !(J isa MatrixOperator && SciMLBase.isconstant(J)) + if gamma == 0 && !(J isa MatrixOperator && isconstant(J)) # Workaround https://github.com/JuliaSparse/SparseArrays.jl/issues/190 # Hopefully `rand()` does not match any value in the array (prob ~ 0, with a check) # Then `one` is required since gamma is zero @@ -451,7 +451,7 @@ function do_newJW(integrator, alg, nlsolver, repeat_step)::NTuple{2, Bool} isfreshJ = isJcurrent(nlsolver, integrator) && !integrator.u_modified iszero(nlsolver.fast_convergence_cutoff) && return isfs && !isfreshJ, isfs mm = integrator.f.mass_matrix - is_varying_mm = isconstant(mm) + is_varying_mm = !isconstant(mm) if isfreshJ jbad = false smallstepchange = true From 309693e7160f273c99e6f12ec78d0e4e17b9d6e2 Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Tue, 9 May 2023 16:19:47 -0400 Subject: [PATCH 24/54] all splitode/ splitfunction tests are passing --- src/derivative_utils.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/derivative_utils.jl b/src/derivative_utils.jl index 4d24901e15..480a81961a 100644 --- a/src/derivative_utils.jl +++ b/src/derivative_utils.jl @@ -273,9 +273,9 @@ mutable struct WOperator{IIP, T, end end function WOperator{IIP}(f, u, gamma; transform = false) where {IIP} - if isa(f, Union{SplitFunction, DynamicalODEFunction}) - error("WOperator does not support $(typeof(f)) yet") - end + # if isa(f, Union{SplitFunction, DynamicalODEFunction}) + # error("WOperator does not support $(typeof(f)) yet") + # end mass_matrix = f.mass_matrix # TODO: does this play nicely with time-state dependent mass matrix? if !isa(mass_matrix, Union{AbstractMatrix, UniformScaling}) From 9132722609a2cb2a0dec37c795b84369ecbc1d5d Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Tue, 9 May 2023 16:21:49 -0400 Subject: [PATCH 25/54] update scimlops compat --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 6a925fbf75..1e33087b03 100644 --- a/Project.toml +++ b/Project.toml @@ -69,7 +69,7 @@ Preferences = "1.3" RecursiveArrayTools = "2.36" Reexport = "0.2, 1.0" SciMLBase = "1.90" -SciMLOperators = "0.2.3" +SciMLOperators = "0.2.4" SciMLNLSolve = "0.1" SimpleNonlinearSolve = "0.1.4" SimpleUnPack = "1" From 5ac7e62a2cb3e17c75c0a1a3ff3a2d654bddf69f Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Tue, 9 May 2023 16:45:15 -0400 Subject: [PATCH 26/54] reverting WOperator split function decision --- src/derivative_utils.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/derivative_utils.jl b/src/derivative_utils.jl index 480a81961a..4d24901e15 100644 --- a/src/derivative_utils.jl +++ b/src/derivative_utils.jl @@ -273,9 +273,9 @@ mutable struct WOperator{IIP, T, end end function WOperator{IIP}(f, u, gamma; transform = false) where {IIP} - # if isa(f, Union{SplitFunction, DynamicalODEFunction}) - # error("WOperator does not support $(typeof(f)) yet") - # end + if isa(f, Union{SplitFunction, DynamicalODEFunction}) + error("WOperator does not support $(typeof(f)) yet") + end mass_matrix = f.mass_matrix # TODO: does this play nicely with time-state dependent mass matrix? if !isa(mass_matrix, Union{AbstractMatrix, UniformScaling}) From 959d4b7938fd3b49022c8722eadf9a8ae725d9df Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Thu, 11 May 2023 12:04:03 -0400 Subject: [PATCH 27/54] fix ETD2Fsal constructor for scimlops --- src/caches/linear_nonlinear_caches.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/caches/linear_nonlinear_caches.jl b/src/caches/linear_nonlinear_caches.jl index ccafa6ac45..1406f2cd61 100644 --- a/src/caches/linear_nonlinear_caches.jl +++ b/src/caches/linear_nonlinear_caches.jl @@ -839,6 +839,16 @@ mutable struct ETD2Fsal{rateType} lin::rateType nl::rateType nlprev::rateType + + function ETD2Fsal(lin, nl, nlprev) + lin = convert(Number, lin) + nl = convert(Number, nl) + nlprev = convert(Number, nlprev) + + T = promote_type(eltype.((lin, nl, nlprev))...) + + new{typeof(lin)}(T(lin), T(nl), T(nlprev)) + end end function ETD2Fsal(rate_prototype) ETD2Fsal(zero(rate_prototype), zero(rate_prototype), zero(rate_prototype)) From e4244fa27857a12758d834c4fc1b0a8a50599932 Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Fri, 12 May 2023 09:06:52 -0400 Subject: [PATCH 28/54] fixed ETDFsal constructor --- src/caches/linear_nonlinear_caches.jl | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/caches/linear_nonlinear_caches.jl b/src/caches/linear_nonlinear_caches.jl index 1406f2cd61..92bd228576 100644 --- a/src/caches/linear_nonlinear_caches.jl +++ b/src/caches/linear_nonlinear_caches.jl @@ -841,13 +841,16 @@ mutable struct ETD2Fsal{rateType} nlprev::rateType function ETD2Fsal(lin, nl, nlprev) - lin = convert(Number, lin) - nl = convert(Number, nl) - nlprev = convert(Number, nlprev) + if size(lin) == () + # convert to same type if Number or AbstractSciMLScalarOperator + T = promote_type(eltype.((lin, nl, nlprev))...) - T = promote_type(eltype.((lin, nl, nlprev))...) + lin = convert(T, lin) + nl = convert(T, nl) + nlprev = convert(T, nlprev) + end - new{typeof(lin)}(T(lin), T(nl), T(nlprev)) + new{typeof(lin)}(lin, nl, nlprev) end end function ETD2Fsal(rate_prototype) From 049f5b295c72e6d61b743e9a9fdbf560b23d2648 Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Fri, 12 May 2023 09:25:03 -0400 Subject: [PATCH 29/54] bump scimlops compat --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 1e33087b03..dd6830e0d5 100644 --- a/Project.toml +++ b/Project.toml @@ -69,7 +69,7 @@ Preferences = "1.3" RecursiveArrayTools = "2.36" Reexport = "0.2, 1.0" SciMLBase = "1.90" -SciMLOperators = "0.2.4" +SciMLOperators = "0.2.5" SciMLNLSolve = "0.1" SimpleNonlinearSolve = "0.1.4" SimpleUnPack = "1" From fd1c6c026ac31a107ee13de0a080115215dabe39 Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Fri, 12 May 2023 09:25:54 -0400 Subject: [PATCH 30/54] pass autodiff tag to jacvec https://github.com/JuliaDiff/SparseDiffTools.jl/pull/241 --- src/derivative_utils.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/derivative_utils.jl b/src/derivative_utils.jl index 4d24901e15..7ad14cc8df 100644 --- a/src/derivative_utils.jl +++ b/src/derivative_utils.jl @@ -851,8 +851,8 @@ function build_J_W(alg, u, uprev, p, t, dt, f::F, ::Type{uEltypeNoUnits}, # be overridden with concrete_jac. _f = islin ? (isode ? f.f : f.f1.f) : f - jacvec = JacVec(UJacobianWrapper(_f, t, p), copy(u), - p, t; autodiff = alg_autodiff(alg)) + jacvec = JacVec(UJacobianWrapper(_f, t, p), copy(u), p, t; + autodiff = alg_autodiff(alg), tag = OrdinaryDiffEqTag()) J = jacvec W = WOperator{IIP}(f.mass_matrix, dt, J, u, jacvec) @@ -867,8 +867,8 @@ function build_J_W(alg, u, uprev, p, t, dt, f::F, ::Type{uEltypeNoUnits}, else deepcopy(f.jac_prototype) end - jacvec = JacVec(UJacobianWrapper(_f, t, p), copy(u), - p, t; autodiff = alg_autodiff(alg)) + jacvec = JacVec(UJacobianWrapper(_f, t, p), copy(u), p, t; + autodiff = alg_autodiff(alg), tag = OrdinaryDiffEqTag()) W = WOperator{IIP}(f.mass_matrix, dt, J, u, jacvec) elseif islin || (!IIP && DiffEqBase.has_jac(f)) From 6393ec1773ef2251add29ebcf345ad4d23524f35 Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Fri, 12 May 2023 09:26:41 -0400 Subject: [PATCH 31/54] unnecessary dots mess with ScalarOperators --- src/nlsolve/newton.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/nlsolve/newton.jl b/src/nlsolve/newton.jl index f9389f1f0a..f00a6f5147 100644 --- a/src/nlsolve/newton.jl +++ b/src/nlsolve/newton.jl @@ -252,10 +252,10 @@ end else ustep = @. tmp + γ * z if mass_matrix === I - ztmp = (dt .* f(ustep, p, tstep) .- z) .* invγdt + ztmp = (dt * f(ustep, p, tstep) - z) * invγdt else update_coefficients!(mass_matrix, ustep, p, tstep) - ztmp = (dt .* f(ustep, p, tstep) .- mass_matrix * z) .* invγdt + ztmp = (dt * f(ustep, p, tstep) - mass_matrix * z) * invγdt end end end From 114644147c3bb74e711bbfe69a0b705490e0cc1e Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Fri, 12 May 2023 09:26:57 -0400 Subject: [PATCH 32/54] rm more unnecessary dots --- src/perform_step/kencarp_kvaerno_perform_step.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/perform_step/kencarp_kvaerno_perform_step.jl b/src/perform_step/kencarp_kvaerno_perform_step.jl index d67824ed60..3087b1b8db 100644 --- a/src/perform_step/kencarp_kvaerno_perform_step.jl +++ b/src/perform_step/kencarp_kvaerno_perform_step.jl @@ -206,7 +206,7 @@ end if typeof(integrator.f) <: SplitFunction # Explicit tableau is not FSAL # Make this not compute on repeat - z₁ = dt .* f(uprev, p, t) + z₁ = dt * f(uprev, p, t) else # FSAL Step 1 z₁ = dt * integrator.fsalfirst From cbfe2ce22c2123d3c33660bf77328d1513a3f1c3 Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Fri, 12 May 2023 11:12:57 -0400 Subject: [PATCH 33/54] fix typo --- test/interface/dae_initialize_integration.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/interface/dae_initialize_integration.jl b/test/interface/dae_initialize_integration.jl index 3315fbdfa4..705d523bee 100644 --- a/test/interface/dae_initialize_integration.jl +++ b/test/interface/dae_initialize_integration.jl @@ -37,5 +37,5 @@ p0 = [ prob = ODEProblem(connected, u0, tspan, p0) sol = solve(prob, Rodas5(), initializealg = BrownFullBasicInit()) @test prob.u0 == sol[1] -sol = solve(prob, Rodas5(), initiailizealg = ShampineCollocationInit()) +sol = solve(prob, Rodas5(), initializealg = ShampineCollocationInit()) @test prob.u0 == sol[1] From 4f9aa2707618d454f212b14577c1deef0832faae Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Fri, 12 May 2023 12:17:24 -0400 Subject: [PATCH 34/54] fix update func in linear_method_tests --- test/algconvergence/linear_method_tests.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/test/algconvergence/linear_method_tests.jl b/test/algconvergence/linear_method_tests.jl index 237c85be62..7658218e03 100644 --- a/test/algconvergence/linear_method_tests.jl +++ b/test/algconvergence/linear_method_tests.jl @@ -222,13 +222,17 @@ function B(y::AbstractMatrix) return b end +function update_func(A, u, p, t) + B(u) +end + function update_func!(A, u, p, t) A .= B(u) return nothing end η = diagm([1.0, 2, 3, 4, 5]) -A = MatrixOperator(Matrix{eltype(η)}(I(size(η, 1))), update_func! = update_func!) +A = MatrixOperator(Matrix{eltype(η)}(I(size(η, 1))), update_func = update_func, update_func! = update_func!) dts = 1 ./ 2 .^ (10:-1:2) tspan = (0.0, 20.0) From d0ffa6c639fd0246b43db56566a2d1c94438ac46 Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Fri, 12 May 2023 12:20:13 -0400 Subject: [PATCH 35/54] split/sdirk tests remain broken --- test/regression/iipvsoop_tests.jl | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/test/regression/iipvsoop_tests.jl b/test/regression/iipvsoop_tests.jl index 6860063f5b..dafc2750b6 100644 --- a/test/regression/iipvsoop_tests.jl +++ b/test/regression/iipvsoop_tests.jl @@ -109,13 +109,8 @@ end sol_ip = solve(prob_ip, alg, dt = 0.0125) sol_scalar = solve(prob_scalar, alg, dt = 0.0125) - if alg isa Kvaerno4 - @test_broken sol_ip(ts, idxs = 1) ≈ sol_scalar(ts) - @test_broken sol_ip.t ≈ sol_scalar.t && sol_ip[1, :] ≈ sol_scalar.u - else - @test sol_ip(ts, idxs = 1) ≈ sol_scalar(ts) - @test sol_ip.t ≈ sol_scalar.t && sol_ip[1, :] ≈ sol_scalar.u - end + @test_broken sol_ip(ts, idxs = 1) ≈ sol_scalar(ts) + @test_broken sol_ip.t ≈ sol_scalar.t && sol_ip[1, :] ≈ sol_scalar.u end working_rosenbrock_algs = [Rosenbrock23(), ROS3P(), Rodas3(), From ddc979172de6573482896669ae13f3b7a74cdb71 Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Mon, 15 May 2023 08:48:44 -0400 Subject: [PATCH 36/54] update scimlop compat --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index dd6830e0d5..8f616ae126 100644 --- a/Project.toml +++ b/Project.toml @@ -69,7 +69,7 @@ Preferences = "1.3" RecursiveArrayTools = "2.36" Reexport = "0.2, 1.0" SciMLBase = "1.90" -SciMLOperators = "0.2.5" +SciMLOperators = "0.2.8" SciMLNLSolve = "0.1" SimpleNonlinearSolve = "0.1.4" SimpleUnPack = "1" From 9c0fa848f231a37db95731d612ede98adc1fdefc Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Tue, 16 May 2023 07:42:40 -0400 Subject: [PATCH 37/54] sparsedifftools compat --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 8f616ae126..5545645144 100644 --- a/Project.toml +++ b/Project.toml @@ -73,7 +73,7 @@ SciMLOperators = "0.2.8" SciMLNLSolve = "0.1" SimpleNonlinearSolve = "0.1.4" SimpleUnPack = "1" -SparseDiffTools = "2" +SparseDiffTools = "2.3" PrecompileTools = "1" StaticArrayInterface = "1.2" StaticArrays = "0.11, 0.12, 1.0" From 0b019055b3e4998d75d13c957d30e7afdbc7d89d Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Tue, 16 May 2023 07:43:33 -0400 Subject: [PATCH 38/54] rm comments from CayleyEuler in linear_perform_step --- src/perform_step/linear_perform_step.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/perform_step/linear_perform_step.jl b/src/perform_step/linear_perform_step.jl index 4619d664a5..c768bc704f 100644 --- a/src/perform_step/linear_perform_step.jl +++ b/src/perform_step/linear_perform_step.jl @@ -820,7 +820,6 @@ function perform_step!(integrator, cache::CayleyEulerConstantCache, repeat_step A = f.f end - # OOP update_coefficients() L = update_coefficients(A, uprev, p, t) L = convert(AbstractMatrix, L) @@ -858,8 +857,6 @@ function perform_step!(integrator, cache::CayleyEulerCache, repeat_step = false) end update_coefficients!(L, uprev, p, t) - - # TODO CayleyEuler cache. even if we store and cache cay(L * dt), we dont know how to apply inv(AddedOperator) L = convert(AbstractMatrix, L) cay!(V, L * dt) From 1ef9d90965fb1056741bbf4c0d473ee645c93034 Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Tue, 16 May 2023 07:45:07 -0400 Subject: [PATCH 39/54] rm Kencarp3, CFNLIRK3 split tests woperator doesn't support split problems --- test/algconvergence/linear_nonlinear_convergence_tests.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/test/algconvergence/linear_nonlinear_convergence_tests.jl b/test/algconvergence/linear_nonlinear_convergence_tests.jl index e6f3864020..edeb2f0e42 100644 --- a/test/algconvergence/linear_nonlinear_convergence_tests.jl +++ b/test/algconvergence/linear_nonlinear_convergence_tests.jl @@ -5,7 +5,7 @@ using OrdinaryDiffEq: alg_order println("Caching Out-of-place") μ = 1.01 linnonlin_f2 = (u, p, t) -> μ * u - linnonlin_f1 = DiffEqScalar(μ) + linnonlin_f1 = ScalarOperator(μ) linnonlin_fun = SplitFunction(linnonlin_f1, linnonlin_f2; analytic = (u0, p, t) -> u0 .* exp.(2μ * t)) prob = SplitODEProblem(linnonlin_fun, 1 / 2, (0.0, 1.0)) @@ -20,9 +20,10 @@ using OrdinaryDiffEq: alg_order ETDRK4, HochOst4, ETD2, - KenCarp3, - CFNLIRK3, + # KenCarp3, # WOperator doesn't support split problems + # CFNLIRK3, # WOperator doesn't support split problems ] + @info "$Alg" sim = test_convergence(dts, prob, Alg()) @test sim.𝒪est[:l2]≈alg_order(Alg()) atol=0.2 end From 0b6594e08b409802e3de6bf425466d8da6011941 Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Tue, 16 May 2023 07:45:34 -0400 Subject: [PATCH 40/54] fix update behavious in mass matrix tests --- test/interface/mass_matrix_tests.jl | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/test/interface/mass_matrix_tests.jl b/test/interface/mass_matrix_tests.jl index 83f8e6f82f..0d2093cf5c 100644 --- a/test/interface/mass_matrix_tests.jl +++ b/test/interface/mass_matrix_tests.jl @@ -47,6 +47,12 @@ function update_func1!(A, u, p, t) A[3, 3] = t * cos(t) + 1 end +function update_func1(A, u, p, t) + A = copy(A) + update_func1!(A, u, p, t) + A +end + function update_func2!(A, u, p, t) A[1, 1] = cos(t) A[2, 1] = sin(t) @@ -59,10 +65,16 @@ function update_func2!(A, u, p, t) A[3, 3] = t * cos(t) + 1 end +function update_func2(A, u, p, t) + A = copy(A) + update_func2!(A, u, p, t) + A +end + almost_I = Matrix{Float64}(1.01I, 3, 3) mm_A = Float64[-2 1 4; 4 -2 1; 2 1 3] -dependent_M1 = MatrixOperator(ones(3, 3), update_func! = update_func1!) -dependent_M2 = MatrixOperator(ones(3, 3), update_func! = update_func2!) +dependent_M1 = MatrixOperator(ones(3, 3), update_func = update_func1, update_func! = update_func1!) +dependent_M2 = MatrixOperator(ones(3, 3), update_func = update_func2, update_func! = update_func2!) @testset "Mass Matrix Accuracy Tests" for mm in (almost_I, mm_A) # test each method for exactness From 5ccc27506ed068b9e24ac8cf57d944adc1b99f88 Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Tue, 16 May 2023 08:17:00 -0400 Subject: [PATCH 41/54] flattening to Nx1 solves the iterativesolvers indexing error --- test/interface/nojac.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/test/interface/nojac.jl b/test/interface/nojac.jl index 86800bf824..9d957f2631 100644 --- a/test/interface/nojac.jl +++ b/test/interface/nojac.jl @@ -5,6 +5,9 @@ const xyd_brusselator = range(0, stop = 1, length = N) brusselator_f(x, y, t) = (((x - 0.3)^2 + (y - 0.6)^2) <= 0.1^2) * (t >= 1.1) * 5.0 limit(a, N) = a == N + 1 ? 1 : a == 0 ? N : a function brusselator_2d_loop(du, u, p, t) + u = reshape(u, N, N, 2) + du = reshape(u, N, N, 2) + A, B, alpha, dx = p alpha = alpha / dx^2 @inbounds for I in CartesianIndices((N, N)) @@ -20,6 +23,7 @@ function brusselator_2d_loop(du, u, p, t) 4u[i, j, 2]) + A * u[i, j, 1] - u[i, j, 1]^2 * u[i, j, 2] end + vec(du) end p = (3.4, 1.0, 10.0, step(xyd_brusselator)) @@ -34,7 +38,7 @@ function init_brusselator_2d(xyd) end u end -u0 = init_brusselator_2d(xyd_brusselator) +u0 = init_brusselator_2d(xyd_brusselator) |> vec prob_ode_brusselator_2d = ODEProblem(brusselator_2d_loop, u0, (0.0, 5.0), p) From 8fd1f1fea0a72a335a7585b9d0c99929d3db04fd Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Tue, 16 May 2023 08:17:33 -0400 Subject: [PATCH 42/54] flattening PDE solve --- test/interface/preconditioners.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/test/interface/preconditioners.jl b/test/interface/preconditioners.jl index 5d6180c56c..a83d29bcfc 100644 --- a/test/interface/preconditioners.jl +++ b/test/interface/preconditioners.jl @@ -13,6 +13,10 @@ limit(a, N) = a == N + 1 ? 1 : a == 0 ? N : a const iter = Ref(0) function brusselator_2d_loop(du, u, p, t) global iter[] += 1 + + u = reshape(u, N, N, 2) + du = reshape(u, N, N, 2) + A, B, alpha, dx = p alpha = alpha / dx^2 @inbounds for I in CartesianIndices((N, N)) @@ -28,6 +32,8 @@ function brusselator_2d_loop(du, u, p, t) 4u[i, j, 2]) + A * u[i, j, 1] - u[i, j, 1]^2 * u[i, j, 2] end + + vec(du) end p = (3.4, 1.0, 10.0, step(xyd_brusselator)) @@ -42,7 +48,7 @@ function init_brusselator_2d(xyd) end u end -u0 = init_brusselator_2d(xyd_brusselator) +u0 = init_brusselator_2d(xyd_brusselator) |> vec prob_ode_brusselator_2d = ODEProblem(brusselator_2d_loop, u0, (0.0, 11.5), p) From 3e7cc18a179e9612206f1d3415bcd386dea62066 Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Tue, 16 May 2023 12:03:27 -0400 Subject: [PATCH 43/54] diffeqbase compat --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 5545645144..50576b60a5 100644 --- a/Project.toml +++ b/Project.toml @@ -47,7 +47,7 @@ TruncatedStacktraces = "781d530d-4396-4725-bb49-402e4bee1e77" Adapt = "1.1, 2.0, 3.0" ArrayInterface = "6, 7" DataStructures = "0.18" -DiffEqBase = "6.122.0" +DiffEqBase = "6.125.0" DocStringExtensions = "0.8, 0.9" ExponentialUtilities = "1.22" FastBroadcast = "0.1.9, 0.2" From a5ae19ec2fa6e73672d68aded538dd8056cdee19 Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Tue, 16 May 2023 17:58:36 -0400 Subject: [PATCH 44/54] precond, nojac tests passing --- test/interface/nojac.jl | 2 +- test/interface/preconditioners.jl | 35 ++++++++++++++++++++++++++++--- 2 files changed, 33 insertions(+), 4 deletions(-) diff --git a/test/interface/nojac.jl b/test/interface/nojac.jl index 9d957f2631..0696e7101f 100644 --- a/test/interface/nojac.jl +++ b/test/interface/nojac.jl @@ -6,7 +6,7 @@ brusselator_f(x, y, t) = (((x - 0.3)^2 + (y - 0.6)^2) <= 0.1^2) * (t >= 1.1) * 5 limit(a, N) = a == N + 1 ? 1 : a == 0 ? N : a function brusselator_2d_loop(du, u, p, t) u = reshape(u, N, N, 2) - du = reshape(u, N, N, 2) + du = reshape(du, N, N, 2) A, B, alpha, dx = p alpha = alpha / dx^2 diff --git a/test/interface/preconditioners.jl b/test/interface/preconditioners.jl index a83d29bcfc..8cabf7de9a 100644 --- a/test/interface/preconditioners.jl +++ b/test/interface/preconditioners.jl @@ -15,7 +15,7 @@ function brusselator_2d_loop(du, u, p, t) global iter[] += 1 u = reshape(u, N, N, 2) - du = reshape(u, N, N, 2) + du = reshape(du, N, N, 2) A, B, alpha, dx = p alpha = alpha / dx^2 @@ -49,8 +49,7 @@ function init_brusselator_2d(xyd) u end u0 = init_brusselator_2d(xyd_brusselator) |> vec -prob_ode_brusselator_2d = ODEProblem(brusselator_2d_loop, - u0, (0.0, 11.5), p) +prob_ode_brusselator_2d = ODEProblem(brusselator_2d_loop, u0, (0.0, 11.5), p) du0 = copy(u0) jac = ModelingToolkit.Symbolics.jacobian_sparsity((du, u) -> brusselator_2d_loop(du, u, p, @@ -115,6 +114,11 @@ sol4 = solve(prob_ode_brusselator_2d_sparse, iter4 = iter[]; iter[] = 0; +@test sol1.retcode === ReturnCode.Success +@test sol2.retcode === ReturnCode.Success +@test sol3.retcode === ReturnCode.Success +@test sol4.retcode === ReturnCode.Success + @test iter2 < iter1 @test iter3 < iter1 @test iter4 < iter1 @@ -139,6 +143,11 @@ sol4 = solve(prob_ode_brusselator_2d_sparse, iter4 = iter[]; iter[] = 0; +@test sol1.retcode === ReturnCode.Success +@test sol2.retcode === ReturnCode.Success +@test sol3.retcode === ReturnCode.Success +@test sol4.retcode === ReturnCode.Success + @test iter2 < iter1 @test iter3 < iter1 @test iter4 < iter1 @@ -163,6 +172,11 @@ sol4 = solve(prob_ode_brusselator_2d_sparse, iter4 = iter[]; iter[] = 0; +@test sol1.retcode === ReturnCode.Success +@test sol2.retcode === ReturnCode.Success +@test sol3.retcode === ReturnCode.Success +@test sol4.retcode === ReturnCode.Success + @test iter2 < iter1 @test iter3 < iter1 @test iter4 < iter1 @@ -187,6 +201,11 @@ sol4 = solve(prob_ode_brusselator_2d_sparse, iter4 = iter[]; iter[] = 0; +@test sol1.retcode === ReturnCode.Success +@test sol2.retcode === ReturnCode.Success +@test sol3.retcode === ReturnCode.Success +@test sol4.retcode === ReturnCode.Success + @test iter2 < iter1 @test iter3 < iter1 @test iter4 < iter1 @@ -211,6 +230,11 @@ sol4 = solve(prob_ode_brusselator_2d_sparse, iter4 = iter[]; iter[] = 0; +@test sol1.retcode === ReturnCode.Success +@test sol2.retcode === ReturnCode.Success +@test sol3.retcode === ReturnCode.Success +@test sol4.retcode === ReturnCode.Success + @test iter2 < iter1 @test iter3 < iter1 @test iter4 < iter1 @@ -235,6 +259,11 @@ sol4 = solve(prob_ode_brusselator_2d_sparse, iter4 = iter[]; iter[] = 0; +@test sol1.retcode === ReturnCode.Success +@test sol2.retcode === ReturnCode.Success +@test sol3.retcode === ReturnCode.Success +@test sol4.retcode === ReturnCode.Success + @test iter2 < iter1 @test iter3 < iter1 @test iter4 < iter1 From 5ff1759ac6f82650dd3c982c4b1e11b507306d8b Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Wed, 17 May 2023 13:39:34 -0400 Subject: [PATCH 45/54] switch out iterativesolvers --> krylovjl to fix inf/nan errors --- test/interface/linear_nonlinear_tests.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/test/interface/linear_nonlinear_tests.jl b/test/interface/linear_nonlinear_tests.jl index 8b773b4a39..86cc328668 100644 --- a/test/interface/linear_nonlinear_tests.jl +++ b/test/interface/linear_nonlinear_tests.jl @@ -36,35 +36,35 @@ end sol = @test_nowarn solve(prob, TRBDF2(autodiff = false)); @test length(sol.t) < 20 sol = @test_nowarn solve(prob, - TRBDF2(autodiff = false, linsolve = IterativeSolversJL_GMRES())); + TRBDF2(autodiff = false, linsolve = KrylovJL_GMRES())); @test length(sol.t) < 20 solref = @test_nowarn solve(prob, - TRBDF2(autodiff = false, linsolve = IterativeSolversJL_GMRES(), + TRBDF2(autodiff = false, linsolve = KrylovJL_GMRES(), smooth_est = false)); @test length(sol.t) < 20 sol = @test_nowarn solve(prob, - TRBDF2(autodiff = false, linsolve = IterativeSolversJL_GMRES(), + TRBDF2(autodiff = false, linsolve = KrylovJL_GMRES(), precs = precsl, smooth_est = false, concrete_jac = true)); @test length(sol.t) < 20 sol = @test_nowarn solve(prob, - TRBDF2(autodiff = false, linsolve = IterativeSolversJL_GMRES(), + TRBDF2(autodiff = false, linsolve = KrylovJL_GMRES(), precs = precsr, smooth_est = false, concrete_jac = true)); @test length(sol.t) < 20 sol = @test_nowarn solve(prob, - TRBDF2(autodiff = false, linsolve = IterativeSolversJL_GMRES(), + TRBDF2(autodiff = false, linsolve = KrylovJL_GMRES(), precs = precslr, smooth_est = false, concrete_jac = true)); @test length(sol.t) < 20 sol = @test_nowarn solve(prob, - QNDF(autodiff = false, linsolve = IterativeSolversJL_GMRES(), + QNDF(autodiff = false, linsolve = KrylovJL_GMRES(), concrete_jac = true)); @test length(sol.t) < 25 sol = @test_nowarn solve(prob, Rosenbrock23(autodiff = false, - linsolve = IterativeSolversJL_GMRES(), + linsolve = KrylovJL_GMRES(), precs = precslr, concrete_jac = true)); @test length(sol.t) < 20 sol = @test_nowarn solve(prob, - Rodas4(autodiff = false, linsolve = IterativeSolversJL_GMRES(), + Rodas4(autodiff = false, linsolve = KrylovJL_GMRES(), precs = precslr, concrete_jac = true)); @test length(sol.t) < 20 From 908cd692f457e5e253f59b700e98a2e1e856fbaf Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Wed, 17 May 2023 13:40:06 -0400 Subject: [PATCH 46/54] add more checks to norecompile --- test/interface/norecompile.jl | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/test/interface/norecompile.jl b/test/interface/norecompile.jl index 1443b76265..48f9ff95cc 100644 --- a/test/interface/norecompile.jl +++ b/test/interface/norecompile.jl @@ -14,14 +14,19 @@ function lorenz(du, u, p, t) end lorenzprob = ODEProblem(lorenz, [1.0; 0.0; 0.0], (0.0, 1.0), Float64[]) -t1 = @elapsed sol = solve(lorenzprob, Rosenbrock23()) -t2 = @elapsed sol = solve(lorenzprob, Rosenbrock23(autodiff = false)) +t1 = @elapsed sol1 = solve(lorenzprob, Rosenbrock23()) +t2 = @elapsed sol2 = solve(lorenzprob, Rosenbrock23(autodiff = false)) lorenzprob2 = ODEProblem{true, SciMLBase.FullSpecialize}(lorenz, [1.0; 0.0; 0.0], (0.0, 1.0), Float64[]) -t3 = @elapsed sol = solve(lorenzprob2, Rosenbrock23()) -t4 = @elapsed sol = solve(lorenzprob2, Rosenbrock23(autodiff = false)) +t3 = @elapsed sol3 = solve(lorenzprob2, Rosenbrock23()) +t4 = @elapsed sol4 = solve(lorenzprob2, Rosenbrock23(autodiff = false)) + +@test sol1.retcode === ReturnCode.Success +@test sol2.retcode === ReturnCode.Success +@test sol3.retcode === ReturnCode.Success +@test sol4.retcode === ReturnCode.Success if VERSION >= v"1.8" @test t1 < t3 From c7e778e2a7e55a5aaa06563db307157dc312f8d1 Mon Sep 17 00:00:00 2001 From: Gaurav Arya Date: Thu, 18 May 2023 16:29:12 -0400 Subject: [PATCH 47/54] Fix UJacobianWrapper update coeffs --- src/derivative_utils.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/derivative_utils.jl b/src/derivative_utils.jl index 7ad14cc8df..ff45fca8de 100644 --- a/src/derivative_utils.jl +++ b/src/derivative_utils.jl @@ -301,10 +301,9 @@ function SciMLOperators.update_coefficients!(W::WOperator, u, p, t) W end -function SciMLOperators.update_coefficients!(J::FunctionOperator{UJacobianWrapper}, u, p, t) - copyto!(J.x, u) - J.f.t = t - J.f.p = p +function SciMLOperators.update_coefficients!(J::UJacobianWrapper, u, p, t) + J.p = p + J.t = t end function Base.convert(::Type{AbstractMatrix}, W::WOperator{IIP}) where {IIP} From 83d38d4769017af3d0c08c327d5b63263580c2c0 Mon Sep 17 00:00:00 2001 From: Gaurav Arya Date: Sat, 20 May 2023 22:30:50 -0400 Subject: [PATCH 48/54] Test update coefficients of jacobian operator --- test/algconvergence/ode_rosenbrock_tests.jl | 23 +++++++++++++++++---- 1 file changed, 19 insertions(+), 4 deletions(-) diff --git a/test/algconvergence/ode_rosenbrock_tests.jl b/test/algconvergence/ode_rosenbrock_tests.jl index 9a130c3fb6..daa76047ae 100644 --- a/test/algconvergence/ode_rosenbrock_tests.jl +++ b/test/algconvergence/ode_rosenbrock_tests.jl @@ -1,10 +1,11 @@ +using OrdinaryDiffEq, DiffEqDevTools, Test, LinearAlgebra, LinearSolve +import ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear, + prob_ode_bigfloatlinear, prob_ode_bigfloat2Dlinear +import LinearSolve + @testset "Rosenbrock Tests" begin ## Breakout these since no other test of their adaptivity - using OrdinaryDiffEq, DiffEqDevTools, Test, LinearAlgebra, LinearSolve - import ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear, - prob_ode_bigfloatlinear, prob_ode_bigfloat2Dlinear - dts = (1 / 2) .^ (6:-1:3) testTol = 0.2 @@ -435,3 +436,17 @@ prob = ODEProblem((u, p, t) -> 0.9u, 0.1, (0.0, 1.0)) @test_nowarn solve(prob, Rosenbrock23(autodiff = false)) end + +@testset "Convergence with time-dependent matrix-free Jacobian" + time_derivative(du, u, p, t) = (du[1] = t * u[1]) + time_derivative_analytic(u0, p, t) = u0 * exp(t^2 / 2) + ff_time_derivative = ODEFunction(time_derivative, analytic=time_derivative_analytic) + prob = ODEProblem(ff_time_derivative, [1.0], (0.0, 1.0)) + + dts = (1 / 2) .^ (6:-1:3) + testTol = 0.2 + # Check convergence of Rodas3 with time-dependent matrix-free Jacobian. + # Primarily to check that the Jacobian is being updated correctly as t changes. + sim = test_convergence(dts, prob, Rodas3(linsolve=LinearSolve.KrylovJL())); + @test sim.𝒪est[:final]≈3 atol=testTol +end \ No newline at end of file From 477a58b2ad93cae448cdc372467a272af4d2d541 Mon Sep 17 00:00:00 2001 From: Gaurav Arya Date: Sat, 20 May 2023 23:01:36 -0400 Subject: [PATCH 49/54] Fix typo --- test/algconvergence/ode_rosenbrock_tests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/algconvergence/ode_rosenbrock_tests.jl b/test/algconvergence/ode_rosenbrock_tests.jl index daa76047ae..25393948bc 100644 --- a/test/algconvergence/ode_rosenbrock_tests.jl +++ b/test/algconvergence/ode_rosenbrock_tests.jl @@ -437,7 +437,7 @@ import LinearSolve @test_nowarn solve(prob, Rosenbrock23(autodiff = false)) end -@testset "Convergence with time-dependent matrix-free Jacobian" +@testset "Convergence with time-dependent matrix-free Jacobian" begin time_derivative(du, u, p, t) = (du[1] = t * u[1]) time_derivative_analytic(u0, p, t) = u0 * exp(t^2 / 2) ff_time_derivative = ODEFunction(time_derivative, analytic=time_derivative_analytic) From af44803eeebb5885b8926856c466c89e2c707b24 Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Mon, 22 May 2023 14:40:41 -0400 Subject: [PATCH 50/54] scalar test fixes --- src/OrdinaryDiffEq.jl | 2 +- test/algconvergence/linear_nonlinear_convergence_tests.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index 9c280c2ef0..fa2ba56f97 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -33,7 +33,7 @@ import DiffEqBase: solve!, step!, initialize!, isadaptive import DiffEqBase: ODE_DEFAULT_NORM, ODE_DEFAULT_ISOUTOFDOMAIN, ODE_DEFAULT_PROG_MESSAGE, ODE_DEFAULT_UNSTABLE_CHECK -import SciMLOperators: SciMLOperators, AbstractSciMLOperator, +import SciMLOperators: SciMLOperators, AbstractSciMLOperator, AbstractSciMLScalarOperator, MatrixOperator, FunctionOperator, update_coefficients, update_coefficients!, DEFAULT_UPDATE_FUNC, isconstant diff --git a/test/algconvergence/linear_nonlinear_convergence_tests.jl b/test/algconvergence/linear_nonlinear_convergence_tests.jl index edeb2f0e42..60ffee48c2 100644 --- a/test/algconvergence/linear_nonlinear_convergence_tests.jl +++ b/test/algconvergence/linear_nonlinear_convergence_tests.jl @@ -20,8 +20,8 @@ using OrdinaryDiffEq: alg_order ETDRK4, HochOst4, ETD2, - # KenCarp3, # WOperator doesn't support split problems - # CFNLIRK3, # WOperator doesn't support split problems + KenCarp3, + CFNLIRK3, ] @info "$Alg" sim = test_convergence(dts, prob, Alg()) From e37f0ef7c3aca4a49f919b0a70700dee8541127d Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Mon, 22 May 2023 14:41:14 -0400 Subject: [PATCH 51/54] rm comment in derv_wrappers --- src/derivative_wrappers.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/derivative_wrappers.jl b/src/derivative_wrappers.jl index e36e224810..e63625739e 100644 --- a/src/derivative_wrappers.jl +++ b/src/derivative_wrappers.jl @@ -280,7 +280,7 @@ function build_jac_config(alg, f::F1, uf::F2, du1, uprev, u, tmp, du2, end sparsity, colorvec = sparsity_colorvec(f, u) - # TODO: more generc, do we need this? + if alg_autodiff(alg) isa AutoForwardDiff _chunksize = get_chunksize(alg) === Val(0) ? nothing : get_chunksize(alg) # SparseDiffEq uses different convection... From f6122eba129b7bc4b1250060a4d5363eb01c7c2d Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Tue, 23 May 2023 09:47:09 -0400 Subject: [PATCH 52/54] retriggering CI with dummy commit --- test/algconvergence/linear_nonlinear_convergence_tests.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/test/algconvergence/linear_nonlinear_convergence_tests.jl b/test/algconvergence/linear_nonlinear_convergence_tests.jl index 60ffee48c2..b053431e16 100644 --- a/test/algconvergence/linear_nonlinear_convergence_tests.jl +++ b/test/algconvergence/linear_nonlinear_convergence_tests.jl @@ -23,7 +23,6 @@ using OrdinaryDiffEq: alg_order KenCarp3, CFNLIRK3, ] - @info "$Alg" sim = test_convergence(dts, prob, Alg()) @test sim.𝒪est[:l2]≈alg_order(Alg()) atol=0.2 end From 1da5423e238bb720ad9b229b532c230f54b6038b Mon Sep 17 00:00:00 2001 From: Gaurav Arya Date: Tue, 23 May 2023 21:42:48 -0400 Subject: [PATCH 53/54] Fix forwarddiffs_model --- src/alg_utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/alg_utils.jl b/src/alg_utils.jl index 97b2b495f0..3f76c0e505 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -17,7 +17,7 @@ function SciMLBase.forwarddiffs_model(alg::Union{OrdinaryDiffEqAdaptiveImplicitA DAEAlgorithm, OrdinaryDiffEqImplicitAlgorithm, ExponentialAlgorithm}) - alg_autodiff(alg) + alg_autodiff(alg) isa AutoForwardDiff end SciMLBase.forwarddiffs_model_time(alg::RosenbrockAlgorithm) = true From 67786682fac9e8ed24edde73f0aed9a5b6582346 Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Wed, 24 May 2023 13:08:24 -0400 Subject: [PATCH 54/54] update downstream compat --- test/downstream/Project.toml | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/test/downstream/Project.toml b/test/downstream/Project.toml index 97aead5433..ac5220e2ba 100644 --- a/test/downstream/Project.toml +++ b/test/downstream/Project.toml @@ -4,14 +4,17 @@ DDEProblemLibrary = "f42792ee-6ffc-4e2a-ae83-8ee2f22de800" DelayDiffEq = "bcd4f6db-9728-5f36-b5f7-82caef46ccdb" ODEInterface = "54ca160b-1b9f-5127-a996-1867f4bc2a2c" ODEInterfaceDiffEq = "09606e27-ecf5-54fc-bb29-004bd9f985bf" +OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1" +StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [compat] -CUDA = "3.12" +CUDA = "4.2" DDEProblemLibrary = "0.1" -DelayDiffEq = "5.39" +DelayDiffEq = "5.42" ODEInterface = "0.5" -ODEInterfaceDiffEq = "3.11" -SciMLSensitivity = "7.11" -Zygote = "0.6" \ No newline at end of file +ODEInterfaceDiffEq = "3.13" +SciMLSensitivity = "7.30" +StochasticDiffEq = "6.60.1" +Zygote = "0.6.61"