This model uses a pressure boundary condition to force an uplift.
This model also utilises scaling our numbers into dimensionless units.
import numpy as np
import underworld as uw
import math
from underworld import function as fn
import glucifer
import os
import unsupported.geodynamics.scaling as sca
u = sca.UnitRegistry
# reference units
KL_meters = 100e3 * u.meter
K_viscosity = (1e16 * u.pascal * u.second)
K_density = (3.3e3 * u.kilogram / (u.meter)**3 )
KM_kilograms = K_density * KL_meters**3
KT_seconds = KM_kilograms / ( KL_meters * K_viscosity )
K_substance = 1. * u.mole
Kt_degrees = 1. * u.kelvin
sca.scaling["[length]"] = KL_meters.to_base_units()
sca.scaling["[temperature]"] = Kt_degrees.to_base_units()
sca.scaling["[time]"] = KT_seconds.to_base_units()
sca.scaling["[mass]"] = KM_kilograms.to_base_units()
# all nondimensional units
gravity = sca.nonDimensionalize(9.81 * u.meter / u.second**2)
density = sca.nonDimensionalize( 3300 * u.kilogram / u.meter**3)
viscosity = sca.nonDimensionalize( 1e22 * u.Pa * u.sec)
bulk_visc = sca.nonDimensionalize( 1e11 * u.Pa *u.sec)
Lx = sca.nonDimensionalize( 100e3 * u.meter)
Ly = sca.nonDimensionalize( 60e3 * u.meter)
dx = sca.nonDimensionalize( 5e3 * u.meter)
dy = sca.nonDimensionalize( 5e3 * u.meter)
center = sca.nonDimensionalize(50e3 * u.meter)
width = sca.nonDimensionalize(3e3*u.meter)
lithostaticPressure = 0.6*Ly*density*gravity
resUnit = 5
boxLength = Lx
boxHeight = Ly
elType = "Q1/dQ0"
resx = 100
resy = 60
minCoord = [0.,0.]
maxCoord = [boxLength,boxHeight]
mesh = uw.mesh.FeMesh_Cartesian( elementType = (elType),
elementRes = (resx, resy),
minCoord = minCoord,
maxCoord = maxCoord )
velocityField = mesh.add_variable( nodeDofCount=mesh.dim )
tractionField = mesh.add_variable( nodeDofCount=2 )
pressureField = mesh.subMesh.add_variable( nodeDofCount=1 )
velocityField.data[:] = [0.,0.]
pressureField.data[:] = 0.
tractionField.data[:] = [0.0,0.0]
# Traction is force per unit area associated with a specific surface
# ie, traction = stress * surface_unit_normal
for ii in mesh.specialSets['MinJ_VertexSet']:
coord = mesh.data[ii]
tractionField.data[ii] = [0.0,lithostaticPressure*(1.+0.2*np.exp((-1/width*(coord[0]-center)**2)))]
# visualise the bottom stress condition
if uw.rank() == 0:
uw.matplotlib_inline()
import matplotlib.pyplot as pyplot
import matplotlib.pylab as pylab
pyplot.ion()
pylab.rcParams[ 'figure.figsize'] = 12, 6
xcoord = mesh.data[mesh.specialSets['MinJ_VertexSet'].data][:,0]
stress = tractionField.data[mesh.specialSets['MinJ_VertexSet'].data][:,1]
pyplot.plot( xcoord, stress, 'o', color = 'black', label='numerical')
pyplot.xlabel('Y coords')
pyplot.ylabel('Temperature')
pyplot.show()
# Initialise a swarm.
swarm = uw.swarm.Swarm( mesh=mesh, particleEscape=True )
advector= uw.systems.SwarmAdvector(velocityField, swarm, order=2)
# Add a data variable which will store an index to determine material.
materialVariable = swarm.add_variable( dataType="double", count=1 )
# Create a layout object that will populate the swarm across the whole domain.
swarmLayout = uw.swarm.layouts.PerCellSpaceFillerLayout( swarm=swarm, particlesPerCell=20 )
# Populate.
swarm.populate_using_layout( layout=swarmLayout )
# material 0 - compressible Lambda=10, density = 0
# material 1 - incompressible Lambda=0, density = 1
materialVariable.data[:]=0
for index,coord in enumerate(swarm.particleCoordinates.data):
if coord[1] < boxHeight*0.6:
materialVariable.data[index]=1
# population control regulars particle creation and deletion
# important for inflow/outflow problems
population_control = uw.swarm.PopulationControl(swarm,
aggressive=True,splitThreshold=0.15, maxDeletions=2,maxSplits=10,
particlesPerCell=20)
# build tracer swarm for fluid level - only 1 particle
mswarm = uw.swarm.Swarm( mesh=mesh, particleEscape=True )
msAdvector= uw.systems.SwarmAdvector(velocityField, mswarm, order=2)
# initial height at 'air' level
particleCoordinates = np.zeros((1,2))
particleCoordinates[:,0] = 0.5*Lx
particleCoordinates[:,1] = 0.6*Ly
ignore=mswarm.add_particles_with_coordinates(particleCoordinates)
fig1 = glucifer.Figure(rulers=True, boundingBox=((0.0, 0.0), (Lx, Ly)))
fig1.append( glucifer.objects.Points(mswarm, colourBar=False ) )
fig1.append( glucifer.objects.Points(swarm, materialVariable, fn_size=2. ) )
fig1.show()
# lambdaFn is created for pseudo compressible air layer
lambdaFn = uw.function.branching.map( fn_key=materialVariable,
mapping={ 0: 1/bulk_visc, 1: 0.0 } )
densityFn = uw.function.branching.map( fn_key=materialVariable,
mapping={ 0: 0.0, 1: density } )
forceFn = densityFn * (0.0,-gravity)
iWalls = mesh.specialSets["MinI_VertexSet"] + mesh.specialSets["MaxI_VertexSet"]
jWalls = mesh.specialSets["MinJ_VertexSet"] + mesh.specialSets["MaxJ_VertexSet"]
bottomWall = mesh.specialSets["MinJ_VertexSet"]
allWalls = iWalls + jWalls
# Now, using these sets, decide which degrees of freedom (on each node) should be considered Dirichlet.
stokesBC = uw.conditions.DirichletCondition( variable = velocityField,
indexSetsPerDof = (iWalls, jWalls-bottomWall) )
# add neumann bcs
nbc = uw.conditions.NeumannCondition( fn_flux=tractionField, variable = velocityField,
indexSetsPerDof = (None, bottomWall ) )
# setup solver
stokesPIC = uw.systems.Stokes( velocityField = velocityField,
pressureField = pressureField,
conditions = [stokesBC, nbc],
fn_viscosity = viscosity,
fn_bodyforce = forceFn,
fn_one_on_lambda = lambdaFn )
solver = uw.systems.Solver( stokesPIC )
# setup analysis function
vdotv = fn.math.dot(velocityField,velocityField)
v2sum_integral = uw.utils.Integral( mesh=mesh, fn=vdotv )
volume_integral = uw.utils.Integral( mesh=mesh, fn=1. )
velmag = fn.math.sqrt(vdotv)
steps = 0
finalTimestep = 3
fieldDict = {'velocity':velocityField, 'pressure':pressureField}
swarmDict = {'material':materialVariable}
prefix='uplift/'
if prefix is None:
prefix=''
else:
try:
import os
if not os.path.exists("./"+prefix+"/"):
os.makedirs("./"+prefix+'/')
except:
raise
outfile = open(prefix+'buildMount.txt', 'w+')
string = "steps, timestep, vrms, peak height"
print(string)
outfile.write( string+"\n")
# initialise loop
dt = -1
h1 = mswarm.particleCoordinates.data[:,1].copy()
while steps<finalTimestep:
# Get solution
solver.solve()
# calculate metrics
v2int = v2sum_integral.evaluate()[0]
vol = volume_integral.evaluate()[0]
# get time step first time around
if dt < 0:
dt = advector.get_max_dt()
h0 = h1.copy() # NOTE the copy()
# Advect particles
advector.integrate(dt)
msAdvector.integrate(dt)
# update peak heigh
h1 = mswarm.particleCoordinates.data[:,1]
diffH = h1-h0
string = "{}, {}, {}, {}".format(steps,
sca.Dimensionalize(dt, u.year),
sca.Dimensionalize(np.sqrt(v2int/vol), u.cm/u.year),
sca.Dimensionalize(diffH, u.metre)[0] )
print(string)
outfile.write(string+"\n")
# population control
population_control.repopulate()
fig1.save(prefix+"particals-"+str(steps)+".png")
steps += 1
outfile.close()
cm_per_year = sca.Dimensionalize(1,u.centimeter/u.year)
# note, scaling by the 'Quantity' object is slow, use the magnitude insead
type(cm_per_year.magnitude)
fig2 = glucifer.Figure()
fig2.append( glucifer.objects.Surface(mesh, cm_per_year.magnitude*velmag) )
fig2.show()
fig1.append( glucifer.objects.VectorArrows(mesh, cm_per_year.magnitude*0.1*velocityField) )
fig1.show()
# for testing purposes
dimensionalise=(diffH*sca.Dimensionalize(1,u.meter))
if np.fabs(dimensionalise.magnitude-245.140) > 0.05*245.140:
raise RuntimeError("Height of passive tracer outside expected 5% tolerance")