notebooks/notebook5.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "notebooks/notebook5.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/notebook5.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 5: Monai Valley Tsunami runup \n",
        "\n",
        "Validation of the AnuGA implementation of the shallow water wave equation.\n",
        "This script sets up Okushiri Island benchmark as published at the\n",
        "\n",
        "The Third International Workshop on Long-Wave Runup Models\n",
        "June 17-18 2004\n",
        "\n",
        "Wrigley Marine Science Center\n",
        "Catalina Island, California\n",
        "http://www.cee.cornell.edu/longwave/\n",
        "\n",
        "\n",
        "The validation data was downloaded and made available in this directory\n",
        "for convenience but the original data is available at\n",
        "http://www.cee.cornell.edu/longwave/index.cfm?page=benchmark&problem=2\n",
        "where a detailed description of the problem is also available.\n",
        "\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",
        "If on github, first follow the link `View in Colaboratory' to start running on google's colab environment. \n",
        "\n",
        "Then 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",
        "The `anuga` and `anuga_clinic-2018` packages 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",
        "\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/validation_tests/experimental_data/okishuri/run_okishuri.py`.\n",
        "\n",
        "First we define the location of our data files, which specify:\n",
        "\n",
        "   (1) the incoming tsunami wave, \n",
        "   \n",
        "   (2) the bathymetry file, \n",
        "   \n",
        "   (3) the measured stage height at 3 gauge locations.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "53vsHGNPAJ36",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Incoming boundary wave\n",
        "boundary_filename = '/content/anuga-clinic-2018/data/Okushiri/Benchmark_2_input.txt'\n",
        "\n",
        "# Digital Elevation Model\n",
        "bathymetry_filename = '/content/anuga-clinic-2018/data/Okushiri/Benchmark_2_Bathymetry.xya'\n",
        "\n",
        "# Observed timeseries\n",
        "gauge_filename = '/content/anuga-clinic-2018/data/Okushiri/output_ch5-7-9.txt'\n",
        "\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "38DaBV-eBJ1x",
        "colab_type": "text"
      },
      "source": [
        "### Load Barthymetry Data\n",
        "\n",
        "Using in the barthymetry data provided from the workshop. We need to reshape the data and form a raster (x,y,Z)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "03R3vG2gbuBQ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "xya = np.loadtxt(bathymetry_filename, skiprows=1, delimiter=',')\n",
        "  \n",
        "X = xya[:,0].reshape(393,244)\n",
        "Y = xya[:,1].reshape(393,244)\n",
        "Z = xya[:,2].reshape(393,244)\n",
        "\n",
        "\n",
        "plt.contourf(X,Y,Z, 20, cmap=plt.get_cmap('gist_earth'));\n",
        "plt.title('Barthymetry')\n",
        "plt.colorbar();\n",
        "\n",
        "# Create raster tuple\n",
        "x = X[:,0]\n",
        "y = Y[0,:]\n",
        "Zr = np.flipud(Z.T)\n",
        "\n",
        "raster = x,y,Zr"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uFqsJONpKLJX",
        "colab_type": "text"
      },
      "source": [
        "### Load Incoming Wave \n",
        "\n",
        "We will apply an incoming wave on the left boundary. So first we load the data from the `boundary_filename` file. \n",
        "\n",
        "From the data we form an interpolation functon called  `wave_function`, which will be used to specify the boundary condition on the left. \n",
        "\n",
        "And we also plot the function.  The units of the data in the file are metres, and the scale of the experimental setup is 1 in 400. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "TrxYp6Xbi5zN",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "bdry = np.loadtxt(boundary_filename, skiprows=1)\n",
        "\n",
        "bdry_t = bdry[:,0]\n",
        "bdry_v = bdry[:,1]\n",
        "\n",
        "\n",
        "import scipy\n",
        "wave_function = scipy.interpolate.interp1d(bdry_t, bdry_v, kind='zero', fill_value='extrapolate')\n",
        "\n",
        "t = np.linspace(0.0,25.0, 100)\n",
        "\n",
        "plt.plot(t,wave_function(t)*400);\n",
        "plt.xlabel('Seconds')\n",
        "plt.ylabel('Metres')\n",
        "plt.title('Incoming Wave');"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "m9Jj58FyL7oe",
        "colab_type": "text"
      },
      "source": [
        "## Setup Domain\n",
        "\n",
        "We define the `domain` for our simulation. This object encapsulates the mesh for our problem, which is defined by setting up a bounding polygon and associated tagged boundary. We use the `base_resolution` variable to set the maximum area of the triangles of our mesh.  \n",
        "\n",
        "At the end we use `matplotlib` to visualise the mesh associated with the `domain`.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "m7kWU_05BT-Q",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "base_resolution = 0.01\n",
        "#base_resolution = 0.0005\n",
        "\n",
        "# Basic geometry and bounding polygon\n",
        "xleft   = 0\n",
        "xright  = 5.448\n",
        "ybottom = 0\n",
        "ytop    = 3.402\n",
        "\n",
        "point_sw = [xleft, ybottom]\n",
        "point_se = [xright, ybottom]\n",
        "point_nw = [xleft, ytop]    \n",
        "point_ne = [xright, ytop]\n",
        "\n",
        "bounding_polygon = [point_se,\n",
        "                    point_ne,\n",
        "                    point_nw,\n",
        "                    point_sw]\n",
        "\n",
        "domain = anuga.create_domain_from_regions(bounding_polygon,\n",
        "                                 boundary_tags={'wall': [0, 1, 3],\n",
        "                                                'wave': [2]},     \n",
        "                                 maximum_triangle_area=base_resolution,\n",
        "                                 use_cache=False,\n",
        "                                 verbose=False)\n",
        "\n",
        "\n",
        "domain.set_name('okushiri')  # Name of output sww file \n",
        "domain.set_minimum_storable_height(0.001) # Don't store w-z < 0.001m\n",
        "domain.set_flow_algorithm('DE1')\n",
        "\n",
        "print 'Number of Elements ', domain.number_of_elements\n",
        "\n",
        "dplotter = anuga.Domain_plotter(domain, min_depth=0.001)  \n",
        "plt.triplot(dplotter.triang, linewidth = 0.4);"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "z6kDzyKSjnlA",
        "colab_type": "text"
      },
      "source": [
        "### Setup Quantities\n",
        "\n",
        "We use the `raster` created earlier to set the `quantity` called `elevation`.  We also set the  `stage` and the Mannings `friction`. \n",
        "\n",
        "We also visualise the `elevation` quantity. \n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "HqXswJ3EjwDi",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "domain.set_quantity('elevation',raster=raster, location='centroids')\n",
        "domain.set_quantity('stage', 0.0)\n",
        "domain.set_quantity('friction', 0.0025)\n",
        "\n",
        "\n",
        "plt.tripcolor(dplotter.triang, \n",
        "              facecolors = dplotter.elev, \n",
        "              edgecolors='k', \n",
        "              cmap='gist_earth')\n",
        "plt.colorbar();"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9BE9r3AEcSbi",
        "colab_type": "text"
      },
      "source": [
        "### Setup Boundary Conditions\n",
        "\n",
        "Excuse the verbose boundary type name `Transmissive_n_momentum_zero_t_momentum_set_stage_boundary`, but we use that to set the incoming wave boundary condition. \n",
        "\n",
        "On the other boundaries we will have just reflective boundaries."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "1l0-N381KKGO",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "Bts = anuga.Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(domain, wave_function)\n",
        "\n",
        "Br = anuga.Reflective_boundary(domain)\n",
        "\n",
        "domain.set_boundary({'wave': Bts, 'wall': Br})"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rO1w6g4pSagt",
        "colab_type": "text"
      },
      "source": [
        "### Setup Interogation Variables\n",
        "\n",
        "We will record the stage at the 3 gauge locations and at the Monai valley. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "XzncoR3CsQim",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "yieldstep = 0.5\n",
        "finaltime = 22.5\n",
        "\n",
        "nt = int(finaltime/yieldstep)+1\n",
        "\n",
        "\n",
        "# area for gulleys\n",
        "x1 = 4.85\n",
        "x2 = 5.25\n",
        "y1 = 2.05\n",
        "y2 = 1.85\n",
        "\n",
        "# indices in gulley area\n",
        "x = domain.centroid_coordinates[:,0]\n",
        "y = domain.centroid_coordinates[:,1]\n",
        "v = np.sqrt( (x-x1)**2 + (y-y1)**2 ) + np.sqrt( (x-x2)**2 + (y-y2)**2 ) < 0.5\n",
        "\n",
        "\n",
        "# Gauge and bounday locations\n",
        "gauge = [[4.521, 1.196],  [4.521, 1.696],  [4.521, 2.196]] #Ch 5-7-9\n",
        "bdyloc = [0.00001, 2.5]\n",
        "g5_id = domain.get_triangle_containing_point(gauge[0])\n",
        "g7_id = domain.get_triangle_containing_point(gauge[1])\n",
        "g9_id = domain.get_triangle_containing_point(gauge[2])\n",
        "bc_id = domain.get_triangle_containing_point(bdyloc)\n",
        "\n",
        "# Arrays to store data\n",
        "meanstage = np.nan*np.ones((nt,))\n",
        "g5 = np.nan*np.ones((nt,))\n",
        "g7 = np.nan*np.ones((nt,))\n",
        "g9 = np.nan*np.ones((nt,))\n",
        "bc = np.nan*np.ones((nt,))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "B6YcfllrTqV8",
        "colab_type": "text"
      },
      "source": [
        "### Evolve"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "RQxmxRGzTaaV",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "Stage = domain.quantities['stage'].centroid_values\n",
        "stage_qoi = Stage[v]\n",
        "\n",
        "k = 0\n",
        "for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):\n",
        "\n",
        "    domain.print_timestepping_statistics()\n",
        "    \n",
        "    # stage\n",
        "    stage_qoi = Stage[v]\n",
        "    meanstage[k] = np.mean(stage_qoi)\n",
        "    g5[k] = Stage[g5_id]\n",
        "    g7[k] = Stage[g7_id]\n",
        "    g9[k] = Stage[g9_id]\n",
        "    \n",
        "    k = k+1\n",
        "     \n",
        "    #dplotter.save_depth_frame()\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": "J5iK5WGl6J-w",
        "colab_type": "text"
      },
      "source": [
        "### Animation Using swwfile\n",
        "\n",
        "Read in the sww file and then iterate through the time slices to produce an animation."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "qXxdE6vgi-2-",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "swwplotter = anuga.SWW_plotter('okushiri.sww', min_depth = 0.001)\n",
        "\n",
        "n = len(swwplotter.time)\n",
        "\n",
        "for k in range(n):\n",
        "  if k%10 == 0:\n",
        "     print ' '\n",
        "  swwplotter.save_stage_frame(frame=k, vmin=-0.02, vmax = 0.1)\n",
        "  print '(', swwplotter.time[k], k,')',\n",
        "  \n",
        "\n",
        "print ' ' \n",
        "swwplotter.make_stage_animation()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VCat3cQwIcsL",
        "colab_type": "text"
      },
      "source": [
        "### View Time  Series"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "AjtOijEcIiD1",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "old_figsize = plt.rcParams['figure.figsize']\n",
        "\n",
        "plt.rcParams['figure.figsize'] = [12, 5]\n",
        "\n",
        "gauge = np.loadtxt(gauge_filename, skiprows=1)\n",
        "\n",
        "gauge_t = gauge[:,0]\n",
        "gauge_5 = gauge[:,1]\n",
        "gauge_7 = gauge[:,2]\n",
        "gauge_9 = gauge[:,3]\n",
        "\n",
        "nt = int(finaltime/yieldstep)+1\n",
        "\n",
        "import scipy\n",
        "gauge_5_f = scipy.interpolate.interp1d(gauge_t, gauge_5, kind='zero', fill_value='extrapolate')\n",
        "gauge_7_f = scipy.interpolate.interp1d(gauge_t, gauge_7, kind='zero', fill_value='extrapolate')\n",
        "gauge_9_f = scipy.interpolate.interp1d(gauge_t, gauge_9, kind='zero', fill_value='extrapolate')\n",
        "\n",
        "t = np.linspace(0.0,22.5, 451)\n",
        "\n",
        "tt= np.linspace(0.0,22.5, nt)\n",
        "\n",
        "plt.subplot(1,4,1)\n",
        "plt.plot(t,gauge_5_f(t)*4)\n",
        "plt.plot(tt,g5*400)\n",
        "plt.title('Gauge 5')\n",
        "\n",
        "plt.subplot(1,4,2)\n",
        "plt.plot(t,gauge_7_f(t)*4)\n",
        "plt.plot(tt,g7*400)\n",
        "plt.title('Gauge 7')\n",
        "\n",
        "plt.subplot(1,4,3)\n",
        "plt.plot(t,gauge_9_f(t)*4)\n",
        "plt.plot(tt,g9*400)\n",
        "plt.title('Gauge 9')\n",
        "         \n",
        "plt.subplot(1,4,4)\n",
        "plt.plot(tt,meanstage*400)\n",
        "plt.title('Runup');\n",
        "         \n",
        "\n",
        "plt.rcParams['figure.figsize'] = old_figsize"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "X_w7J8IUbMUA",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": 0,
      "outputs": []
    }
  ]
}