{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "# Start-to-Finish Example: Unit Testing `GiRaFFE_NRPy`: Interpolating Metric Face-Values\n", "\n", "## Author: Patrick Nelson\n", "\n", "## This module Validates the `FCVAL` routine for `GiRaFFE`.\n", "\n", "**Notebook Status:** Validated\n", "\n", "**Validation Notes:** This module will validate the routines in [Tutorial-GiRaFFE_NRPy-Metric_Face_Values](Tutorial-GiRaFFE_NRPy-Metric_Face_Values.ipynb).\n", "\n", "### NRPy+ Source Code for this module: \n", "* [GiRaFFE_NRPy/GiRaFFE_NRPy_Metric_Face_Values.py](../../edit/in_progress/GiRaFFE_NRPy/GiRaFFE_NRPy_Metric_Face_Values.py) [\\[**tutorial**\\]](Tutorial-GiRaFFE_NRPy-FCVAL.ipynb) Generates the C code to interpolate gridfunctions to cell faces.\n", "\n", "## Introduction:\n", "\n", "This notebook validates the code that will interpolate the metric gridfunctions on cell faces. These values, along with the reconstruction of primitive variables on the faces, are necessary for the Riemann solvers to compute the fluxes through the cell faces.\n", "\n", "It is, in general, good coding practice to unit test functions individually to verify that they produce the expected and intended output. We will generate test data with arbitrarily-chosen analytic functions and calculate gridfunctions at the cell centers on a small numerical grid. We will then compute the values on the cell faces in two ways: first, with our interpolator, then second, we will shift the grid and compute them analytically. Then, we will rerun the function at a finer resolution. Finally, we will compare the results of the two runs to show fourth-order convergence.\n", "\n", "When this notebook is run, the difference between the approximate and exact metric gridfunctions will be output to text files that can be found in the same directory as this notebook. These will be read in in [Step 3](#convergence), and used there to confirm convergence order of the algorithm. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\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 unit testing the FCVAL algorithm\n", " 1. [Step 1.a](#expressions) Write expressions for the metric gridfunctions\n", " 1. [Step 1.b](#ccodekernels) Generate C functions to calculate the gridfunctions\n", " 1. [Step 1.c](#free_parameters) Set free parameters in the code\n", "1. [Step 2](#mainc): `FCVAL_unit_test.c`: The Main C Code\n", " 1. [Step 2.a](#compile_run): Compile and run the 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": [ "\n", "\n", "# Step 1: Set up core functions and parameters for unit testing the FCVAL algorithm \\[Back to [top](#toc)\\]\n", "$$\\label{setup}$$\n", "\n", "We'll start by appending the relevant paths to `sys.path` so that we can access sympy modules in other places. Then, we'll import NRPy+ core functionality and set up a directory in which to carry out our test. " ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "execution": { "iopub.execute_input": "2020-12-17T02:17:20.339920Z", "iopub.status.busy": "2020-12-17T02:17:20.339450Z", "iopub.status.idle": "2020-12-17T02:17:20.525246Z", "shell.execute_reply": "2020-12-17T02:17:20.524741Z" } }, "outputs": [], "source": [ "import os, sys # Standard Python modules for multiplatform OS-level functions\n", "# First, we'll add the parent directory to the list of directories Python will check for modules.\n", "nrpy_dir_path = os.path.join(\"..\")\n", "if nrpy_dir_path not in sys.path:\n", " sys.path.append(nrpy_dir_path)\n", "nrpy_dir_path = os.path.join(\"..\",\"..\")\n", "if nrpy_dir_path not in sys.path:\n", " sys.path.append(nrpy_dir_path)\n", "\n", "from outputC import outCfunction, lhrh # NRPy+: Core C code output module\n", "import sympy as sp # SymPy: The Python computer algebra package upon which NRPy+ depends\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 reference_metric as rfm # NRPy+: Reference metric support\n", "import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface\n", "\n", "out_dir = \"Validation/\"\n", "cmd.mkdir(out_dir)\n", "subdir = \"FCVAL\"\n", "cmd.mkdir(os.path.join(out_dir,subdir))\n", "\n", "thismodule = \"Start_to_Finish_UnitTest-GiRaFFE_NRPy-Metric_Face_Values\"\n", "\n", "# Set the finite-differencing order to 2\n", "par.set_parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\", 2)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 1.a: Write expressions for the metric gridfunctions \\[Back to [top](#toc)\\]\n", "$$\\label{expressions}$$\n", "\n", "Now, we'll choose some functions with arbitrary forms to generate test data. We'll need to set ten gridfunctions, so expressions are being pulled from several previously written unit tests.\n", "\n", "\\begin{align}\n", "\\gamma_{xx} &= ax^3 + by^3 + cz^3 + dy^2 + ez^2 + f \\\\\n", "\\gamma_{yy} &= gx^3 + hy^3 + lz^3 + mx^2 + nz^2 + p \\\\\n", "\\gamma_{zz} &= px^3 + qy^3 + rz^3 + sx^2 + ty^2 + u. \\\\\n", "\\gamma_{xy} &= a \\exp\\left(-\\left((x-b)^2+(y-c)^2+(z-d)^2\\right)\\right) \\\\\n", "\\gamma_{xz} &= f \\exp\\left(-\\left((x-g)^2+(y-h)^2+(z-l)^2\\right)\\right) \\\\\n", "\\gamma_{yz} &= m \\exp\\left(-\\left((x-n)^2+(y-o)^2+(z-p)^2\\right)\\right), \\\\\n", "\\beta^x &= \\frac{2}{\\pi} \\arctan(ax + by + cz) \\\\\n", "\\beta^y &= \\frac{2}{\\pi} \\arctan(bx + cy + az) \\\\\n", "\\beta^z &= \\frac{2}{\\pi} \\arctan(cx + ay + bz) \\\\\n", "\\alpha &= 1 - \\frac{1}{2+x^2+y^2+z^2} \\\\\n", "\\end{align}\n" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "execution": { "iopub.execute_input": "2020-12-17T02:17:20.542121Z", "iopub.status.busy": "2020-12-17T02:17:20.537300Z", "iopub.status.idle": "2020-12-17T02:17:20.599024Z", "shell.execute_reply": "2020-12-17T02:17:20.598582Z" } }, "outputs": [], "source": [ "a,b,c,d,e,f,g,h,l,m,n,o,p,q,r,s,t,u = par.Cparameters(\"REAL\",thismodule,[\"a\",\"b\",\"c\",\"d\",\"e\",\"f\",\"g\",\"h\",\"l\",\"m\",\"n\",\"o\",\"p\",\"q\",\"r\",\"s\",\"t\",\"u\"],1e300)\n", "M_PI = par.Cparameters(\"#define\",thismodule,[\"M_PI\"], \"\")\n", "\n", "gammaDD = ixp.register_gridfunctions_for_single_rank2(\"AUXEVOL\",\"gammaDD\",\"sym01\",DIM=3)\n", "betaU = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\",\"betaU\",DIM=3)\n", "alpha = gri.register_gridfunctions(\"AUXEVOL\",\"alpha\")\n", "\n", "par.set_parval_from_str(\"reference_metric::CoordSystem\",\"Cartesian\")\n", "rfm.reference_metric()\n", "x = rfm.xx_to_Cart[0]\n", "y = rfm.xx_to_Cart[1]\n", "z = rfm.xx_to_Cart[2]\n", "\n", "gammaDD[0][0] = a*x**3 + b*y**3 + c*z**3 + d*y**2 + e*z**2 + f\n", "gammaDD[1][1] = g*x**3 + h*y**3 + l*z**3 + m*x**2 + n*z**2 + o\n", "gammaDD[2][2] = p*x**3 + q*y**3 + r*z**3 + s*x**2 + t*y**2 + u\n", "gammaDD[0][1] = a * sp.exp(-((x-b)**2 + (y-c)**2 + (z-d)**2))\n", "gammaDD[0][2] = f * sp.exp(-((x-g)**2 + (y-h)**2 + (z-l)**2))\n", "gammaDD[1][2] = m * sp.exp(-((x-n)**2 + (y-o)**2 + (z-p)**2))\n", "\n", "betaU[0] = (sp.sympify(2)/M_PI) * sp.atan(a*x + b*y + c*z)\n", "betaU[1] = (sp.sympify(2)/M_PI) * sp.atan(b*x + c*y + a*z)\n", "betaU[2] = (sp.sympify(2)/M_PI) * sp.atan(c*x + a*y + b*z)\n", "\n", "alpha = sp.sympify(1) - sp.sympify(1) / (sp.sympify(2) + x**2 + y**2 + z**2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 1.b: Generate C functions to calculate the gridfunctions \\[Back to [top](#toc)\\]\n", "$$\\label{ccodekernels}$$\n", "\n", "Here, we will use the NRPy+ function `outCfunction()` to generate C code that will calculate our metric gridfunctions over an entire grid. We will also call the function to generate the function we are testing. " ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "execution": { "iopub.execute_input": "2020-12-17T02:17:20.631473Z", "iopub.status.busy": "2020-12-17T02:17:20.622112Z", "iopub.status.idle": "2020-12-17T02:17:20.694966Z", "shell.execute_reply": "2020-12-17T02:17:20.695401Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Output C function calculate_metric_gfs() to file Validation/calculate_metric_gfs.h\n" ] } ], "source": [ "metric_gfs_to_print = [\\\n", " lhrh(lhs=gri.gfaccess(\"aux_gfs\",\"gammaDD00\"),rhs=gammaDD[0][0]),\\\n", " lhrh(lhs=gri.gfaccess(\"aux_gfs\",\"gammaDD01\"),rhs=gammaDD[0][1]),\\\n", " lhrh(lhs=gri.gfaccess(\"aux_gfs\",\"gammaDD02\"),rhs=gammaDD[0][2]),\\\n", " lhrh(lhs=gri.gfaccess(\"aux_gfs\",\"gammaDD11\"),rhs=gammaDD[1][1]),\\\n", " lhrh(lhs=gri.gfaccess(\"aux_gfs\",\"gammaDD12\"),rhs=gammaDD[1][2]),\\\n", " lhrh(lhs=gri.gfaccess(\"aux_gfs\",\"gammaDD22\"),rhs=gammaDD[2][2]),\\\n", " lhrh(lhs=gri.gfaccess(\"aux_gfs\",\"betaU0\"),rhs=betaU[0]),\\\n", " lhrh(lhs=gri.gfaccess(\"aux_gfs\",\"betaU1\"),rhs=betaU[1]),\\\n", " lhrh(lhs=gri.gfaccess(\"aux_gfs\",\"betaU2\"),rhs=betaU[2]),\\\n", " lhrh(lhs=gri.gfaccess(\"aux_gfs\",\"alpha\"),rhs=alpha),\\\n", " ]\n", "\n", "desc = \"Calculate the metric gridfunctions\"\n", "name = \"calculate_metric_gfs\"\n", "outCfunction(\n", " outfile = os.path.join(out_dir,name+\".h\"), desc=desc, name=name,\n", " params =\"const paramstruct *restrict params,REAL *restrict xx[3],REAL *restrict auxevol_gfs\",\n", " body = fin.FD_outputC(\"returnstring\",metric_gfs_to_print,params=\"outCverbose=False\").replace(\"IDX4\",\"IDX4S\"),\n", " loopopts=\"AllPoints,Read_xxs\")\n", "\n", "import GiRaFFE_NRPy.GiRaFFE_NRPy_Metric_Face_Values as FCVAL\n", "FCVAL.GiRaFFE_NRPy_FCVAL(os.path.join(out_dir,subdir))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 1.c: Set free parameters in the code \\[Back to [top](#toc)\\]\n", "$$\\label{free_parameters}$$\n", "\n", "We also need to create the files that interact with NRPy's C parameter interface. " ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "execution": { "iopub.execute_input": "2020-12-17T02:17:20.698542Z", "iopub.status.busy": "2020-12-17T02:17:20.698111Z", "iopub.status.idle": "2020-12-17T02:17:20.700388Z", "shell.execute_reply": "2020-12-17T02:17:20.700013Z" } }, "outputs": [], "source": [ "# Step 3.d.i: Generate declare_Cparameters_struct.h, set_Cparameters_default.h, and set_Cparameters[-SIMD].h\n", "# par.generate_Cparameters_Ccodes(os.path.join(out_dir))\n", "\n", "# Step 3.d.ii: Set free_parameters.h\n", "with open(os.path.join(out_dir,\"free_parameters.h\"),\"w\") as file:\n", " file.write(\"\"\"\n", "// Override parameter defaults with values based on command line arguments and NGHOSTS.\n", "params.Nxx0 = atoi(argv[1]);\n", "params.Nxx1 = atoi(argv[2]);\n", "params.Nxx2 = atoi(argv[3]);\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] = {-1.0,-1.0,-1.0};\n", "const REAL xxmax[3] = { 1.0, 1.0, 1.0};\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", "printf(\"dxx0,dxx1,dxx2 = %.5e,%.5e,%.5e\\\\n\",params.dxx0,params.dxx1,params.dxx2);\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", "par.generate_Cparameters_Ccodes(os.path.join(out_dir))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 2: `FCVAL_unit_test.c`: The Main C Code \\[Back to [top](#toc)\\]\n", "$$\\label{mainc}$$\n", "\n" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "execution": { "iopub.execute_input": "2020-12-17T02:17:20.703482Z", "iopub.status.busy": "2020-12-17T02:17:20.702797Z", "iopub.status.idle": "2020-12-17T02:17:20.705800Z", "shell.execute_reply": "2020-12-17T02:17:20.705430Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Writing Validation//FCVAL_unit_test.c\n" ] } ], "source": [ "%%writefile $out_dir/FCVAL_unit_test.c\n", "// These are common packages that we are likely to need.\n", "#include \"stdio.h\"\n", "#include \"stdlib.h\"\n", "#include \"math.h\"\n", "#include \"string.h\" // Needed for strncmp, etc.\n", "#include \"stdint.h\" // Needed for Windows GCC 6.x compatibility\n", "#include // Needed to set a random seed.\n", "\n", "#define REAL double\n", "#include \"declare_Cparameters_struct.h\"\n", "\n", "const int NGHOSTS = 3;\n", "\n", "REAL a,b,c,d,e,f,g,h,l,m,n,o,p,q,r,s,t,u;\n", "\n", "// Standard NRPy+ memory access:\n", "#define IDX4S(g,i,j,k) \\\n", "( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * ( (k) + Nxx_plus_2NGHOSTS2 * (g) ) ) )\n", "\n", "const int kronecker_delta[4][3] = { { 0,0,0 },\n", " { 1,0,0 },\n", " { 0,1,0 },\n", " { 0,0,1 } };\n", "\n", "// Give gridfunctions their names:\n", "#define GAMMADD00GF 0\n", "#define GAMMADD01GF 1\n", "#define GAMMADD02GF 2\n", "#define GAMMADD11GF 3\n", "#define GAMMADD12GF 4\n", "#define GAMMADD22GF 5\n", "#define BETAU0GF 6\n", "#define BETAU1GF 7\n", "#define BETAU2GF 8\n", "#define ALPHAGF 9\n", "#define GAMMA_FACEDD00GF 10\n", "#define GAMMA_FACEDD01GF 11\n", "#define GAMMA_FACEDD02GF 12\n", "#define GAMMA_FACEDD11GF 13\n", "#define GAMMA_FACEDD12GF 14\n", "#define GAMMA_FACEDD22GF 15\n", "#define BETA_FACEU0GF 16\n", "#define BETA_FACEU1GF 17\n", "#define BETA_FACEU2GF 18\n", "#define ALPHA_FACEGF 19\n", "#define GAMMAUU00GF 20\n", "#define GAMMAUU01GF 21\n", "#define GAMMAUU02GF 22\n", "#define GAMMAUU11GF 23\n", "#define GAMMAUU12GF 24\n", "#define GAMMAUU22GF 25\n", "#define GAMMA_FACEUU00GF 26\n", "#define GAMMA_FACEUU01GF 27\n", "#define GAMMA_FACEUU02GF 28\n", "#define GAMMA_FACEUU11GF 29\n", "#define GAMMA_FACEUU12GF 30\n", "#define GAMMA_FACEUU22GF 31\n", "#define NUM_AUXEVOL_GFS 32\n", "\n", "#include \"calculate_metric_gfs.h\"\n", "#include \"FCVAL/interpolate_metric_gfs_to_cell_faces.h\"\n", "\n", "int main(int argc, const char *argv[]) {\n", " paramstruct params;\n", "#include \"set_Cparameters_default.h\"\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", "#include \"set_Cparameters-nopointer.h\"\n", "\n", " // Step 0e: 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\n", "\n", "## Step 2.a: Compile and run the code \\[Back to [top](#toc)\\]\n", "$$\\label{compile_run}$$\n", "\n", "Now that we have our file, we can compile it and run the executable." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "execution": { "iopub.execute_input": "2020-12-17T02:17:20.711697Z", "iopub.status.busy": "2020-12-17T02:17:20.709852Z", "iopub.status.idle": "2020-12-17T02:17:21.364427Z", "shell.execute_reply": "2020-12-17T02:17:21.364844Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Now compiling, should take ~2 seconds...\n", "\n", "Compiling executable...\n", "(EXEC): Executing `gcc -Ofast -fopenmp -march=native -funroll-loops Validation/FCVAL_unit_test.c -o Validation/FCVAL_unit_test -lm`...\n", "(BENCH): Finished executing in 0.20867443084716797 seconds.\n", "Finished compilation.\n", "Finished in 0.21781635284423828 seconds.\n", "\n", "\n", "Now running...\n", "\n", "(EXEC): Executing `taskset -c 0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15 ./Validation/FCVAL_unit_test 10 10 10`...\n", "(BENCH): Finished executing in 0.20765209197998047 seconds.\n", "(EXEC): Executing `taskset -c 0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15 ./Validation/FCVAL_unit_test 20 20 20`...\n", "(BENCH): Finished executing in 0.20685338973999023 seconds.\n", "Finished in 0.43360447883605957 seconds.\n", "\n", "\n" ] } ], "source": [ "import time\n", "\n", "print(\"Now compiling, should take ~2 seconds...\\n\")\n", "start = time.time()\n", "cmd.C_compile(os.path.join(out_dir,\"FCVAL_unit_test.c\"), os.path.join(out_dir,\"FCVAL_unit_test\"))\n", "end = time.time()\n", "print(\"Finished in \"+str(end-start)+\" seconds.\\n\\n\")\n", "\n", "print(\"Now running...\\n\")\n", "start = time.time()\n", "cmd.Execute(os.path.join(\"Validation\",\"FCVAL_unit_test\"), \"10 10 10\",\"out.txt\")\n", "# To do a convergence test, we'll also need a second grid with twice the resolution.\n", "cmd.Execute(os.path.join(\"Validation\",\"FCVAL_unit_test\"), \"20 20 20\",\"out.txt\")\n", "end = time.time()\n", "print(\"Finished in \"+str(end-start)+\" seconds.\\n\\n\")\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\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}$$\n", "\n", "Here, we import the data at two resolutions and wrote to text files. This data consists of the absolute error of a metric gridfunction at each in the grid. We'll plot a portion of this data along the axis at the lower resolution along with that same data at the higher resolution scaled to demonstrate that this error converges to 0 at the expected rate. Since our algorithm uses a third-order polynomial, we expect fourth-order convergence here." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "execution": { "iopub.execute_input": "2020-12-17T02:17:21.373006Z", "iopub.status.busy": "2020-12-17T02:17:21.372463Z", "iopub.status.idle": "2020-12-17T02:17:21.783223Z", "shell.execute_reply": "2020-12-17T02:17:21.783450Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEZCAYAAAC99aPhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA19ElEQVR4nO3dd3gUVdvH8e+dBEhIKAECSI1BuoCQiGAFaaIIPCIoRUVBBLs+YkGUIhZs7yMWFEGpYgNEBEQREBRRE6UYFJCqQSC0hACBlPP+MUlcQkJ2s2V2k/tzXXOxO3N29pdhN3emnSPGGJRSSqniCrI7gFJKqcCmhUQppZRbtJAopZRyixYSpZRSbtFCopRSyi1aSJRSSrlFC4lSSim3lMpCIiL3iki8iJwSkeleWP8qEUkXkbScaYun30MppfxFqSwkwF5gAvCeF9/jXmNMRM7U2Ivvo5RStiqVhcQYM98Y8xlwKP8yEekhIutF5KiIrBWRlr5PqJRSgaNUFpLCiEhrrL2Uu4CqwDvA5yJSrhire15EDorI9yLSwXMplVLKv2ghOdMw4B1jzI/GmCxjzAzgFNDOxfU8BsQAtYEpwCIRaeDZqEop5R+0kJypPvDfnMNaR0XkKFAXqAV5J9FNIdN3uSvJKUTHjDGncorR98C1dvxASinlbSF2B/AzfwHPGmOeLWihMaZDMddrACluKKWU8melco9EREJEJBQIBoJFJFREQoB3geEicolYwkXkOhGp4MK6K4tIt9x1ishA4ErgS+/8NEopZa9SWUiA0cBJ4HFgUM7j0caYeOBO4A3gCPAnMNjFdZfBurQ4GTgI3Af0NsZs9UhypZTyM6IDWymllHJHad0jUUop5SGl7mR7tWrVTHR0tN0xlFIqoCQkJBw0xkQVtKzUFZLo6Gji4+PtjqGUUgFFRHYXtkwPbSmllHKL3xUSZ3vmFZHBIpLl0MNumnZFopRSvuePh7Zye+btBoQV0fYHY8zl3o+klFKqMH5XSIwx8wFEJA6oY3McpZRSRfC7QuKi1iJyEDgMzAKeN8Zk5m8kIsOwOmSkXr16vk2oSqXU1FQOHDhARkaG3VGUckqZMmWoXr06FStWdPm1gVxIVgMXAruB5sBHQCbwfP6GxpgpWL3wEhcXp3dgKq9KTU1l//791K5dm7CwMES0mzXl34wxnDx5kqSkJACXi4nfnWx3ljFmhzFmpzEm2xizCRgP3Gh3LqUOHDhA7dq1KV++vBYRFRBEhPLly1O7dm0OHDjg8usDeY8kP+1hN5BkZ0HyH5CUAOmpUKM51GwJ4VXtTua2jIwMwsKKuk5EKf8TFhZWrMOxfldIcnrhDcGhZ14gM/+5DxHpDvxijNkvIk2Ap4BPfB5YOSclySoaSfGQ9Avs/RVOp53drmJtqNnizKlyNAQF1s6z7omoQFTcz63fFRKsnnnHODwfBIwTkfeAzUAzY8weoBMwXUQigP3AbOA5X4dVBTh1zCoUf8fnFI8EOPaPtSy4rFUcLhoAteOgThyEVoL9v8G+Tf9O274GkwXA6eBwTlVrRoX6ra29lpotoHpTCCnOCMhKKU/zu0JijBkLjC1kcYRDu0eAR3wQSRVl/2b460drb+PvBOuQFTnXNFSJgegrrIJRO9YqAgUVgJgO1pQrI53fN/7EnIWLaJS1k+b79nDRoTkEZ06xlgeFQLXG1vpiB0P99t79GUuwwYMHU6dOHSZMmGB3FBWg/K6QqADy10+w4hnYudp6HlbFKhjNe1t7G7XbQPkqxVt3mVBWpNbig4wOZJsOBAs8fFVD7mkVfOaey59fw8aPoO0w6PQ0lIsoet1KKY/SQqJc988GWDEBtn0F4VHQ9Vloci1Eng8ePDfQLqYqZUOCyMjMpkxIEO0aREHVSKjawCpWAKfSrGL24zuwdSlcPwkadPRYBqVU0QLrDKay14Hf4aNb4J0rrb2RTmPggQ1w6b3WISwPn2COrR/JnKHteLhrY+YMbUds/cizG5WLgO4T4fal1vmXWb3h8/sgPcWjWXwpYfcR3lz5Jwm7j3hl/b/++itt2rShQoUK3HTTTaSnpwMwceJELrnkEjIzretaJk+eTPPmzfOWK1UYLSSqSL9t+pUtb92Meas9bF8JVz0OD26EKx6GsuFefe/Y+pHc0/GCgouIo/rtYfh3cNkD8OtseLMdbF3m1WzekLD7CAOnruOVr7YwcOo6jxeT06dP07t3b2655RYOHz5M3759mTdvHgAjR46kXLlyTJgwgW3btjFq1Chmz55NaGioRzOokkcLiSrc0b9I/uAumnx6NfX2f8O72dezvs9q6PiEdaWVvykTBl3Gw9DlVr4P+sH8YXDisN3JnLZuxyFOZ2aTbSAjM5t1Ow55dv3r1pGRkcGDDz5ImTJluPHGG7n44osBCAoKYubMmUyaNImePXvy6KOP0rp1a4++vyqZtJCosx3bB0tGwuttqLJtHrOyunDlqf8xMeNmvt+bbXe6otWOhbu+haseg9/mwZttYfNCu1M5Jfe8ULBgnReK8ewNmnv37qV27dpn3C9Qv379vMfR0dF07NiRXbt2cc8993j0vVXJpYVE/ev4IfjqKXjtIoh/Dy4aQOKN3zJRbuewVPbKLzavCSkHHUfBsFVQsRZ8fKs1pbne/YMvOXVeyA3nnXceSUlJGPNvl3N79uzJe7x48WJ++OEHOnXqxMiRIz363qrk0qu2lHVieu0bsO4tOH0cWt4EHR6DKjG0BOYMrcW6HYdoF1PV47/YvK5mCxi6AtZOglXPW5cqXzMRWvbz+MUBnhJbP9Jr27l9+/aEhIQwadIk7r77bhYtWsRPP/1Ex44dOXjwIEOHDmXatGlccskltGjRgl69enHttdd6JYsqOcTxL5PSIC4uzuiY7Q62LYeFd0PafmjWGzo8AdWb2J3KO5K3wMJ74e+foGE36PF/UKm2x9/m999/p2nTph5fr6fEx8dz55138ueff+YViYYNG7J582aqV6/O22+/DcDSpUsZMmQImzZtomrVANkTVW4r7PMrIgnGmLiCXqOFpLTKSIflY+HHyRDVFP4zGWqVghOr2VnWPSffjIfgMtDtWWhzq0ffwt8LiVLnUpxCooe2SqP9iTBvKBzYDJcMh85jrSueSoOgYGh/NzS+Bj6/37rn5FSaNU8pVSx6sr00MQbWvQ1TOsLxgzDwU+tmvtJSRBxViYFbF0KTHrBsFCR+ZncipQKWFpLS4th+mHMjfPmY1YXIiLXQsIvdqewVFAx9pkKdi637TXb/YHcipQKSFpLS4I8lMLk97PoernsF+n8IEVF2p/IPZcKs7VG5Lsy9GZK32p1IqYCjhaQkO30CvngIPuxv3Utx17dw8VC/vezVNuFVrcN8wWVgdh9r700p5TQtJCXVPxtgylXWjYWX3g9Dv4Goxnan8l9VzocBH8OJg/BBX+sEvFLKKVpISprsbPh+ErzbyRqp8NaF0PUZHU3QGbXbQN/p1jgnn9wGWa6PXa1UaaSFpCRJ3QuzesHXT0Hj7tYJdcdRB1XRGnWD616FP5dbhwVL2X1WShWHFpKSYvNCeKu9NdRtzzeg38zij05Y2sXdDleOhF9n8eP0x7w2Loi/GDx4MKNHj7Y7Bnv27CEiIoKsrCxbcyQnJ9OkSRNOnjxpa478Xn/9dR577LFCl0+bNo233nrLh4n+pYUk0GVlwOJHrA4Jq8TA8DXQ5hY9oe6mhJi7+Sz7Si7Z/Q7zpr1Q4ouJK+bOncuAAQPYunUrvXr1IioqiipVqtCtWze2bNni9Hqio6NZvnx53vN69eqRlpZGcHCwN2I77YUXXmDw4MGEhVn3Vz3yyCM0bNiQChUq0KRJE2bOnHlG+/Xr1xMbG0v58uWJjY1l/fr1Z62zcePGbN26lZUrV9KxY0cqVapEdHT0We127dpFx44dKV++PE2aNDlj+9x5553MmTOHAwfO7nj0pZdeYsyYMbzyyivn/KNg27ZthIaGMmjQICe3hnO0kASytGSY0RN+fhfa3wtDvrKGoVVuW7fzMI9lDGV1VgvGyxSS4hfZHclvLF68mGuvvZajR4/Ss2dPtmzZwv79+2nbti29evWyO55bTp06xYwZM874RRseHs6iRYtISUlhxowZPPDAA6xduxawBgrr1asXgwYN4siRI9x222306tWL06dP571++/btZGVl0ahRI8LDw7njjjt46aWXCnz//v3707p1aw4dOsSzzz7LjTfeSHJyMgChoaF07979rEI2Y8YMJk+ezOrVq1m9ejXz5s3jjTfeKHD999xzT974Mx5ljClVU2xsrCkRkn4x5pWmxjxTw5iNn9idpsSJ33XYNB69xLR8/GOz+ekWJvOZ84xJ+tWp127evNm74dz0yy+/mNatW5uIiAjTr18/c9NNN5knn3zSGGPMCy+8YNq2bWsyMjKMMca89dZbplmzZubkyZPGGGOysrJM9erVTXJy8lnrPXTokAHMwYMHi8wwaNAgIyImNDTUhIeHm4kTJ5qdO3caIO+9r7rqKvPkk0+a9u3bm/DwcNOjRw9z8OBBM2DAAFOhQgUTFxdndu7cmbfO33//3XTu3NlERkaaRo0amY8++ihv2eLFi03Tpk1NRESEqVWrlnnppZcKzPXtt9+aBg0anDP79ddfb15++WVjjDHLli0ztWrVMtnZ2XnL69ata5YuXZr3/LXXXjP33XffGev4+uuvTf369c+Yt2XLFlO2bFmTmpqaN+/yyy83kydPzns+e/Zs06FDh7znX3zxhWnatKnZs2dP3rz9+/ebVq1amY8//viM9c+dO9f07dvXjBkzxgwcOLDQn6+wzy8Qbwr5vap9bQWiDR/CogcgvDoMWQbntbI7UYmTOy7Iuh2HyKj5McFL+1gjLg75GiLrF70CR0sft64E86aaLaD7C0U2yx1q98EHH+Tee+9l4cKF9O/fP+/Y+8iRI1m8eDETJkxg4MCBjBo1ihUrVuQNt/vTTz8RExNDtWrVzlr36tWrqVmzplM9Bc+aNYs1a9YwdepUOnfuDFiHdfL78MMPWbZsGdWqVaN9+/a0b9+et956ixkzZnDHHXcwbtw43n//fY4fP06XLl0YP348S5cuZdOmTXTp0oULL7yQZs2aMWTIED7++GOuuOIKjhw5ws6dOwvMtWnTJho3Lvwy+ZMnT/Lzzz9z991W32yJiYm0bNnyjIHCWrZsSWJiItdccw0AS5Ys4aGHHipymyQmJhITE0OFChXy5rVq1YrExMS8502bNmXDhg15z6+77jquu+66M9ZTvXr1sw6vpaam8vTTT7NixQqmTp1aZBZX6aGtQJKVCV8+AQvusrr1GLZKi4gX5Y4X37JpExj0KWSmW93MBNDQvfmda6hdKHq43dzDWvn9/fff3HPPPbz66qsezXv77bfToEEDKlWqRPfu3WnQoAGdO3cmJCSEvn378uuvvwLwxRdfEB0dze23305ISAitW7emT58+fPLJJwCUKVOGzZs3k5qaSmRkJG3atCnw/Y4ePXrGL/L8hg8fTqtWrejWrRsAaWlpVKp05rDTlSpV4tixYwCcOHGCn3/+mQ4dOhT5sxa1LoAKFSqQkpJS5Lrye+qppxgyZAh16tRx+bXO8Ls9EhG5FxgMtADmGmMGn6PtQ8BjQHngU2CEMeaUD2L63Pot26m29C7qHP0ZLhlh3RsSXMbuWKVH9aZw8wcw6z/w4QC45TMoE+rca53YU/CVoobahX+H212yZMlZw+0uWbKEKVOmnDEvOTmZrl27cvfdd9O/f3+P5q1Ro0be47CwsLOep6VZN47u3r2bH3/8kcqVK+ctz8zM5JZbbgFg3rx5TJgwgccff5yWLVvywgsv0L59+7PeLzIy8oxf3I5GjhzJb7/9xsqVK/O2X0REBKmpqWe0S01NzStG33zzDZdeeinlyhV9H1dR6wI4duzYWcWmKOvXr2f58uV5Rdcb/HGPZC8wAXjvXI1EpBvwONAJqA/EAOO8ns4Gm3/5jqgPuhJ1ZD2PZ48godljWkTsEH059J4Me36ABcOsmz8DTFFD7ULhw+3u27ePf/7554y/5o8cOULXrl3p2bMnTz75pEtZxINXFtatW5errrqKo0eP5k1paWlMnjwZgIsvvpiFCxdy4MABevfuTb9+/QpcT8uWLdm69ez+1saMGcPSpUv56quvqFixYt785s2bs3HjxjO258aNG2nevDlgFV5nR5hs3rw5O3bsOKOQbdiwIW9dYI0V0qqVa0chVq1axa5du6hXrx41a9bk5ZdfZt68eYXulRWH3xUSY8x8Y8xnwKEimt4GTDPGJBpjjgDPYO3JlCybPqXhFzcgZNP39NN8knEF63YUtWmU17S4Ebo8Y92385X99164ynGo3YyMDObPn89PP/2Utzx3uN2pU6cyY8YMFi1axJIlSwBrxMRrrrkmrwCkpqbSrVs3LrvsMl544ey9rlWrVp2zWNSoUYMdO3Z45Ofq0aMHW7duZdasWWRkZJCRkcHPP//M77//zunTp5kzZw4pKSmUKVOGihUrEhRU8K++tm3bcvToUZKSkvLmPf/883zwwQcsX778rPM/HTp0IDg4mEmTJnHq1Km8q6WuvvpqwNpmjucwsrOzSU9PJyMjA2MM6enpeVd4NWrUiIsuuohx48aRnp7OggUL2LhxI3369Ml7/bfffkv37t1d2jbDhg1j+/btrF+/nvXr1zN8+HCuu+46li1b5tJ6zsXvCokLmgMbHJ5vAGqIyFln+kRkmIjEi0h87qV0fi87C75+GuYNIT2qJX2znyORBpQJCaJdjA57aqtL74O2d8G6N+HHKUW39yNly5Zl/vz5TJ8+nSpVqvDRRx9xww035C0fNmxY3jjtVatWZdq0aQwdOpRDhw6ddX5kwYIF/Pzzz7z//vtERETkTbl7OH/99ReXXnppoVmeeOIJJkyYQOXKlXn55Zfd+rkqVKjAV199xYcffkitWrWoWbMmjz32GKdOWUe6Z82aRXR0NBUrVuTtt99mzpw5hW6fwYMHM3v27Lx5o0aNYs+ePVxwwQV5P+Nzzz2X1/6zzz5j5syZVK5cmffee4/PPvuMsmXL8ttvvxEREUG9evXy1rV69WrCwsK49tpr2bNnD2FhYXTt2jVv+Ycffkh8fDyRkZE8/vjjfPrpp0RFWT11p6ens2TJEm677TaXtk358uWpWbNm3hQREUFoaGjeej2isMu57J6wDm9NP8fy7cA1Ds/LAAaIPtd6A+Ly3+OHjJnZ25gxFY354mFjMk6Z+F2HzRsrtpn4XYftTqeMMSYr05jZfY0ZH2VM8rYzFvn75b/FkZGRYapWrWpSUlKcfs2QIUPMl19+6cVU3nHgwAHTuHFjc+LECbfWM3HiRDNy5EgPpTJm0qRJHl1fYUrb5b9pQEWH57mPCz5TFij2b7ZO5qb8DddPgljrr4/Y+mWJrR9pcziVJygYek6CN9taw/UOXgyFHC4pCQ4fPswzzzxzxvmBonjjMlNfiIqK4o8//nB7PdHR0Vx//fUeSGS57777PLYuTwvkT34i4HjWqRWw3xgTuCcQNi+EqZ0h4yTcviSviCg/VaEmdJ0Ae9ZCwvt2p/Gq6tWrM2LECLtjBJR+/frRtGlTu2P4hN8VEhEJEZFQIBgIFpFQESloz2kmMEREmolIZWA0MN13ST0oOxu+ecbqL6tGM+v+kLpt7U6lnNH6Fjj/Svh6DKQkFd1eqRLI7woJVkE4iXVp76Ccx6NFpJ6IpIlIPQBjzJfAi8BKYA+wGxhjT2Q3pKdah7LWvGz9Uhq8GCqeZ3cq5SwRuP41yM6Exf/VbudVqeR350iMMWOBsYUsjsjX9lXAs7fS+tKh7VYRObgNrn1Zh8ENVFVi4OonrcuBE+dDcDOys7MLvcRUKX+VXcx7o/STbpftK+HdqyFtP9yyANreqUUkkF0yAmq1hiWPEl4uhKSkJE6fPn3GjWpK+StjDKdPnyYpKYnw8HCXX+93eyQlnjHw49uwbBRENbG63ahyvt2plLuCQ6wBxaZcRZ0N/8fBy8aye/duMjMz7U6mlFNCQkKoVKlSgR1yFvlaL+RRhck8BV88DOtnQ5Me8J+3oVzhHcSpAFPzQrjsQYLWvEz1C2+gesPOdidSyif00JavHNsP03tYReSqx6DfLC0iJdGVI6FqQ/jiQTiVZncapXxCC4kvJP0CUzrA/t+g7wzoOKpE37xWqpUJhZ6vQ8pfsOIZu9Mo5RP628zbNn4C73eHoBBrKNzmve1OpLytfnu4+E748R3466ei2ysV4LSQeEtup4vzh0LtWBi20hrFTpUOncdAxdpW9ymZJXKIHKXyaCHxhvQUmHszfP8axN1hDYIU7vqVECqAlasAPf4Pkv+ANYF7q5NSztBC4mkH/4R3O8H2FXDdq9Yvk5CydqdSdmjUFVr0hTWvwIHf7U6jlNdoIfGkbcutmwxPHoZbF8LFQ+xOpOx2zQvW3snCe63DnUqVQFpIPGXtG/BBX6hcF+5caQ3LqlR4Neg+EZLirZPvSpVAWkg85fRx6ybDIV9BZH270yh/0qIvXNDFuhz4yC670yjlcVpIPOXKkdY9ImVd76dGlXAi1rkyCYJFD2oPwarE0ULiKUFBepOhKlzlutB5LOxYCRvm2p1GKY/S33xK+UrcEKjbDr58AtIO2J1GKY/RQuKkhN1HeHPlnyTsPmJ3FBWogoKs7lMyTsDSR+1Oo5THaO+/TkjYfYSBU9dxOjObsiFBzBnajtj6kXbHUoEoqhFc+SisnGCdhG9ynd2JlHKb7pE4Yd2OQ5zOzCbbQEZmNut2HLI7kgpklz0A1ZtbQ/Omp9idRim3aSFxQruYqpQNCSJYoExIEO1iqtodSQWykLLQ63U4ts+6612pAKeHtpwQWz+SOUPbsW7HIdrFVNXDWsp9tWOh1c2w7m2rp+DKde1OpFSx6R6Jk2LrR3JPxwu0iCjP6fik9e/K5+zNoZSbtJAoZZfKdeGSYdZ9Jfs22Z1GqWLTQqKUna74L4RWguVj7U6iVLFpIVHKTmGRVjH5cznsWGV3GqWKRQuJUnZrOwwq1bVG1MzOtjuNUi7zu0IiIlVEZIGIHBeR3SIyoJB2Y0UkQ0TSHKYYX+dVym1lQuHq0fDPBkicb3capVzmd4UEeBM4DdQABgKTRaR5IW0/MsZEOEw7fJZSKU9q0Q9qtIBvxukY7yrg+FUhEZFwoA/wlDEmzRjzHfA5cIu9yZTysqAg6DIOju6Bn6fZnUYpl/hVIQEaAZnGmK0O8zYAhe2RXC8ih0UkUURGFLZSERkmIvEiEp+cnOzJvEp5zgWdIKYjrH4RTh61O41STvO3QhIBpOablwJUKKDtx0BTIAq4E3haRPoXtFJjzBRjTJwxJi4qKsqTeZXyrC7j4OQR+P5/didRymn+VkjSgIr55lUEjuVvaIzZbIzZa4zJMsasBV4DbvRBRqW857xW1vmSdZMhJcnuNEo5xd8KyVYgREQaOsxrBSQ68VoDiFdSKeVLV48mOzub3+c+ruPfqIDgV4XEGHMcmA+MF5FwEbkM6AXMyt9WRHqJSKRY2gL3Awt9m1gpz0tIrciMzC40+mcR46Z+rMVE+T2/KiQ57gbCgAPAXGCEMSZRRK4QkTSHdjcDf2Id9poJTDTGzPB5WqU8bN2OQ0zK6EUaYTzEBzr+jfJ7fteNvDHmMNC7gPlrsE7G5z4v8MS6UoGuXUxVXg+pyNtZvXgsZC61wrYAF9gdS6lC+eMeiVKlWu74N5U63Mvp8Fo03viSdp2i/JoWEqX8UGz9SIZ3vpCyXZ6Cvb/C5gV2R1KqUFpIlPJnLW+CGhfCN+Mh87TdaZQqkBYSpfxZUDB0HgdHdkH8e3anUapAWkiU8ncXdILzr7S6TklPsTuNUmfRQqKUvxOBLuPhxCH4/jW70yh1Fi0kSgWCWq3hwhvhh7cgda/daZQ6gxYSpQJFp6cgOxNWPW93EqXOoIVEqUARGQ1t74RfZ8OBP+xOo1QeLSRKBZIrR0LZCrB8rN1JlMqjhUSpQFK+Clz+IGxdCru+tzuNUoAWEqUCT7sRUKEWfP0UGGN3GqW0kCgVcMqEQccnICkBtiy1O41SWkiUCkitBkCVGFj5rHboqGznUjfyInIRcAVQDYfRCI0xT3s2llLqnIJDoMMomD8UNn8GF95gdyJVijm9RyIiw4DvgauBx4AWwH/RgRKUsseFN0BUU1j5HGRl2p1GlWKuHNp6FLjGGPMf4GTOvzcCGV5JppQ6t6Bg6DgKDm2DTZ/YnUaVYq4Ukuo5oxQCZItIkDFmKXC9F3IppZzR9Hqo2dK62z1L/6ZT9nClkPwtItE5j7cCvUTkCkAHSVDKLiJw9VNwdLd1x7tSNnClkLwINM15PB6YDawAxnk6lFLKBQ27QN1LYPVLkJFudxpVCjldSIwx03MOZZHzbyQQaYyZ7K1wSikniMDVoyE1CRLetzuNKoVcvo9ERKqLSAxQB8h9rJSy0/lXWtOaV+D0cbvTqFLGlct/rxGRJOAf4E+HaZuXsimlXNFxNBxPhp+m2J1ElTKu7JG8CTwDRBhjghymYC9lU0q5ot4l0LArfPc/HZJX+ZQrhSQSeMcYc9JbYQBEpIqILBCR4yKyW0QGFNJORGSiiBzKmSaKiBTUVqlSo+MoSD8K6/TUpfIdVwrJNOB2bwVx8CbWJcU1gIHAZBFpXkC7YUBvoBXQEut+lrt8kE8p/1WrtXVvyQ9vwonDdqdRpcQ5C4mIrBGR1SKyGmgHvC0iW3PnOSzzCBEJB/oATxlj0owx3wGfA7cU0Pw24BVjzN/GmCTgFWCwp7IoFbA6jIJTx2DtJLuTqFKiqE4bpxbx3NMaAZnGmK0O8zYAVxXQtnnOMsd2Be255PYTNgygXr16nkmqlL+q0Qxa3Ag/vgPt7oaI6nYnUiXcOQuJMWaGr4LkiABS881LASoU0jYlX7sIERFjzhztxxgzBZgCEBcXpyMBqZKvwxPw23xY8yp0f8HuNKqEc+k+EhG5Q0S+FpHEnH+HePgEdxpQMd+8isAxJ9pWBNLyFxGlSqWqDeCiARA/DVL+tjuNKuFcuY/kRazu4+cDI4F5wCPARA/m2QqEiEhDh3mtgMQC2ibmLCuqnVKl01WPWkPxrn7Z7iSqhHNlj2Qw0MkYM9kYs8QY8zbQFQ9eyWWMOY5VqMaLSLiIXAb0AmYV0Hwm8LCI1BaRWlhjo0z3VBalAl7lehA7GH6dBYd32p1GlWCuFJJjnH2I6Rhnn9Nw191AGHAAmAuMMMYkisgVIpLm0O4dYBGwCfgNWJwzTymV64r/QlAIfPui3UlUCSbOnlIQkfuw7tt4AfgbqIt1iGshsCS3nTFmh8dTelBcXJyJj4+3O4ZSvrPsSVj3Ftz9I0Q1sjuNClAikmCMiStwmQuFJNuJZsbfu0zRQqJKneMH4X8toVE36Ku9A6viOVchcaUb+SAnJr8uIkqVSuHVoN0ISJwP+36zO40qgVzuRl4pFYAuvQ9CK8HKZ+1Ookqgc96QKCJrgCKPfRljrvRYIqWU54VVtorJignwdwLUibU7kSpBXO0iRSkVqC4ZbvUKvHIC3LLA7jSqBHGpixQRqQG0BaoB2mW7UoGkXAW4/CH4ajTs+h6iL7M7kSohXLmzvTfWiIjjse7XuC/n34J65lVK+aO4IRBR0zrEpb0JKQ9x5WT7BOAOY0xr4HjOv8OABK8kU0p5XtnycOUjsGctbF9hdxpVQrhSSOoZYz7JN28GcKsH8yilvK3NrVCpHnwzDrKduT1MqXNzpZAcyDlHArBLRNoDDQC9d0SpQBJSDq4eDf9sgN/m2Z1GlQCuFJJ3gctzHv8fsBJrMKm3PB1KKeVlLfpCzRawYjxknrI7jQpwrtzZPtEYMy/n8Uys0QxjjTFPeSucUspLgoKgy3g4ugd+nmZ3GhXgin1nuzFmjzHmd0+GUUr5UIOrIaYjrH4RTh61O40KYNpFilKlWZdxVhH5/n92J1EBTAuJUqXZea2gZT/rjncdklcVkxYSpUq7jk+CyYaVz9udRAUoLSRKlXaR9aHtMNjwAexPtDuNCkBaSJRS1pC85SrA8rF2J1EBSAuJUgrKV7GKybavYOdqu9OoAKOFRCllaXsXVKwDXz+tXacol2ghUUpZyoRaXafs/ZUd387mzZV/krD7iN2pVADQQqKU+lfLfpyIbELIqmeY9FUiA6eu02KiiqSFRCn1r6BgvqlzD/XkAP2DlpORmc26HYfsTqX8nBYSpdQZasX2YK25kPtCFhAZkk67mKp2R1J+zq8KiYhUEZEFInJcRHaLyIBztB0rIhkikuYwxfgyr1IlUWx0FSJ7Pk9VOcai1gnE1o+0O5Lyc35VSIA3gdNADWAgMFlEmp+j/UfGmAiHaYdPUipVwjWNvRJa9OW8ze9B6l674yg/5zeFRETCgT7AU8aYNGPMd8Dn6JjwStnj6tFgsmCVdp2izs1vCgnW+CaZxpitDvM2AOfaI7leRA6LSKKIjCiskYgME5F4EYlPTk72VF6lSrbIaLh4KPw6Gw78YXca5cf8qZBEAKn55qUAFQpp/zHQFIgC7gSeFpH+BTU0xkwxxsQZY+KioqI8lVepku/KkVBWu05R5+azQiIiq0TEFDJ9B6QBFfO9rCJwrKD1GWM2G2P2GmOyjDFrgdeAG737UyhVypSvApc/CFuXwq7v7U6j/JTPCokxpoMxRgqZLge2AiEi0tDhZa0AZ7sjNYB4OrdSpV67EVCxttV1ijF2p1F+yG8ObRljjgPzgfEiEi4ilwG9gFkFtReRXiISKZa2wP3AQt8lVqqUKBMGHUdBUjxs1q+YOpvfFJIcdwNhwAFgLjDCGJMIICJXiEiaQ9ubgT+xDn3NBCYaY2b4OK9SpUOr/lC9GXwzDrIy7E6j/EyI3QEcGWMOA70LWbYG64R87vMCT6wrpbwgKBg6j4MP+kLCdGh7p92JlB/xtz0SpZS/atgFoq+AVS9Aev4LLFVppoVEKeUcEegyDk4chLWv251G+REtJEop59WOheY3wA9vQMrfdqdRfkILiVLKNZ3HWv8u/q9eDqwALSRKKVdF1rf64dr6Jfw2z+40yg9oIVFKue6S4dZhrqWPwnEd+Kq000KilHJdUDD0fB3SU2DZE3anUTbTQqKUKp4azeHyh2HjR7Btud1plI20kCiliu/KR6BaY/jiQThVYP+qqhTQQqKUKr6QctYhrpS/4Ztn7E6jbKKFRCnlnnqXQNth8NMU+Osnu9MoG2ghUUq5r9PTUKkOLLwXMk/ZnUb5mBYSpZT7ykVAj//BwS2w+mW70ygf00KilPKMhp2h5U3w3auw39nx6FRJoIVEKeU53Z6H0Erw+X2QnWV3GuUjWkiUUp4TXhW6vwhJCfDj23anUT6ihUQp5VkX9oGG3WDFBDiyy+40yge0kCilPEsEerwKEgyLHtAegksBLSRKKc+rVAc6j4Edq2D9B3anUV6mhUQp5R1xQ6Bee6tTx2P77U6jvEgLiVLKO4KCrO5TMtJh6Ui70ygv0kKilPKeag3hqkdh80L4/Qu70ygv0UKilPKuyx6AGi2soXlPHrU7jfICLSRKKe8KLgM9J8HxA/D103anUV6ghUQp5X2120D7e+CXGbBzjd1plIf5TSERkXtFJF5ETonIdCfaPyQi+0QkVUTeE5FyPoiplCquDqMg8nxYdD+cPmF3GuVBflNIgL3ABOC9ohqKSDfgcaATUB+IAcZ5NZ1Syj1ly8P1r8HhHbBCB8EqSfymkBhj5htjPgMOOdH8NmCaMSbRGHMEeAYY7MV4SilPiLkKLr4T1r0F8e/bnUZ5iN8UEhc1BzY4PN8A1BCRqgU1FpFhOYfN4pOTk30SUClViGuehwu6wOKHYesyu9MoDwjUQhIBpDg8z31coaDGxpgpxpg4Y0xcVFSU18Mppc4huAz0nQ41W8AngyHpF7sTKTf5pJCIyCoRMYVM3xVjlWlARYfnuY+PuZ9WKeV15SJgwCcQXg0+6AeHd9qdSLnBJ4XEGNPBGCOFTJcXY5WJQCuH562A/cYYZ86vKKX8QYUaMGg+ZGfC7D5wXL++gcpvDm2JSIiIhALBQLCIhIpISCHNZwJDRKSZiFQGRgPTfZNUKeUx1RpC/w8h5W+YezNknLQ7kSoGvykkWMXgJNZlvYNyHo8GEJF6IpImIvUAjDFfAi8CK4E9wG5gjB2hlVJuqtcO+rwLf/8M84aSsPMgb678k4TdR+xOppwkppQNOhMXF2fi4+PtjqGUym/dZPjycWZnd+PpjFspGxLMnKHtiK0faXeywGcMbPoEGl0DoRWLbl8AEUkwxsQVtMyf9kiUUqVZuxGsrzOQQUHLGBK0mIzMbNbt0PMmbjtxGD6+BebfCfHTvPIWhZ2DUEopn8vq/AxL3tvOk2U+4FB2NdrFXGp3pMC2YxUsGA7HD0KX8dD+Pq+8je6RKKX8Rmx0VWreNp29lVrzcsjbxJpEuyMFpsxT8NVTMLM3lI2Aocut7vyDvPMrXwuJUsqvtGlwHrXumk9QlfPhwwFw4He7IwWW5C0wtROsnQRxt8Ndq6HWRV59Sy0kSin/U74KDPoUQkJh9o2Q+o/difyfMfDzVHjnSkjdCzfPhR7/Z3WW6WVaSJRS/qlyPRj4CaQfhTl9IT3V7kT+Ky3Zug9n8X+h/mUwYi00udZnb6+FRCnlv85rBf1mwIHN1pVHmaftTuR/ti2HyZfC9pVwzUQY+ClUqOnTCFpIlFL+7YLO1lC9O1ZZg2KVsnvfCpWRDksfgzl9rD7Lhq2EdsO9dkL9XPTyX6WU/2s9CFKSYNVzUKkOXD3a7kT22p8I84Zae2qXDIfOY6FMmG1xtJAopQLDVY9Cyl+w+iVIT4FOY6xehEuT7Gz46R34egyEVrIOYzXsYncqLSRKqQAhknMVUjj8+A5s/RJ6vg4xHexO5hup/8DCe2D7N1ZXJz3fgAj/GF9Jz5EopQJHcBnoPhFuXwrBZWFmL/j8fmsPpaQ6cRiWj4XX28DutXDdK1aPyX5SRED3SJRSgah+exj+Hax6Hta+Dtu+huv/B4262Z3Mc9JTrY4sf3gDTh2DFjdChyegagO7k51FC4lSKjCVCbP6j2rWCz67xxppseXN1pjw5avYnc5lCbuPsG7HIS6tG0brfZ/A9/+Dk0egSQ/o+CTUaGZ3xEJpIVFKBbbasXDXt7DmFWva/o11+KdZL7uTOS1h9xEGT11Dn+zl1An5DCQFLugCVz8JtVrbHa9Ieo5EKRX4QspBx1EwbBVUrAUf32pNaQfsTla0rAxSv3+XZUEPMrbMDLab2sy7aJrVRUwAFBHQQqKUKklqtoChK6xLg7cshTfbwsaP/fMmxuws2PAhvBFHx63PcoBIbjk9isHmKaJbd7I7nUv00JZSqmQJDoErHoYm18HCe60BnX6bZ106XLGW3emse0F+/xxWPgcHt1jFb8DHZJW9mHY7D/NgTNWAGxVSh9pVSpVc2VnWPSffjLcuHe46Adrcat2T4mtZGfDnclj5LOzbBNUaW4fjmva0pVsTV51rqF3dI1FKlVxBwdD+bmh8jXW/yaL7IWE6RF8GNVtaewNVG1p7Mfx75VQ7d/cKjIEjuyApwZr+jod/NkDWKYiMhv9MsS7nDQr2xE9pO90jUUqVDtnZ8Mt0SJhhDZaVdcqaHxIK1ZuRHNGYt/4IY1NmPXYER/Pu0A7OF5OTR3IKRgIkxVuPT+SMNx8SZg0sVTsW6raFxtdae0cB5lx7JFpIlFKlT1YGHNxmHWLatxH2bSL9r/WEZlp3yGcbIaV8PSJj2lh7Lbl7LxE1IOs07PstZ28j3trbOLw9Z8UCUY2hdhzUibWKR/VmAVk48tNDW0op5Si4jHWDX41m0OomABJ3Hea/0xbTMHsnFwbv4dbqx2Dvr5C44N/Xla9q3WWelTMuSkRNqBNn9U5cO9a6XDe0og0/kL20kCilFBAbXYVXhvbIO0cSmXtYKz3F6rZ93yZrCou0ikftOOsqMDtO3PsZvykkInIvMBhoAcw1xgw+R9vBwDTgpMPsHsaYVd5LqJQq6WLrR559XiS0EtS/1JpUgfymkAB7gQlAN8CZEVp+MMZc7t1ISimliuI3hcQYMx9AROKAOjbHUUop5ST/vwumcK1F5KCIbBWRp0Sk0KIoIsNEJF5E4pOTk32ZUSmlSrxALSSrgQuB6kAfoD8wsrDGxpgpxpg4Y0xcVJT/DAajlFIlgU8KiYisEhFTyPSdq+szxuwwxuw0xmQbYzYB44EbPZ9cKaVUUXxyjsQY08HbbwHoNXhKKWUDvzm0JSIhIhIKBAPBIhJa2HkPEekuIjVyHjcBngIW+i6tUkqpXH5TSIDRWPeFPA4Mynk8GkBE6olImojUy2nbCdgoIseBJcB84DnfR1ZKKVXq+toSkWRgt01vXw04aNN7F0WzFY9mKx7NVjx2ZqtvjCnwaqVSV0jsJCLxhXV6ZjfNVjyarXg0W/H4azZ/OrSllFIqAGkhUUop5RYtJL41xe4A56DZikezFY9mKx6/zKbnSJRSSrlF90iUUkq5RQuJUkopt2ghUUop5RYtJB4kIvfmdFd/SkSmO9H+IRHZJyKpIvKeiJRzWBYtIitF5ISI/CEind3MVkVEFojIcRHZLSIDztF2aU5PArnTaRHZ5LB8l4icdFj+lQ+zjRWRjHz5YhyWXyQiCTnbLUFELvJhtpEi8puIHBORnSIyMt9yt7ebs3nEMlFEDuVME0X+HRPW09vJxWxe305uZPPp58vFbD79XrrEGKOThybgBqA3MBmYXkTbbsB+oDkQCawCXnBY/gPwKtZokX2Ao0CUG9nmAh8BEcDlQArQ3MnXrgKedni+C+jswe3mdDZgLDC7kGVlsXoteAgoB9yf87ysj7I9CrTB6gy1cc573+zJ7eZsHuAuYAvWIHG1gc3AcG9tJxezeX07uZHNp58vVz9j+V7n1e+lSz+DHW9a0iesIYOnF9HmA+A5h+edgH05jxsBp4AKDsvX5P4iKEaecOA00Mhh3iwcCtc5XhsNZAHRDvM89oF1NVsRX/SuQBI5VyPmzNsDXOPr7ZbTdhLwuqe2myt5gLXAMIfnQ4B13thOHviMeXQ7ubndfPb5cme7eft76eqkh7bs0xzY4PB8A1BDRKrmLNthjDmWb3nzYr5XIyDTGLO1GOu7FVhjjNmVb/4cEUkWka9EpFUxcxU32/UiclhEEkVkhMP85sBGk/OtyrGxiHV5OhtgHVoCrgAS8y1yZ7u5kqegz1dzh2We3E6uZsvjpe3kbjZffb6Kky2Xt7+XLtFCYp8IrF3YXLmPKxSwLHd5BTfeK7WY67sVmJ5v3kCsv4jqAyuBZSJS2UfZPgaaAlHAncDTItLfYV3+st3GYn2/3neY5+52cyVPQZ+viJxf3J7eTq5mczQWz28nd7L58vPlajZH3v5eukQLiZPEw6M8AmlARYfnuY+PFbAsd/kxCuBENpfW57Dey4GawKeO840x3xtjThpjThhjnsc6f3OFL7IZYzYbY/YaY7KMMWuB1/h3dEx/2W73Yn3RrzPGnHLI7vR2K4QreQr6fKXl/DVdrJ/Lg9kAr26nYmfz5OfL09lyeeJ76WlaSJxkjOlgjJFCpsuLscpEwHHXsxWw3xhzKGdZjIhUyLc8/+6/s9m2AiEi0tCZ9Tm4DZhvjEkrol2hI1R6MVtB750ItHS8OgloWdi6vJFNRO7AGlOnkzHmbxeyO8OVPAV9vhIdljm9nbyQzdvbya1s53hv27dbDre/lx5nx4mZkjphXYUSCjyPdcIsFAgppO01wD6gGVAZWMGZV22tA17OWcd/cP+qrQ+xrg4JBy6jiCtDsK4WSwGuzje/Xs7ry+ZkGwkkA1V9kQ3ohXWVmwBtsU5+3pazLPeqmgewrqq5F/ev2nIl28Cc/9OmBSzzyHZzNg8wHPgd64qtWli/mPJfteWx7eRiNq9vJzey+fTzVYzPmM++ly79DL54k9IyYR3vNfmmsQ7/0WlAPYf2D2NdApyKdYy4nMOyaKzL+05iXcbp7mWjVYDPgONYV5oMcFh2BdZhD8f2/XO+JJJvfnOsE4zHgUPAN0Ccr7LlfOEO5WzLP4D7862rNZCQs91+AVr7MNtOICMnW+70tie3W2F5CsgiwIvA4ZzpRc682sij28nFbF7fTm5k8+nny5Vsvv5eujJpp41KKaXcoudIlFJKuUULiVJKKbdoIVFKKeUWLSRKKaXcooVEKaWUW7SQKKWUcosWEqWUUm7RQqKUUsotWkiUUkq5RQuJUjYSkQY5Y1+0yXleK2c8iQ72JlPKedpFilI2E5E7sYZvjQMWAJuMMY/Ym0op52khUcoPiMjnwPlYHX1ebBzG51DK3+mhLaX8w7vAhVhjl2sRUQFF90iUspmIRGCN070S6A60MMYctjeVUs7TQqKUzURkGhBhjLlJRKYAlY0x/ezOpZSz9NCWUjYSkV5Yo2WOyJn1MNBGRAbal0op1+geiVJKKbfoHolSSim3aCFRSinlFi0kSiml3KKFRCmllFu0kCillHKLFhKllFJu0UKilFLKLVpIlFJKueX/Ab+NU/SbYhs7AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "plt.rcParams.update({'font.size': 12})\n", "\n", "Data1 = np.loadtxt(\"out10-numer.txt\")\n", "Data2 = np.loadtxt(\"out20-numer.txt\")\n", "\n", "def IDX4(i,j,k,Nxx_plus_2NGHOSTS0,Nxx_plus_2NGHOSTS1,Nxx_plus_2NGHOSTS2):\n", " return (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * ( (k) + Nxx_plus_2NGHOSTS2 * (0) ) )\n", "\n", "x1 = np.zeros(10)\n", "a1 = np.zeros(10)\n", "for i in range(10):\n", " x1[i] = Data1[IDX4(i+3,8,8,16,16,16),1]\n", " a1[i] = Data1[IDX4(i+3,8,8,16,16,16),0]\n", "x2 = np.zeros(20)\n", "a2 = np.zeros(20)\n", "for i in range(20):\n", " x2[i] = Data2[IDX4(i+3,13,13,26,26,26),1]\n", " a2[i] = Data2[IDX4(i+3,13,13,26,26,26),0]\n", "\n", "plt.figure()\n", "a = plt.plot(x1,a1,'.',label=\"dx\")\n", "b = plt.plot(x2,a2*(2**4),label=\"dx/2, times (20/10)^4\")\n", "plt.legend()\n", "plt.xlabel(\"x\")\n", "plt.ylabel(\"alpha\")\n", "plt.show()\n", "\n", "convergence_satisfactory = np.log2(a1/a2[0:-1:2])>3\n", "if not convergence_satisfactory.all():\n", " sys.exit(1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\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_UnitTest-GiRaFFE_NRPy-Metric_Face_Values.pdf](Tutorial-Start_to_Finish_UnitTest-GiRaFFE_NRPy-Metric_Face_Values.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": 8, "metadata": { "execution": { "iopub.execute_input": "2020-12-17T02:17:21.786537Z", "iopub.status.busy": "2020-12-17T02:17:21.786190Z", "iopub.status.idle": "2020-12-17T02:17:24.051215Z", "shell.execute_reply": "2020-12-17T02:17:24.051632Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Created Tutorial-Start_to_Finish_UnitTest-GiRaFFE_NRPy-\n", " Metric_Face_Values.tex, and compiled LaTeX file to PDF file Tutorial-\n", " Start_to_Finish_UnitTest-GiRaFFE_NRPy-Metric_Face_Values.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_UnitTest-GiRaFFE_NRPy-Metric_Face_Values\",location_of_template_file=os.path.join(\"..\"))" ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "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.9.1" } }, "nbformat": 4, "nbformat_minor": 2 }