Tutorial-Start_to_Finish-ScalarWaveCurvilinear.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<script async src=\"https://www.googletagmanager.com/gtag/js?id=UA-59152712-8\"></script>\n",
    "<script>\n",
    "  window.dataLayer = window.dataLayer || [];\n",
    "  function gtag(){dataLayer.push(arguments);}\n",
    "  gtag('js', new Date());\n",
    "\n",
    "  gtag('config', 'UA-59152712-8');\n",
    "</script>\n",
    "\n",
    "# Start-to-Finish Example: Numerical Solution of the Scalar Wave Equation, in Curvilinear Coordinates\n",
    "\n",
    "## Author: Zach Etienne\n",
    "### Formatting improvements courtesy Brandon Clark\n",
    "\n",
    "## This module solves the scalar wave equation in *spherical coordinates* (though other coordinates, including Cartesian, may be chosen).\n",
    "\n",
    "**Notebook Status:** <font color =\"green\"><b> Validated </b></font>\n",
    "\n",
    "**Validation Notes:** This module has been validated to converge at the expected order to the exact solution (see [plot](#convergence) at bottom).\n",
    "\n",
    "### NRPy+ Source Code for this module: \n",
    "* [ScalarWave/ScalarWaveCurvilinear_RHSs.py](../edit/ScalarWave/ScalarWaveCurvilinear_RHSs.py) [\\[**tutorial**\\]](Tutorial-ScalarWaveCurvilinear.ipynb) Generates the right-hand side for the Scalar Wave Equation in curvilinear coordinates\n",
    "* [ScalarWave/InitialData.py](../edit/ScalarWave/InitialData.py) [\\[**tutorial**\\]](Tutorial-ScalarWave.ipynb) Generating C code for either plane wave or spherical Gaussian initial data for the scalar wave equation \n",
    "\n",
    "## Introduction:\n",
    "As outlined in the [previous NRPy+ tutorial notebook](Tutorial-ScalarWaveCurvilinear.ipynb), we first use NRPy+ to generate initial data for the scalar wave equation, and then we use it to generate the RHS expressions for [Method of Lines](https://reference.wolfram.com/language/tutorial/NDSolveMethodOfLines.html) time integration based on the [explicit Runge-Kutta fourth-order scheme](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods) (RK4).\n",
    "\n",
    "The entire algorithm is outlined below, with NRPy+-based components highlighted in <font color='green'>green</font>.\n",
    "\n",
    "1. Allocate memory for gridfunctions, including temporary storage for the RK4 time integration.\n",
    "1. <font color='green'>Set gridfunction values to initial data.</font>\n",
    "1. Evolve the system forward in time using RK4 time integration. At each RK4 substep, do the following:\n",
    "    1. <font color='green'>Evaluate scalar wave RHS expressions.</font>\n",
    "    1. Apply boundary conditions.\n",
    "1. At the end of each iteration in time, output the relative error between numerical and exact solutions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='toc'></a>\n",
    "\n",
    "# Table of Contents\n",
    "$$\\label{toc}$$\n",
    "\n",
    "This notebook is organized as follows\n",
    "\n",
    "1. [Step 1](#writec): Generate C code to solve the scalar wave equation in curvilinear coordinates\n",
    "    1. [Step 1.a](#id_rhss): C code generation: Initial data and scalar wave right-hand-sides\n",
    "    1. [Step 1.b](#boundaryconditions): C code generation: Boundary condition driver\n",
    "    1. [Step 1.c](#cparams_rfm_and_domainsize): Generate Cparameters files; set reference metric parameters, including `domain_size`\n",
    "    1. [Step 1.d](#cfl): C code generation: Finding the minimum proper distance between grid points, needed for [CFL](https://en.wikipedia.org/w/index.php?title=Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition&oldid=806430673)-limited timestep\n",
    "1. [Step 2](#mainc): `ScalarWaveCurvilinear_Playground.c`: The Main C Code\n",
    "1. [Step 3](#compileexec): Compile generated C codes & solve the scalar wave equation\n",
    "1. [Step 4](#convergence): Code validation: Plot the numerical error, and confirm that it converges to zero at expected rate with increasing numerical resolution (sampling)\n",
    "1. [Step 5](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='writec'></a>\n",
    "\n",
    "# Step 1: Using NRPy+ to generate necessary C code to solve the scalar wave equation in curvilinear, singular coordinates \\[Back to [top](#toc)\\]\n",
    "$$\\label{writec}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='id_rhss'></a>\n",
    "\n",
    "## Step 1.a: C code generation: Initial data and scalar wave RHSs \\[Back to [top](#toc)\\]\n",
    "$$\\label{id_rhss}$$\n",
    "\n",
    "\n",
    "We choose simple plane wave initial data, which is documented in the [Cartesian scalar wave module](Tutorial-ScalarWave.ipynb). Specifically, we implement monochromatic (single-wavelength) wave traveling in the $\\hat{k}$ direction with speed $c$\n",
    "$$u(\\vec{x},t) = f(\\hat{k}\\cdot\\vec{x} - c t),$$\n",
    "where $\\hat{k}$ is a unit vector.\n",
    "\n",
    "The scalar wave RHSs in curvilinear coordinates (documented [in the previous module](Tutorial-ScalarWaveCurvilinear.ipynb)) are simply the right-hand sides of the scalar wave equation written in curvilinear coordinates\n",
    "\\begin{align}\n",
    "\\partial_t u &= v \\\\\n",
    "\\partial_t v &= c^2 \\left(\\hat{g}^{ij} \\partial_{i} \\partial_{j} u - \\hat{\\Gamma}^i \\partial_i u\\right),\n",
    "\\end{align}\n",
    "where $\\hat{g}^{ij}$ is the inverse reference 3-metric (i.e., the metric corresponding to the underlying coordinate system we choose$-$spherical coordinates in our example below), and $\\hat{\\Gamma}^i$ is the contracted Christoffel symbol $\\hat{\\Gamma}^\\tau = \\hat{g}^{\\mu\\nu} \\hat{\\Gamma}^\\tau_{\\mu\\nu}$.\n",
    "\n",
    "Below we generate \n",
    "+ the initial data by calling `InitialData(Type=\"PlaneWave\")` inside the NRPy+ [ScalarWave/InitialData.py](../edit/ScalarWave/InitialData.py) module (documented in [this NRPy+ Jupyter notebook](Tutorial-ScalarWave.ipynb)), and \n",
    "+ the RHS expressions by calling `ScalarWaveCurvilinear_RHSs()` inside the NRPy+ [ScalarWave/ScalarWaveCurvilinear_RHSs.py](../edit/ScalarWave/ScalarWaveCurvilinear_RHSs.py) module (documented in [this NRPy+ Jupyter notebook](Tutorial-ScalarWaveCurvilinear.ipynb))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:32:26.686250Z",
     "iopub.status.busy": "2021-03-07T17:32:26.675699Z",
     "iopub.status.idle": "2021-03-07T17:32:30.350287Z",
     "shell.execute_reply": "2021-03-07T17:32:30.350799Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Output C function exact_solution_single_point() to file ScalarWaveCurvilinear_Playground_Ccodes/exact_solution_single_point.h\n",
      "Output C function exact_solution_all_points() to file ScalarWaveCurvilinear_Playground_Ccodes/exact_solution_all_points.h\n",
      "Output C function rhs_eval() to file ScalarWaveCurvilinear_Playground_Ccodes/rhs_eval.h\n"
     ]
    }
   ],
   "source": [
    "# Step P1: Import needed NRPy+ core modules:\n",
    "from outputC import lhrh,outCfunction # NRPy+: Core C code output module\n",
    "import finite_difference as fin       # NRPy+: Finite difference C code generation module\n",
    "import NRPy_param_funcs as par        # NRPy+: Parameter interface\n",
    "import grid as gri                    # NRPy+: Functions having to do with numerical grids\n",
    "import reference_metric as rfm        # NRPy+: Reference metric support\n",
    "import cmdline_helper as cmd          # NRPy+: Multi-platform Python command-line interface\n",
    "import shutil, os, sys                # Standard Python modules for multiplatform OS-level functions\n",
    "\n",
    "# Step P2: Create C code output directory:\n",
    "Ccodesdir = os.path.join(\"ScalarWaveCurvilinear_Playground_Ccodes/\")\n",
    "# First remove C code output directory if it exists\n",
    "# Courtesy https://stackoverflow.com/questions/303200/how-do-i-remove-delete-a-folder-that-is-not-empty\n",
    "shutil.rmtree(Ccodesdir, ignore_errors=True)\n",
    "# Then create a fresh directory\n",
    "cmd.mkdir(Ccodesdir)\n",
    "\n",
    "# Step P3: Create executable output directory:\n",
    "outdir = os.path.join(Ccodesdir,\"output/\")\n",
    "cmd.mkdir(outdir)\n",
    "\n",
    "# Step 1: Set the spatial dimension parameter\n",
    "#         to three this time, and then read\n",
    "#         the parameter as DIM.\n",
    "par.set_parval_from_str(\"grid::DIM\",3)\n",
    "DIM = par.parval_from_str(\"grid::DIM\")\n",
    "\n",
    "# Step 2: Set some core parameters, including CoordSystem, boundary condition,\n",
    "#                                             MoL, timestepping algorithm, FD order,\n",
    "#                                             floating point precision, and CFL factor:\n",
    "\n",
    "# Step 2.a: Set the coordinate system for the numerical grid\n",
    "# Choices are: Spherical, SinhSpherical, SinhSphericalv2, Cylindrical, SinhCylindrical,\n",
    "#              SymTP, SinhSymTP\n",
    "CoordSystem     = \"SinhSpherical\"\n",
    "par.set_parval_from_str(\"reference_metric::CoordSystem\",CoordSystem)\n",
    "rfm.reference_metric()\n",
    "\n",
    "# Step 2.b: Set defaults for Coordinate system parameters.\n",
    "#           These are perhaps the most commonly adjusted parameters,\n",
    "#           so we enable modifications at this high level.\n",
    "\n",
    "# domain_size sets the default value for:\n",
    "#   * Spherical's params.RMAX\n",
    "#   * SinhSpherical*'s params.AMAX\n",
    "#   * Cartesians*'s -params.{x,y,z}min & .{x,y,z}max\n",
    "#   * Cylindrical's -params.ZMIN & .{Z,RHO}MAX\n",
    "#   * SinhCylindrical's params.AMPL{RHO,Z}\n",
    "#   * *SymTP's params.AMAX\n",
    "domain_size     = 10.0 # Needed for all coordinate systems.\n",
    "\n",
    "# sinh_width sets the default value for:\n",
    "#   * SinhSpherical's params.SINHW\n",
    "#   * SinhCylindrical's params.SINHW{RHO,Z}\n",
    "#   * SinhSymTP's params.SINHWAA\n",
    "sinh_width      = 0.4 # If Sinh* coordinates chosen\n",
    "\n",
    "# sinhv2_const_dr sets the default value for:\n",
    "#   * SinhSphericalv2's params.const_dr\n",
    "#   * SinhCylindricalv2's params.const_d{rho,z}\n",
    "sinhv2_const_dr = 0.05# If Sinh*v2 coordinates chosen\n",
    "\n",
    "# SymTP_bScale sets the default value for:\n",
    "#   * SinhSymTP's params.bScale\n",
    "SymTP_bScale    = 1.0 # If SymTP chosen\n",
    "\n",
    "# Step 2.c: Set the order of spatial and temporal derivatives;\n",
    "#           the core data type, and the CFL factor.\n",
    "# RK_method choices include: Euler, \"RK2 Heun\", \"RK2 MP\", \"RK2 Ralston\", RK3, \"RK3 Heun\", \"RK3 Ralston\",\n",
    "#              SSPRK3, RK4, DP5, DP5alt, CK5, DP6, L6, DP8\n",
    "RK_method = \"RK4\"\n",
    "FD_order  = 4        # Finite difference order: even numbers only, starting with 2. 12 is generally unstable\n",
    "REAL      = \"double\" # Best to use double here.\n",
    "CFL_FACTOR= 1.0\n",
    "\n",
    "# Step 3: Generate Runge-Kutta-based (RK-based) timestepping code.\n",
    "#         Each RK substep involves two function calls:\n",
    "#      3.A: Evaluate RHSs (RHS_string)\n",
    "#      3.B: Apply boundary conditions (post_RHS_string)\n",
    "import MoLtimestepping.C_Code_Generation as MoL\n",
    "from MoLtimestepping.RK_Butcher_Table_Dictionary import Butcher_dict\n",
    "RK_order  = Butcher_dict[RK_method][1]\n",
    "cmd.mkdir(os.path.join(Ccodesdir,\"MoLtimestepping/\"))\n",
    "\n",
    "RHS_string = \"rhs_eval(&rfmstruct, &params, RK_INPUT_GFS, RK_OUTPUT_GFS);\"\n",
    "\n",
    "post_RHS_string = \"apply_bcs_curvilinear(&params, &bcstruct, NUM_EVOL_GFS, evol_gf_parity, RK_OUTPUT_GFS);\"\n",
    "\n",
    "MoL.MoL_C_Code_Generation(RK_method, RHS_string = RHS_string, post_RHS_string = post_RHS_string,\n",
    "                          outdir = os.path.join(Ccodesdir,\"MoLtimestepping/\"))\n",
    "\n",
    "# Step 4: Import the ScalarWave.InitialData module.\n",
    "#         This command only declares ScalarWave initial data\n",
    "#         parameters and the InitialData() function.\n",
    "import ScalarWave.InitialData as swid\n",
    "\n",
    "# Step 5: Import ScalarWave_RHSs module.\n",
    "#         This command only declares ScalarWave RHS parameters\n",
    "#         and the ScalarWave_RHSs function (called later)\n",
    "import ScalarWave.ScalarWaveCurvilinear_RHSs as swrhs\n",
    "\n",
    "# Step 6: Set the finite differencing order to FD_order (set above).\n",
    "par.set_parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\",FD_order)\n",
    "\n",
    "# Step 7: Call the InitialData() function to set up initial data.\n",
    "#         Options include:\n",
    "#    \"PlaneWave\": monochromatic (single frequency/wavelength) plane wave\n",
    "#    \"SphericalGaussian\": spherically symmetric Gaussian, with default stdev=3\n",
    "swid.InitialData(CoordSystem=CoordSystem,Type=\"PlaneWave\")\n",
    "\n",
    "# Step 8: Generate SymPy symbolic expressions for\n",
    "#         uu_rhs and vv_rhs; the ScalarWave RHSs.\n",
    "#         This function also declares the uu and vv\n",
    "#         gridfunctions, which need to be declared\n",
    "#         to output even the initial data to C file.\n",
    "cmd.mkdir(os.path.join(Ccodesdir,\"rfm_files/\"))\n",
    "par.set_parval_from_str(\"reference_metric::enable_rfm_precompute\",\"True\")\n",
    "par.set_parval_from_str(\"reference_metric::rfm_precompute_Ccode_outdir\",os.path.join(Ccodesdir,\"rfm_files/\"))\n",
    "swrhs.ScalarWaveCurvilinear_RHSs()\n",
    "# Step 8.a: Now that we are finished with all the rfm hatted\n",
    "#           quantities, let's restore them to their closed-\n",
    "#           form expressions.\n",
    "par.set_parval_from_str(\"reference_metric::enable_rfm_precompute\",\"False\") # Reset to False to disable rfm_precompute.\n",
    "rfm.ref_metric__hatted_quantities()\n",
    "\n",
    "# Step 9: Copy SIMD/SIMD_intrinsics.h to $Ccodesdir/SIMD/SIMD_intrinsics.h\n",
    "cmd.mkdir(os.path.join(Ccodesdir,\"SIMD\"))\n",
    "shutil.copy(os.path.join(\"SIMD/\")+\"SIMD_intrinsics.h\",os.path.join(Ccodesdir,\"SIMD/\"))\n",
    "\n",
    "# Step 10: Generate all needed C functions\n",
    "enable_FD_functions = False\n",
    "par.set_parval_from_str(\"finite_difference::enable_FD_functions\",enable_FD_functions)\n",
    "\n",
    "desc=\"Part P3: Declare the function for the exact solution at a single point. time==0 corresponds to the initial data.\"\n",
    "name=\"exact_solution_single_point\"\n",
    "outCfunction(\n",
    "    outfile  = os.path.join(Ccodesdir,name+\".h\"), desc=desc, name=name,\n",
    "    params   =\"const REAL xx0,const REAL xx1,const REAL xx2,const paramstruct *restrict params,REAL *uu_exact,REAL *vv_exact\",\n",
    "    body     = fin.FD_outputC(\"returnstring\",[lhrh(lhs=\"*uu_exact\",rhs=swid.uu_ID),\n",
    "                                              lhrh(lhs=\"*vv_exact\",rhs=swid.vv_ID)]),\n",
    "    loopopts = \"\")\n",
    "\n",
    "desc=\"Part P4: Declare the function for the exact solution at all points. time==0 corresponds to the initial data.\"\n",
    "name=\"exact_solution_all_points\"\n",
    "outCfunction(\n",
    "    outfile  = os.path.join(Ccodesdir,name+\".h\"), desc=desc, name=name,\n",
    "    params   =\"const paramstruct *restrict params,REAL *restrict xx[3], REAL *restrict in_gfs\",\n",
    "    body     =\"\"\"exact_solution_single_point(xx[0][i0],xx[1][i1],xx[2][i2],params,\n",
    "                            &in_gfs[IDX4S(UUGF,i0,i1,i2)],&in_gfs[IDX4S(VVGF,i0,i1,i2)]);\"\"\",\n",
    "    loopopts = \"AllPoints\")\n",
    "\n",
    "desc=\"Part P5: Declare the function to evaluate the scalar wave RHSs\"\n",
    "includes = None\n",
    "if enable_FD_functions:\n",
    "    includes = [\"finite_difference_functions.h\"]\n",
    "name=\"rhs_eval\"\n",
    "outCfunction(\n",
    "    outfile  = os.path.join(Ccodesdir,name+\".h\"), includes=includes, desc=desc, name=name,\n",
    "    params   =\"\"\"rfm_struct *restrict rfmstruct,const paramstruct *restrict params,\n",
    "                 const REAL *restrict in_gfs, REAL *restrict rhs_gfs\"\"\",\n",
    "    body     =fin.FD_outputC(\"returnstring\",[lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"uu\"),rhs=swrhs.uu_rhs),\n",
    "                                             lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"vv\"),rhs=swrhs.vv_rhs)],\n",
    "                             params=\"enable_SIMD=True\"),\n",
    "    loopopts = \"InteriorPoints,enable_SIMD,enable_rfm_precompute\")\n",
    "\n",
    "# Step 10.b Output functions for computing all finite-difference stencils\n",
    "if enable_FD_functions:\n",
    "    fin.output_finite_difference_functions_h(path=Ccodesdir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='boundaryconditions'></a>\n",
    "\n",
    "## Step 1.b: Output needed C code for boundary condition driver \\[Back to [top](#toc)\\]\n",
    "$$\\label{boundaryconditions}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:32:30.356745Z",
     "iopub.status.busy": "2021-03-07T17:32:30.356078Z",
     "iopub.status.idle": "2021-03-07T17:32:31.543735Z",
     "shell.execute_reply": "2021-03-07T17:32:31.544260Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wrote to file \"ScalarWaveCurvilinear_Playground_Ccodes/boundary_conditions/parity_conditions_symbolic_dot_products.h\"\n",
      "Evolved parity: ( uu:0, vv:0 )\n",
      "\n",
      "\n",
      "Wrote to file \"ScalarWaveCurvilinear_Playground_Ccodes/boundary_conditions/EigenCoord_Cart_to_xx.h\"\n"
     ]
    }
   ],
   "source": [
    "import CurviBoundaryConditions.CurviBoundaryConditions as cbcs\n",
    "cbcs.Set_up_CurviBoundaryConditions(os.path.join(Ccodesdir,\"boundary_conditions/\"),\n",
    "                                    Cparamspath=os.path.join(\"../\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='cparams_rfm_and_domainsize'></a>\n",
    "\n",
    "## Step 1.c: Output C codes needed for declaring and setting Cparameters; also set `free_parameters.h` \\[Back to [top](#toc)\\]\n",
    "$$\\label{cparams_rfm_and_domainsize}$$\n",
    "\n",
    "Based on declared NRPy+ Cparameters, first we generate `declare_Cparameters_struct.h`, `set_Cparameters_default.h`, and `set_Cparameters[-SIMD].h`.\n",
    "\n",
    "Then we output `free_parameters.h`, which sets initial data parameters, as well as grid domain & reference metric parameters, applying `domain_size` and `sinh_width`/`SymTP_bScale` (if applicable) as set above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:32:31.550970Z",
     "iopub.status.busy": "2021-03-07T17:32:31.550257Z",
     "iopub.status.idle": "2021-03-07T17:32:32.120970Z",
     "shell.execute_reply": "2021-03-07T17:32:32.121447Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 1.c.i: Set free_parameters.h\n",
    "with open(os.path.join(Ccodesdir,\"free_parameters.h\"),\"w\") as file:\n",
    "    file.write(\"\"\"\n",
    "// Set free-parameter values.\n",
    "params.time = 0.0; // Initial simulation time time corresponds to exact solution at time=0.\n",
    "params.wavespeed = 1.0;\\n\"\"\")\n",
    "\n",
    "# Append to $Ccodesdir/free_parameters.h reference metric parameters based on generic\n",
    "#    domain_size,sinh_width,sinhv2_const_dr,SymTP_bScale,\n",
    "#    parameters set above.\n",
    "rfm.out_default_free_parameters_for_rfm(os.path.join(Ccodesdir,\"free_parameters.h\"),\n",
    "                                        domain_size,sinh_width,sinhv2_const_dr,SymTP_bScale)\n",
    "\n",
    "# Step 1.c.ii: Generate set_Nxx_dxx_invdx_params__and__xx.h:\n",
    "rfm.set_Nxx_dxx_invdx_params__and__xx_h(os.path.join(Ccodesdir))\n",
    "\n",
    "# Step 1.c.iii: Generate xx_to_Cart.h, which contains xx_to_Cart() for\n",
    "#               (the mapping from xx->Cartesian) for the chosen\n",
    "#               CoordSystem:\n",
    "rfm.xx_to_Cart_h(\"xx_to_Cart\",\"./set_Cparameters.h\",os.path.join(Ccodesdir,\"xx_to_Cart.h\"))\n",
    "\n",
    "# Step 1.c.iv: Generate declare_Cparameters_struct.h, set_Cparameters_default.h, and set_Cparameters[].h\n",
    "par.generate_Cparameters_Ccodes(os.path.join(Ccodesdir))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='cfl'></a>\n",
    "\n",
    "## Step 1.d: Output needed C code for finding the minimum proper distance between grid points, needed for [CFL](https://en.wikipedia.org/w/index.php?title=Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition&oldid=806430673)-limited timestep \\[Back to [top](#toc)\\]\n",
    "$$\\label{cfl}$$\n",
    "\n",
    "In order for our explicit-timestepping numerical solution to the scalar wave equation to be stable, it must satisfy the [CFL](https://en.wikipedia.org/w/index.php?title=Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition&oldid=806430673) condition:\n",
    "$$\n",
    "\\Delta t \\le \\frac{\\min(ds_i)}{c},\n",
    "$$\n",
    "where $c$ is the wavespeed, and\n",
    "$$ds_i = h_i \\Delta x^i$$ \n",
    "is the proper distance between neighboring gridpoints in the $i$th direction (in 3D, there are 3 directions), $h_i$ is the $i$th reference metric scale factor, and $\\Delta x^i$ is the uniform grid spacing in the $i$th direction:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:32:32.143346Z",
     "iopub.status.busy": "2021-03-07T17:32:32.142704Z",
     "iopub.status.idle": "2021-03-07T17:32:32.145242Z",
     "shell.execute_reply": "2021-03-07T17:32:32.145786Z"
    }
   },
   "outputs": [],
   "source": [
    "# Output the find_timestep() function to a C file.\n",
    "rfm.out_timestep_func_to_file(os.path.join(Ccodesdir,\"find_timestep.h\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='mainc'></a>\n",
    "\n",
    "# Step 2: `ScalarWaveCurvilinear_Playground.c`: The Main C Code \\[Back to [top](#toc)\\]\n",
    "$$\\label{mainc}$$\n",
    "\n",
    "Just as in [the start-to-finish, solving the scalar wave equation in Cartesian coordinates module](Tutorial-Start_to_Finish-ScalarWave.ipynb), we will implement the scalar wave equation via the Method of Lines. As discussed above, the critical differences between this code and the Cartesian version are as follows:\n",
    "1. The CFL-constrained timestep depends on the proper distance between neighboring gridpoints\n",
    "1. The boundary conditions must account for the fact that ghost zone points lying in the domain exterior can map either to the interior of the domain, or lie on the outer boundary. In the former case, we simply copy the data from the interior. In the latter case, we apply the usual outer boundary conditions.\n",
    "1. The numerical grids must be staggered to avoid direct evaluation of the equations on coordinate singularities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:32:32.151336Z",
     "iopub.status.busy": "2021-03-07T17:32:32.150666Z",
     "iopub.status.idle": "2021-03-07T17:32:32.152988Z",
     "shell.execute_reply": "2021-03-07T17:32:32.153458Z"
    }
   },
   "outputs": [],
   "source": [
    "# Part P0: Define REAL, set the number of ghost cells NGHOSTS (from NRPy+'s FD_CENTDERIVS_ORDER),\n",
    "#          and set the CFL_FACTOR (which can be overwritten at the command line)\n",
    "\n",
    "with open(os.path.join(Ccodesdir,\"ScalarWaveCurvilinear_Playground_REAL__NGHOSTS__CFL_FACTOR.h\"), \"w\") as file:\n",
    "    file.write(\"\"\"\n",
    "// Part P0.a: Set the number of ghost cells, from NRPy+'s FD_CENTDERIVS_ORDER\n",
    "#define NGHOSTS \"\"\"+str(int(FD_order/2))+\"\"\"\n",
    "// Part P0.b: Set the numerical precision (REAL) to double, ensuring all floating point\n",
    "//            numbers are stored to at least ~16 significant digits\n",
    "#define REAL \"\"\"+REAL+\"\"\"\n",
    "// Part P0.c: Set the number of ghost cells, from NRPy+'s FD_CENTDERIVS_ORDER\n",
    "REAL CFL_FACTOR = \"\"\"+str(CFL_FACTOR)+\";\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:32:32.163854Z",
     "iopub.status.busy": "2021-03-07T17:32:32.160067Z",
     "iopub.status.idle": "2021-03-07T17:32:32.166402Z",
     "shell.execute_reply": "2021-03-07T17:32:32.165694Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing ScalarWaveCurvilinear_Playground_Ccodes//ScalarWaveCurvilinear_Playground.c\n"
     ]
    }
   ],
   "source": [
    "%%writefile $Ccodesdir/ScalarWaveCurvilinear_Playground.c\n",
    "\n",
    "// Step P0: Define REAL and NGHOSTS; and declare CFL_FACTOR. This header is generated in NRPy+.\n",
    "#include \"ScalarWaveCurvilinear_Playground_REAL__NGHOSTS__CFL_FACTOR.h\"\n",
    "\n",
    "#include \"rfm_files/rfm_struct__declare.h\"\n",
    "\n",
    "#include \"declare_Cparameters_struct.h\"\n",
    "\n",
    "// All SIMD intrinsics used in SIMD-enabled C code loops are defined here:\n",
    "#include \"SIMD/SIMD_intrinsics.h\"\n",
    "\n",
    "// Step P1: Import needed header files\n",
    "#include \"stdio.h\"\n",
    "#include \"stdlib.h\"\n",
    "#include \"math.h\"\n",
    "#include \"stdint.h\" // Needed for Windows GCC 6.x compatibility\n",
    "#ifndef M_PI\n",
    "#define M_PI 3.141592653589793238462643383279502884L\n",
    "#endif\n",
    "#ifndef M_SQRT1_2\n",
    "#define M_SQRT1_2 0.707106781186547524400844362104849039L\n",
    "#endif\n",
    "\n",
    "// Step P2: Declare the IDX4S(gf,i,j,k) macro, which enables us to store 4-dimensions of\n",
    "//           data in a 1D array. In this case, consecutive values of \"i\"\n",
    "//           (all other indices held to a fixed value) are consecutive in memory, where\n",
    "//           consecutive values of \"j\" (fixing all other indices) are separated by\n",
    "//           Nxx_plus_2NGHOSTS0 elements in memory. Similarly, consecutive values of\n",
    "//           \"k\" are separated by Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1 in memory, etc.\n",
    "#define IDX4S(g,i,j,k) \\\n",
    "( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * ( (k) + Nxx_plus_2NGHOSTS2 * (g) ) ) )\n",
    "#define IDX3S(i,j,k) ( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * ( (k) ) ) )\n",
    "#define LOOP_REGION(i0min,i0max, i1min,i1max, i2min,i2max) \\\n",
    "  for(int i2=i2min;i2<i2max;i2++) for(int i1=i1min;i1<i1max;i1++) for(int i0=i0min;i0<i0max;i0++)\n",
    "#define LOOP_ALL_GFS_GPS(ii) _Pragma(\"omp parallel for\") \\\n",
    "  for(int (ii)=0;(ii)<Nxx_plus_2NGHOSTS_tot*NUM_EVOL_GFS;(ii)++)\n",
    "\n",
    "// Step P3: Set UUGF and VVGF macros, as well as xx_to_Cart()\n",
    "#include \"boundary_conditions/gridfunction_defines.h\"\n",
    "\n",
    "// Step P4: Set xx_to_Cart(const paramstruct *restrict params,\n",
    "//                     REAL *restrict xx[3],\n",
    "//                     const int i0,const int i1,const int i2,\n",
    "//                     REAL xCart[3]),\n",
    "//           which maps xx->Cartesian via\n",
    "//    {xx[0][i0],xx[1][i1],xx[2][i2]}->{xCart[0],xCart[1],xCart[2]}\n",
    "#include \"xx_to_Cart.h\"\n",
    "\n",
    "// Step P5: Defines set_Nxx_dxx_invdx_params__and__xx(const int EigenCoord, const int Nxx[3],\n",
    "//                                       paramstruct *restrict params, REAL *restrict xx[3]),\n",
    "//          which sets params Nxx,Nxx_plus_2NGHOSTS,dxx,invdx, and xx[] for\n",
    "//          the chosen Eigen-CoordSystem if EigenCoord==1, or\n",
    "//          CoordSystem if EigenCoord==0.\n",
    "#include \"set_Nxx_dxx_invdx_params__and__xx.h\"\n",
    "\n",
    "// Step P6: Include basic functions needed to impose curvilinear\n",
    "//          parity and boundary conditions.\n",
    "#include \"boundary_conditions/CurviBC_include_Cfunctions.h\"\n",
    "\n",
    "// Step P7: Find the CFL-constrained timestep\n",
    "#include \"find_timestep.h\"\n",
    "\n",
    "// Part P8: Declare the function for the exact solution at a single point. time==0 corresponds to the initial data.\n",
    "#include \"exact_solution_single_point.h\"\n",
    "\n",
    "// Part P9: Declare the function for the exact solution at all points. time==0 corresponds to the initial data.\n",
    "#include \"exact_solution_all_points.h\"\n",
    "\n",
    "// Part P10: Declare the function to evaluate the scalar wave RHSs\n",
    "#include \"rhs_eval.h\"\n",
    "\n",
    "// main() function:\n",
    "// Step 0: Read command-line input, set up grid structure, allocate memory for gridfunctions, set up coordinates\n",
    "// Step 1: Set up scalar wave initial data\n",
    "// Step 2: Output relative error between numerical and exact solution.\n",
    "// Step 3: Evolve scalar wave initial data forward in time using Method of Lines with chosen RK-like algorithm,\n",
    "//         applying quadratic extrapolation outer boundary conditions.\n",
    "// Step 4: Free all allocated memory\n",
    "int main(int argc, const char *argv[]) {\n",
    "    paramstruct params;\n",
    "#include \"set_Cparameters_default.h\"\n",
    "\n",
    "    // Step 0a: Read command-line input, error out if nonconformant\n",
    "    if(argc != 4 || atoi(argv[1]) < NGHOSTS || atoi(argv[2]) < NGHOSTS || atoi(argv[3]) < NGHOSTS) {\n",
    "        printf(\"Error: Expected one command-line argument: ./ScalarWaveCurvilinear_Playground Nx0 Nx1 Nx2,\\n\");\n",
    "        printf(\"where Nx[0,1,2] is the number of grid points in the 0, 1, and 2 directions.\\n\");\n",
    "        printf(\"Nx[] MUST BE larger than NGHOSTS (= %d)\\n\",NGHOSTS);\n",
    "        exit(1);\n",
    "    }\n",
    "    // Step 0b: Set up numerical grid structure, first in space...\n",
    "    const int Nxx[3] = { atoi(argv[1]), atoi(argv[2]), atoi(argv[3]) };\n",
    "    if(Nxx[0]%2 != 0 || Nxx[1]%2 != 0 || Nxx[2]%2 != 0) {\n",
    "        printf(\"Error: Cannot guarantee a proper cell-centered grid if number of grid cells not set to even number.\\n\");\n",
    "        printf(\"       For example, in case of angular directions, proper symmetry zones will not exist.\\n\");\n",
    "        exit(1);\n",
    "    }\n",
    "\n",
    "    // Step 0c: Set free parameters, overwriting Cparameters defaults\n",
    "    //          by hand or with command-line input, as desired.\n",
    "#include \"free_parameters.h\"\n",
    "\n",
    "   // Step 0d: Uniform coordinate grids are stored to *xx[3]\n",
    "    REAL *xx[3];\n",
    "    // Step 0d.i: Set bcstruct\n",
    "    bc_struct bcstruct;\n",
    "    {\n",
    "        int EigenCoord = 1;\n",
    "        // Step 0d.ii: Call set_Nxx_dxx_invdx_params__and__xx(), which sets\n",
    "        //             params Nxx,Nxx_plus_2NGHOSTS,dxx,invdx, and xx[] for the\n",
    "        //             chosen Eigen-CoordSystem.\n",
    "        set_Nxx_dxx_invdx_params__and__xx(EigenCoord, Nxx, &params, xx);\n",
    "        // Step 0d.iii: Set Nxx_plus_2NGHOSTS_tot\n",
    "#include \"set_Cparameters-nopointer.h\"\n",
    "        const int Nxx_plus_2NGHOSTS_tot = Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2;\n",
    "        // Step 0e: Find ghostzone mappings; set up bcstruct\n",
    "#include \"boundary_conditions/driver_bcstruct.h\"\n",
    "        // Step 0e.i: Free allocated space for xx[][] array\n",
    "        for(int i=0;i<3;i++) free(xx[i]);\n",
    "    }\n",
    "\n",
    "    // Step 0f: Call set_Nxx_dxx_invdx_params__and__xx(), which sets\n",
    "    //          params Nxx,Nxx_plus_2NGHOSTS,dxx,invdx, and xx[] for the\n",
    "    //          chosen (non-Eigen) CoordSystem.\n",
    "    int EigenCoord = 0;\n",
    "    set_Nxx_dxx_invdx_params__and__xx(EigenCoord, Nxx, &params, xx);\n",
    "\n",
    "    // Step 0g: Set all C parameters \"blah\" for params.blah, including\n",
    "    //          Nxx_plus_2NGHOSTS0 = params.Nxx_plus_2NGHOSTS0, etc.\n",
    "#include \"set_Cparameters-nopointer.h\"\n",
    "    const int Nxx_plus_2NGHOSTS_tot = Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2;\n",
    "\n",
    "    // Step 0h: Time coordinate parameters\n",
    "    const REAL t_final =  0.7*domain_size; /* Final time is set so that at t=t_final,\n",
    "                                            * data at the origin have not been corrupted\n",
    "                                            * by the approximate outer boundary condition */\n",
    "\n",
    "    // Step 0i: Set timestep based on smallest proper distance between gridpoints and CFL factor\n",
    "    REAL dt = find_timestep(&params, xx);\n",
    "    //printf(\"# Timestep set to = %e\\n\",(double)dt);\n",
    "    int N_final = (int)(t_final / dt + 0.5); // The number of points in time.\n",
    "                                             // Add 0.5 to account for C rounding down\n",
    "                                             // typecasts to integers.\n",
    "    int output_every_N = (int)((REAL)N_final/800.0);\n",
    "    if(output_every_N == 0) output_every_N = 1;\n",
    "\n",
    "    // Step 0j: Error out if the number of auxiliary gridfunctions outnumber evolved gridfunctions.\n",
    "    //              This is a limitation of the RK method. You are always welcome to declare & allocate\n",
    "    //              additional gridfunctions by hand.\n",
    "    if(NUM_AUX_GFS > NUM_EVOL_GFS) {\n",
    "        printf(\"Error: NUM_AUX_GFS > NUM_EVOL_GFS. Either reduce the number of auxiliary gridfunctions,\\n\");\n",
    "        printf(\"       or allocate (malloc) by hand storage for *diagnostic_output_gfs. \\n\");\n",
    "        exit(1);\n",
    "    }\n",
    "\n",
    "    // Step 0k: Allocate memory for gridfunctions\n",
    "#include \"MoLtimestepping/RK_Allocate_Memory.h\"\n",
    "\n",
    "    // Step 0l: Set up precomputed reference metric arrays\n",
    "    // Step 0l.i: Allocate space for precomputed reference metric arrays.\n",
    "#include \"rfm_files/rfm_struct__malloc.h\"\n",
    "\n",
    "    // Step 0l.ii: Define precomputed reference metric arrays.\n",
    "    {\n",
    "#include \"set_Cparameters-nopointer.h\"\n",
    "#include \"rfm_files/rfm_struct__define.h\"\n",
    "    }\n",
    "\n",
    "    // Step 1: Set up initial data to be exact solution at time=0:\n",
    "    params.time = 0.0; exact_solution_all_points(&params, xx, y_n_gfs);\n",
    "\n",
    "    for(int n=0;n<=N_final;n++)\n",
    "       { // Main loop to progress forward in time.\n",
    "\n",
    "        // Step 1a: Set current time to correct value & compute exact solution\n",
    "        params.time = ((REAL)n)*dt;\n",
    "\n",
    "        // Step 2: Code validation: Compute log of L2 norm of difference\n",
    "        //         between numerical and exact solutions:\n",
    "        //   log_L2_Norm = log10( sqrt[Integral( [numerical - exact]^2 * dV)] ),\n",
    "        //         where integral is within 30% of the grid outer boundary (domain_size)\n",
    "        if(n%output_every_N == 0) {\n",
    "            REAL integral = 0.0;\n",
    "            REAL numpts   = 0.0;\n",
    "#pragma omp parallel for reduction(+:integral,numpts)\n",
    "            LOOP_REGION(NGHOSTS,Nxx_plus_2NGHOSTS0-NGHOSTS,\n",
    "                        NGHOSTS,Nxx_plus_2NGHOSTS1-NGHOSTS,\n",
    "                        NGHOSTS,Nxx_plus_2NGHOSTS2-NGHOSTS) {\n",
    "                REAL xCart[3]; xx_to_Cart(&params,xx,i0,i1,i2, xCart);\n",
    "                if(sqrt(xCart[0]*xCart[0] + xCart[1]*xCart[1] + xCart[2]*xCart[2]) < domain_size*0.3) {\n",
    "                    REAL uu_exact,vv_exact; exact_solution_single_point(xx[0][i0],xx[1][i1],xx[2][i2],&params,\n",
    "                                                                        &uu_exact,&vv_exact);\n",
    "                    double num   = (double)y_n_gfs[IDX4S(UUGF,i0,i1,i2)];\n",
    "                    double exact = (double)uu_exact;\n",
    "                    integral += (num - exact)*(num - exact);\n",
    "                    numpts   += 1.0;\n",
    "                }\n",
    "            }\n",
    "            // Compute and output the log of the L2 norm.\n",
    "            REAL log_L2_Norm = log10(sqrt(integral/numpts));\n",
    "            printf(\"%e %e\\n\",(double)params.time,log_L2_Norm);\n",
    "        }\n",
    "\n",
    "        // Step 3: Step forward one timestep (t -> t+dt) in time using\n",
    "        //           chosen RK-like MoL timestepping algorithm\n",
    "#include \"MoLtimestepping/RK_MoL.h\"\n",
    "\n",
    "    } // End main loop to progress forward in time.\n",
    "\n",
    "    // Step 4: Free all allocated memory\n",
    "#include \"rfm_files/rfm_struct__freemem.h\"\n",
    "#include \"boundary_conditions/bcstruct_freemem.h\"\n",
    "#include \"MoLtimestepping/RK_Free_Memory.h\"\n",
    "    for(int i=0;i<3;i++) free(xx[i]);\n",
    "    return 0;\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='compileexec'></a>\n",
    "\n",
    "# Step 3: Compile generated C codes & solve the scalar wave equation \\[Back to [top](#toc)\\]\n",
    "$$\\label{compileexec}$$\n",
    "\n",
    "To aid in the cross-platform-compatible (with Windows, MacOS, & Linux) compilation and execution, we make use of `cmdline_helper` [(**Tutorial**)](Tutorial-cmdline_helper.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:32:32.177447Z",
     "iopub.status.busy": "2021-03-07T17:32:32.176623Z",
     "iopub.status.idle": "2021-03-07T17:32:36.302885Z",
     "shell.execute_reply": "2021-03-07T17:32:36.304046Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Compiling executable...\n",
      "(EXEC): Executing `gcc -std=gnu99 -Ofast -fopenmp -march=native -funroll-loops ScalarWaveCurvilinear_Playground_Ccodes/ScalarWaveCurvilinear_Playground.c -o ScalarWaveCurvilinear_Playground_Ccodes/output/ScalarWaveCurvilinear_Playground -lm`...\n",
      "(BENCH): Finished executing in 0.8029422760009766 seconds.\n",
      "Finished compilation.\n",
      "(EXEC): Executing `taskset -c 0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15 ./ScalarWaveCurvilinear_Playground 16  8 16`...\n",
      "(BENCH): Finished executing in 0.20212268829345703 seconds.\n",
      "(EXEC): Executing `taskset -c 0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15 ./ScalarWaveCurvilinear_Playground 24 12 24`...\n",
      "(BENCH): Finished executing in 0.6029894351959229 seconds.\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd\n",
    "\n",
    "cmd.C_compile(os.path.join(Ccodesdir,\"ScalarWaveCurvilinear_Playground.c\"),\n",
    "              os.path.join(outdir,\"ScalarWaveCurvilinear_Playground\"),compile_mode=\"optimized\")\n",
    "# !clang -Ofast -fopenmp -mavx2 -mfma ScalarWave/ScalarWaveCurvilinear_Playground.c -o ScalarWaveCurvilinear_Playground -lm\n",
    "# !icc -align -qopenmp -xHost -O2 -qopt-report=5 -qopt-report-phase ipo -qopt-report-phase vec -vec-threshold1 -qopt-prefetch=4 ScalarWave/ScalarWaveCurvilinear_Playground.c -o ScalarWaveCurvilinear_Playground\n",
    "# !gcc-7 -Ofast -fopenmp -march=native ScalarWave/ScalarWaveCurvilinear_Playground.c -o ScalarWaveCurvilinear_Playground -lm\n",
    "\n",
    "# Change to output directory\n",
    "os.chdir(outdir)\n",
    "# Clean up existing output files\n",
    "cmd.delete_existing_files(\"out-*resolution.txt\")\n",
    "# Run executable\n",
    "if par.parval_from_str(\"reference_metric::CoordSystem\") == \"Cartesian\":\n",
    "    cmd.Execute(\"ScalarWaveCurvilinear_Playground\", \"16 16 16\", \"out-lowresolution.txt\")\n",
    "    cmd.Execute(\"ScalarWaveCurvilinear_Playground\", \"24 24 24\", \"out-medresolution.txt\")\n",
    "else:\n",
    "    cmd.Execute(\"ScalarWaveCurvilinear_Playground\", \"16  8 16\", \"out-lowresolution.txt\")\n",
    "    # 4.28s with icc and FD order = 10.\n",
    "    cmd.Execute(\"ScalarWaveCurvilinear_Playground\", \"24 12 24\", \"out-medresolution.txt\")\n",
    "    ########################################\n",
    "    # BENCHMARK 48x24x48 RUN, FD order = 4. desktop: 17.33s\n",
    "    # laptop: 51.82s on icc. 45.02s on GCC 9, 45.03s on GCC 7, 51.67s on clang\n",
    "    #     cmd.Execute(\"ScalarWaveCurvilinear_Playground\", \"48 24 48\", \"out-hghresolution.txt\")\n",
    "\n",
    "# %timeit cmd.Execute(\"ScalarWaveCurvilinear_Playground\", \"48 24 48\", \"out-hghresolution.txt\", verbose=False)\n",
    "# FD functions disabled:\n",
    "# 16 s ± 702 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
    "# FD functions enabled:\n",
    "# 16.1 s ± 384 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
    "\n",
    "# Return to root directory\n",
    "os.chdir(os.path.join(\"../../\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='convergence'></a>\n",
    "\n",
    "# Step 4: Code validation: Plot the numerical error, and confirm that it converges to zero at expected rate with increasing numerical resolution (sampling) \\[Back to [top](#toc)\\]\n",
    "$$\\label{convergence}$$\n",
    "The numerical solution $u_{\\rm num}(x0,x1,x2,t)$ should converge to the exact solution $u_{\\rm exact}(x0,x1,x2,t)$ at fourth order, which means that\n",
    "$$\n",
    "u_{\\rm num}(x0,x1,x2,t) = u_{\\rm exact}(x0,x1,x2,t) + \\mathcal{O}\\left((\\Delta x0)^4\\right)+ \\mathcal{O}\\left((\\Delta x1)^4\\right)+ \\mathcal{O}\\left((\\Delta x2)^4\\right)+ \\mathcal{O}\\left((\\Delta t)^4\\right).\n",
    "$$\n",
    "\n",
    "Thus the relative error $E_{\\rm rel}$ should satisfy:\n",
    "$$\n",
    "\\left|\\frac{u_{\\rm num}(x0,x1,x2,t) - u_{\\rm exact}(x0,x1,x2,t)}{u_{\\rm exact}(x0,x1,x2,t)}\\right| + \\mathcal{O}\\left((\\Delta x0)^4\\right)+ \\mathcal{O}\\left((\\Delta x1)^4\\right)+ \\mathcal{O}\\left((\\Delta x2)^4\\right)+ \\mathcal{O}\\left((\\Delta t)^4\\right).\n",
    "$$\n",
    "\n",
    "We confirm this convergence behavior by first solving the scalar wave equation at two resolutions: $16\\times 8\\times 16$ (or $16^3$ if `reference_metric::CoordSystem` is set to `Cartesian`), and $24\\times 12\\times 24$ (or $24^3$ if `reference_metric::CoordSystem` is set to `Cartesian`) and evaluating the maximum logarithmic relative error $\\log_{10} E_{\\rm rel,max}$ between numerical and exact solutions within a region $R < 0.1 {\\rm RMAX}$ at all iterations. \n",
    "\n",
    "Since we increase the resolution uniformly over all four coordinates $(x0,x1,x2,t)$, $E_{\\rm rel}$ should drop uniformly as $(\\Delta x0)^4$:\n",
    "$$\n",
    "E_{\\rm rel} \\propto (\\Delta x0)^4.\n",
    "$$\n",
    "\n",
    "So at the two resolutions, we should find that\n",
    "$$\n",
    "\\frac{E_{\\rm rel}(16\\times 8\\times 16)}{E_{\\rm rel}(24\\times 12\\times 24)} = \\frac{E_{\\rm rel}(16^3)}{E_{\\rm rel}(24^3)} \\approx \\left(\\frac{(\\Delta x0)_{16}}{(\\Delta x0)_{24}}\\right)^{4} = \\left(\\frac{24}{16}\\right)^4 \\approx 5.\n",
    "$$\n",
    "\n",
    "Since we're measuring logarithmic relative error, this should be\n",
    "$$\n",
    "\\log_{10}\\left(\\frac{E_{\\rm rel}(16\\times 8\\times 16)}{E_{\\rm rel}(24\\times 12\\times 24)}\\right) = \\log_{10}\\left(\\frac{E_{\\rm rel}(16^3)}{E_{\\rm rel}(24^3)}\\right) \\approx \\log_{10}(5).\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:32:36.327636Z",
     "iopub.status.busy": "2021-03-07T17:32:36.327008Z",
     "iopub.status.idle": "2021-03-07T17:32:36.949759Z",
     "shell.execute_reply": "2021-03-07T17:32:36.949251Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaEAAAEWCAYAAADPZygPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABGQElEQVR4nO3dd3wUdfrA8c+TXiH0DkGaFAGlKIqCqFhOxYKi59kb1vPUs3t2z9Pz9Dz9neDZ9exdFBUUsaHSBDxQiIKJSKgppCf7/P6YSbJpmyVkM7vJ83698sruzuzMM7Mz88x85zvfr6gqxhhjjBeivA7AGGNM22VJyBhjjGcsCRljjPGMJSFjjDGesSRkjDHGM5aEjDHGeMaSUDMSkQNF5Aev4wAQkVtF5LkWnqeKyMCWnKdpfiKyU0T2CHLckPzmIjJZRLKacXrNsm+KyHoRObQ5YmpuIrJARM5zX58mIh96HVMwdjkJuT9CkYjki0iOiHwpIjNFJGISWnNt4LV3QFX9TFWH7O50G5lnRxHZIiKf+33WrDtsuBCRHiLyuIj85m5va0TkNhFJ9jq2SCYiaSLyhIhsctfrjyJyXeVwVU1R1Z+aYT5xInK/iGS5iW29iDy4u9NtipbYNwFEZLyIvOceG7eLyDcicnao51ubqj6vqlObY1qhTrxNTRzHqGoq0A+4B7gWeLzZogoDIhLjdQwN+Buw2ssAmnvd1Dc9EekIfAUkAhPc7e0wIA0Y0Jzz3x0iEu11DE3wAJACDAXaA8cC60Iwn+uBscB4IBWYDCwNwXwCaql9WUQmAB8DnwIDgU7ARcCRzTwfiaST/kap6i79AeuBQ2t9Nh7wASPc9/HA34FfgGzgUSDRHTYZyAKuATYDvwHHAUcBPwLbgRv8ph0PPAhsdP8eBOJrTesqv2md7ffdo4D/AfnAr8DVQDJQ5Ma70/3rCdwKvAo8B+QB57nL9RWQ4077YSDOnfZCQIECdxozKuOpta6uBlYAucBLQILf8Gvc6W5056fAwADrfn83nrOBz93PAi3Py8Az7vJ/D4xt5Lc9H+dgtB14G+jpN0yBS4C1wM/uZ3/2i/8c//iD3AauBTYBz9YTy53ASiCqkfXxrbtuvwX29xu2ALgD+MJd/g+Bzu6w94FLa03rO+AE9/WewEfuevgBONlvvKeAfwPvub/9ocA+wDJ3Pq+4v/Odft85GljubkdfAiN3YRuZ5n43D8gAjnA/b49z4vcbzrZ9JxAd5D68CjguwHD/3/Ep4BFgjrt8XwMDao07090uctxxxR32LnBFI8eS63H20R3Ak5XLTuP79i5tX9TdN/sArwNbgG3Aw+7nA3ASyTZgK/A8kBbo+Oc37HPgkd3Yxxrbnu/C2Z6LcJLcYcAad/yHcZLfee74Z+EeI4L4nRpcZnfd+dx57gSucT/fD2dbzsHZdyb7zess4Cd3e/kZOC3gOglmo61nw6nzI7gbw0Xu6wfcFdwR5wzoHeCvfhtIOfAXINb9UbYA/3XHHe4ucH93/NuBRUBXoIu74HfUmtbt7rSOAgqBDu7w34AD3dcdgH38N9Ja8d8KlOEkxCicM/Ax7sqOAdJxrkCu8PtOjaRRe7ruuvoGJyl0dL8/0x12BM4OMhxIwkl+DSYhIBrnLHIMdTewhpan2F0n0cBfgUUBftcpOBvgPjg7+L+AhbWW9SN3ORLd+LOBETiJ8L/UPHgFsw38zZ1XYj3xLAJuCxBvR5wD1+nu73Oq+76T306bAQx2410A3OMOOwP4wm9aw3B2pnh3WTJxEn0MsLe7Xob5HZRzgQPc7aQdsAH4I842eAJQipuE3O9vBvZ1f4cz3e2i8kQq0DYy3p3XYe68egF7usPeAGa58XZ1p3GhO6yvuzx9G1h3/8E5KTkbGFTP8NpJaJsbSwzOAerFWuO+i3OF2hdnX65MlDfhHBcuBvbCPejV2j9W4SSEjjgH2Mr1NpnA+/YubV/47SPu7/CdO41kIAGY6A6rPLjH4xxvFgIPBnH8SwIqgIObso8R3Pb8C87xIsaNLR+Y7q6fP7nLHCgJNfQ77dIy42yH29zfJMr97jb3u8k4J0xD3HF7AMNbKgktAm4EBOcM0f9saQLVZ8+TcZJMtPs+1V1B+/qNvwT3TA3nQHKU37DDgfW1phXjN3wzsJ/7+hfgQqBdrVgnU/9Be2Ejy34F8EZ9O2t903XX1R/83t8LPOq+fgJ3p/HbEAIloT8B/25gA2toeebVOtAWBVi2x4F7/d6n4CTldL9lneI3/Ancg7r7fnBl/EFuA6X4nfHXE89a3INxA8NPB76p9dlXwFl+O+1NfsMuBub6bXMFQD/3/V3AE+7rGcBntaY7C7jFff0U8IzfsINwrkTE77PPqT6Y/hv3pMlv+A/ApCC2kVnAA/UsezegBL/kjXPQ+iTQ9us3biJwA85+VoZzZn5kfdu1u7z/8Rt2FLCm1rgT/d6/DFznvo7GuXr+wo13I3Bmrf1jZq1pZ/htI/Xu203ZvqiZhCbgHIRjAq0nd9zjgGW1Yq7v+NfLXRd7NmUfI7jt+Xa/YWfgd1LprpMsAiehen+nXV1mnCvMZ2t95wOcE6xknBOgE6nn5LK+v+YsV+yFc4nZBeesYIl7cy4HmOt+Xmmbqla4r4vc/9l+w4twfiBwzhA3+A3b4H7mP61yv/eFft89EWfD3iAin7pltoFk+r8RkcEi8q57AzcPuBvo3Mg0atvUQGw9a82vxrxrxdETuBwnye/OvBNEJMatObPT/XvfL56q9ayqO3HObno1EGPt+P1/o2C2gS2qWhwg9m04Z1ENqb1dVMbgH2+9615V83GKl05xh52Kc4YPzn3OfSvjdmM/DejuN63a6+FXdffEeob3A66qNb0+1NyGG9pG+uCchNXWD+fs9ze/ac7CuSJqlKoWqerdqjoG577Fy8Ar7n24+jQUX8Dhqlqhqo+o6gE4Z+B3AU+IyFC/8WtvQ8Hs27u7ffUBNtSaNgAi0k1EXhSRX919/jmC2+d34BRbBb3N1trHgtmeG9z/3O2vwWOIq97fqQnL3A84qdY2PRHooaoFOCdyM3G2zzkismegoJolCYnIOJyV9TnO5WYRziVYmvvXXlVrb7jB2oiz0JX6up81SlW/VdVpODvnmzg7GzhnBfV+pdb7f+OUuQ5S1XY4Z48SXNiN+g3o7fe+T4Bxx+Ns3P8TkU3AP4HxbnKMrifugNSpOZPi/lXeNK2xnt0aaJ1wzvKrvlorfv+Y+/q9DmYbaCzmecDxAW7A1t4uKmP4tZ5x6/MCcKp7YpIAfOJ+ngl86hd3mrueLmog9t+AXiLiv134r5dM4K5a00tS1ReCiDGT+ithZOJcWXT2m2Y7VR0exDRrUNXKk6tkoP+ufn8X5lOkqo/gHKyH+Q2qvQ0Fs2/v7vaVCfRtoMLC3e5393L3+T8QxD6vqoU4Vy4nBhgt0D4WzPbc4P7nbn+BjiGBNLbMtddlJs6VkP82nayq9wCo6geqehjOMWsN8Figme9WEhKRdiJyNPAi8JyqrlRVnzvTB0SkqzteLxE5vImzeQG4SUS6iEhnnHtJjT7/4lYPPU1E2qtqGU45pc8dnA10EpH2jUwm1f3eTjebX1RreDYQ1PMU9XgZOFtEhopIEnBzgHHfx7lkH+3+/QXnRvho94oy2OUJ5AU3ntEiEo+zYX6tqusDxH+WiAxz47+lckAzbQP/wLnf8rSI9PObxj9EZCROxYDBIvJ79+puBs7B7d0gp/8ezk5/O/CSGzPu9weLyOkiEuv+jat19u7vK5x7AZe6cUzDOWmo9BgwU0T2dWs1JYvI70QkNYgYH8f5TQ4RkSh3+fdU1d9wKlrc7+6DUSIyQEQmBbPgInKzu0xxIpKAcz8rB6eYsNmIyBXiPD6Q6K6bM3H2qWV+o10iIr3dq7AbcSpmBNQM29c3OAfxe9zfI0FEDnCHpeLcgM8VkV44lW+CdQ3OPvFnEenkxjVKRF50hwfax3Z1e54DDBeRE9xkejk1r9Z3RWPLXPs49xxwjIgcLiLR7vqb7P6O3URkmptgS9zp+gigqUnoHRHJx8mIN+IcMPzrwl+LU868yL28mwc0tY7+ncBinNpDK3Fuzt8Z5HdPB9a7MczEKVZBVdfgbBA/uZeTPRv4/tXA73FuAD5G3R3kVpyDZI6InBz0EjkxvA88hHMGvg7nnho4P1ztcUtUdVPlH87N6jL39a4sT6B45uEkwtdwdtABVBdXNRT/gzi1ata5//3t1jagqttxaguVAV+729t8nGVfp6rbcGqdXYVTpHENcLSqbg1y+iU4taMOxalUUfl5PjAVZ9k34hRhVN7grm86pTiVEc7FOZD/AefAUeIOX4xT+eZhnKuAdTjl9cHE+A3OfvWAu9yfUn22fAYQR3XNsldxi4JEpK84Ra1960zUnTROTbSt7jIeBvzOLR5qToXA/TjrcCvO/aETteYzSP/FSag/4RQ9BrtvN3n7ck/cjsG5f/kLzr2UGe7g23AqDuTiHOhfDzIeVPVLnMoHU3D2xe3AbJwEE3Af29Xt2f38JJxHZLYBg3DuvTVFY8v8V5wLgRwRuVpVM3Fqbd6Ac28tEydxRbl/V+JsV9uBSdQ9ea+hsoqe8Zh7pr0Kp9ZUnbJqEzlE5GucygVPeh1LOBOR9Tg30ud5HYvxTut54CkCicjxIhIvIh1wzrbfsQQUeURkkoh09ytyGolzo9wY0whLQt66EKfaaQbOfYWAl60mbA3Bee4kB6c4Zbp738YY0wgrjjPGGOMZuxIyxhjjmXBtpLNenTt31vT0dK/DMMaYiLJkyZKtqtql8TFbXkQlofT0dBYvXux1GMYYE1FEpHZrDGHDiuOMMcZ4xpKQMcYYz1gSMsYY4xlLQsYYYzxjScgYY8KUz+fjscce4w9/+APl5a2zMRVPkpCI3Ccia0RkhYi8ISJpXsRhjDHh6u2336ZTp05ccMEFvPDCCyxbtqzxL0Ugr66EPgJGqOpI4EecfuaNMaZN+/HHH3n99dc5//zzOeGEE8jNzeWQQw4hIyODcePGeR1eSHjynJCqfuj3dhFOP+nGGNPmlJaW8s4773DnnXeyfPlyAOLi4pg5cybXXnstffo0ta+6yBAOD6ueQxAdWRljTGuiqlx55ZU8+uijFBc7PZGnpqZy8cUX86c//Ylu3bp5HGHLCFkSEpF51N/T342q+pY7zo1AOfB8gOlcAFwA0LdvQ310GWNMeKuoqGDhwoXMnj2bbt268cEHH7BmzRpEhPHjx3PllVdywgknEBsb63WoLcqzVrRF5CycrgwOcftnb9TYsWPVmu0xxkSKnJwcZs2axRtvvMGyZcsoLS0FICoqikmTJjFjxgymT59Op06dQhqHiCxR1bEhnUkTeVIcJyJH4HRfOynYBGSMMeFMVVm2bBmvvfYa2dnZ/PLLL3z55ZcUFBQAkJiYyJQpUzj77LM59thjadeunccRhwev7gk9DMQDH4kIwCJVnelRLMYYs0vKy8v54YcfWLZsGbNnz2bNmjVs27YNn89XNc7IkSM5/fTT6dWrF6eccgoDBgzAPd4ZP17VjhvoxXyNMWZXVF7dfPDBB3z99desXr2ajRs3UlhYWJVwRIS4uDj69+/PqFGjOPTQQznppJPo3Lmzx9FHhnCoHWeMMZ4qLS1lxYoVzJs3j0WLFpGZmUn79u1ZuXIlW7durRpPRGjXrh1jx47lsssuY/To0aSnp5OSkuJh9JHNkpAxps3Iz8/n+++/Z+HChfTt25dVq1bx0ksvsW7duhrjiQhjx45l2rRpJCcnk56ezuGHH87QoUOtSK2ZWRIyxrQ627dvJyUlhezsbJ544gn++9//kpWVRWFhzXpQ0dHRdO/enaFDhzJs2DAmTJjAYYcdxogRI4iKsqY1W4IlIWNMRMvKyuKFF17gq6++4vvvvyczM5OioiJSUlLYuXNn1XiJiYkMGjSoKtkceuihDB8+nISEBA+jN5aEjDFhr6ioiO+//57FixezcOFCli1bxogRIygoKGDx4sVs2bKlatzY2FjS09OZOHEiEyZMYMSIEey111506NDBwyUwDbEkZIwJG6rKr7/+yrJly/D5fKgqCxcu5MEHH6T2g/Vr165l2LBhTJkyhT322IP999+fkSNH0qdPH7tvE0EsCRljPFGZVDZt2sSll17Kd999R1ZWFiUlJXXGTUtLY+DAgYwdO5aJEycyatQohgwZ0uaauGmNLAkZY0IuOzubBQsW8Mknn7B06VIyMjKqrla2bdtWNV5CQgKDBw9m5MiRHHTQQYwfP57hw4dbFehWzJKQMabZlJaWsmDBAubNm8ePP/5I165dWb16NYsWLarRM2h0dDS9evVi6tSpjBgxouq+TdeuXT2M3njBkpAxZpdt376dpUuXkp2dTUZGBu+++y6rVq2iqKioxnidOnVi6NChHH744fTr14/Jkyczbtw4+vbta1WgDWBJyBhTD5/Px6ZNm8jIyCAjI4Nvv/2Wr7/+mszMTLZv317jqgagY8eOJCUlMWzYMEaNGsVBBx3E1KlT6dGjh0dLYCKFJSFj2hhVJTc3l6ysLH755RcyMzP56aef+P777/n555/57bffyM3NrdEYp4igqkRHR9OpUyf69evH8OHDOf/889l7771JTEz0cIlMJLMkZEyEq6ioIC8vjx07dpCTk8OOHTvYunUr2dnZVX8bN25k/fr1bNmyhZycnDpXMlFRUTWSTlJSEt26deOyyy7jmGOOoX379hQVFVn1Z9PsLAkZ00JUlZKSEoqKiiguLqagoIDCwkIKCgrqfb1z504KCwurXm/durUq2eTm5pKXl0dRUREFBQV1nqGpJCJ0796dtLQ0Vq9eXfV55RXNxRdfzHnnnUdcXBwLFixg4MCBDBgwwPq6MS3GkpBpkyoqKti5cyfFxcV069at6qHI9evXk5ubS05ODnl5ecTExDBmzBjy8/OZM2cOmzZtoqioiLKyMsrLy0lKSmLw4MEUFRWxdOlSCgoKqKioqPqLjo4mLi6OoqKiOjftg5WSkkJiYmKNVgEqjR49mmOPPZbk5GRuuOEGUlNTSUtLo3PnznTr1o2zzjqL6dOnU1JSwocffkjv3r3p3bs3nTt3rnNFc9JJJzUpPmN2h2fdezeFde9tAlm9ejWrV6/m559/JjMzk40bN1JQUMC0adPYsWMHL774Ij/++COlpaVVxVExMTEkJyezc+dOKioqgp6XiCAixMfH07t3bxISEsjMzKS8vJyoqCiio6OrGsecNGkSiYmJfPjhh5SWlhIbG0tcXBxxcXEMGzaM6dOnk5SUxOOPP46IkJCQQFJSEgkJCey7775MmzYNgKeeeoqEhARSUlJITk4mJSWFXr160bNnT8C50rKiMlOfcO7e25KQCVuqyubNm+natSsFBQW88sorvP/++2RmZpKdnc2OHTsoKCjgwAMPZPPmzaxbt47i4uIGp1d5VZKYmEhycjKpqal07NiRvffem9TUVIqKioiLi6N9+/akpaVVXVH06NGDlJQUUlNTSU5OJj4+3g72JqKEcxKy4jjjqZKSEqKjo9m8eTNvv/02b775Jhs2bGDz5s3k5uZSUVFB+/btyc3NrfPd6OjoqpaSBw4cyIgRI0hNTaVfv37sscce9OrViw4dOtCxY0c6dOhAYmKiJQ9jwowlIdMiKioq+Pbbb3n11VdZsWIFGRkZZGdnU1BQQHx8fJ32wmJjY+nQoQPdunVj3333ZdCgQfTo0YP+/fvTq1cvunfvTnJyskdLY4xpLpaETLMqLi5m3rx5zJ8/n8WLF7Nu3ToSExPZuHFjjUQjIqSkpDB48GAmTpzIPvvsQ3p6OnvssQd9+/a1BGNMG2FJyDRZVlYWr7/+Oj///DMbN25k6dKldbpJjouLY+TIkUyfPp309HTat2/PxIkT6dOnjzXbYoyxJGSCs3nzZr766iuefvppli9fTlZWFmVlZVXD09PT2Wuvvapu9B988MFMnTrVOhIzxgRkScjUUXn/5sUXX2TBggVV7YVVEhG6dOnCsGHDmDx5MtOnT2f48OEeRmyMiVSWhAwVFRUsX76cjz/+mCeffJIffvihRhMuXbt25d5772XChAl07dqVPfbYg5gY23SMMbvPjiRt1LJly3j00UeZP38+69evr3pQs1OnTnTt2pVx48Zx/PHHc/zxx5OWluZtsMaYVsuSUBtRVlbG559/zv/93/8xZ86cqiZkKtsWu+GGGzjxxBOt6X1jTIvyJAmJyB3ANMAHbAbOUtWNXsTSmq1Zs4b77ruPDz74gJycHAoKCoiJiSEpKYkDDzyQM844gxNPPJGEhASvQzXGtFFeXQndp6o3A4jI5cBfgJkexdKqbNy4kcsvv5z58+eTk5MDOM30T548mcsuu4xDDz2UlJQUb4M0xhiXJ0lIVfP83iYDkdOAXZhRVebPn897773HunXrmDt3LmVlZcTFxTFx4kQuuugiZsyYQXR0tNehGmNMHZ7dExKRu4AzgFzg4ADjXQBcANC3b9+WCS7MqSrffvstf/3rX/nwww8pLCwEoHfv3lx66aUcffTRTJ482R4GNcaEvZC1oi0i84Du9Qy6UVXf8hvveiBBVW9pbJptvRVtVeXLL79k5syZrFq1CnBaJJg8eTJXX301hxxyiCUeY0wdbbIVbVU9NMhRnwfeAxpNQm2Rz+fj6aef5r777iM/P5+srCzi4+MZN24cV111FSeeeKI9s2OMiVhe1Y4bpKpr3bfTgDVexBHOMjMzueqqq3j77berGv4cPXo0d911F8cffzypqakeR2iMMbuv0SQkImOBA4GeQBGwCvhIVXfsxnzvEZEhOFW0N2A14wDnqmfBggXMmjWLl19+GYCkpCROO+007rjjDvr37+9xhMYY07waTEIicjZwGfAzsAT4AUgAJgLXisgq4GZV/WVXZ6qqJzYt3NZp48aNPPzwwzz22GNs3bqVtLQ0Jk2axFlnncWZZ55pHbEZY1qtQFdCScABqlpU30ARGQ0MAnY5CRmnksFXX33F3/72N9555x1UFRHh/vvv56KLLiIxMdHrEI0xJuQaTEKq+oiIRIvIn1T1gXqGLw9pZK1UeXk5L7/8Mg888AD+Nf2OPfZY7rvvPgYPHuxhdMYY07IC1udV1Qrg9y0US6tWXl7OU089xZ577slpp51GTk4OycnJHHzwwaxcuZK33nrLEpAxps0Jpnbc5yLyMPASUFD5oaouDVlUrYiq8tJLL3HTTTeRkZFBjx49eO211zjuuOPIysqiT58+ds/HGNNmBZOERrv/b/f7TIEpzR5NK7N06VIuv/xyvvjiC/r370/Hjh3ZvHkzXbt2JSoqylqAMMa0eY0mIVVtsEkdU7/CwkJuuOEGHnroITp16sTkyZNZsGABe+21Fx9++CFjxozxOkRjjAkLwTwn1B6nNYOD3I8+BW5X1dxQBhapvvrqK84880zWrl3LJZdcQmZmJm+//TZXXXUVd999N3FxcV6HaIwxYSOY4rgncB5QPdl9fzrwJHBCqIKKRKrKAw88wDXXXEOvXr2YP38+U6ZMYfHixZx33nkcc8wxXodojDFhJ5gkNKDWw6W3icjyEMUTkQoLCznnnHN46aWXOP7449l///2ZM2cOU6ZMYezYsGwz0BhjwkIwTS4XicjEyjcicgBO8z0G2L59O4cccgivvPIKd999N3379uXPf/4z69ato6yszOvwjDEmrAVzJTQTeMa9NwSwAzgzdCFFjqysLA4//HAyMjJ48cUXmTt3Lk888QR//OMfuf/++60jOWOMaUTAJCQi0cDpqjpKRNpBnV5R26xNmzZx8MEHk52dzdy5c3nmmWd48skn+ctf/sKtt95qz/4YY0wQAiYhVa2oLIqz5FNt+/btTJ06ld9++42PPvqICRMmkJuby5AhQ7j22mu9Ds8YYyJGMMVxy0TkbeAVaraY8HrIogpjpaWlHHfccfzwww+89957tG/vlFJOmzaNadOmeRydMcZElmAqJiQA23BaSDjG/Ts6lEGFs8svv5zPPvuMp556iuXLlzNy5Ei+/PJLr8MyxpiIFMw9oW2qenULxRPWZs2axaxZs7j22mtJTk7m6quv5qSTTmLffff1OjRjjIlIwdwTOqClggln//vf/7jiiis4/PDDOf3005kwYQJjx47lmWeesVpwxhjTRMHcE1re1u8JlZaWctppp5GSksIjjzzC0UcfTUJCAq+//joJCQleh2eMMRErmCTkf0+okgJtJgndddddLF++nDfffJP09HTOPfdcxowZQ58+fbwOzRhjIpqoqtcxBG3s2LHq3xtpS1i7di0jRozgxBNP5LnnniMqKpi6HMYYEz5EZImqhmUbYo0eUUVksIjMF5FV7vuRInJT6EPznqpy2WWXkZCQwHXXXceIESP46KOPvA7LGGNajWBO6x8DrgfKAFR1BXBKKIMKF3PmzOGDDz7g9ttv59577yUjI4OePXt6HZYxxrQawSShJFX9ptZn5aEIJpz4fD5uvPFGBg4cyIgRI3j++ee55pprGD58uNehGWNMqxFMxYStIjIApzICIjId+C2kUYWBV155hRUrVvD000/zxz/+kfT0dK6//nqvwzLGmFYlmCR0CTAb2FNEfgV+Bk4LaVQeq6io4JZbbmH48OHExMTw/fff89Zbb5GUlOR1aMYY06o0moRU9SfgUBFJBqJUNb+5Zi4iVwF/B7qo6tbmmu7uevvtt/nhhx94+eWXmT59OgMHDmT8+PFeh2WMMa1OMFdCAKhqQeNjBU9E+gBTgV+ac7rN4YEHHiA9PZ2jjjoKEbEEZIwxIeLlQy8PANfg3msKF4sXL+azzz7j3HPPZcCAATz77LNeh2SMMa2WJ0lIRKYBv6rqd0GMe4GILBaRxVu2bAl5bA8//DApKSkUFxeTnZ3N6NGjQz5PY4xpqxptMUFEkoCrgL6qer6IDAKGqOq7jXxvHtC9nkE3AjcAU1U1V0TWA2ODuScU6hYT8vPz6d69OzNmzODdd99l3LhxzJkzJ2TzM8aYlhDOLSYEc0/oSWAJMMF9/ytOY6YBk5CqHlrf5yKyF9Af+M7tArs3sFRExqvqpiDjDolXX32VwsJCunbtypYtW7j6auvBwhhjQimYK6HFqjpWRJap6t7uZ9+p6qhmCSCMroQOOuggNm3aRHR0NElJSSxevBg3URpjTMSK9CuhUhFJpPph1QFASUij8sD69ev57LPPuOuuu5g2bRoFBQWWgIwxJsSCSUK3AnOBPiLyPHAAcFZzBaCq6c01rd3xxhtvADBjxgwGDBjgcTTGGNM2NFo7TlU/BE7ASTwv4BSdLQhtWC3v9ddfZ/jw4dx2222sWrXK63CMMaZNCKYrh3dwHipdoKrvhlPLBs0lOzubL774gv79+/Pss89SUNCsz+UaY4xpQDDPCf0dOBD4n4i8KiLTRaRV9Wn97rvvoqps3LiRoUOHWgsJxhjTQoIpjvtUVS8G9gBmAScDm0MdWEv68MMP6datG8uWLWPGjBlWIcEYY1pIUC0muLXjTgRmAuOAp0MZVEvy+XzMnz+f9PR0VJUTTzzR65CMMabNaLR2nIi8DIzHqSH3MPCpqvpCHVhLWb58Odu2bePYY4+lXbt21mmdMca0oGCqaD8OnKqqFaEOxgvz5s0D4K677qJHjx4eR2OMMW1Lg0lIRKao6sdAMjCt9n0SVX09xLG1iE8++YQhQ4bQtWtXr0Mxxpg2J9CV0CTgY+CYeoYpEPFJyOfz8fXXX9O1a1cGDRpERkaGVUowxpgW1GASUtVb3Je3q+rP/sNEpH9Io2oha9euZceOHYgIkyZNsgRkjDEtLJjaca/V89mrzR2IFxYtWgTA9u3bOeKIIzyOxhhj2p5A94T2BIYD7UXkBL9B7YBW8bDqokWLSEhIoLi4mClTpngdjjHGtDmB7gkNAY4G0qh5XygfOD+EMbWYRYsWkZaWhqpao6XGGOOBQPeE3gLeEpEJqvpVC8bUIoqLi1m5ciWnnHIKRx55pN0PMsYYDwTznNAyEbkEp2iuqhhOVc8JWVQtYPXq1VRUVHDsscdy8sknex2OMca0ScFUTHgW6A4cDnyK0x13fiiDagkrV64EoLCwkPLyco+jMcaYtimYJDRQVW8GClT1aeB3wL6hDSv0Vq5cSXR0NBdddBE+X6tphcgYYyJKMMVxZe7/HBEZAWwCIr55gRUrVpCYmMjQoUOJi4vzOhxjjGmTgklCs0WkA3Az8DaQAvwlpFG1gBUrVlBWVsaoUaO8DsUYY9qsRpOQqv7HffkpTp9CES8vL49NmzYBMHr0aG+DMcaYNizQw6pXBvqiqv6j+cNpGRkZGVWv7UrIGGO8E+hKKLXFomhh69atA+C5555jzJgxHkdjjDFtV6CHVW9ryUBaUuWV0LHHHktiYqLH0RhjTNvVaBVtERksIvNFZJX7fqSI3BT60EJn3bp1pKSksHTpUq9DMcaYNi2Y54QeA67HraqtqiuAU0IZVKitXbuWgoIC3n//fa9DMcaYNi2YJJSkqt/U+my3mhgQkVtF5FcRWe7+HbU709tVP/74I6rK4MGDW3K2xhhjagnmOaGtIjIApzdVRGQ68FszzPsBVf17M0xnl5SUlFRVzx40aFBLz94YY4yfYJLQJcBsYE8R+RX4GTgtpFGF0MaNG6te25WQMcZ4K2BxnIhEAxer6qFAF2BPVZ2oqhuaYd6XisgKEXnCbZGhoRguEJHFIrJ4y5Ytuz3TrKwsAJKSkujaNeJbHzLGmIgWMAmpagUw0X1doKpBt54tIvNEZFU9f9OAfwMDgNE4RXv3B4hhtqqOVdWxXbp0CXb2DapMQvPmzbM+hIwxxmPB9if0NvAKUFD5oaq+HuhL7tVTo0TkMeDdYMZtDpVJaPjw4S01S2OMMQ0IpnZcArANmILTzfcxON1+N5mI9PB7ezywanemtysyMzOJjY3l888/b6lZGmOMaUAwDZieHYL53isio3Fq3K0HLgzBPOq1fv16ysrKWLNmDUcd1aI1w40xxtQSTHFcs1PV072YLzhJCKBv375ehWCMMcYVTHFcq7J582bAkpAxxoSDNpeEcnNzAejTp4/HkRhjjAmmAdNnRaS93/t+IjI/tGGFRklJCcXFxSQmJtKtWzevwzHGmDYvmCuhz4GvReQoETkf+Ah4MKRRhcjWrVsBeOCBB4iKanMXgcYYE3aCqR03S0S+Bz4BtgJ7q+qmkEcWApUtLnTu3NnjSIwxxkBwxXGnA08AZwBPAe+JSET2iV2ZhObPj8jSRGOMaXWCqaJ9IjBRVTcDL4jIG8DTOE3uRJTKJFTZirYxxhhvBVMcd1yt99+IyPiQRRRClUmod+/eHkdijDEGgkhCIpIAnAsMx2nCp9I5oQoqVCrbjevXr5/HkRhjjIHgasc9C3QHDgc+BXoDQbemHU5++83pi8+qZxtjTHgIJgkNVNWbgQJVfRr4HbBvaMMKjby8PGJiYqw4zhhjwkQwSajM/Z8jIiOA9kBE9gYXHR3NsGHDmDx5stehGGOMIbjacbPdnk9vBt4GUoC/hDSqEMnLy6Ndu3Zeh2GMMcYVTO24/7gvPwX2CG04oZWRkUF+fj7l5eXExHjSgLgxxhg/DR6JReTKQF9U1X80fzihlZeXR25uriUgY4wJE4GOxn8HlgPvAyWAtERAoVRcXExcXJzXYRhjjHEFSkJ7A6fi1IZbArwAzFdVbYnAQqG0tJTU1FSvwzDGGONqsHacqn6nqtep6mjgcWAa8D8RObalgmtOZWVlVFRUkJSU5HUoxhhjXME0YNoF56poLyAL2BzqoEIhP995vtZaSzDGmPARqGLCOcDJOE31vAqc7DZiGpHy8vIAuOCCCzyOxBhjTKVA94T+A6wCNuA02TNVpLpugqpGVLFcZbfe9pyQMcaEj0BJ6OAWi6IFVF4JzZ07lxNOOMHjaIwxxkCAJKSqn7ZkIKFWeSVUXl7ucSTGGGMqNVgxQUTeEZFjRCS2nmF7iMjt7n2jiLBjxw4AOnTo4HEkxhhjKgUqjjsfuBJ4UES2A1twKimkAxnAw6r6VsgjbCbbt28HLAkZY0w4CVQctwm4BrhGRNKBHkAR8KOqFu7ujEXkMuASoAKYo6rX7O40A8nJyQGgY8eOoZyNMcaYXRBMz6rdgI44Tff81kwJ6GCch19HqWqJiIS8a4iyMqdHigEDBoR6VsYYY4IU6Dmh0cCjOP0H/ep+3FtEcoCLVXXpbsz3IuAeVS0BaInnj1JSUgA46KCDQj0rY4wxQQp0JfQUcKGqfu3/oYjsBzwJjNqN+Q4GDhSRu4Bi4GpV/XY3pteo4uJiABISEkI5G2OMMbsgUBJKrp2AAFR1kYgkNzZhEZkHdK9n0I3ufDsC+wHjgJdFZI/6GkcVkQuACwD69u3b2Gwb9N133wGwfv16+vfv3+TpGGOMaT6BktD7IjIHeAbIdD/rA5wBzG1swqp6aEPDROQi4HU36XwjIj6gM04NvNrTmQ3MBhg7dmyTW/CufFg1Ojq6qZMwxhjTzALVjrtcRI7EqUDQy/34V+ARVX1vN+f7Jk6LDJ+IyGAgDti6m9MMqLI4LjExMZSzMcYYswsC1o5T1fdxOrVrbk8AT4jIKqAUODPU/RQVFRUBloSMMSacNKmfaxGZrapNbo5aVUuBPzT1+01RUlICWBIyxphwEqiKdkNPdQpwVGjCCZ34+Hg6dOhg94SMMSaMBLoS2oLTjYP4fabu+5A/XNrcOnToYL2qGmNMmAmUhH4CDlHVX2oPEJHMesYPa0VFRSQnN1qz3BhjTAsKlIQeBDoAdZIQcG9IogmhH3/8sarpHmOMMeEhUBXtRwIM+1dowgkdS0DGGBN+gmnAtL5uSHOBlS3R5ltzqaioIC4uzuswjDHG+Ammiva5wATgE/f9ZGAJ0F9EblfVZ0MUW7Py+XzExDSpRroxxpgQCeaoHAMMVdVsqOra4RlgX2AhEDFJKCqqwY5kjTHGeCCYo3KfygTk2ux+th2ImBst8fHxu9UAqjHGmOYXzJXQAhF5F3jFfT/d/SwZyAlVYM0tKSmJvffe2+swjDHG+AkmCV0CnABMdN8/DbzmtvV2cKgCa25lZWV2T8gYY8JMo0dlVVUR+RynoVEFvgl1Y6OhsGPHDj766COvwzDGGOOn0XtCInIy8A1OMdzJwNciMj3UgTU3VbV244wxJswEUz51IzCu8pkgEekCzANeDWVgzU1VrTjOGGPCTDC146JqPZS6LcjvhR1LQsYYE16COSrPFZEPgBfc9zOA3e1ZtUVV3sKKjY31OBJjjDH+gqmY8GcRORE4wP1otqq+EdqwmldFRQUAAwYM8DgSY4wx/oIqn1LV14DXQhxLyJSXlwMwcuRIjyMJXz6fj4zVK8grKgOJyNJWY9omVeKknF7pg+jYqZPX0eyyQD2r5uNUya4zCKfmdruQRdXMKpOQiDQyZtv1W9YGJOtb9v76BqJ8pV6HY4wJkk9iKEwbwrryeygp7k+PXn28DmmXNHjKq6qpqtqunr/USEpAAKWlzkF1wYIF3gYSxrZuyabPd/dbAjImwkRpOSk7vmfgouvIzFjDxvVrvQ5pl7SJcpeioiLAKiYEUq4xxBVlNz6iMSYsJeX8QFRSGm8++QAlRYVehxO0NpWErD+hAEQQ9XkdhTGmiaK0HJEoysvKyNux1etwgtamkpBdCRljWjsRweeLnBPKNpGEiouLAbsSimTpD+Yz76dyr8OoQW7LY9326p39g3XlHPdi5BSDAPzr61Ku/ajY6zBMG9YmklBlSwkDBw70OBLjpbPeLCLujjxS7q7+G/Xozmab/o0fF3PdxOoTHbktj73+vROfX3u/N31czFlvFlW9X76pgjGzd5J0Vx5jZu9k+aaKXZ7vfV+UMOL/dpL61zz6/zOf+74oqXe8T9eXI7flcdPH1Unn/DGxPL+yjM0FkXPmbFqXNtGOTWJiIgCDBw/2OBLjtWsOiOPOKQmNjlfuU2Kigq/S/+2vFeSWwH69a+5SG/OVF1eV8/u96hYFl1Yo014s5Ip947h4XByzlpQy7cVC1l6WQlx08PNW4JnjExnZLYqM7T6mPldIn/ZRnDKiep5lFcof5xazb6+ajfgmxAhHDozhme/KuHr/+KDnaUxz8eRKSEReEpHl7t96EVkeyvlV3hOqbDnBRK6ScuWKucX0vD+fnvfnc8XcYkrKq6807v2ihB7usP8sLa1TZNaQ9Tk+5LY8Hl9aSt8H8pnytFOs9sSyUoY+spMOf8vj8OcK2JBT/7TeX1fOpH51W2m/Zv84bllQQrmv7iN3C9ZXUO6DK/aLIz5GuHzfeFTh4593bTu95oB49ukRTUyUMKRzNNOGxPLFLzWLLu//qpSpA2LYs3PdXX5yegxz1oZXUadpOzy5ElLVGZWvReR+IDeU81u/fj0AS5cu5bTTTgvlrFqNK+YWN6loaFeM7h7Ng0c0flXi767PSliUVcHymckIMO3FIu5cWMIdUxKYu66cf3xVyvwzkujfIYoL3tn1ex2fbqhg9SUpRAm8taaMuz8r4Z1TkxjUKYp7Pi/l1NeK+PLc5DrfW7m5gvE96yahE4bG8PL/ynhqeRnn7VPznuT3mysY2S2qxkPUI7tF8/3mCo4YGMM9n5dwz+f1F60B5FxX93E9VeWzX8q5cEz1vDbk+HhiWRlLL0zm0vfqrpOhXaL4LsS/tTEN8bQ4Tpy972RgSijnU3kFZP0JRb7nV5bxryMT6JrsnNHfMimeC98t4o4pCbz8fRlnj45leFfnd751cjzPryyr8f2/f1nKw99UP5A7bc9Ynj4user9rZPjSY5zksKjS0q5fmI8Q7s407vhwDju/qyEDTk++qXVvKLIKVZS4+sWoYkIdxwcz0VzijljVM0iuZ2l0L7Wd9onQH6pc9V03cR4rpu4a0Vkty4owadw9ujqeV0+t5g7Do4nJa7+Ir7UOMhtONcZE1Je3xM6EMhW1QYf8RWRC4ALAPr27dukmVQ22xMV1SbqYTSLXb1CaSkb85V+7at/x35pwsZ8dYf5GNuz+uDbp13dg+7V+we+J+T/nQ05zn2Uqz6svnpQ4Nf8ukmoQ4KQX1J/h8NHDYqld7tSZi2umRBT4iCv1sE/rwRSG0gWjXn4m1KeWVHGZ2cnEx/jTOOdH8rIL1FmjGj48YT8Umhvt4OMR0KWhERkHtC9nkE3qupb7utTqe4iol6qOhuYDTB27NgmdStemYSsP6HI1zNV2JDrq7ra+SVX6ZnqHHB7pEaRlVd9zyYzb9c3F//mBfu0F248MJHTRjb+fNnIbtH8uK3he093TYnn1NeKONUvGQzvGs39X5WiqlVFciuyK7hknFOUdvdnJdz9WcOXKDtvqC6Oe2JZKfd8XsLCs5Pp3a46Qc7/uYLFGyvo/vd8AHJLlGiBlZt9vHVKEgCrt/gY1d1KCYw3QnZUVtVDAw0XkRjgBGBMqGKoVFkcZ1dCke/UEbHcubCUcT2jEYHbPy3hD26SOHlYDOe8XczpI2PplxbFHQt3r4xp5pg4bv6khNHdoxjeNZrcYuXDjHJOGl43KR01KIZTXi2qZyqOyekxjOgaxdPflXHM4Bj3s2iio+Chr0uZOTaOx5Y6V0pT+lcW/8Vzw4GNX6I8v6KMG+aX8MmZSezRoeY2fsfB8TWqjf9xbjE9U6K4eVL1dD/dUM6RA+0EzXjDyy3vUGCNqmaFekapqakADBkyJNSzMiF200Hx5JWUMPLRAgBOGhbLTQc5B9QjB8Vy+XgfBz9dSJTAzQfF88x3ZcT7neTf+0UpDy6qvieUECNsvSa13nkdPzSWnaXKKa8VsSHHR/sE4bA9YupNQvv0iKZ9AnydVc6+vevfre48OJ79Hq9+mDUuWnhzRhLnvVPEdfNLGNo5ijdnJO1S9WyAmz4pZluRMu6xgqrP/jAylkePTiQ1Xmrcq0qMEZLjoGOi81lxufLe2nKWXFC3soUxLUFUm1TCtfszFnkKWKSqjwb7nbFjx+rixYt3eV5z587lyCOP5KuvvmK//fbb5e+3BUuWLGHMOyGtH9LiVm+pYMS/Cyi5KXWXnvlpqg8zyvm/b0t50y3migT/+rqUzDwf9x4WnvcAza5ZcszHfPn6LE6+6Aa69U6v+lxElqjqWO8ia5hnV0KqelZLzavyOaGSEqsC1Nq9sbqMowbFUFgG184r4ZjBMS2SgACmDohh6oDIKta6bF9rysp4q03cJMnMzATgm2++8TgSE2qzlpTS9e/5DHgon+go+Pfv7AzfmHAWWadtTWQVE9qOuX+wexvGRBI7KhtjjPGMJSFjjDGesSRkjDHGM20iCaWlpQEwbNgwbwMxxhhTQ5tIQikpKUDT254zxhgTGm0iCZWWOk/I5+aGtMcIEyLh2LU3tI7uvd/5oYwZr4Ym5mD7cgrGrMWlXDG3ZjcUBzxRwLLfwqsLiuydPoY+srNGH1f+lmys4LYFJWzM37X1UlKu7PnwTra0wh5w20QS2rRpE2DPCbV1rbV7709+Lufgpwtof08e6Q/m1zvOPxeV0P+f+STfncfQR3by4zZnPscMieX7zT5WZIfXwdxfaYVy52cl/Hn/6nX7zg9lpMbB3j2cNplWba7g8OcK6HxvPnJbXr3TeXFVGUMf2Uny3XkMeCifzzbUPLH562cl3DC/mEVZ5Rz2bAEd/5ZHl/vyOemVQn6rJ2mUVihDH9lJ739Ur/NuKVEcnB7N7CVldcZfs7WCI54vZP7P5Rz+XCE5xQ23VnPIMwXIbXlVnSHGxwjn7B3LPZ+XNvidSNUmkpBXTROZ8HPNAXHsvKFd1d93M1PqHa++nlADaax77/pUdu/9h71i2XFtKmeOimXai4WUVuzavJPjhHNGx3JfA03v/GdpKY8vK2PO75PYeX0q756aROek6l3/1BGxzF4S3MHt1gXF3Lpg1zsL3B1vrSlnz85R9PJrHfzRJWWcPrI6KcVGwcnDYnn82PrXwUcZ5Vw7r5gnpyWQf30qC89KrtPY65y15Rw1KIYdRcoF+8Sx/opUNlyRQmqccPZbdRunve+LUrok1W2N47S9YplVa31m5fk48vlC/npIPJ+elcTB6TEc+0IhxfVcMT2/ooyyes4Jfr9XLE9/V9bgVVaksiRkIop1713X+F7RnD4qrs5BFcCnym2flvDA4QkM6xKNiDCgY1RVA6bgtOYd6u69c4uVM94oost9+fR7MJ87F5ZUXSFW+JSrPiim87359P9nPg9/U1rjKsBZt9XJvbRC+fjncialV6/vIZ2jOXefuKouPmq7ZUEJfzkonv16xxAlQq92NZPajiLlx20+JvSO5shBsZw0PJZ28UJSrHDp+Di+yKz5m/y8w8dzK8u4fmLdZo/27R3NTzt8VdvK9iLld/8t5I6D4zlvnzhEhIeOTGD/PtGc8moRFX7bRm6x83vde1jd1tN7t4uiQ6KwKCt8r1qbok20mFDJvxtl07jJTxXU+ezk4bFcPC6OwjLlqOfr3ks4a3QsZ42OY2uhj+kv1z17vGhsHDNGxJKZ66NP+10/B7LuvavV1713bVl5SlaesmpzBWe9VURMFJwxMpZbJscT5c53aJdo1ucoeSVKu3p6h20Ol71fTG6J8tPlKWwrUqY+W0iPFOHcfZwuLN5fV87ymckkxwonvVJzu1q5uaJGVxNrt/mIEmr0mxRIhU9ZvLGCY4fEMPChfIrL4bg9Y7jvsAQSY53l/SCjnEP2iCa6nnYGF24oZ3iXmr/tZe8Xc/eU+Krv+4uJEgZ2jOK77Ar6pTkJv74r7nsOrXvVdsP8Yi4aG0v3lPqXbWhnZ7qT0lvPobtNXAmZ1uP5lWX8ZVIcXZOj6JIcxS2T4nl2hVP+7t+9d1KscOvkumeTf/+ylLR78qr+znyzZqKs7N47MVZqdO8dEyXccGAcyzdV1Hs11Fj33ncsLKlTzBZM994517Vr8C8YlZ38ffhTBSsvSuGTM5N5YVU5jy+tvmeRGle9DKFQ4VNeXFXGXw+JJzVeSE+L4qoJcTV+tz/uG1d1pl+7S/OcYmqs25xi3aXeZ7MLlDIfvPo/p9fZ5TOTWbbJx51+/U3NWVvOUfX0qbQiu4LbPy3hPr8rkzdWl1GhyvFDG+7sMDVednl9Lt5YwReZFQEblU2Nd9ZHa9J60mkAnTp1AmDEiBEeRxJZFpzVcDtsSbEScHjnpKiAw5tyFQTWvfeuSnS7+b5m/zjSEoS0BOHCMbG8t66c88c4B7t89/ZFWkL98z36v4V8/otTXFfsltpV9sk0sW8M7/4+cNcVWwudJFDzd4vi16rfTenTvnretX+3DgnUWLcdEqUqUQejch1cNj6OHqlODFfuF8edn5Vw1yFOkeVHGeX8Y2rN5Lduu3Mf559HJHCgWxxYUKpcM6+E936fGHCe+SXa4Pqsj0+Vi+cU8c8jEtxW3+tfvvwSSGtlbfK2iSSUnOwcDHv27OlxJGZ3Wffe1fy7927IkM5RxEXXXK7apdKrt1SQniYNFsX5J5nKSgm3Tg7+SNg5SYiNgg25PoZ1qfzdfPSq+t2ELL/fqvbvVnvdDuwYhSr8muercV+nIR0Shd7tpMF18O2vFfRLE7okV09rQ46PQ58p4OaD4jl9VPWVydrtPtbn+DjwSafIsLRCyS2B7n/PZ9F5yaSnRVHuU9Zt9zGqW/BdpueVwOKNPma4vfNWXjT3/sdOXjkpsSoJrt7q46oJrav7jTaRhIqLnR1n+/btHkdidpd1712XT5XSCiirUBSnt9QocXpuTYoVZgyP5d4vStm7ezS5JcrsJWU1qjt/uqEipN17R0cJJw+P5caPS3jmuES2Fyn/+KqUq90YTh4Wyz+/LuV3g2JIjhP+9kXN3+2oQTE8uriUG3HWRVy0cOgeMXy6oYLf7+UkDlWlpIKqIs/ickVwqjYDnD06ln99U8oRA2OIjRIeWFTK0YOcZX5vbTm/G1T9m/6a52PKMwVcOj6OmWNrHvBHdI0i80/V93e+zKzg0veKWXphclVNuW9+rSA9LarO1XIg7eNh41XV083MVcb/p4AlFyTTJVmq4tpepOzXO/jkFgnaRBKqfE5o6dKlHHDAAR5HY3aHde9d18INFRz8dPV0E+/KZ1K/6Kri0IePSuCCd4ro+Y980hKE8/eJ45y9q5fhhVVlPHd84OKl3fWvIxO47P1i9nhoJwkx1Ijh/DGx/LjNx8hHC2gXD5ePj2PB+goqV8Mxg2O4Ym4xG/N99HSL0y4cE8vD35by+72caWzIVfr/s/qZr8S78unXXlh/hfPb3nxQPFsLlcH/2klCjJsU3e1mztpyHj26evn/s7SMn3Yoty4o4dYF1Qlx5w3tiIkSuqdU/z4dE51KEv4VCZ5fUcbMsY1fPfsTqTnd4nLnyq9bilR1yvjflWWcOSq2KrG2Fp51790UTe3e+5577uH666/noYce4rLLLgtBZJHPuvfefZHYvfc7P5Tx7IoyXj4pfGJ+f20ZM+cUs+GK6pOD2UtK+d8WHw8eUV0MeMATBTx8ZELVA6tNkb3Tx96zCvj1ypRmqT27ucDHpKcKWXZhMgnNmCxKypVRjxaw8OwkuiY3fIVl3XuHqcpEa1W0Wz/r3nvXHDMklmOG7NpZe3MrKlM+WV/O1AExZO9Ubvu0lOP3rBnTBWPq3gf54pzd78Awt0S5f2pCsx0buiZHsfqS+h+A3h3xMcKaS5t/uuEgsvaY3WRJqPWbtaSUs94qIlpgUnoM/3dUK6tK1AopzsOkM14tIjFG+N3gGG4/uPF7Yc1hcKdoBndqXfdYIk2bSkKm9bPuvSNPUqzw7fmt8yzfNK5NPKzatWtXAEaOHOlxJGFMFZU2sTkY0yr5JAY08lrZbhNHnaQk56ZrZTIydcVIOaWJ3bwOwxjTRIVpQ/AV19+CeDhrE0mosNCpvrplyxaPIwlfnbt045eRf8IX1boehDOmtfNJDDs7DOfHcXey6ZefUFWioiLn0N4m7gllZ2cDsHLlSiZOnOhxNOGpR+9+rN46iqVHzkGi7EatMRFDffiK89i0PoPsDT8QExtLuw6dvY4qaJ4kIREZDTwKJADlwMWqGrIe5yLpWSivREVFMXyf/Vj1zULmv/qM1SQ0JtKIEBMTy3Fn/4n4xPB57qsxXl0J3Qvcpqrvi8hR7vvJHsVi/IwYfxB9Bw5jZ94OfL7Iu8lpTFsVFR1NWqeuJKUE18J6uPAqCSlQuabaAxtbYqZ2dh+cdh07065j5FzOG2Mil1dJ6ArgAxH5O07liP0bGlFELgAuAOjbt2+LBGeMMaZlhCwJicg8oHs9g24EDgH+pKqvicjJwOPAofVNR1VnA7PBaTuuKbF07+6EMXr06KZ83RhjTIh40oCpiOQCaaqq4pSR5apqowWZIrIF2NDE2XYGtjbxu16IpHgjKVaIrHgjKVaIrHjbUqz9VLVLcwXTnLwqjtsITAIWAFOAtcF8aXdWoogsDtdWZOsTSfFGUqwQWfFGUqwQWfFarOHBqyR0PvBPEYkBinHv+RhjjGlbPElCqvo5MMaLeRtjjAkfkdO2w+6b7XUAuyiS4o2kWCGy4o2kWCGy4rVYw0BE9axqjDGmdWlLV0LGGGPCjCUhY4wxnmkTSUhEjhCRH0RknYhc53U8gYjIEyKyWURWeR1LY0Skj4h8IiL/E5HvReSPXsfUEBFJEJFvROQ7N9bbvI6pMSISLSLLRORdr2NpjIisF5GVIrJcRBZ7HU9jRCRNRF4VkTUislpEJngdU31EZIi7Tiv/8kTkCq/jak6t/p6QiEQDPwKHAVnAt8Cpqvo/TwNrgIgcBOwEnlHVEV7HE4iI9AB6qOpSEUkFlgDHheO6dR+KTlbVnSISC3wO/FFVF3kcWoNE5EpgLNBOVY/2Op5ARGQ9MFZVI+LhTxF5GvhMVf8jInFAkqrmeBxWQO6x7FdgX1Vt6kP7YactXAmNB9ap6k+qWgq8CEzzOKYGqepCYLvXcQRDVX9T1aXu63xgNdDL26jqp46d7ttY9y9sz8BEpDfwO+A/XsfS2ohIe+AgnObCUNXScE9ArkOAjNaUgKBtJKFeQKbf+yzC9EAZyUQkHdgb+NrjUBrkFm8tBzYDH6lq2MYKPAhcA0RKfxoKfCgiS9xGh8NZf2AL8KRb3PkfEUn2OqggnAK84HUQza0tJCETYiKSArwGXKGqYdvJvapWqOpooDcwXkTCsrhTRI4GNqvqEq9j2QUTVXUf4EjgErdYOVzFAPsA/1bVvYECINzvFccBxwKveB1Lc2sLSehXoI/f+97uZ6YZuPdXXgOeV9XXvY4nGG7RyyfAER6H0pADgGPd+ywvAlNE5DlvQwpMVX91/28G3sApBg9XWUCW35XwqzhJKZwdCSxV1WyvA2lubSEJfQsMEpH+7tnEKcDbHsfUKrg3+x8HVqvqP7yOJxAR6SIiae7rRJyKKms8DaoBqnq9qvZW1XSc7fVjVf2Dx2E1SESS3YopuMVaU4Gwrd2pqpuATBEZ4n50CBB2lWlqOZVWWBQH3jVg2mJUtVxELgU+AKKBJ1T1e4/DapCIvIDT1XlnEckCblHVx72NqkEHAKcDK917LQA3qOp73oXUoB7A024NoyjgZVUN+6rPEaIb8Ibbc3EM8F9VnettSI26DHjePTH9CTjb43ga5Cb2w4ALvY4lFFp9FW1jjDHhqy0UxxljjAlTloSMMcZ4xpKQMcYYz1gSMsYY4xlLQsYYYzxjSciYBrgtLV/svu4pIq96HZMxrY1V0TamAW57eO+Ge2vmxkSyVv+wqjG74R5ggPsg7lpgqKqOEJGzgOOAZGAQ8HcgDufB3RLgKFXdLiIDgEeALkAhcL6qhmUrDcZ4xYrjjGnYdThN548G/lxr2AjgBGAccBdQ6DaG+RVwhjvObOAyVR0DXA38X0sEbUwksSshY5rmE7cPpXwRyQXecT9fCYx0WxbfH3jFbc4GIL7lwzQmvFkSMqZpSvxe+/ze+3D2qyggx72KMsY0wIrjjGlYPpDalC+6/Sr9LCIngdPiuIiMas7gjGkNLAkZ0wBV3QZ8ISKrgPuaMInTgHNF5Dvge8K4W3ljvGJVtI0xxnjGroSMMcZ4xpKQMcYYz1gSMsYY4xlLQsYYYzxjScgYY4xnLAkZY4zxjCUhY4wxnvl/qa5p94x+np8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import mpmath as mp\n",
    "import csv\n",
    "\n",
    "def file_reader(filename):\n",
    "    with open(filename) as file:\n",
    "        reader = csv.reader(file, delimiter=\" \")\n",
    "        data  = list(zip(*reader))\n",
    "        # data is a tuple of strings. Tuples are immutable, and we need to perform math on\n",
    "        #   the data, so here we convert tuple to lists of floats:\n",
    "        data0 = []\n",
    "        data1 = []\n",
    "        for i in range(len(data[0])):\n",
    "            data0.append(float(data[0][i]))\n",
    "            data1.append(float(data[1][i]))\n",
    "        return data0,data1\n",
    "\n",
    "first_col16,second_col16 = file_reader(os.path.join(outdir,'out-lowresolution.txt'))\n",
    "first_col24,second_col24 = file_reader(os.path.join(outdir,'out-medresolution.txt'))\n",
    "\n",
    "second_col16_rescaled4o = []\n",
    "second_col16_rescaled5o = []\n",
    "for i in range(len(second_col16)):\n",
    "    # data16 = data24*(16/24)**4\n",
    "    # -> log10(data24) = log10(data24) + 4*log10(16/24)\n",
    "    second_col16_rescaled4o.append(second_col16[i] + 4*mp.log10(16./24.))\n",
    "    second_col16_rescaled5o.append(second_col16[i] + 5*mp.log10(16./24.))\n",
    "\n",
    "# https://matplotlib.org/gallery/text_labels_and_annotations/legend.html#sphx-glr-gallery-text-labels-and-annotations-legend-py\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "plt.title(\"Demonstrating 4th-order Convergence: \"+par.parval_from_str(\"reference_metric::CoordSystem\")+\" Coordinates\")\n",
    "plt.xlabel(\"time\")\n",
    "plt.ylabel(\"log10(Max relative error)\")\n",
    "\n",
    "ax.plot(first_col24, second_col24, 'k-', label='logErel(N0=24)')\n",
    "ax.plot(first_col16, second_col16_rescaled4o, 'k--', label='logErel(N0=16) + log((16/24)^4)')\n",
    "ax.set_ylim([-8.05,-1.7]) # Manually set the y-axis range case, since the log10\n",
    "                          # relative error at t=0 could be -inf or about -16,\n",
    "                          # resulting in very different-looking plots\n",
    "                          # despite the data being the same to roundoff.\n",
    "if par.parval_from_str(\"reference_metric::CoordSystem\") == \"Cartesian\":\n",
    "    ax.set_ylim([-2.68,-1.62])\n",
    "if par.parval_from_str(\"reference_metric::CoordSystem\") == \"Cylindrical\":\n",
    "    ax.plot(first_col16, second_col16_rescaled5o, 'k.',  label='(Assuming 5th-order convergence)')\n",
    "legend = ax.legend(loc='lower right', shadow=True, fontsize='large')\n",
    "legend.get_frame().set_facecolor('C1')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 5: Output this notebook to $\\LaTeX$-formatted PDF file \\[Back to [top](#toc)\\]\n",
    "$$\\label{latex_pdf_output}$$\n",
    "\n",
    "The following code cell converts this Jupyter notebook into a proper, clickable $\\LaTeX$-formatted PDF file. After the cell is successfully run, the generated PDF may be found in the root NRPy+ tutorial directory, with filename\n",
    "[Tutorial-Start_to_Finish-ScalarWaveCurvilinear.pdf](Tutorial-Start_to_Finish-ScalarWaveCurvilinear.pdf) (Note that clicking on this link may not work; you may need to open the PDF file through another means.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:32:36.954666Z",
     "iopub.status.busy": "2021-03-07T17:32:36.954008Z",
     "iopub.status.idle": "2021-03-07T17:32:41.056967Z",
     "shell.execute_reply": "2021-03-07T17:32:41.056210Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-Start_to_Finish-ScalarWaveCurvilinear.tex, and compiled\n",
      "    LaTeX file to PDF file Tutorial-Start_to_Finish-\n",
      "    ScalarWaveCurvilinear.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-Start_to_Finish-ScalarWaveCurvilinear\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}