{ "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": 3, "metadata": {}, "outputs": [], "source": [ "%setup_logging --out test_schemes.log logger" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "%cuda_context_handler --no_autotuning my_context " ] }, { "cell_type": "code", "execution_count": 5, "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": 12, "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": 7, "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": 13, "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": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzkAAADhCAYAAADrl3vGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAmq0lEQVR4nO3de5SV1Zmg8eetorgjRLmooHIRpVUiGtRoSAwSBdMqyUhIcMaIRmkVjWPH5SWumU7S7ZKZhGiMRqQVQ1pbcQVsbRbxQjQo2l5KRQUUg4hQBgXEQQW5VNWeP86hUpK6nLqeqlPPb62z6nyXs7+3dnZCvdl7vydSSkiSJElSoSjKdwCSJEmS1JxMciRJkiQVFJMcSZIkSQXFJEeSJElSQTHJkSRJklRQTHIkSZIkFRSTHEmSJEl5ExFzImJjRCyv5XpExC0RsToiXouIY+tr0yRHkiRJUj79FphQx/XTgeHZ1zTg9voaNMmRJEmSlDcppaeALXXcMhH4Xcp4DugTEQfU1Wan5gxQkiRJUscwYcKEtHnz5nrve+mll1YAO6qdmp1Smt2ARw0E1lc7Lsue21DbB0xyJEmSJDXY5s2befHFF+u9r6ioaEdKaXQTHhU1nEt1fcAkR5IkSVKjpFRnrtFcyoCDqh0PAv5S1weatCcnIiZExKpspYNrm9KWJEmSpPYjpURlZWW9r2bwMPD9bJW1LwNbU0q1LlWDJszkREQxcBtwKpns6sWIeDiltLKxbUqSJElqP5pjJici7gO+DvSNiDLgn4CSbPuzgEXAN4HVwHbg/PrabMpyteOB1SmlNdng7idT+cAkR5IkSeoAmiPJSSlNqed6AqY3pM2mJDk1VTk4Ye+bImIamXrW9OjR40sjRoxo9AM/+ugj1qxZU3V88MEH069fv0a3J0mSJHU0a9euZfPmzTVt5m+wZlqO1uyakuTkVOUgWx5uNsDo0aNTaWlpox84b948zjnnnKrOvOaaa7j00ksb3Z4kSZLU0Ywe3ZRCZ3+VUmqtwgMN1pQkp8FVDiRJkiQVjkJMcl4EhkfEEOA94HvAOc0SlSRJkqQ2r+CWq6WUyiPiMuBRoBiYk1Ja0WyRSZIkSWrTCnEmh5TSIjIl3SRJkiR1IIW6J0eSJElSB2aSI0mSJKmgFNyeHEmSJEkdl8vVJEmSJBUckxxJkiRJBcXlapIkSZIKijM5kiRJkgqGe3IkSZIkFRyXq0mSJEkqKM7kSJIkSSooJjmSJEmSCkZKyeVqkiRJkgqLMzmSJEmSCopJjiRJkqSCYpIjSZIkqWC4J0eSJElSwXEmR5IkSVJBMcmRJEmSVDBcriZJkiSp4DiTI0mSJKmgmORIkiRJKiguV5MkSZJUMFJKzuRIkiRJKiwmOZIkSZIKSltdrlaU7wAkSZIktU97lqzV9apPREyIiFURsToirq3heu+I+M+IeDUiVkTE+fW1aZIjSZIkqcFySXDqS3Iiohi4DTgdOAKYEhFH7HXbdGBlSulo4OvAzIjoXFe7JjmSJEmSGqUZZnKOB1anlNaklHYB9wMT934M0CsiAugJbAHK62rUPTmSJEmSGiXHPTl9I6K02vHslNLs7PuBwPpq18qAE/b6/K3Aw8BfgF7Ad1NKdT643iQnIg4CfgfsD1Rmg/pVROwLzAMGA2uBySmlj+prT5IkSVJhyLG62uaU0uharkVNze51PB5YBpwCDAMej4inU0of1/bAXJarlQM/Sin9HfBlYHp2ndy1wB9TSsOBP2aPJUmSJHUAzbEnh8zMzUHVjgeRmbGp7nxgQcpYDbwDjKir0XqTnJTShpTSy9n3nwBvkJlWmgjMzd42F/hWfW1JkiRJKhyVlZX1vurxIjA8IoZkiwl8j8zStOrWAeMAImIAcDiwpq5GG7QnJyIGA8cAzwMDUkobIJMIRUT/hrQlSZIkqX1r6peBppTKI+Iy4FGgGJiTUloRERdnr88C/hn4bUS8TmZ52zUppc11tZtzkhMRPYH5wP9MKX2cKW6Q0+emAdMADj744FwfJ0mSJKmNa2qSk21jEbBor3Ozqr3/C3BaQ9rMqYR0RJSQSXDuTSktyJ7+ICIOyF4/ANhYS9CzU0qjU0qj+/Xr15DYJEmSJLVRKaXmWK7WIupNcrL1qO8C3kgp/bLapYeB87LvzwMeav7wJEmSJLVVzVB4oEXkslztK8C5wOsRsSx77sfADOCBiPgBmc1A32mRCCVJkiS1SflKYupTb5KTUlpKzfWrIVvlQJIkSVLHsme5WlvUoOpqajtKS0t59dVXm9xOjx49+Pu//3t69erVDFFJkiSpI2m3MzlqmxYsWMCNN97Y5HYOOuggjjvuOJMcSZIkNZhJjhrk008/5e6772bz5ppLgD/99NPN8pytW7dy8803s++++9Z4ffz48Zx00knN8ixJkiQVFpMc5WTPQPn444+57bbbWLVqVa335vpdRXX55JNPuPXWW2u93qtXL0488cRmfaYkSZLaP/fkKCd33nknzzzzDACfffYZ77//fq33nn322ZxxxhlNfuaWLVuYMWMGmzZtqvH6vHnzWLFiBQDdunXjRz/6EcOGDWvycyVJktT+OZOjv1FRUcHu3buBzAB54oknuO+++6qud+nSha5du9b42RNOOIGpU6c2OYaysjLuuusuPvnkkxqvv/LKK5SWlgKZWZ3vfve7DBw4EICioiI6d+7c5BgkSZLUPpnk6G889thj/PKXv6waHHtmTAD23XdfZs6cyUEHHVTjZ5trNqVfv37867/+K5999lmN12fNmsXvf/97ALZv385VV11F7969ARgxYgQzZsygZ8+ezRKLJEmS2heXq6nKrl272Lp1K2+++SaLFy+uOr/PPvvQv39/APr378+YMWM49NBDWzSWLl261FlY4Omnn+app54CMpn6smXLKC8vB+DDDz/k/fffp1+/flWJjyRJkjqGlJIzOfqr0tJSpk+fzsaNG6vORQTXXnstZ555JgCdO3fm4IMPzleIVS655BImTZoEZGZypk2bVvX9PKtWreLMM89k/PjxzJw5k+Li4nyGKkmSpFZmktOBVVRUsHbtWnbu3AnAypUrWblyJd27d2fEiBEUFRURERx55JEcddRReY728wYMGMCAAQOATFnrI444omof0fbt23nrrbcYNGgQK1eurEpy9t9//1pLUkuSJKlwuFytA9uyZQtTp07lz3/+MwA7d+5k165dnH322dx0001V97X1JV89evTg1ltvrUpyXnnlFSZPnszTTz/NqaeeWnXfz372M6ZNm5avMCVJktRKnMnpgCorK3n11Vd5++23WbduHR988AGQKSowduxYRo8eXTVL0h5ExOdmaIYMGcLYsWNZu3Ytr7/+etUgX7ZsGU888QQjR46kX79++QpXkiRJLagt78kpyncAhWzXrl1cf/31nHvuuaxfv77q/DHHHMOCBQu4/PLL8xhd0x122GHMmzeP6667jqKivw6lO++8k4kTJ/Jf//VfeYxOkiRJLa2ysrLeVz44k9NCli5dyptvvsm7774LwFlnnUWfPn0AOOqoo+jZsyedOrXv7o8IunbtyvDhw/n+979fNYiff/553nrrLRYvXsy2bds49dRT6du3b56jlSRJUnNrqzM57fuv7DZszpw5/Pa3vyWlxIABA/iXf/kXjjzySCCTHBSSL33pS9x1111Vx1deeSVvvvkmv/71r/m3f/s3Hn30UZMcSZKkAmSS00E89dRTLFmyhFdeeYXi4mK++93vMnLkSAYMGFBwyU111X+38ePH0717dxYsWMDatWuZM2cOL7zwAueee26bL64gSZKk3LTlPTkmOc3siSee4Kc//SkA3bp14/zzz2fcuHF5jqp1nX766Zx22mm8+eabrFq1ijvuuIPBgwdz5plnmuRIkiQVEEtIF7ilS5cyd+5cSktLAZg6dSonn3wyf/d3f5fnyPKjqKiI6dOnM3bsWH75y1+yefNmrrnmGo4++miuvPJKunbtmu8QJUmS1ERtdSbH6mpNlFJi165dvPHGG9x55528/vrrdOnShbFjxzJ16lQOPPDAfIeYFxHBuHHjOO+88xg4cCC7du1i3rx5PPTQQ3z66aeUl5fnO0RJkiQ10Z4la3W98sEkp4leffVVvvOd73DLLbcAMGXKFB566CFOOeWUPEfWNnTv3p1f/OIXzJ49m/79+/Pmm28yZcoUZsyY0WanNyVJklS/lJIlpAvVhx9+yGOPPcaOHTsAOPTQQxk/fnyeo2o7OnXqxJe//GUGDBhAt27d2LhxI4sXL6ZPnz5tdnpTkiRJuWmrf8+Z5EiSJElqFJOcArNr1y7WrVvH+vXrSSnRu3dv9t9/f/bbb798h9YmlZSUMHToUFJKlJWV8cknn7Bq1Sr69etHv3798h2eJEmSGmjPcrW2yD05jfTuu+8yadIkrrrqKnbu3MmECRNYvHgx5513Xr5Da5MOOOAA7r//fm699Va6d+/OkiVL+MY3vsHtt9+e79AkSZLUSBYeKBC7d++mtLSUZ599lnXr1lFRUcGYMWM45phjGDRoEL169cp3iG1ScXEx/fv3Z9iwYXz1q19l6NChbNiwgRUrVrBkyRLKysryHaIkSZIayCSnQGzdupVLLrmEiy++mI8++ohRo0bx0EMPceWVV+Y7tHZhxIgRzJ8/n+uuu46ioiIWLFjAN7/5TRYuXJjv0CRJktQAbbm6mklOA6WU2LFjR1U1teLiYrp3707nzp3zHFn7UFRURLdu3ejatSsRQXl5Odu3b2f37t35Dk2SJEkN1BwzORExISJWRcTqiLi2lnu+HhHLImJFRCypr00LD0iSJElqlKYuR4uIYuA24FSgDHgxIh5OKa2sdk8f4DfAhJTSuojoX1+7OSc52QBKgfdSSmdExL7APGAwsBaYnFL6KOffqJ2pqKjg97//PcuXL2fTpk3st99+nHPOORx99NF06mSu2FCHH344V199Nc899xxPPvkkjz/+ODt27OBb3/oWw4cPz3d4kiRJykEz7Lk5HlidUloDEBH3AxOBldXuOQdYkFJal33mxvoabchf51cAbwD7ZI+vBf6YUpqRnVa6FrimAe21K+Xl5cyZM4fHHnsMgCOOOILrr7+eAQMG5Dmy9mnkyJGMHDmSmTNn8uSTT/Kf//mf/OEPf+Cwww4zyZEkSWoHmqmE9EBgfbXjMuCEve45DCiJiD8BvYBfpZR+V1ejOe3JiYhBwN8Dd1Y7PRGYm30/F/hWLm1JkiRJKgw57snpGxGl1V7TqjURNTW713En4Etk8pHxwP+KiMPqiivXmZybgavJZE57DEgpbcj+chtyWRvXXlVUVLB79+6qTLVTp0506tSJiJr+M1FDFBcXU1JSQkVFBfDXvrZ/JUmS2r4cl6ttTimNruVaGXBQteNBwF9quGdzSmkbsC0ingKOBt6q7YH1zuRExBnAxpTSS/XdW8vnp+3J2jZt2tSYJvLuN7/5DZMnT2bZsmX07t2bm266iVtuuYU+ffrkO7R276yzzqoqI11RUcGNN97Iueeey9tvv53v0CRJklSPZigh/SIwPCKGRERn4HvAw3vd8xDw1YjoFBHdySxne6OuRnNZrvYV4KyIWAvcD5wSEfcAH0TEAQDZnzVuAEopzU4pjU4pje7Xr18Oj2t7Xn31Vf7whz+wefNmOnfuzNe//nVOPvlky0Y3g6FDh3LGGWcwZMgQUkqUlpby6KOPsnXr1nyHJkmSpDrkslStvpmelFI5cBnwKJnE5YGU0oqIuDgiLs7e8wbwCPAa8AJwZ0ppeV3t1rtcLaV0HXAdZOpTA1ellP5HRPwcOA+Ykf35UH1tSZIkSSoczVBdjZTSImDRXudm7XX8c+DnubbZlNrHM4AHIuIHwDrgO01oS5IkSVI70wzV1VpEg5KclNKfgD9l338IjGv+kCRJkiS1B80xk9MS/BZLSZIkSQ2Wy56bfDHJkSRJktQoBbFcraNZuXIly5cvZ82aNRQVFTFmzBgOO+wwevfune/QCs4xxxzDpEmTePbZZ9m6dSuPP/44H3zwAWPHjqVbt275Dk+SJEk1aKszObmUkO6wFixYwJQpU/jTn/5ESUkJ1113HXfccQcHHXRQ/R9Wg0ydOpV77rmH0aNHs23bNq6//nquuOIKPvzww3yHJkmSpFo0tYR0S3Empx7Vp+CKioooKjIvbAkRQXFxMREBZPq9La/zlCRJ6uja8t9qJjmSJEmSGsU9OZIkSZIKijM5kiRJkgqKSY4kSZKkgpFScrmaJEmSpMLiTI4kSZKkgmKSI0mSJKlguFxNkiRJUsFxJkeSJElSQTHJkSRJklRQXK4mSZIkqWCklJzJkSRJklRYTHIkSZIkFRSTHEmSJEkFxT05kiRJkgqGe3LaseLi4qoMtbKykoqKCoqLi/McVeFJKVFeXl7V10VFRRQXFxMReY5MkiRJtWmrSU5RvgNoyyZNmsS8efM45ZRT2L17NzfccAMXXXQR69evz3doBefuu+/mnHPOobS0lB49enDjjTdyyy23sN9+++U7NEmSJNWisrKy3lc+mOTUYcSIEZx99tkMHjyYyspKli5dysKFC9m6dWu+Qys4y5Yt48EHH2TDhg2UlJQwbtw4xo8fT7du3fIdmiRJkmqxZ8laXa98cLmaJEmSpAZzT44kSZKkgmN1NUmSJEkFxZkcSZIkSQXFJEeSJElSwUgpte/lahHRB7gTOApIwAXAKmAeMBhYC0xOKX3UEkHm26hRozjjjDN4/vnn2blzJ08++SQbN27kK1/5Cl26dMl3eO3a22+/zfLly1mzZg0RwXHHHcewYcPo06dPvkOTJElSPdrqTE6uJaR/BTySUhoBHA28AVwL/DGlNBz4Y/a4IF166aXcd999jBo1io8//pgrr7ySK664wlLSzeDhhx9m0qRJLFq0iOLiYq655hp+97vfMWzYsHyHJkmSpHo0RwnpiJgQEasiYnVE1JpTRMRxEVEREZPqa7PeJCci9gG+BtyV/UV2pZT+HzARmJu9bS7wrXp/g3aqqKiIkpISIgKAiooKKisr22zm2p5UVlZSXl5e1ZfFxcV06uQqSkmSpPagqUlORBQDtwGnA0cAUyLiiFru+z/Ao7nElctMzlBgE3B3RLwSEXdGRA9gQEppQ/aX2wD0z+WBkiRJktq/PXty6nvV43hgdUppTUppF3A/mcmUvV0OzAc25hJbLv+XeSfgWODylNLzEfErGrA0LSKmAdMADj744Fw/1uZ06tSJiy66iBNPPJE77riDDz74gJ/+9KeMGjWKCy64wNmHBnrttde49957eeGFFwCYOHEiX/va1zjyyCPzHJkkSZJylePKpr4RUVrteHZKaXb2/UBgfbVrZcAJ1T8cEQOBbwOnAMfl8sBc/jIvA8pSSs9nj39PJsn5ICIOSCltiIgDqCWryv4CswFGjx7dbtd3FRcXM2nSJE4++WTmz5/P8uXLuf322xk3bhzf//73TXIa6K233mLmzJlUVFQAMG7cOC6//PI8RyVJkqSGyDHJ2ZxSGl3Ltaip2b2ObwauSSlV7Nk+Up96/zJPKb0fEesj4vCU0ipgHLAy+zoPmJH9+VBOT5QkSZLU7jVTCeky4KBqx4OAv+x1z2jg/myC0xf4ZkSUp5T+o7ZGc51+uBy4NyI6A2uA88ns53kgIn4ArAO+k2Nb7VpE0L17d7p168aOHTuoqKhg27ZtRITlpHNQWVnJjh072LFjByklSkpK6Ny5MyUlJfkOTZIkSQ3UDIW4XgSGR8QQ4D3ge8A5ez1jyJ73EfFbYGFdCQ7kmOSklJaRyaD2Ni6XzxeS3r17c/vtt7Ny5Up++MMfsmzZMs466yzOOOMMrrvuunyH1+atXLmSq666infffZfKykomT57M9OnTLRktSZLUDjU1yUkplUfEZWSqphUDc1JKKyLi4uz1WY1p140kDVRSUsKxxx5Lr169OOSQQygrK+PZZ59l0KBBlJWVsc8++7DPPvvkO8w2p6Kigk2bNrFmzRqeeeYZKioqGDhwIEcccQRjxozJd3iSJElqhGZYrkZKaRGwaK9zNSY3KaWpubSZ65eBai+HHHII8+fP5xe/+AVdunThkUceYezYscydO7f+D3dAGzZsYMqUKVx22WVs376dk08+mcWLF3PJJZfkOzRJkiQ1Qi7fkZOv75V0JqeROnfuzNChQ3nnnXeICD7++GM+/vhjtmzZku/Q2qTdu3fzzjvvsH59pkJgz549GT58OMXFxXmOTJIkSY2VrySmPiY5kiRJkhqlOZartQSXqzVR3759mTBhAl/84hcBWLVqFYsWLaKsrCzPkbUN5eXlPPvsszz55JN89tln9OnTh9NOO41jjz2WXOucS5IkqW1qq8vVTHKa6Itf/CIPPPBA1RdZzps3j29/+9s88cQTeY6sbdi+fTtXX301F198MRs3buTwww/n3//937n66qspKnL4SZIktVfuySlgEUFJSQlHHnkk//AP/0BpaSkvvfQSTzzxBLt372b8+PEMGjQo32G2upQSjz/+OCtXruS9996ja9euTJo0iVGjRtGjRw/34kiSJBUA9+QUuBNPPJETTzyRn/zkJ7z00kvMnTuXefPmsXDhwg6Z5FRWVjJr1iwefPBBAAYPHsyNN97IIYcckufIJEmS1Fza6p4ck5xm9o1vfIPOnTszf/58XnvtNe666y6ef/55LrzwQvr375/v8FrFwoULWbp0KStWrKBLly5ccMEFjBw5ki984Qv5Dk2SJEnNyJmcDmLMmDGMGTOG1atX88orr3DfffcxYMAAzjrrLPr161d1X6Ftuq8+wBcvXsyvfvUrIoLevXtz/vnnc9xxx+UxOkmSJDW3fO65qY9JTgu56KKL+OpXv8rMmTN5++23+fGPf0zv3r0BGDlyJP/4j/9Ip06F0f2lpaX8+te/rpquLC0tpaioiMsvv5yTTjqJoUOH5jlCSZIktQSXq3UwJ554Iscccwzz589nzZo1LFy4sCrTHTduHBdddBE9evSgc+fOeY608SorK9m5cyerV6/m3nvvpaKiAoCSkhJ69uzJKaecwllnnZXnKCVJktRS2upMjjV8W1Dnzp254YYbuOeeezj44IOrzr/88st8+9vf5pZbbsljdE331ltvMXnyZG644YaqBAdg2rRpPPzww5x00kl5jE6SJEktyRLSHVRRURFHH300Bx54IIMHD2bXrl0A7Ny5kyVLlnDggQeyYcOGqvv79OlDt27d8hVuvVJKbNmyper3WLNmDUuWLGH37t3sv//+VfuMjj76aE4++eR8hipJkqRW4HK1Dmzffffl7rvvrkoOli5dyqWXXsqiRYt4+eWXgUwhghkzZjBx4sR8hlqnbdu2MX36dJYtWwbAjh072LZtG+PGjePmm2+u+u6bjlJFTpIkqaNrq8vVTHJaQXFxMUOGDKk63rJlCyNHjuT9999n1apVQCbJWbFiRdUm/U6dOjF06FC6dOmSl5j32LBhA5s3bwYySc4bb7xRFXO3bt0YMWIERx55JIcffrhf8ClJktTBmOSoyujRo3nkkUe45557uPLKK4HMAJkxYwY33XQTAAMGDOA//uM/OPTQQ/MZKrNmzeI3v/kNkIlx69atVdcOP/xw5s+fT9++fU1wJEmSOpiUksvV9FclJSX07duXESNGcNppp1Wdf/3116v26FRWVrJkyRLWrFlTYxvDhw//3OxQY+3cuZMXX3yR7du313h9xYoVVTM5xcXFfOlLX6JPnz5AJsnZf//96d69e5PjkCRJUvvjTI7+xmmnncbYsWOBzAC58MILuffee4HMkrZLL7201i8N/dnPfsbVV1/d5Bg2bdrERRddxDvvvFPj9fLy8qr33bt35+c//zknnHACkCmsUFJS0uQYJEmS1D6Z5OhvFBUVfW7PzSmnnELXrl0B2L59O4sWLfrc8rDqnnvuOe66664mx7BlyxY+/PBDdu7cWeP1448/npEjRwKZPTiDBg3K+z4hSZIktQ0mOarXBRdcwAUXXABkNvy//PLLtSY5Dz74IA8++GCLxzR58mR+9KMftfhzJEmS1L64J0cN1qtXL374wx9W7YfZ2+OPP87SpUub/Jx99tmH888/n3333bfG636hpyRJkmrjTI4apGfPnlx66aW1Xt+xYwfPPPNMk5/Tp08fLr/8coYNG9bktiRJktSxmOSoWU2aNInDDjusye306NHDL++UJElSg7lcTc3u2GOP5dhjj813GJIkSerAnMmRJEmSVFBMciRJkiQVFJerSZIkSSoYKSVnciRJkiQVlraa5BTlclNEXBkRKyJieUTcFxFdI2LfiHg8Iv6c/fmFlg5WkiRJUttRWVlZ7ysf6k1yImIg8ENgdErpKKAY+B5wLfDHlNJw4I/ZY0mSJEkdxJ4la3W98iGnmRwyy9q6RUQnoDvwF2AiMDd7fS7wrWaPTpIkSVKblEuCk0uSExETImJVRKyOiL+ZOImI/x4Rr2Vfz0bE0fW1WW+Sk1J6D/gFsA7YAGxNKT0GDEgpbcjeswGo8RslI2JaRJRGROmmTZvqe5wkSZKkdqKpSU5EFAO3AacDRwBTIuKIvW57Bzg5pfRF4J+B2fXFlctytS+QmbUZAhwI9IiI/1Hf5/ZIKc1OKY1OKY3u169frh+TJEmS1MY1w56c44HVKaU1KaVdwP1kco8qKaVnU0ofZQ+fAwbV12guy9W+AbyTUtqUUtoNLABOAj6IiAMAsj835tCWJEmSpAKR40xO3z0ru7KvadWaGAisr3Zclj1Xmx8Af6gvrlxKSK8DvhwR3YHPgHFAKbANOA+Ykf35UA5tSZIkSSoADSgssDmlNLqWa1FT0zXeGDGWTJIzpr4H1pvkpJSej4jfAy8D5cArZNbB9QQeiIgfkEmEvlNfW5IkSZIKRzOUiC4DDqp2PIhMkbPPiYgvAncCp6eUPqyv0Zy+DDSl9E/AP+11eieZWR1JkiRJHVAzlIh+ERgeEUOA98h8Vc051W+IiIPJbJk5N6X0Vi6N5pTkSJIkSdLemprkpJTKI+Iy4FEy38c5J6W0IiIuzl6fBfxvYD/gNxEBUF7H8jfAJEeSJElSI6SUmmO5GimlRcCivc7Nqvb+QuDChrRpkiNJkiSpUZphuVqLMMmRJEmS1CgmOZIkSZIKRnMtV2sJJjmSJEmSGsWZHEmSJEkFxSRHkiRJUkExyZEkSZJUMNyTI0mSJKngOJMjSZIkqaCY5EiSJEkqKC5XkyRJklQwUkrO5EiSJEkqLCY5kiRJkgqKy9UkSZIkFRRnciRJkiQVDPfkSJIkSSo4LleTJEmSVFCcyZEkSZJUUExyJEmSJBWMlJLL1SRJkiQVFmdyJEmSJBUUkxxJkiRJBcUkR5IkSVLBcE+OJEmSpILjTI4kSZKkgmKSI0mSJKlguFxNkiRJUsFxJkeSJElSQWmrSU60ZmARsQnYBmxutYcKoC/2eWuyv1uffd667O/WZ5+3Pvu8ddnfreuQlFK/pjZSVFSUOnWqf85k9+7dL6WURjf1eQ3RqjM5KaV+EVHa2r9kR2efty77u/XZ563L/m599nnrs89bl/3dPqWU2uxMTlG+A5AkSZLUPu1JdOp61SciJkTEqohYHRHX1nA9IuKW7PXXIuLY+tp0T44kSZKkRmlqdbWIKAZuA04FyoAXI+LhlNLKaredDgzPvk4Abs/+rFU+ZnJm5+GZHZ193rrs79Znn7cu+7v12eetzz5vXfZ3O9UMMznHA6tTSmtSSruA+4GJe90zEfhdyngO6BMRB9TVaKsWHpAkSZJUGCLiETJFI+rTFdhR7Xh2Sml2to1JwISU0oXZ43OBE1JKl1V7zkJgRkppafb4j8A1KaXS2h7ocjVJkiRJDZZSmtAMzURNTTfins+x8IAkSZKkfCkDDqp2PAj4SyPu+ZxWS3Lqq5qg5hERayPi9YhYFhGl2XP7RsTjEfHn7M8v5DvO9iwi5kTExohYXu1crX0cEddlx/2qiBifn6jbr1r6+ycR8V52nC+LiG9Wu2Z/N1FEHBQRT0bEGxGxIiKuyJ53nLeAOvrbcd5CIqJrRLwQEa9m+/yn2fOO8RZQR387xgXwIjA8IoZERGfge8DDe93zMPD9bJW1LwNbU0ob6mq0VfbkZKsmvEW1qgnAlL2qJqgZRMRaYHRKaXO1c/8X2JJSmpFNML+QUromXzG2dxHxNeBTMhvgjsqeq7GPI+II4D4ym+oOBBYDh6WUKvIUfrtTS3//BPg0pfSLve61v5tBdjPnASmllyOiF/AS8C1gKo7zZldHf0/Gcd4iIiKAHimlTyOiBFgKXAH8Nxzjza6O/p6AY1xANsG9GSgG5qSUboiIiwFSSrOyY+hWMmNmO3B+XftxoPVmcnKpmqCWMxGYm30/l8w/nmqklNJTwJa9TtfWxxOB+1NKO1NK7wCryfz3QTmqpb9rY383g5TShpTSy9n3nwBvAANxnLeIOvq7NvZ3E2UrNH2aPSzJvhKO8RZRR3/Xxv7uYFJKi1JKh6WUhqWUbsiem5VSmpV9n1JK07PXR9aX4EDrJTkDgfXVjsuo+3/A1XgJeCwiXoqIadlzA/ZM6WV/9s9bdIWrtj527LecyyLzhWBzqi0psb+bWUQMBo4Bnsdx3uL26m9wnLeYiCiOiGXARuDxlJJjvAXV0t/gGFcLaa0kp8EVEdRoX0kpHUvmS5OmZ5f6KH8c+y3jdmAYMArYAMzMnre/m1FE9ATmA/8zpfRxXbfWcM5+b6Aa+ttx3oJSShUppVFkNjAfHxFH1XG7fd5EtfS3Y1wtprWSnAZXRFDjpJT+kv25EXiQzPTuB9k133vWfm/MX4QFq7Y+duy3gJTSB9l/MCuBf+Wvyxjs72aSXTc/H7g3pbQge9px3kJq6m/HeetIKf0/4E9k1vo7xltY9f52jKsltVaSk0vVBDVRRPTIblolInoApwHLyfT1ednbzgMeyk+EBa22Pn4Y+F5EdImIIcBw4IU8xFdQ4vPfcvxtMuMc7O9mkd3geRfwRkrpl9UuOc5bQG397ThvORHRLyL6ZN93A74BvIljvEXU1t+OcbWkVvky0JRSeURcBjzKX6smrGiNZ3cwA4AHM/9e0gn495TSIxHxIvBARPwAWAd8J48xtnsRcR/wdaBvRJQB/wTMoIY+TimtiIgHgJVAOTDd6jANU0t/fz0iRpFZvrAW+Aewv5vRV4Bzgdeza+gBfozjvKXU1t9THOct5gBgbrb6axHwQEppYUT8F47xllBbf/+bY1wtpVVKSEuSJElSa2m1LwOVJEmSpNZgkiNJkiSpoJjkSJIkSSooJjmSJEmSCopJjiRJkqSCYpIjSZIkqaCY5EiSJEkqKP8feuGPgBVXll4AAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "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": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Constructed in 0.46072816848754883 seconds\n", "Simulating to t=0.500000. Estimated 488 timesteps (dt=0.001023)\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": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/lustre/storeB/users/martinls/src/ShallowWaterGPU/data/euler_shock-bubble.nc created Tue Mar 22 13:54:57 2022 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 94706225002752'}\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" } ], "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": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzkAAADhCAYAAADrl3vGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAYZElEQVR4nO3df7SlVX3f8fdnLhdRwDjDL4dfgciojUZFR6SlMSSooG0DxmLFRokla5oWU0y0BV2rJV35h0Y0mhplTYEwNAZKBQNhWQglWupSfgxIgGGCEB1xZGQYxoiMAjP3fvvHOaQXvPf8uPecc8997vu11rPOfc6zn/1s9tqXud+19/4+qSokSZIkqSlWLHYDJEmSJGmQDHIkSZIkNYpBjiRJkqRGMciRJEmS1CgGOZIkSZIaxSBHkiRJUqMY5EiSJElaFEmOSPLlJJuTbEpyzixlkuSPkjyU5J4kr+9W717Daa4kSZIkdbUH+HBV3ZVkf+DOJDdV1f0zyrwdWNM+3gR8rv05J2dyJEmSJC2KqtpWVXe1f/4RsBk47HnFTgUur5ZbgZckWd2pXmdyJEmSJPXt5F/etx7fOdW13J33PL0JeGrGV+urav3zyyU5CjgWuO15lw4DvjvjfGv7u21zPdMgR5IkSVLfduyc4rYbD+9abnL13z5VVWs7lUmyH3A18KGqeuL5l2e5pTrVZ5AjSZIkaR6KqZpecC1JJmkFOJ+vqmtmKbIVOGLG+eHAI53qXNCenCSnJHmgnengvIXUJUmSJGnpKGAPU12PTpIEuATYXFWfnKPYdcD721nWjgd+WFVzLlWDBczkJJkA/hh4K63o6o4k1z0vE4IkSZKkBiqKqeq4aqwXJwDvA+5Ncnf7u48BRwJU1UXAl4B3AA8BPwY+0K3ShSxXOw54qKq+BZDkSlqZDwxyJEmSpGVguvPWmK6q6qvMvudmZpkCzu6n3oUEObNlOfipfNVJ1gHrAPZ9Ud7wymP2XsAjJUmSJC3Elu/uZsfOqY6BRS8K2M3C9+QMw0KCnJ6yHLTTw60HWPvafer2G4/4qZskSZIkjcZxJ3+3e6EeFAxiudpQLCTI6TvLgSRJkqTmGM95nIUFOXcAa5IcDXwPeA/w3oG0SpIkSdJYqyqeadpMTlXtSfJB4EZgAri0qjYNrGWSJEmSxlbRzJkcqupLtFK6SZIkSVpWwlTnxGiLZkFBjiRJkqTlqYDp8VytZpAjSZIkqX8FPMOKxW7GrAxyJEmSJM3LdLlcTZIkSVJDFLgnR5IkSVJzFGF3TSx2M2ZlkCNJkiSpb87kSJIkSWqYMFUmHpAkSZLUEAXsxuVqkiRJkhqiypkcSZIkSQ0z7Z4cSZIkSU1RhGdqPMOJ8WyVJEmSpLFWwDQuV5MkSZLUIFPlcjVJkiRJDVGEKWdyJEmSJDVFAbvdkyNJkiSpKYq4XE2SJElSs5h4QJIkSVJjVIXdNbHYzZiVQY4kSZKkvhUwVc7kSJIkSWoQs6tJkiRJaozC5WqSJEmSGqSAaZerSZIkSWqOMIUppCVJkiQ1ROtloOO5XG0855ckSZIkjbWqMF0ruh7dJLk0yfYk981x/WeS/EWSv06yKckHutVpkCNJkiRpXqZqRdejB5cBp3S4fjZwf1W9FjgR+ESSvTtV6HI1SZIkSX0rYHoAe3Kq6pYkR3V51P5JAuwH7AT2dKrTIEeSJElS34qwe7qnPTkHJtk443x9Va3v41GfAa4DHgH2B/5FVU13uqFrkJPkCOBy4KXAdLtRn06yCvgfwFHAFuDdVfWDPhorSZIkaQnr8WWgO6pq7QIeczJwN/ArwMuAm5L836p6Yq4bemnVHuDDVfUPgOOBs5P8PHAecHNVrQFubp9LkiRJWgaKMF3djwH4AHBNtTwEfBt4ZacbugY5VbWtqu5q//wjYDNwGHAqsKFdbANw2vzbLUmSJGkpqWqlkO52DMDDwEkASQ4BXgF8q9MNfe3JaW8IOha4DTikqrZBKxBKcvA8GixJkiRpiRrETE2SK2hlTTswyVbgfGASoKouAn4fuCzJvUCAc6tqR6c6ew5ykuwHXA18qKqeaCU36Om+dcA6gCMPM8+BJEmS1ASt5WoLfyNNVZ3R5fojwNv6qbOnViWZpBXgfL6qrml//WiS1e3rq4HtczRqfVWtraq1Bx0wnm9ElSRJktSfAnbXiq7HYuj61HY+6kuAzVX1yRmXrgPObP98JnDt4JsnSZIkaTy1ZnK6HYuhl/VjJwDvA+5Ncnf7u48BFwBXJTmL1mag04fSQkmSJEljaRAvAx2GrkFOVX0V5mz9SYNtjiRJkqSloIpeXwY6cmYCkNSXh/c8ybvu+Vc8sWufWa//m1ffwodWbhlto6Rl6sfTz/DOB36NLTtW9XzP5OQUF7/uco7fZzz/MJG0dDz7npxxZJAjqS93PHUoB50bVm26Z9brn/6Tk/jQyZeMuFXS8rRj+hme+sNDOer623u+Z+LAA7jyxjdx/OqNQ2yZpOViyS5XG6R7dx7EMVf81igfKWnAJn+4gqO3Pzjn9dU3TnLMTn/P+1GBdW+9mXMPaPXrt3c/ya/81TlM7Jxc5JZp3K3YHY558DGm+rinntzFX/3ZcRxz6NqhtUvS/EztP8W1b/uvvGbv566W+PffP5YvfvlNA3vO93Z+aiD1FIN5T84wjDTIecHWXbzsw7eO8pGShqDTH1T7X3kr+185sqY0Q8KfXvNGPrLqAQAe3L2SV3zmKerOOxe5YVoK+glwAKafeorVn/jaUNoiaWEm1vwcd735SF41+f3nfH/1vceyZoB/Q++oXYOpqMIe9+RIkmZVxapL9+P1t3wQgImfFC/9zjf7/uNVkrTEPbqDz1z4Lv5w/+fOjvzs5t2L1KDOCperSZI62OcvbuelM84NcCRp+Zl64gkOuPjri92MvrhcTZIkSVJjFLBnenFe9tmNQY4kSZKkvplCWpIkSVLjuCdHkiRJUmNUuVxNkiRJUsO4XE2SJElSY7gnR5IkSVLjTJXL1SRJkiQ1RJXL1SRJkiQ1TBnkSJIkSWoO9+RIkiRJapACpkwhLUmSJKkxqrUvZxwZ5EiSJEmal2lcriZJkiSpIYq4XE2SJElSs7hcTZIkSVKjmEJakiRJUmNUmV1NkiRJUsO4XE2SJElSo7hcTZIkSVJjFGF6TIOc8VxEJ0mSJGm8VWsmp9vRTZJLk2xPcl+HMicmuTvJpiT/p1udBjmSJEmS5qd6OLq7DDhlrotJXgJ8FvjVqnoVcHq3CnsOcpJMJPlGkuvb56uS3JTkwfbnyl7rkiRJkrT0DWImp6puAXZ2KPJe4Jqqerhdfnu3OvuZyTkH2Dzj/Dzg5qpaA9zcPpckSZK0DBQwPZ2uxwC8HFiZ5CtJ7kzy/m439BTkJDkc+CfAxTO+PhXY0P55A3Baf22VJEmStGQVUOl+wIFJNs441vX5pL2AN9CKR04G/mOSl3e7oRefAv4DsP+M7w6pqm0AVbUtycF9NlaSJEnSEtbje3J2VNXaBTxma7uOXcCuJLcArwW+OdcNXWdykvxTYHtV3TmfFiVZ92zUtpun51OFJEmSpLETarr7MQDXAr+YZK8kLwLexHO30fyUXmZyTgB+Nck7gH2AFyf5U+DRJKvbszirgVk3AFXVemA9wIuzakzfiSpJkiSpbwP46z7JFcCJtJa1bQXOByYBquqiqtqc5AbgHmAauLiq5kw3DT0EOVX1UeCj7QacCHykqn49yceBM4EL2p/Xzu8/S5IkSdKS035PzoKrqTqjhzIfBz7ea5297smZzQXAVUnOAh6mh3zVkiRJkhpkAEHOMPQV5FTVV4CvtH9+HDhp8E2SJEmStCSM6WaUhczkSJIkSVrODHIkSZIkNUYxqOxpA2eQI0mSJGl+nMmRJEmS1ChNSDwgSZIkSc+KMzmSJEmSGqMC7smRJEmS1CjO5EiSJElqFIMcSZIkSY1RuFxNkiRJUrOYeECSJElSsxjkSJIkSWoSZ3IkSZIkNYsvA5UkSZLUGIXL1SRJkiQ1S6YXuwWzM8iRJEmSND/O5EiSJElqFIMcSZIkSU2RMruaJEmSpKaZNruaJEmSpAZxJkeSJElSsxjkSJIkSWqMMoW0JEmSpKZxJkeSJElSk7gnR5IkSVKzGORIkiRJagzfkyNJkiSpcQxyJEmSJDVFGN/sait6KZTkJUm+kORvkmxO8g+TrEpyU5IH258rh91YSZIkSWOkejgWQU9BDvBp4IaqeiXwWmAzcB5wc1WtAW5un0uSJElaDtp7crod3SS5NMn2JPd1KffGJFNJ/nm3OrsGOUleDLwZuASgqp6pqr8DTgU2tIttAE7rVpckSZKkBhnMTM5lwCmdCiSZAP4LcGMvFfYyk/NzwGPAnyT5RpKLk+wLHFJV2wDanwf38kBJkiRJzZDp7kc3VXULsLNLsd8Grga299KuXoKcvYDXA5+rqmOBXfSxNC3JuiQbk2zczdO93iZJkiRp3PU2k3Pgs/FA+1jXzyOSHAa8E7io13t6ya62FdhaVbe1z79AK8h5NMnqqtqWZDVzRFVVtR5YD/DirBrTJHOSJEmS+tL7crQdVbV2AU/6FHBuVU0l6emGrkFOVX0/yXeTvKKqHgBOAu5vH2cCF7Q/r51vqyVJkiQtPSNKIb0WuLId4BwIvCPJnqr687lu6PU9Ob8NfD7J3sC3gA/QWup2VZKzgIeB0xfQcEmSJElLTC/Z0xaqqo7+++cllwHXdwpwoMcgp6ruphVBPd9JvTdPkiRJUqMMIMhJcgVwIq29O1uB84FJgKrqeR/OTL3O5EiSJEnS3+v1PTjdVNUZfZT9jV7KGeRIkiRJmp8xTStmkCNJkiRpXkaxJ2c+DHIkSZIkzc9osqv1zSBHkiRJUv8GtCdnGAxyJEmSJM2PQY4kSZKkJnEmR5IkSVKjxD05kiRJkhqjcLmaJEmSpIYxyJEkSZLUFMHlapIkSZIaJjWeUzkGOZIkSZL6554cSZIkSU3jcjVJkiRJjeJ7ciRJkiQ1i0GOJEmSpMYol6tJkiRJapDgcjVJkiRJTWMKaWl+9jrqSLaedjjTHUbrxNNw2DVb2PO9R0bXMEmSpGXOmRxpnna96hD+/Hf+gKMn95uzzD3PPMXv3vtvmTDIkSRJGo2CTC12I2Y30iDn6SP25aGPHD/KR6oB9j5sF6smJjqWOXRiiu+sm2b6NMeXpMWX3eHll+1katMDs17fc9Ib2PLPJkfcKmn09v7BCo7+4weY2vH4nGWy1148+lvH8cQxY7qDvYGevvDWwVXmTA78wsrHuP3dF43ykWqMF3a8euDEvjx44mWjaYokdbF9ahfv+vrv8qL7M+v1R9/4Av723Z8dcauk0bv8iQO58soT4fGdc5bJ3nuz+l1b+MYrvjS6hi1zx13y2MDqcrmaJEnLxM+s2JuV53yHh977qlmvv+WoO0bcImlx/OILt/CZT4YnfzL77wLAihXFhUd8YYSt0sAUZHo8oxyDHEmSBuwFmeS6NTfAmsVuibS4jp7cj9uP/Z+L3QwN03jGOAY5kiRJkvrne3IkSZIkNUuVy9UkSZIkNcx4xjgGOZIkSZLmZ1yXq63opVCS30myKcl9Sa5Isk+SVUluSvJg+3PlsBsrSZIkaUwUMFXdj0XQNchJchjw74C1VfVqYAJ4D3AecHNVrQFubp9LkiRJWiZS3Y/F0NNMDq1lbS9MshfwIuAR4FRgQ/v6BuC0gbdOkiRJ0viq6n50keTSJNuT3DfH9X+Z5J728bUkr+1WZ9cgp6q+B1wIPAxsA35YVX8JHFJV29pltgEHz9GodUk2Jtn42ONT3R4nSZIkaYkY0EzOZcApHa5/G/ilqnoN8PvA+m4V9rJcbSWtWZujgUOBfZP8ei+tBaiq9VW1tqrWHnTARK+3SZIkSRpjKch0dT26qapbgJ0drn+tqn7QPr0VOLxbnb1kV3sL8O2qegwgyTXAPwIeTbK6qrYlWQ1s76EuSZIkSU0x3VOpA5NsnHG+vqq6zsbM4Szgf3Ur1EuQ8zBwfJIXAT8BTgI2AruAM4EL2p/XzrOhkiRJkpag9LDnBthRVWsX/Kzkl2kFOf+4W9muQU5V3ZbkC8BdwB7gG7TWwe0HXJXkLFqB0OkLabQkSZKkJaQKeliONghJXgNcDLy9qh7vVr6nl4FW1fnA+c/7+mlaszqSJEmSlqFRpIhOciRwDfC+qvpmL/f0FORIkiRJ0k/pbblaR0muAE6ktXdnK63JlclW9XUR8J+AA4DPJgHY0235m0GOJEmSpP4VZGrhQU5VndHl+m8Cv9lPnQY5kiRJkuZnNFty+maQI0mSJGleesyuNnIGOZIkSZL6V8AAlqsNg0GOJEmSpL6FciZHkiRJUsMY5EiSJElqFIMcSZIkSY0xoBTSw2CQI0mSJGl+nMmRJEmS1BxlkCNJkiSpQUwhLUmSJKlpTCEtSZIkqVkMciRJkiQ1RhVMTS92K2ZlkCNJkiRpfpzJkSRJktQoBjmSJEmSGqMKpqYWuxWzMsiRJEmSND/O5EiSJElqjAKmDXIkSZIkNcm02dUkSZIkNUa5XE2SJElSgxTO5EiSJElqGGdyJEmSJDVHwZQzOZIkSZKaoqDKIEeSJElSk5hCWpIkSVJjVMHU1GK3YlYGOZIkSZLmx8QDkiRJkpqkTCENd97z9I6J1Q/tAnaM8rniQOzzUbK/R88+Hy37e/Ts89Gzz0fL/h6tnx1ILWV2NQCq6qAkG6tq7Sifu9zZ56Nlf4+efT5a9vfo2eejZ5+Plv29hI1pdrUVi90ASZIkSUtPATVdXY9uklyaZHuS++a4niR/lOShJPckeX23Og1yJEmSJPWvipqa6nr04DLglA7X3w6saR/rgM91q3Axgpz1i/DM5c4+Hy37e/Ts89Gyv0fPPh89+3y07O+lqqa7H92qqLoF2NmhyKnA5dVyK/CSJKs71Zka07RvkiRJksZXkhtoJY3oZh/gqRnn66vqOYFtkqOA66vq1bM853rggqr6avv8ZuDcqto41wNNIS1JkiSpb1XVaYnZIGW2x3e6wT05kiRJksbZVuCIGeeHA490umFkQU6SU5I80M6KcN6onrvcJNmS5N4kdyfZ2P5uVZKbkjzY/ly52O1cymbLANKpj5N8tD3uH0hy8uK0eumao79/L8n32uP87iTvmHHN/l6gJEck+XKSzUk2JTmn/b3jfAg69LfjfEiS7JPk9iR/3e7z/9z+3jE+BB362zGuXl0HvL+dZe144IdVta3TDSPZk5NkAvgm8FZakdgdwBlVdf/QH77MJNkCrK2qHTO++wNgZ1Vd0A4wV1bVuYvVxqUuyZuBJ2ltgHt1+7tZ+zjJzwNXAMcBhwL/G3h5VfWUakRz9vfvAU9W1YXPK2t/D0B7M+fqqroryf7AncBpwG/gOB+4Dv39bhznQ5EkwL5V9WSSSeCrwDnAr+EYH7gO/X0KjnEBSa4ATqS1v+dR4HxgEqCqLmqPoc/QGjM/Bj7QaT8OjG4m5zjgoar6VlU9A1xJK0uCRuNUYEP75w20/vHUPM2RAWSuPj4VuLKqnq6qbwMP0fp9UI96yLgyk/09AFW1raruav/8I2AzcBiO86Ho0N9zsb8XqJ2h6cn26WT7KBzjQ9Ghv+dify8zVXVGVa2uqsmqOryqLqmqi6rqovb1qqqzq+plVfUL3QIcGF2Qcxjw3RnnW+n8P3DNXwF/meTOJOva3x3y7JRe+/PgRWtdc83Vx4794flgWi8Eu3TGkhL7e8DSynZzLHAbjvOhe15/g+N8aJJMJLkb2A7cVFWO8SGao7/BMa4hGVWQ03dGBM3bCVX1elovTTq7vdRHi8exPxyfA14GvA7YBnyi/b39PUBJ9gOuBj5UVU90KjrLd/Z7n2bpb8f5EFXVVFW9jtYG5uOS/FTa2hns8wWao78d4xqaUQU5fWdE0PxU1SPtz+3AF2lN7z7aXvP97Nrv7YvXwsaaq48d+0NQVY+2/8GcBv4b/38Zg/09IO1181cDn6+qa9pfO86HZLb+dpyPRlX9HfAVWmv9HeNDNrO/HeMaplEFOXcAa5IcnWRv4D20siRogJLs2960SpJ9gbcB99Hq6zPbxc4Erl2cFjbaXH18HfCeJC9IcjSwBrh9EdrXKHnuW47fSWucg/09EO0NnpcAm6vqkzMuOc6HYK7+dpwPT5KDkryk/fMLgbcAf4NjfCjm6m/HuIZpJC8Drao9ST4I3AhMAJdW1aZRPHuZOQT4YuvfS/YC/qyqbkhyB3BVkrOAh4HTF7GNS97MDCBJttLKAHIBs/RxVW1KchVwP7AHONvsMP2Zo79PTPI6WssXtgD/GuzvAToBeB9wb3sNPcDHcJwPy1z9fYbjfGhWAxva2V9XAFdV1fVJvo5jfBjm6u//7hjXsIwkhbQkSZIkjcrIXgYqSZIkSaNgkCNJkiSpUQxyJEmSJDWKQY4kSZKkRjHIkSRJktQoBjmSJEmSGsUgR5IkSVKj/D/+cAtALsfBNQAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "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": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Constructed in 0.01354837417602539 seconds\n", "Simulating to t=10.000000. Estimated 4346 timesteps (dt=0.002301)\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": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/lustre/storeB/users/martinls/src/ShallowWaterGPU/data/euler_kelvin_helmholtz.nc created Tue Mar 22 14:00:36 2022 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 94706225002752'}\n", "0% [##########====================] 100%. Total: 16s, elapsed: 5s, remaining: 10s\n", "0% [################==============] 100%. Total: 19s, elapsed: 10s, remaining: 9s\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "Animation size has reached 21325363 bytes, exceeding the limit of 20971520.0. If you're sure you want a larger animation embedded, set the animation.embed_limit rc parameter to a larger value (in MB). This and further frames will be dropped.\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" } ], "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": "iVBORw0KGgoAAAANSUhEUgAAA6cAAAJWCAYAAACkpfF4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzt3X2wZHd5J/bv03f0YsSLRuLFsiRKIpYds46xsRbkfXER5IAgrOVUIIE4thYrpXIWOxDsLMKuLZJNagsnqcW44tirAlZQ5ViwmBjFxRrLvKw3tUYGYSwLawVjTNCAbBkQGL8gIc0vf/S5fXvu3Jm5c/vld+6dz6eq1X1On779dN/zaPq5z3NOV2stAAAA0NOkdwAAAACgOAUAAKA7xSkAAADdKU4BAADoTnEKAABAd4pTAAAAujttcVpVb6uqB6vqnrl1/1tV/fuquruq/u+qunDuvtdX1ZGquq+qXji3/rph3ZGqunn5LwXOTnIUxk2OwrjJURiP3XROb01y3bZ1dyT5ztbadyX5VJLXJ0lVPTPJy5P8reEx/2dVbVTVRpJfTPKiJM9M8ophW2Bxt0aOwpjdGjkKY3Zr5CiMwmmL09ba7yT58rZ1v9Vae3RY/EiSy4bb1ye5rbX2cGvtT5IcSfKc4XKktfaZ1tojSW4btgUWJEdh3OQojJschfFYxjGnP5bkXw+3L01y/9x9R4d1J1t/gqq6qao+NlxuWkJ8cLaTozBuchTGTY7Cmhxa5MFV9bNJHk3yK5urdtisZeciuO30M1trtyS5JUmefNFGu/pZ5/+LRWKEMbnr7oe/2Fp7yrqeb9U5uvH4C9p5T79cjnJgPHL/UTkKI3bQcvTcOq89sS6SoxwYX8tDC+XonovTqrohyUuSXNta20y+o0kun9vssiRfGG6fbP1JXXH5Ofm9919+us1g39i45Mj/t67nWkeOHrroonzLT71mCdHCOHz2NT8tR2HEDlqOnp8L8tzJDywhWhiH3z72rxbK0T2N9VbVdUlel+QHW2t/PXfX7UleXlXnVdWVSa5K8ntJPprkqqq6sqrOzfRA8tsXCRw4OTkK4yZHYdzkKPRx2s5pVf1qkucleXJVHU3yhkzPWHZekjuqKkk+0lr78dbaJ6vqXUn+KNMRiFe11h4bfs5PJHl/ko0kb2utfXIFrwfOOnIUxk2Owrh1z9Faxilg4GCorSmF8bn6Wec3Y70cJBuXHLmrtXZ17ziW5bynX96MDHKQfPY1Py1HYcQOWo4+cXJRu+bQC0+/IewTd3zjtoVy1J9qAAAA6G6hs/UCAAB7VcZ6YY7iFAAAepns9O00cHbypxoAAAC60zkFAIAOKslwNmAgilMAAOijKpkYZIRNsgEAAIDudE4BAKAXY70wozgFAIBejPXCjGwAAACgO51TAADooZytF+bpnAIAANCdzikAAHRRyUTnFDYpTgEAoJcyyAibZAMAAADd6ZwCAEAPFWO9MEfnFAAAgO50TgEAoBdfJQMzilMAAOiinBAJ5sgGAAAAutM5BQCAXpwQCWYUpwAA0EMl5ZhTmDHWCwAAQHc6pwAA0IvOKcwoTgEAoItKJgYZYZNsAAAAoDudUwAA6MVYL8woTgEAoIeK4hTmGOsFAACgO51TAADoZaJzCpt0TgEAAOhO5xQAALoox5zCHMUpAAD0ojiFGWO9AAAAdKdzCgAAPfgqGTiOzikAAADd6ZwCAEAnzVfJwIziFAAAejHWCzPGegEAAOhO5xQAALrwPacwT3EKAAA9OFsvHMdYLwAAAN3pnAIAQC/O1gszilMAAOjFWC/MGOsFAACgO51TAADoRecUZhSnAADQQ1Wa4hRmjPUCAADQnc4pAAD0olUEM9IBAACA7nROAQCgF8ecwoziFAAAelGcwoyxXgAAALrTOQUAgA5axVfJwBydUwAAALrTOQUAgF60imBGcQoAAL0Y64UZf6sBAACgO51TAADoonROYY7iFAAAeqjpGXuBKWO9AAAAdKdzCgAAvUy0TmGTzikAAADd6ZwCAEAnzQmRYEZxCgAAvahNYcZYLwAAAN3pnAIAQA8V33MKcxSnAADQie85hS3GegEAAOhO5xQAAHox1gszilMAAOjFHCPMSAcAAAC60zkFAIAeqtKM9cKM4hQAAHpRm8KMsV4AAAC60zkFAIAOWmKsF+bonAIAANCdzikAAPSicQozilMAAOikKU5hxlgvAAAA3emcAgBAD5VkonUKm07bOa2qt1XVg1V1z9y6i6rqjqr69HB9eFhfVfULVXWkqu6uqmfPPeaGYftPV9UNq3k5cPaRozBuchTGTY7CeOxmrPfWJNdtW3dzkg+01q5K8oFhOUlelOSq4XJTkl9Kpgme5A1JnpvkOUnesJnkwMJujRyFMbs1chTG7NZ0zNFWLi4H57Ko0xanrbXfSfLlbauvT/L24fbbk/zQ3Pp3tKmPJLmwqi5J8sIkd7TWvtxaeyjJHTnxfwLAHshRGDc5CuPWPUerXFwOzmVBez0h0tNaaw8kyXD91GH9pUnun9vu6LDuZOuB1ZCjMG5yFMZNjkIHyz4h0k7lcjvF+hN/QNVNmY5J5OmXOl8TLNlSc3TjsKlCWDI5CuO21Bw97/wLlzIKCQfFXjunfzaMMGS4fnBYfzTJ5XPbXZbkC6dYf4LW2i2ttatba1c/5eKNPYYHZ7215OjG4y9YeuBwlpCjMG5rydFzzrtgWta6uByUy4L2WpzenmTzLGQ3JHnv3PofHc5kdk2Srw6jEO9P8oKqOjwcHP6CYR2wGnIUxk2OwrjJUejgtHOzVfWrSZ6X5MlVdTTTM5G9Mcm7qurGJJ9L8rJh8/cleXGSI0n+Oskrk6S19uWq+p+TfHTY7p+21rYfeA7sgRyFcZOjMG7dc3QJ3SY4KE5bnLbWXnGSu67dYduW5FUn+TlvS/K2M4oOOC05CuMmR2HceudoW8IZTuGg2OtYLwAAACyN0+ECAEAvGqcwozgFAIAOWiXNHCPMSAcAAAC60zkFAIBejPXCjM4pAAAA3emcAgBAF+WrZGCO4hQAAHpRm8KMsV4AAAC60zkFAIAeavp1MsCUzikAAADd6ZwCAEAvOqcwozgFAIBOjPXCFmO9AAAAdKdzCgAAvfieU5hRnAIAQCfGemGLsV4AAAC60zkFAIAeKs7WC3MUpwAA0EGLsV6YZ6wXAACA7nROAQCgF51TmFGcAgBAJ8Z6YYuxXgAAALrTOQUAgF50TmFG5xQAAIDudE4BAKCHSlppncImxSkAAPSiNoUZY70AAAB0p3MKAACd+CoZ2KJzCgAAQHc6pwAA0IvOKcwoTgEAoBNjvbDFWC8AAADd6ZwCAEAPFWO9MEfnFAAAgO50TgEAoIMWx5zCPMUpAAD0ojiFGWO9AAAAdKdzCgAAveicwoziFAAAOnHMKWwx1gsAAEB3OqcAANCD7zmF4yhOAQCgE2O9sMVYLwAAAN3pnAIAQC86pzCjOAUAgE5aqU5hk7FeAAAAutM5BQCAXjROYUbnFAAAgO50TgEAoAffcwrHUZwCAEAnvucUthjrBQAAoDudUwAA6EXnFGZ0TgEAAOhO5xQAADpoccwpzFOcAgBAL4pTmDHWCwAAQHc6pwAA0IPvOYXjKE4BAKATx5zCFmO9AAAAdKdzCgAAveicwoziFAAAOjHWC1uM9QIAANCdzikAAPSicwozilMAAOjBV8nAcYz1AgAA0J3OKQAAdOKESLBF5xQAAIDudE4BAKAXnVOYUZwCAEAnxnphi7FeAAAAulOcAgAA0J2xXgAA6MVYL8zonAIAANCdzikAAHTQygmRYJ7OKQAAAN3pnAIAQC86pzCjOAUAgF4UpzBjrBcAAIDudE4BAKATJ0SCLYpTAADoRXEKM8Z6AQAA6G6h4rSq/vuq+mRV3VNVv1pV51fVlVV1Z1V9uqreWVXnDtueNywfGe6/YhkvADg5OQrjJkdh3Faeo+XicsAuC9pzcVpVlyb575Jc3Vr7ziQbSV6e5OeSvKm1dlWSh5LcODzkxiQPtda+Ncmbhu2AFZGjMG5yFMZtXTnaysXl4FwWtehY76Ek31RVh5I8LskDSZ6f5N3D/W9P8kPD7euH5Qz3X1tVS3gJwCnIURg3OQrjJkdhjfZ8QqTW2uer6n9P8rkkf5Pkt5LcleQrrbVHh82OJrl0uH1pkvuHxz5aVV9NcnGSL87/3Kq6KclNSfL0S52vCfZqHTm6cfjwql8GHFhyFMZtHTl66EmHk2qrfimwbywy1ns4078QXZnkW5JckORFO2y6mXE7/eXohGxsrd3SWru6tXb1Uy7e2Gt4cNZbR45uPP6CZYULZx05CuO2lhy94ILuY5guLsu8LGqRsd4fSPInrbU/b619I8l7kvydJBcOow9JclmSLwy3jya5PEmG+5+U5MsLPD9wanIUxk2OwrjJUVizRYrTzyW5pqoeN8zTX5vkj5J8KMlLh21uSPLe4fbtw3KG+z/YWjPHAKsjR2Hc5CiM23pytFxcDtBlQXsuTltrd2Z6sPfHk/zh8LNuSfK6JK+tqiOZztm/dXjIW5NcPKx/bZKbF4gbOA05CuMmR2Hc5Cis30JnHGqtvSHJG7at/kyS5+yw7deTvGyR5wPOjByFcZOjMG4rz9EldZvgoHA6XAAA6GQZJ5GBg2LR7zkFAACAhemcAgBALzqnMKNzCgAAQHc6pwAA0IvOKcwoTgEAoIM2XIApY70AAAB0p3MKAAC9GOuFGZ1TAAAAutM5BQCAHio6pzBHcQoAAJ00xSnMGOsFAACgO51TAADopXyZDGxSnAIAQC/GemHGWC8AAADd6ZwCAEAnTogEWxSnAADQg6+SgeMY6wUAAKA7nVMAAOhF5xRmFKcAANCJY05hi7FeAAAAutM5BQCAXqr1jgBGQ+cUAACA7nROAQCgF8ecwoziFAAAevA9p3AcY70AAAB0p3MKAAC96JzCjM4pAAAA3emcAgBAFy3NV8nAjOIUAAB6MdYLM8Z6AQAA6E7nFAAAetE5hRmdUwAAALpTnAIAANCdsV4AAOihkjhbL8woTgEAoBfHnMKMsV4AAAC60zkFAIBedE5hRnEKAAC9OOYUZoz1AgAA0J3OKQAA9GKsF2Z0TgEAAOhO5xQAAHqo6JzCHMUpAAD04oRIMGOsFwAAgO50TgEAoJMy1gszOqcAAAB0p3MKAAC9OOYUZhSnAADQg7P1wnGM9QIAANCdzikAAHTRUsZ6YUZxCgAAvRjrhRljvQAAAHSncwoAAJ34nlPYojgFAIBeHHMKM8Z6AQAA6E7nFAAAOjHVC1sUpwAA0EFVfJUMzDHWCwAAQHc6pwAA0IvOKczonAIAANCdzikAAHTie05hi+IUAAA6cUIk2GKsFwAAgO50TgEAoBNjvbBFcQoAAD1UM9YLc4z1AgAA0J3OKQAAdFBxQiSYp3MKAABAdzqnAADQiRMiwRbFKQAAdGKsF7YY6wUAAKA7nVMAAOhE5xS26JwCAADQnc4pAAB04nxIsEVxCgAAHVQlE2O9MGOsFwAAgO50TgEAoIvmhEgwR3EKAACdKE5hi7FeAAAAutM5PcBe+C3fvfSf+f4vfGLpPxMA4GxUSSZO17tnq+g6t+YX0tNCxWlVXZjkLUm+M0lL8mNJ7kvyziRXJPlskv+itfZQVVWSNyd5cZK/TvIPW2sfX+T5zyarKDT3Yi9xKGj7kaMwbnIUxm0dOWqsd1z8PvpatHP65iS/2Vp7aVWdm+RxSX4myQdaa2+sqpuT3JzkdUlelOSq4fLcJL80XJ/1xlJ4rspuXp8CdmXkKIybHIVxk6NL4OtykmM6sruy5+K0qp6Y5PuT/MMkaa09kuSRqro+yfOGzd6e5MOZJuz1Sd7RWmtJPlJVF1bVJa21B/Yc/T5y0AvQRZ3q/VG47o0chXGTozBua8lR33N61vB73p1FOqfPSPLnSf5lVT0ryV1JXp3kaZtJ2Fp7oKqeOmx/aZL75x5/dFh3XMJW1U1JbkqSp1867kNiFZzrcSbvs0L2OCvP0Y3Dh1f6AuCAk6MwbivP0XOf+sRRj5E6/nJ8xry/LMMi1d+hJM9O8pOttTur6s2ZjjWczE579wnvbmvtliS3JMnVzzq/67uv+Nx/Tvc7O8uK15Xn6HlPv/xg/x8SVkuOwritPEcv+LZLRp2jB70QYnwWKU6PJjnaWrtzWH53pgn7Z5sjDFV1SZIH57a/fO7xlyX5wgLPvzDF59nnLCte932OwgEnR2HcVp6jlbbScU/HObLd2MeL91ycttb+tKrur6pvb63dl+TaJH80XG5I8sbh+r3DQ25P8hNVdVumB4d/dVnHySgyWZa97ktjLGrHlKPAieQojNuYclSRybLsdV9aV1G76EGdP5nkV4azl30mySuTTJK8q6puTPK5JC8btn1fpqfWPpLp6bVfebof/qm7H6fwZF/Y/X56ZKVx7GClOQosTI7CuI0iR8fe7YJlWag4ba19IsnVO9x17Q7btiSvWuT5gDMjR2Hc5CiMmxyF9Zr0DgAAAADG/V0tAABwgDmaFLbonAIAANCd4hQAAIDujPUCAEAnzsQLW3ROAQAA6E5xCgAAQHeKUwAAALpzzCkAAHRQccwpzNM5BQAAoDvFKQAAAN0Z6wUAgB4qKWO9MKNzCgAAQHeKUwAAALoz1gsAAJ04Wy9s0TkFAACgO8UpAAAA3RnrBQCADirNWC/M0TkFAACgO8UpAAAA3SlOAQAA6M4xpwAA0MkkjjmFTTqnAAAAdKc4BQAAoDtjvQAA0En5KhmY0TkFAACgO8UpAAAA3RnrBQCADirJxFgvzOicAgAA0J3iFAAAgO6M9QIAQA9lrBfm6ZwCAADQneIUAACA7oz1AgBAB5VmrBfm6JwCAADQneIUAACA7oz1AgBAJ5MY64VNOqcAAAB0pzgFAACgO8UpAAAA3TnmFAAAOpnUsd4hwGjonAIAANCd4hQAAIDujPUCAEAHlWRSvkoGNumcAgAA0J3iFAAAgO6M9QIAQBfNWC/M0TkFAACgO8UpAAAA3RnrBQCADqqcrRfm6ZwCAADQneIUAACA7oz1AgBAJ5MY64VNOqcAAAB0pzgFAACgO2O9AADQQcXZemGezikAAADdKU4BAADoTnEKAABAd445BQCALlomdax3EDAaOqcAAAB0pzgFAACgO2O9AADQQSXZ8FUyMKNzCgAAQHeKUwAAALoz1gsAAJ1MYqwXNumcAgAA0J3iFAAAgO6M9QIAQCcTZ+uFGZ1TAAAAulOcAgAA0J2xXgAA6KCqZVLHeocBo6FzCgAAQHeKUwAAALoz1gsAAJ1sOFsvzOicAgAA0J3iFAAAgO4UpwAAAHTnmFMAAOigkkzimFPYpHMKAABAd4pTAAAAult4rLeqNpJ8LMnnW2svqaork9yW5KIkH0/yI621R6rqvCTvSPK9Sb6U5L9srX120ecHTk2OwrjJURi3VefopI6tLHbYb5bROX11knvnln8uyZtaa1cleSjJjcP6G5M81Fr71iRvGrYDVk+OwrjJURg3OQprslBxWlWXJflPk7xlWK4kz0/y7mGTtyf5oeH29cNyhvuvHbYHVkSOwrjJURg3OQrrtehY788n+cdJnjAsX5zkK621R4flo0kuHW5fmuT+JGmtPVpVXx22/+L8D6yqm5LclCTn53ELhgdnvZXm6MbhwysNHs4CchTGbaU5+oRvflwm5Wy9sGnPndOqekmSB1trd82v3mHTtov7tla0dktr7erW2tXn5Ly9hgdnvXXk6MbjL1hCpHB2kqMwbuvI0W867LMuzFukc/p3k/xgVb04yflJnpjpX5curKpDw1+ULkvyhWH7o0kuT3K0qg4leVKSLy/w/MCpyVEYNzkK4yZHYc323Dltrb2+tXZZa+2KJC9P8sHW2g8n+VCSlw6b3ZDkvcPt24flDPd/sLVmjgFWRI7CuMlRGLd15ehGjrm4HJjLolbxPaevS/LaqjqS6Zz9W4f1b01y8bD+tUluXsFzA6cnR2Hc5CiMmxyFFVn4e06TpLX24SQfHm5/Jslzdtjm60letoznA86MHIVxk6MwbnIU1mMpxSkAAHBmKnG2XpizirFeAAAAOCOKUwAAALoz1gsAAD2UsV6Yp3MKAABAd4pTAAAAulOcAgAA0J1jTgEAoINKy0aO9Q4DRkPnFAAAgO4UpwAAAHRnrBcAADrxVTKwRecUAACA7hSnAAAAdGesFwAAOqi0bJSz9cImnVMAAAC6U5wCAADQnbFeAADoZBJn64VNOqcAAAB0pzgFAACgO2O9AADQQSXO1gtzdE4BAADoTnEKAABAd8Z6AQCgE2frhS06pwAAAHSnOAUAAKA7Y70AANBBpTlbL8zROQUAAKA7xSkAAADdKU4BAADozjGnAADQQyUTx5zCjM4pAAAA3SlOAQAA6M5YLwAAdFBJNtJ6hwGjoXMKAABAd4pTAAAAujPWCwAAXTRn64U5OqcAAAB0pzgFAACgO2O9AADQgbP1wvF0TgEAAOhOcQoAAEB3xnoBAKATZ+uFLTqnAAAAdKc4BQAAoDtjvQAA0EGlZSPGemGTzikAAADdKU4BAADoTnEKAABAd445BQCATibVeocAo6FzCgAAQHeKUwAAALoz1gsAAB1U4qtkYI7OKQAAAN0pTgEAAOjOWC8AAHTRslHGemGTzikAAADdKU4BAADoTnEKAABAd445BQCADirJxFfJwIzOKQAAAN0pTgEAAOjOWC8AAPRQyUa13lHAaOicAgAA0J3iFAAAgO6M9QIAQAeVlg1n64UZnVMAAAC6U5wCAADQnbFeAADoZFLGemGTzikAAADdKU4BAADozlgvAAB0UEk20nqHAaOhcwoAAEB3ilMAAAC6M9YLAAAdVFo2nK0XZnROAQAA6E5xCgAAQHeKUwAAALpzzCkAAHQyiWNOYZPOKQAAAN0pTgEAAOjOWC8AAHTgq2TgeHvunFbV5VX1oaq6t6o+WVWvHtZfVFV3VNWnh+vDw/qqql+oqiNVdXdVPXtZLwI4kRyFcZOjMG5yFNZvkbHeR5P8VGvtO5Jck+RVVfXMJDcn+UBr7aokHxiWk+RFSa4aLjcl+aUFnhs4PTkK4yZHYdzkKKzZnsd6W2sPJHlguP21qro3yaVJrk/yvGGztyf5cJLXDevf0VprST5SVRdW1SXDzwGWTI7CuMlRGLd15ehG2mpeAOxDSzkhUlVdkeR7ktyZ5GmbSThcP3XY7NIk98897OiwDlgxOQrjJkdh3OQorMfCxWlVPT7JryV5TWvtL0616Q7rTvhTUVXdVFUfq6qPfSMPLxoenPVWmaOP/eVfLStMOGvJURi3VeboV7/82LLChANhobP1VtU5mSbrr7TW3jOs/rPNEYaquiTJg8P6o0kun3v4ZUm+sP1nttZuSXJLkjyxLjLnAAtYdY6e9/TL5SgsQI7CuK06R7/tP/qmNomz9cKmRc7WW0nemuTe1to/n7vr9iQ3DLdvSPLeufU/OpzJ7JokX3WcDKyOHIVxk6MwbnIU1m+RzunfTfIjSf6wqj4xrPuZJG9M8q6qujHJ55K8bLjvfUlenORIkr9O8soFnhs4PTkK4yZHYdzkKKzZImfr/X+z82x9kly7w/Ytyav2+nzAmZGjMG5yFMZtPTnaslHGemHTUs7WCwAAAItQnAIAANDdQmfrBQAA9qaSbJz4bTNw1tI5BQAAoDvFKQAAAN0pTgEAAOjOMacAANBDJRNfJQMzOqcAAAB0pzgFAACgO2O9AADQQaVlI8Z6YZPOKQAAAN0pTgEAAOjOWC8AAHSykdY7BBgNnVMAAAC6U5wCAADQnbFeAADooJJMytl6YZPOKQAAAN0pTgEAAOjOWC8AAHTibL2wRecUAACA7hSnAAAAdGesFwAAOqg0Y70wR+cUAACA7hSnAAAAdGesFwAAOpmUsV7YpHMKAABAd4pTAAAAulOcAgAA0J1jTgEAoINKfJUMzNE5BQAAoDvFKQAAAN0Z6wUAgE4mxnphRucUAACA7hSnAAAAdGesFwAAOqgkG2WsFzbpnAIAANCd4hQAAIDujPUCAEAHlZYNZ+uFGZ1TAAAAulOcAgAA0J2xXgAA6ESnCLbIBwAAALpTnAIAANCdsV4AAOigkmxU7yhgPHROAQAA6E5xCgAAQHeKUwAAALpzzCkAAHSiUwRb5AMAAADdKU4BAADozlgvAAB0UEk2egcBI6JzCgAAQHeKUwAAALoz1gsAAD1UZaOqdxQwGjqnAAAAdKc4BQAAoDtjvQAA0EFFpwjmyQcAAAC6U5wCAADQnbFeAADoZCPO1gubdE4BAADoTnEKAABAd8Z6AQCgg0oyKWO9sEnnFAAAgO4UpwAAAHSnOAUAAKA7x5wCAEAnvkoGtuicAgAA0J3iFAAAgO6M9QIAQAeVykSvCGZkAwAAAN0pTgEAAOjOWC8AAHSyUc7WC5t0TgEAAOhOcQoAAEB3xnoBAKCDSpytF+bIBgAAALpTnAIAANCdsV4AAOhkEmfrhU06pwAAAHSnOAUAAKA7Y70AANBBpbJRekWwae3ZUFXXVdV9VXWkqm5e9/MDpyZHYdzkKIybHIW9W2txWlUbSX4xyYuSPDPJK6rqmeuMATg5OQrjJkdh3OQoLGbdndPnJDnSWvtMa+2RJLcluX7NMQAnJ0dh3OQojJschQWs+5jTS5PcP7d8NMlz5zeoqpuS3DQs/uVvt3d/KckX1xPe0jw5Yl6H/Rjzt/cO4DTOOEc/+5qflqPrIeb1OGg5+vBnX/PT96wptmXaj/uOmNfjwOXoxiVH5Oh6iHk9FsrRdRenO32RUztuobVbktwye0DVx1prV686sGUS83rs15h7x3AacnSkxLweBy1H9+PvINmfcYt5PeToOOzHuMW8Hovm6LrHeo8muXxu+bIkX1hzDMDJyVEYNzkK4yZHYQHrLk4/muSqqrqyqs5N8vIkt685BuDk5CiMmxyFcZOjsIC1jvW21h6tqp9I8v4kG0ne1lr75Gkedstp7h8jMa+HmJdMjo6amNdj1DHvIUdH/XpOYT/GLeb1GHXMcnTUxLweC8VcrbXTbwUAAAArtO6xXgAAADiB4hQjfPqjAAAQJUlEQVQAAIDuRlucVtV1VXVfVR2pqpt7x7Opqt5WVQ9W1T1z6y6qqjuq6tPD9eFhfVXVLwyv4e6qenanmC+vqg9V1b1V9cmqevXY466q86vq96rqD4aY/6dh/ZVVdecQ8zuHkw2kqs4blo8M91+x7pjnYt+oqt+vqt/YLzHvhRxdasxydL2xy9GO5OjaYpajIydHlxqzHF1v7CvL0VEWp1W1keQXk7woyTOTvKKqntk3qplbk1y3bd3NST7QWrsqyQeG5WQa/1XD5aYkv7SmGLd7NMlPtda+I8k1SV41vJ9jjvvhJM9vrT0ryXcnua6qrknyc0neNMT8UJIbh+1vTPJQa+1bk7xp2K6XVye5d255P8R8RuTo0snR9ZKjfd0aOboOcnTE5OjSydH1Wl2OttZGd0nyfUneP7f8+iSv7x3XXDxXJLlnbvm+JJcMty9Jct9w+18kecVO23WO/71J/pP9EneSxyX5eJLnJvlikkPb95NMz4r3fcPtQ8N21SHWyzL9n9/zk/xGpl/GPeqY9/g65ehq45ejq4tVjo7gIkfXHq8cHdlFjq48fjm6ulhXmqOj7JwmuTTJ/XPLR4d1Y/W01toDSTJcP3VYP7rXMbTTvyfJnRl53MPIwCeSPJjkjiR/nOQrrbVHd4hrFvNw/1eTXLzeiJMkP5/kHyc5NixfnPHHvBej2EfOwKj39XlydOXk6DiNel+fJ0dXTo6O06j39XlydOVWmqNjLU5rh3X78TtvRvU6qurxSX4tyWtaa39xqk13WLf2uFtrj7XWvjvTv9A8J8l37LTZcN095qp6SZIHW2t3za/eYdPRxLyA/Rz7vFG9Djm6WnJ038Q+b1SvQ46ulhzdN7HPG9XrkKOrtY4cHeX3nFbV9yX5H1trLxyWX39xnvbPHskjSdX0VW79J6nt1/P35fj7Use/TXX8dm37/Tv8jDa/7oTbtfWO77DNyR57qsdsaieL63Q/Z9vyKX/OccvtlPFMl9vsdm27b+6dmHubt2/fTrhvfrk2n6ba1q9hp8cM645bTlK1uXzy+ybb18/iaFsxbotnPvZKzcU2vzy9vuvuh7+R5HFt+sXcs327qt4/3P7dqjqU5E+TPKWNMSm32TlHv/mfPZKHUzvm4elydH4HOsMc3WHbseboSfNz27p2iuc4ft3qcrTmnmPvOXriupPn6Lb8m7vvdDm6uTyZf+xc7HL0ZDm67fq4dbP/nDpHT9juVDlac9vkJNtk5xzdRR6eLkdPeM4dtt3Nv6HJaXL0uOU95uj2fzOz+xytbY/fTY7O/8zd5+hWTh/3c4/7eafO0RreHzm6yxzdKS/nbx+0HN0hX8aao8f9b3SXOVrb8uFMcnT2b1x2n6Pbczs58xzdim91OXroZHd09tEkV1XVlUk+n+Tlj+SRXHPoBUlNUhsbyaSm11XJTstVyWSSTGrHdW0ymfaNJ5Npkm1Mt5utr0obHts2arqTTypttm3SJsP64f62+bj525PMLrPl2f3DfcP6zd/6TtukhkSbuz/br7c/JtufJ8mknfCY2bra2iZz1224rs31NSxXm62rSUtVy2RYntTW8qRaJpNj01/D5FgmQ7JsTI5lUm22bnP50LB8aHO5Ntc/lkmmP++c4fahyWPZ2LyulkP1WDZybHb/OZNHs5GWSR3LOZv31WPZGJYnOZZz67FMarjOsZxbj2aSNluebntseOyxbKTlnGrDdbKRyjlVmaRyTibZqMqhbGSjJplkko1LPv2XSV6a5LYkN2R6HESS3D4s/+5w/wf3wz+ogx1y9OFcc+gFQx5OzixHJxvD9WRYf/ocbZPJdP/cXH8mOTrLsdrKj13m6Ak5dQY5uv1nJLvI0UmS2mOODnlZw884kxzdGB67mZdnkqPz+TgZcuV0OXrOkIebeXomOXru8DznZuv55OgucnRjkppM5vLxDHJ0Y5JWlQz5t6sc3cy3jWmenlGObmxbPpMc3fxZm//mnUGO7pjjp8vRWX5mlou7zdHJbHnrQ+duc3Sj5vJzyNnd5ug52653k6Pz/5Zu5NgZ5ei5OZZJJeekZUOO7jpHNz//TvNy89/L3eXo/L+prWpXOTqfj5mcYY5uz6kzyNEd83C3ObqZj9see6oc3XzMZh6eSY7O52VVO6McPVTH/5t6Jjl63GfZ2WfTU+foOfVoNqpNr3PsjHL0nFQ2apqfk+H2qnJ0lGO9w0zyT2R6EO29Sd7VNyJYyNEkr62qI5nO2b91WP/WJBcP61+brTPIjZ4c5YCRozBuchTGbWk5OtbOaVpr70vyvs3lJ9ZF/0vHcGARj7TWnrN9ZWvt60le1iGepZCjHCByFMZNjsK4LS1HR9k5BQAA4OyiOAUAAKA7xSkAAADdKU4BAADoTnEKAABAd7Vfvg6qqu5J8vXecQyenOSLvYMYiGVnY4rl/Nbad/YOYtVGlqO7MaZ9ZDf2U7z7KdZEjo7VftqP9lOsyf6L92zJ0d/M9HczRvttn1kWr3t3vthau24ZTzzar5LZwddba1f3DiJJqupjYjmRWHZWVR/rHcOajCZHd2NM+8hu7Kd491OsiRwdq/20H+2nWJP9GW/vGNZhWR/uV2G/7TPL4nWvn7FeAAAAulOcAgAA0N1+Kk5v6R3AHLHsTCw7G1Msq7TfXqd4V2c/xZrsv3j3ar+9zv0U736KNREvZ+5s/R143Wu2b06IBAAAwMG1nzqnAAAAHFCKUwAAALobXXFaVW+rqgeH72Pb6f6qql+oqiNVdXdVPbtTHD88PP/dVfXvqupZq4hjN7HMbfe3q+qxqnppz1iq6nlV9Ymq+mRV/ZtesVTVk6rq/6mqPxhieeUKY7m8qj5UVfcOz/XqHbZZy767SlV1UVXdUVWfHq4Pn2S736yqr1TVb2xbf2VV3Tk8/p1Vde5I4r1h2ObTVXXD3PoPV9V9w/78iap66gpivG54jiNVdfMO9583vFdHhvfuirn7Xj+sv6+qXrjs2JYZb1VdUVV/M/de/vIIYv3+qvp4VT26/f+bJ9snxk6OylE5yiKq6sKqendV/fuafqb5vm33V+3zzzI72cXrXtvn/nU63eue227lNcZMa21UlyTfn+TZSe45yf0vTvKvk1SSa5Lc2SmOv5Pk8HD7RauKYzexDNtsJPlgkvcleWnH38+FSf4oydOH5ad2jOVnkvzccPspSb6c5NwVxXJJkmcPt5+Q5FNJnrltm7Xsu6u8JPlfk9w83L558/3dYbtrk/yDJL+xbf27krx8uP3LSf7b3vEmuSjJZ4brw8Ptzdz+cJKrVxjfRpI/TvKMJOcm+YMd9pt/lOSXh9svT/LO4fYzh+3PS3Ll8HM2Vvx+LhLvFaf6f1inWK9I8l1J3jH//81T7RNjv8jRUe3zcnTxWA9cjo79kuTtSf6b4fa5SS7cdv++/yyzx9e9ts/9Y3rdw/q11Bibl9F1Tltrv5NpEXEy1yd5R5v6SJILq+qSdcfRWvt3rbWHhsWPJLls2THsNpbBTyb5tSQPriqOXcbyXyV5T2vtc8P2K4tnF7G0JE+oqkry+GHbR1cUywOttY8Pt7+W5N4kl27bbC377opdn+n/yDJc/9BOG7XWPpDka/Prht/D85O8+3SPX6LdxPvCJHe01r485PQdSdb1RejPSXKktfaZ1tojSW7LNOZ586/h3UmuHd7L65Pc1lp7uLX2J0mODD9vrPGu22ljba19trV2d5Jj2x7bc59YlBxdLjm6Omdrjo5WVT0x0z/8vzVJWmuPtNa+sm2zg/BZ5ji7ed3r/Ny/Lrv8fSdrqjE2ja443YVLk9w/t3w0JxYB63Zjpn9F6qKqLk3yn2X6V+7evi3J4WHU6q6q+tGOsfwfSb4jyReS/GGSV7fWtv8Dt3TDiNT3JLlz211j3HfP1NNaaw8k04I8yZmM0F2c5Cuttc0/EKzj9e8m3tP9Xv7lMOL2T1bwAW43+8Rsm+G9+2qm72WP/WmReJPkyqr6/ar6N1X190cQ6yoe25scXS452jfWVTyWk3tGkj/PNKd+v6reUlUXbNvmIL73u3nd87p+7l+i077uHjXGoXU90RLt9A9Pt+/Dqar/ONOd9O/1iiHJzyd5XWvtsT5//DzOoSTfm+nI2Dcl+d2q+khr7VMdYnlhkk9k2gn4D5LcUVX/trX2F6t6wqp6fKZ/XXrNDs8zqn33ZKrqt5N88w53/eyiP3qHdQu//iXEe6q4fri19vmqekKmv9cfyXS8bFl2856cbJse+9Mi8T6Q6bj/l6rqe5P8elX9rRXm4yLvz6hzVY4mkaMnI0dZxKFMD5f6ydbanVX15kzH7f/J3DYH8b3fzetOMprP/cuym9e99hpjPxanR5NcPrd8WaadsbWrqu9K8pYkL2qtfalHDIOrk9w27DRPTvLiqnq0tfbrHWI5muSLrbW/SvJXVfU7SZ6V6TGY6/bKJG9s04H5I1X1J0n+wyS/t4onq6pzMv1w9CuttffssMlo9t1Taa39wMnuq6o/q6pLWmsPDGM8ZzLi8cVMx38ODX+tX8rrX0K8R5M8b275skyPY0tr7fPD9deq6v/KdAxtmR98d7NPbG5ztKoOJXlSpiPqPfanPcc75OHDSdJau6uq/jjTSYuPdYz1VI993rbHfngpUS2BHJWjpyBHWcTRJEdba5uTX+/OtFjZvs3oP8ucod287jF97l+W3bzutdcY+3Gs9/YkP1pT1yT56uZI0DpV1dOTvCfJj3TqCs601q5srV3RWrsi0x3rH3UqTJPkvUn+flUdqqrHJXlupsdf9vC5TDu4qaqnJfn2TE+asHTDKNlbk9zbWvvnJ9lsFPvugm5PsnlWxBsy/X3vyvDB50NJNs/0dkaP36PdxPv+JC+oqsM1PVPoC5K8f9iHn5zM/vDwkiSnPGP2Hnw0yVU1PUPquZmenOT2U7yGlyb54PBe3p7k5TU98+aVSa7Kiv7wsox4q+opVbWRJFX1jCHeleTjGcR6MjvuEyuKc9nk6HLJ0b6xnsx+ztHRaq39aZL7q+rbh1XXZnqSy3kH4bPMcXbzusf0uX9ZdvO6u9QYbQRnipq/JPnVTEdLvpFpRX9jkh9P8uPD/ZXkFzM9w9sfZkVn6dtFHG9J8lCmY6OfSPKxXu/Jtm1vzWrP1nvaWJL8D5nu3PdkOt7aa1/5liS/Newn9yT5r1cYy9/LdKzl7rl94sU99t1VXjI9LukDST49XF80rL86yVvmtvu3mR7H8DfD7+aFw/pnZPrh7EiSf5XkvJHE+2NDTEeSvHJYd0GSu4bf6SeTvDkrONPmsJ98atgvfnZY90+T/OBw+/zhvToyvHfPmHvszw6Puy/Tv+SuYx/YU7xJ/vPhffyDJB9P8g9GEOvfHvbPv0rypSSfPNU+sR8uclSOylGXBX8n351pt/zuJL+e6dmQD9RnmT2+7rV97h/T69627a1Zw9l6a3gyAAAA6GY/jvUCAABwwChOAQAA6E5xCgAAQHeKUwAAALpTnAIAANCd4hQAAIDuFKcAAAB09/8DG2ArreKAbFMAAAAASUVORK5CYII=", "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": [ "Writing output to c:\\Users\\anbro\\Documents\\projects\\GPUSimulators\\data\\euler_rayleigh-taylor_0002.nc\n", "Opening c:\\Users\\anbro\\Documents\\projects\\GPUSimulators\\data\\euler_rayleigh-taylor_0002.nc\n", "Keyword arguments: {'mode': 'w', 'clobber': False}\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Constructed in 0.01562356948852539 seconds\n", "Simulating to t=30.000000. Estimated 30697 timesteps (dt=0.000977)\n", "0% [####==========================] 100%. Total: 1m 26s, elapsed: 10s, remaining: 1m 16s\n", "0% [#######=======================] 100%. Total: 1m 26s, elapsed: 20s, remaining: 1m 6s\n", "0% [##########====================] 100%. Total: 1m 26s, elapsed: 30s, remaining: 56s\n", "0% [##############================] 100%. Total: 1m 27s, elapsed: 40s, remaining: 46s\n", "0% [#################=============] 100%. Total: 1m 27s, elapsed: 50s, remaining: 37s\n", "0% [#####################=========] 100%. Total: 1m 27s, elapsed: 1m, remaining: 27s\n", "0% [########################======] 100%. Total: 1m 28s, elapsed: 1m 10s, remaining: 17s\n", "0% [############################==] 100%. Total: 1m 28s, elapsed: 1m 20s, remaining: 7s\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "Closing c:\\Users\\anbro\\Documents\\projects\\GPUSimulators\\data\\euler_rayleigh-taylor_0002.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": [ "Writing output to c:\\Users\\anbro\\Documents\\projects\\GPUSimulators\\data\\euler_rayleigh-taylor_0002.nc\n", "Opening c:\\Users\\anbro\\Documents\\projects\\GPUSimulators\\data\\euler_rayleigh-taylor_0002.nc\n", "Arguments: ('r',)\n", "Animation.save using \n", "figure size (inches) has been adjusted from 3.4 x 8.0 to 3.388888888888889 x 8.0\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "c:\\Users\\anbro\\Documents\\projects\\GPUSimulators\\data\\euler_rayleigh-taylor_0002.nc created Fri Nov 9 10:31:24 2018 contains 301 timesteps\n", "Simulator arguments: \n", " {'rho': '[...]', 'rho_u': '[...]', 'rho_v': '[...]', 'E': '[...]', 'nx': 151, 'ny': 453, 'dx': 0.0033112582781456954, 'dy': 0.0033112582781456954, 'dt': 0.0009772880151828723, 'g': 0.1, 'gamma': 1.4, 'boundary_conditions': '[north=Type.Reflective, south=Type.Reflective, east=Type.Reflective, west=Type.Reflective]', 'context': 'CudaContext id 346033291232'}\n", "0% [########======================] 100%. Total: 19s, elapsed: 5s, remaining: 14s\n", "0% [##############================] 100%. Total: 21s, elapsed: 10s, remaining: 11s\n", "0% [#####################=========] 100%. Total: 22s, elapsed: 15s, remaining: 7s\n", "0% [###########################===] 100%. Total: 22s, elapsed: 20s, remaining: 2s\n" ] }, { "data": { "text/html": [ "\n", "\n", "\n", "\n", "
\n", " \n", "
\n", " \n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", " Once \n", " Loop \n", " Reflect \n", "
\n", "
\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Closing c:\\Users\\anbro\\Documents\\projects\\GPUSimulators\\data\\euler_rayleigh-taylor_0002.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": { "interpreter": { "hash": "d80e56d67bdb125526bdf12740b058f9c8b2e6eb30981cd0c9aaae49693d1172" }, "kernelspec": { "display_name": "Python [conda env:anaconda3]", "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.7.12" } }, "nbformat": 4, "nbformat_minor": 2 }