From 5b3aa6b7013c8ad708514efce8140072d3e4e5c9 Mon Sep 17 00:00:00 2001 From: Ally Smith <89595486+AllySmith03@users.noreply.github.com> Date: Thu, 30 Apr 2026 03:32:04 -0600 Subject: [PATCH] Jupyter Notebook for ramp A Jupyter notebook containing an example of the ramp problem for the compressible solver --- examples/ramp.ipynb | 343 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 343 insertions(+) create mode 100644 examples/ramp.ipynb diff --git a/examples/ramp.ipynb b/examples/ramp.ipynb new file mode 100644 index 000000000..77346a7c4 --- /dev/null +++ b/examples/ramp.ipynb @@ -0,0 +1,343 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "02953843", + "metadata": {}, + "source": [ + "# Investigating the Ramp problem" + ] + }, + { + "cell_type": "markdown", + "id": "c241bf01", + "metadata": {}, + "source": [ + "Pyro's compressible solver supports the problem ```ramp```, which is based on the Woodward & Colella 1984 test problem Double Mach Reflection of a Strong Shock [1]. This is an extension of the classic shock tube test for the accuracy of hydrodynamics solvers, where a shock, modeled by a high pressure, density, and velocity fluid on the left hand side, collides at an angle with a stationary normal pressure and density fluid on the right. This setup results in discontinuities, i.e. shocks and contact discontinuities. \n", + "\n", + "This problem is valuable for analyzing the effictiveness of how a solver handles discontinuities, since both shocks and contact discontinuities occur frequently in compressible hydrodynamics. \n", + "\n", + "Note: This is not a rigourous investigation of Pyro's performance on this problem. Instead we focus on how this problem is run as well as the purpose and value of running this test. " + ] + }, + { + "cell_type": "markdown", + "id": "a5cfdd70", + "metadata": {}, + "source": [ + "## Setup\n", + "\n", + "Since ```inputs.ramp``` is already set up for the Woodward and Colella [1] Double Mach Reflection problem, we can simply import it. This includes all of the specifications for the mesh, boundary conditions, and initial conditions like density, x velocity, y velocity, and pressure for both the left and right states. " + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "2b9d9fb1", + "metadata": {}, + "outputs": [], + "source": [ + "from pyro import Pyro\n", + "\n", + "# Select solver, problem, and input file\n", + "solver = \"compressible\"\n", + "problem_name = \"ramp\"\n", + "param_file = \"inputs.ramp\"\n", + "\n", + "\n", + "# Initialize object\n", + "pyro_sim = Pyro(solver)\n", + "pyro_sim.initialize_problem(problem_name, inputs_file=param_file)" + ] + }, + { + "cell_type": "markdown", + "id": "b88819a0", + "metadata": {}, + "source": [ + "To see these values, we can print the state of the object." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "ffa33991", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Solver = compressible\n", + "Problem = ramp\n", + "Simulation time = 0.0\n", + "Simulation step number = 0\n", + "\n", + "Runtime Parameters\n", + "------------------\n", + "compressible.cvisc = 0.1\n", + "compressible.delta = 0.33\n", + "compressible.grav = 0.0\n", + "compressible.limiter = 2\n", + "compressible.riemann = HLLC\n", + "compressible.small_dens = -1e+200\n", + "compressible.small_eint = -1e+200\n", + "compressible.use_flattening = 1\n", + "compressible.z0 = 0.75\n", + "compressible.z1 = 0.85\n", + "driver.cfl = 0.8\n", + "driver.fix_dt = -1.0\n", + "driver.init_tstep_factor = 0.01\n", + "driver.max_dt_change = 2.0\n", + "driver.max_steps = 12500\n", + "driver.tmax = 0.25\n", + "driver.verbose = 0\n", + "eos.gamma = 1.4\n", + "io.basename = double_mach_reflection_\n", + "io.do_io = 0\n", + "io.dt_out = 0.1\n", + "io.force_final_output = 0\n", + "io.n_out = 10000\n", + "mesh.grid_type = Cartesian2d\n", + "mesh.nx = 1024\n", + "mesh.ny = 256\n", + "mesh.xlboundary = ramp\n", + "mesh.xmax = 4.0\n", + "mesh.xmin = 0.0\n", + "mesh.xrboundary = outflow\n", + "mesh.ylboundary = ramp\n", + "mesh.ymax = 1.0\n", + "mesh.ymin = 0.0\n", + "mesh.yrboundary = ramp\n", + "particles.do_particles = 0\n", + "particles.n_particles = 100\n", + "particles.particle_generator = grid\n", + "ramp.pl = 116.5\n", + "ramp.pr = 1.0\n", + "ramp.rhol = 8.0\n", + "ramp.rhor = 1.4\n", + "ramp.ul = 7.1447096\n", + "ramp.ur = 0.0\n", + "ramp.vl = -4.125\n", + "ramp.vr = 0.0\n", + "sponge.do_sponge = 0\n", + "sponge.sponge_rho_begin = 0.01\n", + "sponge.sponge_rho_full = 0.001\n", + "sponge.sponge_timescale = 0.01\n", + "vis.dovis = 0\n", + "vis.store_images = 0\n", + "\n" + ] + } + ], + "source": [ + "print(pyro_sim)" + ] + }, + { + "cell_type": "markdown", + "id": "770b7b25", + "metadata": {}, + "source": [ + "For more of a visual, we can plot the initial data. Here we can see the \"at rest\" right hand side of the shock tube as well as the shock on the left hand side hitting this stationary fluid. " + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "f8a483c5", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAHiCAYAAAAOOXDQAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAANwNJREFUeJzt3Xt0VFWeNv5nn7qFWwoSIEWGJKQBbRCMGHg1+SHGC2GCMIO6ENFJg0L3ojswQnR6QHrkMg5hiQP0jA02dr+K2L1EbkrbjCaOSKRpOhDBRryFl2gwVjomdi7QpipVZ//+qEpBkQupIsnO4TyftWqdc3adU3yrEurJPnvXKSGllCAiIlPSVBdARETqMASIiEyMIUBEZGIMASIiE2MIEBGZGEOAiMjEGAJERCbGECAiMjGGABGRiTEEiIhMjCFARGRiDAEiMrSdO3ciLS0Nffv2Rd++fTFjxgy43W7VZRkGQ4CIDOsnP/kJFixYgDlz5uCNN97AypUr8dZbb2HevHmqSzMMwauIEpER/eY3v8EPfvADHDp0CJMnTw61/+AHP8Arr7yCv/71r3A6nQorNAb2BIjIkP7jP/4D9957b1gAAMB1110HKSX+9re/KarMWBgCRGQ4n376KT755BPMnDmz1X1fffUVBgwYgISEBAWVGQ9DgIgM58iRIwCA5OTksHZd1/Hmm29i1qxZ0DS+vXUGXyUiMpxjx44BAMrKysLan332WfzlL3/BkiVLVJRlSFbVBRARRaqkpARJSUlYuXIl7HY7EhISsH//fvzyl7/Ehg0bMGnSJNUlGgZnBxGRoXg8HgwYMABPPvkkBg0ahGeffRbV1dW44YYbsHz5cjzwwAOqSzQUhgARGcqf/vQn3HrrrXjzzTdxzz33qC7H8DgmQESG0jIeMHHiRMWVXBsYAkRkKMeOHUNSUhKngHYRng4iIjIx9gSIiEyMIUBEZGIMASIiE2MIEBGZGEOAiMjEGAJEJlNQUIBJkyZhwIABGDp0KGbNmoXPPvusw2Pee+89CCFa3T799NMeqjrc6tWrW9Xicrk6PObQoUNIT09HTEwMvve97+H555/voWrbNmLEiDZf07y8vDb3766fAa8dRGQyhw4dQl5eHiZNmgSfz4eVK1ciOzsbH3/8Mfr169fhsZ999hliY2ND20OGDOnuctt1ww034J133gltWyyWdvctLy/H9OnT8cMf/hCvvPIK/vCHP+AnP/kJhgwZgvvvv78nym3l2LFj8Pv9oe2PPvoIU6dOxezZszs8rqt/BgwBIpN56623wrZffPFFDB06FKWlpZgyZUqHxw4dOhQDBw7sxuo6z2q1XvGv/xbPP/88kpOTsXnzZgDAmDFjcPz4cTz77LPKQuDyN+/169dj5MiRuP322zs8rqt/BjwdRGRy9fX1AIC4uLgr7jthwgQMGzYMd911Fw4ePNjdpXWorKwMiYmJSE1NxYMPPoizZ8+2u+8f//hHZGdnh7VNmzYNx48fR3Nzc3eXekVerxevvPIKHn30UQghOty3q38GDAEiE5NSIj8/H5MnT8a4cePa3W/YsGHYtm0b9uzZg7179+L666/HXXfdheLi4h6s9qJbbrkFL7/8Mt5++2288MILqKqqQmZmJmpra9vcv6qqqtVlJhISEuDz+VBTU9MTJXfo9ddfR11dHebPn9/uPt31M+BlI4hMLC8vD7///e9x+PBhDB8+PKJjZ86cCSEE9u/f303Vdd6FCxcwcuRI/PSnP0V+fn6r+6+77jo88sgjWLFiRajtD3/4AyZPngy3293p00rdZdq0abDb7fjd734X0XFd8TNgT4DIpJYsWYL9+/fj4MGDEQcAANx6662tvtlLlX79+mH8+PHt1uNyuVBVVRXWVl1dDavVivj4+J4osV1ffvkl3nnnHSxcuDDiY7viZ8AQIDIZKSUWL16MvXv34t1330VqampUj3PixAkMGzasi6uLjsfjwSeffNJuPRkZGSgqKgprKywsxMSJE2Gz2XqixHa1DMxH890IXfIzkERkKj/+8Y+l0+mU7733nnS73aHb3/72t9A+y5cvl7m5uaHtTZs2yX379snPP/9cfvTRR3L58uUSgNyzZ4+KpyAff/xx+d5778mzZ8/Ko0ePyhkzZsgBAwbIL774os36z549K/v27SuXLVsmP/74Y/nrX/9a2mw2uXv3biX1t/D7/TI5OVn+67/+a6v7eupnwCmiRCazdetWAEBWVlZY+4svvhgamHS73aioqAjd5/V68cQTT6CyshJ9+vTBDTfcgN///veYPn16T5Ud5quvvsLcuXNRU1ODIUOG4NZbb8XRo0eRkpLSZv2pqak4cOAAli1bhl/84hdITEzEf/3XfymbHtrinXfeQUVFBR599NFW9/XUz4ADw0REJsYxASIiE2MIEBGZGEOAiMjEGAJERCbGECAiMjGGABGRiTEEiIhMjCFARK14PB6sXr0aHo9HdSlR43PoHH5YjIhaaWhogNPpRH19fdi3WBkJn0PnsCdARGRiDAEiIhPjBeSIDK6pqQler7dLH7OhoSFsaURmew52ux0xMTER/xscEyAysKamJjj7DIIXTapLIcVcLhfKy8sjDgL2BIgMzOv1wosmTLHPgtXigLBYAIsF0LTAuqYFtwVgCa4LDbBokJoWOCGsBdctAlITgBZYyku3g+tSE5ACwW0E24JLCyA1AKKlDWH74NI2cXEZ1t7GfhDh7dBk2L4QEhIX70NoHwnR0iYkhCYhRKBNaHpgW5PQhAy8TJoOTZOwaDoskNA0HRZNh1UEl5oOi5CwBtetQodV88MCHTbND6vwwyYkLJofNhG8hdp12IQfVuG7eJ/wwQYdNuGDRUjYhA9W+GEP3m8VOuxoOVbCCgm7ELAJASs02IQGKyywCQsaGnWkpH8Br9fLECAyIytssAo7hLAAIhgCWjAEQstL1tsNgYth0G4IaB2HwMW2DkJA62QIaAiFwMX72wgB0ckQ0C6GgBYWAjI8BERweVkIBN74w0PAKvSwN3urFniDtgkNttC6HzYhYBMILgPr9pbgEBI2IWEDYBcIrAvADhk8RsJ2SQjYQiGgwSYsV/W7w4FhIiITYwgQEZkYQ4CIyMQYAkREJsaBYSKFiouLsWHDBpSWlsLtdmPfvn2YNWtWxI/jQzMgNQhpAaQF0LXAIDG0wEAxRGBWkLCERlSl1NAyrUYGR2ClFIAMDgDLSwaG5WWzg/Q2BoZ1Y80OksGBYSlk8NhAGzQdgITUdEhNB0RwqenBfQPbUuiQmh86gkvhhxQS/uC6LvzQNT/8wg+/0OEXfviEDz7hD9588CEwaygwMOyHFX40twwwCx126O3MDgJsArBCwCb8aGjUo/4dZAgQKXThwgWkpaXhkUcewf333x/x8Xa7HS6XC8VVr3d9cWQoLpcLdrs94uP4YTGiXkIIEVVPoDs+MUzGE+0nhtkTIDIQj8cTdllhXdfx7bffIj4+HkIIhZWRSlJK1NTUIDExEZoW2VAvQ4DIQAoKCrBmzRrVZVAvde7cOQwfPjyiY3g6iKiX6MzpoMt7AvX19UhOTsaXH4xAbP+un+x3094FSF1e0uWPS13Lh2YcxgHU1dXB6XRGdCx7AkQG4nA44HA4WrXH9tcQO+DqLh/QltErTgDC1uWPS10s+Kd8NKcE+TkBImrTtMQ01SVQD2BPgEih8+fP48yZM6Ht8vJynDx5EnFxcUhOTlZYGZkFQ4BIoePHj+OOO+4Ibefn5wMA5s2bh5deeklRVewFmAlDgEihrKws9La5GfO+nAKgXnUZ1EM4JkBEYaoyGABmwhAgopBROxepLoF6GEOAiEJGLvuj6hKohzEEiAgAMCXvR6pLIAUYAkQEAOiz70+qSyAFGAJExCmhJsYQIDK5wJRQMiuGAFEU5s+fj+LiYtVldAlOCTU3hgBRFBobG5GdnY3Ro0dj3bp1qKysVF1SVDgYTAwBoijs2bMHlZWVWLx4MXbt2oURI0YgJycHu3fvRnNzs+ryOo2DwcQQIIpSfHw8HnvsMZw4cQIlJSUYNWoUcnNzkZiYiGXLlqGsrEx1iR3iYDABDAGiq+Z2u1FYWIjCwkJYLBZMnz4dp0+fxtixY7Fp0ybV5bWJg8HUgiFAFIXm5mbs2bMHM2bMQEpKCnbt2oVly5bB7XZj+/btKCwsxI4dO7B27VrVpbaJg8HUglcRJYrCsGHDoOs65s6di5KSEtx0002t9pk2bRoGDhzY47VdyaidizASvDwEBTAEiKKwadMmzJ49GzExMe3uM2jQIJSXl/dgVZ3D6wPRpXg6iCgKubm5HQZAb8UpoXQ5hgCRiXBKKF2OIUBkEpwSSm1hCBCZAKeEUnsYAkQmwCmh1B6GANE1jl8ZSR1hCBD1Alu2bEFqaipiYmKQnp6O999/v8sem1NCqSMMASLFdu7ciaVLl2LlypU4ceIEbrvtNuTk5KCiouKqH5uDwXQlDAEixTZu3IgFCxZg4cKFGDNmDDZv3oykpCRs3bpVdWlkAvzEMJFCXq8XpaWlWL58eVh7dnY2jhw50mp/j8cDj8cT2q6vDwz4NpzXW+1773XjARjnstYUPV/w5yyljPhYhgCRQjU1NfD7/UhISAhrT0hIQFVVVav9CwoKsGbNmlbtKTd/0cajn+2iKskoamtr4XQ6IzqGIUDUCwghwrallK3aAGDFihXIz88PbdfV1SElJQUVFRUR/+e/VjQ0NCApKQnnzp1DbGys6nKUqK+vR3JyMuLi4iI+liFApNDgwYNhsVha/dVfXV3dqncAAA6HAw6Ho1W70+k07Rtgi9jYWNO/BpoW+TAvB4aJFLLb7UhPT0dRUVFYe1FRETIzMxVVRWbCngCRYvn5+cjNzcXEiRORkZGBbdu2oaKiAosW8UNe1P0YAkSKzZkzB7W1tVi7di3cbjfGjRuHAwcOICUl5YrHOhwOrFq1qs1TRGbB1+DqXgMho5lTRERE1wSOCRARmRhDgIjIxBgCREQmxhAgIjIxhgARkYkxBIiITIwhQERkYgwBIiITYwgQEZkYQ4CIyMQYAkREJsYQICIyMYYAEZGJMQSIiEyMIUBEZGIMASIylNWrV0MIgZqamjbvHzduHLKysnq2KANjCBARmRhDgIjIxBgCREQmxhAgIjIxhgARkYkxBIiITIwhQESGYrVaAQB+v7/N+30+H2w2W0+WZGgMASIylISEBABAZWVlq/uklHC73aF96MoYAkRkKHfeeSeEENi5c2er+9566y00NDTg7rvvVlCZMVlVF0BEFImRI0di8eLF2LBhA+rq6jB9+nT06dMHx44dw/r16zFx4kQ89NBDqss0DCGllKqLICKKhJQSv/zlL/HrX/8aH3/8MXw+H1JSUnDffffhZz/7Gfr376+6RMNgCBARmRjHBIiITIwhQERkYgwBIiITYwgQEZkYQ4CIyMQYAkREJsYQICIyMYYAEZGJMQSIDGzLli1ITU1FTEwM0tPT8f7776suqVdr+ZL6S28ul0t1WUoxBIgMaufOnVi6dClWrlyJEydO4LbbbkNOTg4qKipUl9ar3XDDDXC73aHbqVOnVJekFEOAyKA2btyIBQsWYOHChRgzZgw2b96MpKQkbN26VXVpvZrVaoXL5QrdhgwZorokpRgCRAbk9XpRWlqK7OzssPbs7GwcOXJEUVXGUFZWhsTERKSmpuLBBx/E2bNnVZekFEOAyIBqamrg9/tbfXlKQkICqqqqFFXV+91yyy14+eWX8fbbb+OFF15AVVUVMjMzUVtbq7o0Zfh9AkQGJoQI25ZStmqji3JyckLr48ePR0ZGBkaOHInt27cjPz9fYWXqsCdAZECDBw+GxWJp9Vd/dXU1v1oxAv369cP48eNRVlamuhRlGAJEBmS325Geno6ioqKw9qKiImRmZiqqyng8Hg8++eQTDBs2THUpyvB0EJFB5efnIzc3FxMnTkRGRga2bduGiooKLFq0SHVpvdYTTzyBmTNnIjk5GdXV1Xj66afR0NCAefPmqS5NGYYAkUHNmTMHtbW1WLt2LdxuN8aNG4cDBw4gJSVFdWm91ldffYW5c+eipqYGQ4YMwa233oqjR4+a+jXj10sSEZkYxwSIiEyMIUBEZGIMASIiE2MIEBGZGEOAiMjEGAJERCbGECAyMI/Hg9WrV8Pj8aguxXD42gXwcwJEBtbQ0ACn04n6+nrExsaqLsdQ+NoFsCdARGRiDAEiIhPjtYOIrgENDQ2qSzCcltfsWnnt7HY7YmJiIj6OYwJEBlZfX4/4gYPhh091KaSYy+VCeXl5xEHAngCRgQkh4IcPU+yzYLU4ICwWwGIBNC2wrmnBbQFYgutCAywapKYFTghrwXWLgNQEoAWW8tLt4LrUBKRAcBvBtuDSAkgNgGhpQ9g+uLRNXFyGtbexH0R4OzQZti+EhMTF+xDaR0K0tAkJoUkIEWgTmh7Y1iQ0IQMvk6ZD0yQsmg4LJDRNh0XTYRXBpabDIiSswXWr0GHV/LBAh03zwyr8sAkJi+aHTQRvoXYdNuGHVfgu3id8sEGHJnTYg9tW+IPrfliFDjtajpWwQsIuBGxCwAoNNqHBCgtswoKGRh0p6V/A6/UyBIjMyAobrMIOISyACIaAFgyB0PKS9XZD4GIYtBsCWschcLGtgxDQOhkCGkIhcPH+NkJAdDIEtIshoIWFgAwPARFcXhYCgTf+8BCwCj3szd6qBd6gbUKDLbTuh00I2ASCy8C6XfihCQG7AGxCwgZcXBeAHTJ4jITtkhCwhUJAg01Yrup3hwPDREQmxhAgIjIxhgARkYkxBIiITIwDw0TXAB+aAalBSAsgLYCuBQaJoQUGiiECs4KEJTSiKqWGlmk1MjgCK6UAZHAAWF4yMCwvmx2ktzEwrBtrdpAMDgxLIYPHBtqg6QAkpKZDajoggktND+4b2JZCh9T80BFcCj+kkPAH13Xhh6754Rd++IUOv/DDJ3zwCX/w5oMvODvI1zIjCH40twwwCx126O3MDgJsArBCwCb8aGjUo/7dYQgQKVRcXIwNGzagtLQUbrcb+/btw6xZszp9vN1uh8vlQnHV691WIxmDy+WC3W6P+DiGAJFCFy5cQFpaGh555BHcf//9ER8fExOD8vJyeL3ebqiOjCTaTwwzBIgUysnJQU5OzlU9RkxMTFT/+YkAhgCRoXg8nrDr3+u6jm+//Rbx8fEQQiisjFSSUqKxsRGJiYnQtMjm+zAEiAykoKAAa9asUV0G9VLnzp3D8OHDIzqGF5Aj6iWEEFccGL68J1BfX4/k5GR8+cEIxPbv+hnf9143vssfk7qeD804jAOoq6uD0+mM6Fj2BIgMxOFwwOFwtGqP7a8hdsDVXUOmLVZh6/LHpG4Q/FM+mlOC/LAYEbVpWmKa6hKoB7AnQKTQ+fPncebMmdB2eXk5Tp48ibi4OCQnJyur6/82JCj7t6lnMQSIFDp+/DjuuOOO0HZ+fj4AYN68eXjppZcUVQXs/L5L2b9NPYshQKRQVlYWetvcjIkfPIB4fKa6DOohHBMgojDxMxgAZsIQIKIQDgabD0OAiMjEGAJEBIC9ALNiCBAR/p/vvOoSSBGGABHhJ8n/n+oSSBGGAFEU5s+fj+LiYtVldImJHzygugRSiCFAFIXGxkZkZ2dj9OjRWLduHSorK1WXFDVOCTU3hgBRFPbs2YPKykosXrwYu3btwogRI5CTk4Pdu3ejublZdXmdxsFgYggQRSk+Ph6PPfYYTpw4gZKSEowaNQq5ublITEzEsmXLUFZWprrEDvH6QAQwBIiumtvtRmFhIQoLC2GxWDB9+nScPn0aY8eOxaZNm1SX1y5eH4gAhgBRVJqbm7Fnzx7MmDEDKSkp2LVrF5YtWwa3243t27ejsLAQO3bswNq1a1WX2ib2AqgFLyBHFIVhw4ZB13XMnTsXJSUluOmmm1rtM23aNAwcOLDHa+sM9gKoBUOAKAqbNm3C7NmzERMT0+4+gwYNQnl5eQ9W1TkcDKZLMQSIopCbm6u6BKIuwTEBIhNhL4AuxxAgMgleH4jawhAgMgleH4jawhAgMgFeH4jawxAgMgFeH4jawxAgusZxMJg6whAg6gW2bNmC1NRUxMTEID09He+//77qksgkGAJEiu3cuRNLly7FypUrceLECdx2223IyclBRUXFVT82ewF0JQwBIsU2btyIBQsWYOHChRgzZgw2b96MpKQkbN269aoel1NCqTP4iWEihbxeL0pLS7F8+fKw9uzsbBw5cqTV/h6PBx6PJ7RdX18PAGg4r7fa90fX/R8AxvluA4qeL/hzllJGfCxDgEihmpoa+P1+JCSEX9UzISEBVVVVrfYvKCjAmjVrWrWn3PxFG49+touqJKOora2F0+mM6BiGAFEvIIQI25ZStmoDgBUrViA/Pz+0XVdXh5SUFFRUVET8n/9a0dDQgKSkJJw7dw6xsbGqy1Givr4eycnJiIuLi/hYhgCRQoMHD4bFYmn1V391dXWr3gEAOBwOOByOVu1Op9O0b4AtYmNjTf8aaFrkw7wcGCZSyG63Iz09HUVFRWHtRUVFyMzMVFQVmQl7AkSK5efnIzc3FxMnTkRGRga2bduGiooKLFq0SHVpZAIMASLF5syZg9raWqxduxZutxvjxo3DgQMHkJKScsVjHQ4HVq1a1eYpIrPga3B1r4GQ0cwpIiKiawLHBIiITIwhQERkYgwBIiITYwgQEZkYQ4CIyMQYAkREJsYQICIyMYYAEZGJMQSIiEyMIUBEZGIMASIiE2MIEBGZGEOAiMjEGAJERCbGECAiMjGGABGRiTEEiIhMjCFARIazevVqCCFw4sQJ3HfffYiNjYXT6cQ//dM/4ZtvvlFdnqEwBIjIsO69916MGjUKu3fvxurVq/H6669j2rRpaG5uVl2aYfCL5onIsO677z4888wzAIDs7GwkJCTg4YcfxmuvvYaHH35YcXXGwJ4AERnW5W/0DzzwAKxWKw4ePKioIuNhCBCRYblcrrBtq9WK+Ph41NbWKqrIeBgCRGRYVVVVYds+nw+1tbWIj49XVJHxMASIyLB+85vfhG2/9tpr8Pl8yMrKUlOQAXFgmIgMa+/evbBarZg6dSpOnz6Nf/u3f0NaWhoeeOAB1aUZBnsCRGRYe/fuxaeffor77rsPTz31FGbOnInCwkLY7XbVpRkGewJEZFjJycnYv3+/6jIMjT0BIiITYwgQEZmYkFJK1UUQEZEa7AkQEZkYQ4CIyMQYAkREJsYQICIyMYYAEZGJMQSIiEyMIUBEyrR8TeSlt0svDy2lxOrVq5GYmIg+ffogKysLp0+fVlhxx4qLizFz5kwkJiZCCIHXX3897P7OPB+Px4MlS5Zg8ODB6NevH/7hH/4BX331VbfVzBAgIqVuuOEGuN3u0O3UqVOh+5555hls3LgRzz33HI4dOwaXy4WpU6eisbFRYcXtu3DhAtLS0vDcc8+1eX9nns/SpUuxb98+vPrqqzh8+DDOnz+PGTNmwO/3d0/RkohIkVWrVsm0tLQ279N1XbpcLrl+/fpQW1NTk3Q6nfL555/voQqjB0Du27cvtN2Z51NXVydtNpt89dVXQ/tUVlZKTdPkW2+91S11sidAREqVlZUhMTERqampePDBB3H27FkAQHl5OaqqqpCdnR3a1+Fw4Pbbb8eRI0dUlRu1zjyf0tJSNDc3h+2TmJiIcePGddtzZggQkTK33HILXn75Zbz99tt44YUXUFVVhczMTNTW1oa+NSwhISHsmISEhFbfKGYEnXk+VVVVsNvtGDRoULv7dDVeSpqIlMnJyQmtjx8/HhkZGRg5ciS2b9+OW2+9FQAghAg7RkrZqs1Ionk+3fmc2RMgol6jX79+GD9+PMrKykKzhC7/C7i6urrVX9NG0Jnn43K54PV68de//rXdfboaQ4CIeg2Px4NPPvkEw4YNQ2pqKlwuF4qKikL3e71eHDp0CJmZmQqrjE5nnk96ejpsNlvYPm63Gx999FG3PWeeDiIiZZ544gnMnDkTycnJqK6uxtNPP42GhgbMmzcPQggsXboU69atw+jRozF69GisW7cOffv2xUMPPaS69DadP38eZ86cCW2Xl5fj5MmTiIuLQ3Jy8hWfj9PpxIIFC/D4448jPj4ecXFxeOKJJzB+/Hjcfffd3VN0t8w5IiLqhDlz5shhw4ZJm80mExMT5X333SdPnz4dul/Xdblq1Srpcrmkw+GQU6ZMkadOnVJYcccOHjwoAbS6zZs3T0rZuefz3XffycWLF8u4uDjZp08fOWPGDFlRUdFtNfNLZYiITIxjAkREJsYQICIyMYYAEZGJMQSIiEyMIUBEZGIMASIiE2MIEFGv5PF4sHr1ang8HtWldIve8vz4OQEi6pUaGhrgdDpRX1+P2NhY1eV0ud7y/NgTICIyMYYAEZGJ8QJyRAbX1NQEr9eruowu19DQELa81nT187Pb7YiJiYn4OI4JEBlYU1MTnH0GwYsm1aWQYi6XC+Xl5REHAXsCRAbm9XrhRROm2GfBanFAWCyAxQJoWmBd04LbArAE14UGWDRITQucENaC6xYBqQlACyzlpdvBdakJSIHgNoJtwaUFkBoA0dKGsH1waZu4uAxrb2M/iPB2aDJsXwgZuFRn8D6E9pEQLW1CQmgSQgTahKYHtjUJTcjAy6Tp0DQJi6bDAglN02HRdFhFcKnpsAgJa3DdKnRYNT8s0GHT/LAKP2xCwqL5YRPBW6hdh034YRW+i/cJH2zQYRM+WISETfhghR/24P1WocOOlmMlrJCwCwGbELBCg01osMICm7CgoVFHSvoX8Hq9DAEiM7LCBquwQwgLIIIhoAVDILS8ZL3dELgYBu2GgNZxCFxs6yAEtE6GgIZQCFy8v40QEJ0MAe1iCGhhISDDQ0AEl5eFQOCNPzwErEIPe7O3aoE3aJvQYAut+2ETAjaB4DKwbm8JDiFhExI2AHaBwLoA7JDBYyRsl4SALRQCGmzCclW/OxwYJiIyMYYAEZGJMQSIiEyMIUBEZGIcGCZSqLi4GBs2bEBpaSncbjf27duHWbNmRfw4PjQDUoOQFkBaAF0LDBJDCwwUQwRmBQlLaERVSg0t02pkcARWSgHI4ACwvGRgWF42O0hvY2BYN9bsIBkcGJZCBo8NtEHTAUhITYfUdEAEl5oe3DewLYUOqfmhI7gUfkgh4Q+u68IPXfPDL/zwCx1+4YdP+OAT/uDNBx8Cs4YCA8N+WOFHc8sAs9Bhh97O7CDAJgArBGzCj4ZGPerfQYYAkUIXLlxAWloaHnnkEdx///0RH2+32+FyuVBc9XrXF0eG4nK5YLfbIz6OHxYj6iWEEFH1BK7VTwxTZKL9xDB7AkQG4vF4wi49rOs6vv32W8THx0MIobAyUklKiZqaGiQmJkLTIhvqZQgQGUhBQQHWrFmjugzqpc6dO4fhw4dHdAxPBxH1Ep05HXR5T6C+vh7Jycn48oMRiO3f9ZP9/q06DR9Njn7QkXqGD804jAOoq6uD0+mM6Fj2BIgMxOFwwOFwtGqP7a8hdsDVXT6gLZ9eb4H1Ki9LQD0g+Kd8NKcE+TkBImrT9woXqC6BegB7AkQKnT9/HmfOnAltl5eX4+TJk4iLi0NycrLCyoDR848r/fepZ7AnQKTQ8ePHMWHCBEyYMAEAkJ+fjwkTJuCpp55SWtfE1T9W+u9Tz2FPgEihrKws9Ma5GfHbjqgugXoIewJEFGZaYprqEqgHMQSIKOSnf5mgugTqYQwBIgr5cAI/E2A2DAEiAsApoWbFECAiAJwSalYMASLilFATYwgQEaeEmhhDgCgK8+fPR3FxseoyugSnhJobQ4AoCo2NjcjOzsbo0aOxbt06VFZWqi4pKhwMJoYAURT27NmDyspKLF68GLt27cKIESOQk5OD3bt3o7m5WXV5ncbBYGIIEEUpPj4ejz32GE6cOIGSkhKMGjUKubm5SExMxLJly1BWVqa6xA5xMJgAhgDRVXO73SgsLERhYSEsFgumT5+O06dPY+zYsdi0aZPq8trFwWACGAJEUWlubsaePXswY8YMpKSkYNeuXVi2bBncbje2b9+OwsJC7NixA2vXrlVdaps4GEwteBVRoigMGzYMuq5j7ty5KCkpwU033dRqn2nTpmHgwIE9XtuVBK4PxMtDUABDgCgKmzZtwuzZsxETE9PuPoMGDUJ5eXkPVtU5vD4QXYqng4iikJub22EA9FacEkqXYwgQmQinhNLlGAJEJsEpodQWhgCRSXBKKLWFIUBkApwSSu1hCBBd4/iVkdQRhgDRNY5TQqkjDAGiXmDLli1ITU1FTEwM0tPT8f7773fJ43JKKF0JQ4BIsZ07d2Lp0qVYuXIlTpw4gdtuuw05OTmoqKi46sfmlFC6EoYAkWIbN27EggULsHDhQowZMwabN29GUlIStm7delWPyymh1Bm8bASRQl6vF6WlpVi+fHlYe3Z2No4caT2l0+PxwOPxhLbr6+sBAA3nW5/3d/7yEHxdXC/1Tj4EvsNCShnxsQwBIoVqamrg9/uRkJAQ1p6QkICqqqpW+xcUFGDNmjWt2lNu/qKNRz/bRVWSUdTW1sLpdEZ0DEOAqBcQQoRtSylbtQHAihUrkJ+fH9quq6tDSkoKKioqIv7Pf61oaGhAUlISzp07h9jYWNXlKFFfX4/k5GTExcVFfCxDgEihwYMHw2KxtPqrv7q6ulXvAAAcDgccDkerdqfTado3wBaxsbGmfw00LfJhXg4MEylkt9uRnp6OoqKisPaioiJkZmYqqorMhD0BIsXy8/ORm5uLiRMnIiMjA9u2bUNFRQUWLVqkujQyAYYAkWJz5sxBbW0t1q5dC7fbjXHjxuHAgQNISUm54rEOhwOrVq1q8xSRWfA1uLrXQMho5hQREdE1gWMCREQmxhAgIjIxhgARkYkxBIiITIwhQERkYgwBIiITYwgQEZkYQ4CIyMQYAkREJsYQICIyMYYAEZGJMQSIiEyMIUBEZGIMASIiE2MIEBGZGEOAiMjEGAJERCbGECAiwyorK8NDDz2EoUOHwuFwYMyYMfjFL36huixD4XcME5Ehffzxx8jMzERycjL+8z//Ey6XC2+//Tb++Z//GTU1NVi1apXqEg2B3zFMRIb093//9zh9+jROnz6N2NjYUPuSJUvwq1/9Cl9//TUGDRqksEJj4OkgIjKcpqYm/O///i/uvfde9O3bFz6fL3SbPn06mpqacPToUdVlGgJ7AkRkOJWVlRg+fHiH+7z88svIzc3toYqMi2MCRGQ4gwYNgsViQW5uLvLy8trcJzU1tYerMib2BIjIkKZOnYpvvvkGJSUlsNvtqssxLIYAERnSxx9/jMmTJ2P06NH48Y9/jBEjRqCxsRFnzpzB7373O7z77ruqSzQEng4iIkMaO3YsPvjgA/z7v/87fvazn6G6uhoDBw7E6NGjMX36dNXlGQZ7AkREJsYpokREJsYQICIyMYYAEZGJMQSIiEyMIUBEZGIMASIiE2MIEBGZGEOAiMjEGAJE1KHi4mLMnDkTiYmJEELg9ddfD7tfSonVq1cjMTERffr0QVZWFk6fPt2tNRUUFGDSpEkYMGAAhg4dilmzZuGzzz5TWtfWrVtx4403IjY2FrGxscjIyMD//M//KKunsxgCRNShCxcuIC0tDc8991yb9z/zzDPYuHEjnnvuORw7dgwulwtTp05FY2Njt9V06NAh5OXl4ejRoygqKoLP50N2djYuXLigrK7hw4dj/fr1OH78OI4fP44777wT//iP/xh6o1fxOnWKJCLqJABy3759oW1d16XL5ZLr168PtTU1NUmn0ymff/75HqururpaApCHDh3qVXUNGjRI/upXv+o19bSFPQEiilp5eTmqqqqQnZ0danM4HLj99ttx5MiRHqujvr4eABAXF9cr6vL7/Xj11Vdx4cIFZGRkKK+nIwwBIopaVVUVACAhISGsPSEhIXRfd5NSIj8/H5MnT8a4ceOU1nXq1Cn0798fDocDixYtwr59+zB27Nhe8Tq1h5eSJqKrJoQI25ZStmrrLosXL8af//xnHD58WHld119/PU6ePIm6ujrs2bMH8+bNw6FDh5TV0xnsCRBR1FwuFwC0+mu2urq61V+93WHJkiXYv38/Dh48GPadw6rqstvtGDVqFCZOnIiCggKkpaXh5z//ufLXqSMMASKKWmpqKlwuF4qKikJtXq8Xhw4dQmZmZrf9u1JKLF68GHv37sW7777b6vuEVdXVVp0ej6fX1NMWng4iog6dP38eZ86cCW2Xl5fj5MmTiIuLQ3JyMpYuXYp169Zh9OjRGD16NNatW4e+ffvioYce6raa8vLy8Nvf/hZvvPEGBgwYEPoL2+l0ok+fPhBC9HhdTz75JHJycpCUlITGxka8+uqreO+99/DWW28pqafTVE5NIqLe7+DBgxJAq9u8efOklIHpmKtWrZIul0s6HA45ZcoUeerUqW6tqa16AMgXX3wxtE9P1/Xoo4/KlJQUabfb5ZAhQ+Rdd90lCwsLldXTWfx6SSIiE+OYABGRiTEEiIhMjCFARGRiDAEiIhNjCBARmRhDgIjIxBgCREQmxhAgoqh4PB6sXr0aHo9HdSkhvbEmoPfWBQD8sBgRRaWhoQFOpxP19fWIjY1VXQ6A3lkT0HvrAtgTICIyNYYAEZGJ8SqiRAbX1NQEr9fb4/9uQ0ND2LI36I01AT1Tl91uR0xMTMTHcUyAyMCamprg7DMIXjSpLoUUc7lcKC8vjzgI2BMgMjCv1wsvmjDFPgtWiwPCYgEsFkDTAuuaFtwWgCW4LjTAokFqWuCEsBZctwhITQBaYCkv3Q6uS01ACgS3EWwLLi2A1ACIljaE7YNL28TFZVh7G/tBhLdDk2H7QsjAdaSD9yG0j4RoaRMSQpMQItAmND2wrUloQgZeJk2HpklYNB0WSGiaDoumwyqCS02HRUhYg+tWocOq+WGBDpvmh1X4YRMSFs0PmwjeQu06bMIPq/BdvE/4YIMOm/DBIiRswgcr/LAH77cKHXa0HCthhYRdCNiEgBUabEKDFRbYhAUNjTpS0r+A1+tlCBCZkRU2WIUdQlgAEQwBLRgCoeUl6+2GwMUwaDcEtI5D4GJbByGgdTIENIRC4OL9bYSA6GQIaBdDQAsLARkeAiK4vCwEAm/84SFgFXrYm71VC7xB24QGW2jdD5sQsAkEl4F1e0twCAmbkLABsAsE1gVghwweI2G7JARsoRDQYBOWq/rd4cAwEZGJMQSIiEyMIUBEZGIMASIiE+PAMJFCxcXF2LBhA0pLS+F2u7Fv3z7MmjUr4sfxoRmQGoS0ANIC6FpgkBhaYKAYIjArSFhCI6pSamiZViODI7BSCkAGB4DlJQPD8rLZQXobA8O6sWYHyeDAsBQyeGygDZoOQEJqOqSmAyK41PTgvoFtKXRIzQ8dwaXwQwoJf3BdF37omh9+4Ydf6PALP3zCB5/wB28++BCYNRQYGPbDCj+aWwaYhQ479HZmBwE2AVghYBN+NDTqUf8OMgSIFLpw4QLS0tLwyCOP4P7774/4eLvdDpfLheKq17u+ODIUl8sFu90e8XH8sBhRLyGEiKonoOoTw9S7RPuJYfYEiAzE4/GEXY5Y13V8++23iI+PhxBCYWWkkpQSNTU1SExMhKZFNtTLECAykIKCAqxZs0Z1GdRLnTt3DsOHD4/oGJ4OIuolOnM66PKeQH19PZKTk/HlByMQ27/rJ/s16E2Y9/1JXf641LV8aMZhHEBdXR2cTmdEx7InQGQgDocDDoejVXtsfw2xA67u8gFtmZ2YCSvPMvV+wT/lozklyM8JEFGbpiWmqS6BegB7AkQKnT9/HmfOnAltl5eX4+TJk4iLi0NycrLCysgsGAJECh0/fhx33HFHaDs/Px8AMG/ePLz00kuKqmIvwEwYAkQKZWVlobfNzdAR/adPyXg4JkBEYXISJ6gugXoQQ4CIQl6/0E91CdTDGAJEFLJ19CjVJVAPYwgQEQAOBpsVQ4CIyMQYAkTEXoCJMQSITI5TQs2NIUAUhfnz56O4uFh1GV2CU0LNjSFAFIXGxkZkZ2dj9OjRWLduHSorK1WXFBVOCSWGAFEU9uzZg8rKSixevBi7du3CiBEjkJOTg927d6O5uVl1eZ3GKaHEECCKUnx8PB577DGcOHECJSUlGDVqFHJzc5GYmIhly5ahrKxMdYkd4mAwAQwBoqvmdrtRWFiIwsJCWCwWTJ8+HadPn8bYsWOxadMm1eW1iYPB1IIhQBSF5uZm7NmzBzNmzEBKSgp27dqFZcuWwe12Y/v27SgsLMSOHTuwdu1a1aW2iYPB1IJXESWKwrBhw6DrOubOnYuSkhLcdNNNrfaZNm0aBg4c2OO1XQkHg+lSDAGiKGzatAmzZ89GTExMu/sMGjQI5eXlPVhV53AwmC7FECCKQm5uruoSosLBYLocxwSIiEyMIUBkEuwFUFsYAkQmwCmh1B6GAJEJcEootYchQHSNa9C/U10C9WIMAaJeYMuWLUhNTUVMTAzS09Px/vvvd9ljzx5+a5c9Fl17GAJEiu3cuRNLly7FypUrceLECdx2223IyclBRUXFVT82B4PpShgCRIpt3LgRCxYswMKFCzFmzBhs3rwZSUlJ2Lp1q+rSyAT4YTEihbxeL0pLS7F8+fKw9uzsbBw5cqTV/h6PBx6PJ7RdX18PAGg433r2z73XjQdgnMtaU/R8wZ+zlDLiYxkCRArV1NTA7/cjISEhrD0hIQFVVVWt9i8oKMCaNWtatafc/EUbj362i6oko6itrYXT6YzoGIYAUS8ghAjbllK2agOAFStWID8/P7RdV1eHlJQUVFRURPyf/1rR0NCApKQknDt3DrGxsarLUaK+vh7JycmIi4uL+FiGAJFCgwcPhsViafVXf3V1daveAQA4HA44HI5W7U6n07RvgC1iY2NN/xpoWuTDvBwYJlLIbrcjPT0dRUVFYe1FRUXIzMxUVBWZCXsCRIrl5+cjNzcXEydOREZGBrZt24aKigosWrRIdWlkAgwBIsXmzJmD2tparF27Fm63G+PGjcOBAweQkpJyxWMdDgdWrVrV5ikis+BrcHWvgZDRzCkiIqJrAscEiIhMjCFARGRiDAEiIhNjCBARmRhDgMjAuvMS1L1dQUEBJk2ahAEDBmDo0KGYNWsWPvvsM9VlKVNQUAAhBJYuXRrRcQwBIoPqzktQG8GhQ4eQl5eHo0ePoqioCD6fD9nZ2bhw4YLq0nrcsWPHsG3bNtx4440RH8spokQGdcstt+Dmm28Ou+T0mDFjMGvWLBQUFCisTI1vvvkGQ4cOxaFDhzBlyhTV5fSY8+fP4+abb8aWLVvw9NNP46abbsLmzZs7fTx7AkQG1HIJ6uzs7LD29i5BbQYtl9WO5iJqRpaXl4d77rkHd999d1TH8xPDRAYU6SWor3VSSuTn52Py5MkYN26c6nJ6zKuvvooPPvgAx44di/oxGAJEBtbZS1Bf6xYvXow///nPOHz4sOpSesy5c+fw2GOPobCwEDExMVE/DkOAyIAivQT1tWzJkiXYv38/iouLMXz4cNXl9JjS0lJUV1cjPT091Ob3+1FcXIznnnsOHo8HFovlio/DMQEiA+IlqAO9nsWLF2Pv3r149913kZqaqrqkHnXXXXfh1KlTOHnyZOg2ceJEPPzwwzh58mSnAgBgT4DIsMx+Ceq8vDz89re/xRtvvIEBAwaEekVOpxN9+vRRXF33GzBgQKvxj379+iE+Pj6icRGGAJFBXc0lqK8FLVNjs7KywtpffPFFzJ8/v+cLMih+ToCIyMQ4JkBEZGIMASIiE2MIEBGZGEOAiMjEGAJERCbGECAiMjGGABGRiTEEiIhMjCFARGRiDAEiIhNjCBARmRhDgIgM55tvvoHL5cK6detCbX/6059gt9tRWFiosDLj4QXkiMiQDhw4gFmzZuHIkSP4/ve/jwkTJuCee+6J6EvWiSFARAaWl5eHd955B5MmTcKHH36IY8eOXdVXLZoRQ4CIDOu7777DuHHjcO7cORw/fhw33nij6pIMJ6IxgaysLCxdurSbSjGOU6dO4fbbb0efPn3wd3/3d1i7di2YpUQ97+zZs/j666+h6zq+/PJL1eUYEr9ZLEINDQ2YOnUq7rjjDhw7dgyff/455s+fj379+uHxxx9XXR6RaXi9Xjz88MOYM2cOvv/972PBggU4deoUEhISVJdmKJ0+HTR//nxs3749rK28vBwjRozojrp6ra1bt2LFihX4y1/+AofDAQBYv349/vu//xtfffUVhBCKKyQyh3/5l3/B7t278eGHH6J///644447MGDAALz55puqSzOUTp8O+vnPf46MjAz88Ic/hNvthtvtRlJSUpv7Llq0CP379+/wVlFR0WVPoif98Y9/xO233x4KAACYNm0avv76a3zxxRfqCiMykffeew+bN2/Gjh07EBsbC03TsGPHDhw+fDj03cPUOZ0+HeR0OmG329G3b1+4XK4O9127di2eeOKJDvdJTEzs7D/dq1RVVbXq/bR0P6uqqpCamqqgKiJzycrKQnNzc1hbcnIy6urq1BRkYN0yJjB06FAMHTq0Ox66V7j8lE/LGTWeCiIio+mWTwxfy6eDXC4Xqqqqwtqqq6sBgANSRGQ4EfUE7HY7/H7/Ffe7lk8HZWRk4Mknn4TX64XdbgcAFBYWIjEx0XSD5ERkfBF9WOxHP/oRTp48iddeew39+/dHXFwcNM1clx+qr6/H9ddfjzvvvBNPPvkkysrKMH/+fDz11FOcIkpEhhNRCHz++eeYN28ePvzwQ3z33XemnCIKBD4slpeXh5KSEgwaNAiLFi3CU089xTEBIjIcXjaCiMjEzHUuh4iIwjAEiIhMjCFARGRiDAEiIhNjCBARmRhDgIjIxBgCREQmxhAgIjIxhgARkYkxBIiITIwhQERkYgwBIiIT+/8BppJeCE2S32IAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "pyro_sim.sim.dovis()" + ] + }, + { + "cell_type": "markdown", + "id": "5e6bba71", + "metadata": {}, + "source": [ + "## Run Simulation\n", + "\n", + "We can edit the parameters in a couple useful ways. First, we can make sure that it does not plot every time step by setting ```vis.dovis``` to ```False```. Second, we can decrease the number of time steps with ```driver.max_steps```. This ensures that the simulation runs in an appropriate amount of time at the cost of some accuracy. In this case, 500 is enough to see that the behavior starts to match what is expected from Woodward and Colella [1]. \n", + "\n", + "After changing these we need to re-initialize the problem, and then we can run it. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3d15152d", + "metadata": {}, + "outputs": [], + "source": [ + "extra_parameters = {'vis.dovis': False, 'driver.max_steps': 500}\n", + "pyro_sim.initialize_problem(problem_name, inputs_file=param_file, inputs_dict=extra_parameters)\n", + "pyro_sim.run_sim()" + ] + }, + { + "cell_type": "markdown", + "id": "44e61a89", + "metadata": {}, + "source": [ + "Now, if we plot the solution, it will visualize the density $\\rho$, the magnitude of the velocity $U = [u, v]^T$, the pressure $p$, and the internal energy $e$." + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "id": "de027154", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAHiCAYAAAAOOXDQAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAQfpJREFUeJzt3X98FPWdP/DX5zP7IwGThYBkyZHEVOkVS4tH4BRUBFqjQbQRTzm1KXjoo/QIZ8jDx7dY2uNHe8SDFrgrxUrPIqVaKL+UazmbWDXAUQVjODH+OLzGBmNiJJoEkOxmZz7fP3azZNkEs5tNZod5PR+PfczMZ2Y2b6Z2X/uZz8ysUEopEBGRLUmzCyAiIvMwBIiIbIwhQERkYwwBIiIbYwgQEdkYQ4CIyMYYAkRENsYQICKyMYYAEZGNMQSIiGyMIUBEZGMMASKytB07dmDChAkYMmQIhgwZgtmzZ6OxsdHssiyDIUBElvWP//iPWLBgAebOnYvnnnsOy5Ytw/PPP4958+aZXZplCD5FlIis6Omnn8a3vvUtVFVV4YYbbgi3f+tb38Kvf/1rfPrpp/B4PCZWaA3sCRCRJf3Lv/wL7rzzzogAAIAvfvGLUErhs88+M6kya2EIEJHlvPPOO3j77bdx++23R6374IMPkJaWhszMTBMqsx6GABFZzuHDhwEAOTk5Ee2GYeB3v/sdioqKICU/3vqCR4mILOfo0aMAgBMnTkS0//jHP8ZHH32ExYsXm1GWJTnMLoCIKFZHjhxBdnY2li1bBpfLhczMTOzbtw9PPPEE1q5di8mTJ5tdomXw6iAishSfz4e0tDR873vfw/Dhw/HjH/8Yzc3N+PKXv4ylS5finnvuMbtES2EIEJGlvPrqq7juuuvwu9/9DrfddpvZ5VgexwSIyFK6xgMmTZpkciWXBoYAEVnK0aNHkZ2dzUtAE4Sng4iIbIw9ASIiG2MIEBHZGEOAiMjGGAJERDbGECAisjGGABGZ7sCBA7j99tuRlZUFIQSeffbZiPVKKaxYsQJZWVlITU3F9OnTUVtba06xA6C8vByTJ09GWloaRo0ahaKiIrz77rsR2wzUMWAIEJHpzp49iwkTJmDjxo09rl+zZg3WrVuHjRs34ujRo/B6vbj55ptx+vTpQa50YFRVVWHRokV45ZVXUFlZiUAggIKCApw9eza8zYAdA0VElEQAqL1794aXDcNQXq9XPfbYY+G2jo4O5fF41M9//nMTKhx4zc3NCoCqqqpSSg3sMWBPgIiSWl1dHZqamlBQUBBuc7vduOmmm8K/K3CpaWtrAwBkZGQAGNhjwBAgoqTW1NQEAFGPicjMzAyvu5QopVBWVoYbbrgB48ePBzCwx4C/J0BEliCEiFhWSkW1XQpKSkrwxhtv4NChQ1HrBuIYsCdAREnN6/UCQNQ33ubm5kvuIXKLFy/Gvn378NJLL2HMmDHh9oE8BgwBIkpqeXl58Hq9qKysDLf5/X5UVVVh6tSpJlaWOEoplJSUYM+ePXjxxReRl5cXsX4gjwFPBxGR6c6cOYP33nsvvFxXV4djx44hIyMDOTk5KC0txerVqzF27FiMHTsWq1evxpAhQ3DfffeZWHXiLFq0CM888wyee+45pKWlhb/xezwepKamQggxcMegX9cWERElwEsvvaQARL3mzZunlApeIrl8+XLl9XqV2+1W06ZNU8ePHze36ATq6d8OQG3ZsiW8zUAdA/6eABGRjXFMgIjIxhgCREQ2xhAgIrIxhgARkY0xBIiIbIwhQERkYwwBIiIbYwgQUdLy+XxYsWIFfD6f2aWYZqCPAW8WI6Kk1d7eDo/Hg7a2NqSnp5tdjikG+hiwJ0BEZGMMASIiG+NTRIksrqOjA36/3+wyBkR7e3vE1I76egxcLhdSUlJifn+OCRBZWEdHBzypw+FHh9mlkMm8Xi/q6upiDgL2BIgszO/3w48OTHMVwaG5ITQN0DRAyuC8lKFlAWiheSEBTUJJGTwhLEPzmoCSApDBqeq+HJpXUkAJhJYRagtNNUBJAKKrDRHboHubOD+NaO9hO4jIdkgVsS2ECj53ObQO4W0URFebUBBSQYhgm5BGcFkqSKGCh0kakFJBkwY0KEhpQJMGHCI0lQY0oeAIzTuEAYfUocGAU+pwCB1OoaBJHU4ReoXbDTiFDocInF8nAnDCgFMEoAkFpwjAAR2u0HqHMOBC174KDii4hIBTCDgg4RQSDmhwCg3tpw3k5r8Pv9/PECCyIweccAgXhNAAEQoBGQqB8LTbfK8hcD4Meg0BefEQON92kRCQfQwBiXAInF/fQwiIPoaAPB8CMiIEVGQIiND0ghAIfvBHhoBDGBEf9g4Z/IB2CglneF6HUwg4BULT4LyrKziEglMoOAG4BILzAnBBhfZRcHYLAWc4BCScQuvXfzscGCYisjGGABGRjTEEiIhsjCFARGRjHBgmMtGBAwewdu1aVFdXo7GxEXv37kVRUVHM7xNAJ6AkhNIApQGGDA4SQwYHiiGCVwUJLTyiqpRE12U1KjQCq5QAVGgAWHUbGFYXXB1k9DAwbFjr6iAVGhhWQoX2DbZBGgAUlDSgpAGI0FQaoW2Dy0oYUFKHgdBU6FBCQQ/NG0KHIXXoQocuDOhCR0AEEBB66BVAAMGrhoIDwzoc0NHZNcAsDLhg9HJ1EOAUgAMCTqGj/bQR93+DDAEiE509exYTJkzAAw88gLvuuivm/V0uF7xeLw40PZv44shSvF4vXC5XzPvxZjGiJCGEiKsncCnfMUx9F+8dw+wJEFmIz+eLeKSwYRj45JNPMGLECAghTKyMzKSUwqlTp5CVlQUpYxvqZQgQWUh5eTlWrlxpdhmUpE6ePIkxY8bEtA9PBxElib6cDrqwJ9DW1oacnBx8+f5/huaK/VRAshEGMPxXRwEV/0CnHQXQiUPYj9bWVng8npj2ZU+AyELcbjfcbndUu+ZKsX4IKAACcCD06Avqu9BX+XhOCfI+ASLqv36eTxAGAAGM2Hw4IeVQ37EnQGSiM2fO4L333gsv19XV4dixY8jIyEBOTo6JlcUg9A2+P/srCbjbeArIDAwBIhO99tprmDFjRni5rKwMADBv3jw89dRTJlU1iLoFyGU7XjG1FLtiCBCZaPr06bDytRmi6y7heHQLgMw/NiKQqKIoJhwTIKK4JSIAACDwf3WJKIfiwBAgorh0DebG7IIAuPyp1xJUEcWDIUBEcVESsV8V1MMgssFHXpiKIUBEMRNxXMgjDER94vCSUPMxBIgoJvEMBgs9uufgOGfdAfFLCUOAiGKiYhwHCI8dXLCfZ9ufElUS9QNDgIj6LK7BYBXdC7j8tbYEVkX9wRAgisP8+fNx4MABs8sYXCrO00A9PAbIeL02MTVRvzEEiOJw+vRpFBQUYOzYsVi9ejUaGhrMLmnAiXhO4ffQaxjxC94ZnEwYAkRx2L17NxoaGlBSUoKdO3fiiiuuQGFhIXbt2oXOzk6zy0u8eHsBF+6jwMdEJxmGAFGcRowYgYcffhg1NTU4cuQIrrrqKhQXFyMrKwtLlizBiRMnzC4xYWLuBfQSGiN+wUtCkw1DgKifGhsbUVFRgYqKCmiahlmzZqG2thZXX3011q9fb3Z5/RbXJaE9DCC729kDSEYMAaI4dHZ2Yvfu3Zg9ezZyc3Oxc+dOLFmyBI2Njdi6dSsqKiqwbds2rFq1yuxS+y2eS0J7Ggy+bDvHApIRnyJKFIfRo0fDMAzce++9OHLkCK655pqobW655RYMGzZs0GtLpLieEtrDqaPh/9uRkHoo8RgCRHFYv3497r77bqSk9P6TjsOHD0ddnbWfjpmoS0Lly68npiBKOIYAURyKi4vNLmHA9eu3Arrh84GSG8cEiChagm4M03x8PlCyYwgQUZR4Lgnt6flAw7by+UDJjiFARBHivST0wucDZbx9LqF10cDgmACRjXWMBFKvP4VhKedgKAGlBAKhBFCha0O72g0lcM7vhPHfw+Du9vy33i4JFQdrBuOfQP3EECCymU+m+jH5qvcRUBIBQ8JQAgaCH/KGEtCUASMUAOEgCJ3ncTkCUF/vgG5I6IZAy6eXQTSlwHCe7wIIXWDY2/H87iSZgaeDiJLApk2bkJeXh5SUFOTn5+PgwYMJ/xtffqAWEx/8H9z41yfgkDocwoBL0+GQBhwi9JLnXy5NhzP0cmsBOGVw6tICcDkCcDl0pAzxQ79Mh/QLyE4BZ5vEsLdFjz0DSk7sCRCZbMeOHSgtLcWmTZtw/fXX44knnkBhYSHeeust5OTk9Pv9L7/3L/ib4R/gnO6Ez3BAlwII/dRjwNDgkEa4FxDuFYjgN3l5kRHic51OeIZ0IJChwWi9DEMbAOdZBSVFzD89TOZhT4DIZOvWrcOCBQvw4IMPYty4cdiwYQOys7Px+OOP9/u9v73kWXzt8nfhFDqkUJBCQbvIB3vXNg5hQEKd7yGEXi6pI0ULwCV1DHX6MTL1LAAg7X3A8Rng6AAcHSqu3yAmc7AnQGQiv9+P6upqLF26NKK9oKAAhw9H32Tl8/ng8/nCy21twRFa3R/9WIbyJf+BDz4bDj2gwRfQ4NMlPgsI+AwBHQp+3QED53sBAMLTLnpo2tUj0AEEDAmXpmOE6ywcQkeT5kCH7oAMANAVoAOBgEBAXYKP1E5SAQSPtVKx98EYAkQmOnXqFHRdR2ZmZkR7ZmYmmpqaorYvLy/HypUro9prn45+UN1tWxJXJ1lDS0sLPB5PTPswBIiSgBCR38CVUlFtAPDoo4+irKwsvNza2orc3FzU19fH/H/+S0V7ezuys7Nx8uRJpKenm12OKdra2pCTk4OMjIyY92UIEJlo5MiR0DQt6lt/c3NzVO8AANxuN9xud1S7x+Ox7Qdgl/T0dNsfAyljH+blwDCRiVwuF/Lz81FZWRnRXllZialTp5pUFdkJewJEJisrK0NxcTEmTZqEKVOmYPPmzaivr8fChQvNLo1sgCFAZLK5c+eipaUFq1atQmNjI8aPH4/9+/cjNzf3c/d1u91Yvnx5j6eI7ILHoH/HQKh4rikiIqJLAscEiIhsjCFARGRjDAEiIhtjCBAR2RhDgIjIxhgCREQ2xhAgIrIxhgARkY0xBIiIbIwhQERkYwwBIiIbYwgQEdkYQ4CIyMYYAkRENsYQICKyMYYAEVnKihUrIITAqVOnelw/fvx4TJ8+fXCLsjCGABGRjTEEiIhsjCFARGRjDAEiIhtjCBAR2RhDgIjIxhgCRGQpDocDAKDreo/rA4EAnE7nYJZkaQwBIrKUzMxMAEBDQ0PUOqUUGhsbw9vQ52MIEJGlzJw5E0II7NixI2rd888/j/b2dnz96183oTJrcphdABFRLK688kqUlJRg7dq1aG1txaxZs5CamoqjR4/isccew6RJk3DfffeZXaZlCKWUMrsIIqJYKKXwxBNP4Mknn8Rbb72FQCCA3NxczJkzB9///vdx2WWXmV2iZTAEiIhsjGMCREQ2xhAgIrIxhgARkY0xBIiIbIwhQERkYwwBIiIbYwgQEdkYQ4CIyMYYAkQWtmnTJuTl5SElJQX5+fk4ePCg2SUllQMHDuD2229HVlYWhBB49tlnI9YrpbBixQpkZWUhNTUV06dPR21trTnFmoQhQGRRO3bsQGlpKZYtW4aamhrceOONKCwsRH19vdmlJY2zZ89iwoQJ2LhxY4/r16xZg3Xr1mHjxo04evQovF4vbr75Zpw+fXqQKzUPHxtBZFHXXnstJk6ciMcffzzcNm7cOBQVFaG8vNzEypKTEAJ79+5FUVERgGAvICsrC6Wlpfjud78LAPD5fMjMzMS//uu/4tvf/raJ1Q4e9gSILMjv96O6uhoFBQUR7QUFBTh8+LBJVVlLXV0dmpqaIo6h2+3GTTfdZKtjyBAgsqBTp05B1/WoH0/JzMxEU1OTSVVZS9dxsvsxZAgQWZgQImJZKRXVRhdn92PIECCyoJEjR0LTtKhvrM3NzfxpxT7yer0AYPtjyBAgsiCXy4X8/HxUVlZGtFdWVmLq1KkmVWUteXl58Hq9EcfQ7/ejqqrKVseQPy9JZFFlZWUoLi7GpEmTMGXKFGzevBn19fVYuHCh2aUljTNnzuC9994LL9fV1eHYsWPIyMhATk4OSktLsXr1aowdOxZjx47F6tWrMWTIEHv9PKUiIsv62c9+pnJzc5XL5VITJ05UVVVVZpeUVF566SUFIOo1b948pZRShmGo5cuXK6/Xq9xut5o2bZo6fvy4uUUPMt4nQERkYxwTICKyMYYAEZGNMQSIiGyMIUBEZGMMASIiG2MIEBHZGEOAyMJ8Ph9WrFgBn89ndimWw2MXxPsEiCysvb0dHo8HbW1tSE9PN7scS+GxC2JPgIjIxhgCREQ2xgfIEV0C2tvbzS7BcrqO2aVy7FwuF1JSUmLej2MCRBbW1taGEcNGQkfA7FLIZF6vF3V1dTEHAXsCRBYmhICOAKa5iuDQ3BCaBmgaIGVwXsrQsgC00LyQgCahpAyeEJaheU1ASQHI4FR1Xw7NKymgBELLCLWFphqgJADR1YaIbdC9TZyfRrT3sB1EZDukitgWQgUfDRpah/A2CqKrTSgIqSBEsE1II7gsFaRQwcMkDUipoEkDGhSkNKBJAw4RmkoDmlBwhOYdwoBD6tBgwCl1OIQOp1DQpA6nCL3C7QacQodDBM6vEwE4YUAKA67QsgN6aF6HQxhwoWtfBQcUXELAKQQckHAKCQc0OIWG9tMGcvPfh9/vZwgQ2ZEDTjiEC0JogAiFgAyFQHjabb7XEDgfBr2GgLx4CJxvu0gIyD6GgEQ4BM6v7yEERB9DQJ4PARkRAioyBERoekEIBD/4I0PAIYyID3uHDH5AO4WEMzyvwykEnAKhaXDeJXRIIeASgFMoOIHz8wJwQYX2UXB2CwFnOAQknELr1387HBgmIrIxhgARkY0xBIiIbIwhQERkYxwYJroEBNAJKAmhNEBpgCGDg8SQwYFiiOBVQUILj6gqJdF1WY0KjcAqJQAVGgBW3QaG1QVXBxk9DAwb1ro6SIUGhpVQoX2DbZAGAAUlDShpACI0lUZo2+CyEgaU1GEgNBU6lFDQQ/OG0GFIHbrQoQsDutAREAEEhB56BRAIXR0U6LoiCDo6uwaYhQEXjF6uDgKcAnBAwCl0tJ824v5vhyFAZKIDBw5g7dq1qK6uRmNjI/bu3YuioqI+7+9yueD1enGg6dkBq5Gswev1wuVyxbwfQ4DIRGfPnsWECRPwwAMP4K677op5/5SUFNTV1cHv9w9AdWQl8d4xzBAgMlFhYSEKCwv79R4pKSlx/Z+fCGAIEFmKz+eLeP69YRj45JNPMGLECAghTKyMzKSUwunTp5GVlQUpY7vehyFAZCHl5eVYuXKl2WVQkjp58iTGjBkT0z58gBxRkhBCfO7A8IU9gba2NuTk5CB7xQ8gk+GUkFChS3ril7f0SIKKsY8AOnEI+9Ha2gqPxxPTvuwJEFmI2+2G2+2OapcpKeaHQAICILVRwCGcCSrIRkJf5eM5JcibxYio/xIQAACQteZwAoqhWLAnQGSiM2fO4L333gsv19XV4dixY8jIyEBOTo6JlcUgQQHw1+tP8lcRTMAQIDLRa6+9hhkzZoSXy8rKAADz5s3DU089ZVJVMUhQAABAoP6DhLwPxYYhQGSi6dOnw8rXZigRzIH+unLJn/r/JhQXjgkQUXyEgjAS0AuwbgZeEhgCRBSXBJ0FwpVl7AWYiSFARDFTMjG9AFcr73I2G0OAiGIjVOzjAL3skL2Sl4SajSFARDFRArGfC+ph+7xnfT1sSIONIUBEfZeowWAAsur1hLwP9Q9DgCgO8+fPx4EDB8wuY9AlbDCYl4QmDYYAURxOnz6NgoICjB07FqtXr0ZDQ4PZJQ24RA0GCz0BxVDCMASI4rB79240NDSgpKQEO3fuxBVXXIHCwkLs2rULnZ2dZpeX1L7wCHsByYQhQBSnESNG4OGHH0ZNTQ2OHDmCq666CsXFxcjKysKSJUtw4sQJs0tMmLh6AT1cEZTaxEtCkw1DgKifGhsbUVFRgYqKCmiahlmzZqG2thZXX3011q9fb3Z5/RfPJaFAjwMIWf/KS0KTDUOAKA6dnZ3YvXs3Zs+ejdzcXOzcuRNLlixBY2Mjtm7dioqKCmzbtg2rVq0yu9R+i+uS0B5kvsrnQyQjPkCOKA6jR4+GYRi49957ceTIEVxzzTVR29xyyy0YNmzYoNeWUPFcEtrLk0Uv2/5KgoqiRGIIEMVh/fr1uPvuu5FykV/zGj58OOrq6gaxqsSL6ymhPQQALwlNXgwBojgUFxebXcLAS9CNYdpnHAxOZhwTIKIeJerGsCuWcTA4mTEEiChKom4MG33QSEA1NJAYAkQUKd5LQnswZM+riXkjGjAcEyCiMNkhoPnP9wCUBAwHYLhUzF8ZORhsDQwBIptK+UhAdgY/6CGD065X1zIUIHVA+ERonQIEAAEo7YI3VIAIADIgoDTeE2AVDAEiG0n5SMB59oIP/FiJ8wEhdEA7Fzl24M8MYFT2pwmplwYexwSIksCmTZuQl5eHlJQU5Ofn4+DBgwl9/8v+IjDsbQF3GyCM0LX/Kjjtmu9aBtDzj78LdT4ApAKkgnIqBIaGNhaAPvYzDLn8LHyBC7sJlKwYAkQm27FjB0pLS7Fs2TLU1NTgxhtvRGFhIerr6xPy/iNrAFebCn74Gzj/gY/z890/9C8eBMGXEoDhUDCcCsoRHC8I5J2DMgR8PicCOkPAKhgCRCZbt24dFixYgAcffBDjxo3Dhg0bkJ2djccff7zf7+39bwXZqYLn9Q0FhIKgeyCIHnoEF4q4Z0Co4Dl/hwKG6FCeTmjj2pEyxA+jU4PSBRwafzTAKjgmQGQiv9+P6upqLF26NKK9oKAAhw9H32Tl8/ng853/bd62tjYAgNHREbVtzvN+dA7RYLgljICA4QgO5hoOAaVFjguEX1poGw3hUz9KQ3CgWFMwHAACCoYEVIoOl9uPv85sxh2jjqFTaRjvbkBD5zA82zIRLdABxfsEBkMAwd+wUCr2AXmGAJGJTp06BV3XkZmZGdGemZmJpqamqO3Ly8uxcuXKqPaTK34Y1faXxJV5UX8G8F+R1QA4Pkh/nbpraWmBx+OJaR+GAFESECLyChulVFQbADz66KMoKysLL7e2tiI3Nxf19fUx/5//UtHe3o7s7GycPHkS6enpZpdjira2NuTk5CAjIyPmfRkCRCYaOXIkNE2L+tbf3Nwc1TsAALfbDbfbHdXu8Xhs+wHYJT093fbHQMrYh3k5MExkIpfLhfz8fFRWVka0V1ZWYurUqSZVRXbCngCRycrKylBcXIxJkyZhypQp2Lx5M+rr67Fw4UKzSyMbYAgQmWzu3LloaWnBqlWr0NjYiPHjx2P//v3Izc393H3dbjeWL1/e4ykiu+Ax6N8xECqea4qIiOiSwDEBIiIbYwgQEdkYQ4CIyMYYAkRENsYQICKyMYYAEZGNMQSIiGyMIUBEZGMMASIiG2MIEBHZGEOAiMjGGAJERDbGECAisjGGABGRjTEEiIhsjCFARGRjDAEiIhtjCBCR5axYsQJCCNTU1GDOnDlIT0+Hx+PBN7/5TXz88cdml2cpDAEisqw777wTV111FXbt2oUVK1bg2WefxS233ILOzk6zS7MM/tA8EVnWnDlzsGbNGgBAQUEBMjMzcf/99+O3v/0t7r//fpOrswb2BIjIsi78oL/nnnvgcDjw0ksvmVSR9TAEiMiyvF5vxLLD4cCIESPQ0tJiUkXWwxAgIstqamqKWA4EAmhpacGIESNMqsh6GAJEZFlPP/10xPJvf/tbBAIBTJ8+3ZyCLIgDw0RkWXv27IHD4cDNN9+M2tpa/OAHP8CECRNwzz33mF2aZbAnQESWtWfPHrzzzjuYM2cO/vmf/xm33347Kioq4HK5zC7NMtgTICLLysnJwb59+8wuw9LYEyAisjGGABGRjQmllDK7CCIiMgd7AkRENsYQICKyMYYAEZGNMQSIiGyMIUBEZGMMASIiG2MIEFHClJeXY/LkyUhLS8OoUaNQVFSEd999N2IbpRRWrFiBrKwspKamYvr06aitrY3YxufzYfHixRg5ciSGDh2KO+64Ax988MGg/juEECgtLbVU3fFgCBBRwlRVVWHRokV45ZVXUFlZiUAggIKCApw9eza8zZo1a7Bu3Tps3LgRR48ehdfrxc0334zTp0+HtyktLcXevXuxfft2HDp0CGfOnMHs2bOh6/qA/xuOHj2KzZs346tf/WpEe7LXHTdFRDRAmpubFQBVVVWllFLKMAzl9XrVY489Ft6mo6NDeTwe9fOf/1wppVRra6tyOp1q+/bt4W0aGhqUlFI9//zzA1rv6dOn1dixY1VlZaW66aab1MMPP2yJuvuDPQEiGjBtbW0AgIyMDABAXV0dmpqaUFBQEN7G7XbjpptuwuHDhwEA1dXV6OzsjNgmKysL48ePD28zUBYtWoTbbrsNX//61yPak73u/uBTRIloQCilUFZWhhtuuAHjx48HcP6XwDIzMyO2zczMxF/+8pfwNi6XC8OHD4/a5sJfEkuk7du34/XXX8fRo0ej1iVz3f3FECCiAVFSUoI33ngDhw4dilonhIhYVkpFtV2oL9vE6+TJk3j44YdRUVGBlJSUXrdLtroTgaeDiCjhFi9ejH379uGll17CmDFjwu1dPwx/4Tfj5ubm8Ldsr9cLv9+PTz/9tNdtEq26uhrNzc3Iz8+Hw+GAw+FAVVUV/v3f/x0OhyP8d5Ot7kRgCBBRwiilUFJSgj179uDFF19EXl5exPq8vDx4vV5UVlaG2/x+P6qqqjB16lQAQH5+PpxOZ8Q2jY2NePPNN8PbJNrXvvY1HD9+HMeOHQu/Jk2ahPvvvx/Hjh3DF77whaSsOyHMHJUmokvLd77zHeXxeNTLL7+sGhsbw6/PPvssvM1jjz2mPB6P2rNnjzp+/Li699571ejRo1V7e3t4m4ULF6oxY8aoF154Qb3++utq5syZasKECSoQCAzav6X71UFWqjtWDAEiShgAPb62bNkS3sYwDLV8+XLl9XqV2+1W06ZNU8ePH494n3PnzqmSkhKVkZGhUlNT1ezZs1V9ff2g/lsuDAGr1B0r/qgMEZGNcUyAiMjGGAJERDbGECAisjGGABGRjTEEiIhsjCFARGRjDAEiMoXP58OKFSvg8/nMLiVmVq79QrxPgIhM0d7eDo/Hg7a2NqSnp5tdTkysXPuF2BMgIrIxhgARkY3x9wSILK6jowN+v9/sMmLW3t4eMbWSZKzd5XJd9LcQesMxASIL6+jogCd1OPzoMLsUMpnX60VdXV3MQcCeAJGF+f1++NGBaa4iODQ3hKYBmgZIGZyXMrQsAC00LySgSSgpgyeEZWheE1BSADI4Vd2XQ/NKCiiB0DJCbaGpBigJQHS1IWIbdG8T56cR7T1sBxHZDqkitoVQwUeVhtYhvI2C6GoTCkIqCBFsE9IILksFKVTwMEkDUipo0oAGBSkNaNKAQ4Sm0oAmFByheYcw4JA6NBhwSh0OocMpFDSpwylCr3C7AafQ4RCB8+tEAE4YkMKAK7TsgB6a1+EQBlzo2lfBAQWXEHAKAQcknELCAQ1OoaH9tIHc/Pfh9/sZAkR25IATDuGCEBogQiEgQyEQnnab7zUEzodBryEgLx4C59suEgKyjyEgEQ6B8+t7CAHRxxCQ50NARoSAigwBEZpeEALBD/7IEHAII+LD3iGDH9BOIeEMz+twCgGnQGganHcJHVIIuATgFApO4Py8AFxQoX0UnN1CwBkOAQmn0Pr13w4HhomIbIwhQERkYwwBIiIbYwgQEdkYB4aJTHTgwAGsXbsW1dXVaGxsxN69e1FUVBTz+wTQCSgJoTRAaYAhg4PEkMGBYojgVUFCC4+oKiXRdVmNCo3AKiUAFRoAVt0GhtUFVwcZPQwMG9a6OkiFBoaVUKF9g22QBgAFJQ0oaQAiNJVGaNvgshIGlNRhIDQVOpRQ0EPzhtBhSB260KELA7rQERABBIQeegUQCF0dFOi6Igg6OrsGmIUBF4xerg4CnAJwQMApdLSfNuL+b5AhQGSis2fPYsKECXjggQdw1113xby/y+WC1+vFgaZnE18cWYrX64XL5Yp5P94sRpQkhBBx9QSsescwJVa8dwyzJ0BkIT6fL+LxxYZh4JNPPsGIESMghDCxMjKTUgqnTp1CVlYWpIxtqJchQGQh5eXlWLlypdllUJI6efIkxowZE9M+PB1ElCT6cjrowp5AW1sbcnJycN20pXA4Yj8VkIxS3m5EoLHJ7DIsJYBOHMJ+tLa2wuPxxLQvewJEFuJ2u+F2u6PaHY6USyYE0NQCh3CaXYW1hL7Kx3NKkPcJEFHScL70P2aXYDvsCRCZ6MyZM3jvvffCy3V1dTh27BgyMjKQk5NjYmXmUIFOs0uwHYYAkYlee+01zJgxI7xcVlYGAJg3bx6eeuopk6qKkQKQgAuTHJWv9f9NKGYMASITTZ8+HZa+NiNBASD98d/xSv3DMQEiMp2set3sEmyLIUBE8UlQLyDlzy39fxOKG0OAiEwV+L86s0uwNYYAEcVMGIqDwZcIhgARxUYFHw9NlwaGABGZgr2A5MAQIKK+S9BgsNah9/9NKCEYAkRxmD9/Pg4cOGB2GZYlDtaYXQKFMASI4nD69GkUFBRg7NixWL16NRoaGswuacAlajA4pfbSP1ZWwhAgisPu3bvR0NCAkpIS7Ny5E1dccQUKCwuxa9cudHZems+/SdRgcODDxoS8DyUGQ4AoTiNGjMDDDz+MmpoaHDlyBFdddRWKi4uRlZWFJUuW4MSJE2aXmDDCSMyjLfiU0OTDECDqp8bGRlRUVKCiogKapmHWrFmora3F1VdfjfXr15tdXv8l8JJQPiU0+TAEiOLQ2dmJ3bt3Y/bs2cjNzcXOnTuxZMkSNDY2YuvWraioqMC2bduwatUqs0tNGrwkNDnxKaJEcRg9ejQMw8C9996LI0eO4Jprrona5pZbbsGwYcMGvbaE4lNCL3kMAaI4rF+/HnfffTdSUnr/Scfhw4ejro7PxQH4lNBkxhAgikNxcbHZJQw4YaiEjAU4X6yBhX8x4ZLHMQEi6lHCBoN13h2czBgCRBQlUZeEcjA4+TEEiChSgi4JlT72AKyAYwJENqe7JZr+1gH/cAOQKnhFkBIQCug6mS8U4DgjMeo1Hc6zfftwlwf4fCArYAgQ2VBdkQYM0YPf+g0FqM7gB78BQAkACqpbCCgF+Ifp+OBrAsIIfmxo5ySy/jsAn0eD/zIBwwkoCUAEp95Kk/5xFBOGAJGNnFzQCWVISEOHocvgB79QgCGglOr6/A/2AoBwCARXAFAq+EGvAD3VwIc3aEirC66TgdBm4vzmlPw4JkCUBDZt2oS8vDykpKQgPz8fBw8eTOj7f7L4DJr/8Rw0hwGpGRASoakR/MSXCkIqCAlAKqjwK/itXmkq+HIowKEAZ3DecAFns4IBIAKhKe8LsxSGAJHJduzYgdLSUixbtgw1NTW48cYbUVhYiPr6+n6/d/0tTqj/9zFcmg6HpkOKrg97I3gncPilwlMhEAwDgeAYQbeX0EIh4FBQbgOGK/iJr7uDf092Kkg/g8BKGAJEJlu3bh0WLFiABx98EOPGjcOGDRuQnZ2Nxx9/vF/vm/PD/8XfTn0HTk2Hy6FDCEBKBdH1QS8Qmlc9BEL3bUKvrjDQFKAZkC4duCwAw62CH/oScLcZcJ01EvKoCRocHBMgMpHf70d1dTWWLl0a0V5QUIDDhw9Hbe/z+eDz+cLLbW1tAIBAoCNiO0MKjJcncOxsNkSHRMCnAF2H6nTACEgYuoRhCChdQKngq+uqIADRYwFdgwRdvQWpIBwGNKngT3PCr1xw+IFOw4AOAZ8QCCg+MXSwBBA81krFfn8HQ4DIRKdOnYKu68jMzIxoz8zMRFNTU9T25eXlWLlyZVT7Kwcei2o7/ELi6ozHO+b+eVtqaWmBx+OJaR+GAFESECLy/IlSKqoNAB599FGUlZWFl1tbW5Gbm4v6+vqY/89/qWhvb0d2djZOnjyJ9PR0s8sxRVtbG3JycpCRkRHzvgwBIhONHDkSmqZFfetvbm6O6h0AgNvthtvtjmr3eDy2/QDskp6ebvtjIGXsw7wcGCYykcvlQn5+PiorI++sqqysxNSpU02qiuyEPQEik5WVlaG4uBiTJk3ClClTsHnzZtTX12PhwoVml0Y2wBAgMtncuXPR0tKCVatWobGxEePHj8f+/fuRm5v7ufu63W4sX768x1NEdsFj0L9jIFQ81xQREdElgWMCREQ2xhAgIrIxhgARkY0xBIiIbIwhQERkYwwBIiIbYwgQEdkYQ4CIyMYYAkRENsYQICKyMYYAEZGNMQSIiGyMIUBEZGMMASIiG2MIEBHZGEOAiMjGGAJERDbGECAiyzpx4gTuu+8+jBo1Cm63G+PGjcPPfvYzs8uyFP7GMBFZ0ltvvYWpU6ciJycHP/nJT+D1evGHP/wB//RP/4RTp05h+fLlZpdoCfyNYSKypFtvvRW1tbWora1Fenp6uH3x4sX4j//4D3z44YcYPny4iRVaA08HEZHldHR04I9//CPuvPNODBkyBIFAIPyaNWsWOjo68Morr5hdpiWwJ0BEltPQ0IAxY8ZcdJtf/epXKC4uHqSKrItjAkRkOcOHD4emaSguLsaiRYt63CYvL2+Qq7Im9gSIyJJuvvlmfPzxxzhy5AhcLpfZ5VgWQ4CILOmtt97CDTfcgLFjx+I73/kOrrjiCpw+fRrvvfce/vM//xMvvvii2SVaAk8HEZElXX311Xj99dfxwx/+EN///vfR3NyMYcOGYezYsZg1a5bZ5VkGewJERDbGS0SJiGyMIUBEZGMMASIiG2MIEBHZGEOAiMjGGAJERDbGECAisjGGABGRjTEEiGjAlJeXY/LkyUhLS8OoUaNQVFSEd999N2Kb+fPnQwgR8bruuutMqrhnK1asiKrR6/WG1yulsGLFCmRlZSE1NRXTp09HbW2tiRX3HUOAiAZMVVUVFi1ahFdeeQWVlZUIBAIoKCjA2bNnI7a79dZb0djYGH7t37/fpIp79+UvfzmixuPHj4fXrVmzBuvWrcPGjRtx9OhReL1e3HzzzTh9+rSJFfcNnx1ERAPm+eefj1jesmULRo0aherqakybNi3c7na7I75ZJyOHw9FjjUopbNiwAcuWLcOcOXMAAFu3bkVmZiaeeeYZfPvb3x7sUmPCngARDZq2tjYAQEZGRkT7yy+/jFGjRuGLX/wiHnroITQ3N5tR3kWdOHECWVlZyMvLw9///d/jz3/+MwCgrq4OTU1NKCgoCG/rdrtx00034fDhw2aV22cMASIaFEoplJWV4YYbbsD48ePD7YWFhXj66afx4osv4ic/+QmOHj2KmTNnwufzmVhtpGuvvRa/+tWv8Ic//AG/+MUv0NTUhKlTp6KlpQVNTU0AgMzMzIh9MjMzw+uSGU8HEdGgKCkpwRtvvIFDhw5FtM+dOzc8P378eEyaNAm5ubn4/e9/Hz69YrbCwsLw/Fe+8hVMmTIFV155JbZu3RoexBZCROyjlIpqS0bsCRDRgFu8eDH27duHl1566XN/G3j06NHIzc3FiRMnBqm62A0dOhRf+cpXcOLEifA4wYXf+pubm6N6B8mIIUBEA0YphZKSEuzZswcvvvhin373t6WlBSdPnsTo0aMHocL4+Hw+vP322xg9ejTy8vLg9XpRWVkZXu/3+1FVVYWpU6eaWGXfMASIaMAsWrQIv/71r/HMM88gLS0NTU1NaGpqwrlz5wAAZ86cwSOPPII//elPeP/99/Hyyy/j9ttvx8iRI3HnnXeaXP15jzzyCKqqqlBXV4dXX30Vf/d3f4f29nbMmzcPQgiUlpZi9erV2Lt3L958803Mnz8fQ4YMwX333Wd26Z9PERENEAA9vrZs2aKUUuqzzz5TBQUF6vLLL1dOp1Pl5OSoefPmqfr6enMLv8DcuXPV6NGjldPpVFlZWWrOnDmqtrY2vN4wDLV8+XLl9XqV2+1W06ZNU8ePHzex4r7jz0sSEdkYTwcREdkYQ4CIyMYYAkRENsYQICKyMYYAEZGNMQSIiGyMIUBEZGMMASIyhc/nw4oVK5LqaaHxsPq/gzeLEZEp2tvb4fF40NbWhvT0dLPLiZvV/x3sCRAR2RhDgIjIxvijMkQW19HRAb/fb3YZMWtvb4+YWlWy/DtcLhdSUlJi3o9jAkQW1tHRAU/qcPjRYXYpZDKv14u6urqYg4A9ASIL8/v98KMD01xFcGhuCE0DNA2QMjgvZWhZAFpoXkhAk1BSBk8Iy9C8JqCkAGRwqrovh+aVFFACoWWE2kJTDVASgOhqQ8Q26N4mzk8j2nvYDiKyHVJFbAuhgs+nDq1DeBsF0dUmFIRUECLYJqQRXJYKUqjgYZIGpFTQpAENClIa0KQBhwhNpQFNKDhC8w5hwCF1aDDglDocQodTKGhSh1OEXuF2A06hwyEC59eJAJwwIIUBV2jZAT00r8MhDLjQta+CAwouIeAUAg5IOIWEAxqcQkP7aQO5+e/D7/czBIjsyAEnHMIFITRAhEJAhkIgPO0232sInA+DXkNAXjwEzrddJARkH0NAIhwC59f3EAKijyEgz4eAjAgBFRkCIjS9IASCH/yRIeAQRsSHvUMGP6CdQsIZntfhFAJOgdA0OO8SOqQQcAnAKRScwPl5AbigQvsoOLuFgDMcAhJOofXrvx0ODBMR2RhDgIjIxhgCREQ2xhAgIrIxDgwTmejAgQNYu3Ytqqur0djYiL1796KoqCjm9wmgE1ASQmmA0gBDBgeJIYMDxRDBq4KEFh5RVUqi67IaFRqBVUoAKjQArLoNDKsLrg4yehgYNqx1dZAKDQwroUL7BtsgDQAKShpQ0gBEaCqN0LbBZSUMKKnDQGgqdCihoIfmDaHDkDp0oUMXBnShIyACCAg99AogELo6KNB1RRB0dHYNMAsDLhi9XB0EOAXggIBT6Gg/bcT93yBDgMhEZ8+exYQJE/DAAw/grrvuinl/l8sFr9eLA03PJr44shSv1wuXyxXzfrxZjChJCCHi6glY9Y5hSqx47xhmT4DIQnw+X8Qjiw3DwCeffIIRI0ZACGFiZWQmpRROnTqFrKwsSBnbUC9DgMhCysvLsXLlSrPLoCR18uRJjBkzJqZ9eDqIKEn05XTQhT2BtrY25OTk4K8e+x5kHKcCko3SFMYuOQ4V6DS7FEsJoBOHsB+tra3weDwx7cueAJGFuN1uuN3uqHaZkgKZau0QUJqC1AU0HYBwml2OtYS+ysdzSpD3CRCR6ZSmIHSBsQuPmF2K7bAnQGSiM2fO4L333gsv19XV4dixY8jIyEBOTo6JlQ2ergCQHfxOagaGAJGJXnvtNcyYMSO8XFZWBgCYN28ennrqKZOqGjxKAkIPnsK4svQVk6uxJ4YAkYmmT58OS1+boQQg4qtfSUCEbnS9/E/9exwyxY8hQETxiycARCg7uj3pYNjWPyWuJooJT8IRUXxUHDen9RAAHAw2F0OAiOITYy9AhT5tugdAXEFCCcUQIKLYxdULUOHr2buM/c6riamH4sYQIKLYxDMY7FThq4C6ONo4GJwMGAJEFJtYTwM5FJQe3XPI+y4Hg5MBQ4CI+i6O00DigoFgAEh7hxcmJguGAFEc5s+fjwMHDphdxuCL4zQQAtHB4d1wOEEFUX8xBIjicPr0aRQUFGDs2LFYvXo1GhoazC5p4MXRC1AKUYPBXyypSUw9lBAMAaI47N69Gw0NDSgpKcHOnTtxxRVXoLCwELt27UJn5yX4GOR4B4N76AXwMdHJhSFAFKcRI0bg4YcfRk1NDY4cOYKrrroKxcXFyMrKwpIlS3DixAmzSzSN0hRUD799zhvDkg9DgKifGhsbUVFRgYqKCmiahlmzZqG2thZXX3011q9fb3Z5/RdHL0B0ezAcQhM+JTQ58X8Vojh0dnZi9+7dmD17NnJzc7Fz504sWbIEjY2N2Lp1KyoqKrBt2zasWrXK7FL7L9ZLQrULBoNDu/MpocmJ12kRxWH06NEwDAP33nsvjhw5gmuuuSZqm1tuuQXDhg0b9NoSKs5eAPTItstf4Y1hyYohQBSH9evX4+6770bKRX7Xd/jw4airqxvEqgZAPJeEdkYPBg97ijeGJSuGAFEciouLzS5h4MXaCxDBS0IvjAAOBic3jgkQUbR4Lgl1RF8SKj/jaaBkxxAgon5TEj1eEnplGU8DJTuGABFFimcwWIt+SuiQep5ttgL+r0RkYzKtE18e0wiXFrycxwg9GqK3qd/QcOJkJnDu/Gke1cNpIAD4q9V8PpAVMASIbGbMFz7GlektCCjZ7UM+NN91Y5dQMJSImAKAQxr4yhUNoX0EWs4NRUfAAX9nMBSEAJwOHc9O+A88hBsG/x9HMePpIKIksGnTJuTl5SElJQX5+fk4ePBgwv9GwcQ3cWv+G/jSsGZIYcARekmhIENThzAgoYIvEXw5pBGe774MAJcPOQPDENANCcOQGDbkHH7xlW04GRia8PppYLAnQGSyHTt2oLS0FJs2bcL111+PJ554AoWFhXjrrbeQk5PT7/e/6Zp3cLn7DHQlEDCC39glJDpDc1AI9QAMGEqGP+C77vQ1IOC44AcBDCHgkjraO91IS/HB3+nAHVcex/VpJ/CxfhmGCn+/66bBwZ4AkcnWrVuHBQsW4MEHH8S4ceOwYcMGZGdn4/HHH+/3ez90fRXyhpxCqvRDEwoOqcMpDEhhhKdSdP/mH2wDEA4D2e1Z0F1tLhlAitYJl9ShSQM/nPAcrk87gQ7DiU7lgB51twAlK/YEiEzk9/tRXV2NpUuXRrQXFBTg8OHogVWfzwefzxdebmtrAwAYHR1R2373xt/js3Mu6IYbfsMBGAEYSoOuNChDgwCgDA0KAoYSwXEBiPAyEBwr6E5IHamyE7oSOGc44VRnkdLZiaY2F/7K8Qme/mgSrhz6Me4ZdhQBxUdGD5ZAqF+nVIz3doAhQGSqU6dOQdd1ZGZmRrRnZmaiqakpavvy8nKsXLkyqr1h6eqotpLElfm5KsNz/wsA2AAAeG4QKyAAaGlpgcfjiWkfhgBREhAi8vSJUiqqDQAeffRRlJWVhZdbW1uRm5uL+vr6mP/Pf6lob29HdnY2Tp48ifT0dLPLMUVbWxtycnKQkZER874MASITjRw5EpqmRX3rb25ujuodAIDb7Ybb7Y5q93g8tv0A7JKenm77YyBl7MO8HBgmMpHL5UJ+fj4qKysj2isrKzF16lSTqiI7YU+AyGRlZWUoLi7GpEmTMGXKFGzevBn19fVYuHCh2aWRDTAEiEw2d+5ctLS0YNWqVWhsbMT48eOxf/9+5Obmfu6+brcby5cv7/EUkV3wGPTvGAgVzzVFRER0SeCYABGRjTEEiIhsjCFARGRjDAEiIhtjCBBZ2GA8gjpZlZeXY/LkyUhLS8OoUaNQVFSEd9991+yyTFNeXg4hBEpLS2PajyFAZFFdj6BetmwZampqcOONN6KwsBD19fVmlzYoqqqqsGjRIrzyyiuorKxEIBBAQUEBzp49a3Zpg+7o0aPYvHkzvvrVr8a8Ly8RJbKoa6+9FhMnTox45PS4ceNQVFSE8vJyEyszx8cff4xRo0ahqqoK06ZNM7ucQXPmzBlMnDgRmzZtwo9+9CNcc8012LBhQ5/3Z0+AyIK6HkFdUFAQ0d7bI6jtoOux2vE8RM3KFi1ahNtuuw1f//rX49qfdwwTWVCsj6C+1CmlUFZWhhtuuAHjx483u5xBs337drz++us4evRo3O/BECCysL4+gvpSV1JSgjfeeAOHDh0yu5RBc/LkSTz88MOoqKhASkpK3O/DECCyoFgfQX0pW7x4Mfbt24cDBw5gzJgxZpczaKqrq9Hc3Iz8/Pxwm67rOHDgADZu3AifzwdN0z73fTgmQGRBfAR1sNdTUlKCPXv24MUXX0ReXp7ZJQ2qr33tazh+/DiOHTsWfk2aNAn3338/jh071qcAANgTILIsuz+CetGiRXjmmWfw3HPPIS0tLdwr8ng8SE1NNbm6gZeWlhY1/jF06FCMGDEipnERhgCRRfXnEdSXgq5LY6dPnx7RvmXLFsyfP3/wC7Io3idARGRjHBMgIrIxhgARkY0xBIiIbIwhQERkYwwBIiIbYwgQEdkYQ4CIyMYYAkRENsYQICKyMYYAEZGNMQSIiGyMIUBElvPxxx/D6/Vi9erV4bZXX30VLpcLFRUVJlZmPXyAHBFZ0v79+1FUVITDhw/jS1/6Ev7mb/4Gt912W0w/sk4MASKysEWLFuGFF17A5MmT8T//8z84evRov35q0Y4YAkRkWefOncP48eNx8uRJvPbaa/jqV79qdkmWE9OYwPTp01FaWjpApcTG5/Nh8eLFGDlyJIYOHYo77rgDH3zwwefut2nTJuTl5SElJQX5+fk4ePBgxHohRI+vtWvXAgDef//9XrfZuXNneJsFCxYgLy8PqampuPLKK7F8+XL4/f6Iv/XHP/4RU6dORVpaGkaPHo3vfve7CAQC4fW9/a3nn3++v4eP6JLw5z//GR9++CEMw8Bf/vIXs8uxJMsODJeWlmLv3r3Yvn07Dh06hDNnzmD27NnQdb3XfXbs2IHS0lIsW7YMNTU1uPHGG1FYWIj6+vrwNo2NjRGvX/7ylxBC4K677gIAZGdnR22zcuVKDB06FIWFhQCAd955B4Zh4IknnkBtbS3Wr1+Pn//85/je974X/jtvvPEGZs2ahVtvvRU1NTXYvn079u3bh6VLl0bV/cILL0T8vZkzZybqMBJZlt/vx/3334+5c+fiRz/6ERYsWICPPvrI7LKsR/XRvHnzFICIV11dXV93T6jW1lbldDrV9u3bw20NDQ1KSqmef/75Xvf727/9W7Vw4cKIti996Utq6dKlve7zjW98Q82cOfOi9VxzzTXqH/7hHy66zZo1a1ReXl54+dFHH1WTJk2K2Gbv3r0qJSVFtbe3K6WUqqurUwBUTU3NRd+byI4eeeQRdcUVV6i2tjal67qaNm2auu2228wuy3L63BP4t3/7N0yZMgUPPfRQ+BtpdnZ2j9suXLgQl1122UVf3b99x6q6uhqdnZ0oKCgIt2VlZWH8+PE4fPhwj/v4/X5UV1dH7AMABQUFve7z0Ucf4fe//z0WLFhw0VqOHTt20W0AoK2tDRkZGeFln88XNYCVmpqKjo4OVFdXR7TfcccdGDVqFK6//nrs2rXron+HyA5efvllbNiwAdu2bUN6ejqklNi2bRsOHToU/u1h6ps+/9C8x+OBy+XCkCFD4PV6L7rtqlWr8Mgjj1x0m6ysrL7+6ShNTU1wuVwYPnx4RHtmZiaampp63OfUqVPQdR2ZmZl93mfr1q1IS0vDnDlzeq3lySefxLhx4zB16tRet/m///s//PSnP8VPfvKTcNstt9yCDRs24De/+Q3uueceNDU14Uc/+hGA4CkpALjsssuwbt06XH/99ZBSYt++fZg7dy62bt2Kb37zm73+PaJL3fTp09HZ2RnRlpOTg9bWVnMKsrA+h0AsRo0ahVGjRg3EW1+UUgpCiItuc+H6i+3zy1/+Evfff3+vl5ydO3cOzzzzDH7wgx/0+vc+/PBD3Hrrrbj77rvx4IMPhtsLCgqwdu1aLFy4EMXFxXC73fjBD36AQ4cOQdM0AMDIkSOxZMmS8D6TJk3Cp59+ijVr1jAEiCghBmRgeKBPB3m9Xvj9fnz66acR7c3NzVHf9LuMHDkSmqZFfevvbZ+DBw/i3XffjfjgvtCuXbvw2Wef4Vvf+laP6z/88EPMmDEDU6ZMwebNm6PWl5WVobW1FfX19Th16hS+8Y1vAADy8vJ6/ZvXXXcdTpw40et6IqJYxNQTcLlcF736pstAnw7Kz8+H0+lEZWUl7rnnHgDBUyhvvvkm1qxZ0+M+LpcL+fn5qKysxJ133hlur6ysDH/4dvfkk08iPz8fEyZM6LWOJ598EnfccQcuv/zyqHUNDQ2YMWMG8vPzsWXLFkjZc94KIcLH4je/+Q2ys7MxceLEXv9mTU0NRo8e3et6IqKYxDKK/NBDD6nJkyeruro69fHHHytd1wdmuLoPFi5cqMaMGaNeeOEF9frrr6uZM2eqCRMmqEAgEN5m5syZ6qc//Wl4efv27crpdKonn3xSvfXWW6q0tFQNHTpUvf/++xHv3dbWpoYMGaIef/zxXv/+iRMnlBBC/dd//VfUuoaGBnXVVVepmTNnqg8++EA1NjaGX92tWbNGvfHGG+rNN99Uq1atUk6nU+3duze8/qmnnlJPP/20euutt9Q777yj1q5dq5xOp1q3bl2sh4uIqEcxhcC7776rrrvuOpWammrqJaJKKXXu3DlVUlKiMjIyVGpqqpo9e7aqr6+P2CY3N1ctX748ou1nP/uZys3NVS6XS02cOFFVVVVFvfcTTzyhUlNTVWtra69//9FHH1VjxozpMQi3bNkSdTlt16u7GTNmKI/Ho1JSUtS1116r9u/fH7H+qaeeUuPGjVNDhgxRaWlpKj8/X23btu3zDg0RUZ/xsRFERDZm2TuGiYio/xgCREQ2xhAgIrIxhgARkY0xBIiIbIwhQERkYwwBIiIbYwgQEdkYQ4CIyMYYAkRENsYQICKyMYYAEZGN/X8JYXNPAsPRFAAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "pyro_sim.sim.dovis()" + ] + }, + { + "cell_type": "markdown", + "id": "530159e5", + "metadata": {}, + "source": [ + "Since the equations are solving for the conserved quantities, i.e. density, momentum ($\\rho U$), and total energy density ($\\rho E$), some of these results require solving for other values. Pyro does this implicitly when calling ```dovis()```, but we can also do it on our own via ```cons_to_prim()```.\n", + "\n", + "Here is an example of plotting just one of the outputs. Density ($\\rho$) is shown since it is easy to visually compare with the results from Woodward and Colella [1]. They additionally show the rest of the primitive variables---density $\\rho$, pressure $p$, x velocity $u$, and y velocity $v$---as well as a measure of specific entropy $p/\\rho^\\gamma$. " + ] + }, + { + "cell_type": "code", + "execution_count": 77, + "id": "e23bf7f1", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA5IAAAE6CAYAAABkoW24AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAeMhJREFUeJzt3Xl8FdXdP/DPzNwtLAmbZGEHEVAEKYgsLiAVxL2VSrVlcasIqECtgj6t6ONjqq1KXQBtQeqGVAFxQQV/CmhBH1BQK6j4CCRCQggI2UjuvTPf3x93yd2X5C7Jzef9ek1N5p6ZOXeY3Oab7znfo4iIgIiIiIiIiChGaro7QERERERERM0LA0kiIiIiIiKKCwNJIiIiIiIiigsDSSIiIiIiIooLA0kiIiIiIiKKCwNJIiIiIiIiigsDSSIiIiIiIooLA0kiIiIiIiKKCwNJIiIiIiIiigsDSSKKy6effopf/OIX6N69O6xWK3JzczFy5Ej8/ve/b9D5pk+fjp49ezbo2K1bt2LhwoU4fvx40GtjxozBmDFjGnTeVDl27Bh+/etfo3PnzlAUBVdddRV2796NhQsXYv/+/Snrx/79+6EoSsjtlVdeCWr/ww8/4Je//CXatWuHNm3a4KKLLsLnn38e8tyvvPIKzjrrLNhsNhQUFGDOnDmoqqqK2qdDhw5h4cKF2LVrV9BrCxcuhKIocb/PVFm/fj0WLlyY7m5EtX//flx66aXo0KEDFEXBnDlz0tIPRVGaxP3q2bMnpk+f7v0+0jNIRESAKd0dIKLm4+2338YVV1yBMWPG4JFHHkF+fj5KSkqwY8cOvPLKK3j00UdT2p+tW7fi/vvvx/Tp09GuXTu/1xYvXpzSvjTEf//3f2Pt2rVYvnw5+vTpgw4dOuDLL7/E/fffjzFjxjQ4wG6o2267Ddddd53fvr59+/p9f+TIEZx33nlo3749li9fDpvNhsLCQowZMwbbt29Hv379vG1feukl/Pa3v8VNN92Exx9/HN999x3uvvtu7N69Gxs2bIjYl0OHDuH+++9Hz549cdZZZ/m9dtNNN+Hiiy9u3JtNovXr1+Ppp59uEsFRJHPnzsWnn36K5cuXIy8vD/n5+enuUlqtXbsW2dnZ3u8jPYNERMRAkoji8Mgjj6BXr1547733YDLVf3z8+te/xiOPPJLGngU7/fTT092FqP7zn/+gT58++M1vfuPd9+WXXyblWjU1NWjVqlXENt27d8eIESMitvnLX/6CI0eOYOvWrejRowcA4Nxzz0WfPn3wpz/9CatWrQIA6LqOP/zhDxg/fjz+/ve/AwDGjh2Ltm3b4je/+Q3eeecdTJw4sUHvpWvXrujatWuDjqV6//nPfzB8+HBcddVVCTmfrutwOp2wWq0JOV+qDRkyJN1dICJqVji0lYhidvToUXTq1MkviPRQVf+PE8Mw8Mgjj6B///6wWq3o3Lkzpk6dih9//DHiNTzDLFesWBH0mu8QuIULF+IPf/gDAKBXr17eoZibNm0CEHpo67FjxzBz5kx06dIFFosFvXv3xr333ou6urqg68yePRsvvPACBgwYgFatWmHw4MF46623IvYdAGpra/H73/8eZ511FnJyctChQweMHDkS69atC3qP77//Pvbs2ePt+4oVK/CrX/0KgCvo8t3v8f7772PcuHHIzs5Gq1atMHr0aPy///f//PrgGfr5+eefY9KkSWjfvj369OkTte+xWLt2LS688EJvEAkA2dnZ+OUvf4k333wTTqcTAPDJJ5+gpKQE119/vd/xv/rVr9CmTRusXbs27DU2bdqEs88+GwBw/fXXe++D77994NDWnj174rLLLsNbb72FIUOGICsrCwMGDPD+m61YsQIDBgxA69atMXz4cOzYsSPoujt27MAVV1yBDh06wGazYciQIfjXv/7l16ampgZ33nknevXqBZvNhg4dOmDYsGFYuXIlANdQ7aeffhoA/IYIe4YqiwgWL16Ms846C1lZWWjfvj0mTZqEH374we86Y8aMwcCBA/HRRx9hxIgRyMrKQpcuXfDHP/4Ruq77tV2yZAkGDx6MNm3aoG3btujfvz/uueeeiPdXURR8//33eOedd4L6WFRUhN/+9rfo3LkzrFYrBgwYgEcffRSGYXjP4XmGH3nkETz44IPo1asXrFYrPvzww7DXraiowM0334yOHTuiTZs2uPjii/Hdd9+FbLt3715cd911fn3w3NfA97Fy5Urce++9KCgoQHZ2Nn7+85/j22+/9Wu7c+dOXHbZZd7zFRQU4NJLL/X7PPId2hrpGXzhhRegKAq2bdsW1O8HHngAZrMZhw4dCnsfiIgyhhARxeimm24SAHLbbbfJJ598Ina7PWzb3/3udwJAZs+eLe+++64sXbpUTjnlFOnWrZscOXLE227atGnSo0cP7/f79u0TAPLcc88FnROA3HfffSIiUlxcLLfddpsAkDVr1si2bdtk27ZtcuLECRERueCCC+SCCy7wHnvy5EkZNGiQtG7dWv7617/Khg0b5I9//KOYTCa55JJLgq7Ts2dPGT58uPzrX/+S9evXy5gxY8RkMsn//d//RbxHx48fl+nTp8sLL7wgH3zwgbz77rty5513iqqq8s9//lNERGpra2Xbtm0yZMgQ6d27t7fv+/fvl4ceekgAyNNPP+3dX1ZWJiIiL7zwgiiKIldddZWsWbNG3nzzTbnssstE0zR5//33vX247777BID06NFD7r77btm4caO8/vrrYfvsuecdO3YUs9ksWVlZMnr0aFm3bp1fu5qaGlEURf7whz8EneOpp54SAPLtt9+KiMjSpUsFgHz99ddBbYcNGyYjR44M258TJ07Ic889JwDkv/7rv7z3obi42O/9+erRo4d07dpVBg4cKCtXrpT169fLOeecI2azWf70pz/J6NGjZc2aNbJ27Vo57bTTJDc3V2pqarzHf/DBB2KxWOS8886TVatWybvvvivTp08PehZvueUWadWqlTz22GPy4YcfyltvvSV//vOf5cknnxQRke+//14mTZokALz93rZtm9TW1oqIyM033yxms1l+//vfy7vvvisvv/yy9O/fX3Jzc6W0tNR7nQsuuEA6duwoBQUF8sQTT8h7770nt99+uwCQWbNmedutXLnS+zO5YcMGef/992Xp0qVy++23R7y/27Ztk7y8PBk9erRfH8vKyqRLly5yyimnyNKlS+Xdd9+V2bNnCwC59dZbvefwPDNdunSRsWPHymuvvSYbNmyQffv2hbymYRgyduxYsVqt8j//8z+yYcMGue+++6R3795+P9ciIl9//bXk5OTImWeeKc8//7xs2LBBfv/734uqqrJw4UJvuw8//ND7s/qb3/xG3n77bVm5cqV0795d+vbtK06nU0REqqqqpGPHjjJs2DD517/+JZs3b5ZVq1bJjBkzZPfu3X7P0LRp06I+g3V1dZKXlye/+c1v/N6jw+GQgoIC+dWvfhX23hMRZRIGkkQUs/Lycjn33HMFgAAQs9kso0aNksLCQqmsrPS227NnjwCQmTNn+h3/6aefCgC55557vPsaGkiKiPzlL38RACF/eQ0MJD2Bzb/+9S+/dg8//LAAkA0bNvhdJzc3VyoqKrz7SktLRVVVKSwsDHd7QnI6neJwOOTGG2+UIUOGBPXxjDPO8Nv36quvCgD58MMP/fZXV1dLhw4d5PLLL/fbr+u6DB48WIYPH+7d5wm0/vSnP8XUx0OHDsnNN98s//rXv+Sjjz6Sl156SUaMGCEA5O9//7u33cGDBwVAyHvw8ssvCwDZunWriIj8z//8jwCQkpKSoLbjx4+X0047LWKftm/fHvY5CBdIZmVlyY8//ujdt2vXLgEg+fn5Ul1d7d3/+uuvCwB54403vPv69+8vQ4YMEYfD4Xfeyy67TPLz80XXdRERGThwoFx11VUR+z5r1qyg/omIbNu2TQDIo48+6re/uLhYsrKy5K677vLuu+CCCwRAUDB/8803i6qqcuDAARERmT17trRr1y5if8Lp0aOHXHrppX775s+fLwDk008/9dt/6623iqIo3j8UeH5O+/TpE/EPSh7vvPOOAJC//e1vfvs9z4nvz/WECROka9eu3j8KecyePVtsNpscO3ZMROoDycA/BP3rX//yBvIiIjt27BAAEf+YIuIfSIpEfwYtFoscPnzYu2/VqlUCQDZv3hzxOkREmYJDW4koZh07dsRHH32E7du3489//jOuvPJKfPfdd1iwYAHOPPNMlJeXA4B3eJtvBUQAGD58OAYMGBA0FDMVPvjgA7Ru3RqTJk3y2+/pY2CfPPP5PHJzc9G5c2ccOHAg6rVeffVVjB49Gm3atIHJZILZbMayZcuwZ8+eBvd/69atOHbsGKZNmwan0+ndDMPAxRdfjO3bt6O6utrvmKuvvjqmc+fn5+PZZ5/Fr371K5x77rm47rrrsGXLFgwZMgTz58/3Dlf1iFQxNfC1cG2TUXX1rLPOQpcuXbzfDxgwAIBrmKjv/FDPfs+/5ffff49vvvnGO1fV9/5ecsklKCkp8Q6VHD58ON555x3Mnz8fmzZtwsmTJ2Pu31tvvQVFUfDb3/7W7xp5eXkYPHiwd1i2R9u2bXHFFVf47bvuuutgGAa2bNni7c/x48dx7bXXYt26dd6fwYb64IMPcPrpp2P48OF++6dPnw4RwQcffOC3/4orroDZbI56Xs9ngu98YABBxZ1qa2vx//7f/8MvfvELtGrVKujfora2Fp988klQH3wNGjQIQP2/76mnnor27dvj7rvvxtKlS7F79+6o/Y3m1ltvBQDv/F8AeOqpp3DmmWfi/PPPb/T5iYiaAwaSRBS3YcOG4e6778arr76KQ4cOYe7cudi/f7+34M7Ro0cBIGQVyIKCAu/rqXT06FHk5eUFBTCdO3eGyWQK6lPHjh2DzmG1WqMGDmvWrME111yDLl264MUXX8S2bduwfft23HDDDaitrW1w/w8fPgwAmDRpEsxms9/28MMPQ0Rw7Ngxv2MaU4XTbDZj8uTJOHr0KPbu3QsAaN++PRRFCfnv57l2hw4dANTfv3BtPe0SKfCcFosl4n7Pv4fn3t55551B93bmzJkA4A3QnnjiCdx99914/fXXMXbsWHTo0AFXXXWV9x5FcvjwYYgIcnNzg67zySefBAWBubm5QefIy8sDUH9fp0yZguXLl+PAgQO4+uqr0blzZ5xzzjnYuHFj1P6EcvTo0bA/t77X9Yj1GTt69ChMJlPQz5Xn/fi2czqdePLJJ4Pu0SWXXAIAQfcp8JyeYj+en9WcnBxs3rwZZ511Fu655x6cccYZKCgowH333QeHwxFT/wPl5uZi8uTJeOaZZ6DrOr788kt89NFHmD17doPOR0TUHLFqKxE1itlsxn333YfHH38c//nPfwDU/2JXUlISVF3z0KFD6NSpU9jz2Ww2AAgqgNPY4LNjx4749NNPISJ+wWRZWRmcTmfEPsXjxRdfRK9evbBq1Sq/6wS+n3h5+vfkk0+GrawaGHg0NusnIgDqCyllZWXh1FNPxVdffRXU9quvvkJWVhZ69+4NADjzzDO9+30r6DqdTnzzzTe49tprG9W3RPLc2wULFuCXv/xlyDaeZU1at26N+++/H/fffz8OHz7szU5efvnl+Oabb6JeR1EUfPTRRyErmwbu8wS4vkpLSwH4B0/XX389rr/+elRXV2PLli247777cNlll+G7777zK4oUi44dO6KkpCRov6d4TODPSazPWMeOHeF0OnH06FG/vnvej0f79u2haRqmTJmCWbNmhTxXr169YrqmrzPPPBOvvPIKRARffvklVqxYgQceeABZWVmYP39+3OcDgDvuuAMvvPAC1q1bh3fffRft2rULyrgSEWUyZiSJKGahfsEE4B2y6claXHjhhQBcQZWv7du3Y8+ePRg3blzYa+Tm5sJmswUtg+Fb9dQjMPMQybhx41BVVYXXX3/db//zzz/vfT0RFEWBxWLx+wW7tLQ0ZP9DCfeeRo8ejXbt2mH37t0YNmxYyM2TaUsEh8OBVatWoVOnTjj11FO9+3/xi1/ggw8+QHFxsXdfZWUl1qxZgyuuuMJb0fecc85Bfn5+UPXd1157DVVVVWEDNo94/m0bq1+/fujbty+++OKLsPfWd5izR25uLqZPn45rr70W3377LWpqaiL2/bLLLoOI4ODBgyGv4Qm+PSorK/HGG2/47Xv55ZehqmrI4ZOtW7fGxIkTce+998Jut+Prr7+O+16MGzcOu3fvxueff+63//nnn4eiKBg7dmzc5wTgPe6ll17y2//yyy/7fd+qVSuMHTsWO3fuxKBBg0Lep1CjBWKlKAoGDx6Mxx9/HO3atQt6n76iPYNDhw7FqFGj8PDDD+Oll17C9OnT0bp16wb3jYiouWFGkohiNmHCBHTt2hWXX345+vfvD8MwsGvXLjz66KNo06YN7rjjDgCuX8x/97vf4cknn4Sqqpg4cSL279+PP/7xj+jWrRvmzp0b9hqeOWTLly9Hnz59MHjwYPzv//5v0C+cQH3W629/+xumTZsGs9mMfv36hfylf+rUqXj66acxbdo07N+/H2eeeSY+/vhjPPTQQ7jkkkvw85//PCH36LLLLsOaNWswc+ZMTJo0CcXFxfjv//5v5OfnxzT8ceDAgQCAZ599Fm3btoXNZkOvXr3QsWNHPPnkk5g2bRqOHTuGSZMmoXPnzjhy5Ai++OILHDlyBEuWLGlQn+fNmweHw4HRo0cjLy8PxcXFePLJJ7Fr1y4899xz0DTN2/bOO+/ECy+8gEsvvRQPPPAArFYr/vznP6O2tta7PAcAaJqGRx55BFOmTMEtt9yCa6+9Fnv37sVdd92Fiy66CBdffHHEPvXp0wdZWVl46aWXMGDAALRp0wYFBQXeP1Yk2jPPPIOJEydiwoQJmD59Orp06YJjx45hz549+Pzzz/Hqq68CcAXIl112GQYNGoT27dtjz549eOGFFzBy5EjvPEzPc/nwww9j4sSJ0DQNgwYNwujRo/G73/0O119/PXbs2IHzzz8frVu3RklJCT7++GOceeaZ3rl3gCuLd+utt6KoqAinnXYa1q9fj7///e+49dZb0b17dwDAzTffjKysLIwePRr5+fkoLS1FYWEhcnJyvMtXxGPu3Ll4/vnnvf++PXr0wNtvv43Fixfj1ltvxWmnndag+zt+/Hicf/75uOuuu1BdXY1hw4bh3//+N1544YWgtn/7299w7rnn4rzzzsOtt96Knj17orKyEt9//z3efPPNoHma0bz11ltYvHgxrrrqKvTu3RsigjVr1uD48eO46KKLwh4XyzN4xx13YPLkyVAUxTsMmoioxUhfnR8iam5WrVol1113nfTt21fatGkjZrNZunfvLlOmTPEroy/iqib68MMPy2mnnSZms1k6deokv/3tb71LOHgEVm0VcZXev+mmmyQ3N1dat24tl19+uezfvz+ouqOIyIIFC6SgoEBUVfWrdhpYtVVE5OjRozJjxgzJz88Xk8kkPXr0kAULFniXZvBAwBILHoFVHcP585//LD179hSr1SoDBgyQv//97yErjYaq2ioismjRIunVq5domhZUNXLz5s1y6aWXSocOHcRsNkuXLl3k0ksvlVdffdXbxnMt32VWIlm2bJkMHz5cOnToICaTSdq3by8TJkyQ9957L2T777//Xq666irJzs6WVq1aybhx4+Szzz4L2fbll1+WQYMGicVikby8PLn99tv9KvxGsnLlSunfv7+YzWa/f/twVVsDK5CKhP639FQc/ctf/uK3/4svvpBrrrlGOnfuLGazWfLy8uTCCy+UpUuXetvMnz9fhg0bJu3btxer1Sq9e/eWuXPnSnl5ubdNXV2d3HTTTXLKKaeIoihBlYWXL18u55xzjrRu3VqysrKkT58+MnXqVNmxY4e3jefZ2LRpkwwbNkysVqvk5+fLPffc41dZ9p///KeMHTtWcnNzxWKxSEFBgVxzzTXy5ZdfRr2/4e7ZgQMH5LrrrvMuB9OvXz/5y1/+4q1cG+keRnL8+HG54YYbpF27dtKqVSu56KKL5Jtvvgn5c71v3z654YYbpEuXLmI2m+WUU06RUaNGyYMPPuht46na6vvs+/bN83PzzTffyLXXXit9+vSRrKwsycnJkeHDh8uKFSuC7kfgz3e4Z9Cjrq5OrFarXHzxxTHfByKiTKGIuCfBEBERUZMwZswYlJeXe+cdU9P05ptv4oorrsDbb7/tLQZERNRScGgrERERURx2796NAwcO4Pe//z3OOussTJw4Md1dIiJKORbbISIiIorDzJkzccUVV6B9+/ZYuXJlUtZFJSJq6ji0lYiIiIiIiOKS1ozkli1bcPnll6OgoACKogSV5Q9l8+bNGDp0KGw2G3r37o2lS5cmv6NERERERETkldZAsrq6GoMHD8ZTTz0VU/t9+/bhkksuwXnnnYedO3finnvuwe23347Vq1cnuadERERERETk0WSGtiqKgrVr1+Kqq64K2+buu+/GG2+84V38HABmzJiBL774Atu2bUtBL4mIiIiIiKhZVW3dtm0bxo8f77dvwoQJWLZsGRwOB8xmc9AxdXV1qKur835vGAaOHTuGjh07cnI8EREREVGaiQgqKytRUFAAVW1+tUBra2tht9tjamuxWGCz2ZLco9RoVoFkaWkpcnNz/fbl5ubC6XSivLwc+fn5QccUFhbi/vvvT1UXiYiIiIioAYqLi9G1a9d0dyMutbW16NWjDUrL9Jja5+XlYd++fRkRTDarQBJAUBbRMzI3XHZxwYIFmDdvnvf7EydOoHv37jjjN3+CZmn+/4BE1EIIoAig6ALrCQNZb30OiJHuXhERETWaEw58jPVo27ZtursSN7vdjtIyHfs+64HstpGzqRWVBnoNPQC73c5AMtXy8vJQWlrqt6+srAwmkwkdO3YMeYzVaoXVag3ar1lsDCSJqOlzB5AAoDoFWT/paP1VCZzQAEVLb9+IiIgSwf3/c8152llWG0FWm8ilZxxNozRNwjSrQHLkyJF48803/fZt2LABw4YNCzk/koioWfMJIkUBVAdgO3wSRsnh9PaLiIiI/BgwEG2cUPQWzUtaZ7NWVVVh165d2LVrFwDX8h67du1CUVERANew1KlTp3rbz5gxAwcOHMC8efOwZ88eLF++HMuWLcOdd96Zju4TESVPQBAJAFnlTmgHDsOIcUI/ERERpYYuEtOWSdKakdyxYwfGjh3r/d4zl3HatGlYsWIFSkpKvEElAPTq1Qvr16/H3Llz8fTTT6OgoABPPPEErr766pT3nYgomRQJ/t5WUgX96E/p6RARERGFZUBgIHKgGO315iatgeSYMWMQaRnLFStWBO274IIL8PnnnyexV0REaebzsejJRmp2AYpKIE5HevpEREREYRkQ6AwkiYgobUIEkaouMFcZ0E9UpKdPREREFBEzkkRElD6eqnVSH0QCriI7WUfsXO6DiIioiYplDiTnSBIRUeJJ8LxIKIDqENiO6zAXlcOZlo4RERFRNIZ7i9YmkzCQJCJKt4AAUhQA3rmRXPKDiIioqdNjmCMZ7fXmJq3LfxARtXgBcyJ9g0jFANp/+RPw2W4u+UFERNSEOSS2LR6FhYU4++yz0bZtW3Tu3BlXXXUVvv32W782IoKFCxeioKAAWVlZGDNmDL7++uuo5169ejVOP/10WK1WnH766Vi7dm18nQMDSSKi9An1fyg+cyO9lVp1PWVdIiIiovgZUKBH2Qzf/5OPwebNmzFr1ix88skn2LhxI5xOJ8aPH4/q6mpvm0ceeQSPPfYYnnrqKWzfvh15eXm46KKLUFlZGfa827Ztw+TJkzFlyhR88cUXmDJlCq655hp8+umncfVPkUjrb2SgiooK5OTkYND1D0Gz2NLdHSJqqaIEkYoOZB3VkbVuO4vsEBFRRnOKA5uwDidOnEB2dna6uxMXT2yx4+tctGkbOUdXVWlg2BmHG/w+jxw5gs6dO2Pz5s04//zzISIoKCjAnDlzcPfddwMA6urqkJubi4cffhi33HJLyPNMnjwZFRUVeOedd7z7Lr74YrRv3x4rV66MuT/MSBIRNUGaXWArr2MQSURE1AxEy0Z6NsAVfPpudXV1MV3jxIkTAIAOHToAAPbt24fS0lKMHz/e28ZqteKCCy7A1q1bw55n27ZtfscAwIQJEyIeEwoDSSKiVIuSjQSArGM6zPuPpKQ7RERE1DjxBJLdunVDTk6OdyssLIx6fhHBvHnzcO6552LgwIEAgNLSUgBAbm6uX9vc3Fzva6GUlpbGfUworNpKRJRKodaKDJwyIYCttAZGGQNJIiKi5sAQBYZEngPpeb24uNhvaKvVao16/tmzZ+PLL7/Exx9/HPSaovhfV0SC9iXimEAMJImIUiWWIBKA5hCoRaVwslIrERFRs+CbcYzUBgCys7PjmiN522234Y033sCWLVvQtWtX7/68vDwArgxjfn6+d39ZWVlQxtFXXl5eUPYx2jGhcGgrUUsmCD3MkhLPJ4iMRDEAc5VAP/ZT8vtERERECaFDjWmLh4hg9uzZWLNmDT744AP06tXL7/VevXohLy8PGzdu9O6z2+3YvHkzRo0aFfa8I0eO9DsGADZs2BDxmFCYkSRqqRhApk5AEBkxG2kXZB2xc8kPIiKiZkRiGNoqUV4PNGvWLLz88stYt24d2rZt680i5uTkICsrC4qiYM6cOXjooYfQt29f9O3bFw899BBatWqF6667znueqVOnokuXLt65mHfccQfOP/98PPzww7jyyiuxbt06vP/++yGHzUbCQJKoJfINIuP7TKN4SYgAEgh53801guz/q4H6+bdgrVYiIqLmwy4azBI542iPM5BcsmQJAGDMmDF++5977jlMnz4dAHDXXXfh5MmTmDlzJn766Secc8452LBhA9q2bettX1RUBFWt79uoUaPwyiuv4L/+67/wxz/+EX369MGqVatwzjnnxNU/riNJ1NIE/sQzkEyeODKRANC22IFWXx6E8+ChpHeNiIioqciEdSTf/rI3WrfVIratrtRx6aAfmuX7DIUZSaKWjEFk8gQE7NGCSFel1mroh1mplYiIqLmJp9hOpmAgSdSScEhrWkQNIuGaG6kWHYbT6UhJn4iIiChxdFGhRxnaqmfYQFAGkkREiRbn/08oOmCpEujHjyelO0RERJRcBhQYUf5KH+315oaBJFFLkVl/BGu6Qt3nKP+/odkFtiN1rNRKRETUTBkxLO9hZNgvYwwkiVoCFthJjQYEkapDYD2hw/LjT3AmpVNERESUbBzaSkSZj0Fkk6LZAVtZHYyDpenuChERETWQARUGM5JERBS3BmQjAcB21AnTgcNw1tUmvEtERESUGroo0KOsExnt9eaGgSRRJuOQ1uQL98fFGO61qVaQ9cleOFlkh4iIqFlziAkOibyOpIOBJBERAWhUEKnogKXSgH6iIqFdIiIiotTTYyi2o3NoKxE1C8xGJlcjgkjAXam13A6IkbAuERERUXoYiD50NdP+H5+BJBFRosQRrNt+0mEuKmelViIiogwQW7GdyK83NwwkiVoCZiMTq7EjUwSwlZ6EUXI4Id0hIiKi9Ipt+Q8GkkTU1GXWEPympZFDWgHXsFat+DCcdntCukRERETpZUCBEeWXgWivNzcMJIkyDYPI5ElAEKkYgLlGoJcfTUiXiIiIKP2YkSQiovjE+cdFzS7IOmKH6Hpy+kNEREQpF1vVVgaSRNRUhcqYZdYoivRJQKZXdQosFQYsPx5nkR0iIqIMYogCI1rVVq4jSUTNRmZ9XqVPAoa0AoBmB7LK6mAcLGl0l4iIiKjpMGLISLJqKxE1TQIo7oBHFDCITJREZXkFaLenEvjqOxgsskNERJRRHKJBEy1Km8wqZJFZYTERUSIlcKiwZheoRaUMIomIiDKQIWpMW7y2bNmCyy+/HAUFBVAUBa+//rrf64qihNz+8pe/hD3nihUrQh5TW1sbV9+YkSTKBJn1B66mIYH3VNEBc7VAP/ZT4k5KRERETYYOQI/y1+aGlNmrrq7G4MGDcf311+Pqq68Oer2kxH+6zDvvvIMbb7wxZFtf2dnZ+Pbbb/322Wy2uPqW9ozk4sWL0atXL9hsNgwdOhQfffRRxPYvvfQSBg8ejFatWiE/Px/XX389jh5lGX0iv2Gt1DgJmhPpwUqtREREmS1ZGcmJEyfiwQcfxC9/+cuQr+fl5flt69atw9ixY9G7d++I51UUJejYeKU1kFy1ahXmzJmDe++9Fzt37sR5552HiRMnoqioKGT7jz/+GFOnTsWNN96Ir7/+Gq+++iq2b9+Om266KcU9J2riFJ+N4pPgIFJ1CqwVBiw/MhtJRESUqTzrSEbbAKCiosJvq6urS0gfDh8+jLfffhs33nhj1LZVVVXo0aMHunbtissuuww7d+6M+3ppDSQfe+wx3HjjjbjpppswYMAALFq0CN26dcOSJUtCtv/kk0/Qs2dP3H777ejVqxfOPfdc3HLLLdixY0eKe07UhAQW2aGGS3AQCbgqtdrKamH8yEqtREREmUqgwIiyifsXim7duiEnJ8e7FRYWJqQP//znP9G2bduw2UuP/v37Y8WKFXjjjTewcuVK2Gw2jB49Gnv37o3remmbI2m32/HZZ59h/vz5fvvHjx+PrVu3hjxm1KhRuPfee7F+/XpMnDgRZWVleO2113DppZeGvU5dXZ1flF9RUZGYN0DUVPlmIjl3Mu1sR50wHSiDsy6+CexERETUfPhmHCO1AYDi4mJkZ2d791ut1oT0Yfny5fjNb34Tda7jiBEjMGLECO/3o0ePxs9+9jM8+eSTeOKJJ2K+XtoykuXl5dB1Hbm5uX77c3NzUVpaGvKYUaNG4aWXXsLkyZNhsViQl5eHdu3a4cknnwx7ncLCQr+Iv1u3bgl9H0RpxUAxcZKQjYQAttJq6GXljTgJERERNXWGKDFtgKvQje+WiEDyo48+wrffftugKX+qquLss8+OOyOZ9mI7iuL/W5qIBO3z2L17N26//Xb86U9/wmeffYZ3330X+/btw4wZM8Kef8GCBThx4oR3Ky4uTmj/idJN8Q2AOC+yYZIRRMJVZEc5UAJxOhp3IiIiImrSdKgxbcmybNkyDB06FIMHD477WBHBrl27kJ+fH9dxaRva2qlTJ2iaFpR9LCsrC8pSehQWFmL06NH4wx/+AAAYNGgQWrdujfPOOw8PPvhgyDdvtVoTli4malLEP4iUwCCS2crYJCmINJ0UtD1QyyU/iIiIWgCnaNBEi9LGiPu8VVVV+P77773f79u3D7t27UKHDh3QvXt3AK6pe6+++ioeffTRkOeYOnUqunTp4p2Lef/992PEiBHo27cvKioq8MQTT2DXrl14+umn4+pb2jKSFosFQ4cOxcaNG/32b9y4EaNGjQp5TE1NDVTVv8ua5voHE+FvzdSyKHzkGy9JQSQA2I7rsBQda/yJiIiIqMnTRYlpi9eOHTswZMgQDBkyBAAwb948DBkyBH/605+8bV555RWICK699tqQ5ygqKvJbb/L48eP43e9+hwEDBmD8+PE4ePAgtmzZguHDh8fVN0XSGIGtWrUKU6ZMwdKlSzFy5Eg8++yz+Pvf/46vv/4aPXr0wIIFC3Dw4EE8//zzAIAVK1bg5ptvxhNPPIEJEyagpKQEc+bMgaqq+PTTT2O6ZkVFBXJycjDo+oegWeJbdJOoyQiXjWRGMnZJDCIhQKcvq6Ds+g4Gi+wQERFF5BQHNmEdTpw44VeEpjnwxBa3bLka1jbmiG3rqhx45vzVzfJ9hpK2oa0AMHnyZBw9ehQPPPAASkpKMHDgQKxfvx49evQAAJSUlPitKTl9+nRUVlbiqaeewu9//3u0a9cOF154IR5++OF0vQWipolBZGTJDCIBaA6BVsRKrURERC2FiAojStVWifJ6c5PWjGQ6MCNJzZ77J9YvI6mC2chYRbo3CQgkFQOwVBpou/ozFtkhIiKKQSZkJG/cfA0sUTKS9ioHll3wr2b5PkNJa0aSiBom4vxIBpHhJTmIBFyVWrPKHQwiiYiIWhBD4F3eI1KbTMJAkqiZC6rWSvFL0P1TdFc20vrjCTgTc0oiIiJqBowYhrZGe725YSBJ1JwIq7U2WArum2YX2MrtMIoPJf9iRERE1GQYUGBE+ct0tNebGwaSRM0Y146MUQqGtAKA7Scd5qJyOGtqEndSIiIiavJiWd6jIct/NGUMJImai1iW/KBgSa7Q6qHVCVpv3w/n4bLEnpiIiIiaPKdoUA0taptMwkCSqDkIFQwxgIwuRUGkYgDmGoFefjSxJyYiIqJmQWIY2ioZ9ssbA0miZiLq3EgOa/WXoiAS8MyNdEB0PfEnJyIioibPECWGqq0MJIkozYTzIhsm1Oe3gkbdQ9UpsFQYsP54nJVaiYiIWihWbSWipoeVWuMT771q5L3V7ICtvA7Gj6zUSkRE1FIxI0lETR7XjYwghcNZPWw/6TDvPwLnyZPJuwgRERE1aVz+g4iaFmYimzYBbKU1MMqOpLsnRERElEbMSBJRkxM0rDWzPoMSJ9TyKEBS75dmF6jFh+G025N3ESIiImryGEgSUdMRIhuZYZ8/iRNwr8Lep0YW1vE7lQGYqwX60WOJOSERERE1WwwkiajJCMxEcm5kGPFkIhM4VLjVER1t/lMGJ5f8ICIiavF0UaBEqcqqM5AkoqTj3MjYpGE4q+e6trJaGD+WJPlCRERE1BwwI0lETRKzkSGEWxYlBfdJcwhMRUfgrKtN/sWIiIioyWuJgWRmrYpJlAmYjWy4VHw+C6DVCfTDrNRKRERELp5AMtoWry1btuDyyy9HQUEBFEXB66+/7vf69OnToSiK3zZixIio5129ejVOP/10WK1WnH766Vi7dm3cfWMgSdQEsVJrdN57FOpeJTGDq9kFWUecEKcjORcgIiKiZidZgWR1dTUGDx6Mp556Kmybiy++GCUlJd5t/fr1Ec+5bds2TJ48GVOmTMEXX3yBKVOm4JprrsGnn34aV984tJWoicuwURCJESKIDHmfElilFQAUHbBUCawHK8ASO0REROQhokCi/NIW7fVQJk6ciIkTJ0ZsY7VakZeXF/M5Fy1ahIsuuggLFiwAACxYsACbN2/GokWLsHLlypjPw4wkUVMSbt4fBQt3n5J4/zS7wHakDig+lLyLEBERUbNjQIlpA4CKigq/ra6urlHX3rRpEzp37ozTTjsNN998M8rKyiK237ZtG8aPH++3b8KECdi6dWtc12UgSdSExVxkpyUV44mlyI74bAlkPaHDUnQMelVVYk9MREREzVo8Q1u7deuGnJwc71ZYWNjg606cOBEvvfQSPvjgAzz66KPYvn07LrzwwojBaWlpKXJzc/325ebmorS0NK5rc2grUVMRat3IRhyfkTxBZGPvVQOvnXW4FkbJ4RRcjIiIiJqTeIa2FhcXIzs727vfarU2+LqTJ0/2fj1w4EAMGzYMPXr0wNtvv41f/vKXYY9TFP++ikjQvmgYSBI1IQ0ustMSgshwUpSN1RwCraiMS34QERFRkHiW/8jOzvYLJBMpPz8fPXr0wN69e8O2ycvLC8o+lpWVBWUpo+HQVqKmINRwzZYyVDUOQdlIJTXZSFOtoG2xHc6S+IZ8EBERUctgGCr0KJthJD/0Onr0KIqLi5Gfnx+2zciRI7Fx40a/fRs2bMCoUaPiuhYzkkRNECu1hhAiiEwVc5UB68EKOFN3SSIiImpGBIBEGSHWkAFkVVVV+P77773f79u3D7t27UKHDh3QoUMHLFy4EFdffTXy8/Oxf/9+3HPPPejUqRN+8YtfeI+ZOnUqunTp4p2Leccdd+D888/Hww8/jCuvvBLr1q3D+++/j48//jiuvjGQJEq3RFRqTfAyF01OmHsUczGiRl4764gdRtHBJF+IiIiImisDCpQov5QYDfilZceOHRg7dqz3+3nz5gEApk2bhiVLluCrr77C888/j+PHjyM/Px9jx47FqlWr0LZtW+8xRUVFUNX6bOioUaPwyiuv4L/+67/wxz/+EX369MGqVatwzjnnxNU3BpJETUyDgqNMDiIRusCOfwP3f5NwH1SnwFxUDmdNTeJPTkRERBkhWetIjhkzBhIh1fnee+9FPcemTZuC9k2aNAmTJk2Kuz++GEgSpVOGB4AJE2s2MtGZWQFMdWClViIiIorIEAVKjMV2MgUDSaI0Y5GdKKIFhkkMJjWHwHbUCcNuT8wJiYiIKCOJxDBHMsMSCAwkidIlXJaN/IScP5qCuZGKAZhOCmyHKqEn91JERETUzCVraGtTxkCSKI0aXWQn04ULIqNJQFZSswuyyhxAMZf8ICIiosgYSBJR2qSkAmmGSMXnsKXSgPXH43CeqEj+xYiIiKhZ4xxJIkqNgOUsGEQGUwz3FwFrR3rvVbj7lYAsr+oQtP3sIJxFPzb+ZERERJTxDANQjCiBpBHx5WaHgSRRqnE4a+wC7lXEgDtR99VTqbW0LEEnJCIiokzHoa1ElBKs1BpFqHUjUxFEwjU30lbuYKVWIiIiilm0Ja8Rw+vNjZruDixevBi9evWCzWbD0KFD8dFHH0VsX1dXh3vvvRc9evSA1WpFnz59sHz58hT1lijxMuyPU82aYgDmGlelViIiIqJYeTKS0bZMktaM5KpVqzBnzhwsXrwYo0ePxjPPPIOJEydi9+7d6N69e8hjrrnmGhw+fBjLli3DqaeeirKyMjidzhT3nKiBhJVaowp1j5TUBNyaXZB1xM5KrURERBSfFpiSTGsg+dhjj+HGG2/ETTfdBABYtGgR3nvvPSxZsgSFhYVB7d99911s3rwZP/zwAzp06AAA6NmzZyq7TJRQLLITm5BBZBI+jC0VBiw/HoezghlJIiIiikMsGccMy0imbWir3W7HZ599hvHjx/vtHz9+PLZu3RrymDfeeAPDhg3DI488gi5duuC0007DnXfeiZMnT4a9Tl1dHSoqKvw2orTIsL9CJYUnGxlQqdXvv0m8dlZZHYyDJYBkWFk1IiIiSiqR2LZMkraMZHl5OXRdR25urt/+3NxclJaGHlb2ww8/4OOPP4bNZsPatWtRXl6OmTNn4tixY2HnSRYWFuL+++9PeP+JGoJLfjRQKoa1OgSmA2VwRvjDFBEREVEoLbFqa9qL7SiK/w0VkaB9HoZhQFEUvPTSSxg+fDguueQSPPbYY1ixYkXYrOSCBQtw4sQJ71ZcXJzw90AUVYb9BSpZQmUjU/KZK4BWJ9APH0nBxYiIiCjjiBLblkHSlpHs1KkTNE0Lyj6WlZUFZSk98vPz0aVLF+Tk5Hj3DRgwACKCH3/8EX379g06xmq1wmq1JrbzRA3AbGQEYYa0puo+ZR0z0Hb3MehOR/IvRkRERBlHjOgzYzJt5kzaMpIWiwVDhw7Fxo0b/fZv3LgRo0aNCnnM6NGjcejQIVRVVXn3fffdd1BVFV27dk1qf4kajJVamzxbuR0oPpTubhAREVEz1RKX/0jr0NZ58+bhH//4B5YvX449e/Zg7ty5KCoqwowZMwC4hqVOnTrV2/66665Dx44dcf3112P37t3YsmUL/vCHP+CGG25AVlZWut4GUVyYjQwWlI1E6kZ/qE6BpfgYdJ8/UBERERHFTaJsGSaty39MnjwZR48exQMPPICSkhIMHDgQ69evR48ePQAAJSUlKCoq8rZv06YNNm7ciNtuuw3Dhg1Dx44dcc011+DBBx9M11sgiiwDPzQSLtKHawoqtWp1gHGQ60YSERFRw7HYThrMnDkT+/fvR11dHT777DOcf/753tdWrFiBTZs2+bXv378/Nm7ciJqaGhQXF+PRRx9lNpKatKBhrZn1GZIcKbpHmkNgO+qEUVebmgsSERFRZoqWjWxgVnLLli24/PLLUVBQAEVR8Prrr3tfczgcuPvuu3HmmWeidevWKCgowNSpU3HoUOTpOitWrICiKEFbbW18vw+lPZAkylghPiwy7A9RjRdp/mgKspGmWoGtpDrJFyIiIqLMp8S4xae6uhqDBw/GU089FfRaTU0NPv/8c/zxj3/E559/jjVr1uC7777DFVdcEfW82dnZKCkp8dtsNltcfUvr0FaiTMciOw2QourYml2QdcQBpbgk+RcjIiKizBZLxrEBvxdOnDgREydODPlaTk5OUOHSJ598EsOHD0dRURG6d+8e9ryKoiAvLy/+DvlgRpIoRVhkJ0CoJT9SyFJpwHqwAvpPJ9LTASIiIsoccQxtraio8Nvq6uoS1o0TJ05AURS0a9cuYruqqir06NEDXbt2xWWXXYadO3fGfS0GkkTJwCU/GiaFa0fayu2QH0syb1EnIiIiSj1RYtsAdOvWDTk5Od6tsLAwIV2ora3F/Pnzcd111yE7Oztsu/79+2PFihV44403sHLlSthsNowePRp79+6N63oc2kqUaOHmRjIb6Sfskh8pmBtprhGYPt0DnUV2iIiIKAHEiP63ac/rxcXFfoGe1Wpt9PUdDgd+/etfwzAMLF68OGLbESNGYMSIEd7vR48ejZ/97Gd48skn8cQTT8R8TQaSREnASq1RhBrSmsJKrVnlrNRKRERECSQxFHlwv56dnR0xYxgvh8OBa665Bvv27cMHH3wQ97lVVcXZZ58dd0aSQ1uJkoyVWoOlq1KrYgCmkwJbSVVyL0REREQtiiKxbYnmCSL37t2L999/Hx07doz7HCKCXbt2IT8/P67jmJEkSiTOjYwu1P1JZaXWcgdQxEqtRERElEBJqtpaVVWF77//3vv9vn37sGvXLnTo0AEFBQWYNGkSPv/8c7z11lvQdR2lpaUAgA4dOsBisQAApk6dii5dunjnYt5///0YMWIE+vbti4qKCjzxxBPYtWsXnn766bj6xkCSKIk4N7JpsVQasP54As4TFenuChEREWWSOIa2xmPHjh0YO3as9/t58+YBAKZNm4aFCxfijTfeAACcddZZfsd9+OGHGDNmDACgqKgIqlo/EPX48eP43e9+h9LSUuTk5GDIkCHYsmULhg8fHlffGEgSJQozkdGFWfIjVUV2bOV2GMWHWKmViIiIEitJGckxY8ZAJPyBkV7z2LRpk9/3jz/+OB5//PH4OxOAgSRRAvkOa2U2MljIdSNTdJ9Up8BcVA5nTU3yL0ZEREQtS5ICyaaMgSRRImTYB0PGEcBUBxglh9PdEyIiIspEDCSJqKG45EcUYbKR3ukCSkDbBNLsAlu5A4bdntgTExEREQFJmyPZlDGQJEqCDPucaLw0V7M117iW/NDT1wUiIiLKYLEs75Fplf3jXkdy+vTp2LJlSzL6QtQ8CedGRhPqg1NU1+a9V4LYhoXESXUKcnYdgf6f7xJ7YiIiIiIPiXHLIHEHkpWVlRg/fjz69u2Lhx56CAcPHkxGv4iahwz7QEiqdNwrATQ7YPzISq1ERESUPArqs5Jht3R3MsHiDiRXr16NgwcPYvbs2Xj11VfRs2dPTJw4Ea+99hocDkcy+kjUpHFuZAOkKGurOQS2n3QYJ08m/2JERETUcnnmSEbbMkjcgSQAdOzYEXfccQd27tyJ//3f/8Wpp56KKVOmoKCgAHPnzsXevXsT3U+ipifUcM3M+nxIjMChHKm6R+5spK2Uy30QERFRknFoa3xKSkqwYcMGbNiwAZqm4ZJLLsHXX3+N008/PSGLXBI1dZk2aTrhwkw8T0XArdkFWUccUIu55AcRERElGQPJ6BwOB1avXo3LLrsMPXr0wKuvvoq5c+eipKQE//znP7Fhwwa88MILeOCBB5LRX6Imi0V2AoQKIpXU3SdztcBaUgn96LHkX4yIiIhatKjzI9NcwT4Z4l7+Iz8/H4Zh4Nprr8X//u//4qyzzgpqM2HCBLRr1y4B3SNqojLwwyBp0nSfssrtwIFDEJ2LfhAREVGSxZJxzLDfHeMOJB9//HH86le/gs1mC9umffv22LdvX6M6RtScMBsZTAnxgZqq+6Q6BZbin+Csqk7+xYiIiIgYSEY3ZcqUZPSDqPnIsA+BpEjnPADvkh8lXPKDiIiIUiKWoauZNpot7kCSiPw/CJiNDBDugzQV90kAS5WBnG8rYdTVJvliRERERG6G4tqitckgDCSJ4sG5kQ2ToqWTXJVanVAPlIK5SCIiIkoVZiSJKC7MRgaQ0HMjU8VSJbAerIB+/Hh6OkBEREQtE+dIElFYGfbDn0qpCrht5XagmJVaiYiIKMViGbWWYb9LMpAkigPnRkYWMhuZonvkqtR6DM6qqtRckIiIiMijBWYk1XR3gKhZyLAf/JRLQZEdrQ4wDpYm+UJEREREIUiMW5y2bNmCyy+/HAUFBVAUBa+//rr/ZUWwcOFCFBQUICsrC2PGjMHXX38d9byrV6/G6aefDqvVitNPPx1r166Nu28MJIliFDRcgdlIf2GykSkpsuMQ2I46WamViIiI0sJTbCfaFq/q6moMHjwYTz31VMjXH3nkETz22GN46qmnsH37duTl5eGiiy5CZWVl2HNu27YNkydPxpQpU/DFF19gypQpuOaaa/Dpp5/G1TcObSVqgFQERxlNQUKzvKZaga2kmpVaiYiIKKNMnDgREydODPmaiGDRokW499578ctf/hIA8M9//hO5ubl4+eWXccstt4Q8btGiRbjooouwYMECAMCCBQuwefNmLFq0CCtXroy5b8xIEkXDJT+iC3WPlNTNI8064oBSXJL8CxERERGFEsfQ1oqKCr+trq6uQZfct28fSktLMX78eO8+q9WKCy64AFu3bg173LZt2/yOAYAJEyZEPCYUBpJEkYQIIFlkJ1ioIjsR71OiAnMBTCcF5m17oB/7KUEnJSIiIoqPIoBiRNncv/9069YNOTk53q2wsLBB1ywtddWGyM3N9dufm5vrfS3ccfEeEwqHthJFwbmRUYSr1JqC+6Q6BbbjOoyamuRfjIiIiCicOKq2FhcXIzs727vbarU26tKK4v9Ll4gE7UvEMYEYSBKFEy4bSZEpYb5O9PBgATQ7YDt8kkV1iYiIKK1iKabjeT07O9svkGyovLw8AK4MY35+vnd/WVlZUMYx8LjA7GO0Y0JJ+9DWxYsXo1evXrDZbBg6dCg++uijmI7797//DZPJhLPOOiu5HaQWjXMjowjzoRky4E5wEK45BFnlTmgHDif2xERERETxStLyH5H06tULeXl52Lhxo3ef3W7H5s2bMWrUqLDHjRw50u8YANiwYUPEY0JJa0Zy1apVmDNnDhYvXozRo0fjmWeewcSJE7F7925079497HEnTpzA1KlTMW7cOBw+zF8iKTU4NzKAhJ4bmSqmkwLboUro5UfT0wEiIiIit3gykvGoqqrC999/7/1+37592LVrFzp06IDu3btjzpw5eOihh9C3b1/07dsXDz30EFq1aoXrrrvOe8zUqVPRpUsX71zMO+64A+effz4efvhhXHnllVi3bh3ef/99fPzxx3H1La0Zycceeww33ngjbrrpJgwYMACLFi1Ct27dsGTJkojH3XLLLbjuuuswcuTIFPWUWhxWam2YVFZqLXMARSUQXU/+xYiIiIgiSVJGcseOHRgyZAiGDBkCAJg3bx6GDBmCP/3pTwCAu+66C3PmzMHMmTMxbNgwHDx4EBs2bEDbtm295ygqKkJJSX11+1GjRuGVV17Bc889h0GDBmHFihVYtWoVzjnnnLj6lraMpN1ux2effYb58+f77R8/fnzE0rPPPfcc/u///g8vvvgiHnzwwajXqaur8yupW1FR0fBOU4vFbGSwmIa0JikYV3TA+uNxOCvCL7ZLRERElDJxFNuJx5gxYyAS/kBFUbBw4UIsXLgwbJtNmzYF7Zs0aRImTZoUf4d8pC0jWV5eDl3X4yo9u3fvXsyfPx8vvfQSTKbYYuDCwkK/8rrdunVrdN8pwzETGbtQ9yrZAbcAml1g/HgIECPJFyMiIiKKzjO0NdqWSdJebCfW0rO6ruO6667D/fffj9NOOy3m8y9YsAAnTpzwbsXFxY3uM2U+LvnRAClc8iPrmA7j5MnkX4yIiIgoFmkotpNuaRva2qlTJ2iaFnPp2crKSuzYsQM7d+7E7NmzAQCGYUBEYDKZsGHDBlx44YVBx1mt1kavzUItSIi/FnHJjwChiuwoqblPWp2gzSEHbF8cgDP5lyMiIiKKiWK4tmhtMknaAkmLxYKhQ4di48aN+MUvfuHdv3HjRlx55ZVB7bOzs/HVV1/57Vu8eDE++OADvPbaa+jVq1fS+0wtD+dGBoi03EcK7pO5WmAtqYR+9FjyL0ZEREQUqyTNkWzK0rr8x7x58zBlyhQMGzYMI0eOxLPPPouioiLMmDEDgGtY6sGDB/H8889DVVUMHDjQ7/jOnTvDZrMF7SdqkAz74U6qNN2rrHI7cOAQK7USERFRk5Ks5T+asrQGkpMnT8bRo0fxwAMPoKSkBAMHDsT69evRo0cPAEBJSQmKiorS2UVqYXx/wJmNDBZy3cgUzo20/Hgczqrq5F+MiIiIKB4tMCOpSKR6shmooqICOTk5GHT9Q9AstnR3h5oK908BA8kIJMzYfgWQZJftEsB8UtDu1Z0w6mqTfDEiIiJKJac4sAnrcOLECWRnZ6e7O3HxxBYDZj4EzRo5ttDrarFn8T3N8n2GktaMJFFTwkqtDRBYZMfzdYL/PKU5BLajTgaRRERE1CTFkn/ItF8tGUgShcBKrQHSvPaRViewHa5BhhU7IyIiokzRAoe2MpAkysAFYlMi1JIfgvo/ySXwnmYdcUItPsxAkoiIiJokFtshIs6NDCFkkR3viwHfJziYVHTAWlIB/dhPjT8ZERERUTIwI0nUwmTYD3QqRQy4JcJrcdLsAhSVcMkPIiIiatpa2O+VDCSp5WKl1tg0dMmPxn6YCmCuEWT/Xw30yspGnoyIiIgoeRQjTHX7gDaZhIEktWiZNlY94UKN9w81NzIJPJVaTUVlcCb/ckREREQNxjmSRC0Ys5HBYv7A871vCfqQNNUKbCXV0MvKE3NCIiIiomThHEmiFoKVWqOTgP8C9dnIFATcWUccUIpLYDgdyb8YERERUSMwI0nUQjEb2UgJLK4DuCu1HqyA86cTiTspERERUbIwI0nUAmTYD3GyhFryI2ql1kQQV6VWo+ggIBk2K52IiIgyUwsMJNV0d4AoHUIVkCEfDa3UmgCqU2A7rsOoqUn+xYiIiIgSwDO0NdoWj549e0JRlKBt1qxZIdtv2rQpZPtvvvkmAe8wGDOS1OKlogIpxUgAUx1gO3wy0/5oR0RERJksCRnJ7du3Q/dZR/s///kPLrroIvzqV7+KeNy3336L7Oxs7/ennHJKfBeOEQNJanEybaJzwkVb8iMJFVo9PEt+aAcOc8kPIiIiajYUESgS+RejaK8HCgwA//znP6NPnz644IILIh7XuXNntGvXLq5rNQSHtlLLEe+cv5ZIYpgbGaqaa4K0PWCH7ZPv4DxclviTExERESWJYsS2AUBFRYXfVldXF/X8drsdL774Im644QYoSuRfXocMGYL8/HyMGzcOH374YSLeXkgMJKlF8c20MYiMUaj7lIQgUnUKrD8eh15RmfiTExERESWTxLgB6NatG3JycrxbYWFh1NO//vrrOH78OKZPnx62TX5+Pp599lmsXr0aa9asQb9+/TBu3Dhs2bKlce8tDA5tpZaBw1ljEiobmRICqA7A+PEQK7USERFRsxPPOpLFxcV+cxitVmvU8y9btgwTJ05EQUFB2Db9+vVDv379vN+PHDkSxcXF+Otf/4rzzz8/6jXixUCSWgxWam0AJTXFiFSnIOuYDuPkyeRfjIiIiCjR4ii2k52d7RdIRnPgwAG8//77WLNmTdzdGjFiBF588cW4j4sFA0lqkVipNYTAD8BQxXWSxFQH2EprwFwkERERNUfxZCTj9dxzz6Fz58649NJL4z52586dyM/Pb9iFo2AgSZmvAev2tDhh7lGqAu6sIw6oxYcZSBIREVHzlITlPwDAMAw899xzmDZtGkwm/9BtwYIFOHjwIJ5//nkAwKJFi9CzZ0+cccYZ3uI8q1evxurVq+O/cAwYSFKLwyI7wcIGkSm4T4oBWEsqoR89lvyLERERESVJMhIX77//PoqKinDDDTcEvVZSUoKioiLv93a7HXfeeScOHjyIrKwsnHHGGXj77bdxySWXJL5jYCBJmY6ZyNiFGtaaAppdgAOHID4L7hIRERE1KyKuLVqbOI0fPx4S5rgVK1b4fX/XXXfhrrvuivsaDcVAkjIei+xEEW4oRuAcyWQt+XHCgF5VnfiTExEREaVIMudINlUMJClzhfiBZpGdAKE+9HwrtSr++xMZTGp2QavDTrT68iCcXPKDiIiImjFFBxQ1eptMwkCSWgzOjQzgCSLjCQ4TGExqdQJbaTX0w0cSc0IiIiKidElSsZ2mjIEkZSZWam2YWALtBAWTWUecUIsOw+l0NP5kRERERGnEoa1EGYrZyGDhspF+90qQlPum6ID1YAX048cTf3IiIiKiVEtSsZ2mjIEkZZ7M+hlNqZABd2AwmYD7q9kFKGalViIiIsoMzEgSZQhWao0iVDYy0j1K4Aef6hBYT+jQq6oSd1IiIiKidOIcSSJq0ZIdcAug2YGsw7WZ9llKRERELRgzkkQZgEt+RBFtyY8k0hwC21EntKIyOJN/OSIiIqLU4BxJomYuVBDJQLJeQ5b8SCBTrcBWUg39yNH0dICIiIgoCZiRJGrOAn6AGUQGCxlEKqm5V4oOtNtZDud3PwBiJPdiRERERCmkGK4tWptMwkCSMkOG/YUn5VIyN1JgFB9iEElERESZxxDXFq1NBmEgSRkjZKVWd4AUbf5f0LGZ9XPuEuo9pWhupOoU2H7SYdTUJP9iRERERKnGqq1EzYQ7ABLV9bVhcgdEKqBbfP6r1O8PDJi8waPhHtduABBAdbr+qzjrhymoPl97hyU0xw8D3z6nKIgEAFMdYCs92SxvGREREVE0CmKYI5mSnqRO2gPJxYsX4y9/+QtKSkpwxhlnYNGiRTjvvPNCtl2zZg2WLFmCXbt2oa6uDmeccQYWLlyICRMmpLjXlBI+waJhdv/XBBgWQEyAbgVEE/c+AUwC0QSKVYeqGdDMrohPUQRQwvzwun/gdUOBiAIxFBhOV9QpdhWKrgA6oNpVKDqgOhRota7AUq1zzftTdUB1pLeITVSRJoCn4FPNVu6AVnyYlVqJiIgoM7Fqa2qtWrUKc+bMweLFizF69Gg888wzmDhxInbv3o3u3bsHtd+yZQsuuugiPPTQQ2jXrh2ee+45XH755fj0008xZMiQNLwDSigFMDRXJtEwA3qW67+GReBsJYDZgJgEms0J1WTAZDKgqgYUBdBUA4oi0FSBogjUgIhOCREs+f4sGz7RlIgCQxT354ErwHTqqne/w6HBcCqQGhPgUKE6FZiqFKgOQLUDplpXNlNzeE6YhHsVp1BBZKqKESkGYCupgl7OSq1ERESUmVi1NcUee+wx3HjjjbjpppsAAIsWLcJ7772HJUuWoLCwMKj9okWL/L5/6KGHsG7dOrz55psMJJsh0dwZRhPgbO0KIHWbQG9rABYD5jZ2aJoBs2agjUmHqriCRE0R79cAoMbwU2mEGcMp7v1aiGjP9xhP0GlA8QaURo4C3VBhGAqcTg12pwrDoQGVJqh2BZYKBardHVyedGUuXSeL4yYlUpquq9kFKCqB6Hr0xkRERETNEedIpo7dbsdnn32G+fPn++0fP348tm7dGtM5DMNAZWUlOnToELZNXV0d6urqvN9XVFQ0rMPUeO5hqo5WruGp9nYC3SYQmwEt2w6zWYfVpMNidkJTBJpqRAwYQwWHEiZg9BwbeIwScE7f4/2u597tG3B6spZAfYCpGwoc2RoMQ4XDboJepwG1GkwVKrRaBdpJwFztnnOZqrgqzJIfqchGqk6BpcKAXlGZ/IsRERERpYkiAiXK0NVorzc3arouXF5eDl3XkZub67c/NzcXpaWlMZ3j0UcfRXV1Na655pqwbQoLC5GTk+PdunXr1qh+U3zEXfSmtgNQnQ9U9BJUDnCg9oyTyOp3Au37/IRTuv6EvA4V6JxTiY6tq9HWUofWljrYTA6YVR0mxYApxMI7qjsz6SswMAx3jO8WeHyoLdy5XENpAU0RmFQDVpOONjY72mbVIie7Bu06VaFtfiW0PlXQ+1Wjpr8dlX0MVHYH6nIAZytXRjZpQV0MwyySRgDNDtjK67jkBxEREWU2I8YtDgsXLoSiKH5bXl5exGM2b96MoUOHwmazoXfv3li6dGn87yVGaS+2owRMXhORoH2hrFy5EgsXLsS6devQuXPnsO0WLFiAefPmeb+vqKhgMJlknqI4jmxXoORsY0A62mHJcqCNxYEssxMmTYfZPdbTN5gLlTH0ZAkD5z165jUGHh8tmPT2U4KPD8X3nKEynqEypYoCb/ArqivLKqJAb1UHe1sTdF1FVUcb1JOuTKW1XIFmdwVeSc9UpnjJD/P+IyyyQ0RERBktWRnJM844A++//773e03Twrbdt28fLrnkEtx888148cUX8e9//xszZ87EKaecgquvvjrua0eTtkCyU6dO0DQtKPtYVlYWlKUMtGrVKtx444149dVX8fOf/zxiW6vVCqvV2uj+Umiiuaqn1nYAHDkGpI2O7E5VsJmdaG+2e4enRgoWffeFahcuMNQgwYFdHAFSuD8KBZ7Tf4irhJ1v6W0fuEPxGRKrAVlmp+scbU5CxBUQO5wanLqKWl2D/agNaq0K22HVNbfS7l6SJF6ebGQilvzwXY8zYOkV33N6l0jRgbY/OmD7+iCch0oacEEiIiKiZsQQ1xatTZxMJlPULKTH0qVL0b17d29dmQEDBmDHjh3461//mlmBpMViwdChQ7Fx40b84he/8O7fuHEjrrzyyrDHrVy5EjfccANWrlyJSy+9NBVdpUDudRudNuBknkDP1tE2twrtLA5YTE60NttdwaM7gvFkDgODRSPMfMRIQWVQVwIyhb5ZSyNKVKmGCQp9A9dYMpCB/YuUEQ08n6IAZsWA2eIKa3VDQZ3FAYdTQ007K1BlgnZSRdZhBVqde5mRGIdFRFyOJPBteQJF1SdY1NwBo+JabkVUn/+q7tc963hq7nOIq6iQVgtYttdCKqti6ywRERFRMxZP1dbAmi2REl979+5FQUEBrFYrzjnnHDz00EPo3bt3yLbbtm3D+PHj/fZNmDABy5Ytg8PhgNlsju3NxCitQ1vnzZuHKVOmYNiwYRg5ciSeffZZFBUVYcaMGQBcw1IPHjyI559/HoAriJw6dSr+9re/YcSIEd5sZlZWFnJyctL2PloMBdDNQG0nwN7egGQ70bFzBWwmJ9pY6mBSjNBBIgQGlKDXwgWKsQaafl0LMfQ0cCis93iEPm+oc8c6TDbaZGPf7KqI4ldECKjvt6IANrMTVpMTFrMTzhxXprKqbRbUkxpMVQps5a6AsiFZSk/gJ2p98OcbFBpmn69NrjdmaICYxHWsCkCR+oyk6lqjE5pAPF8bCrQaFYpDZ6VWIiIiahniWEcycJrdfffdh4ULFwY1P+ecc/D888/jtNNOw+HDh/Hggw9i1KhR+Prrr9GxY8eg9qWlpSHrzzidTpSXlyM/Pz++9xRFWgPJyZMn4+jRo3jggQdQUlKCgQMHYv369ejRowcAoKSkBEVFRd72zzzzDJxOJ2bNmoVZs2Z590+bNg0rVqxIdfdbDgWoywacbYC6DgbadqtAB5urGE4bd/bRl+EOlAKDSQAhA0ogckYyMLgLtx/wn1MZTmCmNNS5Q4k0pDVS9ddwbcNlRF1tXENg4R4GezJPh1NX4XRqqMzOcmUpy+qHvQZlKX0yit6Moer6Q4CY3EuvmOuzi66vBVB9g02fQNcTKXv2Kf5fiwIohuIaTmtXYapQ4cy2wmSxADU1Ye8FERERUSbwTO+J1gYAiouLkZ2d7d0fLhs5ceJE79dnnnkmRo4ciT59+uCf//ynXw0Yv2uEqD8Tan8ipL3YzsyZMzFz5syQrwUGh5s2bUp+h6ieAtjbAo62gL13LbJzanCKrQ45llqY1NDZR1+hAiVPdtLzeqjjg4LQKJnLwNdCZSdjzVhGGg4bbUirr6Aht7FmNhXPD7t/QGwzOwAzoIuKk1116LqKyraugFKrUWA9Xv/h5A0O3RnF+kAxdMDoyTC6LlzfF7+35um/520IAKN+cqRqAIquwHJCganaVTTI3s4Cs9US0/smIiIiatbiyEhmZ2f7BZKxat26Nc4880zs3bs35Ot5eXkh68+YTKaQGczGSnsgSU2PaIAzC3C0Aer61iKnXQ16tqmERdVDBpChAr9wrwH+wWTYNmGydZGGuAb2K9Sw1KjzH2NYKTbS0NjAfkWbL+ntc7hr+Qx5BVyVYNvY6qCLitp8HU5dg+7QUHnMAtWpQK1V6wNKz1qRitRnJzUJXWwnxD7FN2gME2C75mEqUByAVqfA+lP9PE7DpABt2wCHy8LeAyIiIqKMEKk2hW+bRqirq8OePXtw3nnnhXx95MiRePPNN/32bdiwAcOGDUv4/EggjetIUhOkALoNOHGqoHKQHeafHcdpXQ+jR85PaGepRSuTw289R99qrIEBle/3oYKtwIAt1JqO4dZ6jPZ6rNm/WNeLDNX3SAFnpH7F8l4C71249S5tFida2exo1aoO5lNqgU51cLbRodsEhllc8xo1qS+KEy6IBOo//EJtYSjurKSiA+ZqBZbjgFbnXr7EHUiKzQwo/JghIiKizOZZ/iPaFo8777wTmzdvxr59+/Dpp59i0qRJqKiowLRp0wC46slMnTrV237GjBk4cOAA5s2bhz179mD58uVYtmwZ7rzzzoS+Vw9mJAkAYM92b7lO5Pc4imxrLdqY7bCoThhSHwioiniDC98hqoaEL6bj+5rf6yEqrEabFxl4jsDXI2VCYxEtY+l3zTDrWobqV7Rr+Qp33cD3pCgCFYBqchW0sWsaaq0miO4acuodjeozn1GRxq0hWZ+ldGdKDUCtU2A+AZhqARj11zDMgFhMUFQFwpo7RERElMniGNoaqx9//BHXXnstysvLccopp2DEiBH45JNPwtaT6dWrF9avX4+5c+fi6aefRkFBAZ544omkLP0BMJBs8XSbawhr7ekn0bF9FdplnUQH60mofplHwy+Y9O4PmO8YrkhOpKGvoc4V6hxB144w9DVUVdhQogWacQd6UcYrhAuWwxUf8h3uGrS2ZcD9UhRxrdlp1SGGAuWkCjjcrytwV8OpDyYbRdxFdZyu4aymGsB8Et7xuZ6uiqZAzJo7I8lIkoiIiDKYIPxcJd82cXjllVcivh6q2OgFF1yAzz//PL4LNRADyRbInu1awsOZV4cBPUvQxmRHK5MdgH8wY0DxCXbqg8lYgrhQr/meP94Mpa+GLB8SqU/xinEZR69Qa1wGXt9TzdYQxRvAaj4FewzUz6n0BpBSX5THs2kWHYYmMCwa1J/M9dXD3HMZo4o09BXugj4GoDpcy5CYTtYvQ1I/J9P9nkyAnmWCyWKGOB3Rr01ERETUTCmGQIlStlUxGvvX/KaFgWQLYphdVVjtZ1Ujt30lOmVVo53lpF8bv4wepD4oEMWbpTRErR/OGrCURtRCO7Hui7JER6Shr7FmJH3FMww21kDVI/R6kaGXKYlUTMh3DUpDFGiKQMRVPLW+DaBqBhQroFs0qHY1tsnfnnOEayfuhXZ111BW00nAXI0QhX18MpIKYJhVKEmY3E1ERETUpCRhaGtTx0CyhajtBNTkG7DmV2NwXinamOpgVl3DDXVRQg5dDceTnQwMOmMZ5hr4eqR9nvMCwQFltOqt8c6PjGdpj2h9iLw+pPh9HS5b6WoQPDw4MJgMdX5FAUQxoLfSIboCxamEHc4q0QJIz3kNVxvtpGt5D3NNwFpJAUGkt/+aAsXEjxkiIiLKcAbCj+zybZNB+BtehtMtgCMbUM4+gX7tjqODtQZtTXXe1w0o7iGUoedBRhOUwYT/MM1Q7UJ979kHhA7gAucEhutDpH3xiifz2JDKsb7DWCMV9vG08fzreD6DVLj+CCCieM+lqK5AU1EE5tYO2A0FWpVWP1/S78QxzpcUd1EdhwLLCZ/hrIL64DFM98WkABZmJImIiCizxVKVNd6qrU0dA8kMpluAE2foyOlyAoNOKUFrU13EYDFapjFU20iiFZWJVqE13uGugcdE2ucr0XMpG1stFvCZV+k3j9LdH/c1PNfxq+Tqey5NYGrlhO5UoDm0+qxjpCykErDf/bV3OGuVTyYyWhCpKK5lRxhIEhERUabj0FbKBIYJsOcANd10nHn6AXSw1iBLcxXTCQwAg6ulhg4mfedKBmps9i/e4j3RspOB/YwWCCYiexmqD+FEm8PpyS6qPi97sryur93ZSPjPTdVFgaq6ojxPIGqyOGG0VoBKrf5k0T7DPIV1PP/VFZir/OdEApGDSNeBrsqt4BxJIiIiynQMJKm5ExU40d9A2x4ncGa7n9C7TTkMUaH7/cZfXzTHE4ho7kAEiJyZ9J//GH1ZEFe7hg9rbUgwGem4SO0DJTK4jHatUNfzzVJ6gkrPew61NIniV3RHYBgKFABmqxMOiwVqXfi5kl6+rxsKVB1Q6wBLBaA6XEMyRFH8qrP69cFnnUpxz5GESQtuSERERJRJGEhSs6W4KrLW5AsGDt6P3KxKtNZccyFVxXBFmD5tvVVYfYLJaEIFb6GCznBt/Y+LPQsZbyEe3+OAhgeEjR3SGqkgULjrGaKEDmrhXv7DnR0OFXRqPse5lgxxB52qwGitQ9FNUKKtwuH7b+gAtDrAVA1o9gZ88CmAmOBaS5KIiIgok7HYDjVLClDVFZDTqtHnlKM4te0Rb0ASaV3HwGyjb1YyVLtw+3yDSSD24a/xBpNhzx1DQBnuOh6JyD5GCjzjGe4abvkUb9EdQVDGN1r3tSwnDKcC9YQWeT1Jn/UiPUt8mGrq94uqBFdn9RwXuM+dkWQgSURERJmOxXaoeVGAY2fq6NL3CH7e+QcYUGFSdFfwoQAGfIedhhmG6gkCvRGEqxJoLMV0ogm1xiQQvSCOZ3+8wV/gEhqBwVbEvsZYZbWhx0Y7n+85Qp3P9z7p4n8/BQJNc71fTwVXEffwVlEghgLNZEDJsUOvs0Grhc/40/preNeJtLsCSOtx8SvQ4x3S6m7r10tPkR7Ff3irYQL01pYYc95EREREzZRuIGrKUc+slCQDyWZKtwEnOwODB+1Hz9ZHYVZ0ADp0+BbSMQBFDcoyBgoMGj2ZyVjnSsYr1mOjDU2N+nqU4bXxamjA2NDzhRvmqgd875kbqUEgEABqUMbQMFw7xGbA0DWoDgEMn2ynABBXFtJU4xrOqhi+8x19Ak/fYDJSctM9j9KwqGBOkoiIiDIa50hSc1DbEajpa0e/XiUY2q4IdYYJOlRoAX8F8c9IhqjYGkcwmOi1GmMpwBPrdRq6/mRTFxhEBt6vwHmRAKAbKhQIFCXgPbuzlLAaMJwKFEP1q8AKca0Naa50BZKK4TNcVqkfiiGK4peljCWYNEwMJImIiCjTxRBIRi2d37wwkGxmnFlA1rnlOLtjKbrafgIAmFUdqufBFQNGiF/b6wu51Bfe0cMEWNEymKH4BqlB60XGGMwlbRmRWK8fZ9AZOJQWCD9PM9Ix4RhQYlo7E3ANYxVRIIpAJDhrKYZnSRCBWFzBpOZ0BYWKAJod0GoAc7V494X6p/BWbQ0UMEfSe7wCGObmGcgTERERxYwZSWqyFKC2A1DTrw6/7vItckw1MESFQ/yDRlURV0AkClQYMKCGLKIDuIZChqoAGo9YlgNp6Pli2R/z8ZHWnvQJ7ny/DqpQG0MQGE+gGMu5vMGiO1DzfW+aGjDOXhHohn8w75m6KO65kp7oUDSBqAJFV6DogMm9TqS3qI7vW/fJQLou45mbqQRnJQMCSlEU11qSRERERJnMcM8TitomczCQbCYquwPZZ5Xj3FN+RFut1ru/vkiOAd0dxGkwoClwzZcUI2iIaySe5SWA2IruxDw0NuTSIQ1b6zHaMbEGk9GCvkQGhSHP34g5l56++a4nKaJAUdzDWsV/XUlXAwWiu4NJd7ZQ1V1LfJhjWeLD9+VQgWO4yq0mBpJERESU4cRwbdHaZBAGkk2d4lofctAFe3Fm9kG00WqDspC+QgdQBnRoQQFWuKGtHuEymeGvHTlgTdQQ11gC0GhrT6ZCowLFgOq0vvfNNxMZ2E41BOKeP+nH/UcyRVcAXQEMBapDgeq7TmSE7oY6nYL6zKSC0FlJUQHDDEBRM+7Dk4iIiMiLQ1upqTl5CoAhFRjWbj9sitM/Q+iThQQATTEAqK4lIgKykPEGhZH4BoSNndcY9hoJCiaBxKwRGa+GBJHxrDUZ+F/Pa4o7iPQMizUUd8gn7iI5ugLFqUDRFWi1gLkKMNe4LxCUwXTvDjFSw31ad3JTgoNJb0NXRlLRNIiTgSQRERFlKA5tpaZEtwB5Y37EOZ0OwKzorgyiqK6COQAgKoyAYDIUA2rA2oPBw1VDDWENDDxDLQMS7RzJ1Nh5k4nsRyKPjViF1id49JvbqQicUGFSDIiiQFEQIivpykYqBqDWKTBXuoJIxRm+f54gMlQwKWr9HMz6yZgICkgNDe5A0hH+QkRERETNWQvMSKbut36KWV074MTYGoye/hl+3WUHelmPwKY43GtF+jMrOsyq0xtUaO4g03cpENXna0MUb7CnwT+zFYpvYOgJbAzUD7UMXhIkvqxTY9dmjPV4VZGY2nraxbPFc/1Qx/p+b1KNoOGqvucwKYZ3s2g6bCYnbCYnLJoOi6pDUw2YVR1mTYeqClS1vnqO6Aq0OgXm4yqyygBLJaDocGUWVcW1+awXqRiu11XdtTSI67/i+q9nn3u/Yoi3CA8Cgk7RADWnbUz3iYiIiKhZMgQwjChbfL/3FhYW4uyzz0bbtm3RuXNnXHXVVfj2228jHrNp0yYoihK0ffPNN415dyExI9nUKID5vGO4qtvX6GU74v+aJ9sjqnd+o+e/mntZj3BzED1DW0NlJCMPIY0vyxixME8TWM8x2nqVqe5LqK8jze0MCka940/hDUQB173WVAO6Up82VAxAO6nAXAWYTtafW5T6OZCK1C//EfJ2hFn6w7OqTOCyIaK4157MssV0T4iIiIiaJcMAECWhYsSXcNm8eTNmzZqFs88+G06nE/feey/Gjx+P3bt3o3Xr1hGP/fbbb5Gdne39/pRTTonr2rFgINmUKEBNLnBDrx3oZKqE5h626gnAzIoOCCJWYfXOjfP5bT6Wqq2JGJLa0HMke+hprH1I1XnDBZCxnMcve6n4Zp0FUAGTYQCqu/AOxDtn0jM/0lTtCiIVPfQ6keEYnvpOiv9xvt33XCboNQUMJImIiCizJWFo67vvvuv3/XPPPYfOnTvjs88+w/nnnx/x2M6dO6Ndu3ZxXS9eHNraVChAdQEw4KLv0dVyDDbV4R2SqkK8Q1b9g5D6QMIzT9K3CI7v8FbPfEffNR9D8Rwfb2GeaEFkpGxkqoPIVGQeYxn2Gvi6a9ix/xZ4Pt/hrapiwKTqMKm662vFgEmt/6+mugNJABBXkR3TSddw1JBzwUNMqRTVFUSKZ3PvEw3QLQqcVgW6RfFfAiTUOaxmV+VWIiIiokzkCSSjbQAqKir8trq6upguceLECQBAhw4dorYdMmQI8vPzMW7cOHz44YcNf18R8De7JqKqCzBs/G5c1XknzIrTvRZkfQDpS1UM7/xGX77Bml8Q4g4ovXMcRY05e+hpF2+wF0v7wGAp3uMb0i+PZASTscybDPea4RPoRwoi6zd3MAnX8FaTN6g0vJuiiDeYFAFUuwLVgdAVWP121AeLvgGkYgCaAzDVwXUeBRCTqyiUblW860WGPJ/NBEULv2wNERERUbNmSGwbgG7duiEnJ8e7FRYWRj29iGDevHk499xzMXDgwLDt8vPz8eyzz2L16tVYs2YN+vXrh3HjxmHLli0Je6seHNraBIgKjBz/H1zS4Uu0UuvgEM0959HkmgOpAIbP2pEajJiGqwYuAQK4AkM9RMDpmQsZba3HWOY5NjbDmOwg0iOVcyIDrxW0TmRAxjjU8Z5hyyZV937tySwbogIC13BX99BWs6bDqauugjuGe9kPA/7LevhSXGtC+n7vodkFml1cASQA1R1oOm0KxAw4ze7CO3b3AQHVW3WbCSaziZVbiYiIKCOJGJAoa2Z7Xi8uLvabv2i1WqOef/bs2fjyyy/x8ccfR2zXr18/9OvXz/v9yJEjUVxcjL/+9a9Rh8PGi4FkurmHtF7Z8XO0VuugQ/Uu1+HKPKohCzypigFdQmd4QgUjmjdqMLzFevzmUUYqkuNTrTXc+SMJFXw2NhPZnIQKWP3vffjMrG9RHU8Q6RnWqvmcV1N06OKupmsA0FwFmEyaAZOmu6vnuAvi6GEK6QD1cyB9uyPuY5yoH9svCrQ619cOTYFuAwwTYNZdmUt3E+9/dasKs9UKnPSp8kNERESUKaQ+4xixDYDs7Gy/QDKa2267DW+88Qa2bNmCrl27xt21ESNG4MUXX4z7uGgYSKZZTR5w9kW70VatBQCo0F0VOKFCg7gCPMU9B9JdfEdH/dBUDeLNMEbLJgZSFYHTqA9Gw82PjFQMJ965kYkIIJtKsBnvHMhAgcNZA48NNZzVFVC6gkg1aIyqKyPpykCaYNF06KLCpGlQTALDLDAs7qxkmGAyKIj0ObdrU1zzHk2ur1UnoNUB9rau1w2HKysZtISlWYVitUS8H0RERETNloRYdDtkm3hOKbjtttuwdu1abNq0Cb169WpQ13bu3In8/PwGHRsJA8l0UoA+Y/fhio67YFMc3uGqhniGs7q+D1f4RoN4y+l4Krx6eIa1es4ZeA7fANTVzj+IDLU8SGCgGm3+ZLg5m5H2hTxPCgLHwDmJybxGLMNYPf+tz0Ya3u81n31+x7n/x/NvaNOcMESBQ9OgmnUYNg2ONhoMzTXXUXG614r0/dwL+Az0vCYKYJgVKHr9i6K6507qAtEU11DXLMBkuAv6+L53kwLYog/bICIiImqWDKN+WFY4UYa+Bpo1axZefvllrFu3Dm3btkVpaSkAICcnB1lZWQCABQsW4ODBg3j++ecBAIsWLULPnj1xxhlnwG6348UXX8Tq1auxevXq+N9TFAwk08jRCrg673MUmH6CRdGhoz4Y1MU1vLU+IFT8AsV413cE/AMXPUTaKTCIDBzSGtM1YhzGGi04bCpZx8YIN3zVQJgMpE8EFzikNfB1T5tAqriGsELVYVF1OFUVVs0Jk1WHvbUOe7YKk1mB0wlodngDR09AqTpcQaDiBDSHeF9TneJqo7uu4ynA4/3nVgDD7NrvDVB9PisNDYDZ3LAbSURERNTUJSEjuWTJEgDAmDFj/PY/99xzmD59OgCgpKQERUVF3tfsdjvuvPNOHDx4EFlZWTjjjDPw9ttv45JLLonr2rFgIJkOClDXDhj2i//gvFY/wBCgVjTY4RpmqoqrWqsBFZoYgKJCFYH7d3hoMOBwF+JxzXVU4RDNteakO/g0RPHOtXQFMfVFdozAgDGG4jqh5lP6r1UZnH0MDJaaQgYy2nxF3+/DZSZjHbLqdP97AIDTCA76TWr9ki6eNSFd8yBdQ1hdy3qIX2GdwGxkYPVeTXVChwKnocGm1Re2yWl1EtWajhoB9DoVil2Fald8Akl3xVVPAGgAqtM1dFVxAKaTrq89mUZRXdlHw+QKIJ2tAN1av3CkqruGvHrviUWB2BhIEhERUWYSXYd4/uIero1Efj24ffTAc8WKFX7f33XXXbjrrrviuk5DMZBMA1GA6p46zsn5AaeoGmpFh24Y0EWBa/VIBQ4xwS4a7KLBISY4RHMFi3AFjZ7A0GGYXO0NLWIAGRgIBg5l9cyVDAwcww1fjTb3MVQwGdc8yDiyoOEEzSGMsw/hRJvn6ck4eoLHSPMfffsZy9qTkXiy1BoEUHWYRIUFCpyiuovuGFA0wzXM1BCIXh9IQnVlH11TcRXvH9REBRT3NFo1oJCOp8COmADDIhCTQBG4hrkGxM2iwFXqlYiIiCgTGRKhmqFbnBnJpo6BZDooQKvcKpzQs7DXqQHQcFzPQq2YUStmnNBbwe4OHmsMCxyGCXWGyRsougJJTyZS9QaGunuf4a7e6WkDBAc5nn2eYMfp064hwZ+EyE76HisNCOBCZQYVn/mDHoEBoxLwQxxLcBbtfQYu1+EhPvfV81qo9+rbbxGB4S5g5DeEVfyL68Cov65J1aGKAsOTmYQCE3QYUOr75hkW7e6LpvgPjVXclVsV73BUqS+t6tdZ935vgR0AqntYv89QVu8ak6r70p64VA1xWkWBKPE/A0RERETNggiAaHMkMyuQTHuKYPHixejVqxdsNhuGDh2Kjz76KGL7zZs3Y+jQobDZbOjduzeWLl2aop4mjqIDtndz8OoTFzXpuYBNoW+pXOsxnQLnR/q9FjBxW4vxnujiykY6DRUOXYNT1yC64spEGnBFe54/nvmuLSmKa7UQw3+oq+9lDZ85kqICorleVHQF2klXMR8vAazHnFAPlcXUbyIiIqLmRgyJacskac1Irlq1CnPmzMHixYsxevRoPPPMM5g4cSJ2796N7t27B7Xft28fLrnkEtx888148cUX8e9//xszZ87EKaecgquvvjoN76Dx+pmdqBUdrRUHakVDrZiQrZ5ErVhgFw21hhkOMaFWzN75j75DW+sMc9DcSId7fclYh7YaPnP5YhnaGq2gTixfR5Osoa1xHR9nEOs7NxJA2KGtnrmRnrmQvvt850v6zpE0qa4x9ZoiMLnH36uKeOdIeoJNz9BWz9xZwx1Iinj+CKZEHXURxCfY9BzrF2AaruGwokr9/lDXyLAPTyIiIiIvcf9iFLVN5khrIPnYY4/hxhtvxE033QTAVa72vffew5IlS1BYWBjUfunSpejevTsWLVoEABgwYAB27NiBv/71r802kDxi6DAEcLjXhwRcAZBZcVU1MRTVlYEyAAc0b+BhQKC6Fpz0FttRRbxDHQ1RYCgqdBHXkEioPpVaXYuCGKK62rqjA7+qrYo7EHW/Fq3YTqi4zzN0M16GKI0OAsOJpdhOQ87hu18V97BVTbxDhQMDau+/kc+/l99w3biD2PphrZ4/IDh9Nofhzkg6FShOBYquQHGiPhPpnhfpm4VU3UuEqA7AdNI9R1J3JzI119eGGVDNgG4FYPKcG0GBpGIIRI9vgjkRERFRcyGGQKL8/hZL8ZzmJG2BpN1ux2effYb58+f77R8/fjy2bt0a8pht27Zh/PjxfvsmTJiAZcuWweFwwBxieYG6ujrU1dWXjzxx4gQAQLfXNvYtJMT//WRDjlYLhwicIrALUCMGHKLCIQZqRWBAh0MUGHAtMO8QDQZUdwBpQAfgMEwQ1BfOCVV4R7zZyfrr+1ZuDVr+w30+CTHnsaGVW0OJ1KYxw2tDBWPhQplQBXHCtY01yBNRXFMMRfFOGqzP/PqfS/eZQ2lSDW9m0rN+pB311VrN3v3BVVsdPnNnT+oaapwqahxO1FQ4YT+pQI6rUBwKFIcCOBAyy6gY7uU/3ENUVYdAq/VZ/sMTSFoV6FYFugWABsAkkDoFqASMOp9pAALodjuczpPQxQEiIiIiX064fj9ozoGWU+qiZhw97zNTpC2QLC8vh67ryM3N9dufm5vrXWwzUGlpacj2TqcT5eXlyM/PDzqmsLAQ999/f9D+r196oBG9T5yJz6W7B0RERERE6Xf06FHk5OSkuxtxsVgsyMvLw8el62Nqn5eXB4vFkuRepUbaq7YqAZUcRSRoX7T2ofZ7LFiwAPPmzfN+f/z4cfTo0QNFRUXN7kFtrioqKtCtWzcUFxcjOzs73d1pEXjPU4/3PPV4z1OP9zz1eM9Tj/c89U6cOIHu3bujQ4cO6e5K3Gw2G/bt2we73R5Te4vFApvNluRepUbaAslOnTpB07Sg7GNZWVlQ1tEjLy8vZHuTyYSOHTuGPMZqtcJqtQbtz8nJ4YdDimVnZ/Oepxjveerxnqce73nq8Z6nHu956vGep57aTNecttlsGRMcxiNt/1oWiwVDhw7Fxo0b/fZv3LgRo0aNCnnMyJEjg9pv2LABw4YNCzk/koiIiIiIiBIvrWH/vHnz8I9//APLly/Hnj17MHfuXBQVFWHGjBkAXMNSp06d6m0/Y8YMHDhwAPPmzcOePXuwfPlyLFu2DHfeeWe63gIREREREVGLk9Y5kpMnT8bRo0fxwAMPoKSkBAMHDsT69evRo0cPAEBJSQmKioq87Xv16oX169dj7ty5ePrpp1FQUIAnnngirqU/rFYr7rvvvpDDXSk5eM9Tj/c89XjPU4/3PPV4z1OP9zz1eM9Tj/e8eVKkOdfZJSIiIiIiopRrnjNaiYiIiIiIKG0YSBIREREREVFcGEgSERERERFRXBhIEhERERERUVwyMpBcvHgxevXqBZvNhqFDh+Kjjz6K2H7z5s0YOnQobDYbevfujaVLl6aop5kjnnu+adMmKIoStH3zzTcp7HHztmXLFlx++eUoKCiAoih4/fXXox7D57xx4r3nfM4bp7CwEGeffTbatm2Lzp0746qrrsK3334b9Tg+5w3XkHvO57xxlixZgkGDBnkXvh85ciTeeeediMfwGW+ceO85n/HEKywshKIomDNnTsR2fNabvowLJFetWoU5c+bg3nvvxc6dO3Heeedh4sSJfsuI+Nq3bx8uueQSnHfeedi5cyfuuece3H777Vi9enWKe958xXvPPb799luUlJR4t759+6aox81fdXU1Bg8ejKeeeiqm9nzOGy/ee+7B57xhNm/ejFmzZuGTTz7Bxo0b4XQ6MX78eFRXV4c9hs954zTknnvwOW+Yrl274s9//jN27NiBHTt24MILL8SVV16Jr7/+OmR7PuONF+899+Aznhjbt2/Hs88+i0GDBkVsx2e9mZAMM3z4cJkxY4bfvv79+8v8+fNDtr/rrrukf//+fvtuueUWGTFiRNL6mGnivecffvihAJCffvopBb3LfABk7dq1EdvwOU+sWO45n/PEKisrEwCyefPmsG34nCdWLPecz3nitW/fXv7xj3+EfI3PeHJEuud8xhOnsrJS+vbtKxs3bpQLLrhA7rjjjrBt+aw3DxmVkbTb7fjss88wfvx4v/3jx4/H1q1bQx6zbdu2oPYTJkzAjh074HA4ktbXTNGQe+4xZMgQ5OfnY9y4cfjwww+T2c0Wj895+vA5T4wTJ04AADp06BC2DZ/zxIrlnnvwOW88XdfxyiuvoLq6GiNHjgzZhs94YsVyzz34jDferFmzcOmll+LnP/951LZ81puHjAoky8vLoes6cnNz/fbn5uaitLQ05DGlpaUh2zudTpSXlyetr5miIfc8Pz8fzz77LFavXo01a9agX79+GDduHLZs2ZKKLrdIfM5Tj8954ogI5s2bh3PPPRcDBw4M247PeeLEes/5nDfeV199hTZt2sBqtWLGjBlYu3YtTj/99JBt+YwnRjz3nM94Yrzyyiv4/PPPUVhYGFN7PuvNgyndHUgGRVH8vheRoH3R2ofaT+HFc8/79euHfv36eb8fOXIkiouL8de//hXnn39+UvvZkvE5Ty0+54kze/ZsfPnll/j444+jtuVznhix3nM+543Xr18/7Nq1C8ePH8fq1asxbdo0bN68OWxgw2e88eK553zGG6+4uBh33HEHNmzYAJvNFvNxfNabvozKSHbq1AmapgVlwsrKyoL+quGRl5cXsr3JZELHjh2T1tdM0ZB7HsqIESOwd+/eRHeP3PicNw18zuN322234Y033sCHH36Irl27RmzL5zwx4rnnofA5j4/FYsGpp56KYcOGobCwEIMHD8bf/va3kG35jCdGPPc8FD7j8fnss89QVlaGoUOHwmQywWQyYfPmzXjiiSdgMpmg63rQMXzWm4eMCiQtFguGDh2KjRs3+u3fuHEjRo0aFfKYkSNHBrXfsGEDhg0bBrPZnLS+ZoqG3PNQdu7cifz8/ER3j9z4nDcNfM5jJyKYPXs21qxZgw8++AC9evWKegyf88ZpyD0Phc9544gI6urqQr7GZzw5It3zUPiMx2fcuHH46quvsGvXLu82bNgw/OY3v8GuXbugaVrQMXzWm4m0lPhJoldeeUXMZrMsW7ZMdu/eLXPmzJHWrVvL/v37RURk/vz5MmXKFG/7H374QVq1aiVz586V3bt3y7Jly8RsNstrr72WrrfQ7MR7zx9//HFZu3atfPfdd/Kf//xH5s+fLwBk9erV6XoLzU5lZaXs3LlTdu7cKQDksccek507d8qBAwdEhM95MsR7z/mcN86tt94qOTk5smnTJikpKfFuNTU13jZ8zhOrIfecz3njLFiwQLZs2SL79u2TL7/8Uu655x5RVVU2bNggInzGkyHee85nPDkCq7byWW+eMi6QFBF5+umnpUePHmKxWORnP/uZX+nyadOmyQUXXODXftOmTTJkyBCxWCzSs2dPWbJkSYp73PzFc88ffvhh6dOnj9hsNmnfvr2ce+658vbbb6eh182Xpxx54DZt2jQR4XOeDPHecz7njRPqXgOQ5557ztuGz3liNeSe8zlvnBtuuMH7/52nnHKKjBs3zhvQiPAZT4Z47zmf8eQIDCT5rDdPioh75ioRERERERFRDDJqjiQRERERERElHwNJIiIiIiIiigsDSSIiIiIiIooLA0kiIiIiIiKKCwNJIiIiIiIiigsDSSIiIiIiIooLA0kiIiIiIiKKCwNJIiIiIiIiigsDSSIiIiIiIooLA0kiIiIiIiKKCwNJIiIiIiIiigsDSSIialaOHDmCvLw8PPTQQ959n376KSwWCzZs2JDGnhEREbUciohIujtBREQUj/Xr1+Oqq67C1q1b0b9/fwwZMgSXXnopFi1alO6uERERtQgMJImIqFmaNWsW3n//fZx99tn44osvsH37dthstnR3i4iIqEVgIElERM3SyZMnMXDgQBQXF2PHjh0YNGhQurtERETUYnCOJBERNUs//PADDh06BMMwcODAgXR3h4iIqEVhRpKIiJodu92O4cOH46yzzkL//v3x2GOP4auvvkJubm66u0ZERNQiMJAkIqJm5w9/+ANee+01fPHFF2jTpg3Gjh2Ltm3b4q233kp314iIiFoEDm0lIqJmZdOmTVi0aBFeeOEFZGdnQ1VVvPDCC/j444+xZMmSdHePiIioRWBGkoiIiIiIiOLCjCQRERERERHFhYEkERERERERxYWBJBEREREREcWFgSQRERERERHFhYEkERERERERxYWBJBEREREREcWFgSQRERERERHFhYEkERERERERxYWBJBEREREREcWFgSQRERERERHFhYEkERERERERxeX/A1LutvL0PQ2YAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "from pyro.compressible import simulation as s\n", + "from matplotlib import pyplot as plt\n", + "\n", + "# Get outputs for plotting. Change into primitive variables for ease of use\n", + "conserved_data = pyro_sim.sim.cc_data\n", + "q = s.cons_to_prim(conserved_data.data, conserved_data.get_aux(\"gamma\"), s.Variables(conserved_data), conserved_data.grid)\n", + "\n", + "# Isolate variable to plot. If choosing another, replace name in pcolormesh and title\n", + "# p = q[:, :, s.Variables(conserved_data).ip] # Pressure\n", + "# u = q[:, :, s.Variables(conserved_data).iu] # x Velocity\n", + "# v = q[:, :, s.Variables(conserved_data).iv] # y Velocity\n", + "rho = q[:, :, s.Variables(conserved_data).irho] # Density\n", + "\n", + "# Set up for plotting - necessary for plt.pcolormesh, relates to labeling axes correctly\n", + "myg = conserved_data.grid\n", + "x = myg.scratch_array()\n", + "y = myg.scratch_array()\n", + "x.v()[:, :] = myg.x2d.v()[:, :]\n", + "y.v()[:, :] = myg.y2d.v()[:, :]\n", + "\n", + "# Plot\n", + "plt.figure(figsize=(12, 3))\n", + "plt.pcolormesh(x.v(), y.v(), rho[0:len(x.v()),0:len(x.v()[0])], shading=\"nearest\", cmap=pyro_sim.sim.cm)\n", + "plt.xlabel(\"x\")\n", + "plt.ylabel(\"y\")\n", + "plt.title(\"Solution after 500 timesteps for density\")\n", + "plt.colorbar()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "ae86e8cc", + "metadata": {}, + "source": [ + "## Results\n", + "\n", + "While 500 time steps is not enough to evolve the solution completely, it does resemble the results shown in Woodward and Colella [1]. The same sections within the shock can be identified, including the turbulent jet that forms near the bottom wall. This is an important feature, since it requires more resolution to accurately model. \n", + "\n", + "Pyro does not use one of the flux limiting schemes explicitly tested in Woodward and Colella [1]; instead, it uses a hybrid approach proposed by Colella in \"Multidimensional upwind methods for hyperbolic conservation laws\" [2]. However, the results are comparable, indicating that Pyro is adept at handling discontinuities.\n", + "\n", + "For more information about the Double Mach Reflection problem, as well as an explanation of three different types of flux limiters and their performance on the problem, see \"The Numerical Simulation of Two-Dimensional Fluid Flow with Strong Shocks\" [1]. " + ] + }, + { + "cell_type": "markdown", + "id": "dcca0f23", + "metadata": {}, + "source": [ + "## References\n", + "\n", + "[1] P.R. Woodward & P. Colella. The Numerical Simulation of Two-Dimensional Fluid Flow with Strong Shocks, *Journal of Computational Physics*, 54:115-173, April 1984. [10.1016/0021-9991(84)90142-6](10.1016/0021-9991(84)90142-6).\n", + "\n", + "[2] P. Colella. Multidimensional upwind methods for hyperbolic conservation laws. *Journal of Computational Physics*, 87:171–200, March 1990. [doi:10.1016/0021-9991(90)90233-Q](doi:10.1016/0021-9991(90)90233-Q)." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.9" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}