physical units

Stream: Natural sciences

[Tomas Skrivan]

I got interested in physical units again and wanted to see if I can state Navier-Stokes equation with units. Hopefully someone will find this interesting :) so here is what I got:

Let's start with a simple kinematics, let's introduce time, velocity and acceleration:

variable
  (t : ℝ⟦second⟧)
  (v : ℝ⟦meter*second^(-1:Int)⟧)
  (a : ℝ⟦meter*second^(-2:Int)⟧)

now some simple formulas from kinematics do typecheck

#check (v*t) + ((1/2:ℝ)*a*t*t)  -- ℝ⟦meter⟧

To do Navier-Stoke equation we need differentiation. The type of differentiation respecting units would look like this:

def uderiv {X Y : Type} {u v w : PhysicalUnit} (f : X⟦u⟧ → Y⟦v⟧) (x : X⟦u⟧) (dx : X⟦w⟧) : Y⟦(v*w)/u⟧ := sorry

now we can state Navier-Stokes equation(in one dimension) as

variable
  (u : ℝ⟦second⟧ → ℝ⟦meter⟧ → ℝ⟦meter*second^(-1:Int)⟧)
  (p : ℝ⟦second⟧ → ℝ⟦meter⟧ → ℝ⟦pascal⟧)
  (density : ℝ⟦second⟧ → ℝ⟦meter⟧ → ℝ⟦kilogram*meter^(-3:Int)⟧)
  (ν : ℝ⟦(meter^(2:Int)) * second^(-1:Int)⟧)

#check fun t x => (uderiv (u · x) t (1:ℝ⟦1⟧))                              -- ∂u/∂t
                + (uderiv (u t ·) x (u t x))                               -- (u⬝∇)u
                + (((1:ℝ⟦1⟧) / density t x) * uderiv (p t ·) x (1:ℝ⟦1⟧))   -- 1/ρ ∇p
                + ν * uderiv ((uderiv (u t ·)) · (1:ℝ⟦1⟧)) x (1:ℝ⟦1⟧)      -- ν Δu
                =
                0

---

How does this work? The key piece is to have a custom elaborator of x * y and x + y which calls norm_num on the types of x and y and then checks if the units are the same. Have a look at the code here. The unfortunate thing is that the elaborators adds lots of cast to the expressions which will probably be quite painful to work with later on but I was unable to figure out how to do it without casts.

[Max Bobbin]

is X⟦u⟧ an actual Lean thing or did you define a notation to get this? I was having a discussion in another channel about defining physical units. Can you send how you define PhysicalUnit?

[Max Bobbin]

The basic idea that was recommended was def PhysicalUnits (α) := monoid_algebra ℝ (dimension α) and this makes use of the dimensional analysis already defined and just adds a value to a dimension.

[Eric Wieser]

Can you send how you define

It's line 10 of the link above

[Eric Wieser]

I prefer Tomas Skrivan's approach of discarding nonsense like 1m + 2s in the type system

[Tomas Skrivan]

Max Bobbin said:

The basic idea that was recommended was def PhysicalUnits (α) := monoid_algebra ℝ (dimension α) and this makes use of the dimensional analysis already defined and just adds a value to a dimension.

Yes I will generalize it like that at some point. Right now, it is just a proof of concept. I was mainly interested in how to do some complicated reduction in the types to check that two terms have the same units.

[Tomas Skrivan]

The notation X⟦u⟧ is defined on the line 131

[Eric Wieser]

Can you outline what makes the reduction complicated? Isn't everything just defeq in your example?

[Tomas Skrivan]

I explicitly call norm_num and I was having a hard time replying purely on defeq. I will have to think for some time to give you a concrete example where it is obvious that defeq is not enough.

[Eric Wieser]

It stops being defeq as soon as your powers are free variables

[Tomas Skrivan]

Yes that is the exact problem I got when I tried to define an adjoint of differential with units. I had to prove a * b * c / b / a = c.

[Tomas Skrivan]

But even with concrete units it didn't work very well.

[Tomas Skrivan]

Ohh I know, synthesizing HAdd A[u] A[v] ?? is kind of problematic if u and v are defeq.

[Eric Wieser]

HAdd strikes again

[Tomas Skrivan]

But that can be solved with a custom elaborator, however there is still the issue with writing functions polymorphic in units. There units appear as free variables and you can't rely on defeq.

[Tomas Skrivan]

Ohh the main problem is that operations over rationals are not defeq, 3+2 is not defeq 5 or 3*2/3 is not defeq to 2 over rationals.

You can probably rely on defeq when checking dimensions where you are dealing only with integer powers. If you want do to units, you have to deal with the scaling which has to be at least rational number.

[Eric Wieser]

This is because docs#Rat.ofScientific is an irreducible_def because apparently you don't want to unfold it

[Eric Wieser]

In fact pretty much all of the rationals are built this way

[Tomas Skrivan]

1+1 =?= 2 does not contain any ofScientific and still fails

open Lean Meta Qq in
#eval show MetaM Unit from do
  let a := q((1:ℚ) + (1:ℚ))
  let b := q((2:ℚ))
  IO.println s!"Is {← ppExpr a} defeq to {b}: {← isDefEq a b}"

[Eric Wieser]

Right, https://github.com/leanprover/std4/commit/04b3c9831e0c141713a70e68af7e40973ec9a1ff is what's really to blame

[Tomas Skrivan]

Hmm, anyway I can't rely on defeq because I need to deal with units as free variables.

[Tomas Skrivan]

Btw. I got interested in physical units again because of the latest talk at Topos Institute.

[Tomas Skrivan]

Max Bobbin said:

The basic idea that was recommended was def PhysicalUnits (α) := monoid_algebra ℝ (dimension α) and this makes use of the dimensional analysis already defined and just adds a value to a dimension.

I do not like the idea of using monoid_algebra. As Eric mentioned, I don't want the type system to allow nonsense like 1m + 2s. I even want to restrict stuff like 1m + 1cm and probably require an explicit cast to meters or centimeters.

[Adomas Baliuka]

I'd like to add that many people use units like √Hz. That's also not in monoid_algebra and I think should be allowed, if only to be inclusive.

[Eric Wieser]

That works just fine with MonoidAlgebra ℚ Real