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