CreatingANewThorn-WaveEqn.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Creating a New Thorn -- WaveEqn\n",
    "---\n",
    "\n",
    "In this tutorial we will walk through the process of creating a new thorn from scratch. As a new user, you may be content at first to run simulations that have already been assembled by someone else (such as the gallery examples). Beyond that you may want to assemble your own simulation using existing thorns. However at some point you may have an application that involves a step, or calculation, or algorithm, that no existing thorn performs. At that point, you need to write your own thorn to carry out that subroutine. The case study we will consider here is to create a thorn to evolve the standard wave equation.\n",
    "\n",
    "Let's start with a brief overview of what a thorn needs to include:\n",
    "\n",
    "* a Cactus thorn must have a name.\n",
    "* a Cactus thorn must live in an arrangement.\n",
    "* a Cactus thorn must have four ccl (Cactus Configuration Language) files:\n",
    " * `interface.ccl`\n",
    " * `schedule.ccl`\n",
    " * `param.ccl`\n",
    " * `configuration.ccl`\n",
    "* a Cactus thorn must have a `src` directory\n",
    "* a Cactus thorn must have a `make.code.defn` file in that source directory\n",
    "\n",
    "The above thorn files and directories should reside in `Cactus/arrangements/ArrangementName/ThornName/`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0. Preliminaries\n",
    "\n",
    "First let's create a variable for the Cactus directory on your machine and then move to that directory (in case it is not already the working directory)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Linux example:\n",
    "#cactus_dir='/home/ejwest/ETK/Cactus'\n",
    "\n",
    "## Windows example:\n",
    "#cactus_dir='C:\\\\Users\\\\EJWest\\\\ETK\\\\Cactus'\n",
    "import os\n",
    "home_dir = os.environ[\"HOME\"]\n",
    "\n",
    "## fill in the path to your Cactus directory\n",
    "cactus_dir=os.path.join(home_dir,'Cactus') #<-- COMPLETE THIS LINE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%cd {cactus_dir}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next let's specify the name and arrangement of the thorn we want to make. Here we are creating a thorn to evolve the wave equation, so we will name it `WaveEqn` and place it in the arrangment `Tests`. As a result, the files we create below will be located in `Cactus/arrangements/Tests/WaveEqn`. The following block creates some environment variables that will be used below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Define some basic parameters describing a new thorn\n",
    "thorn_pars = {\n",
    "   \"thorn_name\" : \"WaveEqn\",\n",
    "   \"arrangement_name\" : \"Tests\",\n",
    " }\n",
    "\n",
    "## Define environment variables ARR and THORN to be used\n",
    "## throughout the rest of the tutorial\n",
    "import os\n",
    "os.environ[\"ARR\"]=thorn_pars[\"arrangement_name\"]\n",
    "os.environ[\"THORN\"]=thorn_pars[\"thorn_name\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This next command is only needed if you want to start over. Uncomment and run this line to delete the thorn you created and start over."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%rm -rf arrangements/$ARR/$THORN"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we create the thorn directory and sub-directories."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## The equivalent of mkdir -p\n",
    "def create_dir(dir):\n",
    "    print(\"Ensuring directory '\"+dir+\"'\")\n",
    "    os.makedirs(dir, exist_ok=True)\n",
    "        \n",
    "## Create the thorn directory inside the Cactus source tree\n",
    "arrangement_dir = os.path.join(cactus_dir, \"arrangements\", thorn_pars[\"arrangement_name\"])\n",
    "thorn_dir = os.path.join(arrangement_dir, thorn_pars[\"thorn_name\"])\n",
    "create_dir(thorn_dir)\n",
    "\n",
    "## Create the source directory\n",
    "src_dir = os.path.join(thorn_dir, \"src\")\n",
    "create_dir(src_dir)\n",
    "\n",
    "## Other dirs and files not strictly needed\n",
    "## for compiling and running Cactus.\n",
    "test_dir = os.path.join(thorn_dir, \"test\")\n",
    "create_dir(test_dir)\n",
    "par_dir = os.path.join(thorn_dir, \"par\")\n",
    "create_dir(par_dir)\n",
    "doc_dir = os.path.join(thorn_dir, \"doc\")\n",
    "create_dir(doc_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now move to the thorn directory that we just created and examine its contents. All that appears so far are the sub-directories that we created above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print()\n",
    "print('Directory contents:')\n",
    "%ls {thorn_dir}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're finally ready to start creating the files that make up the thorn. We can write the contents of a notebook cell to a file by using the `%%writefile` command. A command that begins with `%%` is called a \"cell-magic\", since it operates on the entire content of a cell (as opposed to a \"line-magic\", prefixed by `%`, which only operates on the same line). We will use the `%%writefile` cell-magic repeatedly throughout this tutorial to create files. Of course you could also create the files outside of the notebook, using your favorite editor. Working in the notebook is just a shortcut. \n",
    "\n",
    "As a first application of the `%%writefile` magic, let's create a README file for our thorn. Below is the standard template for that. Change the names and email addresses to your own and then execute the cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile {thorn_dir}/README\n",
    "Author(s)    : Eric J. West <ewest@d.umn.edu>\n",
    "Maintainer(s): Eric J. West <ewest@d.umn.edu>\n",
    "Licence      : BSD\n",
    "--------------------------------------------------------------------------\n",
    "\n",
    "1. Purpose\n",
    "\n",
    "Basic implementation evolving the standard 3D wave equation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now re-examine the contents of the thorn directory. The README file we just created should be listed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%ls {thorn_dir}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Source code and make file\n",
    "\n",
    "In this section we will create the contents of the `src` subdirectory. This includes any source code (written in either C++ or Fortran) and a make file."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The equation we want to solve is the standard (linear) 3D wave equation. Let's call the dependent variable $U(t,x,y,z)$. In terms of Cartesian coordinates, the wave equation is written as\n",
    "\n",
    "\\begin{align}\n",
    "  \\partial_{t}^{2}U\n",
    "  = \\partial_{x}^{2}U + \\partial_{y}^{2}U + \\partial_{z}^{2}U,\n",
    "\\end{align}\n",
    "\n",
    "where units have been chosen so that the wave speed $c=1$. In this simulation, we will use the method of lines (implemented by the MoL thorn). This requires that the equations are reduced to first-order as follows\n",
    "\n",
    "\\begin{align}\n",
    "  & \\partial_{t}U = K \\\\\n",
    "  & \\partial_{t}K \n",
    "  = \\partial_{x}^{2}U + \\partial_{y}^{2}U + \\partial_{z}^{2}U.\n",
    "\\end{align}\n",
    "\n",
    "What is passed to the MoL thorn is the right-hand side (RHS) for each variable, which here means the RHS for $U$ and $K$. Using centered finite-differences gives \n",
    "\n",
    "\\begin{align}\n",
    "  & (\\text{RHS for }U) = K \\\\\n",
    "  & \\\\\n",
    "  & (\\text{RHS for }K) \n",
    "  = \\frac{U^{n}_{i+1,j,k} + U^{n}_{i-1,j,k} - 2U^{n}_{i,j,k}}{(\\Delta{x})} \n",
    "  + \\frac{U^{n}_{i,j+1,k} + U^{n}_{i,j-1,k} - 2U^{n}_{i,j,k}}{(\\Delta{y})} \n",
    "  + \\frac{U^{n}_{i,j,k+1} + U^{n}_{i,j,k-1} - 2U^{n}_{i,j,k}}{(\\Delta{z})}\n",
    "\\end{align}\n",
    "\n",
    "For simplicity we will seek solutions with boundary conditions\n",
    "\n",
    "\\begin{align}\n",
    "  & \\partial_{t}K(t,0,y,z) = \\partial_{t}K(t,L_{x},y,z) = 0 \\\\\n",
    "  & \\partial_{t}K(t,x,0,z) = \\partial_{t}K(t,x,L_{y},z) = 0 \\\\\n",
    "  & \\partial_{t}K(t,x,y,0) = \\partial_{t}K(t,x,y,L_{z}) = 0\n",
    "\\end{align}\n",
    "\n",
    "We will take the initial configuration to be a Gaussian\n",
    "\n",
    "\\begin{align}\n",
    "  U(0,x,y,z) = e^{-(x^2 + y^2)}\n",
    "\\end{align}\n",
    "\n",
    "The equations, boundary conditions, and initial condition are implemented in the following source code files."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following file registers our variables with the MoL thorn."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile {thorn_dir}/src/mol_register.c\n",
    "#include <cctk.h>\n",
    "#include <cctk_Arguments.h>\n",
    "#include <cctk_Parameters.h>\n",
    "\n",
    "// Register variables with MoL\n",
    "void WaveEqn_MoLRegister(CCTK_ARGUMENTS)\n",
    "{\n",
    "  DECLARE_CCTK_ARGUMENTS;  // Declare all grid functions from interface.ccl\n",
    "  CCTK_INT ierr = 0, group, rhs;\n",
    "\n",
    "  // Look up the group index of the evol_group\n",
    "  group = CCTK_GroupIndex (\"WaveEqn::evol_group\");\n",
    "  \n",
    "  // Look up the group index of the rhs_group\n",
    "  rhs = CCTK_GroupIndex (\"WaveEqn::rhs_group\");\n",
    "\n",
    "  if (CCTK_IsFunctionAliased (\"MoLRegisterEvolvedGroup\"))\n",
    "  {\n",
    "    // *** KEY STEP! *** Tell MoL what we variables\n",
    "    // we are evolving and what the rhs variables are.\n",
    "    ierr += MoLRegisterEvolvedGroup(group, rhs);\n",
    "  }\n",
    "  else\n",
    "  {\n",
    "    CCTK_WARN (0, \"MoL function not aliased !\");\n",
    "    ierr++;\n",
    "  }\n",
    "\n",
    "  if (ierr)\n",
    "  {\n",
    "    CCTK_WARN(CCTK_WARN_ABORT, \"Problems registering with MoL !\");\n",
    "  }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following file initializes the solution, implementing the initial conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile {thorn_dir}/src/init.c\n",
    "#include <cctk.h>\n",
    "#include <cctk_Arguments.h>\n",
    "#include <cctk_Parameters.h>\n",
    "\n",
    "// Initialize grid functions\n",
    "void WaveEqn_Initialize(CCTK_ARGUMENTS)\n",
    "{\n",
    "  DECLARE_CCTK_ARGUMENTS;  // Declare all grid functions from interface.ccl\n",
    "  DECLARE_CCTK_PARAMETERS; // Declare all parameters from param.ccl\n",
    "\n",
    "  for (int k = 0; k < cctk_lsh[2]; k++) // loop over the z direction\n",
    "  {\n",
    "    for (int j = 0; j < cctk_lsh[1]; j++) // loop over the y direction\n",
    "    {\n",
    "      for (int i = 0; i < cctk_lsh[0]; i++) // loop over the x direction\n",
    "      {\n",
    "        const size_t ijk = CCTK_GFINDEX3D(cctkGH, i, j, k);\n",
    "        U[ijk] = exp(-x[ijk] * x[ijk] - y[ijk] * y[ijk]);\n",
    "        K[ijk] = 0.0;\n",
    "      }\n",
    "    }\n",
    "  }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following file updates the RHS's at each timestep."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile {thorn_dir}/src/evolve.c\n",
    "#include <cctk.h>\n",
    "#include <cctk_Arguments.h>\n",
    "#include <cctk_Parameters.h>\n",
    "\n",
    "// Update RHS functions\n",
    "void WaveEqn_RHS(CCTK_ARGUMENTS)\n",
    "{\n",
    "  DECLARE_CCTK_ARGUMENTS;  // Declare all grid functions from interface.ccl\n",
    "  DECLARE_CCTK_PARAMETERS; // Declare all parameters from param.ccl\n",
    "\n",
    "  const int gz = cctk_nghostzones[2];\n",
    "  const int gy = cctk_nghostzones[1];\n",
    "  const int gx = cctk_nghostzones[0];\n",
    "\n",
    "  const CCTK_REAL dt = CCTK_DELTA_TIME;\n",
    "\n",
    "  const CCTK_REAL dx = CCTK_DELTA_SPACE(0);\n",
    "  const CCTK_REAL dy = CCTK_DELTA_SPACE(1);\n",
    "  const CCTK_REAL dz = CCTK_DELTA_SPACE(2);\n",
    "\n",
    "  for (int k = gz; k < cctk_lsh[2] - gz; k++) // loop over the z direction\n",
    "  {\n",
    "    for (int j = gy; j < cctk_lsh[1] - gy; j++) // loop over the y direction\n",
    "    {\n",
    "      for (int i = gx; i < cctk_lsh[0] - gx; i++) // loop over the x direction\n",
    "      {\n",
    "        const size_t ijk = CCTK_GFINDEX3D(cctkGH, i, j, k);\n",
    "        const size_t ip1jk = CCTK_GFINDEX3D(cctkGH, i + 1, j, k);\n",
    "        const size_t im1jk = CCTK_GFINDEX3D(cctkGH, i - 1, j, k);\n",
    "        const size_t ijp1k = CCTK_GFINDEX3D(cctkGH, i, j + 1, k);\n",
    "        const size_t ijm1k = CCTK_GFINDEX3D(cctkGH, i, j - 1, k);\n",
    "        const size_t ijkp1 = CCTK_GFINDEX3D(cctkGH, i, j, k + 1);\n",
    "        const size_t ijkm1 = CCTK_GFINDEX3D(cctkGH, i, j, k - 1);\n",
    "\n",
    "        const CCTK_REAL laplacian_U = (U[ip1jk] + U[im1jk] - 2 * U[ijk]) / (dx * dx)\n",
    "                                    + (U[ijp1k] + U[ijm1k] - 2 * U[ijk]) / (dy * dy)\n",
    "                                    + (U[ijkp1] + U[ijkm1] - 2 * U[ijk]) / (dz * dz);\n",
    "\n",
    "        U_rhs[ijk] = K[ijk];\n",
    "        K_rhs[ijk] = laplacian_U;\n",
    "      }\n",
    "    }\n",
    "  }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following file imposes the boundary conditions at each timestep. This should be done after the interior points are updated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile {thorn_dir}/src/boundary.c\n",
    "#include <cctk.h>\n",
    "#include <cctk_Arguments.h>\n",
    "#include <cctk_Parameters.h>\n",
    "\n",
    "// Update grid functions at boundary points\n",
    "void WaveEqn_Boundary(CCTK_ARGUMENTS)\n",
    "{\n",
    "  DECLARE_CCTK_ARGUMENTS;  // Declare all grid functions from interface.ccl\n",
    "  DECLARE_CCTK_PARAMETERS; // Declare all parameters from param.ccl\n",
    "\n",
    "  const int gz = cctk_nghostzones[2];\n",
    "  const int gy = cctk_nghostzones[1];\n",
    "  const int gx = cctk_nghostzones[0];\n",
    "\n",
    "  const CCTK_REAL dt = CCTK_DELTA_TIME;\n",
    "\n",
    "  // Lower X Boundary\n",
    "  if (cctk_bbox[0])\n",
    "  {\n",
    "    for (int k = 0; k < cctk_lsh[2]; k++) // loop over the z direction\n",
    "    {\n",
    "      for (int j = 0; j < cctk_lsh[1]; j++) // loop over the y direction\n",
    "      {\n",
    "        for (int i = 0; i < gx; i++) // loop over the x direction\n",
    "        {\n",
    "          const size_t ijk = CCTK_GFINDEX3D(cctkGH, i, j, k);\n",
    "          U_rhs[ijk] = K[ijk];\n",
    "          K_rhs[ijk] = 0.0;\n",
    "        }\n",
    "      }\n",
    "    }\n",
    "  }\n",
    "\n",
    "  // Upper X Boundary\n",
    "  if (cctk_bbox[1])\n",
    "  {\n",
    "    for (int k = 0; k < cctk_lsh[2]; k++) // loop over the z direction\n",
    "    {\n",
    "      for (int j = 0; j < cctk_lsh[1]; j++) // loop over the y direction\n",
    "      {\n",
    "        for (int i = cctk_lsh[0] - gx; i < cctk_lsh[0]; i++) // loop over the x direction\n",
    "        {\n",
    "          const size_t ijk = CCTK_GFINDEX3D(cctkGH, i, j, k);\n",
    "          U_rhs[ijk] = K[ijk];\n",
    "          K_rhs[ijk] = 0.0;\n",
    "        }\n",
    "      }\n",
    "    }\n",
    "  }\n",
    "\n",
    "  // Lower Y Boundary\n",
    "  if (cctk_bbox[2])\n",
    "  {\n",
    "    for (int k = 0; k < cctk_lsh[2]; k++) // loop over the z direction\n",
    "    {\n",
    "      for (int j = 0; j < gy; j++) // loop over the y direction\n",
    "      {\n",
    "        for (int i = 0; i < cctk_lsh[0]; i++) // loop over the x direction\n",
    "        {\n",
    "          const size_t ijk = CCTK_GFINDEX3D(cctkGH, i, j, k);\n",
    "          U_rhs[ijk] = K[ijk];\n",
    "          K_rhs[ijk] = 0.0;\n",
    "        }\n",
    "      }\n",
    "    }\n",
    "  }\n",
    "\n",
    "  // Upper Y Boundary\n",
    "  if (cctk_bbox[3])\n",
    "  {\n",
    "    for (int k = 0; k < cctk_lsh[2]; k++) // loop over the z direction\n",
    "    {\n",
    "      for (int j = cctk_lsh[1] - gy; j < cctk_lsh[1]; j++) // loop over the y direction\n",
    "      {\n",
    "        for (int i = 0; i < cctk_lsh[0]; i++) // loop over the x direction\n",
    "        {\n",
    "          const size_t ijk = CCTK_GFINDEX3D(cctkGH, i, j, k);\n",
    "          U_rhs[ijk] = K[ijk];\n",
    "          K_rhs[ijk] = 0.0;\n",
    "        }\n",
    "      }\n",
    "    }\n",
    "  }\n",
    "\n",
    "  // Lower Z Boundary\n",
    "  if (cctk_bbox[4])\n",
    "  {\n",
    "    for (int k = 0; k < gz; k++) // loop over the z direction\n",
    "    {\n",
    "      for (int j = 0; j < cctk_lsh[1]; j++) // loop over the y direction\n",
    "      {\n",
    "        for (int i = 0; i < cctk_lsh[0]; i++) // loop over the x direction\n",
    "        {\n",
    "          const size_t ijk = CCTK_GFINDEX3D(cctkGH, i, j, k);\n",
    "          U_rhs[ijk] = K[ijk];\n",
    "          K_rhs[ijk] = 0.0;\n",
    "        }\n",
    "      }\n",
    "    }\n",
    "  }\n",
    "    \n",
    "  // Upper Z Boundary\n",
    "  if (cctk_bbox[5])\n",
    "  {\n",
    "    for (int k = cctk_lsh[2] - gz; k < cctk_lsh[2]; k++) // loop over the z direction\n",
    "    {\n",
    "      for (int j = 0; j < cctk_lsh[1]; j++) // loop over the y direction\n",
    "      {\n",
    "        for (int i = 0; i < cctk_lsh[0]; i++) // loop over the x direction\n",
    "        {\n",
    "          const size_t ijk = CCTK_GFINDEX3D(cctkGH, i, j, k);\n",
    "          U_rhs[ijk] = K[ijk];\n",
    "          K_rhs[ijk] = 0.0;\n",
    "        }\n",
    "      }\n",
    "    }\n",
    "  }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above source code contains four functions:\n",
    "\n",
    "- WaveEqn_MoLRegister(CCTK_ARGUMENTS)\n",
    "- WaveEqn_Initialize(CCTK_ARGUMENTS)\n",
    "- WaveEqn_RHS(CCTK_ARGUMENTS)\n",
    "- WaveEqn_Boundary(CCTK_ARGUMENTS)\n",
    "\n",
    "We need to tell Cactus when to execute each of these functions. We do that in `schedule.ccl`, below.\n",
    "\n",
    "The source code contains four variables that we want to track:\n",
    "\n",
    "- U\n",
    "- K\n",
    "- U_rhs\n",
    "- K_rhs\n",
    "\n",
    "We need to tell Cactus about these variables and their properties. We do that in `interface.ccl`, below. In addition, we need to allocate storage for these variables, which we do in `schedule.ccl`.\n",
    "\n",
    "As is, the source code does not contain any parameters. If it did, we would declare those in `param.ccl`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also need a make file for our thorn, which is always called `make.code.defn`. This file tells Cactus the names of the source code files to compile. The following cell provides the content for the make file. Notice in our case we have four source code files to list: `init.c`, `evolve.c`, `boundary.c`, and `mol_register.c`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile {thorn_dir}/src/make.code.defn\n",
    "## Main make.code.defn file for thorn WaveEqn\n",
    "\n",
    "## Source files in this directory\n",
    "SRCS = init.c evolve.c boundary.c mol_register.c\n",
    "\n",
    "## Subdirectories containing source files\n",
    "SUBDIRS ="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. CCL files\n",
    "\n",
    "As mentioned above, every thorn needs four ccl (Cactus Configuration Language) files. These are `interface.ccl`, `schedule.ccl`, `param.ccl`, and `configure.ccl`. Here we will create these files with explanatory comments."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the `interface.ccl` file, you specify which thorns your thorn inherits variables from, as well as declaring groups of variables that other thorns may access. Recall that we have four variables that we need to declare here: U, K, U_rhs, K_rhs. We put them in different groups, because they have different storage allocation needs. In particular, we need to keep data for U and K for three time levels, whereas we only need to keep one time level's worth of data for the RHS functions U_rhs and K_rhs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile {thorn_dir}/interface.ccl\n",
    "## Interface definitions for thorn WaveEqn\n",
    "inherits: Grid\n",
    "\n",
    "## An implementation name is required for all thorns. No\n",
    "## two thorns in a configuration can implement the same\n",
    "## interface.\n",
    "implements: WaveEqn\n",
    "\n",
    "## The groups declared below can be public, private, or protected.\n",
    "public:\n",
    "\n",
    "## A group defines a set of variables that are allocated together\n",
    "## and share common properties, i.e. TYPE, TIMELEVELS, and TAGS \n",
    "## such as the Prolongation=None tag. \n",
    "## Note that the number of timelevels can be an integer parameter.\n",
    "## The TYPE tag can take on the values GF, SCALAR, or ARRAY.\n",
    "## GF stands for \"Grid Function\" and refers to a distributed array\n",
    "## data structure.\n",
    "CCTK_REAL evol_group TYPE=GF TIMELEVELS=3 TAGS='tensortypealias=\"Scalar\"'\n",
    "{\n",
    "  U, K\n",
    "} \"Wave equation fields\"\n",
    "\n",
    "CCTK_REAL rhs_group TYPE=GF TIMELEVELS=1 TAGS='tensortypealias=\"Scalar\"'\n",
    "{\n",
    "  U_rhs, K_rhs\n",
    "} \"RHS functions for wave equation fields\"\n",
    "\n",
    "## Scalars are single variables that are available on all processors.\n",
    "#CCTK_REAL scalar_group TYPE=SCALAR \n",
    "#{\n",
    "#  scalar1, scalar2\n",
    "#}\n",
    "\n",
    "##----------------------------------------------------\n",
    "##   MoL functions to register variables\n",
    "##----------------------------------------------------\n",
    "\n",
    "CCTK_INT FUNCTION MoLRegisterEvolvedGroup(CCTK_INT IN EvolvedIndex, CCTK_INT IN RHSIndex)\n",
    "REQUIRES FUNCTION MoLRegisterEvolvedGroup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `schedule.ccl` file is where you declare storage allocations as well as telling Cactus when to execute thorn functions. \n",
    "\n",
    "Cactus schedules its work in several bins. While there are many of them, here we only consider a handful:\n",
    "\n",
    "1) CCTK_INITIAL - This runs once at the beginning. Initialize your grid functions here.\n",
    "\n",
    "2) CCTK_EVOL - Evolve a single timestep forward. This step will run repeatedly until the simulation finishes. Actually, this is rarely used directly. Usually, MoL_CalcRHS is used instead. More on that in a little bit.\n",
    "\n",
    "3) CCTK_ANALYSIS - Called right before output routines. This ensures values are filled in sensibly. Call it to calculate values that you want to look at, but which aren't necessarily needed for the evolution.\n",
    "\n",
    "4) MOL_REGISTER - This registers variables with the MoL thorn. (Note: events are scheduled IN, not AT, MoL_Register)  \n",
    "\n",
    "5) MOL_CALCRHS - This calculates the RHS's for variables registered with the MoL thorn. (Note: events are scheduled IN, not AT, MoL_CalcRHS)\n",
    "\n",
    "Recall that we have four functions that we need to schedule:\n",
    "\n",
    "- WaveEqn_MoLRegister\n",
    "- WaveEqn_Initialize\n",
    "- WaveEqn_RHS\n",
    "- WaveEqn_Boundary\n",
    "\n",
    "and two groups to allocate storage for:\n",
    "\n",
    "- evol_group\n",
    "- rhs_group\n",
    "\n",
    "The next cell generates the `schedule.ccl` file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile {thorn_dir}/schedule.ccl\n",
    "## Schedule definitions for thorn WaveEqn\n",
    "\n",
    "## There won't be any storage allocated for a group\n",
    "## unless a corresponding storage declaration exists\n",
    "## for it in the schedule file. In square brackets,\n",
    "## we specify the number of storage levels to allocate.\n",
    "STORAGE: evol_group[2]\n",
    "STORAGE: rhs_group\n",
    "\n",
    "## Schedule a function defined in this thorn to run at various stages\n",
    "## of the simulation. The minimum you need to specify for a schedule\n",
    "## item is what language it's written in. Choices are: C (which includes\n",
    "## C++) and Fortran (which means Fortran90).\n",
    "\n",
    "SCHEDULE WaveEqn_MoLRegister IN MoL_Register\n",
    "{\n",
    "  LANG: C\n",
    "  OPTIONS: GLOBAL\n",
    "} \"Register WaveEqn variables\"\n",
    "\n",
    "SCHEDULE WaveEqn_Initialize AT CCTK_INITIAL\n",
    "{\n",
    "   LANG: C\n",
    "} \"Initialize evolved variables\"\n",
    "\n",
    "SCHEDULE WaveEqn_RHS IN MoL_CalcRHS\n",
    "{\n",
    "   LANG: C\n",
    "} \"Evolve the scalar wave equation\"\n",
    "\n",
    "SCHEDULE WaveEqn_Boundary IN MoL_CalcRHS AFTER WaveEqn_RHS\n",
    "{\n",
    "   LANG: C\n",
    "} \"Scalar wave equation BC\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `param.ccl` file is where you declare parameters. In doing so, you specify the allowed values, types, and defaults of those parameters. The parameters declared here are assignable at runtime from a parfile. In our example, we have no parameters, so we leave this file empty (except for comments that you may find useful for future reference)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile {thorn_dir}/param.ccl\n",
    "## Parameter definitions for thorn WaveEqn\n",
    "\n",
    "## There are five types of parameters: INT, REAL, KEYWORD, STRING, \n",
    "## and BOOLEAN. The following comments provide prototypes of each.\n",
    "#\n",
    "#CCTK_INT one_to_five \"This integer parameter goes from 1 to 5\"\n",
    "#{\n",
    "#  1:5 :: \"Another comment\"\n",
    "#} 3 # This is the default value\n",
    "#\n",
    "#CCTK_REAL from_2p5_to_3p8e4 \"This integer parameter goes from 2.5 to 3.8e4\"\n",
    "#{\n",
    "#  2.5:3.8e4 :: \"Another comment\"\n",
    "#} 4.4e3 # This is the default value\n",
    "#\n",
    "## This keyword example defines the parameter wavemaker_type and 8 possible values.\n",
    "#CCTK_KEYWORD wavemaker_type \"types of wave makers\"\n",
    "#{\n",
    "#  \"ini_rec\" :: \"initial rectangular hump, need xc,yc and wid\"\n",
    "#  \"lef_sol\" :: \"initial solitary wave, WKN B solution, need amp, dep\"\n",
    "#  \"ini_oth\" :: \"other initial distribution specified by users\"\n",
    "#  \"ini_gau\" :: \"initial Gaussian hump, need amp, xc, yc, and wid\"\n",
    "#  \"ini_sol\" :: \"initial solitary wave, xwavemaker\"\n",
    "#} \"ini_gau\" # This is the default value\n",
    "#\n",
    "#CCTK_STRING a_string_par \"a comment\"\n",
    "#{\n",
    "#  .* :: \"This is a perl 5 regular expression defining what the string may contain\"\n",
    "#} \"blah blah blah\" # This is the default value\n",
    "#\n",
    "#BOOLEAN a_boolean_par \"a comment\"\n",
    "#{\n",
    "#} true"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `configure.ccl` file... [add description here]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile {thorn_dir}/configure.ccl\n",
    "## Configuration definitions for thorn WaveEqn\n",
    "\n",
    "## You should not need to include \"mpi.h\", but if you\n",
    "## do, you will need this next line.\n",
    "# REQUIRES MPI\n",
    "# REQUIRES HDF5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Update thornlist\n",
    "\n",
    "Before we can use the thorn, we need to add it to the thornlist and then rebuild Cactus."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%cd {cactus_dir}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thorn_list=os.path.join(os.getcwd(), 'configs', 'sim', 'ThornList')\n",
    "thorn_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_thorn_path=os.path.join(os.environ['ARR'], os.environ['THORN'])\n",
    "new_thorn_path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make sure the new thorn is in our ThornList\n",
    "contents = None\n",
    "with open(thorn_list, 'r') as f:\n",
    "    contents = f.read()\n",
    "if new_thorn_path not in contents:\n",
    "    with open(thorn_list, 'a') as f:\n",
    "        f.write(new_thorn_path + '\\n')\n",
    "!tail {thorn_list}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Rebuild Cactus\n",
    "\n",
    "Now rebuild cactus. If all goes well, you will notice Cactus compiling your new thorn."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%cd {cactus_dir}\n",
    "!time ./simfactory/bin/sim build --mdbkey make 'make -j2'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Create parfile\n",
    "\n",
    "We're almost done! The last thing to do before running a simulation is to create a parfile. In the next cell, indicate the path to the directory that will contain the parfile for this simulation. This can be stored anywhere, even outside of your Cactus directory. You may want to store your customized parfiles in a directory where they will not be overwritten. It is up to you."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the next cell we create the parfile, `wave_eqn.par`, in the parfile directory that you just specified."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile {par_dir}/wave_eqn.par\n",
    "Cactus::cctk_run_title = \"scalar\"\n",
    "Cactus::cctk_full_warnings         = yes\n",
    "Cactus::highlight_warning_messages = no\n",
    "Cactus::terminate       = \"time\"\n",
    "Cactus::cctk_final_time = 50\n",
    "\n",
    "ActiveThorns = \"Carpet CarpetLib CarpetInterp CarpetReduce CarpetSlab\"\n",
    "Carpet::verbose           = no\n",
    "Carpet::veryverbose       = no\n",
    "Carpet::schedule_barriers = no\n",
    "Carpet::storage_verbose   = no\n",
    "#Carpet::timers_verbose    = no\n",
    "CarpetLib::output_bboxes  = no\n",
    "\n",
    "Carpet::domain_from_coordbase = yes\n",
    "Carpet::max_refinement_levels = 1\n",
    "driver::ghost_size            = 1\n",
    "Carpet::init_fill_timelevels  = yes\n",
    "\n",
    "ActiveThorns = \"NaNChecker\"\n",
    "NaNChecker::check_every     = 1 # 512\n",
    "#NaNChecker::verbose         = \"all\"\n",
    "#NaNChecker::action_if_found = \"just warn\"\n",
    "NaNChecker::action_if_found = \"terminate\"\n",
    "NaNChecker::check_vars      = \"\n",
    "        WaveEqn::U\n",
    "\"\n",
    "\n",
    "ActiveThorns = \"MoL\"\n",
    "MoL::ode_method                        = \"RK4\"\n",
    "MoL::mol_intermediate_steps            = 4\n",
    "MoL::mol_num_scratch_levels            = 1\n",
    "\n",
    "ActiveThorns = \"Boundary CartGrid3D CoordBase SymBase\"\n",
    "\n",
    "CartGrid3D::type = \"coordbase\"\n",
    "\n",
    "CoordBase::domainsize = \"minmax\"\n",
    "CoordBase::spacing    = \"numcells\"\n",
    "\n",
    "CoordBase::xmin = -10.0\n",
    "CoordBase::ymin = -10.0\n",
    "CoordBase::zmin = -10.0\n",
    "CoordBase::xmax = +10.0\n",
    "CoordBase::ymax = +10.0\n",
    "CoordBase::zmax = +10.0\n",
    "CoordBase::ncells_x  = 50\n",
    "CoordBase::ncells_y  = 50\n",
    "CoordBase::ncells_z  = 50\n",
    "\n",
    "CoordBase::boundary_size_x_lower     = 1\n",
    "CoordBase::boundary_size_y_lower     = 1\n",
    "CoordBase::boundary_size_z_lower     = 1\n",
    "CoordBase::boundary_size_x_upper     = 1\n",
    "CoordBase::boundary_size_y_upper     = 1\n",
    "CoordBase::boundary_size_z_upper     = 1\n",
    "\n",
    "ActiveThorns = \"Time\"\n",
    "Time::dtfac = 0.25\n",
    "\n",
    "ActiveThorns = \"InitBase\"\n",
    "\n",
    "ActiveThorns = \"WaveEqn\"\n",
    "\n",
    "ActiveThorns = \"IOUtil\"\n",
    "IO::out_dir = $parfile\n",
    "IO::checkpoint_dir                  = $parfile\n",
    "IO::checkpoint_ID                   = no\n",
    "IO::checkpoint_every_walltime_hours = 6.0\n",
    "IO::checkpoint_on_terminate         = yes\n",
    "IO::recover     = \"autoprobe\"\n",
    "IO::recover_dir = $parfile\n",
    "\n",
    "ActiveThorns=\"CarpetIOHDF5\"\n",
    "IOHDF5::out2D_every            = 16\n",
    "IOHDF5::out2D_xz               = no\n",
    "IOHDF5::out2D_yz               = no\n",
    "IOHDF5::output_buffer_points   = yes\n",
    "#IOHDF5::one_file_per_group     = yes\n",
    "IOHDF5::output_symmetry_points = no\n",
    "IOHDF5::compression_level      = 1\n",
    "IOHDF5::use_checksums          = yes\n",
    "IOHDF5::out2D_vars             = \"\n",
    "        WaveEqn::U \n",
    "        WaveEqn::K\n",
    "        Grid::Coordinates{out_every=1000000000} #we only want this to print once\n",
    "\"\n",
    "IOHDF5::checkpoint             = no"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Run simulation\n",
    "\n",
    "Finally, we are ready to run the simulation using our new thorn. In the next cell, indicate the path to the directory where simulation output is written on your machine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Linux example:\n",
    "#sims_dir = '/home/ejwest/ETK/simulations'\n",
    "\n",
    "## Windows example:\n",
    "#sims_dir='C:\\\\Users\\\\EJWest\\\\ETK\\\\Simulations'\n",
    "\n",
    "## fill in the path to your simulations directory\n",
    "sims_dir=os.path.join(home_dir, \"simulations\") #<-- COMPLETE THIS LINE!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remove any old versions of the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%rm -rf {sims_dir}/wave_eqn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now run the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%cd {cactus_dir}\n",
    "!time ./simfactory/bin/sim create-run wave_eqn --parfile={par_dir}/wave_eqn.par --procs=2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. Visualize output\n",
    "\n",
    "The output of this simulation are `.h5` files, which are HDF5 files. 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 normally use for visualizing large simulations. At this point you could switch over to VisIt in order to visualize the output. However it is also possible to use a Jupyter notebook using the H5Py library, which is a Python library that provides methods for interfacing with HDF5 files. Once the data are extracted from the `.h5` files, they can be visualized as normal using matplotlib."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extract data from hdf5 files\n",
    "\n",
    "Let's first move to the output directory. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%cd {sims_dir}/wave_eqn/output-0000/wave_eqn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And then list the available `.h5` files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%ls *.h5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will extract data from `U.xy.h5` into the notebook, specifically into NumPy arrays. To do this, we need to import NumPy and H5Py."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## import libraries\n",
    "import numpy as np\n",
    "import h5py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## create file handle\n",
    "u_h5 = h5py.File(\"U.xy.h5\")"
   ]
  },
  {
   "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 nm in u_h5:\n",
    "    if hasattr(u_h5[nm], 'shape'):\n",
    "        ts = u_h5[nm].attrs['timestep']\n",
    "        time_steps = np.append(time_steps, ts)\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. This information is not stored in the `U.xy.h5` file, but instead is found in `x.xy.h5` and `y.xy.h5`. The following cell retrieves this information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## get xdata, ydata\n",
    "x_h5 = h5py.File(\"x.xy.h5\")\n",
    "y_h5 = h5py.File(\"y.xy.h5\")\n",
    "\n",
    "for nm in x_h5:\n",
    "    if hasattr(x_h5[nm], 'shape'):\n",
    "        xdata = x_h5[nm][...]\n",
    "for nm in y_h5:\n",
    "    if hasattr(y_h5[nm], 'shape'):\n",
    "        ydata = y_h5[nm][...]\n",
    "\n",
    "## print to screen (optional)        \n",
    "print(\"xdata = \", xdata)\n",
    "print(\"ydata = \", ydata)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell extracts the values of U over the xy-grid at each time step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## get zdata\n",
    "zdata = np.zeros((len(time_steps), xdata.shape[0], xdata.shape[1]))\n",
    "for i in range(len(time_steps)):\n",
    "    tstep = time_steps[i]\n",
    "    ##find dataset with matching time step\n",
    "    for nm in u_h5:\n",
    "        if hasattr(u_h5[nm], 'shape'):\n",
    "            if u_h5[nm].attrs['timestep'] == tstep: \n",
    "                ## save dataset as element of numpy array\n",
    "                zdata[i,:,:] = u_h5[nm][...]\n",
    "\n",
    "## print to screen (optional)                \n",
    "#print('zdata =',zdata)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot using matplotlib's plot_surface( )\n",
    "\n",
    "Now that the data we need are extracted from the HDF5 files, we can plot them as we normally would. Here we use matplotlib's `plot_surface()` command."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot a snapshot of the solution using `plot_surface()`."
   ]
  },
  {
   "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",
    "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(zdata)\n",
    "maxval = np.max(zdata)\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 = xdata\n",
    "yi = ydata\n",
    "zi = zdata[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": [
    "Put all time steps together into an animation using matplotlib's `animation.FuncAnimation()` command."
   ]
  },
  {
   "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",
    "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(zdata)\n",
    "maxval = np.max(zdata)\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 = xdata\n",
    "yi = ydata\n",
    "zi = zdata[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, xdata, ydata, zdata, scene):\n",
    "    time = time_steps[frame_number]\n",
    "    fig.suptitle('time = %.3f' % time)\n",
    "    scene[0].remove()\n",
    "    xi = xdata\n",
    "    yi = ydata\n",
    "    zi = zdata[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, xdata, ydata, zdata, 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": [
    "Save the animation to a file, if you like. In order to save as an `.mp4` file, you need to have the ffmpeg library installed on your machine. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## save animation to file\n",
    "file = 'wave_eqn'\n",
    "#anim.save(file + '.mp4', writer='ffmpeg', fps=fps)      # requires ffmpeg\n",
    "anim.save(file + '.gif', writer='imagemagick', fps=fps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8. Exercises\n",
    "\n",
    "1. Change the time stepsize and termination time of the simulation.\n",
    "2. Change the grid size and grid spacing of the simulation.\n",
    "3. Change the height/width/location of the initial pulse.\n",
    "4. Add a wave speed parameter and run for different values.\n",
    "5. Write a new function that calculates the energy of the configuration.\n",
    "6. Change the boundary conditions: Dirichlet, Neumann, mixed, periodic \n",
    "7. Change the underlying equation: inviscid Burger's equation, Navier-Stokes equation \n",
    "8. Add adaptive mesh refinement using Carpet. (advanced!)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Acknowledgments\n",
    "\n",
    "This notebook grew out of a series of notebook-based tutorials delivered by Steve Brandt at the 2017 Einstein Toolkit school and workshop at NCSA. Those tutorials were based on his \"Funwave\" example--a simulation that models coastal waves. The current notebook is an attempt at a condensed version of the \"create a new thorn\" tutorial for new user's. To that end, the Funwave example has been replaced by a simpler case study--the standard wave equation--in order to fit everything into a single notebook. The aim has been to provide a more transparent example of how to create a new thorn from scratch. For new users, this is a daunting task, which the author of this notebook fully recognizes.\n",
    "\n",
    "A big thanks goes to Steve Brandt for his outstanding series of Jupyter notebooks on the in's and out's of the Einstein Toolkit. His tutorials go far beyond this one, both in depth and breadth. They cover many common gotcha issues that this notebook glosses over. I would encourage any new user to work through them; it is well worth the effort to get them to work on your machine. They can be found in the Google drive [here](https://drive.google.com/drive/folders/0B4gNfWainf-5R3BldGJuRlo4dHc) and video of the 2017 workshop sessions can be found on YouTube [here](https://www.youtube.com/watch?v=jqgEaUjl23I), starting at the 4:08:37 mark.\n",
    "\n",
    "Thanks also goes to Yosef Zlochower, whose wave equation example has been used in this notebook.\n",
    "\n",
    "My role has basically been to merge the pedagogical clarity and thoroughness of Steve's tutorials (including his fully commented ccl file templates) with the simplicity of Yosef's example. I am fully responsible for any bugs or confusions that remain in the final result. --Eric J West"
   ]
  }
 ],
 "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
}