Support Unitful

[Eben60]

using DiffEqOperators, Unitful

nknots = 100
h = 1.0u"mm"/(nknots+1)

Δ = CenteredDifference(2, 2, h, nknots)

results in an error

[xtalax]

What is the error?

[Eben60]

LoadError: MethodError: no method matching CenteredDifference{1}(::Int64, ::Int64, ::Quantity{Float64, 𝐋 , Unitful.FreeUnits{(mm,), 𝐋 , nothing}}, ::Int64)
Closest candidates are:
  CenteredDifference{N}(::Int64, ::Int64, ::AbstractVector{T}, ::Int64) where {T<:Real, N} at C:\Users\Eben60\.julia\packages\DiffEqOperators\tW3ze\src\derivative_operators\derivative_operator.jl:120
  CenteredDifference{N}(::Int64, ::Int64, ::AbstractVector{T}, ::Int64, ::Any) where {T<:Real, N} at C:\Users\Eben60\.julia\packages\DiffEqOperators\tW3ze\src\derivative_operators\derivative_operator.jl:120
  CenteredDifference{N}(::Int64, ::Int64, ::T, ::Int64) where {T<:Real, N} at C:\Users\Eben60\.julia\packages\DiffEqOperators\tW3ze\src\derivative_operators\derivative_operator.jl:65
  ...
Stacktrace:
 [1] CenteredDifference(::Int64, ::Vararg{Any})
   @ DiffEqOperators C:\Users\Eben60\.julia\packages\DiffEqOperators\tW3ze\src\derivative_operators\derivative_operator.jl:290
 [2] top-level scope
   @ C:\_MyOwnDocs\SoftwareDevelopment\MFPDE\MFPDE.jl\src\Heat-1D\broken-unitful.jl:6
 [3] eval
   @ .\boot.jl:373 [inlined]
 [4] include_string(mapexpr::typeof(identity), mod::Module, code::String, filename::String)
   @ Base .\loading.jl:1196
in expression starting at C:\_MyOwnDocs\SoftwareDevelopment\MFPDE\MFPDE.jl\src\Heat-1D\broken-unitful.jl:6

---

As I see, CenteredDifference is defined like this:
CenteredDifference{N}(::Int64, ::Int64, ::T, ::Int64) where {T<:Real, N}

whereas:

julia> 1.0u"mm" isa Real
false

Could probably this help? -
Union{Real, Unitful.AbstractQuantity{T}} where {T<:Real}

[xtalax]

We could allow T<:Number and warn when not <:Real, I'll take a look and see what this would take.

[Eben60]

I must admit in the meanwhile I went a bit down the rabbit hole before giving up. The parameters passed to calculate_weights() do not have compatible units. I have tried a couple of things, but without a good understanding of the whole picture it makes apparently little sense.

Still, in the file fornberg.jl the line 24 should probably be C[1,1] = one(T) or C[1,1] = oneunit(T) instead of C[1,1] = 1 . Also, line 21 - for Unitful types, one() and oneunit() are not the same - which ist the proper one here? In the derivative_operator.jl file, it is always oneunit().

[xtalax]

I really appreciate your input here, I'm not sure what's better.
I see that:

julia> using Unitful

julia> T = typeof(1000u"g")
Quantity{Int64, 𝐌, Unitful.FreeUnits{(g,), 𝐌, nothing}}

julia> one(T)
1

julia> oneunit(T)
1 g

julia> typeof(one(T))
Int64

julia> typeof(oneunit(T))
Quantity{Int64, 𝐌, Unitful.FreeUnits{(g,), 𝐌, nothing}}

If usual math methods promote from usual Real types to the type of the quantity , one could work as long as your u0 and bcs are defined with Unitful quantities.

As the following is the case:

julia> 1*1u"g"
1 g

julia> 1+1u"g"
ERROR: DimensionError: 1 and 1 g are not dimensionally compatible.
Stacktrace:
 [1] +(x::Quantity{Int64, NoDims, Unitful.FreeUnits{(), NoDims, nothing}}, y::Quantity{Int64, 𝐌, Unitful.FreeUnits{(g,), 𝐌, nothing}})
   @ Unitful ~/.julia/packages/Unitful/woO6b/src/quantities.jl:132
 [2] +(x::Int64, y::Quantity{Int64, 𝐌, Unitful.FreeUnits{(g,), 𝐌, nothing}})
   @ Base ./promotion.jl:379
 [3] top-level scope
   @ REPL[9]:1

I will need to check for cases where values are added and make sure these are using oneunit. I think the operators are mostly applied with dot so this may not be an issue. The only place I can think of that needs a change here is in the affine part of the BC operators.

[Eben60]

The proper way would be check the unit type of each and every variable - as a physist I fully agree here with Tim Holy (see his following posts, too).

It appears, there are two type of "x" variables here: those with the same units as dx (which could typically be some length unit, let's assume it to be in cm), and dimensionless ("unitless"). So, low_boundary_x appears to be unitful, interior_x unitless. Is that the case? - then it would help to rename them so you can distinguish between them on the first glance.

Further, the coefficients as returned by calculate_weights() could probably have some other units (dimensionless? cm? 1/cm?). How do they scale with dx? Are they all the same units?

x0 and x , as passed to calculate_weights(), are both in cm . What are the units of the coefficients c1..c5? c3, c4, c5 are apparently in cm (line 22, 29, 33). What is c2? cm or cm^2 or cm^n? Here, at the line 34, I've lost the trail.