Creating `BVProblem` from `ODESystem` #3251

vyudu · 2024-12-01T21:25:49Z

WIP PR for the BVProblem conversion.

…ition

vyudu · 2024-12-02T02:57:56Z

for some reason trying to solve the BVProblem constructed from the Lorenz IVP with MIRK4() is not respecting the boundary conditions (the u[1] is something like [0.911, -0.073, 0.001] instead of [1., 0., 0.]). Not sure if this bug is in the way I'm constructing the problem

ChrisRackauckas · 2024-12-03T01:33:57Z

Lorenz will be a very difficult one to solve because the BVP solver is going to correct global error rather than local, and getting good global error on a chaotic system is almost impossible 😅. Switch to something like Lotka Volterra.

src/systems/diffeqs/abstractodesystem.jl

ChrisRackauckas · 2024-12-03T19:05:59Z

This looks right. How does it do?

vyudu · 2024-12-03T20:07:28Z

Getting there, but I'm getting this error when I try to solve the pendulum out-of-place BVP (it doesn't happen for the Lotka-Volterra out-of-place BVP).

julia> sol2 = solve(bvp2, MIRK4(), dt = 0.01);
ERROR: MethodError: no method matching create_array(::Type{…}, ::Nothing, ::Val{…}, ::Val{…}, ::ForwardDiff.Dual{…}, ::ForwardDiff.Dual{…})
The function `create_array` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  create_array(::Type{<:SubArray{T, N, P, I, L}}, ::Any, ::Val, ::Val, ::Any...) where {T, N, P, I, L}
   @ SymbolicUtils ~/.julia/packages/SymbolicUtils/jf8aQ/src/code.jl:546
  create_array(::Type{<:Array}, ::Any, ::Val, ::Val, ::Any...)
   @ SymbolicUtils ~/.julia/packages/SymbolicUtils/jf8aQ/src/code.jl:498
  create_array(::Type{<:Array}, ::Nothing, ::Val, ::Val{dims}, ::Any...) where dims
   @ SymbolicUtils ~/.julia/packages/SymbolicUtils/jf8aQ/src/code.jl:502
  ...

