I have written a small python module which leverages Pint, a python package used to define and manipulate physical units. Pint uses a system of unit registries to define a large set of units and most of their prefixes. A Pint Quantity object is a physical value with a specific unit attached to itself. Quantity objects can be used in mathematical expressions and return Quantities in appropriate units.
Romain Beucher [mailto:romain.beucher@unimelb.edu.au]
First of all make sure that you have Pint installed on your system. The easiest is to use pip in a terminal:
[sudo] pip install pintThe scaling submodule can then be imported as follow:
import unsupported.geodynamics.scaling as sca
The module gives you access to a Pint UnitRegistry which is basically a registry of all units available through Pint.
Lets create an alias to make the following steps easier to read:
u = sca.UnitRegistry
You can now define quantities as follow:
length = 1.0 * u.nautical_mile
length = 1.0 * u.parsec
length = 1.0 * u.kilometer
width = 2.0 * u.kilometer
The length and width objects are Pint Quantities. I strongly encourage you to explore their different attributes and methods. Use tab + tab after the . to get a list of all the options available.
print length.dimensionality
print length.magnitude
print length.units
Units are not important here. You can use whatever you want, you can even define your own:
u.define('distance_to_coffee_machine = 57.34 * meters = dtcm')
length = 3.2 * u.distance_to_coffee_machine
print length
And you can of course convert between units:
length = length.to(u.meter)
print length
Pint Quantities can be used to calculate derived quantities:
area = length * width
print "area = {0}".format(area)
Note Underworld does not accept pint Quantities as input (yet). That is where the scaling comes to play.
Lets assume that you want a build a simple thermal convection model but with some physical units instead of dimensionless values. The example solves 2D dimensionless isoviscous thermal convection with a Rayleigh number of $10^4$, see Blankenbach et al. 1989 for details.
Here are the physical dimension of our model.
# Set simulation box size.
boxHeight = 1000. * u.kilometer
boxLength = 2000. * u.kilometer
# Set min/max temperatures.
tempMin = 273.15 * u.degK
tempMax = 1273.15 * u.degK
viscosity = 2.5e19 * u.pascal * u.second
gravity = 9.81 * u.metre / u.second**2
By default the scaling module gives you a set of scaling coefficients to normalize your values. It can be accessed as follow:
sca.scaling
You can see that 5 fundamental dimensions have been defined in terms of [Length], [Mass], [Temperature], [Time] and [Substance], following Pint convention.
The scaling coefficients are all equal to 1 by default. Now you can redefine the models physical dimensions to be used by Underworld. The nonDimensinalize function ensures all units are converted into the SI unit system.
boxHeight = sca.nonDimensionalize(boxHeight)
boxLength = sca.nonDimensionalize(boxLength)
tempMin = sca.nonDimensionalize(tempMin)
tempMax = sca.nonDimensionalize(tempMax)
viscosity = sca.nonDimensionalize(viscosity)
gravity = sca.nonDimensionalize(gravity)
print "boxHeight = {0}".format(boxHeight)
print "boxLength = {0}".format(boxLength)
print "tempMin = {0}".format(tempMin)
print "tempMax = {0}".format(tempMax)
print "viscosity = {0}".format(viscosity)
print "gravity = {0}".format(gravity)
All the values above are dimensionless. They have been normalized using the values from the scaling dictionary and uniformised into the SI system.
We can use this same process to define a different scaling, with length, the temperature gradient and the viscosity are all equal to one. Let's start by redefining the model parameters.
# re-defining the parameters as above
boxHeight = 1000. * u.kilometer
boxLength = 2000. * u.kilometer
tempMin = 273.15 * u.degK
tempMax = 1273.15 * u.degK
viscosity = 2.5e19 * u.pascal * u.second
gravity = 9.81 * u.metre / u.second**2
We can know define some scaling cooefficients
# define scaling coefficients
KL = boxHeight
KT = tempMax - tempMin
Kt = (gravity**-1 * KL)**0.5
KM = viscosity * KL * Kt
Let's check the scaling coefficients
sca.scaling["[length]"] = KL.to_base_units()
sca.scaling["[temperature]"] = KT.to_base_units()
sca.scaling["[time]"] = Kt.to_base_units()
sca.scaling["[mass]"] = KM.to_base_units()
sca.scaling
We can now nondimensionalize the model parameters
boxHeight = sca.nonDimensionalize(boxHeight)
boxLength = sca.nonDimensionalize(boxLength)
tempMin = sca.nonDimensionalize(tempMin)
tempMax = sca.nonDimensionalize(tempMax)
viscosity = sca.nonDimensionalize(viscosity)
gravity = sca.nonDimensionalize(gravity)
print "boxHeight = {0}".format(boxHeight)
print "boxLength = {0}".format(boxLength)
print "tempMin = {0}".format(tempMin)
print "tempMax = {0}".format(tempMax)
print "viscosity = {0}".format(viscosity)
print "gravity = {0}".format(gravity)
To go back to dimension units, call the Dimensionalize function:
print "boxHeight in Kilometers: {0}".format(sca.Dimensionalize(boxHeight, units=u.kilometer))
print "boxHeight in dtcm: {0}".format(sca.Dimensionalize(boxHeight, units=u.dtcm))
From now on, any new quantity must be scaled using a call to sca.nonDimenionalize, this makes sure that the value is properly scaled.
thermal_conductivity = sca.nonDimensionalize(3.0 * u.watt / u.meter / u.degK)
print "Thermal Conductivity = {0}".format(thermal_conductivity)
Back to our thermal convection example, using the previously defined model parameters
import underworld as uw
from underworld import function as fn
import glucifer
import math
nd = sca.nonDimensionalize
conductivity = nd(3. * u.watt / u.meter / u.delta_degC)
density = nd(3300. * u.kilogram * u.meter**3)
capacity = nd(914. * u.joules / u.kelvin / u.kilogram)
diffusivity = conductivity / (density * capacity)
# Set the resolution.
resx = 32
resy = 16
We create the mesh
mesh = uw.mesh.FeMesh_Cartesian( elementType = ("Q1/dQ0"),
elementRes = (resx, resy),
minCoord = (0., 0.),
maxCoord = (boxLength, boxHeight))
velocityField = mesh.add_variable( nodeDofCount=2 )
pressureField = mesh.subMesh.add_variable( nodeDofCount=1 )
temperatureField = mesh.add_variable( nodeDofCount=1 )
temperatureDotField = mesh.add_variable( nodeDofCount=1 )
# Initialise values
velocityField.data[...] = nd(0. * u.meter / u.second)
pressureField.data[...] = nd(0. * u.megapascal)
temperatureDotField.data[...] = nd(0.)
# Rayleigh number.
Ra = 1.0e4
# Construct our density function.
densityFn = Ra * temperatureField
# Define our vertical unit vector using a python tuple (this will be automatically converted to a function).
z_hat = ( 0.0, gravity )
# Now create a buoyancy force vector using the density and the vertical unit vector.
buoyancyFn = densityFn * z_hat
Set a sinusoidal perturbation in the temperature field to seed the onset of convection.
pertStrength = 0.2
deltaTemp = tempMax - tempMin
for index, coord in enumerate(mesh.data):
pertCoeff = math.cos( math.pi * coord[0] ) * math.sin( math.pi * coord[1] )
temperatureField.data[index] = tempMin + deltaTemp*(boxHeight - coord[1]) + pertStrength * pertCoeff
temperatureField.data[index] = max(tempMin, min(tempMax, temperatureField.data[index]))
Set top and bottom wall temperature boundary values.
for index in mesh.specialSets["MinJ_VertexSet"]:
temperatureField.data[index] = tempMax
for index in mesh.specialSets["MaxJ_VertexSet"]:
temperatureField.data[index] = tempMin
Construct sets for I (vertical) and J (horizontal) walls.
iWalls = mesh.specialSets["MinI_VertexSet"] + mesh.specialSets["MaxI_VertexSet"]
jWalls = mesh.specialSets["MinJ_VertexSet"] + mesh.specialSets["MaxJ_VertexSet"]
Create Direchlet, or fixed value, boundary conditions. More information on setting boundary conditions can be found in the Systems section of the user guide.
# 2D velocity vector can have two Dirichlet conditions on each vertex,
# v_x is fixed on the iWalls (vertical), v_y is fixed on the jWalls (horizontal)
velBC = uw.conditions.DirichletCondition( variable = velocityField,
indexSetsPerDof = (iWalls, jWalls) )
# Temperature is held constant on the jWalls
tempBC = uw.conditions.DirichletCondition( variable = temperatureField,
indexSetsPerDof = (jWalls,) )
Render initial conditions for temperature
figtemp = glucifer.Figure( figsize=(800,400) )
figtemp.append( glucifer.objects.Surface(mesh, sca.Dimensionalize(temperatureField, u.degC), colours="blue white red") )
figtemp.append( glucifer.objects.Mesh(mesh) )
figtemp.show()
Setup a Stokes system
Underworld uses the Stokes system to solve the incompressible Stokes equations.
stokes = uw.systems.Stokes( velocityField = velocityField,
pressureField = pressureField,
conditions = velBC,
fn_viscosity = viscosity,
fn_bodyforce = buoyancyFn )
# get the default stokes equation solver
solver = uw.systems.Solver( stokes )
Set up the advective diffusive system
Underworld uses the AdvectionDiffusion system to solve the temperature field given heat transport through the velocity field. More information on the advection diffusion equation can be found here.
advDiff = uw.systems.AdvectionDiffusion( phiField = temperatureField,
phiDotField = temperatureDotField,
velocityField = velocityField,
fn_diffusivity = diffusivity,
conditions = tempBC )
surfaceArea = uw.utils.Integral(fn=1.0,mesh=mesh, integrationType='surface', surfaceIndexSet=mesh.specialSets["MaxJ_VertexSet"])
surfacePressureIntegral = uw.utils.Integral(fn=pressureField, mesh=mesh, integrationType='surface', surfaceIndexSet=mesh.specialSets["MaxJ_VertexSet"])
# define an update function
def update():
# Retrieve the maximum possible timestep for the advection-diffusion system.
dt = advDiff.get_max_dt()
# Advect using this timestep size.
advDiff.integrate(dt)
# pressure correction
(area,) = surfaceArea.evaluate()
(p0,) = surfacePressureIntegral.evaluate()
pressureField.data[:] = pressureField.data[:] -(p0 / area)
return time+dt, step+1
# init these guys
time = 0.
step = 0
steps_end = 20
# perform timestepping
while step < steps_end:
# Solve for the velocity field given the current temperature field.
solver.solve()
time, step = update()
Plot final temperature and velocity field
# plot figure
figtemp = glucifer.Figure( figsize=(800,400) )
figtemp.append( glucifer.objects.Surface(mesh, sca.Dimensionalize(temperatureField, u.degC), colours="blue white red") )
figtemp.append( glucifer.objects.VectorArrows(mesh, velocityField/20.0, arrowHead=0.2, scaling=0.1) )
figtemp.show()