Closed mateuszbaran closed 6 months ago
ChrisRackauckas
commented
6 months ago I could write an MWE but at this point I think our use case should be just added to your test suite? It keeps failing so often and it's not too complex.
That would be good yes. It sounds like we're missing something.
mateuszbaran
commented
6 months ago Here is a full independent example:
using BoundaryValueDiffEq
struct EmbeddedTorus
R::Float64
r::Float64
end
function affine_connection!(M::EmbeddedTorus, Zc, i, a, Xc, Yc)
θ = a[1] .+ i[1]
sinθ, cosθ = sincos(θ)
Γ¹₂₂ = (M.R + M.r * cosθ) * sinθ / M.r
Γ²₁₂ = -M.r * sinθ / (M.R + M.r * cosθ)
Zc[1] = Xc[2] * Γ¹₂₂ * Yc[2]
Zc[2] = Γ²₁₂ * (Xc[1] * Yc[2] + Xc[2] * Yc[1])
return Zc
end
M = EmbeddedTorus(3, 2)
a1 = [0.5, -1.2]
a2 = [-0.5, 0.3]
i = (0, 0)
solver=MIRK4()
dt=0.05
tspan = (0.0, 1.0)
function bc1!(residual, u, p, t)
mid = div(length(u[1]), 2)
residual[1:mid] = u[1][1:mid] - a1
return residual[(mid + 1):end] = u[end][1:mid] - a2
end
function chart_log_problem!(du, u, params, t)
M, i = params
mid = div(length(u), 2)
a = u[1:mid]
dx = u[(mid + 1):end]
ddx = similar(dx)
affine_connection!(M, ddx, i, a, dx, dx)
ddx .*= -1
du[1:mid] .= dx
du[(mid + 1):end] .= ddx
return du
end
u0 = [vcat(a1, zero(a1)), vcat(a2, zero(a1))]
bvp1 = BVProblem(chart_log_problem!, bc1!, u0, tspan, (M, i))
sol1 = solve(bvp1, solver, dt=dt)I get the error locally with
[4fba245c] ArrayInterface v7.6.1
[764a87c0] BoundaryValueDiffEq v5.5.0
[2b5f629d] DiffEqBase v6.143.0
[459566f4] DiffEqCallbacks v2.35.0
⌃ [1dea7af3] OrdinaryDiffEq v6.59.3
⌅ [731186ca] RecursiveArrayTools v2.38.10
avik-pal
commented
6 months ago @mateuszbaran I will add this in our test suite. What is ddx supposed to be in affine_connection!(M, ddx, i, a, dx, dx)?
mateuszbaran
commented
6 months ago Thanks! This is essentially the geodesic equation on a torus in R^3, see https://en.wikipedia.org/wiki/Geodesic#Affine_geodesics . ddx is the second derivative in that equation. I'm using BoundaryValueDiffEq.jl to find the shortest geodesic between two points.
avik-pal
commented
6 months ago I meant is ddx = zeros(2)? it is missing in the code
mateuszbaran
commented
6 months ago ddx is filled by affine_connection!(M::EmbeddedTorus, Zc, i, a, Xc, Yc), where it's called Zc.
mateuszbaran
commented
6 months ago Ah, sorry, I understand now what you mean. It was supposed to be ddx = similar(dx). I've fixed that and tried re-running the example and now I get
Arrays with non-number element types, such as
`Array{Array{Float64}}`, are not supported by the
solvers.which I guess is a breaking change introduced at some point recently without making a breaking release?
avik-pal
commented
6 months ago I am fixing that. Basically the order of fixes will be:
VectorOfArray(...)mateuszbaran
commented
6 months ago I see, thanks!
avik-pal
commented
6 months ago @mateuszbaran I will also add a downstream test on Manifolds.jl. Is there a particular group I should test https://github.com/JuliaManifolds/Manifolds.jl/blob/b2d83b337e214158f64d22a39776b9a80ffd2981/test/utils.jl#L4 for boundaryvaluediffeq?
mateuszbaran
commented
6 months ago test_integration group has the test for BoundaryValueDiffEq but currently the example I've pasted above is the only use case we have for BVP and we don't have any plans to add new functionality relying on BVP. Adding the test case I've posted above should be enough IMO.
BoundaryValueDiffEq.jl is again causing test failures in Manifolds.jl:
See here: https://github.com/JuliaManifolds/Manifolds.jl/actions/runs/7211150273/job/19646434007?pr=692 . I could write an MWE but at this point I think our use case should be just added to your test suite? It keeps failing so often and it's not too complex.