This example solves 2D dimensionless isoviscous thermal convection with a Rayleigh number of $10^4$, see Blankenbach et al. 1989 for details.
This example introduces:
Keywords: material parameters, Stokes system, advective diffusive systems
References
B. Blankenbach, F. Busse, U. Christensen, L. Cserepes, D. Gunkel, U. Hansen, H. Harder, G. Jarvis, M. Koch, G. Marquart, D. Moore, P. Olson, H. Schmeling and T. Schnaubelt. A benchmark comparison for mantle convection codes. Geophysical Journal International, 98, 1, 23–38, 1989 http://onlinelibrary.wiley.com/doi/10.1111/j.1365-246X.1989.tb05511.x/abstract
import underworld as uw
from underworld import function as fn
import glucifer
import math
# Set simulation box size.
boxHeight = 1.0
boxLength = 2.0
# Set the resolution.
res = 16
# Set min/max temperatures.
tempMin = 0.0
tempMax = 1.0
The mesh object has both a primary and sub mesh. "Q1/dQ0" produces a primary mesh with element type Q1 and a sub-mesh with elements type dQ0. Q1 elements have nodes at the element corners, dQ0 elements have a single node at the elements centre.
mesh = uw.mesh.FeMesh_Cartesian( elementType = ("Q1/dQ0"),
elementRes = (2*res, res),
minCoord = (0., 0.),
maxCoord = (boxLength, boxHeight))
Create mesh variables. Note the pressure field uses the sub-mesh.
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[:] = [0.,0.]
pressureField.data[:] = 0.
temperatureDotField.data[:] = 0.
Set functions for viscosity, density and buoyancy force. These functions and variables only need to be defined at the beginning of the simulation, not each timestep.
# Set viscosity to be a constant.
viscosity = 1.
# 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, 1.0 )
# 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, temperatureField, 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 = 1.0,
conditions = tempBC )
# 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)
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, temperatureField, colours="blue white red") )
figtemp.append( glucifer.objects.VectorArrows(mesh, velocityField/100.0, arrowHead=0.2, scaling=0.1) )
figtemp.show()