Rayleigh-Taylor instability¶
This notebook models the Rayleigh-Taylor instability outlined in van Keken et al. (1997).
The system of equations is given by
$$ \nabla \cdot \left( \eta \nabla \dot\varepsilon \right) - \nabla p = -\rho g \mathbf{\hat z} $$$$ \nabla \cdot \mathbf{v} = 0 $$Keywords: particle swarms, Stokes system, advective diffusive systems
References
- van Keken, P.E., S.D. King, H. Schmeling, U.R. Christensen, D.Neumeister and M.-P. Doin. A comparison of methods for the modeling of thermochemical convection. Journal of Geophysical Research, 102, 22,477-22,495, 1997.
http://onlinelibrary.wiley.com/doi/10.1029/97JB01353/abstract
In [1]:
import underworld as uw
from underworld import function as fn
import glucifer
import math
import numpy as np
Setup parameters¶
In [2]:
res = 64
boxLength = 0.9142
boxHeight = 1.0
# light material viscosity / dense material viscosity
viscosityRatio = 1.0
In [3]:
inputPath = 'input/1_06_Rayleigh_Taylor/'
outputPath = 'output/'
# Make output directory if necessary.
if uw.rank()==0:
import os
if not os.path.exists(outputPath):
os.makedirs(outputPath)
Create mesh and finite element variables¶
In [4]:
mesh = uw.mesh.FeMesh_Cartesian( elementType = ("Q1/dQ0"),
elementRes = (res, res),
minCoord = (0., 0.),
maxCoord = (boxLength, boxHeight))
velocityField = mesh.add_variable( nodeDofCount=2 )
pressureField = mesh.subMesh.add_variable( nodeDofCount=1 )
# initialise
velocityField.data[:] = [0.,0.]
pressureField.data[:] = 0.
Create a particle swarm¶
In [5]:
# Create a swarm.
swarm = uw.swarm.Swarm( mesh=mesh )
# Create a data variable. It will be used to store the material index of each particle.
materialIndex = swarm.add_variable( dataType="int", count=1 )
# Create a layout object, populate the swarm with particles.
swarmLayout = uw.swarm.layouts.GlobalSpaceFillerLayout( swarm=swarm, particlesPerCell=20 )
swarm.populate_using_layout( layout=swarmLayout )
Initialise each particle's material index¶
In [6]:
# define these for convience.
denseIndex = 0
lightIndex = 1
# material perturbation from van Keken et al. 1997
wavelength = 2.0*boxLength
amplitude = 0.02
offset = 0.2
k = 2. * math.pi / wavelength
# Create function to return particle's coordinate
coord = fn.coord()
# Define the material perturbation, a function of the x coordinate (accessed by `coord[0]`).
perturbationFn = offset + amplitude*fn.math.cos( k*coord[0] )
# Setup the conditions list.
# If z is less than the perturbation, set to lightIndex.
conditions = [ ( perturbationFn > coord[1] , lightIndex ),
( True , denseIndex ) ]
# The swarm is passed as an argument to the evaluation, providing evaluation on each particle.
# Results are written to the materialIndex swarm variable.
materialIndex.data[:] = fn.branching.conditional( conditions ).evaluate(swarm)
Plot the particles by material
In [7]:
fig1 = glucifer.Figure()
fig1.append( glucifer.objects.Points(swarm, materialIndex, pointSize=2, colourBar=False) )
fig1.show()
Map properties to material index¶
The Map function allows us to create 'per material' type behaviour. Again we use the branching function to set up a (condition, action) command.
In [8]:
# Set a density of '0.' for light material, '1.' for dense material.
densityMap = { lightIndex:0., denseIndex:1. }
densityFn = fn.branching.map( fn_key = materialIndex, mapping = densityMap )
# Set a viscosity value of '1.' for both materials.
viscosityMap = { lightIndex:viscosityRatio, denseIndex:1. }
fn_viscosity = fn.branching.map( fn_key = materialIndex, mapping = viscosityMap )
# Define a vertical unit vector using a python tuple.
z_hat = ( 0.0, 1.0 )
# Create buoyancy force vector
buoyancyFn = -densityFn*z_hat
Boundary conditions¶
Create free-slip condition on the vertical boundaries, and a no-slip condition on the horizontal boundaries.
In [9]:
# Construct node sets using the mesh specialSets
iWalls = mesh.specialSets["MinI_VertexSet"] + mesh.specialSets["MaxI_VertexSet"]
jWalls = mesh.specialSets["MinJ_VertexSet"] + mesh.specialSets["MaxJ_VertexSet"]
allWalls = iWalls + jWalls
# Prescribe degrees of freedom on each node to be considered Dirichlet conditions.
# In the x direction on allWalls flag as Dirichlet
# In the y direction on jWalls (horizontal) flag as Dirichlet
stokesBC = uw.conditions.DirichletCondition( variable = velocityField,
indexSetsPerDof = (allWalls, jWalls) )
Create systems¶
In [10]:
stokes = uw.systems.Stokes( velocityField = velocityField,
pressureField = pressureField,
voronoi_swarm = swarm,
conditions = stokesBC,
fn_viscosity = fn_viscosity,
fn_bodyforce = buoyancyFn,
fn_source=0.)
solver = uw.systems.Solver( stokes )
# Optional solver settings
if(uw.nProcs()==1):
solver.set_inner_method("lu")
# solver.set_penalty(1.0e6)
solver.set_penalty(0.)
solver.options.scr.ksp_rtol = 1.0e-3
# Create a system to advect the swarm
advector = uw.systems.SwarmAdvector( swarm=swarm, velocityField=velocityField, order=2 )
In [11]:
if 0.!=None:
print("hi")
In [24]:
stokes.fn_source = 0.
# stokes._cforceVecTerm.fn = 0.
# stokes._fn_source
Time stepping¶
In [25]:
# Initialise time and timestep.
time = 0.
step = 0
# We set timeEnd=3 so that the simulation completes quickly,
# but generally you will want to set timeEnd~300 to capture
# the peak V_rms.
timeEnd = 3.
outputEvery = 20
timeVal = []
vrmsVal = []
# functions for calculating RMS velocity
vdotv = fn.math.dot(velocityField,velocityField)
# Save mesh and retain file handle for future xdmf creation
meshFileHandle = mesh.save(outputPath+"Mesh.h5")
In [26]:
# define an update function
def update():
# Retrieve the maximum possible timestep for the advection system.
dt = advector.get_max_dt()
# Advect using this timestep size.
advector.integrate(dt)
return time+dt, step+1
In [27]:
while time < timeEnd:
# Get instantaneous Stokes solution
solver.solve()
# Calculate the RMS velocity.
vrms = math.sqrt( mesh.integrate(vdotv)[0] / mesh.integrate(1.)[0] )
# Record values into arrays
if(uw.rank()==0):
vrmsVal.append(vrms)
timeVal.append(time)
# Output to disk
if step%outputEvery == 0:
if(uw.rank()==0):
print 'step = {0:6d}; time = {1:.3e}; v_rms = {2:.3e}'.format(step,time,vrms)
velFilename = outputPath+"/velocityField."+str(step)
vFH = velocityField.save(velFilename+".h5")
velocityField.xdmf( velFilename, vFH, "velocity", meshFileHandle, "Mesh", time )
pressureFilename = outputPath+"/pressureField."+str(step)
pFH = pressureField.save(pressureFilename+".h5")
pressureField.xdmf(pressureFilename, pFH, "pressure", meshFileHandle, "Mesh", time )
outputFilename = outputPath+"image"+str(step).zfill(4)
fig1.save_image(outputFilename)
# We are finished with current timestep, update.
time, step = update()
if(uw.rank()==0):
print 'step = {0:6d}; time = {1:.3e}; v_rms = {2:.3e}'.format(step,time,vrms)
Post simulation analysis¶
In [23]:
fig1.append( glucifer.objects.VectorArrows( mesh, velocityField, scaling=1.0e2))
fig1.show()
In [17]:
if(uw.rank()==0):
maxIndex = np.argmax(vrmsVal)
print('Viscosity ratio = {0:.2f}'.format(viscosityRatio))
print(' t(max vrms) = {0:.2f}'.format(timeVal[maxIndex]))
print(' vrms = {0:.4f}'.format(vrmsVal[maxIndex]))
# output a summary file with benchmark values (useful for parallel runs)
np.savetxt(outputPath+'summary.txt', [viscosityRatio, timeVal[maxIndex], vrmsVal[maxIndex]])
The benchmark values from van Keken et al. 1997 are approximately:
| Viscosity ratio ($\frac{\eta}{\eta_r}$) | t (max $v_{rms}$) | max $v_{rms}$ |
|---|---|---|
| 1.00 | 208.99 | 0.0031 |
| 0.10 | 71.92 | 0.0095 |
| 0.01 | 49.57 | 0.0146 |
Plot RMS velocity
In [18]:
if uw.nProcs() == 1:
if viscosityRatio == 1.0 :
data = np.loadtxt(inputPath+'VrmsCaseA.txt', unpack=True )
elif viscosityRatio == 0.1 :
data = np.loadtxt(inputPath+'VrmsCaseB.txt', unpack=True )
elif viscosityRatio == 0.01 :
data = np.loadtxt(inputPath+'VrmsCaseC.txt', unpack=True )
else :
print('No specific data found - default to Case A')
data = np.loadtxt(inputPath+'VrmsCaseA.txt', unpack=True )
# Load into data arrays to compare with timevals and vrmsvals from above.
timeCompare, vrmsCompare = data[0], data[1]
# These can be copied onto timevals and vrmsvals to skip running the time loop.
uw.matplotlib_inline()
import matplotlib.pyplot as pyplot
fig = pyplot.figure()
fig.set_size_inches(12, 6)
ax = fig.add_subplot(1,1,1)
ax.plot(timeCompare, vrmsCompare, color = 'black')
ax.plot(timeVal, vrmsVal, color = 'red', marker=".", markersize=10)
ax.set_xlabel('Time')
ax.set_ylabel('RMS velocity')
ax.set_xlim([0.0,1000.0])
fig.show()
