From bdf7b4292c97a592ec749aaa770a3a37690d284d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20R=2E=20Brodtkorb?= Date: Thu, 8 Nov 2018 22:07:31 +0100 Subject: [PATCH] Implemented proper boundary conditions handling --- EulerTesting.ipynb | 401433 +++++++++++++++++++- GPUSimulators/EE2D_KP07_dimsplit.py | 10 +- GPUSimulators/Simulator.py | 9 + GPUSimulators/cuda/EE2D_KP07_dimsplit.cu | 44 +- GPUSimulators/cuda/common.h | 189 +- 5 files changed, 401229 insertions(+), 456 deletions(-) diff --git a/EulerTesting.ipynb b/EulerTesting.ipynb index 1c43390..2a5af7c 100644 --- a/EulerTesting.ipynb +++ b/EulerTesting.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ @@ -17,7 +17,6 @@ "from matplotlib import animation, rc\n", "from matplotlib import pyplot as plt\n", "from matplotlib.colorbar import ColorbarBase\n", - "from matplotlib.colors import Normalize\n", "\n", "from IPython.display import display, HTML\n", "\n", @@ -38,30 +37,56 @@ "except ImportError:\n", " from io import StringIO\n", "\n", - "from GPUSimulators import Common, IPythonMagic" + "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": null, + "execution_count": 2, "metadata": {}, - "outputs": [], + "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.6.7 |Anaconda custom (64-bit)| (default, Oct 28 2018, 19:44:12) [MSC v.1915 64 bit (AMD64)]\n" + ] + } + ], "source": [ "%setup_logging --out test_schemes.log" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Registering my_context in user workspace\n", + "PyCUDA version 2018.1.1\n", + "CUDA version (10, 0, 0)\n", + "Driver version 10000\n", + "Using 'GeForce 840M' GPU\n", + "Created context handle <716326090688>\n", + "Using CUDA cache dir c:\\Users\\anbro\\Documents\\projects\\GPUSimulators\\GPUSimulators\\cuda_cache\n" + ] + } + ], "source": [ "%cuda_context_handler --no_autotuning my_context " ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ @@ -73,263 +98,66 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ - "def genShockBubble(nx, ny, gamma):\n", - " width = 4.0\n", - " height = 1.0\n", - " dx = width / float(nx)\n", - " dy = height / float(ny)\n", + "def plotVars(rho, rho_u, rho_v, E):\n", + " plt.figure()\n", + " plt.subplot(1,4,1)\n", + " plt.imshow(rho, origin='bottom')\n", + " plt.colorbar(orientation='horizontal', pad=0.02, shrink=0.9)\n", "\n", - " x_center = 0.5 #dx*nx/2.0\n", - " y_center = 0.5 #dy*ny/2.0\n", + " plt.subplot(1,4,2)\n", + " plt.imshow(rho_u, origin='bottom')\n", + " plt.colorbar(orientation='horizontal', pad=0.02, shrink=0.9)\n", "\n", + " plt.subplot(1,4,3)\n", + " plt.imshow(rho_v, origin='bottom')\n", + " plt.colorbar(orientation='horizontal', pad=0.02, shrink=0.9)\n", "\n", - " rho = np.ones((ny, nx), dtype=np.float32)\n", - " u = np.zeros((ny, nx), dtype=np.float32)\n", - " v = np.zeros((ny, nx), dtype=np.float32)\n", - " E = np.zeros((ny, nx), dtype=np.float32)\n", - " p = np.ones((ny, nx), dtype=np.float32)\n", - "\n", - " x = (dx*(np.arange(0, nx, dtype=np.float32)+0.5) - x_center)\n", - " y = (dy*(np.arange(0, ny, dtype=np.float32)+0.5) - y_center)\n", - " xv, yv = np.meshgrid(x, y, sparse=False, indexing='xy')\n", - " \n", - " #Bubble\n", - " radius = 0.25\n", - " bubble = np.sqrt(xv**2+yv**2) <= radius\n", - " rho = np.where(bubble, 0.1, rho)\n", - " \n", - " #Left boundary\n", - " left = (xv - xv.min() < 0.1)\n", - " rho = np.where(left, 3.81250, rho)\n", - " u = np.where(left, 2.57669, u)\n", - " \n", - " #Energy\n", - " p = np.where(left, 10.0, p)\n", - " E = 0.5*rho*(u**2+v**2) + p/(gamma-1.0)\n", - " \n", - " #Estimate dt\n", - " scale = 0.9\n", - " max_rho_estimate = 5.0\n", - " max_u_estimate = 5.0\n", - " dx = width/nx\n", - " dy = height/ny\n", - " dt = scale * min(dx, dy) / (max_u_estimate + np.sqrt(gamma*max_rho_estimate))\n", - "\n", - " return rho, rho*u, rho*v, E, dx, dy, dt\n", - "\n", - "nx = 400\n", - "rho0, rho_u0, rho_v0, E0, dx, dy, dt = genShockBubble(nx, nx//4, 1.4)\n", - "plt.figure()\n", - "plt.imshow(rho0)\n", - "plt.colorbar(orientation='vertical', aspect=20, pad=0.02, shrink=0.3)" + " plt.subplot(1,4,4)\n", + " plt.imshow(E, origin='bottom')\n", + " plt.colorbar(orientation='horizontal', pad=0.02, shrink=0.9)" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 6, "metadata": {}, "outputs": [], "source": [ - "def genKelvinHelmholtz(nx, ny, gamma, roughness=0.5):\n", - " \"\"\"\n", - " Roughness parameter in (0, 1.0] determines how squiggly the interface betweeen the zones is\n", - " \"\"\"\n", - " \n", - " def genZones(nx, ny, n):\n", - " \"\"\"\n", - " Generates the zones of the two fluids of K-H\n", - " \"\"\"\n", - " zone = np.zeros((ny, nx), dtype=np.int32)\n", - "\n", - " dx = 1.0 / nx\n", - " dy = 1.0 / ny\n", - "\n", - " def genSmoothRandom(nx, n):\n", - " assert (n <= nx), \"Number of generated points nx must be larger than n\"\n", - " \n", - " if n == nx:\n", - " return np.random.random(nx)-0.5\n", - " else:\n", - " from scipy.interpolate import interp1d\n", - "\n", - " #Control points and interpolator\n", - " xp = np.linspace(0.0, 1.0, n)\n", - " yp = np.random.random(n) - 0.5\n", - " f = interp1d(xp, yp, kind='cubic')\n", - "\n", - " #Interpolation points\n", - " x = np.linspace(0.0, 1.0, nx)\n", - " return f(x)\n", - "\n", - "\n", - "\n", - " x = np.linspace(0, 1, nx)\n", - " y = np.linspace(0, 1, ny)\n", - "\n", - " _, y = np.meshgrid(x, y)\n", - "\n", - " #print(y+a[0])\n", - "\n", - " a = genSmoothRandom(nx, n)*dy\n", - " zone = np.where(y > 0.25+a, zone, 1)\n", - "\n", - " a = genSmoothRandom(nx, n)*dy\n", - " zone = np.where(y < 0.75+a, zone, 1)\n", - " \n", - " return zone\n", - " \n", - " width = 2.0\n", - " height = 1.0\n", - " dx = width / float(nx)\n", - " dy = height / float(ny)\n", - "\n", - " rho = np.empty((ny, nx), dtype=np.float32)\n", - " u = np.empty((ny, nx), dtype=np.float32)\n", - " v = np.zeros((ny, nx), dtype=np.float32)\n", - " p = 2.5*np.ones((ny, nx), dtype=np.float32)\n", - "\n", - " #Generate the different zones \n", - " zones = genZones(nx, ny, max(1, min(nx, int(nx*roughness))))\n", - " \n", - " #Zone 0\n", - " zone0 = zones == 0\n", - " rho = np.where(zone0, 1.0, rho)\n", - " u = np.where(zone0, 0.5, u)\n", - " \n", - " #Zone 1\n", - " zone1 = zones == 1\n", - " rho = np.where(zone1, 2.0, rho)\n", - " u = np.where(zone1, -0.5, u)\n", - " \n", - " E = 0.5*rho*(u**2+v**2) + p/(gamma-1.0)\n", - " \n", - " #Estimate dt\n", - " scale = 0.9\n", - " max_rho_estimate = 3.0\n", - " max_u_estimate = 2.0\n", - " dx = width/nx\n", - " dy = height/ny\n", - " dt = scale * min(dx, dy) / (max_u_estimate + np.sqrt(gamma*max_rho_estimate))\n", - "\n", - " return rho, rho*u, rho*v, E, dx, dy, dt\n", - "\n", - "nx = 400\n", - "rho0, rho_u0, rho_v0, E0, dx, dy, dt = genKelvinHelmholtz(nx, nx//4, 1.4, 0.25)\n", - "plt.figure()\n", - "plt.imshow(rho0)\n", - "plt.colorbar(orientation='vertical', aspect=20, pad=0.02, shrink=0.3)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "def genRayleighTaylor(nx, ny, gamma):\n", - " width = 0.5\n", - " height = 1.5\n", - " dx = width / float(nx)\n", - " dy = height / float(ny)\n", - " g = 0.1\n", - "\n", - " rho = np.zeros((ny, nx), dtype=np.float32)\n", - " u = np.zeros((ny, nx), dtype=np.float32)\n", - " v = np.zeros((ny, nx), dtype=np.float32)\n", - " p = np.zeros((ny, nx), dtype=np.float32)\n", - " \n", - " x = dx*np.arange(0, nx, dtype=np.float32)-width*0.5\n", - " y = dy*np.arange(0, ny, dtype=np.float32)-height*0.5\n", - " xv, yv = np.meshgrid(x, y, sparse=False, indexing='xy')\n", - " \n", - " #y_threshold = -0.01*np.cos(2*np.pi*np.abs(x)/0.5)\n", - " #rho = np.where(yv <= y_threshold, 1.0, rho)\n", - " #rho = np.where(yv > y_threshold, 2.0, rho)\n", - " \n", - " rho = np.where(yv <= 0.0, 1.0, rho)\n", - " rho = np.where(yv > 0.0, 2.0, rho)\n", - " \n", - " v = 0.01*(1.0 + np.cos(2*np.pi*xv/0.5))/4\n", - " \n", - " p = 2.5 - rho*g*yv\n", - " E = 0.5*rho*(u**2+v**2) + p/(gamma-1.0)\n", - " \n", - " #Estimate dt\n", - " scale = 0.9\n", - " max_rho_estimate = 3.0\n", - " max_u_estimate = 1.0\n", - " dx = width/nx\n", - " dy = height/ny\n", - " dt = scale * min(dx, dy) / (max_u_estimate + np.sqrt(gamma*max_rho_estimate))\n", - "\n", - " return rho, rho*u, rho*v, E, dx, dy, dt\n", - "\n", - "nx = 400\n", - "rho0, rho_u0, rho_v0, E0, dx, dy, dt = genRayleighTaylor(nx, nx*3, 1.4)\n", - "plt.figure()\n", - "plt.subplot(1,4,1)\n", - "plt.imshow(rho0, origin='bottom')\n", - "plt.colorbar(orientation='horizontal', aspect=20, pad=0.02, shrink=0.9)\n", - "\n", - "plt.subplot(1,4,2)\n", - "plt.imshow(rho_u0, origin='bottom')\n", - "plt.colorbar(orientation='horizontal', aspect=20, pad=0.02, shrink=0.9)\n", - "\n", - "plt.subplot(1,4,3)\n", - "plt.imshow(rho_v0, origin='bottom')\n", - "plt.colorbar(orientation='horizontal', aspect=20, pad=0.02, shrink=0.9)\n", - "\n", - "plt.subplot(1,4,4)\n", - "plt.imshow(E0, origin='bottom')\n", - "plt.colorbar(orientation='horizontal', aspect=20, pad=0.02, shrink=0.9)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "def makeSimulation(outfile, t_end, rho0, rho_u0, rho_v0, E0, nx, ny, dx, dy, dt, gamma):\n", - " #Construct simulator\n", - " arguments = {\n", - " 'context': my_context,\n", - " 'rho': rho0, 'rho_u': rho_u0, 'rho_v': rho_v0, 'E': E0,\n", - " 'nx': nx, 'ny': ny,\n", - " 'dx': dx, 'dy': dy, 'dt': 0.45*dt,\n", - " 'gamma': gamma\n", - " } \n", + "def runSimulation(outfile, t_end, sim_args):\n", " with Common.Timer(\"construct\") as t:\n", - " from GPUSimulators import EE2D_KP07_dimsplit\n", - " importlib.reload(EE2D_KP07_dimsplit)\n", - " sim = EE2D_KP07_dimsplit.EE2D_KP07_dimsplit(**arguments)\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', nx)\n", - " outdata.ncfile.createDimension('y', ny)\n", + " outdata.ncfile.createDimension('x', sim.nx)\n", + " outdata.ncfile.createDimension('y', sim.ny)\n", "\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)\n", - " outdata.rho_u_var = outdata.ncfile.createVariable('rho_u', np.dtype('float32').char, ('time', 'y', 'x'), zlib=True)\n", - " outdata.rho_v_var = outdata.ncfile.createVariable('rho_v', np.dtype('float32').char, ('time', 'y', 'x'), zlib=True)\n", - " outdata.E_var = outdata.ncfile.createVariable('E', np.dtype('float32').char, ('time', 'y', 'x'), zlib=True)\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", - " outdata.x_var[:] = np.linspace(0, nx*dx, nx)\n", - " outdata.y_var[:] = np.linspace(0, ny*dy, ny)\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", - " sim.check()\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", @@ -354,46 +182,7 @@ }, { "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "def genSchlieren(rho):\n", - " #Compute length of z-component of normalized gradient vector \n", - " normal = np.gradient(rho) #[x, y, 1]\n", - " length = 1.0 / np.sqrt(normal[0]**2 + normal[1]**2 + 1.0**2)\n", - " schlieren = np.power(length, 128)\n", - " return schlieren\n", - "\n", - "\n", - "def genCurl(rho, rho_u, rho_v):\n", - " u = rho_u / rho\n", - " v = rho_v / rho\n", - " u = np.sqrt(u**2 + v**2)\n", - " u_max = u.max()\n", - " \n", - " du_dy, du_dx = np.gradient(u)\n", - " dv_dy, dv_dx = np.gradient(v)\n", - " \n", - " curl = dv_dx - du_dy\n", - " return curl\n", - "\n", - "\n", - "def genColors(rho, rho_u, rho_v, cmap, vmax, vmin):\n", - " schlieren = genSchlieren(rho)\n", - " curl = genCurl(rho, rho_u, rho_v)\n", - " \n", - " colors = Normalize(vmin, vmax, clip=True)(curl)\n", - " colors = cmap(colors)\n", - " for k in range(3):\n", - " colors[:,:,k] = colors[:,:,k]*schlieren\n", - " \n", - " return colors" - ] - }, - { - "cell_type": "code", - "execution_count": null, + "execution_count": 7, "metadata": {}, "outputs": [], "source": [ @@ -401,7 +190,7 @@ " Schlieren = 0\n", " Density = 1\n", "\n", - "def makeAnimation(infile, vis_type, vmin, vmax, save_anim=False, cmap=plt.cm.coolwarm, fig_args=None):\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", @@ -417,7 +206,7 @@ " extents = [0, x.max(), 0, y.max()]\n", " \n", " if (vis_type == VisType.Schlieren):\n", - " im = ax1.imshow(genColors(rho, rho_u, rho_v, cmap, vmax, vmin), origin='bottom', extent=extents, cmap='gray', vmin=0.0, vmax=1.0)\n", + " im = ax1.imshow(Visualization.genColors(rho, rho_u, rho_v, cmap, vmax, vmin), origin='bottom', 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='bottom', extent=extents, cmap=cmap, vmin=vmin, vmax=vmax)\n", @@ -426,6 +215,7 @@ " 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", @@ -447,7 +237,7 @@ " rho_v = indata.ncfile.variables['rho_v'][i]\n", " \n", " if (vis_type == VisType.Schlieren):\n", - " im.set_data(genColors(rho, rho_u, rho_v, cmap, vmax, vmin))\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", @@ -490,6 +280,236882 @@ "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": "iVBORw0KGgoAAAANSUhEUgAAAzkAAADhCAYAAADrl3vGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAGFVJREFUeJzt3XGsZOV53/Hvj8XECXaC7V1bFHCyTdbF1GrAvsW0SCkB0i60YqmEK0jrEIt2WwkSJzipcVrhlDZSnLahjUTcbgyBRI4xIU5ZWVsTi4CiVoXsrqG2dwlis0awZuvdjYGkcTHZ7NM/5lyYHWbuzNw7d+7Mud+PdHTnnDnznndeXmbPM+/zvpOqQpIkSZLa4pS1roAkSZIkTZJBjiRJkqRWMciRJEmS1CoGOZIkSZJaxSBHkiRJUqsY5EiSJElqFYMcSZIkSWsmyV1JjiT56oDnk+RXkhxI8uUk7x1WpkGOJEmSpLV0N7B1ieevALY023bgk8MKNMiRJEmStGaq6g+Aby5xyjbgN6rjUeCMJGcuVeapk6ygJEmSpPVh69atdezYsSXP2bt37z7g5a5DO6pqx5iXOgt4rmv/UHPs8KAXGORIkiRJGtuxY8fYs2fPkuckebmqFlZ4qfQ5Vku9wCBHkiRJ0rKcOHFiGpc5BJzTtX828PxSL1jRnJwkW5M81ax0cMtKypIkSZI0X6pqyW1CdgI/1qyydhHwUlUNTFWDFYzkJNkA3AH8CJ3oaneSnVW1f7llSpIkSZoPkwpkknwGuATYmOQQ8HHgDc01/guwC7gSOAB8C/jQsDJXkq52IXCgqg42lbuXzsoHBjmSJEnSOjCJdLWqum7I8wXcOE6ZKwly+q1y8P7ek5Jsp7OeNaeffvr7zj333BVc8jV79+599fH73ve+iZQpSZIktd0zzzzDsWPH+k3mH9sEU9ImaiVBzkirHDRLxO0AWFhYqGErMIx88bx2+UmVKUmSJLXdwsJKFzt7TRuDnLFXOZAkSZLUDlU1rdXVxraS1dV2A1uSbE5yGnAtnZUPJEmSJK0DU1pdbWzLHsmpquNJbgIeBDYAd1XVvonVTJIkSdJMa2O6GlW1i86SbpIkSZLWmVYGOZIkSZLWp1mek2OQI0mSJGlZHMmRJEmS1CoGOZIkSZJaw3Q1SZIkSa3jSI4kSZKkVjHIkSRJktQapqtJkiRJah1HciRJkiS1ikGOJEmSpFYxXU2SJElSa1SVIzmSJEmS2sUgR5IkSVKrGORIkiRJahXn5EiSJElqDefkSJIkSWodgxxJkiRJrWK6miRJkqRWcSRHkiRJUms4J0eSJElS65iuJkmSJKlVHMmRJEmS1CoGOZIkSZJao6pMV5MkSZLULrM6knPKWldAkiRJ0nxaXGFt0DaKJFuTPJXkQJJb+jz/ziQPJ3k8yZeTXDmsTIMcSZIkScuy0iAnyQbgDuAK4DzguiTn9Zz2r4H7quoC4FrgV4eVa7qaJEmSpLFNaE7OhcCBqjoIkOReYBuwv/tSwHc3j78HeH5YoQY5kiRJkpZlhNGajUn2dO3vqKodXftnAc917R8C3t9Txs8Dv5fkJ4DTgcuHXXRoulqSc5ocuCeT7Evy4eb4W5N8McnTzd+3DCtLkiRJUnuMkK52rKoWurYdPUWkX7E9+9cBd1fV2cCVwG8mWTKOGWVOznHgI1X1buAi4MYmT+4W4KGq2gI81OxLkiRJWgcW09WW2kZwCDina/9sXp+OdgNwX3PN/wW8Edi4VKFDg5yqOlxVX2oe/xnwJJ1hpW3APc1p9wBXD30LkiRJklpjAqur7Qa2JNmc5DQ6Cwvs7DnnWeAygCTvphPkHF2q0LHm5CT5PuAC4DHgHVV1GDqBUJK3j1OWJEmSpPm20t/JqarjSW4CHgQ2AHdV1b4ktwF7qmon8BHg15L8NJ1Uth+vIRceOchJ8ibgd4Cfqqo/Tfqlz/V93XZgO8A73/nOUS8nSZIkacZNYHU1qmoXsKvn2K1dj/cDF49T5ki/k5PkDXQCnE9X1eeaw99Icmbz/JnAkQGV3rE40WjTpk3j1E2SJEnSjBqWqrbSUZ6VGGV1tQB3Ak9W1S93PbUTuL55fD3wwOSrJ0mSJGlWzWqQM0q62sXAB4GvJHmiOfZzwC8C9yW5gc5koA+sThUlSZIkzaJJpKuthqFBTlX9D/qvXw3NKgeSJEmS1p+1HK1Zylirq0mSJEkSsOYpaUsxyJlTo65uN4pZ7ZySJEmabbN6H2mQM8MmGcgs9zqz2nElSZK09uZ2To4kSZIk9TOrX4gb5MyQcUZuJtmhlrpu73Oz2pElSZI0Xc7JkSRJktQ6pqupr1mYD7PUdXrr58iOJEmSFs3qvaBBjiRJkqRlMcjRSfqN4MxiJ+mt06CRnVmsuyRJklZPVZmuJkmSJKldZvWLboOcKRk092ZWO8Ygg0Z2ut/fvL0nSZIkLc+s3vcZ5EiSJEkam+lq61hbRnAGWXwf3e/TeTqSJEnrw6ze7xnkSJIkSVoWg5x1Zr39nkz3++udp9P29y5JkrRezep9nkGOJEmSpLE5J2edWe8rjfXO03FER5IkqZ1m9f7OIGeC1ntw08tgR5Ikqd1m9b7OIEeSJEnSspiu1mLrbZGBcTmiI0mS1D5VNbP3cwY5kiRJkpbFIEeSJElSq5iuJkmSJKlVHMlpIefijMe5OZIkSe3hnBxJkiRJrWO6Wos4grMy/UZ0bENJkqT5M6v3cKesdQUkSZIkzZ/FdLWltlEk2ZrkqSQHktwy4Jx/lGR/kn1JfmtYmY7kSJIkSVqWlY7kJNkA3AH8CHAI2J1kZ1Xt7zpnC/Ax4OKqeiHJ24eVO/JITpINSR5P8vlmf3OSx5I8neSzSU4b901JkiRJml8nTpxYchvBhcCBqjpYVa8A9wLbes75Z8AdVfUCQFUdGVboOOlqHwae7Nr/BHB7VW0BXgBuGKOsuZTkpPk4s7yixDzobr/etpUkSdLsm0C62lnAc137h5pj3d4FvCvJ/0zyaJKtwwodKchJcjbw94FPNfsBLgXub065B7h6lLIkSZIkzb8R5+RsTLKna9veU0y/b7l7o6NTgS3AJcB1wKeSnLFU3Uadk/OfgH8JvLnZfxvwYlUdb/b7RVySJEmSWmyElLRjVbWwxPOHgHO69s8Gnu9zzqNV9RfA15I8RSfo2T2o0KEjOUn+AXCkqvZ2H+5zat/xqCTbFyO3o0ePDrvcTDJNbXV1t6Vpa5IkSfNjAulqu4EtzXz/04BrgZ095/w34IcBkmykk752cKlCR0lXuxi4KskzdCYCXUpnZOeMJIsjQf0iLgCqakdVLVTVwqZNm0a4nCRJkqR5sNIgp8kMuwl4kM78//uqal+S25Jc1Zz2IPAnSfYDDwM/W1V/slS5Q9PVqupjdJZsI8klwM9U1T9O8tvANXQCn+uBB4a+C0mSJEmtUFWjrqA2rJxdwK6eY7d2PS7g5mYbyUp+DPSjwM1JDtCZo3PnCsqSJEmSNGcm8WOgq2GsHwOtqkeAR5rHB+msay1JkiRpHZrVeepjBTmSJEmSBJNLV1sNBjmSJEmSlsWRHEmSJEmtYpAzp3p/H0erY7FtF9t78a9tLkmSNLtMV5MkSZLUGmu9gtpSDHIkSZIkLYtBjiRJkqRWMciRJEmS1CrOyZEkSZLUGs7JkSRJktQ6BjmSJEmSWsV0NUmSJEmt4kiOJEmSpNZwTo4kSZKk1jFdTZIkSVKrOJIjSZIkqVUMciRJkiS1RlWZriZJkiSpXRzJkSRJktQqBjmSJEmSWsUgR5IkSVJrOCdnjlUVSQBe/TurEes8W2zbRbaxJEnS7JvVezaDHEmSJEnLYpAjSZIkqTVMV5MkSZLUOo7kSJIkSWoVgxxJkiRJrWK6miRJkqTWqKqZHck5ZZSTkpyR5P4kf5TkySR/K8lbk3wxydPN37esdmUlSZIkzY7FQGfQtlZGCnKA/wx8oarOBX4QeBK4BXioqrYADzX7rdT7HynJ637XRcvX3ZZr/T+EJEmSRnfixIklt1Ek2ZrkqSQHkgyMKZJck6SSLAwrc2iQk+S7gR8C7gSoqleq6kVgG3BPc9o9wNWjvAlJkiRJ7bDSkZwkG4A7gCuA84DrkpzX57w3Az8JPDZKvUYZyfmrwFHg15M8nuRTSU4H3lFVh5s3dxh4+ygXlCRJkjT/hgU4I2bnXAgcqKqDVfUKcC+dwZRe/xb4JeDlUQodJcg5FXgv8MmqugD4c8ZITUuyPcmeJHuOHj066stmkmlrk9XdfqapSZIkzZ8RgpyNi7FAs23vKeIs4Lmu/UPNsVcluQA4p6o+P2q9Rlld7RBwqKoWh4bupxPkfCPJmVV1OMmZwJF+L66qHcAOgIWFBe9iJUmSpJYYYd7Nsapaag5NvxGDV2OGJKcAtwM/Pk69ho7kVNX/AZ5L8teaQ5cB+4GdwPXNseuBB8a5sCRJkqT5NoF0tUPAOV37ZwPPd+2/GXgP8EiSZ4CLgJ3DFh8Y9XdyfgL4dJLTgIPAh+gESPcluQF4FvjAiGVJkiRJmnMTmm6wG9iSZDPwdeBa4Ee7rvESsHFxP8kjwM9U1Z6lCh0pyKmqJ4B+0dJlo7y+bRb/Yy7OJ+meV6Lheucx2W6SJEnzadRlogepquNJbgIeBDYAd1XVviS3AXuqaudyyh11JEeSJEmSTjKJL6urahewq+fYrQPOvWSUMg1yVsARnfE4giNJktQus3o/Z5AjSZIkaWxVteJ0tdVikCNJkiRpWRzJkSRJktQqBjkt1m9uzqz+B18LzsWRJElqH9PVJEmSJLXOrH55bZAzQd0jOq605giOJElS283q/Z1BziqoqnW9rLTBjSRJ0vowq/d5BjmSJEmSxuacnHVo0A+F9j7fFr3vD9r3HiVJknSyWb3fM8iRJEmStCwGOetU74jOorbM1XEER5IkaX0yXU2SJElS68zql9sGOVPS2wH6zdWZ1U7Srd/IDcxH3SVJkjRZs3oPaJAjSZIkaVlMV9NJ+s3VmcUV2AaN3CyahTpKkiRp+qpqZu8FDXIkSZIkLYtBjvrq7hiDVmAb9rqVGjZas1rXlSRJ0nwzXU2SJElSq8zqF+AGOTNk0Aps/Ywz+rISs9pxJUmStLack6NlWarTTDLImdXOKUmSpNk2q/eRBjmSJEmSlsU5OZqoWY2aJUmStH7M6j2pQY4kSZKksTknR5IkSVLrmK4mSZIkqVUcyZEkSZLUKgY5kiRJklqjqmY2Xe2UUU5K8tNJ9iX5apLPJHljks1JHkvydJLPJjlttSsrSZIkaXYsLj4waFsrQ4OcJGcBPwksVNV7gA3AtcAngNuragvwAnDDalZUkiRJ0myZ2yCncSrwnUlOBb4LOAxcCtzfPH8PcPXkqydJkiRpFi2mqy21jSLJ1iRPJTmQ5JY+z9+cZH+SLyd5KMn3DitzaJBTVV8H/gPwLJ3g5iVgL/BiVR1vTjsEnDWg0tuT7Emy5+jRo8MuJ0mSJGlOrHQkJ8kG4A7gCuA84Lok5/Wc9jidrLK/QWeQ5ZeGlTtKutpbgG3AZuCvAKc3lejV911U1Y6qWqiqhU2bNg27nCRJkqQ5MYF0tQuBA1V1sKpeAe6lE3t0X+PhqvpWs/socPawQkdJV7sc+FpVHa2qvwA+B/xt4IwmfY3mQs+P8i4kSZIktcMIQc7GxayuZtveU8RZwHNd+wMzxBo3AP99WL1GWUL6WeCiJN8F/D/gMmAP8DBwDZ1o63rggRHKkiRJktQCIy4hfayqFpZ4Pv2K7nti8k+ABeDvDLvoKHNyHqOT+/Yl4CvNa3YAHwVuTnIAeBtw57CyJEmSJLXHBNLVDgHndO33zRBLcjnwr4Crqurbwwod6cdAq+rjwMd7Dh+kk0MnSZIkaR2awDLRu4EtSTYDX6fzUzU/2n1CkguA/wpsraojoxQ6UpAjSZIkSb1GXSZ6kKo6nuQm4EE6v8d5V1XtS3IbsKeqdgL/HngT8NtJAJ6tqquWKtcgR5IkSdLYJvWDn1W1C9jVc+zWrseXj1umQY4kSZKkZZlEkLMaDHIkSZIkLctK09VWi0GOJEmSpGVxJEeSJElSa0xqTs5qMMiRJEmStCymq0mSJElqFUdyJEmSJLWKQY4kSZKk1nBOjiRJkqTWcU6OJEmSpFZxJEeSJElSqxjkSJIkSWqNqjJdTZIkSVK7OJIjSZIkqVUMciRJkiS1hulqkiRJklrHkRxJkiRJrWKQI0mSJKlVTFeTJEmS1BpV5UiOJEmSpHYxyJEkSZLUKgY5kiRJklrFOTmSJEmSWsM5OZIkSZJaxyBHkiRJUquYriZJkiSpVRzJkSRJktQazsmRJEmS1DqmqwF79+49luTPgWOTLDfJJItrm41MuL01lG0+fbb5dNne02ebT59tPl2293R976QKciQHqKpNSfZU1cI0r7ue2d7TZ5tPn20+Xbb39Nnm02ebT5ftPb9mNcg5Za0rIEmSJGn+VBUnTpxYchtFkq1JnkpyIMktfZ7/jiSfbZ5/LMn3DSvTIEeSJEnSsiwuPjBoGybJBuAO4ArgPOC6JOf1nHYD8EJV/QBwO/CJYeWuRZCzYw2uuZ7Z3tNnm0+fbT5dtvf02ebTZ5tPl+09p1Ya5AAXAgeq6mBVvQLcC2zrOWcbcE/z+H7gsgyZlJ9ZzaOTJEmSNLuSfIHOohFLeSPwctf+jqp6NahNcg2wtar+abP/QeD9VXVT1zlfbc451Oz/cXPOwMUqXEJakiRJ0tiqausEiuk3ItM7CjPKOSdxTo4kSZKktXIIOKdr/2zg+UHnJDkV+B7gm0sVOrUgZ9iqCZqMJM8k+UqSJ5LsaY69NckXkzzd/H3LWtdzniW5K8mRZuh08VjfNk7HrzT9/stJ3rt2NZ9PA9r755N8vennTyS5suu5jzXt/VSSv7c2tZ5vSc5J8nCSJ5PsS/Lh5rj9fBUs0d7281WS5I1J/jDJ/27a/N80xzc3Kzc93azkdFpzfOyVnfSaJdr77iRf6+rj5zfH/UxZX3YDW5r//04DrgV29pyzE7i+eXwN8Ps1ZM7NVIKcEVdN0OT8cFWd37Xe/C3AQ1W1BXio2dfy3Q30Ds8OauMrgC3Nth345JTq2CZ38/r2Bri96efnV9UugOZz5Vrgrzev+dXm80fjOQ58pKreDVwE3Ni0rf18dQxqb7Cfr5ZvA5dW1Q8C5wNbk1xEZ8Wm25s+/gKdFZ1gGSs76SSD2hvgZ7v6+BPNMT9T1pGqOg7cBDwIPAncV1X7ktyW5KrmtDuBtyU5ANzMCPey0xrJGWXVBK2e7hUp7gGuXsO6zL2q+gNeP0Q6qI23Ab9RHY8CZyQ5czo1bYcB7T3INuDeqvp2VX0NOEDn80djqKrDVfWl5vGf0flH5yzs56tiifYexH6+Qk1f/b/N7huarYBL6azcBK/v42Ot7KTXLNHeg/iZss5U1a6qeldVfX9V/UJz7Naq2tk8frmqPlBVP1BVF1bVwWFlTivIOQt4rmv/EEt/gGv5Cvi9JHuTbG+OvaOqDkPnH1Pg7WtWu/Ya1Mb2/dVzU5PGcFdXCqbtPWFNWs4FwGPYz1ddT3uD/XzVJNmQ5AngCPBF4I+BF5tvleHkdn21zZvnXwLeNt0az7fe9q6qxT7+C00fvz3JdzTH7ONasWkFOWOviKBlu7iq3ktnqPfGJD+01hVa5+z7q+OTwPfTSXs4DPzH5rjtPUFJ3gT8DvBTVfWnS53a55jtPqY+7W0/X0VV9ZdVdT6dSc4XAu/ud1rz1zZfod72TvIe4GPAucDfBN4KfLQ53fbWik0ryBll1QRNQFU93/w9AvwunQ/ubywO8zZ/j6xdDVtrUBvb91dBVX2j+QfzBPBrvJaqY3tPSJI30Lnh/nRVfa45bD9fJf3a234+HVX1IvAInflQZ6SzchOc3K5jr+yk/rrae2uTqllV9W3g17GPa4KmFeSMsmqCVijJ6UnevPgY+LvAVzl5RYrrgQfWpoatNqiNdwI/1qwUcxHw0mK6j5avJzf7H9Lp59Bp72ublZA205m0+ofTrt+8a+Ya3Ak8WVW/3PWU/XwVDGpv+/nqSbIpyRnN4+8ELqczF+phOis3wev7+FgrO+k1A9r7j7q+NAmd+U/dfdzPFK3IVH4MtKqOJ1lcNWEDcFdV7ZvGtdeZdwC/28yFPBX4rar6QpLdwH1JbgCeBT6whnWce0k+A1wCbExyCPg48Iv0b+NdwJV0JgZ/C/jQ1Cs85wa09yXNUqMFPAP8c4BmNZb7gP10Vqy6sar+ci3qPecuBj4IfKXJoQf4Oeznq2VQe19nP181ZwL3NKvSnUJnNafPJ9kP3Jvk3wGP0wk+af7+ZrOy0zfpfFmr0Q1q799PsolOetoTwL9ozvczRSsWv4iQJEmS1CZT+zFQSZIkSZoGgxxJkiRJrWKQI0mSJKlVDHIkSZIktYpBjiRJkqRWMciRJEmS1CoGOZIkSZJa5f8DL+rH19igU4QAAAAASUVORK5CYII=\n", + "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": [ + "Writing output to c:\\Users\\anbro\\Documents\\projects\\GPUSimulators\\data\\euler_shock-bubble_0001.nc\n", + "Opening c:\\Users\\anbro\\Documents\\projects\\GPUSimulators\\data\\euler_shock-bubble_0001.nc\n", + "Keyword arguments: {'mode': 'w', 'clobber': False}\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Constructed in 2.126896381378174 seconds\n", + "Simulating to t=0.500000. Estimated 1699 timesteps (dt=0.000294)\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Closing c:\\Users\\anbro\\Documents\\projects\\GPUSimulators\\data\\euler_shock-bubble_0001.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": [ + "Writing output to c:\\Users\\anbro\\Documents\\projects\\GPUSimulators\\data\\euler_shock-bubble_0001.nc\n", + "Opening c:\\Users\\anbro\\Documents\\projects\\GPUSimulators\\data\\euler_shock-bubble_0001.nc\n", + "Arguments: ('r',)\n", + "Animation.save using \n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0% [###########################===] 100%. Total: 5s, 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", + " 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_shock-bubble_0001.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.RdYlBu, save_anim=False, fig_args={'figsize':(16, 5), 'tight_layout':True})" + ] + }, + { + "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": "iVBORw0KGgoAAAANSUhEUgAAAzkAAADhCAYAAADrl3vGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAF8BJREFUeJzt3X+w5XV93/HniwVC+BVcQLouGKjZiMYKmg3S0qYISUDaCWSCLTSjlCGzzRRTbUwLOtPSTtoZnTSSZBJltkIWOxakSAphqJQhWuokgguu/NooW3Bgw9YVFsSaiuy97/5xvmuvyz2/zz3n3O8+HzPfOfec7+f7+X7uh88e9r2fz+f9TVUhSZIkSW1x0KwbIEmSJEmTZJAjSZIkqVUMciRJkiS1ikGOJEmSpFYxyJEkSZLUKgY5kiRJklrFIEeSJEnSTCQ5Kcnnk2xP8liS9y9TJkl+P8mOJA8neXu/eg9emeZKkiRJUl97gQ9W1UNJjgIeTHJPVT2+pMy7gA3N8Q7gE81rV87kSJIkSZqJqtpVVQ81P38H2A6s36/YhcCnquNLwDFJ1vWq15kcSZIkSUM7751H1PN7FnqWefDhlx8Dvrfko81VtXm5sklOBt4G3L/fqfXAM0ve72w+29XtvgY5kiRJkob2/J4FHrj79T3LrFn3xPeqamO/upIcCXwW+EBVvbT/6WUuqV71GeRIkiRJGlpRvFJ7x64nySF0ApxPV9VtyxTZCZy05P2JwLO96hxrT06S85N8rcl0cPU4dUmSJElaPQpYpHoe/SQJcD2wvao+1qXYHcB7myxrZwLfrqquS9VgjJmcJGuAPwR+nk509eUkd+yXCUGSJElSSy2yOG4VZwHvAR5Jsq357MPA6wGq6jrgLuACYAfwV8Dl/SodZ7naGcCOqnoSIMnNdDIfGORIkiRJLddZrjZekFNVX2T5PTdLyxRw5TD1jhPkLJfl4FX5qpNsAjYBHHF4fvrUnzh0jFtKkiRJGsc3nnmF5/Ys9AwsBlHAwgBL0mZhnCBnoCwHTYq4zQAbTzusHrj7pFddJEmSJGk6zjjvmf6FBjTIvptZGCfIGTrLgSRJkqR2KOCVms8gZ5zsal8GNiQ5JcmhwCV0Mh9IkiRJarmiWOhzzMrIMzlVtTfJ+4C7gTXADVX12MRaJkmSJGl+FSzM50TOeA8Draq76KR0kyRJknQA6TwnZz6NFeRIkiRJOjAV4ZUaO0nbijDIkSRJkjSShd6PuJkZgxxJkiRJQ+s8J8cgR5IkSVJLdFJIj5OseeUY5EiSJEkaWhEWxnoizcoxyJEkSZI0kkUTD0iSJElqiyJ8v9bMuhnLMsiRJEmSNLTOc3JcriZJkiSpRcyuJkmSJKk1qsIrLleTJEmS1Bad5+S4XE2SJElSa4QFn5MjSZIkqS1MPCBJkiSpVUwhLUmSJKl1Fl2uJkmSJKktTDwgSZIkqVUKU0hLkiRJapEqzK4mSZIkqU3CIpl1I5ZlkCNJkiRpaAV8v+YznJjPVkmSJEmaa0VYLGdyJEmSJLWI2dUkSZIktYbZ1SRJkiS1SjG/DwOdz1ZJkiRJmnsLpOcxiCQ3JNmd5NEu538syZ8k+WqSx5Jc3q9OgxxJkiRJQ6sKi3VQz2NAW4Dze5y/Eni8qk4DzgZ+J8mhvSp0uZokSZKkoRVMZE9OVd2X5OQ+tzoqSYAjgT3A3l51GuRIkiRJGkFY6D9bc1ySrUveb66qzUPe6A+AO4BngaOAf1hVi70u6BvkJDkJ+BTw14DFpmG/l2Qt8BngZOAbwD+oqheGbLAkSZKkVaiTeKDvvpvnqmrjmLc6D9gGnAO8Abgnyf+sqpe6XTDIQrm9wAer6k3AmcCVSd4MXA3cW1UbgHub95IkSZIOAPtSSPc6JuRy4Lbq2AE8BZza64K+QU5V7aqqh5qfvwNsB9YDFwI3NsVuBC4ao+GSJEmSVplFDup5TMjTwLkASU4A3gg82euCofbkNBuC3gbcD5xQVbugEwglee3w7ZUkSZK0GlXBQv/lan0luYlO1rTjkuwErgEO6dyjrgN+C9iS5BEgwFVV9VyvOgcOcpIcCXwW+EBVvdRJbjDQdZuATQCvX2+eA0mSJKkNirB3cSLZ1S7tc/5Z4BeGqXOgOaQkh9AJcD5dVbc1H38zybrm/Dpgd5dGba6qjVW18fhjJ7YuT5IkSdKMTeJhoCuhb5DT5KO+HtheVR9bcuoO4LLm58uA2yffPEmSJEnzaF92tV7HrAyyfuws4D3AI0m2NZ99GPgIcEuSK+hsBnr3yjRRkiRJ0vwJeyeXQW2i+gY5VfVF6DrXdO5kmyNJkiRpNZhU4oGVYCYASZIkSSNZrImliZ4ogxxJkiRJQytmu++ml6kGOV9/+HDOe93p07zlVN397Lb+hfYz6f4YpQ37m8f/RuP8XvPw+wzT/pVu7yTGSC/Dtn8e/tzsb1J9tH87J933sxrbg/we8/DnThrXvrE+rfG80t/P3RwIf15Xsm/nfXws176v1/PjNgfoJB7Y60yOJEmSpDaZ1+Vqqaqp3ezorK13xFwFkiRJ0qzcX/fyUu0Ze53Z2lNfW+fe8Ms9y9x61nUPVtXGce81LGdyJEmSJA3N5WqSJEmSWmXfw0DnkUGOJEmSpJEY5EiSJElqjSIuV5MkSZLUIuVMjiRJkqQWcU+OJEmSpFYpwt5Fl6tJkiRJapFyJkeSJElSmyxikCNJkiSpJcrEA5IkSZLaJSy4J0eSJElSm7gnR5IkSVJrmEJakiRJUrsULBjkSJIkSWqLwuVqkiRJklolLleTJEmS1C6LiwY5kiRJklqiyuVqkiRJklrG5WqSJEmSWsXlapIkSZJao8jcLlc7aNYNkCRJkrQ6VZ9jEEluSLI7yaM9ypydZFuSx5L8j351GuRIkiRJGl6TeKDXMaAtwPndTiY5Bvg48ItV9VPAu/tVOHCQk2RNkq8kubN5f0qS+5M8keQzSQ4dtC5JkiRJq18tpucxUB1V9wF7ehT5R8BtVfV0U353vzqHmcl5P7B9yfuPAtdW1QbgBeCKIeqSJEmStMp10kh3PybkJ4HXJPlCkgeTvLffBQMFOUlOBP4e8MnmfYBzgFubIjcCF43UZEmSJEmrTjHQcrXjkmxdcmwa4VYHAz9NJx45D/hXSX6y3wWD+F3gXwJHNe+PBV6sqr3N+53A+qGbK0mSJGl1KgZZkvZcVW0c8047m3q+C3w3yX3AacDXu13QdyYnyd8HdlfVg0s/XqboshNSSTbti9xe4eV+t5MkSZK0WkwivVp/twN/J8nBSQ4H3sEPb6N5lUFmcs4CfjHJBcBhwNF0ZnaOSXJwM5tzIvDschdX1WZgM8DRWTu5X1WSJEnSDE3mOTlJbgLOprO0bSdwDXAIQFVdV1Xbk3wOeBhYBD5ZVV3TTcMAQU5VfQj4UNOAs4HfrKpfSfJfgIuBm4HL6ERYkiRJkg4Egy1X619N1aUDlPlt4LcHrXOc5+RcBfxGkh109uhcP0ZdkiRJklab6SxXG9qgiQcAqKovAF9ofn4SOGPyTZIkSZK0Oow/k7MShgpyJEmSJOkHFmfdgOUZ5EiSJEkaXgETSDywEgxyJEmSJI2k5jR3skGOJEmSpNFMILvaSjDIkSRJkjSSOJMjSZIkqTVmnCa6F4McSZIkSSOIiQckSZIktYwppCVJkiS1isvVJEmSJLWGz8mRJEmS1DaZ0+VqB826AZIkSZI0Sc7kSJIkSRqJz8mRJEmS1B4FLLonR5IkSVKbOJMjSZIkqU1criZJkiSpXeY0u5pBjiRJkqShpZzJkSRJktQ2PgxUkiRJUqs4kyNJkiSpTeKeHEmSJEmt4Z4cSZIkSa1jkCNJkiSpTeZ1udpBs26AJEmSJE2SMzmSJEmSRuNyNUmSJEmtUfO7XM0gR5IkSdJo5nQmZ6A9OUmOSXJrkr9Isj3J30yyNsk9SZ5oXl+z0o2VJEmSNB9CJ4V0r2NWBk088HvA56rqVOA0YDtwNXBvVW0A7m3eS5IkSToQNMvVeh2DSHJDkt1JHu1T7meSLCS5uF+dfYOcJEcDPwtcD1BV36+qF4ELgRubYjcCF/WrS5IkSVKLVJ9jMFuA83sVSLIG+Chw9yAVDjKT89eBbwF/lOQrST6Z5AjghKraBdC8vnaQG0qSJElqiQkEOVV1H7CnT7FfBz4L7B6kzkGCnIOBtwOfqKq3Ad9liKVpSTYl2Zpk6yu8POhlkiRJkubcAHtyjtsXCzTHpqHvkawHfgm4btBrBsmuthPYWVX3N+9vpRPkfDPJuqralWQdXaKqqtoMbAY4OmvnNP+CJEmSpKEU0H/fzXNVtXHMO/0ucFVVLSQZ6IK+QU5V/e8kzyR5Y1V9DTgXeLw5LgM+0rzePnKzJUmSJK06U8qgthG4uQlwjgMuSLK3qv5rtwsGfU7OrwOfTnIo8CRwOZ2lbrckuQJ4Gnj3OC2XJEmStMpMIcipqlP2/ZxkC3BnrwAHBgxyqmobnQhqf+cO00BJkiRJ7TFomuiedSQ3AWfT2b+zE7gGOASgqgbeh7PUoDM5kiRJkvT/DZcmuns1VZcOUfYfD1LOIEeSJEnS0NIc88ggR5IkSdJIJrFcbSUY5EiSJEkazZw+IMYgR5IkSdJoDHIkSZIktUa5XE2SJElSy0zpYaBDM8iRJEmSNBqDHEmSJElt4kyOJEmSpPYowD05kiRJktoiOJMjSZIkqW0MciRJkiS1RkEW5zPKMciRJEmSNBKXq0mSJElqF4McSZIkSW0Ss6tJkiRJao1yuZokSZKktjHIkSRJktQWwexqkiRJklrG5WpTcPez20a+9rzXnT7Blhw4hulz+3h2Rv2zMex/s2ndZx6M832zz2r8vZda6T//q7WPR2n3ah8L/XTrk7b/3quFf39aPebu+6WY2+VqqZpeyzaedlg9cPdJU7ufJEmSpB92xnnPsPWr38u49Ry59qR6689/oGeZP7/lNx+sqo3j3mtYrZrJkSRJkjQ9ppCWJEmS1B4FTHFV2DAMciRJkiSNxMQDkiRJklqjk0J61q1YnkGOJEmSpOFVuVxNkiRJUru4XE2SJElSq8zrcrWDBimU5J8neSzJo0luSnJYklOS3J/kiSSfSXLoSjdWkiRJ0pwoYLF6HzPSN8hJsh74Z8DGqnoLsAa4BPgocG1VbQBeAK5YyYZKkiRJmjPV55iRgWZy6Cxr+9EkBwOHA7uAc4Bbm/M3AhdNvnmSJEmS5lUWq+cxUB3JDUl2J3m0y/lfSfJwc/xZktP61dk3yKmqvwT+A/A0neDm28CDwItVtbcpthNY36VRm5JsTbL1W88v9LudJEmSpFUi1fsY0Bbg/B7nnwL+blW9FfgtYHO/CgdZrvYa4ELgFOB1wBHAu5YpuuyvUVWbq2pjVW08/tg1/W4nSZIkaTXot1RtwCCnqu4D9vQ4/2dV9ULz9kvAif3qHCS72s8BT1XVtwCS3Ab8LeCYJAc3szknAs8OUJckSZKkFgiQ/s/JOS7J1iXvN1dV35mYHq4A/lu/QoMEOU8DZyY5HPi/wLnAVuDzwMXAzcBlwO0jN1WSJEnSqpOFvkHOc1W1cSL3St5JJ8j52/3KDrIn5346CQYeAh5prtkMXAX8RpIdwLHA9WO0WZIkSdJqMqHlaoNI8lbgk8CFVfV8v/IDPQy0qq4Brtnv4yeBM4ZuoSRJkqQWKOi/XG1sSV4P3Aa8p6q+Psg1AwU5kiRJkrS/QdNE96wjuQk4m87+nZ10JlcOAaiq64B/TWfl2MeTAOzttwTOIEeSJEnS8AqyOIFqqi7tc/5XgV8dpk6DHEmSJEmjmcJytVEY5EiSJEkaySSWq60EgxxJkiRJo3EmR5IkSVJrFDCBPTkrwSBHkiRJ0tBCkcX5jHIMciRJkiSNxuVqkiRJklrD5WqSJEmS2ibO5EiSJElqjwL35EiSJElqjcI9OZIkSZJaZj4ncgxyJEmSJI3GFNKSJEmS2qOARZerSZIkSWqNck+OJEmSpJZxuZokSZKk1nC5miRJkqR2KShnciRJkiS1RQELBjmSJEmS2sTEA5IkSZJaxSBHkiRJUnuYQlqSJElSmxSwsDDrVizLIEeSJEnSaJzJkSRJktQe5XNyJEmSJLVIQblcTZIkSVKruFxNkiRJUmtUwaIPA5UkSZLUIi5XAx58+OXn1qzb8V3guWne9wB3HPb3tNnn02efT5f9PX32+fTZ59Nlf0/Xj0+mGp+TA0BVHZ9ka1VtnOZ9D2T29/TZ59Nnn0+X/T199vn02efTZX+vUsXcZlc7aNYNkCRJkrT6FJ3lar2OQSS5IcnuJI92OZ8kv59kR5KHk7y9X50GOZIkSZKGVwW12PsYzBbg/B7n3wVsaI5NwCf6VTiLIGfzDO55ILO/p88+nz77fLrs7+mzz6fPPp8u+3uVqsXqeQxUR9V9wJ4eRS4EPlUdXwKOSbKuV52pOd0sJEmSJGl+JfkcnaQRvRwGfG/J+81V9aqgNsnJwJ1V9ZZlzt0JfKSqvti8vxe4qqq2drupKaQlSZIkDa2qei0xm6Qsd/teF7gnR5IkSdI82wmctOT9icCzvS6YWpCT5PwkX2uyIlw9rfseaJJ8I8kjSbYl2dp8tjbJPUmeaF5fM+t2rmbLZQDp1sejZAPRD+vS3/8myV8243xbkguWnPtQ099fS3LebFq9uiU5Kcnnk2xP8liS9zefO85XQI/+dpyvkCSHJXkgyVebPv+3zeenJLm/GeOfSXJo8/mPNO93NOdPnmX7V5se/b0lyVNLxvjpzed+p2h/dwDvbcbGmcC3q2pXrwumEuQkWQP8IZ3MCG8GLk3y5mnc+wD1zqo6fUm++auBe6tqA3Bv816j28KrM4B06+Ohs4HoVbawfMaVa5txfnpV3QXQfK9cAvxUc83Hm+8fDWcv8MGqehNwJnBl07eO85XRrb/Bcb5SXgbOqarTgNOB85u/OH2UTp9vAF4ArmjKXwG8UFU/AVzblNPguvU3wL9YMsa3NZ/5nXKASXIT8OfAG5PsTHJFkl9L8mtNkbuAJ4EdwH8E/mm/Oqc1k3MGsKOqnqyq7wM308mSoOm4ELix+flG4KIZtmXV65IBpFsfD50NRD9sgIwrS10I3FxVL1fVU3S+DM9Ysca1VFXtqqqHmp+/A2wH1uM4XxE9+rsbx/mYmrH6f5q3hzRHAecAtzaf7z/G9439W4Fzkyy3R0DL6NHf3fidcoCpqkural1VHVJVJ1bV9VV1XVVd15yvqrqyqt5QVX+jV8KBfaYV5KwHnlnyfie9v8A1ugL+e5IHk2xqPjth35Re8/rambWuvbr1sWN/5byvWcZww5IlmPb3hDXLct4G3I/jfMXt19/gOF8xSdYk2QbsBu4B/hfwYlXtbYos7dcf9Hlz/tvAsdNt8eq2f39X1b4x/u+bMX5tkh9pPnOMa2zTCnKGzoigkZ1VVW+nM9V7ZZKfnXWDDnCO/ZXxCeANdJY97AJ+p/nc/p6gJEcCnwU+UFUv9Sq6zGf2+5CW6W/H+QqqqoWqOp3OBuYzgDctV6x5tc/HtH9/J3kL8CHgVOBngLXAVU1x+1tjm1aQM3RGBI2mqp5tXncDf0zni/ub+6Z5m9fds2tha3XrY8f+Cqiqbzb/w1ykszZ331Id+3tCkhxC5y/cn66q25qPHecrZLn+dpxPR1W9CHyBzn6oY5Lse7zG0n79QZ8353+MwZfRaokl/X1+s1Szqupl4I9wjGuCphXkfBnY0GQtOZTOhsk7pnTvA0aSI5Icte9n4BeAR+n09WVNscuA22fTwlbr1sdDZwNRf/utzf4lOuMcOv19SZMJ6RQ6m1YfmHb7Vrtmr8H1wPaq+tiSU47zFdCtvx3nKyfJ8UmOaX7+UeDn6OyF+jxwcVNs/zG+b+xfDPxp+TT1gXXp779Y8o8mobP/aekY9ztFY5nKw0Cram+S9wF3A2uAG6rqsWnc+wBzAvDHzV7Ig4H/XFWfS/Jl4JYkVwBPA++eYRtXvSYDyNnAcUl2AtcAH2H5Pr4LuIDOxuC/Ai6feoNXuS79fXaTarSAbwD/BKCqHktyC/A4nYxVV1bVwizavcqdBbwHeKRZQw/wYRznK6Vbf1/qOF8x64Abm6x0BwG3VNWdSR4Hbk7y74Cv0Ak+aV7/U5IddGZwLplFo1exbv39p0mOp7M8bRuwNJOW3ykaS/yHCEmSJEltMrWHgUqSJEnSNBjkSJIkSWoVgxxJkiRJrWKQI0mSJKlVDHIkSZIktYpBjiRJkqRWMciRJEmS1Cr/D3HjNPA0mhxEAAAAAElFTkSuQmCC\n", + "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": [ + "Writing output to c:\\Users\\anbro\\Documents\\projects\\GPUSimulators\\data\\euler_kelvin_helmholtz_0001.nc\n", + "Opening c:\\Users\\anbro\\Documents\\projects\\GPUSimulators\\data\\euler_kelvin_helmholtz_0001.nc\n", + "Keyword arguments: {'mode': 'w', 'clobber': False}\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Constructed in 0.008003711700439453 seconds\n", + "Simulating to t=10.000000. Estimated 8998 timesteps (dt=0.001111)\n", + "0% [#####################=========] 100%. Total: 14s, elapsed: 10s, remaining: 4s\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Closing c:\\Users\\anbro\\Documents\\projects\\GPUSimulators\\data\\euler_kelvin_helmholtz_0001.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": [ + "Writing output to c:\\Users\\anbro\\Documents\\projects\\GPUSimulators\\data\\euler_kelvin_helmholtz_0001.nc\n", + "Opening c:\\Users\\anbro\\Documents\\projects\\GPUSimulators\\data\\euler_kelvin_helmholtz_0001.nc\n", + "Arguments: ('r',)\n", + "Animation.save using \n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0% [#####=========================] 100%. Total: 30s, elapsed: 5s, remaining: 25s\n", + "0% [###########===================] 100%. Total: 28s, elapsed: 10s, remaining: 18s\n", + "0% [################==============] 100%. Total: 28s, elapsed: 15s, remaining: 12s\n", + "0% [#####################=========] 100%. Total: 29s, elapsed: 20s, remaining: 8s\n", + "0% [###########################===] 100%. Total: 29s, elapsed: 26s, remaining: 3s\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_kelvin_helmholtz_0001.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), 'tight_layout':True})" + ] + }, + { + "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", @@ -498,87 +237164,164160 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 15, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6cAAAJWCAYAAACkpfF4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzt3X2QbHd5J/bv03P1YsSLXnixLImSiGXHrGNsrAV5X1wEOSAIazkVSCCOrcVKqZzFDgQ7i7Bri2ST2sJJajGuOPaqgBVUORYsJkZxscYyL+tNrZFBGMvCWsE1JuiCbBkQGL8gId1f/ujTPX1Hc++dmZ7u35m5n09Vq/ucPj39dM95dPuZ5zmnq7UWAAAA6GnSOwAAAABQnAIAANCd4hQAAIDuFKcAAAB0pzgFAACgO8UpAAAA3Z22OK2qt1XVA1V198K6/62q/n1V3VVV/3dVnb9w3+ur6mhV3VtVL1xYf+2w7mhV3bT/LwXOTHIUxk2OwrjJURiPnXROb0ly7ZZ1tyf5ztbadyX5VJLXJ0lVPTPJy5P8reEx/2dVbVTVRpJfTPKiJM9M8ophW2B5t0SOwpjdEjkKY3ZL5CiMwmmL09ba7yT58pZ1v9Vae2RY/EiSS4fb1yW5tbX2UGvtT5IcTfKc4XK0tfaZ1trDSW4dtgWWJEdh3OQojJschfHYj2NOfyzJvx5uX5LkvoX7jg3rTrb+Marqxqr62HC5cR/igzOdHIVxk6MwbnIU1uTIMg+uqp9N8kiSX5mt2mazlu2L4Lbdz2yt3Zzk5iR58oUb7apnnfsvlokRxuTOux76YmvtKet6vlXn6Mbjz2vnPP0yOcqh8fB9x+QojNhhy9Gz65z2xLpQjnJofC0PLpWjey5Oq+r6JC9Jck1rbZZ8x5JctrDZpUm+MNw+2fqTuvyys/J777/sdJvBgbFx8dH/b13PtY4cPXLhhfmWn3rNPkQL4/DZ1/y0HIURO2w5em7Oy3MnP7AP0cI4/Pbxf7VUju5prLeqrk3yuiQ/2Fr764W7bkvy8qo6p6quSHJlkt9L8tEkV1bVFVV1dqYHkt+2TODAyclRGDc5CuMmR6GP03ZOq+pXkzwvyZOr6liSN2R6xrJzktxeVUnykdbaj7fWPllV70ryR5mOQLyqtfbo8HN+Isn7k2wkeVtr7ZMreD1wxpGjMG5yFMate47WfpwCBg6H2pxSGJ+rnnVuM9bLYbJx8dE7W2tX9Y5jv5zz9MuakUEOk8++5qflKIzYYcvRJ04ubFcfeeHpN4QD4vZv3LpUjvpTDQAAAN0tdbZeAABgr8pYLyxQnAIAQC+T7b6dBs5M/lQDAABAdzqnAADQQSUZzgYMRHEKAAB9VCUTg4wwIxsAAADoTucUAAB6MdYLc4pTAADoxVgvzMkGAAAAutM5BQCAHsrZemGRzikAAADd6ZwCAEAXlUx0TmFGcQoAAL2UQUaYkQ0AAAB0p3MKAAA9VIz1wgKdUwAAALrTOQUAgF58lQzMKU4BAKCLckIkWCAbAAAA6E7nFAAAenFCJJhTnAIAQA+VlGNOYc5YLwAAAN3pnAIAQC86pzCnOAUAgC4qmRhkhBnZAAAAQHc6pwAA0IuxXphTnAIAQA8VxSksMNYLAABAdzqnAADQy0TnFGZ0TgEAAOhO5xQAALoox5zCAsUpAAD0ojiFOWO9AAAAdKdzCgAAPfgqGTiBzikAAADd6ZwCAEAnzVfJwJziFAAAejHWC3PGegEAAOhO5xQAALrwPaewSHEKAAA9OFsvnMBYLwAAAN3pnAIAQC/O1gtzilMAAOjFWC/MGesFAACgO51TAADoRecU5hSnAADQQ1Wa4hTmjPUCAADQnc4pAAD0olUEc9IBAACA7nROAQCgF8ecwpziFAAAelGcwpyxXgAAALrTOQUAgA5axVfJwAKdUwAAALrTOQUAgF60imBOcQoAAL0Y64U5f6sBAACgO51TAADoonROYYHiFAAAeqjpGXuBKWO9AAAAdKdzCgAAvUy0TmFG5xQAAIDudE4BAKCT5oRIMKc4BQCAXtSmMGesFwAAgO50TgEAoIeK7zmFBYpTAADoxPecwiZjvQAAAHSncwoAAL0Y64U5xSkAAPRijhHmpAMAAADd6ZwCAEAPVWnGemFOcQoAAL2oTWHOWC8AAADd6ZwCAEAHLTHWCwt0TgEAAOhO5xQAAHrROIU5xSkAAHTSFKcwZ6wXAACA7nROAQCgh0oy0TqFmdN2TqvqbVX1QFXdvbDuwqq6vao+PVxfMKyvqvqFqjpaVXdV1bMXHnP9sP2nq+r61bwcOPPIURg3OQrjJkdhPHYy1ntLkmu3rLspyQdaa1cm+cCwnCQvSnLlcLkxyS8l0wRP8oYkz03ynCRvmCU5sLRbIkdhzG6JHIUxuyUdc7SVi8vhuSzrtMVpa+13knx5y+rrkrx9uP32JD+0sP4dbeojSc6vqouTvDDJ7a21L7fWHkxyex77PwFgD+QojJschXHrnqNVLi6H57KkvZ4Q6WmttfuTZLh+6rD+kiT3LWx3bFh3svXAashRGDc5CuMmR6GD/T4h0nblcjvF+sf+gKobMx2TyNMvcb4m2Gf7mqMbF5gqhH0mR2Hc9jVHzzn3/H0ZhYTDYq+d0z8bRhgyXD8wrD+W5LKF7S5N8oVTrH+M1trNrbWrWmtXPeWijT2GB2e8teToxuPP2/fA4QwhR2Hc1pKjZ51z3rSsdXE5LJcl7bU4vS3J7Cxk1yd578L6Hx3OZHZ1kq8OoxDvT/KCqrpgODj8BcM6YDXkKIybHIVxk6PQwWnnZqvqV5M8L8mTq+pYpmcie2OSd1XVDUk+l+Rlw+bvS/LiJEeT/HWSVyZJa+3LVfU/J/nosN0/ba1tPfAc2AM5CuMmR2HcuufoPnSb4LA4bXHaWnvFSe66ZpttW5JXneTnvC3J23YVHXBachTGTY7CuPXO0bYPZziFw2KvY70AAACwb5wOFwAAetE4hTnFKQAAdNAqaeYYYU46AAAA0J3OKQAA9GKsF+Z0TgEAAOhO5xQAALooXyUDCxSnAADQi9oU5oz1AgAA0J3OKQAA9FDTr5MBpnROAQAA6E7nFAAAetE5hTnFKQAAdGKsFzYZ6wUAAKA7nVMAAOjF95zCnOIUAAA6MdYLm4z1AgAA0J3OKQAA9FBxtl5YoDgFAIAOWoz1wiJjvQAAAHSncwoAAL3onMKc4hQAADox1gubjPUCAADQnc4pAAD0onMKczqnAAAAdKdzCgAAPVTSSusUZhSnAADQi9oU5oz1AgAA0J3OKQAAdOKrZGCTzikAAADd6ZwCAEAvOqcwpzgFAIBOjPXCJmO9AAAAdKdzCgAAPVSM9cICxSkAAHTQYqwXFhnrBQAAoDudUwAA6EXnFOYUpwAA0IviFOaM9QIAANCdzikAAHTihEiwSecUAACA7nROAQCgB99zCidQnAIAQCfGemGTsV4AAAC60zkFAIBedE5hTnEKAACdtFKdwoyxXgAAALrTOQUAgF40TmFO5xQAAIDudE4BAKAH33MKJ1CcAgBAJ77nFDYZ6wUAAKA7nVMAAOhF5xTmdE4BAADoTucUAAA6aHHMKSxSnAIAQC+KU5gz1gsAAEB3OqcAANCD7zmFEyhOAQCgE8ecwiZjvQAAAHSncwoAAL3onMKc4hQAADox1gubjPUCAADQnc4pAAD0onMKc4pTAADowVfJwAmM9QIAANCdzikAAHTihEiwSecUAACA7nROAQCgF51TmFOcAgBAJ8Z6YZOxXgAAALpTnAIAANCdsV4AAOjFWC/M6ZwCAADQnc4pAAB00MoJkWCRzikAAADd6ZwCAEAvOqcwpzgFAIBeFKcwZ6wXAACA7nROAQCgEydEgk2KUwAA6EVxCnPGegEAAOhuqeK0qv77qvpkVd1dVb9aVedW1RVVdUdVfbqq3llVZw/bnjMsHx3uv3w/XgBwcnIUxk2OwritPEfLxeWQXZa05+K0qi5J8t8luaq19p1JNpK8PMnPJXlTa+3KJA8muWF4yA1JHmytfWuSNw3bASsiR2Hc5CiM27pytJWLy+G5LGvZsd4jSb6pqo4keVyS+5M8P8m7h/vfnuSHhtvXDcsZ7r+mqvbhJQCnIEdh3OQojJschTXa8wmRWmufr6r/PcnnkvxNkt9KcmeSr7TWHhk2O5bkkuH2JUnuGx77SFV9NclFSb64+HOr6sYkNybJ0y9xvibYq3Xk6MYFF6z6ZcChJUdh3NaRo0eedEFSbdUvBQ6MZcZ6L8j0L0RXJPmWJOcledE2m84ybru/HD0mG1trN7fWrmqtXfWUizb2Gh6c8daRoxuPP2+/woUzjhyFcVtLjp53XvcxTBeX/bwsa5mx3h9I8iettT9vrX0jyXuS/J0k5w+jD0lyaZIvDLePJbksSYb7n5Tky0s8P3BqchTGTY7CuMlRWLNlitPPJbm6qh43zNNfk+SPknwoyUuHba5P8t7h9m3Dcob7P9haM8cAqyNHYdzkKIzbenK0XFwO0WVJey5OW2t3ZHqw98eT/OHws25O8rokr62qo5nO2b91eMhbk1w0rH9tkpuWiBs4DTkK4yZHYdzkKKzfUmccaq29Ickbtqz+TJLnbLPt15O8bJnnA3ZHjsK4yVEYt5Xn6D51m+CwcDpcAADoZD9OIgOHxbLfcwoAAABL0zkFAIBedE5hTucUAACA7nROAQCgF51TmFOcAgBAB224AFPGegEAAOhO5xQAAHox1gtzilMAAOihojiFBcZ6AQAA6E7nFAAAOmk6pzCnOAUAgF7K+XphxlgvAAAA3emcAgBAL8Z6YU5xCgAAnTjmFDYZ6wUAAKA7nVMAAOjB95zCCXROAQAA6E7nFAAAetE5hTnFKQAAdOKESLDJWC8AAADd6ZwCAEAv1XpHAKOhcwoAAEB3OqcAANCLY05hTnEKAAA9+J5TOIGxXgAAALrTOQUAgF50TmFO5xQAAIDudE4BAKCLluarZGBOcQoAAL0Y64U5Y70AAAB0p3MKAAC96JzCnM4pAAAA3SlOAQAA6M5YLwAA9FBJnK0X5hSnAADQi2NOYc5YLwAAAN3pnAIAQC86pzCnOAUAgF4ccwpzxnoBAADoTucUAAB6MdYLczqnAAAAdKdzCgAAPVR0TmGB4hQAAHpxQiSYM9YLAABAdzqnAADQSRnrhTmdUwAAALrTOQUAgF4ccwpzilMAAOjB2XrhBMZ6AQAA6E7nFAAAumgpY70wpzgFAIBejPXCnLFeAAAAutM5BQCATnzPKWxSnAIAQC+OOYU5Y70AAAB0p3MKAACdmOqFTYpTAADooCq+SgYWGOsFAACgO51TAADoRecU5nROAQAA6E7nFAAAOvE9p7BJcQoAAJ04IRJsMtYLAABAdzqnAADQibFe2KQ4BQCAHqoZ64UFxnoBAADoTucUAAA6qDghEizSOQUAAKA7nVMAAOjECZFgk+IUAAA6MdYLm4z1AgAA0J3OKQAAdKJzCpsUpwAA0IlDTmGTsV4AAAC60zkFAIAOqpKJsV6YU5wCAEAXzTGnsMBYLwAAAN3pnAIAQCc6p7BJcQoAAB1UkonT9cKc4vQQeeG3fPfan/P9X/jE2p8TAIAzT48uc2v+erBOSxWnVXV+krck+c4kLcmPJbk3yTuTXJ7ks0n+i9bag1VVSd6c5MVJ/jrJP2ytfXyZ5z8T9Cg4d2O38Slm10uOwrjJURi3deSosd5x8/tZr2U7p29O8puttZdW1dlJHpfkZ5J8oLX2xqq6KclNSV6X5EVJrhwuz03yS8P1GW3sxed+O93rVbzuOzkK4yZHYdzk6JJ8Vc6JjuvEntKei9OqemKS70/yD5OktfZwkoer6rokzxs2e3uSD2easNcleUdrrSX5SFWdX1UXt9bu33P0B8CZVnwuS/G6f+QojJschXFbS476ntMzjt/3qS3TOX1Gkj9P8i+r6llJ7kzy6iRPmyVha+3+qnrqsP0lSe5bePyxYd0JCVtVNya5MUmefsk4D4lVcPaz0/deEZtkDTm6ccEFK30BcMjJURi3lefo2U994ijHRh1nOX5j3G/2wzLV35Ekz07yk621O6rqzZmONZzMdnv5Y97V1trNSW5OkquedW6Xd13xefDpwCZZQ46e8/TLDuf/GWE95CiM28pz9Lxvu3iUOXpYCx/Gb5ni9FiSY621O4bld2easH82G2GoqouTPLCw/WULj780yReWeP49UXiS7Gw/OAQF7IHMUTiDyFEYt5XnaKXt+5inYxrZqTGOGO+5OG2t/WlV3VdV395auzfJNUn+aLhcn+SNw/V7h4fcluQnqurWTA8O/+rpjpP51F2PU0zSzWr2vaMr+JnbW0eOAnsnR2Hc1pGjLaWYpJsx7nvLHtT5k0l+ZTh72WeSvDLJJMm7quqGJJ9L8rJh2/dlemrto5meXvuVSz43cHpyFMZNjsK4yVFYo6WK09baJ5Jctc1d12yzbUvyqmWeD9gdOQrjJkdh3Fado6sY64WDbNI7AAAAAFCcAgAA0N04v0gUAADOAOM7JQ30o3MKAABAd4pTAAAAujPWCwAAnThbL2zSOQUAAKA7xSkAAADdGesFAIAOKsZ6YZHOKQAAAN0pTgEAAOjOWC8AAPRQSRnrhTmdUwAAALpTnAIAANCdsV4AAOjE2Xphk84pAAAA3SlOAQAA6M5YLwAAdFBpxnphgc4pAAAA3SlOAQAA6E5xCgAAQHeOOQUAgE4mccwpzOicAgAA0J3iFAAAgO6M9QIAQCflq2RgTucUAACA7hSnAAAAdGesFwAAOqgkE2O9MKdzCgAAQHeKUwAAALoz1gsAAD2UsV5YpHMKAABAd4pTAAAAujPWCwAAHVSasV5YoHMKAABAd4pTAAAAujPWCwAAnUxirBdmdE4BAADoTnEKAABAd4pTAAAAunPMKQAAdDKp471DgNHQOQUAAKA7xSkAAADdGesFAIAOKsmkfJUMzOicAgAA0J3iFAAAgO6M9QIAQBfNWC8s0DkFAACgO8UpAAAA3RnrBQCADqqcrRcW6ZwCAADQneIUAACA7oz1AgBAJ5MY64UZnVMAAAC6U5wCAADQnbFeAADooOJsvbBI5xQAAIDuFKcAAAB0pzgFAACgO8ecAgBAFy2TOt47CBgNnVMAAAC6U5wCAADQnbFeAADooJJs+CoZmNM5BQAAoDvFKQAAAN0Z6wUAgE4mMdYLMzqnAAAAdKc4BQAAoDtjvQAA0MnE2XphTucUAACA7hSnAAAAdGesFwAAOqhqmdTx3mHAaOicAgAA0J3iFAAAgO6M9QIAQCcbztYLczqnAAAAdKc4BQAAoDvFKQAAAN055hQAADqoJJM45hRmdE4BAADoTnEKAABAd0uP9VbVRpKPJfl8a+0lVXVFkluTXJjk40l+pLX2cFWdk+QdSb43yZeS/Jettc8u+/zAqclRGDc5CuO26hyd1PGVxQ4HzX50Tl+d5J6F5Z9L8qbW2pVJHkxyw7D+hiQPtta+Ncmbhu2A1ZOjMG5yFMZNjsKaLFWcVtWlSf7TJG8ZlivJ85O8e9jk7Ul+aLh93bCc4f5rhu2BFZGjMG5yFMZNjsJ6LTvW+/NJ/nGSJwzLFyX5SmvtkWH5WJJLhtuXJLkvSVprj1TVV4ftv7j4A6vqxiQ3Jsm5edyS4cEZb6U5unHBBSsNHs4AchTGbaU5+oRvflwm5Wy9MLPnzmlVvSTJA621OxdXb7Np28F9mytau7m1dlVr7aqzcs5ew4Mz3jpydOPx5+1DpHBmkqMwbuvI0W+6wGddWLRM5/TvJvnBqnpxknOTPDHTvy6dX1VHhr8oXZrkC8P2x5JcluRYVR1J8qQkX17i+YFTk6MwbnIUxk2OwprtuXPaWnt9a+3S1trlSV6e5IOttR9O8qEkLx02uz7Je4fbtw3LGe7/YGvNHAOsiByFcZOjMG7rytGNHHdxOTSXZa3ie05fl+S1VXU00zn7tw7r35rkomH9a5PctILnBk5PjsK4yVEYNzkKK7L095wmSWvtw0k+PNz+TJLnbLPN15O8bD+eD9gdOQrjJkdh3OQorMe+FKcAAMDuVOJsvbBgFWO9AAAAsCuKUwAAALoz1gsAAD2UsV5YpHMKAABAd4pTAAAAulOcAgAA0J1jTgEAoINKy0aO9w4DRkPnFAAAgO4UpwAAAHRnrBcAADrxVTKwSecUAACA7hSnAAAAdGesFwAAOqi0bJSz9cKMzikAAADdKU4BAADozlgvAAB0Momz9cKMzikAAADdKU4BAADozlgvAAB0UImz9cICnVMAAAC6U5wCAADQnbFeAADoxNl6YZPOKQAAAN0pTgEAAOjOWC8AAHRQac7WCwt0TgEAAOhOcQoAAEB3ilMAAAC6c8wpAAD0UMnEMacwp3MKAABAd4pTAAAAujPWCwAAHVSSjbTeYcBo6JwCAADQneIUAACA7oz1AgBAF83ZemGBzikAAADdKU4BAADozlgvAAB04Gy9cCKdUwAAALpTnAIAANCd4hQAAIDuHHMKAACd+CoZ2KRzCgAAQHeKUwAAALoz1gsAAB1UWjZirBdmdE4BAADoTnEKAABAd8Z6AQCgk0m13iHAaOicAgAA0J3iFAAAgO6M9QIAQAeVOFsvLNA5BQAAoDvFKQAAAN0Z6wUAgC5aNspYL8zonAIAANCd4hQAAIDuFKcAAAB055hTAADooJJMfJUMzOmcAgAA0J3iFAAAgO6M9QIAQA+VbFTrHQWMhs4pAAAA3SlOAQAA6M5YLwAAdFBp2XC2XpjTOQUAAKA7xSkAAADdGesFAIBOJmWsF2Z0TgEAAOhOcQoAAEB3xnoBAKCDSrKR1jsMGA2dUwAAALpTnAIAANCdsV4AAOig0rLhbL0wp3MKAABAd4pTAAAAulOcAgAA0J1jTgEAoJNJHHMKMzqnAAAAdKc4BQAAoDtjvQAA0IGvkoET7blzWlWXVdWHquqeqvpkVb16WH9hVd1eVZ8eri8Y1ldV/UJVHa2qu6rq2fv1IoDHkqMwbnIUxk2OwvotM9b7SJKfaq19R5Krk7yqqp6Z5KYkH2itXZnkA8NykrwoyZXD5cYkv7TEcwOnJ0dh3OQojJschTXb81hva+3+JPcPt79WVfckuSTJdUmeN2z29iQfTvK6Yf07WmstyUeq6vyqunj4OcA+k6MwbnIUxm1dObqRtpoXAAfQvpwQqaouT/I9Se5I8rRZEg7XTx02uyTJfQsPOzasA1ZMjsK4yVEYNzkK67F0cVpVj0/ya0le01r7i1Ntus26x/ypqKpurKqPVdXHvpGHlg0PznirzNFH//Kv9itMOGPJURi3VeboV7/86H6FCYfCUmfrraqzMk3WX2mtvWdY/WezEYaqujjJA8P6Y0kuW3j4pUm+sPVnttZuTnJzkjyxLjTnAEtYdY6e8/TL5CgsQY7CuK06R7/tP/qmNomz9cLMMmfrrSRvTXJPa+2fL9x1W5Lrh9vXJ3nvwvofHc5kdnWSrzpOBlZHjsK4yVEYNzkK67dM5/TvJvmRJH9YVZ8Y1v1MkjcmeVdV3ZDkc0leNtz3viQvTnI0yV8neeUSzw2cnhyFcZOjMG5yFNZsmbP1/r/ZfrY+Sa7ZZvuW5FV7fT5gd+QojJschXFbT462bJSxXpjZl7P1AgAAwDIUpwAAAHS31Nl6AQCAvakkG4/9thk4Y+mcAgAA0J3iFAAAgO4UpwAAAHTnmFMAAOihkomvkoE5nVMAAAC6U5wCAADQnbFeAADooNKyEWO9MKNzCgAAQHeKUwAAALoz1gsAAJ1spPUOAUZD5xQAAIDuFKcAAAB0Z6wXAAA6qCSTcrZemNE5BQAAoDvFKQAAAN0Z6wUAgE6crRc26ZwCAADQneIUAACA7oz1AgBAB5VmrBcW6JwCAADQneIUAACA7oz1AgBAJ5My1gszOqcAAAB0pzgFAACgO8UpAAAA3TnmFAAAOqjEV8nAAp1TAAAAulOcAgAA0J2xXgAA6GRirBfmdE4BAADoTnEKAABAd8Z6AQCgg0qyUcZ6YUbnFAAAgO4UpwAAAHRnrBcAADqotGw4Wy/M6ZwCAADQneIUAACA7oz1AgBAJzpFsEk+AAAA0J3iFAAAgO6M9QIAQAeVZKN6RwHjoXMKAABAd4pTAAAAulOcAgAA0J1jTgEAoBOdItgkHwAAAOhOcQoAAEB3xnoBAKCDSrLROwgYEZ1TAAAAulOcAgAA0J2xXgAA6KEqG1W9o4DR0DkFAACgO8UpAAAA3RnrBQCADio6RbBIPgAAANCd4hQAAIDujPUCAEAnG3G2XpjROQUAAKA7xSkAAADdGesFAIAOKsmkjPXCjM4pAAAA3SlOAQAA6E5xCgAAQHeOOQUAgE58lQxs0jkFAACgO8UpAAAA3RnrBQCADiqViV4RzMkGAAAAulOcAgAA0J2xXgAA6GSjnK0XZnROAQAA6E5xCgAAQHfGegEAoINKnK0XFsgGAAAAulOcAgAA0J2xXgAA6GQSZ+uFGZ1TAAAAulOcAgAA0J2xXgAA6KBS2Si9IphZezZU1bVVdW9VHa2qm9b9/MCpyVEYNzkK4yZHYe/WWpxW1UaSX0zyoiTPTPKKqnrmOmMATk6OwrjJURg3OQrLWXfn9DlJjrbWPtNaezjJrUmuW3MMwMnJURg3OQrjJkdhCes+5vSSJPctLB9L8tzFDarqxiQ3Dot/+dvt3V9K8sX1hLdvnhwxr8NBjPnbewdwGrvO0c++5qfl6HqIeT0OW44+9NnX/PTda4ptPx3EfUfM63HocnTj4qNydD3EvB5L5ei6i9PtvsipnbDQ2s1Jbp4/oOpjrbWrVh3YfhLzehzUmHvHcBpydKTEvB6HLUcP4u8gOZhxi3k95Og4HMS4xbwey+bousd6jyW5bGH50iRfWHMMwMnJURg3OQrjJkdhCesuTj+a5MqquqKqzk7y8iS3rTkG4OTkKIybHIVxk6OwhLWO9bbWHqmqn0jy/iQbSd7WWvvkaR5282nuHyMxr4eY95kcHTUxr8eoY95Djo769ZzCQYxbzOsx6pjl6KiJeT2Wirlaa6ffCgAAAFZo3WO9AAAA8BiKUwAAALobbXFaVddW1b2wm2HiAAAQFklEQVRVdbSqbuodz0xVva2qHqiquxfWXVhVt1fVp4frC4b1VVW/MLyGu6rq2Z1ivqyqPlRV91TVJ6vq1WOPu6rOrarfq6o/GGL+n4b1V1TVHUPM7xxONpCqOmdYPjrcf/m6Y16IfaOqfr+qfuOgxLwXcnRfY5aj641djnYkR9cWsxwdOTm6rzHL0fXGvrIcHWVxWlUbSX4xyYuSPDPJK6rqmX2jmrslybVb1t2U5AOttSuTfGBYTqbxXzlcbkzyS2uKcatHkvxUa+07klyd5FXD+znmuB9K8vzW2rOSfHeSa6vq6iQ/l+RNQ8wPJrlh2P6GJA+21r41yZuG7Xp5dZJ7FpYPQsy7Ikf3nRxdLzna1y2Ro+sgR0dMju47Obpeq8vR1troLkm+L8n7F5Zfn+T1veNaiOfyJHcvLN+b5OLh9sVJ7h1u/4skr9huu87xvzfJf3JQ4k7yuCQfT/LcJF9McmTrfpLpWfG+b7h9ZNiuOsR6aab/83t+kt/I9Mu4Rx3zHl+nHF1t/HJ0dbHK0RFc5Oja45WjI7vI0ZXHL0dXF+tKc3SUndMklyS5b2H52LBurJ7WWrs/SYbrpw7rR/c6hnb69yS5IyOPexgZ+ESSB5LcnuSPk3yltfbINnHNYx7u/2qSi9YbcZLk55P84yTHh+WLMv6Y92IU+8gujHpfXyRHV06OjtOo9/VFcnTl5Og4jXpfXyRHV26lOTrW4rS2WXcQv/NmVK+jqh6f5NeSvKa19hen2nSbdWuPu7X2aGvtuzP9C81zknzHdpsN191jrqqXJHmgtXbn4uptNh1NzEs4yLEvGtXrkKOrJUcPTOyLRvU65OhqydEDE/uiUb0OObpa68jRUX7PaVV9X5L/sbX2wmH59Rflaf/s4TycVE1f5eZ/ktp6vXhfTrwvdeLbVCdu17bev83PaIvrHnO7Nt/xbbY52WNP9ZiZdrK4Tvdztiyf8uecsNxOGc90uc1v15b7Ft6Jhbd56/btMfctLtfsaapt/hq2e8yw7oTlJFWz5ZPfN9m6fh5H24xxSzyLsVdqIbbF5en1nXc99I0kj2vTL+ae79tV9f7h9u9W1ZEkf5rkKW2MSbnF9jn6zf/s4TyU2jYPT5ejizvQLnN0m23HmqMnzc8t69opnuPEdavL0Vp4jr3n6GPXnTxHt+Tfwn2ny9HZ8mTxsQuxy9GT5eiW6xPWzf9z6hx9zHanytFa2CYn2Sbb5+gO8vB0OfqY59xm2538G5qcJkdPWN5jjm79NzM7z9Ha8vid5Ojiz9x5jm7m9Ak/94Sfd+ocreH9kaM7zNHt8nLx9mHL0W3yZaw5esL/RneYo7UlH3aTo/N/47LzHN2a28nuc3QzvtXl6JGT3dHZR5NcWVVXJPl8kpc/nIdz9ZEXJDVJbWwkk5peVyXbLVclk0kyqW3Xtclk2jeeTKZJtjHdbr6+Km14bNuo6U4+qbT5tkmbDOuH+9vscYu3J5lf5svz+4f7hvWz3/p226SGRFu4P1uvtz4mW58nyaQ95jHzdbW5TRau23Bds/U1LFebr6tJS1XLZFie1ObypFomk+PTX8PkeCZDsmxMjmdSbb5utnxkWD4yW67Z+kczyfTnnTXcPjJ5NBuz62o5Uo9mI8fn9581eSQbaZnU8Zw1u68ezcawPMnxnF2PZlLDdY7n7Hokk7T58nTb48Njj2cjLWdVG66TjVTOqsoklbMyyUZVjmQjGzXJJJNsXPzpv0zy0iS3Jrk+0+MgkuS2Yfl3h/s/eBD+QR1sk6MP5eojLxjycLK7HJ1sDNeTYf3pc7RNJtP9c7Z+Nzk6z7HazI8d5uhjcmoXObr1ZyQ7yNFJktpjjg55WcPP2E2ObgyPneXlbnJ0MR8nQ66cLkfPGvJwlqe7ydGzh+c5O5vPJ0d3kKMbk9RkspCPu8jRjUlaVTLk345ydJZvG9M83VWObmxZ3k2Ozn7W7N+8XeTotjl+uhyd52fmubjTHJ3Mlzc/dO40RzdqIT+HnN1pjp615XonObr4b+lGju8qR8/O8UwqOSstG3J0xzk6+/w7zcvZv5c7y9HFf1Nb1Y5ydDEfM9lljm7NqV3k6LZ5uNMcneXjlseeKkdnj5nl4W5ydDEvq9qucvRInfhv6m5y9ITPsvPPpqfO0bPqkWxUm17n+K5y9KxUNmqan5Ph9qpydJRjvcNM8k9kehDtPUne1TciWMqxJK+tqqOZztm/dVj/1iQXDetfm80zyI2eHOWQkaMwbnIUxm3fcnSsndO01t6X5H2z5SfWhf9Lx3BgGQ+31p6zdWVr7etJXtYhnn0hRzlE5CiMmxyFcdu3HB1l5xQAAIAzi+IUAACA7hSnAAAAdKc4BQAAoDvFKQAAAN3VQfk6qKq6O8nXe8cxeHKSL/YOYiCW7Y0plnNba9/ZO4hVG1mO7sSY9pGdOEjxHqRYEzk6VgdpPzpIsSYHL94zJUd/M9PfzRgdtH1mv3jdO/PF1tq1+/HEo/0qmW18vbV2Ve8gkqSqPiaWxxLL9qrqY71jWJPR5OhOjGkf2YmDFO9BijWRo2N1kPajgxRrcjDj7R3DOuzXh/tVOGj7zH7xutfPWC8AAADdKU4BAADo7iAVpzf3DmCBWLYnlu2NKZZVOmivU7yrc5BiTQ5evHt10F7nQYr3IMWaiJfdO1N/B173mh2YEyIBAABweB2kzikAAACHlOIUAACA7kZXnFbV26rqgeH72La7v6rqF6rqaFXdVVXP7hTHDw/Pf1dV/buqetYq4thJLAvb/e2qerSqXtozlqp6XlV9oqo+WVX/plcsVfWkqvp/quoPhlheucJYLquqD1XVPcNzvXqbbday765SVV1YVbdX1aeH6wtOst1vVtVXquo3tqy/oqruGB7/zqo6eyTxXj9s8+mqun5h/Yer6t5hf/5EVT11BTFeOzzH0aq6aZv7zxneq6PDe3f5wn2vH9bfW1Uv3O/Y9jPeqrq8qv5m4b385RHE+v1V9fGqemTr/zdPtk+MnRyVo3KUZVTV+VX17qr69zX9TPN9W+6vOuCfZbazg9e9ts/963S6172w3cprjLnW2qguSb4/ybOT3H2S+1+c5F8nqSRXJ7mjUxx/J8kFw+0XrSqOncQybLOR5INJ3pfkpR1/P+cn+aMkTx+Wn9oxlp9J8nPD7ack+XKSs1cUy8VJnj3cfkKSTyV55pZt1rLvrvKS5H9NctNw+6bZ+7vNdtck+QdJfmPL+ncleflw+5eT/Le9401yYZLPDNcXDLdnuf3hJFetML6NJH+c5BlJzk7yB9vsN/8oyS8Pt1+e5J3D7WcO25+T5Irh52ys+P1cJt7LT/X/sE6xXp7ku5K8Y/H/m6faJ8Z+kaOj2ufl6PKxHrocHfslyduT/DfD7bOTnL/l/gP/WWaPr3ttn/vH9LqH9WupMWaX0XVOW2u/k2kRcTLXJXlHm/pIkvOr6uJ1x9Fa+3ettQeHxY8kuXS/Y9hpLIOfTPJrSR5YVRw7jOW/SvKe1trnhu1XFs8OYmlJnlBVleTxw7aPrCiW+1trHx9ufy3JPUku2bLZWvbdFbsu0/+RZbj+oe02aq19IMnXFtcNv4fnJ3n36R6/j3YS7wuT3N5a+/KQ07cnWdcXoT8nydHW2mdaaw8nuTXTmBctvoZ3J7lmeC+vS3Jra+2h1tqfJDk6/Lyxxrtup421tfbZ1tpdSY5veWzPfWJZcnR/ydHVOVNzdLSq6omZ/uH/rUnSWnu4tfaVLZsdhs8yJ9jJ617n5/512eHvO1lTjTEzuuJ0By5Jct/C8rE8tghYtxsy/StSF1V1SZL/LNO/cvf2bUkuGEat7qyqH+0Yy/+R5DuSfCHJHyZ5dWtt6z9w+24YkfqeJHdsuWuM++5uPa21dn8yLciT7GaE7qIkX2mtzf5AsI7Xv5N4T/d7+ZfDiNs/WcEHuJ3sE/Nthvfuq5m+lz32p2XiTZIrqur3q+rfVNXfH0Gsq3hsb3J0f8nRvrGu4rGc3DOS/HmmOfX7VfWWqjpvyzaH8b3fyete1PVz/z467evuUWMcWdcT7aPt/uHp9n04VfUfZ7qT/r1eMST5+SSva6092uePnyc4kuR7Mx0Z+6Ykv1tVH2mtfapDLC9M8olMOwH/QZLbq+rfttb+YlVPWFWPz/SvS6/Z5nlGte+eTFX9dpJv3uaun132R2+zbunXvw/xniquH26tfb6qnpDp7/VHMh0v2y87eU9Otk2P/WmZeO/PdNz/S1X1vUl+var+1grzcZn3Z9S5KkeTyNGTkaMs40imh0v9ZGvtjqp6c6bj9v9kYZvD+N7v5HUnGc3n/v2yk9e99hrjIBanx5JctrB8aaadsbWrqu9K8pYkL2qtfalHDIOrktw67DRPTvLiqnqktfbrHWI5luSLrbW/SvJXVfU7SZ6V6TGY6/bKJG9s04H5o1X1J0n+wyS/t4onq6qzMv1w9Cuttfdss8lo9t1Taa39wMnuq6o/q6qLW2v3D2M8uxnx+GKm4z9Hhr/W78vr34d4jyV53sLypZkex5bW2ueH669V1f+V6Rjafn7w3ck+MdvmWFUdSfKkTEfUe+xPe453yMOHkqS1dmdV/XGmkxYf6xjrqR77vC2P/fC+RLUP5KgcPQU5yjKOJTnWWptNfr0702Jl6zaj/yyzSzt53WP63L9fdvK6115jHMSx3tuS/GhNXZ3kq7ORoHWqqqcneU+SH+nUFZxrrV3RWru8tXZ5pjvWP+pUmCbJe5P8/ao6UlWPS/LcTI+/7OFzmXZwU1VPS/LtmZ40Yd8No2RvTXJPa+2fn2SzUey7S7otyeysiNdn+vvekeGDz4eSzM70tqvH79FO4n1/khdU1QU1PVPoC5K8f9iHn5zM//DwkiSnPGP2Hnw0yZU1PUPq2ZmenOS2U7yGlyb54PBe3pbk5TU98+YVSa7Miv7wsh/xVtVTqmojSarqGUO8K8nHXcR6MtvuEyuKc7/J0f0lR/vGejIHOUdHq7X2p0nuq6pvH1Zdk+lJLhcdhs8yJ9jJ6x7T5/79spPX3aXGaCM4U9TiJcmvZjpa8o1MK/obkvx4kh8f7q8kv5jpGd7+MCs6S98O4nhLkgczHRv9RJKP9XpPtmx7S1Z7tt7TxpLkf8h057470/HWXvvKtyT5rWE/uTvJf73CWP5epmMtdy3sEy/use+u8pLpcUkfSPLp4frCYf1VSd6ysN2/zfQ4hr8ZfjcvHNY/I9MPZ0eT/Ksk54wk3h8bYjqa5JXDuvOS3Dn8Tj+Z5M1ZwZk2h/3kU8N+8bPDun+a5AeH2+cO79XR4b17xsJjf3Z43L2Z/iV3HfvAnuJN8p8P7+MfJPl4kn8wglj/9rB//lWSLyX55Kn2iYNwkaNyVI66LPk7+e5Mu+V3Jfn1TM+GfKg+y+zxda/tc/+YXveWbW/JGs7WW8OTAQAAQDcHcawXAACAQ0ZxCgAAQHeKUwAAALpTnAIAANCd4hQAAIDuFKcAAAB0pzgFAACgu/8fFIkpHMSWI+cAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ - "nx = 301\n", + "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": 16, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Writing output to c:\\Users\\anbro\\Documents\\projects\\GPUSimulators\\data\\euler_rayleigh-taylor.nc\n", + "Opening c:\\Users\\anbro\\Documents\\projects\\GPUSimulators\\data\\euler_rayleigh-taylor.nc\n", + "Keyword arguments: {'mode': 'w', 'clobber': False}\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Constructed in 0.014011383056640625 seconds\n", + "Simulating to t=30.000000. Estimated 30697 timesteps (dt=0.000977)\n", + "0% [#######=======================] 100%. Total: 45s, elapsed: 10s, remaining: 35s\n", + "0% [#############=================] 100%. Total: 45s, elapsed: 20s, remaining: 25s\n", + "0% [####################==========] 100%. Total: 46s, elapsed: 30s, remaining: 15s\n", + "0% [##########################====] 100%. Total: 46s, elapsed: 40s, remaining: 5s\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Closing c:\\Users\\anbro\\Documents\\projects\\GPUSimulators\\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", - "rho0, rho_u0, rho_v0, E0, dx, dy, dt = genRayleighTaylor(nx, ny, gamma)\n", - "outfile = makeSimulation(outfile, t_end, rho0, rho_u0, rho_v0, E0, nx, ny, dx, dy, dt, gamma)" + "arguments = InitialConditions.genRayleighTaylor(nx, ny, gamma)\n", + "arguments['context'] = my_context\n", + "outfile = runSimulation(outfile, t_end, arguments)" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 17, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Writing output to c:\\Users\\anbro\\Documents\\projects\\GPUSimulators\\data\\euler_rayleigh-taylor.nc\n", + "Opening c:\\Users\\anbro\\Documents\\projects\\GPUSimulators\\data\\euler_rayleigh-taylor.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": [ + "0% [####==========================] 100%. Total: 37s, elapsed: 5s, remaining: 32s\n", + "0% [########======================] 100%. Total: 37s, elapsed: 10s, remaining: 27s\n", + "0% [############==================] 100%. Total: 41s, elapsed: 15s, remaining: 25s\n", + "0% [###############===============] 100%. Total: 39s, elapsed: 20s, remaining: 19s\n", + "0% [###################===========] 100%. Total: 39s, elapsed: 25s, remaining: 14s\n", + "0% [#######################=======] 100%. Total: 40s, elapsed: 30s, remaining: 9s\n", + "0% [###########################===] 100%. Total: 39s, elapsed: 35s, remaining: 4s\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.nc\n" + ] + } + ], "source": [ "#outfile = 'data/euler_rayleigh-taylor_0007.nc'\n", - "makeAnimation(outfile, vis_type=VisType.Density, vmin=1, vmax=2, cmap=plt.cm.RdBu, save_anim=True, fig_args={'figsize':(3.4, 8), 'tight_layout':True})" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "nx = 1600\n", - "ny = nx//4\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", - "rho0, rho_u0, rho_v0, E0, dx, dy, dt = genShockBubble(nx, ny, gamma)\n", - "outfile = makeSimulation(outfile, t_end, rho0, rho_u0, rho_v0, E0, nx, ny, dx, dy, dt, gamma)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "#outfile = 'data/euler_shock-bubble_0011.nc'\n", - "makeAnimation(outfile, vis_type=VisType.Schlieren, vmin=-0.2, vmax=0.2, cmap=plt.cm.RdBu, save_anim=True, fig_args={'figsize':(16, 5), 'tight_layout':True})" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "nx = 1600\n", - "ny = nx//2\n", - "gamma = 1.4\n", - "roughness = 0.125\n", - "t_end = 10.0\n", - "n_save = 1000\n", - "outfile = \"data/euler_kelvin_helmholtz.nc\"\n", - "outdata = None\n", - "\n", - "rho0, rho_u0, rho_v0, E0, dx, dy, dt = genKelvinHelmholtz(nx, ny, gamma, roughness)\n", - "outfile = makeSimulation(outfile, t_end, rho0, rho_u0, rho_v0, E0, nx, ny, dx, dy, dt, gamma)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "#outfile='data/euler_kelvin_helmholtz_0012.nc'\n", - "makeAnimation(outfile, vis_type=VisType.Density, vmin=1.0, vmax=2.0, save_anim=True, fig_args={'figsize':(16, 9), 'tight_layout':True})" + "createAnimation(outfile, vis_type=VisType.Density, vmin=1, vmax=2, cmap=plt.cm.RdBu, save_anim=False, fig_args={'figsize':(3.4, 8), 'tight_layout':True})" ] }, { diff --git a/GPUSimulators/EE2D_KP07_dimsplit.py b/GPUSimulators/EE2D_KP07_dimsplit.py index def0c05..a7bad02 100644 --- a/GPUSimulators/EE2D_KP07_dimsplit.py +++ b/GPUSimulators/EE2D_KP07_dimsplit.py @@ -61,7 +61,7 @@ class EE2D_KP07_dimsplit (BaseSimulator): gamma, \ theta=1.3, \ order=2, \ - boundaryConditions=BoundaryCondition(), \ + boundary_conditions=BoundaryCondition(), \ block_width=16, block_height=8): # Call super constructor @@ -73,7 +73,7 @@ class EE2D_KP07_dimsplit (BaseSimulator): self.gamma = np.float32(gamma) self.theta = np.float32(theta) self.order = np.int32(order) - self.boundaryConditions = boundaryConditions.asCodedInt() + self.boundary_conditions = boundary_conditions.asCodedInt() #Get kernels self.kernel = context.get_prepared_kernel("cuda/EE2D_KP07_dimsplit.cu", "KP07DimsplitKernel", \ @@ -112,7 +112,7 @@ class EE2D_KP07_dimsplit (BaseSimulator): self.gamma, \ self.theta, \ Simulator.stepOrderToCodedInt(step=0, order=self.order), \ - self.boundaryConditions, \ + self.boundary_conditions, \ self.u0[0].data.gpudata, self.u0[0].data.strides[0], \ self.u0[1].data.gpudata, self.u0[1].data.strides[0], \ self.u0[2].data.gpudata, self.u0[2].data.strides[0], \ @@ -132,8 +132,8 @@ class EE2D_KP07_dimsplit (BaseSimulator): self.g, \ self.gamma, \ self.theta, \ - Simulator.stepOrderToCodedInt(step=0, order=self.order), \ - self.boundaryConditions, \ + Simulator.stepOrderToCodedInt(step=1, order=self.order), \ + self.boundary_conditions, \ self.u0[0].data.gpudata, self.u0[0].data.strides[0], \ self.u0[1].data.gpudata, self.u0[1].data.strides[0], \ self.u0[2].data.gpudata, self.u0[2].data.strides[0], \ diff --git a/GPUSimulators/Simulator.py b/GPUSimulators/Simulator.py index a4702f7..5ae232f 100644 --- a/GPUSimulators/Simulator.py +++ b/GPUSimulators/Simulator.py @@ -66,6 +66,12 @@ class BoundaryCondition(object): self.south = types['south'] self.east = types['east'] self.west = types['west'] + + if (self.north == BoundaryCondition.Type.Neumann \ + or self.south == BoundaryCondition.Type.Neumann \ + or self.east == BoundaryCondition.Type.Neumann \ + or self.west == BoundaryCondition.Type.Neumann): + raise(NotImplementedError("Neumann boundary condition not supported")) def asCodedInt(self): @@ -87,6 +93,9 @@ class BoundaryCondition(object): + + + class BaseSimulator(object): def __init__(self, \ diff --git a/GPUSimulators/cuda/EE2D_KP07_dimsplit.cu b/GPUSimulators/cuda/EE2D_KP07_dimsplit.cu index 487e09b..2e1a4ef 100644 --- a/GPUSimulators/cuda/EE2D_KP07_dimsplit.cu +++ b/GPUSimulators/cuda/EE2D_KP07_dimsplit.cu @@ -156,17 +156,11 @@ __global__ void KP07DimsplitKernel( __shared__ float F[4][h+4][w+4]; //Read into shared memory - readBlock( rho0_ptr_, rho0_pitch_, Q[0], nx_, ny_); - readBlock(rho_u0_ptr_, rho_u0_pitch_, Q[1], nx_, ny_); - readBlock(rho_v0_ptr_, rho_v0_pitch_, Q[2], nx_, ny_); - readBlock( E0_ptr_, E0_pitch_, Q[3], nx_, ny_); + readBlock( rho0_ptr_, rho0_pitch_, Q[0], nx_, ny_, boundary_conditions_); + readBlock(rho_u0_ptr_, rho_u0_pitch_, Q[1], nx_, ny_, boundary_conditions_); + readBlock(rho_v0_ptr_, rho_v0_pitch_, Q[2], nx_, ny_, boundary_conditions_); + readBlock( E0_ptr_, E0_pitch_, Q[3], nx_, ny_, boundary_conditions_); __syncthreads(); - - //Fix boundary conditions - noFlowBoundary(Q[0], nx_, ny_); - noFlowBoundary(Q[1], nx_, ny_); - noFlowBoundary(Q[2], nx_, ny_); - noFlowBoundary(Q[3], nx_, ny_); //Step 0 => evolve x first, then y @@ -185,7 +179,6 @@ __global__ void KP07DimsplitKernel( minmodSlopeY(Q, Qx, theta_); __syncthreads(); - computeFluxG(Q, Qx, F, gamma_, dy_, dt_); __syncthreads(); @@ -236,24 +229,21 @@ __global__ void KP07DimsplitKernel( //This is the RK2-part if (getOrder(step_order_) == 2) { - const int tx = threadIdx.x + gc; - const int ty = threadIdx.y + gc; - const float q1 = Q[0][ty][tx]; - const float q2 = Q[1][ty][tx]; - const float q3 = Q[2][ty][tx]; - const float q4 = Q[3][ty][tx]; - __syncthreads(); + const int i = threadIdx.x + gc; + const int j = threadIdx.y + gc; + const int tx = blockDim.x*blockIdx.x + i; + const int ty = blockDim.y*blockIdx.y + j; - readBlock( rho1_ptr_, rho1_pitch_, Q[0], nx_, ny_); - readBlock(rho_u1_ptr_, rho_u1_pitch_, Q[1], nx_, ny_); - readBlock(rho_v1_ptr_, rho_v1_pitch_, Q[2], nx_, ny_); - readBlock( E1_ptr_, E1_pitch_, Q[3], nx_, ny_); - __syncthreads(); + const float q1 = ((float*) ((char*) rho1_ptr_ + rho1_pitch_*ty))[tx]; + const float q2 = ((float*) ((char*) rho_u1_ptr_ + rho_u1_pitch_*ty))[tx]; + const float q3 = ((float*) ((char*) rho_v1_ptr_ + rho_v1_pitch_*ty))[tx]; + const float q4 = ((float*) ((char*) E1_ptr_ + E1_pitch_*ty))[tx]; - Q[0][ty][tx] = 0.5f*( Q[0][ty][tx] + q1 ); - Q[1][ty][tx] = 0.5f*( Q[1][ty][tx] + q2 ); - Q[2][ty][tx] = 0.5f*( Q[2][ty][tx] + q3 ); - Q[3][ty][tx] = 0.5f*( Q[3][ty][tx] + q4 ); + Q[0][j][i] = 0.5f*( Q[0][j][i] + q1 ); + Q[1][j][i] = 0.5f*( Q[1][j][i] + q2 ); + Q[2][j][i] = 0.5f*( Q[2][j][i] + q3 ); + Q[3][j][i] = 0.5f*( Q[3][j][i] + q4 ); + __syncthreads(); } } diff --git a/GPUSimulators/cuda/common.h b/GPUSimulators/cuda/common.h index 5729ee7..3bf4400 100644 --- a/GPUSimulators/cuda/common.h +++ b/GPUSimulators/cuda/common.h @@ -86,45 +86,135 @@ __device__ float desingularize(float x_, float eps_) { + + + +/** + * Returns the step stored in the leftmost 16 bits + * of the 32 bit step-order integer + */ +inline __device__ int getStep(int step_order_) { + return step_order_ >> 16; +} + +/** + * Returns the order stored in the rightmost 16 bits + * of the 32 bit step-order integer + */ +inline __device__ int getOrder(int step_order_) { + return step_order_ & 0x0000FFFF; +} + + +enum BoundaryCondition { + Dirichlet = 0, + Neumann = 1, + Periodic = 2, + Reflective = 3 +}; + +inline __device__ BoundaryCondition getBCNorth(int bc_) { + return static_cast(bc_ & 0x000F); +} + +inline __device__ BoundaryCondition getBCSouth(int bc_) { + return static_cast((bc_ >> 8) & 0x000F); +} + +inline __device__ BoundaryCondition getBCEast(int bc_) { + return static_cast((bc_ >> 16) & 0x000F); +} + +inline __device__ BoundaryCondition getBCWest(int bc_) { + return static_cast(bc_ >> 24); +} + + + + + + +template +inline __device__ int handlePeriodicBoundaryX(int k, int nx_, int boundary_conditions_) { + const int gc_pad = 2*ghost_cells; + + if ((k < gc_pad) + && getBCWest(boundary_conditions_) == Periodic) { + k += (nx_+2*ghost_cells - 2*gc_pad); + } + else if ((k >= nx_+2*ghost_cells-gc_pad) + && getBCEast(boundary_conditions_) == Periodic) { + k -= (nx_+2*ghost_cells - 2*gc_pad); + } + + return k; +} + +template +inline __device__ int handlePeriodicBoundaryY(int l, int ny_, int boundary_conditions_) { + const int gc_pad = 2*ghost_cells; + + if ((l < gc_pad) + && getBCSouth(boundary_conditions_) == Periodic) { + l += (ny_+2*ghost_cells - 2*gc_pad); + } + else if ((l >= ny_+2*ghost_cells-gc_pad) + && getBCNorth(boundary_conditions_) == Periodic) { + l -= (ny_+2*ghost_cells - 2*gc_pad); + } + + return l; +} + + /** * Reads a block of data with ghost cells */ -template +template inline __device__ void readBlock(float* ptr_, int pitch_, - float shmem[block_height+2*ghost_cells][block_width+2*ghost_cells], - const int nx_, const int ny_) { + float Q[block_height+2*ghost_cells][block_width+2*ghost_cells], + const int nx_, const int ny_, + const int boundary_conditions_) { //Index of block within domain const int bx = blockDim.x * blockIdx.x; const int by = blockDim.y * blockIdx.y; - const int gc_pad = 4; - //Read into shared memory //Loop over all variables for (int j=threadIdx.y; j= ny_+2*ghost_cells-gc_pad) ) * (ny_+2*ghost_cells - 2*gc_pad); - const int l = min(y + y_offset, ny_+2*ghost_cells-1); - */ - + //Handle periodic boundary conditions here + int l = handlePeriodicBoundaryY(by + j, ny_, boundary_conditions_); + l = min(l, ny_+2*ghost_cells-1); float* row = (float*) ((char*) ptr_ + pitch_*l); for (int i=threadIdx.x; i(bx + i, nx_, boundary_conditions_); + k = min(k, nx_+2*ghost_cells-1); - /* - const int x = bx + i; - const int gc_pad = 4; - const int x_offset = ( (int) (x < gc_pad) - (int) (x >= nx_+2*ghost_cells-gc_pad) ) * (nx_+2*ghost_cells - 2*gc_pad); - const int k = min(x + x_offset, nx_+2*ghost_cells-1); - */ - - shmem[j][i] = row[k]; + //Read from global memory + Q[j][i] = row[k]; } } + __syncthreads(); + + //Handle reflective boundary conditions + if (getBCNorth(boundary_conditions_) == Reflective) { + bcNorthReflective(Q, nx_, ny_); + __syncthreads(); + } + if (getBCSouth(boundary_conditions_) == Reflective) { + bcSouthReflective(Q, nx_, ny_); + __syncthreads(); + } + if (getBCEast(boundary_conditions_) == Reflective) { + bcEastReflective(Q, nx_, ny_); + __syncthreads(); + } + if (getBCWest(boundary_conditions_) == Reflective) { + bcWestReflective(Q, nx_, ny_); + __syncthreads(); + } } @@ -165,17 +255,6 @@ inline __device__ void writeBlock(float* ptr_, int pitch_, -template -__device__ void noFlowBoundary(float Q[block_height+2*ghost_cells][block_width+2*ghost_cells], const int nx_, const int ny_) { - bcEastReflective(Q, nx_, ny_); - bcWestReflective(Q, nx_, ny_); - __syncthreads(); - bcNorthReflective(Q, nx_, ny_); - bcSouthReflective(Q, nx_, ny_); - __syncthreads(); -} - - // West boundary template __device__ void bcWestReflective(float Q[block_height+2*ghost_cells][block_width+2*ghost_cells], const int nx_, const int ny_) { @@ -355,50 +434,6 @@ __device__ void memset(float Q[vars][shmem_height][shmem_width], float value) { } -/** - * Returns the step stored in the leftmost 16 bits - * of the 32 bit step-order integer - */ -inline __device__ int getStep(int step_order_) { - return step_order_ >> 16; -} - -/** - * Returns the order stored in the rightmost 16 bits - * of the 32 bit step-order integer - */ -inline __device__ int getOrder(int step_order_) { - return step_order_ & 0x0000FFFF; -} - - -enum BoundaryCondition { - Dirichlet = 0, - Neumann = 1, - Periodic = 2, - Reflective = 3 -}; - -inline __device__ BoundaryCondition getBCNorth(int bc_) { - return static_cast(bc_ & 0x000F); -} - -inline __device__ BoundaryCondition getBCSouth(int bc_) { - return static_cast((bc_ >> 8) & 0x000F); -} - -inline __device__ BoundaryCondition getBCEast(int bc_) { - return static_cast((bc_ >> 16) & 0x000F); -} - -inline __device__ BoundaryCondition getBCWest(int bc_) { - return static_cast(bc_ >> 24); -} - - - - -