Deprecated/Tutorial-Start_to_Finish-ScalarWave.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 Cartesian Coordinates\n",
    "\n",
    "## Author: Zach Etienne\n",
    "### Formatting improvements courtesy Brandon Clark\n",
    "\n",
    "## This module solves the scalar wave equation in Cartesian coordinates, using the [Method of Lines](Tutorial-Method_of_Lines-C_Code_Generation.ipynb).\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/ScalarWave_RHSs.py](../edit/ScalarWave/ScalarWave_RHSs.py) [\\[**tutorial**\\]](Tutorial-ScalarWave.ipynb) Generates the right-hand side for the Scalar Wave Equation in cartesian coordinates\n",
    "* [ScalarWave/InitialData.py](../edit/ScalarWave/InitialData.py) [\\[**tutorial**\\]](Tutorial-ScalarWave.ipynb) Generates C code for plane wave or spherical Gaussian initial data for the scalar wave equation\n",
    "\n",
    "## Introduction:\n",
    "\n",
    "As outlined in the [previous NRPy+ tutorial notebook](Tutorial-ScalarWave.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 as follows, with links to the relevant NRPy+ tutorial notebooks listed at each step:\n",
    "\n",
    "1. Allocate memory for gridfunctions, including temporary storage for the Method of Lines time integration\n",
    "    * [**NRPy+ tutorial notebook on Method of Lines algorithm**](Tutorial-Method_of_Lines-C_Code_Generation.ipynb).\n",
    "1. Set gridfunction values to initial data \n",
    "    * [**NRPy+ tutorial notebook section on plane-wave solution to scalar wave equation**](Tutorial-ScalarWave.ipynb#planewavesoln)\n",
    "1. Next, integrate the initial data forward in time using the Method of Lines coupled to a Runge-Kutta explicit timestepping algorithm:\n",
    "    1. At the start of each iteration in time, output the difference between the numerical and exact solution\n",
    "        * [**NRPy+ tutorial notebook section on plane-wave solution to scalar wave equation**](Tutorial-ScalarWave.ipynb#planewavesoln).\n",
    "    1. At each RK time substep, do the following:\n",
    "        1. Evaluate scalar wave RHS expressions \n",
    "            * [**NRPy+ tutorial notebook section on right-hand sides of scalar wave equation, in 3 spatial dimensions**](Tutorial-ScalarWave.ipynb#rhss3d)\n",
    "        1. Apply boundary conditions [*a la* the SENR/NRPy+ paper](https://arxiv.org/abs/1712.07658)\n",
    "1. Repeat above steps at two numerical resolutions to confirm convergence to zero."
   ]
  },
  {
   "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](#setup): Set up core functions and parameters for solving scalar wave equation\n",
    "    1. [Step 1.a](#applybcs) `apply_bcs()`: outer boundary condition driver function\n",
    "    1. [Step 1.b](#mol) Generate Runge-Kutta-based Method of Lines timestepping code\n",
    "    1. [Step 1.c](#freeparams) Output C codes needed for declaring and setting Cparameters; also set `free_parameters.h`\n",
    "1. [Step 2](#mainc): `ScalarWave_Playground.c`: The Main C Code\n",
    "1. [Step 3](#convergence): Code validation: Verify that relative error in numerical solution converges to zero at the expected order\n",
    "1. [Step 4](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='setup'></a>\n",
    "\n",
    "# Step 1: Set up core functions and parameters for solving scalar wave equation \\[Back to [top](#toc)\\]\n",
    "$$\\label{setup}$$\n",
    "\n",
    "Let's pick up where we left off in the [previous module](Tutorial-ScalarWave.ipynb), interfacing with the [ScalarWave/InitialData](../edit/ScalarWave/InitialData.py) and [ScalarWave/ScalarWave_RHSs](../edit/ScalarWave/ScalarWave_RHSs.py) NRPy+ modules to generate\n",
    "* monochromatic (single-wavelength) plane wave scalar wave initial data, and\n",
    "* the scalar wave equation RHSs at **4th** finite difference order in **3 spatial dimensions**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:32:43.423384Z",
     "iopub.status.busy": "2021-03-07T17:32:43.416415Z",
     "iopub.status.idle": "2021-03-07T17:32:44.006636Z",
     "shell.execute_reply": "2021-03-07T17:32:44.007147Z"
    }
   },
   "outputs": [],
   "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 indexedexp as ixp              # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) 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(\"ScalarWave_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",
    "# !rm -r ScalarWaveCurvilinear_Playground_Ccodes\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 P4: Set domain_size, the physical extent of numerical grid;\n",
    "#          in Cartesian coordinates xmin=ymin=zmin=-domain_size,\n",
    "#          and xmax=ymax=zmax=+domain_size\n",
    "domain_size = 10.0\n",
    "\n",
    "# Step P5: Set timestepping algorithm (we adopt the Method of Lines)\n",
    "RK_method = \"RK4\"\n",
    "\n",
    "# Step P6: Set the finite differencing order to 4.\n",
    "par.set_parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\",4)\n",
    "\n",
    "# Step 1: 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 2: Import ScalarWave_RHSs module.\n",
    "#         This command only declares ScalarWave RHS parameters\n",
    "#         and the ScalarWave_RHSs function (called later)\n",
    "import ScalarWave.ScalarWave_RHSs as swrhs\n",
    "\n",
    "# Step 3: Set the spatial dimension parameter\n",
    "#         to 3, and then read the parameter as DIM.\n",
    "par.set_parval_from_str(\"grid::DIM\",3)\n",
    "DIM = par.parval_from_str(\"grid::DIM\")\n",
    "\n",
    "# Step 4: 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(Type=\"PlaneWave\")\n",
    "\n",
    "# Step 5: 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",
    "swrhs.ScalarWave_RHSs()\n",
    "\n",
    "# Step 6: Enable \"FD functions\". In other words, all finite-difference stencils\n",
    "#         will be output as inlined static functions. This is essential for\n",
    "#         compiling highly complex FD kernels with using certain versions of GCC;\n",
    "#         GCC 10-ish will choke on BSSN FD kernels at high FD order, sometimes\n",
    "#         taking *hours* to compile. Unaffected GCC versions compile these kernels\n",
    "#         in seconds. FD functions do not slow the code performance, but do add\n",
    "#         another header file to the C source tree.\n",
    "enable_FD_functions = True\n",
    "par.set_parval_from_str(\"finite_difference::enable_FD_functions\", enable_FD_functions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:32:43.423384Z",
     "iopub.status.busy": "2021-03-07T17:32:43.416415Z",
     "iopub.status.idle": "2021-03-07T17:32:44.006636Z",
     "shell.execute_reply": "2021-03-07T17:32:44.007147Z"
    },
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Output C function exact_solution_single_point() to file ScalarWave_Ccodes/exact_solution_single_point.h\n",
      "Output C function exact_solution_all_points() to file ScalarWave_Ccodes/exact_solution_all_points.h\n",
      "Output C function rhs_eval() to file ScalarWave_Ccodes/rhs_eval.h\n"
     ]
    }
   ],
   "source": [
    "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",
    "\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(xx0,xx1,xx2,params,&in_gfs[IDX4S(UUGF,i0,i1,i2)],&in_gfs[IDX4S(VVGF,i0,i1,i2)]);\",\n",
    "    loopopts = \"AllPoints,Read_xxs\")\n",
    "\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   =\"const paramstruct *restrict params, 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\")\n",
    "\n",
    "# Step 6.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='applybcs'></a>\n",
    "\n",
    "## Step 1.a: `apply_bcs()`: outer boundary condition driver function\n",
    "$$\\label{applybcs}$$\n",
    "\n",
    "When solving the wave equation on a 3D Cartesian numerical grid cube (or, if you like, rectangular prism), at each step in time, we first evaluate the right-hand sides (RHSs) of the $\\partial_t u$ and $\\partial_t v$ equations. \n",
    "\n",
    "These RHSs generally contain spatial derivatives, which we evaluate using finite-difference differentiation ([**tutorial**](Tutorial-Finite_Difference_Derivatives.ipynb)). Each finite-difference derivative depends on neighboring points on the left and right, so the RHSs can only be evaluated in the grid interior. For example, a standard fourth-order centered finite difference derivative depends on two points to the left and right of the point at which the derivative is being evaluated. In order for the same interior to be filled at the next time step, we need to fill in the data at the boundaries; i.e., we need to apply boundary conditions.\n",
    "\n",
    "Here we quadratically extrapolate data to the outer boundary using the `FACE_UPDATE()` C macro defined below. The C code function `apply_bcs()` below updates all 6 faces of the cube. To ensure that all gridpoints on the outer boundary (also known as \"ghost cells\") are filled, the following algorithm is implemented, starting at the innermost ghost cells (i.e., the ghost cells closest to the grid interior):\n",
    "\n",
    "1. The lower $x$ face is updated on only the interior points of the face.\n",
    "1. The upper $x$ face is updated on only the interior points of the face.\n",
    "1. The lower $y$ face is updated on the interior points of that face, plus the lower and upper $x$ boundary points\n",
    "1. The upper $y$ face is updated on the interior points of that face, plus the lower and upper $x$ boundary points\n",
    "1. The lower $z$ face is updated on the interior points of that face, plus the lower and upper $x$ boundary points, plus the lower and upper $y$ boundary points\n",
    "1. The upper $z$ face is updated on the interior points of that face, plus the lower and upper $x$ boundary points, plus the lower and upper $y$ boundary points\n",
    "1. The above is repeated on the next outer ghost cell, until all outer boundary points are filled."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:32:44.013381Z",
     "iopub.status.busy": "2021-03-07T17:32:44.012289Z",
     "iopub.status.idle": "2021-03-07T17:32:44.016350Z",
     "shell.execute_reply": "2021-03-07T17:32:44.016939Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing ScalarWave_Ccodes//apply_bcs.h\n"
     ]
    }
   ],
   "source": [
    "%%writefile $Ccodesdir/apply_bcs.h\n",
    "\n",
    "// Declare boundary condition FACE_UPDATE macro,\n",
    "//          which updates a single face of the 3D grid cube\n",
    "//          using quadratic polynomial extrapolation.\n",
    "const int MAXFACE = -1;\n",
    "const int NUL     = +0;\n",
    "const int MINFACE = +1;\n",
    "#define  FACE_UPDATE(which_gf, i0min,i0max, i1min,i1max, i2min,i2max, FACEX0,FACEX1,FACEX2) \\\n",
    "  for(int i2=i2min;i2<i2max;i2++) for(int i1=i1min;i1<i1max;i1++) for(int i0=i0min;i0<i0max;i0++) { \\\n",
    "        gfs[IDX4S(which_gf,i0,i1,i2)] =                                  \\\n",
    "          +3.0*gfs[IDX4S(which_gf,i0+1*FACEX0,i1+1*FACEX1,i2+1*FACEX2)]  \\\n",
    "          -3.0*gfs[IDX4S(which_gf,i0+2*FACEX0,i1+2*FACEX1,i2+2*FACEX2)]  \\\n",
    "          +1.0*gfs[IDX4S(which_gf,i0+3*FACEX0,i1+3*FACEX1,i2+3*FACEX2)]; \\\n",
    "      }\n",
    "\n",
    "void apply_bcs(const paramstruct *restrict params,REAL *restrict gfs) {\n",
    "#pragma omp parallel for\n",
    "    for(int which_gf=0;which_gf<NUM_EVOL_GFS;which_gf++) {\n",
    "#include \"set_Cparameters.h\"\n",
    "        int imin[3] = { NGHOSTS, NGHOSTS, NGHOSTS };\n",
    "        int imax[3] = { Nxx_plus_2NGHOSTS0-NGHOSTS, Nxx_plus_2NGHOSTS1-NGHOSTS, Nxx_plus_2NGHOSTS2-NGHOSTS };\n",
    "        for(int which_gz = 0; which_gz < NGHOSTS; which_gz++) {\n",
    "            // After updating each face, adjust imin[] and imax[]\n",
    "            //   to reflect the newly-updated face extents.\n",
    "            FACE_UPDATE(which_gf, imin[0]-1,imin[0], imin[1],imax[1], imin[2],imax[2], MINFACE,NUL,NUL); imin[0]--;\n",
    "            FACE_UPDATE(which_gf, imax[0],imax[0]+1, imin[1],imax[1], imin[2],imax[2], MAXFACE,NUL,NUL); imax[0]++;\n",
    "\n",
    "            FACE_UPDATE(which_gf, imin[0],imax[0], imin[1]-1,imin[1], imin[2],imax[2], NUL,MINFACE,NUL); imin[1]--;\n",
    "            FACE_UPDATE(which_gf, imin[0],imax[0], imax[1],imax[1]+1, imin[2],imax[2], NUL,MAXFACE,NUL); imax[1]++;\n",
    "\n",
    "            FACE_UPDATE(which_gf, imin[0],imax[0], imin[1],imax[1], imin[2]-1,imin[2], NUL,NUL,MINFACE); imin[2]--;\n",
    "            FACE_UPDATE(which_gf, imin[0],imax[0], imin[1],imax[1], imax[2],imax[2]+1, NUL,NUL,MAXFACE); imax[2]++;\n",
    "        }\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='mol'></a>\n",
    "\n",
    "## Step 1.b: Generate Runge-Kutta-based Method of Lines timestepping code \\[Back to [top](#toc)\\]\n",
    "$$\\label{mol}$$\n",
    "\n",
    "The Method of Lines algorithm is described in detail in the [**NRPy+ tutorial notebook on Method of Lines algorithm**](Tutorial-Method_of_Lines-C_Code_Generation.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:32:44.023378Z",
     "iopub.status.busy": "2021-03-07T17:32:44.022377Z",
     "iopub.status.idle": "2021-03-07T17:32:44.037707Z",
     "shell.execute_reply": "2021-03-07T17:32:44.036956Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 1.b: Generate Runge-Kutta-based (RK-based) timestepping code.\n",
    "#       As described above the Table of Contents, this is a 2-step process:\n",
    "#       1.b.A: Evaluate RHSs (RHS_string)\n",
    "#       1.b.B: Apply boundary conditions (post_RHS_string, pt 1)\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",
    "MoL.MoL_C_Code_Generation(RK_method,\n",
    "    RHS_string      = \"rhs_eval(&params,   RK_INPUT_GFS, RK_OUTPUT_GFS);\",\n",
    "    post_RHS_string = \"apply_bcs(&params,  RK_OUTPUT_GFS);\",\n",
    "    outdir = os.path.join(Ccodesdir,\"MoLtimestepping/\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='freeparams'></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{freeparams}$$\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": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:32:44.044871Z",
     "iopub.status.busy": "2021-03-07T17:32:44.044110Z",
     "iopub.status.idle": "2021-03-07T17:32:44.048339Z",
     "shell.execute_reply": "2021-03-07T17:32:44.048859Z"
    }
   },
   "outputs": [],
   "source": [
    "domain_size_str=str(domain_size)\n",
    "# Step 3.d.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",
    "\n",
    "// Set free-parameter values for the initial data.\n",
    "params.time = 0.0; params.wavespeed = 1.0;\n",
    "//params.kk0 = 1.0; params.kk1 = 1.0; params.kk2 = 1.0;\n",
    "\n",
    "const REAL domain_size    = \"\"\"+str(domain_size)+\"\"\";\n",
    "\n",
    "// Override parameter defaults with values based on command line arguments and NGHOSTS.\n",
    "const int Nx0x1x2 = atoi(argv[1]);\n",
    "params.Nxx0 = Nx0x1x2;\n",
    "params.Nxx1 = Nx0x1x2;\n",
    "params.Nxx2 = Nx0x1x2;\n",
    "params.Nxx_plus_2NGHOSTS0 = params.Nxx0 + 2*NGHOSTS;\n",
    "params.Nxx_plus_2NGHOSTS1 = params.Nxx1 + 2*NGHOSTS;\n",
    "params.Nxx_plus_2NGHOSTS2 = params.Nxx2 + 2*NGHOSTS;\n",
    "// Step 0d: Set up space and time coordinates\n",
    "// Step 0d.i: Declare \\Delta x^i=dxx{0,1,2} and invdxx{0,1,2}, as well as xxmin[3] and xxmax[3]:\n",
    "const REAL xxmin[3] = {-\"\"\"+domain_size_str+\"\"\",-\"\"\"+domain_size_str+\"\"\",-\"\"\"+domain_size_str+\"\"\" };\n",
    "const REAL xxmax[3] = {+\"\"\"+domain_size_str+\"\"\",+\"\"\"+domain_size_str+\"\"\",+\"\"\"+domain_size_str+\"\"\" };\n",
    "\n",
    "params.dxx0 = (xxmax[0] - xxmin[0]) / ((REAL)params.Nxx0);\n",
    "params.dxx1 = (xxmax[1] - xxmin[1]) / ((REAL)params.Nxx1);\n",
    "params.dxx2 = (xxmax[2] - xxmin[2]) / ((REAL)params.Nxx2);\n",
    "params.invdx0 = 1.0 / params.dxx0;\n",
    "params.invdx1 = 1.0 / params.dxx1;\n",
    "params.invdx2 = 1.0 / params.dxx2;\n",
    "\\n\"\"\")\n",
    "\n",
    "# Generates declare_Cparameters_struct.h, set_Cparameters_default.h, and set_Cparameters[-SIMD].h\n",
    "import deprecated_NRPy_param_funcs as evil_par\n",
    "evil_par.generate_Cparameters_Ccodes(os.path.join(Ccodesdir))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='mainc'></a>\n",
    "\n",
    "# Step 2: `ScalarWave_Playground.c`: The Main C Code \\[Back to [top](#toc)\\]\n",
    "$$\\label{mainc}$$\n",
    "\n",
    "Next we will write the C code infrastructure necessary to make use of the above NRPy+-generated codes. Again, we'll be using RK4 time integration via the Method of Lines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:32:44.054481Z",
     "iopub.status.busy": "2021-03-07T17:32:44.053716Z",
     "iopub.status.idle": "2021-03-07T17:32:44.056916Z",
     "shell.execute_reply": "2021-03-07T17:32:44.056340Z"
    }
   },
   "outputs": [],
   "source": [
    "# Part P0: Set the number of ghost cells, from NRPy+'s FD_CENTDERIVS_ORDER\n",
    "with open(os.path.join(Ccodesdir,\"ScalarWave_NGHOSTS.h\"), \"w\") as file:\n",
    "    file.write(\"// Part P0: Set the number of ghost cells, from NRPy+'s FD_CENTDERIVS_ORDER\\n\")\n",
    "    file.write(\"#define NGHOSTS \"+str(int(par.parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\")/2))+\"\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:32:44.063546Z",
     "iopub.status.busy": "2021-03-07T17:32:44.062595Z",
     "iopub.status.idle": "2021-03-07T17:32:44.066404Z",
     "shell.execute_reply": "2021-03-07T17:32:44.065774Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing ScalarWave_Ccodes//ScalarWave_Playground.c\n"
     ]
    }
   ],
   "source": [
    "%%writefile $Ccodesdir/ScalarWave_Playground.c\n",
    "\n",
    "// Part P0: Import NGHOSTS, which is based on FD_CENTDERIVS_ORDER\n",
    "#include \"ScalarWave_NGHOSTS.h\"\n",
    "// Part P0a: set REAL=double, so that all floating point numbers are stored to at least ~16 significant digits.\n",
    "#define REAL double\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",
    "const int NSKIP_0D_OUTPUT = 1;\n",
    "const int NSKIP_2D_OUTPUT = 10;\n",
    "\n",
    "// Part P1: Import needed header files\n",
    "#include \"stdio.h\"\n",
    "#include \"stdlib.h\"\n",
    "#include \"math.h\"\n",
    "\n",
    "// Part P2: Add needed #define's to set data type, the IDX4S() macro, and the gridfunctions\n",
    "// Part P2a: 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 LOOP_ALL_GFS_GPS(ii) _Pragma(\"omp parallel for\") \\\n",
    "  for(int (ii)=0;(ii)<Nxx_plus_2NGHOSTS_tot*NUM_EVOL_GFS;(ii)++)\n",
    "// Part P2b: Set UUGF and VVGF macros\n",
    "#define NUM_EVOL_GFS 2\n",
    "#define UUGF 0\n",
    "#define VVGF 1\n",
    "\n",
    "// Part P3: 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 P4: 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 P5: Declare the function to evaluate the scalar wave RHSs\n",
    "#include \"rhs_eval.h\"\n",
    "\n",
    "// Part P6: Boundary condition driver routine: Apply BCs to all six\n",
    "//          boundary faces of the cube, filling in the innermost\n",
    "//          ghost zone first, and moving outward.\n",
    "#include \"apply_bcs.h\"\n",
    "\n",
    "// Part P7: 2D File output routine, for comparing numerical results to exact solution\n",
    "void output_2D(const paramstruct *restrict params, const int iter, REAL *xx[3],\n",
    "               const REAL *numerical_gridfunction_data,REAL *gridfunction_to_store_exact) {\n",
    "#include \"set_Cparameters.h\"\n",
    "    // Step 2a: Validation: Output to 2D data files numerical and exact solutions\n",
    "    char filename[100];\n",
    "    sprintf(filename,\"out2D__resolution_%dx%dx%d__iter_%d.txt\",Nxx0,Nxx1,Nxx2,iter);\n",
    "    FILE *out2D = fopen(filename, \"w\");\n",
    "    const int i2  = (int)((Nxx2+ 2*NGHOSTS)*0.5);\n",
    "    const REAL xx2 = xx[2][i2];\n",
    "    for(int i0=0;i0<Nxx0+2*NGHOSTS;i0++) {\n",
    "        for(int i1=0;i1<Nxx1+2*NGHOSTS;i1++) {\n",
    "            if(i0> (Nxx0+2*NGHOSTS)*.25 && i0< (Nxx0+2*NGHOSTS)*.75 &&\n",
    "               i1> (Nxx1+2*NGHOSTS)*.25 && i1< (Nxx1+2*NGHOSTS)*.75) {\n",
    "                const REAL xx0 = xx[0][i0];\n",
    "                const REAL xx1 = xx[1][i1];\n",
    "                REAL uu_exact,vv_exact; exact_solution_single_point(xx0,xx1,xx2,params, &uu_exact,&vv_exact);\n",
    "                fprintf(out2D,\"%e %e %e %e\\n\", xx0,xx1,\n",
    "                        numerical_gridfunction_data[IDX4S(0,i0,i1, (int)((Nxx2+ 2*NGHOSTS)*0.5))], uu_exact);\n",
    "            }\n",
    "        }\n",
    "    }\n",
    "    fclose(out2D);\n",
    "}\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: Evolve scalar wave initial data forward in time using Method of Lines with RK4 algorithm,\n",
    "//         applying quadratic extrapolation outer boundary conditions.\n",
    "// Step 3: Output relative error between numerical and exact solution.\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",
    "    // Step 0a: Read command-line input, error out if nonconformant\n",
    "    if(argc != 2 || atoi(argv[1]) < NGHOSTS) {\n",
    "      printf(\"Error: Expected one command-line argument: ./ScalarWave_Playground [Nx(=Ny=Nz)],\\n\");\n",
    "      printf(\"where Nx is the number of grid points in the x,y, and z directions.\\n\");\n",
    "      printf(\"Nx MUST BE larger than NGHOSTS (= %d)\\n\",NGHOSTS);\n",
    "      exit(1);\n",
    "    }\n",
    "    if(atoi(argv[1])%2 != 0) {\n",
    "        printf(\"Error: Algorithm for setting up cell-centered grids here requires Nx, Ny, and Nz to be a multiple of 2 .\\n\");\n",
    "        exit(1);\n",
    "    }\n",
    "\n",
    "    // Step 0b: Set free parameters, overwriting Cparameters defaults\n",
    "    //          by hand or with command-line input, as desired.\n",
    "#include \"free_parameters.h\"\n",
    "    //          ... and then set up the numerical grid structure in time:\n",
    "    const REAL CFL_FACTOR = 0.5; // Set the CFL Factor\n",
    "    #define MIN(A, B) ( ((A) < (B)) ? (A) : (B) )\n",
    "    REAL dt = CFL_FACTOR * MIN(params.dxx0,MIN(params.dxx1,params.dxx2)); // CFL condition\n",
    "\n",
    "    // Now that params struct has been properly set up, create\n",
    "    //   list of const's containing each parameter. E.g.,\n",
    "    //  const REAL dxx0 = params.dxx0;\n",
    "#include \"set_Cparameters-nopointer.h\"\n",
    "\n",
    "    // Step 0c: Allocate memory for gridfunctions\n",
    "    const int Nxx_plus_2NGHOSTS_tot = Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2;\n",
    "    // Step 0d: Allocate memory for gridfunctions\n",
    "#include \"MoLtimestepping/RK_Allocate_Memory.h\"\n",
    "\n",
    "    // Step 0e: Set t_final, and number of timesteps based on t_final\n",
    "    const REAL t_final = xxmax[0]*0.8; /* 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",
    "    int Nt = (int)(t_final / dt + 0.5); // The number of points in time.\n",
    "                                        //Add 0.5 to account for C rounding down integers.\n",
    "\n",
    "    // Step 0f: Set up cell-centered Cartesian coordinate grids\n",
    "    REAL *xx[3];\n",
    "    xx[0] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS0);\n",
    "    xx[1] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS1);\n",
    "    xx[2] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS2);\n",
    "    for(int j=0;j<Nxx_plus_2NGHOSTS0;j++) xx[0][j] = xxmin[0] + ((REAL)(j-NGHOSTS) + (1.0/2.0))*params.dxx0;\n",
    "    for(int j=0;j<Nxx_plus_2NGHOSTS1;j++) xx[1][j] = xxmin[1] + ((REAL)(j-NGHOSTS) + (1.0/2.0))*params.dxx1;\n",
    "    for(int j=0;j<Nxx_plus_2NGHOSTS2;j++) xx[2][j] = xxmin[2] + ((REAL)(j-NGHOSTS) + (1.0/2.0))*params.dxx2;\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<=Nt;n++) { // Main loop to progress forward in time.\n",
    "        // Step 1a: Set current time to correct value & compute exact solution\n",
    "        params.time = ((REAL)n)*dt;\n",
    "\n",
    "        /* Step 2: Validation: Output relative error between numerical and exact solution, */\n",
    "        // Step 2b: Output to 2D grid (region of x-y plane near origin)\n",
    "        //          every NSKIP_2D_OUTPUT iterations.\n",
    "        if((n)%NSKIP_2D_OUTPUT ==0) {\n",
    "            // Step 2a: Output data on gridpoints closest to xy-plane, ignoring points \"too close\" to outer boundary\n",
    "            output_2D(&params, n, xx, y_n_gfs, diagnostic_output_gfs);\n",
    "        }\n",
    "\n",
    "        if((n)%NSKIP_0D_OUTPUT ==0) {\n",
    "            // Step 2c: Output relative error between exact & numerical at center of grid.\n",
    "            const int i0mid=Nxx_plus_2NGHOSTS0/2;\n",
    "            const int i1mid=Nxx_plus_2NGHOSTS1/2;\n",
    "            const int i2mid=Nxx_plus_2NGHOSTS2/2;\n",
    "            REAL uu_exact,vv_exact; exact_solution_single_point(xx[0][i0mid],xx[1][i1mid],xx[2][i2mid],&params,\n",
    "                                                                &uu_exact,&vv_exact);\n",
    "\n",
    "            const REAL numerical = y_n_gfs[IDX4S(UUGF,i0mid,i1mid,i2mid)];\n",
    "            const REAL relative_error = fabs((uu_exact-numerical)/uu_exact) + 1e-16; // Add 1e-16 in case relerror\n",
    "                                                                                     // is exactly zero.\n",
    "            printf(\"%e %e || %e %e %e: %e %e\\n\",params.time, log10(relative_error),\n",
    "                   xx[0][i0mid],xx[1][i1mid],xx[2][i2mid], numerical,uu_exact);\n",
    "        }\n",
    "\n",
    "        // Step 3: Evolve scalar wave initial data forward in time using Method of Lines with RK4 algorithm,\n",
    "        //         applying quadratic extrapolation outer boundary conditions.\n",
    "        // Step 3.b: Step forward one timestep (t -> t+dt) in time using\n",
    "        //           chosen RK-like MoL timestepping algorithm\n",
    "#include \"MoLtimestepping/RK_MoL.h\"\n",
    "    } // End main loop to progress forward in time.\n",
    "\n",
    "    // Step 4: Free all allocated memory\n",
    "#include \"MoLtimestepping/RK_Free_Memory.h\"\n",
    "    for(int i=0;i<3;i++) free(xx[i]);\n",
    "    return 0;\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:32:44.075262Z",
     "iopub.status.busy": "2021-03-07T17:32:44.072124Z",
     "iopub.status.idle": "2021-03-07T17:32:49.224157Z",
     "shell.execute_reply": "2021-03-07T17:32:49.223317Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Compiling executable...\n",
      "(EXEC): Executing `gcc -std=gnu99 -Ofast -fopenmp -march=native -funroll-loops ScalarWave_Ccodes/ScalarWave_Playground.c -o ScalarWave_Ccodes/output/ScalarWave_Playground -lm`...\n",
      "(BENCH): Finished executing in 1.0047273635864258 seconds.\n",
      "Finished compilation.\n",
      "(EXEC): Executing `taskset -c 1,3,5,7,9,11,13,15 ./ScalarWave_Playground 48`...\n",
      "(BENCH): Finished executing in 0.20269250869750977 seconds.\n",
      "(EXEC): Executing `taskset -c 1,3,5,7,9,11,13,15 ./ScalarWave_Playground 64`...\n",
      "(BENCH): Finished executing in 0.4042799472808838 seconds.\n",
      "(EXEC): Executing `taskset -c 1,3,5,7,9,11,13,15 ./ScalarWave_Playground 96`...\n",
      "(BENCH): Finished executing in 1.6105124950408936 seconds.\n"
     ]
    }
   ],
   "source": [
    "cmd.C_compile(os.path.join(Ccodesdir,\"ScalarWave_Playground.c\"),\n",
    "              os.path.join(outdir,\"ScalarWave_Playground\"))\n",
    "#!icc -align -qopenmp -xHost -O2 -qopt-report=5 -qopt-report-phase ipo -qopt-report-phase vec -vec-threshold1 -qopt-prefetch=4 ScalarWave/ScalarWave_Playground.c -o ScalarWave_Playground\n",
    "\n",
    "# 10o FD testing:\n",
    "# 4.46s\n",
    "# !icc -align -qopenmp -xHost -O2 -qopt-report=5 -qopt-report-phase ipo -qopt-report-phase vec -vec-threshold1 -qopt-prefetch=4 ScalarWave/ScalarWave_Playground.c -o ScalarWave_Playground\n",
    "# 4.65s\n",
    "# !gcc -Ofast -fopenmp -march=native ScalarWave/ScalarWave_Playground.c -fopt-info-vec-optimized-missed -o ScalarWave_Playground -lm 2>&1 |grep RHS\n",
    "# 5.45s\n",
    "# !clang -Ofast -fopenmp -mavx2 -mfma ScalarWave/ScalarWave_Playground.c -o ScalarWave_Playground -lm\n",
    "\n",
    "# Change to output directory\n",
    "os.chdir(outdir)\n",
    "# Clean up existing output files\n",
    "cmd.delete_existing_files(\"out??.txt\")\n",
    "cmd.Execute(\"ScalarWave_Playground\", \"48\", \"out48.txt\")\n",
    "cmd.Execute(\"ScalarWave_Playground\", \"64\", \"out64.txt\")\n",
    "cmd.Execute(\"ScalarWave_Playground\", \"96\", \"out96.txt\")\n",
    "# for benchmarking:\n",
    "# %timeit cmd.Execute(\"ScalarWave_Playground\", \"148\", \"out148.txt\", verbose=False)\n",
    "# FD functions enabled\n",
    "# 7.06 s ± 2.27 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
    "# disabled:\n",
    "# 7.06 s ± 2.58 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 3: Code Validation: Verify that relative error in numerical solution converges to zero at the expected order \\[Back to [top](#toc)\\]\n",
    "$$\\label{convergence}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:32:49.240042Z",
     "iopub.status.busy": "2021-03-07T17:32:49.238337Z",
     "iopub.status.idle": "2021-03-07T17:32:49.808041Z",
     "shell.execute_reply": "2021-03-07T17:32:49.808549Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEWCAYAAACXGLsWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAAsTAAALEwEAmpwYAABVtklEQVR4nO3dd3xT9frA8c+TNG1aWnYBGQWRpaCgIoiioOIeIE5ccN3zqtd1XfeKiuM68ecEFQfiRNwbRWUIgog4cLJxMISOpE2aPL8/zklJ96BtAjzv1yuvJjnJOU/Sk/Oc7zjfr6gqxhhjTFU8iQ7AGGNM8rNkYYwxplqWLIwxxlTLkoUxxphqWbIwxhhTLUsWxhhjqmXJIomIyAwROTvRcWwPROQ7ERmaBHF0EREVkZRG3OZTInJrY23PbBssWTQyEVkmIkERyReRP90fbmYt11HtAUZEbhKRsIjkubefRORBEdlhyz9F43G/r2FbuI5yB0dV7a2qM7YouOq3+6T7f+oW99wWf55kI45/isi3IlIgIqtE5GUR2TXRsZn6Y8kiMY5W1UxgD6A/cEMDbedFVc0CWgLHAu2ABVtbwqhKY56R14aIDAZ2SoI4vPW4rsq+6/HApcA/cfa1HsBrwJH1te0tlaz7yVZFVe3WiDdgGTAs7vFdwFvu/RnA2e59D04SWQ78BTwDNHOXrQAUyHdvgyrYzk3A5DLPeYFFwN1xzx0FfA1sBGYDu5WJ9SrgG6AAeAJoC7wL5AEfAS3iXn8M8J27rhnAzmXWdaW7rk3Ai4DfXdYaeMt93wbgc/fzPwtEgaD7Oa8Gurif/Sz3e/jMXcfLwB/uuj8DervPnwuEgZC7jjfL/h/c7+ol9zvOcz9D/7jY9wAWustedmO/tYr/cYr7+t3cWLu5z1f1eUa7n2cdcH01+9DO7ve70Y31mLhlTwGPAO+4/7NhwO7AV278LwIvxMdfg33gGvf/VgSklImlOxABBlQRbzP3u12Lsz/fAHjcZWOAmcDdwN/AUuBwd9lJwPwy67oceMO9n+a+bwXwJ/AokO4uGwqscmP/w/3u04Gn3e384H7/q+LW3R6Y6sa5FPhnmd9TVftIJ+BV973rgQfjlp3pbu9v4H2gc6KPQ3U6diU6gO3tRumDVCd3p7vFfTyDzcniTOAXoCuQ6e6Iz7rLuuAcYFKq2M5NlEkW7vM3A3Pd+7vjJKKBOIlktBtfWlysX+AkiA7ua79y3+cHPgb+6762B87B6WDA5/4QfwFS49Y1z/1BtnR/POe7y253f+g+97YfIGW/rzKf/RmgSdzB4UwgC+cAcj/wddx7nqLMwZ3yyaIQOML9Hm4HvnCXpeIc4C51YxuJk3iqShZXAePd+yXJoprPMxHnYNYX56C8cyXr9rnf63VubAfiHLx6xn3WTcC+OAm3qRv/5e57j8dJnrfWYh/4GmdfTa8gnvOB5dXs888Ar7v/ny7AT8BZ7rIxbjznuNu/AFgDCJDhfrbucev6EjjZvX8f8AbO/pQFvAnc7i4bChQDd7r7RDpwB/Ap0ALoiJMAV7mv9wALgP+432tX4Dfg0BrsI7GTsPtw9kk/MNhdNtz9f+2McxJxAzA70cehOh27Eh3A9nZzf3z5OGdxy4GH2XzAm8HmZDEduDDufT3dH1UKW5Yszgd+du8/gpuo4pb/CAyJi/XUuGVTgUfiHl8CvObevxF4KW6ZB1gNDI1b12lxy/8HPOrevxnnYNKtgniXUfHBtWsVn725+5pm7uOnqD5ZfBS3bBcg6N7f3/0cErd8Ztn1xS3r5B4cYtuuabLoGPfcPNwDYgXr3w/nTNkT99zzwE1xn/WZuGX74x58456bzeZkUZN94MwqvuvrcQ+alSz34iTXXeKeOw+Y4d4fA/wStyzD/T7auY8nA/9x73fHSR4ZOMmkANgp7r2DgKXu/aHudv1xy0sO/u7js9mcLAYCK8rEfi0wqQb7yCCcEkW53yNOKfysMr+LAFth6cLaLBJjhKo2V9XOqnqhqgYreE17nGQSsxwnUbTdwm13wKnqAegMXCEiG2M3nINd+7jX/xl3P1jB41jjfKl4VTUKrHS3F/NH3P1A3HvvwjnAfiAiv4nIv2vwOVbG7oiIV0TuEJFfRSQX5wAHTvVWTZWNze/Wc7cHVqv7Sy+77QrcD9ysqptqse2Ktp8J4HaEiN1y3HhWut9vzHJKf8/x8VUUf/x+VZN9oKrPux6oqg2sNU6Jpuy+XOF+oaoB925s35gCjHLvn4JzchIAsnGSxoK4uN9zn49Zq6qFcY/bl/ks8fc7A+3LfA/XUfr3Vtk+0gmndFVc/uPTGRgft84NOImuQwWvTWqWLJLXGpwdLSYHp1j9J86ZV62JiAc4GqdNAJwfyzg3ccVuGar6/JbGKyKC8yNaXd0bVTVPVa9Q1a447R7/EpGDYosre1vc/VNwivvDcOrHu8TCqGYdNfE70MH9PDGdqnj9QcBdIvKHiMQOLnNE5JS6xKKqmXG3FTjfcyf3fxmTQ+nvOX4bFcWfE3e/JvtAVTFPBzqKSP9Klq/DKRGX3Zer3S9cHwLZItIPJ2lMiVtvEKdtKhZ3M3U6jlQW9+841U8x8f/HlTilkvjvIUtVj6hBjCuBnEoa0VcC55VZb7qqzq7BepOKJYvk9TxwuYjs6HatvQ2nd1MxTpE3ilOvWi0RSRGRnd11tgPudRdNBM4XkYFu98cmInKkiGTVId6XgCNF5CAR8QFX4NS9V/ujEJGjRKSbe0DbhNNgGjtz/pPqP2eWu631OGebt5VZXpN1VGaOG8/F7vc4HBhQxet74LQ79HNv4CToafUQC8BcnLPaq0XE514rcjROo3VF5uCcZPzTff1ISse/RfuAqv6MU5X6vIgMFZFUEfGLyMki8m9VjeDsG+NEJEtEOgP/wqleqsn6wzidCu7CaZv40H0+6sZ+n4i0ARCRDiJyaBWrewm4VkRaiEgH4OK4ZfOAPBG5RkTS3dJqHxHZqwZhzsNJRHe4359fRPZ1lz3qbrO3G2MzETmhJp892ViySF5P4vTg+AynZ0YhThtBrKg+DpjlFm/3rmQdJ4lIPs4B+A2cg+meqrrGXc98nIbFB3F6avyCU4dca6r6I3Aa8H84Z31H43QRDtXg7d1xelbl4xzcHlbVT9xltwM3uJ/zykre/wxO1cZq4HucRvl4TwC7uOt4reafCtz4R+L0vtqI8xnfwklOFb3+L1X9I3Zzn14XV9VYk89TXTxHA4fjfM8PA2eo6pJq4h+DUwVyEk5nidjy+tgH/um+/yGc7+hXnK7ab7rLL8FpX/gNp71nCs7+XVNTcEqNL5ep6rnGjfcLt/rxI5y2vcrcjNNDaqn72ldw/49uUjsKJ8EvxfluH8cpqVbJfe/RQDecnlmrcL5nVHUaTiP7C26M3+L877Y6sR4nxpgaEpG5OI3zkxIdi6k7EbkApyPBkETHsjWwkoUx1RCRISLSzq2GGo1z/cR7iY7L1I6I7CAi+4qIR0R64lSVTqvufcZhVzUaU72eOPXdTXCqUo5X1d8TG5Kpg1TgMWBHnOqyF3Cq8UwNWDWUMcaYalk1lDHGmGptk9VQrVu31i5duiQ6DGOM2aosWLBgnapmV7Rsm0wWXbp0Yf78+YkOwxhjtioisryyZVYNZYwxplqWLIwxxlTLkoUxxphqWbIwxhhTLUsWxhhjqmXJwhhjTLUsWRhjjKnWVnOdhYhcgTM5e7aqrkt0PMaYbcOcOXOYMWMGQ4cOBSi5P2jQoFLLBg4cyMcff8xHH31E9+7dadWqFV9++SWLFi2iZ8+eHH744SxevJjXX3+d1q1b065dO9asWcPvv//OwIEDufHGG/nxxx954IEH6NatG7179yYtLY2mTZvSo0cPcnJyqg40wbaKsaFEpBPO2PK9cOZjqDJZ9O/fX+2iPGMMlE4G8Qlg//33Z+PGjRx33HGEw2FikwlGIhFEhJYtW7J+/fpGiXG33XZjxowZtGjRgiOPPJKsrCw6d+5M586d2X333Rk4cCAeT8NXBInIAlWtcNbDraVkcR9wNfB6ogMxxiSXypLB4MGD+fPPPznttNMIhUJ4vV5233135s+fT3UnyapaLlF07dqVpUuXoqqICJ06dWLlypWoKh6Phz333JMFCxYQjUbxeDzsuuuuLF68uORx9+7d+emnn0q2nZWVRV5eHgDffPMNLVu2ZIcddiA3NxePx0MgECASiQBw6aWXcv/995e8t/QsuY1EVZP6hjO38nj3/jKgdXXv2XPPPdUYs22YPXu23nbbbTp79uxyj2fNmqV+v189Ho/6fD4dNmyYejwexZl/u9zN7/eXepyZmakej0dFRFNTUzUtLU29Xq/6/X596KGHND09Xb1er6anp+tjjz1Wr4/PP/989Xq9CqjH49GBAweWPI7d2rVrp4cddphOmjRJo9GoLly4ULt27ar//ve/deHChRqNRuv1uwbmayXH1aSohhKRj3Dmhi7reuA64BBV3SQiy4D+WkE1lIicC5wLkJOTs+fy5ZUOcWKMSWJl2xAOOuggQqEQqampXHzxxdx///0UFxcjIni9XsLhcMl7PR4P0Wi05LHP56O42JmJNS0tjfHjx3PppZcSDodJTU1l+vTpADVqs6jvx2U/2+jRo5k4cSKRSASv18vxxx/P6tWrmT17NtFolK5duzJkyBB+/vln5syZQyQSoWfPnsyZM4cWLVrUy3dfVTVUUiSLyojIrsB0nAnqAToCa4ABunl+43KszWLbV9sfqUlelR1AU1JS6NevH3Pnzq3RelJTU3nggQe47LLLCIVCpKWllUsGybZvVJU87r///pLP4vF42H333VmwYAFer5fLLruMTp06sWjRIiZOnFhv8Wy1bRaquhhoE3tcVcnCbN1qc/CPRqMcfPDBFBUVkZqayhVXXMHdd99NOBzG6/UycuRIXn31VYqLi/F6vdxwww2kpKTwxRdf0KxZM3bccUf++OMPVqxYQffu3Rk9ejTLly/nk08+oVevXuy777789NNPLFy4kAMOOIBDDz2UL774ImkOMFu7ig6QRUVFeL1edtxxR4LBIOA0NMcnCq/XywUXXMDEiRNLksk999xDjx49+Oqrr0r+N7HG4vj/Vfz/bNCgQUnzPywby/Tp00tinzFjBqFQqKTdYsSIEUybNo0rr7yS//3vf/Tq1YsJEyY0WqxJXbIoq6bJwkoWyamyLoqw+YzK5/NxwQUX8OCDD1JcXIzH46FHjx4sWbKk2kbJhtK1a1eWLVtWUr3RqlUrMjMzUVW6dOnCs88+y7Jly3jxxRfp1q0bw4YNY4cddqBVq1aJaYhMMlUlh/bt2xNfZez3+wmFQiXfdfPmzdl9993x+/2cfvrpjBo1KqlKBg1pzpw5pUoa8aWktLQ0HnjgAZYvX865557LnXfeSfPmzbd4m1WVLBLegN0QN2vgTg7xDZGzZ88uadxLTU3V1NRU9Xg8mpKSol26dKm0QRLQjIyMUo9TUlJK7ns8Hh01apSmpaWVrO/YY49Vn89X0uh57LHHljR6ejwe7dOnj4pIyeNdd9215LGIaNu2bUttr0WLFlXGV9kttq6vvvpKp0+frieddJJeccUV+t577+kvv/yi4XA40f+iBlHV/71Hjx6lviOfz1fqfzlr1iydPXu2nnzyyfrCCy/UewPu1qay7zI9PV2nT5+uV1xxhXo8Hm3Xrp2+9NJLW7w9qmjgTviBvSFuliwSo6ode+TIkSUH5LK3Jk2aqNfrLVkefz8lJUVvuummUuuKrb+yHjLVxbKlvVjOPvvsUslnl112KfXZyia3qpLJHnvsUfL9vfPOOzp37lwNBAIJ+f/VVWXftd/v10GDBpX6zLHeRvHf1bBhw/T444/XqVOnJvqjJLXbbrut5Lvzer162223qarqggULdI899tDRo0dv8TYsWZgGUdUB+Ywzzqg0OZS9nXXWWTp79my95ZZbdOrUqRqJRKpNBlsS65Y+ri75xHeJ9Hq92rdv31LfRfzZNKCtW7fWffbZR1NTU0sSUO/evfX000+vl7PFhlST/3v84w8//FBnz56tY8aM0XvvvVcLCwsT/RG2GmW/6/jfQjgc1ry8vC3ehiULU+/K7rjxB8iyN6/Xq2PHjtW33npLmzdvrjk5OTpy5Eg95ZRT9N133030R6mT+komIqLdu3cvd4D1+/3q8/l04MCBumbNGg2FQjp06FC98cYbddasWQmtwor/rPFnu2VvHo9Hs7KyFNBWrVrpaaedpuvWrUtY3NuCLT1pqo4lC1MvKjtIeL1e7dOnT6kDxXHHHafHHXecZmZmalZWlubn56uqam5uboI/ReOoTTKp6OKssheWdejQQZs2bVqqHeXEE0/UL7/8slE/S9nYzzzzzFKJ7vTTT9fJkyfrpZdeqoAeeeSROn369O2+7aGh1HfysGRhtlhFZ8vxV8OmpKTogQceqCNHjtTDDz9c09LSFNDDDjtM3333XY1EIon+CEmlpsnD6/WWSx7Z2dm6++67a3Z2tn766aeqqvr999/rE088oUVFRfUeZ3UlyM6dO+u+++6re+21l44aNarkvStWrKjXWExpVVVL1ZUlC1MnVZUkBg8erF6vV30+n5511lm6fv16VVWdO3euZmRk6AUXXKA//PBDgj/B1qM2JY/43mM777yzTps2Ta+99loFtGPHjjp+/HgtKCiol7jK/t+POuqoUj3JbrrpJj3ggAMUnKEp7rjjjnrZrqleZQ3eW8KSham1yurd489wzzzzTF29erXOnj1bx44dW/LeTZs2JTDybUN1bSB+v19FpOTAPWjQIL3rrrt0v/32K2k0v/fee+t127Euyzk5OXreeefplClTtEuXLpqWlqYPPvigNVY3MitZWLJICmXPWq688krt3LmzAtq/f3/9+uuvVVV18uTJmpaWpjvttJP+/fffiQ16G1ZVKe/YY4/V9u3bl7QRTJo0SQ8//HC9/PLLVVU1Go3q2rVra7ydsgegxx57TFu2bKmAXnbZZSVde//880/de++9de7cuQ32uU3VrM3CkkVCVNVfvnXr1tqyZUt94403NBqNaiQS0euvv14BHTJkiPVyaUQVlTRuuukmveCCC7R58+YqIjpmzJiSzgQfffSRZmZmliT4qpRNRPvvv39J1dcnn3yikUhEJ02apKFQSFXVGq63MZYsthGbNm3Sn376SX/44Qf97rvvdMOGDaqqmpeXp3PnztXZs2fr8uXL67Tuis4oZ8+erRdddJE2bdpU27Vrp4sXLy55/emnn66Ann322fXeqGqqF0vsZbvlvvfee3rVVVepx+PRo446SsPhsP7000/avn177d69e4VVhJWdJMSSxjnnnKO5ubmal5enxx57rAI6efLkBHxq09AsWWxl/v77b501a5b+8ssvqur0dOnYsWO5fuxPPPGEqjo/9vjnu3TpomPGjNFvv/22xtusqLHs008/1aysLO3SpUtJLDEvvvii3nvvvXZmmWCVNXI+8sgjJck8Go3qp59+ql6vV0844YRS/7PKThLGjBmjgN5www2qqvrrr79qnz591OPx6H333Wf/922UJYskV1RUpE8//bQOHjy41FARnTp10o4dO2rXrl21WbNm2qpVK23durU2b95c27dvr4cccojecMMNOnHiRH3wwQf19ddf1/Hjx+vIkSO1VatW+tVXX6mqM4zE6NGjddKkSbpx48YKYyh70LjnnnvU7/drr169dOXKlaqqOm/ePJ0yZUqjfS+melU1ct5www0K6H/+8x9VVb3jjjsU0AcffLDkNRUlm9zcXM3JydGdd95ZCwsL9dNPP9WWLVtqixYt9MMPP2z0z2gaT1XJIqmHKN9WqSrfffcdU6ZM4ZtvvmHDhg3MmTOn1Gs8Hg+rVq1C1Zmycf/998fn87F+/Xo6dOhASkoK3333HR999FHJCJ1er5dhw4bxv//9jz59+pSMeLpy5Ureeustnn76aR544AFmz56N3+8vN3pnbHjkaDTKv//9b/r06cP7779PdnY2ixYtYujQoXTs2JHjjjuO1NTURv/eTHnx/7eyw7kffvjhrFmzhptvvpn27dtz1VVXMXv2bAKBQMn7hw4dSmpqasnIpkOHDuWGG25g5cqVzJw5k7S0NDIzM+nWrRvPPfcc3bp1S+CnNQlVWRbZmm/JWrJ4+OGHtVu3bpqenl5u7JyuXbuWGkjv4IMPLnXGd/7551dYXXDzzTfr008/rVOmTNFrrrlGW7RooR6PR8855xz9448/SrYdiUT0pZdeKqmaqOyMNBAIlAwAFyuFbNiwQbt27art27fXNWvWJOS7MzVT9v/62Wef6RFHHKEej0dfe+21Ci+OLNtmISJ60UUX6SeffFLyGqt22j5g1VCJMXv2bB03bpxOmDBBu3fvXm7wuPjkUDYZVDc4XUXJQ1V1/fr1eumll2pKSopmZWXp7bffrsFgsCSm6667Tj0ej1566aUV1nXH2j9ee+01VXWSzBFHHKE+n6/BxqMx9aeiaqX8/HwdMGCA+v1+nTVrlqqqvvfee3rmmWeWSgJFRUW6yy67aKdOnfScc85RQN95551EfRSTAJYsEuDdd98tNSxC2fuVlRTqMp5Q7KAQ//olS5bo0UcfXdLgPXPmTFVVLS4u1i+//LLSksV9992nQEkJ4s0331RAH3roocR8kaZWKvu//vXXX9q9e3dt1qyZ3nnnnSXtF/fcc0/Je2+66SYFdO+991ZAr7rqKhumZTtjyaKRzJ49W6+//nodOHBgudE3zzzzzGqTQ03WX9ORTWPr/Oijj7R9+/Y6ZMiQcuu7//779cYbbyy1/VGjRmnHjh1LvW7GjBlWDbEVqWy/Wrp0qR5yyCEKzhwiXbt2Va/Xq7NmzdJvv/1WfT6ftmrVSkVEH3jggQRFb7bEll6kZ8miEcycObPUDG4ZGRnq8/m2KDlUp6qreuPHibnsssvU7/eXuh5i9erVmpqaqsccc0yps8eddtpJjzvuOP31119rdBGX2brMnj27ZJDH2IlMs2bNdK+99tKsrCz1+/02CdFWqj6G/7Bk0UBiB+knnnhCc3JySv0Ax40b1+Bjz5eNpbId5ZVXXlFA58yZU+o948ePV6Bk8Le1a9cqoLfccov27dtXO3ToYOP9bGPiTyo8Ho926NChZL999tln9ffff090iKaO6mNgwaqShXWdraPYZOrBYBCArKwsfD4f0WiU1NRUDjjgAAYNGtRoE8pX1IUyZt999wVg5syZ7L333iXPX3LJJcyePZvrrruOAQMGUFhYCMCnn37KN998wzvvvENaWlqjxG8aR9musi+//DJLliyhsLCQU089taS7tdn6VNQNul5VlkW25ltjlCxiFzzFbrG6/8YqSdRWt27ddPjw4eWez8vL0169emmbNm306quvLumhdfPNNzd+kKZRJPN+arZMQ7ZZiLN829K/f3+dP39+va83drHTbrvtxpgxY1i3bh0Afr+fjz/+uNFKEXXxj3/8g1dffZVrrrmmpNQT8/333zNjxgymTJnCrFmzOOqoo3j99dfxeDwJjNgY09hEZIGq9q9wmSWLmolVOxUVFZVcMd2nTx9GjhzJYYcdlpSJIv4K7TfffJPbb78dj8dDWloa06dPLxWzqpKdnU3nzp2ZPn06zZs3T1zgpl6VvVLfmMpUlSyszaKGZsyYQSgUKkkUe+21F1988UXSnn3Hklus/nL48OEARKNRQqEQM2bMKHXgWLp0KevXr2fcuHGWKLYhZfeDsicJxtRUch7pyhCRS0RkiYh8JyL/S0QMPXv2LEkUaWlpjB8/PmkTBWxObpFIhFAoRLNmzQAQkQobv+bOnQvAwIEDGztU04DK7gczZsxIdEhmK5W8RzuXiBwADAf6qmpv4O7GjmH9+vVcfvnlpKens88++/DWW28l/dlZrGeE1+slNTWV0aNHs//++9O8efMKzy7nzZtHeno6vXv3TlDEpiGU3Q/qvYeM2W5sDdVQFwB3qGoRgKr+1dgBzJkzhxUrVnDBBRfwyCOP0Ldv38YOodYq6ko7YsQIPvvsMzp37lzu9fPmzWOPPfbA5/MlIFrTUKrqUm1MbWwNyaIHsJ+IjAMKgStV9cuyLxKRc4FzAXJycuplw7GGwUgkAji9ngAyMzPrZf0Nrex1HoMHDwac6y1OPPHEkufD4TBfffUVF154YaPHaBpeY17vY7ZdSZEsROQjoF0Fi67HibElsDewF/CSiHTVMt24VHUCMAGc3lBbGlN8w6DX6wUgFArh8XhKksbWpl+/fmRkZJRLFosXL6awsJABAwYkMDpjTDJLimShqsMqWyYiFwCvuslhnohEgdbA2oaMKb5hMJaXQqEQmZmZW+1Vrj6fj7333puZM2eWej7WuG3JwhhTmaRv4AZeAw4AEJEeQCqwrqE3Gt8wGCtZFBYWkpWV1dCbblD77bcfixYtIjc3t+S5efPmkZ2dTZcuXRIXmDEmqSVFyaIaTwJPisi3QAgYXbYKqiHENwz+9ttvPP7449x///0UFRU19KYb1ODBg4lGo3zxxRcccsghgJMsBgwYsNWWmIwxDS/pSxaqGlLV01S1j6ruoaofN9a2Bw0axLXXXkvr1q3x+Xy0bNmSHXbYobE23yAGDhyI1+vl888/ByA3N5cffvjBqqCMMVVK+mSRDILBIBkZGTz22GM8++yziQ5ni2RlZdGvX7+Sdov58+ejqnYxnjGmSpYsaiAQCJCens7EiRN58cUXEx3OFttvv/2YO3cuoVCIefPmAc7wJcYYUxlLFjUQCATIyMggLy9vq7nGoiqDBw8mGAyycOFC5s6dS7du3WjZsmWiwzLGJDFLFjUQSxb5+fnbRLKInwxp3rx5VgVljKmWJYsaiLVZbCvJol27dnTr1o0XXniBNWvWWOO2MaZalixqINZmkZ+fv9VfZxEzePBgYnN+WLIwxlRna7jOIuECgQDZ2dkEAgG2lcmi9ttvP5566il8Ph/9+vVLdDjGmCRnyaIGYm0WaWlpiQ6l3sQGFezbt+9WO9aVMabxWDVUDQSDQUSESy65hK+++irR4dSL7t27s9NOOzFsWKXDchljTIlalyxEpAlQqKqRBognKcWqnx588EEOOOAA9thjj0SHtMVEhG+++YbU1NREh2KM2QpUW7IQEY+InCIib4vIX8AS4HcR+V5E7hKRbg0fZmIFAoGSwQS3hd5QMRkZGaSkWE2kMaZ6NamG+gTYCbgWaKeqnVS1DTAY+AK4U0ROa8AYEy4YDJbMt70tJQtjjKmpmpxWDgOKgY6qGo09qaobgKnAVBHZZufiDIfDFBcXlzy2ZGGM2R5VW7JQ1bA7JPg7Vb2mXqNKIoFAAKBkXgtLFsaY7VFtekN9JSLb3WhzsWQxePBgwuEwO+64Y4IjMsaYxleb1s2BwKkishwoAARQVd2tQSJLEsFgEHAag21yIGPM9qo2yeLQBosiicVKFgsXLmTmzJlMmDAhwREZY0zjq3E1lKouB5oDR7u35u5z27RYsvjtt9946aWXEhyNMcYkRo2ThYhcCjwHtHFvk0XkkoYKLFnEkkU4HLbGbWPMdqs21VBnAQNVtQBARO4E5gD/1xCBJYtYm4UlC2PM9qw2vaEEiB/iI+I+t02LlSwKCwstWRhjtlu1KVlMAuaKyDT38QjgiXqPKMnEkkVGRgZNmzZNcDTGGJMYNUoW4vQZfRmYgTPMB8A/VHVhA8WVNGLJYtKkSbRv3z7B0RhjTGLUKFmoqorIO6q6K7BtjNFdQ/HXWRhjzPYq6a/gFpF+IvKFiHwtIvNFpFHnAI2VLC666CK7xsIYs93aGq7g/h8wVlXfFZEj3MdDG3ibJQKBAB6PhzfeeIM2bdo01maNMSap1KbN4lwgERfhKRBrWW4GrGnMjcemVC0oKLDeUMaY7VZt2iwectssGttlwPsicjdOtdk+Fb1IRM7FSWjk5OTU28aDwSDp6enk5+eTlZVVb+s1xpitSVK0WYjIRyLybQW34cAFwOWq2gm4nEq666rqBFXtr6r9s7Oz6y22QCCA3+8HbC4LY8z2q7ZtFqeJyDLquc1CVYdVtkxEngEudR++DDy+pdurjViy6N69O/WZhIwxZmuyNYw6uwYYgnONx4HAz4258UAgQNOmTZk/f35jbtYYY5JKbZLFCuBUoKuq3iwiOUA7Gr7R+xxgvIikAIW47RKNJRgM2jUWxpjtXm3aLB4GBgGj3Md5wEP1HlEZqjpTVfdU1b6qOlBVFzT0NuMFAgFCoRAHHnggS5YsacxNG2NM0qhNshioqhfhnN2jqn8DqQ0SVRIJBAKoKp988gnh8DY71bgxxlSpNskiLCJenOseEJFsINogUSWR2EV5YL2hjDHbr9okiweAaUAbERkHzARua5CokkgwGCyZe9uShTFme1XjBm5VfU5EFgAH4XSbHaGqPzRYZEkiNjYUWLIwxmy/atMbClVdAmxXrbyxrrP9+vUruTjPGGO2N7WphtruRCIRQqEQ++yzDwsXLiypjjLGmO2NJYsq2FwWxhjjsGRRhVh7xbvvvsuJJ56Y4GiMMSZxatxm4Q5TXu4KblWd12DRJVgsWaxdu5ZNmzYlOBpjjEmcpL+CO5FiySIcDltPKGPMdq1Wo86q6h4ishCcK7hFZJu+gjvWZhEKhSxZGGO2a3YFdxViJYuioiKb+MgYs12zK7irEEsWO++8M3379k1wNMYYkzh2BXcVYsnirrvuol+/fokNxhhjEqg2vaH+Bbyoqtt0o3Y8u87CGGMctamGygI+EJHPReRiEWnbUEEli1jJ4qCDDuKJJyqc+tsYY7YLNU4WqjpWVXsDFwE7AJ+KyEcNFlkSiCWLVatWUVRUlOBojDEmcepyBfdfwB/AeqBN/YaTXGzEWWOMcdQ4WYjIhSIyA/gYaAWco6q7NVRgySDWZgGWLIwx27faXJTXEbhUVRc1VDDJJhAIkJaWRlFRkSULY8x2rdqShYjMdO/+E/hcRHLdW56I5DZseIkVCARIT09nxIgRdOjQIdHhGGNMwlRbslDVwe7f7e7UOhAIkJWVxbRp0xIdijHGJFRt2izurMlz25JgMGjXWBhjDLXrDXVwBc8dXh9BiMgJIvKdiERFpH+ZZdeKyC8i8qOIHFof26upQCBAIBCgTZs2rF+/vjE3bYwxSaXaaigRuQC4EOgqIt/ELcoCZtVTHN8CI4HHymx7F+BkoDfQHvhIRHqoaqSetlulQCCAx+Nh7dq1VsIwxmzXatIbagrwLnA78O+45/NUdUN9BBEbY6qCOa6HAy+oahGwVER+AQYAc+pju9UJBAKICF6vF7/f3xibNMaYpFSTBu5NwCY2T3rUmDoAX8Q9XuU+V46InAucC5CTk1MvGw8Gg4gImZmZFSUyY4zZbtTmOgtEpAXQHSg5zVbVz2r43o+AdhUsul5VX69NHBVR1QnABID+/fvrlq4PNl/BbddYGGO2d7UZdfZs4FKci/O+BvbGqQ46sCbvV9VhdYhvNdAp7nFH97lGEQgEyMnJYdCgQY21SWOMSUq16Q11KbAXsFxVDwB2BzY2RFBx3gBOFpE0EdkRp1Qzr4G3WSIQCLD77rtz9913N9YmjTEmKdUmWRSqaiGAiKSp6hKgZ30EISLHisgqYBDwtoi8D6Cq3wEvAd8D7wEXNVZPKHDaLNLT0xtrc8YYk7Rq02axSkSaA6/hdGHdACyvjyBUdRrOlK0VLRsHjKuP7dRGNBolGAwyefJkfvvtN6ZOndrYIRhjTNKozbSqx7p3bxKRT4CmOGf726TCwkIAwuEwqampCY7GGGMSqyYX5eUBsd5FUsH9pg0TWmLFekKFw2HrDWWM2e7V5DqLrMYIJNnE5rIoKioiK2u7/AqMMaZEbQYSFBE5TURudB93EpEBDRdaYsVKFjaXhTHG1K431MM4vZVOcR/nAw/Ve0RJIpYsjj76aPbZZ58ER2OMMYlVm95QA1V1DxFZCKCqf4vINtvyG0sWF198MYccckiCozHGmMSqTckiLCJe3AZuEckGog0SVRKIn39btV5GDzHGmK1WbZLFAzjXQrQRkXHATOC2BokqCcRKFoceeqjNlGeM2e7V5jqL50RkAXAQTrfZEbGhxbdFsWQBWG8oY8x2r1ajzrpDfCwBEBGPiJyqqs81SGQJFp8srDeUMWZ7V201lIg0dac2fVBEDnG70F4M/Aac2PAhJkZ8m4UlC2PM9q4mJYtngb9xhiM/G7iOzdVQXzdcaIll1VDGGLNZTZJFV1XdFUBEHgd+B3JiI9Buq2LJ4tprr6VVq1YJjsYYYxKrJskiHLujqhERWbWtJwpwkkVaWhq33bbNdvgyxpgaq0my6Csiue59AdLdxwKoqm6TAwkGg0H8fj9r164lOzs70eEYY0xC1WQgQW9jBJJsAoEAxcXFdO/enY0bNyY6HGOMSaia9IaS+njN1iYQCOD1eq0nlDHGULMruD8RkUtEJCf+SRFJFZEDReRpYHTDhJc4gUAAEbGeUMYYQ83aLA4DzgSeF5EdgY2AH/ACHwD3q+rCBoswQYLBICJiJQtjjKFmbRaFOMOTPywiPqA1EFTVjQ0cW0IFAgFU1ZKFMcZQ++E+wjjXWWzzAoEAXbt25cILL0x0KMYYk3C1GXV2uxIIBOjWrRsnnHBCokMxxpiEs2RRiWAwSCgUYu3atYkOxRhjEq7WyUJEThGRF0TkORGZIiKjGiKwRAsEArz99tvcfvvtiQ7FGGMSri4liyGqerKqnqqqpwCDtzQIETlBRL4TkaiI9I97/mARWSAii92/B27ptmoqdlGeNXAbY0wtG7hdaSJyJLAS6Aik10Mc3wIjgcfKPL8OOFpV14hIH+B9oEM9bK9KqloykKBdZ2GMMXVLFhfiHNh3BVYBF21pELEZ98peCF7m+o3vcMalSlPVoi3dZlVCoVDJvNtWsjDGmDokC1UNAJNjj0XkGuDO+gyqEscBX1WWKETkXOBcgJycnIpeUmPb2yx5G9avZ/WynwlpCmx7I7cYY1RJlWI6dOlOyzpOuVDrZCEiL8U/BPpRg2QhIh8B7SpYdL2qvl7Ne3u72ziksteo6gRgAkD//v21uniqEksWZ5xxBoMGDdqSVSW9DevXs/LX79lpzjVkbPwRjxYnOiRjTD2LSgqB5j35qWgcv0Wj7LHPAXg8tWuyrks1VK6qnh17ICKP1ORNqjqsDttCRDoC04AzVPXXuqyjtmJTqh588MF069atMTaZMKuX/cxOc64h8+/vEh2KMaaBeLSYzL+/o8eX1/Ndv5v49M3nGXL0qFoljLr0hhpX5vH1dVhHjYhIc+Bt4N+qOquhtlNWrGSxcuXKUlVS26KQppCx8cdEh2GMaQQZG38kpUkLvp49nT9XLa3Ve2udLFR1aZnHG2q7jrJE5FgRWQUMAt4WkffdRRcD3YD/iMjX7q3Nlm6vOrEEcd111/Hbb7819OYSS8SqnozZTni0GMSDx+Mhf9PftXpvjauhRORfFTy9CVigql/XaqtlqOo0nKqmss/fCty6Jeuui+2tgdsYs51RiEYjtXpLbUoW/YHzca5z6ACchzN8+UQRubpWW01ysTYLsGRhjDFQu2TREdhDVa9Q1SuAPYE2wP7AmAaILWGsZLF1GPNaEBmby9UfFpZ6flVuFBmby4xlyVO9FggrvR/OR8bmMnNF6bgW/h7h0MkFtLwzl6a357LvkwV8+GvyxG4M1C5ZtAHir3EIA21VNVjm+a1eLFl4vV7S0tISHI2pij8FHpgbYvnGaKJDqdKFbxeyU4vyP7dAWDlkcoBmacLn/2jC/HOb0K+th6OfD7AsyT+T2b7UJlk8B8wVkf+KyE3ALGCKiDQBvm+I4BIlliweeeSRcleVm+SyTycvfdt5uO7jwkpf89J3YVJvyWXe6s11tM8sCpE+Lpdv/qxdvW1dPP11iK//iHDXweVPPH5aH2VdQLlx/zR6t/HSo5WXO4b5KYrAoj8aPjZjaqrGDdyqeouIvAvs6z51vqrOd++fWu+RJVCszeKkk05KcCSJM/SpgnLPndjbx4V7pRIIK0c8V75L8Zh+Psb0S2VdIMrxLwXLLb+gfyon9fGxclOU06eVXz5jTJNaxynA3Qf7GfJUgMv3jtC/vbfCuD/8tZhRUwMsPC+TP/KjXPROIfcc4me3tuVfH3P4cwV8vrzqA/a7p2awX+fKf0Y/rI1w1YdFfPaPDNJSyp94dG/poW0TYdLXYcYd6CHFA4/OD9EyXRjUqfLYjGlstb0oLwxEAXXvb5NiJYvFixez7777VvNqk2j7dU5heK8UrvygsNKEM/5wP3tNLODsN4L8vCHKsK4pXLhXapXrffzodILVNB10yKq85BkIKye8HOSOYWn0au2tsFqpSarw2T8yOO6lIPd/EcIj0KaJ8P5pGbRpYtPNmORRm66zlwLnAFNxTugmi8gEVf2/hgouUQKBACLCZZddxpdffpnocBKiqrP8DJ9Uubx1hqfK5Z2aVb28Lu4clkbvhwt448cwe+xQ/ow8wye8eHw6/R4toG2mMP2M6gdL7tB0yw7W/3y3kF3bejhz98qTUjCsnPl6ITu39jDxaD8+jzBhQYijnw8w9+wm5DSzhGGSQ21KFmcBA1W1AEBE7gTmANtksvB4PDY8+VakRysv5+3p45qPinj31IwKXzNzhVOltKlQWVsQpWV61dU8W1oN9dFvxazMVV7+LrfU80OfCnBQVy/vn9aE578N882fET4ZnYXP65RSHmufzvQH8piwIMStB/qr3L4xjaU2yUKA+F9OxH1umxMMBhER6za7lfnvkDSe/SafCQtC5ZZ9+1eEf71fyOPH+HltSTEnTw3yxVlNKmxHiNnSaqgPTs8gFPeLWZOnHDo5wKTh/pIEUxByBvr1lFmN1yPoFg2HaUz9qk2ymITTG2oaTpIYATzZEEElWqzNwpLF1iW7iYd/75vGLZ+V7sldWKyMmhpkRK8UxvRL5ZiePvo+ms/VHxYx/vDKz9y3tBqqR6vSJZfMVKfNYscWHro0d9Z9aDcvV38EZ71RyFX7pJLigccWhPl1Q5RjetZlnE9jGkaNfw2qei/wD2A9zgx2o1X1voYKLJEsWWy9Lh+USuuM0qfpl79XSEFIefQop52iZbowZWQ6D88P8fZPie2n0aOVl3dPzWDZxiiDJxUw4PECZq+MMO2kdAZ2tGRhkke1e6OI5OH0fip5Km6ZqmrThggskQKBAN26dePSSy9NdCimCk+NKN9I7U8RVlxeuq3pkaPKv26/zimEb2zcXbdLcw/63/LbHNolhRljLDGY5FbtHqqq210rbzAYpG3btvTu3TvRoRhjTFKwfnkVKCgoYOPGjfz6a6PMtWSMMUnPkkUF8vLyWLRoER988EGiQzHGmKRgyaICBQXOUBfWwG2MMQ5LFhUoLHQGpbOL8owxxmHJogKxgQStZGGMMQ5LFhWIlSwsWRhjjMOSRRnhcJhIJMJ5551Hnz59Eh2OMcYkBUsWZcSqoHr06GElC2OMcVmyKCM21MfChQsJhcoPSGfM9qjL/Xnc+lntZ08e81qQYc+Un0irvhzzfIC7Z29dszo/vzjMXhPz0SpGinznZ2d2x7Z35/HjuuSYMdGSRRmxZDF58uQq/5km8ca8FkTG5nL1h6WnVF2VG0XG5jJjWTVDxjaSiQtC9H00H/+tubS8M5ejny8/yyBAVJWDnilAxuYy+ZvkPlHp9kAeN82ofCrbxjD9t2LmrY5w8YDS84Us3xjl1FcDtP5fHv5bc+n5YD5v/ljxGGAfLy3Ge3Mu3R7Iq3D5pIUhdn8sv1bve/jLELs8lE/GuFx2uCeP0a8F+TN/88RXJ/dJIRCG5xZXHNMnS4s57qUgVwxKZUAHL8OerX4+9qe+DiFjcxs0MSdFshCRE0TkOxGJikj/CpbniEi+iFzZ0LHEkoXH4yE1teqZ1Ezi+VPggbkhllfzY0qUGz8u5IZPivjX3qksvqAJM89swqm7+ip87c2fFtHEt02O+t8g7v0ixBl9ffjjhplfnRtl7ycKUIW3TklnycWZTDjKT8cKRhD+Iz/K6NeCHLJT5fOaTFtSzMhepf9fVb3v5e/CXPpeIf8alMr3F2Xy8gnpLFgT4YzXNk8jLCKctbuP+78of0Iwd1Uxw18IcNuBadw+zM+0k9I5oEsKBz1TwJq8ivfx79dGuG56Eft3bthpeJMiWQDfAiOBzypZfi/wbmMEEmuzSE9PR8R+uMlun05e+rbzcN3HlZ/lvvSdU6Sft3pzcf6ZRSHSx+XyzZ8NV8T/dUOU22aGeGZEOqP7pdK9lZddsr2c3Kd8svh4aTFPLgwzaXjdJzuSsbn839wQJ70SoMltueTcl8cr34fZVKic+mqArNtz6To+j6nfbz6jXbbRKYXNXFG6FFZVyWHoUwX8+rcy9lPnbFbG5lZ75nvfnCI63JtHxrhcTng5wIagU2qfscw5Q1+5qfT7n1kUotkduRSEKi7drw9Eee+XYkb0Kj283XUfF9GluYcpx2Wwd8cUujT3MKRLCruXmT0xqspprwa5aK9UBnao+CCbH1I+/K2YY3dOqfH7Zq2MsFtbD2fvkUqX5h4G56Rw3p6ppfY9gGN7+Vjwe5QlcVVMi/6IcPTzQR48ws/lg9IASPEIT4/wM7ynj4OfDbAuUPp7CoSVE18Ocu+hfnZs3rCH86RIFqr6g6r+WNEyERkBLAW+a4xYYiWL9PTqp900iSfA3Qf7eX5xMfPXVHzgP7G3j9F9fYyaGiC3SPlpfYSL3inknkP87Na28rOxw58rIPO23Cpvny+vvKpr2pIwPg+sCyi9H86n/T15HDklwLd/lY7zz/woZ0wL8syx6bTK2LKf5LjPiziiWwqLzs/kqB4pnD4tyMlTAxzcNYWF52VyZPcUzngtyPpA3Utir56UQZfmwhWDUvn9ikx+vyKTTk0rP7GatzrCJ8sivHdqBu+cmsHXf0Q56w3npGxolxS6t/Tw5MLSVTITvwpzSh8fTVIrXu/MFREESk2hG1XltSVhBnX0MmpqgDZ35bHrI/nc/nkRxdHSSeeWT0OIwDX7Vl578O7PxXRq6qFPG2+N3zc4x8sPa6PMWFaMqvJHfpRXfghzZPfSJwg7tvDQponwydLN+0Lfdl7+uiqLM/qWXreIcO+hfr67MJPWZfaPi94pZGCHik9A6ltSj4ssIpnANcDBQJVVUCJyLnAuQE5OTp23GUsWTZrU7xzRW5PL3ivk6z8av1GtXzsv9x9W+zPr/TqnMLxXCld+UFjp3N7jD/ez18QCzn4jyM8bogzrmsKFe1VdzbilM+X9uiFKVOG/MwoZf5ifNk083DW7iP0nFbDk4kzaNPEQVeXUV4OcubuPoV22/Od4ch8fo/s5n2vs0DQemR+mWwsPY9znbj7Az4NfhpmzKsJRPeqWmFqmC16BzFRol1n9OqIKzx6bTjO/8109dISfQycH+GVDlG4tPZy7p4/xc0PcOCQVjwhL1kWYuSLCA1XsC0s3RmmVIaWqoNYWKLlFTpvBRXul8v5paXy/Nsol7wbJDynjDnLW98nSYh5dEGLheU2qrD2YtiTMsXEll5q87/hdfGwqVI54LkA4CsVROLJ7Ck8cU/6zdMgSfvu77kn7mUUhvlgVYf45jXOsarSShYh8JCLfVnAbXsXbbgLuU9XyLUxlqOoEVe2vqv2zs7PrHGcsWTz00EN1XodpfHcOS2PWyghvVNKQmeETXjw+nVd/KOavAuWJY6ovOXZo6qFby6pv6VW0MUQVwlEYf5ifI3v42KuDl2eOTccjwuRvnDhv+zxEUcSZErY+9G27+Sed3cSDVyhVemqRLqR64a+Cxuu8sUu2pyRRAOzbyYnn+7XOCcnovj7+KlDe/8XJzI9/FWbPHTzlqo7iBcNOe1W8WOFh17Ye7jnUz+47eDl1Nx/X7ZfGg1867QPrAlFOmxZk0vD0KhNdKKK8/fPmKqiavu/z5cVc93ERdx3sZ8G5TXjnlHSWbYxy5hvBcq/1pwjB4rr9H35cF+Hy94t44bj0Sktf9a3RShaqOqwObxsIHC8i/wOaA1ERKVTVB+s1uDixNouePXs21CaSXl3O7hOtRysv5+3p45qPinj31IwKXzNzhXNw2lSorC2I0jK96gbBw58r4PPlVZew3j01o2Q+7bJ2cEsdveOqMfwpwk4theUbnYPER78VM3tlhLRbS/eqGf1aIbd+FmLJxbW71sdXwUcq+5yw+cAam/u7bMe/cCP2F2iV4eH4XXxM/CrMQV1TeGZRmFsPrDp5ZjeRknaPmNYZgs8DvbNLf+De2R5yi+DvoPLtX1HW5ClHTdncIy2qzuxuKTfn8syx6Zyyq4/pvxWTmSol7RI1fd/1HxcxslcKF7k9tHZr6yUzVdj/qQBjhzolqZgNQSW7jtWOc1ZF2BBU9pywufdT7H+acnMun47JYN+c+j28J3U1lKruF7svIjcB+Q2ZKGBzyWLWrFl069atITdl6tl/h6Tx7Df5TFhQvpfJt39F+Nf7hTx+jJ/XlhRz8tQgX5zVhLSUys/KtrQaar+cFCDEknXRkjm3QxFl6d/Kyb2d900ank5BuPRBb9dHChh3YBrH7dzwP89sdwraNXmbY/irIMrq3KrPeFO9QqSGCeWHdVFyi5Smac62Zq90EvAucQf18/b0ccDTAR6bHyZYrIyqpg5+jx285IdgxaYoOc2c79bnFQZ29LJkXenAflwfpVmaU6raq72XxReUrrZ5+MsQb/1UzDunZtDJ7TX16g/FjOiZUlLdVNP3FYS1JAHHeN18EN8VPxBWfv07Sv/2dUsWI3r56N++dFK84eMi/ixQJh7tp2uL+q80SopkISLHAv8HZANvi8jXqnpoImKJJYu33nqL0aNHJyIEU0fZTTz8e980bilz8VhhsTJqapARvVIY0y+VY3r66PtoPld/WMT4wysvRXWooLtlbRy4o5e9O3q57L1CJhztp00T4Y6ZIaIKp+3mHAx3rORH3bGp0L1Vw3aFBEj3Cft28vK/2UX0au2hOArXf1xIWjVHhh1beJi1MsKKTVEyfE47hqeSenwBzpgW5NYD09gQVC56p5BjeqaUOssenJNCz1YervywkDN285GVVnXVSr92HnbIFD5dVszpcQ3C1w5O5agpQf77SSGn7ebjh3VRxn0e4tKBzmuapEqpBmuANk2cqrnY81FV3vipmCkjN1dV1uR9ACN6+rh9ZhEDOnjZv3MKq3KjXPZ+Ibu19bBT3OedtSJCmheG1LGdqrlfaO73lnsuP6Tl4qwvydIbapqqdlTVNFVtW1GiUNWbVPXuho4lliyaNWvW0JsyDeDyQam0zih9oLn8vUIKQsqj7lzcLdOFKSPTeXh+iLd/qriNoz6ICG+cnM5eHbwc83yAQU8U8Ed+lE/HZJDdpHY/vaFPFTD0qYa54OrJ4X4yU4V9nizg5KlBzt0zlR0yqz5Yjx2axsZCpeeD+WTflc+KTZWXRAZ08DI4x8vBzwY4bHKAXdt6eLKCBt9z9vARisC5e1Z/fZNHhPP2TOXZb0r//47o7uP549J55Ydidn2kgCs/KOLKQancWIs2oVkrIoQjypAutT/oXrdfKjfsn8ZtM0Ps8nA+J70SpFdrD2+OyiiVTCcvDnPqrj4yG6m9oT7ItniVcv/+/XX+/Pl1eu+NN97IrbfeyuWXX869995bz5ElnwULFrDnmwcmOgxTjZz78rigfyrX7lc/DeHJ6OoPC/nwt2IWnlezdpq/g06yev+0jCobw2vr8vcK2VCoPD2iYbrPr9wUZbdH8/n6vEw6N/C1ERVZcPTHzHz5EQ4bdQ49+w4stUxEFqhquQujIUmqoZJJbJY8m/jIJItv/ozgTxGu2GfbHFFgU6Hy0/ooExaEeKCKasGyWqQLk0emsyYvWq/JYudsT7n2gPq0bGOUiUenJyRRbAlLFmXk5uYCNpeFSR67tfXy0yXb7v44/IUAc1dHOLmPr6Qtp6YO2an+D2E1qQbbEpX1nkt2W2fUDSgUCtGpUyfOP//8RIdizHahsgspTXLZuspBjaCwsJDMzEyrhjLGmDiWLMpYunQpK1eu5MUXX0x0KMYYkzQsWcSZM2cOCxYsID8/n9GjRzNnzpxEh2SMMUnBkkWcGTNmlFxlGQ6HmTFjRmIDMsaYJGHJIs7QoUNLLu9PTU1l6NChiQ3IGGOShCWLOIMGDaJ58+YAPPfccwwaNCixARljTJKwZFFGNOoMQjZkyJAER2KMMcnDkkUZIsL5559Py5YtEx2KMUmjy/153FpmgMaaGPNakGHPNMyYVgDHPB/g7tm1jyuRnl8cZq+J+VQ11NI7PztTAbe9O48f1zX+RGQVsWRRRjAYpFmzZjb/9lZgzGtBZGwuV39Yeq7oVbnOvNIzllUzvngjeG1JmAET88m8zfnh//PdQoLh8geJ15aEGfh4Punjcml2Ry77Tyogv5L5p5NBVXN0N5bpvxUzb3WEiwdsvuJ6XSDKeW8G6XRfHunjchkwMb/C/WBdIMoFbwVpf08eabfmsuP4PCZWMLT9pIUhdn+s/NxrHy915g7v9kBeuWUPfxlil4fyyRiXyw735DH6tSB/5m8eNv3kPikEwvDc4ooHsfxkaTHHvRTkikGpDOjgZdizgWrnOH/qa2c+9IZMzJYs4kQiEYqKiqzL7FbEnwIPzA2xvJofUyJ88Kvzoz+5j49F52fywnHpvP9rMWe/WXrWtCe+CnHGtCCn7urjq3ObMO/sJlwyIBWvna9U6d4vQpzR11cytaqqcuyLQRb8HuHF49P55vwmHLJTCodNLj3veX5I2X9SgF/+jvL8cen8eHEmU0ams3N2+cPhtCXFjOxVegiSP/KjjH4tyCE7lR8/6uXvwlz6XiH/GpTK9xdl8vIJ6SxYE+GM1zb/z0WEs3b3cf8X5ZPT3FXFDH8hwG0HpnH7MD/TTkrngC4pHPRMAWvyKt7Hv18b4brpRezfuWGHtLdkESc2S96iRYsSHMnWac7KYm7/vIg5KxvvjH6fTl76tvNw3ceVn+W+9J1TpJ+3evMB45lFIdLH5fLNnw1XxH9mUZiDu3r516A0dmrp4YAdU7jjoDSmLC5mqTv3cm6Rcvn7hdx1sJ9/Dkxj52wvPVt7OaG3r8opWysiY3P5v7khTnolQJPbcsm5L49Xvg+zqVA59dUAWbfn0nV8HlO/33xGu2yjUwqbuaL0/6yqksPQpwr49W9l7KfO2ayMza32zPe+OUV0uDePjHG5nPByoGSWuxnLnDP0lZtKv/+ZRSGa3ZFLQSWlq/WBKO/9UsyIuDmyf/1bmbkiwoNH+NmnUwrdW3m59UA/PVt7uGv25gPzXbOKCISVt0ZlMKRLCl2aexjUKYXBZWaWyw8pH/62eWpVcOa6OO3VIBftlVoyi168WSsj7NbWw9l7pNKluYfBOSmct2dqqX0P4NhePhb8HmVJXBXToj8iHP18kAeP8HP5IGd04RSP8PQIP8N7+jj42QDrAqW/p0BYOfHlIPce6mfHBh6Y0MaGihOby8Lv3/qmFU20OSuLOeiZAKEIpHph+hkZDOrU8LuXAHcf7GfIUwEu3ztS4WihJ/b28eGvxYyaGmDheZn8kR/loncKuecQf6n5qcva0mlVC4u15Kw3JpYAPltezI4tUvng12LyQpDug/4T8lmZq+zc2sMtB6TVacC5cZ8XceewNMYd6OfeOUWcPi3I0C5eTurtY+xQP+O/KOKM15znWtVxSs9XT8pgzwn5HLezjyvdkXCzMypPbPNWR8jwCe+dmsH6oHLOm4Wc9UaQaSdlMLRLCt1benhyYZj/Dt08/PrEr8Kc0sdX6fzSM1dEEJwZ82IK3fmsy33nKcJnyzcnw6k/FDM4J4XL3y9k2pJimqUJR/dIYewBaWTEJeh3fy6mU1NPqcmEbvk0hAhcs28qYz8t31YyOMfLhAUhZiwrZkhnL38WKK/8EObI7qVLJzu28NCmifDJ0gi9Wjvr79vOy19XlR9mSES499CKj0kXvVPIwA5eTu7j471fGvYkzZJFnFjJIj29Ycax35bNWBYhFIGIQijiPG6MZAHOKJ7De6Vw5QeFlQ5KN/5wP3tNLODsN4L8vCHKsK4pXLhX1aOLbum0qod3S+GCtwt548cwR/VIYXWulsziF5vG9NcNzpniddOLuPsQPz1aeZi0MMRBzwT4+vwmpaYerYmT+/gY3c/5XGOHpvHI/DDdWngY4z538wF+HvwyzJxVEY7qUbdk0TJd8ApkpkK7zOrXEVV49th0mvmd7+qhI/wcOjnALxucOanP3dPH+LkhbhySikeEJesizFwR4YEq5oJfujFKqwwplRh6tfawY3Ph+o+LmDTcT8t0YfI3YeatjpSah/zXv6P8siHK8bv4eHNUBmvyolz8TiFr8qM8N3Lz/O3TloQ5Nq7k8snSYh5dEGLheU0qbdM8fhcfmwqVI54LEI5CcRSO7J7CExVM9tQhS/jt77pXnz6zKMQXqyLMP6dxBmK0ZBEnVrLIyMio5pWmrKFdvKR6KSlZDK3DLGNb4s5hafR+uIA3fgyXOtuMyfAJLx6fTr9HC2ibKUw/o/oTgi2dVvXM3X0s2xhl1NQgRcVO+8p/h6Qxe2WkZJ7mqFvLcu3gNE52553eY4d0ZiyP8Oj8MA8cXrvvsW/bzTFnN/HgFUqVnlqkO1OB/lXQeI3nu2R7ShIFwL6dnHi+XxuhW0sPo/v6uP7jIt7/pZjDu/t4/Kswe+7gqXKOimDY+T7jpXiEaSdlcPabQdrenY9XnFn6Tt3NxytxVW9RhVbpwqThfnxeAbyEInDCy0H+73ClZboQiihv/1zM+6c5x4J1gSinTQsyaXh6lQny8+XFXPdxEXcd7Ge/zl5W50a56sMiznwjWCoRgVMCChbX7f/w47oIl79fxMdnZFRa+qpvlizixJKFzWVRe4M6pTD9jAxmLIswtIu30UoVMT1aeTlvTx/XfFTEu6dWnOxnrnCqlDYVKmsLorRMr/pAvKXVUCLCLQf6GXtAGr/nKa0ynDPJqz8qKpmPeQe3ZNK7TekD0C7ZHpZvqv1Zp6+Cj1T2OWFzkoolrbK9OMON2F+gVYaH43fxMfGrMAd1TeGZRWFuPbDqGQGzm0hJu0e8vu28fHlOJnlFSiCstM30cOLLAXaKm+t8h0yhS3OPmygcvd3G7eUbnf1i+m/FZKZKSbvEt39FWZOnHDUlUPKeqIICKTfn8syx6Zyyq5P0RvZK4SK3h9Zubb1kpgr7PxVg7NBoqXnHNwSV7DpWBc5ZFWFDUNlzwubeT7H/acrNuXw6JoN9c+r3N2jJIk52djbXXHMNo0ePTnQoW6VBnVIaPUnE+++QNJ79Jp8JFXSB/PavCP96v5DHj/Hz2pJiTp4a5IuzmpCWUvlZ2ZZWQ8V4ROjQ1HndlMVhslI3T9qzn/uDXrIuytAum9/z47ooQxphkpxYW0OsWgzgr4Ioq3OrPuNN9QqRGiaUH9ZFyS1SmqY525q90knA8VVs5+3p44CnAzw2P0ywWBnVp+pJkPbYwUt+CFZsipLTrPwBNytNyEoT1geivP9rMZcN3FzluF/nFD5ZWkxxVElxs+WP650P08VtJH71h2JG9EwpqW7aq72XxReUru55+MsQb/1UzDunZtDJLYUWhLUkAcd43fDir6sIhJVf/47Sv33dksWIXr5y7XM3fFzEnwXKxKP9dG1R/43dlizi5OTkcMcddyQ6DFNH2U08/HvftJJ2gZjCYmXU1CAjeqUwpl8qx/T00ffRfK7+sIjxVUzjuaXVUH8HlSmLwxywo5dIFF75PswdM0M8dpS/5MC5U0sPJ/ZOYeynRXRqKvRo5TT2LlkX5YXjazdrXF2k+4R9O3n53+wierX2UByF6z8uJK2aI8OOLTzMWhlhxaYoGT6nHcNTST2+AGdMC3LrgWlsCCoXvVPIMT1TSp1lD85JoWcrD1d+WMgZu/nISqs6Cfdr52GHTOHTZcWc3ndzIpj6fZjmfqFrCw8/ro9y1YeFdMjycNW+m0sqVw5K5aXvwlz4ttPF9fc85coPCjmjr48W6UJUlTd+KmbKyM1VlU1SpVRDN0CbJk6VXvzzI3r6uH1mEQM6eNm/cwqrcqNc9n4hu7X1lJQmAWatiJDmhSFd6nYIbu4Xmvu95Z7LD2m5OOuLdZ0t48wzz+Tll19OdBimji4flErrMj1zLn+vkIKQ8uhRzo+/ZbowZWQ6D88P8fZPFV8YVV+mfBtm78cLGPB4Ae/9WswrJ6Zz1h6lG9YnDU/nuJ1T+MfrhewxoYCZKyNMPyOjpJcMON1Vhz7VMBdcPTncT2aqsM+TBZw8Nci5e6ayQ2bVB+uxQ9PYWKj0fDCf7LvyWbGp8pLIgA5eBud4OfjZAIdNDrBrWw9PVtDge84ePkKRmk1r6hHhvD1Tefab0v+/P/KVM98I0vPBfP7xepDBnbx89o8MMuPq9fu28/LOKRks/CNCv0cL+MfrQY7t5eORI52YZq2IEI4oQ+rQ7nbdfqncsH8at80MscvD+Zz0SpBerT28OSqjVDKdvDjMqbv6SsWV7KSqS863Vv3799f58+fX+n3RaBSv18t//vMfxo4d2wCRJZ8FCxaw55sHJjoMU42c+/K4oH8q1+5XdV3+1uzqDwv58LdiFp5XszbDv4NOsnr/tIwqG8Nr6/L3CtlQqDw9omF6Ra7cFGW3R/P5+rxMOjfwtREVWXD0x8x8+REOG3UOPfsOLLVMRBaoav+K3mfVUHGsgdsko2/+jOBPEa7Yp/oz7q3RpkLlp/VRJiwI8UAV1YJltUgXJo9MZ01etF6Txc7Zngqv16kvyzZGmXh0ekISxZZIimQhIicANwE7AwNUdX7cst2Ax4CmQBTYS1UbZFCa/HxnDBhLFiaZ7NbWy0+XbLv75PAXAsxdHeHkPj5O26127TSxjgL1qSbVYFuiLhdbJoNkifpbYCROUighIinAZOB0VV0kIq2ABqtkzstzBgXLyip/FaUxpmFUdiGlSS5JkSxU9QegoqsiDwG+UdVF7uvWN2QcoVCI7OxsmjVr1pCbMcaYrU5SJIsq9ABURN4HsoEXVPV/Fb1QRM4FzgWnC2xd9O7dm7/++quOoW6lVIlKCh5N/HDexpiGFZUU0LpdcdloLSwi8pGIfFvBbXgVb0sBBgOnun+PFZGDKnqhqk5Q1f6q2j87O7sBPsG2KVWKCTTvmegwjDGNINC8J9HCXKDCmpwqNVrJQlWH1eFtq4DPVHUdgIi8A+wBTK/P2LZnHbp055fw7XSbey0ZG3+0EoYx26CopBBo3pOf9rqVP5b9iqKkN6ld22yyV0O9D1wtIhlACBgC3JfYkLYtLVu1ori4Bz8U34L4myKydXXnM8bUgEaJFuby+9JfWP79Atp06EybDl1qtYqkSBYicizwfzjtEm+LyNeqeqiq/i0i9wJf4ozZ9Y6qvp3IWLdFbdruQFZmFm8/9zCrl/5kU8oas42KRqN02LEHR512EWn+2l10aFdwm1LCoRDhUPlJXYwxWz9fahq+1MqvI7EruE2N+VJTq9yZjDHbJ6ugNsYYUy1LFsYYY6plycIYY0y1tskGbhFZCyzfglW0BtbVUzj1JRljAourtiyu2rG4amdL4+qsqhVe1bxNJostJSLzK+sRkCjJGBNYXLVlcdWOxVU7DRmXVUMZY4ypliULY4wx1bJkUbEJiQ6gAskYE1hctWVx1Y7FVTsNFpe1WRhjjKmWlSyMMcZUy5KFMcaYalmyiCMih4nIjyLyi4j8O9HxAIjIkyLyl4h8m+hY4olIJxH5RES+F5HvROTSRMcEICJ+EZknIovcuMYmOqZ4IuIVkYUi8laiY4kRkWUislhEvhaRpBmBU0Sai8grIrJERH4QkUFJEFNP93uK3XJF5LJExwUgIpe7+/y3IvK8iPjrdf3WZuEQES/wE3AwzqRLXwKjVPX7BMe1P5APPKOqfRIZSzwR2QHYQVW/EpEsYAEwIgm+LwGaqGq+iPiAmcClqvpFIuOKEZF/Af2Bpqp6VKLjASdZAP1jk4wlCxF5GvhcVR8XkVQgQ1U3JjisEu4xYzUwUFW35CLg+oilA86+vouqBkXkJZwpHZ6qr21YyWKzAcAvqvqbqoaAF4CqpnxtFKr6GbAh0XGUpaq/q+pX7v084AegQ2KjAnXkuw997i0pzohEpCNwJPB4omNJdiLSDNgfeAJAVUPJlChcBwG/JjpRxEkB0kUkBcgA1tTnyi1ZbNYBWBn3eBVJcPDbGohIF2B3YG6CQwFKqnq+Bv4CPlTVpIgLuB+4GogmOI6yFPhARBaIyLmJDsa1I7AWmORW2z0uIk0SHVQZJwPPJzoIAFVdDdwNrAB+Bzap6gf1uQ1LFmaLiEgmMBW4TFVzEx0PgKpGVLUf0BEYICIJr74TkaOAv1R1QaJjqcBgVd0DOBy4yK36TLQUYA/gEVXdHSgAkqIdEcCtFjsGeDnRsQCISAucmpAdgfZAExE5rT63Yclis9VAp7jHHd3nTCXcNoGpwHOq+mqi4ynLrbb4BDgswaEA7Asc47YPvAAcKCKTExuSwz0rRVX/AqbhVMkm2ipgVVyp8BWc5JEsDge+UtU/Ex2IaxiwVFXXqmoYeBXYpz43YMlisy+B7iKyo3vWcDLwRoJjSlpuQ/ITwA+qem+i44kRkWwRae7eT8fpsLAkoUEBqnqtqnZU1S44+9bHqlqvZ351ISJN3A4KuNU8hwAJ73mnqn8AK0Wkp/vUQUBCO0+UMYokqYJyrQD2FpEM97d5EE47Yr2xaVVdqlosIhcD7wNe4ElV/S7BYSEizwNDgdYisgr4r6o+kdioAOdM+XRgsds+AHCdqr6TuJAA2AF42u2p4gFeUtWk6aaahNoC05zjCynAFFV9L7EhlbgEeM49efsN+EeC4wFKkurBwHmJjiVGVeeKyCvAV0AxsJB6HvrDus4aY4ypllVDGWOMqZYlC2OMMdWyZGGMMaZaliyMMcZUy5KFMcaYalmyMKYeuCOkXujeb+92YzRmm2FdZ42pB+74WG8l08jAxtQnuyjPmPpxB7CTe4Hiz8DOqtpHRMYAI4AmQHecwd5ScS5oLAKOUNUNIrIT8BCQDQSAc1Q14VeeGxNj1VDG1I9/4wxX3Q+4qsyyPsBIYC9gHBBwB8ebA5zhvmYCcImq7glcCTzcGEEbU1NWsjCm4X3izvmRJyKbgDfd5xcDu7kj9+4DvOwOuwGQ1vhhGlM5SxbGNLyiuPvRuMdRnN+gB9jolkqMSUpWDWVM/cgDsuryRncekKUicgI4I/qKSN/6DM6YLWXJwph6oKrrgVki8i1wVx1WcSpwlogsAr4jCab0NSaedZ01xhhTLStZGGOMqZYlC2OMMdWyZGGMMaZaliyMMcZUy5KFMcaYalmyMMYYUy1LFsYYY6r1/4CDcUm71Ym4AAAAAElFTkSuQmCC\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_col48,second_col48 = file_reader(os.path.join(outdir,\"out48.txt\"))\n",
    "first_col64,second_col64 = file_reader(os.path.join(outdir,\"out64.txt\"))\n",
    "first_col96,second_col96 = file_reader(os.path.join(outdir,\"out96.txt\"))\n",
    "\n",
    "for i in range(len(second_col64)):\n",
    "    # data64 = data48*(64/48)**4\n",
    "    # -> log10(data64) = log10(data48) + 4*log(64/48)\n",
    "    second_col64[i] += 4*mp.log10(64./48.)\n",
    "for i in range(len(second_col96)):\n",
    "    # data96 = data48*(96/48)**4\n",
    "    # -> log10(data96) = log10(data48) + 4*log(96/48)\n",
    "    second_col96[i] += 4*mp.log10(96./48.)\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(\"Plot Demonstrating 4th-order Convergence\")\n",
    "plt.xlabel(\"time\")\n",
    "plt.ylabel(r\"$\\log_{10}$(Relative error)\")\n",
    "\n",
    "ax.plot(first_col48, second_col48, 'k--', label='Nx = 48')\n",
    "ax.plot(first_col64, second_col64, 'k-', label='Nx = 64, mult by (64/48)^4')\n",
    "ax.plot(first_col96, second_col96, 'k.', label='Nx = 96, mult by (96/48)^4')\n",
    "\n",
    "legend = ax.legend(loc='lower right', shadow=True, fontsize='x-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 4: 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-ScalarWave.pdf](Tutorial-Start_to_Finish-ScalarWave.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": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:32:49.815973Z",
     "iopub.status.busy": "2021-03-07T17:32:49.812835Z",
     "iopub.status.idle": "2021-03-07T17:32:53.918708Z",
     "shell.execute_reply": "2021-03-07T17:32:53.917654Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-Start_to_Finish-ScalarWave.tex, and compiled LaTeX file to\n",
      "    PDF file Tutorial-Start_to_Finish-ScalarWave.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-ScalarWave\")"
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}