FiniteVolumeGPU/MPITest.ipynb
André R. Brodtkorb c51afef9fc MPI prototype
2018-11-21 07:49:39 +01:00

762 lines
39 KiB
Plaintext

{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Starting cluster\n",
"\n",
"## Prerequisites\n",
"First, you need to install MPI, on windows use MS-MPI:\n",
"https://msdn.microsoft.com/en-us/library/bb524831(v=vs.85).aspx\n",
"\n",
"\n",
"## With a profile (not working)\n",
"In theory, you should be able to create a profile using\n",
"```\n",
"ipython profile create --parallel --profile=myprofile\n",
"```\n",
"and then set\n",
"```\n",
"c.IPClusterEngines.engine_launcher_class = 'MPIEngineSetLauncher'\n",
"```\n",
"in ```<IPYTHON-DIR>/profile_myprofile/ipcluster_config.py```. This should then enable you to start a cluster using\n",
"```\n",
"ipcluster start --profile=myprofile\n",
"```\n",
"or alternatively through the Clusters tab in Jupyter\n",
"\n",
"\n",
"## Without a profile (not working)\n",
"An alternative is to run\n",
"```\n",
"ipcluster start --engines=MPI\n",
"```\n",
"\n",
"\n",
"## Manual start (working)\n",
"This, however, does *not* work for me on Windows. What does work is the following:\n",
"\n",
"Start a controller using\n",
"```\n",
"ipcontroller --ip='*'\n",
"```\n",
"and then start several engines using mpiexec:\n",
"```\n",
"mpiexec -n 4 ipengine --mpi\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"from GPUSimulators import IPythonMagic"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Console logger using level INFO\n",
"File logger using level DEBUG to mpi.log\n",
"Python version 3.6.7 |Anaconda custom (64-bit)| (default, Oct 28 2018, 19:44:12) [MSC v.1915 64 bit (AMD64)]\n"
]
}
],
"source": [
"%setup_logging --out mpi.log"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Starting IPController\n",
"Starting IPEngines\n",
"Waiting for cluster...\n",
"Waiting for cluster...\n",
"Waiting for cluster...\n",
"Waiting for cluster...\n",
"Waiting for cluster...\n",
"Waiting for cluster...\n",
"Waiting for cluster...\n",
"Waiting for cluster...\n",
"Waiting for cluster...\n",
"Waiting for cluster...\n",
"Waiting for cluster...\n",
"Waiting for cluster...\n",
"Done\n"
]
}
],
"source": [
"%setup_mpi --num_engines 4"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"profile: default\n",
"Number of ids: 4\n",
"IDs: [0, 1, 2, 3]\n"
]
}
],
"source": [
"import ipyparallel\n",
"\n",
"# attach to a running cluster\n",
"cluster = ipyparallel.Client()#profile='mpi')\n",
"\n",
"print('profile:', cluster.profile)\n",
"print('Number of ids:', len(cluster.ids))\n",
"print(\"IDs:\", cluster.ids) # Print process id numbers"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[stdout:0] MPI rank 0 of 4 OK\n",
"[stdout:1] MPI rank 1 of 4 OK\n",
"[stdout:2] MPI rank 2 of 4 OK\n",
"[stdout:3] MPI rank 3 of 4 OK\n"
]
}
],
"source": [
"%%px\n",
"\n",
"from mpi4py import MPI\n",
"comm = MPI.COMM_WORLD\n",
"print(\"MPI rank {:d} of {:d} OK\".format(comm.rank, comm.size))\n",
"comm.Barrier() "
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[stdout:0] \n",
"Starting\n",
"0: sent data to 1: [91 81 57 75 60 73 2 68 11 50 71 17 65 20 4 64 48 34 66 37 5 85 42 80\n",
" 25 78 83 94 38 72 27 93 82 56 92 51 88 22 29 79 14 13 47 35 0 30 63 3\n",
" 12 99 24 23 15 32 44 1 62 67 19 8 87 45 16 61 84 7 49 10 21 69 89 70\n",
" 39 86 98 90 36 9 53 54 33 43 26 58 77 40 95 6 41 74 46 28 59 76 31 97\n",
" 55 18 96 52]\n",
"\n",
"[stdout:1] \n",
"Starting\n",
"1: received data from 0: [91 81 57 75 60 73 2 68 11 50 71 17 65 20 4 64 48 34 66 37 5 85 42 80\n",
" 25 78 83 94 38 72 27 93 82 56 92 51 88 22 29 79 14 13 47 35 0 30 63 3\n",
" 12 99 24 23 15 32 44 1 62 67 19 8 87 45 16 61 84 7 49 10 21 69 89 70\n",
" 39 86 98 90 36 9 53 54 33 43 26 58 77 40 95 6 41 74 46 28 59 76 31 97\n",
" 55 18 96 52]\n",
"\n",
"[stdout:2] \n",
"Starting\n",
"2: idle\n",
"\n",
"[stdout:3] \n",
"Starting\n",
"3: idle\n",
"\n"
]
}
],
"source": [
"%%px\n",
"\n",
"from mpi4py import MPI\n",
"import numpy\n",
"\n",
"comm = MPI.COMM_WORLD\n",
"rank = comm.Get_rank()\n",
"\n",
"print(\"Starting\")\n",
"# passing MPI datatypes explicitly\n",
"if rank == 0:\n",
" data = numpy.arange(100, dtype='i')\n",
" numpy.random.shuffle(data)\n",
" comm.Send([data, MPI.INT], dest=1, tag=77)\n",
" print(\"{0}: sent data to 1: {1}\".format(rank, data))\n",
"elif rank == 1:\n",
" data = numpy.empty(100, dtype='i')\n",
" comm.Recv([data, MPI.INT], source=0, tag=77)\n",
" print(\"{0}: received data from 0: {1}\".format(rank, data))\n",
"else:\n",
" print(\"{0}: idle\".format(rank))\n",
" \n",
"print()"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"%%px\n",
"from GPUSimulators import IPythonMagic"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"[stderr:0] \n",
"Console logger using level INFO\n",
"File logger using level DEBUG to mpi_0.log\n",
"Python version 3.6.7 |Anaconda custom (64-bit)| (default, Oct 28 2018, 19:44:12) [MSC v.1915 64 bit (AMD64)]\n",
"[stderr:1] \n",
"Console logger using level INFO\n",
"File logger using level DEBUG to mpi_1.log\n",
"Python version 3.6.7 |Anaconda custom (64-bit)| (default, Oct 28 2018, 19:44:12) [MSC v.1915 64 bit (AMD64)]\n",
"[stderr:2] \n",
"Console logger using level INFO\n",
"File logger using level DEBUG to mpi_2.log\n",
"Python version 3.6.7 |Anaconda custom (64-bit)| (default, Oct 28 2018, 19:44:12) [MSC v.1915 64 bit (AMD64)]\n",
"[stderr:3] \n",
"Console logger using level INFO\n",
"File logger using level DEBUG to mpi_3.log\n",
"Python version 3.6.7 |Anaconda custom (64-bit)| (default, Oct 28 2018, 19:44:12) [MSC v.1915 64 bit (AMD64)]\n"
]
}
],
"source": [
"%%px\n",
"%setup_logging --out \"'mpi_' + str(MPI.COMM_WORLD.rank) + '.log'\""
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"[stderr:0] \n",
"Registering my_context in user workspace\n",
"PyCUDA version 2018.1.1\n",
"CUDA version (10, 0, 0)\n",
"Driver version 10000\n",
"Using 'GeForce 840M' GPU\n",
"Created context handle <863101499808>\n",
"Using CUDA cache dir c:\\Users\\anbro\\Documents\\projects\\GPUSimulators\\GPUSimulators\\cuda_cache\n",
"Autotuning enabled. It may take several minutes to run the code the first time: have patience\n",
"[stderr:1] \n",
"Registering my_context in user workspace\n",
"PyCUDA version 2018.1.1\n",
"CUDA version (10, 0, 0)\n",
"Driver version 10000\n",
"Using 'GeForce 840M' GPU\n",
"Created context handle <912376276640>\n",
"Using CUDA cache dir c:\\Users\\anbro\\Documents\\projects\\GPUSimulators\\GPUSimulators\\cuda_cache\n",
"Autotuning enabled. It may take several minutes to run the code the first time: have patience\n",
"[stderr:2] \n",
"Registering my_context in user workspace\n",
"PyCUDA version 2018.1.1\n",
"CUDA version (10, 0, 0)\n",
"Driver version 10000\n",
"Using 'GeForce 840M' GPU\n",
"Created context handle <390889602416>\n",
"Using CUDA cache dir c:\\Users\\anbro\\Documents\\projects\\GPUSimulators\\GPUSimulators\\cuda_cache\n",
"Autotuning enabled. It may take several minutes to run the code the first time: have patience\n",
"[stderr:3] \n",
"Registering my_context in user workspace\n",
"PyCUDA version 2018.1.1\n",
"CUDA version (10, 0, 0)\n",
"Driver version 10000\n",
"Using 'GeForce 840M' GPU\n",
"Created context handle <550832168304>\n",
"Using CUDA cache dir c:\\Users\\anbro\\Documents\\projects\\GPUSimulators\\GPUSimulators\\cuda_cache\n",
"Autotuning enabled. It may take several minutes to run the code the first time: have patience\n"
]
}
],
"source": [
"%%px\n",
"%cuda_context_handler my_context"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"%%px\n",
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"%%px\n",
"import numpy as np\n",
"from matplotlib import pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"%%px\n",
"\n",
"\n",
"def getFactors(number, num_factors):\n",
" # Adapted from https://stackoverflow.com/questions/28057307/factoring-a-number-into-roughly-equal-factors\n",
" # Original code by https://stackoverflow.com/users/3928385/ishamael\n",
" \n",
" #Dictionary to remember already computed permutations\n",
" memo = {}\n",
" def dp(n, left): # returns tuple (cost, [factors])\n",
" \"\"\"\n",
" Recursively searches through all factorizations\n",
" \"\"\"\n",
" \n",
" #Already tried: return existing result\n",
" if (n, left) in memo: \n",
" return memo[(n, left)]\n",
"\n",
" #Spent all factors: return number itself\n",
" if left == 1:\n",
" return (n, [n])\n",
"\n",
" #Find new factor\n",
" i = 2\n",
" best = n\n",
" bestTuple = [n]\n",
" while i * i < n:\n",
" #If factor found\n",
" if n % i == 0:\n",
" #Factorize remainder\n",
" rem = dp(n // i, left - 1)\n",
" \n",
" #If new permutation better, save it\n",
" if rem[0] + i < best:\n",
" best = rem[0] + i\n",
" bestTuple = [i] + rem[1]\n",
" i += 1\n",
"\n",
" #Store calculation\n",
" memo[(n, left)] = (best, bestTuple)\n",
" return memo[(n, left)]\n",
" \n",
" assert(isinstance(number, int))\n",
" assert(isinstance(num_factors, int))\n",
" \n",
" factors = dp(number, num_factors)[1]\n",
" \n",
" if (len(factors) < num_factors):\n",
" #Split problematic 4\n",
" if (4 in factors):\n",
" factors.remove(4)\n",
" factors.append(2)\n",
" factors.append(2)\n",
" \n",
" factors = factors + [1]*(num_factors - len(factors))\n",
" return factors"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"%%px\n",
"\n",
"def getCoordinate(rank, grid):\n",
" i = (rank % grid[0])\n",
" j = (rank // grid[0])\n",
" \n",
" return i, j\n",
"\n",
"\n",
"def getRank(i, j, grid):\n",
" return j*grid[0] + i\n",
"\n",
"\n",
"def getEast(rank, grid):\n",
" i, j = getCoordinate(rank, grid)\n",
" i = (i+1) % grid[0]\n",
" return getRank(i, j, grid)\n",
"\n",
"def getWest(rank, grid):\n",
" i, j = getCoordinate(rank, grid)\n",
" i = (i+grid[0]-1) % grid[0]\n",
" return getRank(i, j, grid)\n",
"\n",
"def getNorth(rank, grid):\n",
" i, j = getCoordinate(rank, grid)\n",
" j = (j+1) % grid[1]\n",
" return getRank(i, j, grid)\n",
"\n",
"def getSouth(rank, grid):\n",
" i, j = getCoordinate(rank, grid)\n",
" j = (j+grid[1]-1) % grid[1]\n",
" return getRank(i, j, grid)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[stdout:0] \n",
"[3, 4]\n",
"0 1 2 \n",
"3 4 5 \n",
"6 7 8 \n",
"9 10 11 \n",
"0 (0, 0) 3 9 1 2\n",
"1 (1, 0) 4 10 2 0\n",
"2 (2, 0) 5 11 0 1\n",
"3 (0, 1) 6 0 4 5\n",
"4 (1, 1) 7 1 5 3\n",
"5 (2, 1) 8 2 3 4\n",
"6 (0, 2) 9 3 7 8\n",
"7 (1, 2) 10 4 8 6\n",
"8 (2, 2) 11 5 6 7\n",
"9 (0, 3) 0 6 10 11\n",
"10 (1, 3) 1 7 11 9\n",
"11 (2, 3) 2 8 9 10\n"
]
}
],
"source": [
"%%px \n",
"\n",
"if (MPI.COMM_WORLD.rank == 0):\n",
" n_procs = 12\n",
" grid = getFactors(n_procs, 2)\n",
" \n",
" print(grid)\n",
" for j in range(grid[1]):\n",
" for i in range(grid[0]):\n",
" print(j*grid[0]+i, end=\" \")\n",
" print()\n",
"\n",
" for i in range(n_procs):\n",
" print(getRank(*getCoordinate(i, grid), grid), getCoordinate(i, grid), getNorth(i, grid), getSouth(i, grid), getEast(i, grid), getWest(i, grid))"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"%%px\n",
"\n",
"def plotSolution(data, grid):\n",
" ny, nx = data.shape\n",
" \n",
" if (MPI.COMM_WORLD.rank != 0):\n",
" mpi_request = MPI.COMM_WORLD.Isend(data, dest=0, tag=MPI.COMM_WORLD.rank)\n",
" mpi_request.wait()\n",
" else:\n",
" def my_imshow(data, idx):\n",
" i, j = getCoordinate(idx, grid)\n",
" x = i * width\n",
" y = j * height \n",
" extent=[x, x+width, y, y+height]\n",
" plt.imshow(data, extent=extent, vmin=0.4, vmax=0.6)\n",
" \n",
" mpi_requests = []\n",
" for k in range(1, MPI.COMM_WORLD.size):\n",
" buffer = np.empty((ny, nx), dtype=np.float32)\n",
" mpi_requests += [(buffer, MPI.COMM_WORLD.Irecv(buffer, source=k, tag=k))]\n",
"\n",
" plt.figure()\n",
" my_imshow(data, 0)\n",
" idx = 1\n",
" for buffer, request in mpi_requests:\n",
" request.wait()\n",
" my_imshow(buffer, idx)\n",
" idx += 1\n",
" plt.axis('tight')\n",
" plt.colorbar()"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"%%px\n",
"\n",
"from GPUSimulators.helpers import InitialConditions\n",
"from GPUSimulators.Simulator import BoundaryCondition\n",
"\n",
"nx = 16\n",
"ny = 15\n",
"g = 9.81\n",
"dt = 0.05\n",
"width = 100\n",
"height = 100\n",
"\n",
"\n",
"if (MPI.COMM_WORLD.rank == 0):\n",
" h0, hu0, hv0, dx, dy = InitialConditions.bump(nx, ny, width, height, g, x_center=0.75, y_center=0.75)\n",
"else:\n",
" h0, hu0, hv0, dx, dy = InitialConditions.bump(nx, ny, width, height, g, h_ref=0.5, h_amp=0.0, u_amp=0.0, v_amp=0.0)\n",
" \n",
"bc = BoundaryCondition({\n",
" 'north': BoundaryCondition.Type.Dirichlet,\n",
" 'south': BoundaryCondition.Type.Dirichlet,\n",
" 'east': BoundaryCondition.Type.Dirichlet,\n",
" 'west': BoundaryCondition.Type.Dirichlet\n",
"})\n",
"\n",
"arguments = {\n",
" 'context': my_context,\n",
" 'h0': h0, 'hu0': hu0, 'hv0': hv0,\n",
" 'nx': nx, 'ny': ny,\n",
" 'dx': dx, 'dy': dy, \n",
" 'g': g,\n",
" 'boundary_conditions': bc\n",
"} \n"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[output:0]"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAACzCAYAAACZ+efrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAACYZJREFUeJzt3M+LVfcdxvHnqcZIJATFKpKEElRCJaWTMthFoBhCgslGEyjElYvABKl/gLtk6SaELoJgWhk3MXQRGxeSRNy4KSUTIukEKmOsbSaK0zAgIWB09NPFXMs4znhnzjn3/PjM+wVy7z3eO+fJ9eOTc4/fcx0RAgB038+aDgAAqAaFDgBJUOgAkASFDgBJUOgAkASFDgBJUOgAkASFDgBJUOgAkMTqOne2xg/HWq0r9NpNz9yoOA3aamp8baHX3dCPuhk/ueI4S8JsYykGPdulCt32bkl/lLRK0p8i4vCDnr9W6/Rbv1BoXwf+erHQ69A9R7ZvK/S6v8fZyjIw2xiEQc924VMutldJek/Sy5J2SNpne0fRnwe0BbONripzDn2npIsRcSkibkr6UNKeamIBjWK20UllCv1xSd/OeTzZ23YP2yO2x2yP3dJPJXYH1IbZRieVKfSFTtDf9128EXE0IoYjYvghPVxid0BtmG10UplCn5T05JzHT0i6Ui4O0ArMNjqpTKF/Lmm77adsr5H0uqRT1cQCGsVso5MKL1uMiBnbByV9qtmlXcci4uvKkgENYbbRVaXWoUfEaUmnK8oCtAazjS7i0n8ASIJCB4AkKHQASIJCB4AkKHQASKLWr8/d9MwNvlkOfR2YKDYj/9rb3NfQlpnt17Z+WXEatNbEs4VettTZ5ggdAJKg0AEgCQodAJKg0AEgCQodAJKg0AEgCQodAJKodR06gOp99E2xtc1YvrZfM8AROgAkQaEDQBIUOgAkQaEDQBIUOgAkQaEDQBIUOgAkwTp0oOVYZ94e/f4sml6nzhE6ACRBoQNAEhQ6ACRBoQNAEhQ6ACRBoQNAEhQ6ACRRah267cuSfpB0W9JMRAxXEQpoGrONLqriwqLnI+L7Cn4O0DbMNjqFUy4AkETZQg9Jn9n+wvbIQk+wPWJ7zPbY9emZkrsDasNso3PKnnJ5LiKu2N4k6Yztf0bEublPiIijko5K0rZfPRIl9wfUhdlG55Q6Qo+IK73bKUknJe2sIhTQNGYbXVS40G2vs/3o3fuSXpI0XlUwoCnMNrqqzCmXzZJO2r77cz6IiE8qSQU0i9lGJxUu9Ii4JOnXFWYBWoHZRlexbBEAkqDQASAJCh0AkqDQASAJCh0AkqDQASAJCh0AkqDQASAJCh0AkqDQASAJCh0AkqDQASAJCh0AkqDQASAJCh0AkqDQASAJCh0AkqDQASAJCh0AkqDQASAJCh0AkqDQASAJCh0AkqDQASAJCh0AkqDQASAJCh0AkqDQASAJCh0Akuhb6LaP2Z6yPT5n2wbbZ2xP9G7XDzYmUD1mG9ks5Qh9VNLuedsOSTobEdslne09BrpmVMw2Eulb6BFxTtL0vM17JB3v3T8uaW/FuYCBY7aRTdFz6Jsj4qok9W43LfZE2yO2x2yPXZ+eKbg7oDbMNjpr4P8oGhFHI2I4IoYf27B60LsDasNso22KFvo121skqXc7VV0koFHMNjqraKGfkrS/d3+/pI+riQM0jtlGZy1l2eIJSX+T9LTtSdtvSDos6UXbE5Je7D0GOoXZRjZ9T/xFxL5FfuuFirMAtWK2kQ1XigJAEhQ6ACRBoQNAEhQ6ACRBoQNAEhQ6ACTB9cor1NMPlbsA8sKtRb/iBGjUSp5tjtABIAkKHQCSoNABIAkKHQCSoNABIAkKHQCSYNniCvXLNY+Uev2FWxUFASpWdrYnbt1Z9PfutPwYuN3pAABLRqEDQBIUOgAkQaEDQBIUOgAkQaEDQBIUOgAkwTr0FerUj+XW6gJtVXa2277W/EG6mxwAcA8KHQCSoNABIAkKHQCSoNABIAkKHQCSYNniCtXlpVkrzWtbv3zg73/0zbM1JemGQc52vz+LpvX9L7d9zPaU7fE52962/Z3t871frww2JlA9ZhvZLOV/ZaOSdi+w/d2IGOr9Ol1tLKAWo2K2kUjfQo+Ic5Kma8gC1IrZRjZlTjYdtP1V72Pr+sWeZHvE9pjtsevTMyV2B9SG2UYnFS30I5K2ShqSdFXSO4s9MSKORsRwRAw/toF/g0XrMdvorEKFHhHXIuJ2RNyR9L6kndXGAprBbKPLChW67S1zHr4qaXyx5wJdwmyjy/p+TrR9QtIuSRttT0p6S9Iu20OSQtJlSW8OMCMwEFlmu+1ro1GfvoUeEfsW2PznAWQBasVsIxsuFwSAJCh0AEiCQgeAJCh0AEiCQgeAJCh0AEii1uuVp8bX6sj2bYVee2DiYsVp0FZFZ2Qq/l1xkmXsu8Rsa4LvM18pBj3bHKEDQBIUOgAkQaEDQBIUOgAkQaEDQBIUOgAkQaEDQBIUOgAkQaEDQBIUOgAkQaEDQBIUOgAkQaEDQBIUOgAk4Yiob2f2fyXN/R7IjZK+ry3A0pFredqS6xcR8fMmdjxvttvyfsxHruVrS7YlzXathX7fzu2xiBhuLMAiyLU8bc3VlLa+H+RavjZnWwinXAAgCQodAJJoutCPNrz/xZBredqaqyltfT/ItXxtznafRs+hAwCq0/QROgCgIo0Uuu3dti/Yvmj7UBMZFmL7su1/2D5ve6zhLMdsT9ken7Ntg+0ztid6t+tbkutt29/13rfztl+pO1dbMNt9czDXA1R7odteJek9SS9L2iFpn+0dded4gOcjYqgFS5VGJe2et+2QpLMRsV3S2d7juo3q/lyS9G7vfRuKiNM1Z2oFZntJRsVcD0wTR+g7JV2MiEsRcVPSh5L2NJCj1SLinKTpeZv3SDreu39c0t5aQ2nRXJjFbPfBXA9WE4X+uKRv5zye7G1rg5D0me0vbI80HWYBmyPiqiT1bjc1nGeug7a/6n10rf0jc0sw28Uw1xVpotC9wLa2LLV5LiJ+o9mPzH+w/bumA3XEEUlbJQ1JuirpnWbjNIbZzqVzc91EoU9KenLO4yckXWkgx30i4krvdkrSSc1+hG6Ta7a3SFLvdqrhPJKkiLgWEbcj4o6k99W+960uzHYxzHVFmij0zyVtt/2U7TWSXpd0qoEc97C9zvajd+9LeknS+INfVbtTkvb37u+X9HGDWf7v7l/GnlfVvvetLsx2Mcx1RVbXvcOImLF9UNKnklZJOhYRX9edYwGbJZ20Lc2+Lx9ExCdNhbF9QtIuSRttT0p6S9JhSX+x/Yak/0j6fUty7bI9pNnTC5clvVl3rjZgtvtjrgeLK0UBIAmuFAWAJCh0AEiCQgeAJCh0AEiCQgeAJCh0AEiCQgeAJCh0AEjif1vVgWLjT56TAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"engine": 0
},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAACzCAYAAACZ+efrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAACjdJREFUeJzt3M+PVfUZx/HPBxRbibYSkRo1jZ2iLWnDUG9wYdJgDAbdgE1MZMXCBBflD3CnSzfGdGFMsCXDRk1NAFkQlbBh0xgvUSyaWhxC6whhMBOlMSmKPF1waUaYmXvnnnPPj2fer4Sce869Z75PDs98OPfwPccRIQBA+y2ruwAAQDkIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCRuqHKwFb4pfqSVQ+173wO/KDT2Vxc/KbQ/BvfTm9YV2v+fx04Ntd9/9Y2+jYsuNPiQ6O2loem97SK3/tveIulPkpZL+nNEvLDQ52/1qnjQjww11uHLbw6131X7JjcU2h+D+8PYB4X237zsyaH2ey+O6ELMlBLobertA5PrC+2PwW0bO15o/1H39tCXXGwvl/SypMckrZO03Xaxf76ABqC30VZFrqFvlPRZRJyKiG8lvSFpazllAbWit9FKRQL9Lkmfz1qf6m37Ads7bXdtd7/TxQLDAZWht9FKRQJ9rus5112Qj4jdEdGJiM6NuqnAcEBl6G20UpFAn5J0z6z1uyWdKVYO0Aj0NlqpSKC/L2mt7Xttr5D0lKSD5ZQF1IreRisNPQ89Ii7Z3iXpHV2Z2rUnIj4urbKSFZ1Kh6Wjbb1ddCod8ih0Y1FEHJJ0qKRagMagt9FG3PoPAEkQ6ACQBIEOAEkQ6ACQBIEOAEkUetriYnU6neh2u5WNh6XF9rGI6NQxNr2NURq0tzlDB4AkCHQASIJAB4AkCHQASIJAB4AkCHQASIJAB4AkCj1tsU32TW6ou4Qlg0cVA/XgDB0AkiDQASAJAh0AkiDQASAJAh0AkiDQASAJAh0Akqh0HvpXFz9ZcD54kfnLzDNvjn5/F6P6ex77zY8fGPoHF/TVxU90YHL9vO9vGzteYTVoo4X6Z9De5gwdAJIg0AEgCQIdAJIg0AEgCQIdAJIg0AEgiSXz+FwA11toqlxWmaeQFgp026cl/UfS95IuRUSnjKKAutHbaKMyztAfjogvS/g5QNPQ22gVrqEDQBJFAz0kvWv7mO2dc33A9k7bXdvdr2cuFRwOqMyievsCvY0GKHrJ5aGIOGP7DkmHbf8jIo7O/kBE7Ja0W5J++dubo+B4QFXobbROoTP0iDjTW05L2i9pYxlFAXWjt9FGQwe67ZW2b7n6WtKjkk6UVRhQF3obbVXkkssaSfttX/05r0XE26VUBdQrTW8vxXnm/fQ7Jm2epz50oEfEKUl0C9Kht9FWTFsEgCQIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCT6BrrtPbanbZ+YtW2V7cO2T/aWt422TKB89DayGeQMfULSlmu2PSvpSESslXSktw60zYTobSTSN9Aj4qikmWs2b5W0t/d6r6RtJdcFjBy9jWyGvYa+JiLOSlJvecd8H7S903bXdvfrmUtDDgdUZqjevkBvowFG/p+iEbE7IjoR0fnJqhtGPRxQmdm9fSu9jQYYNtDP2b5TknrL6fJKAmpFb6O1hg30g5J29F7vkPRWOeUAtaO30VqDTFt8XdLfJN1ve8r205JekLTZ9klJm3vrQKvQ28im74W/iNg+z1uPlFwLUCl6G9lwpygAJEGgA0ASBDoAJEGgA0ASBDoAJEGgA0AS3K+8RN1/Y7EbID/9bt5HnAC1+tWK83WXUBvO0AEgCQIdAJIg0AEgCQIdAJIg0AEgCQIdAJJg2uIS9esVNxfa/+R3l+d97zLnCY2xbez4gu8fmFxfUSXVue/GlQu+v+xnJyuqpHr85gFAEgQ6ACRBoANAEgQ6ACRBoANAEgQ6ACRBoANAEsxDX6IOflNsHjpzzXPoN08d7cJvJQAkQaADQBIEOgAkQaADQBIEOgAkQaADQBKOiMoG63Q60e12Kxtvtn2TG2oZdyn6w9gHtYxr+1hEdOoYu87eRn6D9nbfM3Tbe2xP2z4xa9vztr+w/WHvz+NFCwaqRm8jm0EuuUxI2jLH9pciYrz351C5ZQGVmBC9jUT6BnpEHJU0U0EtQKXobWRT5D9Fd9n+qPe19bb5PmR7p+2u7e758+cLDAdUht5GKw0b6K9IGpM0LumspBfn+2BE7I6ITkR0Vq9ePeRwQGXobbTWUIEeEeci4vuIuCzpVUkbyy0LqAe9jTYbKtBt3zlr9QlJJ+b7LNAm9DbarO/jc22/LmmTpNttT0l6TtIm2+OSQtJpSc+MsMZS1DU3Gs2VpbeBq/oGekRsn2PzX0ZQC1ApehvZcOs/ACRBoANAEgQ6ACRBoANAEgQ6ACRBoANAEpU+D/1Wr4oH/chQ+x6+/GbJ1aCpNi97cqj93osjuhAzLrmcgdDbGMSoe5szdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCT6Pm0xi32TG+ouYcngUcXVOjC5vu4SloxtY8frLmFBnKEDQBIEOgAkQaADQBIEOgAkQaADQBIEOgAkQaADQBKVPj7X9nlJ/5q16XZJX1ZWwOCoa3GaUtfPI2J1HQNf09tNOR7Xoq7Fa0ptA/V2pYF+3eB2NyI6tRUwD+panKbWVZemHg/qWrwm1zYXLrkAQBIEOgAkUXeg7655/PlQ1+I0ta66NPV4UNfiNbm269R6DR0AUJ66z9ABACWpJdBtb7H9qe3PbD9bRw1zsX3a9t9tf2i7W3Mte2xP2z4xa9sq24dtn+wtb2tIXc/b/qJ33D60/XjVdTUFvd23Dvp6hCoPdNvLJb0s6TFJ6yRtt72u6joW8HBEjDdgqtKEpC3XbHtW0pGIWCvpSG+9ahO6vi5Jeql33MYj4lDFNTUCvT2QCdHXI1PHGfpGSZ9FxKmI+FbSG5K21lBHo0XEUUkz12zeKmlv7/VeSdsqLUrz1oUr6O0+6OvRqiPQ75L0+az1qd62JghJ79o+Zntn3cXMYU1EnJWk3vKOmuuZbZftj3pfXSv/ytwQ9PZw6OuS1BHonmNbU6baPBQRv9OVr8x/tP37ugtqiVckjUkal3RW0ov1llMbejuX1vV1HYE+JemeWet3SzpTQx3XiYgzveW0pP268hW6Sc7ZvlOSesvpmuuRJEXEuYj4PiIuS3pVzTtuVaG3h0Nfl6SOQH9f0lrb99peIekpSQdrqOMHbK+0fcvV15IelXRi4b0qd1DSjt7rHZLeqrGW/7v6y9jzhJp33KpCbw+Hvi7JDVUPGBGXbO+S9I6k5ZL2RMTHVdcxhzWS9tuWrhyX1yLi7bqKsf26pE2Sbrc9Jek5SS9I+qvtpyX9W9KTDalrk+1xXbm8cFrSM1XX1QT0dn/09WhxpygAJMGdogCQBIEOAEkQ6ACQBIEOAEkQ6ACQBIEOAEkQ6ACQBIEOAEn8D/m8021OLK8PAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"engine": 0
},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAACzCAYAAACZ+efrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAACXJJREFUeJzt3M2PVHUaxfFzANGEOAmEFomaiUMwEzZ2nAouTAyGaMANuCCRFQuTdiF/ADtcsjFmFsYEZ0izEaMLlAVRSW/YGYsEnXYxAxJGW5Buw2Imk4wIPLPoYtL0W3XfunVfHr6fpFNVt6v6d3LzcHLrcqscEQIAtN+augMAAMpBoQNAEhQ6ACRBoQNAEhQ6ACRBoQNAEhQ6ACRBoQNAEhQ6ACSxrsrF1vvheEQbCr32mT/9oeQ0aKp/XLhS6HX/1X90K351yXFWhNnGSgx7tgcqdNt7JP1Z0lpJf4mIY8s9/xFt0PPeXWitc91PCr0O7fPymgOFXvdVTJSWgdnGMAx7tgufcrG9VtJ7kvZK2iHpoO0dRf8e0BTMNtpqkHPoOyVdjogrEXFL0keS9pUTC6gVs41WGqTQn5D045zHU71t97E9Zrtru/ubfh1gOaAyzDZaaZBCX+wE/YLv4o2I4xHRiYjOQ3p4gOWAyjDbaKVBCn1K0lNzHj8p6dpgcYBGYLbRSoMU+teSttt+2vZ6Sa9LOlNOLKBWzDZaqfBlixFx2/ZhSV9o9tKuExHxXWnJgJow22irga5Dj4izks6WlAVoDGYbbcRH/wEgCQodAJKg0AEgCQodAJKg0AEgCUcs+ADc0HQ6neh2u5WthweL7QsR0aljbWYbw7TS2eYIHQCSoNABIAkKHQCSoNABIAkKHQCSoNABIAkKHQCSoNABIAkKHQCSoNABIAkKHQCSoNABIAkKHQCSoNABIAkKHQCSoNABIAkKHQCSoNABIAkKHQCSoNABIAkKHQCSoNABIIl1dQdog0+/f7buCJXbv+2buiMAWKWBCt32VUn/lnRH0u2I6JQRCqgbs402KuMI/aWI+KWEvwM0DbONVuEcOgAkMWihh6QvbV+wPbbYE2yP2e7a7s7MzAy4HFAZZhutM2ihvxARz0naK+kt2y/Of0JEHI+ITkR0RkZGBlwOqAyzjdYZqNAj4lrvdlrSaUk7ywgF1I3ZRhsVLnTbG2w/eu++pFckTZYVDKgLs422GuQqly2STtu+93c+jIjPS0lVsQfxOvN++u2T5Nepp5ltPFgKF3pEXJFEEyIdZhttxWWLAJAEhQ4ASVDoAJAEhQ4ASVDoAJAEX58LoJC7P2+vbe01j1+qbe0m4wgdAJKg0AEgCQodAJKg0AEgCQodAJKg0AEgCQodAJLgOnQAS6rzWvPlLJfrQb5GnSN0AEiCQgeAJCh0AEiCQgeAJCh0AEiCQgeAJCh0AEiCQgeAJCh0AEiCQgeAJCh0AEiCQgeAJCh0AEiCQgeAJCh0AEiib6HbPmF72vbknG2bbJ+zfal3u3G4MYHyMdvIZiVH6OOS9szbdkTSRERslzTRewy0zbiYbSTSt9Aj4rykm/M275N0snf/pKT9JecCho7ZRjZFz6FviYjrktS7fWypJ9oes9213Z2ZmSm4HFAZZhutNfT/FI2I4xHRiYjOyMjIsJcDKsNso2mKFvoN21slqXc7XV4koFbMNlqraKGfkXSod/+QpM/KiQPUjtlGa63r9wTbpyTtkrTZ9pSko5KOSfrY9huSfpB0YJghgWFgtvtb8/ilJX939+ftFSa533K5HmR9Cz0iDi7xq90lZwEqxWwjGz4pCgBJUOgAkASFDgBJUOgAkASFDgBJUOgAkETfyxaR0x/X890jGAzXgjcPR+gAkASFDgBJUOgAkASFDgBJUOgAkASFDgBJcNmipP3bvln2959+/2xFSarzzEMblv09l6QB7cMROgAkQaEDQBIUOgAkQaEDQBIUOgAkQaEDQBIUOgAkwXXoK9DvOnUAaAKO0AEgCQodAJKg0AEgCQodAJKg0AEgCQodAJKg0AEgib6FbvuE7Wnbk3O2vW37J9sXez+vDjcmUD5mG9ms5Ah9XNKeRba/GxGjvZ+z5cYCKjEuZhuJ9C30iDgv6WYFWYBKMdvIZpBz6Idtf9t727pxqSfZHrPdtd2dmZkZYDmgMsw2Wqloob8vaZukUUnXJb2z1BMj4nhEdCKiMzIyUnA5oDLMNlqrUKFHxI2IuBMRdyV9IGlnubGAejDbaLNChW5765yHr0maXOq5QJsw22izvl+fa/uUpF2SNtueknRU0i7bo5JC0lVJbw4xIzAUzDay6VvoEXFwkc1/HUIWoFLMNrLhk6IAkASFDgBJUOgAkASFDgBJUOgAkASFDgBJOCIqW+x33hTPe3eh1567+0nJadBUL685UOh1X8WE/hU3XXKcFWG2sRLDnm2O0AEgCQodAJKg0AEgCQodAJKg0AEgCQodAJKg0AEgCQodAJKg0AEgCQodAJKg0AEgCQodAJKg0AEgCQodAJKo9Otzbc9I+uecTZsl/VJZgJUj1+o0JdfvI2KkjoXnzXZT9sd85Fq9pmRb0WxXWugLFre7EdGpLcASyLU6Tc1Vl6buD3KtXpOzLYZTLgCQBIUOAEnUXejHa15/KeRanabmqktT9we5Vq/J2Rao9Rw6AKA8dR+hAwBKUkuh295j+++2L9s+UkeGxdi+avtvti/a7tac5YTtaduTc7Ztsn3O9qXe7caG5Hrb9k+9/XbR9qtV52oKZrtvDuZ6iCovdNtrJb0naa+kHZIO2t5RdY5lvBQRow24VGlc0p55245ImoiI7ZImeo+rNq6FuSTp3d5+G42IsxVnagRme0XGxVwPTR1H6DslXY6IKxFxS9JHkvbVkKPRIuK8pJvzNu+TdLJ3/6Sk/ZWG0pK5MIvZ7oO5Hq46Cv0JST/OeTzV29YEIelL2xdsj9UdZhFbIuK6JPVuH6s5z1yHbX/be+ta+VvmhmC2i2GuS1JHoXuRbU251OaFiHhOs2+Z37L9Yt2BWuJ9SdskjUq6LumdeuPUhtnOpXVzXUehT0l6as7jJyVdqyHHAhFxrXc7Lem0Zt9CN8kN21slqXc7XXMeSVJE3IiIOxFxV9IHat5+qwqzXQxzXZI6Cv1rSdttP217vaTXJZ2pIcd9bG+w/ei9+5JekTS5/Ksqd0bSod79Q5I+qzHL/937x9jzmpq336rCbBfDXJdkXdULRsRt24clfSFpraQTEfFd1TkWsUXSadvS7H75MCI+ryuM7VOSdknabHtK0lFJxyR9bPsNST9IOtCQXLtsj2r29MJVSW9WnasJmO3+mOvh4pOiAJAEnxQFgCQodABIgkIHgCQodABIgkIHgCQodABIgkIHgCQodABI4n+7BIZRkAQARgAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"engine": 0
},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAACzCAYAAACZ+efrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAACTNJREFUeJzt3M+LVfcdxvHnUWsCkoLixEgSShqmCzcd0otZBIJBEjQbk4UQVy4Ck0X8A9yZZTYhdBECppVxE0OykLiQJDIbdyFXkHSyaBWxyUTjTHDRUmis+ulijmWcX3fm3HPPj0/fL5B77/Gc+T5cPj6cezx3HBECAHTfpqYDAACqQaEDQBIUOgAkQaEDQBIUOgAkQaEDQBIUOgAkQaEDQBIUOgAksaXOxbb6kXhU20od+7s//LbiNGirv126Vuq4f+tfuhO/uOI468JsYz1GPdtDFbrtA5L+KGmzpD9FxLtr7f+otul57y+11oX+Z6WOQ/e8vOlwqeO+junKMjDbGIVRz3bpSy62N0v6QNJBSXskHbG9p+zPA9qC2UZXDXMNfa+kqxFxLSLuSPpE0qFqYgGNYrbRScMU+pOSflj0erbY9hDbk7b7tvv/0S9DLAfUhtlGJw1T6CtdoF/2u3gj4mRE9CKi9ys9MsRyQG2YbXTSMIU+K+npRa+fknRjuDhAKzDb6KRhCv0bSeO2n7G9VdIbks5VEwtoFLONTip922JE3LV9TNKXWri161REfFdZMqAhzDa6aqj70CPivKTzFWUBWoPZRhfx1X8ASIJCB4AkKHQASIJCB4AkKHQASMIRy74ANzK9Xi/6/X5t6+H/i+1LEdFrYm1mG6O03tnmDB0AkqDQASAJCh0AkqDQASAJCh0AkqDQASAJCh0AkqDQASAJCh0AkqDQASAJCh0AkqDQASAJCh0AkqDQASAJCh0AkqDQASAJCh0AkqDQASAJCh0AkqDQASAJCh0AkqDQASAJCh0AktgyzMG2r0v6p6R7ku5GRK+KUEDTmG100VCFXngpIn6u4OcAbcNso1O45AIASQxb6CHpK9uXbE+utIPtSdt92/35+fkhlwNqw2yjc4Yt9Bci4jlJByW9bfvFpTtExMmI6EVEb2xsbMjlgNow2+icoQo9Im4Uj3OSzkraW0UooGnMNrqodKHb3mb7sQfPJb0iaaaqYEBTmG101TB3ueySdNb2g5/zcUR8UUkqoFnMNjqpdKFHxDVJv68wC9AKzDa6itsWASAJCh0AkqDQASAJCh0AkqDQASCJKn45F9Zw/6fxxtbe9MSVxtYGUD/O0AEgCQodAJKg0AEgCQodAJKg0AEgCQodAJKg0AEgCe5Dr0CT95qvZa1c3KMO5MMZOgAkQaEDQBIUOgAkQaEDQBIUOgAkQaEDQBIUOgAkQaEDQBIUOgAkQaEDQBIUOgAkQaEDQBIUOgAkQaEDQBIUOgAkMbDQbZ+yPWd7ZtG2HbYv2L5SPG4fbUygesw2slnPGfqUpANLth2XNB0R45Kmi9dA10yJ2UYiAws9Ii5Kur1k8yFJp4vnpyW9VnEuYOSYbWRT9hr6roi4KUnF4+Or7Wh70nbfdn9+fr7kckBtmG101sj/UzQiTkZELyJ6Y2Njo14OqA2zjbYpW+i3bO+WpOJxrrpIQKOYbXRW2UI/J+lo8fyopM+riQM0jtlGZ20ZtIPtM5L2Sdppe1bSCUnvSvrU9puSvpd0eJQh227TE1dW/bv7P43XmORha+UCs418BhZ6RBxZ5a/2V5wFqBWzjWz4pigAJEGhA0ASFDoAJEGhA0ASFDoAJEGhA0ASA29bxHC4FxxAXThDB4AkKHQASIJCB4AkKHQASIJCB4AkKHQASIJCB4AkKHQASIJCB4AkKHQASIJCB4AkKHQASIJCB4AkKHQASIJCB4AkKHQASIJCB4AkKHQASIJCB4AkKHQASIJCB4AkKHQASGJgods+ZXvO9syibe/Y/tH25eLPq6ONCVSP2UY26zlDn5J0YIXt70fERPHnfLWxgFpMidlGIgMLPSIuSrpdQxagVsw2shnmGvox298WH1u3r7aT7Unbfdv9+fn5IZYDasNso5PKFvqHkp6VNCHppqT3VtsxIk5GRC8iemNjYyWXA2rDbKOzShV6RNyKiHsRcV/SR5L2VhsLaAazjS4rVei2dy96+bqkmdX2BbqE2UaXbRm0g+0zkvZJ2ml7VtIJSftsT0gKSdclvTXCjMBIMNvIZmChR8SRFTb/eQRZgFox28iGb4oCQBIUOgAkQaEDQBIUOgAkQaEDQBIUOgAk4YiobbFfe0c87/2ljr1w/7OK06CtXt50uNRxX8e0/hG3XXGcdWG2sR6jnm3O0AEgCQodAJKg0AEgCQodAJKg0AEgCQodAJKg0AEgCQodAJKg0AEgCQodAJKg0AEgCQodAJKg0AEgCQodAJKo9dfn2p6X9PdFm3ZK+rm2AOtHro1pS67fRMRYEwsvme22vB9LkWvj2pJtXbNda6EvW9zuR0SvsQCrINfGtDVXU9r6fpBr49qcbSVccgGAJCh0AEii6UI/2fD6qyHXxrQ1V1Pa+n6Qa+PanG2ZRq+hAwCq0/QZOgCgIo0Uuu0Dtv9q+6rt401kWInt67b/Yvuy7X7DWU7ZnrM9s2jbDtsXbF8pHre3JNc7tn8s3rfLtl+tO1dbMNsDczDXI1R7odveLOkDSQcl7ZF0xPaeunOs4aWImGjBrUpTkg4s2XZc0nREjEuaLl7XbUrLc0nS+8X7NhER52vO1ArM9rpMibkemSbO0PdKuhoR1yLijqRPJB1qIEerRcRFSbeXbD4k6XTx/LSk12oNpVVzYQGzPQBzPVpNFPqTkn5Y9Hq22NYGIekr25dsTzYdZgW7IuKmJBWPjzecZ7Fjtr8tPrrW/pG5JZjtcpjrijRR6F5hW1tutXkhIp7Twkfmt22/2HSgjvhQ0rOSJiTdlPRes3Eaw2zn0rm5bqLQZyU9vej1U5JuNJBjmYi4UTzOSTqrhY/QbXLL9m5JKh7nGs4jSYqIWxFxLyLuS/pI7Xvf6sJsl8NcV6SJQv9G0rjtZ2xvlfSGpHMN5HiI7W22H3vwXNIrkmbWPqp25yQdLZ4flfR5g1n+58E/xsLrat/7VhdmuxzmuiJb6l4wIu7aPibpS0mbJZ2KiO/qzrGCXZLO2pYW3pePI+KLpsLYPiNpn6SdtmclnZD0rqRPbb8p6XtJh1uSa5/tCS1cXrgu6a26c7UBsz0Ycz1afFMUAJLgm6IAkASFDgBJUOgAkASFDgBJUOgAkASFDgBJUOgAkASFDgBJ/BdGyHgH7VyOJQAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"engine": 0
},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAACzCAYAAACZ+efrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAACPxJREFUeJzt3M+LlIcdx/HPR7cmICkobkSSUNJgC166pIM5BIJBDJqLyUGIJw+BzSH+Ad6So5cQegiBTbuslxiag8SDJJG9eCkhI0i6ObSK2GajuBs8tBQaq3572Mey2R/u7DPPPD++fb9AZubJzD5fhm/ezI4zOiIEAOi+LU0PAACoBkEHgCQIOgAkQdABIAmCDgBJEHQASIKgA0ASBB0AkiDoAJDEWJ0n2+bH4nFtL/XYX/32lxVPg7b66+XrpR73b/1Ld+NHVzzOQNhtDGLUuz1U0G0flvQ7SVsl/T4iTj/q/o9ru17wwVLnutj/tNTj0D2Hthwr9bivYrayGdhtjMKod7v0Wy62t0r6QNIRSfskHbe9r+zPA9qC3UZXDfMe+n5J1yLiekTclfSJpKPVjAU0it1GJw0T9Kckfbfs9nxx7CdsT9ru2+7/Rz8OcTqgNuw2OmmYoK/1Bv2qf4s3IqYiohcRvZ/psSFOB9SG3UYnDRP0eUnPLLv9tKSbw40DtAK7jU4aJuhfS9pr+1nb2yS9Iel8NWMBjWK30UmlP7YYEfdsn5T0hZY+2jUdEd9WNhnQEHYbXTXU59Aj4oKkCxXNArQGu40u4qv/AJAEQQeAJAg6ACRB0AEgCYIOAEk4YtUX4Eam1+tFv9+v7Xz4/2L7ckT0mjg3u41RGnS3eYUOAEkQdABIgqADQBIEHQCSIOgAkARBB4AkCDoAJEHQASAJgg4ASRB0AEiCoANAEgQdAJIg6ACQBEEHgCQIOgAkQdABIAmCDgBJEHQASIKgA0ASBB0AkiDoAJAEQQeAJAg6ACQxNsyDbd+Q9E9J9yXdi4heFUMBTWO30UVDBb3wckT8UMHPAdqG3Uan8JYLACQxbNBD0pe2L9ueXOsOtidt9233FxcXhzwdUBt2G50zbNBfjIjnJR2R9Lbtl1beISKmIqIXEb3x8fEhTwfUht1G5wwV9Ii4WVwuSDonaX8VQwFNY7fRRaWDbnu77SceXpf0iqS5qgYDmsJuo6uG+ZTLbknnbD/8OR9HxOeVTAU0i91GJ5UOekRcl/SbCmcBWoHdRlfxsUUASIKgA0ASBB0AkiDoAJAEQQeAJAg6ACRB0AEgCYIOAEkQdABIgqADQBIEHQCSIOgAkARBB4AkCDoAJEHQASAJgg4ASRB0AEiCoANAEgQdAJIg6ACQBEEHgCQIOgAkQdABIAmCDgBJEHQASIKgA0ASBB0AkiDoAJAEQQeAJDYMuu1p2wu255Yd22n7ou2rxeWO0Y4JVI/dRjaDvEKfkXR4xbFTkmYjYq+k2eI20DUzYreRyIZBj4hLku6sOHxU0pni+hlJr1U8FzBy7DayKfse+u6IuCVJxeWT693R9qTtvu3+4uJiydMBtWG30Vkj/0vRiJiKiF5E9MbHx0d9OqA27DbapmzQb9veI0nF5UJ1IwGNYrfRWWWDfl7SieL6CUmfVTMO0Dh2G501yMcWz0r6k6Rf2563/aak05IO2b4q6VBxG+gUdhvZjG10h4g4vs5/OljxLECt2G1kwzdFASAJgg4ASRB0AEiCoANAEgQdAJIg6ACQBEEHgCQIOgAkQdABIAmCDgBJEHQASIKgA0ASBB0AkiDoAJAEQQeAJAg6ACRB0AEgCYIOAEkQdABIgqADQBIEHQCSIOgAkARBB4AkCDoAJEHQASAJgg4ASRB0AEiCoANAEgQdAJLYMOi2p20v2J5bduxd29/bvlL8eXW0YwLVY7eRzSCv0GckHV7j+PsRMVH8uVDtWEAtZsRuI5ENgx4RlyTdqWEWoFbsNrIZ5j30k7a/KX5t3bHenWxP2u7b7i8uLg5xOqA27DY6qWzQP5T0nKQJSbckvbfeHSNiKiJ6EdEbHx8veTqgNuw2OqtU0CPidkTcj4gHkj6StL/asYBmsNvoslJBt71n2c3XJc2td1+gS9htdNnYRnewfVbSAUm7bM9LekfSAdsTkkLSDUlvjXBGYCTYbWSzYdAj4vgah/8wglmAWrHbyIZvigJAEgQdAJIg6ACQBEEHgCQIOgAkQdABIAlHRG0n+7l3xgs+WOqxFx98WvE0aKtDW46VetxXMat/xB1XPM5A2G0MYtS7zSt0AEiCoANAEgQdAJIg6ACQBEEHgCQIOgAkQdABIAmCDgBJEHQASIKgA0ASBB0AkiDoAJAEQQeAJAg6ACRR6z+fa3tR0t+WHdol6YfaBhgcc21OW+b6RUSMN3HiFbvdludjJebavLbMNtBu1xr0VSe3+xHRa2yAdTDX5rR1rqa09flgrs1r82xr4S0XAEiCoANAEk0Hfarh86+HuTanrXM1pa3PB3NtXptnW6XR99ABANVp+hU6AKAijQTd9mHbf7F9zfapJmZYi+0btv9s+4rtfsOzTNtesD237NhO2xdtXy0ud7Rkrndtf188b1dsv1r3XG3Bbm84B3s9QrUH3fZWSR9IOiJpn6TjtvfVPccjvBwREy34qNKMpMMrjp2SNBsReyXNFrfrNqPVc0nS+8XzNhERF2qeqRXY7YHMiL0emSZeoe+XdC0irkfEXUmfSDrawBytFhGXJN1ZcfiopDPF9TOSXqt1KK07F5aw2xtgr0eriaA/Jem7Zbfni2NtEJK+tH3Z9mTTw6xhd0TckqTi8smG51nupO1vil9da/+VuSXY7XLY64o0EXSvcawtH7V5MSKe19KvzG/bfqnpgTriQ0nPSZqQdEvSe82O0xh2O5fO7XUTQZ+X9Myy209LutnAHKtExM3ickHSOS39Ct0mt23vkaTicqHheSRJEXE7Iu5HxANJH6l9z1td2O1y2OuKNBH0ryXttf2s7W2S3pB0voE5fsL2dttPPLwu6RVJc49+VO3OSzpRXD8h6bMGZ/mfh/8zFl5X+563urDb5bDXFRmr+4QRcc/2SUlfSNoqaToivq17jjXslnTOtrT0vHwcEZ83NYzts5IOSNple17SO5JOS/qj7Tcl/V3SsZbMdcD2hJbeXrgh6a2652oDdntj7PVo8U1RAEiCb4oCQBIEHQCSIOgAkARBB4AkCDoAJEHQASAJgg4ASRB0AEjiv/hca15VZOTsAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"engine": 0
},
"output_type": "display_data"
}
],
"source": [
"%%px\n",
"from GPUSimulators import KP07_dimsplit\n",
"\n",
"sim = KP07_dimsplit.KP07_dimsplit(**arguments)\n",
"\n",
"grid = getFactors(MPI.COMM_WORLD.size, 2)\n",
"rank = MPI.COMM_WORLD.rank\n",
" \n",
"gc_x = int(sim.u0[0].x_halo)\n",
"gc_y = int(sim.u0[0].y_halo)\n",
"nx = int(sim.nx)\n",
"ny = int(sim.ny)\n",
"\n",
"n_steps = 1\n",
"t_end = 5*dt\n",
"for i in range(n_steps):\n",
" t_this = t_end/n_steps\n",
" n_steps = int(t_this/dt)\n",
" for j in range(n_steps):\n",
" #Swap east-west boundary\n",
" read_e = [nx , gc_y, gc_x, ny]\n",
" write_e = [nx+gc_x, gc_y, gc_x, ny]\n",
" out_e = sim.u0[0].download(sim.stream, async=False, extent=read_e)\n",
" in_e = np.empty((ny, gc_x), dtype=np.float32)\n",
" send_e = MPI.COMM_WORLD.Isend(out_e, dest=getEast(rank, grid), tag=1)\n",
" recv_e = MPI.COMM_WORLD.Irecv(in_e, source=getEast(rank, grid), tag=0)\n",
" #print(\"Send 1 to \", getEast(rank, grid), \" receive 0 from \", getEast(rank, grid))\n",
"\n",
" read_w = [gc_x, gc_y, gc_x, ny]\n",
" write_w = [ 0, gc_y, gc_x, ny]\n",
" out_w = sim.u0[0].download(sim.stream, async=False, extent=read_w)\n",
" in_w = np.empty((ny, gc_x), dtype=np.float32)\n",
" send_w = MPI.COMM_WORLD.Isend(out_w, dest=getWest(rank, grid), tag=0)\n",
" recv_w = MPI.COMM_WORLD.Irecv(in_w, source=getWest(rank, grid), tag=1)\n",
" #print(\"Send 0 to \", getWest(rank, grid), \" receive 1 from \", getWest(rank, grid))\n",
"\n",
" #Swap north-south boundary\n",
" read_n = [gc_x, ny , nx, gc_y]\n",
" write_n = [gc_x, ny+gc_y, nx, gc_y]\n",
" out_n = sim.u0[0].download(sim.stream, async=False, extent=read_n)\n",
" in_n = np.empty((gc_y, nx), dtype=np.float32)\n",
" send_n = MPI.COMM_WORLD.Isend(out_n, dest=getNorth(rank, grid), tag=3)\n",
" recv_n = MPI.COMM_WORLD.Irecv(in_n, source=getNorth(rank, grid), tag=2)\n",
"\n",
" read_s = [gc_x, gc_y, nx, gc_y]\n",
" write_s = [gc_x, 0, nx, gc_y]\n",
" out_s = sim.u0[0].download(sim.stream, async=False, extent=read_s)\n",
" in_s = np.empty((gc_y, nx), dtype=np.float32)\n",
" send_s = MPI.COMM_WORLD.Isend(out_s, dest=getSouth(rank, grid), tag=2)\n",
" recv_s = MPI.COMM_WORLD.Irecv(in_s, source=getSouth(rank, grid), tag=3)\n",
"\n",
" \n",
" recv_e.wait()\n",
" recv_w.wait()\n",
" recv_n.wait()\n",
" recv_s.wait()\n",
" \n",
" sim.u0[0].upload(in_e, sim.stream, extent=write_e)\n",
" sim.u0[0].upload(in_w, sim.stream, extent=write_w)\n",
" sim.u0[0].upload(in_n, sim.stream, extent=write_n)\n",
" sim.u0[0].upload(in_s, sim.stream, extent=write_s)\n",
" \n",
" send_e.wait()\n",
" send_w.wait()\n",
" send_n.wait()\n",
" send_s.wait()\n",
" \n",
" if (rank == 0):\n",
" h1 = sim.u0[0].download(sim.stream, extent=[0, 0, nx+2*gc_x, ny+2*gc_y])\n",
" plt.figure()\n",
" plt.subplot(1,2,1)\n",
" plt.imshow(h1)\n",
" \n",
" #print(sim.simTime())\n",
" sim.simulate(dt, dt=dt)\n",
" \n",
" if (rank == 0):\n",
" h1 = sim.u0[0].download(sim.stream, extent=[0, 0, nx+2*gc_x, ny+2*gc_y])\n",
" plt.subplot(1,2,2)\n",
" plt.imshow(h1)\n",
" \n",
" \n",
" h1 = sim.u0[0].download(sim.stream)\n",
" #plotSolution(h1, grid)\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python [default]",
"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.6.7"
}
},
"nbformat": 4,
"nbformat_minor": 2
}