{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "#Lets have matplotlib \"inline\"\n", "%matplotlib inline\n", "\n", "# Add line profiler\n", "%load_ext line_profiler\n", "\n", "#Import packages we need\n", "import numpy as np\n", "from matplotlib import animation, rc\n", "from matplotlib import pyplot as plt\n", "from matplotlib.colorbar import ColorbarBase\n", "\n", "from IPython.display import display, HTML\n", "\n", "from enum import Enum\n", "import subprocess\n", "import time\n", "import os\n", "import gc\n", "import datetime\n", "import importlib\n", "import logging\n", "import netCDF4\n", "import json\n", "\n", "import pycuda.driver as cuda\n", "import pycuda.compiler\n", "\n", "try:\n", " from StringIO import StringIO\n", "except ImportError:\n", " from io import StringIO\n", "\n", "from GPUSimulators import Common, IPythonMagic\n", "from GPUSimulators.EE2D_KP07_dimsplit import EE2D_KP07_dimsplit\n", "from GPUSimulators.helpers import InitialConditions, Visualization" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Console logger using level INFO\n", "File logger using level DEBUG to test_schemes.log\n", "Python version 3.9.19 | packaged by conda-forge | (main, Mar 20 2024, 12:50:21) \n", "[GCC 12.3.0]\n" ] } ], "source": [ "%setup_logging --out test_schemes.log logger" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Registering my_context in user workspace\n", "PyCUDA version 2024.1\n", "CUDA version (11, 8, 0)\n", "Driver version 12070\n", "Using device 0/1 'NVIDIA GeForce GTX 970' (0000:07:00.0) GPU\n", "Created context handle <96726128804080>\n", "Using CUDA cache dir /home/smyalygames/PycharmProjects/FiniteVolumeGPU/GPUSimulators/cuda_cache\n" ] } ], "source": [ "%cuda_context_handler --no_autotuning my_context " ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "#Set large figure sizes as default\n", "rc('figure', figsize=(16.0, 12.0))\n", "rc('animation', html='html5')\n", "rc('animation', bitrate=1800)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def plotVars(rho, rho_u, rho_v, E):\n", " plt.figure()\n", " plt.subplot(1,4,1)\n", " plt.imshow(rho, origin='lower')\n", " plt.colorbar(orientation='horizontal', pad=0.02, shrink=0.9)\n", "\n", " plt.subplot(1,4,2)\n", " plt.imshow(rho_u, origin='lower')\n", " plt.colorbar(orientation='horizontal', pad=0.02, shrink=0.9)\n", "\n", " plt.subplot(1,4,3)\n", " plt.imshow(rho_v, origin='lower')\n", " plt.colorbar(orientation='horizontal', pad=0.02, shrink=0.9)\n", "\n", " plt.subplot(1,4,4)\n", " plt.imshow(E, origin='lower')\n", " plt.colorbar(orientation='horizontal', pad=0.02, shrink=0.9)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def runSimulation(outfile, t_end, sim_args):\n", " with Common.Timer(\"construct\") as t:\n", " sim = EE2D_KP07_dimsplit(**sim_args)\n", " print(\"Constructed in \" + str(t.secs) + \" seconds\")\n", "\n", " #Create netcdf file and simulate\n", " with Common.DataDumper(outfile, mode='w', clobber=False) as outdata:\n", " outdata.ncfile.createDimension('time', None)\n", " outdata.ncfile.createDimension('x', sim.nx)\n", " outdata.ncfile.createDimension('y', sim.ny)\n", "\n", " #Create variables\n", " outdata.time_var = outdata.ncfile.createVariable('time', np.dtype('float32').char, 'time')\n", " outdata.x_var = outdata.ncfile.createVariable('x', np.dtype('float32').char, 'x')\n", " outdata.y_var = outdata.ncfile.createVariable('y', np.dtype('float32').char, 'y')\n", " outdata.rho_var = outdata.ncfile.createVariable('rho', np.dtype('float32').char, ('time', 'y', 'x'), zlib=True, least_significant_digit=3)\n", " outdata.rho_u_var = outdata.ncfile.createVariable('rho_u', np.dtype('float32').char, ('time', 'y', 'x'), zlib=True, least_significant_digit=3)\n", " outdata.rho_v_var = outdata.ncfile.createVariable('rho_v', np.dtype('float32').char, ('time', 'y', 'x'), zlib=True, least_significant_digit=3)\n", " outdata.E_var = outdata.ncfile.createVariable('E', np.dtype('float32').char, ('time', 'y', 'x'), zlib=True, least_significant_digit=3)\n", " \n", " #Create attributes\n", " def toJson(in_dict):\n", " out_dict = in_dict.copy()\n", "\n", " for key in out_dict:\n", " if isinstance(out_dict[key], np.ndarray):\n", " out_dict[key] = out_dict[key].tolist()\n", " else:\n", " try:\n", " json.dumps(out_dict[key])\n", " except:\n", " out_dict[key] = str(out_dict[key])\n", "\n", " return json.dumps(out_dict)\n", " outdata.ncfile.created = time.ctime(time.time())\n", " outdata.ncfile.sim_args = toJson(sim_args)\n", "\n", " outdata.x_var[:] = np.linspace(0, sim.nx*sim.dx, sim.nx)\n", " outdata.y_var[:] = np.linspace(0, sim.ny*sim.dy, sim.ny)\n", "\n", " progress_printer = Common.ProgressPrinter(n_save, print_every=10)\n", " print(\"Simulating to t={:f}. Estimated {:d} timesteps (dt={:f})\".format(t_end, int(t_end / sim.dt), sim.dt))\n", " for i in range(n_save+1):\n", " #Sanity check simulator\n", " try:\n", " sim.check()\n", " except AssertionError as e:\n", " print(\"Error after {:d} steps (t={:f}: {:s}\".format(sim.simSteps(), sim.simTime(), str(e)))\n", " return outdata.filename\n", "\n", " #Simulate\n", " if (i > 0):\n", " sim.simulate(t_end/n_save)\n", "\n", " #Save to file\n", " #rho = sim.u0[0].download(sim.stream)\n", " rho, rho_u, rho_v, E = sim.download()\n", " outdata.time_var[i] = sim.simTime()\n", " outdata.rho_var[i, :] = rho\n", " outdata.rho_u_var[i, :] = rho_u\n", " outdata.rho_v_var[i, :] = rho_v\n", " outdata.E_var[i, :] = E\n", "\n", " #Write progress to screen\n", " print_string = progress_printer.getPrintString(i)\n", " if (print_string):\n", " print(print_string)\n", "\n", " return outdata.filename " ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "class VisType(Enum):\n", " Schlieren = 0\n", " Density = 1\n", "\n", "def createAnimation(infile, vis_type, vmin, vmax, save_anim=False, cmap=plt.cm.coolwarm, fig_args=None):\n", " fig = plt.figure(**fig_args)\n", " \n", " with Common.DataDumper(infile, 'r') as indata:\n", " time = indata.ncfile.variables['time'][:]\n", " x = indata.ncfile.variables['x'][:]\n", " y = indata.ncfile.variables['y'][:]\n", " rho = indata.ncfile.variables['rho'][0]\n", " rho_u = indata.ncfile.variables['rho_u'][0]\n", " rho_v = indata.ncfile.variables['rho_v'][0]\n", " \n", " created = indata.ncfile.created\n", " sim_args = json.loads(indata.ncfile.sim_args)\n", " for key in sim_args:\n", " if isinstance(sim_args[key], list):\n", " sim_args[key] = \"[...]\"\n", " num_frames = len(time)\n", " print(\"{:s} created {:s} contains {:d} timesteps\".format(infile, created, num_frames))\n", " print(\"Simulator arguments: \\n\", sim_args)\n", "\n", " ax1 = fig.gca()\n", " extents = [0, x.max(), 0, y.max()]\n", " \n", " if (vis_type == VisType.Schlieren):\n", " im = ax1.imshow(Visualization.genColors(rho, rho_u, rho_v, cmap, vmax, vmin), origin='lower', extent=extents, cmap='gray', vmin=0.0, vmax=1.0)\n", " fig.suptitle(\"Schlieren / vorticity at t={:.2f}\".format(time[0]))\n", " elif (vis_type == VisType.Density):\n", " im = ax1.imshow(rho, origin='lower', extent=extents, cmap=cmap, vmin=vmin, vmax=vmax)\n", " fig.suptitle(\"Density at t={:.2f}\".format(time[0]))\n", " else:\n", " assert False, \"Wrong vis_type\"\n", "\n", " #Create colorbar\n", " from matplotlib.colors import Normalize\n", " from mpl_toolkits.axes_grid1 import make_axes_locatable\n", " divider = make_axes_locatable(ax1)\n", " ax2 = divider.append_axes(\"right\", size=0.1, pad=0.05)\n", " cb1 = ColorbarBase(ax2, cmap=cmap,\n", " norm=Normalize(vmin=vmin, vmax=vmax),\n", " orientation='vertical')\n", " \n", " #Label colorbar\n", " if (vis_type == VisType.Schlieren):\n", " cb1.set_label('Vorticity')\n", " elif (vis_type == VisType.Density):\n", " cb1.set_label('Density')\n", "\n", " progress_printer = Common.ProgressPrinter(num_frames, print_every=5)\n", " \n", " def animate(i):\n", " rho = indata.ncfile.variables['rho'][i]\n", " rho_u = indata.ncfile.variables['rho_u'][i]\n", " rho_v = indata.ncfile.variables['rho_v'][i]\n", " \n", " if (vis_type == VisType.Schlieren):\n", " im.set_data(Visualization.genColors(rho, rho_u, rho_v, cmap, vmax, vmin))\n", " fig.suptitle(\"Schlieren / vorticity at t={:.2f}\".format(time[i]))\n", " elif (vis_type == VisType.Density):\n", " im.set_data(rho) \n", " fig.suptitle(\"Density at t={:.2f}\".format(time[i]))\n", " \n", " #Write progress to screen\n", " print_string = progress_printer.getPrintString(i)\n", " if (print_string):\n", " print(print_string)\n", "\n", " anim = animation.FuncAnimation(fig, animate, interval=50, frames=range(num_frames))\n", " plt.close()\n", "\n", " if (save_anim):\n", " root, _ = os.path.splitext(infile)\n", " movie_outpath = os.path.abspath(root + \".mp4\")\n", " if (os.path.isfile(movie_outpath)):\n", " print(\"Reusing previously created movie \" + movie_outpath)\n", " else:\n", " print(\"Creating movie \" + movie_outpath)\n", " #from matplotlib.animation import FFMpegFileWriter\n", " #writer = FFMpegFileWriter(fps=25)\n", " from matplotlib.animation import FFMpegWriter\n", " writer = FFMpegWriter(fps=25)\n", " anim.save(movie_outpath, writer=writer)\n", " display(HTML(\"\"\"\n", "
\n", " \n", "
\n", " \"\"\".format(movie_outpath)))\n", " else:\n", " #plt.rcParams[\"animation.html\"] = \"html5\"\n", " plt.rcParams[\"animation.html\"] = \"jshtml\"\n", " display(anim)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Shock-bubble" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAABHsAAAE8CAYAAABQNL2dAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAABFAUlEQVR4nO3deXxU9b3/8ffMJJnsIQvZSAhhEyQVISwCIm5EEbcfVrG0ohZ7TQNapG7IbY20BWvv5UcrBS+tKF6Lci2g9HdRSVUWBSxLEFkElJAFMoSEbJA9c35/eJnbkQCZQDLJmdfz8ZjHQ77ne+Z8hvmeecj78T3fr8UwDEMAAAAAAAAwBau3CwAAAAAAAMDlQ9gDAAAAAABgIoQ9AAAAAAAAJkLYAwAAAAAAYCKEPQAAAAAAACZC2AMAAAAAAGAihD0AAAAAAAAmQtgDAAAAAABgIoQ9AAAAAAAAJkLYAwAAAAAAYCKEPQAAAAAAAO1g06ZNuuOOO5SYmCiLxaJ33333ouds3LhR6enpCgwMVO/evfXKK694fF3CHgAAAAAAgHZw5swZDR48WIsWLWpV/7y8PN12220aO3ascnNz9dxzz+nxxx/XqlWrPLquxTAMoy0FAwAAAAAAoHUsFovWrFmju++++7x9nnnmGa1du1YHDhxwtWVmZuqLL77Q1q1bW30tv0spFAAAAAAAoKupq6tTQ0ODx+cZhiGLxeLWZrfbZbfbL0tdW7duVUZGhlvbLbfcoldffVWNjY3y9/dv1fsQ9gAAAAAAAJ9RV1en1NRUORwOj88NDQ3V6dOn3dqef/55ZWdnX5baHA6H4uLi3Nri4uLU1NSk0tJSJSQktOp9CHsAAAAAAIDPaGhokMPhUEFBgcLDw1t9XlVVlXr27KnCwkK38y7XrJ6zvjtz6OzqO99tvxDCHgAAAAAA4HPCwsIUFhbW6v5nQ5fw8HCPQiJPxMfHnzPjqKSkRH5+foqOjm71+xD2AAAAAAAAn2MYhjzZs6oj9rcaNWqU/va3v7m1rV+/XsOGDWv1ej0SW68DAAAAAAAfdDbs8eTlqdOnT2v37t3avXu3pG+3Vt+9e7cKCgokSbNnz9bUqVNd/TMzM5Wfn69Zs2bpwIEDWrZsmV599VU9+eSTHl2XmT0AAAAAAMDndMTMnh07duiGG25w/XnWrFmSpAcffFCvv/66iouLXcGPJKWmpmrdunV64okn9Mc//lGJiYn6wx/+oHvuucej61qMjpiHBAAAAAAA0AlUVVUpIiJCpaWlHi/QHBMTo8rKynZbs+dyabfHuBYvXqzU1FQFBgYqPT1dmzdvbq9LAQAAAAAAeKQjHuPylnYJe1auXKmZM2dqzpw5ys3N1dixYzVhwgS3qUkAAAAAAADeYuawp10e4xo5cqSGDh2qJUuWuNoGDhyou+++W/Pnz7/guU6nU8ePH1dYWJhHe8gDAAAAAICuzzAMVVdXKzExUVbr5Z+jcvYxrhMnTnj8GFdcXFyXeIzrsi/Q3NDQoJ07d+rZZ591a8/IyNCWLVvO6V9fX6/6+nrXn48dO6Yrr7zycpcFAAAAAAC6kMLCQiUlJbXb+3fGrdcvl8se9pSWlqq5uVlxcXFu7XFxcXI4HOf0nz9/vl544YVz2gsLCy97UnbkyBE9++yz+vDDD11tUVFRmjt3rh544IHLei0AAAAAAOC5qqoqJScnKywsrF2vQ9jTBt99BMswjBYfy5o9e7Zr6zHpf7/U8PDwyx72hIWFyd/f363NarUqKCio00/BAgAAAADAl7T30i6EPR6IiYmRzWY7ZxZPSUnJObN9JMlut8tut1/uMgAAAAAAAM7LzGHPZV/pKCAgQOnp6crJyXFrz8nJ0ejRoy/35QAAAAAAADxm5t242uUxrlmzZumBBx7QsGHDNGrUKC1dulQFBQXKzMxsj8sBAAAAAAB4xMwze9ol7Jk8ebLKyso0d+5cFRcXKy0tTevWrVNKSkp7XA4AAAAAAMAjhD1tkJWVpaysrPZ6ewAAAAAAgDYj7AEAAAAAADARM4c9l32BZgAAAAAAAHgPM3sAAAAAAIDPMfPMHsIeAAAAAADgk7pSgOMJwh4AAAAAAOBzmNkDAAAAAABgIoQ9AAAAAAAAJkLYAwAAAAAAYCKEPQAAAAAAACZC2AMAAAAAAGAihD0AAAAAAAAmQtgDAAAAAABgIoQ9AAAAAAAAJkLYAwAAAAAAYCJmDnus3i4AAAAAAAAAlw8zewAAAAAAgM8x88wewh4AAAAAAOBzCHsAAAAAAABMhLAHAAAAAADARAh7AAAAAAAATISwBwAAAAAAwEQIewAAAAAAAEyEsAcAAAAAAMBkulKA4wnCHgAAAAAA4HOY2QMAAAAAAGAihD0AAAAAAAAmQtgDAAAAAABgImYOe6zeLgAAAAAAAACXDzN7AAAAAACAzzHzzB7CHgAAAAAA4HMIewAAAAAAAEyEsAcAAAAAAMBECHsAAAAAAABMxMxhD7txAQAAAAAAn3M27PHk1RaLFy9WamqqAgMDlZ6ers2bN1+w/1/+8hcNHjxYwcHBSkhI0MMPP6yysjKPrknYAwAAAAAAfE5HhD0rV67UzJkzNWfOHOXm5mrs2LGaMGGCCgoKWuz/6aefaurUqZo2bZr27dund955R9u3b9cjjzzi0XUJewAAAAAAgM/piLBnwYIFmjZtmh555BENHDhQCxcuVHJyspYsWdJi/23btqlXr156/PHHlZqaqmuvvVaPPvqoduzY4dF1CXsAAAAAAIDPaWvYU1VV5faqr69v8f0bGhq0c+dOZWRkuLVnZGRoy5YtLZ4zevRoFRUVad26dTIMQydOnNBf//pXTZw40aPPRtgDAAAAAAB8TlvDnuTkZEVERLhe8+fPb/H9S0tL1dzcrLi4OLf2uLg4ORyOFs8ZPXq0/vKXv2jy5MkKCAhQfHy8unXrppdfftmjz0bYAwAAAAAAfE5bw57CwkJVVla6XrNnz77gdSwWyznX/W7bWfv379fjjz+uX/7yl9q5c6c++OAD5eXlKTMz06PPxtbrAAAAAADA57R16/Xw8HCFh4dftH9MTIxsNts5s3hKSkrOme1z1vz58zVmzBg99dRTkqSrrrpKISEhGjt2rH79618rISGhVbUyswcAAAAAAPic9l6gOSAgQOnp6crJyXFrz8nJ0ejRo1s8p6amRlare1Rjs9lc9baWR2HP/PnzNXz4cIWFhSk2NlZ33323Dh486NbHMAxlZ2crMTFRQUFBuv7667Vv3z5PLgMAAAAAANDlzZo1S3/+85+1bNkyHThwQE888YQKCgpcj2XNnj1bU6dOdfW/4447tHr1ai1ZskRHjhzRZ599pscff1wjRoxQYmJiq6/rUdizceNGTZ8+Xdu2bVNOTo6ampqUkZGhM2fOuPq89NJLWrBggRYtWqTt27crPj5e48ePV3V1tSeXAgAAAAAAaFftue26JE2ePFkLFy7U3LlzdfXVV2vTpk1at26dUlJSJEnFxcUqKChw9X/ooYdcmUpaWpruvfdeXXHFFVq9erVH1/VozZ4PPvjA7c+vvfaaYmNjtXPnTl133XUyDEMLFy7UnDlzNGnSJEnS8uXLFRcXpxUrVujRRx/1qDgAAAAAAID20NY1ezyVlZWlrKysFo+9/vrr57Q99thjeuyxx9p0rbMuac2eyspKSVJUVJQkKS8vTw6Hw20PebvdrnHjxp13D/n6+vpz9qgHAAAAAABoT+29Zo83tTnsMQxDs2bN0rXXXqu0tDRJcq0w7cke8vPnz3fbnz45ObmtJQEAAAAAALQKYU8LZsyYoT179uitt94655gne8jPnj3bbX/6wsLCtpYEAAAAAADQKmYOezxas+esxx57TGvXrtWmTZuUlJTkao+Pj5f07Qyff977/UJ7yNvtdtnt9raUAQAAAAAA0CYdtWaPN3g0s8cwDM2YMUOrV6/Wxx9/rNTUVLfjqampio+Pd9tDvqGhQRs3bjzvHvIAAAAAAAAdjZk9/2P69OlasWKF3nvvPYWFhbnW4YmIiFBQUJAsFotmzpypefPmqV+/furXr5/mzZun4OBgTZkypV0+AAAAAAAAgKfMPLPHo7BnyZIlkqTrr7/erf21117TQw89JEl6+umnVVtbq6ysLJWXl2vkyJFav369wsLCLkvBAAAAAAAAl4qw53+05oNZLBZlZ2crOzu7rTUBAAAAAAC0K8IeoI1OnTqlQ4cOqaampt2uYbFYlJSUpN69e8tms7XbdQAAAAAA5kHYA7TR/v379eKLL+ro0aPtdg2LxaIf/vCHmjFjhkJDQ9vtOgAAAAAA8yDsAb6jublZjY2NF+1XVlamgwcP6uuvv263WqxWq4qKinTmzBn5+V14SFutVvn7+8tisbRbPQAAAAAAeBNhDzzmdDq1detWbdiw4aKBz+HDh1VeXt6u9RiGoe3bt2vBggUKDAy8YN/U1FTddtttio2NbdeaAAAAAACdGzN7gH/idDr1+eef6//+3/+rM2fOXLRvU1NTu9ZjGIZ27typL7744qJ9r7/+eg0bNoywBwAAAAB8HGEPfJbT6dSJEydUVlbmGthNTU0qKipSbW2t6uvrL3h+WFiYEhISFBAQ0K51lpWV6cSJE3I6nRfsV15erkOHDrndpIGBgUpISGC9HwAAAADwIYQ98Fn19fVas2aN/vrXv7qCFMMwVFRUdNGgR5KuvPJKzZgxQ8nJye1Wo9Pp1Jo1a7Rs2bKLzjQ6dOiQ5s+fr5CQEFdb79699bOf/UyDBw9utxoBAAAAAJ0LYQ98xncHb2Njo77++mtt2rRJzc3NbscsFstFFzqOiorSiBEj1L9//8te61lOp1N79+5t1cLLlZWV2rFjh1tbeXm5ysvLz/nsLOIMAAAAAOZF2AOfUFNTo61bt+rw4cOutrq6Ou3Zs8dtUFssFg0ePFjp6eny9/e/4HsOGjRI3bp1a6+SXfWkpaXp4YcfVm1t7QX7Hjt2TJ999plOnTrlaisrK9PatWv11VdfudpCQ0M1evRo9e7du93qBgAAAAB4D2EPfEJ1dbX+67/+S3/9619dbYZhqLa21m0tHJvNpuuuu05PPfWUgoODL/ieAQEBF+1zqSwWi0aPHq2rr776ojff5s2bdfjwYbewx+Fw6M9//rNbcJWUlKRf//rXhD0AAAAAYGJdKcDxBGEPdPr0aVVVVam4uFglJSVuQYjFYlFERIRiYmJcjzXZbDbFx8crKiqq3YOc1rLb7bLb7Rft1717dyUlJbmt7dPQ0KCKigpVV1e72gIDA1VcXKzCwkIFBwerW7dustls7VI7AAAAAKDjMbMHprZ9+3atWLFCJ06cOGf7crvdrjvvvFMTJkyQ1WqVJFmtVvXv379V4Upn07dvXz355JOqrKx0teXn5+v111/X/v37XW0VFRVavny5PvroI40cOVIPPvigoqOjvVEyAAAAAKAdEPbA1PLy8vTuu++qtLT0nGN+fn66+uqrde+995piZktsbKxuueUWt7Y9e/bov//7v93aampqtGXLFklSc3Ozvv/97xP2AAAAAICJEPbAFAzD0OHDh3Xw4EHXzlqGYSg3N1cNDQ2y2+1KS0tz2yY9KChIffr0MfXOVOHh4Ro7dqwiIyNdbbW1tfryyy91/PhxHTt2TB9++KG6d+/uOh4TE6OrrrpK4eHh3igZAAAAAHCJCHtgCk1NTfroo4/0hz/8QXV1da7206dP6/Tp04qOjtaPfvQj3XXXXa5jVqtVkZGRrke4zCgxMVHTp09328nr2LFj+s1vfqPjx49rz549+tWvfuU2s2nkyJHKzs4m7AEAAACALoqwB11ac3OzamtrVVdXJ4fDofz8fLdgw9/fXyEhIYqMjFSPHj2UmprqxWo7XkBAgOLj493a/P391b17d0VERKixsVFFRUVuN3ZiYqJKS0tVVVXV6sWhAQAAAADoCIQ9PqCgoECrVq1Sfn6+du3apcbGRrfjQ4cO1Z133qm4uDhdddVVXqqyc4mIiND999+vYcOGaffu3Xr33XdVXl7uOn706FEtXrxY8fHxGj9+vMaPHy8/P24nAAAAAOgqmNmDLq24uFgrV65Ubm6unE7nOQP0yiuv1LRp09S9e3dTr83jibCwMN16660yDEOrV6/W3//+d7ew5/jx41q5cqUCAwMVERGhG264gbAHAAAAALoQwh50Oc3NzSooKJDD4dCePXtUXV0tp9OpxMREJSUlua0/069fPwUGBpp6XZ62OPv30b17dw0bNsxt4eqqqiodOXJEjY2Nys/P17Zt29StWzelpqYqIiLCWyUDAAAAAFqJsAddTm1trVatWqWVK1eqqqpKhYWFstlsysjI0COPPKLAwEBX3+joaIWEhHix2s7t6quv1ty5c9XQ0OBq27Vrl37729/qm2++0fvvv6/c3Fz1799fTz75pIYPH+7FagEAAAAArUHYgy7DMAw1Nzerrq7OtUaPYRiy2Wyy2+3q0aOHhgwZoqCgIG+X2mV069ZN3bp1c2urqalRRESE/Pz8VFJSIofDocbGRpWXl6uxsVFWq9Vt9hQAAAAAoHMh7EGX8fXXX+ujjz6Sw+FwBT2JiYkaP368evTooWuvvZa1ZS6DpKQkPfDAAyoqKtK2bdu0bds2lZaW6p133lFubq6uvvpqjRs3zm0GFQAAAACg8yDsQZfx1Vdf6eWXX9bRo0fV0NAgwzCUnJysadOmaejQofL395e/v7+3y+zyevbsqUcffVS1tbVasGCBduzYoZKSEr355pvy8/PTww8/rOHDhxP2AAAAAEAnRdiDTq2pqUknT55UdXW1CgsLVVVVpbq6OkVHRysyMlKpqamKiopScHCwt0s1DavVqsDAQNlsNsXFxalv376qrq5WSUmJzpw5o5KSEn399deKiYlR9+7dFRYW5u2SAQAAAAD/hLAHnVpFRYWWL1+uDRs26MSJEyorK1NgYKDuuece3XnnnYqKilJSUpK3yzQlPz8/TZgwQVdccYWOHDmiJUuW6IsvvtBnn32mkpIS9ejRQ4888ojGjRvn7VIBAAAAAP+EsAedWl1dnXbv3q0PP/zQ1RYREaGBAwfqlltuYUv1dmSxWNSnTx/16dNHe/fu1cqVKyVJRUVFKioqUlJSkiZOnOjlKgEAAAAA30XYg07p0KFD2rVrlxwOh44ePSpJSk5O1ogRI9S9e3cNHDhQFovFu0X6kG7duumWW25RcnKyDh06pNzcXJ05c0abNm1SfX29evbsqeHDhys0NNTbpQIAAACAzyPsQadjGIY+//xzzZs3T6WlpTp9+rQkKS0tTc8884x69eqlkJAQwp4OFB8fr3/5l39RfX293nzzTR08eFCVlZVasWKFVq9erQkTJqhPnz6EPQAAAADQCRD2oFOqqanRyZMnVVZW5mqz2+2uRYHRsfz8/BQZGSnDMBQWFiar1Sqn06nKykpVVlaqoqJCzc3N3i4TAAAAAGByhD0AAAAAAMDnMLMHAAAAAADAZLpSgOMJwp4u5vTp0zpw4IDKysq0f/9+NTQ0yG63a8CAAYqPj9eQIUMUHBzs7TJ9msViUXJysm688UaVlJTo4MGDcjgccjgc2rRpk44cOaK+ffsqJSXF26UCAAAAgM9iZg86jeLiYi1atEjbtm1TZWWlzpw5o+7du+vBBx/UxIkTFRoaqujoaG+X6fPGjBmjK664QseOHdOLL76o999/X3v37lV2draio6P1+OOP60c/+pGsVqu3SwUAAAAAn0TYA69rbGxUY2OjKioqlJ+fr0OHDsnPz0+BgYEKCwtTcnKy+vfv7+0y8T8iIiIUERGhoKAgxcbGKiQkRI2NjcrPz9epU6d04sQJnTlzRv7+/goICCD0AQAAAIAORtgDr2pubtamTZu0YcMGFRcX6+jRo5KkoUOH6pZbblF8fLwGDRrk3SLRorCwME2aNEkDBgzQ/v379f/+3/9TTU2N1q9fr9LSUqWmpurOO+9UYmKit0sFAAAAAJ9C2AOvam5u1tatW/X73/9etbW1am5ulsVi0VVXXaWsrCzFxMTIZrN5u0y0ICwsTBMnTtSECRP03nvvafPmzSovL9cnn3yijRs3avTo0Ro5ciRhDwAAAAB0MMIeeJ3T6VRjY6OampokfbsIsM1mk7+/v/z8+Bo7M5vNJpvN5vY9NTc3q7m5WU1NTV3qBwMAAAAAzMLMYQ8LhQAAAAAAAJ9zNuzx5NUWixcvVmpqqgIDA5Wenq7NmzdfsH99fb3mzJmjlJQU2e129enTR8uWLfPomkwJ6eScTqeam5vdBpXVapXVapXFYpHFYvFidfDE2dlYVqvV7cfi7Hd89jsFAAAAALS/jpjZs3LlSs2cOVOLFy/WmDFj9B//8R+aMGGC9u/fr549e7Z4zn333acTJ07o1VdfVd++fVVSUuJ6yqe1CHs6sePHj2vDhg1yOBz6/PPP1dTUpOjoaF133XXq2bOnRo8ercDAQG+XiVbq3bu3pk6dquLiYm3dulV79uzR8ePHtXLlSm3btk1Dhw7ViBEjeCwPAAAAADpAR4Q9CxYs0LRp0/TII49IkhYuXKgPP/xQS5Ys0fz588/p/8EHH2jjxo06cuSIoqKiJEm9evXy+LqX9BjX/PnzZbFYNHPmTFebYRjKzs5WYmKigoKCdP3112vfvn2XchmfdfToUS1evFgvvPCCPv74YzU1Nal79+56+OGH9cILL+iOO+5QUFCQt8tEKw0cOFBPPPGEnnvuOY0aNUoWi0WFhYV65ZVX9MILL+jvf/+7GhsbvV0mAAAAAPiEtj7GVVVV5faqr69v8f0bGhq0c+dOZWRkuLVnZGRoy5YtLZ6zdu1aDRs2TC+99JJ69Oih/v3768knn1Rtba1Hn63NYc/27du1dOlSXXXVVW7tL730khYsWKBFixZp+/btio+P1/jx41VdXd3WS/mspqYmVVdXuw0ePz8/hYSEKCIiQkFBQTz204X4+fkpLCxM4eHhCggIkPTtQs2nT59WZWWl6urqvFwhAAAAAOBikpOTFRER4Xq1NENHkkpLS9Xc3Ky4uDi39ri4ODkcjhbPOXLkiD799FPt3btXa9as0cKFC/XXv/5V06dP96jGNoU9p0+f1g9/+EP96U9/UmRkpKvdMAwtXLhQc+bM0aRJk5SWlqbly5erpqZGK1asaMulAAAAAAAALru2zuwpLCxUZWWl6zV79uwLXue7kzQMwzjvxA2n0ymLxaK//OUvGjFihG677TYtWLBAr7/+ukeze9oU9kyfPl0TJ07UzTff7Nael5cnh8PhNkXJbrdr3Lhx552iVF9ff84UKAAAAAAAgPbU1rAnPDzc7WW321t8/5iYGNlstnNm8ZSUlJwz2+eshIQE9ejRQxEREa62gQMHyjAMFRUVtfqzebwS7Ntvv61du3Zp+/bt5xw7+wFamqKUn5/f4vvNnz9fL7zwgqdlmFZjY6O+/PJLffPNN9q/f78qKytltVo1YMAADRw4UL169VJsbKy3y8Ql8Pf315AhQ1wrrO/evVtVVVU6cOCAVq1apZiYGA0ZMuS8Nz8AAAAA4NK19wLNAQEBSk9PV05Ojv7P//k/rvacnBzdddddLZ4zZswYvfPOOzp9+rRCQ0MlSYcOHZLValVSUlKrr+1R2FNYWKif/exnWr9+/QV3gfJkitLs2bM1a9Ys15+rqqqUnJzsSVmmUltbq9WrV+s///M/VVtbq4qKCvn5+SkjI0PTp09XaGio26Nz6HqCgoJ099136+abb9b27dv1y1/+UuXl5froo4+0fft29e/fX9nZ2YQ9AAAAANCOOmI3rlmzZumBBx7QsGHDNGrUKC1dulQFBQXKzMyU9G0mcuzYMb3xxhuSpClTpuhXv/qVa2Om0tJSPfXUU/rxj3/s0QZNHoU9O3fuVElJidLT011tzc3N2rRpkxYtWqSDBw9K+naGT0JCgqvPhaYo2e3280558kVOp1OnTp1SUVGRnE6npG/TwPDwcCUlJbHVuglYLBZ169ZN3bp109GjR12LNVdXV6u6ulphYWEs1gwAAAAA7awjwp7JkyerrKxMc+fOVXFxsdLS0rRu3TqlpKRIkoqLi1VQUODqHxoaqpycHD322GMaNmyYoqOjdd999+nXv/61R9f1KOy56aab9OWXX7q1PfzwwxowYICeeeYZ9e7dW/Hx8crJydGQIUMkfbvV2MaNG/Xb3/7Wo8IAAAAAAADaS0eEPZKUlZWlrKysFo+9/vrr57QNGDBAOTk5bbrWWR6FPWFhYUpLS3NrCwkJUXR0tKt95syZmjdvnvr166d+/fpp3rx5Cg4O1pQpUy6pUAAAAAAAgMulo8Ieb/B4geaLefrpp1VbW6usrCyVl5dr5MiRWr9+vcLCwi73pQAAAAAAANqsKwU4nrjksGfDhg1uf7ZYLMrOzlZ2dvalvjUAAAAAAEC7YGYPAAAAAACAiRD2AAAAAAAAmAhhDwAAAAAAgIkQ9qDdVVVVyeFwqKysTOXl5TIMQ8HBwUpISFBYWJi6d+8ui8Xi7TJxmYWEhKhPnz6qq6tTaWmpTp48qbq6OuXn52vv3r2KiopSXFycbDabt0sFAAAAAHQRhD2dxBdffKHFixerqKhIeXl5MgxDffv21WOPPaa+ffsqJSVF/v7+3i4Tl1nv3r311FNPqby8XP/1X/+lN998UyUlJVq8eLFWrlypW2+9VT/5yU8UERHh7VIBAAAAwFSY2YN2V1JSoi1btqigoMDV1q1bN40YMUJXXXWVFytDe4qIiNDw4cPV2NionTt3ymq1qqamRrt375Yk9ezZU42Njd4tEgAAAABMiLAHAAAAAADARAh7AAAAAAAATISwBwAAAAAAwEQIewAAAAAAAEyEsAcAAAAAAMBECHsAAAAAAABMhLAHAAAAAADARAh7AAAAAAAATISwBwAAAAAAwEQIewAAAAAAAEzEzGGP1dsFAAAAAAAA4PJhZg8AAAAAAPA5Zp7ZQ9gDAAAAAAB8UlcKcDxB2AMAAAAAAHwOM3sAAAAAAABMhLAHAAAAAADARAh7AAAAAAAATISwBwAAAAAAwEQIewAAAAAAAEyEsAcAAAAAAMBECHsAAAAAAABMhLAHAAAAAADARAh7AAAAAAAATISwB+0uLi5OY8eO1bFjx/T111+rqKhIp06d0tatW3Xy5EmlpqaqV69eslqt3i4Vl1FFRYUOHjyoU6dO6fDhw3I6nQoJCdEVV1yh6OhoDRo0SP7+/t4uEwAAAADQhRD2dBKDBw9Wdna2Tp06pYULF+rtt9/WN998oxdffFHh4eH6yU9+op/85Cey2+3eLhWX0ZEjR/S73/1Oe/fu1alTp9TQ0KAePXropz/9qa699lpFREQoNDTU22UCAAAAgOkwswftLiwsTGFhYaqoqFC3bt1ksVhUW1uro0ePKiAgQCdPnuxSAwutU1NTo7y8PB08eNDVFhAQoJ49e2rAgAFerAwAAAAAzI2wBwAAAAAAwEQIewAAAAAAAEyEsAcAAAAAAMBECHsAAAAAAABMhLAHAAAAAADARAh7AAAAAAAATKYrBTieIOwBAAAAAAA+h5k96DBWq1XR0dHq2bOnamtrderUKUlSZWWlCgoKFBYWpsjISAUGBnq5UrSVYRiqqKhQdXW1Tpw4oYaGBklSeHi4IiIi1KNHDwUFBXm5SgAAAAAwN8IedJigoCDdc889GjJkiPbt26dXX31VhYWF+vDDD5Wfn6/U1FT9+Mc/VlpamrdLRRvV1tbq3Xff1fvvv6+TJ0+qqKhIfn5+uummm3TvvfcqOjpaV1xxhbfLBAAAAABTM3PYY/X0hGPHjulHP/qRoqOjFRwcrKuvvlo7d+50HTcMQ9nZ2UpMTFRQUJCuv/567du377IWbWb+/v66+uqrNWnSJI0bN04RERFyOp366quv9O677yonJ0clJSXeLhOXoLGxUbt379aqVau0YcMGVVRUyGKxaMCAAbr77ruVkZGh2NhYb5cJAAAAAKZ2Nuzx5NUWixcvVmpqqgIDA5Wenq7Nmze36rzPPvtMfn5+uvrqqz2+pkdhT3l5ucaMGSN/f3+9//772r9/v/793/9d3bp1c/V56aWXtGDBAi1atEjbt29XfHy8xo8fr+rqao+LAwAAAAAA6KpWrlypmTNnas6cOcrNzdXYsWM1YcIEFRQUXPC8yspKTZ06VTfddFObrutR2PPb3/5WycnJeu211zRixAj16tVLN910k/r06SPp21Rs4cKFmjNnjiZNmqS0tDQtX75cNTU1WrFiRZsKBAAAAAAAuNw6YmbPggULNG3aND3yyCMaOHCgFi5cqOTkZC1ZsuSC5z366KOaMmWKRo0a1abP5lHYs3btWg0bNkz33nuvYmNjNWTIEP3pT39yHc/Ly5PD4VBGRoarzW63a9y4cdqyZUuL71lfX6+qqiq3F77l5+en8PBwdevWzbUgc1NTk06fPq2KigrV1tZ2qWcGfV1jY6OqqqpUWVmp+vp6SZLNZnMtus2izAAAAADQcdoa9nw3wzj777vvamho0M6dO90yEknKyMg4b0YiSa+99pq++eYbPf/8823+bB4t0HzkyBEtWbJEs2bN0nPPPad//OMfevzxx2W32zV16lQ5HA5JUlxcnNt5cXFxys/Pb/E958+frxdeeKGN5Ztbr169NGPGDDkcDuXk5OjDDz9USUmJli1bppycHI0ZM0a33367QkNDvV0qWuHAgQNatWqViouL9fnnn8swDPXq1Uv33nuvUlJSNGTIEPn7+3u7TAAAAADwCW1doDk5Odmt/fnnn1d2dvY5/UtLS9Xc3NxiRnI2P/muw4cP69lnn9XmzZvl59f2PbU8OtPpdGrYsGGaN2+eJLl2jFqyZImmTp3q6mexWNzOMwzjnLazZs+erVmzZrn+XFVVdc5fnK9KTEzUfffdp7q6OlVUVOijjz7SqVOntHbtWlmtVjU1NZ2TEKLzysvL01/+8hfl5eW5flTi4+P1/e9/X0OHDpXFYjnvfQIAAAAAuLzaGvYUFhYqPDzc1W632y94XmszkubmZk2ZMkUvvPCC+vfv3+q6WuJR2JOQkKArr7zSrW3gwIFatWqVJCk+Pl6S5HA4lJCQ4OpTUlJyTpJ1lt1uv+hfjC+zWCyy2WxuA8EwDDmdzktaDRwdzzAMNTc3y+l0utosFousVqusVo83xgMAAAAAXIK2hj3h4eFuYc/5xMTEyGaznTOL53wZSXV1tXbs2KHc3FzNmDFDklz/9vfz89P69et14403tqpWj/6FOWbMGB08eNCt7dChQ0pJSZEkpaamKj4+Xjk5Oa7jDQ0N2rhxo0aPHu3JpQAAAAAAANpNey/QHBAQoPT0dLeMRJJycnJazEjCw8P15Zdfavfu3a5XZmamrrjiCu3evVsjR45s9bU9mtnzxBNPaPTo0Zo3b57uu+8+/eMf/9DSpUu1dOlSSd/OUpg5c6bmzZunfv36qV+/fpo3b56Cg4M1ZcoUTy6F77BarQoICFBzc7OampokfTvFq7GxUY2NjbLZbMwO6aSamprkdDpd35v07cLMNptN/v7+PLoFAAAAAF7Q1pk9npg1a5YeeOABDRs2TKNGjdLSpUtVUFCgzMxMSd8ubXPs2DG98cYbslqtSktLczs/NjZWgYGB57RfjEdhz/Dhw7VmzRrNnj1bc+fOVWpqqhYuXKgf/vCHrj5PP/20amtrlZWVpfLyco0cOVLr169XWFiYR4Xhf9lsNo0ZM0ZNTU06fvy41q9fr/z8fO3evVt/+MMflJCQoJtuuumcR+zgfdXV1froo4+0d+9effXVV6qsrJTdbtd1112nYcOGKTU19byPOAIAAAAA2k9HhD2TJ09WWVmZ5s6dq+LiYqWlpWndunWuJ6SKi4tVUFDg8ftejMdLO99+++26/fbbz3vcYrEoOzu7xZWo0TY2m01jx47VNddco3379unQoUPKz89Xbm6u9u3bpx49erS4nhK8r7q6WmvWrNE777zjmokVFhamjIwMZWZmyt/fXwEBAd4uEwAAAAB8TkeEPZKUlZWlrKysFo+9/vrrFzy3rflK2/fxQofy8/OTn5+fIiIilJqaqhMnTqiiokInT55UdXW1CgoK9NVXXyk0NFRxcXFs4e1lZ7+bY8eO6eTJk6qtrVVoaKiSk5MVGRmpuLg4BQcH8+gdAAAAAHhJR4U93kDY08UkJCRoxowZ+sEPfqB169Zp2bJlqqio0GuvvaYPPvhAo0ePVmZmpmtnNHQ8wzD06aefatmyZTp58qQOHz4sSbrqqquUmZmp5ORk9enTh6AHAAAAALyIsAedRmhoqNLT02UYhvLy8hQQEKDq6mrt3btXe/fuVUhIiGpra71dps87duyYNm3apLKyMldbbGysrr32WqWmpnqxMgAAAACAZO6wh6kFAAAAAAAAJsLMHgAAAAAA4HPMPLOHsKcLCwkJUVxcnGw2m6qqqlRXV6e6ujqVlJQoODhYYWFhCg4O9naZPqOpqcn1PVRWVsrpdMpqtSo8PFxBQUGKjIyUnx+3HAAAAAB0Fl0pwPEE//LsoiwWi0aOHKns7GwVFxfrrbfe0rZt27Rv3z69+OKLio2N1T333KPx48fLYrF4u1yf4HA49J//+Z/at2+fvv76a505c0YRERGaMmWKrrnmGtdOXAAAAAAA72NmDzqlfv36qV+/fioqKtKWLVu0bds2FRYWqrCwUBERERo0aJBuvvlmwp4OUlFRoZycHH3yySeutrOLMt9///1erAwAAAAA8F2EPejUgoKCNGTIEFVVVamkpET79+9XY2Oj9u3bp3Xr1ik6OlqDBg1SeHi4t0s1HafTqSNHjujrr7/WkSNHXLtv9ezZU1dccYUSEhKUmJjo5SoBAAAAAN9F2INOrVu3bnrooYd0zz33KCcnR/PmzdPx48e1evVqffzxx0pPT9e//uu/Ki0tzdulmk5zc7M++OADvfLKK6qurtbJkydlsVg0evRo/fznP1dMTIxiYmK8XSYAAAAA4DsIe9Cp2Ww2xcXFKS4uTl999ZUiIiJUXl6uyspKlZaWKiYmRhUVFTp9+rQCAgIUEBDg7ZK7PKfT6VoQ2+Fw6JtvvlFDQ4PsdrtCQ0MVGxur3r17KyoqytulAgAAAABaQNiDLmPAgAH62c9+JofDofXr1+uzzz5TYWGhli5dqh49emjcuHG66aab5O/v7+1Su7T8/Hy9++67Kiws1D/+8Q81NjYqNjZWd911l/r27avBgwezExoAAAAAdGKEPegy+vTpo169eqmyslIlJSXasmWLjh8/rhUrVshut8vPz0/jxo0j7LlERUVFevPNN7Vnzx45nU45nU5FR0frnnvu0Q033CCr1Sqr1ertMgEAAAAA50HYgy7DYrHIz89PdrtdvXr10vDhw1VZWamCggI1NDSoqKhIO3fulN1ud/WPiYlRUlKS/PwYDi0pLy9XYWGh6uvrXW0HDhxQVVWVmpublZCQoISEBPXr10+RkZH8PQIAAABAF0DYgy4nKChI3//+9zV27Fjt3r1bCxYs0KFDh7R+/Xrt3btXNptN0rdhz8SJEzVjxgx169bNu0V3Urt379bvf/97ORwOV1t1dbWOHz8uPz8/3XLLLXrooYcUGRmplJQUL1YKAAAAAGgtwh50OTabTSkpKUpJSVFTU5PCwsJktVrlcDh0/Phxt74DBgxQXV2dmpubZbVaZbFYvFR15+J0OmUYhkpLS7Vr1y4VFha6HbfZbLLb7UpJSdGIESMUGBjopUoBAAAAAJ4i7EGXlpCQoB/84AcaM2aMduzYoc8//1xNTU2u4/v379ef/vQnxcXF6brrrtOAAQO8WG3nUF1drY0bN+rw4cP64osvVF1d7Xa8R48euvHGGxUfH69rrrmGR7cAAAAAAJ0G/0L1AT179lRmZqZqa2v18ssva9euXW5hT25urvbv36/4+HhFRkYS9kiqrKzUypUr9d5776mpqUl1dXVux1NSUvTTn/5UaWlproWvAQAAAABdBzN70KXZbDYFBwfL399f8fHx6t27t1t4UVVVpbKyMlVWVqqoqEjffPON65jValVUVJQiIiK8UXqHaGhoUGlpqWpra11tx48f18mTJ1VdXa2QkBD17NnTLdDp2bOnIiMjFRYW5o2SAQAAAACXiLAHpuDn56ebb75ZPXv2VHNzs6RvB+vf//53LV++XFVVVXrzzTf1ySefuM4JCgrSlClTdMcdd5h2K/Hjx49r6dKl2rt3r6uttrZW+/btkyQNHjxY06ZNU/fu3V3Ho6OjlZiY2OG1AgAAAAAuD8IemILFYlHfvn3Vt29ft/bS0lK99dZbqqqq0q5du7Rr1y7XsdDQUF1zzTVdalB7qqqqSlu2bNHGjRtbPJ6QkOAKyQAAAAAA5kDYA1Pr3bu3Jk2apJKSEuXm5io/P991rKmpSbt27dJbb73l2q7darVqwIABGjRoUJdbq+bEiRPatWuXKioqXG0FBQU6ceKEW7/g4GANHTpUycnJGjFihEJCQjq4UgAAAABAeyLsgamNGDFCAwYMUHFxsX71q1+5hT319fX629/+po8//ti1Jbufn5/+5V/+Rf369etyYc/hw4f1u9/9TocOHXK1NTY2uoU/ktStWzdNnTpVEyZMUHBwsKnXLAIAAAAAX0TYA1MLDg5WcHCwLBaLYmNjFRMT4xrEhmGopqZGlZWVrv5+fn5yOBwqLS296IwXf39/hYaGtvt6P3V1daqpqbnozVdaWqri4mIdO3bM1Waz2RQSEuK22HL37t2VkJCgpKSkdqsZAAAAAOA9hD3wCeHh4Zo8ebLS09NdbXV1dVq7dq0+/vhjOZ1OSVJzc7M2bdqk6upq+fv7X/A909LSNHnyZMXFxbVb3YZhaMuWLXrvvffcdtRqyfHjx895ZCshIUH333+/+vfv72oLDQ3VoEGD2qVeAAAAAID3EfbAJwQFBemGG27Q9ddf72qrrq5WXl6e2w5dhmHoiy++0J49ey76nrfeeqtuvfXWdg979u7dqzfeeMNtBtKF+v+zqKgoTZw4UePGjXNrP/vYGgAAAADAnLpSgOMJwh6c459DDn9/f/Xr10833nijmpqaXO0FBQU6evSoawv38ykrK9O2bdtUVFTUbvU6nU4dOnRIjY2NF71RIyMj1a9fP7fHz/r06aPIyEjCHQAAAADwIczsgc+y2+2aNGmSxo0b5/YY1/Lly/XKK69c9LGp/fv361e/+pXsdnu71WgYhk6dOnXRWiSpf//+evbZZ922nw8KClJCQkK71QcAAAAA6HwIe+CzrFar4uLi3B7DampqUlJSkoKDgy862BsaGvTNN9+0+01hs9kUEBBw0X6RkZHq37+/rrzyynatBwAAAADQuRH2AP/EarXqmmuu0ZNPPqmGhoYL9j18+LDef/99lZWVtVs9FotFw4cP1w033KDAwMAL9k1NTVX37t3brRYAAAAAALyNsAces1qtGjlypNuuXefz4Ycfatu2be0e9qSnp2vmzJmKiIi4YF+r1So/P4Y9AAAAAPg6ZvYA32Gz2WSz2S7aLyoqSgMHDrzojJtLrSU5OVnBwcHtujYQAAAAAMA8CHuANrryyiv1i1/8QjU1Ne12DYvFosTERAUFBbXbNQAAAAAA5kLYA7RRVFSUoqKivF0GAAAAAABuCHsAAAAAAABMhLAHAAAAAADARAh7AAAAAAAATISwBwAAAAAAwEQIewAAAAAAAEyEsAcAAAAAAMBEzBz2WD3p3NTUpH/9139VamqqgoKC1Lt3b82dO1dOp9PVxzAMZWdnKzExUUFBQbr++uu1b9++y144AAAAAABAW50Nezx5dRUehT2//e1v9corr2jRokU6cOCAXnrpJf3ud7/Tyy+/7Orz0ksvacGCBVq0aJG2b9+u+Ph4jR8/XtXV1Ze9eAAAAAAAgLYg7PkfW7du1V133aWJEyeqV69e+v73v6+MjAzt2LFD0rd/UQsXLtScOXM0adIkpaWlafny5aqpqdGKFSva5QMAAAAAAADgf3kU9lx77bX66KOPdOjQIUnSF198oU8//VS33XabJCkvL08Oh0MZGRmuc+x2u8aNG6ctW7a0+J719fWqqqpyewEAAAAAALQnM8/s8WiB5meeeUaVlZUaMGCAbDabmpub9Zvf/EY/+MEPJEkOh0OSFBcX53ZeXFyc8vPzW3zP+fPn64UXXmhL7QAAAAAAAG3CAs3/Y+XKlXrzzTe1YsUK7dq1S8uXL9e//du/afny5W79LBaL258Nwzin7azZs2ersrLS9SosLPTwIwAAAAAAAHiuI2b1LF68WKmpqQoMDFR6ero2b9583r6rV6/W+PHj1b17d4WHh2vUqFH68MMPPb6mR2HPU089pWeffVb333+/vve97+mBBx7QE088ofnz50uS4uPjJf3vDJ+zSkpKzpntc5bdbld4eLjbCwAAAAAAoD11xGNcK1eu1MyZMzVnzhzl5uZq7NixmjBhggoKClrsv2nTJo0fP17r1q3Tzp07dcMNN+iOO+5Qbm6uR9f1KOypqamR1ep+is1mc229npqaqvj4eOXk5LiONzQ0aOPGjRo9erRHhQEAAAAAALSXjgh7FixYoGnTpumRRx7RwIEDtXDhQiUnJ2vJkiUt9l+4cKGefvppDR8+XP369dO8efPUr18//e1vf/Pouh6t2XPHHXfoN7/5jXr27KlBgwYpNzdXCxYs0I9//GNJ3z6+NXPmTFcxZwsLDg7WlClTPCoMAAAAAACgvbR1zZ7vbixlt9tlt9vP6d/Q0KCdO3fq2WefdWvPyMg47yZW3+V0OlVdXa2oqKhW1yl5GPa8/PLL+sUvfqGsrCyVlJQoMTFRjz76qH75y1+6+jz99NOqra1VVlaWysvLNXLkSK1fv15hYWEeFQYAAAAAANBe2hr2JCcnu7U///zzys7OPqd/aWmpmpubW9zE6rvL35zPv//7v+vMmTO67777Wl2n5GHYExYWpoULF2rhwoXn7WOxWJSdnd3iBwUAAAAAAOgM2hr2FBYWuq033NKsnn/mySZW/+ytt95Sdna23nvvPcXGxra6TsnDsAcAAAAAAMAM2hr2tHZzqZiYGNlsNo82sTpr5cqVmjZtmt555x3dfPPNra7xLI8WaAYAAAAAADCD9l6gOSAgQOnp6W6bWElSTk7OBTexeuutt/TQQw9pxYoVmjhxYps+GzN7AAAAAACAz2nrzB5PzJo1Sw888ICGDRumUaNGaenSpSooKFBmZqYkafbs2Tp27JjeeOMNSd8GPVOnTtXvf/97XXPNNa5ZQUFBQYqIiGj1dQl7AAAAAACAz+mIsGfy5MkqKyvT3LlzVVxcrLS0NK1bt04pKSmSpOLiYhUUFLj6/8d//Ieampo0ffp0TZ8+3dX+4IMP6vXXX2/1dQl7AAAAAAAA2klWVpaysrJaPPbdAGfDhg2X5ZqEPQAAAAAAwOd0xMwebyHsAQAAAAAAPoewBwAAAAAAwEQIewAAAAAAAEyEsAcAAAAAAMBECHsAAAAAAABMhLAHAAAAAADARAh7AAAAAAAATISwBwAAAAAAwGS6UoDjCcIeAAAAAADgc5jZAwAAAAAAYCKEPQAAAAAAACZC2AMAAAAAAGAiZg57rN4uAAAAAAAAAJcPM3sAAAAAAIDPMfPMHsIeAAAAAADgcwh7AAAAAAAATISwBwAAAAAAwEQIewAAAAAAAEyEsAcAAAAAAMBECHsAAAAAAABMhLAHAAAAAADARAh7AAAAAAAATISwBwAAAAAAwEQIewAAAAAAAEyEsAcAAAAAAMBEzBz2WL1dAAAAAAAAAC4fZvYAAAAAAACfY+aZPYQ9AAAAAADAJ3WlAMcThD0AAAAAAMDneBr0dKVgiLAHAAAAAAD4HMIeAAAAAAAAEyHsAQAAAAAAMBHCHgAAAAAAABMh7AEAAAAAADARwh4AAAAAAAATIewBAAAAAAAwEcKeDnT2L6+qquqyv3d1dbUaGxvd2pxOp2pra9vlegAAAAAAwDNn/33e3uGKmcMei9HJqi0qKlJycrK3ywAAAAAAAF5UWFiopKSky/6+VVVVioiIkN1ul8ViafV5hmGovr5elZWVCg8Pv+x1XU6dbmZPYmKiCgsLZRiGevbsqcLCwk7/l4jLr6qqSsnJyXz/PowxAMaAb+P7B2MAjAHfxvfv2wzDUHV1tRITE71dSpfV6cIeq9WqpKQk17St8PBwbm4fxvcPxgAYA76N7x+MATAGfBvfv++KiIho92uY+TGuThf2AAAAAAAAtDfCHgAAAAAAABMh7PECu92u559/Xna73dulwAv4/sEYAGPAt/H9gzEAxoBv4/tHRzBz2NPpduMCAAAAAABoL2d347JarR7vxuV0Oj3ejWvx4sX63e9+p+LiYg0aNEgLFy7U2LFjz9t/48aNmjVrlvbt26fExEQ9/fTTyszMbPX1JMnqUW8AAAAAAAATMAzD45enVq5cqZkzZ2rOnDnKzc3V2LFjNWHCBBUUFLTYPy8vT7fddpvGjh2r3NxcPffcc3r88ce1atUqj67LzB4AAAAAAOAzzs7skeTxzB5JHs3sGTlypIYOHaolS5a42gYOHKi7775b8+fPP6f/M888o7Vr1+rAgQOutszMTH3xxRfaunVrq2tlZg8AAAAAAPBJbZnVU1VV5faqr69v8b0bGhq0c+dOZWRkuLVnZGRoy5YtLZ6zdevWc/rfcsst2rFjhxobG1v9uQh7AAAAAACAzwgICFB8fHybzg0NDVVycrIiIiJcr5Zm6EhSaWmpmpubFRcX59YeFxcnh8PR4jkOh6PF/k1NTSotLW11nZ0y7Fm8eLFSU1MVGBio9PR0bd682dsloZ1kZ2fLYrG4vf75pjMMQ9nZ2UpMTFRQUJCuv/567du3z4sV41Js2rRJd9xxhxITE2WxWPTuu++6HW/N911fX6/HHntMMTExCgkJ0Z133qmioqIO/BS4FBcbAw899NA5vwnXXHONWx/GQNc1f/58DR8+XGFhYYqNjdXdd9+tgwcPuvXhd8DcWjMG+B0wryVLluiqq65SeHi4wsPDNWrUKL3//vuu49z/5nexMcD9j44SGBiovLw8VVZWevwqKio6p2327NkXvN53HxUzDOOCj4+11L+l9gvpdGGPp4sXoesbNGiQiouLXa8vv/zSdeyll17SggULtGjRIm3fvl3x8fEaP368qqurvVgx2urMmTMaPHiwFi1a1OLx1nzfM2fO1Jo1a/T222/r008/1enTp3X77berubm5oz4GLsHFxoAk3XrrrW6/CevWrXM7zhjoujZu3Kjp06dr27ZtysnJUVNTkzIyMnTmzBlXH34HzK01Y0Did8CskpKS9OKLL2rHjh3asWOHbrzxRt11112uQIf73/wuNgYk7n90nMDAQFfw6MkrIiLinDa73d7iNWJiYmSz2c6ZxVNSUnLO7J2z4uPjW+zv5+en6Ojo1n9Ao5MZMWKEkZmZ6dY2YMAA49lnn/VSRWhPzz//vDF48OAWjzmdTiM+Pt548cUXXW11dXVGRESE8corr3RQhWgvkow1a9a4/tya77uiosLw9/c33n77bVefY8eOGVar1fjggw86rHZcHt8dA4ZhGA8++KBx1113nfccxoC5lJSUGJKMjRs3GobB74Av+u4YMAx+B3xNZGSk8ec//5n734edHQOGwf0PcxoxYoTx05/+1K1t4MCB5804nn76aWPgwIFubZmZmcY111zj0XU71cyetixehK7v8OHDSkxMVGpqqu6//34dOXJE0rdbzjkcDrfxYLfbNW7cOMaDCbXm+965c6caGxvd+iQmJiotLY0xYSIbNmxQbGys+vfvr5/85CcqKSlxHWMMmEtlZaUkKSoqShK/A77ou2PgLH4HzK+5uVlvv/22zpw5o1GjRnH/+6DvjoGzuP9hNrNmzdKf//xnLVu2TAcOHNATTzyhgoICZWZmSpJmz56tqVOnuvpnZmYqPz9fs2bN0oEDB7Rs2TK9+uqrevLJJz26rt9l/RSXqC2LF6FrGzlypN544w31799fJ06c0K9//WuNHj1a+/btc33nLY2H/Px8b5SLdtSa79vhcCggIECRkZHn9OE3whwmTJige++9VykpKcrLy9MvfvEL3Xjjjdq5c6fsdjtjwEQMw9CsWbN07bXXKi0tTRK/A76mpTEg8Ttgdl9++aVGjRqluro6hYaGas2aNbryyitd/1Dn/je/840Bifsf5jR58mSVlZVp7ty5Ki4uVlpamtatW6eUlBRJUnFxsduyNampqVq3bp2eeOIJ/fGPf1RiYqL+8Ic/6J577vHoup0q7DnL08WL0HVNmDDB9d/f+973NGrUKPXp00fLly93LcbGePAtbfm+GRPmMXnyZNd/p6WladiwYUpJSdF///d/a9KkSec9jzHQ9cyYMUN79uzRp59+es4xfgd8w/nGAL8D5nbFFVdo9+7dqqio0KpVq/Tggw9q48aNruPc/+Z3vjFw5ZVXcv/DtLKyspSVldXisddff/2ctnHjxmnXrl2XdM1O9RhXWxYvgrmEhIToe9/7ng4fPuzalYvx4Bta833Hx8eroaFB5eXl5+0Dc0lISFBKSooOHz4siTFgFo899pjWrl2rTz75RElJSa52fgd8x/nGQEv4HTCXgIAA9e3bV8OGDdP8+fM1ePBg/f73v+f+9yHnGwMt4f4H2q5ThT0BAQFKT09XTk6OW3tOTo5Gjx7tparQkerr63XgwAElJCQoNTVV8fHxbuOhoaFBGzduZDyYUGu+7/T0dPn7+7v1KS4u1t69exkTJlVWVqbCwkIlJCRIYgx0dYZhaMaMGVq9erU+/vhjpaamuh3nd8D8LjYGWsLvgLkZhqH6+nrufx92dgy0hPsfuAQeLefcAd5++23D39/fePXVV439+/cbM2fONEJCQoyjR496uzS0g5///OfGhg0bjCNHjhjbtm0zbr/9diMsLMz1fb/44otGRESEsXr1auPLL780fvCDHxgJCQlGVVWVlytHW1RXVxu5ublGbm6uIclYsGCBkZuba+Tn5xuG0brvOzMz00hKSjL+/ve/G7t27TJuvPFGY/DgwUZTU5O3PhY8cKExUF1dbfz85z83tmzZYuTl5RmffPKJMWrUKKNHjx6MAZP46U9/akRERBgbNmwwiouLXa+amhpXH34HzO1iY4DfAXObPXu2sWnTJiMvL8/Ys2eP8dxzzxlWq9VYv369YRjc/77gQmOA+x+4vDpd2GMYhvHHP/7RSElJMQICAoyhQ4e6bccJc5k8ebKRkJBg+Pv7G4mJicakSZOMffv2uY47nU7j+eefN+Lj4w273W5cd911xpdffunFinEpPvnkE0PSOa8HH3zQMIzWfd+1tbXGjBkzjKioKCMoKMi4/fbbjYKCAi98GrTFhcZATU2NkZGRYXTv3t3w9/c3evbsaTz44IPnfL+Mga6rpe9ekvHaa6+5+vA7YG4XGwP8Dpjbj3/8Y9f/43fv3t246aabXEGPYXD/+4ILjQHuf+DyshiGYXTcPCIAAAAAAAC0p061Zg8AAAAAAAAuDWEPAAAAAACAiRD2AAAAAAAAmAhhDwAAAAAAgIkQ9gAAAAAAAJgIYQ8AAAAAAICJEPYAAAAAAACYCGEPAAAAAACAiRD2AAAAAAAAmAhhDwAAAAAAgIkQ9gAAAAAAAJgIYQ8AAAAAAICJ/H+hem21z7rFGgAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "nx = 400\n", "arguments = InitialConditions.genShockBubble(nx, nx//4, 1.4)\n", "plt.figure()\n", "plt.imshow(Visualization.genSchlieren(arguments['rho']), cmap='gray')\n", "plt.colorbar(orientation='vertical', aspect=20, pad=0.02, shrink=0.3)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Creating directory /home/smyalygames/PycharmProjects/FiniteVolumeGPU/data\n", "Initialized /home/smyalygames/PycharmProjects/FiniteVolumeGPU/data/euler_shock-bubble.nc\n", "Opening /home/smyalygames/PycharmProjects/FiniteVolumeGPU/data/euler_shock-bubble.nc\n", "Keyword arguments: {'mode': 'w', 'clobber': False}\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Constructed in 0.07140874862670898 seconds\n", "Simulating to t=0.500000. Estimated 488 timesteps (dt=0.001023)\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "Closing /home/smyalygames/PycharmProjects/FiniteVolumeGPU/data/euler_shock-bubble.nc\n" ] } ], "source": [ "nx = 800#1600\n", "ny = nx//4\n", "g = 0.0\n", "gamma = 1.4\n", "t_end = 0.5#3.0\n", "n_save = 20#500\n", "outfile = \"data/euler_shock-bubble.nc\"\n", "outdata = None\n", "\n", "arguments = InitialConditions.genShockBubble(nx, ny, gamma)\n", "arguments['context'] = my_context\n", "outfile = runSimulation(outfile, t_end, arguments)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Initialized /home/smyalygames/PycharmProjects/FiniteVolumeGPU/data/euler_shock-bubble.nc\n", "Opening /home/smyalygames/PycharmProjects/FiniteVolumeGPU/data/euler_shock-bubble.nc\n", "Arguments: ('r',)\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "/home/smyalygames/PycharmProjects/FiniteVolumeGPU/data/euler_shock-bubble.nc created Wed Feb 12 16:40:00 2025 contains 21 timesteps\n", "Simulator arguments: \n", " {'rho': '[...]', 'rho_u': '[...]', 'rho_v': '[...]', 'E': '[...]', 'nx': 800, 'ny': 200, 'dx': 0.005, 'dy': 0.005, 'g': 0.0, 'gamma': 1.4, 'boundary_conditions': '[north=Type.Reflective, south=Type.Reflective, east=Type.Periodic, west=Type.Periodic]', 'context': 'CudaContext id 96726128804080'}\n" ] }, { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "
\n", " \n", "
\n", " \n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "
\n", "
\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Closing /home/smyalygames/PycharmProjects/FiniteVolumeGPU/data/euler_shock-bubble.nc\n" ] } ], "source": [ "#outfile = 'data/euler_shock-bubble_0008.nc'\n", "createAnimation(outfile, vis_type=VisType.Schlieren, vmin=-0.2, vmax=0.2, cmap=plt.cm.RdBu, save_anim=False, fig_args={'figsize':(16, 5)})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Kelvin-Helmholtz" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAABHsAAAE8CAYAAABQNL2dAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAvZUlEQVR4nO3de3SdVZk4/uekl9NbEinSpKGlBCwCVkptEVtuBW2Zgh1ZOIjiWC7ODP1SwFq5VUYIjjZStatiFQYHuSy5LYUCrkFtZ7SpCPymLY1g4Qs4VFqwscIXm/SWkJz390ftmQm95Q1J05zz+bj2WmSf/b7vk559zmofn713JkmSJAAAAAAoCCU9HQAAAAAAXUeyBwAAAKCASPYAAAAAFBDJHgAAAIACItkDAAAAUEAkewAAAAAKiGQPAAAAQAGR7AEAAAAoIJI9AAAAAAVEsgcAAACggEj2AAAAAHSx2traOOGEE6K0tDSGDRsW55xzTrz44ov7vK6uri7Gjx8fAwYMiCOOOCJuu+221M+W7AEAAADoYnV1dTFr1qx4+umnY+nSpdHa2hpTp06NLVu27PGatWvXxllnnRWnnHJKrF69Or785S/HlVdeGQ899FCqZ2eSJEne7S8AAAAAwJ79+c9/jmHDhkVdXV2ceuqpux1z7bXXxmOPPRYvvPBCvm/mzJnx29/+Np566qkOP6vvu44WAAAAoBfZvn17tLS0pL4uSZLIZDLt+rLZbGSz2X1eu2nTpoiIGDp06B7HPPXUUzF16tR2fWeeeWbccccd8fbbb0e/fv06FKdkDwAAAFA0tm/fHtWjhkTDxrbU1w4ZMiQ2b97cru/GG2+MmpqavV6XJEnMmTMnTj755BgzZswexzU0NERFRUW7voqKimhtbY033ngjhg8f3qE4JXsAAACAotHS0hING9ti7apRUVba8a2MG5tyUT3+1Vi/fn2UlZXl+ztS1XP55ZfHs88+G0888cQ+x76zcmjn7jvv7N8byR4AAACg6AwesqN1VNtfdzwuKytrl+zZlyuuuCIee+yxWL58eYwYMWKvYysrK6OhoaFd38aNG6Nv375x8MEHd/iZkj0AAABA0clFErno+JlVacZG7KjIueKKK2Lx4sWxbNmyqK6u3uc1EydOjJ/+9Kft+pYsWRITJkzo8H49EY5eBwAAAIpQrhP/S2PWrFnxox/9KO67774oLS2NhoaGaGhoiG3btuXHzJ07N2bMmJH/eebMmfHqq6/GnDlz4oUXXogf/vCHcccdd8RVV12V6tmSPQAAAEDRaUuS1C2NW2+9NTZt2hSTJ0+O4cOH59uDDz6YH7Nhw4ZYt25d/ufq6up4/PHHY9myZXH88cfHv/zLv8Qtt9wSn/zkJ1M9O5MkKaMFAAAA6KUaGxujvLw81v/fQ1Nv0Dzy6Ndj06ZNqfbs6QndVtnz/e9/P6qrq2PAgAExfvz4+PWvf91djwIAAABIZeeePWlab9EtyZ4HH3wwZs+eHddff32sXr06TjnllJg2bVq70iQAAACAnpKLJNpStN6U7OmWZVwnnnhifOhDH4pbb70133fMMcfEOeecE7W1tXu9NpfLxR//+McoLS1NdYY8AAAA0PslSRJNTU1RVVUVJSVdX6OycxnXf//fyihNsYyrqSkXRx7d0CuWcXX50estLS2xatWquO6669r1T506NZ588sldxjc3N0dzc3P+59dffz2OPfbYrg4LAAAA6EXWr18fI0aM6Lb7p910Oe0GzT2py5M9b7zxRrS1tUVFRUW7/oqKimhoaNhlfG1tbdx000279L/6zOFRNsRhYQAAAFBMGjfnYtSH/hClpaXd+pzcX1ua8b1Flyd7dnrnEqwkSXa7LGvu3LkxZ86c/M+NjY0xcuTIKBtSkmpXbAAAAKBwdPfWLjv34kkzvrfo8mTPe9/73ujTp88uVTwbN27cpdonIiKbzUY2m+3qMAAAAAD2qC3Z0dKM7y26vHSmf//+MX78+Fi6dGm7/qVLl8akSZO6+nEAAAAAqeU60XqLblnGNWfOnPjc5z4XEyZMiIkTJ8btt98e69ati5kzZ3bH4wAAAABSyUUm2qLjS8VyKcb2tG5J9px//vnx5ptvxle/+tXYsGFDjBkzJh5//PEYNWpUdzwOAAAAIJVcsqOlGd9bdNsGzZdddllcdtll3XV7AAAAgE5rS1nZk2ZsT+u2ZA8AAADAgaqQkz3ONgcAAAAoICp7AAAAgKKTSzKRS1Js0JxibE+T7AEAAACKTiEv45LsAQAAAIpOW5REW4rdbdq6MZauJtkDAAAAFJ0k5TKuxDIuAAAAgAOXZVwAAAAABaQtKYm2JMUyrqQbg+likj0AAABA0clFJnIp9uzJRe/J9kj2AAAAAEXHMi4AAACAApJ+GZfKHgAAAIAD1o5lXB2v1kkztqdJ9gAAAABFJxcl0WbPHgAAAIDCUMjLuDr+WwEAAABwwFPZAwAAABSdXJQ4eh0AAACgULQlmWhLUhy9nmJsT5PsAQAAAIpOW8oNmttU9gAAAAAcuHJJSeRSbNCc60UbNEv2AAAAAEVHZQ8AAABAAclFun14ct0XSpeT7AEAAACKTvrTuDo+tqdJ9gAAAABFpy0pibYUe/akGdvTJHsAAACAopOLTOQizTIuR68DAAAAHLBU9gAAAAAUkPSncUn2AAAAABywckkmcmlO40oxtqf1nrQUAAAAAPuksgcAAAAoOrmUy7gcvQ4AAABwAMslJZFLselymrE9TbIHAAAAKDptkYm2FMeppxnb0yR7AAAAgKKjsgcAAACggLRFumqdtu4Lpcv1nrQUAAAAQBfZWdmTpqW1fPnymD59elRVVUUmk4lHHnlkn9fce++9MXbs2Bg0aFAMHz48Lr744njzzTdTPVeyBwAAACg6bUlJ6pbWli1bYuzYsbFo0aIOjX/iiSdixowZ8fnPfz7WrFkTP/7xj2PFihXxD//wD6meaxkXAAAAUHSSyEQuxTKupBMbNE+bNi2mTZvW4fFPP/10HH744XHllVdGRER1dXVceumlMX/+/FTPVdkDAAAAFJ3OVvY0Nja2a83NzV0W06RJk+K1116Lxx9/PJIkiT/96U/xk5/8JM4+++xU95HsAQAAAIpOLsmkbhERI0eOjPLy8nyrra3tspgmTZoU9957b5x//vnRv3//qKysjPe85z3x3e9+N9V9LOMCAAAAik5blERbihqYnWPXr18fZWVl+f5sNttlMT3//PNx5ZVXxg033BBnnnlmbNiwIa6++uqYOXNm3HHHHR2+j2QPAAAAUHT+d7VOR8dHRJSVlbVL9nSl2traOOmkk+Lqq6+OiIjjjjsuBg8eHKecckp87Wtfi+HDh3foPpI9AAAAQNHJRUnkUlT2pBnbWVu3bo2+fdunavr06RMREUmSdPg+qSKtra2NE044IUpLS2PYsGFxzjnnxIsvvthuTJIkUVNTE1VVVTFw4MCYPHlyrFmzJs1jAAAAAHq9zZs3R319fdTX10dExNq1a6O+vj7WrVsXERFz586NGTNm5MdPnz49Hn744bj11lvjlVdeid/85jdx5ZVXxoc//OGoqqrq8HNTJXvq6upi1qxZ8fTTT8fSpUujtbU1pk6dGlu2bMmPmT9/fixYsCAWLVoUK1asiMrKypgyZUo0NTWleRQAAABAt2lLMqlbWitXroxx48bFuHHjIiJizpw5MW7cuLjhhhsiImLDhg35xE9ExEUXXZTPqYwZMybOO++8eP/73x8PP/xwqudmkjR1QO/w5z//OYYNGxZ1dXVx6qmnRpIkUVVVFbNnz45rr702IiKam5ujoqIibr755rj00kv3ec/GxsYoLy+Pt146IspKHRYGAAAAxaSxKRcHHfVKbNq0qVv2xtmZd7h0+ScjO6Rfh69r3vx2/OupD3VbXF3pXWVTNm3aFBERQ4cOjYgd5UgNDQ0xderU/JhsNhunnXZaPPnkk7u9R3Nz8y5n1AMAAAB0pyQpiVyKliS9pyCl05EmSRJz5syJk08+OcaMGRMREQ0NDRERUVFR0W5sRUVF/rV3qq2tbXc+/ciRIzsbEgAAAECHtEUmdestOp3sufzyy+PZZ5+N+++/f5fXMpn2fwBJkuzSt9PcuXNj06ZN+bZ+/frOhgQAAADQIbnkf45f71jr6Yg7rlNHr19xxRXx2GOPxfLly2PEiBH5/srKyojYUeHzv89+37hx4y7VPjtls9nIZrOdCQMAAACgU3Yuz0ozvrdIFWmSJHH55ZfHww8/HL/85S+jurq63evV1dVRWVkZS5cuzfe1tLREXV1dTJo0qWsiBgAAAHiXcpFJ3XqLVJU9s2bNivvuuy8effTRKC0tze/DU15eHgMHDoxMJhOzZ8+OefPmxejRo2P06NExb968GDRoUFxwwQXd8gsAAAAApJX2OPXOHL3eU1Ile2699daIiJg8eXK7/jvvvDMuuuiiiIi45pprYtu2bXHZZZfFW2+9FSeeeGIsWbIkSktLuyRgAAAAgHerkJdxpUr2JMm+dyPKZDJRU1MTNTU1nY0JAAAAoFvlYsfGy2nG9xad2qAZAKAQrWnZFj9886TYsL18l9dKMkl8bOjz8ZnS1yOb6dcD0QEAXSlJuQ9PItkDAND7/Osbp8b/d8uEKH95664vlmTips8eFWd8fEEc1leyBwB6u51HqqcZ31scsMmeN9u2REtb71kPBwD0fi81DouDnm+KZOXvdn2xpE8MnHxirG8dFIMyW/Z/cABQJJracj0dQq93wCZ7Prr4iigZMKCnwwAAisiQP5TEoQ2vRuvuXkxyUbmiJS4ePCty2X3vYwgAdE5u+/aI+Ofuf44Nmve/I7/2u+ib6d/TYQAARSRpa4vWlpY9vJhE/1/Wx5FP+PsJAHSn1qQlXt0Pz7GMqwfktm6LXGa3/78aAECPSFpbI2n19xMA6E655O3985yUGzQ7jQsAAADgAKayBwAAAKCASPYAAAAAFBDJHgAAAIACItkDAAAAUECSSLfpctJ9oXQ5yR4AAACg6KjsAQAAACggkj0AAAAABUSyBwAAAKCAFHKyp6SnAwAAAACg66jsAQAAAIpOkmQiSVGtk2ZsT5PsAQAAAIpOLjKpjl5PM7anSfYAAAAARaeQ9+yR7AEAAACKjmVcAAAAAAVEZQ8AAABAAVHZAwAAAFBAkpSVPZI9AAAAAAewJCKSJN343kKyBwAAACg6uchExtHrAAAAAIXBnj0AAAAABSSXZCLjNC4AAACAwpAkKffs6UWb9pT0dAAAAAAAdB2VPQAAAEDRsWcPAAAAQAGR7AEAAAAoIDZoBgAAACgghbxBs2QPAAAAUHR2JHvSLOPqxmC6mGQPAAAAUHTs2QMAAABQQJK/tjTjewvJHgAAAKDoFHJlT0lPBwAAAACw3yWdaCktX748pk+fHlVVVZHJZOKRRx7Z5zXNzc1x/fXXx6hRoyKbzcaRRx4ZP/zhD1M9V2UPAAAAUHxSVvZEJyp7tmzZEmPHjo2LL744PvnJT3bomk996lPxpz/9Ke6444543/veFxs3bozW1tZUz5XsAQAAAIrO/jh6fdq0aTFt2rQOj//5z38edXV18corr8TQoUMjIuLwww9P/dx3tYyrtrY2MplMzJ49O9+XJEnU1NREVVVVDBw4MCZPnhxr1qx5N48BAAAA6FI79+xJ0yIiGhsb27Xm5uYui+mxxx6LCRMmxPz58+PQQw+No446Kq666qrYtm1bqvt0OtmzYsWKuP322+O4445r1z9//vxYsGBBLFq0KFasWBGVlZUxZcqUaGpq6uyjAAAAAA4II0eOjPLy8nyrra3tsnu/8sor8cQTT8Tvfve7WLx4cSxcuDB+8pOfxKxZs1Ldp1PJns2bN8dnP/vZ+MEPfhAHHXRQvj9Jkli4cGFcf/31ce6558aYMWPi7rvvjq1bt8Z9993XmUcBAAAAdL0kk75FxPr162PTpk35Nnfu3C4LKZfLRSaTiXvvvTc+/OEPx1lnnRULFiyIu+66K1V1T6eSPbNmzYqzzz47Pvaxj7XrX7t2bTQ0NMTUqVPzfdlsNk477bR48sknd3uv5ubmXUqgAAAAALrTzj170rSIiLKysnYtm812WUzDhw+PQw89NMrLy/N9xxxzTCRJEq+99lqH75M62fPAAw/EM888s9sypYaGhoiIqKioaNdfUVGRf+2damtr25U/jRw5Mm1IAAAAAOnsh6PX0zrppJPij3/8Y2zevDnf99JLL0VJSUmMGDGiw/dJlexZv359fOELX4gf/ehHMWDAgD2Oy2TaH0eWJMkufTvNnTu3XfnT+vXr04QEAAAAkFpnN2hOY/PmzVFfXx/19fURsWNFVH19faxbty4iduREZsyYkR9/wQUXxMEHHxwXX3xxPP/887F8+fK4+uqr45JLLomBAwd2+Lmpkj2rVq2KjRs3xvjx46Nv377Rt2/fqKuri1tuuSX69u2br+h5ZxXPxo0bd6n22Smbze5SAgUAAADQ7bq5qmflypUxbty4GDduXEREzJkzJ8aNGxc33HBDRERs2LAhn/iJiBgyZEgsXbo0/vKXv8SECRPis5/9bEyfPj1uueWWVM/tm2bwRz/60Xjuuefa9V188cVx9NFHx7XXXhtHHHFEVFZWxtKlS/O/SEtLS9TV1cXNN9+cKjAAAACA7pK2WqczlT2TJ0+OJNlzpuiuu+7ape/oo4+OpUuXpn7W/5Yq2VNaWhpjxoxp1zd48OA4+OCD8/2zZ8+OefPmxejRo2P06NExb968GDRoUFxwwQXvKlAAAACALpO2Ymc/7NnTVVIlezrimmuuiW3btsVll10Wb731Vpx44omxZMmSKC0t7epHAQAAAHRS5q8tzfje4V0ne5YtW9bu50wmEzU1NVFTU/Nubw0AAADQPVT2AAAAABQQyR4AAACAApJkdrQ043sJyR4AAACg6CTJjpZmfG9R0tMBAAAAANB1VPYAAAAAxceePQAAAAAFxJ49AAAAAIUjk+xoacb3FpI9AAAAQPGxjAsAAACggFjGBQAAAFBAVPYAAAAAFBDJHgAAAIACItkDAAAAUEDs2QMAAABQOBy9DgAAAFBICngZV0lPBwAAAABA15HsAQAAACgglnEBAAAARScTKffs6bZIup5kDwAAAFB8nMYFAAAAUEAKeINmyR4AAACg+Ej2AAAAABSOTJJyzx7JHgAAAIADmMoeAAAAgAIi2QMAAABQOCzjAgAAACgkjl4HAAAAKCCWcQEAAAAUDsu4AAAAAApJAVf2lPR0AAAAAAB0HZU9AAAAQPFJuYyrN1X2SPYAAAAAxaeAl3FJ9gAAAADFR7IHAAAAoHAU8mlcNmgGAAAAKCAqewAAAIDiYxkXAAAAQOEo5GVckj0AAABAcepFCZw0JHsAAACA4mMZFwAAAEDhsIwLAAAAoJAUcGVP6qPXX3/99fj7v//7OPjgg2PQoEFx/PHHx6pVq/KvJ0kSNTU1UVVVFQMHDozJkyfHmjVrujRoAAAAgHdjZ2VPmpbW8uXLY/r06VFVVRWZTCYeeeSRDl/7m9/8Jvr27RvHH3986uemSva89dZbcdJJJ0W/fv3iZz/7WTz//PPx7W9/O97znvfkx8yfPz8WLFgQixYtihUrVkRlZWVMmTIlmpqaUgcHAAAA0Ftt2bIlxo4dG4sWLUp13aZNm2LGjBnx0Y9+tFPPTbWM6+abb46RI0fGnXfeme87/PDD8/+dJEksXLgwrr/++jj33HMjIuLuu++OioqKuO++++LSSy/tVJAAAAAAXWo/LOOaNm1aTJs2LfV1l156aVxwwQXRp0+fVNVAO6Wq7HnsscdiwoQJcd5558WwYcNi3Lhx8YMf/CD/+tq1a6OhoSGmTp2a78tms3HaaafFk08+udt7Njc3R2NjY7sGAAAA0K2STrSIXXIYzc3NXRrWnXfeGf/93/8dN954Y6fvkSrZ88orr8Stt94ao0ePjl/84hcxc+bMuPLKK+Oee+6JiIiGhoaIiKioqGh3XUVFRf61d6qtrY3y8vJ8GzlyZGd+DwAAAIAO6+yePSNHjmyXx6itre2ymF5++eW47rrr4t57742+fTt/plaqK3O5XEyYMCHmzZsXERHjxo2LNWvWxK233hozZszIj8tkMu2uS5Jkl76d5s6dG3PmzMn/3NjYKOEDAAAAdK9OLuNav359lJWV5buz2WyXhNPW1hYXXHBB3HTTTXHUUUe9q3ulSvYMHz48jj322HZ9xxxzTDz00EMREVFZWRkROyp8hg8fnh+zcePGXap9dspms132BwMAAADQIZ1M9pSVlbVL9nSVpqamWLlyZaxevTouv/zyiNhRdJMkSfTt2zeWLFkSZ5xxRofulWoZ10knnRQvvvhiu76XXnopRo0aFRER1dXVUVlZGUuXLs2/3tLSEnV1dTFp0qQ0jwIAAADoNvvj6PU0ysrK4rnnnov6+vp8mzlzZrz//e+P+vr6OPHEEzt8r1SVPV/84hdj0qRJMW/evPjUpz4V//Vf/xW333573H777RGxY/nW7NmzY968eTF69OgYPXp0zJs3LwYNGhQXXHBBut8SAAAAoLvsh9O4Nm/eHL///e/zP69duzbq6+tj6NChcdhhh8XcuXPj9ddfj3vuuSdKSkpizJgx7a4fNmxYDBgwYJf+fUmV7DnhhBNi8eLFMXfu3PjqV78a1dXVsXDhwvjsZz+bH3PNNdfEtm3b4rLLLou33norTjzxxFiyZEmUlpamCgwAAACgu6St1ulMZc/KlSvj9NNPz/+8c8/iCy+8MO66667YsGFDrFu3Lv2N9yGTJEk3FyKl09jYGOXl5TE5PhF9M/16OhwAAABgP2pN3o5l8Whs2rSpW/bG2Zl3OGbWvOiTHdDh69qat8cL3/tyt8XVlTp/jhcAAABAb7UflnH1FMkeAAAAoOhk/trSjO8tJHsAAACA4lPAlT2pjl4HAAAA4MCmsgcAAAAoOvvjNK6eItkDAAAAFJ8CXsYl2QMAAAAUp16UwElDsgcAAAAoOpZxAQAAABQSy7gAAAAACofKHgAAAIBCorIHAAAAoHCo7AEAAAAoJCp7AAAAAAqIZA8AAABA4bCMCwAAAKCQFHBlT0lPBwAAAABA11HZAwAAABSdTJJEJul4uU6asT1NsgcAAAAoPgW8jEuyBwAAACg6NmgGAAAAKCQqewAAAAAKh8oeAAAAgEKisgcAAACgcKjsAQAAACgkKnsAAAAACktvqtZJQ7IHAACAd61k0KBoO350bB0+IAZt2B596l+O3NatPR0W7FmS7GhpxvcSJT0dAAAAAL1fScUh8fLn+8XfffUX8fLn+0VJxSE9HRLs1c49e9K03uKArewpGTwoSjL99zwgSSJpaYmktXWf98r07RuZbLYLo4P/kbS8HcnbLfseWNInSgZkIzKZ7g8KAAD2s7b3DIlRI9+I2Qf9IR4Z+Ua0vWdIlAwe3NNhsR8kzc0d/7d5//77/DdRSdISsaWrotsLe/bsf7+/4QNRMmDAHl/vuyUThy3ZHiV1q/d+o0wmtk8ZF6+f3jdy/XvRO0OvkGmLGP6bJAb/dPU+Ez6Z8cfGH84ujZb35PZTdAAAsP/kylrjqhHLIyLivBGr4ltXnBkljR/s4ajobiXbMzFiWWv0/8XKfS5zaj35uFh3ZjbaBu59XG779ohruzLK4nPAJnvq/nZRlJbueZXZr7e/N25ouCgq6vZ+n0yfPvGnE/rFA+ctjCP67jvTCGn8v1wuzkyuitG/6LfPZM//O3ZIfPH8R+K8Ib/fT9EBAMD+U5LJxKBM/4joE/9U/of4+6nfjVwv2uOEznnx7Wxc0nRFjFxSEpG07XXsn48fEIvO+0GckN2013FNTbmo3g/JnkxuR0szvrc4YJM97+kzKMr67DnZc1jft2LLyCSSk47f631yfTOx/dC3o6pPSxzUZ0gXR0mx65fZHlG1PVo+cnT02b73L7amwzNxeL8/x0F9Bu2n6AAAoGf0y/SJ8szAng6D/aAq2RzbRrRG7qTjIpPbc3IvyURsGZmLQ/s27vPfRH367KesimVcB5739Uvi6umPxlOnHrnXcSWRxEUHvRBD+9izh643MNM/vnXCT+Lfjxwbrcne9zs/q3R9fCj7l4iwbhkAACgMFX2ycdPpD0fduKMjF3vfi+e88v+OUX0PnD1M0266bIPm/WBIyYD4p/I/xj+V/7GDV/Tr1ngoTn0yJXHO4M1xzuDfdPAKiR4AAKBwZDP9YkbZGzGj7IkOXrHnvXn3uwI+er3XJnsAAAAAOktlDwAAAEAhsWcPAAAAQOFQ2QMAAABQSOzZAwAAAFA4VPYAAAAAFJIC3rOnJM3g1tbW+Od//ueorq6OgQMHxhFHHBFf/epXI5fL5cckSRI1NTVRVVUVAwcOjMmTJ8eaNWu6PHAAAACAztpZ2ZOm9Rapkj0333xz3HbbbbFo0aJ44YUXYv78+fHNb34zvvvd7+bHzJ8/PxYsWBCLFi2KFStWRGVlZUyZMiWampq6PHgAAACATskl6VsvkSrZ89RTT8UnPvGJOPvss+Pwww+Pv/u7v4upU6fGypUrI2JHVc/ChQvj+uuvj3PPPTfGjBkTd999d2zdujXuu+++bvkFAAAAAPgfqZI9J598cvznf/5nvPTSSxER8dvf/jaeeOKJOOussyIiYu3atdHQ0BBTp07NX5PNZuO0006LJ598crf3bG5ujsbGxnYNAAAAoFslnWi9RKoNmq+99trYtGlTHH300dGnT59oa2uLr3/96/GZz3wmIiIaGhoiIqKioqLddRUVFfHqq6/u9p61tbVx0003dSZ2AAAAgE7JRMrTuLotkq6XqrLnwQcfjB/96Edx3333xTPPPBN33313fOtb34q777673bhMpv0fQZIku/TtNHfu3Ni0aVO+rV+/PuWvAAAAAJBSkqRvKS1fvjymT58eVVVVkclk4pFHHtnr+IcffjimTJkShxxySJSVlcXEiRPjF7/4Rernpkr2XH311XHdddfFpz/96fjgBz8Yn/vc5+KLX/xi1NbWRkREZWVlRPxPhc9OGzdu3KXaZ6dsNhtlZWXtGgAAAEB32h+ncW3ZsiXGjh0bixYt6tD45cuXx5QpU+Lxxx+PVatWxemnnx7Tp0+P1atXp3puqmVcW7dujZKS9vmhPn365I9er66ujsrKyli6dGmMGzcuIiJaWlqirq4ubr755lSBAQAAAHSbtPvwdCLZM23atJg2bVqHxy9cuLDdz/PmzYtHH300fvrTn+bzLB2RKtkzffr0+PrXvx6HHXZYfOADH4jVq1fHggUL4pJLLomIHcu3Zs+eHfPmzYvRo0fH6NGjY968eTFo0KC44IIL0jwKAAAAoNtkkiQyKZZm7Rz7zoOlstlsZLPZLo1tp1wuF01NTTF06NBU16VK9nz3u9+Nr3zlK3HZZZfFxo0bo6qqKi699NK44YYb8mOuueaa2LZtW1x22WXx1ltvxYknnhhLliyJ0tLSVIEBAAAAdJvcX1ua8RExcuTIdt033nhj1NTUdFVU7Xz729+OLVu2xKc+9alU16VK9pSWlsbChQt3KSv63zKZTNTU1HTbLwoAAADwbnW2smf9+vXt9hvurqqe+++/P2pqauLRRx+NYcOGpbo2VbIHAAAAoCB0cs+e/XG41IMPPhif//zn48c//nF87GMfS329ZA8AAABQfNIep96Jo9c74/77749LLrkk7r///jj77LM7dQ/JHgAAAKDopD1OvTNHr2/evDl+//vf539eu3Zt1NfXx9ChQ+Owww6LuXPnxuuvvx733HNPROxI9MyYMSO+853vxEc+8pFoaGiIiIiBAwdGeXl5h59bsu8hAAAAAAVmZ2VPmpbSypUrY9y4cflj0+fMmRPjxo3LH3S1YcOGWLduXX78v/7rv0Zra2vMmjUrhg8fnm9f+MIXUj1XZQ8AAABAN5g8eXIke0kS3XXXXe1+XrZsWZc8V7IHAAAAKDqZ3I6WZnxvIdkDAAAAFJ8DdIPmriDZAwAAABSfTh693htI9gAAAABFJ5MkkUlRrZNmbE+T7AEAAACKj2VcAAAAAAUkiYg0my73nlyPZA8AAABQfCzjAgAAACgkSaRcxtVtkXQ5yR4AAACg+NizBwAAAKCA5CIik3J8LyHZAwAAABQde/YAAAAAFBLLuAAAAAAKSAEne0p6OgAAAAAAuo7KHgAAAKD4FHBlj2QPAAAAUHycxgUAAABQOJzGBQAAAFBILOMCAAAAKCC5JCKTIoGTk+wBAAAAOHCp7AEAAAAoJCmTPSHZAwAAAHDgUtkDAAAAUEBySaSq1rFnDwAAAMABLMntaGnG9xKSPQAAAEDxsYwLAAAAoIAU8DKukp4OAAAAAICuo7IHAAAAKD6WcQEAAAAUkCRSJnu6LZIuJ9kDAAAAFB+VPQAAAAAFJJeLiBTHqeccvQ4AAABw4FLZAwAAAFBAJHsAAAAACkguiVS7LuckewAAAAAOWEmSiyTp+D48acb2NMkeAAAAoPgkSbpqHcu4AAAAAA5gScplXJI9nZf89Q+vcXPvKY8CAAAAusbOfEDS3cmVXC4ikyL3YBlX5zU1NUVExKgP/aFnAwEAAAB6TFNTU5SXl3ffA1T27D9VVVWxfv36SJIkDjvssFi/fn2UlZX1dFjsZ42NjTFy5EjvfxEzBzAHipv3H3MAc6C4ef+LW5Ik0dTUFFVVVT0dSq91wCV7SkpKYsSIEdHY2BgREWVlZT7cRcz7jzmAOVDcvP+YA5gDxc37X7y6taLnr5JcLpIUy7icxgUAAABwILOMCwAAAKCA5JKIjGTPfpXNZuPGG2+MbDbb06HQA7z/mAOYA8XN+485gDlQ3Lz/7BdJEhFpTuPqPcmeTNLtZ5kBAAAAHBgaGxujvLw8Tu/7d9E306/D17Umb8evWn8SmzZt6vBeUsuXL49vfvObsWrVqtiwYUMsXrw4zjnnnL1eU1dXF3PmzIk1a9ZEVVVVXHPNNTFz5swOxxkRUZJqNAAAAEAhSHLpW0pbtmyJsWPHxqJFizo0fu3atXHWWWfFKaecEqtXr44vf/nLceWVV8ZDDz2U6rkH7DIuAAAAgO6S5JJIUuzZ05mFUdOmTYtp06Z1ePxtt90Whx12WCxcuDAiIo455phYuXJlfOtb34pPfvKTHb6PZA8AAABQdFqT5lTVOq3xdkTsWAb2v2Wz2S7bX+qpp56KqVOntus788wz44477oi33347+vXr2LIzyR4AAACgaPTv3z8qKyvjiYbHU187ZMiQGDlyZLu+G2+8MWpqaroktoaGhqioqGjXV1FREa2trfHGG2/E8OHDO3SfA3LPnu9///tRXV0dAwYMiPHjx8evf/3rng6JblJTUxOZTKZdq6yszL+eJEnU1NREVVVVDBw4MCZPnhxr1qzpwYh5N5YvXx7Tp0+PqqqqyGQy8cgjj7R7vSPvd3Nzc1xxxRXx3ve+NwYPHhx/+7d/G6+99tp+/C14N/Y1By666KJdvhM+8pGPtBtjDvRetbW1ccIJJ0RpaWkMGzYszjnnnHjxxRfbjfE9UNg6Mgd8DxSuW2+9NY477rgoKyuLsrKymDhxYvzsZz/Lv+7zX/j2NQd8/tlfBgwYEGvXro1Nmzalbq+99toufXPnzu3S+DKZTLufdy4fe2f/3hxwyZ4HH3wwZs+eHddff32sXr06TjnllJg2bVqsW7eup0Ojm3zgAx+IDRs25Ntzzz2Xf23+/PmxYMGCWLRoUaxYsSIqKytjypQp0dTU1IMR01n72pysI+/37NmzY/HixfHAAw/EE088EZs3b46Pf/zj0dbWtr9+Dd6FjmxQ9zd/8zftvhMef7z9/+NiDvRedXV1MWvWrHj66adj6dKl0draGlOnTo0tW7bkx/geKGwdmQMRvgcK1YgRI+Ib3/hGrFy5MlauXBlnnHFGfOITn8gndHz+C9++5kCEzz/7z4ABA/KJxzStvLx8l76uWsIVEVFZWRkNDQ3t+jZu3Bh9+/aNgw8+uOM3Sg4wH/7wh5OZM2e26zv66KOT6667rociojvdeOONydixY3f7Wi6XSyorK5NvfOMb+b7t27cn5eXlyW233bafIqS7RESyePHi/M8deb//8pe/JP369UseeOCB/JjXX389KSkpSX7+85/vt9jpGu+cA0mSJBdeeGHyiU98Yo/XmAOFZePGjUlEJHV1dUmS+B4oRu+cA0nie6DYHHTQQcm//du/+fwXsZ1zIEl8/il8u/v77ztdc801yTHHHNOub+bMmclHPvKRVM86oCp7WlpaYtWqVbtsRjR16tR48skneygqutvLL78cVVVVUV1dHZ/+9KfjlVdeiYgdR841NDS0mw/ZbDZOO+0086EAdeT9XrVqVbz99tvtxlRVVcWYMWPMiQKybNmyGDZsWBx11FHxj//4j7Fx48b8a+ZAYdm0aVNERAwdOjQifA8Uo3fOgZ18DxS+tra2eOCBB2LLli0xceJEn/8i9M45sJPPP4Vm8+bNUV9fH/X19RGx4+879fX1+dVLc+fOjRkzZuTHz5w5M1599dWYM2dOvPDCC/HDH/4w7rjjjrjqqqtSPfeA2qD5jTfeiLa2tt1uRvTOMiYKw4knnhj33HNPHHXUUfGnP/0pvva1r8WkSZNizZo1+fd8d/Ph1Vdf7Ylw6UYdeb8bGhqif//+cdBBB+0yxndEYZg2bVqcd955MWrUqFi7dm185StfiTPOOCNWrVoV2WzWHCggSZLEnDlz4uSTT44xY8ZEhO+BYrO7ORDhe6DQPffcczFx4sTYvn17DBkyJBYvXhzHHnts/h/qPv+Fb09zIMLnn8K0cuXKOP300/M/z5kzJyIiLrzwwrjrrrtiw4YN7batqa6ujscffzy++MUvxve+972oqqqKW265JdWx6xEHWLJnp91tRpRmIyJ6j2nTpuX/+4Mf/GBMnDgxjjzyyLj77rvzm7GZD8WlM++3OVE4zj///Px/jxkzJiZMmBCjRo2Kf//3f49zzz13j9eZA73P5ZdfHs8++2w88cQTu7zme6A47GkO+B4obO9///ujvr4+/vKXv8RDDz0UF154YdTV1eVf9/kvfHuaA8cee6zPPwVp8uTJ+Q2Wd+euu+7ape+0006LZ5555l0994BaxvXe9743+vTps9vNiN6Z5acwDR48OD74wQ/Gyy+/nD+Vy3woDh15vysrK6OlpSXeeuutPY6hsAwfPjxGjRoVL7/8ckSYA4XiiiuuiMceeyx+9atfxYgRI/L9vgeKx57mwO74Higs/fv3j/e9730xYcKEqK2tjbFjx8Z3vvMdn/8isqc5sDs+/9B5B1Syp3///jF+/PhYunRpu/6lS5fGpEmTeigq9qfm5uZ44YUXYvjw4VFdXR2VlZXt5kNLS0vU1dWZDwWoI+/3+PHjo1+/fu3GbNiwIX73u9+ZEwXqzTffjPXr18fw4cMjwhzo7ZIkicsvvzwefvjh+OUvfxnV1dXtXvc9UPj2NQd2x/dAYUuSJJqbm33+i9jOObA7Pv/wLqTaznk/eOCBB5J+/fold9xxR/L8888ns2fPTgYPHpz84Q9/6OnQ6AZf+tKXkmXLliWvvPJK8vTTTycf//jHk9LS0vz7/Y1vfCMpLy9PHn744eS5555LPvOZzyTDhw9PGhsbezhyOqOpqSlZvXp1snr16iQikgULFiSrV69OXn311SRJOvZ+z5w5MxkxYkTyH//xH8kzzzyTnHHGGcnYsWOT1tbWnvq1SGFvc6CpqSn50pe+lDz55JPJ2rVrk1/96lfJxIkTk0MPPdQcKBD/5//8n6S8vDxZtmxZsmHDhnzbunVrfozvgcK2rznge6CwzZ07N1m+fHmydu3a5Nlnn02+/OUvJyUlJcmSJUuSJPH5LwZ7mwM+/9C1DrhkT5Ikyfe+971k1KhRSf/+/ZMPfehD7Y7jpLCcf/75yfDhw5N+/folVVVVybnnnpusWbMm/3oul0tuvPHGpLKyMslms8mpp56aPPfccz0YMe/Gr371qyQidmkXXnhhkiQde7+3bduWXH755cnQoUOTgQMHJh//+MeTdevW9cBvQ2fsbQ5s3bo1mTp1anLIIYck/fr1Sw477LDkwgsv3OX9NQd6r9299xGR3HnnnfkxvgcK277mgO+BwnbJJZfk/45/yCGHJB/96EfziZ4k8fkvBnubAz7/0LUySbKXnYIAAAAA6FUOqD17AAAAAHh3JHsAAAAACohkDwAAAEABkewBAAAAKCCSPQAAAAAFRLIHAAAAoIBI9gAAAAAUEMkeAAAAgAIi2QMAAABQQCR7AAAAAAqIZA8AAABAAZHsAQAAACgg/z/z5aQ5BMKuTwAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "nx = 400\n", "arguments = InitialConditions.genKelvinHelmholtz(nx, nx//4, 1.4)\n", "\n", "plt.figure()\n", "plt.imshow(arguments['rho'])\n", "plt.colorbar(orientation='vertical', aspect=20, pad=0.02, shrink=0.3)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Initialized /home/smyalygames/PycharmProjects/FiniteVolumeGPU/data/euler_kelvin_helmholtz.nc\n", "Opening /home/smyalygames/PycharmProjects/FiniteVolumeGPU/data/euler_kelvin_helmholtz.nc\n", "Keyword arguments: {'mode': 'w', 'clobber': False}\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Constructed in 0.001973867416381836 seconds\n", "Simulating to t=10.000000. Estimated 4346 timesteps (dt=0.002301)\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "Closing /home/smyalygames/PycharmProjects/FiniteVolumeGPU/data/euler_kelvin_helmholtz.nc\n" ] } ], "source": [ "nx = 400\n", "ny = nx//2\n", "roughness = 0.125\n", "t_end = 10.0\n", "n_save = 100#1000\n", "outfile = \"data/euler_kelvin_helmholtz.nc\"\n", "outdata = None\n", "\n", "arguments = InitialConditions.genKelvinHelmholtz(nx, ny, gamma, roughness)\n", "arguments['context'] = my_context\n", "outfile = runSimulation(outfile, t_end, arguments)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Initialized /home/smyalygames/PycharmProjects/FiniteVolumeGPU/data/euler_kelvin_helmholtz.nc\n", "Opening /home/smyalygames/PycharmProjects/FiniteVolumeGPU/data/euler_kelvin_helmholtz.nc\n", "Arguments: ('r',)\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "/home/smyalygames/PycharmProjects/FiniteVolumeGPU/data/euler_kelvin_helmholtz.nc created Wed Feb 12 16:40:04 2025 contains 101 timesteps\n", "Simulator arguments: \n", " {'rho': '[...]', 'rho_u': '[...]', 'rho_v': '[...]', 'E': '[...]', 'nx': 400, 'ny': 200, 'dx': 0.005, 'dy': 0.005, 'g': 0.0, 'gamma': 1.4, 'boundary_conditions': '[north=Type.Periodic, south=Type.Periodic, east=Type.Periodic, west=Type.Periodic]', 'context': 'CudaContext id 96726128804080'}\n", "0% [######################========] 100%. Total: 6s, elapsed: 5s, remaining: 1s\n" ] }, { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "
\n", " \n", "
\n", " \n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "
\n", "
\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Closing /home/smyalygames/PycharmProjects/FiniteVolumeGPU/data/euler_kelvin_helmholtz.nc\n" ] } ], "source": [ "#outfile='data/euler_kelvin_helmholtz_0012.nc'\n", "createAnimation(outfile, vis_type=VisType.Density, vmin=1.0, vmax=2.0, save_anim=False, fig_args={'figsize':(16, 9)})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Rayleigh-Taylor\n", "\n", "* Liska and Wendroff, COMPARISON OF SEVERAL DIFFERENCE SCHEMES ON 1D AND 2D TEST PROBLEMS FOR THE EULER EQUATIONS, http://www-troja.fjfi.cvut.cz/~liska/CompareEuler/compare8.pdf\n", "* https://www.astro.princeton.edu/~jstone/Athena/tests/rt/rt.html\n", "* " ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAABQsAAANJCAYAAACrg2qeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAABbGklEQVR4nO3de3xdZZ0v/s/aSZsCthEobahcRE8FmSJiUQRUcIDizCA6N3DgoDOH8cBBcDqIAuqM4E9bxBF9KYIyhwOOg9SfIuqcQaRzsYqgYm1HLg5e6HBtrfIraQu90azfH9lJm7Zps9MkO8l6v3mtV5K1n737PDt5vqx8813PU5RlWQYAAAAAqLxaszsAAAAAAIwOkoUAAAAAQBLJQgAAAACgTrIQAAAAAEgiWQgAAAAA1EkWAgAAAABJJAsBAAAAgDrJQgAAAAAgSdLa7A4Ml66urjz11FOZPHlyiqJodneAEVSWZdasWZMZM2akVvM3EfEQqkks7EsshOoSD/sSD6GaGomF4zZZ+NRTT+XAAw9sdjeAJnr88cdzwAEHNLsbTSceQrWJhd3EQkA87CYeQrUNJBaO22Th5MmTkySP/uTFmfICfz2CKlm9tisHv+q/euNA1fW8Dwdc8cHUJk1qcm+AkdK1fn2euOIjYmGdWAjVJR721fM+vC6/n9ZMaHJvgJHyfDbl7twxoFg4bpOFPeXUU15Qy5TJkoVQRW6r6NbzPtQmTfILMlSQWNhNLATEw24970NrJqS1kCyEyii7PwwkFo7bZCEAAADQj6LoPoCKKHoThrui5A4AAAAASKKyEAAAAKqnqHUfQEXUVBYCAAAAAI1RWQgAAAAVU9QKm75AhRRlkXQNrK3KQgAAAAAgicpCAAAAqB5rFkLFDHy+iwwAAAAAQBLJQgAAAACgzm3IAAAAUDW1IrHBCVRHOfD5rrIQAAAAAEiishAAAAAqpyiSQmUhVEYjs11lIQAAAACQRGUhAAAAVE+tlhTqh6AyyoHPd5EBAAAAAEiishAAAACqp7AbMlSL3ZABAAAAgAZJFgIAAAAASdyGDAAAANVjgxOoFhucAAAAAACNUlkIAAAAFVMURQobnEBlFDY4AQAAAAAapbIQAAAAqqZWJCoLoTpKlYUAAAAAQINUFgIAAEDVFHZDhmqxGzIAAAAA0CDJQgAAAAAgiduQAQAAoHqKonuTE6AaumxwAgAAAAA0SGUhAAAAVE1RdB9ANTQw31UWAgAAAABJVBYCAABA9RS17gOohgbmu8gAAAAAACRRWQgAAADVU7MbMlSLNQsBAAAAgAZJFgIAAAAASdyGDAAAAJVTFEWKwm3IUBWNzHeVhQAAAABAEpWFAAAAUD1F0X0A1aCyEAAAAABolMpCAAAAqJpakdTUD0F1qCwEAAAAABqkshAAAACqxpqFUC3WLAQAAAAAGqWyEAAAAKpGZSFUi8pCAAAAAKBRkoUAAAAAQBK3IQMAAED11IruA6gItyEDAAAAAA1SWQgAAABVY4MTqBYbnAAAAAAAjVJZCAAAAFWjshCqRWUhAAAAANAolYUAAABQNSoLoVpUFgIAAAAAjZIsBAAAAACSuA0ZAAAAKqcsipQ1tyFDVZSl25ABAAAAgAapLAQAAICqscEJVIsNTgAAAACARqksBAAAgKpRWQjVorIQAAAAAGiUykIAAACoGpWFUC0qCwEAAACARkkWAgAAAABJ3IYMAAAA1VMrug+gGkq3IQMAAAAADVJZCAAAAFVjgxOoFhucAAAAAACNUlkIAAAAVaOyEKpFZSEAAAAA0CiVhQAAAFAxZVGkVFkIldHIfFdZCAAAAAAkkSwEAAAAAOrchgwAAABVU4vyIaiScuBNhQYAAAAAIInKQgAAAKieoug+gGqwwQkAAAAA0CiVhQAAAFA1KguhWlQWAgAAAACNUlkIAAAAFVMWRUqVhVAZjcx3lYUAAAAAQBLJQgAAAACgzm3IAAAAUDW1KB+CKikH3lRoAAAAAACSqCwEAACA6imS2OAEqqOB6a6yEAAAAABIorIQAAAAqqcoVBZClTQw31UWAgAAAABJVBYCAABA5ZRJSoWFUBkNbIasshAAAAAA6CZZCAAAAAAkcRsyAAAAVE+t6D6AamhgvqssBAAAAACSqCwEAACAyimLImWhshCqopH5rrIQAAAAAEiishAAAACqp6gfQDU0MN9VFgIAAAAASVQWAgAAQPUURfcBVIM1CwEAAACARkkWAgAAAABJ3IYMAAAAlVMW3QdQDY3Md5WFAAAAAEASlYUAAABQPTY4gWqxwQkAAAAA0CiVhQAAAFA1tSgfgippYL4LDQAAAABAEpWFAAAAUDllUaS0ZiFURiPzXWUhAAAAAJBEshAAAAAAqHMbMgAAAFRNUT+AamhgvqssBAAAAACSqCwEAACAyrHBCVSLDU4AAAAAgIapLAQAAICqsWYhVIs1CwEAAACARqksBAAAgIopi+4DqIZG5rvKQgAAAAAgiWQhAAAAAFDnNmQAAAComlrRfQDV0MB8b7iy8Lvf/W7e/OY3Z8aMGSmKIl//+tf7PF6WZa644orMmDEje+yxR0488cQ8+OCDfdps2LAhF110UaZOnZq99torp59+ep544ok+bVatWpVzzjkn7e3taW9vzznnnJNnnnmm0e4CDAuxEKCbeAggFgLjS8PJwmeffTZHHnlkrr322h0+fvXVV+eaa67Jtddem/vuuy8dHR055ZRTsmbNmt42c+fOze23354FCxbk7rvvztq1a3Paaadl8+bNvW3OOuusLF26NHfeeWfuvPPOLF26NOecc84ghggw9MRCgG7iIcDYjIU9G5w4HI7qHANVlGVZDiqyJCmKIrfffnve+ta3dgebssyMGTMyd+7cXHrppUm6/zoyffr0fOxjH8t5552Xzs7O7LfffvniF7+YM888M0ny1FNP5cADD8wdd9yRU089NT/72c9y+OGH5wc/+EGOOeaYJMkPfvCDHHvssfnP//zPHHroobvs2+rVq9Pe3p5VP39Jpky2NCNUyeo1Xdn7ZY+ks7MzU6ZMGfZ/bzTHwmRLPDzoqo+kNmnS0L8BwKjUtX59HrvsgyMWC5PRHQ/FQqiukY6HozkWJlvi4XEnX5HWVvEQquL559fnnn+5YkCxcEizaMuWLcuKFSsyZ86c3nNtbW054YQTcs899yRJFi9enE2bNvVpM2PGjMyaNau3zb333pv29vbeAJgkr33ta9Pe3t7bZlsbNmzI6tWr+xwAzdDMWJiIh8Do4doQYDRfGxZJ4XA4KnOkGHDcGtJk4YoVK5Ik06dP73N++vTpvY+tWLEiEydOzN57773TNtOmTdvu9adNm9bbZlvz58/vXbehvb09Bx544G6PB2AwmhkLE/EQGD1cGwK4NgTGnmG5P7co+mYry7Lc7ty2tm2zo/Y7e53LL788nZ2dvcfjjz8+iJ4DDJ1mxMJEPARGH9eGAKPv2rDZa6c5HI6RPwZqSJOFHR0dSbLdXzVWrlzZ+1eUjo6ObNy4MatWrdppm1//+tfbvf5vfvOb7f4a06OtrS1TpkzpcwA0QzNjYSIeAqOHa0MA14bA2DOkycJDDjkkHR0dWbhwYe+5jRs3ZtGiRTnuuOOSJLNnz86ECRP6tFm+fHkeeOCB3jbHHntsOjs786Mf/ai3zQ9/+MN0dnb2tgEYrcRCgG7iIYBYCIw9rY0+Ye3atfnlL3/Z+/WyZcuydOnS7LPPPjnooIMyd+7czJs3LzNnzszMmTMzb9687LnnnjnrrLOSJO3t7Tn33HPznve8J/vuu2/22WefXHLJJTniiCNy8sknJ0le/vKX501velPe+c535vOf/3yS5H/+z/+Z0047bcA7PAEMJ7EQoJt4CDBGY2GRRvY7AMa6BuZ7w8nCH//4x3njG9/Y+/XFF1+cJHnHO96Rm2++Oe973/uybt26XHDBBVm1alWOOeaY3HXXXZk8eXLvcz75yU+mtbU1Z5xxRtatW5eTTjopN998c1paWnrb3HLLLXn3u9/duxvU6aefnmuvvbbR7gIMC7EQoJt4CCAWAuNLUZZl2exODIfVq1envb09q37+kkyZPCz7uACj1Oo1Xdn7ZY+ks7PTmizZEg8PuuojqU2a1OzuACOka/36PHbZB8XCOrEQqks87KsnHh576pVpnSAeQlU8v2l97v32hwYUC2XRAAAAAIAkg7gNGQAAABjbyqJIWVi0EKqikfmushAAAAAASKKyEAAAAKrHbshQLQ3Md5WFAAAAAEASyUIAAAAAoM5tyAAAAFAxZa37AKqhkfkuNAAAAAAASVQWAgAAQPXY4ASqxQYnAAAAAECjVBYCAABAxZRFkbJQWghV0ch8V1kIAAAAACRRWQgAAADVY81CqBZrFgIAAAAAjZIsBAAAAACSuA0ZAAAAKqcsug+gGhqZ7yoLAQAAAIAkKgsBAACgemxwAtWishAAAAAAaJTKQgAAAKgYaxZCtVizEAAAAABomMpCAAAAqJqi6D6AamhgvqssBAAAAACSqCwEAACAyrFmIVSLNQsBAAAAgIZJFgIAAAAASdyGDAAAANVT1A+gGtyGDAAAAAA0SmUhAAAAVI0NTqBaVBYCAAAAAI1SWQgAAABVY81CqBaVhQAAAABAo1QWAgAAQMWUsWYhVEnZQFuVhQAAAABAEslCAAAAAKDObcgAAABQNTY4gWqxwQkAAAAA0CiVhQAAAFAxZVGkLJQWQlU0Mt9VFgIAAAAASVQWAgAAQPVYsxCqxZqFAAAAAECjVBYCAABAxZRF9wFUQyPzXWUhAAAAAJBEshAAAAAAqHMbMgAAAFSNDU6gWtyGDAAAAAA0SmUhAAAAVIwNTqBabHACAAAAADRMZSEAAABUjTULoVpUFgIAAAAAjVJZCAAAABVjzUKoFmsWAgAAAAANkywEAAAAAJK4DRkAAACqxwYnUC1uQwYAAAAAGqWyEAAAAKpGZSFUi8pCAAAAAKBRKgsBAACgYsqi+wCqoZH5rrIQAAAAAEiishAAAACqx5qFUC0qCwEAAACARkkWAgAAAABJ3IYMAAAAlWODE6gWG5wAAAAAAA1TWQgAAABVY4MTqBaVhQAAAABAo1QWAgAAQMWURZGyUFoIVdHIfFdZCAAAAAAkUVkIAAAA1WPNQqgWaxYCAAAAAI2SLAQAAAAAkrgNGQAAAKrHbchQLW5DBgAAAAAapbIQAAAAKqZMUqoshMooG2irshAAAAAASKKyEAAAAKrHmoVQLdYsBAAAAAAapbIQAAAAKqYsrFkIVdLIfFdZCAAAAAAkkSwEAAAAAOrchgwAAABVY4MTqBa3IQMAAAAAjVJZCAAAAFWjshCqRWUhAAAAANAolYUAAABQMWXRfQDV0Mh8V1kIAAAAACRRWQgAAADVY81CqBaVhQAAAABAoyQLAQAAAIAkbkMGAACAyrHBCVSLDU4AAAAAgIapLAQAAICqscEJVIvKQgAAAACgUSoLAQAAoGpUFkK1qCwEAAAAABqlshAAAAAqxm7IUC12QwYAAAAAGiZZCAAAAAAkcRsyAAAAVI8NTqBa3IYMAAAAADRKZSEAAABUjA1OoFpscAIAAAAANEyyEAAAAABIIlkIAAAAANRZsxAAAACqxm7IUC3WLAQAAAAAGiVZCAAAAAAkcRsyAAAAVE5ZdB9ANTQy31UWAgAAAABJVBYCAABA9djgBKpFZSEAAAAA0CiVhQAAAFA1KguhWlQWAgAAAACNUlkIAAAAFWM3ZKgWuyEDAAAAAA2TLAQAAAAAkrgNGQAAAKrJbcjADgx5ZeHzzz+fD37wgznkkEOyxx575CUveUk+/OEPp6urq7dNWZa54oorMmPGjOyxxx458cQT8+CDD/Z5nQ0bNuSiiy7K1KlTs9dee+X000/PE088MdTdBRgWYiFAN/EQQCwExpYhTxZ+7GMfy+c+97lce+21+dnPfparr746H//4x/OZz3ymt83VV1+da665Jtdee23uu+++dHR05JRTTsmaNWt628ydOze33357FixYkLvvvjtr167Naaedls2bNw91lwGGnFgI0E08BBilsbBwOByVOwaoKMuyHHjzXTvttNMyffr03Hjjjb3n/viP/zh77rlnvvjFL6Ysy8yYMSNz587NpZdemqT7ryPTp0/Pxz72sZx33nnp7OzMfvvtly9+8Ys588wzkyRPPfVUDjzwwNxxxx059dRTd9mP1atXp729Pat+/pJMmWxpRqiS1Wu6svfLHklnZ2emTJnSlD6MlliYbImHB131kdQmTRr6wQKjUtf69Xnssg82NRYmoyceioVQXaMhHo6WWJhsiYcvvXxeWsRDqIzN69fnV/PfP6BYOORZtNe97nX513/91/z85z9PkvzHf/xH7r777vz+7/9+kmTZsmVZsWJF5syZ0/uctra2nHDCCbnnnnuSJIsXL86mTZv6tJkxY0ZmzZrV22ZbGzZsyOrVq/scAM3SrFiYiIfA6OLaEGB0XhuWhcPhqNoxUEO+wcmll16azs7OHHbYYWlpacnmzZvz0Y9+NH/2Z3+WJFmxYkWSZPr06X2eN3369Dz66KO9bSZOnJi99957uzY9z9/W/Pnzc+WVVw71cAAGpVmxMBEPgdHFtSGAa0NgbBnyysIvf/nL+cd//Md86Utfyk9+8pN84QtfyN/93d/lC1/4Qp92RdE3pVmW5XbntrWzNpdffnk6Ozt7j8cff3z3BgKwG5oVCxPxEBhdXBsCjNJrw6J0OBxVOwZoyCsL3/ve9+ayyy7L2972tiTJEUcckUcffTTz58/PO97xjnR0dCTp/qvI/vvv3/u8lStX9v4VpaOjIxs3bsyqVav6/NVk5cqVOe6443b477a1taWtrW2ohwMwKM2KhYl4CIwurg0BXBsCY8uQVxY+99xzqdX6vmxLS0vvlvCHHHJIOjo6snDhwt7HN27cmEWLFvUGuNmzZ2fChAl92ixfvjwPPPDAToMgwGghFgJ0Ew8BxEJgbBnyysI3v/nN+ehHP5qDDjoov/M7v5MlS5bkmmuuyf/4H/8jSVIURebOnZt58+Zl5syZmTlzZubNm5c999wzZ511VpKkvb095557bt7znvdk3333zT777JNLLrkkRxxxRE4++eSh7jLAkBMLAbqJhwCjNBY2uOEBMMY1c4OTz3zmM/mbv/mbXHDBBVm5cmVmzJiR8847L3/7t3/b2+Z973tf1q1blwsuuCCrVq3KMccck7vuuiuTJ0/ubfPJT34yra2tOeOMM7Ju3bqcdNJJufnmm9PS0jLUXQYYcmIhQDfxEEAsBMaWoizLga9wOIasXr067e3tWfXzl2TK5CG/2xoYxVav6creL3sknZ2dmTJlSrO703Q98fCgqz6S2qRJze4OMEK61q/PY5d9UCysEwuhusTDvnri4Uv+5qPiIVRI1/r1eeT/+cCAYqEsGgAAAACQZBhuQwYAAABGuSINrWEGjHENzHeVhQAAAABAEpWFAAAAUDml3ZChUhqZ7yoLAQAAAIAkKgsBAACgeqxZCNWishAAAAAAaJRkIQAAAACQxG3IAAAAUD1uQ4ZqcRsyAAAAANAolYUAAABQMWX9AKqhkfmushAAAAAASKKyEAAAAKrHmoVQLdYsBAAAAAAapbIQAAAAqkZlIVSLykIAAAAAoFGShQAAAABAErchAwAAQOWURfcBVEMj811lIQAAAACQRGUhAAAAVE9Rdh9ANTQw31UWAgAAAABJVBYCAABA9RT1A6gGaxYCAAAAAI1SWQgAAAAVYzdkqBa7IQMAAAAADZMsBAAAAACSuA0ZAAAAqscGJ1AtbkMGAAAAABqlshAAAACqRmUhVIvKQgAAAACgUSoLAQAAoGLKJKXKQqiMsoG2KgsBAAAAgCQqCwEAAKB6irL7AKqhgfmushAAAAAASCJZCAAAAADUuQ0ZAAAAqqaoH0A1NDDfVRYCAAAAAElUFgIAAED1qCyEalFZCAAAAAA0SmUhAAAAVI3KQqgWlYUAAAAAQKNUFgIAAEDFlEWZsiib3Q1ghDQy31UWAgAAAABJJAsBAAAAgDq3IQMAAEDV2OAEqsUGJwAAAABAo1QWAgAAQNWoLIRqUVkIAAAAADRKshAAAAAASCJZCAAAAADUWbMQAAAAqqYouw+gGhqY7yoLAQAAAIAkkoUAAAAAQJ3bkAEAAKBqivoBVEMD811lIQAAAACQRGUhAAAAVI/KQqgWlYUAAAAAQKNUFgIAAEDVFGX3AVRDA/NdZSEAAAAAkERlIQAAAFSPNQuhWqxZCAAAAAA0SrIQAAAAAEjiNmQAAACoHrchQ7W4DRkAAAAAaJTKQgAAAKiaouw+gGpoYL6rLAQAAAAAkqgsBAAAgMopiu4DqIZG5rvKQgAAAAAgicpCAAAAqB5rFkK1WLMQAAAAAGiUZCEAAAAAkMRtyAAAAFA9Rf0AqsEGJwAAAABAo1QWAgAAQMUUKVPY4AQqo4gNTgAAAACABqksBAAAgKqxZiFUizULAQAAAIBGqSwEAACAiimK7gOohkbmu8pCAAAAACCJZCEAAAAAUOc2ZAAAAKiaouw+gGpoYL6rLAQAAAAAkqgsBAAAgMop6gdQDY3Md5WFAAAAAEASlYUAAABQOUVRprBmIVRGI/NdZSEAAAAAkERlIQAAAFSP3ZChWlQWAgAAAACNkiwEAAAAAJK4DRkAAAAqpyi6D6AaGpnvKgsBAAAAgCQqCwEAAKB6ijKFDU6gOmxwAgAAAAA0SmUhAAAAVIw1C6FarFkIAAAAADRMZSEAAABUTGHNQqiURua7ykIAAAAAIIlkIQAAAABQ5zZkAAAAqJjuDU7chgxVYYMTAAAAAKBhKgsBAACgYrorC5vdC2CkqCwEAAAAABqmshAAAAAqpihKaxZChTQy31UWAgAAAABJVBYCAABA5agshGpRWQgAAAAANEyyEAAAAABI4jZkAAAAqJyifgDV0Mh8V1kIAAAAACRRWQgAAACVUyvK1GxwApVR2uAEAAAAAGiUykIAAAComKIoU6gshMpoZL6rLAQAAAAAkqgsBAAAgMpRWQjVorIQAAAAAGiYZCEAAAAAkMRtyDTB2q71ea7cPKC2x9wxNy/4xYRh7lFfaw/dlB++6VMDartn0ZIX1CYNb4cAAACGWK3oPhi7Rvo28rL0AzOWNfLtG5Zk4ZNPPplLL7003/rWt7Ju3bq87GUvy4033pjZs2d3d7Asc+WVV+aGG27IqlWrcswxx+Szn/1sfud3fqf3NTZs2JBLLrkkt956a9atW5eTTjop1113XQ444IDh6DJD6Lvrk0sf/uN+H++6dVr2/favBvRah/5//5Fy08ah6tqAFBMm5s/3+cMBtf3t7700LWeu7Pfxvzv0Kzl+kgLeqhILAbqJhwBiId3G8jqRQ9l3icfRbciThatWrcrxxx+fN77xjfnWt76VadOm5Ve/+lVe+MIX9ra5+uqrc8011+Tmm2/Oy172snzkIx/JKaeckocffjiTJ09OksydOzf/9E//lAULFmTffffNe97znpx22mlZvHhxWlpahrrb7MTKzc/mv56fuN35P7v9okxduv0E3/PXmzLl2z/eySv+KgOrK2yOctPGbP51/wnAre1988rk5v4f/+Cbzstz07afZr89qsytb/3Mdudf3Lox01r2GmhXGcXEQoBu4iHA6IyFNjihmfzsjbxG3vOiLMsh/Q5ddtll+f73v5/vfe97O3y8LMvMmDEjc+fOzaWXXpqk+68j06dPz8c+9rGcd9556ezszH777ZcvfvGLOfPMM5MkTz31VA488MDccccdOfXUU3fZj9WrV6e9vT2rfv6STJmssmsgfrRhU972vfO2O7/P99oy9R9+st35cuPGZGh/fKqjKFJM3D4B+5s/f1VWHb9hu/MLXv/5vKZtZG/HHstWr+nK3i97JJ2dnZkyZUpT+jBaYmGyJR4edNVHUpvktnmoiq716/PYZR9saixMRk88FAuhukZDPBwtsTDZEg9f+dWL07Jn29AMsKJqEl4jqks14m7Z/NyGLP2TawYUC4e8svCb3/xmTj311Pzpn/5pFi1alBe96EW54IIL8s53vjNJsmzZsqxYsSJz5szpfU5bW1tOOOGE3HPPPTnvvPOyePHibNq0qU+bGTNmZNasWbnnnnt2GAQ3bNiQDRu2JFlWr1491EMbs3628bn8cP2L+5z79M/fmPbr+v5wTFizKTO/v31SMEmEwCFWlik3bJ8UnPr5ezP189s3/+Dx78ymyX2Thavf1ZmLZn6nz7ljJv1XXj5xz6HsKYPUrFiYiIfA6OLaEGB0XhvWUkp2Mab4ed09ZQOZnSFPFj7yyCO5/vrrc/HFF+f9739/fvSjH+Xd73532tra8va3vz0rVqxIkkyfPr3P86ZPn55HH300SbJixYpMnDgxe++993Ztep6/rfnz5+fKK68c6uGMOaf/4k158Ccv7nNu+g+TyQt+0Ofcfnl4BHvF7iq+vzTb1iFOvTO5NTP6nPv8n/1RVr6mb7sjXrUsX5/57eHtINtpVixMxENgdHFtCODaEBhbhjxZ2NXVlaOPPjrz5s1Lkhx11FF58MEHc/311+ftb397b7ui6Fs+Wpbldue2tbM2l19+eS6++OLer1evXp0DDzxwsMMYdT7+/700T27o+z+F7/7vV2faD/r+lbz21G/y337dNzFIdUy59QeZcmvfcxunT8ubZpzd59yvj23PCef+qM+5F7Wtynv3GdjGM+xas2JhMv7jITC2uDYEGJ3XhuN5zUKbZ7A7xuu8aGRcQ54s3H///XP44Yf3Offyl788t912W5Kko6MjSfdfRfbff//eNitXruz9K0pHR0c2btyYVatW9fmrycqVK3Pcccft8N9ta2tLW9vYWW9hc9mVdeX2u/z+7n+ck7X37Lfd+UP+4bE8//gTfc7tl3u3KyIdzRuH0Bybf70y2WbDlmlLkp9d17fdLw58eW56+/a3LrS/7tdZeMQt253fo5iYlsJ6oP1pVixMxl48BMY314YArg1H2nhN9sBIGfJk4fHHH5+HH+57i+vPf/7zHHzwwUmSQw45JB0dHVm4cGGOOuqoJMnGjRuzaNGifOxjH0uSzJ49OxMmTMjChQtzxhlnJEmWL1+eBx54IFdfffVQd3lYzV1+dB58Zv/tzj/yHy/KYZ98Yrvz+z6zInuv+cV2558flt7BFs8//kQO/Oj2P5O1yZNz5gv/ZLvz/3nxAXnJK57c7vwrXvhkPrH/jte+rBKxEKCbeAgwOmNhUTR3zUKbVVBFTZ1zzaws/Ou//uscd9xxmTdvXs4444z86Ec/yg033JAbbrghSXdZ9dy5czNv3rzMnDkzM2fOzLx587LnnnvmrLPOSpK0t7fn3HPPzXve857su+++2WeffXLJJZfkiCOOyMknnzzUXd7Oc10b8+jzA0vPXfjLt+WZ//dF/T7ecdeTqf3XY9ud/295XAKQMaFrzZp0rVmz3fn/9tfbJxaT5D8POTizTzmm39fb58wn8umXfnlA//bBra3Zs7b9rtFjwXiIhQBDQTwEqG4slBCEvnY2J0bTBi5Dnix89atfndtvvz2XX355PvzhD+eQQw7Jpz71qZx99pY10973vvdl3bp1ueCCC7Jq1aocc8wxueuuuzJ58uTeNp/85CfT2tqaM844I+vWrctJJ52Um2++OS0tLQ31Z/Y9b0ttz0kNPad4eK8c8vEHBtR2wsZfZ+qG7ZOBPSQEqZrnlz2aqTc82u/jxRfacvHEHe/Utq1HLp2VzHy24T50Pbc+ybyGnzeURlssBGgW8RBgfMVCCUAYHo3MreFOLBZlWY6e1OUQWr16ddrb23Ni3pLWYkKzuwOMoOfLTflOvpHOzs5MmTKl2d1pup54eNBVH0ltUmN/PAHGrq716/PYZR8UC+vEQqgu8bCvnng4+7a5admrWmsZQpVtfnZDFv/xpwYUC4e8shAAAAAY3Zq9ZiEwshpZs9BWpgAAAABAEslCAAAAAKBOshAAAAAASGLNQgAAAKicon4A1dDIfFdZCAAAAAAkkSwEAAAAAOokCwEAAACAJNYsBAAAgMqpFWVqRdnsbgAjpJH5rrIQAAAAAEgiWQgAAAAA1EkWAgAAAABJJAsBAAAAgDobnAAAAEDF2OAEqsUGJwAAAABAwyQLAQAAAIAkkoUAAAAAQJ01CwEAAKBiiqJMYc1CqIxG5rvKQgAAAAAgiWQhAAAAAFAnWQgAAAAAJJEsBAAAAADqbHACAAAAFVMrytRscAKV0ch8V1kIAAAAACSRLAQAAAAA6iQLAQAAAIAk1iwEAACAyrFmIVSLNQsBAAAAgIZJFgIAAAAASSQLAQAAAIA6yUIAAAAAIIkNTgAAAKByailTiw1OoCoame8qCwEAAACAJJKFAAAAAECdZCEAAAAAkMSahQAAAFA5RVGmKKxZCFXRyHxXWQgAAAAAJJEsBAAAAADqJAsBAAAAgCSShQAAAABAnQ1OAAAAoGJqRZmaDU6gMhqZ7yoLAQAAAIAkkoUAAAAAQJ1kIQAAAACQxJqFAAAAUDlF0dgaZsDYVhQDb6uyEAAAAABIIlkIAAAAANRJFgIAAAAASSQLAQAAAIA6G5wAAABAxdSK0gYnUCGNzHeVhQAAAABAEslCAAAAAKBOshAAAAAASGLNQgAAAKicWsrUYs1CqIpG5rvKQgAAAAAgiWQhAAAAAFAnWQgAAAAAJJEsBAAAAADqbHACAAAAFVMrulIruprdDWCENDLfVRYCAAAAAEkkCwEAAACAOslCAAAAACCJNQsBAACgcmpFmVpRNrsbwAhpZL6rLAQAAAAAkkgWAgAAAAB1koUAAAAAQBLJQgAAAACgzgYnAAAAUDE2OIFqscEJAAAAANAwyUIAAAAAIIlkIQAAAABQZ81CAAAAqBhrFkK1WLMQAAAAAGiYZCEAAAAAkESyEAAAAACokywEAAAAAJLY4AQAAAAqp5YytdjgBKqikfmushAAAAAASCJZCAAAAADUSRYCAAAAAEmsWQgAAACVUyvK1AprFkJVNDLfVRYCAAAAAEkkCwEAAACAOslCAAAAACCJZCEAAAAAUGeDEwAAAKiYWtGVWtHV7G4AI6SR+a6yEAAAAABIIlkIAAAAANRJFgIAAAAASaxZCAAAAJVTS5mWomx2N4ARsjkDn+8qCwEAAACAJJKFAAAAAECdZCEAAAAAkESyEAAAAACos8EJAAAAVEwtZWoNbHgAjG2NzHeVhQAAAABAEslCAAAAAKBOshAAAAAASGLNQgAAAKicWlGmVlizEKqikfmushAAAAAASCJZCAAAAADUSRYCAAAAAEkkCwEAAACAOhucAAAAQMXUiq7Uiq5mdwMYIY3Md5WFAAAAAEASyUIAAAAAoE6yEAAAAABIYs1CAAAAqJyWokxLUTa7G8AIaWS+qywEAAAAAJJIFgIAAAAAdZKFAAAAAEASyUIAAAAAoM4GJwAAAFAxtSS12OAEqqKRakGVhQAAAABAEslCAAAAAKBOshAAAAAASDICycL58+enKIrMnTu391xZlrniiisyY8aM7LHHHjnxxBPz4IMP9nnehg0bctFFF2Xq1KnZa6+9cvrpp+eJJ54Y7u4CDAuxEKCbeAgwOmJhrehyOBwVOwYcHwYVVQbovvvuyw033JBXvOIVfc5fffXVueaaa3LttdfmvvvuS0dHR0455ZSsWbOmt83cuXNz++23Z8GCBbn77ruzdu3anHbaadm8efNwdhlgyImFAN3EQwCxEBj9hi1ZuHbt2px99tn5+7//++y9996958uyzKc+9al84AMfyB/90R9l1qxZ+cIXvpDnnnsuX/rSl5IknZ2dufHGG/OJT3wiJ598co466qj84z/+Y+6///78y7/8y3B1GWDIiYUA3cRDALEQGBuGLVn4rne9K3/wB3+Qk08+uc/5ZcuWZcWKFZkzZ07vuba2tpxwwgm55557kiSLFy/Opk2b+rSZMWNGZs2a1dsGYCwQCwG6iYcAYiEwNrQOx4suWLAgP/nJT3Lfffdt99iKFSuSJNOnT+9zfvr06Xn00Ud720ycOLHPX1p62vQ8f1sbNmzIhg0ber9evXr1bo0BYHc1IxYm4iEw+rg2BHBtCIwdQ54sfPzxx/NXf/VXueuuuzJp0qR+2xVF0efrsiy3O7etnbWZP39+rrzyysY7DDAMmhULE/EQGF1cGwKMzmvDWlGmVpS76DkwXjQy34f8NuTFixdn5cqVmT17dlpbW9Pa2ppFixbl05/+dFpbW3v/UrLtXz5WrlzZ+1hHR0c2btyYVatW9dtmW5dffnk6Ozt7j8cff3yohwYwYM2KhYl4CIwurg0BXBsCY8uQJwtPOumk3H///Vm6dGnvcfTRR+fss8/O0qVL85KXvCQdHR1ZuHBh73M2btyYRYsW5bjjjkuSzJ49OxMmTOjTZvny5XnggQd622yrra0tU6ZM6XMANEuzYmEiHgKji2tDANeGwNgy5LchT548ObNmzepzbq+99sq+++7be37u3LmZN29eZs6cmZkzZ2bevHnZc889c9ZZZyVJ2tvbc+655+Y973lP9t133+yzzz655JJLcsQRR2y3ECzAaCQWAnQTDwHEQmBsGZYNTnblfe97X9atW5cLLrggq1atyjHHHJO77rorkydP7m3zyU9+Mq2trTnjjDOybt26nHTSSbn55pvT0tLSjC4DDDmxEKCbeAgw8rGwlq60pGsohwCMYrUG5ntRluW4XNF09erVaW9vz4l5S1qLCc3uDjCCni835Tv5Rjo7O91mkS3x8KCrPpLaThbUBsaXrvXr89hlHxQL68RCqC7xsK+eePi/vvuHaXuB35WhKjas3ZTr33D7gGLhkK9ZCAAAAACMTZKFAAAAAEASyUIAAAAAoK4pG5wAAAAAzVMrytSKcbmFAbADjcx3lYUAAAAAQBLJQgAAAACgTrIQAAAAAEhizUIAAACoHGsWQrVYsxAAAAAAaJhkIQAAAACQRLIQAAAAAKiTLAQAAAAAktjgBAAAACqnJV1pSVezuwGMkEbmu8pCAAAAACCJZCEAAAAAUCdZCAAAAAAksWYhAAAAVE6tKFMrymZ3Axghjcx3lYUAAAAAQBLJQgAAAACgTrIQAAAAAEgiWQgAAAAA1NngBAAAACqmpehKS9HV7G4AI6SR+a6yEAAAAABIIlkIAAAAANRJFgIAAAAASaxZCAAAAJVTS5laymZ3Axghjcx3lYUAAAAAQBLJQgAAAACgTrIQAAAAAEgiWQgAAAAA1NngBAAAACqmVnSlpehqdjeAEVJrYL6rLAQAAAAAkkgWAgAAAAB1koUAAAAAQBJrFgIAAEDl1FKmlrLZ3QBGSCPzXWUhAAAAAJBEshAAAAAAqJMsBAAAAACSSBYCAAAAAHU2OAEAAICKaSm60lJ0NbsbwAhpZL6rLAQAAAAAkkgWAgAAAAB1koUAAAAAQBJrFgIAAEDl1IoyNWsWQmXUinLgbYexHwAAAADAGCJZCAAAAAAkkSwEAAAAAOokCwEAAACAJDY4AQAAgMppSZmWDHzDA2Bsa2S+qywEAAAAAJJIFgIAAAAAdZKFAAAAAEASaxYCAABA5dSKrtSKrmZ3Axghjcx3lYUAAAAAQBLJQgAAAACgTrIQAAAAAEgiWQgAAAAA1NngBAAAACqmJWVaUja7G8AIaWS+qywEAAAAAJJIFgIAAAAAdZKFAAAAAEASaxYCAABA5RRFV2pFV7O7AYyQooH5rrIQAAAAAEgiWQgAAAAA1EkWAgAAAABJJAsBAAAAgDobnAAAAEDFtKQrLbHBCVRFI/NdZSEAAAAAkESyEAAAAACokywEAAAAAJJYsxAAAAAqp1aUqRVls7sBjJBG5rvKQgAAAAAgiWQhAAAAAFAnWQgAAAAAJJEsBAAAAADqbHACAAAAFdOSrrSkq9ndAEZII/NdZSEAAAAAkESyEAAAAACokywEAAAAAJJYsxAAAAAqp6XoSkthzUKoikbmu8pCAAAAACCJZCEAAAAAUCdZCAAAAAAkkSwEAAAAAOpscAIAAAAVU0tXarHBCVRFI/NdZSEAAAAAkESyEAAAAACokywEAAAAAJJYsxAAAAAqp6Uo01KUze4GMEIame8qCwEAAACAJJKFAAAAAECdZCEAAAAAkESyEAAAAACos8EJAAAAVEwtZVrS1exuACOkFhucAAAAAAANkiwEAAAAAJJIFgIAAAAAddYsBAAAgIqppSu1omh2N4ARUmtgjVKVhQAAAABAEslCAAAAAKBOshAAAAAASGLNQgAAAKiclpRpSdnsbgAjpJH5rrIQAAAAAEgiWQgAAAAA1EkWAgAAAABJJAsBAAAAgDobnAAAAEDFtBRdaSmKZncDGCEtRdeA26osBAAAAACSSBYCAAAAAHWShQAAAABAEmsWAgAAQOXU0pVarFkIVVGLNQsBAAAAgAZJFgIAAAAASSQLAQAAAIA6yUIAAAAAIIkNTgAAAKByWoqutBQ2OIGqaClscAIAAAAANGjIk4Xz58/Pq1/96kyePDnTpk3LW9/61jz88MN92pRlmSuuuCIzZszIHnvskRNPPDEPPvhgnzYbNmzIRRddlKlTp2avvfbK6aefnieeeGKouwswLMRCgG7iIYBYCIwtQ54sXLRoUd71rnflBz/4QRYuXJjnn38+c+bMybPPPtvb5uqrr84111yTa6+9Nvfdd186OjpyyimnZM2aNb1t5s6dm9tvvz0LFizI3XffnbVr1+a0007L5s2bh7rLAENOLAToJh4CiIXA2FKUZVkO5z/wm9/8JtOmTcuiRYvyhje8IWVZZsaMGZk7d24uvfTSJN1/HZk+fXo+9rGP5bzzzktnZ2f222+/fPGLX8yZZ56ZJHnqqady4IEH5o477sipp566y3939erVaW9vz4l5S1qLCcM5RGCUeb7clO/kG+ns7MyUKVOa3Z0kzYuFyZZ4eNBVH0lt0qRhGyMwunStX5/HLvvgqIqFSfOvDcVCqJ7RGA9Hw7XhV5Yelj0ntwzbGIHR5bk1m/Onr/zPAcXCYV+zsLOzM0myzz77JEmWLVuWFStWZM6cOb1t2tracsIJJ+See+5JkixevDibNm3q02bGjBmZNWtWbxuAsUQsBOgmHgKIhcDoNqy7IZdlmYsvvjive93rMmvWrCTJihUrkiTTp0/v03b69Ol59NFHe9tMnDgxe++993Ztep6/rQ0bNmTDhg29X69evXrIxgGwO0YyFibiITB6uTYEcG0IjH7DWll44YUX5qc//WluvfXW7R4rttmivSzL7c5ta2dt5s+fn/b29t7jwAMPHHzHAYbQSMbCRDwERi/XhgCuDYHRb9iShRdddFG++c1v5t///d9zwAEH9J7v6OhIku3+8rFy5crev6J0dHRk48aNWbVqVb9ttnX55Zens7Oz93j88ceHcjgAgzLSsTARD4HRybUhgGtDYGwY8mRhWZa58MIL87WvfS3/9m//lkMOOaTP44ccckg6OjqycOHC3nMbN27MokWLctxxxyVJZs+enQkTJvRps3z58jzwwAO9bbbV1taWKVOm9DkAmqVZsTARD4HRxbUhwOi8NizSlZrD4ajMUaRrwDFryNcsfNe73pUvfelL+cY3vpHJkyf3/mWkvb09e+yxR4qiyNy5czNv3rzMnDkzM2fOzLx587LnnnvmrLPO6m177rnn5j3veU/23Xff7LPPPrnkkktyxBFH5OSTTx7qLgMMObEQoJt4CCAWAmPLkCcLr7/++iTJiSee2Of8TTfdlD//8z9Pkrzvfe/LunXrcsEFF2TVqlU55phjctddd2Xy5Mm97T/5yU+mtbU1Z5xxRtatW5eTTjopN998c1pabO0OjH5iIUA38RBALATGlqIsy7LZnRgOq1evTnt7e07MW9JaTGh2d4AR9Hy5Kd/JN9LZ2em2s2yJhwdd9ZHUJk1qdneAEdK1fn0eu+yDYmGdWAjVJR721RMPv7r0ZdlrsiQjVMWzazbnT1758wHFwiGvLAQAAABGt5aiKy272GkZGD9aioGvWThsuyEDAAAAAGOLZCEAAAAAkESyEAAAAACokywEAAAAAJLY4AQAAAAqpyVlWlI2uxvACGlkvqssBAAAAACSSBYCAAAAAHWShQAAAABAEmsWAgAAQOXUijK1oqvZ3QBGSK2wZiEAAAAA0CDJQgAAAAAgiWQhAAAAAFAnWQgAAAAAJLHBCQAAAFROS7rS0uxOACOmJQPf0EhlIQAAAACQRLIQAAAAAKiTLAQAAAAAklizEAAAACqnJWVaUja7G8AIaWS+qywEAAAAAJJIFgIAAAAAdZKFAAAAAEASyUIAAAAAoM4GJwAAAFAxtaIrtaLZvQBGSq3oGnjbYewHAAAAADCGSBYCAAAAAEkkCwEAAACAOmsWAgAAQMXUUqYlZbO7AYyQWgPzXWUhAAAAAJBEshAAAAAAqJMsBAAAAACSSBYCAAAAAHU2OAEAAICKabHBCVRKI/NdZSEAAAAAkESyEAAAAACokywEAAAAAJJYsxAAAAAqp1aUqRXWLISqaGS+qywEAAAAAJJIFgIAAAAAdZKFAAAAAEASyUIAAAAAoM4GJwAAAFAxLSnTEhucQFU0Mt9VFgIAAAAASSQLAQAAAIA6yUIAAAAAIIk1CwEAAKByailTs2YhVEYj811lIQAAAACQRLIQAAAAAKiTLAQAAAAAkkgWAgAAAAB1NjgBAACAimkpyrQUNjiBqmhkvqssBAAAAACSSBYCAAAAAHWShQAAAABAEmsWAgAAQOW0pExLrFkIVdHIfFdZCAAAAAAkkSwEAAAAAOokCwEAAACAJJKFAAAAAECdDU4AAACgYoqoHoIqKRpoKzYAAAAAAEkkCwEAAACAOslCAAAAACCJNQsBAACgclqK7gOohkbmu8pCAAAAACCJZCEAAAAAUCdZCAAAAAAkkSwEAAAAAOpscAIAAAAVU4vqIaiSRua72AAAAAAAJJEsBAAAAADqJAsBAAAAgCTWLAQAAIDKaakfQDU0Mt9VFgIAAAAASSQLAQAAAIA6yUIAAAAAIIlkIQAAAABQZ4MTAAAAqJiWokhLUTS7G8AIaWS+qywEAAAAAJJIFgIAAAAAdZKFAAAAAEASaxYCAABA5dSiegiqpJH5LjYAAAAAAEkkCwEAAACAOslCAAAAACCJZCEAAAAAUGeDEwAAAKiYlhRpSdHsbgAjpJH5rrIQAAAAAEgiWQgAAAAA1EkWAgAAAABJrFkIAAAAlVMURWqFNQuhKooG5rvKQgAAAAAgiWQhAAAAAFAnWQgAAAAAJJEsBAAAAADqbHACAAAAFdOSIi2xwQlURSPzXWUhAAAAAJBEshAAAAAAqJMsBAAAAACSWLMQAAAAKqdW/w+ohkZmu8gAAAAAACSRLAQAAAAA6iQLAQAAAIAkkoUAAAAAQJ0NTgAAAKBiWooiLUXR7G4AI6SR+a6yEAAAAABIIlkIAAAAANRJFgIAAAAASaxZCAAAAJVTq/8HVEMjs11kAAAAAACSSBYCAAAAAHWShQAAAABAEslCAAAAAKDOBicAAABQMbUUqaVodjeAEdLIfFdZCAAAAAAkkSwEAAAAAOokCwEAAACAJGMgWXjdddflkEMOyaRJkzJ79ux873vfa3aXAEacWAggFgL0GIp42FLUHA5HxY6BGtXJwi9/+cuZO3duPvCBD2TJkiV5/etfn9/7vd/LY4891uyuAYwYsRBALAToIR4Cw21UJwuvueaanHvuufnLv/zLvPzlL8+nPvWpHHjggbn++uub3TWAESMWAoiFAD3EQ2C4jdpk4caNG7N48eLMmTOnz/k5c+bknnvu2a79hg0bsnr16j4HwFjXaCxMxENg/BELAbqJh8BIaG12B/rz29/+Nps3b8706dP7nJ8+fXpWrFixXfv58+fnyiuv3O7889mUlMPWTWAUej6bkiRlOfYnf6OxMOk/HnatXz8sfQRGp545Lxb2JRZC9YiHO46Hq9d2DUsfgdGpZ84PJBaO2mRhj6Io+nxdluV255Lk8ssvz8UXX9z79ZNPPpnDDz88d+eOYe8jMDqtWbMm7e3tze7GkBhoLEy2j4fLli3LK1/5yjxxxUeGtY/A6CQWdhMLAfGwW088PPhV/zWcXQRGqYHEwlGbLJw6dWpaWlq2++vIypUrt/srSpK0tbWlra2t9+sXvOAFeeihh3L44Yfn8ccfz5QpU4a9z6PJ6tWrc+CBB1Zy7InxV338ZVlmzZo1mTFjRrO7stsajYXJ9vHw4IMPTpI89thj4+YCuRFVnw/GX93xi4Vi4daqPBcS46/6+MVD8XBrVZ8Pxl/d8TcSC0dtsnDixImZPXt2Fi5cmD/8wz/sPb9w4cK85S1v2eXza7VaXvSiFyVJpkyZUrkfgh5VHnti/FUe/3i58NndWJh0x8Ok+z2p6s9DUu35kBh/VccvFm4hFnar6lzoYfzVHb94uIV42K3K8yEx/qqOf6CxcNQmC5Pk4osvzjnnnJOjjz46xx57bG644YY89thjOf/885vdNYARIxYCiIUAPcRDYLiN6mThmWeemaeffjof/vCHs3z58syaNSt33HFHb9k0QBWIhQBiIUAP8RAYbqM6WZgkF1xwQS644IJBPbetrS0f+tCH+qzPUBVVHnti/FUf/3gkFg6e8Rt/lcc/3oiFg2f8xl/l8Y9H4uHgGb/xV3n8A1WU42H/eAAAAABgt9Wa3QEAAAAAYHSQLAQAAAAAkkgWAgAAAAB1koUAAAAAQJJxnCy87rrrcsghh2TSpEmZPXt2vve97zW7S0Piu9/9bt785jdnxowZKYoiX//61/s8XpZlrrjiisyYMSN77LFHTjzxxDz44IN92mzYsCEXXXRRpk6dmr322iunn356nnjiiREcxeDMnz8/r371qzN58uRMmzYtb33rW/Pwww/3aTOex3/99dfnFa94RaZMmZIpU6bk2GOPzbe+9a3ex8fz2Bk8sXB8zgfxUDykceLh+JsPYqFYSOPEwvE3H8RCsXBYlOPQggULygkTJpR///d/Xz700EPlX/3VX5V77bVX+eijjza7a7vtjjvuKD/wgQ+Ut912W5mkvP322/s8ftVVV5WTJ08ub7vttvL+++8vzzzzzHL//fcvV69e3dvm/PPPL1/0oheVCxcuLH/yk5+Ub3zjG8sjjzyyfP7550d4NI059dRTy5tuuql84IEHyqVLl5Z/8Ad/UB500EHl2rVre9uM5/F/85vfLP/5n/+5fPjhh8uHH364fP/7319OmDChfOCBB8qyHN9jZ3DEwvE7H8RD8ZDGiIfjcz6IhWIhjRELx+d8EAvFwuEwLpOFr3nNa8rzzz+/z7nDDjusvOyyy5rUo+GxbRDs6uoqOzo6yquuuqr33Pr168v29vbyc5/7XFmWZfnMM8+UEyZMKBcsWNDb5sknnyxrtVp55513jljfh8LKlSvLJOWiRYvKsqze+MuyLPfee+/yf//v/13JsbNrYmF15oN4KB6yc+JhNeaDWCgWsnNiYTXmg1goFg6FcXcb8saNG7N48eLMmTOnz/k5c+bknnvuaVKvRsayZcuyYsWKPmNva2vLCSec0Dv2xYsXZ9OmTX3azJgxI7NmzRpz709nZ2eSZJ999klSrfFv3rw5CxYsyLPPPptjjz22UmNnYMTCas0H8VA8pH/iYXXmg1goFtI/sbA680EsFAuHwrhLFv72t7/N5s2bM3369D7np0+fnhUrVjSpVyOjZ3w7G/uKFSsyceLE7L333v22GQvKsszFF1+c173udZk1a1aSaoz//vvvzwte8IK0tbXl/PPPz+23357DDz+8EmOnMWJhdeaDeCgesnPiYTXmg1goFrJzYmE15oNYKBYOldZmd2C4FEXR5+uyLLc7N14NZuxj7f258MIL89Of/jR33333do+N5/EfeuihWbp0aZ555pncdtttecc73pFFixb1Pj6ex87giIVbjNf5IB6KhwyMeLjFeJwPYqFYyMCIhVuMx/kgFoqFQ2XcJQunTp2alpaW7TLAK1eu7JNNXr9+fTZu3DjS3Rtyzz33XFavXp0k2WuvvZIkv/zlL3s/T5Inn3wye++9d1avXp3Jkydn48aNefTRR/tkzpcvX57Zs2f3vtZo9t73vjf/9//+33zrW9/KlClTKjf+adOmZdq0abn88stz77335uMf/3jmzp2bZPyPfVsTJ07MpEmTmt2NUUksrMZ8EA/Fwx7iYf/Ew/E/H8RCsbCHWNg/sXD8zwexUCzc2u7Gw6Isy3II+zMqHHPMMZk9e3auu+663nOHH3543vKWt2T+/PlZv3592vfYOxuzvom9BHZXR0dHli1b5qKwH2IhVId4uHPiIVSDWLhzYiFUx+7Gw3FXWZgkF198cc4555wcffTROfbYY3PDDTfksccey/nnn5+ke3HXjVmf1xWnZUJtYlLUUrTUklqRFEWKWi0piqTnY1GkKLY83vex1L+un29p6f6Yos/5sii6V4gstrxG33Pp/rr+emXP8+qPl8U27Wpb2vc8Vqb7n02tp31Sbv15sfVrpf7aW84nW3+9g49bvVa/7bZ9jeygfa3v1zt6raSfz7d+Ti1JUe7wdbqPMr2Z8N7Hym1ep+x9nZ5vWZ826Tnf/UrFNu2Kotzyea2sf/vKLT8mvW3KPl/XsnXbcsv5okxL0ZVa/fOeo/vb2vdcLT2fd9Xb9JzrSktRprXYnFqRtKRMUW/TUpSppf55z8eiK0V6/t2utCT1j92v1VpsTku6uh/vaZeu1JL6a3TVH+9+7S2Pd7dtqfdrYrE5Rf3fban3vSVl/fN09z1JS5G0pKh/XqSWLUf317VMKFqydk2Zg2f/VzZu3OiCsB+DioW13h/e/uNhT5zbWSzcOmb2Fw/rz91pLOwTF3vO9W1btgwgFm4b//qLh7UBxMI+sXdXbdN/LOwndvX3Gkk/sXDreNbva5bbjGv7WLjtYzuKh0Vt17GwJ2buLB621HYdC2vZ6vN+4mFrrWuXsbDn853FwwnF5l3Gwq1jZ3/xcELx/C5jYc/nu4qHEwYQC7s/r4mHA9BoPOyNc7t7bVhr2XUsHMC1YVmr7TIWjvS1Yfe16gDi17avkb7te/7N3b023PI6u3ltWCt3GQsHdG1Yj3M7i4UDuTZsLbp2GQsHcm04oR7bxuO1YWtaxMIBGo5rw2LrOLc714Ytu46FI31t2PP6ydZtGr827N1WdiexcEDXhvXX2e1rw9rAfk/e5bVhPU4lO4mFA7g2rNXKXcbCgVwbtta6BvR78li8NmxNOaDfk4fy2nBcJgvPPPPMPP300/nwhz+c5cuXZ9asWbnjjjty8MEH92nXmglpLSZ0B8Fiyy+/RdET4HoC4c4uCPsGutRadnh+1xeEfb8eVLKwd/IP7oJwl7/sDiRZuO257ODcQJKFRf9BsO/r7OKCsP78nQXBAScL+3xMnyBYFNsmC/teAO7wgnDrr7f6uKNk4ZaA1/8F4dZfb0kWlluC2U4uCGt9gmDZ5+OEothhAGwpyu5gVRRpSVF/7e5A1R3Eyvrn3f/mxJ7xJTsMgC07uSBsybbJwlpqW/4XRz8GFwu3in39xcPajpKF28TCbWPmjuLhDpOF28TCwVwQ7igW7jDO7ThZuMtYuIsLwh3FuUEnC2sDiIXbxrwdvWZtCJOFu4iF218Qbh8P+yQL+4mF/V0Qbh0Pd5Qs3DYW1raOgf3EwwlFbZexsGUnF4Q98XBC7+v2Hwtbiq1jYP/xcEfJwm1jYc8FoXi4a43Gw+7rwl3EwoFcG277S/Qgrw0HnSzcUTwcomvDAScLtz2XvucGnCwsdh4PB5wsrD+/33jYSLJwZ/Fwu2Th4K4N+0sWNnptuCVZOP6uDVtTEwsHaDiuDftPFvaNebu8Nuw3Wdi8a8O+r7+bycJdxMIBXRtu+3t0BnltWBvY78m7vDasx6ndvTbcLlm4g1g4kGvD/pKF4+HacEI/ycLhvDYcl8nCJLngggtywQUXNLsbAE0lFgJ0Ew8BxEJgYGq7bgIAAAAAVIFkIQAAAACQRLIQAAAAAKiTLAQAAAAAkkgWAgAAAAB1koUAAAAAQBLJQgAAAACgrrXZHWim57MpRVkkqaUoNydlkaRIUdaSFEnPx64iRVH/vKgfqdU/T/3rnsda+n5ddn8siyIps+X5251L99dd3V+XtXqbWv3rYpt2tfrHIr2Plen+Z1PraZ+UW39ebP1aW4bX+9rZ+usdfNzqtfptt+1rZAfta32/3tFrJf18vvVzakmKcoev032U3f9+nz6V27xO2fs6Pd+yPm3Sc777lYpt2hVFueXzWln/cSh7v81lb5uyz9dlutt21b+ubfWxLLpSq3/ec3R/W/ueq6Xn8656m55zXWkpynQVm1MrkpaUKeptWooytdQ/7/lYdKVI98da0ZWWpP6x+7WeLzanJV3dj/e0S1dqSf01uuqPd7/2lse727bU+7Wp6EpR/3db6n1vSVn/PN19T9JSJC0p6p8XqWXL0f11MqEosnZN73eX3dQ3FvadTDuOhy0DiIVbxcz+4mFZ23Us3Oq1+4uHZTmAWLht/OsvHtYGEAv7xN5dtU3/sbCf2NXfayT9xMKt41m/r1luM67tY+G2j+0oHha1XcfCnpi5s3iY2q5jYS1bfd5PPCxrXbuMhT2f7ywedhWbdxkLt46d/cXDCcXmXcbCns93FQ8nDCAWdn8e8XAI9cTDouzKrmPhAK4N05KtJsuOY+EArg3LWm2XsXCkrw27r1UHEL+2fY30bd/zb+7qtZKdx8Mtr7Ob14a1MruKhT1xLtlJPKzHuZ3FwoFcG6Z+/bS714ZlPbaNx2vDVrFwyDVybVh0bRXndufacKvXGC3Xhj2vn/R5Cxq+Nix7SrR2Egv7e36fx+uvs9vXhrWB/Z687WPbxcOeOJWdxMIBXBuWtXJAvyfv6tqwrHUN6PfksXht2JpyQL8nD+W1YSWThWVZ5gUveEHuXvt/k831k883tUvAILzgBS9IWbowHKwdxkJgTBIPd8928VBMhDFJLNw9EydOTEdHR+5e4doQxrqOjo5MnDhx0M+vZLKwKIqsXbs2jz/+eKZMmdLs7jRs9erVOfDAA8ds/5OxP4ax3v9k7I+hp/9Fz1/7aNhYj4WDMdZ/7gejamOu2ngT8XAoVC0eVnmeGPP4JRbuvkmTJmXZsmXZuHFjs7sypKo2F0aK93X4DMV7O3HixEyaNGnQfahksrDHlClTxvQP9VjvfzL2xzDW+5+MjzGwe6r4M2DM41/VxsvQqNrPTdXGmxgz7MqkSZN2K8EwmpkLw8P7Onya+d7a4AQAAAAASCJZCAAAAADUVTJZ2NbWlg996ENpa2trdlcGZaz3Pxn7Yxjr/U/G/hjGev9Hgyq+h8Y8/lVtvEk1xzzUqvYeVm28iTFXQdXGy8D52Rge3tfhMxre26K0XRQAAAAAkIpWFgIAAAAA25MsBAAAAACSSBYCAAAAAHWShQAAAABAknGYLPzud7+bN7/5zZkxY0aKosjXv/71XT5n0aJFmT17diZNmpSXvOQl+dznPjf8Hd2JRsfwta99Laecckr222+/TJkyJccee2y+/e1vj0xnd2Aw34Me3//+99Pa2ppXvvKVw9a/gRjMGDZs2JAPfOADOfjgg9PW1paXvvSl+T//5/8Mf2d3YDD9v+WWW3LkkUdmzz33zP7775+/+Iu/yNNPPz38nd2B+fPn59WvfnUmT56cadOm5a1vfWsefvjhXT5vtM3lZlu1alXOOeectLe3p729Peecc06eeeaZnT7na1/7Wk499dRMnTo1RVFk6dKl27XZsGFDLrrookydOjV77bVXTj/99DzxxBPDM4gGDWbMZVnmiiuuyIwZM7LHHnvkxBNPzIMPPtinzYknnpiiKPocb3vb24ZxJP277rrrcsghh2TSpEmZPXt2vve97+20/UDmxW233ZbDDz88bW1tOfzww3P77bcPV/cHZajHfPPNN2/3/SyKIuvXrx/OYQxYI+Ndvnx5zjrrrBx66KGp1WqZO3fuDtuN9u/xcBMPx188FAvHfyxMxEMG5sknn8x//+//Pfvuu2/23HPPvPKVr8zixYt3+hy/N+xao+/raMtLjGaD+ZntMVI5k3GXLHz22Wdz5JFH5tprrx1Q+2XLluX3f//38/rXvz5LlizJ+9///rz73e/ObbfdNsw97V+jY/jud7+bU045JXfccUcWL16cN77xjXnzm9+cJUuWDHNPd6zR/vfo7OzM29/+9px00knD1LOBG8wYzjjjjPzrv/5rbrzxxjz88MO59dZbc9hhhw1jL/vXaP/vvvvuvP3tb8+5556bBx98MF/5yldy33335S//8i+Huac7tmjRorzrXe/KD37wgyxcuDDPP/985syZk2effbbf54zGudxsZ511VpYuXZo777wzd955Z5YuXZpzzjlnp8959tlnc/zxx+eqq67qt83cuXNz++23Z8GCBbn77ruzdu3anHbaadm8efNQD6Fhgxnz1VdfnWuuuSbXXntt7rvvvnR0dOSUU07JmjVr+rR75zvfmeXLl/cen//854dzKDv05S9/OXPnzs0HPvCBLFmyJK9//evze7/3e3nsscd22H4g8+Lee+/NmWeemXPOOSf/8R//kXPOOSdnnHFGfvjDH47UsHZqOMacJFOmTOnz/Vy+fHkmTZo0EkPaqUbHu2HDhuy33375wAc+kCOPPHKHbUb793gkiIfjKx6KheM/FibiIQOzatWqHH/88ZkwYUK+9a1v5aGHHsonPvGJvPCFL+z3OX5v2LXBvK+jLS8xWg3mve0xojmTchxLUt5+++07bfO+972vPOyww/qcO++888rXvva1w9izgRvIGHbk8MMPL6+88sqh71CDGun/mWeeWX7wgx8sP/ShD5VHHnnksParEQMZw7e+9a2yvb29fPrpp0emUw0YSP8//vGPly95yUv6nPv0pz9dHnDAAcPYs4FbuXJlmaRctGhRv21G+1weaQ899FCZpPzBD37Qe+7ee+8tk5T/+Z//ucvnL1u2rExSLlmypM/5Z555ppwwYUK5YMGC3nNPPvlkWavVyjvvvHPI+j8YgxlzV1dX2dHRUV511VW959avX1+2t7eXn/vc53rPnXDCCeVf/dVfDVvfB+o1r3lNef755/c5d9hhh5WXXXbZDtsPZF6cccYZ5Zve9KY+bU499dTybW972xD1evcMx5hvuummsr29fcj7OhQaHe/W+vs5He3f4+EmHnYbT/FQLOw2nmNhWYqHDMyll15avu51r2voOX5v2LXBvK87MlryEqPJ7ry3I5kzGXeVhY269957M2fOnD7nTj311Pz4xz/Opk2bmtSr3dPV1ZU1a9Zkn332aXZXBuymm27Kr371q3zoQx9qdlcG5Zvf/GaOPvroXH311XnRi16Ul73sZbnkkkuybt26ZndtQI477rg88cQTueOOO1KWZX7961/nq1/9av7gD/6g2V1L0v0XlCQ7/Zkej3N5d9x7771pb2/PMccc03vuta99bdrb23PPPfcM+nUXL16cTZs29XmvZ8yYkVmzZu3W6w6FwYx52bJlWbFiRZ/xtLW15YQTTtjuObfcckumTp2a3/md38kll1yyXaXNcNu4cWMWL1683c/5nDlz+h3fQOZFf22a/f1Mhm/MSbJ27docfPDBOeCAA3LaaaeNir96D2a8AzGav8cjQTzsNl7ioVi4xXiNhYl4yMD1/B72p3/6p5k2bVqOOuqo/P3f//1On+P3hl0bzPu6rbGYlxgJg31vRzpnUvlk4YoVKzJ9+vQ+56ZPn57nn38+v/3tb5vUq93ziU98Is8++2zOOOOMZndlQH7xi1/ksssuyy233JLW1tZmd2dQHnnkkdx999154IEHcvvtt+dTn/pUvvrVr+Zd73pXs7s2IMcdd1xuueWWnHnmmZk4cWI6Ojrywhe+MJ/5zGea3bWUZZmLL744r3vd6zJr1qx+243Hubw7VqxYkWnTpm13ftq0aVmxYsVuve7EiROz99579zk/ffr03XrdoTCYMfec39HPztbPOfvss3PrrbfmO9/5Tv7mb/4mt912W/7oj/5oCHu/a7/97W+zefPmXfZ1awOZF/21afb3Mxm+MR922GG5+eab881vfjO33nprJk2alOOPPz6/+MUvhmcgAzSY8Q7EaP4ejwTxcIvxEA/Fwi3GayxMxEMG7pFHHsn111+fmTNn5tvf/nbOP//8vPvd784//MM/9Pscvzfs2mDe122NtbzESBnMe9uMnMnYzMwMsaIo+nxdluUOz48Ft956a6644op84xvf2OFF4mizefPmnHXWWbnyyivzspe9rNndGbSurq4URZFbbrkl7e3tSZJrrrkmf/Inf5LPfvaz2WOPPZrcw5176KGH8u53vzt/+7d/m1NPPTXLly/Pe9/73px//vm58cYbm9q3Cy+8MD/96U9z991377LteJrL/bniiity5ZVX7rTNfffdl2TH4y7Lcljej+F63WRkxryjn52tz73zne/s/XzWrFmZOXNmjj766PzkJz/Jq171ql2OYSjtqq8Dab/t+UZfc6QN9Zhf+9rX5rWvfW3v48cff3xe9apX5TOf+Uw+/elPD1W3B204vh+j/Xs8GOLhjlUlHoqF4z8WJuIhu9bV1ZWjjz468+bNS5IcddRRefDBB3P99dfn7W9/e7/Pq8LvDbtjsO9rj7GWlxhJjb63zcqZVD5Z2NHRsd1fklauXJnW1tbsu+++TerV4Hz5y1/Oueeem6985Ss5+eSTm92dAVmzZk1+/OMfZ8mSJbnwwguTdE+esizT2tqau+66K7/7u7/b5F7u2v77758XvehFvYnCJHn5y1+esizzxBNPZObMmU3s3a7Nnz8/xx9/fN773vcmSV7xildkr732yutf//p85CMfyf7779+Ufl100UX55je/me9+97s54IADdtp2PM3lnbnwwgt3uevki1/84vz0pz/Nr3/96+0e+81vfrPdX1Ib0dHRkY0bN2bVqlV9qmlWrlyZ4447btCvuzPDOeaOjo4k3X9h3vrnfOXKlTt9n171qldlwoQJ+cUvfjFivxxPnTo1LS0tO/w539n4djUv+muzOz8nQ2W4xrytWq2WV7/61U2vphnMeAdiNH+Pd4d4uGPjPR6KhVuM11iYiIcM3P7775/DDz+8z7mXv/zlO92spCq/N+yOwbyvPcZiXmIkNfreNitnUvnbkI899tgsXLiwz7m77rorRx99dCZMmNCkXjXu1ltvzZ//+Z/nS1/60qhZZ24gpkyZkvvvvz9Lly7tPc4///wceuihWbp0aZ+1dkaz448/Pk899VTWrl3be+7nP/95arXaLpNco8Fzzz2XWq1vOGhpaUmy5a9sI6ksy1x44YX52te+ln/7t3/LIYccssvnjJe5vCtTp07NYYcdttNj0qRJOfbYY9PZ2Zkf/ehHvc/94Q9/mM7Ozt36JXb27NmZMGFCn/d6+fLleeCBB4btl+PhHPMhhxySjo6OPuPZuHFjFi1atNPxPPjgg9m0adOIJtInTpyY2bNnb/dzvnDhwn77OpB50V+b4fp+NmK4xrytsiyzdOnSpv1hpMdgxjsQo/l7vDvEw2rGQ7Fwi/EaCxPxkIE7/vjj8/DDD/c59/Of/zwHH3xwv8+pyu8Nu2Mw72sydvMSI6nR97ZpOZNh3T6lCdasWVMuWbKkXLJkSZmkvOaaa8olS5aUjz76aFmWZXnZZZeV55xzTm/7Rx55pNxzzz3Lv/7rvy4feuih8sYbbywnTJhQfvWrX23WEBoew5e+9KWytbW1/OxnP1suX76893jmmWfGRP+3NRp2Q250DGvWrCkPOOCA8k/+5E/KBx98sFy0aFE5c+bM8i//8i/HRP9vuummsrW1tbzuuuvKX/3qV+Xdd99dHn300eVrXvOapvT/f/2v/1W2t7eX3/nOd/r8TD/33HO9bcbCXG62N73pTeUrXvGK8t577y3vvffe8ogjjihPO+20Pm0OPfTQ8mtf+1rv108//XS5ZMmS8p//+Z/LJOWCBQvKJUuWlMuXL+9tc/7555cHHHBA+S//8i/lT37yk/J3f/d3yyOPPLJ8/vnnR2xs/RnMmK+66qqyvb29/NrXvlbef//95Z/92Z+V+++/f7l69eqyLMvyl7/8ZXnllVeW9913X7ls2bLyn//5n8vDDjusPOqoo0Z8zAsWLCgnTJhQ3njjjeVDDz1Uzp07t9xrr73K//qv/yrLcnDz4vvf/37Z0tJSXnXVVeXPfvaz8qqrripbW1v77KLaTMMx5iuuuKK88847y1/96lflkiVLyr/4i78oW1tbyx/+8IcjPr5tNTresix74/3s2bPLs846q1yyZEn54IMP9j4+2r/HI0E8HF/xUCwc/7GwLMVDBuZHP/pR2draWn70ox8tf/GLX5S33HJLueeee5b/+I//2NvG7w2NG8z7OtryEqPVYN7bbY1EzmTcJQv//d//vUyy3fGOd7yjLMuyfMc73lGecMIJfZ7zne98pzzqqKPKiRMnli9+8YvL66+/fuQ7vpVGx3DCCSfstP1o7/+2RkOycDBj+NnPflaefPLJ5R577FEecMAB5cUXX9wnuTWSBtP/T3/60+Xhhx9e7rHHHuX+++9fnn322eUTTzwx8p0vyx32PUl500039bYZC3O52Z5++uny7LPPLidPnlxOnjy5PPvss8tVq1b1abPt+3rTTTft8L3/0Ic+1Ntm3bp15YUXXljus88+5R577FGedtpp5WOPPTYyg9qFwYy5q6ur/NCHPlR2dHSUbW1t5Rve8Iby/vvv7338scceK9/whjeU++yzTzlx4sTypS99afnud7+7fPrpp0doVH199rOfLQ8++OBy4sSJ5ate9apy0aJFvY8Ndl585StfKQ899NBywoQJ5WGHHVbedtttwz2Mhgz1mOfOnVsedNBB5cSJE8v99tuvnDNnTnnPPfeMxFAGpNHx7mjOHnzwwX3ajPbv8XATD8dfPBQLx38sLEvxkIH5p3/6p3LWrFllW1tbedhhh5U33HBDn8f93jA4jb6voy0vMZoN5md2ayORMynKsgn3GAIAAAAAo07l1ywEAAAAALpJFgIAAAAASSQLAQAAAIA6yUIAAAAAIIlkIQAAAABQJ1kIAAAAACSRLAQAAAAA6iQLAQAAAIAkkoUAAAAAQJ1kIQAAAACQRLIQAAAAAKiTLAQAAAAAkiT/P+N7o/FZC40pAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "nx = 400\n", "arguments = InitialConditions.genRayleighTaylor(nx, nx*3, 1.4, version=0)\n", "plotVars(arguments['rho'], arguments['rho_u'], arguments['rho_v'], arguments['E'])" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Initialized /home/smyalygames/PycharmProjects/FiniteVolumeGPU/data/euler_rayleigh-taylor.nc\n", "Opening /home/smyalygames/PycharmProjects/FiniteVolumeGPU/data/euler_rayleigh-taylor.nc\n", "Keyword arguments: {'mode': 'w', 'clobber': False}\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Constructed in 0.0018711090087890625 seconds\n", "Simulating to t=30.000000. Estimated 15160 timesteps (dt=0.001979)\n", "0% [#######################=======] 100%. Total: 13s, elapsed: 10s, remaining: 3s\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "Closing /home/smyalygames/PycharmProjects/FiniteVolumeGPU/data/euler_rayleigh-taylor.nc\n" ] } ], "source": [ "nx = 151\n", "ny = nx*3\n", "g = 0.1\n", "gamma = 1.4\n", "t_end = 30\n", "n_save = 300\n", "outfile = \"data/euler_rayleigh-taylor.nc\"\n", "outdata = None\n", "\n", "arguments = InitialConditions.genRayleighTaylor(nx, ny, gamma)\n", "arguments['context'] = my_context\n", "outfile = runSimulation(outfile, t_end, arguments)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Initialized /home/smyalygames/PycharmProjects/FiniteVolumeGPU/data/euler_rayleigh-taylor.nc\n", "Opening /home/smyalygames/PycharmProjects/FiniteVolumeGPU/data/euler_rayleigh-taylor.nc\n", "Arguments: ('r',)\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "/home/smyalygames/PycharmProjects/FiniteVolumeGPU/data/euler_rayleigh-taylor.nc created Wed Feb 12 16:40:15 2025 contains 301 timesteps\n", "Simulator arguments: \n", " {'rho': '[...]', 'rho_u': '[...]', 'rho_v': '[...]', 'E': '[...]', 'nx': 151, 'ny': 453, 'dx': 0.0033112582781456954, 'dy': 0.0033112582781456954, 'g': 0.1, 'gamma': 1.4, 'boundary_conditions': '[north=Type.Reflective, south=Type.Reflective, east=Type.Reflective, west=Type.Reflective]', 'context': 'CudaContext id 96726128804080'}\n", "0% [##############================] 100%. Total: 10s, elapsed: 5s, remaining: 5s\n", "0% [#############################=] 100%. Total: 10s, elapsed: 10s, remaining: 1s\n" ] }, { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "
\n", " \n", "
\n", " \n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "
\n", "
\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Closing /home/smyalygames/PycharmProjects/FiniteVolumeGPU/data/euler_rayleigh-taylor.nc\n" ] } ], "source": [ "#outfile = 'data/euler_rayleigh-taylor_0007.nc'\n", "createAnimation(outfile, vis_type=VisType.Density, vmin=1, vmax=2, cmap=plt.cm.coolwarm, save_anim=False, fig_args={'figsize':(3.4, 8)})" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "ShallowWaterGPU", "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.19" } }, "nbformat": 4, "nbformat_minor": 2 }