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 2 years 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

Running into the same issue - unit validation passes on the PDESystem but not on the discretized system. I can post an MWE for my case if it'd help

aditya-sengupta avatar Jul 23 '23 01:07 aditya-sengupta

What would be the successful case, that the result array has the right units?

xtalax avatar Aug 09 '23 15:08 xtalax

Potentially it'd be running unit validation on the PDESystem and trusting that that correctness is preserved on the discretization. The units in the result array being correct is ensured by the user declaring it as a variable with those units and it passing the PDESystem validation - it seems like it'd run just as correctly and with fewer hurdles if the discretized system didn't have to be unit checked (or we figure out autopromotion rules for constants to unit aware constants, but that seems harder)

aditya-sengupta avatar Aug 09 '23 15:08 aditya-sengupta

I guess the immediate successful case would be that the example above would run without error. As far as I can tell it's just an issue of the constants that are added being units (e.g. when the upwind scheme checks for u > 0, the 0 is unitless but u is not).

ctessum avatar Aug 14 '23 22:08 ctessum