Notebook 12: Working with Physical Units

Run this notebook on Binder

Underworld3 has built-in support for physical units throughout the modeling workflow. This makes your models easier to understand and helps catch dimensional errors early.

In this notebook you’ll learn:

  • Creating physical quantities (temperatures, velocities, viscosities)

  • Converting between units

  • Working with unit-aware arrays and coordinates

  • Automatic unit tracking through derivatives

import nest_asyncio
nest_asyncio.apply()

import underworld3 as uw
import numpy as np
import sympy

Creating Physical Quantities

Create quantities using uw.units, which provides access to the Pint units library. You can write units explicitly or use strings:

# Explicit units
plate_velocity = 5 * uw.units.cm / uw.units.year
mantle_temperature = 1500 * uw.units.K

# String notation (Pint can parse these)
mantle_depth = 2900 * uw.units("km")
mantle_viscosity = 1e21 * uw.units("Pa*s")

# Let's see !
mantle_temperature

# View all the quantities
plate_velocity, mantle_depth, mantle_viscosity

Unit Conversions

Converting between compatible units is straightforward using the .to() method:

# Convert velocity to different units
plate_velocity.to("mm/year")

# Convert to SI units
plate_velocity.to("m/s")

# Convert depth
mantle_depth.to("m")

Optional: Setting Up Reference Quantities

Note: You can use units in Underworld3 without setting reference quantities. However, setting them is strongly recommended for better numerical conditioning in solvers.

Reference quantities define the characteristic scales for your problem using any domain-specific terminology. The system automatically derives fundamental scales [L], [M], [T], [θ] through dimensional analysis. For example:

  • Geoscience: domain_depth, plate_velocity, mantle_viscosity

  • Aerodynamics: wingspan, airspeed, air_density

  • Any domain: Use whatever makes sense for YOUR problem

When set, reference quantities ensure your variables have proper scaling coefficients (e.g., scaling_coefficient=1/1500 for a 1500 K temperature reference scale). This improves solver performance and conditioning.

CRITICAL TIMING: Reference quantities must be set BEFORE creating any mesh, because the model locks units once a mesh is created. The mesh domain size is just geometry - it does NOT define the physical length scale. You must provide that explicitly.

Without reference quantities: Variables will use scaling_coefficient=1.0, which can lead to poorly conditioned systems if your physical values span many orders of magnitude.

With reference quantities: Variables automatically get appropriate scaling factors that normalize values to order-1 numbers, improving numerical stability.

# Set up reference quantities BEFORE creating the mesh
# Note: These define the CHARACTERISTIC SCALES for your problem,
# NOT the mesh size!
# 
# Parameter names don't matter for dimensional analysis (system is domain-agnostic),
# but validation expects specific names for recommended use.
model = uw.get_default_model()
model.set_reference_quantities(
    # Standard parameter names (for solver validation and conditioning hints)
    temperature_difference=uw.quantity(1500, "K"),      # Temperature scale
    plate_velocity=uw.quantity(5, "cm/year"),           # Velocity scale
    domain_depth=uw.quantity(1000, "km"),               # Length scale
    mantle_viscosity=uw.quantity(1e21, "Pa*s"),         # Viscosity scale
)

# Get the derived reference scales
scales = model._fundamental_scales
print("Derived reference scales (from dimensional analysis):")
print(f"  Length:      {scales['length']}")
print(f"  Time:        {scales.get('time')}")
print(f"  Temperature: {scales.get('temperature')}")
print()

# Now create a mesh with EXPLICIT UNITS
# These are SEPARATE from the reference scales!
# Example: Earth-like dimensions (6370 km radius)
mesh = uw.meshing.UnstructuredSimplexBox(
    minCoords=(0.0 * uw.units.km, 0.0 * uw.units.km),
    maxCoords=(6370.0 * uw.units.km, 3185.0 * uw.units.km),  # Earth radius
    cellSize=200.0 * uw.units.km,                            # 200 km cells
    qdegree=2,
)

print(f"Mesh domain: (0-6370 km) x (0-3185 km)")
print(f"  Mesh coordinates have units: {uw.get_units(mesh.X[0])}")
print()
print(f"Key relationship:")
print(f"  Reference length scale: 1000 km")
print(f"  Mesh maximum: 6370 km")
print(f"  Ratio: 6370/1000 = 6.37 reference units")
print()


