ODE with heterogeneous input through RecursiveArrayTools

christianhvejsel · March 6, 2023, 5:53pm

Hi all
I am trying to solve a first-order ODE with initial conditions input as a heterogeneous array using of RecursiveArrayTools.jl https://docs.sciml.ai/DiffEqDocs/stable/features/diffeq_arrays/#ArrayPartitions.

I cannot figure out what causes this.

I have made an MWE


using OrdinaryDiffEq, LinearAlgebra, ForwardDiff, RecursiveArrayTools

A    = 0.225
σ_oe = 490e6 
C    = 1.25 
K    = 2.65e-5 
L    = 14.4

      
# Inline functions
s(x) = x - 1.0/3.0*tr(x)*diagm(ones(size(x,1)))  
σ_eff(x1, x2) = sqrt(3/2)*norm(s(x1)-x2,2)
β(x1, x2) = 1/σ_oe * (σ_eff(x1, x2) + A*tr(x1)-σ_oe)

function CD_ODE!(du,u,p,t)
    α .= u.x[1]
    D .= u.x[2]
    σ .= p[1](t)
    dσ .= p[2](t)

    dβ = 1/(σ_oe + C*σ_eff(σ, α))*(3/2*tr((s(σ).-α)*dσ')/σ_eff(σ, α) + A*tr(dσ))

    if  dβ > 0.0 && β(σ, α) > 0.0
        du.x[1] .= dβ*C*(s(σ).-α) 
        du.x[2] .= dβ*K*exp(L*β(σ, α)) 
    else
        du.x[1] .= 0.0
        du.x[2] .= 0.0
    end
end

α₀ = zeros(Float64, 3,3)
D₀ = 0.0::Float64
#u₀ = [α₀ D₀]
u₀ = ArrayPartition(α₀, D₀)


t_end = 100.
tspan = (0.0,t_end) 

f = 1.0

signal  = t -> 500e6*cos(2*pi*f*t)*diagm(ones(size(α₀,1)))
dsignal = t -> ForwardDiff.derivative(signal, t)

# p = [signal]
p = [signal dsignal]

prob = ODEProblem(CD_ODE!,u₀,tspan,p)
sol = solve(prob, Tsit5())

which gives me the following error:

ERROR: MethodError: no method matching similar(::Float64, ::Type{Float64})
Closest candidates are:
  similar(::Union{Adjoint{T, var"#s886"}, Transpose{T, var"#s886"}} where {T, var"#s886"<:(AbstractVector)}, ::Type{T}) where T at ~/.julia/juliaup/julia-1.8.5+0.aarch64.apple.darwin14/share/julia/stdlib/v1.8/LinearAlgebra/src/adjtrans.jl:207
  similar(::Union{Adjoint{T, S}, Transpose{T, S}} where {T, S}, ::Type{T}) where T at ~/.julia/juliaup/julia-1.8.5+0.aarch64.apple.darwin14/share/julia/stdlib/v1.8/LinearAlgebra/src/adjtrans.jl:211
  similar(::Union{Adjoint{T, S}, Transpose{T, S}} where {T, S}, ::Type{T}, ::Tuple{Vararg{Int64, N}}) where {T, N} at ~/.julia/juliaup/julia-1.8.5+0.aarch64.apple.darwin14/share/julia/stdlib/v1.8/LinearAlgebra/src/adjtrans.jl:212
  ...
Stacktrace:
  [1] (::RecursiveArrayTools.var"#52#53"{Float64, ArrayPartition{Float64, Tuple{Matrix{Float64}, Float64}}})(i::Int64)
    @ RecursiveArrayTools ~/.julia/packages/RecursiveArrayTools/9mtCv/src/array_partition.jl:59
  [2] ntuple
    @ ./ntuple.jl:49 [inlined]
  [3] ArrayPartition
    @ ~/.julia/packages/RecursiveArrayTools/9mtCv/src/array_partition.jl:31 [inlined]
  [4] similar
    @ ~/.julia/packages/RecursiveArrayTools/9mtCv/src/array_partition.jl:59 [inlined]
  [5] alg_cache(alg::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, u::ArrayPartition{Float64, Tuple{Matrix{Float64}, Float64}}, rate_prototype::ArrayPartition{Float64, Tuple{Matrix{Float64}, Float64}}, #unused#::Type{Float64}, #unused#::Type{Float64}, #unused#::Type{Float64}, uprev::ArrayPartition{Float64, Tuple{Matrix{Float64}, Float64}}, uprev2::ArrayPartition{Float64, Tuple{Matrix{Float64}, Float64}}, f::ODEFunction{true, SciMLBase.AutoSpecialize, typeof(CD_ODE!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, t::Float64, dt::Float64, reltol::Float64, p::Matrix{Function}, calck::Bool, #unused#::Val{true})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/vIqLB/src/caches/low_order_rk_caches.jl:522
  [6] __init(prob::ODEProblem{ArrayPartition{Float64, Tuple{Matrix{Float64}, Float64}}, Tuple{Float64, Float64}, true, Matrix{Function}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(CD_ODE!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Nothing, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Nothing, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{Int64}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/vIqLB/src/solve.jl:322
  [7] __init (repeats 5 times)
    @ ~/.julia/packages/OrdinaryDiffEq/vIqLB/src/solve.jl:10 [inlined]
  [8] __solve(::ODEProblem{ArrayPartition{Float64, Tuple{Matrix{Float64}, Float64}}, Tuple{Float64, Float64}, true, Matrix{Function}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(CD_ODE!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, ::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/vIqLB/src/solve.jl:5
  [9] __solve
    @ ~/.julia/packages/OrdinaryDiffEq/vIqLB/src/solve.jl:1 [inlined]
 [10] #solve_call#22
    @ ~/.julia/packages/DiffEqBase/jb48J/src/solve.jl:496 [inlined]
 [11] solve_call(_prob::ODEProblem{ArrayPartition{Float64, Tuple{Matrix{Float64}, Float64}}, Tuple{Float64, Float64}, true, Matrix{Function}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(CD_ODE!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/jb48J/src/solve.jl:466
 [12] solve_up(prob::ODEProblem{ArrayPartition{Float64, Tuple{Matrix{Float64}, Float64}}, Tuple{Float64, Float64}, true, Matrix{Function}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(CD_ODE!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::ArrayPartition{Float64, Tuple{Matrix{Float64}, Float64}}, p::Matrix{Function}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/jb48J/src/solve.jl:917
 [13] solve_up
    @ ~/.julia/packages/DiffEqBase/jb48J/src/solve.jl:890 [inlined]
 [14] solve(prob::ODEProblem{ArrayPartition{Float64, Tuple{Matrix{Float64}, Float64}}, Tuple{Float64, Float64}, true, Matrix{Function}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(CD_ODE!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{true}, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/jb48J/src/solve.jl:827
 [15] solve(prob::ODEProblem{ArrayPartition{Float64, Tuple{Matrix{Float64}, Float64}}, Tuple{Float64, Float64}, true, Matrix{Function}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(CD_ODE!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/jb48J/src/solve.jl:817
 [16] top-level scope
    @ ~/UserDir/MWE.jl:50

Any ideas for what’s going wrong here? Thanks!

ChrisRackauckas · March 6, 2023, 7:57pm

Try D₀ = [0.0]? Having one part just be a scalar might be the issue. I think they all need to be arrays if you want to mutate it.

christianhvejsel · March 6, 2023, 8:26pm

Thanks @ChrisRackauckas! Now it solves (…but doesn’t enter the non-trivial part of the conditions, but that’s for me to debug ).

I also modified these lines within the ODE function in the MWE:

    α = u.x[1]
    D = u.x[2]
    σ = p[1](t)
    dσ = p[2](t)