Array could not broadcast to match destination in solve

General Usage Performance
1

@ChrisRackauckas

Here is a MWE solving a simple 3x3 equation, where I use the SizedVector from StaticArrays.jl:

using  DifferentialEquations
using StaticArrays

# Set up the Giesekus model
function dudt!(du, σ, p, t)
    du .= -σ 
end

η0 = 1
τ = 1
α = 0.8
p_giesekus = [η0,τ,α]
σ0 = SizedMatrix{3,3}([1. 0. 0. ; 0. 2. 0. ; 3. 0. 0.])

fct(t) = 1*cos(t)

prob_giesekus = ODEProblem(dudt!, σ0, (0., 2.), p_giesekus)
solve_giesekus = solve(prob_giesekus, Rodas4(), saveat=0.2)  # Broadcast error
solve_giesekus = solve(prob_giesekus, Tsit5(), saveat=0.2)  # Runs fine

Notice the line that says Broadcast error. Why do I get an error?

I feel like I am not using StaticArrays correctly. The MWE is a greatly reduced version of a code which has been giving strange errors. I am solving a UODE and was using static arrays. When using Rodas4, I think that autodifferentiatpion is used, and therefore the types stored in the static arrays were no long bit-compatible variables (ints, floats, etc.) and I got an error. So I changed to mutable arrays SizedArray since my array sizes are almost `9x9. Thanks for any help!

Here is the stack trace:

ERROR: DimensionMismatch: array could not be broadcast to match destination
Stacktrace:
  [1] check_broadcast_shape(shp::Tuple{SOneTo{3}, SOneTo{3}}, Ashp::Tuple{SOneTo{9}}) @ Base.Broadcast ./broadcast.jl:553
  [2] check_broadcast_axes @ ./broadcast.jl:556 [inlined]
  [3] instantiate @ ./broadcast.jl:297 [inlined]
  [4] instantiate @ ~/.julia/packages/StaticArrays/VLqRb/src/broadcast.jl:32 [inlined]
  [5] materialize! @ ./broadcast.jl:884 [inlined]
  [6] materialize! @ ./broadcast.jl:881 [inlined]
  [7] fast_materialize! @ ~/.julia/packages/FastBroadcast/xrPat/src/FastBroadcast.jl:47 [inlined]
  [8] perform_step!(integrator::OrdinaryDiffEq.ODEIntegrator{Rodas4{9, true, LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, true, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Nothing, Float64, Vector{Float64}, Float64, Float64, Float64, Float64, Vector{SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}}, ODESolution{Float64, 3, Vector{SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}}}, ODEProblem{SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(dudt!), LinearAlgebra.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}, Rodas4{9, true, LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.AutoSpecialize, typeof(dudt!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}}, Vector{Float64}, Vector{Vector{SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}}}, OrdinaryDiffEq.Rodas4Cache{SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.RodasTableau{Float64, Float64}, SciMLBase.TimeGradientWrapper{ODEFunction{true, SciMLBase.AutoSpecialize, typeof(dudt!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Vector{Float64}}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.AutoSpecialize, typeof(dudt!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, Vector{Float64}}, LinearSolve.LinearCache{Matrix{Float64}, SizedVector{9, Float64, Vector{Float64}}, SizedVector{9, Float64, Vector{Float64}}, SciMLBase.NullParameters, LUFactorization{LinearAlgebra.RowMaximum}, LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, SizedVector{9, Float64, Vector{Float64}}}}, LinearAlgebra.Diagonal{Float64, SizedVector{9, Float64, Vector{Float64}}}, Float64, true, LinearSolve.OperatorCondition.IllConditioned}, SparseDiffTools.ForwardColorJacCache{SizedMatrix{3, 3, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}, 2, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}}}, SizedMatrix{3, 3, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}, 2, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}}}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Vector{Vector{NTuple{9, Float64}}}, UnitRange{Int64}, Nothing}, SizedMatrix{3, 3, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}, 2, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, Float64, Rodas4{9, true, LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}}, DiffEqBase.Stats, Nothing}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(dudt!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, OrdinaryDiffEq.Rodas4Cache{SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.RodasTableau{Float64, Float64}, SciMLBase.TimeGradientWrapper{ODEFunction{true, SciMLBase.AutoSpecialize, typeof(dudt!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Vector{Float64}}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.AutoSpecialize, typeof(dudt!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, Vector{Float64}}, LinearSolve.LinearCache{Matrix{Float64}, SizedVector{9, Float64, Vector{Float64}}, SizedVector{9, Float64, Vector{Float64}}, SciMLBase.NullParameters, LUFactorization{LinearAlgebra.RowMaximum}, LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, SizedVector{9, Float64, Vector{Float64}}}}, LinearAlgebra.Diagonal{Float64, SizedVector{9, Float64, Vector{Float64}}}, Float64, true, LinearSolve.OperatorCondition.IllConditioned}, SparseDiffTools.ForwardColorJacCache{SizedMatrix{3, 3, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}, 2, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}}}, SizedMatrix{3, 3, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}, 2, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}}}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Vector{Vector{NTuple{9, Float64}}}, UnitRange{Int64}, Nothing}, SizedMatrix{3, 3, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}, 2, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, Float64, Rodas4{9, true, LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Float64, Tuple{}}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.Rodas4Cache{SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.RodasTableau{Float64, Float64}, SciMLBase.TimeGradientWrapper{ODEFunction{true, SciMLBase.AutoSpecialize, typeof(dudt!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Vector{Float64}}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.AutoSpecialize, typeof(dudt!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, Vector{Float64}}, LinearSolve.LinearCache{Matrix{Float64}, SizedVector{9, Float64, Vector{Float64}}, SizedVector{9, Float64, Vector{Float64}}, SciMLBase.NullParameters, LUFactorization{LinearAlgebra.RowMaximum}, LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, SizedVector{9, Float64, Vector{Float64}}}}, LinearAlgebra.Diagonal{Float64, SizedVector{9, Float64, Vector{Float64}}}, Float64, true, LinearSolve.OperatorCondition.IllConditioned}, SparseDiffTools.ForwardColorJacCache{SizedMatrix{3, 3, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}, 2, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}}}, SizedMatrix{3, 3, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}, 2, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}}}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Vector{Vector{NTuple{9, Float64}}}, UnitRange{Int64}, Nothing}, SizedMatrix{3, 3, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}, 2, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, Float64, Rodas4{9, true, LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, repeat_step::Bool) @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/xs0Fk/src/perform_step/rosenbrock_perform_step.jl:1084
  [9] perform_step!(integrator::OrdinaryDiffEq.ODEIntegrator{Rodas4{9, true, LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, true, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Nothing, Float64, Vector{Float64}, Float64, Float64, Float64, Float64, Vector{SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}}, ODESolution{Float64, 3, Vector{SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}}}, ODEProblem{SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(dudt!), LinearAlgebra.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}, Rodas4{9, true, LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.AutoSpecialize, typeof(dudt!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}}, Vector{Float64}, Vector{Vector{SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}}}, OrdinaryDiffEq.Rodas4Cache{SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.RodasTableau{Float64, Float64}, SciMLBase.TimeGradientWrapper{ODEFunction{true, SciMLBase.AutoSpecialize, typeof(dudt!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Vector{Float64}}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.AutoSpecialize, typeof(dudt!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, Vector{Float64}}, LinearSolve.LinearCache{Matrix{Float64}, SizedVector{9, Float64, Vector{Float64}}, SizedVector{9, Float64, Vector{Float64}}, SciMLBase.NullParameters, LUFactorization{LinearAlgebra.RowMaximum}, LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, SizedVector{9, Float64, Vector{Float64}}}}, LinearAlgebra.Diagonal{Float64, SizedVector{9, Float64, Vector{Float64}}}, Float64, true, LinearSolve.OperatorCondition.IllConditioned}, SparseDiffTools.ForwardColorJacCache{SizedMatrix{3, 3, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}, 2, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}}}, SizedMatrix{3, 3, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}, 2, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}}}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Vector{Vector{NTuple{9, Float64}}}, UnitRange{Int64}, Nothing}, SizedMatrix{3, 3, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}, 2, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, Float64, Rodas4{9, true, LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}}, DiffEqBase.Stats, Nothing}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(dudt!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, OrdinaryDiffEq.Rodas4Cache{SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.RodasTableau{Float64, Float64}, SciMLBase.TimeGradientWrapper{ODEFunction{true, SciMLBase.AutoSpecialize, typeof(dudt!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Vector{Float64}}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.AutoSpecialize, typeof(dudt!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, Vector{Float64}}, LinearSolve.LinearCache{Matrix{Float64}, SizedVector{9, Float64, Vector{Float64}}, SizedVector{9, Float64, Vector{Float64}}, SciMLBase.NullParameters, LUFactorization{LinearAlgebra.RowMaximum}, LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, SizedVector{9, Float64, Vector{Float64}}}}, LinearAlgebra.Diagonal{Float64, SizedVector{9, Float64, Vector{Float64}}}, Float64, true, LinearSolve.OperatorCondition.IllConditioned}, SparseDiffTools.ForwardColorJacCache{SizedMatrix{3, 3, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}, 2, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}}}, SizedMatrix{3, 3, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}, 2, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}}}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Vector{Vector{NTuple{9, Float64}}}, UnitRange{Int64}, Nothing}, SizedMatrix{3, 3, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}, 2, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, Float64, Rodas4{9, true, LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Float64, Tuple{}}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.Rodas4Cache{SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.RodasTableau{Float64, Float64}, SciMLBase.TimeGradientWrapper{ODEFunction{true, SciMLBase.AutoSpecialize, typeof(dudt!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Vector{Float64}}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.AutoSpecialize, typeof(dudt!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, Vector{Float64}}, LinearSolve.LinearCache{Matrix{Float64}, SizedVector{9, Float64, Vector{Float64}}, SizedVector{9, Float64, Vector{Float64}}, SciMLBase.NullParameters, LUFactorization{LinearAlgebra.RowMaximum}, LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, SizedVector{9, Float64, Vector{Float64}}}}, LinearAlgebra.Diagonal{Float64, SizedVector{9, Float64, Vector{Float64}}}, Float64, true, LinearSolve.OperatorCondition.IllConditioned}, SparseDiffTools.ForwardColorJacCache{SizedMatrix{3, 3, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}, 2, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}}}, SizedMatrix{3, 3, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}, 2, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}}}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Vector{Vector{NTuple{9, Float64}}}, UnitRange{Int64}, Nothing}, SizedMatrix{3, 3, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}, 2, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, Float64, Rodas4{9, true, LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}) @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/xs0Fk/src/perform_step/rosenbrock_perform_step.jl:1033
 [10] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{Rodas4{9, true, LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, true, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Nothing, Float64, Vector{Float64}, Float64, Float64, Float64, Float64, Vector{SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}}, ODESolution{Float64, 3, Vector{SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}}}, ODEProblem{SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(dudt!), LinearAlgebra.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}, Rodas4{9, true, LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.AutoSpecialize, typeof(dudt!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}}, Vector{Float64}, Vector{Vector{SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}}}, OrdinaryDiffEq.Rodas4Cache{SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.RodasTableau{Float64, Float64}, SciMLBase.TimeGradientWrapper{ODEFunction{true, SciMLBase.AutoSpecialize, typeof(dudt!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Vector{Float64}}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.AutoSpecialize, typeof(dudt!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, Vector{Float64}}, LinearSolve.LinearCache{Matrix{Float64}, SizedVector{9, Float64, Vector{Float64}}, SizedVector{9, Float64, Vector{Float64}}, SciMLBase.NullParameters, LUFactorization{LinearAlgebra.RowMaximum}, LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, SizedVector{9, Float64, Vector{Float64}}}}, LinearAlgebra.Diagonal{Float64, SizedVector{9, Float64, Vector{Float64}}}, Float64, true, LinearSolve.OperatorCondition.IllConditioned}, SparseDiffTools.ForwardColorJacCache{SizedMatrix{3, 3, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}, 2, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}}}, SizedMatrix{3, 3, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}, 2, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}}}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Vector{Vector{NTuple{9, Float64}}}, UnitRange{Int64}, Nothing}, SizedMatrix{3, 3, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}, 2, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, Float64, Rodas4{9, true, LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}}, DiffEqBase.Stats, Nothing}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(dudt!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, OrdinaryDiffEq.Rodas4Cache{SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.RodasTableau{Float64, Float64}, SciMLBase.TimeGradientWrapper{ODEFunction{true, SciMLBase.AutoSpecialize, typeof(dudt!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Vector{Float64}}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.AutoSpecialize, typeof(dudt!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, Vector{Float64}}, LinearSolve.LinearCache{Matrix{Float64}, SizedVector{9, Float64, Vector{Float64}}, SizedVector{9, Float64, Vector{Float64}}, SciMLBase.NullParameters, LUFactorization{LinearAlgebra.RowMaximum}, LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, SizedVector{9, Float64, Vector{Float64}}}}, LinearAlgebra.Diagonal{Float64, SizedVector{9, Float64, Vector{Float64}}}, Float64, true, LinearSolve.OperatorCondition.IllConditioned}, SparseDiffTools.ForwardColorJacCache{SizedMatrix{3, 3, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}, 2, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}}}, SizedMatrix{3, 3, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}, 2, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}}}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Vector{Vector{NTuple{9, Float64}}}, UnitRange{Int64}, Nothing}, SizedMatrix{3, 3, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}, 2, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, Float64, Rodas4{9, true, LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Float64, Tuple{}}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}) @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/xs0Fk/src/solve.jl:520
 [11] #__solve#640 @ ~/.julia/packages/OrdinaryDiffEq/xs0Fk/src/solve.jl:6 [inlined]
 [12] __solve @ ~/.julia/packages/OrdinaryDiffEq/xs0Fk/src/solve.jl:1 [inlined]
 [13] #solve_call#22 @ ~/.julia/packages/DiffEqBase/ihYDa/src/solve.jl:509 [inlined]
 [14] solve_call @ ~/.julia/packages/DiffEqBase/ihYDa/src/solve.jl:479 [inlined]
 [15] #solve_up#29 @ ~/.julia/packages/DiffEqBase/ihYDa/src/solve.jl:932 [inlined]
 [16] solve_up @ ~/.julia/packages/DiffEqBase/ihYDa/src/solve.jl:905 [inlined]
 [17] #solve#27 @ ~/.julia/packages/DiffEqBase/ihYDa/src/solve.jl:842 [inlined]
 [18] top-level scope @ ~/src/2022/rude/giesekus/GE_rude.jl/alex_report_code_2023-03-24/julia19/giesekus_opt_3x3:24

I am running Julia 1.9. The error does not occur when using Tsit5.
More specifically, the error occurs in :

 perform_step!(integrator::OrdinaryDiffEq.ODEIntegrator{Rodas4{9, true, LUFactor ...
2

Hey, not your fault. It seems there was a dispatch that was missing a vec so it was throwing a bounds error out of safety because the sizes didn’t align internally. It seems like it only shows up for Rodas4, and I tested a bunch of cases so it was just that one dispatch.

I’m just going to copy the text of the PR which fixes this:

This fixes a case where you have an AbstractMatrix which was not vec resized so it was throwing an error in broadcast. Normally this is covered by the convergence tests, since most of the convergence tests use matrices for u0. However, for a subset of algorithms there is a separate <:Array dispatch for compile time performance. For those cases, it was found that in one fallback <:AbstractArray dispatch that the _vecs were missing, making a size issue between the linear solver’s _vec(u) vs k1.

This is fixed, and a general test on all of the algorithms with such a compile time improvement fallback was added:

using OrdinaryDiffEq, StaticArrays, Test

function dudt!(du, σ, p, t)
    du .= -σ 
end

η0 = 1
τ = 1
α = 0.8
p_giesekus = [η0,τ,α]

## Exercise code paths that are in-place but not `<:Array` but not AbstractVector
σ0 = SizedMatrix{3,3}([1. 0. 0. ; 0. 2. 0. ; 3. 0. 0.])

prob_giesekus = ODEProblem(dudt!, σ0, (0., 2.), p_giesekus)
solve_giesekus = solve(prob_giesekus, Rodas4(), saveat=0.2, abstol=1e-14, reltol=1e-14) 

for alg in [Rosenbrock23(), Rodas4(), Rodas4P(), Rodas5(), Rodas5P(), Tsit5(), Vern6(), Vern7(), Vern8(), Vern9(), DP5()]
    sol = solve(prob_giesekus, alg, saveat=0.2, abstol=1e-14, reltol=1e-14)
    @test Array(sol) ≈ Array(solve_giesekus)
end

It should be merged in a few hours and be in a patch by the end of the day.

3

Thank you, Chris. Now I continue to debug. Unless my issues are fully solved. :slight_smile:

4

Hi @ChrisRackauckas,

Hopefully you are around :slight_smile:

I return to a previous error. I just put a few lines from the stack trace, because you might be able to see what is happening. I will put my code below the stack trace snippet.

σ: [10.844265263548873 0.0 0.0; 0.0 11.393423775193504 0.0; 1.647475534933883 0.0 10.29510675190425]
tr(σ): 32.53279579064663
==> g: Float32[1.3600748, -0.08879593]
==> F: [0.3971481376901842 0.0 0.0; 0.0 0.34838509685046315 0.0; -0.1462891225191623 0.0 0.4459111785299046]
σ: ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 2, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}, SizedMatrix{3, 3, Float64, 2, Matrix{Float64}}}}[TrackedReal<KGh>(10.844265263548873, 0.0, E3M, 1, LoY) TrackedReal<7nK>(0.0, 0.0, E3M, 4, LoY) TrackedReal<DX9>(0.0, 0.0, E3M, 7, LoY); TrackedReal<6dI>(0.0, 0.0, E3M, 2, LoY) TrackedReal<3Z3>(11.393423775193504, 0.0, E3M, 5, LoY) TrackedReal<K5b>(0.0, 0.0, E3M, 8, LoY); TrackedReal<Ika>(1.647475534933883, 0.0, E3M, 3, LoY) TrackedReal<14S>(0.0, 0.0, E3M, 6, LoY) TrackedReal<H54>(10.29510675190425, 0.0, E3M, 9, LoY)]
tr(σ): TrackedReal<Amk>(32.53279579064663, 0.0, E3M, ---)
ERROR: ArgumentError: Converting an instance of ReverseDiff.TrackedReal{Float64, Float64, Nothing} to Float64 is not defined. Please use `ReverseDiff.value` instead.
Stacktrace:
  [1] convert(#unused#::Type{Float64}, t::ReverseDiff.TrackedReal{Float64, Float64, Nothing})
    @ ReverseDiff ~/.julia/packages/ReverseDiff/Zu4v6/src/tracked.jl:261
  [2] setindex!(A::Vector{Float64}, x::ReverseDiff.TrackedReal{Float64, Float64, Nothing}, i1::Int64)
    @ Base ./array.jl:969
  [3] tbnn_opt(σ::ReverseDiff.TrackedArray{Float64, Float64, 2, SizedMatrix{3, 3,
...

function dudt_univ_opt!(du, σ, p, t, gradv, model_univ, model_weights)
    # the parameters are [NN parameters, ODE parameters)
    η0, τ = @view p[end-1 : end]
    # SizeMatrix wraps an array of known size to increase efficiency. Not static.
    ∇v = SizedMatrix{3,3}([0. 0. 0. ; gradv(t) 0. 0. ; 0. 0. 0.])  # ERROR with gradv ???
    D = 0.5 .* (∇v .+ transpose(∇v))
    T1 = (η0/τ) .* D    
    T2 = (transpose(∇v) * σ) .+ (σ * ∇v)

    # Run stress/strain through a Tensor-Base Neural Network (TBNN)
    # Change tbnn to read D
    F = tbnn_opt(σ, D, model_weights, model_univ, t)
    du .= F  # TEMPORARY to SIMPLIFY the code
end

function tbnn_opt(σ, D, model_weights, model_univ, t)
    # Compute elements of the tensor basis
    I  = SA[1. 0. 0.; 0. 1. 0.; 0. 0. 1.]

    # Compute the integrity basis from scalar invariants. Traces. 
    λ = zeros(2) # must be mutale
    println("σ: ", σ)
    println("tr(σ): ", tr(σ))
    λ[1] = tr(σ)
    λ[2] = 3
    g = re(model_weights)(λ)
    println("==> g: ", g) # should be a Vector
    F = g[1] .* I .+ g[2] .* σ
    println("==> F: ", F)
    return F
end

# Neural network

model_univ = Flux.Chain(Flux.Dense(2, 8, tanh),
                       Flux.Dense(8, 2))
p_model, re = Flux.destructure(model_univ)

and the processing that leads to the error, namely the last line. The code is to complete. I can add more if necessary. It would be nice to simply use an attachment.

# Continutation training loop
adtype = Optimization.AutoZygote()
    klist = [1,1]
    trajectory = klist[1]
	loss_fn(θ) = loss_univ([θ; p_system], protocols[klist], tspans[klist], σ0, σ12_all, trajectory, model_univ, model_weights) 
	cb_fun(θ, l) = callback(θ, l)
	optf = Optimization.OptimizationFunction((x,p) -> loss_fn(x), adtype)
	optprob = Optimization.OptimizationProblem(optf, θi)
	result_univ = Optimization.solve(optprob, Optimisers.AMSGrad(), callback=cb_fun, maxiters=max_nb_iter)
5

λ = zeros(2) # must be mutale will make it Float64 when it shouldn’t be.

λ = similar(σ, (2,)) is what you want.

6

Thank you, Chris. I made the change you suggested (which I never would have guessed), and got a setIndex error. I fixed that by changing the definition of the identity matrix from a 3x3 StaticArray to a 3x3 SizedArray. Why would I have to do this? Basically, g[2]*I created a problem when I was a static array. Why would that be?

Here is the code:

function tbnn_opt(σ, D, model_weights, model_univ, t)
    # Compute elements of the tensor basis
    # I  = SA[1. 0. 0.; 0. 1. 0.; 0. 0. 1.]  # Should work   but doesn't
    I  = SizedMatrix{3,3}([1. 0. 0.; 0. 1. 0.; 0. 0. 1.])

    # Compute the integrity basis from scalar invariants. Traces. 
    # λ = zeros(2) # must be mutale
    # Important to use similar() to ensure the same type of σ, which can change
    λ = similar(σ, (2,)) # must be mutable and the same type of σ
    λ[1] = tr(σ)
    λ[2] = 3
    g = re(model_weights)(λ)
    F = g[1] .* I .+ g[2] .* σ
    return F
end
7

Is this on the new patch?

8

No. The code ran with Tsit5 and had the issue I posted.

I will now try the new patch.

9

Hi @ChrisRackauckas ,

I ran Rodas4 with the new patch and it work. However (I could post another topic), look at this src:

function tbnn_opt(σ, D, model_weights, model_univ, t)
    # Compute elements of the tensor basis
    I  = SA[1. 0. 0.; 0. 1. 0.; 0. 0. 1.]  # Should work
    # I  = SizedMatrix{3,3}([1. 0. 0.; 0. 1. 0.; 0. 0. 1.])

    λ = similar(σ, (2,)) # must be mutable and the same type of σ
    λ[1] = tr(σ)
    λ[2] = 3
    g = re(model_weights)(λ)
    F = g[1] .* I .+ g[2] .* σ
    return F
end

Using Rodas4 and a definition of the identity matrix I as a mutable array using SizedMatrix, the code works. However, if the identity matrix is defined as a static array, I get an error (first few lines of an extremely long trace):

ERROR: setindex!(::SMatrix{3, 3, Float64, 9}, value, ::Int) is not defined.
 Hint: Use `MArray` or `SizedArray` to create a mutable static array
Stacktrace:
   [1] error(s::String)
     @ Base ./error.jl:35
   [2] setindex!(a::SMatrix{3, 3, Float64, 9}, value::Float64, i::Int64)
     @ StaticArrays ~/.julia/packages/StaticArrays/VLqRb/src/indexing.jl:3
   [3] macro expansion
     @ ~/.julia/packages/StaticArrays/VLqRb/src/broadcast.jl:160 [inlined]
   [4] _broadcast!

I don’t understand why that should because the identity never changes, and the expression g[1] .* I should work regardless of the type of g[1]. When g[1] is of type ForwardDiff, there is no issue. the issue occurs when g[1] is of type ReverseDiff, see below:

==> g: ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 9}[Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}}(-0.9629427516671772,-0.003012227512831078,0.0,0.0,0.0,-0.003012227512831078,0.0,0.0,0.0,-0.003012227512831078), Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}}(0.5011892925377939,0.1570962334407297,0.0,0.0,0.0,0.1570962334407297,0.0,0.0,0.0,0.1570962334407297)]
==> g: Float32[-0.9629427, 0.50118923]
==> g: Float32[-0.9629815, 0.5026771]
==> g: Float32[-0.96296215, 0.50202286]
==> g: Float32[-0.96303296, 0.5039894]
==> g: Float32[-0.9632064, 0.5069435]
==> g: Float32[-0.96320677, 0.50695294]
==> g: ReverseDiff.TrackedReal{Float64, Float64, Nothing}[TrackedReal<6SM>(-0.9632086547532805, 0.0, KFo, ---), TrackedReal<Buo>(0.5069759643583189, 0.0, KFo, ---)]
ERROR: setindex!(::SMatrix{3, 3, Float64, 9}, value, ::Int) is not defined.
 Hint: Use `MArray` or `SizedArray` to create a mutable static array
Stacktrace:

I feel that the type structure of Julia gets so complex that as packages keep extending, it becomes harder and harder to ensure that all type combinations are taken into account, especially if people do not define their packages with the highest level of Array abstractions. But that is a discussion for another time :slight_smile: .