x,y = mesh.X

mesh.dm.getCoordinates().array

mesh.X.coords.to("km")

# Create mesh variables with units
# Mesh coordinates automatically inherited units from mesh creation
temperature = uw.discretisation.MeshVariable("T", mesh, 1, degree=2, units="K")
velocity = uw.discretisation.MeshVariable("u", mesh, 2, degree=2, units="m/s")

print("Variable units:")
print(f"  temperature.units = {temperature.units}")
print(f"  velocity.units = {velocity.units}")
print()

print("Coordinate units (automatically set by mesh):")
print(f"  mesh.X[0] units = {uw.get_units(mesh.X[0])}")
print(f"  mesh.X[1] units = {uw.get_units(mesh.X[1])}")
print()

print("Derivatives automatically get correct units:")
dT_dy = temperature.sym.diff(mesh.X[1])
print(f"  dT/dy units = {uw.get_units(dT_dy)}")
print()

print("This is K / km because:")
print(f"  Temperature is in K (from variable definition)")
print(f"  y-coordinate is in km (from mesh)")
print(f"  Dimensional analysis: K / km ✅")

mesh.X.coords

Understanding Reference Scales vs. Mesh Coordinates

Critical distinctions:

  1. Reference scales (temperature_difference=1500 K):

    • Define the characteristic physical size of your PROBLEM

    • Used to normalize variables for numerical conditioning

    • Independent of mesh size

    • Example: “The mantle convection problem I’m studying has temperature variations of ~1500 K”

  2. Mesh size (maxCoords=6370 km):

    • The actual geometric extent of your computational domain

    • Can be larger, smaller, or equal to reference scales

    • Expressed in explicit units (km, m, etc.)

    • Example: “My mesh extends to Earth’s radius: 6370 km”

The system automatically handles the conversion:

  • Mesh maxCoords = 6370 km = 6.37 × (reference scale of 1000 km)

  • Mesh coordinates automatically get units from mesh creation (km in this case)

  • Derivatives automatically have correct units through dimensional analysis (K/km for temperature gradients)

  • All this happens automatically! ✅

About parameter names:

  • The dimensional analysis is truly domain-agnostic: You could use wingspan, airspeed, experiment_duration, etc.

  • However, the validation system recognizes specific parameter names to provide helpful warnings and scaling hints

  • Best practice: Use standard names (temperature_difference, plate_velocity, etc.) to get proper solver conditioning

  • If you use custom names, set the reference quantities and ignore validation warnings if they don’t apply to your domain

Why separate reference scales and mesh size?

Imagine two different problems with the same mesh (6370 km):

  1. Small mantle plume: Reference scale = 100 km → Mesh = 63.7 reference units

  2. Global convection: Reference scale = 2000 km → Mesh = 3.18 reference units

The solver sees different problem sizes through reference scales, even with identical mesh geometry. This flexibility lets you reuse the same mesh for different problem scales.

Initializing Fields

Set up a simple temperature field:

# Get coordinate symbols
x, y = mesh.X

# Initialize temperature: T = 300 + 2.6*y (K)
temperature.array[...] = uw.function.evaluate(
    300.0 + 2.6 * y, temperature.coords
)

velocity.array[...] = uw.function.evaluate(
    sympy.Matrix([5.0, 0.0]), velocity.coords
).reshape(velocity.array.shape)

temperature.min(), temperature.max()

uw.get_units(temperature * velocity[0])

temperature.max()

model.get_fundamental_scales()

1/1500

uw.get_units(model.to_model_magnitude(velocity))

x,y = mesh.X.with_units
uw.get_units(x)

uw.get_units(temperature*temperature)

type(x.sym)

L0 = model.get_fundamental_scales()["length"]

uw.units_of(x)

temperature.diff(x).units

uw.units_of(temperature.diff(x))

Unit-Aware Operations

Variables with units support mathematical operations that preserve dimensional consistency:

# Get units from expressions
uw.get_units(temperature)

%%timeit

# Check dimensionality
uw.get_dimensionality(temperature / velocity[0])

