From 10f8547cf76ac5f4a865299b9e673cbbe9cf6948 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Mon, 22 Jul 2024 11:23:36 -0400 Subject: [PATCH 01/11] fix out of place WOperator --- src/derivative_utils.jl | 37 ++++++++++++++++--------------------- 1 file changed, 16 insertions(+), 21 deletions(-) diff --git a/src/derivative_utils.jl b/src/derivative_utils.jl index 23c4ebe485..749d5e9d84 100644 --- a/src/derivative_utils.jl +++ b/src/derivative_utils.jl @@ -90,13 +90,12 @@ function calc_J(integrator, cache, next_step::Bool = false) J = jacobian(uf, uprev, integrator) end - integrator.stats.njacs += 1 - if alg isa CompositeAlgorithm integrator.eigen_est = constvalue(opnorm(J, Inf)) end end + integrator.stats.njacs += 1 J end @@ -144,12 +143,11 @@ function calc_J!(J, integrator, cache, next_step::Bool = false) end end - integrator.stats.njacs += 1 - if alg isa CompositeAlgorithm integrator.eigen_est = constvalue(opnorm(J, Inf)) end + integrator.stats.njacs += 1 return nothing end @@ -604,21 +602,21 @@ function jacobian2W!(W::Matrix, mass_matrix, dtgamma::Number, J::Matrix, return nothing end -function jacobian2W(mass_matrix::MT, dtgamma::Number, J::AbstractMatrix, - W_transform::Bool)::Nothing where {MT} +function jacobian2W(mass_matrix, dtgamma::Number, J::AbstractMatrix, + W_transform::Bool) # check size and dimension mass_matrix isa UniformScaling || @boundscheck axes(mass_matrix) == axes(J) || _throwJMerror(J, mass_matrix) @inbounds if W_transform invdtgamma = inv(dtgamma) - if MT <: UniformScaling + if mass_matrix isa UniformScaling λ = -mass_matrix.λ W = J + (λ * invdtgamma) * I else W = muladd(-mass_matrix, invdtgamma, J) end else - if MT <: UniformScaling + if mass_matrix isa UniformScaling λ = -mass_matrix.λ W = dtgamma * J + λ * I else @@ -742,15 +740,16 @@ end W = cache.W if isnewton(nlsolver) # we will call `update_coefficients!` for u/p/t in NLNewton - update_coefficients!(W; transform = W_transform, dtgamma) + W = update_coefficients(W; transform = W_transform, dtgamma) else - update_coefficients!(W, uprev, p, t; transform = W_transform, dtgamma) + W = update_coefficients(W, uprev, p, t; transform = W_transform, dtgamma) end - if W.J !== nothing && !(W.J isa AbstractSciMLOperator) - islin, isode = islinearfunction(integrator) + if W.J !== nothing J = islin ? (isode ? f.f : f.f1.f) : calc_J(integrator, cache, next_step) - !isdae && - jacobian2W!(W._concrete_form, mass_matrix, dtgamma, J, W_transform) + if !isdae + integrator.stats.nw += 1 + W = jacobian2W(mass_matrix, dtgamma, J, W_transform) + end end elseif cache.W isa AbstractSciMLOperator && !(cache.W isa StaticWOperator) J = update_coefficients(cache.J, uprev, p, t) @@ -760,11 +759,9 @@ end W = WOperator{false}(mass_matrix, dtgamma, J, uprev; transform = W_transform) elseif DiffEqBase.has_jac(f) J = f.jac(uprev, p, t) + integrator.stats.njacs += 1 if J isa StaticArray && - integrator.alg isa - Union{ - Rosenbrock23, Rodas23W, Rodas3P, Rodas4, Rodas4P, Rodas4P2, - Rodas5, Rodas5P, Rodas5Pe, Rodas5Pr} + integrator.alg isa OrdinaryDiffEqRosenbrockAdaptiveAlgorithm W = W_transform ? J - mass_matrix * inv(dtgamma) : dtgamma * J - mass_matrix else @@ -788,9 +785,7 @@ end W = if W_full isa Number W_full elseif len !== nothing && - integrator.alg isa - Union{Rosenbrock23, Rodas23W, Rodas3P, Rodas4, Rodas4P, - Rodas4P2, Rodas5, Rodas5P, Rodas5Pe, Rodas5Pr} + integrator.alg isa OrdinaryDiffEqRosenbrockAdaptiveAlgorithm StaticWOperator(W_full) else DiffEqBase.default_factorize(W_full) From c76ead00559e12819d1eae936a205708c7d3c6d0 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Mon, 22 Jul 2024 11:38:56 -0400 Subject: [PATCH 02/11] try the ! version of update_coefficients --- 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 749d5e9d84..9366f22e4f 100644 --- a/src/derivative_utils.jl +++ b/src/derivative_utils.jl @@ -740,9 +740,9 @@ end W = cache.W if isnewton(nlsolver) # we will call `update_coefficients!` for u/p/t in NLNewton - W = update_coefficients(W; transform = W_transform, dtgamma) + W = update_coefficients!(W; transform = W_transform, dtgamma) else - W = update_coefficients(W, uprev, p, t; transform = W_transform, dtgamma) + W = update_coefficients!(W, uprev, p, t; transform = W_transform, dtgamma) end if W.J !== nothing J = islin ? (isode ? f.f : f.f1.f) : calc_J(integrator, cache, next_step) From 67d3a4e078a36da61249463cbaa3ec84d1ccd080 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Mon, 22 Jul 2024 12:33:51 -0400 Subject: [PATCH 03/11] don't use WOperator for Array in OOP --- src/derivative_utils.jl | 53 +++++++++++++++++++++-------------------- 1 file changed, 27 insertions(+), 26 deletions(-) diff --git a/src/derivative_utils.jl b/src/derivative_utils.jl index 9366f22e4f..dcb6e52ce2 100644 --- a/src/derivative_utils.jl +++ b/src/derivative_utils.jl @@ -744,7 +744,7 @@ end else W = update_coefficients!(W, uprev, p, t; transform = W_transform, dtgamma) end - if W.J !== nothing + if W.J !== nothing && !(W.J isa AbstractSciMLOperator) J = islin ? (isode ? f.f : f.f1.f) : calc_J(integrator, cache, next_step) if !isdae integrator.stats.nw += 1 @@ -757,27 +757,15 @@ end elseif islin J = isode ? f.f : f.f1.f # unwrap the Jacobian accordingly W = WOperator{false}(mass_matrix, dtgamma, J, uprev; transform = W_transform) - elseif DiffEqBase.has_jac(f) - J = f.jac(uprev, p, t) - integrator.stats.njacs += 1 - if J isa StaticArray && - integrator.alg isa OrdinaryDiffEqRosenbrockAdaptiveAlgorithm - W = W_transform ? J - mass_matrix * inv(dtgamma) : - dtgamma * J - mass_matrix - else - if !isa(J, AbstractSciMLOperator) && (!isnewton(nlsolver) || - nlsolver.cache.W.J isa AbstractSciMLOperator) - J = MatrixOperator(J) - end - W = WOperator{false}(mass_matrix, dtgamma, J, uprev, cache.W.jacvec; - transform = W_transform) - end - integrator.stats.nw += 1 else integrator.stats.nw += 1 J = calc_J(integrator, cache, next_step) if isdae W = J + elseif J isa StaticArray && + integrator.alg isa OrdinaryDiffEqRosenbrockAdaptiveAlgorithm + W = W_transform ? J - mass_matrix * inv(dtgamma) : + dtgamma * J - mass_matrix else W_full = W_transform ? J - mass_matrix * inv(dtgamma) : dtgamma * J - mass_matrix @@ -912,29 +900,42 @@ function build_J_W(alg, u, uprev, p, t, dt, f::F, ::Type{uEltypeNoUnits}, autodiff = alg_autodiff(alg), tag = OrdinaryDiffEqTag()) W = WOperator{IIP}(f.mass_matrix, dt, J, u, jacvec) - elseif islin || (!IIP && DiffEqBase.has_jac(f)) + elseif islin + J = isode ? f.f : f.f1.f # unwrap the Jacobian accordingly + W = WOperator{IIP}(f.mass_matrix, dt, J, u) + elseif 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 = MatrixOperator(J) + if J isa AbstractSciMLOperator + W = WOperator{IIP}(f.mass_matrix, dt, J, u) + else + W = if alg isa DAEAlgorithm + J + elseif IIP + similar(J) + else + len = StaticArrayInterface.known_length(typeof(J)) + if len !== nothing && + alg isa OrdinaryDiffEqRosenbrockAdaptiveAlgorithm + StaticWOperator(J, false) + else + ArrayInterface.lu_instance(J) + end + end end - W = WOperator{IIP}(f.mass_matrix, dt, J, u) else J = if f.jac_prototype === nothing ArrayInterface.undefmatrix(u) else deepcopy(f.jac_prototype) end - isdae = alg isa DAEAlgorithm - W = if isdae + W = if alg isa DAEAlgorithm J elseif IIP similar(J) else len = StaticArrayInterface.known_length(typeof(J)) if len !== nothing && - alg isa - Union{Rosenbrock23, Rodas23W, Rodas3P, Rodas4, Rodas4P, - Rodas4P2, Rodas5, Rodas5P, Rodas5Pe, Rodas5Pr} + alg isa OrdinaryDiffEqRosenbrockAdaptiveAlgorithm StaticWOperator(J, false) else ArrayInterface.lu_instance(J) From 044542983b1830afcda8aa1d10f637413095878a Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Tue, 23 Jul 2024 17:28:21 -0400 Subject: [PATCH 04/11] still broken --- src/derivative_utils.jl | 77 ++++++++++++----------------------------- 1 file changed, 23 insertions(+), 54 deletions(-) diff --git a/src/derivative_utils.jl b/src/derivative_utils.jl index dcb6e52ce2..a8bfbb6d1c 100644 --- a/src/derivative_utils.jl +++ b/src/derivative_utils.jl @@ -736,45 +736,35 @@ end islin, isode = islinearfunction(integrator) !isdae && update_coefficients!(mass_matrix, uprev, p, t) - if cache.W isa WOperator - W = cache.W - if isnewton(nlsolver) - # we will call `update_coefficients!` for u/p/t in NLNewton - W = update_coefficients!(W; transform = W_transform, dtgamma) - else - W = update_coefficients!(W, uprev, p, t; transform = W_transform, dtgamma) - end - if W.J !== nothing && !(W.J isa AbstractSciMLOperator) - J = islin ? (isode ? f.f : f.f1.f) : calc_J(integrator, cache, next_step) - if !isdae - integrator.stats.nw += 1 - W = jacobian2W(mass_matrix, dtgamma, J, W_transform) - end + if cache.W isa StaticWOperator + integrator.stats.nw += 1 + J = calc_J(integrator, cache, next_step) + W = StaticWOperator(W_transform ? J - mass_matrix * inv(dtgamma) : dtgamma * J - mass_matrix) + elseif cache.W isa WOperator + integrator.stats.nw += 1 + J = if islin + isode ? f.f : f.f1.f + elseif cache.W.J === nothing + calc_J(integrator, cache, next_step) + else + J = update_coefficients(cache.W.J, uprev, p, t) end - elseif cache.W isa AbstractSciMLOperator && !(cache.W isa StaticWOperator) - J = update_coefficients(cache.J, uprev, p, t) + W = WOperator{false}(mass_matrix, dtgamma, J, uprev, cache.W.jacvec; transform = W_transform) + elseif cache.W isa AbstractSciMLOperator W = update_coefficients(cache.W, uprev, p, t; dtgamma, transform = W_transform) elseif islin - J = isode ? f.f : f.f1.f # unwrap the Jacobian accordingly + J = isode ? f.f : f.f1.f W = WOperator{false}(mass_matrix, dtgamma, J, uprev; transform = W_transform) else integrator.stats.nw += 1 - J = calc_J(integrator, cache, next_step) + J = islin ? isode ? f.f : f.f1.f : calc_J(integrator, cache, next_step) if isdae W = J - elseif J isa StaticArray && - integrator.alg isa OrdinaryDiffEqRosenbrockAdaptiveAlgorithm - W = W_transform ? J - mass_matrix * inv(dtgamma) : - dtgamma * J - mass_matrix else W_full = W_transform ? J - mass_matrix * inv(dtgamma) : dtgamma * J - mass_matrix - len = StaticArrayInterface.known_length(typeof(W_full)) W = if W_full isa Number W_full - elseif len !== nothing && - integrator.alg isa OrdinaryDiffEqRosenbrockAdaptiveAlgorithm - StaticWOperator(W_full) else DiffEqBase.default_factorize(W_full) end @@ -859,10 +849,11 @@ function build_J_W(alg, u, uprev, p, t, dt, f::F, ::Type{uEltypeNoUnits}, elseif f.jac_prototype isa AbstractSciMLOperator W = WOperator{IIP}(f, u, dt) J = W.J + elseif islin + J = isode ? f.f : f.f1.f # unwrap the Jacobian accordingly + W = WOperator{IIP}(f.mass_matrix, dt, J, u) elseif IIP && f.jac_prototype !== nothing && concrete_jac(alg) === nothing && - (alg.linsolve === nothing || - alg.linsolve !== nothing && - LinearSolve.needs_concrete_A(alg.linsolve)) + (alg.linsolve === nothing || LinearSolve.needs_concrete_A(alg.linsolve)) # If factorization, then just use the jac_prototype J = similar(f.jac_prototype) @@ -879,7 +870,6 @@ function build_J_W(alg, u, uprev, p, t, dt, f::F, ::Type{uEltypeNoUnits}, autodiff = alg_autodiff(alg), tag = OrdinaryDiffEqTag()) J = jacvec W = WOperator{IIP}(f.mass_matrix, dt, J, u, jacvec) - elseif alg.linsolve !== nothing && !LinearSolve.needs_concrete_A(alg.linsolve) || concrete_jac(alg) !== nothing && concrete_jac(alg) # The linear solver does not need a concrete Jacobian, but the user has @@ -899,31 +889,10 @@ function build_J_W(alg, u, uprev, p, t, dt, f::F, ::Type{uEltypeNoUnits}, jacvec = JacVec(__f, copy(u), p, t; autodiff = alg_autodiff(alg), tag = OrdinaryDiffEqTag()) W = WOperator{IIP}(f.mass_matrix, dt, J, u, jacvec) - - elseif islin - J = isode ? f.f : f.f1.f # unwrap the Jacobian accordingly - W = WOperator{IIP}(f.mass_matrix, dt, J, u) - elseif DiffEqBase.has_jac(f) - J = islin ? (isode ? f.f : f.f1.f) : f.jac(uprev, p, t) # unwrap the Jacobian accordingly - if J isa AbstractSciMLOperator - W = WOperator{IIP}(f.mass_matrix, dt, J, u) - else - W = if alg isa DAEAlgorithm - J - elseif IIP - similar(J) - else - len = StaticArrayInterface.known_length(typeof(J)) - if len !== nothing && - alg isa OrdinaryDiffEqRosenbrockAdaptiveAlgorithm - StaticWOperator(J, false) - else - ArrayInterface.lu_instance(J) - end - end - end else - J = if f.jac_prototype === nothing + J = if !IIP && DiffEqBase.has_jac(f) + f.jac(uprev, p, t) + elseif f.jac_prototype === nothing ArrayInterface.undefmatrix(u) else deepcopy(f.jac_prototype) From 04733643e989c67c6081093eb69ccd0025766154 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Wed, 24 Jul 2024 11:50:36 -0400 Subject: [PATCH 05/11] improve static --- src/derivative_utils.jl | 33 ++++++++++++++++++--------------- 1 file changed, 18 insertions(+), 15 deletions(-) diff --git a/src/derivative_utils.jl b/src/derivative_utils.jl index a8bfbb6d1c..4dfdae3ba9 100644 --- a/src/derivative_utils.jl +++ b/src/derivative_utils.jl @@ -739,6 +739,7 @@ end if cache.W isa StaticWOperator integrator.stats.nw += 1 J = calc_J(integrator, cache, next_step) + @assert J isa StaticArray W = StaticWOperator(W_transform ? J - mass_matrix * inv(dtgamma) : dtgamma * J - mass_matrix) elseif cache.W isa WOperator integrator.stats.nw += 1 @@ -770,8 +771,6 @@ end end end end - (W isa WOperator && unwrap_alg(integrator, true) isa NewtonAlgorithm) && - (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 @@ -833,6 +832,7 @@ function update_W!(nlsolver::AbstractNLSolver, nothing end +import StaticArrays: StaticMatrix function build_J_W(alg, u, uprev, p, t, dt, f::F, ::Type{uEltypeNoUnits}, ::Val{IIP}) where {IIP, uEltypeNoUnits, F} # TODO - make J, W AbstractSciMLOperators (lazily defined with scimlops functionality) @@ -881,14 +881,20 @@ function build_J_W(alg, u, uprev, p, t, dt, f::F, ::Type{uEltypeNoUnits}, else deepcopy(f.jac_prototype) end - __f = if IIP - (du, u, p, t) -> _f(du, u, p, t) + W = if J isa StaticMatrix && alg isa OrdinaryDiffEqRosenbrockAdaptiveAlgorithm + StaticWOperator(J, false) + elseif J isa StaticMatrix + ArrayInterface.lu_instance(J) else - (u, p, t) -> _f(u, p, t) + __f = if IIP + (du, u, p, t) -> _f(du, u, p, t) + else + (u, p, t) -> _f(u, p, t) + end + jacvec = JacVec(__f, copy(u), p, t; + autodiff = alg_autodiff(alg), tag = OrdinaryDiffEqTag()) + WOperator{IIP}(f.mass_matrix, dt, J, u, jacvec) end - jacvec = JacVec(__f, copy(u), p, t; - autodiff = alg_autodiff(alg), tag = OrdinaryDiffEqTag()) - W = WOperator{IIP}(f.mass_matrix, dt, J, u, jacvec) else J = if !IIP && DiffEqBase.has_jac(f) f.jac(uprev, p, t) @@ -901,16 +907,13 @@ function build_J_W(alg, u, uprev, p, t, dt, f::F, ::Type{uEltypeNoUnits}, J elseif IIP similar(J) + elseif J isa StaticMatrix && alg isa OrdinaryDiffEqRosenbrockAdaptiveAlgorithm + StaticWOperator(J, false) else - len = StaticArrayInterface.known_length(typeof(J)) - if len !== nothing && - alg isa OrdinaryDiffEqRosenbrockAdaptiveAlgorithm - StaticWOperator(J, false) - else - ArrayInterface.lu_instance(J) - end + ArrayInterface.lu_instance(J) end end + @show W return J, W end From 64d8bae40dad185d9b2ecd6ece7af4188390f7fe Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Wed, 24 Jul 2024 12:26:56 -0400 Subject: [PATCH 06/11] it works --- src/derivative_utils.jl | 21 ++++++--------------- 1 file changed, 6 insertions(+), 15 deletions(-) diff --git a/src/derivative_utils.jl b/src/derivative_utils.jl index 4dfdae3ba9..5613ba896a 100644 --- a/src/derivative_utils.jl +++ b/src/derivative_utils.jl @@ -739,35 +739,27 @@ end if cache.W isa StaticWOperator integrator.stats.nw += 1 J = calc_J(integrator, cache, next_step) - @assert J isa StaticArray W = StaticWOperator(W_transform ? J - mass_matrix * inv(dtgamma) : dtgamma * J - mass_matrix) elseif cache.W isa WOperator integrator.stats.nw += 1 J = if islin - isode ? f.f : f.f1.f - elseif cache.W.J === nothing - calc_J(integrator, cache, next_step) - else - J = update_coefficients(cache.W.J, uprev, p, t) + isode ? f.f : f.f1.f + else + calc_J(integrator, cache, next_step) end W = WOperator{false}(mass_matrix, dtgamma, J, uprev, cache.W.jacvec; transform = W_transform) elseif cache.W isa AbstractSciMLOperator W = update_coefficients(cache.W, uprev, p, t; dtgamma, transform = W_transform) - elseif islin - J = isode ? f.f : f.f1.f - W = WOperator{false}(mass_matrix, dtgamma, J, uprev; transform = W_transform) else integrator.stats.nw += 1 J = islin ? isode ? f.f : f.f1.f : calc_J(integrator, cache, next_step) if isdae W = J else - W_full = W_transform ? J - mass_matrix * inv(dtgamma) : + W = W_transform ? J - mass_matrix * inv(dtgamma) : dtgamma * J - mass_matrix - W = if W_full isa Number - W_full - else - DiffEqBase.default_factorize(W_full) + if !isa(W, Number) + W = DiffEqBase.default_factorize(W) end end end @@ -913,7 +905,6 @@ function build_J_W(alg, u, uprev, p, t, dt, f::F, ::Type{uEltypeNoUnits}, ArrayInterface.lu_instance(J) end end - @show W return J, W end From d421621659ea3fd29137d8734d703c5470fab838 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Wed, 24 Jul 2024 15:53:29 -0400 Subject: [PATCH 07/11] add ParameterizedFunctions to test deps --- Project.toml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 57da92e269..c6d21868f0 100644 --- a/Project.toml +++ b/Project.toml @@ -101,6 +101,7 @@ InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" ODEProblemLibrary = "fdc4e326-1af4-4b90-96e7-779fcce2daa5" +ParameterizedFunctions = "65888b18-ceab-5e60-b2b9-181511a3b968" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" PoissonRandom = "e409e4f3-bfea-5376-8464-e040bb5c01ab" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" @@ -114,4 +115,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", "NLsolve"] +test = ["Calculus", "ComponentArrays", "Symbolics", "AlgebraicMultigrid", "IncompleteLU", "DiffEqCallbacks", "DiffEqDevTools", "ODEProblemLibrary", "ElasticArrays", "InteractiveUtils", "ParameterizedFunctions", "PoissonRandom", "Printf", "Random", "ReverseDiff", "SafeTestsets", "SparseArrays", "Statistics", "Test", "Unitful", "ModelingToolkit", "Pkg", "NLsolve"] From 0772416939300f475b553a621ff9eb88c803224f Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Thu, 25 Jul 2024 09:41:30 -0400 Subject: [PATCH 08/11] add Hires test --- test/interface/linear_solver_test.jl | 40 ++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) diff --git a/test/interface/linear_solver_test.jl b/test/interface/linear_solver_test.jl index 5fb1124665..ae759b14f1 100644 --- a/test/interface/linear_solver_test.jl +++ b/test/interface/linear_solver_test.jl @@ -157,3 +157,43 @@ end atol = 1e-1, rtol = 1e-1) @test isapprox(exp.(p), g_helper(p; alg = KenCarp47(linsolve = KrylovJL_GMRES())); atol = 1e-1, rtol = 1e-1) + +using OrdinaryDiffEq, StaticArrays, LinearSolve, ParameterizedFunctions + +hires = @ode_def Hires begin + dy1 = -1.71*y1 + 0.43*y2 + 8.32*y3 + 0.0007 + dy2 = 1.71*y1 - 8.75*y2 + dy3 = -10.03*y3 + 0.43*y4 + 0.035*y5 + dy4 = 8.32*y2 + 1.71*y3 - 1.12*y4 + dy5 = -1.745*y5 + 0.43*y6 + 0.43*y7 + dy6 = -280.0*y6*y8 + 0.69*y4 + 1.71*y5 - 0.43*y6 + 0.69*y7 + dy7 = 280.0*y6*y8 - 1.81*y7 + dy8 = -280.0*y6*y8 + 1.81*y7 +end + +u0 = zeros(8) +u0[1] = 1 +u0[8] = 0.0057 + +probiip = ODEProblem{true}(hires, u0, (0.0,10.0)) +proboop = ODEProblem{false}(hires, u0, (0.0,10.0)) +probstatic = ODEProblem{false}(hires, SVector{8}(u0), (0.0,10.0)) +probs = (;probiip, proboop, probstatic) +qndf = QNDF() +krylov_qndf = QNDF(linsolve=KrylovJL_GMRES()) +fbdf = FBDF() +krylov_fbdf = FBDF(linsolve=KrylovJL_GMRES()) +rodas = Rodas5P() +krylov_rodas = Rodas5P(linsolve=KrylovJL_GMRES()) +solvers = (;qndf, krylov_qndf, rodas, krylov_rodas, fbdf, krylov_fbdf, ) + +refsol = solve(probiip, FBDF(), abstol=1e-12, reltol=1e-12) +@testset "Hires" begin +@testset "$probname" for (probname, prob) in pairs(probs) + @testset "$solname" for (solname, solver) in pairs(solvers) + sol = solve(prob, solver, abstol=1e-12, reltol=1e-12, maxiters=2e4) + @test sol.retcode == ReturnCode.Success + @test isapprox(sol.u[end], refsol.u[end], rtol=1e-10, atol=1e-10) + end +end +end From 8dbca6eeeba82873f1f6fe1be3b4ae2cef51fc56 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Thu, 25 Jul 2024 11:05:48 -0400 Subject: [PATCH 09/11] loosen test --- test/interface/linear_solver_test.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/test/interface/linear_solver_test.jl b/test/interface/linear_solver_test.jl index ae759b14f1..32e8770020 100644 --- a/test/interface/linear_solver_test.jl +++ b/test/interface/linear_solver_test.jl @@ -188,12 +188,12 @@ krylov_rodas = Rodas5P(linsolve=KrylovJL_GMRES()) solvers = (;qndf, krylov_qndf, rodas, krylov_rodas, fbdf, krylov_fbdf, ) refsol = solve(probiip, FBDF(), abstol=1e-12, reltol=1e-12) -@testset "Hires" begin -@testset "$probname" for (probname, prob) in pairs(probs) - @testset "$solname" for (solname, solver) in pairs(solvers) - sol = solve(prob, solver, abstol=1e-12, reltol=1e-12, maxiters=2e4) - @test sol.retcode == ReturnCode.Success - @test isapprox(sol.u[end], refsol.u[end], rtol=1e-10, atol=1e-10) +@testset "Hires calc_W tests" begin + @testset "$probname" for (probname, prob) in pairs(probs) + @testset "$solname" for (solname, solver) in pairs(solvers) + sol = solve(prob, solver, abstol=1e-12, reltol=1e-12, maxiters=2e4) + @test sol.retcode == ReturnCode.Success + @test isapprox(sol.u[end], refsol.u[end], rtol=1e-8, atol=1e-10) + end end end -end From 7e69639b7d7f393d3d83233a51814a106363c20d Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Fri, 26 Jul 2024 08:48:33 -0400 Subject: [PATCH 10/11] Update derivative_utils.jl --- src/derivative_utils.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/derivative_utils.jl b/src/derivative_utils.jl index 5613ba896a..1f23127bf2 100644 --- a/src/derivative_utils.jl +++ b/src/derivative_utils.jl @@ -824,7 +824,6 @@ function update_W!(nlsolver::AbstractNLSolver, nothing end -import StaticArrays: StaticMatrix function build_J_W(alg, u, uprev, p, t, dt, f::F, ::Type{uEltypeNoUnits}, ::Val{IIP}) where {IIP, uEltypeNoUnits, F} # TODO - make J, W AbstractSciMLOperators (lazily defined with scimlops functionality) From 9c667c4d5f0e8f7c3d034a6f1e26f0f1775d73d3 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Fri, 26 Jul 2024 08:49:14 -0400 Subject: [PATCH 11/11] Update OrdinaryDiffEq.jl --- src/OrdinaryDiffEq.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index 6214e87808..67cd3b2b4a 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -66,7 +66,7 @@ using NonlinearSolve # Required by temporary fix in not in-place methods with 12+ broadcasts # `MVector` is used by Nordsieck forms -import StaticArrays: SArray, MVector, SVector, @SVector, StaticArray, MMatrix, SA +import StaticArrays: SArray, MVector, SVector, @SVector, StaticArray, MMatrix, SA, StaticMatrix # Integrator Interface import DiffEqBase: resize!, deleteat!, addat!, full_cache, user_cache, u_cache, du_cache,