From a9847685255841a4198b0839470882c73e64a302 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Sun, 23 Mar 2025 15:55:07 +0100 Subject: [PATCH 01/26] allow LinearSolve v3 --- Project.toml | 2 +- docs/Project.toml | 2 +- src/mprk.jl | 4 +++- test/Project.toml | 4 +++- test/runtests.jl | 28 ++++++++++++++++++---------- 5 files changed, 26 insertions(+), 14 deletions(-) diff --git a/Project.toml b/Project.toml index 66a17068..200b85d8 100644 --- a/Project.toml +++ b/Project.toml @@ -21,7 +21,7 @@ SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" [compat] FastBroadcast = "0.3.5" LinearAlgebra = "1" -LinearSolve = "2.32" +LinearSolve = "2.32, 3" MuladdMacro = "0.2.4" OrdinaryDiffEqCore = "1.16" RecipesBase = "1.3.4" diff --git a/docs/Project.toml b/docs/Project.toml index bf7d0553..e17a197e 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -27,7 +27,7 @@ DiffEqDevTools = "2.45.1" Documenter = "1" InteractiveUtils = "1" LinearAlgebra = "1" -LinearSolve = "2.32" +LinearSolve = "2.32, 3" OrdinaryDiffEqFIRK = "1.7" OrdinaryDiffEqLowOrderRK = "1.2" OrdinaryDiffEqRosenbrock = "1.4" diff --git a/src/mprk.jl b/src/mprk.jl index cd369eea..07a1cdbc 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -429,7 +429,9 @@ end build_mprk_matrix!(P, P, σ, dt) # Same as linres = P \ uprev + display(linsolve.A) linsolve.A = P + display(linsolve.A) linres = solve!(linsolve) u .= linres @@ -987,7 +989,7 @@ You can optionally choose the linear solver to be used by passing an algorithm from [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl) as keyword argument `linsolve`. You can also choose the parameter `small_constant` which is added to all Patankar-weight denominators -to avoid divisions by zero. To display the default value for data type `type` evaluate +to avoid divisions by zero. To display the default value for data type `type` evaluate `MPRK43II(gamma).small_constant_function(type)`, where `type` can be, e.g., `Float64`. diff --git a/test/Project.toml b/test/Project.toml index 464f1e9f..978b479d 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -11,6 +11,7 @@ OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" OrdinaryDiffEqVerner = "79d7bb75-1356-48c1-b8c0-6832512096c2" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" +RecursiveFactorization = "f2c3362d-daeb-58d1-803e-2bc74f2840b4" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" @@ -22,7 +23,7 @@ ADTypes = "1.12" Aqua = "0.8" ExplicitImports = "1.0.1" LinearAlgebra = "1" -LinearSolve = "2.32" +LinearSolve = "2.32, 3" OrdinaryDiffEqLowOrderRK = "1.2" OrdinaryDiffEqRosenbrock = "1.4" OrdinaryDiffEqSDIRK = "1.2" @@ -30,6 +31,7 @@ OrdinaryDiffEqTsit5 = "1.1" OrdinaryDiffEqVerner = "1.1" Plots = "1.25.11" RecipesBase = "1.3.4" +RecursiveFactorization = "0.2" SparseArrays = "1" StaticArrays = "1.9.7" Statistics = "1" diff --git a/test/runtests.jl b/test/runtests.jl index 02706689..754bf6a6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -15,6 +15,8 @@ using OrdinaryDiffEqTsit5: Tsit5 using OrdinaryDiffEqVerner: Vern7, Vern9 using PositiveIntegrators +# load RecursiveFactorization to get RFLUFactorization +using RecursiveFactorization: RecursiveFactorization using LinearSolve: RFLUFactorization, LUFactorization, KrylovJL_GMRES using Aqua: Aqua @@ -1265,10 +1267,16 @@ end # deliver the same results @testset "Different matrix types (conservative)" begin prod_1! = (P, u, p, t) -> begin - fill!(P, zero(eltype(P))) - for i in 1:(length(u) - 1) - P[i, i + 1] = i * u[i] + # fill!(P, zero(eltype(P))) + for j in axes(P, 2) + for idx in nzrange(P, j) + i = rowvals(P)[idx] + nonzeros(P)[idx] = i * u[i] + end end + # for i in 1:(length(u) - 1) + # P[i, i + 1] = i * u[i] + # end return nothing end @@ -2274,7 +2282,7 @@ end end # Here we run a single scheme multiple times and check - # that the errors are identical and that the computing times + # that the errors are identical and that the computing times # differ only slightly. @testset "work-precision fixed" begin prob = prob_pds_nonlinmod @@ -2295,7 +2303,7 @@ end @test all(y -> y == v[1], v) end - # check that computing times are close enough + # check that computing times are close enough for (i, _) in enumerate(dts) v = [value[i][2] for (key, value) in wp] m1 = mean(v) @@ -2307,7 +2315,7 @@ end end # Here we run a single scheme multiple times and check - # that the errors are identical and that the computing times + # that the errors are identical and that the computing times # differ only slightly. @testset "work-precision adaptive" begin @testset "adatpive_ref = false" begin @@ -2332,13 +2340,13 @@ end @test all(y -> y == v[1], v) end - # check that computing times are close enough + # check that computing times are close enough for (i, _) in enumerate(abstols) v = [value[i][2] for (key, value) in wp] m1 = mean(v) # This test allows computing times that are # 2.5 times the mean value. In a loglog plot these - # differences won't be significant. + # differences won't be significant. @test maximum((v .- m1) ./ m1) < 1.5 end end @@ -2365,13 +2373,13 @@ end @test all(y -> y == v[1], v) end - # check that computing times are close enough + # check that computing times are close enough for (i, _) in enumerate(abstols) v = [value[i][2] for (key, value) in wp] m1 = mean(v) # This test allows computing times that are # 2.5 times the mean value. In a loglog plot these - # differences won't be significant. + # differences won't be significant. @test maximum((v .- m1) ./ m1) < 1.5 end end From ce17d63e636db0eb2d2bfc88200106f397f0376d Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Mon, 24 Mar 2025 11:00:26 +0100 Subject: [PATCH 02/26] use KLUFactorization() instead of LUFactorization() --- src/PositiveIntegrators.jl | 2 +- src/mprk.jl | 8 ++++---- src/sspmprk.jl | 6 +++--- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/PositiveIntegrators.jl b/src/PositiveIntegrators.jl index 51306dec..ba0ec4cd 100644 --- a/src/PositiveIntegrators.jl +++ b/src/PositiveIntegrators.jl @@ -24,7 +24,7 @@ using OrdinaryDiffEqCore: OrdinaryDiffEqCore, OrdinaryDiffEqAlgorithm, ODESoluti using SymbolicIndexingInterface: SymbolicIndexingInterface -using LinearSolve: LinearSolve, LinearProblem, LUFactorization, solve! +using LinearSolve: LinearSolve, LinearProblem, KLUFactorization, solve! import SciMLBase: interp_summary diff --git a/src/mprk.jl b/src/mprk.jl index 07a1cdbc..5d002ffd 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -260,7 +260,7 @@ struct MPE{F, T} <: OrdinaryDiffEqAlgorithm small_constant_function::T end -function MPE(; linsolve = LUFactorization(), small_constant = nothing) +function MPE(; linsolve = KLUFactorization(), small_constant = nothing) if isnothing(small_constant) small_constant_function = floatmin elseif small_constant isa Number @@ -494,7 +494,7 @@ struct MPRK22{T, F, T2} <: OrdinaryDiffEqAdaptiveAlgorithm small_constant_function::T2 end -function MPRK22(alpha; linsolve = LUFactorization(), +function MPRK22(alpha; linsolve = KLUFactorization(), small_constant = nothing) if isnothing(small_constant) small_constant_function = floatmin @@ -904,7 +904,7 @@ struct MPRK43I{T, F, T2} <: OrdinaryDiffEqAdaptiveAlgorithm small_constant_function::T2 end -function MPRK43I(alpha, beta; linsolve = LUFactorization(), +function MPRK43I(alpha, beta; linsolve = KLUFactorization(), small_constant = nothing) if isnothing(small_constant) small_constant_function = floatmin @@ -1016,7 +1016,7 @@ function small_constant_function_MPRK43II(type) return small_constant end -function MPRK43II(gamma; linsolve = LUFactorization(), +function MPRK43II(gamma; linsolve = KLUFactorization(), small_constant = small_constant_function_MPRK43II) if small_constant isa Number small_constant_function = Returns(small_constant) diff --git a/src/sspmprk.jl b/src/sspmprk.jl index 447eb0d4..cc087225 100644 --- a/src/sspmprk.jl +++ b/src/sspmprk.jl @@ -48,7 +48,7 @@ struct SSPMPRK22{T, F, T2} <: OrdinaryDiffEqAdaptiveAlgorithm small_constant_function::T2 end -function SSPMPRK22(alpha, beta; linsolve = LUFactorization(), +function SSPMPRK22(alpha, beta; linsolve = KLUFactorization(), small_constant = nothing) if isnothing(small_constant) small_constant_function = floatmin @@ -466,7 +466,7 @@ You can optionally choose the linear solver to be used by passing an algorithm from [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl) as keyword argument `linsolve`. You can also choose the parameter `small_constant` which is added to all Patankar-weight denominators -to avoid divisions by zero. To display the default value for data type `type` evaluate +to avoid divisions by zero. To display the default value for data type `type` evaluate `SSPMPRK43. small_constant_function(type)`, where `type` can be, e.g., `Float64`. @@ -506,7 +506,7 @@ function small_constant_function_SSPMPRK43(type) return small_constant end -function SSPMPRK43(; linsolve = LUFactorization(), +function SSPMPRK43(; linsolve = KLUFactorization(), small_constant = small_constant_function_SSPMPRK43) if small_constant isa Number small_constant_function = Returns(small_constant) From 04aaeb050929ee4d5bc8e9ea257892a4cb5cf861 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Mon, 24 Mar 2025 11:01:12 +0100 Subject: [PATCH 03/26] remove bandaid --- .github/workflows/Downgrade.yml | 2 +- Project.toml | 3 +-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/.github/workflows/Downgrade.yml b/.github/workflows/Downgrade.yml index 8e5cb87c..d5422872 100644 --- a/.github/workflows/Downgrade.yml +++ b/.github/workflows/Downgrade.yml @@ -38,7 +38,7 @@ jobs: - uses: julia-actions/cache@v2 - uses: julia-actions/julia-downgrade-compat@v1 with: - skip: LinearAlgebra,SparseArrays,Statistics,Test,SciMLBase,SciMLOperators + skip: LinearAlgebra,SparseArrays,Statistics,Test projects: ., test - uses: julia-actions/julia-buildpkg@v1 env: diff --git a/Project.toml b/Project.toml index 6d976786..de92fb39 100644 --- a/Project.toml +++ b/Project.toml @@ -27,8 +27,7 @@ MuladdMacro = "0.2.4" OrdinaryDiffEqCore = "1.16" RecipesBase = "1.3.4" Reexport = "1.2.2" -SciMLBase = "2.68 - 2.77" -SciMLOperators = "0.3 - 0.3.12" +SciMLBase = "2.68" SimpleUnPack = "1" SparseArrays = "1" StaticArrays = "1.9.7" From 3df55a04924a6681488c2b6f82ee11ade6843742 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Mon, 24 Mar 2025 11:07:48 +0100 Subject: [PATCH 04/26] remove debugging statements --- src/mprk.jl | 2 -- test/runtests.jl | 12 +++--------- 2 files changed, 3 insertions(+), 11 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index 5d002ffd..1faef363 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -429,9 +429,7 @@ end build_mprk_matrix!(P, P, σ, dt) # Same as linres = P \ uprev - display(linsolve.A) linsolve.A = P - display(linsolve.A) linres = solve!(linsolve) u .= linres diff --git a/test/runtests.jl b/test/runtests.jl index b0cf1e1a..d248ee83 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1276,16 +1276,10 @@ end # deliver the same results @testset "Different matrix types (conservative)" begin prod_1! = (P, u, p, t) -> begin - # fill!(P, zero(eltype(P))) - for j in axes(P, 2) - for idx in nzrange(P, j) - i = rowvals(P)[idx] - nonzeros(P)[idx] = i * u[i] - end + fill!(P, zero(eltype(P))) + for i in 1:(length(u) - 1) + P[i, i + 1] = i * u[i] end - # for i in 1:(length(u) - 1) - # P[i, i + 1] = i * u[i] - # end return nothing end From de15d15f7559ca3ad48fcdf46e9b0d6589107f41 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Mon, 24 Mar 2025 11:11:41 +0100 Subject: [PATCH 05/26] fix warnings related to alias --- src/mprk.jl | 24 ++++++++++++++++++------ src/sspmprk.jl | 16 ++++++++++++---- 2 files changed, 30 insertions(+), 10 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index 1faef363..08c5fbcc 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -357,7 +357,9 @@ function alg_cache(alg::MPE, u, rate_prototype, ::Type{uEltypeNoUnits}, # as well as to store the system matrix of the linear system # Right hand side of linear system is always uprev linprob = LinearProblem(P, _vec(uprev)) - linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, + linsolve = init(linprob, alg.linsolve, + alias = LinearSolve.LinearAliasSpecifier(; alias_A = true, + alias_b = true), assumptions = LinearSolve.OperatorAssumptions(true)) MPEConservativeCache(P, σ, tab, linsolve) @@ -366,7 +368,9 @@ function alg_cache(alg::MPE, u, rate_prototype, ::Type{uEltypeNoUnits}, # We use P to store the evaluation of the PDS # as well as to store the system matrix of the linear system linprob = LinearProblem(P, _vec(linsolve_rhs)) - linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, + linsolve = init(linprob, alg.linsolve, + alias = LinearSolve.LinearAliasSpecifier(; alias_A = true, + alias_b = true), assumptions = LinearSolve.OperatorAssumptions(true)) MPECache(P, similar(u), σ, tab, linsolve_rhs, linsolve) @@ -664,7 +668,9 @@ function alg_cache(alg::MPRK22, u, rate_prototype, ::Type{uEltypeNoUnits}, # not be altered, since it is needed to compute the adaptive time step # size. linprob = LinearProblem(P2, _vec(tmp)) - linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, + linsolve = init(linprob, alg.linsolve, + alias = LinearSolve.LinearAliasSpecifier(; alias_A = true, + alias_b = true), assumptions = LinearSolve.OperatorAssumptions(true)) MPRK22ConservativeCache(tmp, P, P2, σ, @@ -672,7 +678,9 @@ function alg_cache(alg::MPRK22, u, rate_prototype, ::Type{uEltypeNoUnits}, linsolve) elseif f isa PDSFunction linprob = LinearProblem(P2, _vec(tmp)) - linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, + linsolve = init(linprob, alg.linsolve, + alias = LinearSolve.LinearAliasSpecifier(; alias_A = true, + alias_b = true), assumptions = LinearSolve.OperatorAssumptions(true)) MPRK22Cache(tmp, P, P2, @@ -1262,7 +1270,9 @@ function alg_cache(alg::Union{MPRK43I, MPRK43II}, u, rate_prototype, ::Type{uElt # not be altered, since it is needed to compute the adaptive time step # size. linprob = LinearProblem(P3, _vec(tmp)) - linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, + linsolve = init(linprob, alg.linsolve, + alias = LinearSolve.LinearAliasSpecifier(; alias_A = true, + alias_b = true), assumptions = LinearSolve.OperatorAssumptions(true)) MPRK43ConservativeCache(tmp, tmp2, P, P2, P3, σ, tab, linsolve) elseif f isa PDSFunction @@ -1271,7 +1281,9 @@ function alg_cache(alg::Union{MPRK43I, MPRK43II}, u, rate_prototype, ::Type{uElt D3 = similar(u) linprob = LinearProblem(P3, _vec(tmp)) - linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, + linsolve = init(linprob, alg.linsolve, + alias = LinearSolve.LinearAliasSpecifier(; alias_A = true, + alias_b = true), assumptions = LinearSolve.OperatorAssumptions(true)) MPRK43Cache(tmp, tmp2, P, P2, P3, D, D2, D3, σ, tab, linsolve) diff --git a/src/sspmprk.jl b/src/sspmprk.jl index cc087225..9212acb0 100644 --- a/src/sspmprk.jl +++ b/src/sspmprk.jl @@ -240,7 +240,9 @@ function alg_cache(alg::SSPMPRK22, u, rate_prototype, ::Type{uEltypeNoUnits}, if f isa ConservativePDSFunction linprob = LinearProblem(P2, _vec(tmp)) - linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, + linsolve = init(linprob, alg.linsolve, + alias = LinearSolve.LinearAliasSpecifier(; alias_A = true, + alias_b = true), assumptions = LinearSolve.OperatorAssumptions(true)) SSPMPRK22ConservativeCache(tmp, P, P2, σ, @@ -248,7 +250,9 @@ function alg_cache(alg::SSPMPRK22, u, rate_prototype, ::Type{uEltypeNoUnits}, linsolve) elseif f isa PDSFunction linprob = LinearProblem(P2, _vec(tmp)) - linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, + linsolve = init(linprob, alg.linsolve, + alias = LinearSolve.LinearAliasSpecifier(; alias_A = true, + alias_b = true), assumptions = LinearSolve.OperatorAssumptions(true)) SSPMPRK22Cache(tmp, P, P2, @@ -775,7 +779,9 @@ function alg_cache(alg::SSPMPRK43, u, rate_prototype, ::Type{uEltypeNoUnits}, if f isa ConservativePDSFunction linprob = LinearProblem(P3, _vec(tmp)) - linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, + linsolve = init(linprob, alg.linsolve, + alias = LinearSolve.LinearAliasSpecifier(; alias_A = true, + alias_b = true), assumptions = LinearSolve.OperatorAssumptions(true)) SSPMPRK43ConservativeCache(tmp, tmp2, P, P2, P3, σ, ρ, tab, linsolve) elseif f isa PDSFunction @@ -784,7 +790,9 @@ function alg_cache(alg::SSPMPRK43, u, rate_prototype, ::Type{uEltypeNoUnits}, D3 = similar(u) linprob = LinearProblem(P3, _vec(tmp)) - linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, + linsolve = init(linprob, alg.linsolve, + alias = LinearSolve.LinearAliasSpecifier(; alias_A = true, + alias_b = true), assumptions = LinearSolve.OperatorAssumptions(true)) SSPMPRK43Cache(tmp, tmp2, P, P2, P3, D, D2, D3, σ, ρ, tab, linsolve) From d6504196e6e70048ae10fcd0d2cebe85fb6ca0c7 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Mon, 24 Mar 2025 11:21:29 +0100 Subject: [PATCH 06/26] bump compat of RecursiveFactorization to v0.2.14 --- test/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/Project.toml b/test/Project.toml index 978b479d..7c2b6a3d 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -31,7 +31,7 @@ OrdinaryDiffEqTsit5 = "1.1" OrdinaryDiffEqVerner = "1.1" Plots = "1.25.11" RecipesBase = "1.3.4" -RecursiveFactorization = "0.2" +RecursiveFactorization = "0.2.14" SparseArrays = "1" StaticArrays = "1.9.7" Statistics = "1" From d0f10d351727cbb2a17e9c6c5f605ad91184c6d0 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Mon, 24 Mar 2025 12:15:16 +0100 Subject: [PATCH 07/26] use KLUFactorization only in test --- src/PositiveIntegrators.jl | 2 +- src/mprk.jl | 8 ++++---- src/sspmprk.jl | 4 ++-- test/runtests.jl | 30 ++++++++++++++++++------------ 4 files changed, 25 insertions(+), 19 deletions(-) diff --git a/src/PositiveIntegrators.jl b/src/PositiveIntegrators.jl index ba0ec4cd..51306dec 100644 --- a/src/PositiveIntegrators.jl +++ b/src/PositiveIntegrators.jl @@ -24,7 +24,7 @@ using OrdinaryDiffEqCore: OrdinaryDiffEqCore, OrdinaryDiffEqAlgorithm, ODESoluti using SymbolicIndexingInterface: SymbolicIndexingInterface -using LinearSolve: LinearSolve, LinearProblem, KLUFactorization, solve! +using LinearSolve: LinearSolve, LinearProblem, LUFactorization, solve! import SciMLBase: interp_summary diff --git a/src/mprk.jl b/src/mprk.jl index 08c5fbcc..1328d86d 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -260,7 +260,7 @@ struct MPE{F, T} <: OrdinaryDiffEqAlgorithm small_constant_function::T end -function MPE(; linsolve = KLUFactorization(), small_constant = nothing) +function MPE(; linsolve = LUFactorization(), small_constant = nothing) if isnothing(small_constant) small_constant_function = floatmin elseif small_constant isa Number @@ -496,7 +496,7 @@ struct MPRK22{T, F, T2} <: OrdinaryDiffEqAdaptiveAlgorithm small_constant_function::T2 end -function MPRK22(alpha; linsolve = KLUFactorization(), +function MPRK22(alpha; linsolve = LUFactorization(), small_constant = nothing) if isnothing(small_constant) small_constant_function = floatmin @@ -910,7 +910,7 @@ struct MPRK43I{T, F, T2} <: OrdinaryDiffEqAdaptiveAlgorithm small_constant_function::T2 end -function MPRK43I(alpha, beta; linsolve = KLUFactorization(), +function MPRK43I(alpha, beta; linsolve = LUFactorization(), small_constant = nothing) if isnothing(small_constant) small_constant_function = floatmin @@ -1022,7 +1022,7 @@ function small_constant_function_MPRK43II(type) return small_constant end -function MPRK43II(gamma; linsolve = KLUFactorization(), +function MPRK43II(gamma; linsolve = LUFactorization(), small_constant = small_constant_function_MPRK43II) if small_constant isa Number small_constant_function = Returns(small_constant) diff --git a/src/sspmprk.jl b/src/sspmprk.jl index 9212acb0..f3a325a5 100644 --- a/src/sspmprk.jl +++ b/src/sspmprk.jl @@ -48,7 +48,7 @@ struct SSPMPRK22{T, F, T2} <: OrdinaryDiffEqAdaptiveAlgorithm small_constant_function::T2 end -function SSPMPRK22(alpha, beta; linsolve = KLUFactorization(), +function SSPMPRK22(alpha, beta; linsolve = LUFactorization(), small_constant = nothing) if isnothing(small_constant) small_constant_function = floatmin @@ -510,7 +510,7 @@ function small_constant_function_SSPMPRK43(type) return small_constant end -function SSPMPRK43(; linsolve = KLUFactorization(), +function SSPMPRK43(; linsolve = LUFactorization(), small_constant = small_constant_function_SSPMPRK43) if small_constant isa Number small_constant_function = Returns(small_constant) diff --git a/test/runtests.jl b/test/runtests.jl index d248ee83..fd457e9f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -17,7 +17,7 @@ using PositiveIntegrators # load RecursiveFactorization to get RFLUFactorization using RecursiveFactorization: RecursiveFactorization -using LinearSolve: RFLUFactorization, LUFactorization, KrylovJL_GMRES +using LinearSolve: RFLUFactorization, LUFactorization, KLUFactorization, KrylovJL_GMRES using Aqua: Aqua using RecipesBase: RecipesBase # only for Aqua tests @@ -1310,11 +1310,15 @@ end tspan = (0.0, 1.0) dt = 0.25 - @testset "$alg" for alg in (MPE(), - MPRK22(0.5), MPRK22(1.0), - MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), - MPRK43II(2.0 / 3.0), MPRK43II(0.5), - SSPMPRK22(0.5, 1.0), SSPMPRK43()) + @testset "$alg" for alg in (linsolve -> MPE(; linsolve), + linsolve -> MPRK22(0.5; linsolve), + linsolve -> MPRK22(1.0; linsolve), + linsolve -> MPRK43I(1.0, 0.5; linsolve), + linsolve -> MPRK43I(0.5, 0.75; linsolve), + linsolve -> MPRK43II(2.0 / 3.0; linsolve), + linsolve -> MPRK43II(0.5; linsolve), + linsolve -> SSPMPRK22(0.5, 1.0; linsolve), + linsolve -> SSPMPRK43(; linsolve)) for prod! in (prod_1!, prod_2!, prod_3!) prod = (u, p, t) -> begin P = similar(u, (length(u), length(u))) @@ -1334,14 +1338,16 @@ end prob_sparse_op = ConservativePDSProblem(prod, u0, tspan; p_prototype = P_sparse) - sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg; dt, + alg_lu = alg(LUFactorization()) + alg_klu = alg(KLUFactorization()) + sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg_lu; dt, adaptive = false) - sol_tridiagonal_op = solve(prob_tridiagonal_op, alg; dt, + sol_tridiagonal_op = solve(prob_tridiagonal_op, alg_lu; dt, adaptive = false) - sol_dense_ip = solve(prob_dense_ip, alg; dt, adaptive = false) - sol_dense_op = solve(prob_dense_op, alg; dt, adaptive = false) - sol_sparse_ip = solve(prob_sparse_ip, alg; dt, adaptive = false) - sol_sparse_op = solve(prob_sparse_op, alg; dt, adaptive = false) + sol_dense_ip = solve(prob_dense_ip, alg_lu; dt, adaptive = false) + sol_dense_op = solve(prob_dense_op, alg_lu; dt, adaptive = false) + sol_sparse_ip = solve(prob_sparse_ip, alg_klu; dt, adaptive = false) + sol_sparse_op = solve(prob_sparse_op, alg_lu; dt, adaptive = false) @test sol_tridiagonal_ip.t ≈ sol_tridiagonal_op.t @test sol_dense_ip.t ≈ sol_dense_op.t From 583306b72e50b18c45d3832c21d129ad35b53227 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Mon, 24 Mar 2025 12:25:25 +0100 Subject: [PATCH 08/26] bump compat of RecursiveFactorization.jl to v0.2.23 --- test/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/Project.toml b/test/Project.toml index 7c2b6a3d..7d33dc27 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -31,7 +31,7 @@ OrdinaryDiffEqTsit5 = "1.1" OrdinaryDiffEqVerner = "1.1" Plots = "1.25.11" RecipesBase = "1.3.4" -RecursiveFactorization = "0.2.14" +RecursiveFactorization = "0.2.23" SparseArrays = "1" StaticArrays = "1.9.7" Statistics = "1" From 50a0da4919760578dcbf97375d297d1cda02e943 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Mon, 24 Mar 2025 12:32:01 +0100 Subject: [PATCH 09/26] remove SciMLOperators again --- Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/Project.toml b/Project.toml index de92fb39..6de1efc3 100644 --- a/Project.toml +++ b/Project.toml @@ -12,7 +12,6 @@ OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" -SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961" SimpleUnPack = "ce78b400-467f-4804-87d8-8f486da07d0a" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" From 58f873d1a2db2af22175b040180d0aedb987e5f5 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Mon, 24 Mar 2025 12:56:34 +0100 Subject: [PATCH 10/26] fix change of sparsity pattern also in other tests --- test/runtests.jl | 102 ++++++++++++++++++++++++++++++----------------- 1 file changed, 65 insertions(+), 37 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index fd457e9f..3e2ddf86 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1402,11 +1402,15 @@ end rtol = sqrt(eps(Float32)) - @testset "$alg" for alg in (MPE(), - MPRK22(0.5), MPRK22(1.0), - MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), - MPRK43II(2.0 / 3.0), MPRK43II(0.5), - SSPMPRK22(0.5, 1.0), SSPMPRK43()) + @testset "$alg" for alg in (linsolve -> MPE(; linsolve), + linsolve -> MPRK22(0.5; linsolve), + linsolve -> MPRK22(1.0; linsolve), + linsolve -> MPRK43I(1.0, 0.5; linsolve), + linsolve -> MPRK43I(0.5, 0.75; linsolve), + linsolve -> MPRK43II(2.0 / 3.0; linsolve), + linsolve -> MPRK43II(0.5; linsolve), + linsolve -> SSPMPRK22(0.5, 1.0; linsolve), + linsolve -> SSPMPRK43(; linsolve)) for prod! in (prod_1!, prod_2!, prod_3!) prod = (u, p, t) -> begin P = similar(u, (length(u), length(u))) @@ -1426,12 +1430,14 @@ end prob_sparse_op = ConservativePDSProblem(prod, u0, tspan; p_prototype = P_sparse) - sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg; dt) - sol_tridiagonal_op = solve(prob_tridiagonal_op, alg; dt) - sol_dense_ip = solve(prob_dense_ip, alg; dt) - sol_dense_op = solve(prob_dense_op, alg; dt) - sol_sparse_ip = solve(prob_sparse_ip, alg; dt) - sol_sparse_op = solve(prob_sparse_op, alg; dt) + alg_lu = alg(LUFactorization()) + alg_klu = alg(KLUFactorization()) + sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg_lu; dt) + sol_tridiagonal_op = solve(prob_tridiagonal_op, alg_lu; dt) + sol_dense_ip = solve(prob_dense_ip, alg_lu; dt) + sol_dense_op = solve(prob_dense_op, alg_lu; dt) + sol_sparse_ip = solve(prob_sparse_ip, alg_klu; dt) + sol_sparse_op = solve(prob_sparse_op, alg_lu; dt) @test isapprox(sol_tridiagonal_ip.t, sol_tridiagonal_op.t; rtol) @test isapprox(sol_dense_ip.t, sol_dense_op.t; rtol) @@ -1517,11 +1523,15 @@ end tspan = (0.0, 1.0) dt = 0.25 - @testset "$alg" for alg in (MPE(), - MPRK22(0.5), MPRK22(1.0), - MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), - MPRK43II(2.0 / 3.0), MPRK43II(0.5), - SSPMPRK22(0.5, 1.0), SSPMPRK43()) + @testset "$alg" for alg in (linsolve -> MPE(; linsolve), + linsolve -> MPRK22(0.5; linsolve), + linsolve -> MPRK22(1.0; linsolve), + linsolve -> MPRK43I(1.0, 0.5; linsolve), + linsolve -> MPRK43I(0.5, 0.75; linsolve), + linsolve -> MPRK43II(2.0 / 3.0; linsolve), + linsolve -> MPRK43II(0.5; linsolve), + linsolve -> SSPMPRK22(0.5, 1.0; linsolve), + linsolve -> SSPMPRK43(; linsolve)) for (prod!, dest!) in zip((prod_1!, prod_2!, prod_3!), (dest_1!, dest_2!, dest_3!)) prod = (u, p, t) -> begin @@ -1547,17 +1557,19 @@ end prob_sparse_op = PDSProblem(prod, dest, u0, tspan; p_prototype = P_sparse) - sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg; + alg_lu = alg(LUFactorization()) + alg_klu = alg(KLUFactorization()) + sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg_lu; dt, adaptive = false) - sol_tridiagonal_op = solve(prob_tridiagonal_op, alg; + sol_tridiagonal_op = solve(prob_tridiagonal_op, alg_lu; dt, adaptive = false) - sol_dense_ip = solve(prob_dense_ip, alg; + sol_dense_ip = solve(prob_dense_ip, alg_lu; dt, adaptive = false) - sol_dense_op = solve(prob_dense_op, alg; + sol_dense_op = solve(prob_dense_op, alg_lu; dt, adaptive = false) - sol_sparse_ip = solve(prob_sparse_ip, alg; + sol_sparse_ip = solve(prob_sparse_ip, alg_klu; dt, adaptive = false) - sol_sparse_op = solve(prob_sparse_op, alg; + sol_sparse_op = solve(prob_sparse_op, alg_lu; dt, adaptive = false) @test sol_tridiagonal_ip.t ≈ sol_tridiagonal_op.t @@ -1643,11 +1655,15 @@ end dt = 0.25 rtol = sqrt(eps(Float32)) - @testset "$alg" for alg in (MPE(), - MPRK22(0.5), MPRK22(1.0), - MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), - MPRK43II(2.0 / 3.0), MPRK43II(0.5), - SSPMPRK22(0.5, 1.0), SSPMPRK43()) + @testset "$alg" for alg in (linsolve -> MPE(; linsolve), + linsolve -> MPRK22(0.5; linsolve), + linsolve -> MPRK22(1.0; linsolve), + linsolve -> MPRK43I(1.0, 0.5; linsolve), + linsolve -> MPRK43I(0.5, 0.75; linsolve), + linsolve -> MPRK43II(2.0 / 3.0; linsolve), + linsolve -> MPRK43II(0.5; linsolve), + linsolve -> SSPMPRK22(0.5, 1.0; linsolve), + linsolve -> SSPMPRK43(; linsolve)) for (prod!, dest!) in zip((prod_1!, prod_2!, prod_3!), (dest_1!, dest_2!, dest_3!)) prod! = prod_3! @@ -1675,17 +1691,19 @@ end prob_sparse_op = PDSProblem(prod, dest, u0, tspan; p_prototype = P_sparse) - sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg; + alg_lu = alg(LUFactorization()) + alg_klu = alg(KLUFactorization()) + sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg_lu; dt) - sol_tridiagonal_op = solve(prob_tridiagonal_op, alg; + sol_tridiagonal_op = solve(prob_tridiagonal_op, alg_lu; dt) - sol_dense_ip = solve(prob_dense_ip, alg; + sol_dense_ip = solve(prob_dense_ip, alg_lu; dt) - sol_dense_op = solve(prob_dense_op, alg; + sol_dense_op = solve(prob_dense_op, alg_lu; dt) - sol_sparse_ip = solve(prob_sparse_ip, alg; + sol_sparse_ip = solve(prob_sparse_ip, alg_klu; dt) - sol_sparse_op = solve(prob_sparse_op, alg; + sol_sparse_op = solve(prob_sparse_op, alg_lu; dt) @test isapprox(sol_tridiagonal_ip.t, sol_tridiagonal_op.t; rtol) @@ -1760,13 +1778,23 @@ end prob_sparse2 = PDSProblem(prod_sparse!, dest!, u0, tspan; p_prototype = P_sparse) #solve and test - for alg in (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK43I(1.0, 0.5), - MPRK43I(0.5, 0.75), - MPRK43II(2.0 / 3.0), MPRK43II(0.5), SSPMPRK22(0.5, 1.0), - SSPMPRK43()) + for alg in (linsolve -> MPE(; linsolve), + linsolve -> MPRK22(0.5; linsolve), + linsolve -> MPRK22(1.0; linsolve), + linsolve -> MPRK43I(1.0, 0.5; linsolve), + linsolve -> MPRK43I(0.5, 0.75; linsolve), + linsolve -> MPRK43II(2.0 / 3.0; linsolve), + linsolve -> MPRK43II(0.5; linsolve), + linsolve -> SSPMPRK22(0.5, 1.0; linsolve), + linsolve -> SSPMPRK43(; linsolve)) for prob in (prob_default, prob_tridiagonal, prob_dense, prob_sparse, prob_default2, prob_tridiagonal2, prob_dense2, prob_sparse2) + if prob == prob_sparse || prob == prob_sparse2 + alg = alg(KLUFactorization()) + else + alg = alg(LUFactorization()) + end sol1 = solve(prob, alg; dt, adaptive = false) # test get_tmp_cache and integrator interface - modifying From e23318934470cf17bc0e98ef83e05f4336fafc43 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Mon, 24 Mar 2025 13:33:00 +0100 Subject: [PATCH 11/26] fix --- test/runtests.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 3e2ddf86..86fcfb94 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1791,11 +1791,11 @@ end prob_default2, prob_tridiagonal2, prob_dense2, prob_sparse2) if prob == prob_sparse || prob == prob_sparse2 - alg = alg(KLUFactorization()) + alg_ = alg(KLUFactorization()) else - alg = alg(LUFactorization()) + alg_ = alg(LUFactorization()) end - sol1 = solve(prob, alg; dt, adaptive = false) + sol1 = solve(prob, alg_; dt, adaptive = false) # test get_tmp_cache and integrator interface - modifying # values from the cache should not change the final results @@ -2062,7 +2062,7 @@ end # Check that approximations, and thus the Patankar weights, # remain positive to avoid division by zero. - @testset "Positvity check" begin + @testset "Positivity check" begin # For this problem u[1] decreases montonically to 0 very fast. # We perform 10^5 steps and check that u[end] does not contain any NaNs u0 = [0.9, 0.1] From de6100ed2a22cb20e54a22638682a89550a0a997 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Mon, 24 Mar 2025 13:49:08 +0100 Subject: [PATCH 12/26] finally fixing it --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 86fcfb94..bc134921 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1799,7 +1799,7 @@ end # test get_tmp_cache and integrator interface - modifying # values from the cache should not change the final results - integrator = init(prob, alg; dt, adaptive = false) + integrator = init(prob, alg_; dt, adaptive = false) PositiveIntegrators.OrdinaryDiffEqCore.step!(integrator) cache = @inferred PositiveIntegrators.get_tmp_cache(integrator) @test !isempty(cache) From 4226dd62b9d4037a2583fcaec6247b0f76807d6c Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Mon, 24 Mar 2025 13:58:27 +0100 Subject: [PATCH 13/26] use same pattern as elsewhere --- test/runtests.jl | 97 +++++++++++++++++++++++++----------------------- 1 file changed, 51 insertions(+), 46 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index bc134921..d65f14bb 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1310,15 +1310,16 @@ end tspan = (0.0, 1.0) dt = 0.25 - @testset "$alg" for alg in (linsolve -> MPE(; linsolve), - linsolve -> MPRK22(0.5; linsolve), - linsolve -> MPRK22(1.0; linsolve), - linsolve -> MPRK43I(1.0, 0.5; linsolve), - linsolve -> MPRK43I(0.5, 0.75; linsolve), - linsolve -> MPRK43II(2.0 / 3.0; linsolve), - linsolve -> MPRK43II(0.5; linsolve), - linsolve -> SSPMPRK22(0.5, 1.0; linsolve), - linsolve -> SSPMPRK43(; linsolve)) + algs = (MPE, (; kwargs...) -> MPRK22(1.0; kwargs...), + (; kwargs...) -> MPRK22(0.5; kwargs...), + (; kwargs...) -> MPRK22(1.0; kwargs...), + (; kwargs...) -> MPRK43I(1.0, 0.5; kwargs...), + (; kwargs...) -> MPRK43I(0.5, 0.75; kwargs...), + (; kwargs...) -> MPRK43II(0.5; kwargs...), + (; kwargs...) -> MPRK43II(2.0 / 3.0; kwargs...), + (; kwargs...) -> SSPMPRK22(0.5, 1.0; kwargs...), + (; kwargs...) -> SSPMPRK43(; kwargs...)) + @testset "$alg" for alg in algs for prod! in (prod_1!, prod_2!, prod_3!) prod = (u, p, t) -> begin P = similar(u, (length(u), length(u))) @@ -1402,15 +1403,16 @@ end rtol = sqrt(eps(Float32)) - @testset "$alg" for alg in (linsolve -> MPE(; linsolve), - linsolve -> MPRK22(0.5; linsolve), - linsolve -> MPRK22(1.0; linsolve), - linsolve -> MPRK43I(1.0, 0.5; linsolve), - linsolve -> MPRK43I(0.5, 0.75; linsolve), - linsolve -> MPRK43II(2.0 / 3.0; linsolve), - linsolve -> MPRK43II(0.5; linsolve), - linsolve -> SSPMPRK22(0.5, 1.0; linsolve), - linsolve -> SSPMPRK43(; linsolve)) + algs = (MPE, (; kwargs...) -> MPRK22(1.0; kwargs...), + (; kwargs...) -> MPRK22(0.5; kwargs...), + (; kwargs...) -> MPRK22(1.0; kwargs...), + (; kwargs...) -> MPRK43I(1.0, 0.5; kwargs...), + (; kwargs...) -> MPRK43I(0.5, 0.75; kwargs...), + (; kwargs...) -> MPRK43II(0.5; kwargs...), + (; kwargs...) -> MPRK43II(2.0 / 3.0; kwargs...), + (; kwargs...) -> SSPMPRK22(0.5, 1.0; kwargs...), + (; kwargs...) -> SSPMPRK43(; kwargs...)) + @testset "$alg" for alg in algs for prod! in (prod_1!, prod_2!, prod_3!) prod = (u, p, t) -> begin P = similar(u, (length(u), length(u))) @@ -1523,15 +1525,16 @@ end tspan = (0.0, 1.0) dt = 0.25 - @testset "$alg" for alg in (linsolve -> MPE(; linsolve), - linsolve -> MPRK22(0.5; linsolve), - linsolve -> MPRK22(1.0; linsolve), - linsolve -> MPRK43I(1.0, 0.5; linsolve), - linsolve -> MPRK43I(0.5, 0.75; linsolve), - linsolve -> MPRK43II(2.0 / 3.0; linsolve), - linsolve -> MPRK43II(0.5; linsolve), - linsolve -> SSPMPRK22(0.5, 1.0; linsolve), - linsolve -> SSPMPRK43(; linsolve)) + algs = (MPE, (; kwargs...) -> MPRK22(1.0; kwargs...), + (; kwargs...) -> MPRK22(0.5; kwargs...), + (; kwargs...) -> MPRK22(1.0; kwargs...), + (; kwargs...) -> MPRK43I(1.0, 0.5; kwargs...), + (; kwargs...) -> MPRK43I(0.5, 0.75; kwargs...), + (; kwargs...) -> MPRK43II(0.5; kwargs...), + (; kwargs...) -> MPRK43II(2.0 / 3.0; kwargs...), + (; kwargs...) -> SSPMPRK22(0.5, 1.0; kwargs...), + (; kwargs...) -> SSPMPRK43(; kwargs...)) + @testset "$alg" for alg in algs for (prod!, dest!) in zip((prod_1!, prod_2!, prod_3!), (dest_1!, dest_2!, dest_3!)) prod = (u, p, t) -> begin @@ -1655,15 +1658,16 @@ end dt = 0.25 rtol = sqrt(eps(Float32)) - @testset "$alg" for alg in (linsolve -> MPE(; linsolve), - linsolve -> MPRK22(0.5; linsolve), - linsolve -> MPRK22(1.0; linsolve), - linsolve -> MPRK43I(1.0, 0.5; linsolve), - linsolve -> MPRK43I(0.5, 0.75; linsolve), - linsolve -> MPRK43II(2.0 / 3.0; linsolve), - linsolve -> MPRK43II(0.5; linsolve), - linsolve -> SSPMPRK22(0.5, 1.0; linsolve), - linsolve -> SSPMPRK43(; linsolve)) + algs = (MPE, (; kwargs...) -> MPRK22(1.0; kwargs...), + (; kwargs...) -> MPRK22(0.5; kwargs...), + (; kwargs...) -> MPRK22(1.0; kwargs...), + (; kwargs...) -> MPRK43I(1.0, 0.5; kwargs...), + (; kwargs...) -> MPRK43I(0.5, 0.75; kwargs...), + (; kwargs...) -> MPRK43II(0.5; kwargs...), + (; kwargs...) -> MPRK43II(2.0 / 3.0; kwargs...), + (; kwargs...) -> SSPMPRK22(0.5, 1.0; kwargs...), + (; kwargs...) -> SSPMPRK43(; kwargs...)) + @testset "$alg" for alg in algs for (prod!, dest!) in zip((prod_1!, prod_2!, prod_3!), (dest_1!, dest_2!, dest_3!)) prod! = prod_3! @@ -1777,16 +1781,17 @@ end p_prototype = P_dense) prob_sparse2 = PDSProblem(prod_sparse!, dest!, u0, tspan; p_prototype = P_sparse) - #solve and test - for alg in (linsolve -> MPE(; linsolve), - linsolve -> MPRK22(0.5; linsolve), - linsolve -> MPRK22(1.0; linsolve), - linsolve -> MPRK43I(1.0, 0.5; linsolve), - linsolve -> MPRK43I(0.5, 0.75; linsolve), - linsolve -> MPRK43II(2.0 / 3.0; linsolve), - linsolve -> MPRK43II(0.5; linsolve), - linsolve -> SSPMPRK22(0.5, 1.0; linsolve), - linsolve -> SSPMPRK43(; linsolve)) + # solve and test + algs = (MPE, (; kwargs...) -> MPRK22(1.0; kwargs...), + (; kwargs...) -> MPRK22(0.5; kwargs...), + (; kwargs...) -> MPRK22(1.0; kwargs...), + (; kwargs...) -> MPRK43I(1.0, 0.5; kwargs...), + (; kwargs...) -> MPRK43I(0.5, 0.75; kwargs...), + (; kwargs...) -> MPRK43II(0.5; kwargs...), + (; kwargs...) -> MPRK43II(2.0 / 3.0; kwargs...), + (; kwargs...) -> SSPMPRK22(0.5, 1.0; kwargs...), + (; kwargs...) -> SSPMPRK43(; kwargs...)) + for alg in algs for prob in (prob_default, prob_tridiagonal, prob_dense, prob_sparse, prob_default2, prob_tridiagonal2, prob_dense2, prob_sparse2) From 642565f39ed1bfc0b500d7067e0e22b6d8898fb9 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Mon, 24 Mar 2025 14:13:11 +0100 Subject: [PATCH 14/26] fix keywords --- test/runtests.jl | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index d65f14bb..0416bfa3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1339,8 +1339,8 @@ end prob_sparse_op = ConservativePDSProblem(prod, u0, tspan; p_prototype = P_sparse) - alg_lu = alg(LUFactorization()) - alg_klu = alg(KLUFactorization()) + alg_lu = alg(; linsolve = LUFactorization()) + alg_klu = alg(; linsolve = KLUFactorization()) sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg_lu; dt, adaptive = false) sol_tridiagonal_op = solve(prob_tridiagonal_op, alg_lu; dt, @@ -1432,8 +1432,8 @@ end prob_sparse_op = ConservativePDSProblem(prod, u0, tspan; p_prototype = P_sparse) - alg_lu = alg(LUFactorization()) - alg_klu = alg(KLUFactorization()) + alg_lu = alg(; linsolve = LUFactorization()) + alg_klu = alg(; linsolve = KLUFactorization()) sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg_lu; dt) sol_tridiagonal_op = solve(prob_tridiagonal_op, alg_lu; dt) sol_dense_ip = solve(prob_dense_ip, alg_lu; dt) @@ -1560,8 +1560,8 @@ end prob_sparse_op = PDSProblem(prod, dest, u0, tspan; p_prototype = P_sparse) - alg_lu = alg(LUFactorization()) - alg_klu = alg(KLUFactorization()) + alg_lu = alg(; linsolve = LUFactorization()) + alg_klu = alg(; linsolve = KLUFactorization()) sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg_lu; dt, adaptive = false) sol_tridiagonal_op = solve(prob_tridiagonal_op, alg_lu; @@ -1695,8 +1695,8 @@ end prob_sparse_op = PDSProblem(prod, dest, u0, tspan; p_prototype = P_sparse) - alg_lu = alg(LUFactorization()) - alg_klu = alg(KLUFactorization()) + alg_lu = alg(; linsolve = LUFactorization()) + alg_klu = alg(; linsolve = KLUFactorization()) sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg_lu; dt) sol_tridiagonal_op = solve(prob_tridiagonal_op, alg_lu; @@ -1796,9 +1796,9 @@ end prob_default2, prob_tridiagonal2, prob_dense2, prob_sparse2) if prob == prob_sparse || prob == prob_sparse2 - alg_ = alg(KLUFactorization()) + alg_ = alg(; linsolve = KLUFactorization()) else - alg_ = alg(LUFactorization()) + alg_ = alg(; linsolve = LUFactorization()) end sol1 = solve(prob, alg_; dt, adaptive = false) From 949aabac9025803735a831d70cb392ea2f870767 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Tue, 25 Mar 2025 11:23:38 +0100 Subject: [PATCH 15/26] bump compat bounds to fix downgrade --- Project.toml | 4 ++-- docs/Project.toml | 2 +- test/Project.toml | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index 6de1efc3..3c8ac3c4 100644 --- a/Project.toml +++ b/Project.toml @@ -21,12 +21,12 @@ SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" [compat] FastBroadcast = "0.3.5" LinearAlgebra = "1" -LinearSolve = "2.32, 3" +LinearSolve = "3.7" MuladdMacro = "0.2.4" OrdinaryDiffEqCore = "1.16" RecipesBase = "1.3.4" Reexport = "1.2.2" -SciMLBase = "2.68" +SciMLBase = "2.78" SimpleUnPack = "1" SparseArrays = "1" StaticArrays = "1.9.7" diff --git a/docs/Project.toml b/docs/Project.toml index e17a197e..a806f8e4 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -27,7 +27,7 @@ DiffEqDevTools = "2.45.1" Documenter = "1" InteractiveUtils = "1" LinearAlgebra = "1" -LinearSolve = "2.32, 3" +LinearSolve = "3.7" OrdinaryDiffEqFIRK = "1.7" OrdinaryDiffEqLowOrderRK = "1.2" OrdinaryDiffEqRosenbrock = "1.4" diff --git a/test/Project.toml b/test/Project.toml index 7d33dc27..0f8b8f3a 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -23,7 +23,7 @@ ADTypes = "1.12" Aqua = "0.8" ExplicitImports = "1.0.1" LinearAlgebra = "1" -LinearSolve = "2.32, 3" +LinearSolve = "3.7" OrdinaryDiffEqLowOrderRK = "1.2" OrdinaryDiffEqRosenbrock = "1.4" OrdinaryDiffEqSDIRK = "1.2" From b38354f53edf4b51d2e75b9649f1a35392e665c4 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Tue, 25 Mar 2025 11:45:34 +0100 Subject: [PATCH 16/26] use algs --- test/runtests.jl | 24 ++++++++++-------------- 1 file changed, 10 insertions(+), 14 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 9277d520..69b1640a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2334,9 +2334,9 @@ end labelsB = ["B_1"; "B_2"; "B_3"; "B_4"; "B_5"] dts = (last(prob.tspan) - first(prob.tspan)) / 10.0 * 0.5 .^ (4:13) alg_ref = Vern7() - wp = work_precision_fixed(prob, [alg; alg; alg; alg; alg], labelsA, dts, + wp = work_precision_fixed(prob, algs, labelsA, dts, alg_ref) - work_precision_fixed!(wp, prob, [alg; alg; alg; alg; alg], labelsB, dts, + work_precision_fixed!(wp, prob, algs, labelsB, dts, alg_ref) # check that errors agree @@ -2369,12 +2369,10 @@ end abstols = 1 ./ 10 .^ (4:8) reltols = 1 ./ 10 .^ (3:7) alg_ref = Vern7() - wp = work_precision_adaptive(prob, [alg; alg; alg; alg; alg], labelsA, - abstols, - reltols, alg_ref) - work_precision_adaptive!(wp, prob, [alg; alg; alg; alg; alg], labelsB, - abstols, - reltols, alg_ref) + wp = work_precision_adaptive(prob, algs, labelsA, + abstols, reltols, alg_ref) + work_precision_adaptive!(wp, prob, algs, labelsB, + abstols, reltols, alg_ref) # check that errors agree for (i, _) in enumerate(abstols) @@ -2402,12 +2400,10 @@ end abstols = 1 ./ 10 .^ (4:8) reltols = 1 ./ 10 .^ (3:7) alg_ref = Rodas4P() - wp = work_precision_adaptive(prob, [alg; alg; alg; alg; alg], labelsA, - abstols, - reltols, alg_ref; adaptive_ref = true) - work_precision_adaptive!(wp, prob, [alg; alg; alg; alg; alg], labelsB, - abstols, - reltols, alg_ref; adaptive_ref = true) + wp = work_precision_adaptive(prob, algs, labelsA, + abstols, reltols, alg_ref; adaptive_ref = true) + work_precision_adaptive!(wp, prob, algs, labelsB, + abstols, reltols, alg_ref; adaptive_ref = true) # check that errors agree for (i, _) in enumerate(abstols) From f12c93bd26204946083947c29adb4c63e311f97a Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Tue, 25 Mar 2025 11:47:50 +0100 Subject: [PATCH 17/26] bump compat of OrdinaryDiffEqRosenbrock.jl to v1.5 --- debug.jl | 33 +++++++++++++++++++++++++++++++++ test/Project.toml | 2 +- 2 files changed, 34 insertions(+), 1 deletion(-) create mode 100644 debug.jl diff --git a/debug.jl b/debug.jl new file mode 100644 index 00000000..e01bd24d --- /dev/null +++ b/debug.jl @@ -0,0 +1,33 @@ +using LinearAlgebra +using SparseArrays +using PositiveIntegrators +using LinearSolve + +prod! = (P, u, p, t) -> begin + fill!(P, zero(eltype(P))) + # for j in axes(P, 2) + # for idx in nzrange(P, j) + # i = rowvals(P)[idx] + # nonzeros(P)[idx] = i * u[i] + # end + # end + for i in 1:(length(u) - 1) + P[i, i + 1] = i * u[i] + end + return nothing +end + +n = 4 +P_tridiagonal = Tridiagonal([0.1, 0.2, 0.3], + zeros(n), + [0.4, 0.5, 0.6]) +P_sparse = sparse(P_tridiagonal) +u0 = [1.0, 1.5, 2.0, 2.5] +tspan = (0.0, 1.0) +dt = 0.25 + +alg = MPE(; linsolve = KLUFactorization()) +display(P_sparse) +prob_sparse_ip = ConservativePDSProblem(prod!, u0, tspan; + p_prototype = P_sparse) +sol_sparse_ip = solve(prob_sparse_ip, alg; dt, adaptive = false) diff --git a/test/Project.toml b/test/Project.toml index 0f8b8f3a..20b39432 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -25,7 +25,7 @@ ExplicitImports = "1.0.1" LinearAlgebra = "1" LinearSolve = "3.7" OrdinaryDiffEqLowOrderRK = "1.2" -OrdinaryDiffEqRosenbrock = "1.4" +OrdinaryDiffEqRosenbrock = "1.5" OrdinaryDiffEqSDIRK = "1.2" OrdinaryDiffEqTsit5 = "1.1" OrdinaryDiffEqVerner = "1.1" From ac45d5ba16580c48a48af1e3fa156e020b4e4deb Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Tue, 25 Mar 2025 12:26:03 +0100 Subject: [PATCH 18/26] remove accidently pushed file again --- debug.jl | 33 --------------------------------- 1 file changed, 33 deletions(-) delete mode 100644 debug.jl diff --git a/debug.jl b/debug.jl deleted file mode 100644 index e01bd24d..00000000 --- a/debug.jl +++ /dev/null @@ -1,33 +0,0 @@ -using LinearAlgebra -using SparseArrays -using PositiveIntegrators -using LinearSolve - -prod! = (P, u, p, t) -> begin - fill!(P, zero(eltype(P))) - # for j in axes(P, 2) - # for idx in nzrange(P, j) - # i = rowvals(P)[idx] - # nonzeros(P)[idx] = i * u[i] - # end - # end - for i in 1:(length(u) - 1) - P[i, i + 1] = i * u[i] - end - return nothing -end - -n = 4 -P_tridiagonal = Tridiagonal([0.1, 0.2, 0.3], - zeros(n), - [0.4, 0.5, 0.6]) -P_sparse = sparse(P_tridiagonal) -u0 = [1.0, 1.5, 2.0, 2.5] -tspan = (0.0, 1.0) -dt = 0.25 - -alg = MPE(; linsolve = KLUFactorization()) -display(P_sparse) -prob_sparse_ip = ConservativePDSProblem(prod!, u0, tspan; - p_prototype = P_sparse) -sol_sparse_ip = solve(prob_sparse_ip, alg; dt, adaptive = false) From 5c50df61b08a68b3879f92dbf61d21a926271ad7 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Thu, 27 Mar 2025 17:49:14 +0100 Subject: [PATCH 19/26] try Downgrade with backported fix of SimpleNonlinearSolve.jl --- Project.toml | 4 ++-- docs/Project.toml | 2 +- test/Project.toml | 4 ++-- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/Project.toml b/Project.toml index 3c8ac3c4..6de1efc3 100644 --- a/Project.toml +++ b/Project.toml @@ -21,12 +21,12 @@ SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" [compat] FastBroadcast = "0.3.5" LinearAlgebra = "1" -LinearSolve = "3.7" +LinearSolve = "2.32, 3" MuladdMacro = "0.2.4" OrdinaryDiffEqCore = "1.16" RecipesBase = "1.3.4" Reexport = "1.2.2" -SciMLBase = "2.78" +SciMLBase = "2.68" SimpleUnPack = "1" SparseArrays = "1" StaticArrays = "1.9.7" diff --git a/docs/Project.toml b/docs/Project.toml index a806f8e4..e17a197e 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -27,7 +27,7 @@ DiffEqDevTools = "2.45.1" Documenter = "1" InteractiveUtils = "1" LinearAlgebra = "1" -LinearSolve = "3.7" +LinearSolve = "2.32, 3" OrdinaryDiffEqFIRK = "1.7" OrdinaryDiffEqLowOrderRK = "1.2" OrdinaryDiffEqRosenbrock = "1.4" diff --git a/test/Project.toml b/test/Project.toml index 20b39432..7d33dc27 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -23,9 +23,9 @@ ADTypes = "1.12" Aqua = "0.8" ExplicitImports = "1.0.1" LinearAlgebra = "1" -LinearSolve = "3.7" +LinearSolve = "2.32, 3" OrdinaryDiffEqLowOrderRK = "1.2" -OrdinaryDiffEqRosenbrock = "1.5" +OrdinaryDiffEqRosenbrock = "1.4" OrdinaryDiffEqSDIRK = "1.2" OrdinaryDiffEqTsit5 = "1.1" OrdinaryDiffEqVerner = "1.1" From 77159bb3b89323fe439d93e561966b7b8f5885f8 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Thu, 27 Mar 2025 20:17:05 +0100 Subject: [PATCH 20/26] bump compats again Downgrade still broken --- Project.toml | 4 ++-- docs/Project.toml | 2 +- test/Project.toml | 4 ++-- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/Project.toml b/Project.toml index 6de1efc3..3c8ac3c4 100644 --- a/Project.toml +++ b/Project.toml @@ -21,12 +21,12 @@ SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" [compat] FastBroadcast = "0.3.5" LinearAlgebra = "1" -LinearSolve = "2.32, 3" +LinearSolve = "3.7" MuladdMacro = "0.2.4" OrdinaryDiffEqCore = "1.16" RecipesBase = "1.3.4" Reexport = "1.2.2" -SciMLBase = "2.68" +SciMLBase = "2.78" SimpleUnPack = "1" SparseArrays = "1" StaticArrays = "1.9.7" diff --git a/docs/Project.toml b/docs/Project.toml index e17a197e..a806f8e4 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -27,7 +27,7 @@ DiffEqDevTools = "2.45.1" Documenter = "1" InteractiveUtils = "1" LinearAlgebra = "1" -LinearSolve = "2.32, 3" +LinearSolve = "3.7" OrdinaryDiffEqFIRK = "1.7" OrdinaryDiffEqLowOrderRK = "1.2" OrdinaryDiffEqRosenbrock = "1.4" diff --git a/test/Project.toml b/test/Project.toml index 7d33dc27..20b39432 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -23,9 +23,9 @@ ADTypes = "1.12" Aqua = "0.8" ExplicitImports = "1.0.1" LinearAlgebra = "1" -LinearSolve = "2.32, 3" +LinearSolve = "3.7" OrdinaryDiffEqLowOrderRK = "1.2" -OrdinaryDiffEqRosenbrock = "1.4" +OrdinaryDiffEqRosenbrock = "1.5" OrdinaryDiffEqSDIRK = "1.2" OrdinaryDiffEqTsit5 = "1.1" OrdinaryDiffEqVerner = "1.1" From 048f96c48156934d4008f49f23823012a135e377 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Fri, 28 Mar 2025 11:13:02 +0100 Subject: [PATCH 21/26] allow LinearSolve v2.39.1 --- Project.toml | 2 +- docs/Project.toml | 2 +- test/Project.toml | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 3c8ac3c4..f1b9749c 100644 --- a/Project.toml +++ b/Project.toml @@ -21,7 +21,7 @@ SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" [compat] FastBroadcast = "0.3.5" LinearAlgebra = "1" -LinearSolve = "3.7" +LinearSolve = "2.39.1, 3.7" MuladdMacro = "0.2.4" OrdinaryDiffEqCore = "1.16" RecipesBase = "1.3.4" diff --git a/docs/Project.toml b/docs/Project.toml index a806f8e4..c860a614 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -27,7 +27,7 @@ DiffEqDevTools = "2.45.1" Documenter = "1" InteractiveUtils = "1" LinearAlgebra = "1" -LinearSolve = "3.7" +LinearSolve = "2.39.1, 3.7" OrdinaryDiffEqFIRK = "1.7" OrdinaryDiffEqLowOrderRK = "1.2" OrdinaryDiffEqRosenbrock = "1.4" diff --git a/test/Project.toml b/test/Project.toml index 20b39432..aafaedd6 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -23,7 +23,7 @@ ADTypes = "1.12" Aqua = "0.8" ExplicitImports = "1.0.1" LinearAlgebra = "1" -LinearSolve = "3.7" +LinearSolve = "2.39.1, 3.7" OrdinaryDiffEqLowOrderRK = "1.2" OrdinaryDiffEqRosenbrock = "1.5" OrdinaryDiffEqSDIRK = "1.2" From 0ca900c799673dacf15778dab4ed3220debb8077 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Fri, 28 Mar 2025 11:31:06 +0100 Subject: [PATCH 22/26] allow SciMLBase.jl v2.68 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index f1b9749c..02a10e35 100644 --- a/Project.toml +++ b/Project.toml @@ -26,7 +26,7 @@ MuladdMacro = "0.2.4" OrdinaryDiffEqCore = "1.16" RecipesBase = "1.3.4" Reexport = "1.2.2" -SciMLBase = "2.78" +SciMLBase = "2.68" SimpleUnPack = "1" SparseArrays = "1" StaticArrays = "1.9.7" From 0c7234e0d98bb86e31f04764a9de9acdaf33f873 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Fri, 28 Mar 2025 11:36:06 +0100 Subject: [PATCH 23/26] require SciMLBase.jl v2.78 again --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 02a10e35..f1b9749c 100644 --- a/Project.toml +++ b/Project.toml @@ -26,7 +26,7 @@ MuladdMacro = "0.2.4" OrdinaryDiffEqCore = "1.16" RecipesBase = "1.3.4" Reexport = "1.2.2" -SciMLBase = "2.68" +SciMLBase = "2.78" SimpleUnPack = "1" SparseArrays = "1" StaticArrays = "1.9.7" From b8bcba6cb127ef5a2579cf8ce880604f62fa5d12 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Fri, 28 Mar 2025 11:59:54 +0100 Subject: [PATCH 24/26] require LinearSolve.jl v3.7.1 --- Project.toml | 2 +- docs/Project.toml | 2 +- test/Project.toml | 2 +- test/runtests.jl | 143 ++++++++++++++++------------------------------ 4 files changed, 52 insertions(+), 97 deletions(-) diff --git a/Project.toml b/Project.toml index f1b9749c..5e3fd71e 100644 --- a/Project.toml +++ b/Project.toml @@ -21,7 +21,7 @@ SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" [compat] FastBroadcast = "0.3.5" LinearAlgebra = "1" -LinearSolve = "2.39.1, 3.7" +LinearSolve = "2.39.1, 3.7.1" MuladdMacro = "0.2.4" OrdinaryDiffEqCore = "1.16" RecipesBase = "1.3.4" diff --git a/docs/Project.toml b/docs/Project.toml index c860a614..c413e644 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -27,7 +27,7 @@ DiffEqDevTools = "2.45.1" Documenter = "1" InteractiveUtils = "1" LinearAlgebra = "1" -LinearSolve = "2.39.1, 3.7" +LinearSolve = "2.39.1, 3.7.1" OrdinaryDiffEqFIRK = "1.7" OrdinaryDiffEqLowOrderRK = "1.2" OrdinaryDiffEqRosenbrock = "1.4" diff --git a/test/Project.toml b/test/Project.toml index aafaedd6..832d4d94 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -23,7 +23,7 @@ ADTypes = "1.12" Aqua = "0.8" ExplicitImports = "1.0.1" LinearAlgebra = "1" -LinearSolve = "2.39.1, 3.7" +LinearSolve = "2.39.1, 3.7.1" OrdinaryDiffEqLowOrderRK = "1.2" OrdinaryDiffEqRosenbrock = "1.5" OrdinaryDiffEqSDIRK = "1.2" diff --git a/test/runtests.jl b/test/runtests.jl index 69b1640a..0cd16ab2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1310,16 +1310,11 @@ end tspan = (0.0, 1.0) dt = 0.25 - algs = (MPE, (; kwargs...) -> MPRK22(1.0; kwargs...), - (; kwargs...) -> MPRK22(0.5; kwargs...), - (; kwargs...) -> MPRK22(1.0; kwargs...), - (; kwargs...) -> MPRK43I(1.0, 0.5; kwargs...), - (; kwargs...) -> MPRK43I(0.5, 0.75; kwargs...), - (; kwargs...) -> MPRK43II(0.5; kwargs...), - (; kwargs...) -> MPRK43II(2.0 / 3.0; kwargs...), - (; kwargs...) -> SSPMPRK22(0.5, 1.0; kwargs...), - (; kwargs...) -> SSPMPRK43(; kwargs...)) - @testset "$alg" for alg in algs + @testset "$alg" for alg in (MPE(), + MPRK22(0.5), MPRK22(1.0), + MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), + MPRK43II(2.0 / 3.0), MPRK43II(0.5), + SSPMPRK22(0.5, 1.0), SSPMPRK43()) for prod! in (prod_1!, prod_2!, prod_3!) prod = (u, p, t) -> begin P = similar(u, (length(u), length(u))) @@ -1339,16 +1334,14 @@ end prob_sparse_op = ConservativePDSProblem(prod, u0, tspan; p_prototype = P_sparse) - alg_lu = alg(; linsolve = LUFactorization()) - alg_klu = alg(; linsolve = KLUFactorization()) - sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg_lu; dt, + sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg; dt, adaptive = false) - sol_tridiagonal_op = solve(prob_tridiagonal_op, alg_lu; dt, + sol_tridiagonal_op = solve(prob_tridiagonal_op, alg; dt, adaptive = false) - sol_dense_ip = solve(prob_dense_ip, alg_lu; dt, adaptive = false) - sol_dense_op = solve(prob_dense_op, alg_lu; dt, adaptive = false) - sol_sparse_ip = solve(prob_sparse_ip, alg_klu; dt, adaptive = false) - sol_sparse_op = solve(prob_sparse_op, alg_lu; dt, adaptive = false) + sol_dense_ip = solve(prob_dense_ip, alg; dt, adaptive = false) + sol_dense_op = solve(prob_dense_op, alg; dt, adaptive = false) + sol_sparse_ip = solve(prob_sparse_ip, alg; dt, adaptive = false) + sol_sparse_op = solve(prob_sparse_op, alg; dt, adaptive = false) @test sol_tridiagonal_ip.t ≈ sol_tridiagonal_op.t @test sol_dense_ip.t ≈ sol_dense_op.t @@ -1403,16 +1396,11 @@ end rtol = sqrt(eps(Float32)) - algs = (MPE, (; kwargs...) -> MPRK22(1.0; kwargs...), - (; kwargs...) -> MPRK22(0.5; kwargs...), - (; kwargs...) -> MPRK22(1.0; kwargs...), - (; kwargs...) -> MPRK43I(1.0, 0.5; kwargs...), - (; kwargs...) -> MPRK43I(0.5, 0.75; kwargs...), - (; kwargs...) -> MPRK43II(0.5; kwargs...), - (; kwargs...) -> MPRK43II(2.0 / 3.0; kwargs...), - (; kwargs...) -> SSPMPRK22(0.5, 1.0; kwargs...), - (; kwargs...) -> SSPMPRK43(; kwargs...)) - @testset "$alg" for alg in algs + @testset "$alg" for alg in (MPE(), + MPRK22(0.5), MPRK22(1.0), + MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), + MPRK43II(2.0 / 3.0), MPRK43II(0.5), + SSPMPRK22(0.5, 1.0), SSPMPRK43()) for prod! in (prod_1!, prod_2!, prod_3!) prod = (u, p, t) -> begin P = similar(u, (length(u), length(u))) @@ -1432,14 +1420,12 @@ end prob_sparse_op = ConservativePDSProblem(prod, u0, tspan; p_prototype = P_sparse) - alg_lu = alg(; linsolve = LUFactorization()) - alg_klu = alg(; linsolve = KLUFactorization()) - sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg_lu; dt) - sol_tridiagonal_op = solve(prob_tridiagonal_op, alg_lu; dt) - sol_dense_ip = solve(prob_dense_ip, alg_lu; dt) - sol_dense_op = solve(prob_dense_op, alg_lu; dt) - sol_sparse_ip = solve(prob_sparse_ip, alg_klu; dt) - sol_sparse_op = solve(prob_sparse_op, alg_lu; dt) + sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg; dt) + sol_tridiagonal_op = solve(prob_tridiagonal_op, alg; dt) + sol_dense_ip = solve(prob_dense_ip, alg; dt) + sol_dense_op = solve(prob_dense_op, alg; dt) + sol_sparse_ip = solve(prob_sparse_ip, alg; dt) + sol_sparse_op = solve(prob_sparse_op, alg; dt) @test isapprox(sol_tridiagonal_ip.t, sol_tridiagonal_op.t; rtol) @test isapprox(sol_dense_ip.t, sol_dense_op.t; rtol) @@ -1525,16 +1511,11 @@ end tspan = (0.0, 1.0) dt = 0.25 - algs = (MPE, (; kwargs...) -> MPRK22(1.0; kwargs...), - (; kwargs...) -> MPRK22(0.5; kwargs...), - (; kwargs...) -> MPRK22(1.0; kwargs...), - (; kwargs...) -> MPRK43I(1.0, 0.5; kwargs...), - (; kwargs...) -> MPRK43I(0.5, 0.75; kwargs...), - (; kwargs...) -> MPRK43II(0.5; kwargs...), - (; kwargs...) -> MPRK43II(2.0 / 3.0; kwargs...), - (; kwargs...) -> SSPMPRK22(0.5, 1.0; kwargs...), - (; kwargs...) -> SSPMPRK43(; kwargs...)) - @testset "$alg" for alg in algs + @testset "$alg" for alg in (MPE(), + MPRK22(0.5), MPRK22(1.0), + MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), + MPRK43II(2.0 / 3.0), MPRK43II(0.5), + SSPMPRK22(0.5, 1.0), SSPMPRK43()) for (prod!, dest!) in zip((prod_1!, prod_2!, prod_3!), (dest_1!, dest_2!, dest_3!)) prod = (u, p, t) -> begin @@ -1560,19 +1541,17 @@ end prob_sparse_op = PDSProblem(prod, dest, u0, tspan; p_prototype = P_sparse) - alg_lu = alg(; linsolve = LUFactorization()) - alg_klu = alg(; linsolve = KLUFactorization()) - sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg_lu; + sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg; dt, adaptive = false) - sol_tridiagonal_op = solve(prob_tridiagonal_op, alg_lu; + sol_tridiagonal_op = solve(prob_tridiagonal_op, alg; dt, adaptive = false) - sol_dense_ip = solve(prob_dense_ip, alg_lu; + sol_dense_ip = solve(prob_dense_ip, alg; dt, adaptive = false) - sol_dense_op = solve(prob_dense_op, alg_lu; + sol_dense_op = solve(prob_dense_op, alg; dt, adaptive = false) - sol_sparse_ip = solve(prob_sparse_ip, alg_klu; + sol_sparse_ip = solve(prob_sparse_ip, alg; dt, adaptive = false) - sol_sparse_op = solve(prob_sparse_op, alg_lu; + sol_sparse_op = solve(prob_sparse_op, alg; dt, adaptive = false) @test sol_tridiagonal_ip.t ≈ sol_tridiagonal_op.t @@ -1658,16 +1637,11 @@ end dt = 0.25 rtol = sqrt(eps(Float32)) - algs = (MPE, (; kwargs...) -> MPRK22(1.0; kwargs...), - (; kwargs...) -> MPRK22(0.5; kwargs...), - (; kwargs...) -> MPRK22(1.0; kwargs...), - (; kwargs...) -> MPRK43I(1.0, 0.5; kwargs...), - (; kwargs...) -> MPRK43I(0.5, 0.75; kwargs...), - (; kwargs...) -> MPRK43II(0.5; kwargs...), - (; kwargs...) -> MPRK43II(2.0 / 3.0; kwargs...), - (; kwargs...) -> SSPMPRK22(0.5, 1.0; kwargs...), - (; kwargs...) -> SSPMPRK43(; kwargs...)) - @testset "$alg" for alg in algs + @testset "$alg" for alg in (MPE(), + MPRK22(0.5), MPRK22(1.0), + MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), + MPRK43II(2.0 / 3.0), MPRK43II(0.5), + SSPMPRK22(0.5, 1.0), SSPMPRK43()) for (prod!, dest!) in zip((prod_1!, prod_2!, prod_3!), (dest_1!, dest_2!, dest_3!)) prod! = prod_3! @@ -1695,20 +1669,12 @@ end prob_sparse_op = PDSProblem(prod, dest, u0, tspan; p_prototype = P_sparse) - alg_lu = alg(; linsolve = LUFactorization()) - alg_klu = alg(; linsolve = KLUFactorization()) - sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg_lu; - dt) - sol_tridiagonal_op = solve(prob_tridiagonal_op, alg_lu; - dt) - sol_dense_ip = solve(prob_dense_ip, alg_lu; - dt) - sol_dense_op = solve(prob_dense_op, alg_lu; - dt) - sol_sparse_ip = solve(prob_sparse_ip, alg_klu; - dt) - sol_sparse_op = solve(prob_sparse_op, alg_lu; - dt) + sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg; dt) + sol_tridiagonal_op = solve(prob_tridiagonal_op, alg; dt) + sol_dense_ip = solve(prob_dense_ip, alg; dt) + sol_dense_op = solve(prob_dense_op, alg; dt) + sol_sparse_ip = solve(prob_sparse_ip, alg; dt) + sol_sparse_op = solve(prob_sparse_op, alg; dt) @test isapprox(sol_tridiagonal_ip.t, sol_tridiagonal_op.t; rtol) @test isapprox(sol_dense_ip.t, sol_dense_op.t; rtol) @@ -1782,25 +1748,14 @@ end prob_sparse2 = PDSProblem(prod_sparse!, dest!, u0, tspan; p_prototype = P_sparse) # solve and test - algs = (MPE, (; kwargs...) -> MPRK22(1.0; kwargs...), - (; kwargs...) -> MPRK22(0.5; kwargs...), - (; kwargs...) -> MPRK22(1.0; kwargs...), - (; kwargs...) -> MPRK43I(1.0, 0.5; kwargs...), - (; kwargs...) -> MPRK43I(0.5, 0.75; kwargs...), - (; kwargs...) -> MPRK43II(0.5; kwargs...), - (; kwargs...) -> MPRK43II(2.0 / 3.0; kwargs...), - (; kwargs...) -> SSPMPRK22(0.5, 1.0; kwargs...), - (; kwargs...) -> SSPMPRK43(; kwargs...)) - for alg in algs + for alg in (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK43I(1.0, 0.5), + MPRK43I(0.5, 0.75), + MPRK43II(2.0 / 3.0), MPRK43II(0.5), SSPMPRK22(0.5, 1.0), + SSPMPRK43()) for prob in (prob_default, prob_tridiagonal, prob_dense, prob_sparse, prob_default2, prob_tridiagonal2, prob_dense2, prob_sparse2) - if prob == prob_sparse || prob == prob_sparse2 - alg_ = alg(; linsolve = KLUFactorization()) - else - alg_ = alg(; linsolve = LUFactorization()) - end - sol1 = solve(prob, alg_; dt, adaptive = false) + sol1 = solve(prob, alg; dt, adaptive = false) # test get_tmp_cache and integrator interface - modifying # values from the cache should not change the final results From e3d1ca9f46c43b1a99d8a9334ebb352f2a19b2bb Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Fri, 28 Mar 2025 12:01:03 +0100 Subject: [PATCH 25/26] fix --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 0cd16ab2..48be033e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1759,7 +1759,7 @@ end # test get_tmp_cache and integrator interface - modifying # values from the cache should not change the final results - integrator = init(prob, alg_; dt, adaptive = false) + integrator = init(prob, alg; dt, adaptive = false) PositiveIntegrators.OrdinaryDiffEqCore.step!(integrator) cache = @inferred PositiveIntegrators.get_tmp_cache(integrator) @test !isempty(cache) From fa4b682a3db0c62d26aed0554197fd10537e2e39 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Fri, 28 Mar 2025 12:42:16 +0100 Subject: [PATCH 26/26] drop support for LinearSolve.jl v2 --- Project.toml | 2 +- docs/Project.toml | 2 +- test/Project.toml | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 5e3fd71e..69008841 100644 --- a/Project.toml +++ b/Project.toml @@ -21,7 +21,7 @@ SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" [compat] FastBroadcast = "0.3.5" LinearAlgebra = "1" -LinearSolve = "2.39.1, 3.7.1" +LinearSolve = "3.7.1" MuladdMacro = "0.2.4" OrdinaryDiffEqCore = "1.16" RecipesBase = "1.3.4" diff --git a/docs/Project.toml b/docs/Project.toml index c413e644..d7fa0d15 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -27,7 +27,7 @@ DiffEqDevTools = "2.45.1" Documenter = "1" InteractiveUtils = "1" LinearAlgebra = "1" -LinearSolve = "2.39.1, 3.7.1" +LinearSolve = "3.7.1" OrdinaryDiffEqFIRK = "1.7" OrdinaryDiffEqLowOrderRK = "1.2" OrdinaryDiffEqRosenbrock = "1.4" diff --git a/test/Project.toml b/test/Project.toml index 832d4d94..65a43072 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -23,7 +23,7 @@ ADTypes = "1.12" Aqua = "0.8" ExplicitImports = "1.0.1" LinearAlgebra = "1" -LinearSolve = "2.39.1, 3.7.1" +LinearSolve = "3.7.1" OrdinaryDiffEqLowOrderRK = "1.2" OrdinaryDiffEqRosenbrock = "1.5" OrdinaryDiffEqSDIRK = "1.2"