Stacktrace:
  [1] macro expansion
    @ ~/.julia/packages/SymbolicUtils/jf8aQ/src/code.jl:480 [inlined]
  [2] macro expansion
    @ ~/.julia/packages/RuntimeGeneratedFunctions/M9ZX8/src/RuntimeGeneratedFunctions.jl:163 [inlined]
  [3] macro expansion
    @ ./none:0 [inlined]
  [4] generated_callfunc
    @ ./none:0 [inlined]
  [5] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{…})(::Base.ReinterpretArray{…}, ::Vector{…}, ::Float64)
    @ RuntimeGeneratedFunctions ~/.julia/packages/RuntimeGeneratedFunctions/M9ZX8/src/RuntimeGeneratedFunctions.jl:150
  [6] (::ModelingToolkit.var"#f#946"{…})(u::Base.ReinterpretArray{…}, p::MTKParameters{…}, t::Float64)
    @ ModelingToolkit ~/.julia/dev/ModelingToolkit/src/systems/diffeqs/abstractodesystem.jl:381
  [7] (::ODEFunction{…})(::Base.ReinterpretArray{…}, ::Vararg{…})
    @ SciMLBase ~/.julia/packages/SciMLBase/ZyZAV/src/scimlfunctions.jl:2358
  [8] (::BVPFunction{…})(::Base.ReinterpretArray{…}, ::Vararg{…})
    @ SciMLBase ~/.julia/packages/SciMLBase/ZyZAV/src/scimlfunctions.jl:2411
  [9] Φ(fᵢ_cache::PreallocationTools.DiffCache{…}, k_discrete::Vector{…}, f::BVPFunction{…}, TU::BoundaryValueDiffEqMIRK.MIRKTableau{…}, y::Vector{…}, u::Vector{…}, p::MTKParameters{…}, mesh::Vector{…}, mesh_dt::Vector{…}, stage::Int64)
    @ BoundaryValueDiffEqMIRK ~/.julia/packages/BoundaryValueDiffEqMIRK/wGmgf/src/collocation.jl:54
 [10] Φ
    @ ~/.julia/packages/BoundaryValueDiffEqMIRK/wGmgf/src/collocation.jl:33 [inlined]
 [11] __mirk_loss_collocation
    @ ~/.julia/packages/BoundaryValueDiffEqMIRK/wGmgf/src/mirk.jl:302 [inlined]
 [12] #43
    @ ~/.julia/packages/BoundaryValueDiffEqMIRK/wGmgf/src/mirk.jl:223 [inlined]
 [13] Fix2
    @ ./operators.jl:1144 [inlined]
 [14] forwarddiff_color_jacobian(J::BandedMatrices.BandedMatrix{…}, f::Base.Fix2{…}, x::Vector{…}, jac_cache::ForwardColorJacCache{…})
    @ SparseDiffTools ~/.julia/packages/SparseDiffTools/m4eae/src/differentiation/compute_jacobian_ad.jl:205
 [15] sparse_jacobian!
    @ ~/.julia/packages/SparseDiffTools/m4eae/src/highlevel/forward_mode.jl:53 [inlined]
 [16] __mirk_mpoint_jacobian(J::FastAlmostBandedMatrices.AlmostBandedMatrix{…}, J_c::BandedMatrices.BandedMatrix{…}, x::Vector{…}, bc_diffmode::AutoForwardDiff{…}, nonbc_diffmode::AutoSparse{…}, bc_diffcache::SparseDiffTools.ForwardDiffJacobianCache{…}, nonbc_diffcache::SparseDiffTools.ForwardDiffJacobianCache{…}, loss_bc::Base.Fix2{…}, loss_collocation::Base.Fix2{…}, L::Int64)
    @ BoundaryValueDiffEqMIRK ~/.julia/packages/BoundaryValueDiffEqMIRK/wGmgf/src/mirk.jl:406
 [17] #59
    @ ~/.julia/packages/BoundaryValueDiffEqMIRK/wGmgf/src/mirk.jl:359 [inlined]
 [18] JacobianCache
    @ ~/.julia/packages/NonlinearSolveBase/sduQ7/src/jacobian.jl:170 [inlined]
 [19] JacobianCache
    @ ~/.julia/packages/NonlinearSolveBase/sduQ7/src/jacobian.jl:130 [inlined]
 [20] step!(cache::NonlinearSolveFirstOrder.GeneralizedFirstOrderAlgorithmCache{…}; recompute_jacobian::Nothing)
    @ NonlinearSolveFirstOrder ~/.julia/packages/NonlinearSolveFirstOrder/XQs9h/src/solve.jl:224
 [21] step!
    @ ~/.julia/packages/NonlinearSolveFirstOrder/XQs9h/src/solve.jl:218 [inlined]
 [22] #step!#161
    @ ~/.julia/packages/NonlinearSolveBase/sduQ7/src/solve.jl:219 [inlined]
 [23] step!
    @ ~/.julia/packages/NonlinearSolveBase/sduQ7/src/solve.jl:213 [inlined]
 [24] solve!(cache::NonlinearSolveFirstOrder.GeneralizedFirstOrderAlgorithmCache{…})
    @ NonlinearSolveBase ~/.julia/packages/NonlinearSolveBase/sduQ7/src/solve.jl:11
 [25] #__solve#144
    @ ~/.julia/packages/NonlinearSolveBase/sduQ7/src/solve.jl:6 [inlined]
 [26] __solve
    @ ~/.julia/packages/NonlinearSolveBase/sduQ7/src/solve.jl:1 [inlined]
 [27] macro expansion
    @ ~/.julia/packages/NonlinearSolveBase/sduQ7/src/solve.jl:139 [inlined]
 [28] __generated_polysolve(::NonlinearProblem{…}, ::NonlinearSolveBase.NonlinearSolvePolyAlgorithm{…}; stats::SciMLBase.NLStats, alias_u0::Bool, verbose::Bool, kwargs::@Kwargs{…})
    @ NonlinearSolveBase ~/.julia/packages/NonlinearSolveBase/sduQ7/src/solve.jl:111
 [29] __generated_polysolve
    @ ~/.julia/packages/NonlinearSolveBase/sduQ7/src/solve.jl:111 [inlined]
 [30] #__solve#152
    @ ~/.julia/packages/NonlinearSolveBase/sduQ7/src/solve.jl:108 [inlined]
 [31] __perform_mirk_iteration(cache::BoundaryValueDiffEqMIRK.MIRKCache{…}, abstol::Float64, adaptive::Bool; nlsolve_kwargs::@NamedTuple{}, kwargs::@Kwargs{…})
    @ BoundaryValueDiffEqMIRK ~/.julia/packages/BoundaryValueDiffEqMIRK/wGmgf/src/mirk.jl:163
 [32] __perform_mirk_iteration
    @ ~/.julia/packages/BoundaryValueDiffEqMIRK/wGmgf/src/mirk.jl:159 [inlined]
 [33] solve!(cache::BoundaryValueDiffEqMIRK.MIRKCache{…})
    @ BoundaryValueDiffEqMIRK ~/.julia/packages/BoundaryValueDiffEqMIRK/wGmgf/src/mirk.jl:140
 [34] #__solve#29
    @ ~/.julia/packages/BoundaryValueDiffEqCore/SKdo9/src/BoundaryValueDiffEqCore.jl:34 [inlined]
 [35] __solve
    @ ~/.julia/packages/BoundaryValueDiffEqCore/SKdo9/src/BoundaryValueDiffEqCore.jl:32 [inlined]
 [36] #solve_call#44
    @ ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:632 [inlined]
 [37] solve_call
    @ ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:589 [inlined]
 [38] solve_up(prob::BVProblem{…}, sensealg::Nothing, u0::Vector{…}, p::MTKParameters{…}, args::MIRK4{…}; kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:1120
 [39] solve_up
    @ ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:1099 [inlined]
 [40] solve(prob::BVProblem{…}, args::MIRK4{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:1036
 [41] top-level scope
    @ REPL[147]:1
Some type information was truncated. Use `show(err)` to see complete types.

vyudu · 2024-12-03T20:59:00Z

Seems like it's reinterpreting the array but I can't figure out why

ChrisRackauckas · 2024-12-03T23:16:41Z

Reinterpreting which array?

vyudu · 2024-12-04T00:10:38Z

I'm not sure which array this corresponds to, but it looks like this line is returning tmp as a ReinterpretArray, which is causing issues when f is called on it later.

vyudu · 2024-12-16T16:00:57Z

@ChrisRackauckas Think this should be ready now. I was looking totally in the wrong place, the bug was just that create_array in SymbolicUtils.Code has no method for a ReinterpretArray. I added the method here but if it would make sense to change create_array to work for anything <: AbstractArray instead of just <: Array or something like that I can make a PR to SymbolicUtils as well.

test/bvproblem.jl

ChrisRackauckas · 2024-12-18T15:23:10Z

src/systems/diffeqs/abstractodesystem.jl

+@inline function SymbolicUtils.Code.create_array(
+        ::Type{<:Base.ReinterpretArray}, T, ::Val{1}, ::Val{dims}, elems...) where {dims}
+    T[elems...]
+end


This is piracy. Why is there a ReinterpretArray in the first place?

I haven't understood the internals well enough to really know what's going on but fetching from the DiffCache is returning a ReinterpretArray of Dual numbers in the line I linked above. Should that not happen?

Oh that would cause a reinterpret array. Okay, upstream this.

SciML/PreallocationTools.jl#114

I meant this bit of code should go to SymbolicUtils

src/systems/diffeqs/abstractodesystem.jl

vyudu added 9 commits October 28, 2024 17:14

make Observable definition comply with new record_from_solution defin…

add7059

…ition

add test

eb681c1

init

9733460

Merge remote-tracking branch 'origin/master' into MTK

d95e4a7

up

b3da813

Merge remote-tracking branch 'origin/master' into MTK

86c82ce

up

4affeac

up

a3429ea

up

f751fbb

up

a9fdfd6

ChrisRackauckas reviewed Dec 3, 2024

View reviewed changes

src/systems/diffeqs/abstractodesystem.jl Outdated Show resolved Hide resolved

ChrisRackauckas reviewed Dec 3, 2024

View reviewed changes

src/systems/diffeqs/abstractodesystem.jl Outdated Show resolved Hide resolved

vyudu added 5 commits December 3, 2024 19:11

up

a9f2106

up

18fdd5f

fixing create_array

9d65a33

revert Project.toml

999ec30

Up

9226ad6

vyudu marked this pull request as ready for review December 16, 2024 15:56

vyudu added 4 commits December 17, 2024 00:03

Merge remote-tracking branch 'origin/master' into MTK

0cb4893

formatting

67d8164

up

25988f3

up

bb28d4f

fix

b2bf7c0

ChrisRackauckas reviewed Dec 18, 2024

View reviewed changes

test/bvproblem.jl Outdated Show resolved Hide resolved

ChrisRackauckas reviewed Dec 18, 2024

View reviewed changes

src/systems/diffeqs/abstractodesystem.jl Outdated Show resolved Hide resolved

ChrisRackauckas reviewed Dec 18, 2024

View reviewed changes

src/systems/diffeqs/abstractodesystem.jl Outdated Show resolved Hide resolved

up

3751c2a

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Creating `BVProblem` from `ODESystem` #3251

Creating `BVProblem` from `ODESystem` #3251

vyudu commented Dec 1, 2024

vyudu commented Dec 2, 2024

ChrisRackauckas commented Dec 3, 2024

ChrisRackauckas commented Dec 3, 2024

vyudu commented Dec 3, 2024

vyudu commented Dec 3, 2024

ChrisRackauckas commented Dec 3, 2024

vyudu commented Dec 4, 2024

vyudu commented Dec 16, 2024

ChrisRackauckas Dec 18, 2024

vyudu Dec 20, 2024

ChrisRackauckas Dec 20, 2024

vyudu Dec 20, 2024

ChrisRackauckas Dec 20, 2024

vyudu Dec 20, 2024

Creating BVProblem from ODESystem #3251

Are you sure you want to change the base?

Creating BVProblem from ODESystem #3251

Conversation

vyudu commented Dec 1, 2024

vyudu commented Dec 2, 2024

ChrisRackauckas commented Dec 3, 2024

ChrisRackauckas commented Dec 3, 2024

vyudu commented Dec 3, 2024

vyudu commented Dec 3, 2024

ChrisRackauckas commented Dec 3, 2024

vyudu commented Dec 4, 2024

vyudu commented Dec 16, 2024

ChrisRackauckas Dec 18, 2024

Choose a reason for hiding this comment

vyudu Dec 20, 2024

Choose a reason for hiding this comment

ChrisRackauckas Dec 20, 2024

Choose a reason for hiding this comment

vyudu Dec 20, 2024

Choose a reason for hiding this comment

ChrisRackauckas Dec 20, 2024

Choose a reason for hiding this comment

vyudu Dec 20, 2024

Choose a reason for hiding this comment

Creating `BVProblem` from `ODESystem` #3251

Creating `BVProblem` from `ODESystem` #3251