notebooks/notebook3.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "notebooks/notebook3.ipynb",
      "version": "0.3.2",
      "provenance": [],
      "collapsed_sections": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python2",
      "display_name": "Python 2"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/stoiver/anuga-clinic-2018/blob/master/notebooks/notebook3.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MFkDMm10DPK7",
        "colab_type": "text"
      },
      "source": [
        "# CSDMS ANUGA Clinic 2018\n",
        "\n",
        "## Notebook 3: Setting up a simple Erosion Operator \n",
        "\n",
        "Here we go through the process of creating an operator (fractional step operator) which implements a simple sand dune erosion operator. \n",
        "\n",
        "### Installation\n",
        "\n",
        "These notebooks have been designed to run in the google `colaboratory` environment, which provides a jupyter notebook environment running on a virtual machine on the cloud. To use this environment you need a google account so that your notebook can be saved on google drive. \n",
        "\n",
        "To start interacting with the notebook follow the \n",
        "`View in Colaboratory` link above. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Dt-Ti7kdQX3C",
        "colab_type": "text"
      },
      "source": [
        "## Setup Environment\n",
        "\n",
        "If on github, first follow the link `View in Colaboratory' to start running on google's colab environment. Then ....\n",
        "\n",
        "Run the following cell to install the dependencies and some extra code for visualising on Colaboratory.\n",
        "\n",
        "Wait until you see the comment *(5) Ready to go* before proceeding to subsequent commands. \n",
        "\n",
        "The install should take less than a minute (and quicker if you have already run this earlier)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "xdCgMgk7_ZcW",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "try:\n",
        "  # On colab we can install all the packages we need via the notebook\n",
        "  #\n",
        "  # First download the clinic repository\n",
        "  import os\n",
        "  os.chdir('/content')\n",
        "  !git clone https://github.com/stoiver/anuga-clinic-2018.git\n",
        "\n",
        "  # Now install environment using tool\n",
        "  !/bin/bash /content/anuga-clinic-2018/anuga_tools/install_anuga_colab.sh\n",
        " \n",
        "except:\n",
        "  pass\n",
        "\n",
        "# Make inline animate code available\n",
        "if not 'workbookDir' in globals():\n",
        "    workbookDir = os.getcwd()\n",
        "\n",
        "import sys\n",
        "sys.path.append(os.path.join(workbookDir,\"anuga-clinic-2018\"))\n",
        "                \n",
        "%matplotlib inline\n",
        "\n",
        "# Allow inline jshtml animations\n",
        "from matplotlib import rc\n",
        "rc('animation', html='jshtml')\n",
        "\n",
        "import anuga"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "n6CIHtFYYbkx",
        "colab_type": "text"
      },
      "source": [
        "## The Code \n",
        "\n",
        "Here we use the `Sanddune_erosion_operator` which was contibuteed to the anuga code by Ted Rigby. \n",
        "\n",
        "[The code is available from anuga github reposiory](https://github.com/stoiver/anuga_core/blob/master/anuga/operators/sanddune_erosion_operator.py)\n",
        "\n",
        "We need to setup a quantity to represent the bedrock (can erode down to the bedrock) and also the fractional step operator that does the erosion."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "1l0-N381KKGO",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "\"\"\"sand dune erosion test                          run_sanddune_erosion.py\n",
        "\n",
        "This tests the sanddune_erosion operator confirming that; \n",
        "    1. flow creates erosion when bed shear > critical and that erosion rates \n",
        "\tare higher in higher bed shear zones (typically higher velocity areas).\n",
        "\t2. that erosion is augmented by collapse of the sand face whenever\n",
        "\terosion creates slopes > the angle of repose.  \n",
        "\tthis process leads to widening of the notch (laterally)and head\n",
        "    like recession of \n",
        "\tthe main scour zone  as scour acts to steepen the longitudinal \n",
        "\tbed, triggering collapse of the steep bed and reshaping back \n",
        "\tto the angle of repose).\n",
        "\t3. that the operator can handle multiple erosion polygons with\n",
        "    different base levels\n",
        "\t \n",
        "\"\"\"\n",
        "\n",
        "#------------------------------------------------------------------------------\n",
        "# Import necessary modules\n",
        "#------------------------------------------------------------------------------\n",
        "import anuga\n",
        "from anuga import Region, Geo_reference\n",
        "import numpy as num\n",
        "\n",
        "x0 = 314036.58727982\n",
        "y0 = 6224951.2960092\n",
        "geo = Geo_reference(56, x0, y0)\n",
        "\n",
        "#------------------------------------------------------------------------------\n",
        "# Function to describe topography as function of X and y\n",
        "#------------------------------------------------------------------------------\n",
        "def topography(x,y):\n",
        "    \"\"\"Complex topography defined by a function of vectors x and y.\"\"\"\n",
        "    print ' Creating topography....'\n",
        "    \n",
        "    z = 0.0*(x)                             # horizontal plane \n",
        "\n",
        "    N = len(x)\n",
        "   \n",
        "    for i in range(N):\n",
        "        # First notched sand dune across Channel\n",
        "        if 110 < x[i] <= 120:\n",
        "            z[i] +=  1.0*(x[i]-110)       # Sloping U/S Face 1:1\n",
        "        if 120 < x[i] < 130 :\n",
        "            z[i] +=  10                   # Crest of Embankment at +10\n",
        "        if 130 <= x[i] < 140:\n",
        "            z[i] +=  10-1.0*(x[i] -130.0)  # Sloping D/S Face at 1:1\n",
        "\n",
        "        # add notch in crest 1m wide by nom 300 deep\n",
        "        # note sides are near vertical so will collapse back to repose even without erosion\n",
        "        \n",
        "        if 117 <= x[i] <= 133 and 2 <= y[i] <= 3:\n",
        "            z[i] =  9                   # add 4m Notch in Embankment crest \n",
        "            \n",
        "        # second lower plain sand dune across Channel\n",
        "        if 230.0 < x[i] <= 240.0:\n",
        "            z[i] +=  0.2*(x[i] -230)      # Sloping U/S Face 1:5        \n",
        "        if 240 < x[i] < 250 :\n",
        "            z[i] +=  2                    # Crest of Embankment at +2.0\n",
        "        if 250 <= x[i] < 260:\n",
        "            z[i] +=  2-0.2*(x[i] -250)   # Sloping D/S Face at 1:5      \n",
        "    return z\n",
        "\n",
        "#------------------------------------------------------------------------------\n",
        "# build the check points and erosion polygons now so available to all processors\n",
        "#------------------------------------------------------------------------------\n",
        "\n",
        "polygon1    = num.array([ [106, 0.0], [144, 0.0], [144, 5.0], [106, 5.0] ])\n",
        "polygon2    = num.array([ [230, 0.0], [260, 0.0], [260, 5.0], [230, 5.0] ])\n",
        "poly1nsbase = 5  # set poly no scour base level in m\n",
        "poly2nsbase = 3\n",
        "\n",
        "polygon1 += num.array([x0,y0])\n",
        "polygon2 += num.array([x0,y0])\n",
        "\n",
        "#------------------------------------------------------------------------------\n",
        "# Setup computational domain\n",
        "#------------------------------------------------------------------------------\n",
        "\n",
        "print '>>>>> DUNE EROSION TEST SCRIPT'\n",
        "print '>>>>> Setting  up Domain.'\n",
        "length = 360.\n",
        "width = 5.\n",
        "dx = dy = 1.0            # Resolution: Length of subdivisions on both axes\n",
        "print '>>>>> Domain has L = %f, W = %f, dx=dy= %f' %(length,width,dx) \n",
        "points, vertices, boundary = anuga.rectangular_cross(int(length/dx), int(width/dy), len1=length, len2=width)\n",
        "\n",
        "\n",
        "\n",
        "domain = anuga.Domain(points, vertices, boundary, geo_reference = geo)\n",
        "domain.set_flow_algorithm('DE0')\n",
        "domain.set_name('sanddune') # Output name\n",
        "\n",
        "domain.set_store_vertices_uniquely(False)\n",
        "domain.set_quantities_to_be_stored({'elevation': 2,'stage': 2,'xmomentum': 2,'ymomentum': 2})\n",
        "\n",
        "domain.set_quantity('elevation', topography, location='centroids')           # elevation is a function\n",
        "domain.set_quantity('friction', 0.01, location='centroids')                  # Constant friction\n",
        "domain.set_quantity('stage', expression='elevation', location='centroids')   # Dry initial condition\n",
        "\n",
        "print domain.statistics()\n",
        "\n",
        "\n",
        "#------------------------------------------------------------------------------\n",
        "# Set up a new quantity to represent the bedrock\n",
        "#------------------------------------------------------------------------------\n",
        "\n",
        "# get the indices of triangles in each erosion poly so can setup nsbase_ in domain\n",
        "poly1ind = (Region(domain, polygon=polygon1)).indices\n",
        "poly2ind = (Region(domain, polygon=polygon2)).indices\n",
        "print '>>>>> poly1ind is of length ', len(poly1ind), ' and contains triangles', poly1ind[0], ' to ' , poly1ind[-1]\n",
        "print '>>>>> poly2ind is of length ', len(poly2ind), ' and contains triangles', poly2ind[0], ' to ' , poly2ind[-1]\n",
        "\n",
        "# get the initial model surface elevation\n",
        "nsbase_elev_c = domain.get_quantity('elevation').get_values(location='centroids')\n",
        "print '>>>>> nsbase_elev_c is of length ',len(nsbase_elev_c) \n",
        "\n",
        "\n",
        "# build the no scour base surface  by combining initial elev where < nsbase and nsbase in each scour poly\n",
        "nsbase_elev_c[poly1ind] = num.minimum(nsbase_elev_c[poly1ind], poly1nsbase)\n",
        "nsbase_elev_c[poly2ind] = num.minimum(nsbase_elev_c[poly2ind], poly2nsbase)\n",
        "\n",
        "\n",
        "# create new Anuga quantity and assign nsbase_elev_c values to domain so can distribute later\n",
        "anuga.Quantity(domain, name='nsbase_elevation', register=True)\n",
        "domain.set_quantity('nsbase_elevation',nsbase_elev_c, location='centroids')\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "#------------------------------------------------------------------------------\n",
        "# Print out polygon points and ids as check\n",
        "#------------------------------------------------------------------------------\n",
        "print '>>>>> erosion polygon1 contains ', polygon1\n",
        "print '>>>>> erosion polygon2 contains ', polygon2\n",
        "\n",
        "\n",
        "\n",
        "#------------------------------------------------------------------------------\n",
        "# Setup boundary conditions on the distributed domain\n",
        "#------------------------------------------------------------------------------\n",
        "Bi = anuga.Dirichlet_boundary([12.0,0.0,0.0]) # Inflow at depth\n",
        "Br = anuga.Reflective_boundary(domain)              # Solid reflective side walls\n",
        "Bo = anuga.Dirichlet_boundary([-5, 0, 0])           # uncontrolled outflow \n",
        "\n",
        "domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})\n",
        "\n",
        "#------------------------------------------------------------------------------\n",
        "# Setup sanddune erosion operator \n",
        "#------------------------------------------------------------------------------\n",
        "\n",
        "print '>>>>> Setting up Erosion Area(s) to test...'\n",
        "\n",
        "# power up the erosion operator\n",
        "from anuga import Sanddune_erosion_operator\n",
        "\n",
        "# assign ns base elevations to partitioned domains\n",
        "nsbase_elev_c = domain.get_quantity('nsbase_elevation').get_values(location='centroids')\n",
        "\n",
        "poly1ind = (Region(domain, polygon=polygon1)).indices\n",
        "poly2ind = (Region(domain, polygon=polygon2)).indices\n",
        "indices_union = list(set(poly1ind) | set(poly2ind))\n",
        "\n",
        "# setup and create operator within polys setting the scour base elevations for each poly\n",
        "op0 = Sanddune_erosion_operator(domain, base=nsbase_elev_c, indices=indices_union, Ra=45)   # both dunes\n",
        "\n",
        "\n",
        "#------------------------------------------------------------------------------\n",
        "# Evolve simulation through time  \n",
        "#------------------------------------------------------------------------------\n",
        "for t in domain.evolve(yieldstep=2.0, duration=100.0):\n",
        "    domain.print_timestepping_statistics()\n",
        "    \n",
        "\n",
        "#  run completed - tidy up \n",
        "print ' >>>>> Simulation completed successfully '\n",
        "\n",
        " \n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QpGVChWjwi5G",
        "colab_type": "text"
      },
      "source": [
        "## Visualise Flow\n",
        "\n",
        "Let's open up the `swwfile` and look at the evolution of the `depth` the `stage` and the `elev`."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "nSJhl9y3K_aJ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Create a wrapper for contents of sww file\n",
        "swwfile = 'sanddune.sww'\n",
        "splotter = anuga.SWW_plotter(swwfile)\n",
        "\n",
        "# plot stage\n",
        "splotter.triang.set_mask(None)\n",
        "for i,time in enumerate(splotter.time):\n",
        "  #print time\n",
        "  splotter.save_stage_frame(frame=i)\n",
        "  \n",
        "splotter.make_stage_animation()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "JVEM8ZO0ai_m",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# plot depth\n",
        "splotter.triang.set_mask(None)\n",
        "for i,time in enumerate(splotter.time):\n",
        "  print time\n",
        "  splotter.save_depth_frame(frame=i)\n",
        "  \n",
        "splotter.make_depth_animation()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LHcS1-qcxfVW",
        "colab_type": "text"
      },
      "source": [
        "## 3D visualisation\n",
        "\n",
        "We can use the `mpl` module to produce some 3D plots. We need to create vertex values given the centroid values.\n",
        "\n",
        "This code produes an extra extraneous plot (any suggestion to fix that welcome).\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_6bT1VYTZ0u0",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "X = np.vstack((splotter.xc, splotter.yc)).transpose()\n",
        "E = splotter.elev\n",
        "\n",
        "splotter.nodes = np.vstack((splotter.x,splotter.y)).transpose()\n",
        "import scipy\n",
        "\n",
        "stage_v = []\n",
        "for i,S in enumerate(splotter.stage):\n",
        "  #print i, S.shape, X.shape\n",
        "  Interp = scipy.interpolate.NearestNDInterpolator(X,S)\n",
        "  stage_v.append(Interp(splotter.nodes))\n",
        "  \n",
        "stage_v = np.array(stage_v)\n",
        "\n",
        "elev_v = []\n",
        "for i,E in enumerate(splotter.elev):\n",
        "  #print i, E.shape, X.shape\n",
        "  Interp = scipy.interpolate.NearestNDInterpolator(X,E)\n",
        "  elev_v.append(Interp(splotter.nodes))  \n",
        "\n",
        "elev_v = np.array(elev_v)\n",
        "\n",
        "from matplotlib import animation\n",
        "from mpl_toolkits.mplot3d import Axes3D\n",
        "\n",
        "\n",
        "fig = plt.figure()\n",
        "#fig, ax = plt.subplots()\n",
        "ax = plt.gca(projection='3d')\n",
        "\n",
        "def update_image(i):\n",
        "    ax.clear()\n",
        "    ax.plot_trisurf(splotter.triang, elev_v[i*step], color='navajowhite', edgecolors='none', antialiased=False, shade=True)\n",
        "    ax.plot_trisurf(splotter.triang, stage_v[i*step], color='dodgerblue', edgecolors='none', alpha=1.0, linewidth=0, antialiased=False, shade=True)\n",
        "    ax.set_zlim(elev_v.min(), np.max([stage_v.max(), elev_v.max()]))\n",
        "    \n",
        "\n",
        "# increase step to skip frames; increase interval to increase movie speed\n",
        "step = 4\n",
        "anim = animation.FuncAnimation(fig, update_image, frames=len(splotter.depth)/step, interval=20);\n",
        "\n",
        "anim"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-RhWS6HTye-3",
        "colab_type": "text"
      },
      "source": [
        "## Examine elevation at points\n",
        "\n",
        "As in the previous notebook, it is instructive to examine the value of the quantity of interest at specific locations. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "sk0rK5COca-G",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Alias\n",
        "xc = splotter.xc\n",
        "yc = splotter.yc\n",
        "\n",
        "# Observation points across the first step\n",
        "point_observations = [[125.3, 0.0], \n",
        "                      [125.3, 0.5], \n",
        "                      [125.3, 1.0], \n",
        "                      [125.3, 1.5],\n",
        "                      [125.3, 2.0],\n",
        "                      [125.3, 2.5],\n",
        "                     ]\n",
        "\n",
        "# Find nearest centroid to observation points\n",
        "nearest_points = []\n",
        "for row in point_observations:\n",
        "    nearest_points.append(np.argmin( (xc-row[0])**2 + (yc-row[1])**2 ))\n",
        "    \n",
        "fig, ax = plt.subplots()\n",
        "\n",
        "for i, npt in enumerate(nearest_points):\n",
        "  ax.plot(splotter.time, splotter.elev[:,npt], label=str(point_observations[i][1]) )\n",
        "\n",
        "ax.legend(loc = 'center right')\n",
        "plt.title('Elevation')\n",
        "plt.xlabel('time')\n",
        "plt.ylabel('height')\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "print nearest_points\n",
        " "
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "beEW6rWnibak",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "95v97G0jvGXy",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": 0,
      "outputs": []
    }
  ]
}