Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Update StaticArray tests for eltype requirements #2087

Merged
merged 4 commits into from
Dec 24, 2023
Merged

Update StaticArray tests for eltype requirements #2087

Changes from 2 commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
393a506
Update StaticArray tests for eltype requirements
ChrisRackauckas Dec 17, 2023
f35f470
typo
ChrisRackauckas Dec 17, 2023
1ed5ccc
Update test/interface/static_array_tests.jl
ChrisRackauckas Dec 18, 2023
8296bc3
Update test/interface/static_array_tests.jl
ChrisRackauckas Dec 18, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 4 additions & 4 deletions test/interface/static_array_tests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@ using StaticArrays, Test
using OrdinaryDiffEq
using RecursiveArrayTools

u0 = [fill(2, MVector{2, Float64}), ones(MVector{2, Float64})]
u0 = VectorOfArray([fill(2, MVector{2, Float64}), ones(MVector{2, Float64})])
g(u, p, t) = SA[u[1] + u[2], u[1]]
f = (du, u, p, t) -> begin
for i in 1:2
Expand All @@ -17,7 +17,7 @@ sol = solve(ode, Tsit5())
sol = solve(ode, Vern9())
@test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0)

u0 = [fill(2, SVector{2, Float64}), ones(SVector{2, Float64})]
u0 = VectorOfArray([fill(2, SVector{2, Float64}), ones(SVector{2, Float64})])
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does this mean that plain Arrays of SVectors will not be supported anymore and we have to switch to VectorOfArray?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

See SciML/SciMLBase.jl#561. Basically I am trying to get things to a well-defined interface. This little bit sticks out as sore area of DiffEq where we don't throw any errors about arrays of arrays since 10 out of like 300 methods supports arrays of arrays. But we should be able to support such types of targets more broadly without overhead by just using a small wrapper that makes everything conform to the interface. If this is allowed, then we can throw informative errors as to what's going on and we can vastly clean up a lot of the code to not require hacks for specific methods to avoid assumptions about the lack of Number eltype. Since it's an undocumented bit only used downstream by you and @Datseris, it's a small discussion to be had with many good consequences to most users, and would make things like transitioning a code to implicit methods more readily available since the Vector of Vector form will likely never be able to support linear algebra appropriately.

As you can see from the PR, there is a keyword argument being made to bypass the interface erroring, but I hope that is only used during development and is not a crutch that is actually relied on for these kinds of tasks given that there should be zero overhead to using a more structured and interface-compliant form.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We use arrays of SVectors at least for some solvers in Trixi.jl with explicit Runge-Kutta methods. If OrdinaryDiffEq.jl transitions to another interface requirement, it would be great to get a breaking release describing this change.

We would also need more interfaces such as resize! and pointer for VectorOfArray.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As you can see from the PR, there is a keyword argument being made to bypass the interface erroring, but I hope that is only used during development and is not a crutch that is actually relied on for these kinds of tasks given that there should be zero overhead to using a more structured and interface-compliant form.

Which PR do you refer to? I could only see changes to the tests here and new docs at SciML/SciMLBase.jl#561

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure undocumented side pieces really constitute breaking. We are trying to impose such interfaces so that there is a strong sense of what is actually breaking.

We would also need more interfaces such as resize! and pointer for VectorOfArray.

Add suggestions in the interface PR and I'll happily accept that.

Copy link
Member

@avik-pal avik-pal Dec 19, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In BVDiffEq we were using a AbstractArray{<:AbstractArray} for specifying the initial guess on a mesh. We can upperbound SciMLBase there and cut a breaking release.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is bounding SciMLBase compat enough to address the impact of this change?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not SciMLBase my bad. DiffEqBase needs a upperbound. See SciML/BoundaryValueDiffEq.jl#152

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In BVDiffEq we were using a AbstractArray{<:AbstractArray} for specifying the initial guess on a mess. We can upperbound SciMLBase there and cut a breaking release.

BVDiffEq we need to specialize there because that's a legitimate use case, that's quite different.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How do you plan to do that? We may use it as inspiration for us in Trixi.jl.

ode = ODEProblem(f, u0, (0.0, 1.0))
sol = solve(ode, Euler(), dt = 1e-2)
@test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0)
Expand All @@ -28,14 +28,14 @@ sol = solve(ode, SSPRK22(), dt = 1e-2)
sol = solve(ode, ROCK4())
@test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0)

u0 = ones(MVector{2, Float64})
u0 = VectorOfArray(ones(MVector{2, Float64}))
ranocha marked this conversation as resolved.
Show resolved Hide resolved
ChrisRackauckas marked this conversation as resolved.
Show resolved Hide resolved
ode = ODEProblem(g, u0, (0.0, 1.0))
sol = solve(ode, Euler(), dt = 1e-2)
@test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0)
sol = solve(ode, Tsit5(), dt = 1e-2)
@test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0)

u0 = ones(SVector{2, Float64})
u0 = VectorOfArray(ones(SVector{2, Float64}))
ranocha marked this conversation as resolved.
Show resolved Hide resolved
ChrisRackauckas marked this conversation as resolved.
Show resolved Hide resolved
f = (u, p, t) -> u
ode = ODEProblem(f, u0, (0.0, 1.0))
sol = solve(ode, Euler(), dt = 1e-2)
Expand Down
Loading