# 'BoundsError' on ODEProblem accelerated with ModelingToolkit.jl

**URL:** <https://discourse.julialang.org/t/boundserror-on-odeproblem-accelerated-with-modelingtoolkit-jl/86287>\
**Category:** General Usage\
**Tags:** modelingtoolkit, ordinarydiffeq\
**Created:** [August 24, 2022, 7:18pm UTC](https://discourse.julialang.org/t/boundserror-on-odeproblem-accelerated-with-modelingtoolkit-jl/86287 "2022-08-24T19:18:05Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![luismunozheinen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luismunozheinen/32/39070_2.png) [@luismunozheinen](https://discourse.julialang.org/u/luismunozheinen)\
**Post date:** [August 24, 2022, 7:18pm UTC](https://discourse.julialang.org/t/boundserror-on-odeproblem-accelerated-with-modelingtoolkit-jl/86287/1 "2022-08-24T19:18:05Z")

</div>

Hi,

I have tried to update my ODE problem to run on the lastest versions of OrdinaryDiffEq and ModelingToolkit, however I now encounter a ‘ERROR: BoundsError’ within OrdinaryDiffEq…

I have included a reduced version of my code that reproduces the error.

```julia
using OrdinaryDiffEq
using ModelingToolkit

##Constants
const Nx = 20;
const Ny = 10;
const ρ = 2.5;
const kp = 60.0;
const kc = 200.0;
const νd = 1e-9;
const ν0 = 1e-3;
const σn = 100.0;
const τ0 = 0.6;
const Dc = 1e-3;
const a = 0.04;
const b = 0.06;
##ODE Function
function f(u,p,t)
          @inbounds for i in 1:Nx, j in 1:Ny
              du[j,i,1] = @. u[j,i,2] - νd
              du[j,i,3] = @. 1.0 - u[j,i,2]*u[j,i,3]/Dc
          end
           @inbounds for i in 2:Nx-1, j in 2:Ny-1
                du[j,i,2] = @. 1/ρ*kc*(u[j,i+1,1]+u[j,i-1,1]+u[j+1,i,1]+u[j-1,i,1]-4*u[j,i,1]) - 1/ρ*kp*u[j,i,1] - 1/ρ*σn*a*asinh(u[j,i,2]/2/ν0*exp((τ0+b*log(ν0*u[j,i,3]/Dc))/a))
           end
           @inbounds for j in 2:Ny-1
               du[j,1,2] = @. 1/ρ*kc*(u[j,2,1]+u[j+1,1,1]+u[j-1,1,1]-3*u[j,1,1]) - 1/ρ*kp*u[j,1,1] - 1/ρ*σn*a*asinh(u[j,1,2]/2/ν0*exp((τ0+b*log(ν0*u[j,1,3]/Dc))/a))
               du[j,Nx,2] = @. 1/ρ*kc*(u[j,Nx-1,1]+u[j+1,Nx,1]+u[j-1,Nx,1]-3*u[j,Nx,1]) - 1/ρ*kp*u[j,Nx,1] - 1/ρ*σn*a*asinh(u[j,Nx,2]/2/ν0*exp((τ0+b*log(ν0*u[j,Nx,3]/Dc))/a))
           end
           @inbounds for i in 2:Nx-1
               du[1,i,2] = @. 1/ρ*kc*(u[1,i+1,1]+u[1,i-1,1]+u[2,i,1]-3*u[1,i,1]) - 1/ρ*kp*u[1,i,1] - 1/ρ*σn*a*asinh(u[1,i,2]/2/ν0*exp((τ0+b*log(ν0*u[1,i,3]/Dc))/a))
               du[Ny,i,2] = @. 1/ρ*kc*(u[Ny,i+1,1]+u[Ny,i-1,1]+u[Ny-1,i,1]-3*u[Ny,i,1]) - 1/ρ*kp*u[Ny,i,1] - 1/ρ*σn*a*asinh(u[Ny,i,2]/2/ν0*exp((τ0+b*log(ν0*u[Ny,i,3]/Dc))/a))
           end
           @inbounds begin
                du[1,1,2] = @. 1/ρ*kc*(u[1,2,1]+u[2,1,1]-2*u[1,1,1]) - 1/ρ*kp*u[1,1,1] - 1/ρ*σn*a*asinh(u[1,1,2]/2/ν0*exp((τ0+b*log(ν0*u[1,1,3]/Dc))/a))
                du[1,Nx,2] = @. 1/ρ*kc*(u[1,Nx-1,1]+u[2,Nx,1]-2*u[1,Nx,1]) - 1/ρ*kp*u[1,Nx,1] - 1/ρ*σn*a*asinh(u[1,Nx,2]/2/ν0*exp((τ0+b*log(ν0*u[1,Nx,3]/Dc))/a))
                du[Ny,1,2] = @. 1/ρ*kc*(u[Ny,2,1]+u[Ny-1,1,1]-2*u[Ny,1,1]) - 1/ρ*kp*u[Ny,1,1] - 1/ρ*σn*a*asinh(u[Ny,1,2]/2/ν0*exp((τ0+b*log(ν0*u[Ny,1,3]/Dc))/a))
                du[Ny,Nx,2] = @. 1/ρ*kc*(u[Ny,Nx-1,1]+u[Ny-1,Nx,1]-2*u[Ny,Nx,1]) - 1/ρ*kp*u[Ny,Nx,1] - 1/ρ*σn*a*asinh(u[Ny,Nx,2]/2/ν0*exp((τ0+b*log(ν0*u[Ny,Nx,3]/Dc))/a))
            end
end
##ODE Speed-up
@variables u[1:Ny, 1:Nx, 1:3]
du = similar(u);
f(collect(u), nothing, 0.0);
du = simplify.(du);
fastf = eval(ModelingToolkit.build_function(du,u,parallel=Symbolics.MultithreadedForm())[2]);
jac = Symbolics.sparsejacobian(vec(du), vec(collect(u)));
fjac = eval(Symbolics.build_function(jac,u,parallel=Symbolics.MultithreadedForm())[2]);
##Initial Conditions
u0 = Array{Float64, 3}(undef,Ny,Nx,3);
u0[:,:,1] = - σn * a /kp* asinh(νd/(2*ν0)*exp((τ0+b*log(ν0/νd))/a)).*ones(Ny,Nx) -1e-6.*rand(Ny,Nx);
u0[:,:,2] = νd.*ones(Ny,Nx);
u0[:,:,3] = Dc/νd.*ones(Ny,Nx);
##solver
tspan = (0.0,1e1)
prob_jac = ODEProblem(ODEFunction((du, u, p, t) -> fastf(du, u),
                                   jac = (du, u, p, t) -> fjac(du, u),
                                   jac_prototype = similar(jac, Float64)),
                                   u0, tspan);
#solver algorithm which produces error
solve(prob_jac,Rodas5(),reltol = 1e-8,abstol = 1e-8,maxiters = 1e10,save_everystep = false,dtmin = 1e-20,dt = 1e-10,force_dtmin =true);
#other solver algorithms don't produce errors
solve(prob_jac,Tsit5(),reltol = 1e-8,abstol = 1e-8,maxiters = 1e10,save_everystep = false,dtmin = 1e-20,dt = 1e-10,force_dtmin =true);

```

The error message when using the ‘Rodas5’ solver algorithm is

```julia
ERROR: BoundsError
Stacktrace:
  [1] _copyto_impl!(dest::Vector{Float64}, doffs::Int64, src::Vector{Float64}, soffs::Int64, n::Int64)
    @ Base ./array.jl:329
  [2] copyto!
    @ ./array.jl:322 [inlined]
  [3] copyto!(dest::Vector{Float64}, src::Vector{Float64})
    @ Base ./array.jl:346
  [4] solve(cache::LinearSolve.LinearCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, LinearSolve.KLUFactorization, KLU.KLUFactorization{Float64, Int64}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64}, alg::LinearSolve.KLUFactorization; kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:reltol,), Tuple{Float64}}})
    @ LinearSolve ~/.julia/packages/LinearSolve/FDu3j/src/factorization.jl:296
  [5] #solve#5
    @ ~/.julia/packages/LinearSolve/FDu3j/src/common.jl:143 [inlined]
  [6] #dolinsolve#3
    @ ~/.julia/packages/OrdinaryDiffEq/SmImO/src/misc_utils.jl:105 [inlined]
  [7] perform_step!(integrator::OrdinaryDiffEq.ODEIntegrator{Rodas5{12, true, LinearSolve.KLUFactorization, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, true, Array{Float64, 3}, Nothing, Float64, SciMLBase.NullParameters, Float64, Float64, Float64, Float64, Vector{Array{Float64, 3}}, ODESolution{Float64, 4, Vector{Array{Float64, 3}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Array{Float64, 3}}}, ODEProblem{Array{Float64, 3}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, var"#17#19",LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, var"#18#20", Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Rodas5{12, true, LinearSolve.KLUFactorization, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, var"#17#19", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, var"#18#20", Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Array{Float64, 3}}, Vector{Float64}, Vector{Vector{Array{Float64, 3}}}, OrdinaryDiffEq.Rosenbrock5Cache{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, OrdinaryDiffEq.Rodas5Tableau{Float64, Float64}, SciMLBase.TimeGradientWrapper{ODEFunction{true, var"#17#19", LinearAlgebra.UniformScaling{Bool},Nothing, Nothing, var"#18#20", Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Array{Float64, 3}, SciMLBase.NullParameters}, SciMLBase.UJacobianWrapper{ODEFunction{true, var"#17#19", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, var"#18#20", Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, SciMLBase.NullParameters}, LinearSolve.LinearCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, LinearSolve.KLUFactorization, KLU.KLUFactorization{Float64, Int64}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64}, Nothing, Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}, 3}, Float64, Rodas5{12, true, LinearSolve.KLUFactorization, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}}, DiffEqBase.DEStats}, ODEFunction{true, var"#17#19", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, var"#18#20", Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, OrdinaryDiffEq.Rosenbrock5Cache{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, OrdinaryDiffEq.Rodas5Tableau{Float64, Float64}, SciMLBase.TimeGradientWrapper{ODEFunction{true, var"#17#19", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, var"#18#20", Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Array{Float64, 3}, SciMLBase.NullParameters}, SciMLBase.UJacobianWrapper{ODEFunction{true, var"#17#19", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, var"#18#20", Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, SciMLBase.NullParameters}, LinearSolve.LinearCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, LinearSolve.KLUFactorization, KLU.KLUFactorization{Float64, Int64},LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64}, Nothing, Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}, 3}, Float64, Rodas5{12, true, LinearSolve.KLUFactorization, 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, Float64, Tuple{}, Tuple{}, Tuple{}}, Array{Float64, 3}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.Rosenbrock5Cache{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}, SparseArrays.SparseMatrixCSC{Float64,Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, OrdinaryDiffEq.Rodas5Tableau{Float64, Float64}, SciMLBase.TimeGradientWrapper{ODEFunction{true, var"#17#19", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, var"#18#20", Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Array{Float64, 3}, SciMLBase.NullParameters}, SciMLBase.UJacobianWrapper{ODEFunction{true, var"#17#19", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, var"#18#20", Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, SciMLBase.NullParameters}, LinearSolve.LinearCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, LinearSolve.KLUFactorization, KLU.KLUFactorization{Float64, Int64}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64}, Nothing, Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}, 3}, Float64, Rodas5{12, true, LinearSolve.KLUFactorization, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, repeat_step::Bool)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/SmImO/src/perform_step/rosenbrock_perform_step.jl:1800
  [8] perform_step!
    @ ~/.julia/packages/OrdinaryDiffEq/SmImO/src/perform_step/rosenbrock_perform_step.jl:1742 [inlined]
  [9] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{Rodas5{12, true, LinearSolve.KLUFactorization, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, true, Array{Float64, 3}, Nothing, Float64, SciMLBase.NullParameters, Float64, Float64, Float64, Float64, Vector{Array{Float64, 3}}, ODESolution{Float64, 4, Vector{Array{Float64, 3}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Array{Float64, 3}}}, ODEProblem{Array{Float64, 3}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, var"#17#19", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, var"#18#20", Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Rodas5{12, true, LinearSolve.KLUFactorization,typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, var"#17#19", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, var"#18#20", Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Array{Float64, 3}}, Vector{Float64}, Vector{Vector{Array{Float64, 3}}}, OrdinaryDiffEq.Rosenbrock5Cache{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, OrdinaryDiffEq.Rodas5Tableau{Float64, Float64}, SciMLBase.TimeGradientWrapper{ODEFunction{true, var"#17#19", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, var"#18#20", Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Array{Float64, 3}, SciMLBase.NullParameters}, SciMLBase.UJacobianWrapper{ODEFunction{true, var"#17#19", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, var"#18#20", Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, SciMLBase.NullParameters}, LinearSolve.LinearCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, LinearSolve.KLUFactorization, KLU.KLUFactorization{Float64, Int64}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64}, Nothing, Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}, 3}, Float64, Rodas5{12, true, LinearSolve.KLUFactorization, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}}, DiffEqBase.DEStats}, ODEFunction{true, var"#17#19", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, var"#18#20", Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, OrdinaryDiffEq.Rosenbrock5Cache{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, OrdinaryDiffEq.Rodas5Tableau{Float64, Float64}, SciMLBase.TimeGradientWrapper{ODEFunction{true, var"#17#19", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, var"#18#20", Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Array{Float64, 3}, SciMLBase.NullParameters}, SciMLBase.UJacobianWrapper{ODEFunction{true, var"#17#19", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, var"#18#20", Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, SciMLBase.NullParameters}, LinearSolve.LinearCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, LinearSolve.KLUFactorization, KLU.KLUFactorization{Float64, Int64}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64}, Nothing, Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}, 3}, Float64, Rodas5{12, true, LinearSolve.KLUFactorization, 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, Float64, Tuple{}, Tuple{}, Tuple{}}, Array{Float64, 3}, Float64, Nothing, OrdinaryDiffEq.DefaultInit})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/SmImO/src/solve.jl:475
 [10] #__solve#562
    @ ~/.julia/packages/OrdinaryDiffEq/SmImO/src/solve.jl:5 [inlined]
 [11] #solve_call#24
    @ ~/.julia/packages/DiffEqBase/x9wgy/src/solve.jl:451 [inlined]
 [12] solve_up(prob::ODEProblem{Array{Float64, 3}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, var"#17#19", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, var"#18#20", Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Array{Float64,3}, p::SciMLBase.NullParameters, args::Rodas5{0, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}; kwargs::Base.Pairs{Symbol, Real, NTuple{7, Symbol}, NamedTuple{(:reltol, :abstol, :maxiters, :save_everystep, :dtmin, :dt, :force_dtmin), Tuple{Float64, Float64, Float64, Bool, Float64, Float64, Bool}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/x9wgy/src/solve.jl:814
 [13] #solve#29
    @ ~/.julia/packages/DiffEqBase/x9wgy/src/solve.jl:784 [inlined]
 [14] top-level scope
    @ REPL[54263]:1

```

My versions are

```julia
(@v1.8) pkg> st
Status `~/.julia/environments/v1.8/Project.toml`
  [6e4b80f9] BenchmarkTools v1.3.1
  [336ed68f] CSV v0.10.4
  [35d6a980] ColorSchemes v3.19.0
  [a93c6f00] DataFrames v1.3.4
  [0c46a032] DifferentialEquations v7.2.0
  [5789e2e9] FileIO v1.15.0
  [033835bb] JLD2 v0.4.22
  [961ee093] ModelingToolkit v8.20.0
  [1dea7af3] OrdinaryDiffEq v6.23.0
  [91a5bcdd] Plots v1.31.7
  [10745b16] Statistics

```

I appreciate any help to fix this issue. The objective here is to _start_ solving the ode error-free with the ‘Rodas5’-solver, and not to completely solve the problem for the full time-span (I’ve simplified the actual physical problem/setup expected to be solved).

Regards,  
Luis

---

<div class="post-metadata">

**Author:** ![SteffenPL](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/steffenpl/32/206270_2.png) [@SteffenPL](https://discourse.julialang.org/u/SteffenPL)\
**Post date:** [August 24, 2022, 7:29pm UTC](https://discourse.julialang.org/t/boundserror-on-odeproblem-accelerated-with-modelingtoolkit-jl/86287/2 "2022-08-24T19:29:54Z")

</div>

Could it be that the right-hand side expects input of dimension `Ny x Nx x 3` whereas the Jacobian is for a vector with length `Nx * Ny * 3`.

Does it work if you add this after the line `du = simplify.(du)`

```julia
u = vec(u)
du = vec(du)

```

and accordingly `u0 = vec(u0)`…  
That way, both functions (right-hand side and jacobian) would use the vectorized version.

The alternative could be to add `reshape` commands around the call of `fjac`.

---

<div class="post-metadata">

**Author:** ![luismunozheinen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luismunozheinen/32/39070_2.png) [@luismunozheinen](https://discourse.julialang.org/u/luismunozheinen)\
**Post date:** [August 25, 2022, 4:16pm UTC](https://discourse.julialang.org/t/boundserror-on-odeproblem-accelerated-with-modelingtoolkit-jl/86287/3 "2022-08-25T16:16:57Z")

</div>

Thank you @SteffenPL for your reply. I provide few more precisions: My optimization is ‘copied’ from the example in the documentation of an earlier version of ModelingToolkit.jl, see this [link](https://mtk.sciml.ai/v5.8/tutorials/auto_parallel/), and used to work with older versions of Julia, OrdinaryDiffEq and ModelingToolkit. Still, I have tried your changes but the vectorized version resulted in exactly the same error…

Using the current ModelingToolkit documention ([link](https://mtk.sciml.ai/stable/mtkitize_tutorials/sparse_jacobians/)), I am able to run the vectorized form via `sys = modelingtoolkitize(prob)` and `sparseprob = ODEProblem(sys,u0,(0.,1e1),jac=true,sparse=true)`. However, in previous benchmarks, I determined the `Ny x Nx x 3` form to be significantly faster than the vectorized `Nx * Ny * 3` for my specific problem (also the new multithreaded version lead to a segmentation fault…). Hence I would like the `Ny x Nx x 3` form to work…

I wonder why the example only fails for certain solvers?

---

<div class="post-metadata">

**Author:** ![SteffenPL](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/steffenpl/32/206270_2.png) [@SteffenPL](https://discourse.julialang.org/u/SteffenPL)\
**Post date:** [August 25, 2022, 4:23pm UTC](https://discourse.julialang.org/t/boundserror-on-odeproblem-accelerated-with-modelingtoolkit-jl/86287/4 "2022-08-25T16:23:53Z")

</div>

A quick answer to the last question is that some solvers do not use the Jacobian. In particular most explicit solvers like `Tsit5` do not need this information so they will never call it. Most likely the provided Jacobian function has some dimension error somewhere. (Which is of course always tricky to find 😉 )

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 26, 2022, 1:45am UTC](https://discourse.julialang.org/t/boundserror-on-odeproblem-accelerated-with-modelingtoolkit-jl/86287/5 "2022-08-26T01:45:04Z")

</div>

What if you just let ModelingToolkit do its thing instead of manually using the Symbolics to generate the sparse Jacobian and such? If you just follow the ModelingToolkit documented form, i.e.

```julia
prob_jac = ODEProblem(sys, u0, tspan, jac=true,sparse=true)

```

is it fine? That would narrow it down.

---

<div class="post-metadata">

**Author:** ![luismunozheinen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luismunozheinen/32/39070_2.png) [@luismunozheinen](https://discourse.julialang.org/u/luismunozheinen)\
**Post date:** [August 26, 2022, 5:31pm UTC](https://discourse.julialang.org/t/boundserror-on-odeproblem-accelerated-with-modelingtoolkit-jl/86287/6 "2022-08-26T17:31:09Z")

</div>

Hi @ChrisRackauckas, thank you for the reply. Using `prob_jac = ODEProblem(sys, u0, tspan, jac=true,sparse=true)` works, however using `prob_jac_Mthread = ODEProblem(sys, u0, tspan, jac=true,sparse=true,parallel = ModelingToolkit.MultithreadedForm())` only worked for the model with (`Nx = 20;Ny = 10`); solving (`Nx=60,Ny=20`) caused different `TaskFailedException` errors such as:

```julia
julia> @time solve(prob_jacMthread,Rodas5(),reltol = 1e-8,abstol = 1e-8,maxiters = 1e10,save_everystep = false,dtmin = 1e-20,dt = 1e-10,force_dtmin =true);
ERROR: TaskFailedException
Stacktrace:
  [1] wait
    @ ./task.jl:345 [inlined]
  [2] fetch
    @ ./task.jl:360 [inlined]
  [3] map
    @ ./tuple.jl:224 [inlined]
  [4] macro expansion
    @ ~/.julia/packages/Symbolics/J8IHJ/src/build_function.jl:351 [inlined]
  [5] macro expansion
    @ ~/.julia/packages/SymbolicUtils/qulQp/src/code.jl:351 [inlined]
  [6] macro expansion
    @ ~/.julia/packages/RuntimeGeneratedFunctions/KrkGo/src/RuntimeGeneratedFunctions.jl:129 [inlined]
  [7] macro expansion
    @ ./none:0 [inlined]
  [8] generated_callfunc
    @ ./none:0 [inlined]
  [9] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag",ModelingToolkit.var"#_RGF_ModTag", (0x8e58222a, 0x6dccd51c, 0xaf0903ef, 0x37e53944, 0xbe3d02b3)})(::SparseArrays.SparseMatrixCSC{Float64, Int64}, ::Vector{Float64}, ::Nothing, ::Float64)
    @ RuntimeGeneratedFunctions ~/.julia/packages/RuntimeGeneratedFunctions/KrkGo/src/RuntimeGeneratedFunctions.jl:117
 [10] _jac
    @ ~/.julia/packages/ModelingToolkit/5qzZn/src/systems/diffeqs/abstractodesystem.jl:295 [inlined]
 [11] calc_J!
    @ ~/.julia/packages/OrdinaryDiffEq/ccwlJ/src/derivative_utils.jl:148 [inlined]
 [12] calc_W!
    @ ~/.julia/packages/OrdinaryDiffEq/ccwlJ/src/derivative_utils.jl:652 [inlined]
 [13] calc_W!
    @ ~/.julia/packages/OrdinaryDiffEq/ccwlJ/src/derivative_utils.jl:599 [inlined]
 [14] ... 

```

Besides my inability to run the multithreaded version, I have two more concerns when following the current ModelingToolkit documentation:

- I assume that `modelingtoolkitize` automatically vectorizes the problem and changes the type of variable `u` from `Array{Float64, 3}` to `Vector{Float64}`; I previously ran benchmarks in which `Array{Float64, 3}` performed better for my problem and would hence prefer that option (or at least be able to compare both options)
- I often resubmit the same problem, just changing `u0`. Running `prob_jac = ODEProblem(sys, u0, tspan, jac=true,sparse=true)` for larger problems is more time-consuming than computing a sparse `jac` once and simply providing it to `ODEFunction`

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 26, 2022, 5:44pm UTC](https://discourse.julialang.org/t/boundserror-on-odeproblem-accelerated-with-modelingtoolkit-jl/86287/7 "2022-08-26T17:44:17Z")

</div>

Can you add `checkbounds=true` to the ODEProblem? It might be indexing out of bounds in the multithreaded sharded form? That’s my guess.

> [@luismunozheinen](#):
>
> I assume that `modelingtoolkitize` automatically vectorizes the problem and changes the type of variable `u` from `Array{Float64, 3}` to `Vector{Float64}`; I previously ran benchmarks in which `Array{Float64, 3}` performed better for my problem and would hence prefer that option (or at least be able to compare both options)

That’s because of the SIMD you’d get from the Cartesian iteration. The MTK generated code will lose that and will be slower, changing to `Array{Float64, 3}` will not help that in the current code gen since it would be more about the SLP vectorization pass. For this case, it would be most optimal to build the sparse Jacobian code and add that to the original `f`, though I wouldn’t be surprised if the difference in the `f`’s don’t make a major difference in the final ODE time since it’s probably dominated by Jacobian time.

> [@luismunozheinen](#):
>
> I often resubmit the same problem, just changing `u0`. Running `prob_jac = ODEProblem(sys, u0, tspan, jac=true,sparse=true)` for larger problems is more time-consuming than computing a sparse `jac` once and simply providing it to `ODEFunction`

That makes sense. I’m mostly pushing you down this route to reduce the kinds of possible bugs to hone in on what’s actually going on. For your real case, that split approach is a good idea.

It sounds like the issue is the multithreaded sharded form probably “completes” the indexing, so if it splits to 4 parts with 4 threads on 15 ODEs, it indexes 1:16 with even parts in all 4, and errors out with a bounds error (because it should’ve cut the last one early). Right now that’s just a guess though given your last post.

---

<div class="post-metadata">

**Author:** ![luismunozheinen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luismunozheinen/32/39070_2.png) [@luismunozheinen](https://discourse.julialang.org/u/luismunozheinen)\
**Post date:** [August 29, 2022, 5:58pm UTC](https://discourse.julialang.org/t/boundserror-on-odeproblem-accelerated-with-modelingtoolkit-jl/86287/8 "2022-08-29T17:58:16Z")

</div>

Thanks @ChrisRackauckas for providing the explanations on how different implementations impact the ode solving performance. I’ve tested `checkbounds=true` in `ODEProblem(sys,nothing,(0.,1e1),jac=true,sparse=true,parallel = ModelingToolkit.MultithreadedForm(),checkbounds=true)`, however it doesn’t solve the issue, and both setups (with or without `checkbounds`) will error with `ERROR: LoadError: TaskFailedException` for larger model sizes (i.e. `Nx>50;Ny=20` for a single-node job submission with ncpus=8 and memory=96gb on a cluster computer). To determine if the issue is specific to my ODE system, I have run a multi-threaded version of the Brusselator PDE from the ModelingToolkit documentation; it works for `N=32`, but fails with segfault for `N=40` (with or without `checkbounds`).

Regarding the original problem (using the split approach), I have determined that same setup without multi-threading also results in the original error message `ERROR: BoundsError` detailed in the first post. Hence, the issue should be to the novel generation of `fjac` and `jac` (in the current version of ModelingTookit) rather than the multi-threading part.

> [@ChrisRackauckas](#):
>
> For this case, it would be most optimal to build the sparse Jacobian code and add that to the original `f`, though I wouldn’t be surprised if the difference in the `f`’s don’t make a major difference in the final ODE time since it’s probably dominated by Jacobian time.

I run my simulations for an extended amount of time, so that small gains will make certainly make a difference in longer simulations.

> [@ChrisRackauckas](#):
>
> It sounds like the issue is the multithreaded sharded form probably “completes” the indexing, so if it splits to 4 parts with 4 threads on 15 ODEs, it indexes 1:16 with even parts in all 4, and errors out with a bounds error (because it should’ve cut the last one early).

The number of ODEs should be the system size times the degree of freedom per node, i.e. `Nx*Ny*3=20*60*3` which then shouldn’t be an issue using `export JULIA_NUM_THREADS` equal to 2, 4 or 8…

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 24, 2022, 6:38am UTC](https://discourse.julialang.org/t/boundserror-on-odeproblem-accelerated-with-modelingtoolkit-jl/86287/9 "2022-12-24T06:38:40Z")

</div>

These PDE cases now have a lot more benchmarks, see [Brusselator Work-Precision Diagrams · The SciML Benchmarks](https://docs.sciml.ai/SciMLBenchmarksOutput/stable/StiffODE/Bruss/) . It would be nice to capture this one as well.

---

<div class="post-metadata">

**Author:** ![luismunozheinen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luismunozheinen/32/39070_2.png) [@luismunozheinen](https://discourse.julialang.org/u/luismunozheinen)\
**Post date:** [January 14, 2023, 5:28pm UTC](https://discourse.julialang.org/t/boundserror-on-odeproblem-accelerated-with-modelingtoolkit-jl/86287/10 "2023-01-14T17:28:57Z")

</div>

Hi @ChrisRackauckas, thank you for the update. I can confirm that I am now able to run my ODE system on the latest versions of OrdinaryDiffEq and ModelingToolkit, using an implementation similar to the [Brusselator Work-Precision Diagrams](https://docs.sciml.ai/SciMLBenchmarksOutput/stable/StiffODE/Bruss/). Despite the ‘loss of SIMD’, I observe a speed-up with respect to my previous implementation (as presented in the first post).

Regarding the work-precision diagrams, I first provide information on the type of PDE system:

- Each element at `j,i` is characterised by its position `u[j,i,1]` (-1e1,1e1), velocity `u[j,i,2]` (1e-15,1e2) and state `u[j,i,3]` (1e0,1e9). The approximate ranges of the variables are indicated in parentheses.
- The system specifies the linear interactions between elements, as well as a non-linear resistance term (which in turn describes a velocity weakening effect)
- For suitable parameter ranges, the model describes a chaotic dynamical system and ‘self-organizes’ into critical states which result in a fractal distribution of (stick and) slip motion over time. (one larger, ‘system-wide’ slip is used for the work-precision benchmark).

Regarding the work-precision benchmark, I have ‘tuned’ it to my needs, to establish which solver is most suitable for my studies.

- To avoid excessive benchmark times, I focus on a (random) single large slip observed during the chaotic cycle simulations. (If needed, one might also be able to approximate the initial conditions by spatial random fields, so that this benchmark can be reproduced without the full cycle simulations).
- I specify an error metric based on the maximal difference in the cumulative slip rate, sampled at 1-sec time intervals, `error = maximum(abs((sum(u[:,:,2])-sum(u_ref[:,:,2]))/sum(u_ref[:,:,2])))`, where `u_ref` corresponds to the reference solution or `appxsol` in `WorkPrecisionSet` (here using `CVODE_BDF MTK KLU`, ` abstol = 1e-20`, `reltol = 1e-15`). I then compute the average error from 4 simulation runs.
- solver tolerances based on any combination of `abstols = 1.0 ./ 10.0 .^ (6:2:20)` and `reltols = 1.0 ./ 10.0 .^ (4:20)`

So this is what the benchmark simulation looks like (total spatial slip and velocity in time respectively), for a system using` j = 1:18` and `i = 1:66` elements:

 ![Fig_Simulation](https://global.discourse-cdn.com/julialang/original/3X/6/2/62799d5244d2dd0ea1dd4a0fe92841dbb0f1a1ea.png)

And here is the Work-Precision Diagram, for the error metric defined previously and a maximum elapsed time of 4h. (I was relying on `Rodas5()` before this excercise…)

 ![Fig_WorkPrecision](https://global.discourse-cdn.com/julialang/original/3X/2/8/285dc5ccd2338b8f1db8842694d30d08db2fd1f8.png)

Few additional observations:

- KLU factorization enables more solvers to complete within the limited time, though no significant speed-ups have been observed between e.g. `Rodas5()` and `Rodas5(linsolve=KLUFactorization())`
- no setup using `GMRES` succeeded with the simulation
- no gain from `ModelingToolkit.MultithreadedForm()`
- no completion for `FBDF`,`ESERK5`,`ROCK4`,`GRK4A`,`Kvaerno4`,`Cash4`,`Hairer4`  
(initially tested for abstol = 1.0 ./ 10.0 .^ (12:2:14), reltol = 1.0 ./10.0 .^ (4:2:10)) and `CVODE_BDF{:Newton, :GMRES)` (initially tested for abstol 1.0 ./ 10.0 .^ (10:2:18), reltol 1.0 ./ 10.0 .^(10:2:18))
- not enough precision from `QNDF`

Computer information:

```julia
Julia Version 1.8.4
Commit 00177ebc4fc (2022-12-23 21:32 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 256 × AMD EPYC 7742 64-Core Processor
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, znver2)
  Threads: 1 on 256 virtual cores
Environment:
  JULIA_NUM_THREADS = 1

```

Current versions:

```julia
(@v1.8) pkg> st
Status `~/.julia/environments/v1.8/Project.toml`
 [2169fc97] AlgebraicMultigrid v0.5.1
 [fbb218c0] BSON v0.3.6
 [6e4b80f9] BenchmarkTools v1.3.2
⌃ [336ed68f] CSV v0.10.8
 [35d6a980] ColorSchemes v3.20.0
 [a93c6f00] DataFrames v1.4.4
 [f3b72e0c] DiffEqDevTools v2.33.0
 [5789e2e9] FileIO v1.16.0
 [40713840] IncompleteLU v0.2.1
 [033835bb] JLD2 v0.4.29
 [682c06a0] JSON v0.21.3
 [7f56f5a3] LSODA v0.7.2
⌃ [7ed4a6bd] LinearSolve v1.32.1
⌃ [961ee093] ModelingToolkit v8.38.0
 [09606e27] ODEInterfaceDiffEq v3.12.0
⌃ [1dea7af3] OrdinaryDiffEq v6.36.2
 [65888b18] ParameterizedFunctions v5.15.0
⌃ [91a5bcdd] Plots v1.38.0
⌃ [438e738f] PyCall v1.94.1
 [d330b81b] PyPlot v2.11.0
 [c3572dad] Sundials v4.11.4
 [0c5d862f] Symbolics v4.14.0
 [b8865327] UnicodePlots v3.3.2

```

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [January 14, 2023, 5:50pm UTC](https://discourse.julialang.org/t/boundserror-on-odeproblem-accelerated-with-modelingtoolkit-jl/86287/11 "2023-01-14T17:50:18Z")

</div>

Just for the record, the pileup where almost all of your solvers are showing 1e-5 (ish) accuracy strongly suggests that your reference solution isn’t accurate enough for the precision you are testing at.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 16, 2023, 3:46pm UTC](https://discourse.julialang.org/t/boundserror-on-odeproblem-accelerated-with-modelingtoolkit-jl/86287/12 "2023-01-16T15:46:57Z")

</div>

> [@luismunozheinen](#):
>
> Hi @ChrisRackauckas, thank you for the update. I can confirm that I am now able to run my ODE system on the latest versions of OrdinaryDiffEq and ModelingToolkit, using an implementation similar to the [Brusselator Work-Precision Diagrams](https://docs.sciml.ai/SciMLBenchmarksOutput/stable/StiffODE/Bruss/). Despite the ‘loss of SIMD’, I observe a speed-up with respect to my previous implementation (as presented in the first post).

This looks pretty great! Could you do a PR to add this to the benchmarking platform so we can play with it some more?

> [@Oscar\_Smith](#):
>
> Just for the record, the pileup where almost all of your solvers are showing 1e-5 (ish) accuracy strongly suggests that your reference solution isn’t accurate enough for the precision you are testing at.

This. Whatever was chosen as the reference solution is definitely not the most accurate solution.

> [@luismunozheinen](#):
>
> no completion for `FBDF`,`ESERK5`,`ROCK4`,`GRK4A`,`Kvaerno4`,`Cash4`,`Hairer4`  
> (initially tested for abstol = 1.0 ./ 10.0 .^ (12:2:14), reltol = 1.0 ./10.0 .^ (4:2:10)) and `CVODE_BDF{:Newton, :GMRES)` (initially tested for abstol 1.0 ./ 10.0 .^ (10:2:18), reltol 1.0 ./ 10.0 .^(10:2:18))

You cannot get a relative tolerance lower than 1e-16 with 64-bit floats. Anything lower is just going to increase the numerical error.

> [@luismunozheinen](#):
>
> KLU factorization enables more solvers to complete within the limited time, though no significant speed-ups have been observed between e.g. `Rodas5()` and `Rodas5(linsolve=KLUFactorization())`

It uses KLU by default, so it shouldn’t be close to the same, it should be the same 😅

> [@luismunozheinen](#):
>
> no setup using `GMRES` succeeded with the simulation

You likely need to use a preconditioner. Unpreconditioned GMRES is never that great.

> [@luismunozheinen](#):
>
> no gain from `ModelingToolkit.MultithreadedForm()`

That’s interesting. This should be a good case for this. Probably a good test case for optimizing then.

> [@luismunozheinen](#):
>
> - no completion for `FBDF`,`ESERK5`,`ROCK4`,`GRK4A`,`Kvaerno4`,`Cash4`,`Hairer4`  
> (initially tested for abstol = 1.0 ./ 10.0 .^ (12:2:14), reltol = 1.0 ./10.0 .^ (4:2:10)) and `CVODE_BDF{:Newton, :GMRES)` (initially tested for abstol 1.0 ./ 10.0 .^ (10:2:18), reltol 1.0 ./ 10.0 .^(10:2:18))
> - not enough precision from `QNDF`

Only FBDF is unexpected there.

---

<div class="post-metadata">

**Author:** ![luismunozheinen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luismunozheinen/32/39070_2.png) [@luismunozheinen](https://discourse.julialang.org/u/luismunozheinen)\
**Post date:** [January 26, 2023, 7:13pm UTC](https://discourse.julialang.org/t/boundserror-on-odeproblem-accelerated-with-modelingtoolkit-jl/86287/13 "2023-01-26T19:13:35Z")

</div>

Hi, I have created a file describing the model and simulation, similar to what I used for my work-precision benchmark. I’ve added it to this [PR](https://github.com/SciML/SciMLBenchmarks.jl/pull/571) (I have little experience with PRs, so I hope I did it correctly). A suitable Work-precision benchmark then needs to be defined and added to the file…

---

<div class="post-metadata">

**Author:** ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)\
**Post date:** [January 29, 2023, 5:26pm UTC](https://discourse.julialang.org/t/boundserror-on-odeproblem-accelerated-with-modelingtoolkit-jl/86287/14 "2023-01-29T17:26:41Z")

</div>

Hey @luismunozheinen , do you have a link to the mathematical description of the system to ease me also benchmarking MethodOfLines against this?

---

<div class="post-metadata">

**Author:** ![luismunozheinen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luismunozheinen/32/39070_2.png) [@luismunozheinen](https://discourse.julialang.org/u/luismunozheinen)\
**Post date:** [January 31, 2023, 9:27pm UTC](https://discourse.julialang.org/t/boundserror-on-odeproblem-accelerated-with-modelingtoolkit-jl/86287/15 "2023-01-31T21:27:48Z")

</div>

Hi @xtalax, a description of the model (except few modifications) can be found, for example, [here](https://royalsocietypublishing.org/doi/10.1098/rsta.2017.0391). The latter should facilitate the understanding of my implementation.

---

<div class="post-metadata">

**Author:** ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)\
**Post date:** [February 1, 2023, 3:54pm UTC](https://discourse.julialang.org/t/boundserror-on-odeproblem-accelerated-with-modelingtoolkit-jl/86287/16 "2023-02-01T15:54:48Z")

</div>

Thanks this is useful, I don’t recognize these boundary conditions though, what are they - neumann/diriclet/robin?

Edit: Ah, looks like neumann

---

<div class="post-metadata">

**Author:** ![luismunozheinen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luismunozheinen/32/39070_2.png) [@luismunozheinen](https://discourse.julialang.org/u/luismunozheinen)\
**Post date:** [February 1, 2023, 5:37pm UTC](https://discourse.julialang.org/t/boundserror-on-odeproblem-accelerated-with-modelingtoolkit-jl/86287/17 "2023-02-01T17:37:36Z")

</div>

Actually there are no constraints applied to the boundary blocks and they are free to move along their one spatial degree of freedom (as shown in figure 1 of the linked publication). Fyi, the ‘boundary conditions’ of the actual physical problem are replicated by applying a heterogeneous parameter distribution.