# %%timeit

# Units work naturally with arithmetic
nd = (uw.non_dimensionalise(temperature / velocity[0]))

(temperature/velocity[0]).units

uw.non_dimensionalise(temperature)

uw.get_units(2 * temperature * velocity[0])

Automatic Unit Tracking Through Derivatives

When you take derivatives, units are tracked automatically:

# Derivative automatically has correct units
dTdy = temperature.diff(y)

# The derivative has units!
uw.get_units(dTdy)

uw.get_units(dTdy[0,0])

The derivative has units that make physical sense - temperature units divided by coordinate units.

Computing Gradients with Projection

To get numerical gradient values at mesh nodes, use the Projection system:

# Create a variable to hold the gradient
gradT = uw.discretisation.MeshVariable(
    "gradT",
    mesh,
    1,
    degree=1,
    units="K/m",  # Specify units directly
)

# Project the derivative onto the mesh
proj = uw.systems.Projection(mesh, gradT)
proj.uw_function = temperature.diff(y)[0]
proj.solve()

gradT.min(), gradT.max()

# Gradient array is unit-aware
gradT.array[100, ...]

# Unit conversion example
# Note: Array views don't have .to() method, but we can extract units info
print("Gradient values:")
print(f"  gradT[0,0,0] = {gradT.array[0, 0, 0]}")
print(f"  Units: {uw.get_units(gradT)}")
print()

# To convert units, use a Pint quantity instead
gradT_value_km_per_K = 2.6 * uw.units.K / uw.units.km  # Direct Pint quantity

# Convert to K/mm
gradT_value_mm_per_K = gradT_value_km_per_K.to("K/mm")
print(f"Gradient conversion:")
print(f"  2.6 K/km = {gradT_value_mm_per_K}")
print()
print("(MeshArrayView doesn't support .to() directly, but you can work with")
print(" extracted values using Pint quantities as shown above)")

# SimpleMeshArrayView doesn't support .to() method directly
# But we can work with Pint quantities for unit conversion

print("Gradient array values:")
print(f"  gradT.array[0, 0, 0] = {gradT.array[0, 0, 0]}")
print(f"  Units: {uw.get_units(gradT)}")
print()

# Create a Pint quantity with the correct value and units
gradT_value_in_K_per_m = 2.6 * uw.units.K / uw.units.m

# Convert to K/km
gradT_in_K_per_km = gradT_value_in_K_per_m.to("K/km")
print(f"Unit conversion example:")
print(f"  2.6 K/m = {gradT_in_K_per_km}")
print()
print("Note: To convert array values, extract individual values and use")
print("      Pint quantities as shown above.")

gradT.to("K/km")

uw.get_dimensionality(uw.units("K/cm"))

Dimensional Analysis

Units help catch errors. For example, trying to add quantities with incompatible dimensions will fail:

# This works - same dimensions
total_temperature = mantle_temperature + 100 * uw.units.K
total_temperature

# This fails - incompatible dimensions
plantle_temperocity = mantle_temperature + plate_velocity

Summary

The units system in Underworld3:

  • Creation: Use uw.units to create physical quantities

  • Conversion: Use .to(target_units) to convert between compatible units

  • Derivatives: Automatically get correct units (e.g., temperature.diff(y) has units K/m)

  • Type checking: uw.get_units() and uw.get_dimensionality() let you inspect units

  • Error prevention: Incompatible unit operations raise errors

Units make your code clearer and help catch errors early!

What’s Next?

Notebook 13 introduces non-dimensional scaling and shows how to:

  • Set up reference quantities for automatic scaling

  • Solve problems in non-dimensional form for better numerical conditioning

  • Convert easily between dimensional and non-dimensional representations

This makes working with multi-scale physical problems much more robust.

Try It Yourself

Exercises to explore:

# 1. Create different quantities
density = 3300 * uw.units("kg/m^3")
gravity = 9.81 * uw.units("m/s^2")
stress = density * gravity * mantle_depth

# 2. Check the units
uw.get_units(stress)

# 3. Convert to different units
stress.to("GPa")

# 4. Create a vector derivative
div_velocity = velocity[0].diff(x) + velocity[1].diff(y)
uw.get_units(div_velocity)