VisualizingOutput-WaveMoL.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualizing Output Data -- WaveMoL\n",
    "---\n",
    "\n",
    "Cactus output files can be confusing to a new user. The aim of this tutorial is to provide you with hands-on practice working with output files and visualizing their data. To do this, we will run a simple simulation and then visualize its output.\n",
    "\n",
    "**Simulation:** we will use the WaveMol example as our test case. The source code is located in `Cactus/arrangements/CactusExamples/WaveMoL`. We will not be writing any source code, creating new thorns, re-compiling Cactus, or anything like that. We will simply modify the WaveMoL parfile (below) so that the simulation generates examples of many of the kinds of output files that you will encounter. If you have a working installation of the Einstein Toolkit, and if you have the necessary Python libraries installed, everything in this notebook should \"just work\".\n",
    "\n",
    "**Visualization:** we will use Python's matplotlib graphics library to visualize our simulation data. This is not necessarily the best choice, but it is convenient, especially for a Jupyter notebook-based tutorial such as this. At the very least, visualization in Python can be a useful preliminary diagnostic, before moving to a more advanced tool such as VisIt.\n",
    "\n",
    "**Python Libraries:** parts of this notebook require NumPy, Matplotlib, and H5Py. Verify that these are installed on your machine before proceding.\n",
    "\n",
    "**Further Reading:** there is a lot of excellent I/O documentation for Cactus. For example, see the documentation for the \n",
    "[IOUtil](http://einsteintoolkit.org/thornguide/CactusBase/IOUtil/documentation.html), \n",
    "[IOBasic](http://einsteintoolkit.org/thornguide/CactusBase/IOBasic/documentation.html), \n",
    "[IOASCII](http://einsteintoolkit.org/thornguide/CactusBase/IOASCII/documentation.html), \n",
    "and \n",
    "[IOHDF5](http://einsteintoolkit.org/thornguide/CactusPUGHIO/IOHDF5/documentation.html)\n",
    "thorns. A condensed version of that material is presented here.\n",
    "\n",
    "**What this tutorial is NOT:** there are two important topics that this tutorial does not get into.\n",
    "\n",
    "1. Carpet output. The WaveMoL thorn uses a uniform grid and the corresponding PUGH driver. It is very likely that the applications you are interested in will require the Carpet driver instead of PUGH, which allows for mesh refinement. The output files generated by Carpet are slightly more complicated. Unfortunately this tutorial will not examine those output files.\n",
    "\n",
    "\n",
    "2. Output from multiple processors. This tutorial is meant as an introduction to Cactus output. The imagined audience is someone working on a laptop or desktop workstation. As such, it is assumed here that you are running this notebook on a single processor, not a cluster. How to handle output from multiple processors adds a layer of complexity that we will not get into here.\n",
    "\n",
    "Hopefully both of these important topics will be addressed in a future tutorial."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0. Preliminaries\n",
    "\n",
    "First, specify the paths on your machine to the Cactus directory, simulations directory, and the directory where you want to store the edited copy of the WaveMoL parfile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Linux example:\n",
    "#cactus_dir = '/home/ejwest/ETK/Cactus'\n",
    "#sims_dir   = '/home/ejwest/ETK/simulations'\n",
    "#par_dir    = '/home/ejwest/ETK/parfiles'\n",
    "\n",
    "## Windows example:\n",
    "#cactus_dir = 'C:\\\\Users\\\\EJWest\\\\ETK\\\\Cactus'\n",
    "#sims_dir   = 'C:\\\\Users\\\\EJWest\\\\ETK\\\\Simulations'\n",
    "#par_dir    = 'C:\\\\Users\\\\EJWest\\\\ETK\\\\Parfiles'\n",
    "\n",
    "## Fill in paths to directories:\n",
    "import os\n",
    "home_dir = os.environ[\"HOME\"]\n",
    "cactus_dir = os.path.join(home_dir, \"Cactus\") #<-- Path to Cactus directory\n",
    "sims_dir   = os.path.join(home_dir, \"simulations\") #<-- Path to simulations directory\n",
    "par_dir    = os.path.join(cactus_dir, \"parfiles\") #<-- Path to parfile directory\n",
    "\n",
    "os.makedirs(par_dir, exist_ok=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Edit Parfile\n",
    "\n",
    "Next, we create a modified copy of the parfile in the directory that you specified above. The parts of the file that have to do with controlling output have commented headers. For the sake of illustation, here we are generating all three kinds of output files: xgraph-style ASCII, gnuplot-style ASCII, and HDF5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile {par_dir}/wave_mol.par\n",
    "Cactus::cctk_itlast = 100\n",
    "\n",
    "ActiveThorns = \"CoordBase CartGrid3D SymBase Boundary\"\n",
    "ActiveThorns = \"PUGHReduce PUGH PUGHSlab\"\n",
    "Grid::domain = \"full\"\n",
    "Grid::type = \"byspacing\"\n",
    "Grid::avoid_origin = \"no\"\n",
    "driver::global_nx = 51\n",
    "driver::global_ny = 51\n",
    "driver::global_nz = 51\n",
    "Grid::dxyz = 0.02\n",
    "driver::ghost_size = 1\n",
    "\n",
    "ActiveThorns = \"NaNChecker\"\n",
    "ActiveThorns = \"LocalReduce\"\n",
    "\n",
    "ActiveThorns = \"MoL\"\n",
    "MethodOfLines::ode_method = \"icn\" #default: icn\n",
    "        \n",
    "ActiveThorns = \"IDWaveMoL\"\n",
    "IDWaveMoL::initial_data = \"gaussian\"\n",
    "IDWaveMoL::amplitude = 1.0 #default: 1.0\n",
    "IDWaveMoL::radius = 0.0 #default: 0\n",
    "IDWaveMoL::sigma = 0.1 #default: 0.1\n",
    "IDWaveMoL::centrex = 0.0 #default: 0\n",
    "IDWaveMoL::centrey = 0.0 #default: 0\n",
    "IDWaveMoL::centrez = 0.0 #default: 0\n",
    "IDWaveMoL::kx = 0.0 #default: 0\n",
    "IDWaveMoL::ky = 0.0 #default: 0\n",
    "IDWaveMoL::kz = 0.0 #default: 0\n",
    "IDWaveMoL::offsett = 0.0 #default: 0\n",
    "IDWaveMoL::offsetx = 0.0 #default: 0\n",
    "IDWaveMoL::offsety = 0.0 #default: 0\n",
    "IDWaveMoL::offsetz = 0.0 #default: 0\n",
    "IDWaveMoL::slopet = 0.0 #default: 0\n",
    "IDWaveMoL::slopex = 0.0 #default: 0\n",
    "IDWaveMoL::slopey = 0.0 #default: 0\n",
    "IDWaveMoL::slopez = 0.0 #default: 0\n",
    "\n",
    "ActiveThorns = \"WaveMoL\"\n",
    "WaveMoL::bound = \"radiation\"\n",
    "\n",
    "ActiveThorns = \"Time\"\n",
    "Time::timestep_method = \"courant_static\"\n",
    "Time::dtfac = 0.5 #default: 0.5\n",
    "\n",
    "##--------------------------------------------------------##\n",
    "## IOUtil                                                 ##\n",
    "##--------------------------------------------------------##\n",
    "ActiveThorns = \"IOUtil\"\n",
    "IO::out_mode = \"proc\"\n",
    "IO::out_dir = $parfile\n",
    "\n",
    "##--------------------------------------------------------##\n",
    "## IOBasic - provides I/O methods for outputting scalar   ##\n",
    "## values in ASCII format into files and for printing     ## \n",
    "## them as runtime information to screen                  ##\n",
    "##--------------------------------------------------------##\n",
    "ActiveThorns = \"IOBasic\"\n",
    "\n",
    "## stdout\n",
    "IOBasic::outInfo_every = 1\n",
    "IOBasic::outInfo_vars  = \"\n",
    "        WaveMoL::phi \n",
    "        #WaveMoL::energy\n",
    "\"\n",
    "IOBasic::outInfo_reductions = \"minimum maximum\"\n",
    "\n",
    "## ASCII output for scalars\n",
    "IOBasic::outScalar_every = 1\n",
    "IOBasic::outScalar_vars  = \"\n",
    "        WaveMoL::phi\n",
    "        #WaveMoL::energy\n",
    "\"\n",
    "IOBasic::outScalar_style = \"xgraph\" #default\n",
    "#IOBasic::outScalar_style = \"gnuplot\"\n",
    "\n",
    "##--------------------------------------------------------##\n",
    "## IOASCII - provides I/O methods for 1D, 2D, and 3D      ## \n",
    "## output of grid arrays and grid functions into files in ##\n",
    "## ASCII format. The precise format is designed for       ##\n",
    "## visualisation using the clients xgraph or gnuplot.     ##\n",
    "##--------------------------------------------------------##\n",
    "ActiveThorns = \"IOASCII\"\n",
    "\n",
    "## ASCII output for 1D slices\n",
    "IOASCII::out1D_every = 2\n",
    "IOASCII::out1D_x     = yes\n",
    "IOASCII::out1D_y     = yes\n",
    "IOASCII::out1D_z     = yes\n",
    "IOASCII::out1D_d     = yes       \n",
    "IOASCII::out1D_vars  = \"\n",
    "        WaveMoL::phi\n",
    "        #WaveMoL::energy\n",
    "\"\n",
    "#IOASCII::out1D_style = \"xgraph\" #default\n",
    "#IOASCII::out1D_style = \"gnuplot f(x)\"\n",
    "IOASCII::out1D_style = \"gnuplot f(t,x)\"\n",
    "\n",
    "## ASCII output for 2D slices\n",
    "IOASCII::out2D_every = 2\n",
    "IOASCII::out2D_xyplane_z = 0.0\n",
    "IOASCII::out2D_yzplane_x = 0.0\n",
    "IOASCII::out2D_xzplane_y = 0.0\n",
    "IOASCII::out2D_vars  = \"\n",
    "        WaveMoL::phi\n",
    "        #WaveMoL::energy\n",
    "\"\n",
    "#IOASCII::out2D_style = \"gnuplot f(x,y)\" #default\n",
    "IOASCII::out2D_style = \"gnuplot f(t,x,y)\"\n",
    "\n",
    "##--------------------------------------------------------##\n",
    "## IOHDF5 - provides an I/O method to output variables in ##\n",
    "## HDF5 file format. It also implements checkpointing and ##\n",
    "## recovery functionality using HDF5.                     ##\n",
    "##--------------------------------------------------------##\n",
    "ActiveThorns = \"IOHDF5\"\n",
    "IOHDF5::out_every            = 5\n",
    "IOHDF5::out_vars             = \"\n",
    "        WaveMoL::phi\n",
    "        #WaveMoL::energy\n",
    "        Grid::Coordinates{out_every=1000000000} #we only want this to print once\n",
    "\"\n",
    "IOHDF5::checkpoint                  = no"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Run Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%rm -rf {sims_dir}/wave_mol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%cd {cactus_dir}\n",
    "!./simfactory/bin/sim create-run wave_mol --parfile={par_dir}/wave_mol.par --procs=2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Analyze Output\n",
    "\n",
    "Now move to the output directory for the WaveMoL simulation that we just ran. (This assumes you entered the correct path at the beginning of the notebook.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%cd {sims_dir}/wave_mol/output-0000/wave_mol"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take a look at the contents of this directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%ls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are three output file format options in Cactus: xgraph-style ASCII (file extension .xg), gnuplot-style ASCII (file extension .asc), or HDF5 (file extension .h5). Notice that all three types of files are present in this example, by design. \n",
    "\n",
    "In the parfile above, we selected \"scalars\" to be written in xgraph-style ASCII format (see the `IOBasic::outScalar` lines). Scalars are quantities that have one value per time step. There are four scalar output files here:\n",
    "\n",
    "- `phi_maximum.xg` (maximum value of phi over the entire grid, at each time step)\n",
    "- `phi_minimum.xg` (minimum value of phi over the entire grid, at each time step)\n",
    "- `phi_norm1.xg` (L1 norm of phi over the entire grid, at each time step)\n",
    "- `phi_norm2.xg` (L2 norm of phi over the entire grid, at each time step)\n",
    "\n",
    "In the parfile we selected 1D ASCII output to be written in gnuplot-style format (see the `IOASCII::out1D` lines). 1D output is usually the value of some function along a 1D slice of the computational domain. There are four instances of 1D ASCI output files here:\n",
    "\n",
    "- `phi_3D_diagonal.asc` (phi along a 3D diagonal of the grid, at each time step)\n",
    "- `phi_x_[25][25].asc` (phi along the x-axis of the grid, at each time step)\n",
    "- `phi_y_[25][25].asc` (phi along the y-axis of the grid, at each time step))\n",
    "- `phi_z_[25][25].asc` (phi along the z-axis of the grid, at each time step)\n",
    "\n",
    "In the parfile we selected 2D ASCII output to be written in gnuplot-style format (see the `IOASCII::out2D` lines). 2D output is usually the value of some function along a 2D slice of the computational domain. There are three instances of 2D ASCI output files here:\n",
    "\n",
    "- `phi_xy_[25].asc` (phi along the xy-plane of the grid, at each time step)\n",
    "- `phi_xz_[25].asc` (phi along the xz-plane of the grid, at each time step)\n",
    "- `phi_yz_[25].asc` (phi along the yz-plane of the grid, at each time step)\n",
    "\n",
    "Finally, we selected HDF5 output to also be written (see the `IOHDF5` lines). There are five instances of HDF5 output files here:\n",
    "\n",
    "- `phi.h5` (phi over the full 3D grid, at each time step)\n",
    "- `r.h5` (???, at each time step)\n",
    "- `x.h5` (x coordinates, in 3D meshgrid form, at each time step)\n",
    "- `y.h5` (y coordinates, in 3D meshgrid form, at each time step)\n",
    "- `z.h5` (z coordinates, in 3D meshgrid form, at each time step)\n",
    "\n",
    "We will examine each of these types of output files in turn, starting with the simplest case of scalar output and finishing with the most difficult case of HDF5 output."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Plotting ASCII Scalar Output Data\n",
    "\n",
    "Scalar output is the simplest to understand, and so this serves as a warm-up case. Scalars usually refer to some global property over the grid at each time step. These are refered to as \"reductions\" in the Cactus documentation (because information over the entire grid is reduced to a single number at each time step). An example is the average value of a grid function, over the grid, at a given time step. Other examples are the various norms (L1, L2, etc) of a grid function, over the grid, at a given time step.\n",
    "\n",
    "Scalar output is very simple. It contains two columns: time step and the scalar value at that time step. This is true regardless of whether the output is xgraph-style ASCII (.xg) or gnuplot-style ASCII (.asc). \n",
    "\n",
    "Unless you make changes, the parfile above is setup to print xgraph-style ASCII files for scalars. You can repeat the following steps for gnuplot-style ASCII by modifying the parfile above, re-running the simulation, and then toggle the commented lines below to switch from xgraph to gnuplot style output (i.e., comment out the lines that immediately follow \"## xgraph\" and uncomment lines that immediately follow \"## gnuplot\")."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's list the xgraph-style ASCII (.xg) files that we have to choose from."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## xgraph\n",
    "%ls *.xg\n",
    "\n",
    "## gnuplot\n",
    "#%ls *.asc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Choose one of these files to work with. Here we will choose the `phi_norm2` file. To see the structure of the contents of the file, we print the first 10 lines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "### choose a file ###\n",
    "\n",
    "## xgraph\n",
    "file1 = 'phi_norm2.xg'\n",
    "\n",
    "## gnuplot\n",
    "#file1 = 'phi_norm2.asc'\n",
    "\n",
    "## print first 10 lines of file\n",
    "!head -n 10 {file1}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that xgraph-style files use quotation marks (\") to begin comments. Here we see the structure as advertised above. After a few lines of comments, the file contains two columns, with the structure\n",
    "\n",
    "- Column 1: time step\n",
    "- Column 2: scalar data\n",
    "\n",
    "Next we assign those data to a numpy array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## import numpy library\n",
    "import numpy as np\n",
    "\n",
    "## xgraph\n",
    "data1 = np.loadtxt(file1, comments='\"')\n",
    "\n",
    "## gnuplot\n",
    "#data1 = np.loadtxt(file1, comments='#')\n",
    "\n",
    "## print data to screen (optional)\n",
    "print(data1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then for the sake of code-readability, we can assign meaningful names to different parts of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## get time steps\n",
    "time_steps = data1[:,0]  #every row, first col\n",
    "\n",
    "## get scalar values\n",
    "phi_norm2  = data1[:,1]  #every row, second col\n",
    "\n",
    "## print data to screen (optional)\n",
    "print('time_steps = ', time_steps)\n",
    "print('phi_norm2 = ', phi_norm2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And finally we can plot the scalar data versus time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## set graphics backend\n",
    "%matplotlib notebook\n",
    "## import graphics libraries\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "## plot\n",
    "plt.plot(time_steps, phi_norm2, color='red', linewidth=2)\n",
    "## NOTE: you can use LaTeX in your labels by prefixing with r.\n",
    "## The r tells matplotlib to \"render\" the string as a LaTeX string.\n",
    "plt.title(r'$L_{2}$ norm of the grid function $\\phi(x,y,z)$')  \n",
    "plt.xlabel('time')\n",
    "plt.ylabel(r'$||\\phi||_{2}$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises\n",
    "\n",
    "1) Plot the maximum value of phi as a function of time.\n",
    "\n",
    "2) Modify the parfile above to print gnuplot-style ASCII output for scalars. Then plot the L2 norm of phi as a function of time.\n",
    "\n",
    "3) Modify the parfile above to print scalar output for energy. Then plot the L2 norm for energy as a function of time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Plotting ASCII 1D Output Data\n",
    "\n",
    "A 1D slice returns the value of some grid function along a ray in the computational domain. An example would be some evolved variable along the x-axis at each time step."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1D output data is slightly more complicated than the scalar case. The form of the output depends on the style that you specify in the parfile:\n",
    "\n",
    "- `IOASCII::out1D_style = \"xgraph\"`$\\rightarrow$ two columns: x, f(x)\n",
    "- `IOASCII::out1D_style = \"gnuplot f(x)\"` $\\rightarrow$ two columns: x, f(x)\n",
    "- `IOASCII::out1D_style = \"gnuplot f(t,x)\"` $\\rightarrow$ three columns: t, x, f(t,x)\n",
    "\n",
    "Unless you make changes, the parfile above is setup to print gnuplot-style ASCII files for 1D slices using the \"f(t,x)\" format."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's list the gnuplot-style ASCII (.asc) files that we have to choose from."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## xgraph\n",
    "#%ls *.xg\n",
    "\n",
    "## gnuplot\n",
    "%ls *.asc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Choose one of these files to work with. Here we will choose the `phi_x` file. To see the structure of the contents of the file, we print the first 10 lines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "### choose a file ###\n",
    "\n",
    "## xgraph\n",
    "#file2 = 'phi_x_[25][25].xg'\n",
    "\n",
    "## gnuplot\n",
    "file2 = 'phi_x_[25][25].asc'\n",
    "\n",
    "## print first 10 lines of file\n",
    "!head -n 10 {file2}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that gnuplot-style files use the hash symbol (#) to begin comments. After a few lines of comments, the file contains three columns, with the structure\n",
    "\n",
    "- Column 1: time step \n",
    "- Column 2: grid values (slice coordinate)\n",
    "- Column 3: grid function values\n",
    "\n",
    "Time steps are printed in the first column because \"gnuplot f(t,x)\" was selected in the parfile. The other options would have been to use \"gnuplot f(x)\" or \"xgraph\", in which case only two columns would appear (which appear here as columns 2 and 3).\n",
    "\n",
    "Next we assign the data to a numpy array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## import numpy library\n",
    "import numpy as np\n",
    "\n",
    "## xgraph\n",
    "#data2 = np.loadtxt(file2, comments='\"')\n",
    "\n",
    "## gnuplot\n",
    "data2 = np.loadtxt(file2, comments='#')\n",
    "\n",
    "## print data to screen (optional)\n",
    "print(data2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we would like to assign meaningful names to different parts of the data. However, this time things are slightly more complicated by the block structure of the data, which schematically looks something like this\n",
    "\n",
    "\\begin{align}\n",
    "\\text{block 1}\\left\\{\n",
    "  \\begin{array}{ccc}\n",
    "  \\tt{t1} & \\tt{x1} & \\tt{f(t1, x1)} \\\\\n",
    "  \\tt{t1} & \\tt{x2} & \\tt{f(t1, x2)} \\\\\n",
    "  \\tt{t1} & \\tt{x3} & \\tt{f(t1, x3)} \\\\\n",
    "  \\end{array}\\right.\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "\\text{block 2}\\left\\{\n",
    "  \\begin{array}{ccc}\n",
    "  \\tt{t2} & \\tt{x1} & \\tt{f(t2, x1)} \\\\\n",
    "  \\tt{t2} & \\tt{x2} & \\tt{f(t2, x2)} \\\\\n",
    "  \\tt{t2} & \\tt{x3} & \\tt{f(t2, x3)} \\\\\n",
    "  \\end{array}\\right.\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "\\text{block 3}\\left\\{\n",
    "  \\begin{array}{ccc}\n",
    "  \\tt{t3} & \\tt{x1} & \\tt{f(t3, x1)} \\\\\n",
    "  \\tt{t3} & \\tt{x2} & \\tt{f(t3, x2)} \\\\\n",
    "  \\tt{t3} & \\tt{x3} & \\tt{f(t3, x3)} \\\\\n",
    "  \\end{array}\\right.\n",
    "\\end{align}\n",
    "\n",
    "One way to handle this is as follows.\n",
    "\n",
    "1. Create a 1D array of time steps by scanning the entire first column and keeping the unique elements.\n",
    "2. Create a 1D array of spatial grid values by scanning the entire second column and keeping the unique elements. \n",
    "3. Create a 2D array of grid function values: each row corresponding to a different time step.\n",
    "\n",
    "NOTE: step 2 is a compromise of efficiency for convenience. Since we are dealing with a fixed grid, we don't need to scan the entire second column; we could just scan the second column of the first block. But it's easier to just scan the entire second column. Doing things in a more clever way is left as an exercise to the interested reader! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "### gnuplot, f(t,x) ###\n",
    "\n",
    "## identify columns\n",
    "t_col = 0\n",
    "x_col = 1\n",
    "f_col = 2\n",
    "\n",
    "## step 1: get unique time steps\n",
    "time_steps = np.unique(data2[:,t_col])\n",
    "tdim = len(time_steps)\n",
    "\n",
    "## step 2: get unique spatial grid values\n",
    "xvals = np.unique(data2[:,x_col])\n",
    "xdim = len(xvals)\n",
    "\n",
    "## step 3: get data\n",
    "phi_x = np.zeros((tdim, xdim))\n",
    "for i in range(tdim):\n",
    "    ## get time step\n",
    "    tstep = time_steps[i]\n",
    "    ## fetch corresponding data\n",
    "    tsdata = data2[data2[:,t_col]==tstep, :] \n",
    "    ## save data to arrays\n",
    "    phi_x[i,:] = tsdata[:, f_col]\n",
    "\n",
    "## print data to screen (optional)\n",
    "print('time_steps = ', time_steps)\n",
    "print('xvals =', xvals)\n",
    "print('phi_x = ', phi_x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell contains the analogous commands if you use \"gnuplot f(x)\" style output instead. (Try it!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### gnuplot, f(x) ###\n",
    "\n",
    "## identify columns\n",
    "#x_col = 0\n",
    "#f_col = 1\n",
    "\n",
    "## step 1: get unique time steps\n",
    "#import re\n",
    "#fh = open(file2).read()\n",
    "#time_steps = []\n",
    "#for line in re.findall('#Time =.*', fh):\n",
    "#    tstep = float(line.split()[2])\n",
    "#    time_steps.append(tstep)\n",
    "#f.close()\n",
    "#time_steps = np.array(time_steps)\n",
    "#time_steps\n",
    "#tdim = len(time_steps)\n",
    "\n",
    "## step 2: get unique spatial grid values\n",
    "#xvals = np.unique(data2[:,x_col])\n",
    "#xdim = len(xvals)\n",
    "\n",
    "## step 3: get data\n",
    "#phi_x = np.zeros((tdim, xdim))\n",
    "#for i in range(tdim):\n",
    "#    ## get grid value\n",
    "#    xstep = xvals[i]\n",
    "#    ## fetch corresponding data\n",
    "#    xsdata = data2[data2[:,x_col]==xstep, :] \n",
    "#    ## save data to arrays\n",
    "#    phi_x[:,i] = xsdata[:, f_col]\n",
    "\n",
    "## print data to screen (optional)\n",
    "#print('time_steps = ', time_steps)\n",
    "#print('xvals =', xvals)\n",
    "#print('phi_x = ', phi_x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are ready to plot the data in whatever way we like."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot snapshot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## set graphics backend\n",
    "%matplotlib notebook\n",
    "## import graphics libraries\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "## figure properties\n",
    "fig_width = 8\n",
    "fig_height = 4\n",
    "minval = np.min(phi_x)\n",
    "maxval = np.max(phi_x)\n",
    "\n",
    "## choose a time step\n",
    "tstep = 1\n",
    "\n",
    "## plot\n",
    "fig = plt.figure(figsize=(fig_width, fig_height))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot(xvals[:], phi_x[tstep,:], color='red', linewidth=2)\n",
    "ax.set_ylim(minval, 1.1*maxval)\n",
    "ax.set_title(r'$\\phi$ along $x$ axis')\n",
    "ax.set_xlabel('x')\n",
    "ax.set_ylabel(r'$\\phi$')\n",
    "ax.text(0.35, 0.95, 'time = %.3f' % time_steps[tstep])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot animation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## import graphics libraries\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.animation as animation\n",
    "\n",
    "## figure properties\n",
    "fig_width = 8\n",
    "fig_height = 4\n",
    "minval = np.min(phi_x)\n",
    "maxval = np.max(phi_x)\n",
    "\n",
    "### initialize plot\n",
    "plt.ioff()\n",
    "fig = plt.figure(figsize=(fig_width, fig_height))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.set_ylim(minval, 1.1*maxval)\n",
    "time = time_steps[tstep]\n",
    "ax.text(0.35, 0.95, 'time = %.3f' % time)\n",
    "scene = ax.plot(xvals[:], phi_x[0,:], color='red', linewidth=2)\n",
    "\n",
    "## update plot\n",
    "def update_scene(frame_number, xvals, phi_x, scene):\n",
    "    ax.cla()\n",
    "    ax.set_ylim(minval, 1.1*maxval)\n",
    "    time = time_steps[frame_number]\n",
    "    ax.text(0.35, 0.95, 'time = %.3f' % time)\n",
    "    scene = ax.plot(xvals[:], phi_x[frame_number,:], color='red', linewidth=2)\n",
    "\n",
    "## make animation\n",
    "fps=5\n",
    "frames=np.arange(0, len(time_steps), 1)\n",
    "anim = animation.FuncAnimation(\n",
    "    fig, update_scene, frames, fargs=(xvals, phi_x, scene), \n",
    "    interval=1000/fps, blit=False, repeat=False)\n",
    "plt.close()\n",
    "\n",
    "## play animation\n",
    "from IPython.display import HTML\n",
    "HTML(anim.to_html5_video())  #playback option 1\n",
    "#HTML(anim.to_jshtml())       #playback option 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises\n",
    "\n",
    "1) Plot phi_z as a function of time.\n",
    "\n",
    "2) Modify the parfile to use the \"gnuplot f(x)\" style for 1D output and then replot phi_x as a function of time.\n",
    "\n",
    "3) Modify the parfile to use the \"xgraph\" style for 1D output and then replot phi_x as a function of time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Plotting ASCII 2D Output Data\n",
    "\n",
    "A 2D slice returns the value of some grid function along a plane in the computational domain. An example would be some evolved variable over the xy-plane at each time step."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output for 2D slices is more complicated still. In this case there is no xgraph-style ASCII option. But as with 1D slices, there are two gnuplot-style options to choose from in the parfile:\n",
    "\n",
    "- `IOASCII::out2D_style = \"gnuplot f(x,y)\"` $\\rightarrow$ three columns: x, y, f(x,y)\n",
    "- `IOASCII::out2D_style = \"gnuplot f(t,x,y)\"` $\\rightarrow$ four columns: t, x, y, f(t,x,y)\n",
    "\n",
    "Unless you make changes, the parfile above is setup to print gnuplot-style ASCII files for 2D slices using the \"f(t,x,y)\" format."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Re-list the gnuplot-style ASCII files that we have to choose from."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%ls *.asc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Choose a file to work with. Here we will take the `phi_xy` file, which gives the grid function phi on the xy-plane, at each time step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### choose a file ###\n",
    "file3 = 'phi_xy_[25].asc'\n",
    "\n",
    "## print first 10 lines\n",
    "!head -n 10 {file3}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After a few lines of comments, the file contains four columns, with the structure\n",
    "\n",
    "- Column 1: time step \n",
    "- Column 2: grid values (first slice coordinate)\n",
    "- Column 3: grid values (second slice coordinate)\n",
    "- Column 4: grid function values\n",
    "\n",
    "Time steps are printed in the first column because \"gnuplot f(t,x,y)\" was selected in the parfile. The other option would have been to use \"gnuplot f(x,y)\", in which case only three columns would appear (which appear here as columns 2-4). (Try it!)\n",
    "\n",
    "Next we assign the data to a numpy array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## import numpy library\n",
    "import numpy as np\n",
    "\n",
    "## get data\n",
    "data3 = np.loadtxt(file3, comments='#')\n",
    "\n",
    "## print data to screen (optional)\n",
    "print(data3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we would like to assign meaningful names to different parts of the data. This time the structure schematically looks something like this\n",
    "\n",
    "\\begin{align}\n",
    "\\text{block 1}\\left\\{\n",
    "  \\begin{array}{cccc}\n",
    "  \\tt{t1} & \\tt{x1} & \\tt{y1} & \\tt{f(t1, x1, y1)} \\\\\n",
    "  \\tt{t1} & \\tt{x1} & \\tt{y2} & \\tt{f(t1, x1, y2)} \\\\\n",
    "  \\tt{t1} & \\tt{x1} & \\tt{y3} & \\tt{f(t1, x1, y3)} \\\\\n",
    "  \\end{array}\\right.\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "\\text{block 2}\\left\\{\n",
    "  \\begin{array}{cccc}\n",
    "  \\tt{t1} & \\tt{x2} & \\tt{y1} & \\tt{f(t1, x2, y1)} \\\\\n",
    "  \\tt{t1} & \\tt{x2} & \\tt{y2} & \\tt{f(t1, x2, y2)} \\\\\n",
    "  \\tt{t1} & \\tt{x2} & \\tt{y3} & \\tt{f(t1, x2, y3)} \\\\\n",
    "  \\end{array}\\right.\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "\\text{block 3}\\left\\{\n",
    "  \\begin{array}{cccc}\n",
    "  \\tt{t1} & \\tt{x3} & \\tt{y1} & \\tt{f(t1, x3, y1)} \\\\\n",
    "  \\tt{t1} & \\tt{x3} & \\tt{y2} & \\tt{f(t1, x3, y2)} \\\\\n",
    "  \\tt{t1} & \\tt{x3} & \\tt{y3} & \\tt{f(t1, x3, y3)} \\\\\n",
    "  \\end{array}\\right.\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "\\text{block 4}\\left\\{\n",
    "  \\begin{array}{cccc}\n",
    "  \\tt{t2} & \\tt{x1} & \\tt{y1} & \\tt{f(t2, x1, y1)} \\\\\n",
    "  \\tt{t2} & \\tt{x1} & \\tt{y2} & \\tt{f(t2, x1, y2)} \\\\\n",
    "  \\tt{t2} & \\tt{x1} & \\tt{y3} & \\tt{f(t2, x1, y3)} \\\\\n",
    "  \\end{array}\\right.\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "\\text{block 5}\\left\\{\n",
    "  \\begin{array}{cccc}\n",
    "  \\tt{t2} & \\tt{x2} & \\tt{y1} & \\tt{f(t2, x2, y1)} \\\\\n",
    "  \\tt{t2} & \\tt{x2} & \\tt{y2} & \\tt{f(t2, x2, y2)} \\\\\n",
    "  \\tt{t2} & \\tt{x2} & \\tt{y3} & \\tt{f(t2, x2, y3)} \\\\\n",
    "  \\end{array}\\right.\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "\\text{block 6}\\left\\{\n",
    "  \\begin{array}{cccc}\n",
    "  \\tt{t2} & \\tt{x3} & \\tt{y1} & \\tt{f(t2, x3, y1)} \\\\\n",
    "  \\tt{t2} & \\tt{x3} & \\tt{y2} & \\tt{f(t2, x3, y2)} \\\\\n",
    "  \\tt{t2} & \\tt{x3} & \\tt{y3} & \\tt{f(t2, x3, y3)} \\\\\n",
    "  \\end{array}\\right.\n",
    "\\end{align}\n",
    "\n",
    "We can handle this is as follows.\n",
    "\n",
    "1. Create a 1D array of time steps by scanning the entire first column and keeping the unique elements.\n",
    "2. Create a 1D array of spatial grid values by scanning the entire second column and keeping the unique elements.\n",
    "3. Create a 1D array of spatial grid values by scanning the entire third column and keeping the unique elements.\n",
    "4. Create a 3D array of grid function values: each row corresponding to a different time step.\n",
    "\n",
    "NOTE: as before, we are trading efficiency for convenience by scanning the entire second and third columns (steps 2 and 3). The interested reader is challenged to come up with a better approach!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "### gnuplot, f(t,x,y) ###\n",
    "\n",
    "## identify columns\n",
    "t_col = 0\n",
    "x_col = 1\n",
    "y_col = 2\n",
    "f_col = 3\n",
    "\n",
    "## step 1: get unique time steps\n",
    "time_steps = np.unique(data3[:,t_col])\n",
    "tdim = len(time_steps)\n",
    "\n",
    "## step 2: get unique spatial grid values\n",
    "xvals = np.unique(data3[:,x_col])\n",
    "xdim = len(xvals)\n",
    "\n",
    "## step 3: get unique spatial grid values\n",
    "yvals = np.unique(data3[:,y_col])\n",
    "ydim = len(yvals)\n",
    "\n",
    "## step 4: get data\n",
    "phi_xy = np.zeros((tdim, xdim, ydim))\n",
    "for n in range(tdim):\n",
    "    ## get time step\n",
    "    tstep = time_steps[n]\n",
    "    ## fetch corresponding data\n",
    "    tsdata = data3[data3[:,t_col]==tstep, :]\n",
    "    for i in range(xdim):\n",
    "        ## get grid value\n",
    "        xval = xvals[i]\n",
    "        ## fetch corresponding data\n",
    "        xsdata = tsdata[tsdata[:,x_col]==xval, :]\n",
    "        ## save data to arrays\n",
    "        phi_xy[n,i,:] = xsdata[:, f_col]\n",
    "\n",
    "## print data to screen (optional)\n",
    "print('time_steps = ', time_steps)\n",
    "print('xvals =', xvals)\n",
    "print('yvals =', yvals)\n",
    "#print('phi_xy = ', phi_xy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that the data is extracted into numpy arrays, we can visualize them in various ways. Here we will illustrate a couple of examples: we will make a color map using imshow() and a surface plot using plot_surface(). In addition, as above, we may either create a snapshot of the solution or an animation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Color map snapshot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "### plot snapshot using imshow() ###\n",
    "\n",
    "## set graphics backend\n",
    "%matplotlib notebook\n",
    "## import graphics libraries\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "\n",
    "## figure properties\n",
    "cmap = cm.gist_rainbow\n",
    "fig_width = 8\n",
    "fig_height = 4\n",
    "minval = np.min(phi_xy)\n",
    "maxval = np.max(phi_xy)\n",
    "\n",
    "## choose a time step\n",
    "tstep = 3\n",
    "\n",
    "## plot\n",
    "time = time_steps[tstep]\n",
    "fig = plt.figure(figsize=(fig_width, fig_height))\n",
    "fig.suptitle('time = %.3f' % time)\n",
    "ax = fig.add_subplot(111)\n",
    "zi = phi_xy[tstep,:,:]\n",
    "ax.imshow(zi, cmap=cmap, clim=(minval, maxval), origin='lower')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Color map animation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### animate using imshow() ###\n",
    "\n",
    "## set graphics backend\n",
    "%matplotlib notebook\n",
    "## import graphics libraries\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "import matplotlib.animation as animation\n",
    "\n",
    "## figure properties\n",
    "cmap = cm.gist_rainbow\n",
    "fig_width = 8\n",
    "fig_height = 4\n",
    "minval = np.min(phi_xy)\n",
    "maxval = np.max(phi_xy)\n",
    "\n",
    "## initialize plot\n",
    "plt.ioff()\n",
    "time = time_steps[tstep]\n",
    "fig = plt.figure(figsize=(fig_width, fig_height))\n",
    "fig.suptitle('time = %.3f' % time)\n",
    "ax = fig.add_subplot(111)\n",
    "zi = phi_xy[0,:,:]\n",
    "scene = [ax.imshow(zi, cmap=cmap, clim=(minval, maxval), origin='lower')]\n",
    "\n",
    "## update plot\n",
    "def update_plot(frame_number, time_steps, phi_xy, scene):\n",
    "    time = time_steps[frame_number]\n",
    "    fig.suptitle('time = %.3f' % time)\n",
    "    zi = phi_xy[frame_number,:,:]\n",
    "    scene[0] = ax.imshow(zi, cmap=cmap, clim=(minval, maxval), origin='lower')\n",
    "\n",
    "## make animation\n",
    "fps=5\n",
    "frames=np.arange(0, len(time_steps), 1)\n",
    "anim = animation.FuncAnimation(\n",
    "    fig, update_plot, frames, fargs=(time_steps, phi_xy, scene), \n",
    "    interval=1000/fps, blit=False, repeat=False)\n",
    "plt.close()\n",
    "\n",
    "## play animation\n",
    "from IPython.display import HTML\n",
    "HTML(anim.to_html5_video())  #playback option 1\n",
    "#HTML(anim.to_jshtml())       #playback option 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Surface plot snapshot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### plot snapshot using plot_surface() ###\n",
    "\n",
    "## set graphics backend\n",
    "%matplotlib notebook\n",
    "## import graphics libraries\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import matplotlib.cm as cm\n",
    "\n",
    "## figure properties\n",
    "cmap = cm.gist_rainbow\n",
    "fig_width = 8\n",
    "fig_height = 4\n",
    "minval = np.min(phi_xy)\n",
    "maxval = np.max(phi_xy)\n",
    "\n",
    "## choose a time step\n",
    "tstep = 1\n",
    "\n",
    "## plot\n",
    "fig = plt.figure(figsize=(fig_width, fig_height))\n",
    "time = time_steps[tstep]\n",
    "fig.suptitle('time = %.3f' % time)\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "ax.set_zlim(minval, maxval)\n",
    "xi, yi = np.meshgrid(xvals, yvals)\n",
    "zi = phi_xy[tstep,:,:]\n",
    "ax.plot_surface(xi, yi, zi, cstride=1, rstride=1, cmap=cmap, clim=(minval, maxval))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Surface plot animation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### animate using plot_surface() ###\n",
    "\n",
    "## set graphics backend\n",
    "%matplotlib notebook\n",
    "## import graphics libraries\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import matplotlib.cm as cm\n",
    "import matplotlib.animation as animation\n",
    "\n",
    "## figure properties\n",
    "cmap = cm.gist_rainbow\n",
    "fig_width = 8\n",
    "fig_height = 4\n",
    "minval = np.min(phi_xy)\n",
    "maxval = np.max(phi_xy)\n",
    "\n",
    "## initialize plot\n",
    "plt.ioff()\n",
    "fig = plt.figure(figsize=(fig_width, fig_height))\n",
    "time = time_steps[0]\n",
    "fig.suptitle('time = %.3f' % time)\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "ax.set_zlim(minval, maxval)\n",
    "xi, yi = np.meshgrid(xvals, yvals)\n",
    "zi = phi_xy[0,:,:]\n",
    "scene = [ax.plot_surface(xi, yi, zi, cstride=1, rstride=1, cmap=cmap, clim=(minval, maxval))]\n",
    "\n",
    "## update plot\n",
    "def update_plot(frame_number, time_steps, phi_xy, scene):\n",
    "    time = time_steps[frame_number]\n",
    "    fig.suptitle('time = %.3f' % time)\n",
    "    scene[0].remove()\n",
    "    zi = phi_xy[frame_number,:,:]\n",
    "    scene[0] = ax.plot_surface(xi, yi, zi, cstride=1, rstride=1, cmap=cmap, clim=(minval, maxval))\n",
    "\n",
    "## make animation\n",
    "fps=5\n",
    "frames=np.arange(0, len(time_steps), 1)\n",
    "anim = animation.FuncAnimation(\n",
    "    fig, update_plot, frames, fargs=(time_steps, phi_xy, scene), \n",
    "    interval=1000/fps, blit=False, repeat=False)\n",
    "plt.close()\n",
    "\n",
    "## play animation\n",
    "from IPython.display import HTML\n",
    "HTML(anim.to_html5_video())  #playback option 1\n",
    "#HTML(anim.to_jshtml())       #playback option 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## save animation to file\n",
    "anim_file = 'wave_mol'\n",
    "#anim.save(anim_file + '.mp4', writer='ffmpeg', fps=fps) # requires ffmpeg\n",
    "anim.save(anim_file + '.gif', writer='imagemagick', fps=fps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises\n",
    "\n",
    "1) Plot phi_yz as a function of time.\n",
    "\n",
    "2) Modify the parfile to use the \"gnuplot f(x,y)\" style for 2D output and then replot phi_xy as a function of time.\n",
    "\n",
    "3) Plot energy as a function of space and time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. Extracting Data from HDF5 Files \n",
    "\n",
    "HDF5 (Hierarchical Data Format 5) is a portable binary data format. As such, it is far more efficient to read and write than ASCII formats, and it is probably what you should use for visualizing large simulations. A nice tool for visualizing simulation data in HDF5 format is VisIt. However, the H5Py library also provides methods for interfacing with HDF5 files using Python. This is the approach we will take here. Once the data are extracted from the .h5 files, they can be visualized as normal using matplotlib.\n",
    "\n",
    "In order to extract information from an HDF5 file, it helps to have at least some familiarity with how they are structured. This section will take a brief detour into exploring the basics of HDF5 so that the remainder of the tutorial doesn't seem like \"black magic\". If you just want to get to visualization and are ok with a little black magic, you can skip to the next section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before doing anything else, let's import the H5Py library."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## import h5py library\n",
    "import h5py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now list the available HDF5 files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%ls *.h5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then choose a file to work with. Here we will examine the `phi.h5` file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## create file handle\n",
    "file4 = h5py.File('phi.h5', mode='r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "HDF5 files are not plain text files, so you can't read them using a text editor (try!). The basic elements of an HDF5 file are \"groups\" and \"datasets\". You can think of groups like directories or folders, and you can think of datasets like files. In addition, both groups and datasets have \"attributes\" that contain metadata about that object. Let's start by exploring the primary groups in our file. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## get (name, value) tuples for members of the root group\n",
    "## (the root group, called '/', is the top-level group in the file)\n",
    "list(file4.items())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we see that the root group contains two subgroups, `Cactus Parameters` and `Global Attributes`, in addition to several datasets with names of the form `WAVEMOL::phi timelevel 0 at iteration NN`. In the parfile above, we asked for HDF5 output every 5 time steps (see the `IOHDF5::out_every` line). Correspondingly, here we see that there is a dataset for every fifth time step (iteration). \n",
    "\n",
    "In general an HDF5 file could have a complicated group hierarchy with subgroups nested within subgroups, nested within subgroups, etc. Of the two groups in our file, one is empty (`Global Attributes` is listed as having zero \"members\", i.e., nothing in it) and one is not (`Cactus Parameters` is listed as having one member, which may be a dataset or another group). The following cell gives us a quick way to preview the entire hierarchy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## traverse the group hierarchy\n",
    "file4.visit(lambda x: print(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the above, we learn that the hierarchy in this case is not very deep or complicated. It should be clear that our data are in the `WAVEMOL::phi` datasets; there are no further datasets burried elsewhere in the hierarchy. The next cell extracts the first of these datasets in order to examine it further."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## dictionary-like lookup (works for datasets as well as groups)\n",
    "dset1 = file4['/WAVEMOL::phi timelevel 0 at iteration 0']\n",
    "print(dset1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next cell prints some of the main attributes of the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## get dataset attributes\n",
    "print('dataset name is', dset1.name)\n",
    "print('dataset type is', dset1.dtype)\n",
    "print('dataset shape is', dset1.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next cell lists any additional attributes (metadata) that are stored with the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## retreive (name, value) tuples of dataset attributes \n",
    "list(dset1.attrs.items())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, the next cell illustrates how to retreive data in a dataset and store it in a NumPy array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## extract data from dataset using the '[...]' syntax\n",
    "npdset1 = dset1[...]\n",
    "print(npdset1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell verifies that this resulting object is indeed a NumPy array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(npdset1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point, we treat the data object as we would treat any other NumPy array. In the next section, we use what we have learned in this section to extract our data and then visualize them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8. Plotting HDF5 2D Output Data\n",
    "\n",
    "Now we return to the task of visualizing the data in an HDF5 (.h5) file. First import H5Py and NumPy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## import libraries\n",
    "import numpy as np\n",
    "import h5py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "List the available HDF5 files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%ls *.h5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Choose a file to work with. Here we will examine the `phi.h5` file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## create file handle\n",
    "phi_h5 = h5py.File('phi.h5', mode='r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell extracts the time steps from the metadata of the HDF5 file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## get time steps\n",
    "time_steps = np.array([])\n",
    "for member in phi_h5:\n",
    "    ## find members with shape attributes\n",
    "    ## (only datasets have shape attribute)\n",
    "    if hasattr(phi_h5[member], 'shape'):\n",
    "        ## get time of dataset\n",
    "        ts = phi_h5[member].attrs['time']\n",
    "        ## append time to time_steps\n",
    "        time_steps = np.append(time_steps, ts)\n",
    "    ## sort\n",
    "    time_steps = np.sort(time_steps)\n",
    "\n",
    "## print to screen (optional)\n",
    "print(\"time_steps = \", time_steps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to plot using `plot_surface()`, we need to construct an xy-meshgrid. The grid information is not stored in the `phi.h5` file, but instead is found in `x.h5`, `y.h5`, and `z.h5`. The following cell retrieves this information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## get 3D xdata, ydata, zdata\n",
    "x_h5 = h5py.File(\"x.h5\")\n",
    "y_h5 = h5py.File(\"y.h5\")\n",
    "z_h5 = h5py.File(\"z.h5\")\n",
    "\n",
    "for member in x_h5:\n",
    "    if hasattr(x_h5[member], 'shape'):\n",
    "        xdata = x_h5[member][...]\n",
    "        ## Stop! We only need to do this for one time \n",
    "        ## step, since we are dealing with a fixed grid!\n",
    "        break \n",
    "for member in y_h5:\n",
    "    if hasattr(y_h5[member], 'shape'):\n",
    "        ydata = y_h5[member][...]\n",
    "        ## Stop! We only need to do this for one time \n",
    "        ## step, since we are dealing with a fixed grid!\n",
    "        break \n",
    "for member in z_h5:\n",
    "    if hasattr(z_h5[member], 'shape'):\n",
    "        zdata = z_h5[member][...]\n",
    "        ## Stop! We only need to do this for one time \n",
    "        ## step, since we are dealing with a fixed grid!\n",
    "        break \n",
    "\n",
    "## print to screen (optional)\n",
    "print(\"xdata = \", xdata)\n",
    "print(\"ydata = \", ydata)\n",
    "print(\"zdata = \", zdata)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next cell extracts the values of phi over the xy-grid at each time step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## get dimensions of time and space grids\n",
    "tdim = len(time_steps)\n",
    "xdim = len(xdata[0,0,:]) # indices: [z,y,x]\n",
    "ydim = len(ydata[0,:,0]) # indices: [z,y,x]\n",
    "zdim = len(ydata[:,0,0]) # indices: [z,y,x]\n",
    "\n",
    "## get 3D phidata\n",
    "phidata = np.zeros((tdim, xdim, ydim, zdim))\n",
    "for i in range(tdim):\n",
    "    ## get time step\n",
    "    tstep = time_steps[i]\n",
    "    ## find dataset with matching time step\n",
    "    for member in phi_h5:\n",
    "        if hasattr(phi_h5[member], 'shape'):\n",
    "            if phi_h5[member].attrs['time'] == tstep: \n",
    "                ## save dataset as element of numpy array\n",
    "                phidata[i,:,:,:] = phi_h5[member][...]\n",
    "\n",
    "## get zstep for xy-plane\n",
    "for i in range(zdim):\n",
    "    if zdata[i,0,0]==0:\n",
    "        z_xyplane = i\n",
    "\n",
    "## get phi_xy_data\n",
    "phi_xy = phidata[:, z_xyplane, :, :]\n",
    "\n",
    "## print to screen (optional)\n",
    "print('phi_xy =', phi_xy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next cell gets 2D meshgrid slices of xdata and ydata."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## get 2D meshgrid data\n",
    "xvals = xdata[z_xyplane,:,:]\n",
    "yvals = ydata[z_xyplane,:,:]\n",
    "\n",
    "## print to screen (optional)\n",
    "print('xvals =', xvals)\n",
    "print('yvals =', yvals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now plot exactly as before. (The cells used to make the color map plots are exactly as before, copied and pasted from above. The cells used to make the surface plots are slightly modifed, due to the fact that our `xvals` and `yvals` are already in meshgrid form. See comments below.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Color map snapshot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "### plot snapshot using imshow() ###\n",
    "\n",
    "## set graphics backend\n",
    "%matplotlib notebook\n",
    "## import graphics libraries\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "\n",
    "## figure properties\n",
    "cmap = cm.gist_rainbow\n",
    "fig_width = 8\n",
    "fig_height = 4\n",
    "minval = np.min(phi_xy)\n",
    "maxval = np.max(phi_xy)\n",
    "\n",
    "## choose a time step\n",
    "tstep = 3\n",
    "\n",
    "## plot\n",
    "time = time_steps[tstep]\n",
    "fig = plt.figure(figsize=(fig_width, fig_height))\n",
    "fig.suptitle('time = %.3f' % time)\n",
    "ax = fig.add_subplot(111)\n",
    "zi = phi_xy[tstep,:,:]\n",
    "ax.imshow(zi, cmap=cmap, clim=(minval, maxval), origin='lower')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Color map animation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### animate using imshow() ###\n",
    "\n",
    "## set graphics backend\n",
    "%matplotlib notebook\n",
    "## import graphics libraries\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "import matplotlib.animation as animation\n",
    "\n",
    "## figure properties\n",
    "cmap = cm.gist_rainbow\n",
    "fig_width = 8\n",
    "fig_height = 4\n",
    "minval = np.min(phi_xy)\n",
    "maxval = np.max(phi_xy)\n",
    "\n",
    "## initialize plot\n",
    "plt.ioff()\n",
    "time = time_steps[tstep]\n",
    "fig = plt.figure(figsize=(fig_width, fig_height))\n",
    "fig.suptitle('time = %.3f' % time)\n",
    "ax = fig.add_subplot(111)\n",
    "zi = phi_xy[0,:,:]\n",
    "scene = [ax.imshow(zi, cmap=cmap, clim=(minval, maxval), origin='lower')]\n",
    "\n",
    "## update plot\n",
    "def update_plot(frame_number, time_steps, phi_xy, scene):\n",
    "    time = time_steps[frame_number]\n",
    "    fig.suptitle('time = %.3f' % time)\n",
    "    zi = phi_xy[frame_number,:,:]\n",
    "    scene[0] = ax.imshow(zi, cmap=cmap, clim=(minval, maxval), origin='lower')\n",
    "\n",
    "## make animation\n",
    "fps=5\n",
    "frames=np.arange(0, len(time_steps), 1)\n",
    "anim = animation.FuncAnimation(\n",
    "    fig, update_plot, frames, fargs=(time_steps, phi_xy, scene), \n",
    "    interval=1000/fps, blit=False, repeat=False)\n",
    "plt.close()\n",
    "\n",
    "## play animation\n",
    "from IPython.display import HTML\n",
    "HTML(anim.to_html5_video())  #playback option 1\n",
    "#HTML(anim.to_jshtml())       #playback option 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Surface plot snapshot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### plot snapshot using plot_surface() ###\n",
    "\n",
    "## set graphics backend\n",
    "%matplotlib notebook\n",
    "## import graphics libraries\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import matplotlib.cm as cm\n",
    "\n",
    "## figure properties\n",
    "cmap = cm.gist_rainbow\n",
    "fig_width = 8\n",
    "fig_height = 4\n",
    "minval = np.min(phi_xy)\n",
    "maxval = np.max(phi_xy)\n",
    "\n",
    "## choose a time step\n",
    "tstep = 1\n",
    "\n",
    "## plot\n",
    "fig = plt.figure(figsize=(fig_width, fig_height))\n",
    "time = time_steps[tstep]\n",
    "fig.suptitle('time = %.3f' % time)\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "ax.set_zlim(minval, maxval)\n",
    "## this was the old line which no longer works here\n",
    "#xi, yi = np.meshgrid(xvals, yvals)\n",
    "## this is the new way of doing things\n",
    "xi = xvals #xvals is already in mesh grid form\n",
    "yi = yvals #yvals is already in mesh grid form\n",
    "zi = phi_xy[tstep,:,:]\n",
    "ax.plot_surface(xi, yi, zi, cstride=1, rstride=1, cmap=cmap, clim=(minval, maxval))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Surface plot animation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "### animate using plot_surface() ###\n",
    "\n",
    "## set graphics backend\n",
    "%matplotlib notebook\n",
    "## import graphics libraries\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import matplotlib.cm as cm\n",
    "import matplotlib.animation as animation\n",
    "\n",
    "## figure properties\n",
    "cmap = cm.gist_rainbow\n",
    "fig_width = 8\n",
    "fig_height = 4\n",
    "minval = np.min(phi_xy)\n",
    "maxval = np.max(phi_xy)\n",
    "\n",
    "## initialize plot\n",
    "plt.ioff()\n",
    "fig = plt.figure(figsize=(fig_width, fig_height))\n",
    "time = time_steps[0]\n",
    "fig.suptitle('time = %.3f' % time)\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "ax.set_zlim(minval, maxval)\n",
    "## this was the old line which no longer works here\n",
    "#xi, yi = np.meshgrid(xvals, yvals)\n",
    "## this is the new way of doing things\n",
    "xi = xvals #xvals is already in mesh grid form\n",
    "yi = yvals #yvals is already in mesh grid form\n",
    "zi = phi_xy[0,:,:]\n",
    "scene = [ax.plot_surface(xi, yi, zi, cstride=1, rstride=1, cmap=cmap, clim=(minval, maxval))]\n",
    "\n",
    "## update plot\n",
    "def update_plot(frame_number, time_steps, phi_xy, scene):\n",
    "    time = time_steps[frame_number]\n",
    "    fig.suptitle('time = %.3f' % time)\n",
    "    scene[0].remove()\n",
    "    zi = phi_xy[frame_number,:,:]\n",
    "    scene[0] = ax.plot_surface(xi, yi, zi, cstride=1, rstride=1, cmap=cmap, clim=(minval, maxval))\n",
    "\n",
    "## make animation\n",
    "fps=5\n",
    "frames=np.arange(0, len(time_steps), 1)\n",
    "anim = animation.FuncAnimation(\n",
    "    fig, update_plot, frames, fargs=(time_steps, phi_xy, scene), \n",
    "    interval=1000/fps, blit=False, repeat=False)\n",
    "plt.close()\n",
    "\n",
    "## play animation\n",
    "from IPython.display import HTML\n",
    "HTML(anim.to_html5_video())  #playback option 1\n",
    "#HTML(anim.to_jshtml())       #playback option 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## save animation to file\n",
    "anim_file = 'wave_mol'\n",
    "#anim.save(anim_file + '.mp4', writer='ffmpeg', fps=fps) # requires ffmpeg\n",
    "anim.save(anim_file + '.gif', writer='imagemagick', fps=fps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises\n",
    "\n",
    "1) Plot phi_yz as a function of time using the data in `phi.h5`.\n",
    "\n",
    "2) Plot phi_xz as a function of time using the data in `phi.h5`.\n",
    "\n",
    "3) Plot phi_xy on the upper/lower z=const boundary.\n",
    "\n",
    "4) Change the frequency of hdf5 output and replot phi_xy.\n",
    "\n",
    "5) Plot energy in the xy plane as a function of time."
   ]
  },
  {
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}