moresi-2014-material-layout.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Material Layout**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import underworld as uw\n",
    "from underworld import function as fn\n",
    "from underworld.scaling import units as u\n",
    "from underworld.scaling import dimensionalise, non_dimensionalise\n",
    "from underworld.scaling import non_dimensionalise as nd\n",
    "\n",
    "import math\n",
    "import glucifer\n",
    "import numpy as np\n",
    "import os\n",
    "import scipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import model_properties as modprop\n",
    "import material_properties as matprop\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "%run moresi-2014-scaling.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# I assume here the origin is a the top, front, middle\n",
    "# 'middle' being the slab hinge at top, front\n",
    "pert = 0.3  # nondimensional pert \n",
    "slab_xStart = 2500. * u.kilometer\n",
    "slab_dx = 3000.0 * u.kilometer  # was 7000 km in Moresi 2014\n",
    "slab_dy =  100.0 * u.kilometer\n",
    "slab_dz = 3000.0 * u.kilometer # this is the entire domain width\n",
    "slab_layers = 4\n",
    "\n",
    "slab_crust = 7.0 * u.kilometer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "backarc_dx = 1200. * u.kilometer\n",
    "backarc_dy =  100. * u.kilometer\n",
    "backarc_xStart = slab_xStart - backarc_dx\n",
    "backarc_layers = 2\n",
    "\n",
    "trans_dx =  350. * u.kilometer\n",
    "trans_dy =  100. * u.kilometer\n",
    "trans_xStart = slab_xStart - backarc_dx - trans_dx\n",
    "trans_layers = 2\n",
    "\n",
    "craton_dx = 750. * u.kilometer\n",
    "craton_dy = 150. * u.kilometer\n",
    "craton_xStart = slab_xStart - backarc_dx - trans_dx - craton_dx\n",
    "craton_layers = 2\n",
    "\n",
    "ribbon_dx =  500. * u.kilometer\n",
    "ribbon_dy =   50. * u.kilometer\n",
    "ribbon_dz = 1500. * u.kilometer \n",
    "ribbon_xStart = slab_xStart + 500. * u.kilometer\n",
    "\n",
    "bouyStrip_dx = 500. * u.kilometer\n",
    "bouyStrip_dy =  50. * u.kilometer\n",
    "bouyStrip_xStart = slab_xStart + slab_dx - bouyStrip_dx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def slabGeo(x, y, dx, dy,):\n",
    "    slabShape = np.array([ (x,y), (x+dx,y), (x+dx,y-dy), (x,y-dy), (x-pert,y-dy-pert), (x-pert,y-pert) ])\n",
    "    return fn.shape.Polygon(slabShape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def backArcGeo(x, y, dx, dy,):\n",
    "    backArcShape = np.array([ (x,y), (x+dx,y), (x+dx-pert/2.,y-dy), (x,y-dy)])\n",
    "    return fn.shape.Polygon(backArcShape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def boxGeo(x, y, dx, dy):\n",
    "    boxShape = np.array([(x,y), (x+dx,y), (x+dx,y-dy), (x,y-dy)]) \n",
    "    return (fn.shape.Polygon(boxShape))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "top50km = fn.input()[1] > nd(-50.*u.kilometer)\n",
    "fn_top50kmGeo = fn.branching.conditional( ( (top50km,  matprop.subplate1['index']),\n",
    "                                            (   True,  matprop.subplate2['index']),  \n",
    "                                          )\n",
    "                                        )\n",
    "\n",
    "outside_ribbon = fn.input()[2] < nd(ribbon_dz)\n",
    "fn_ribbonGeo = fn.branching.conditional( ( (outside_ribbon, fn_top50kmGeo  ),\n",
    "                                           (          True, matprop.ribbon['index']), \n",
    "                                         )\n",
    "                                       )  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set materials\n",
    "conditions_2d = [ \n",
    "               ( fn.input()[1] < nd( -600.0 * 10**3 * u.meter ), matprop.lm['index']) ,\n",
    "\n",
    "                # define ribbon gemetry - this is not 2D/3D auto safe\n",
    "               ( boxGeo(nd(ribbon_xStart), 0., nd(ribbon_dx), nd(ribbon_dy)), matprop.ribbon['index']),\n",
    "\n",
    "                # define buoyant strip gemetry\n",
    "               ( boxGeo(nd(bouyStrip_xStart), 0., nd(bouyStrip_dx), nd(bouyStrip_dy)), matprop.buoyStrip['index']),\n",
    "\n",
    "                # define slab geometery\n",
    "               ( slabGeo(nd(slab_xStart), -0.*nd(slab_crust)         , nd(slab_dx), nd(slab_dy)/slab_layers), matprop.subplate1['index']),\n",
    "               ( slabGeo(nd(slab_xStart), -1.*nd(slab_dy)/slab_layers, nd(slab_dx), nd(slab_dy)/slab_layers), matprop.subplate2['index']),\n",
    "               ( slabGeo(nd(slab_xStart), -2.*nd(slab_dy)/slab_layers, nd(slab_dx), nd(slab_dy)/slab_layers), matprop.subplate3['index']),\n",
    "               ( slabGeo(nd(slab_xStart), -3.*nd(slab_dy)/slab_layers, nd(slab_dx), nd(slab_dy)/slab_layers), matprop.subplate4['index']),\n",
    "                # dedine back arc geometry\n",
    "               ( backArcGeo(nd(backarc_xStart),  -0.*nd(backarc_dy)/backarc_layers, nd(backarc_dx), nd(backarc_dy)), matprop.backArc1['index']),\n",
    "               ( backArcGeo(nd(backarc_xStart),  -1.*nd(backarc_dy)/backarc_layers, nd(backarc_dx), nd(backarc_dy)), matprop.backArc2['index']),\n",
    "                # define transition gemetry\n",
    "               ( boxGeo(nd(trans_xStart),  -0.*nd(trans_dy)/trans_layers, nd(trans_dx), nd(trans_dy)), matprop.trans1['index']),\n",
    "               ( boxGeo(nd(trans_xStart),  -1.*nd(trans_dy)/trans_layers, nd(trans_dx), nd(trans_dy)), matprop.trans2['index']),\n",
    "                # define craton gemetry\n",
    "               ( boxGeo(nd(craton_xStart), -0.*nd(craton_dy)/craton_layers, nd(craton_dx), nd(craton_dy)), matprop.craton1['index']),\n",
    "               ( boxGeo(nd(craton_xStart), -1.*nd(craton_dy)/craton_layers, nd(craton_dx), nd(craton_dy)), matprop.craton2['index']),\n",
    "                # otherwise upper mantle!\n",
    "                ( True, matprop.um['index']),\n",
    "             ] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "conditions_3d = [ \n",
    "               ( fn.input()[1] < nd( -600.0 * 10**3 * u.meter ), matprop.lm['index']) ,\n",
    "\n",
    "                # define ribbon gemetry - this is not 2D/3D auto safe\n",
    "               ( boxGeo(nd(ribbon_xStart), 0., nd(ribbon_dx), nd(ribbon_dy)), fn_ribbonGeo),\n",
    "\n",
    "                # define buoyant strip gemetry\n",
    "               ( boxGeo(nd(bouyStrip_xStart), 0., nd(bouyStrip_dx), nd(bouyStrip_dy)), matprop.buoyStrip['index']),\n",
    "\n",
    "                # define slab geometery\n",
    "               ( slabGeo(nd(slab_xStart), -0.*nd(slab_crust)         , nd(slab_dx), nd(slab_dy)/slab_layers), matprop.subplate1['index']),\n",
    "               ( slabGeo(nd(slab_xStart), -1.*nd(slab_dy)/slab_layers, nd(slab_dx), nd(slab_dy)/slab_layers), matprop.subplate2['index']),\n",
    "               ( slabGeo(nd(slab_xStart), -2.*nd(slab_dy)/slab_layers, nd(slab_dx), nd(slab_dy)/slab_layers), matprop.subplate3['index']),\n",
    "               ( slabGeo(nd(slab_xStart), -3.*nd(slab_dy)/slab_layers, nd(slab_dx), nd(slab_dy)/slab_layers), matprop.subplate4['index']),\n",
    "                # dedine back arc geometry\n",
    "               ( backArcGeo(nd(backarc_xStart),  -0.*nd(backarc_dy)/backarc_layers, nd(backarc_dx), nd(backarc_dy)), matprop.backArc1['index']),\n",
    "               ( backArcGeo(nd(backarc_xStart),  -1.*nd(backarc_dy)/backarc_layers, nd(backarc_dx), nd(backarc_dy)), matprop.backArc2['index']),\n",
    "                # define transition gemetry\n",
    "               ( boxGeo(nd(trans_xStart),  -0.*nd(trans_dy)/trans_layers, nd(trans_dx), nd(trans_dy)), matprop.trans1['index']),\n",
    "               ( boxGeo(nd(trans_xStart),  -1.*nd(trans_dy)/trans_layers, nd(trans_dx), nd(trans_dy)), matprop.trans2['index']),\n",
    "                # define craton gemetry\n",
    "               ( boxGeo(nd(craton_xStart), -0.*nd(craton_dy)/craton_layers, nd(craton_dx), nd(craton_dy)), matprop.craton1['index']),\n",
    "               ( boxGeo(nd(craton_xStart), -1.*nd(craton_dy)/craton_layers, nd(craton_dx), nd(craton_dy)), matprop.craton2['index']),\n",
    "                # otherwise upper mantle!\n",
    "                ( True, matprop.um['index']),\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
}