Incompatibility with units

[lamorton]

Here's what I tried to do:

using ModelingToolkit,Unitful,NonlinearSolve
vars = @variables τ_lawson 
pars = @parameters τ_ratio τ_slow
eqs = [
    τ_ratio ~ τ_slow/τ_lawson,
]

ns = NonlinearSystem(eqs, vars, pars)
ps =[τ_ratio => 1.0,
     τ_slow=>1.0u"s"]
guess = [τ_lawson=>0.1u"s"]
prob = NonlinearProblem(ns,guess,ps)
sol = solve(prob,NewtonRaphson())

and then this happened:

julia> sol = solve(prob,NewtonRaphson())
ERROR: DimensionError: s and 10.0 are not dimensionally compatible.
Stacktrace:
  [1] convert(#unused#::Type{Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}, x::Float64)
    @ Unitful ~/.julia/packages/Unitful/yI5QN/src/conversion.jl:112
  [2] setindex!(A::Vector{Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}, x::Float64, i1::Int64)
    @ Base /Applications/Julia-1.6.app/Contents/Resources/julia/share/julia/base/array.jl:839
  [3] macro expansion
    @ ~/.julia/packages/SymbolicUtils/SXYR9/src/code.jl:325 [inlined]
  [4] macro expansion
    @ ~/.julia/packages/Symbolics/H8dtg/src/build_function.jl:359 [inlined]
  [5] macro expansion
    @ ~/.julia/packages/SymbolicUtils/SXYR9/src/code.jl:283 [inlined]
  [6] macro expansion
    @ ~/.julia/packages/RuntimeGeneratedFunctions/3SZ1T/src/RuntimeGeneratedFunctions.jl:124 [inlined]
  [7] macro expansion
    @ ./none:0 [inlined]
  [8] generated_callfunc
    @ ./none:0 [inlined]
  [9] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x73c43a35, 0x06c47ba7, 0xb1d644a0, 0xa2947d04, 0xc57c3159)})(::Vector{Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}, ::Vector{Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}, ::Vector{Quantity{Float64, D, U} where {D, U}})
    @ RuntimeGeneratedFunctions ~/.julia/packages/RuntimeGeneratedFunctions/3SZ1T/src/RuntimeGeneratedFunctions.jl:112
 [10] f
    @ ~/.julia/packages/ModelingToolkit/VigaA/src/systems/nonlinear/nonlinearsystem.jl:156 [inlined]
 [11] NonlinearFunction
    @ ~/.julia/packages/SciMLBase/1aTqd/src/scimlfunctions.jl:342 [inlined]
 [12] init(::NonlinearProblem{Vector{Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}, true, Vector{Quantity{Float64, D, U} where {D, U}}, NonlinearFunction{true, ModelingToolkit.var"#f#302"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xd260deed, 0x7dc90363, 0x08fe8abd, 0xd742811b, 0x0b2c6f74)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x73c43a35, 0x06c47ba7, 0xb1d644a0, 0xa2947d04, 0xc57c3159)}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, ModelingToolkit.var"#generated_observed#305"{NonlinearSystem, Dict{Any, Any}}, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::NewtonRaphson{12, true, DataType, NonlinearSolve.DefaultLinSolve}; alias_u0::Bool, maxiters::Int64, tol::Float64, internalnorm::Function, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ NonlinearSolve ~/.julia/packages/NonlinearSolve/WjsHK/src/solve.jl:52
 [13] init
    @ ~/.julia/packages/NonlinearSolve/WjsHK/src/solve.jl:43 [inlined]
 [14] solve(::NonlinearProblem{Vector{Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}, true, Vector{Quantity{Float64, D, U} where {D, U}}, NonlinearFunction{true, ModelingToolkit.var"#f#302"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xd260deed, 0x7dc90363, 0x08fe8abd, 0xd742811b, 0x0b2c6f74)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x73c43a35, 0x06c47ba7, 0xb1d644a0, 0xa2947d04, 0xc57c3159)}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, ModelingToolkit.var"#generated_observed#305"{NonlinearSystem, Dict{Any, Any}}, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::NewtonRaphson{12, true, DataType, NonlinearSolve.DefaultLinSolve}; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ NonlinearSolve ~/.julia/packages/NonlinearSolve/WjsHK/src/solve.jl:4
 [15] solve(::NonlinearProblem{Vector{Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}, true, Vector{Quantity{Float64, D, U} where {D, U}}, NonlinearFunction{true, ModelingToolkit.var"#f#302"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xd260deed, 0x7dc90363, 0x08fe8abd, 0xd742811b, 0x0b2c6f74)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x73c43a35, 0x06c47ba7, 0xb1d644a0, 0xa2947d04, 0xc57c3159)}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, ModelingToolkit.var"#generated_observed#305"{NonlinearSystem, Dict{Any, Any}}, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::NewtonRaphson{12, true, DataType, NonlinearSolve.DefaultLinSolve})
    @ NonlinearSolve ~/.julia/packages/NonlinearSolve/WjsHK/src/solve.jl:4
 [16] top-level scope
    @ REPL[87]:1

It looks like the culprit is in #init#10(alias_u0, maxiters, tol, internalnorm, kwargs, , prob, alg, args) at /Users/lmorton/.julia/packages/NonlinearSolve/WjsHK/src/solve.jl:35

   50    if iip
> 51      fu = zero(u)
   52      f(fu, u, p)

Here fu is:

Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}[0.0 s]

which propagates all the way down to the setindex! call that's trying to store the dimensionless result there.

[ChrisRackauckas]

Thanks! Yes, I don't think we have looked into units here so far. I think it might be really hard with some solvers because units don't play too well with linear algebra, but we could try.

[lamorton]

I'll close this as it seems to be not directly related to solvers.