MethodOfLines.jl icon indicating copy to clipboard operation
MethodOfLines.jl copied to clipboard

Units don't seem to work correctly

Open ctessum opened this issue 1 year ago • 4 comments

Hello,

Here is a small example with units:

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets
using Unitful

function ModelingToolkit.get_unit(a::Pair{Num, Float64})
    ModelingToolkit.get_unit(a[1])
end

@parameters x [unit=u"m"] 
@parameters t [unit=u"s"]
@variables u(..) [unit=u"kg"]
Dt = Differential(t)
Dx = Differential(x)

x_min = t_min = 0.0
x_max = t_max = 1.0

@parameters α = 10. [unit=u"m/s"]

eq = [Dt(u(x,t)) ~ α * Dx(u(x,t))]

@assert ModelingToolkit.get_unit(eq[1].lhs) == u"kg/s"
@assert ModelingToolkit.get_unit(eq[1].rhs) == u"kg/s"

domains = [x ∈ Interval(x_min, x_max),
              t ∈ Interval(t_min, t_max)]


bcs = [u(x,t_min) ~ 0.0,
       u(x_min,t) ~ 0.0,
       u(x_max,t) ~ sin(1.0),
] 

@named pdesys = PDESystem(eq,bcs,domains,[x,t],[u(x,t)], [α => 10.])

discretization = MOLFiniteDifference([x=>6], t)
@time prob = discretize(pdesys,discretization)

It works if there aren't any units, and as you can see the units on the two sides of the equation match. However, running the code to the end gives this error: ERROR: ModelingToolkit.ValidationError("Some equations had invalid units. See warnings for details.")

Here are the warnings in question:

┌ Warning:  in eq. #1left, in comparison >, units [m s⁻¹] and [] do not match.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/rrbUl/src/systems/validation.jl:152
┌ Warning:  in eq. #2left, in comparison >, units [m s⁻¹] and [] do not match.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/rrbUl/src/systems/validation.jl:152
┌ Warning:  in eq. #3left, in comparison >, units [m s⁻¹] and [] do not match.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/rrbUl/src/systems/validation.jl:152
┌ Warning:  in eq. #4left, in comparison >, units [m s⁻¹] and [] do not match.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/rrbUl/src/systems/validation.jl:152

It looks like what is happening is that the upwind differencing is comparing the velocity to zero, but the zero doesn't have any units, which causes the problem. I also tried advection_scheme=WENOScheme(), but it gave a similar error.

ctessum avatar Jun 25 '23 05:06 ctessum