3D Stokes sinker for UWGeo - 2.10.2

3DStokesSinker.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stokes Sinker\n",
    "\n",
    "Demonstration example for setting up particle swarms with different material properties. This system consists of a dense, high viscosity sphere falling through a background lower density and viscosity fluid.\n",
    "\n",
    "![Stokes 2D](./images/Stokes2D.gif)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "loaded rc file /opt/venv/lib/python3.7/site-packages/UWGeodynamics/uwgeo-data/uwgeodynamicsrc\n"
     ]
    }
   ],
   "source": [
    "import UWGeodynamics as GEO\n",
    "from UWGeodynamics import visualisation as vis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "u = GEO.UnitRegistry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "velocity = 1.0 * u.centimeter / u.hour\n",
    "model_length = 2. * u.meter\n",
    "model_height = 1. * u.meter\n",
    "refViscosity = 1e6 * u.pascal * u.second\n",
    "bodyforce = 200 * u.kilogram / u.metre**3 * 9.81 * u.meter / u.second**2\n",
    "\n",
    "KL = model_height\n",
    "Kt = KL / velocity\n",
    "KM = bodyforce * KL**2 * Kt**2\n",
    "\n",
    "GEO.scaling_coefficients[\"[length]\"] = KL\n",
    "GEO.scaling_coefficients[\"[time]\"] = Kt\n",
    "GEO.scaling_coefficients[\"[mass]\"]= KM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Initialise the Model\n",
    "\n",
    "This will set up a 2 meters x 1 meter box, the resolution is 64 x 64."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "Model = GEO.Model(elementRes=(16, 16, 16), \n",
    "                  minCoord=(-1. * u.meter, -50. * u.centimeter, -1. * u.meter), \n",
    "                  maxCoord=(1. * u.meter, 50. * u.centimeter, 1. * u.meter),\n",
    "                  gravity = (0.0, -9.8 * u.meter / u.second**2, 0.))\n",
    "\n",
    "Model.outputDir= \"3DStokesSinker\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define the Materials\n",
    "\n",
    "We define a `heavyMaterial` which will represent the background medium in which the ball will fall.\n",
    "The Ball itself is defined using a `lightMaterial` with an initial disk shape."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "lightMaterial = Model.add_material(name=\"Light\", shape=GEO.shapes.Layer3D(top=Model.top, bottom=Model.bottom))\n",
    "heavyMaterial = Model.add_material(name=\"Heavy\", shape=GEO.shapes.Sphere(center=(0.,30.*u.centimetre,0.), radius=10. * u.centimetre))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Material properties\n",
    "\n",
    "The materials have the same viscosity but their density differs, the `heavyMaterial` is 50 times heavier than the\n",
    "surrounding materials."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "lightMaterial.density = 10 * u.kilogram / u.metre**3\n",
    "heavyMaterial.density = 500 * u.kilogram / u.metre**3\n",
    "\n",
    "lightMaterial.viscosity = 1e6 * u.pascal * u.second\n",
    "heavyMaterial.viscosity = 1e6 * u.pascal * u.second"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define Boundary Conditions\n",
    "\n",
    "The boundary conditions are freeslip everywhere (zero shear stress)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<underworld.conditions._conditions.DirichletCondition at 0x7f6f5eda69e8>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Model.set_velocityBCs(left=[0,None,0], right=[0,None,0], \n",
    "                      top=[None,0,0], bottom=[None,0,0],\n",
    "                      front=[0,0,None], back=[0,0,None])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualise Initial Set up"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Analysis tools\n",
    "-----\n",
    "\n",
    "We define a set of metrics to monitor the evolution of the model through time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import underworld as uw\n",
    "import underworld.function as fn\n",
    "import math"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**RMS velocity**\n",
    "\n",
    "The root mean squared velocity is defined by intergrating over the entire simulation domain via\n",
    "\n",
    "\\\\[\n",
    "\\begin{aligned}\n",
    "v_{rms}  =  \\sqrt{ \\frac{ \\int_V (\\mathbf{v}.\\mathbf{v}) dV } {\\int_V dV} }\n",
    "\\end{aligned}\n",
    "\\\\]\n",
    "\n",
    "where $V$ denotes the volume of the box."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "vdotv = fn.math.dot(Model.velocityField, Model.velocityField)\n",
    "v2sum_integral  = uw.utils.Integral(mesh=Model.mesh, fn=vdotv )\n",
    "volume_integral = uw.utils.Integral(mesh=Model.mesh, fn=1. )\n",
    "vrms = math.sqrt(v2sum_integral.evaluate()[0] / volume_integral.evaluate()[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Position of the bottom of the Ball**\n",
    "\n",
    "We will use a passive tracers, initially located at the bottom of the disk.\n",
    "Note that because the ball is going to deform, the position of the passive tracers may change lateraly.\n",
    "In practice this is very minimal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "coords = np.ndarray((1, 3))\n",
    "coords[...] = 0.\n",
    "coords[:, 1] = GEO.nd(20. * u.centimetre)\n",
    "\n",
    "pt = Model.add_passive_tracers(name=\"tip\", vertices=coords)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the tools have been defined, we need a way to get them to execute after each time step.\n",
    "This can be done using the `Model.postSolveHook` entry point. \n",
    "\n",
    "We need to define a python function that will process the output of Model and extract information: In the following we define 2 containers in the form of python list objects. The position of the passive tracers and the vrms will be appended to their respective list after each timestep using the `post_solve_hook` python function. Note that the lists must be defined as global inside the python function so that we can retrieve them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# parallel safe way of finding the particles vertical coordinate.\n",
    "fn_y = fn.coord()[1]\n",
    "fn_y_minmax = fn.view.min_max(fn_y)\n",
    "\n",
    "tSinker = [0.]\n",
    "fn_y_minmax.reset()\n",
    "fn_y_minmax.evaluate(pt)\n",
    "ypos = fn_y_minmax.max_global()\n",
    "\n",
    "ypos = GEO.dimensionalise(ypos, u.centimetre)\n",
    "ypos = ypos.magnitude\n",
    "ySinker = [ypos]\n",
    "vrms = [math.sqrt(v2sum_integral.evaluate()[0] / volume_integral.evaluate()[0])]\n",
    "\n",
    "def post_solve_hook():\n",
    "    global tSinker\n",
    "    global ySinker\n",
    "    global vrms\n",
    "    \n",
    "    time = Model.time.to(u.hour).magnitude\n",
    "    fn_y_minmax.reset()\n",
    "    fn_y_minmax.evaluate(pt)\n",
    "    ypos = fn_y_minmax.max_global()\n",
    "    \n",
    "    ypos = GEO.dimensionalise(ypos, u.centimetre)\n",
    "    ypos = ypos.magnitude\n",
    "    \n",
    "    tSinker.append(time)\n",
    "    ySinker.append(ypos)\n",
    "    vrms.append(math.sqrt(v2sum_integral.evaluate()[0] / volume_integral.evaluate()[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `post_solve_hook` function is \"attached\" to `Model.postSolveHook`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "Model.post_solve_functions[\"Measurements\"] = post_solve_hook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "Model.init_model()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running with UWGeodynamics version 2.10.2\n",
      "Options:  -Q22_pc_type uw -ksp_type bsscr -pc_type none -ksp_k2_type NULL -rescale_equations False -remove_constant_pressure_null_space False -change_backsolve False -change_A11rhspresolve False -restore_K False -A11_ksp_type fgmres -A11_ksp_rtol 1e-06 -scr_ksp_type fgmres -scr_ksp_rtol 1e-05\n",
      "Step:     1 Model Time: 37.1 minute dt: 37.1 minute (2020-12-15 22:13:13)\n",
      "Step:     2 Model Time: 1.2 hour dt: 36.1 minute (2020-12-15 22:13:24)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Model.run_for(nstep=2, checkpoint_interval=1)#, restartStep=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualisation of the results\n",
    "\n",
    "**Material Field**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quick Analysis\n",
    "\n",
    "The position of the Sinker through time can be plotted using the `tSinker` and `ySinker` lists.\n",
    "\n",
    "Here we use Matplotlib to make the plot. Matplotlib is not parallel safe and will return message errors when attempting to run\n",
    "this Model on multiple processors. To avoid this, the user will need to run the Matplotlib function on one CPUs.\n",
    "This can be achieved using a condition on the processor `rank`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dfZyNdf7H8ddnhhEiZPKLtDZbWt1RQyKJLLqjkFXRnchN97tabW27bVtbqW1rQwnRbomNUBGKUnLTKIWSVCpUVMTWlJv5/P64LrvTdGbmGOec68zM+/l4nMecc53rOud9ZfTxvb7X9/s1d0dERKSwjKgDiIhIelKBEBGRmFQgREQkJhUIERGJSQVCRERiqhR1gESqW7euN2rUKOoYIiJlxrJly7509+xY75WrAtGoUSNyc3OjjiEiUmaY2cdFvadLTCIiEpMKhIiIxKQCISIiMalAiIhITCoQIiISU9IKhJk1NLP5Zvauma0ys2vC7XXMbK6ZvR/+rF3E8ReH+7xvZhcnK+e0NzfQ5s55/HzYc7S5cx7T3tyQrK8SESlTktmC2AX8xt1/CbQChphZU2AY8KK7Hw68GL7+ETOrA/wROBFoCfyxqEKyL6a9uYEbp65gw9Y8HNiwNY8bp65QkRARIYkFwt0/c/c3wufbgXeBBkA3YEK42wTgnBiHdwbmuvvX7r4FmAt0SXTG4bPfI2/nbq7KnMqx9gEAeTt3M3z2e4n+KhGRMiclfRBm1ghoDiwB6rn7ZxAUEeCgGIc0AD4t8Hp9uC3WZw8ws1wzy928efNe5dq4NY+a/IcLKs3j6axbGFZpIlXYwcateXv1OSIi5VHSC4SZ7Q9MAa51923xHhZjW8yVjdx9tLvnuHtOdnbM0eJFql+rKtvYn04/3M2k3acysNIzzMy6kc411+3V54iIlEdJLRBmVpmgODzu7lPDzV+Y2cHh+wcDm2Icuh5oWOD1IcDGROcb2rkJVStnsp1q/H5Xfy7ccSP72S5G7bgJZg2DHd8m+itFRMqMZN7FZMBY4F13/1uBt2YAe+5KuhiYHuPw2UAnM6sddk53Crcl1DnNG/DX7sfQoFZVDFhXsyVvnDUTa9kfloyCkSfBRwsS/bUiImWCJWtNajM7GXgFWAHkh5t/T9APMRk4FPgEOM/dvzazHGCgu18eHn9ZuD/A7e7+aEnfmZOT4wmbrG/dQphxJXz9IeRcBh1vhf1qJuazRUTShJktc/ecmO8lq0BEIaEFAmDHdzD/dlg8EmrUh7Pvh8M7Ju7zRUQiVlyB0Ejq4mRVg863w2VzIKs6PN4Dpg2GvC1RJxMRSToViHg0bAFXLIC2v4G3noQRrWD1zKhTiYgklQpEvCrvB6fdAv3nQfW68OT58FQ/+ParqJOJiCSFCsTeqt8M+s+H9jfBO9NhREtYORXKUV+OiAioQJROpSxodwNc8TLUaghPXQqT+sD2L6JOJiKSMCoQ+6LeUdDvheAW2PfnBq2Jt55Ua0JEygUViH2VWQlOvhYGLYTsJvD0FfBEL/hGM8KKSNmmApEodQ+HS2dBlzvho1dgZCtYNkGtCREps1QgEikjE1oNgsGvwcHHwTNXw2PdYMvHUScTEdlrKhDJUOcwuGgGnHUfbHgjmNNpyWjIzy/5WBGRNKECkSwZGcEcToMXwaGtYNZQGH8mfPVB1MlEROKiApFstRpCnynQbSRsWgWjWsNr/4D83VEnExEplgpEKphB8wth8BJo3AHm3AxjO8Gm1VEnExEpkgpEKtU8GHo/AT3GBtOIP9wWFtwDu3dGnUxE5CdUIFLNDI7pCUOWQpMzYN5t8EgH+OztqJOJiPxIMleUG2dmm8xsZYFtx5nZIjNbYWbPmFnMFXjMbF24z3IzS+ACD2lk/2zoNQF6PQbbP4dH2sO822HXjqiTiYgAyW1BjAe6FNo2Bhjm7scATwNDizm+vbs3K2ohi3KjaTcYsgSO7gkL7obR7WDDsqhTiYgkr0C4+wLg60KbmwB7FnmeC/RI1veXKdXqQPeH4YLJkLcVxnSEubfAzryok4lIBZbqPoiVQNfw+XlAwyL2c2COmS0zswHFfaCZDTCzXDPL3bx5cwKjRuCIzjBkMTTvCwvvh4dOhk8WR51KRCqoVBeIy4AhZrYMqAEUdcG9jbsfD5we7n9KUR/o7qPdPcfdc7KzsxOfONX2OwC6PgB9pwX9EeO6wKxhsOPbqJOJSAWT0gLh7qvdvZO7nwBMBGIOK3b3jeHPTQR9FS1TlzJNNG4fjMJu2R+WjAqm6/hoQcnHiYgkSEoLhJkdFP7MAG4GHoqxT3Uzq7HnOdCJ4NJUxVNlfzhjOFwyM5gIcMLZ8Ox18P22qJOJSAWQzNtcJwKLgCZmtt7M+gHnm9kaYDWwEXg03Le+mc0MD60HvGpmbwFLgefc/flk5SwTGrWBgQvhpCth2figNfH+C1GnEpFyzrwcrVeQk5Pjubnlc9jEf336OkwfAl++B80uhM63Q9XaUacSkTLKzJYVNZxAI6nLmoYtYOAr0PY3wfKmI1rB6pklHycispdUIMqiSlXgtFug/zyoXheePB+e6gfffhV1MhEpR1QgyrL6zaD/fGh/E7wzHUa0hJVTtcypiCSECkRZVykL2t0AV7wcrD3x1KUwqQ9s/yLqZCJSxqlAlBf1joJ+L0DHW+H9uUFr4q0n1ZoQkVJTgShPMivBydfCoIWQ3QSevgKe6AXfbIg6mYiUQSoQ5VHdw+HSWdDlLlj3KoxsBcsmqDUhIntFBaK8ysiEVgOD1sTBx8EzV8Nj3WDLx1EnE5EyQgWivKtzGFw0A866Dza8EYzCXjIa8vOjTiYiaU4FoiLIyICcy4LJ/w5tBbOGwvgz4auYcyWKiAAqEBVLrYbQZwp0GwmbVsGo1vDaPyB/d9TJRCQNqUBUNGbQ/EIYvAQad4A5N8PYTrBpddTJRCTNqEBUVDUPht5PQI+x8PWH8HBbWHAP7N4ZdTIRSRMqEBWZGRzTE4YshSPPhHm3wSMd4LO3o04mImlABUJg/2w4bzz0+ids/xweaQ/zbg+WPBWRCksFQv6naVcYsgSO7gkL7obR7WDDsqhTiUhEkrmi3Dgz22RmKwtsO87MFpnZCjN7xsxqFnFsFzN7z8zWmtmwZGWUGKrVge4PwwWTIW8rjOkIc2+BnXlRJxORFEtmC2I80KXQtjHAMHc/BngaGFr4IDPLBEYApwNNCZYpbZrEnBLLEZ1hyGJo3hcW3g8PnQyfLI46lYikUNIKhLsvAL4utLkJsCB8PhfoEePQlsBad//Q3XcATwLdkpVTirHfAdD1Aeg7LeiPGNcFZg2DHd9GnUxEUiDVfRArga7h8/OAhjH2aQB8WuD1+nBbTGY2wMxyzSx38+bNCQsqBTRuH4zCbtkflowKpuv4aEHJx4lImZbqAnEZMMTMlgE1gFi3yViMbUVOQ+ruo909x91zsrOzExRTfqLK/nDGcLhkZjAR4ISz4dnr4PttUScTkSRJaYFw99Xu3sndTwAmArEmA1rPj1sWhwAbU5FP4tCoDQxcCCddCcvGB62J91+IOpWIJEFKC4SZHRT+zABuBh6KsdvrwOFm9nMzywJ6AzNSl1JKlFUNOt8Ol82BrOrweA+YNhjytkSdTEQSKJm3uU4EFgFNzGy9mfUjuCNpDbCaoFXwaLhvfTObCeDuu4ArgdnAu8Bkd1+VrJyyDxq2gIGvQNvfBMubjmgFq2dGnUpEEsS8HK0ylpOT47m5uVHHqJg2LofpQ+CLlcFAu9PvhuoHRp1KREpgZsvcPSfWexpJLYlRvxn0nw/tb4J3psOIlrByqpY5FSnDVCAkcSplQbsb4IqXg7UnnroUJvWB7V9EnUxESkEFQhKv3lHQ7wXoeCu8PzdoTbz1pFoTImWMCoQkR2YlOPlaGLQQspvA01fAE73gmw1RJxOROKlASHLVPRwunQVd7oJ1r8LIVrBsgloTImWACoQkX0YmtBoYtCYOPg6euRoe6wZbPo46mYgUQwVCUqfOYXDRDDjrPtjwRjAKe8loyM+POpmIxKACIamVkQE5lwWT/x3aCmYNhfFnwlexZl0RkSipQEg0ajWEPlOg20jYtApGtYbX/gH5u6NOJiKhuAqEmTUws9ZmdsqeR7KDSQVgBs0vhMFLoHEHmHMzjO0Em1ZHnUxEgEol7WBmdwG/Bt4B9vzzzvnfwj8i+6bmwdD7CVg5BWYOhYfbQrvfQZtrILNy1OlEKqwSCwRwDtDE3X9IdhipwMzgmJ7w83ZBv8S824IpO7qNgIOPjTqdSIUUzyWmDwH9M05SY/9sOG889PonbP8cHmkP824PljwVkZSKpwXxHbDczF4E/tuKcPerk5ZKpGlXaHQyPH8jLLgbVj8L3R6EBidEnUykwoinBTEDuA14DVhW4CGSXNXqQPeH4YLJkLcVxnSEubfAzryok4lUCCUWCHefQLA86J7C8ES4rVhmNs7MNpnZygLbmpnZYjNbbma5ZtayiGN3h/ssNzOtJlfRHdEZhiyG5n1h4f3w0MnwyeKoU4mUeyUWCDM7FXgfGAGMBNbEeZvreKBLoW13A7e6ezPglvB1LHnu3ix8dI3ju6S82+8A6PoA9J0Gu3fAuC4waxjs+DbqZCLlVjyXmO4FOrl7O3c/BegM3FfSQe6+APi68GagZvj8AIJlR0Xi17g9DFoELfvDklHBdB0f6Y5rkWSIp0BUdvf39rxw9zWU/q6ma4HhZvYpcA9wYxH77RdeglpsZucU94FmNiDcN3fz5s2ljCVlSpX94YzhcMnMYCLACWfDs9fB99uiTiZSrsRTIHLNbKyZnRo+HqH0ndSDgOvcvSFwHTC2iP0ODddIvQD4u5k1LuoD3X20u+e4e052dnYpY0mZ1KgNDFwIJ10Jy8YHrYn3X4g6lUi5EU+BGASsAq4GriEYUT2wlN93MTA1fP5vIGYntbtvDH9+CLwENC/l90l5l1UNOt8Ol82BrOrweA+YNhjytkSdTKTMi+cuph/c/W/u3t3dz3X3+/ZhVPVGoF34vANB5/ePmFltM6sSPq8LtCEoSiJFa9gCBr4CbX8TLG86ohWsnhl1KpEyrcgCYWaTw58rzOztwo+SPtjMJgKLgCZmtt7M+gH9gXvN7C3gDmBAuG+OmY0JD/0lwWWtt4D5wJ3urgIhJatUBU67BfrPg+p14cnz4al+8O1XUScTKZPMi1j60cwOdvfPzOxnsd5397RbDiwnJ8dzc3OjjiHpYNcOWPh3ePnu4BbZM4bDUecGcz6JyH+Z2bKwz/cnimxBuPtn4dPB7v5xwQcwOBlBRRKmUha0uwGuWBCsPfHUpTCpD2z/IupkImVGPJ3Uv4qx7fREBxFJinpNod8L0PFWeH8ujGgZ9FEU0XIWkf8prg9ikJmtIOhDKNj/8BFQYh+ESNrIrAQnXwuDFkJ2E3j6CniiF3yzIepkImmtuD6IA4DawF+BYQXe2u7uhUdIpwX1QUiJ8nfD0kfgxVshoxJ0+gscf5H6JqTCKlUfBODuvg4YAmwv8MDM6iQ6pEhKZGRCq4FBa+Lg4+CZq+GxbrAl7e65EIlccQXiifDnMiCXH0/1rX+mS9lW5zC4aAacdR9seCMYhb1kNOTnR51MJG0UeYmpLNIlJimVrZ/CM9fABy/Coa2DhYkOLHJ2F5FypbSXmPYc3MbMqofP+5jZ38zs0ESHFIlMrYbQZwp0GwmbVsGo1vDaP4L+CpEKLJ7bXEcB35nZccANwMfAP5OaSiTVzKD5hTB4CTTuAHNuhrGdYNPqqJOJRCaeArHLg+tQ3YD73f1+oEZyY4lEpObB0PsJ6DEWvv4QHm4LC+6B3TujTiaScvEUiO1mdiPQF3jOzDIp/XoQIunPDI7pCUOWwpFnwrzb4JEO8JmG/0jFEk+B+DXwA3CZu38ONACGJzWVSDrYPxvOGw+9/gnbP4dH2sO824N5nkQqgHim+/4ceBw4wMzOAr5398eSnkwkXTTtCkOWwNE9YcHdMLodbCjtmlkiZUc8dzH1ApYC5wG9gCVm1jPZwUTSSrU60P1huGAy5G2FMR1h7i2wMy/qZCJJE88lppuAFu5+sbtfRLAK3B+SG0skTR3RGYYshuZ9YeH98NDJ8MniqFOJJEU8BSLD3TcVeP1VnMdhZuPMbJOZrSywrZmZLTaz5WaWa2Yxlx01s4vN7P3wcXE83yeSEvsdAF0fgL7TYPcOGNcFZg2DHd9GnUwkoeL5H/3zZjbbzC4xs0uA54B413IcD3QptO1u4FZ3bwbcEr7+kXCupz8CJxK0WP5oZrXj/E6R1GjcHgYtgpb9YcmoYLqOjxZEnUokYeLppB4KPAwcCxwHjHb338Xz4e6+ACg886sDNcPnBxCsU11YZ2Cuu3/t7luAufy00IhEr8r+wWp1l8wMJgKccDY8ex18vy3qZCL7rFKc+70G7Abygdf38TuvBWab2T0EBap1jH0aAJ8WeL0+3PYTZjaAcG3rQw/VDCASkUZtYOBCmH87LB4Ja+bA2ffD4R2jTiZSavHcxXQ5wV1M5wI9gcVmdtk+fOcg4Dp3bwhcB4yN9bUxtsWcVdDdR7t7jrvnZGdn70MskX2UVQ063w6XzYGs6vB4D5g2GPK2RJ1MpFTiaUEMBZq7+1cAZnYgQYtiXCm/82LgmvD5v4ExMfZZD5xa4PUhwEul/D6R1GrYAga+Ai/fBa/+Hda+yOKmN/ObtxuwcWse9WtVZWjnJpzTPGajWCRtxNNJvZ5woaDQdn58+WdvbQTahc87AO/H2Gc20MnMaoed053CbSJlQ6UqcNot0H8e32QcQKulV3LDt8OpxTY2bM3jxqkrmPamljyV9BZPgdhAMDjuT2b2R2AxsNbMrjez64s70MwmAosI1rVeb2b9gP7AvWb2FnAHYf+BmeWY2RiAcEnT2wj6O14H/pyuy5yKFKt+M7r+cBv37uzJ6RlLeKHKULpmLCRv5y6Gz34v6nQixYrnEtMH4WOP6eHPEmd0dffzi3jrhBj75gKXF3g9jtJfxhJJG598s4t/0J3Z+S24q/IjPJA1gnN3v8oftu5LV55I8pVYINz91lQEESmv6teqyoateazxhvTY8ScuypzD0EqTmLPf72Dx98E4iozMqGOK/ERcI6JFpPSGdm5C1cpBAcgng/G7u9A1/162H5QDz/8OxnWGTe9GnFLkp1QgRJLsnOYN+Gv3Y2hQqyoGNKhVlau6d6DeoGeh+yPw1QfwUFuYfwfs+iHquCL/ZcFicUW8GSwOdLW735e6SKWXk5Pjubm5UccQ2Tvffgmzfw9vT4K6TaDrP+DQE6NOJRWEmS1z95xY7xXbgnD33QRLjYpIslSvC91Hw4VTYOd3wSWn536r6TokcvFcYlpoZg+aWVszO37PI+nJRCqawzvC4MVw4kB4fQyMbAXvPR91KqnAir3EBGBm82NsdnfvkJxIpadLTFJurM+FGVfBpnfg6B7Q5a5gCVSRBCvuElM8t7m2T3wkESnWITkw4GVY+HdYMBw+mAed74DjzgeLNVWZSOLFM1lfPTMba2azwtdNwxHRIpJMlbKg3Q0w8NWg83raIPjnubBlXdTJpIKIpw9iPME8SPXD12sIpuwWkVTIbgKXzoIz7w0uPY08CV57EHbvijqZlHPxFIi67j6ZYC0I3H0XwdoQIpIqGRnQ4vJgPeyfnwJzboKxHeHzFVEnk3IsngLxbTjFtwOYWSvgm6SmEpHYDjgEzn8Sej4K36yH0afCi3+Gnd9HnUzKoXgKxPXADKCxmS0EHgOuSmoqESmaGRzdHYYshWN/Da/cCw+1gXWvRp1Mypl41qR+g2D9htbAFcBRgOYpFolatTpwzkjoOw1274TxZ8Iz18D3auBLYsRzF9M4d9/l7qvcfSWQBcxMfjQRiUvj9jB4EbS+Ct54DB5sCe8+G3UqKQfiWjDIzEYBhKu7zQX+VdJBZjbOzDaZ2coC2yaZ2fLwsc7Mlhdx7DozWxHup5FvIiXJqg6d/gKXvwjVs2HShTCpL2z/POpkUobFc4npD8A2M3sImAPc6+6PxvHZ44EuhT7r1+7ezN2bAVOAqcUc3z7cN+YIPxGJocHxMGA+nPZHWDMbRrQMWhUlzJggEkuRBcLMuu95AEuBVsCbgIfbiuXuC4CYy4SamQG9gImlSi0iRcusDG2vh0GvQb1jgik7JpwdTCsusheKa0GcXeBxFkFxqFzg9b5oC3zh7u8X8b4Dc8xsmZkNKO6DzGyAmeWaWe7mzZv3MZZIOVL3F3DxM3D2/fDZ2zCqNbx6nwbYSdxKnKxvnz7crBHwrLsfXWj7KGCtu99bxHH13X2jmR1E0OdxVdgiKZYm6xMpwrbPYOZvYfWz8H/HBmtO1G8WdSpJA/s0WZ+ZZQP9gUYF93f3Uq24bmaVgO7ACUXt4+4bw5+bzOxpoCVQYoEQkSLUPBh6Pw7vzAgKxSMd4KQhcOqNkFUt6nSSpuK5i2k6cADwAvBcgUdpdQRWu/v6WG+aWXUzq7HnOdAJWBlrXxHZS027BgPsmveB1x4ILjt9+HLUqSRNxVMgqrn779x9srtP2fMo6SAzmwgsApqY2foCM8D2plDntJnVN7M9YyvqAa+a2VsEnePPubtWTRFJlKq1oOsDQf+EGTzWFaYPgbwtUSeTNBPPgkF/AV5z97QfHKc+CJG9tDMPXr4LFj4A1Q6EM+6GpudozYkKpNRrUoeuAZ41szwz22Zm281Mi+WKlAeVq0LHP8GAl4J+in9fAk9eCNs2RhpL0kM8A+VquHuGu1d195rh65qpCCciKXLwsXD5PPjVbcHqdSNOhNfHQn5+1MkkQsUNlDsy/Hl8rEfqIopISmRWgjZXw+DXoH5zeO76YALAL4sariTlXXG3uV4PDABijVVwoENSEolItOocBhdNh+WPw+ybgjud2t0Ara8JlkGVCiOpA+VSTZ3UIgm2/Qt4/new6mk46KhggN0hRQ5hkjKoVJ3UZtbCzP6vwOuLzGy6mT1gZnWSEVRE0kyNenDeeOg9MbgNdmxHeP73sOPbqJNJChTXSf0wsAPAzE4B7iRYTe4bYHTyo4lI2jjyjGA97BMuhcUjYGQrWPtC1KkkyYorEJnuvmc21l8Do8NBcn8AfpH8aCKSVvY7AM76G1z6PGRWgX/1gKlXwHcxJ22WcqDYAhHOmwRwGjCvwHslzuEkIuXUz06Cga/CKUNh5VPwYAtY8ZTWnCiHiisQE4GXzWw6kAe8AmBmvyC4zCQiFVXl/aDDzXDFAqj9M5jSD57oBVs/jTqZJFCRBcLdbwd+Q7Ay3Mn+v9udMoCrkh9NRNJevaOg31zociesezXom1gyGvJ3R51MEqDYkdTuvtjdn3b3bwtsW+PubyQ/moiUCRmZ0GoQDF4MDU+EWUNhXBfY9G7UyWQfxTMXk4hIyWr/DPpMgXNHw1dr4aG2MP+vsOuHqJNJKalAiEjimMFxv4YrX4ejzoGX74SHT4FPl0adTEpBBUJEEq96XegxBi74N/zwHxjbCWYOhR+2R51M9oIKhIgkzxGdggF2J14BSx+BEa1gzeyoU0mcklYgzGycmW0ys5UFtk0ys+XhY52ZLS/i2C5m9p6ZrTWzYcnKKCIpUKUGnH4X9JsDVfYPbod9qh/8Z3PUyaQEyWxBjAe6FNzg7r9292bu3gyYAkwtfJCZZQIjgNOBpsD5ZtY0iTlFJBUatoQrXoFTfw/vTIcRLWD5RA2wS2NJKxDuvgCIOQbfzAzoRaG1qUMtgbXu/qG77wCeBLolK6eIpFClLDj1d8FI7LpHwLSB8K/usGVd1Mkkhqj6INoCX7h7rJVIGgAFh2OuD7fFZGYDzCzXzHI3b1aTVaRMOOjIYE6nM+4J7nAaeRIsGqEBdmkmqgJxPrFbDwCxVksvsg3q7qPdPcfdc7KzsxMSTkRSICMDWvaHIUugUVuY/XsY0xE+X1nysZISKS8Q4QSA3YFJReyyHmhY4PUhgFZQFymvDjgELpgEPcbC1k9gdDt48TbY+X3UySq8KFoQHYHV7r6+iPdfBw43s5+bWRbQG5iRsnQiknpmcEzPYIDdMb3glXvgoTawbmHUySq0ZN7mOhFYBDQxs/Vm1i98qzeFLi+ZWX0zmwng7ruAK4HZwLvAZHdflaycIpJGqtWBc0dB36dh9w4YfwY8ex18rwmko6A1qUUkPe34FubfAYtHwv714Mx74cgzo05V7pRqTWoRkUhlVYfOt8PlL0C1A+HJC2DyRbD9i6iTVRgqECKS3hqcAANegg5/gPeeDwbYvfFPDbBLARUIEUl/mZXhlN/CoIVQ72iYcSU81hW++iDqZOWaCoSIlB11D4eLn4Wz/g4bl8Oo1rDwfti9K+pk5ZIKhIiULRkZkHNpMMCu8Wkw9xZ4pD189lbUycodFQgRKZtq1ofej0Ovx+A/X8Do9kGx2JkXdbJyQwVCRMouM2jaLWhNNLsguNw08iT4aEHUycoFFQgRKfuq1oZuD8JF4aQLE86G6VdC3pZoc5VxKhAiUn4c1g4GL4I218LyJ2DEicHaE7oltlRUIESkfKlcFX51K/SfF4zAnnwRTOoD2zTn595SgRCR8ql+M+g/H371Z1j7QtCayB0H+flRJyszVCBEpPzKrARtroFBrwUF49nrYMJZ8GWstcqkMBUIESn/DmwcdGB3fRC+WAmj2sCCe2D3zqiTpTUVCBGpGMzg+L4w5HVocjrMuw1GnwoblkWdLG2pQIhIxVKjHvSaAL2fgO++CpY5nX1TML24/EgyFwwaZ2abzGxloe1Xmdl7ZrbKzO4u4th1ZrbCzJabmRZ4EJHEO/LMYIDdCZfAogdhZCtY+2LUqdJKMlsQ44EuBTeYWXugG3Csux8F3FPM8e3dvVlRC1mIiOyz/Q6As+6DS2dBZhb8qzs8PQi++zrqZGkhaQXC3RcAhf8rDwLudPcfwn02Jev7RUTi9rPWMHAhtP0trJgMD7aAFU9V+AF2qe6DOAJoa2ZLzOxlM2tRxH4OzDGzZWY2oLgPNLMBZpZrZrmbN29OeGARqSAq7wen/QEGvAy1DoUp/WBib/hmfdTJIpPqAlEJqA20AoYCk83MYuzXxt2PB04HhiXld+MAAAqWSURBVJjZKUV9oLuPdvccd8/Jzs5OSmgRqUD+7+hgmdPOdwST/o04EZY+UiEH2KW6QKwHpnpgKZAP1C28k7tvDH9uAp4GWqY0pYhUbBmZcNKQYF6nhi1h5m/h0S6waXXUyVIq1QViGtABwMyOALKALwvuYGbVzazGnudAJ2AlIiKpVrsR9JkK5z4MX66Bh9vCS3fBrh1RJ0uJZN7mOhFYBDQxs/Vm1g8YBxwW3vr6JHCxu7uZ1TezmeGh9YBXzewtYCnwnLs/n6ycIiLFMoPjegcD7H7ZFV66Ax4+BT5dGnWypDMvR730OTk5npurYRMikkRrZsOz18O2DdByQNCxXaVG1KlKzcyWFTWcQCOpRUT2xhGdYchiaNkflo4OVrB7f27UqZJCBUJEZG9VqQFnDIfLZkPlavB4T5hyOXz7ZcnHliEqECIipXXoiTDwFTj1Rlg1LRhg99akcjPATgVCRGRfVKoCpw4LCsWBjeHpAfCvHrDl46iT7TMVCBGRRDjol8Elp9OHw6dLgsn/Fo2E/N1RJys1FQgRkUTJyIQTB8DgxdDoZJh9I4z9FXyxKupkpaICISKSaLUawgWTocdY2LIuGDcx7y+w8/uok+0VFQgRkWQwg2N6BgPsjjkPFgwPRmJ/vCjqZHFTgRARSabqB8K5D0GfKUEL4tEuwUC777dFnaxEKhAiIqnwi47B5H+thsCyR4NZYlfPLPm4CKlAiIikSpX9ocsd0O8FqFobnjwf/n0J/Cc9105TgRARSbVDToABL0GHm2H1c8EAuzf/lXYD7FQgRESiUCkLThkaLHV6UFOYPgQe6wZffxh1sv9SgRARiVL2EXDJc3Dm32DDGzCyNSx8AHbvijqZCoSISOQyMqBFP7hyKTRuD3P/AGM6wGdvRxsrWR9sZuPMbFO4OFDB7VeZ2XtmtsrM7i7i2C7hPmvNbFiyMoqIpJWa9aH3E3DeBNj2GYw+FV74E+zMiyROMlsQ44EuBTeYWXugG3Csux8F3FP4IDPLBEYApwNNgfPNrGkSc4qIpA8zOOocGLIEmp0Pr94Ho1rDR6+kPErSCoS7LwC+LrR5EHCnu/8Q7hPr3q6WwFp3/9DddxAsTdotWTlFRNJStTrQbQRcNB08HyacBTOugrytKYuQ6j6II4C2ZrbEzF42sxYx9mkAfFrg9fpwW0xmNsDMcs0sd/PmzQmOKyISscNOhUGLoPXVwa2wI1rCOzNS8tWVUvItP/6+2kAroAUw2cwO8x8vjG0xjivy5mB3Hw2MhmBN6gRmFRFJD1nVoNNtcHQPmHElTO4LR57F84f+ltsWbGHj1jzq16rK0M5NOKd5kf+e3mupbkGsB6Z6YCmQD9SNsU/DAq8PATamKJ+ISPqq3wz6z4eOt7J7zVxazz6DttufA/LZsDWPG6euYNqbGxL2dakuENOADgBmdgSQBRRexPV14HAz+7mZZQG9gdS0p0RE0l1mZTj5Wi6o9DdW5TfizspjmJh1O1X5nryduxk++72EfVUyb3OdCCwCmpjZejPrB4wDDgtvfX0SuNjd3czqm9lMAHffBVwJzAbeBSa7e9lcbUNEJEmWbqvN+Ttv4oad/VmXX4889gNg49bE3RKbtD4Idz+/iLf6xNh3I3BGgdczgfSe5lBEJEL1a1Vlw9Y8Ju9uz2Ta/2h7omgktYhIGTS0cxOqVs780baqlTMZ2rlJwr4j1XcxiYhIAuy5W2n47PeSdheTCoSISBl1TvMGCS0IhekSk4iIxKQCISIiMalAiIhITCoQIiISkwqEiIjEZJ5mi2TvCzPbDHxcysPr8tNpP8oinUd6KQ/nUR7OAXQeRfmZu2fHeqNcFYh9YWa57p4TdY59pfNIL+XhPMrDOYDOozR0iUlERGJSgRARkZhUIP5ndNQBEkTnkV7Kw3mUh3MAncdeUx+EiIjEpBaEiIjEpAIhIiIxVbgCYWZdzOw9M1trZsNivF/FzCaF7y8xs0apT1m8OM7hejN7x8zeNrMXzexnUeQsSUnnUWC/nmbmZpaWtyjGcx5m1iv8M1llZk+kOmM84vi9OtTM5pvZm+Hv1hmxPidqZjbOzDaFK1fGet/M7IHwPN82s+NTnbEkcZzDhWH2t83sNTM7LilB3L3CPIBM4APgMIL1sN8CmhbaZzDwUPi8NzAp6tylOIf2QLXw+aB0O4d4zyPcrwawAFgM5ESdu5R/HocDbwK1w9cHRZ27lOcxGhgUPm8KrIs6dxHncgpwPLCyiPfPAGYBBrQClkSduRTn0LrA79PpyTqHitaCaAmsdfcP3X0HwbrY3Qrt0w2YED5/CjjNzCyFGUtS4jm4+3x3/y58uRg4JMUZ4xHPnwXAbcDdwPepDLcX4jmP/sAId98C4O6bUpwxHvGchwM1w+cHABtTmC9u7r4A+LqYXboBj3lgMVDLzA5OTbr4lHQO7v7ant8nkvh3vKIViAbApwVerw+3xdzH3XcB3wAHpiRdfOI5h4L6EfxrKd2UeB5m1hxo6O7PpjLYXornz+MI4AgzW2hmi82sS8rSxS+e8/gT0MfM1hOsGX9VaqIl3N7+HUp3Sfs7XtFWlIvVEih8n288+0Qp7nxm1gfIAdolNVHpFHseZpYB3AdckqpApRTPn0clgstMpxL8S+8VMzva3bcmOdveiOc8zgfGu/u9ZnYS8M/wPPKTHy+h0v3veNzMrD1BgTg5GZ9f0VoQ64GGBV4fwk+byf/dx8wqETSli2uuplo854CZdQRuArq6+w8pyrY3SjqPGsDRwEtmto7gWvGMNOyojvd3arq773T3j4D3CApGOonnPPoBkwHcfRGwH8HEcWVNXH+H0p2ZHQuMAbq5+1fJ+I6KViBeBw43s5+bWRZBJ/SMQvvMAC4On/cE5nnYE5QmSjyH8NLMwwTFIR2vd0MJ5+Hu37h7XXdv5O6NCK6zdnX33GjiFime36lpBDcOYGZ1CS45fZjSlCWL5zw+AU4DMLNfEhSIzSlNmRgzgIvCu5laAd+4+2dRh9obZnYoMBXo6+5rkvZFUffWp/pBcAfDGoI7Nm4Kt/2Z4H8+EPzS/xtYCywFDos6cynO4QXgC2B5+JgRdebSnEehfV8iDe9iivPPw4C/Ae8AK4DeUWcu5Xk0BRYS3OG0HOgUdeYizmMi8Bmwk6C10A8YCAws8OcxIjzPFen4exXHOYwBthT4O56bjByaakNERGKqaJeYREQkTioQIiISkwqEiIjEpAIhIiIxqUCIiEhMFW0ktUhCmNmBwIvhy/8DdvO/MQHfuXvrSIKJJJBucxXZR2b2J+A/7n5P1FlEEkmXmEQSzMz+E/481cxeNrPJZrbGzO4M5/FfamYrzKxxuF+2mU0xs9fDR5toz0AkoAIhklzHAdcAxwB9gSPcvSXBSNg9s6HeD9zn7i2AHuF7IpFTH4RIcr3u4Tw/ZvYBMCfcvoJwfiagI9C0wLIjNc2shrtvT2lSkUJUIESSq+BMuvkFXufzv79/GcBJ7p6XymAiJdElJpHozQGu3PPCzJpFmEXkv1QgRKJ3NZATLkD/DsGsnSKR022uIiISk1oQIiISkwqEiIjEpAIhIiIxqUCIiEhMKhAiIhKTCoSIiMSkAiEiIjH9P41SxPIt8kmxAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "if(GEO.rank==0):\n",
    "    import matplotlib.pyplot as plt\n",
    "    plt.plot(tSinker, ySinker, \"o\") \n",
    "    plt.plot(tSinker, ySinker) \n",
    "    plt.xlabel('Time')\n",
    "    plt.ylabel('Sinker position')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}