Physical Units: Dimensional Safety in Simulation

println("Julia ", VERSION)

Julia 1.12.6

Overview

Engineers lose spacecraft and mix up inches with millimetres when numbers have no units. Rasmah ships a native physical-units system so a length never silently becomes a pressure: every quantity carries its physical dimensions, and the system raises an error the moment you try to do something dimensionally nonsensical (add a length to a pressure, compare a force to a temperature).

The key idea is separation of concerns:

Dimensions are the physics — which of the seven SI base quantities (length, mass, time, electric current, thermodynamic temperature, amount of substance, luminous intensity), plus plane angle, a value measures.
The unit system is the presentation — how bare numbers are interpreted at the boundary and how quantities are printed.

Internally everything is stored in SI base units as a plain Float64, so a Quantity never enters the solver kernel or the automatic-differentiation path. That is what makes units “free” — you get dimensional safety without paying for it in the hot loop.

Quantities and the unit string macro

A Quantity is a number tagged with a Dimensions value. The u"…" string macro builds units, and juxtaposing a number with a unit builds a quantity:

using Rasmah

l = 10u"mm"        # a length
F = 5u"kN"         # a force
T = 25u"°C"        # an *affine* temperature (offset + scale)
println(l)         # printed in the active unit system
println(F)
println(T)

0.01 m
5000.0 N
298.15 K

A quantity is always stored in SI. ustrip gives the raw SI number, and uconvert re-expresses it in any other unit:

println("10 mm in metres = ", ustrip(l))             # 0.01
println("10 mm in inches = ", uconvert(u"inch", l))   # 0.3937…

10 mm in metres = 0.01
10 mm in inches = 0.3937007874015748

Dimensions and arithmetic

Dimensions are the eight base exponents (the seven SI base quantities plus plane angle), so arithmetic is dimension-aware: multiplying lengths divides out, dividing a force by an area gives a pressure, and comparing quantities compares SI magnitudes.

d_force    = udim(1u"N")                                  # kg·m/s²
d_mass_acc = udim(1u"kg") * udim(1u"m") / udim(1u"s")^2
println("N ≡ kg·m/s² ? ", d_force == d_mass_acc)         # true

p = 5u"kN" / (2u"m" * 2u"m")                             # a pressure
println("p = ", p)
println("dimension of p = ", si_unit_string(udim(p)))    # "kg/(m·s^2)" = pascal

N ≡ kg·m/s² ? true
p = 1250.0 Pa
dimension of p = kg/(m·s^2)

Trying to combine inconsistent dimensions is an error, not a silent wrong answer:

try
    1u"m" + 1u"s"
catch e
    println("caught: ", typeof(e))                       # DimensionError
end

caught: DimensionError

Parsing unit strings

uparse turns a unit expression into its (dimensions, factor, offset) triple, and unit_name / unit_info describe a unit in words:

println(uparse("kg/m^3"))          # (Dimensions, 1.0, 0.0)
println(unit_name(u"GPa"))         # "gigapascal"
println(unit_info(u"MPa"))         # "MPa — megapascal — dimension: kg/(m·s^2)"

