notebooks/notebook4.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "notebooks/notebook4.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/notebook4.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 4: Setting up a Pumping Operator separated by a River Wall\n",
        "\n",
        "Here we go through the process of creating an operator (fractional step operator) which implements a simple pump, and separate the inlet from teh outlet with a River wall.  \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": [
        "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\"))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uCN-Ah6pIg2k",
        "colab_type": "text"
      },
      "source": [
        "### Import ANUGA and Tools\n",
        "\n",
        "Once `anuga` and `anuga_clinic-2018` have been installed we can now `import` our standard libraries."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ZOHiEet7IYhe",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\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",
        "This code is taken from `anuga_core/examples/structures/run_pump.py`.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "m9Jj58FyL7oe",
        "colab_type": "text"
      },
      "source": [
        "### Parameters to Setup Domain\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "1l0-N381KKGO",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "boundaryPolygon = [ [0., 0.], [0., 100.], [100.0, 100.0], [100.0, 0.0]]\n",
        "wallLoc = 50.\n",
        "# The boundary polygon + riverwall breaks the mesh into multiple regions\n",
        "# Must define the resolution in these areas with an xy point + maximum area\n",
        "# Otherwise triangle.c gets confused\n",
        "length = 2.0\n",
        "res = length*length*0.5\n",
        "regionPtAreas = [ [99., 99., res],\n",
        "                  [1., 1., res],\n",
        "                  [45, 50, res],\n",
        "                  [55, 50, res]]\n",
        "\n",
        "wallHeight=10.\n",
        "InitialOceanStage=2.\n",
        "InitialLandStage=6.\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "asQCz8wjQnbf",
        "colab_type": "text"
      },
      "source": [
        "### RiverWall Setup and Utility Procedure\n",
        "\n",
        "We will be setting up a `riverwall` which implicitly has an associated polyline. Setup a procedure to pull out this polyline so that it can be used in the construction of the `domain`.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "mSt6iBYLQnuo",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "riverWall = { 'centralWall':\n",
        "                           [ [wallLoc, 0.0, wallHeight],\n",
        "                             [wallLoc, 100.0, wallHeight]] \n",
        "                        }\n",
        "\n",
        "riverWall_Par = {'centralWall':{'Qfactor':1.0}}\n",
        "\n",
        "def toPolyline(riverWall):\n",
        "  return [np.array(riverWall.values()[0])[:,0:2].tolist()]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uNFTBGMLRepO",
        "colab_type": "text"
      },
      "source": [
        "### Setup Pump geometry\n",
        "\n",
        "First setup the inlet and outlet polygonal regions of the pump. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7yMvW8gywa_F",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "end_point0 = [49.0,50.0]\n",
        "end_point1 = [51.0,50.0]\n",
        "\n",
        "end_points = [end_point0, end_point1]\n",
        "\n",
        "inlet1_poly = [[[end_point0[0]-10, 45.0],[end_point0[0]-10,55],\n",
        "                [end_point0[0],55],[end_point0[0],45],[end_point0[0]-10, 45.0]]]\n",
        "              \n",
        "inlet2_poly = [[[end_point1[0], 45.0],[end_point1[0],55],\n",
        "                [end_point1[0]+10,55],[end_point1[0]+10,45],[end_point1[0], 45.0]]]\n",
        "              "
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EnFwhpcfQFID",
        "colab_type": "text"
      },
      "source": [
        "### Create and Setup Domain\n",
        "\n",
        "Create  a domain with a breakline down the middle (where the riverWall will occur) and two polygonal regions where the inlet and oulet polygons will sit. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kYnnXxk-QFax",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "domain = anuga.create_domain_from_regions(boundaryPolygon, \n",
        "                           boundary_tags={'left': [0],\n",
        "                                        'top': [1],\n",
        "                                        'right': [2],\n",
        "                                        'bottom': [3]},\n",
        "                           maximum_triangle_area = 10.0,\n",
        "                           minimum_triangle_angle = 28.0,\n",
        "                           interior_regions =[ ], #[ [higherResPolygon, 1.*1.*0.5],\n",
        "                                                  #  [midResPolygon, 3.0*3.0*0.5]],\n",
        "                           breaklines=toPolyline(riverWall)+inlet1_poly+inlet2_poly,\n",
        "                           regionPtArea=regionPtAreas,\n",
        "                           use_cache=False,\n",
        "                           verbose=False\n",
        "                           )\n",
        "\n",
        "\n",
        "#domain.set_flow_algorithm('DE1')\n",
        "domain.set_name('run_pump')\n",
        "domain.set_store_vertices_uniquely(True)\n",
        "\n",
        "d_plotter = anuga.Domain_plotter(domain)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EmtgsITKO_A4",
        "colab_type": "text"
      },
      "source": [
        "### Setup Initial Conditions\n",
        "\n",
        "We will have a slight slope down towards the right, with water height `InitialOceanStage` on the  right, and `InitialLandStage` on the left."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "83zqqhutO_Vt",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def topography(x,y):\n",
        "    return -x/150. \n",
        "\n",
        "def stagefun(x,y):\n",
        "    stg = InitialOceanStage*(x>=wallLoc) + InitialLandStage*(x<wallLoc)\n",
        "    return stg \n",
        "\n",
        "\n",
        "# NOTE: Setting quantities at centroids is important for exactness of tests\n",
        "domain.set_quantity('elevation',topography,location='centroids')     \n",
        "domain.set_quantity('stage', stagefun,location='centroids')    "
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "K_h6MlQFNsbD",
        "colab_type": "text"
      },
      "source": [
        "### Setup RiverWall \n",
        "\n",
        "From the `riverWall` polyline and elevation data create a riverWall  `operator`. Looking back at the definition of `riverWall` we see that this will create a wall of height `wallHeight` down the middle of the domain. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "z3C0XzD4Ns8r",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "domain.riverwallData.create_riverwalls(riverWall,riverWall_Par,verbose=False) "
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RQih8Hn9N7Hk",
        "colab_type": "text"
      },
      "source": [
        "### Boundary conditions\n",
        "\n",
        "Simple reflective BC all around\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Ipai69CAN4Gv",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "Br = anuga.Reflective_boundary(domain)\n",
        "domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom':Br})"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "o2ew3M_qNKw3",
        "colab_type": "text"
      },
      "source": [
        "### Setup pump Operator\n",
        "\n",
        "(1) First setup the pump characteristics\n",
        "\n",
        "(2) Then locate and create the pump using the operator"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "QnzR4HquNLgL",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "pump_function = anuga.pumping_station_function(\n",
        "            domain=domain,\n",
        "            pump_capacity=100.0,\n",
        "            hw_to_start_pumping=0.0,\n",
        "            hw_to_stop_pumping=-1.0,\n",
        "            initial_pump_rate=100.0, \n",
        "            pump_rate_of_increase = 50.0, \n",
        "            pump_rate_of_decrease = 50.0, \n",
        "            verbose=True)\n",
        "\n",
        "\n",
        "end_points = [[49.0,50.0],[51.0,50.0]]\n",
        "pump = anuga.Internal_boundary_operator(domain, pump_function,\n",
        "                                        width = 10.0,\n",
        "                                        height = 1.0,\n",
        "                                        apron = 10.0,\n",
        "                                        end_points=end_points,\n",
        "                                        verbose=True)\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "v7AZXi7uMv2v",
        "colab_type": "text"
      },
      "source": [
        "### Evolve\n",
        "\n",
        "Monitor the amount of water on each side of the wall. The sum should remain constant, and the change should be match the pump capacity."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "EblPWS4CMyAQ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "region1 = anuga.Region(domain, polygon=[[0.0,0.0], [50.0,0.0], [50.0, 100.0], [0.0,100.0]])\n",
        "region2 = anuga.Region(domain, polygon=[[50.0,0.0], [100.0,0.0], [100.0, 100.0], [50.0,100.0]])\n",
        "\n",
        "for t in domain.evolve(yieldstep=2, duration=120):\n",
        "    domain.print_timestepping_statistics()\n",
        "    stage = domain.get_quantity('stage')\n",
        "    elev  = domain.get_quantity('elevation')\n",
        "    height = stage - elev\n",
        "    \n",
        "    #d_plotter.plot_depth_frame(vmin=2.0, vmax=6.0)\n",
        "\n",
        "    print anuga.indent + 'Integral1 = ', height.get_integral(region=region1)\n",
        "    print anuga.indent + 'Integral2 = ', height.get_integral(region=region2)\n",
        "    print anuga.indent + 'Total Integral = ', height.get_integral()\n",
        "    #pump.print_timestepping_statistics()"
      ],
      "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": "OCRY9K9R-DDX",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Create a wrapper for contents of sww file\n",
        "swwfile = 'run_pump.sww'\n",
        "splotter = anuga.SWW_plotter(swwfile);\n",
        "\n",
        "plt.triplot(splotter.triang, linewidth = 0.4);"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "nSJhl9y3K_aJ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# 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, vmin=1.0, vmax=7.0)\n",
        "  \n",
        "splotter.make_stage_animation()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LHcS1-qcxfVW",
        "colab_type": "text"
      },
      "source": [
        "## Animation\n",
        "\n",
        "We ca use the `FuncAnimation` function to animate a sequence of timeslices (in this case  `stage`)\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "cFciK00M1zUF",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def update_frame(i):\n",
        "  \n",
        "  ax.clear()\n",
        "  j = step*i\n",
        "  artist = ax.tripcolor(splotter.triang, \n",
        "              facecolors = splotter.stage[j], \n",
        "              edgecolors='k', \n",
        "              cmap='Spectral',\n",
        "              vmax = 6.0,\n",
        "              vmin = 2.0) \n",
        "  ax.set_title('Stage (time %05.2f)' % splotter.time[j],loc='right')\n",
        "  fig.colorbar(artist, cax=cax)\n",
        "  \n",
        "\n",
        "from matplotlib import animation\n",
        "  \n",
        "fig = plt.figure();\n",
        "ax = plt.gca();\n",
        "cax = fig.add_axes([0.91, 0.125, 0.03, 0.75])\n",
        "\n",
        "# increase step to skip frames; increase interval to increase movie speed\n",
        "step = 4\n",
        "anim = animation.FuncAnimation(fig, update_frame, frames=len(splotter.depth)/step, interval=200);\n",
        "\n",
        "\n",
        "plt.close()\n",
        "\n",
        "anim\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2Mf_ALyZxo72",
        "colab_type": "text"
      },
      "source": [
        "###  3D Surface Plots\n",
        "\n",
        "We can use the mpl module to produce some 3D plots. We need to create vertex values given the centroid values."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_6bT1VYTZ0u0",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "\n",
        "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",
        "Interp = scipy.interpolate.NearestNDInterpolator(X,E)\n",
        "elev_v = Interp(splotter.nodes)\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "from matplotlib.animation import FuncAnimation\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, 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 = FuncAnimation(fig, update_image, frames=len(splotter.depth)/step, interval=200);\n",
        "\n",
        "plt.close()\n",
        "\n",
        "anim"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "llOycSk7UlIA",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": 0,
      "outputs": []
    }
  ]
}