Using-Cactus/WaveEquations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1>A Toy Wave Equation</h1>\n",
    "Here we build a very small Cactus application to evolve the wave equation. The main new thing we are learning\n",
    "about here is the MoL (Method of Lines) thorn. This thorn abstracts the time integration of your simulation\n",
    "from the physics code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%cd ~/CactusFW2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're going through the boilerplate of making a new thorn, but we'll do some special\n",
    "things along the way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define some basic parameters describing a new thorn\n",
    "thorn_pars = {\n",
    "   \"thorn_name\" : \"WaveToyET18\",\n",
    "   \"arrangement_name\" : \"FunwaveUtils\",\n",
    "   \"author\" : \"Steven R. Brandt\",\n",
    "   \"email\" : \"sbrandt@cct.lsu.edu\",\n",
    "   \"license\" : \"BSD\"\n",
    " }\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 delete the thorn you are\n",
    "going to create below and start over."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!rm -fr arrangements/$ARR/$THORN"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we define some functions to help us create files automatically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import re\n",
    "# This function does substitutes all occurances of \"{name}\" in input_str with values[\"name\"]\n",
    "# and returns the new string.\n",
    "def replace_values(input_str,values):\n",
    "    while True:\n",
    "        g = re.search(r'{(\\w+)}',input_str)\n",
    "        if g:\n",
    "            var = g.group(1)\n",
    "            if var in values:\n",
    "                val = values[var]\n",
    "            else:\n",
    "                raise Exception(\"Undefined: <<\"+var+\">>\")\n",
    "            start = g.start(0)\n",
    "            end = start+len(g.group(0))\n",
    "            input_str = input_str[0:start]+val+input_str[end:]\n",
    "            continue\n",
    "        break\n",
    "    return input_str\n",
    "\n",
    "# Createas a file with the given name, uses replace_values()\n",
    "# to update file_contents with file_values, and writes it out.\n",
    "def create_file(file_name,file_contents,file_values):\n",
    "    fd = open(file_name,\"w\")\n",
    "    file_contents = replace_values(file_contents,file_values)\n",
    "    print(\"Over-writing file '\"+file_name+\"'\")\n",
    "    file_contents = re.sub(r'^\\s+','',file_contents)\n",
    "    fd.write(file_contents)\n",
    "    fd.close()\n",
    "    \n",
    "# The equivalent of mkdir -p\n",
    "def create_dir(dir):\n",
    "    print(\"Ensuring directory '\"+dir+\"'\")\n",
    "    if not os.path.exists(dir):\n",
    "        os.makedirs(dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3>Something new!</h3>\n",
    "In the interface.ccl file we not only declare a grid function, but we give it scratch space of three time levels.\n",
    "This gives MoL something to work with. Different time integrators will need different numbers of time levels.\n",
    "\n",
    "Perhaps also unexpected are the rhs variables. The phi_rhs will hold the time derivative of phi, the psi_rhs will\n",
    "hold the time derivative of psi. MoL will call the evolution routine repeatedly, storing the intermediate steps as\n",
    "necessary for its time integration.\n",
    "\n",
    "But wait! There's more! We also declare that we require a method called \"MoLRegisterEvolvedGroup\"\n",
    "be provided by some thorn."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "interface_ccl_contents = \"\"\"\n",
    "## Interface definitions for thorn {thorn_name}\n",
    "inherits: grid\n",
    "## An implementation name is required for all thorns. No\n",
    "## two thorns in a configuration can implement the same\n",
    "## interface.\n",
    "implements: {thorn_name}\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. timelevels, tags such as the\n",
    "## Prolongation=None tag. The type tag can take on the values\n",
    "## GF, Scalar, or Array.\n",
    "\n",
    "## Note that the number of timelevels can be an integer parameter\n",
    "## GF stands for \"Grid Function\" and refers to a distributed array\n",
    "## data structure.\n",
    "#cctk_real force_group type=GF timelevels=3 tags='Prolongation=\"None\"'\n",
    "#{\n",
    "#  force1, force2\n",
    "#}\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",
    "cctk_real evol_group type=GF timelevels=3\n",
    "{\n",
    "  psi, phi\n",
    "}\n",
    "cctk_real rhs_group type=GF \n",
    "{\n",
    "  psi_rhs, phi_rhs\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\n",
    "\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3>Something new!</h3>\n",
    "We schedule three routines. The initial and evolution routines are unsurprising. The register method, however, is\n",
    "new. This is the place where we will tell MoL what variables we are evolving."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "schedule_ccl_contents = \"\"\"\n",
    "## Schedule definitions for thorn {thorn_name}\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",
    "## These commented examples correspond to the commented\n",
    "## examples in the interface file above.\n",
    "# storage: force_group[3]\n",
    "# storage: scalar_group\n",
    "\n",
    "## Schedule a function defined in this thorn to run at the beginning\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",
    "#SCHEDULE init_function at CCTK_INIT\n",
    "#{\n",
    "#  LANG: C\n",
    "#}\"Do some initial stuff\"\n",
    "\n",
    "storage: evol_group[3], rhs_group\n",
    "\n",
    "SCHEDULE init_wave_toy in CCTK_INITIAL\n",
    "{\n",
    "  LANG: C\n",
    "} \"initialize data\"\n",
    "\n",
    "SCHEDULE evolve_wave_toy in MoL_CalcRHS\n",
    "{\n",
    "  LANG: C\n",
    "} \"evolve data\"\n",
    "\n",
    "schedule wave_toy_registervars in MoL_Register\n",
    "{\n",
    "  LANG: C\n",
    "}\"Register funwave variables for MoL\"\n",
    "\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "param_ccl_contents = \"\"\"\n",
    "## Parameter definitions for thorn {thorn_name}\n",
    "## There are five types of parameters: int, real, keyword, string, and boolean.\n",
    "## The 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",
    "#  \"wk_reg\" :: \"Wei and Kirby 1999 internal wave maker, need xc_wk, tperiod, amp_wk, dep_wk, theta_wk, and time_ramp (factor of period)\"\n",
    "#  \"wk_irr\" :: \"Wei and Kirby 1999 TMA spectrum wavemaker, need xc_wk, dep_wk, time_ramp, delta_wk, freqpeak, freqmin, freqmax, hmo, gammatma, theta-peak\"\n",
    "#  \"wk_time_series\" :: \"fft a time series to get each wave component and then use Wei and Kirby's ( 1999) wavemaker. Need input wavecompfile (including 3 columns:   per,amp,pha) and numwavecomp, peakperiod, dep_wk, xc_wk, ywidth_wk\"\n",
    "#  \"ini_gau\" :: \"initial Gaussian hump, need amp, xc, yc, and wid\"\n",
    "#  \"ini_sol\" :: \"initial solitary wave, xwavemaker\"\n",
    "#}\"wk_reg\" # 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\n",
    "\n",
    "\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "configuration_ccl_contents = \"\"\"\n",
    "# Configuration definitions for thorn {thorn_name}\n",
    "## You should not need include \"mpi.h\", but if you\n",
    "## do, you will need this next line.\n",
    "# REQUIRES MPI\n",
    "# REQUIRES HDF5\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "makefile_contents = \"\"\"\n",
    "# Main make.code.defn file for thorn {thorn_name}\n",
    "\n",
    "# Source files in this directory\n",
    "SRCS = init.cc evolve.cc molregister.cc\n",
    "\n",
    "# Subdirectories containing source files\n",
    "SUBDIRS =\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "readme_contents = \"\"\"\n",
    "Author(s)    : {author} <{email}>\n",
    "Maintainer(s): {author} <{email}>\n",
    "Licence      : {license}\n",
    "--------------------------------------------------------------------------\n",
    "\n",
    "1. Purpose\n",
    "\n",
    "not documented\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "# This function will create a complete thorn\n",
    "def create_thorn():\n",
    "  # Create the thorn directory inside the Cactus source tree\n",
    "  arrangement_dir = \"arrangements/\"+thorn_pars[\"arrangement_name\"]\n",
    "  thorn_dir = arrangement_dir + \"/\" + thorn_pars[\"thorn_name\"]\n",
    "  create_dir(thorn_dir)\n",
    "\n",
    "  # Create basic ccl files needed by Cactus\n",
    "  schedule_ccl = thorn_dir+\"/schedule.ccl\"\n",
    "  interface_ccl = thorn_dir+\"/interface.ccl\"\n",
    "  param_ccl = thorn_dir+\"/param.ccl\"\n",
    "  configuration_ccl = thorn_dir+\"/configuration.ccl\"\n",
    "  create_file(schedule_ccl,schedule_ccl_contents,thorn_pars)\n",
    "  create_file(interface_ccl,interface_ccl_contents,thorn_pars)\n",
    "  create_file(param_ccl,param_ccl_contents,thorn_pars)\n",
    "  create_file(configuration_ccl,configuration_ccl_contents,thorn_pars)\n",
    "\n",
    "  # Create the source directory and first makefile\n",
    "  src_dir = thorn_dir+\"/src\"\n",
    "  create_dir(src_dir)\n",
    "  create_file(src_dir+\"/make.code.defn\",makefile_contents,thorn_pars)\n",
    "\n",
    "  # Other dirs and files not strictly needed\n",
    "  # for compiling and running Cactus.\n",
    "  create_file(thorn_dir+\"/README\",readme_contents,thorn_pars)\n",
    "  \n",
    "  test_dir = thorn_dir+\"/test\"\n",
    "  create_dir(test_dir)\n",
    "  par_dir = thorn_dir+\"/par\"\n",
    "  create_dir(par_dir)\n",
    "  doc_dir = thorn_dir+\"/doc\"\n",
    "  create_dir(doc_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "create_thorn()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_thorns_contents=\"\"\"\n",
    "# ./configs/sim/ThornList\n",
    "# This file was automatically generated using the GetComponents script.\n",
    "\n",
    "!CRL_VERSION = 2.0\n",
    "\n",
    "\n",
    "# Component list: funwave.th\n",
    "\n",
    "!DEFINE ROOT = CactusFW2\n",
    "!DEFINE ARR = $ROOT/arrangements\n",
    "!DEFINE ET_RELEASE = trunk\n",
    "!DEFINE FW_RELEASE = FW_2014_05\n",
    "\n",
    "#Cactus Flesh\n",
    "!TARGET   = $ROOT\n",
    "!TYPE     = git\n",
    "!URL      = https://bitbucket.org/cactuscode/cactus.git\n",
    "!NAME     = flesh\n",
    "!CHECKOUT = CONTRIBUTORS COPYRIGHT doc lib Makefile src\n",
    "\n",
    "!TARGET   = $ARR\n",
    "!TYPE     = git\n",
    "!URL      = https://bitbucket.org/stevenrbrandt/cajunwave.git\n",
    "!REPO_PATH= $2\n",
    "# Old version\n",
    "#!AUTH_URL = https://svn.cct.lsu.edu/repos/projects/ngchc/code/branches/$FW_RELEASE/$1/$2\n",
    "#!URL = https://svn.cct.lsu.edu/repos/projects/ngchc/code/branches/$FW_RELEASE/$2\n",
    "#!URL = https://svn.cct.lsu.edu/repos/projects/ngchc/code/CactusCoastal/$2\n",
    "!CHECKOUT =\n",
    "CactusCoastal/Funwave\n",
    "CactusCoastal/FunwaveMesh\n",
    "CactusCoastal/FunwaveCoord\n",
    "CactusCoastal/Tridiagonal\n",
    "CactusCoastal/Tridiagonal2\n",
    "\n",
    "# CactusBase thorns\n",
    "!TARGET   = $ARR\n",
    "!TYPE     = git\n",
    "!URL      = https://bitbucket.org/cactuscode/cactusbase.git\n",
    "!REPO_PATH= $2\n",
    "!CHECKOUT =\n",
    "CactusBase/Boundary\n",
    "CactusBase/CartGrid3D\n",
    "CactusBase/CoordBase\n",
    "CactusBase/Fortran\n",
    "CactusBase/InitBase\n",
    "CactusBase/IOASCII\n",
    "CactusBase/IOBasic\n",
    "CactusBase/IOUtil\n",
    "CactusBase/SymBase\n",
    "CactusBase/Time\n",
    "#\n",
    "# CactusNumerical thorns\n",
    "!TARGET   = $ARR\n",
    "!TYPE     = git\n",
    "!URL      = https://bitbucket.org/cactuscode/cactusnumerical.git\n",
    "!REPO_PATH= $2\n",
    "!CHECKOUT =\n",
    "!CHECKOUT =\n",
    "CactusNumerical/MoL\n",
    "CactusNumerical/LocalInterp\n",
    "\n",
    "CactusNumerical/Dissipation\n",
    "CactusNumerical/SpaceMask\n",
    "CactusNumerical/SphericalSurface\n",
    "CactusNumerical/LocalReduce\n",
    "CactusNumerical/InterpToArray\n",
    "\n",
    "!TARGET   = $ARR\n",
    "!TYPE     = git\n",
    "!URL      = https://bitbucket.org/cactuscode/cactusutils.git\n",
    "!REPO_PATH= $2\n",
    "!CHECKOUT = CactusUtils/Accelerator CactusUtils/OpenCLRunTime\n",
    "CactusUtils/NaNChecker\n",
    "CactusUtils/Vectors\n",
    "CactusUtils/SystemTopology\n",
    "\n",
    "# Carpet, the AMR driver\n",
    "!TARGET   = $ARR\n",
    "!TYPE     = git\n",
    "!URL      = https://bitbucket.org/eschnett/carpet.git\n",
    "!REPO_PATH= $2\n",
    "!CHECKOUT = Carpet/doc\n",
    "Carpet/Carpet\n",
    "Carpet/CarpetEvolutionMask\n",
    "Carpet/CarpetIOASCII\n",
    "Carpet/CarpetIOBasic\n",
    "Carpet/CarpetIOHDF5\n",
    "Carpet/CarpetIOScalar\n",
    "#Carpet/CarpetIntegrateTest\n",
    "Carpet/CarpetInterp\n",
    "Carpet/CarpetInterp2\n",
    "Carpet/CarpetLib\n",
    "Carpet/CarpetMask\n",
    "#Carpet/CarpetProlongateTest\n",
    "Carpet/CarpetReduce\n",
    "Carpet/CarpetRegrid\n",
    "Carpet/CarpetRegrid2\n",
    "#Carpet/CarpetRegridTest\n",
    "Carpet/CarpetSlab\n",
    "Carpet/CarpetTracker\n",
    "Carpet/CycleClock\n",
    "#Carpet/HighOrderWaveTest\n",
    "Carpet/LoopControl\n",
    "#Carpet/ReductionTest\n",
    "#Carpet/ReductionTest2\n",
    "#Carpet/ReductionTest3\n",
    "#Carpet/RegridSyncTest\n",
    "Carpet/TestCarpetGridInfo\n",
    "Carpet/TestLoopControl\n",
    "Carpet/Timers\n",
    "\n",
    "# Additional Cactus thorns\n",
    "!TARGET   = $ARR\n",
    "!TYPE     = svn\n",
    "!URL      = https://svn.cactuscode.org/projects/$1/$2/trunk\n",
    "!CHECKOUT = ExternalLibraries/OpenBLAS ExternalLibraries/OpenCL ExternalLibraries/pciutils ExternalLibraries/PETSc\n",
    "ExternalLibraries/MPI\n",
    "ExternalLibraries/HDF5\n",
    "ExternalLibraries/zlib\n",
    "ExternalLibraries/hwloc\n",
    "\n",
    "# Simulation Factory\n",
    "!TARGET   = $ROOT/simfactory\n",
    "!TYPE     = git\n",
    "!URL      = https://bitbucket.org/simfactory/simfactory2.git\n",
    "!NAME     = simfactory2\n",
    "!CHECKOUT = README.md README_FIRST.txt bin doc etc lib mdb\n",
    "\n",
    "# Various thorns from LSU\n",
    "#!TARGET   = $ARR\n",
    "#!TYPE     = git\n",
    "#!URL      = https://bitbucket.org/einsteintoolkit/archivedthorns-vectors.git\n",
    "#!REPO_PATH= $2\n",
    "#!CHECKOUT =\n",
    "#LSUThorns/Vectors\n",
    "#LSUThorns/QuasiLocalMeasures\n",
    "#LSUThorns/SummationByParts\n",
    "#LSUThorns/Prolong\n",
    "\n",
    "#Roland/MapPoints\n",
    "#Tutorial/BadWaveMoL\n",
    "#Tutorial/BasicWave\n",
    "#Tutorial/BasicWave2\n",
    "#Tutorial/BasicWave3\n",
    "\n",
    "# Various thorns from the AEI\n",
    "# Numerical\n",
    "!TARGET   = $ARR\n",
    "!TYPE     = git\n",
    "!URL      = https://bitbucket.org/cactuscode/numerical.git\n",
    "!REPO_PATH= $2\n",
    "!CHECKOUT =\n",
    "#AEIThorns/ADMMass\n",
    "AEIThorns/AEILocalInterp\n",
    "#AEIThorns/PunctureTracker\n",
    "#AEIThorns/SystemStatistics\n",
    "#AEIThorns/Trigger\n",
    "{arrangement_name}/{thorn_name}\"\"\"\n",
    "create_file(\"my_thorns.th\",my_thorns_contents,thorn_pars)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3>Something new!</h3>\n",
    "Here we create the source code. The initialization, evolution, and molregister functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "src_dir = \"arrangements/\"+thorn_pars[\"arrangement_name\"]+\"/\"+thorn_pars[\"thorn_name\"]+\"/src\"\n",
    "\n",
    "init_cc_src = \"\"\"\n",
    "#include <cctk.h>\n",
    "#include <cctk_Arguments.h>\n",
    "#include <cctk_Parameters.h>\n",
    "#include <math.h>\n",
    "\n",
    "void init_wave_toy(CCTK_ARGUMENTS) {\n",
    "  DECLARE_CCTK_ARGUMENTS;\n",
    "  DECLARE_CCTK_PARAMETERS;\n",
    "  for(int k=0;k<cctk_lsh[2];k++) {\n",
    "      for(int j=0;j<cctk_lsh[1];j++) {\n",
    "          for(int i=0;i<cctk_lsh[0];i++) {\n",
    "              int cc = CCTK_GFINDEX3D(cctkGH,i,j,k);\n",
    "              CCTK_REAL xx = x[cc] - 5.0;\n",
    "              CCTK_REAL yy = y[cc] - 5.0;\n",
    "              psi[cc] = exp(-xx*xx-yy*yy);\n",
    "              phi[cc] = 0;\n",
    "          }\n",
    "      }\n",
    "  }\n",
    "}\n",
    "\"\"\"\n",
    "create_file(src_dir+\"/init.cc\",init_cc_src,{})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "evolve_cc_src = \"\"\"\n",
    "#include <cctk.h>\n",
    "#include <cctk_Arguments.h>\n",
    "#include <cctk_Parameters.h>\n",
    "#include <math.h>\n",
    "#include <assert.h>\n",
    "\n",
    "void evolve_wave_toy(CCTK_ARGUMENTS) {\n",
    "  DECLARE_CCTK_ARGUMENTS;\n",
    "  DECLARE_CCTK_PARAMETERS;\n",
    "  for(int k=0;k<cctk_lsh[2];k++) {\n",
    "      for(int j=1;j<cctk_lsh[1]-1;j++) {\n",
    "          for(int i=1;i<cctk_lsh[0]-1;i++) {\n",
    "              int cc = CCTK_GFINDEX3D(cctkGH,i,j,k);\n",
    "              int im1 = CCTK_GFINDEX3D(cctkGH,i-1,j,k);\n",
    "              int ip1 = CCTK_GFINDEX3D(cctkGH,i+1,j,k);\n",
    "              int jm1 = CCTK_GFINDEX3D(cctkGH,i,j-1,k);\n",
    "              int jp1 = CCTK_GFINDEX3D(cctkGH,i,j+1,k);\n",
    "              CCTK_REAL dx = CCTK_DELTA_SPACE(0);\n",
    "              CCTK_REAL dy = CCTK_DELTA_SPACE(1);\n",
    "              psi_rhs[cc] = phi[cc];\n",
    "              CCTK_REAL divx = (psi[ip1] + psi[im1] - 2*psi[cc])/(dx*dx);\n",
    "              CCTK_REAL divy = (psi[jp1] + psi[jm1] - 2*psi[cc])/(dy*dy);\n",
    "              phi_rhs[cc] = divx + divy;\n",
    "          }\n",
    "      }\n",
    "  }\n",
    "}\n",
    "\"\"\"\n",
    "create_file(src_dir+\"/evolve.cc\",evolve_cc_src,{})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "molregister_cc_src = \"\"\"\n",
    "#include <cctk.h>\n",
    "#include <cctk_Arguments.h>\n",
    "#include <cctk_Parameters.h>\n",
    "\n",
    "extern \"C\" void wave_toy_registervars (CCTK_ARGUMENTS)\n",
    "{\n",
    "  DECLARE_CCTK_ARGUMENTS;\n",
    "  DECLARE_CCTK_PARAMETERS;\n",
    "  CCTK_INT ierr = 0, group, rhs;\n",
    "\n",
    "  // Look up the group index of the evol_group\n",
    "  group = CCTK_GroupIndex (\"{thorn_name}::evol_group\");\n",
    "  \n",
    "  // Look up the group index of the rhs_group\n",
    "  rhs = CCTK_GroupIndex (\"{thorn_name}::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",
    "    CCTK_WARN (0, \"Problems registering with MoL !\");\n",
    "}\n",
    "\n",
    "\"\"\"\n",
    "create_file(src_dir+\"/molregister.cc\",molregister_cc_src,thorn_pars)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!time ./simfactory/bin/sim build --mdbkey make 'make -j2' --thornlist=./my_thorns.th | cat -"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A minimal par file....\n",
    "\n",
    "wave_toy_par_src = \"\"\"\n",
    "ActiveThorns = \"{thorn_name} coordbase carpet cartgrid3d MoL\"\n",
    "\n",
    "Cactus::cctk_itlast = 250\n",
    "\n",
    "# the output dir will be named after the parameter file name\n",
    "IO::out_dir = $parfile\n",
    "\n",
    "# The set of parameters required to specify the grid...\n",
    "Carpet::domain_from_coordbase = \"yes\"\n",
    "CartGrid3D::type = \"coordbase\"\n",
    "CartGrid3D::avoid_origin = \"no\"\n",
    "CoordBase::domainsize = \"minmax\"\n",
    "CoordBase::spacing    = \"gridspacing\"\n",
    "CoordBase::xmin =  0\n",
    "CoordBase::xmax =  30\n",
    "CoordBase::ymin =  0\n",
    "CoordBase::ymax =  30\n",
    "CoordBase::zmin =  0.0\n",
    "CoordBase::zmax =  0.0\n",
    "CoordBase::dx   =  0.25\n",
    "CoordBase::dy   =  0.25\n",
    "\n",
    "# Because we're doing a 2D code...\n",
    "CoordBase::boundary_size_z_lower     = 0\n",
    "CoordBase::boundary_size_z_upper     = 0\n",
    "CoordBase::boundary_shiftout_z_lower = 1\n",
    "CoordBase::boundary_shiftout_z_upper = 1\n",
    "\n",
    "# The time integration scheme...\n",
    "MoL::ODE_Method = \"RK3\"\n",
    "\n",
    "# The next section is required if you want some kind of confirmation\n",
    "# on the console that your code is actually doing anything...\n",
    "ActiveThorns = \"CarpetIOBasic IOUtil CarpetReduce SymBase\"\n",
    "IO::out_fileinfo=\"none\"\n",
    "IOBasic::outInfo_every = 1\n",
    "IOBasic::outInfo_vars = \"{thorn_name}::phi {thorn_name}::psi\"\n",
    "\n",
    "# Get some data!\n",
    "ActiveThorns = \"CarpetIOASCII\"\n",
    "CarpetIOASCII::compact_format = false\n",
    "IOASCII::out2D_every = 25\n",
    "IOASCII::out2D_xyplane_z = 0\n",
    "IOASCII::out2D_vars = \"{thorn_name}::psi {thorn_name}::phi\"\n",
    "IOASCII::out2D_xz = \"no\"\n",
    "IOASCII::out2D_yz = \"no\"\n",
    "IOASCII::output_ghost_points = \"no\"\n",
    "\n",
    "# Taking too big a time step can lead to unstable evolution.\n",
    "# You can control that with dtfac and thorn Time.\n",
    "ActiveThorns = \"Time\"\n",
    "Time::dtfac = .2\n",
    "\"\"\"\n",
    "create_file(\"wave_toy.par\",wave_toy_par_src,thorn_pars)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!rm -fr ../simulations/wave_toy\n",
    "!./simfactory/bin/sim create-run wave_toy.par --procs=1 --num-threads=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This cell enables inline plotting in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "import matplotlib\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "# https://matplotlib.org/examples/color/colormaps_reference.html\n",
    "cmap = cm.prism"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file_data = np.genfromtxt(\"../simulations/wave_toy/output-0000/wave_toy/psi.xy.asc\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sets = np.unique(file_data[:,0])\n",
    "width = 8\n",
    "height = 4\n",
    "print(\"sets=\",sets)\n",
    "mn, mx = np.min(file_data[:,12]),np.max(file_data[:,12])\n",
    "#mn, mx = 0,1\n",
    "for which in sets:\n",
    "    g = file_data[file_data[:,0]==which,:]\n",
    "    print(\"time=\",which)\n",
    "    x = g[:,5]\n",
    "    y = g[:,6]\n",
    "    z = g[:,12]\n",
    "    zi = z.reshape(len(np.unique(y)),len(np.unique(x)))\n",
    "    print('min/max=',np.min(zi),np.max(zi))\n",
    "    plt.figure(figsize=(width, height))\n",
    "    plt.imshow(zi[::-1,:],cmap,clim=(mn,mx))\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3>Question</h3>: Does the above simulation work properly if we use 2 procs? Why or why not? Try it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "schedule_ccl_contents = \"\"\"\n",
    "## Schedule definitions for thorn {thorn_name}\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",
    "## These commented examples correspond to the commented\n",
    "## examples in the interface file above.\n",
    "# storage: force_group[3]\n",
    "# storage: scalar_group\n",
    "\n",
    "## Schedule a function defined in this thorn to run at the beginning\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",
    "#SCHEDULE init_function at CCTK_INIT\n",
    "#{\n",
    "#  LANG: C\n",
    "#}\"Do some initial stuff\"\n",
    "\n",
    "storage: evol_group[3], rhs_group\n",
    "\n",
    "SCHEDULE init_wave_toy in CCTK_INITIAL\n",
    "{\n",
    "  LANG: C\n",
    "} \"initialize data\"\n",
    "\n",
    "SCHEDULE evolve_wave_toy in MoL_CalcRHS\n",
    "{\n",
    "  LANG: C\n",
    "  SYNC: {thorn_name}::evol_group\n",
    "  SYNC: {thorn_name}::rhs_group\n",
    "} \"evolve data\"\n",
    "\n",
    "schedule wave_toy_registervars in MoL_Register\n",
    "{\n",
    "  LANG: C\n",
    "}\"Register funwave variables for MoL\"\n",
    "\n",
    "\"\"\"\n",
    "thorn_dir = \"arrangements/\"+thorn_pars[\"arrangement_name\"]+\"/\"+thorn_pars[\"thorn_name\"]\n",
    "create_file(thorn_dir+\"/schedule.ccl\",schedule_ccl_contents,thorn_pars)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!time ./simfactory/bin/sim build --mdbkey make 'make -j2' --thornlist=./my_thorns.th | cat -"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!rm -fr ../simulations/wave_toy\n",
    "!./simfactory/bin/sim create-run wave_toy.par --procs=2 --num-threads=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file_data = np.genfromtxt(\"../simulations/wave_toy/output-0000/wave_toy/psi.xy.asc\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sets = np.unique(file_data[:,0])\n",
    "width = 8\n",
    "height = 4\n",
    "print(\"sets=\",sets)\n",
    "mn, mx = np.min(file_data[:,12]),np.max(file_data[:,12])\n",
    "#mn, mx = 0,1\n",
    "for which in sets:\n",
    "    g = file_data[file_data[:,0]==which,:]\n",
    "    print(\"time=\",which)\n",
    "    x = g[:,5]\n",
    "    y = g[:,6]\n",
    "    z = g[:,12]\n",
    "    zi = z.reshape(len(np.unique(y)),len(np.unique(x)))\n",
    "    print('min/max=',np.min(zi),np.max(zi))\n",
    "    plt.figure(figsize=(width, height))\n",
    "    plt.imshow(zi[::-1,:],cmap,clim=(mn,mx))\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<h3>Exercises:</h3>\n",
    "<ol>\n",
    "<li>Add a parameter to param.ccl and use it to change the position of the Gaussian wave hump.</li>\n",
    "<li>Advanced: Write a Conway's Game of Life thorn. For this thorn it does not make sense to use MoL. Schedule your routine in CCTK_EVOL.</li>\n",
    "</ol>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>Adaptive Mesh Refinement</h2>\n",
    "Sometimes you care more about what is happening in one part of a grid than another. Wouldn't it be nice if you could\n",
    "spend less time computing in regions that are less important? That way you could explore larger problems with less compute cycles.\n",
    "\n",
    "Well, dream no more. Adaptive Mesh Refinement is supported in Cactus/Carpet. And, at least for the toy example we are\n",
    "considering in this file, all you have to change to make use of it is your parameter file. (Actually, what we're going\n",
    "to consider here is technically Fixed Mesh Refinement (FMR))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3>How AMR works</h3>\n",
    "For reasons of simplicity, we'll focus on a 2 refinement level system. Thus, we have a coarse (refinement level 0) and a fine level (refinement level 1). We start from a time t0, at which data on all refinement levels is consistent and valid.\n",
    "\n",
    "* All data is initialized at time t = t0\n",
    "* The coarse level is evolved forward to time to t = t0 + dt\n",
    "* The fine level is evolved forward to time to t = t0 + dt/2\n",
    "* Boundary data for the fine grid at time to t = t0+dt/2 is interpolated from the coarse grids at times t0 and t0+dt\n",
    "* The fine level is evolved forward in time to t = t0 + dt\n",
    "* Data from the fine level at t=t0+dt is copied back into the coarse grid at time t=t0+dt\n",
    "* If not done, t0 += dt and repeat."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wave_toy2_par_src = \"\"\"\n",
    "ActiveThorns = \"{thorn_name} coordbase carpet cartgrid3d MoL\"\n",
    "\n",
    "Cactus::cctk_itlast = 250\n",
    "\n",
    "# the output dir will be named after the parameter file name\n",
    "IO::out_dir = $parfile\n",
    "\n",
    "# The set of parameters required to specify the grid...\n",
    "Carpet::domain_from_coordbase = \"yes\"\n",
    "CartGrid3D::type = \"coordbase\"\n",
    "CartGrid3D::avoid_origin = \"no\"\n",
    "CoordBase::domainsize = \"minmax\"\n",
    "CoordBase::spacing    = \"gridspacing\"\n",
    "CoordBase::xmin =  0\n",
    "CoordBase::xmax =  30\n",
    "CoordBase::ymin =  0\n",
    "CoordBase::ymax =  30\n",
    "CoordBase::zmin =  0.0\n",
    "CoordBase::zmax =  0.0\n",
    "CoordBase::dx   =  0.25\n",
    "CoordBase::dy   =  0.25\n",
    "\n",
    "# Because we're doing a 2D code...\n",
    "CoordBase::boundary_size_z_lower     = 0\n",
    "CoordBase::boundary_size_z_upper     = 0\n",
    "CoordBase::boundary_shiftout_z_lower = 1\n",
    "CoordBase::boundary_shiftout_z_upper = 1\n",
    "\n",
    "# The time integration scheme...\n",
    "MoL::ODE_Method = \"RK3\"\n",
    "\n",
    "# The next section is required if you want some kind of confirmation\n",
    "# on the console that your code is actually doing anything...\n",
    "ActiveThorns = \"CarpetIOBasic IOUtil CarpetReduce SymBase\"\n",
    "IO::out_fileinfo=\"all\"\n",
    "IOBasic::outInfo_every = 1\n",
    "IOBasic::outInfo_vars = \"{thorn_name}::phi {thorn_name}::psi\"\n",
    "\n",
    "# Get some data!\n",
    "ActiveThorns = \"CarpetIOASCII\"\n",
    "CarpetIOASCII::compact_format = false\n",
    "IOASCII::out2D_every = 25\n",
    "IOASCII::out2D_xyplane_z = 0\n",
    "IOASCII::out2D_vars = \"{thorn_name}::psi {thorn_name}::phi\"\n",
    "IOASCII::out2D_xz = \"no\"\n",
    "IOASCII::out2D_yz = \"no\"\n",
    "IOASCII::output_ghost_points = \"no\"\n",
    "\n",
    "# Taking too big a time step can lead to unstable evolution.\n",
    "# You can control that with dtfac and thorn Time.\n",
    "ActiveThorns = \"Time\"\n",
    "Time::dtfac = .2\n",
    "\n",
    "### MESH REFINEMENT ####\n",
    "ActiveThorns = \"CarpetRegrid2\"\n",
    "Carpet::init_fill_timelevels = \"true\"\n",
    "Carpet::max_refinement_levels = 2\n",
    "CarpetRegrid2::num_centres = 1\n",
    "CarpetRegrid2::num_levels_1 = 2\n",
    "CarpetRegrid2::position_x_1 = 8.0\n",
    "CarpetRegrid2::position_y_1 = 8.0\n",
    "CarpetRegrid2::position_z_1 = 0.0\n",
    "CarpetRegrid2::radius_1[1] = 8.0\n",
    "\n",
    "#MoL::disable_prolongation        = \"no\"\n",
    "\n",
    "### Dissipation is numerical noise added to mitigate\n",
    "### grid boundary effects seen in some AMR simulations.\n",
    "#ActiveThorns = \"Dissipation Spacemask SPHERICALSURFACE\"\n",
    "#Dissipation::vars = \"{thorn_name}::psi {thorn_name}::phi\"\n",
    "#Dissipation::order = 1\n",
    "\"\"\"\n",
    "create_file(\"wave_toy2.par\",wave_toy2_par_src,thorn_pars)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!rm -fr ../simulations/wave_toy2\n",
    "!./simfactory/bin/sim create-run wave_toy2.par --procs=2 --num-threads=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file_data = np.genfromtxt(\"../simulations/wave_toy2/output-0000/wave_toy2/psi.xy.asc\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sets = np.unique(file_data[:,0])\n",
    "width = 8\n",
    "height = 4\n",
    "# NOTE: I changed the colorscheme here\n",
    "cmap = cm.prism\n",
    "print(\"sets=\",sets)\n",
    "\n",
    "mn, mx = np.min(file_data[:,12]),np.max(file_data[:,12])\n",
    "#mn, mx = 0,1\n",
    "for which in sets:\n",
    "    gg = file_data[file_data[:,0]==which,:]\n",
    "    print(\"time=\",which)\n",
    "    for rl in range(10): # refinement levels\n",
    "      g = gg[gg[:,2]==rl,:]\n",
    "      x = g[:,5]\n",
    "      y = g[:,6]\n",
    "      z = g[:,12]\n",
    "      if z.shape[0] == 0:\n",
    "        continue\n",
    "      print(\"refinement level =\",rl)\n",
    "      plt.figure(figsize=(width, height))\n",
    "      zi = z.reshape(len(np.unique(y)),len(np.unique(x)))\n",
    "      plt.imshow(zi[::-1,:],cmap,clim=(mn,mx))\n",
    "      plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3>Excercise:</h3> Play around with the CarpetRegrid2 parameters. Attempt to add more levels of refinement and/or more centers of refinement."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<table><tr><td>This work sponsored by NSF grants <a href=\"https://www.nsf.gov/awardsearch/showAward?AWD_ID=1550551\"> OAC 1550551</a> and <a href=\"https://www.nsf.gov/awardsearch/showAward?AWD_ID=1539567\"> CCF 1539567</a></td><td><img src=\"https://www.nsf.gov/awardsearch/images/common/nsf_logo_bottom.png\"></tr></table>"
   ]
  },
  {
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}