notebooks/notebook2.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "notebooks/notebook2.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/notebook2.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 2: Merewether Flood Case Study\n",
        "\n",
        "Here we look at a case study of a flood in the community of Merewether.  We will add a flow using an `Inlet_operator` and extract flow details at various points by interagating the `sww` file which is produced from each ANUGA run. \n",
        "\n",
        "This notebook is based on the the validation test [merewether](https://github.com/GeoscienceAustralia/anuga_core/tree/master/validation_tests/case_studies/merewether), provided in the ANUGA distribution. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Dt-Ti7kdQX3C",
        "colab_type": "text"
      },
      "source": [
        "## Setup the 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",
        "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": "bWOxaRmEDg4p",
        "colab_type": "text"
      },
      "source": [
        "## Setup Inline graphics and import anuga"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "hUT6s7gXDA95",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "%matplotlib inline\n",
        "\n",
        "# Allow inline jshtml animations\n",
        "from matplotlib import rc\n",
        "rc('animation', html='jshtml')"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "15yi1XYR9wpG",
        "colab_type": "text"
      },
      "source": [
        "## Import ANUGA"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "OSlhVDrx-Bay",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import anuga"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "p_QOy0bRZ37C",
        "colab_type": "text"
      },
      "source": [
        "## Read in Data\n",
        "\n",
        "We have included some topography data and extent data in our `anuga-clinic-2018` download. \n",
        "\n",
        "Let's read that in and create a mesh associated with it. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "nCKdVO1BPAhx",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "data_dir = '/content/anuga-clinic-2018/data/merewether'\n",
        "\n",
        "# Polygon defining broad area of interest\n",
        "bounding_polygon = anuga.read_polygon(os.path.join(data_dir,'extent.csv'))\n",
        "\n",
        "\n",
        "# Polygon defining particular area of interest\n",
        "merewether_polygon = anuga.read_polygon(os.path.join(data_dir,'merewether.csv'))\n",
        "\n",
        "\n",
        "# Elevation Data\n",
        "topography_file = os.path.join(data_dir,'topography1.asc')\n",
        "\n",
        "\n",
        "# Resolution for most of the mesh\n",
        "base_resolution = 80.0  # m^2\n",
        "\n",
        "# Resolution in particular area of interest \n",
        "merewether_resolution = 25.0 # m^2\n",
        "\n",
        "interior_regions = [[merewether_polygon, merewether_resolution]]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mfgw_HkDarZY",
        "colab_type": "text"
      },
      "source": [
        "## Create and View Domain\n",
        "\n",
        "\n",
        "Note that we use a base_resolution to ensure a reasonable refinement overthe whole region, and we use interior_regions to refine the mesh in the area of interest. In this case we pass a list of polygon, resolution pairs. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Bb2oWmRraSU6",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "domain = anuga.create_domain_from_regions(\n",
        "            bounding_polygon,\n",
        "            boundary_tags={'bottom': [0],\n",
        "                           'right':  [1],\n",
        "                           'top':    [2],\n",
        "                           'left':   [3]},\n",
        "            maximum_triangle_area=base_resolution,\n",
        "            interior_regions=interior_regions)\n",
        "\n",
        "\n",
        "domain.set_name('merewether1') # Name of sww file\n",
        "dplotter = anuga.Domain_plotter(domain)  \n",
        "plt.triplot(dplotter.triang, linewidth = 0.4);"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZGGxHd8EEvea",
        "colab_type": "text"
      },
      "source": [
        "## Setup Initial Conditions\n",
        "\n",
        "We have to setup the values of various quantities associated with the domain. In particular we need to setup the `elevation` the elevation of the bed or the bathymetry. In this case we will do this using the DEM file `topography1.asc` .\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "xLn6Mozrnedv",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "domain.set_quantity('elevation', filename=topography_file, location='centroids') # Use function for elevation\n",
        "domain.set_quantity('friction', 0.01, location='centroids')                        # Constant friction \n",
        "domain.set_quantity('stage', expression='elevation', location='centroids')         # Dry Bed \n",
        "\n",
        "plt.tripcolor(dplotter.triang, \n",
        "              facecolors = dplotter.elev, \n",
        "              cmap='Greys_r')\n",
        "plt.colorbar();\n",
        "plt.title(\"Elevation\");"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jlgKD23QvDhF",
        "colab_type": "text"
      },
      "source": [
        "## Setup Boundary Conditions\n",
        "\n",
        "The rectangular domain has 4 tagged boundaries, left, top, right and bottom. We need to set boundary conditons for each of these tagged boundaries. We can set Transmissive type BC on the outflow boundaries and reflective on the others. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "rF7q2bzOn64K",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "Br = anuga.Reflective_boundary(domain)\n",
        "Bt = anuga.Transmissive_boundary(domain)\n",
        "\n",
        "domain.set_boundary({'bottom':   Br,\n",
        "                     'right':    Bt, # outflow\n",
        "                     'top':      Bt, # outflow\n",
        "                     'left':     Br})"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kSoqioIwYQcX",
        "colab_type": "text"
      },
      "source": [
        "## Setup Inflow \n",
        "\n",
        "We need some water to flow. The easiest way to input a specified amount of water is via an `Inlet_operator` where we can specify a discharge `Q`."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "KxLx34vhZIRC",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Setup inlet flow\n",
        "center = (382270.0, 6354285.0)\n",
        "radius = 10.0\n",
        "region0 = anuga.Region(domain, center=center, radius=radius)\n",
        "fixed_inflow = anuga.Inlet_operator(domain, region0 , Q=19.7)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YRJs8HbRwLfM",
        "colab_type": "text"
      },
      "source": [
        "## Run the Evolution\n",
        "\n",
        "We evolve using a `for` statement, which evolves the quantities using the shallow water wave solver. The calculation `yields` every `yieldstep` seconds, up to a given `duration`."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "3yxIDXQiH3XN",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "for t in domain.evolve(yieldstep=20, duration=300):\n",
        "  \n",
        "    #dplotter.plot_depth_frame()\n",
        "    dplotter.save_depth_frame(vmin=0.0, vmax=1.0)\n",
        "    \n",
        "    domain.print_timestepping_statistics()\n",
        "\n",
        "    \n",
        "# Read in the png files stored during the evolve loop\n",
        "dplotter.make_depth_animation() "
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "iYqp9_LkU2Ve",
        "colab_type": "text"
      },
      "source": [
        "## SWW File\n",
        "\n",
        "The evolve loop saves the quantites at the end of each yield step to an `sww` file, with name domain name + extension sww. In this case the sww file is `merewether1.sww`. \n",
        "\n",
        "An sww file can be viewed via our 3D anuga-viewer application, via the crayfish plugin for QGIS, or simply read back into python using netcdf commands. \n",
        "\n",
        "For this clinic we have provided a wrapper called an SWW_plotter to provide easy acces to the saved quantities, stage, elev, depth, xmom, ymom, xvel, yvel, speed which are all time slices of centroid values, and a time variable."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "FnkApd1GPe7t",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Create a wrapper for contents of sww file\n",
        "swwfile = 'merewether1.sww'\n",
        "splotter = anuga.SWW_plotter(swwfile)\n",
        "\n",
        "\n",
        "# Plot Depth and Speed at the last time slice\n",
        "plt.subplot(121)\n",
        "splotter.triang.set_mask(None)\n",
        "plt.tripcolor(splotter.triang, \n",
        "              facecolors = splotter.depth[-1,:], \n",
        "              cmap='viridis')\n",
        "\n",
        "plt.title(\"Depth\")\n",
        "\n",
        "\n",
        "plt.subplot(122)\n",
        "splotter.triang.set_mask(None)\n",
        "plt.tripcolor(splotter.triang, \n",
        "              facecolors = splotter.speed[-1,:], \n",
        "              cmap='viridis')\n",
        "\n",
        "plt.title(\"Speed\");\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YQIzXk6PYvgd",
        "colab_type": "text"
      },
      "source": [
        "### Comparison\n",
        "\n",
        "The data file ObservationPoints.csv contains some comparison depth data. Let's plot the depth for our simulation against the comparison data. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "3_v8qFMxWWZq",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "point_observations = np.genfromtxt(\n",
        "    os.path.join(data_dir,'ObservationPoints.csv'),\n",
        "    delimiter=\",\",skip_header=1)\n",
        "\n",
        "# Convert to absolute corrdinates\n",
        "xc = splotter.xc + splotter.xllcorner\n",
        "yc = splotter.yc + splotter.yllcorner\n",
        "\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",
        "loc_id = point_observations[:,2]\n",
        "\n",
        "fig, ax = plt.subplots()\n",
        "ax.plot(loc_id, point_observations[:,4], '*r', label='ARR')\n",
        "ax.plot(loc_id, point_observations[:,5], '*b', label='Tuflow')\n",
        "ax.plot(loc_id, splotter.stage[-1,nearest_points], '*g', label='Anuga')\n",
        "\n",
        "plt.xticks(range(0,5))\n",
        "plt.xlabel('ID')\n",
        "plt.ylabel('Stage')\n",
        "ax.legend()\n",
        "\n",
        "plt.show()\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rlT-HfB0dt0R",
        "colab_type": "text"
      },
      "source": [
        "## Flow with Houses\n",
        "\n",
        "We have polygonal data which specifies the location of a number of structures (homes) in our study. We can consider the flow in which those houses are cut out of the simulation. \n",
        "\n",
        "First we read in the house polygonal data. To maintain a small mesh size we will only read in structures with an area grester than 60 m^2. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "5UJf_z9JAhsZ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Read in house polygons from data directory and retain those of area > 60 m^2\n",
        "\n",
        "import glob\n",
        "house_files = glob.glob(os.path.join(data_dir,'house*.csv'))\n",
        "\n",
        "house_polygons = []\n",
        "for hf in house_files:\n",
        "  house_poly = anuga.read_polygon(hf)\n",
        "  poly_area = anuga.polygon_area(house_poly)\n",
        "  \n",
        "  # Leave out some small houses\n",
        "  if poly_area > 60:\n",
        "    house_polygons.append(house_poly)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1HaKkAcgiRc5",
        "colab_type": "text"
      },
      "source": [
        "### Create Domain \n",
        "\n",
        "To incorporate the housing information, we will cutout the polygons representing the houses. This is done by passing the list of house polygons to the `interior_holes` argument of the `anuga.create_domain_from_regions` procedure. \n",
        "\n",
        "This will produce a new tagged boundary region called `interior`. We will have to assign a boundsry condition to this new boundary region."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "MxM92-6iiXYD",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Resolution for most of the mesh\n",
        "base_resolution = 20.0  # m^2\n",
        "\n",
        "# Resolution in particular area of interest \n",
        "merewether_resolution = 10.0 # m^2\n",
        "\n",
        "domain = anuga.create_domain_from_regions(\n",
        "            bounding_polygon,\n",
        "            boundary_tags={'bottom': [0],\n",
        "                           'right':  [1],\n",
        "                           'top':    [2],\n",
        "                           'left':   [3]},\n",
        "            maximum_triangle_area=base_resolution,\n",
        "            interior_holes=house_polygons,\n",
        "            interior_regions=interior_regions)\n",
        "\n",
        "\n",
        "domain.set_name('merewether2') # Name of sww file\n",
        "domain.set_low_froude(1)\n",
        "\n",
        "# Setup Initial Conditions\n",
        "domain.set_quantity('elevation', filename=topography_file, location='centroids') # Use function for elevation\n",
        "domain.set_quantity('friction', 0.01, location='centroids')                        # Constant friction \n",
        "domain.set_quantity('stage', expression='elevation', location='centroids')         # Dry Bed \n",
        "\n",
        "# Setup BC\n",
        "Br = anuga.Reflective_boundary(domain)\n",
        "Bt = anuga.Transmissive_boundary(domain)\n",
        "\n",
        "\n",
        "# NOTE: We need to assign a BC to the interior boundary region.\n",
        "domain.set_boundary({'bottom':   Br,\n",
        "                     'right':    Bt, # outflow\n",
        "                     'top':      Bt, # outflow\n",
        "                     'left':     Br,\n",
        "                     'interior': Br})\n",
        "\n",
        "\n",
        "# Setup inlet flow\n",
        "center = (382270.0, 6354285.0)\n",
        "radius = 10.0\n",
        "region0 = anuga.Region(domain, center=center, radius=radius)\n",
        "fixed_inflow = anuga.Inlet_operator(domain, region0 , Q=19.7)\n",
        "\n",
        "\n",
        "dplotter = anuga.Domain_plotter(domain)  \n",
        "plt.triplot(dplotter.triang, linewidth = 0.4);"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vmI7UTWAfcCL",
        "colab_type": "text"
      },
      "source": [
        "### Evolve"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Ktd2SqLL46WR",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "for t in domain.evolve(yieldstep=20, duration=300):\n",
        "  \n",
        "    #dplotter.plot_depth_frame()\n",
        "    dplotter.save_depth_frame(vmin=0.0, vmax=1.0)\n",
        "    \n",
        "    domain.print_timestepping_statistics()\n",
        "\n",
        "    \n",
        "# Read in the png files stored during the evolve loop\n",
        "dplotter.make_depth_animation() "
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Jhc2maOJlSRf",
        "colab_type": "text"
      },
      "source": [
        "### Read in SWW File and Compare\n",
        "\n",
        "Perhaps not conclusive, but with the houses the `anuga` results, especially for id point 0, are much closer to the comparison  results. Note that we are running with a very coarse mesh for this case study. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "PeZ0LxtkJPLR",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Create a wrapper for contents of sww file\n",
        "swwfile2 = 'merewether2.sww'\n",
        "splotter2 = anuga.SWW_plotter(swwfile2)\n",
        "\n",
        "# Convert to absolute corrdinates\n",
        "xc = splotter2.xc + splotter2.xllcorner\n",
        "yc = splotter2.yc + splotter2.yllcorner\n",
        "\n",
        "nearest_points_2 = []\n",
        "for row in point_observations:\n",
        "    nearest_points_2.append(np.argmin( (xc-row[0])**2 + (yc-row[1])**2 ))\n",
        " \n",
        "loc_id = point_observations[:,2]\n",
        "\n",
        "fig, ax = plt.subplots()\n",
        "ax.plot(loc_id, point_observations[:,4], '*r', label='ARR')\n",
        "ax.plot(loc_id, point_observations[:,5], '*b', label='Tuflow')\n",
        "ax.plot(loc_id, splotter2.stage[-1,nearest_points_2], '*g', label='Anuga1')\n",
        "ax.plot(loc_id, splotter.stage[-1,nearest_points], '*k', label='Anuga0')\n",
        "\n",
        "\n",
        "plt.xticks(range(0,5))\n",
        "plt.xlabel('ID')\n",
        "plt.ylabel('Stage')\n",
        "ax.legend()\n",
        "\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "M8FvRSxIfvtS",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "U58VBVfFnaRZ",
        "colab_type": "text"
      },
      "source": [
        "[Next Notebook](https://colab.research.google.com/github/stoiver/anuga-clinic-2018/blob/master/notebooks/notebook3.ipynb)"
      ]
    }
  ]
}