{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"#Lets have matplotlib \"inline\"\n",
"%matplotlib inline\n",
"\n",
"# Add line profiler\n",
"%load_ext line_profiler\n",
"\n",
"#Import packages we need\n",
"import numpy as np\n",
"from matplotlib import animation, rc\n",
"from matplotlib import pyplot as plt\n",
"from matplotlib.colorbar import ColorbarBase\n",
"\n",
"from IPython.display import display, HTML\n",
"\n",
"from enum import Enum\n",
"import subprocess\n",
"import time\n",
"import os\n",
"import gc\n",
"import datetime\n",
"import importlib\n",
"import logging\n",
"import netCDF4\n",
"import json\n",
"\n",
"import pycuda.driver as cuda\n",
"import pycuda.compiler\n",
"\n",
"try:\n",
" from StringIO import StringIO\n",
"except ImportError:\n",
" from io import StringIO\n",
"\n",
"from GPUSimulators import Common, IPythonMagic\n",
"from GPUSimulators.EE2D_KP07_dimsplit import EE2D_KP07_dimsplit\n",
"from GPUSimulators.helpers import InitialConditions, Visualization"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Console logger using level INFO\n",
"File logger using level DEBUG to test_schemes.log\n",
"Python version 3.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": 3,
"metadata": {},
"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 <306858187648>\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": 4,
"metadata": {},
"outputs": [],
"source": [
"#Set large figure sizes as default\n",
"rc('figure', figsize=(16.0, 12.0))\n",
"rc('animation', html='html5')\n",
"rc('animation', bitrate=1800)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"def plotVars(rho, rho_u, rho_v, E):\n",
" plt.figure()\n",
" plt.subplot(1,4,1)\n",
" plt.imshow(rho, origin='bottom')\n",
" plt.colorbar(orientation='horizontal', pad=0.02, shrink=0.9)\n",
"\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",
" 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": 6,
"metadata": {},
"outputs": [],
"source": [
"def runSimulation(outfile, t_end, sim_args):\n",
" with Common.Timer(\"construct\") as t:\n",
" sim = EE2D_KP07_dimsplit(**sim_args)\n",
" print(\"Constructed in \" + str(t.secs) + \" seconds\")\n",
"\n",
" #Create netcdf file and simulate\n",
" with Common.DataDumper(outfile, mode='w', clobber=False) as outdata:\n",
" outdata.ncfile.createDimension('time', None)\n",
" outdata.ncfile.createDimension('x', sim.nx)\n",
" outdata.ncfile.createDimension('y', sim.ny)\n",
"\n",
" #Create variables\n",
" outdata.time_var = outdata.ncfile.createVariable('time', np.dtype('float32').char, 'time')\n",
" outdata.x_var = outdata.ncfile.createVariable('x', np.dtype('float32').char, 'x')\n",
" outdata.y_var = outdata.ncfile.createVariable('y', np.dtype('float32').char, 'y')\n",
" outdata.rho_var = outdata.ncfile.createVariable('rho', np.dtype('float32').char, ('time', 'y', 'x'), zlib=True, least_significant_digit=3)\n",
" outdata.rho_u_var = outdata.ncfile.createVariable('rho_u', np.dtype('float32').char, ('time', 'y', 'x'), zlib=True, least_significant_digit=3)\n",
" outdata.rho_v_var = outdata.ncfile.createVariable('rho_v', np.dtype('float32').char, ('time', 'y', 'x'), zlib=True, least_significant_digit=3)\n",
" outdata.E_var = outdata.ncfile.createVariable('E', np.dtype('float32').char, ('time', 'y', 'x'), zlib=True, least_significant_digit=3)\n",
" \n",
" #Create attributes\n",
" def toJson(in_dict):\n",
" out_dict = in_dict.copy()\n",
"\n",
" for key in out_dict:\n",
" if isinstance(out_dict[key], np.ndarray):\n",
" out_dict[key] = out_dict[key].tolist()\n",
" else:\n",
" try:\n",
" json.dumps(out_dict[key])\n",
" except:\n",
" out_dict[key] = str(out_dict[key])\n",
"\n",
" return json.dumps(out_dict)\n",
" outdata.ncfile.created = time.ctime(time.time())\n",
" outdata.ncfile.sim_args = toJson(sim_args)\n",
"\n",
" outdata.x_var[:] = np.linspace(0, sim.nx*sim.dx, sim.nx)\n",
" outdata.y_var[:] = np.linspace(0, sim.ny*sim.dy, sim.ny)\n",
"\n",
" progress_printer = Common.ProgressPrinter(n_save, print_every=10)\n",
" print(\"Simulating to t={:f}. Estimated {:d} timesteps (dt={:f})\".format(t_end, int(t_end / sim.dt), sim.dt))\n",
" for i in range(n_save+1):\n",
" #Sanity check simulator\n",
" try:\n",
" sim.check()\n",
" except AssertionError as e:\n",
" print(\"Error after {:d} steps (t={:f}: {:s}\".format(sim.simSteps(), sim.simTime(), str(e)))\n",
" return outdata.filename\n",
"\n",
" #Simulate\n",
" if (i > 0):\n",
" sim.simulate(t_end/n_save)\n",
"\n",
" #Save to file\n",
" #rho = sim.u0[0].download(sim.stream)\n",
" rho, rho_u, rho_v, E = sim.download()\n",
" outdata.time_var[i] = sim.simTime()\n",
" outdata.rho_var[i, :] = rho\n",
" outdata.rho_u_var[i, :] = rho_u\n",
" outdata.rho_v_var[i, :] = rho_v\n",
" outdata.E_var[i, :] = E\n",
"\n",
" #Write progress to screen\n",
" print_string = progress_printer.getPrintString(i)\n",
" if (print_string):\n",
" print(print_string)\n",
"\n",
" return outdata.filename "
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"class VisType(Enum):\n",
" Schlieren = 0\n",
" Density = 1\n",
"\n",
"def createAnimation(infile, vis_type, vmin, vmax, save_anim=False, cmap=plt.cm.coolwarm, fig_args=None):\n",
" fig = plt.figure(**fig_args)\n",
" \n",
" with Common.DataDumper(infile, 'r') as indata:\n",
" time = indata.ncfile.variables['time'][:]\n",
" x = indata.ncfile.variables['x'][:]\n",
" y = indata.ncfile.variables['y'][:]\n",
" rho = indata.ncfile.variables['rho'][0]\n",
" rho_u = indata.ncfile.variables['rho_u'][0]\n",
" rho_v = indata.ncfile.variables['rho_v'][0]\n",
" \n",
" created = indata.ncfile.created\n",
" sim_args = json.loads(indata.ncfile.sim_args)\n",
" for key in sim_args:\n",
" if isinstance(sim_args[key], list):\n",
" sim_args[key] = \"[...]\"\n",
" num_frames = len(time)\n",
" print(\"{:s} created {:s} contains {:d} timesteps\".format(infile, created, num_frames))\n",
" print(\"Simulator arguments: \\n\", sim_args)\n",
"\n",
" ax1 = fig.gca()\n",
" extents = [0, x.max(), 0, y.max()]\n",
" \n",
" if (vis_type == VisType.Schlieren):\n",
" im = ax1.imshow(Visualization.genColors(rho, rho_u, rho_v, cmap, vmax, vmin), origin='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",
" fig.suptitle(\"Density at t={:.2f}\".format(time[0]))\n",
" else:\n",
" assert False, \"Wrong vis_type\"\n",
"\n",
" #Create colorbar\n",
" from matplotlib.colors import Normalize\n",
" from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
" divider = make_axes_locatable(ax1)\n",
" ax2 = divider.append_axes(\"right\", size=0.1, pad=0.05)\n",
" cb1 = ColorbarBase(ax2, cmap=cmap,\n",
" norm=Normalize(vmin=vmin, vmax=vmax),\n",
" orientation='vertical')\n",
" \n",
" #Label colorbar\n",
" if (vis_type == VisType.Schlieren):\n",
" cb1.set_label('Vorticity')\n",
" elif (vis_type == VisType.Density):\n",
" cb1.set_label('Density')\n",
"\n",
" progress_printer = Common.ProgressPrinter(num_frames, print_every=5)\n",
" \n",
" def animate(i):\n",
" rho = indata.ncfile.variables['rho'][i]\n",
" rho_u = indata.ncfile.variables['rho_u'][i]\n",
" rho_v = indata.ncfile.variables['rho_v'][i]\n",
" \n",
" if (vis_type == VisType.Schlieren):\n",
" im.set_data(Visualization.genColors(rho, rho_u, rho_v, cmap, vmax, vmin))\n",
" fig.suptitle(\"Schlieren / vorticity at t={:.2f}\".format(time[i]))\n",
" elif (vis_type == VisType.Density):\n",
" im.set_data(rho) \n",
" fig.suptitle(\"Density at t={:.2f}\".format(time[i]))\n",
" \n",
" #Write progress to screen\n",
" print_string = progress_printer.getPrintString(i)\n",
" if (print_string):\n",
" print(print_string)\n",
"\n",
" anim = animation.FuncAnimation(fig, animate, interval=50, frames=range(num_frames))\n",
" plt.close()\n",
"\n",
" if (save_anim):\n",
" root, _ = os.path.splitext(infile)\n",
" movie_outpath = os.path.abspath(root + \".mp4\")\n",
" if (os.path.isfile(movie_outpath)):\n",
" print(\"Reusing previously created movie \" + movie_outpath)\n",
" else:\n",
" print(\"Creating movie \" + movie_outpath)\n",
" #from matplotlib.animation import FFMpegFileWriter\n",
" #writer = FFMpegFileWriter(fps=25)\n",
" from matplotlib.animation import FFMpegWriter\n",
" writer = FFMpegWriter(fps=25)\n",
" anim.save(movie_outpath, writer=writer)\n",
" display(HTML(\"\"\"\n",
" \n",
" \n",
"
\n",
" \"\"\".format(movie_outpath)))\n",
" else:\n",
" #plt.rcParams[\"animation.html\"] = \"html5\"\n",
" plt.rcParams[\"animation.html\"] = \"jshtml\"\n",
" display(anim)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Shock-bubble"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "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_0006.nc\n",
"Opening c:\\Users\\anbro\\Documents\\projects\\GPUSimulators\\data\\euler_shock-bubble_0006.nc\n",
"Keyword arguments: {'mode': 'w', 'clobber': False}\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Constructed in 1.875434160232544 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_0006.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_0006.nc\n",
"Opening c:\\Users\\anbro\\Documents\\projects\\GPUSimulators\\data\\euler_shock-bubble_0006.nc\n",
"Arguments: ('r',)\n",
"Closing c:\\Users\\anbro\\Documents\\projects\\GPUSimulators\\data\\euler_shock-bubble_0006.nc\n",
"Exception caught: Resetting to CUDA context my_context\n",
"Traceback (most recent call last):\n",
" File \"c:\\Users\\anbro\\AppData\\Local\\Continuum\\anaconda3\\lib\\site-packages\\IPython\\core\\interactiveshell.py\", line 3267, in run_code\n",
" exec(code_obj, self.user_global_ns, self.user_ns)\n",
" File \"\", line 2, in \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)})\n",
" File \"\", line 18, in createAnimation\n",
" for key in sim_args:\n",
"RuntimeError: dictionary changed size during iteration\n",
"Popping <306858187648>\n",
"Pushing <306858187648>\n"
]
},
{
"ename": "RuntimeError",
"evalue": "dictionary changed size during iteration",
"output_type": "error",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[1;31mRuntimeError\u001b[0m Traceback (most recent call last)",
"\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m\u001b[0m\n\u001b[0;32m 1\u001b[0m \u001b[1;31m#outfile = 'data/euler_shock-bubble_0008.nc'\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[0mcreateAnimation\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0moutfile\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mvis_type\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mVisType\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mSchlieren\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mvmin\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m0.2\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mvmax\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m0.2\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcmap\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcm\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mRdBu\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msave_anim\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mFalse\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfig_args\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m{\u001b[0m\u001b[1;34m'figsize'\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m16\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m5\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m}\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[1;32m\u001b[0m in \u001b[0;36mcreateAnimation\u001b[1;34m(infile, vis_type, vmin, vmax, save_anim, cmap, fig_args)\u001b[0m\n\u001b[0;32m 16\u001b[0m \u001b[0mcreated\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mindata\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mncfile\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcreated\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 17\u001b[0m \u001b[0msim_args\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mjson\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mloads\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mindata\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mncfile\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msim_args\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 18\u001b[1;33m \u001b[1;32mfor\u001b[0m \u001b[0mkey\u001b[0m \u001b[1;32min\u001b[0m \u001b[0msim_args\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 19\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msim_args\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mkey\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mlist\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 20\u001b[0m \u001b[0msim_args\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpop\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;31mRuntimeError\u001b[0m: dictionary changed size during iteration"
]
},
{
"data": {
"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": null,
"metadata": {},
"outputs": [],
"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": null,
"metadata": {},
"outputs": [],
"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": null,
"metadata": {},
"outputs": [],
"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": null,
"metadata": {},
"outputs": [],
"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": null,
"metadata": {},
"outputs": [],
"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": null,
"metadata": {},
"outputs": [],
"source": [
"#outfile = 'data/euler_rayleigh-taylor_0007.nc'\n",
"createAnimation(outfile, vis_type=VisType.Density, vmin=1, vmax=2, cmap=plt.cm.coolwarm, save_anim=False, fig_args={'figsize':(3.4, 8)})"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python [conda env:anaconda3]",
"language": "python",
"name": "conda-env-anaconda3-py"
},
"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.6.7"
}
},
"nbformat": 4,
"nbformat_minor": 2
}