(Dimensions(-3//1, 1//1, 0//1, 0//1, 0//1, 0//1, 0//1, 0//1), 1.0, 0.0)
gigapascal
MPa — megapascal — dimension: kg/(m·s^2)

Any scale: SI prefixes and simplification

The whole SI prefix family is available as a modifier on any base unit — from yocto (y, 1e-24) up to yotta (Y, 1e24) — so you can describe any magnitude directly, from a femtometre to a yottametre:

println(ustrip(1u"fm"))         # 1.0e-15  (femtometre)
println(ustrip(1u"pm"))         # 1.0e-12  (picometre)
println(ustrip(1u"nm"))         # 1.0e-9   (nanometre)
println(ustrip(1u"µm"))         # 1.0e-6   (micrometre)
println(ustrip(1u"km"))         # 1000.0   (kilometre)
println(unit_name(u"Em"))       # "exameter"
println(ustrip(1u"Ym"))         # 1.0e24   (yottametre)

1.0e-15
1.0e-12
1.0e-9
1.0e-6
1000.0
exameter
1.0e24

When a quantity lands at an awkward scale, simplify re-expresses it at the nearest power-of-1000 prefix so the number reads naturally — 1000 micro is 1 milli, and 1500 metres is 1.5 kilometres:

println(simplify(1000u"µm"))    # "1.0 mm"
println(simplify(1500u"m"))     # "1.5 km"
println(simplify(0.001u"s"))    # "1.0 ms"
println(simplify(2500u"N"))     # "2.5 kN"
println(simplify(1000u"kg"))    # "1.0 Mg"  (mass prefixes attach to the gram)

1.0 mm
1.5 km
1.0 ms
2.5 kN
1.0 Mg

simplify returns a ScaledQuantity — a value plus an explicit display unit. The stored value is still SI underneath, so it round-trips back to a plain Quantity with nothing lost:

s = simplify(1000u"µm")
println(s)                         # "1.0 mm"
println(ustrip(s))                 # 0.001
println(Quantity(s) ≈ 1000u"µm")   # true

1.0 mm
0.001
true

The active unit system

The unit system controls how bare numbers are read and how quantities are rendered. unitsystem() reports it, and set_unitsystem! swaps it (SI is the canonical internal system regardless):

println("default length unit = ", unitsystem().length)   # "m"
set_unitsystem!(; length="mm", mass="tonne")
println(unitsystem().length, " / ", unitsystem().mass)   # "mm" / "tonne"
println(1000u"mm")                                       # now prints "1000 mm"
set_unitsystem!(; length="m", mass="kg")                 # restore SI

default length unit = m
mm / tonne
1000.0 mm

UnitSystem(m, kg, s, A, K, mol, cd, deg, N, Pa, Pa, J, W, m/s, kg/m^3, kg/(m·s), Hz)

Unit-aware material properties

The materials library is unit-aware: every property knows its physical dimension, so a material can be constructed with a Quantity and inspected with one.

steel = material("steel")
println("E of steel = ", quantity(steel, :youngs_modulus))   # e.g. 200 GPa
println(property_dimension(:youngs_modulus) == udim(1u"Pa")) # true (a pressure)

E of steel = 2.0e11 Pa
true

The verification is the point: E (a modulus) is dimensionally a pressure, and the system proves it rather than assuming it.

Alternative and imperial units

The active unit system picks a display unit for each base and derived quantity, so the same SI value can be shown in kelvin, Celsius, or Fahrenheit — or a length in metres, millimetres, or inches. Swap just the quantities you care about:

set_unitsystem!(; temperature="°C")
println(293.15u"K")                                  # "20.0 °C" (room temperature)
println(round(uconvert(u"°F", 273.15u"K"); digits=2)) # 32.0 (the freezing point)
println(round(uconvert(u"°R", 273.15u"K"); digits=2)) # 491.67 (Rankine)
set_unitsystem!(; temperature="K")                   # back to kelvin

20.0 °C
32.0
491.67

UnitSystem(m, kg, s, A, K, mol, cd, deg, N, Pa, Pa, J, W, m/s, kg/m^3, kg/(m·s), Hz)

Temperature scales are affine: °C and °F carry an offset, so 25u"°C" is 298.15 kelvin — a shift, not just a rescale. °R (Rankine) is the absolute Fahrenheit-scaled counterpart of kelvin.

Plane angle is a tracked dimension, so radians and degrees are interchangeable — pick whichever you want to see, while the value stays in radians internally:

θ = (pi / 2)u"rad"
set_unitsystem!(; angle="deg")
println(θ)                                    # "90.0 deg"
set_unitsystem!(; angle="rad")
println(θ)                                    # "1.5707963267948966 rad"
set_unitsystem!(; angle="deg")
println(ustrip(180u"deg"))                    # 3.141592653589793 (π in radians)

90.0 deg
1.5707963267948966 rad
3.141592653589793

Imperial and US-customary units are built in — inch, foot, yard, mile, ounce, pound, kip (1000 lbf), psi/ksi, horsepower, BTU, gallon, mph, rpm, and more:

println(ustrip(1u"inch"))        # 0.0254
println(ustrip(1u"ksi"))         # 6.894757293168361e6
println(ustrip(1u"hp"))          # 745.6998715822702
println(unit_name(u"BTU"))       # "British thermal unit"

0.0254
6.894757293168361e6
745.6998715822702
British thermal unit

For a full switch to the inch–pound–second system, use the ready-made unit system:

set_unitsystem!(IN_LBF_S)
println(1u"inch")                            # "1.0 inch"
println(1u"lbf")                             # "1.0 lbf"
println(1u"psi")                             # "1.0 psi"
set_unitsystem!(SI_UNITS)                    # back to SI

1.0 inch
1.0 lbf
1.0 psi

UnitSystem(m, kg, s, A, K, mol, cd, deg, N, Pa, Pa, J, W, m/s, kg/m^3, kg/(m·s), Hz)

Nothing changes internally — the value is still stored in SI. Only the display moves between unit systems, so a solve never depends on which units you chose to look at it in.

Display rounding

By default a Quantity prints at full Float64 precision. Imperial and affine temperature conversions can carry floating-point noise (273.15 K → 31.999999999999943 °F); set_sigfigs! rounds the display to a chosen number of significant digits — the stored value is never touched:

set_unitsystem!(; temperature="°F")
set_sigfigs!(6)
println(273.15u"K")      # "32.0 °F" (clean)
println(sigfigs())       # 6
set_sigfigs!(nothing)    # back to full precision
println(273.15u"K")      # "31.999999999999943 °F"
set_unitsystem!(; temperature="K")

32.0 °F
6
31.999999999999943 °F

UnitSystem(m, kg, s, A, K, mol, cd, deg, N, Pa, Pa, J, W, m/s, kg/m^3, kg/(m·s), Hz)

set_sigfigs! also governs the show of any quantity and the colorbar labels that vtk(...) renders for a Quantity field.

Units flow through the workflow

Unit awareness threads through the modeling → simulate → visualise pipeline. Geometry constructors accept Quantity (stripped to SI internally, so the solve is unchanged), and vtk derives axis/colorbar units from a Quantity field:

part = box(10u"mm", 20u"mm", 30u"mm")        # a 10×20×30 mm box
mesh = tetrahedralize(part, 4)

# a field carrying displacement dimensions renders in the active length unit
set_unitsystem!(; length="mm")
u = QuantityVector(collect(1.0:length(mesh.nodes)), udim(1u"m"))  # (demo)
# vtk(mesh; field=u, fieldname="displacement")   # colorbar "displacement [mm]"
set_unitsystem!(SI_UNITS)

UnitSystem(m, kg, s, A, K, mol, cd, deg, N, Pa, Pa, J, W, m/s, kg/m^3, kg/(m·s), Hz)

Units flow out of a solve

A solve(...) result is unit-tagged too — each solved field is a QuantityVector carrying its physical dimension, so a displacement prints in metres, a temperature in kelvin, and a potential in volts. vtk(...) reads that dimension and labels the colorbar automatically:

part = box(10u"mm", 10u"mm", 40u"mm")
b = brep(part)
left  = faces_in_box(b, [-6e-3, -1.0, -1.0], [-4e-3, 1.0, 1.0])
right = faces_in_box(b, [ 4e-3, -1.0, -1.0], [ 6e-3, 1.0, 1.0])
mesh = tetrahedralize(b, 4)

u = solve(Elasticity(), mesh, material("steel");
          fixed = left => [0, 0, 0], loads = right => [100u"kN", 0, 0])
println(typeof(u))                          # QuantityVector
println("max |u| = ", round(maximum(abs, u); sigdigits = 2))   # a length
# vtk(mesh; field = u, fieldname = "displacement")   # colorbar "displacement [m]"

QuantityVector{Float64}
max |u| = 5.4e-9 m

Multi-field solves return a named tuple of QuantityVector fields, each with its own dimension:

r = solve(ThermoElasticity(), mesh, material("steel");
          fixed_u = left => [0, 0, 0], loads_u = right => [1e6, 0, 0], fixed_T = [1])
println(udim(r.u))   # length (metres)
println(udim(r.T))   # temperature (kelvin)

Dimensions(1//1, 0//1, 0//1, 0//1, 0//1, 0//1, 0//1, 0//1)
Dimensions(0//1, 0//1, 0//1, 0//1, 1//1, 0//1, 0//1, 0//1)

Verified inputs, normalized solves, and result round-trips

prepare_analysis checks that loads and boundary values carry the right dimensions before anything is solved — a length passed where a force is expected is an error, not a silent wrong answer:

part = box(10u"mm", 10u"mm", 40u"mm")
b = brep(part)
left  = faces_in_box(b, [-6e-3, -1.0, -1.0], [-4e-3, 1.0, 1.0])
right = faces_in_box(b, [ 4e-3, -1.0, -1.0], [ 6e-3, 1.0, 1.0])
mesh = tetrahedralize(b, 4)

try
    prepare_analysis(Elasticity(), mesh, material("steel");
                     fixed = left => [0, 0, 0], loads = right => [1u"mm", 0, 0])
catch e
    println("caught: ", typeof(e))        # DimensionError (a length, not a force)
end

caught: DimensionError

The prepared analysis also computes reference scales and, when solved, runs the linear system in a non-dimensionalized form — matrix entries are O(1) whether you model in metres+pascals or millimetres+megapascals — while returning the same physical answer:

p = prepare_analysis(Elasticity(), mesh, material("steel");
                     fixed = left => [0, 0, 0], loads = right => [100u"kN", 0, 0])
u  = solve(p)                                    # normalized solve
u0 = solve(Elasticity(), mesh, material("steel");
           fixed = left => [0, 0, 0], loads = right => [100u"kN", 0, 0])
println("normalized == plain ? ", u ≈ u0)        # true
println("reference length L0 = ", round(p.L0; sigdigits = 3), " m")

normalized == plain ? true
reference length L0 = 0.0424 m

extract_results re-expresses a solve’s fields in the unit system you choose — here the same displacement printed in millimetres and in metres:

println("max |u| = ", round(maximum(abs, extract_results(u, UnitSystem(length = "mm"))); sigdigits = 3), " mm")
println("max |u| = ", round(maximum(abs, u); sigdigits = 3))    # a length, printed in the active unit

max |u| = 5.35e-6 mm
max |u| = 5.35e-9 m

Units at the file boundary

Formats that carry no unit of their own (STL, OBJ, PLY, MSH, INP, VTK) take a unit keyword on read and write; the value is converted to SI on the way in and back on the way out. Omit unit and the active unit system’s length unit is used (metres in SI):

io = IOBuffer()
write_stl(surface_mesh(part), io; unit = :mm)   # written in millimetres
seekstart(io)
m = read_stl(io; unit = :mm)                    # read back as millimetres → SI
println("max node = ", round(maximum(m.vertices) * 1000; digits = 1), " mm")

max node = 20.0 mm

Formats that do declare units write them: STEP records its header length unit (write_step(b; length_unit = :mm), read_step_length_unit), and 3MF writes its unit attribute (write_3mf(m; unit = "millimeter")), rescaling to match.

Next steps

Materials and multiphysics properties — the material library that the unit system attaches to.
Getting started: differentiable geometry — the features and fields that the units annotate.