{
"cells": [
{
"cell_type": "markdown",
"id": "25dd0180",
"metadata": {},
"source": [
"# Successive Over-relaxation (SOR)\n",
"\n",
"In this section, we want to examine another parallel algorithm called Successive Over-relaxation (SOR). \n",
"\n",
"The SOR algorithm is an iterative method used to solve Laplace equations. The underlying data structure of the SOR algorithm is a two-dimensional grid, whose elements are updated iteratively through some weighted function that considers the old value as well as the values of neighbouring cells. \n",
"\n",
"\n",
"\n",
"This algorithm is applied, for instance, in physics to simulate the climate or temperature of some object.\n",
"\n",
"\n",
"Climate Simulation | Temperature of a metal plate\n",
":------:|:----:\n",
" | \n",
"\n",
"\n",
"## Sequential algorithm\n",
"\n",
"The sequential SOR algorithm is as follows (where `f` is some update function and `N`,`M` are the grid sizes):\n",
"\n",
"```julia\n",
"grid = zeros(N, M)\n",
"grid_new = zeros(N, M)\n",
"for step in 1:NSTEPS\n",
" for i in 2:N #update grid\n",
" for j in 2:M\n",
" grid_new[i,j] = f(grid[i,j], grid[i-1,j], grid[i+1,j], grid[i,j-1], grid[i,j+1])\n",
" end\n",
" end\n",
"end\n",
"```\n",
"\n",
"## Diffusion equation on grid \n",
"\n",
"Consider the diffusion of a chemical substance on a two-dimensional grid. The concentration of the chemical is given as $c(x,y)$, a function of the coordinates $x$ and $y$. We will consider a square grid with $0 \\leq x,y \\leq 1$ and the boundary conditions $c(x,y=1) = 1$ and $c(x,y=0) = 0$. That is, the concentration at the top of the grid is always 1 and the concentration at the very bottom is always 0. Furthermore, in the x-direction we will assume periodic boundary conditions, i.e. $c(x=0, y) = c(x=1,y)$. \n",
"\n",
"\n",
"\n",
"We will take the initial condition $c(x,y) = 0$ for $0 \\leq x \\leq 1, 0 \\leq y < 1$.\n",
"\n",
"The stable state of the diffusion, that is, when the concentraion does not change anymore, can be described by the Laplace equation\n",
"\n",
"$$\n",
"\\nabla^2 c = 0.\n",
"$$\n",
"\n",
"Numerically, we can approximate the solution with the Jacobi iteration \n",
"\n",
"$$\n",
"\\frac{1}{4}(c^k_{i,j+1} + c^k_{i,j-1} + c^k_{i+1,j} + c^k_{i-1,j}) = c^{k+1}_{i,j}.\n",
"$$\n",
"\n",
"The superscript $k$ denotes the $k$-th iteration. The algorithm stepwise updates the cells of the grid, until a steady state is reached. To determine the end of the algorithm, we use the stopping condition \n",
"\n",
"$$\n",
"\\max_{i,j} \\lvert c^{k+1}_{i,j} - c^{k}_{i,j} \\rvert < \\epsilon.\n",
"$$\n",
"\n",
"That is, we stop when all changes to cell values are smaller than some small number, say $\\epsilon = 10^{-5}$.\n",
"\n",
"Furthermore, for this set of initial and boundary conditions, there exists an analytical solution for the stable state, namely \n",
"$$\n",
"c(x,y) = y.\n",
"$$\n",
"That is, the concentration profile is the identity function of the y-coordinate. "
]
},
{
"cell_type": "markdown",
"id": "84167bc6",
"metadata": {},
"source": [
"## Serial Program"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "b082bce8",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"SOR! (generic function with 1 method)"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Update grid until threshold reached\n",
"function update_grid!(new_grid, grid, M, N)\n",
" @assert grid[1,:] == grid[M-1,:]\n",
" @assert grid[M,:] == grid[2,:] \n",
" @assert grid[:,1] == fill(0,M) == new_grid[:,1]\n",
" @assert grid[:,N] == fill(1,M) == new_grid[:,N]\n",
" # Remove ghost cells\n",
" g_left = grid[1:M-2, 2:N-1]\n",
" g_right = grid[3:M, 2:N-1]\n",
" g_up = grid[2:M-1, 3:N]\n",
" g_down = grid[2:M-1, 1:N-2]\n",
" # Jacobi iteration\n",
" new_grid[2:M-1,2:N-1] = 0.25 * (g_up + g_down + g_left + g_right)\n",
" # Update ghost cells\n",
" new_grid[1,:] = new_grid[M-1,:]\n",
" new_grid[M,:] = new_grid[2,:]\n",
" return new_grid\n",
"end\n",
"\n",
"function SOR!(new_grid, grid, ϵ, M, N)\n",
" grid_old = true\n",
" while true\n",
" if grid_old\n",
" update_grid!(new_grid, grid, M, N)\n",
" else \n",
" update_grid!(grid, new_grid, M, N)\n",
" end\n",
" # Flip boolean value\n",
" grid_old = !grid_old\n",
" diffs = abs.(new_grid[2:M-1, :] - grid[2:M-1, :])\n",
" if maximum(diffs) < ϵ\n",
" break\n",
" end\n",
" end \n",
"end\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "e280865a",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"\u001b[32m\u001b[1mTest Passed\u001b[22m\u001b[39m"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Test, Plots\n",
"N = 100\n",
"M = N + 2\n",
"ϵ = 10^-5\n",
"\n",
"# initialize grid\n",
"grid = zeros(M,N)\n",
"grid[:,N] .= 1\n",
"new_grid = zeros(M,N)\n",
"new_grid[:,N] .= 1\n",
"\n",
"update_grid!(new_grid, grid, M, N)\n",
"\n",
"# Test if first iteration successful\n",
"@test all(new_grid[:,N-1] .≈ fill(0.25, M))\n",
"@test all(new_grid[:,N-2] .≈ fill(0,M))\n",
"@test all(new_grid[:,N] == fill(1,M))\n",
"@test all(new_grid[:,1] == fill(0,M))"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "f2f49b93",
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3deXxN19748XVOxiYigwwiQUhMiSgNV2MuidQQhMRMW8NFr7Z6taWtlF5tUUq1fTpcJbQNVZSqSy/VGi8tLR6hpTQhSIKMSqZz9v79sXvPL09EZJ9sdpLzeb/OHysr65y9nOv2a33XZJBlWQAAYKuMencAAAA9EQgBADaNQAgAsGkEQgCATSMQAgBsmr3eHQAA4I7y8vKOHj36+++/d+3aNSwsrMI2KSkpn3/+ub29/bhx44KDg9U+ghEhAKDm6tu376xZsxITE7/77rsKGxw/frxr166Ojo7FxcWdOnVKTU1V+whGhACAmuvQoUN2dnbR0dF3avDWW29NmTLllVdeEUJkZWW9++67S5cuVfWIKo0Ir127tnTp0jFjxjz66KNl64uKimbOnBkaGtq7d+99+/ZZ6g8fPhwVFdWmTZunnnrq1q1bqjoEAICFnZ1d5Q327t3bt29fpRwTE7Nnzx61j6hSILx8+fLJkycDAwN3795dtn7u3Lk//fTTV199NWHChEGDBl27dk0IkZ+f379//5EjR3799dfnzp174YUX1PYJAICqkGU5KyvL19dX+dHX1zcjI0Pth1QpNdq+ffukpKSTJ0+WHW+WlpZ+/PHH27Zta9GiRYsWLZKTkz/55JOZM2euXbs2LCxs0qRJQog333yzS5cub775pouLy+0fm5GRsW3b55MmP66207ZAkiQhhNHIJG5lzGbzXf+1aONMJpO9PTMgleEruivLV2QQnlVoXiALs6rPP3Pm3CuJi8tVRkdHT548+a7vNRgMdnZ2JpPJ0lUHBwdVTxfVmSO8cuVKTk5ORESE8mPHjh1TUlKEECkpKZbK8PBwk8mUmppa4VKfCxcutG1XJIntVvehLjMKIYSkdy9qOIMdX9FdGO35iu6Cr+iuLF+RnRhz18aSfEAWuao+//Tpn9LS0p5//vmylW3btq3i2xs1anTlypWHHnpICHH58uVGjRqperqoTiC8du2ai4uLo6Oj8qOHh8fPP/+s1IeGhlqaeXh4XL16tcJAmJWV5exSYHUHAAD3U2lp6V3HW5JsloVJ1cfKwhwQEJCQkFD1t2RnZ589ezYyMlIIMXDgwE2bNg0cOFAIsWnTptjYWFVPF9UJhB4eHkVFRZbc1I0bNzw9PZX6mzdvWppZ6m/n5+dnlutZ3QEAwP1UlayjLMyyrC41KsuVjckXLFjw3Xff/fzzz+np6Vu3bk1MTOzRo8fBgwcnTZp09epVIcTf//73Ll26xMfHFxYWnjt3bvXq1aqeLqoTCBs1auTg4PDbb7+1bt1aCHH27FllG2OzZs0OHTqktLl06VJRUVGTJk3u9CGyXGoysawUAGoBOz3mUkeOHGlZFCqEaN68uRCie/fu27f/Oa3WpEmTU6dO7dq1y97ePjo62tXVVe0jqvTHkmU5Ly+voKBACJGbm2s0Gt3d3V1cXIYOHbps2bKPPvro7Nmz//rXv3744QchxOjRo19//fWTJ0+Gh4cvXbp0wIABXl5earsFAKiNJNksyypTo5WOIJs1a9asWbNylZ6enh07drT86O7uHh8fr+qhZVUpEGZnZ7ds2VII4ebmFhwc7O/vf+rUKSHE4sWL4+Pj/fz8SkpK5s+frwwNmzZtunTp0p49ezo6OgYGBm7cuNHqzgEAahkrAqHKVaaaq1Ig9Pb2zsnJub0+ICDg0KFDBQUFDzzwQNnc8V//+tcnnniisLCwfv36mvUUAIB7QIOMb4XRzsHBoSrTqmapyGS+eddmAADdOVUhYsiSSfWIUNJ5AwvbSAEAmpHramoUAICqkIXqQChUbrfQnM6BUDYXSaY/9O0DAKBKnPTuwL3BiBAAoBlZMsmS2jlC2x4RAgDqFNks1KZG9T7tlcsNAAA2Tfc5wkKZOUIAqDO0PlnmPiA1CgDQjmwWKucIbX3VKACgTpFMQipV+RbbDoQGc6GhVPVtwgAAaIURIQBAO7JJ/apRsxB296QzVUMgBABox4o5QknnQMj2CQCATdN7RGgqEqV6/kMAAKAliVWjAABbZsXJMjLXMAEA6gzZZFA5IjToPSJkjhAAYNP03kdoKjKW6NsFAIB2JEn1BnluqAcA1BkGyUxqFACA2kT31GiJoVTWtw8AAM3IZtWpUVaNAgDqDiv2Edr4odsAgLrEIJkNWo8Ir1y5snLlytzc3CFDhvTo0eP2Bnl5eUlJSRcvXuzcufOIESMMBoOq5zNHCACoufLy8v7yl79kZGQEBQUNHTr066+/Ltfgxo0bHTt2PHXqVHh4+NKlS1988UW1j9B9jrDYUKL2nHIAQI1lxckylY0gP/nkk5YtW77//vtCCBcXlwULFsTGxpZtsG7dOg8Pj48//lgI0adPn9DQ0FmzZnl6elb9+YwIAQCaMUiSkh1V8ao0Nbpv377o6GilHBUVdfjw4ZKS/7P9XBksKuUmTZqUlpb++OOPqvpMIAQA1FwZGRk+Pj5K2dfXV5blzMzMsg3atWv3ww8/3Lx5UwixZ8+e0tLSK1euqHqE7qnREmOpvl0AAGhHsmb7xMGDBy3DPsWQIUP+9re/CSEcHBxMpj9zrUrB0dGxbMvBgwevX78+PDy8Xbt26enpTZo0cXFxUfV8Vo0CADRjkCW1q0YNktSiRYtZs2aVrWzWrJlSCAgIuHz5slK+dOmSg4ODZYCoMBqNn3/++blz53JycsLCwgICAlq0aKGqAwRCAIB2rBoR+vr6RkVFVfjLQYMGzZ8/PzEx0dHR8YsvvhgwYICdnZ0Q4siRIw0bNmzcuLHSLCQkRAixdOnSwMDADh06qHo+gRAAUHPFxcW9//773bt3b968+e7du3fu3KnUT58+fdSoUTNmzBBCtG7dukOHDpcvX05NTd2+fbvafYS6zxGWGkp0PlwHAKAZyZrUaCW/dXR03L179549e/Lz8999911vb2+lfs2aNV5eXkp5y5YtJ0+edHV17dWrl9oJQqF7IAQA1CmypPpaJfkuJ07b29vfnjht3bp12XLZH9Vi+wQAwKYxIgQAaMZgxcW8tn77RGkJR6wBQN1hxTVM3FAPAKg7auGIkDlCAIBN03v7hLnUUFpy93YAgNrAIKleNVr5odv3AalRAIB2ZPWpUb3nCEmNAgBsGiNCAIB21KdG77qh/l7TOxCWmkQJc4QAUEcYJFn1HCGpUQAAdKT3iBAAUJdYsaHexleNGkylBlKjAFBnqE+N2vwcIQCgLrHi9glJ50DIHCEAwKYxIgQAaMYgMUeoltkkOGINAOoMa+YI2T4BAIB+9B4RAgDqklq4WIZACADQTi1MjeodCEtNoqRU5z4AALRixVmjeo8ImSMEANg0vUeEAIA6RVKf6rTxOUKzWZSadO4DAEArVswRcvsEAAA60ntECACoS+7Bxbzp6en//Oc/8/PzhwwZ0rt379sbXLt2LSkp6eLFi0FBQRMmTPDy8lL1fEaEAADtKPsIVb0qDYS5ubmdO3cuKCho06bNiBEjtmzZUq7BzZs3O3fu/Ouvv3bt2vXYsWORkZHFxcWquqzziFA2yXKpztlhAIBmJPXbISptv2bNmtDQ0OXLlwshnJycFi5cOGTIkLINTpw4kZubu3LlSoPBMHz4cDc3t9OnT3fo0KHqz6/WiPDcuXPDhg0LCQnp0KHDO++8Y6nfs2dPZGRks2bNJk6ceOPGjeo8AgBgyw4cOBAVFaWU+/Tp8+OPP5Yb8AUHB0uSdPr0aSHE8ePHnZycgoKCVD2iWoFw5MiRgYGBR44cWbly5WuvvbZt2zYhRE5OzpAhQ6ZPn37w4MHr16/PnDmzOo8AANQmavOid0uNZmRkeHt7K2VfX19ZljMzM8s28PPzS05O7tKlS7Nmzfr06fPFF194enqq6nK1UqOnT5/+6KOPPD09PT09u3TpkpKSMnDgwLVr17Zv337MmDFCiAULFnTq1GnZsmWurq4Vf4RJFiU67yABAGhGllWnRmX54MGD0dHRZetiY2OffvppIYSjo6PJ9Ocuu5KSEiGEk5NT2ZYXLlyYOHHi//zP//To0ePbb78dP378zz//7O/vX/XnVysQPv7440uWLHn55ZfPnz9/7NixRYsWCSFOnToVERGhNAgNDZVlOTU1tW3bttV5EACgDmvRosWsWbPK1gQHByuFgICAS5cuKeVLly45Ojr6+PiUbfnVV1+FhoaOHTtWCDFhwoQVK1Zs27Zt8uTJVX96tQLhlClT4uPjR44cefXq1fHjx4eEhAghsrOz/fz8LG3c3d2vXbtW4duzsrIc8/NFI0N1+gAAuD9MJpO9/d2ihlUb6n19fS0TgeUMHjx43rx5iYmJTk5O69evHzhwoJ2dnRDi0KFDjRo1atq0qY+PT1paWklJiaOjY2Fh4cWLF8vGoKqwfo7w5s2bUVFRS5cuTUlJuXjx4g8//LBkyRIhhKen5x9//GFpVlBQcKctHX5+fvXr17e6AwCA++nuUVD8NzWq6lVpJnXIkCENGzbs0qVLfHz8Rx99NG/ePKV+xowZmzdvFkIMHTq0SZMmEREREydOjIiIaNu2bf/+/dX9uVS1ListLS0nJ6dfv35CCGdn56ioqKNHjwohmjdvvn//fqXNxYsXS0pKmjZtesdPMRvlUkaEAFBXSEKo3RNXaXsHB4edO3ceOHAgPz9/xYoVloUwycnJHh4eQggnJ6c9e/YcPXr00qVLTz75pGVuruqsHxE2a9bMzc0tOTlZCJGXl/f111+3b99eCDFmzBilT5IkLVq0aNCgQUpfAQCwgp2dXc+ePQcNGlR2OWhISIhlNanBYOjUqVNcXJwVUVBUJxC6uLisX79+0aJF/v7+wcHB4eHhyk6JwMDADz74oF+/fl5eXidPniy7vxAAUMfJsupXrb59IiYm5vTp08oUZdn6cePGjRs3rqioyNnZuXrdAwDUKlqnRu8DDY5YKxcFLaoUBU0GuZTzTgEAuuH2CQCAdmT1Izy9T1UhEAIANCPXvgvq9Q6EstlOLrXTtw8AAM1YMSLUe46Q+TkAgE0jNQoA0I4Vq0ZtPDUKAKhTJIOQVJ4XJut8vpjec4QmO7mUYAwA0A1BCACgGVk2yCpHeJXey3s/EAgBANqphatG9Q6EZiOpUQCoOyRR6+YI2T4BALBpjMYAAJqRZaMsqRtiyXpPEhIIAQDasWL7hI3PEUomO4kj1gAA+mFECADQjmxQv/jFtjfUAwDqElkSsurUKKtGAQDQj95HrJmNsolRKQDUFbJRqFw1qvs+QoIQAEA7kkFtalTtkWyaIxACALRjxWIZGz9rVDbbSRyxBgDQD0EIAKAZWTaoPllG7Zyi1giEAADtWHMx711+n5aW9sEHH9y4cWPw4MExMTHlfvvbb799+eWXZWtGjx7duHHjqj+f7RMAgJorJyfn4YcfNpvNHTp0GD9+/KZNm8o1KCkpyf2v48ePv/LKK87OzqoewRFrAADNyLLGq0ZXr1794IMPLlmyRAjh4ODw5ptvDhs2rGyDsLCwhQsXKuXnnntu0KBBPj4+qjrAiBAAoBlZNqp9Vb7K9ODBg71791bKvXv3PnLkSHFxcYUtTSZTcnLyhAkT1PaZQAgA0I7032nCqr8qDYQZGRne3t5K2cfHR5blzMzMClt+/fXXdnZ2ffv2Vdtl3bdPGCVOlgEA23bw4MHo6OiyNQMGDJgxY4YQwsnJqaSkRKlUCneaAly1atUTTzxhZ6d6uo0gBADQjGzFyTKSoUWLFrNmzSpb2aJFC6UQGBh46dIlpZyenu7o6FjhFGBmZubOnTuXLVtmRZ8JhAAADRmFrHbSzeDr6xsVFVXh74YMGTJnzpzExERnZ+d169YNGjTIaDQKIfbv3x8QENC8eXOlWVJSUteuXUNCQqzqMQAANdXgwYObNGnSuXPnwYMHr1q1at68eUr9c889t3XrVkuzTz75xIplMgq9t0+Y7cwmSd8+AAA0Y1VqtJJBmb29/Y4dOw4fPpyXl/fJJ5+4u7sr9evXr69fv75SNpvN//rXv1Rtov8/j7DubQAA3M6ai3nvdki30Wjs0qVLucqgoCBL2c7OzpIjtQKBEACgodp3+wRzhAAAm6b7HCH7CAGg7pAlo+rbJ1SvMtUYQQgAoBlZNqi9cV4mNQoAgI50T43amU36dgEAoBkrTpZRvbhGa6RGAQDakQ1C7Y3zagOn1giEAADtWDFHKHQOhMwRAgBsWk24hknfLgAANGPd7RP3qDNVRGoUAKAZWahOjeqO1CgAwKbpvX1CMppNtezfDgCAO5Elg+qTZdSuMtUaqVEAgGasOVlG71WjBEIAgGas2VCv96W0zBECAGyaziNCs9nObNa3CwAADZEaBQDYMsmg+sg0vfcRkhoFANg0RoQAAM3IskHtRbu6b8CvAfsIzewjBIA6gu0TAACbJstsnwAAoFbRffuE0WQiGANAHUFqFABg02RZSGoDoY0vlgEAoHLnzp378MMPc3Nz4+LiBg4cWGGb/fv3f/HFFzdv3uzYseOTTz6p6vNJSwIANCPLRllS+ap0RJidnd2lSxcHB4cePXpMnjz5iy++uL3NypUrhw0bFhgY2KtXr7Nnz6rtc03YPmGnbx8AAFqxZo6w0vZJSUkRERELFiwQQhgMhjfffHP48OFlG+Tl5c2YMWPXrl0PP/ywEGL8+PFq+1zdEeFnn33WoUMHb2/vzp07nz9/Xqn897//3b59ez8/v1GjRuXm5lbzEQAAm3Xo0KFevXop5V69ev3888/FxcXlGjRs2NDe3v7VV19999138/Pz1T6iWoFw7dq1L7744uLFi3/99de3337b1dVVCHH9+vWEhIR58+adOXPGbDb//e9/r84jAAC1iDIiVPUScmUfmJGR4e3trZR9fHxkWc7MzCzbIC0tLT8//9lnn/X19f3Pf/7TqVOnW7duqepztVKjc+fOXbx4cVRUlBDC0tHk5OROnToNGTJECDF//vz27du/8847bm5uFX6CJNmRGgWAOkOWhNoN9bJs2LNnT8eOHctWjhgx4vnnnxdCODs7l5SUKJXKWNDZ2blsS0dHx+vXr584ccLf33/q1KlhYWGbNm0aN25c1TtgfSAsKCg4d+5cUVFRr169SkpKxo4dO23aNIPB8Msvv7Rv315p06pVK4PBkJaWFh4ebvWDAAC1hlVzhO3bt1+yZEnZykaNGimFwMDA9PR0pZyenu7k5OTj41O2ZWBgoIuLi7+/vxDCYDAEBwdnZGSo6oD1gfDKlStCiKSkpBUrVmRnZyckJLi5uY0bNy4nJ6dhw4aWZu7u7tevX6/wE7KysvLy8kSDBlb3AQBw35SWljo4ONyLT/bw8IiIiKjwV3FxcbNnz05MTHzggQeSk5MHDx5sNBqFEN9//33jxo1DQkJ69uzp4uJy5MiRTp063bhx48iRI9OnT1f1dOvnCD09PYUQL730UsuWLSMjI6dOnbpp0yYhhJeX1x9//GFplp+f3+AOoc7Pz8/Dw8PqDgAA7qeqREFZqJ8jrFRsbGxISEinTp369+//6aefzps3T6mfPXv2tm3bhBDOzs7Lli0bOHDg8OHD27dvP2DAgJiYGFV/LutHhD4+Pp6enpZcrbOzc2lpqRAiJCRk9+7dSmVqaqrZbG7atOmdPsQsG00Sc4QAUEfIslH9NUyVtbe3t9+2bdvRo0fz8vIiIyPr1aun1G/cuNGy+mTMmDGPPPLIyZMnX3vttZYtW6rts/UjQqPROGHChHfffbekpOT69eurV6/u16+f0qH9+/fv37+/tLT0tddei4uLc3d3t/opAAAbZzAYOnXqFB0dbYmCQojGjRuXzSk2atQoJibGiigoqrl94tVXX3VycvL19W3Xrl1sbOy0adOEEP7+/klJSWPGjPH09Lx8+fLy5cur8wgAQC0iSQa1L9XXNmmtWtsnXF1dk5OTb69PSEhISEiozicDAGolWfUh2pVuI7wfdL+GyWBSmU0GANRYymIZdW/R+/YJghAAwKZxDRMAQDOaH7p9H9SA2ydktk8AQB0hyQbVF/Peo65UGalRAIBNIzUKANAMqVEAgG1Tv31C2HggNMtGE3OEAFBXWHPEmmD7BAAA+iE1CgDQjCQbJLUX89bqI9aqzywZTYxKAaCusGaxzD3qSpURhAAANo3UKABAO2yfAADYMkkWqk+WsfFAyPYJAKhLrLl94h51pcqYIwQA2DRSowAAzcjCoH6DvG2nRgEAdYms/vYJte01p/c1TLLBJJGeBQDohhEhAEAzVmyot/VDtwEAdYkV2yckvZeN6r19QjKaSY0CQF1h1RFr3D4BAIB+SI0CADRjxarRu44gz549+9577+Xl5cXFxcXFxd32dvnFF1+0/BgZGTl48GBVHWBECADQjHKyjMpXZR947dq1rl27enh4PProo9OnT1+3bt3tbRYtWuTs7Ozp6enp6eni4qK2z2yfAADUXElJSZ06dfrHP/4hhJAkafHixaNGjbq92TPPPOPp6WndI0iNAgA0Y832iUoXyxw+fLhXr15KuWfPnuPHjy8qKnJ2di7X7I033nBwcOjevXu/fv1UPp3UKABAO7IwSOpflXxgZmZmgwYNlLK3t7csy5mZmeXajBo1KjAw0NnZedKkSTNnzlTbZ/23T5AaBYA6Q5bVX6skiz179nTs2LFsXXx8/OzZs4UQzs7OJSUlSmVxcbEQotwsoMFgWLt2rVIePHhwhw4dXn75ZS8vr6o/n9QoAEBn7du3X7JkSdmagIAApRAYGHjx4kWlfPHiRWdnZ29v7zt9Ttu2be3s7DIyMgiEAAB9SFYduu3h4REREVHhb4cNG/b8888nJia6uLh8+umncXFxRqNRCLFr166mTZu2bNkyNzfXzc3N3t5eCJGcnOzq6hocHKyqAwRCAIBmZFlUvh1CrYEDByYlJT300EOBgYFnzpzZvXu3Uj9nzpxRo0a1bNly+/btzz33XGhoaEFBQWpq6urVq29fSlM5vecIZYNJZo4QAFAxOzu7LVu2HD9+PD8//y9/+csDDzyg1G/ZssXV1VUIMWbMmK5du6amptarV69Nmzb16tVT+whGhAAAzcjintxH2L59+3I1/v7+lnJQUFBQUJCqh5ZFIAQAaMaqG+p1RloSAGDTdJ8jZB8hANQd1q0avUedqSJSowAAzVixoV73VCqBEACgGfluR6ZV9Bad6Z4aZfsEAEBPjAgBAJqxYkO9thvwrUAgBABoxorFMuqvbdIYaUkAgE2rAXOEUi3begkAuBNZqL6Yl1WjAIC6w6pVowRCAEBdIcmqN8jb+mIZs2ww6T1NCgCwZYwIAQCaYfsEAMCm1cY5QrZPAABsWg2YI+T2CQCoK2SZ7RMAABsmCSGpf4u+CIQAAO2oHxHqfv0EaUkAgE2rAXOE7CMEgLqCG+oBADbNijlCvTOjpEYBALZN5xGhJAtunwCAOsOa2ydIjQIA6gyZ1CgAwJbJ8p976lW87vaZp0+fnjZt2ujRo7/44otKmp05c2b27NnHjh1T22cCIQCg5rp69Wr37t0bNWo0ZMiQmTNnfvbZZxU2M5vNEydOXLFixalTp9Q+gu0TAADNWLF9ovI5wlWrVkVGRiYmJgohTCbT4sWLx44de3uzZcuWdenSpbCwUNWjFRqMCNPT0x9++OF58+ZZar766quWLVu6u7sPHjz4+vXr1X8EAKC2kNW/KvHjjz/27NlTKffo0ePEiRNFRUXl2qSmpq5evXru3LnWdViDQDh16lRJktLS0pQfs7Kyxo4d+95772VmZrq7u8+YMaP6jwAA2KaMjAwvLy+l7O3tLctyZmZm2QaSJD3xxBPLli1zdXW17hHVTY2uWbPG19c3PDzc0rPk5OTIyMi+ffsKIV555ZW2bdu+//779evXr/DtZsnA9gkAqDOsS43u2bOnY8eOZSuHDRv24osvCiFcXFyKi4uVSmUs6OLiUrblBx980KxZs+joaKv7XK1AmJmZ+cYbbxw4cOCtt96yVJ45c6Zdu3ZKOSQkxGg0pqWlWWoAAHWYdSfLtG/ffsmSJWUrGzdurBQCAwMvXryolC9cuODs7Ozt7V225Y4dO/bu3fvVV18JIQoKCqZOnXr8+PFyn1a5agXCv/3tb6+++qqPj0/ZytzcXH9/f8uP7u7u2dnZFb49KysrLy/P+4GKB4sAgBrFZDLZ298lalh1H6Hw8PCIiIio8Lfx8fHPPvtsYmKiq6vrmjVrhg0bZjQahRA7duwICgpq06bN559/XlpaqjTu1avXk08+OW7cOFUdsD4Qbtu27fz58y1atPjpp58yMzOzs7NTUlLatm3r5eV148YNS7P8/Pxy0dvCz8/Pw8NDFFvdBQDA/XPXKHgv9O/ff82aNQ8++GBgYGBqauru3buV+nnz5o0aNapNmzb16tUr20NXV9dyudO7sv5P9ccff9jb20+ZMkUIcfny5dLS0jlz5mzZsqVVq1Y7duxQ2pw7d06SpKZNm97pQ8xCmHQ/VAAAoBErTpapvL2dnd3GjRtPnTqVl5fXsWNHJycnpX7btm23B7ytW7e6u7urfH41AuHIkSNHjhyplGfPnp2Zmbl69WohxOjRo1955ZVvvvmmR48er776akJCwp1WygAA6hgrzhoVVWgfFhZWrqbcrJwiICBA3aOFEFqdLBMQEBAUFKSU/fz81q5dO2PGjMDAwFu3br399tuaPAIAgHtBm4TvU089VfbH2NjY2NhYTT4ZAFCLSLKQVE54qW2vOY5YAwBopjbePsE1TAAAzSi3T6h7i9B5OMTtEwAAm6Z3alTihnoAqDusO1lGX6RGAQCakWUhq4xsattrjtQoAMCmMSIEAGhGFgZJ5eIXvQeEegdCjlgDgLrEin2EpEYBANATqVEAgGZk9alOvQeEegdCs0xqFADqDln9DfVq5xQ1x4gQAKAZ2Yo5P72HQ8wRAgBsGiNCAIBmrLjNTy8AABWrSURBVDhZRm17zekcCCVZmHX/DgAAWpHvycW89xSpUQCATSM1CgDQDKlRAIBNs+ZkmXvTk6pjHyEAQDO1cUM9c4QAAJtGahQAoBnZqkO39V02qndqVBIm3edJAQAase5iXn0DIalRAIBNIzUKANCMdA8u5j158uTy5cvz8/OHDBkyZsyYcr+9devWsmXLTp8+XVJS0q5du+nTp3t6eqrqACNCAIBm5P/uoFD1qkRWVlbPnj1DQkLGjBnz0ksvrVmzplyDmzdv3rhxY8iQIaNHjz548OCAAQPU9lnvI9a4oR4A6hArtk9UbuXKld26dZs9e7YQoqioaMGCBY899ljZBj4+PgsXLlTKERERTZs2zc/Pd3d3r/ojGBECAGquI0eOdO/eXSl369btf//3f4uKiipsWVBQ8Omnn3bs2LF+/fqqHsEcIQBAM1acLFP51oHMzMwGDRooZW9vb6UmKCioXLPw8PAzZ854enpu377dYFA3Scn2CQCAZqxLje7Zs6djx45la4YOHfrSSy8JIVxcXCxDwMLCQiGEq6vr7Z9w8uRJs9m8Zs2amJiYs2fPenl5Vf3pjAgBADpr3779kiVLytZYxnyNGze+cOGCUr5w4cIDDzxgGSCWY2dnN2HChJdffvnYsWN9+vSp+tMJhAAAzVhxsowkCw8Pj4iIiAp/m5CQ8NRTT82ZM8fNzS0pKSk+Pt5oNAohvv766+bNm4eFhV29etXd3d3JyUkI8cMPP+Tk5LRs2VJVBwiEAADNWHP7RKXt+/Xr17lz53bt2jVq1CgjI+Pbb79V6l977bVRo0aFhYXt2bNn2rRprVq1Ki4uTk1Nfe+99xo3bqyqA3rPEbJ9AgDqEM1vnzAajevWrfvtt99ycnI6dOjg6Oio1H/zzTfOzs5CiOHDh/ft2/fMmTMPPPBAcHBwhTOIlWNECACo6Vq0aFGupuzxMR4eHp07d7b6wwmEAADNWHP7xL3pSdURCAEAmqmNF/PqfcSaLMxqb+wAAEA7jAgBAJrRfNXofUAgBABoRlYf2Gw9EJol2aT2Hw8AgJpKvtvZoRW+RV/cPgEAsGmkRgEAmtH89on7gEAIANCMLGRZ/2SnOnrPEQrZpPs8KQDAhjEiBABoxrrbJ/RFIAQAaEaW1W+HsPFAaJZlA6lRAKgrJPWLX3RfLMP2CQCATSM1CgDQDLdPAABsmhXbJ3TfbqH7HGENOGYOAGDDGBECADQjqd8OoftoiEAIANBMbbyYl1WjAACbpvccoZBk/feQAAC0IcuypDLXKeudGyU1CgDQTG3cUE8gBABoRpZl3Ud4aumdGpVlSdb9XwMAANvFiBAAoBmZ1CgAwJZJstB8scyxY8eWL1+em5sbFxf3+OOPl/vt9evXk5OTDx8+XFpa2q1bt2nTpjk5OanqAIEQAFBzZWRk9O7de86cOa1bt54+fbokSRMmTCjb4Lvvvjty5MigQYOcnJzmz59/7NixNWvWqHqE3nOEQhLCrG8fAABa0XxD/apVq3r27Dlz5kwhxB9//PHaa6+VC4TDhw8fPny4Uvbz8+vbt+/q1asNBkPVO2D9hvpbt26NHTu2YcOGTk5Obdu23b59u+VX69evb9Sokaura1RUVGZmptWPAADULrKQJZWvyg/dPnr0aLdu3ZRyt27dUlJSCgsL79Q4NTU1ICBAVRQU1QmEJSUl4eHhR48eLSwsfOGFF+Lj47OysoQQV65cmThx4vr16/Pz84OCgp555hmrHwEAqF0kWVb7qnyKMDMz08vLSyk3aNBAqblTyxdeeOH1119X22frU6MeHh6zZs1SyuPHj3/mmWfOnDnj5+eXnJzcvXv37t27CyFefvnlVq1a5efnu7u7V/ghkkGSDKRGAcCm7dq1Kzg4uGzNuHHj5s2bJ4RwdXUtKipSKpWxYL169W7/hOzs7L59+06ZMiUhIUHt07WZI/zPf/4jSVJ4eLgQ4rfffmvbtq1S36xZMwcHh7S0tAcffLDCN9a6fZcAgEpYdx9hly5dPvzww7KVyuBPCNG4ceO0tDSlnJqa6uLi4u3tXe4T8vLyYmJi+vfvP3fuXCv6rMGh25mZmWPHjl26dKmnp6fSIVdXV8tv69evn5OTU+Ebr169WlBQUP0OAADuA5PJdNc28n9PWav6SxbC1dW1+f9lySMOHz58w4YNSrBYtWpVQkKCMgW4efPmkydPCiEKCgpiYmK6deu2cOFC6/5c1R0RXr9+PSoq6oknnpg4caJS06BBg7LhLT8/38fHp8L3+vr6uru7Szeq2QUAwP1gb6/DRoOYmJgePXqEhob6+/vn5eXt2rVLqV+4cOGoUaPCw8M/++yzH3/88ezZs5988onyq/PnzysDsyqq1p8qLy/v0Ucf7devX2JioqWydevWW7duVcpnzpwRQjRt2vROn2AWZjPbJwCgrlAWgqp9SyW/NRqNa9asSUtLKygoCAsLs7OzU+q//fZbR0dHIcTEiRNHjRpV9i0eHh6qOmB9ILx582ZUVJSfn9/IkSN/+uknIUTz5s09PT3HjBmTmJj45Zdf9urVa86cOSNHjnRzc7P6KQCAWkQWQv0c4d0FBQWVq7FEFicnJ7VHyZRj/RyhMvOXlZU15b+OHj0qhPD29t60adMbb7wRGhrq7Oy8dOnS6vQPAFCLSEJS+9L9VlrrR4SNGzdWIt/toqOjo6Ojrf5kAADuG72PWDOYTYa7L0MCANQKsvo5Qt130XHoNgBAM5ovlrkPNNhHCABA7aV3alRIJkFqFADqCCtOltE9OUpqFACgGVnIkspVoLqnRgmEAADNSAZZMqhcLKOyveaYIwQA2DS95wgNJpOBYAwAdYQsJLWpUfVzihojNQoA0IwsZLUnxeh+sgyjMQCATdM7NSqqcL0VAKCWkEiNAgBsmWSQJYPK7RN6rxolEAIANGPVYhnmCAEA0I/OI0JJmMx6Z4cBAFqR1d8vyBwhAKDukNQfsaZ7ICQ1CgCwaYwIAQCasWKxjNr2mqsJ+wh1/goAAFqRDJLa7RO6p0YZEQIANCMLSRZmtW+5R52pIuYIAQA2Te/tE7LJzO4JAKgrauOGelKjAADNWHX7BNsnAAC4s6NHj44ePXrAgAErVqyosMHOnTvfeOONKVOmnDlzxorPJxACADQjyWa1r8pHkFeuXImKinr44YeffvrphQsXVhgLX3311UuXLq1bt+7KlStW9FnvOUJhNuudHQYAaEXz1OjKlSt79+799NNPCyEWLlw4b968yZMnl2tz8OBBIcTmzZtVdvZPjAgBADXXTz/91KVLF6XcpUuX06dPFxYWavsIAiEAQDOyMEsqX5XvO8zKyvLy8lLKDRo0EEJkZmZq22e9U6NSqVlSt/USAFBjWXf7xK5du4KDg8tWjhkz5h//+IcQwtXV1TIEVApubm4adfZPbJ8AAGhGEpKk/mSZLl26fPjhh2Urvb29lUKTJk1SU1OV8u+//+7q6qqMCzVEahQAoDNXV9fm/1f9+vWVX40YMWLDhg15eXlCiI8//nj48OEGg0EIsWHDhhMnTmjydEaEAADNyLIkyypTo3Jlq0b79u0bHR3dpk0bPz+/4uLinTt3KvVLliwZNWrUgw8+KITo0aNHSkpKfn7+4MGD7e3tDx061KpVq6p3QPfbJ8xmuVTfPgAAtGPF9onK2hsMho8//vjy5ct5eXlt2rQxGv9MZO7du9fe/s8Qtm3bNrP5/+djLaPJKmJECACo6QICAgICAsrWODs7W8pqI185BEIAgGaUw2JUvUVtKlVzBEIAgGZq46HbOgdCWTZJzBECQF0hC7OsckQo9B4Rsn0CAGDTSI0CADSk8arR+0DvI9Zkk1ku0bcPAACtSLIkqd1HqPccIalRAIBNIzUKANCMLKTKb5Oo8C33qDNVRCAEAGhG8yPW7gO95wgF2ycAoC6x4homtk8AAKAfUqMAAM3Isqz6yDQbT43KklmSTPr2AQCgFWvmCNk+AQCAjkiNAgA0Y9X2CZVnk2qNQAgA0IwVqVGbnyOUzTLbJwCg7pCE6u0QzBECAKAfUqMAAM3IQvX2CW6oBwDUJVacLGPjc4TCxBwhAEBHjAgBANqRJaE61WnbI0IAQF1ixT5C9atMNaZ7IDTLgiPWAKDOsGL7BLdPAACgHwIhAEA7svznNKGK113mCA8fPpyQkBATE/P+++/fi1t8dU+NAgDqDlnI6rdDVNb+0qVLMTExCxcubNWq1dSpU41G49SpU6vTw9vpvX1ClmRZ5+NWAQA11sqVK6Ojo6dNmyaEWLBgQWJiouaBkNQoAEBDkvpXZSPCY8eORUZGKuXIyMhffvnl1q1b2vaYQAgA0JBs1euOMjMzPT09lbKXl5cQIisrS9se65kazcvLKy6+1bx5Ex37UGPl5+ebzWblf3VUqLCwMD8/v2HDhnp3pOYym81Xrlxp3Lix3h2p0S5evBgYGGg0Miq4o4yMDE9PT2dn58TExPnz51fe+Pjxn9V+/o4dOxISEoKDg8tWjh49WnmWm5tbYWGhUqmMBd3c3NQ+onJ6BsJHH330xIkTdnZ2OvahxjKbzbIs29uzmqkyxcXFTk5OeveiRuMruiu+oruyfEX+/v734vMfffTRlJQUSfo/uwl9fX2VQpMmTX7//Xel/Pvvv9erV69BgwbadsBwL5aiAgCgiV27dk2cOPH48eNeXl6TJ0+WZfnjjz/W9hEMOAAANVdUVFT//v3btGnj7e1tNBq/+eYbzR/BiBAAUNNlZWXl5eW1bNnSYDBo/uEEQgCATWOhFADApjFHWCPk5ORs27YtJSXFzc0tLi6ubdu2ll/t2rVrx44dfn5+kyZN0nytVG107dq1zZs3d+3aNSwsTKlJTU1dvXp1UVHRiBEjHnroIX27p7vz58+vW7cuOzs7LCzssccec3BwEEJcunQpKSmpoKBg6NChlr3JtkmSpA0bNhw9etTNzS0hIaFNmzZK/R9//PHPf/7z0qVL3bt3j4uL07eT919RUdGJEydOnTrVsGHD/v37W+pNJtPq1atPnz6t/HWyrGM/e/bsp59+ajKZRo8eHR4erlOvNcOIsEaYNWvWli1bfHx8CgoKOnfu/O9//1up//TTTx977LFmzZqdOnWqa9euxcXF+vazJpg+ffrzzz///fffKz9eunSpY8eOt27datCgwSOPPHLo0CF9u6evffv2RUREXLt2LSgo6MCBA8ruq+zs7E6dOl29etXf33/AgAG7du3Su5t6mj59+uuvv96mTRtJkjp16vTzzz8LIWRZjoqK2rdvX0hIyAsvvLBkyRK9u3m/LVy4cPz48cuWLVu+fHnZ+kmTJq1atapFixYrV66cMmWKUnn+/PnOnTvLslyvXr1u3bodP35cjy5rSkYNUFhYaCnPmjVr6NChsixLktSmTZsNGzYo5Q4dOnz22We6dbFm2Lp1a3x8fO/evd99912l5uWXXx4xYoRSnj9/flxcnH6905nZbA4JCUlKSipXv2jRon79+inld955p3fv3ve7ZzWJr6/v7t27lXJCQsLcuXNlWd61a1dgYGBJSYksy/v27fP19S0uLtaxk/efsnF5+fLlffv2tVReuHDByckpKytLluWsrCwnJ6f09HRZlp955pmJEycqbWbNmjV27Fg9uqwlRoQ1grOzs6VcVFTk6uoqhLh69eovv/wSFRUlhDAYDFFRUXv37tWtizVAXl7erFmz3n777bKV+/bt69u3r1KOjo625a/o7Nmz6enpPXv2/OCDD1auXJmXl6fU79u3Lzo6WilHR0fv37+/3M5lmxIaGnrs2DEhxK1bt3755Rclwb53795HHnlESSN37dr15s2bv/76q84dvb8qPFjnwIEDbdu2VTa2+/r6hoWFHTx4UNTF/9MRCGuWkydPJiUlzZw5UwiRkZHh6Ojo4eGh/MrPz+/KlSu69k5nzzzzzLPPPhsQEFC2MiMjw8fHRyn7+vrm5uYWFRXp0Tv9paWlOTs7x8fHZ2dn7969u127dtnZ2eK2r6i0tPT69eu69lRPycnJK1asaNGiRWBgYP/+/RMSEoQQmZmZlq/IaDR6e3vb+P/XFGW/FlHmP0Hl/kZlZmbKtXz3AYGwBrl48eKgQYPeeuutBx98UAjh4OBgNpst/3gvLS215YOgduzYceHChUmTJpWrt7e3N5lMStlkMhkMBps9l85oNObn5y9fvnzOnDlr164NDg5WDuAo9xUJIRwdHfXsqK4mT54cGRm5devWjRs3btiwYfPmzUIIe3t7s/n/3wdXWlpqy1+RxZ2+lnJ/o+zt7e/F3r77yUb/k1EDpaenP/LII88//7zlv/WNGjUym81ZWVnK+X6XL1++Rwf91QrJyckXLlzo1KmTEOLs2bO//fZbRkbG66+/HhAQYPnH++XLl/38/Gw2ECpj5dDQUOXH0NDQCxcuKPVlvyIXFxdLmsHWXLx4cfv27bm5uR4eHm3atPnrX/+6cuXKuLi4gICAEydOKG2Ki4uzs7MbNWqkb1drgoCAgMuXL1t+vHz5svK1lPsbVS5JUxsxIqwRsrKy+vbtO2XKlCeffNJS6enp2b17940bNwohCgsLt23bNmjQIP36qLP58+dv3Ljxo48++uijj1q3bj1q1CjlXwyxsbEbN25UMjMbNmyIjY3Vu6e6CQ0Nbdmy5eHDh4UQsiz/8MMPSlCMjY398ssvlX/a2/hX5OnpaTQaz58/r/x4/vx5b29vIURsbOy3336bk5MjhNi6dWvTpk1btWqlZ0drhj59+ly4cCElJUUIkZKSkp6e3rt3byFEbGzshg0blDZ15G+Uzot1IMuyLI8cOdLZ2Tniv8aMGaPU7927t0GDBuPGjevQocOAAQOUlV0ou2o0Pz+/Xbt2jzzySHx8vL+///nz5/Xtm762bt3q6+s7adKkrl27RkZGKquRb9261blz565du44cOdLHxyclJUXvbupp7ty5fn5+U6ZMGThwoK+vr+XbePzxx1u3bv344497e3tv3rxZ307efzt37oyIiGjcuHH9+vUjIiISExOV+kWLFgUEBEyYMCEgIGDx4sVK5fXr11u1ahUTEzNkyJAmTZooS0lrNY5YqxHOnz9vWeMnhHBxcbHs871y5YqynrtXr15cmaY4e/ash4eH5ZaWoqKi3bt3FxYWRkVF2WzSzyI9Pf3gwYP+/v7dunWz3HFWUlLy/fff5+fn9+nTh2MZzpw5c/LkyXr16nXt2rXszXYHDhxIT0+PjIwMCgrSr3f6yM3NtVx1JIRo0KCB5Us4ceKEsqG+Xbt2lga3bt369ttvzWZzVFSU5rcD3n8EQgCATWOEAQCwaQRCAIBNIxACAGwagRAAYNMIhAAAm0YgBADYNAIhAMCmEQgBADaNQAgAsGkEQgCATSMQAgBs2v8DTyGyCcgrPFYAAAAASUVORK5CYII=",
"image/svg+xml": [
"\n",
"\n"
],
"text/html": [
"\n",
"\n"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"N = 100\n",
"M = N + 2\n",
"ϵ = 10^-5\n",
"\n",
"# initialize grid\n",
"grid = zeros(M,N)\n",
"grid[:,N] .= 1\n",
"new_grid = zeros(M,N)\n",
"new_grid[:,N] .= 1\n",
"\n",
"SOR!(new_grid, grid, ϵ, M, N)\n",
"plt = heatmap(transpose(new_grid[2:M-1,:]))\n",
"display(plt)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "4aefd11d",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"\u001b[32m\u001b[1mTest Passed\u001b[22m\u001b[39m"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# analytical solution\n",
"function analytical_solution(N)\n",
" # Returns the analytical solution as a square grid of size N\n",
" grid = zeros(N,N)\n",
" for i in 2:N\n",
" grid[:,i] .= (i-1)/(N-1)\n",
" end\n",
" return grid\n",
"end\n",
"\n",
"# Test if solution is identical with analytical solution\n",
"sol = analytical_solution(N)\n",
"@test maximum(abs.(sol - new_grid[2:M-1,:])) < 0.01 * N"
]
},
{
"cell_type": "markdown",
"id": "f9f090f1",
"metadata": {},
"source": [
"## Parallel Program with MPI"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "49c4e9b4",
"metadata": {},
"outputs": [],
"source": [
"#] add MPIClusterManagers DelimitedFiles"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "6ec8e2eb",
"metadata": {
"code_folding": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" From worker 5:\tProc 4 gets rows 17:33 and columns 2:17\n",
" From worker 8:\tProc 7 gets rows 34:50 and columns 2:17\n",
" From worker 3:\tProc 2 gets rows 1:16 and columns 18:34\n",
" From worker 6:\tProc 5 gets rows 17:33 and columns 18:34\n",
" From worker 9:\tProc 8 gets rows 34:50 and columns 18:34\n",
" From worker 4:\tProc 3 gets rows 1:16 and columns 35:51\n",
" From worker 10:\tProc 9 gets rows 34:50 and columns 35:51\n",
" From worker 7:\tProc 6 gets rows 17:33 and columns 35:51\n",
" From worker 2:\tProc 1 gets rows 1:16 and columns 2:17\n",
" From worker 2:\tTHRESHOLD SUBCEEDED - TERMINATE SOR\n"
]
}
],
"source": [
"# to import MPIManager\n",
"using MPIClusterManagers\n",
"\n",
"# need to also import Distributed to use addprocs()\n",
"using Distributed\n",
"\n",
"# specify, number of mpi workers, launch cmd, etc.\n",
"manager=MPIWorkerManager(9)\n",
"\n",
"# start mpi workers and add them as julia workers too.\n",
"addprocs(manager)\n",
"\n",
"@mpi_do manager begin\n",
" \n",
"function calculate_partition(p, N, nrows, ncols)\n",
" # Calculates the row and column indices of this processor\n",
" # Get row and column number for processor p \n",
" if mod(p,ncols) == 0\n",
" i = div(p,ncols)\n",
" else\n",
" i = floor(div(p,ncols)) + 1\n",
" end\n",
" j = p - (i-1)*ncols\n",
" # Rows\n",
" if mod(N,nrows) == 0\n",
" prows = div(N, nrows)\n",
" row_range =((i-1) * prows + 1) : (i*prows)\n",
" else\n",
" # nlower processors get the smaller partition\n",
" nlower = nrows - (mod(N,nrows)) \n",
" n_floor = floor(div(N,nrows))\n",
" if i <= nlower\n",
" row_range =((i-1) * n_floor + 1) : (i*n_floor)\n",
" else\n",
" row_range = ((i-1) * n_floor + (i - nlower)) : (i*n_floor + (i-nlower))\n",
" end\n",
" end\n",
" # Columns\n",
" if mod(N,ncols) == 0\n",
" prows = div(N, ncols)\n",
" col_range =((j-1) * prows + 1) : (j*prows)\n",
" else\n",
" nlower = ncols - (mod(N,ncols)) \n",
" n_floor = floor(div(N,ncols))\n",
" if j <= nlower\n",
" col_range =((j-1) * n_floor + 1) : (j*n_floor)\n",
" else\n",
" col_range = ((j-1) * n_floor + (j - nlower)) : (j*n_floor + (j-nlower))\n",
" end\n",
" end\n",
" # Add 1 to each column index because of ghost cells\n",
" col_range = col_range .+ 1\n",
" return row_range, col_range\n",
"end\n",
"\n",
" \n",
"function update_grid(grid)\n",
" # Returns the updated grid as M-2 x N-2 matrix where M and N are sizes of grid\n",
" M = size(grid, 1)\n",
" N = size(grid, 2) \n",
" # Remove ghost cells\n",
" g_left = grid[2:M - 1,1:N-2]\n",
" g_right = grid[2:M - 1, 3:N]\n",
" g_up = grid[1:M-2, 2:N-1]\n",
" g_down = grid[3:M, 2:N-1]\n",
" # Jacobi iteration\n",
" return 0.25 * (g_up + g_down + g_left + g_right)\n",
"end\n",
"\n",
" using MPI\n",
" comm=MPI.COMM_WORLD\n",
" id = MPI.Comm_rank(comm) + 1\n",
" \n",
"\n",
" M = 50\n",
" N = M + 2\n",
" ϵ = 10^-5 # Stopping threshold\n",
" nrows = 3 # Number of grid rows\n",
" ncols = 3 # Number of grid columns\n",
" n_procs = nrows * ncols\n",
" @assert n_procs == MPI.Comm_size(comm)\n",
" max_diffs = ones(n_procs) # Differences between iterations\n",
" max_diff_buf = MPI.UBuffer(max_diffs,1) # Buffer to store maximum differences\n",
" \n",
"\n",
" # initialize grid\n",
" if id == 1\n",
" grid_a = zeros(M,N)\n",
" grid_a[1,:] .= 1\n",
" grid_b = zeros(M,N)\n",
" grid_b[1,:] .= 1\n",
" else\n",
" grid_a = nothing\n",
" grid_b = nothing\n",
" end\n",
" \n",
" # Broadcast matrix to other processors\n",
" grid_a = MPI.bcast(grid_a, 0, comm)\n",
" grid_b = MPI.bcast(grid_b, 0, comm)\n",
"\n",
" # Determine if processor is in top or bottom row of grid\n",
" top_pos = id <= ncols\n",
" bottom_pos = id > ((nrows-1) * ncols)\n",
" local grid_a_old = false # Grid a is the source grid for the first update\n",
"\n",
" # Get local partition\n",
" ind_rows, ind_cols = calculate_partition(id, M, nrows, ncols)\n",
" println(\"Proc $(id) gets rows $(ind_rows) and columns $(ind_cols)\")\n",
"\n",
" # Determine neighbors\n",
" n_left = id - 1\n",
" n_right = id + 1\n",
" n_down = id + ncols\n",
" n_up = id - ncols\n",
" if mod(id, ncols) == 1\n",
" # Left neighbor is last in row\n",
" n_left = id + ncols - 1\n",
" end \n",
" if mod(id, ncols) == 0\n",
" # Right neighbor is first in row\n",
" n_right = id - ncols + 1\n",
" end\n",
"\n",
" #println(\"Proc $(id) has neighbors left $(n_left) and right $(n_right) and up $(n_up) and down $(n_down)\")\n",
"\n",
" local finished = false\n",
"\n",
" #Perform SOR\n",
" while !finished\n",
" # Flip old and new grid\n",
" grid_a_old = !grid_a_old \n",
" \n",
" # Determine which grid is updated\n",
" if grid_a_old \n",
" old_grid = grid_a\n",
" new_grid = grid_b\n",
" else\n",
" old_grid = grid_b\n",
" new_grid = grid_a\n",
" end \n",
" \n",
" # send left and right columns\n",
" left_ind = first(ind_cols)\n",
" right_ind = last(ind_cols)\n",
" left_col = old_grid[ind_rows, left_ind]\n",
" right_col = old_grid[ind_rows, right_ind]\n",
" slreq = MPI.Isend(left_col, comm; dest=n_left-1)\n",
" srreq = MPI.Isend(right_col, comm; dest=n_right-1)\n",
" \n",
" # Send bottom row if not bottom\n",
" bottom_ind = last(ind_rows)\n",
" if !bottom_pos\n",
" bottom_col = old_grid[bottom_ind, ind_cols]\n",
" sbreq = MPI.Isend(bottom_col, comm; dest=n_down -1)\n",
" end\n",
"\n",
" # Send top row if not at the top \n",
" top_ind = first(ind_rows)\n",
" if !top_pos\n",
" top_row = old_grid[top_ind, ind_cols]\n",
" streq = MPI.Isend(top_row, comm; dest = n_up -1)\n",
" end\n",
"\n",
" # Receive left and right column\n",
" left_buf = Array{Float64,1}(undef, length(ind_rows))\n",
" right_buf = Array{Float64,1}(undef, length(ind_rows))\n",
" rlreq = MPI.Irecv!(left_buf, comm; source=n_left -1)\n",
" rrreq = MPI.Irecv!(right_buf, comm; source=n_right -1)\n",
"\n",
" # Receive top row if not at the top\n",
" if !top_pos\n",
" top_buf = Array{Float64,1}(undef, length(ind_cols))\n",
" rtreq = MPI.Irecv!(top_buf, comm; source=n_up -1)\n",
" end\n",
"\n",
" # Receive bottom row if not at the bottom \n",
" if !bottom_pos\n",
" bottom_buf = Array{Float64,1}(undef, length(ind_cols))\n",
" rbreq = MPI.Irecv!(bottom_buf, comm; source=n_down -1)\n",
" end\n",
"\n",
" # Wait for results\n",
" statlr = MPI.Waitall([rlreq, rrreq], MPI.Status)\n",
" old_grid[ind_rows, left_ind - 1] = left_buf\n",
" old_grid[ind_rows, right_ind + 1] = right_buf\n",
" #println(\"Proc $(id) received left $(old_grid[ind_rows, left_ind - 1]) and right $(old_grid[ind_rows, right_ind + 1])\")\n",
" \n",
" if !top_pos\n",
" statt = MPI.Wait(rtreq)\n",
" old_grid[top_ind - 1, ind_cols] = top_buf\n",
" #println(\"Proc $(id) received top $(old_grid[top_ind - 1, ind_cols])\")\n",
" end\n",
"\n",
" if !bottom_pos\n",
" statb = MPI.Wait(rbreq)\n",
" old_grid[bottom_ind + 1, ind_cols] = bottom_buf\n",
" #println(\"Proc $(id) received bottom $(old_grid[bottom_ind + 1, ind_cols])\")\n",
" end \n",
"\n",
" # Get local subgrid\n",
" if !top_pos & !bottom_pos\n",
" local_with_ghosts = old_grid[top_ind - 1 : bottom_ind + 1, left_ind - 1 : right_ind + 1]\n",
" lb_row = top_ind\n",
" ub_row = bottom_ind\n",
" elseif top_pos \n",
" local_with_ghosts = old_grid[top_ind : bottom_ind + 1, left_ind - 1 : right_ind + 1]\n",
" lb_row = top_ind + 1\n",
" ub_row = bottom_ind\n",
" elseif bottom_pos\n",
" local_with_ghosts = old_grid[top_ind - 1 : bottom_ind, left_ind - 1 : right_ind + 1]\n",
" lb_row = top_ind\n",
" ub_row = bottom_ind - 1\n",
" end\n",
"\n",
" # Perform one step of Jacobi iteration\n",
" new_grid[lb_row : ub_row, left_ind : right_ind] = update_grid(local_with_ghosts)\n",
"\n",
" # Calculate max difference\n",
" diffs = abs.(new_grid[lb_row : ub_row, left_ind : right_ind] - old_grid[lb_row : ub_row, left_ind : right_ind])\n",
" maxdiff = maximum(diffs)\n",
" \n",
" # Gather maxdiffs in processor 1 \n",
" MPI.Gather!(maxdiff, max_diff_buf, comm; root=0)\n",
"\n",
" # First processor determines if threshold is exeeded\n",
" if id == 1\n",
" if all(max_diffs .< ϵ)\n",
" finished = true\n",
" println(\"THRESHOLD SUBCEEDED - TERMINATE SOR\") \n",
" end\n",
" end\n",
" \n",
" finished = MPI.bcast(finished, 0, comm)\n",
"\n",
" if finished\n",
" # Set ghost cells to zero again so MPI.Reduce gives correct output\n",
" new_grid[ind_rows, left_ind - 1] .= 0.0\n",
" new_grid[ind_rows, right_ind + 1] .= 0.0\n",
" if !bottom_pos\n",
" new_grid[bottom_ind + 1, ind_cols] .= 0.0\n",
" end\n",
" if !top_pos\n",
" new_grid[top_ind - 1, ind_cols] .= 0.0\n",
" end\n",
" end\n",
" end\n",
"\n",
" using DelimitedFiles\n",
"\n",
" # Reduce matrix & store result\n",
" if !grid_a_old\n",
" sor_result = grid_a\n",
" else \n",
" sor_result = grid_b\n",
" end\n",
"\n",
" MPI.Reduce!(sor_result, +, comm, root = 0)\n",
" sor_result[1,:] .= 1.0\n",
" if id == 1\n",
" writedlm(\"SOR_result.txt\", sor_result)\n",
" end \n",
"\n",
" MPI.Finalize()\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "9465e78a",
"metadata": {},
"outputs": [],
"source": [
"rmprocs(manager);"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "d0bbbb61",
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3deVyVZf74/+s+57CDAgIqIC6IhuZCgKS4fVTEccPdj2upWdo0ad+szNGppiadcvJTmVPu1KgpmVZOi0tKxaS5Z9hgKu5gpYCEyHZfvz/uOj9Cw3Ofbs4Bzuv5OH/c5+I657oePKR372tVpJQCAABXZXJ2BwAAcCYCIQDApREIAQAujUAIAHBpBEIAgEuzOLsDAAD8pvz8/AMHDpw+fToxMbF9+/a3rPPNN9+8/fbbFotl0qRJkZGRepsgIwQA1F79+/d/4oknFixY8Omnn96ywpEjRxITE93d3UtKSuLj47Ozs/U2QUYIAKi9vvzyS7PZnJSU9FsV/vGPfzzwwAN/+ctfhBCXL19+9dVXX3rpJV1NkBECAGovs9lcfYX09PT+/ftrz8nJyXv27NHbhDGBsKioyJDvAQDAdlLKy5cvh4SEaG9DQkJycnL0fomtQ6MLFix47bXXrG8vXrzo5eUlhEhPT7/nnnsKCwsDAwPXrVvXpUuXW35cVdU//vGPDzwwQ2//6j0ppZTSZCI1v43y8nKLhZH825BSqqp62/+DBv+cbFTlF9W5cycbPnRNigpdrWRlnfzLgherFCYlJU2fPv22n1UUxWw2l5eXa2/Ly8vd3Nx0tS5sD4TFxcVTpkyZP3++9tbT01NrcsKECS+88ML48eOXL18+ceLErKwsRVFu/nhJScnKlatff32l3v4BAGoJKctuW0eVX0iRp+trjx8/eObMmccee6xy4Z133mnjx0NDQy9dunTXXXcJIS5evBgaGqqrdaFrsYyXl1dAQEDlkp07d5pMpnHjxgkhpk2bNm/evIyMjO7du+vtBACgflBlhRTluj4iRUVYWNjo0aNt/8iVK1dOnDjRtWtXIcTgwYM3b948ePBgIcTmzZuHDBmiq3Wha47wtdde8/Pzu/POO9euXauVnDx5Mjo6WksBzWZzmzZtTp06pbcHAIB6Q4oKKfW+1Gq+cOHChUlJSYcOHXrttdeSkpI+++wzIURGRkZKSopW4f/9v/+3Y8eOUaNGDRo06MSJEzNm6J6DszUjnDRp0qxZs4KDg3fs2DF+/PgmTZoMGDCgoKDAx8fHWsfPzy8v79YZsZSyvLycRaoAUHeVlZXZMQP3O/3v//6vdVGoEKJVq1ZCiB49enz44YdaSURERGZm5o4dOywWS1JSUuWoZCNbA2GnTj/PkQ4ZMmTy5MnvvffegAEDgoODr127Zq2Tl5dnXbpThaIoFoulvLy6sA8AqM1siYKqrJBS59CorG5xTcuWLVu2bFmlMCAgIC4uzvq2YcOGo0aN0tVoZfasmyopKdGWjLZr1+7o0aMVFRVms/nGjRvffvttu3bt7O4KAKDOsyMQ6lxlajhbxypfeeWVI0eOZGdnr1q1av369WPGjBFCJCYmNm7c+JlnnsnJyZk/f3779u07d+5ck70FAMBgtmaEWVlZK1euvHHjRosWLTZv3tyjRw8hhKIoW7dunT17dvfu3Tt16rRp06ZqvkFRLOyWA4D6TarlujNC1cmzZrYGwsq76Str3br1tm3bjOsPAKAOk3VwaJSzFQAAhpFCdyAU1S6WcQDGKgEALo2MEABgGKmWS1XvHCFDowCAekNWCL1Do8LJi2UYGgUAuDTHZYQmxV2YPBzWHADACYw+WcYBGBoFABhHVgidc4ROXzVKIAQAGEctF+rtry389UfYPgEAgPOQEQIAjCPL9a8arRDCXCOdsQ2BEABgHDvmCFUnB0KGRgEALs2B2ydM7mZF98XBAIC6RGXVKADAldlxsoysI9cwAQBwe7Jc0ZkRKs7OCJkjBAC4NDJCAIBxVFX3Bvm6ckM9AAC3pagVDI0CAFCXOC4jtCgeZsXXYc0BAJxAVugeGmXVKACg/rBjH6GzD90mEAIADKOoFYrRGeGlS5dWrVqVl5c3bNiwnj173lwhPz9/zZo1586dS0hIGDt2rKIoutpnjhAAUHvl5+d36dIlJyenRYsWI0aM+OCDD6pUKCwsjIuLy8zM7NChw0svvfTkk0/qbYKMEABgIDtOlqkug3zzzTfbtGmzbNkyIYS3t/fChQuHDBlSucKGDRv8/f1XrlwphOjbt2+7du2eeOKJgIAA29snIwQAGEZRVW10VMer2qHRzz77LCkpSXvu16/f3r17S0tLK1fQkkXtOSIioqys7KuvvtLVZwIhAKD2ysnJCQ4O1p5DQkKklLm5uZUrdOzYcd++fUVFRUKIPXv2lJWVXbp0SVcTjhsaNSvuHmyfAID6TbVn+0RGRoY17dMMGzbsj3/8oxDCzc2tvPznsVbtwd3dvXLNlJSUjRs3dujQoWPHjufPn4+IiPD29tbVPnOEAADDKFLVu2pUUdWoqKgnnniicmHLli21h7CwsIsXL2rPFy5ccHNzsyaIGpPJ9Pbbb588efLq1avt27cPCwuLiorS1QECIQDAOHZlhCEhIf369bvlD4cOHfrss88uWLDA3d1906ZNgwYNMpvNQoj9+/c3adKkWbNmWrXWrVsLIV566aXw8PCYmBhd7RMIAQC11/Dhw5ctW9ajR49WrVrt2rVr+/btWvlDDz00bty42bNnCyHuuOOOmJiYixcvZmdnf/jhh3r3ERIIAQDGUe0ZGq3mp+7u7rt27dqzZ09BQcGrr74aFBSklaempgYGBmrPW7duPXbsmI+PT+/evfVOEAoCIQDASFLVfa2SlNX/3GKx3Dxwescdd1R+rvxWL7ZPAABcGhkhAMAwih0X87rO7RMW4cE+QgCo5+y4hokb6gEA9UcdzAiZIwQAuDQyQgCAYRRV96rR6g/ddgACIQDAOFL/0Kiz5wgZGgUAuDQyQgCAcfQPjd52Q31Nc1wgdBPuXoLtEwBQnymq1D1HyNAoAABOxNAoAMA4dmyoZ9UoAKD+0D806kJzhACA+s+O2ydUJwdC5ggBAC6NjBAAYBhFZY7wt7lJNy/Vy2HNAQCcwJ45QrZPAADgPAyNAgCMUwcXyxAIAQDGqYNDowRCAIBx7Dhr1NkZIXOEAACXRkYIADCQqn+o02XmCC3C4ik9HNYcAMAJ7Jgj5PYJAACciKFRAIBxauBi3vPnzy9fvrygoGDYsGF9+vS5ucIPP/ywZs2ac+fOtWjRYurUqYGBgbraJyMEABhH20eo61VtIMzLy0tISLh27Vp0dPTYsWO3bt1apUJRUVFCQsJ///vfxMTEw4cPd+3ataSkRFeXyQgBAMZR9W+HqLZ+ampqu3btXn75ZSGEh4fHokWLhg0bVrnC0aNH8/LyVq1apSjKmDFj/Pz8jh8/HhMTY3v7ZIQAgNrriy++6Nevn/bct2/fr776qkrCFxkZqarq8ePHhRBHjhzx8PBo0aKFriYMyAjLy8sLCgr8/f3NZvPv/zYAQB1m9BxhTk5OUFCQ9hwSEiKlzM3Nbd68ubVC48aN161b161bt8DAwLy8vLS0tICAAF3t6wuE+fn5CQkJTZo0SU9P10q2b99+zz33WCwWIcT69et79OjxW591EyYvxU1XcwCAOkZK3UOjUmZkZCQlJVUuGzJkyMMPPyyEcHd3Ly8v1wpLS0uFEB4ev9qJd/bs2WnTpr322ms9e/bcuXPn5MmTDx061LRpU9vb1xcI58yZ07Jly8uXL1v7NHny5GXLlo0YMeLNN9+cPHnyqVOnTCaGWwEAOkRFRT3xxBOVSyIjI7WHsLCwCxcuaM8XLlxwd3cPDg6uXPO9995r167dxIkThRBTp05dsWLFtm3bpk+fbnvrOoLWp59+mp2dPWnSJGvJjh07vLy8RowYIYSYOHFiUVHRF198YfsXAgDqG21Dvc5XSEhIv19r2bKl9n0pKSlbtmzR5gU3btw4ePBgbRruyy+/PHv2rBAiODj4zJkzWrJYXFx87ty5xo0b6+qyrYGwqKho1qxZy5cvVxTFWnj69Om2bdv+/EUmU1RU1OnTp3U1DwCoV7ShUV2vakdShw0b1qRJk27duo0aNeqNN954+umntfLZs2dv2bJFCDFixIiIiIjY2Nhp06bFxsbeeeedAwcO1NVlW4dGH3/88XvvvTcyMnLfvn3WwmvXrnl7e1vf+vr65ufn3/LjUsqKigo2awBA3VVeXq6tCKmOKoTeE9Oqre/m5rZ9+/YvvviioKBgxYoV1oUw69at8/f3F0J4eHjs2bPnwIEDFy5cePDBB2NjY3U2b1toyszM3Lhx49KlS9PS0vbt25efn5+WljZ8+PDg4ODKkS8vLy8kJOSW36AoCmtKAaBOu30UrBlms7lXr15VClu3bm19VhQlPj4+Pj7evu+3dWi0T58+7777buVAWFFR0aFDh6NHj2rreYqLizMzMzt27GhfPwAA9YGUul914vaJ9u3bb9q0SXtev379iy++qL3t2rVrRETE3LlzZ8yYsWTJEm1wtgY7CwCo5YweGnUA3XluRERE5d0eW7dufeyxx4YPH96pU6eNGzdW80E3xeSlMEkIAKhddEem7t27d+/e3fq2efPm1mQRAODqpP4Mz8kjo6zjBAAYR9a9C+oJhAAAA9mRETp7jpDj0AAALo2MEABgHDtWjTI0CgCoP1RFqMrtq1UmddY3muMCocVk8uRiCgBALUNGCAAwjJSK1JnhVXsvryMQCAEAxqmDq0YJhAAA46iizs0RMmkHAHBpZIQAAMNIaZKqvhRLOnuSkEAIADCOHdsnXGeO0KIIT7OTB4IBAKiCjBAAYByp6F/84jIb6gEA9Z5UhdQ9NMqqUQAAnIeMEABgHGkSOleNOn0fIYEQAGAcVdE7NKr3SDbDEQgBAMaxY7GM65w16mYSnmaHtQYAgE3ICAEAhpFS0X2yjN45RaMRCAEAxrHnYt7b/PzMmTP//Oc/CwsLU1JSkpOTq/z0u+++e/fddyuXjB8/vlmzZra3z/YJAEDtdfXq1bvvvruioiImJmby5MmbN2+uUqG0tDTvF0eOHPnLX/7i6empqwkyQgCAYaQ0eNXo2rVrO3XqtHjxYiGEm5vbCy+8MHLkyMoV2rdvv2jRIu15zpw5Q4cODQ4O1tUBMkIAgGGkNOl9Vb/KNCMjo0+fPtpznz599u/fX1JScsua5eXl69atmzp1qt4+EwgBAMZRf5kmtP1VbSDMyckJCgrSnoODg6WUubm5t6z5wQcfmM3m/v376+2yg2+fcFhrAIA6IyMjIykpqXLJoEGDZs+eLYTw8PAoLS3VCrWH35oCXL169ZQpU8xm3ZGGOUIAgGGkHSfLqEpUVNQTTzxRuTAqKkp7CA8Pv3DhgvZ8/vx5d3f3W04B5ubmbt++fcmSJXb0mUAIADCQSUi9k25KSEhIv379bvmzYcOGzZ8/f8GCBZ6enhs2bBg6dKjJZBJCfP7552FhYa1atdKqrVmzJjExsXXr1nb1GACA2iolJSUiIiIhISElJWX16tVPP/20Vj5nzpz333/fWu3NN9+0Y5mMhowQAGAcu4ZGq0nKLBbLRx99tHfv3vz8/DfffLNhw4Za+caNGxs0aKA9V1RU/Pvf/9a1if5XTdj3MQAAbmbPxby3O6TbZDJ169atSmGLFi2sz2az2TpGagcCIQDAQHXv9gnmCAEALs2B+wi5hgkA6jupmnTfPqF7lanBGBoFABhGSkXvjfOSoVEAAJyIjBAAYBg7TpbRvbjGaARCAIBxpCL03jivN3AajUAIADCOHXOEwsmBkDlCAIBLc1xG6GYSnmZnrw0CANQk+26fqKHO2IihUQCAYaTQPTTqdAyNAgBcGhkhAMAwUlV0nyyjd5Wp0QiEAADD2HOyjLNXjRIIAQCGsWdDvVozXbEZc4QAAJfmwNsnFOnF9gkAqOcYGgUAuDJV0X1kmrP3ETI0CgBwaWSEAADDSKnovWjX6RvwCYQAAMOwfQIA4NKkZPsEAAB1iiNvn5CeZmfHfQBATWJoFADg0qQUqt5AyGIZAACqcfLkyddffz0vL2/48OGDBw++ZZ3PP/9806ZNRUVFcXFxDz74oK7vZ44QAGAYKU1S1fmqNiO8cuVKt27d3NzcevbsOX369E2bNt1cZ9WqVSNHjgwPD+/du/eJEyf09pmMEABgGHvmCKutv2bNmtjY2IULFwohFEV54YUXxowZU7lCfn7+7Nmzd+zYcffddwshJk+erLfPZIQAgNrryy+/7N27t/bcu3fvQ4cOlZSUVKnQpEkTi8XyzDPPvPrqqwUFBXqb0BEIi4uLz5w5U1RUVKW8tLQ0Ozu7rKxMb9sAgHpGywh1vUS11zHk5OQEBQVpz8HBwVLK3NzcyhXOnDlTUFDwyCOPhISE/Oc//4mPj79+/bquPts6NDpr1qw1a9YEBwdfvnx50KBBqampnp6eQogPPvhg6tSpjRo1ys/PX7duXd++fX+zJUX1NFfo6hwAoG6RqtC7oV5KZc+ePXFxcZULx44d+9hjjwkhPD09S0tLtUItF9Sij5W7u/uPP/549OjRpk2bzpgxo3379ps3b540aZLtHbA1ED766KMvvfSS2WwuKCjo1q3bihUr/vSnP5WUlEydOjU1NXXgwIGbNm2aMmVKdna22Wy2vXkAQL1i1xxh586dFy9eXLkwNDRUewgPDz9//rz2fP78eQ8Pj+Dg4Mo1w8PDvb29mzZtKoRQFCUyMjInJ0dXB2wdGo2IiNAiXMOGDaOioq5duyaE2L59e4MGDQYOHCiEGDVqVGlp6WeffaareQAA/P39Y39NC2xCiOHDh2/evLm4uFgIsW7dupSUFJPJJITYvXv3yZMnhRC9evXy9vbev3+/EKKwsHD//v2dOnXS1bqOVaMHDhzYtWvXf//734KCgvvvv18IkZ2d3aZNG+2nJpMpMjLyzJkzupoHANQnUv/FvNUbMmTIypUr4+PjIyIijh49unPnTq187ty548aNmz17tqen55IlSwYPHtyrV6+DBw8OGjQoOTlZVxM6AmFJSUleXt6VK1eKioqKioqCg4MLCwu9vLysFXx8fLRM8WZSyooKJggBoA4rKytzc3Orvo6UJv3XMFVX32KxbNu27cCBA/n5+V27dvX19dXK33nnHT8/P+15woQJ//M//3Ps2LHnnnvOmp7ZTkcgTExMTExMFELMnDnzqaeeSk1NDQkJycvLs1bIy8tr3LjxLT+rKApzhwBQp902CtYQRVHi4+OrFDZr1qzy29DQUOu0ol727CMMCwvTMr9OnTodPnxYW89TVFSUmZmpd2QWAFCfqKqi96X72iaj2ZoRLliwoGfPnv7+/seOHVuyZMmyZcuEEF26dGnTps0jjzzywAMP/N///V/Xrl2jo6NrsrcAgNpN6j5Eu9pthI5gayD08vJavHhxQUFBeHh4amqq9djTrVu3zp8//7777ouJidm4cWN1LZkk+wgBoH6zY7FMnbl9Yt68efPmzbu5PDQ0dPXq1YZ2CQAAx+HQbQCAYQw/dNsBCIQAAMOoUtF9MW8NdcVm3D4BAHBpZIQAAMMwNAoAcG36t08I1wmEbiauYQKAes6eI9aEkwMhc4QAAJfG0CgAwDCqVFS9F/PWlSPWAAC4LXsWy9RQV2zG0CgAwKWREQIAjMP2CQCAK1Ol0H2yjOsEQouieprLHdYcAMDx7Ll9ooa6YjPmCAEALo2hUQCAYaRQ9G+Qd5mhUQBAvSf13z6ht77hGBoFALg0MkIAgGHs2FDvQoduAwDqPTu2T6jOXjbqwO0TJtXTwvYJAKjP7DpijTlCAACch6FRAIBh7Fg1etsM8sSJE0uXLs3Pzx8+fPjw4cNv+rh88sknrW+7du2akpKiqwNkhAAAw2gny+h8VfeFP/zwQ2Jior+//4ABAx566KENGzbcXOfvf/+7p6dnQEBAQECAt7e33j6TEQIAaq81a9bEx8f/9a9/FUKoqvriiy+OGzfu5mqzZs0KCAiwrwkCIQDAMPZsn6h2sczevXt79+6tPffq1Wvy5Mk3btzw9PSsUu355593c3Pr0aPHH/7wB52tMzQKADCOFIqq/1XNF+bm5jZq1Eh7DgoKklLm5uZWqTNu3Ljw8HBPT8/77rvv0Ucf1dtnB26fMLN9AgDqOSn1X6skxZ49e+Li4iqXjRo1au7cuUIIT0/P0tJSrbCkpEQIUWUWUFGU9evXa88pKSkxMTF//vOfAwMDbW+foVEAgJN17tx58eLFlUvCwsK0h/Dw8HPnzmnP586d8/T0DAoK+q3vufPOO81mc05ODoEQAOAcql2Hbvv7+8fGxt7ypyNHjnzssccWLFjg7e391ltvDR8+3GQyCSF27NjRvHnzNm3a5OXl+fn5WSwWIcS6det8fHwiIyN1dYBACAAwjJSi+u0Qeg0ePHjNmjV33XVXeHh4VlbWrl27tPL58+ePGzeuTZs2H3744Zw5c9q1a3ft2rXs7Oy1a9fevJSmegRCAEDtZTabt27deuTIkYKCgi5dunh5eWnlW7du9fHxEUJMmDAhMTExOzvb19c3Ojra19dXbxMEQgCAYaSokfsIO3fuXKWkadOm1ucWLVq0aNFCV6OVEQgBAIax64Z6J2MfIQDApTkuI3QzVXhayhzWHADA8exbNVpDnbERQ6MAAMPYsaHe6UOpBEIAgGHk7Y5Mu9VHnIw5QgCASyMjBAAYxo4N9cZuwLcDgRAAYBg7Fsvov7bJYAyNAgBcmgOvYTJVeLiVOqw5AIDjSaH7Yl5WjQIA6g+7Vo0SCAEA9YUqdW+Qd/piGeYIAQAujYwQAGAYtk8AAFxaXZwjZGgUAODSHJcRms0VHm7cPgEA9ZmUbJ8AALgwVQhV/0eci0AIADCO/ozQ6ddPMEcIAHBpZIQAAMNwQz0AwKXZMUfo7JFRhkYBAK7NgbdPmCvc3bl9AgDqM3tun2BoFABQb0iGRgEArkzKn/fU63jd7juPHz8+c+bM8ePHb9q0qZpqWVlZc+fOPXz4sN4+EwgBALXX999/36NHj9DQ0GHDhj366KP/+te/blmtoqJi2rRpK1asyMzM1NsEQ6MAAMPYsX2i+jnC1atXd+3adcGCBUKI8vLyF198ceLEiTdXW7JkSbdu3YqLi3U1rSEjBAAYSep/VeOrr77q1auX9tyzZ8+jR4/euHGjSp3s7Oy1a9c+9dRT9nVYRyAsLCzMysq6du1alfKffvrp22+/vX79un09AADgt+Tk5AQGBmrPQUFBUsrc3NzKFVRVnTJlypIlS3x8fOxrwtah0WHDhu3atSsiIuLChQv33HPPyy+/rCiKECItLW3GjBnNmzc/f/782rVrBw0a9FvfYLZUuHuwfQIA6jP7hkb37NkTFxdXuXDkyJFPPvmkEMLb27ukpEQr1HJBb2/vyjX/+c9/tmzZMikpye4+2xoIJ06c+Pbbb3t6el64cCEmJiYpKWnIkCHFxcUzZszYtGlT375933vvvfvvv//s2bMWC/OOAOCi7DtZpnPnzosXL65c2KxZM+0hPDz83Llz2vPZs2c9PT2DgoIq1/zoo4/S09Pfe+89IcS1a9dmzJhx5MiRKt9WPVuD1qhRo6x9uuOOO86fPy+E+OSTTxo1atS3b18hxNChQ2fOnJmenq69BQC4ILvuIxT+/v6xsbG3/OmoUaMeeeSRBQsW+Pj4pKamjhw50mQyCSE++uijFi1aREdHv/3222VlP19227t37wcffHDSpEm6OqA7e8vMzDx69Ojq1auFEGfPnm3durVWrihKq1atzp49q/cLAQD4LQMHDkxNTe3UqVN4eHh2dvauXbu08qeffnrcuHHR0dG+vr7WyhaLxcfHp8rY6W3pC4Q//PDDqFGjnnnmmaioKCFEUVGRh4eH9afe3t6FhYW3/KCUsqKiQldbAIBapby8/LaTX3acLFN9fbPZ/M4772RmZubn58fFxVmDzrZt224OeO+//37Dhg11tq8nEObl5fXv33/kyJGPPPKIVhISEpKXl2etcPXq1SZNmtzys4qimM1mvZ0DANQetiwBseOsUWFD/fbt21cpCQ4OvrlaWFiYvqaFELZvnygoKBgwYEDv3r2fe+45a2FMTMyhQ4e09TyFhYWZmZmdO3e2oxMAADiLrRnhoEGD8vLyoqOjly9fLoSIiYmJj4+PjY3t2LHjzJkz77///ldeeaV3795t27atyd4CAGo1VQpV5ynaeusbztZAGBMTU1paevDgQe1tSEhIfHy8EGLLli1PPfXUvHnzOnTo8Nprr1XzDWZLubtHye/sLgCgNquLt0/YGghfffXVW5YHBwcvW7bMuP4AAOow7fYJfR8RTr6PkLNGAQAujVNgAACGse9kGeciEAIADCOlkDojm976hmNoFADg0sgIAQCGkUJRdS5+cXZC6MBAaLZUuHENEwDUa3bsI2RoFAAAZ2JoFABgGKl/qNPZCSGBEABgHKn/hnq9c4qGIxACAAwj7Zjzc3ZKyBwhAMClkRECAAxjx8kyeusbznGB0GQpd/Nk+wQA1GuyRi7mrVEMjQIAXBpDowAAwzA0CgBwafacLFMzPbEdgRAAYJi6uKGeOUIAgEsjIwQAGEbadei2c5eNOnD7hFuFxavEYc0BABzPvot5nRsIGRoFALg0hkYBAIZRa+Bi3mPHjr388ssFBQXDhg2bMGFClZ9ev359yZIlx48fLy0t7dix40MPPRQQEKCrA2SEAADDyF92UOh6VePy5cu9evVq3br1hAkT5s2bl5qaWqVCUVFRYWHhsGHDxo8fn5GRMWjQIL19JiMEABjGju0T1Vu1alX37t3nzp0rhLhx48bChQvvueeeyhWCg4MXLVqkPcfGxjZv3rygoKBhw4a2N0FGCACovfbv39+jRw/tuXv37l9//fWNGzduWfPatWtvvfVWXFxcgwYNdDVBRggAMIwdJ8tUf8Rabm5uo0aNtOegoCCtpEWLFlWqdejQISsrKyAg4MMPP1QUfZOUjtw+UW5m+wQA1Gv2DY3u2bMnLi6ucsmIESPmzcwLQuMAABK2SURBVJsnhPD29ramgMXFxUIIHx+fm7/h2LFjFRUVqampycnJJ06cCAwMtL11MkIAgJN17tx58eLFlUusOV+zZs3Onj2rPZ89e9bLy8uaIFZhNpunTp365z//+fDhw3379rW9dQIhAMAwdpwso0rh7+8fGxt7y5+OHj36T3/60/z58/38/NasWTNq1CiTySSE+OCDD1q1atW+ffvvv/++YcOGHh4eQoh9+/ZdvXq1TZs2ujpAIAQAGMae2yeqrf+HP/whISGhY8eOoaGhOTk5O3fu1Mqfe+65cePGtW/ffs+ePTNnzmzbtm1JSUl2dvbSpUubNWumqwMEQgCAYQy/fcJkMm3YsOG77767evVqTEyMu7u7Vv7xxx97enoKIcaMGdO/f/+srCwvL6/IyMhbziBWj0AIAKjtoqKiqpRUPj7G398/ISHB7i8nEAIADGPP7RM10xPbEQgBAIapixfzOnAfoaXc7OWw1gAAsAkZIQDAMIavGnUAAiEAwDBSf2AjEAIA6g95u7NDb/kR5+L2CQCASyMjBAAYxvDbJxyAQAgAMIwUUjp/sFMfxwVCxa3c5FXhsOYAALAFGSEAwDD23T7hXARCAIBhpNS/HYJACACoN1T9i1+cvliG7RMAAJdGRggAMAy3TwAAXJod2yecvt3Cgdsn3FWTt9MDPwAAv0JGCAAwjKp/OwSHbgMA6o+6eDEvq0YBAC6NjBAAYBgppapzrFM6e2yUQAgAMExd3FBPIAQAGEZK6fQMTy8HBkI3RfF2XGsAANiCjBAAYBjJ0CgAwJWpUhi+WObw4cMvv/xyXl7e8OHD77333io//fHHH9etW7d3796ysrLu3bvPnDnTw8NDVwcIhACA2isnJ6dPnz7z58+/4447HnroIVVVp06dWrnCp59+un///qFDh3p4eDz77LOHDx9OTU3V1QSBEABgGMM31K9evbpXr16PPvqoEOKnn3567rnnqgTCMWPGjBkzRntu3Lhx//79165dqyiK7R0gEAIADCOFVA09dPvAgQPdu3fXnrt37/7NN98UFxd7eXndsnJ2dnZYWJiuKCj0nixTVFT0008/VSm8evXqwYMHCwoKdH0VAKD+UaXU+6p+ijA3NzcwMFB7btSokVbyWzUff/zxv/3tb3r7bGtGmJqa+txzz508eXLQoEHbtm2zlr/11luzZs2Kjo7Oyspavnz5iBEjfvMr3CzC26y3fwCAem/Hjh2RkZGVSyZNmvT0008LIXx8fG7cuKEVFhcXCyF8fX1v/oYrV67079//gQceGD16tN7WbQ2EHTt23LBhw44dOzIyMqyFRUVFDz/88Pvvv9+jR4+PPvpo6tSpgwcPdnd319sJAED9YN99hN26dXv99dcrF2rJnxCiWbNmZ86c0Z6zs7O9vb2DgoKqfEN+fn5ycvLAgQOfeuopO/ps69BoTExMXFxclSD3ySefNG7cuEePHkKIAQMGmM3m9PR0OzoBAKgf5C+nrNn+kkL4+Pi0+rWGDRtqXzhmzJi0tLRr164JIVavXj169GhtCnDLli3Hjh0TQly7di05Obl79+6LFi2yr8+/a7HMuXPnWrVqpT0ritKyZctz5879VuU6d+gOAMDpkpOTe/bs2a5du6ZNm+bn5+/YsUMrX7Ro0bhx4zp06PCvf/3rq6++OnHixJtvvqn96NSpUwEBAbY38bsC4fXr1yvvW/Ty8rp5KY1GSqmqTj89AABgv7KyMjc3t+rrqPpXjVZf32Qypaamnjlz5tq1a+3btzebf15rsnPnTm2Qctq0aePGjav8EX9/f10d+F2BsHHjxlevXrW+vXLlStOmTW9ZU1EUa+8BAHXRbaOg+Hkfod45wttr0aJFlRI/Pz/twcPDQ+9RMlX8rot54+LiDh06pC3jKSgoyMzMvOuuu37PFwIA6jRVqHpf0tmnjdoaCL/77rvly5fv3bv33Llzy5cv1xbFdOrUKT4+ftq0aZ9++umUKVOSk5Nbt25dk70FAMBgtg6NarvmAwMDu3btevDgQU9PT6188+bNzz///D/+8Y9OnTo9+eST1X2Fu5vw9vyd3QUA1GZ2nSzjZLYGwoSEhISEhJvLAwICXnzxRUO7BACoqwxfLOMAv2uOEACAuo5DtwEAhrHjZBmnD44SCAEAhpFCqjpXgTp9aJRACAAwjKpIVdG5WEZnfcMxRwgAcGkOzAjd3aWPj+OaAwA4nBSq3qFR/XOKBmNoFABgGCmk3pNi6szJMgAA1EtkhAAAw6gMjQIAXJmqSFXRuX3C2atGCYQAAMPYtViGOUIAAJzHcRmhdPNUvRs4rDkAgONJ/fcLMkcIAKg/VP1HrDk9EDI0CgBwaWSEAADD2LFYRm99wxEIAQCGURVV7/YJpw+NEggBAIaRQpWiQu9HaqgzNmKOEADg0hy4fcLdQ/Vh+wQA1Gd1cUM9Q6MAAMPYdfsE2ycAAPhtBw4cGD9+/KBBg1asWHHLCtu3b3/++ecfeOCBrKwsO76fQAgAMIwqK/S+qs8gL1261K9fv7vvvvvhhx9etGjRLWPhM888c+HChQ0bNly6dMmOPjM0CgAwjOFDo6tWrerTp8/DDz8shFi0aNHTTz89ffr0KnUyMjKEEFu2bNHZ2Z+REQIAaq+DBw9269ZNe+7Wrdvx48eLi4uNbYJACAAwjBQVqs5X9fsOL1++HBgYqD03atRICJGbm2tsnx03NKq6eVd4BzisOQCA49l3+8SOHTsiIyMrF06YMOGvf/2rEMLHx8eaAmoPfn5+BnX2Z8wRAgAMowpV1X+yTLdu3V5//fXKhUFBQdpDREREdna29nz69GkfHx8tLzQQQ6MAACfz8fFp9WsNGvx8AMvYsWPT0tLy8/OFECtXrhwzZoyiKEKItLS0o0ePGtI6GSEAwDBSqlLqHBqV1a0a7d+/f1JSUnR0dOPGjUtKSrZv366VL168eNy4cZ06dRJC9OzZ85tvvikoKEhJSbFYLF9++WXbtm1t7wCBEABgIDu2T1RXX1GUlStXXrx4MT8/Pzo62mT6eSAzPT3dYvk5hG3btq2i4v8fj7VmkzYiEAIAaruwsLCwsLDKJZ6entZnvZGvCgIhAMAw2mExuj6idyjVcARCAIBh6uKh2468hslb9TZ4zSsAoFaRokLqzAiFszNCtk8AAFwaQ6MAAAMZvGrUAQiEAADDqFJV9e4jdPYcIUOjAACXRkYIADCMFGr1t0nc8iM11BkbEQgBAIYx/Ig1B3BgILT4Cq8QxzUHAHACO65hYvsEAADOw9AoAMAwUkrdR6a50NAoAKC+s2eOkO0TAAA4ERkhAMAwdm2f0Hk2qdEIhAAAw9gxNOpCc4Qmi4/Zg+0TAFC/qUL3dgjmCAEAcB6GRgEAhpFC9/YJbqgHANQndpwsw9AoAADOQ0YIADCOVIXuoU6XWTUKAKj37NhHqH+VqcEcFwjNFh8PjyCHNQcAcAY7tk9w+wQAAM5DIAQAGEfKn6cJdbxuM0e4d+/e0aNHJycnL1u2rCZu8WWOEABgGCmk/u0Q1dW/cOFCcnLyokWL2rZtO2PGDJPJNGPGjN/Tw5uREQIAaq9Vq1YlJSXNnDmzT58+CxcufOWVVwxvgkAIADCQqv9VXUZ4+PDhrl27as9du3b99ttvr1+/bmyPCYQAAANJu16/KTc3NyAgQHsODAwUQly+fNnYHjtojtBkMj35xNv/fO0/jmmuDiksLCwtLW3UqJGzO1LbnT9/PjQ01Gw2O7sjtdr169d/+umnkBCuebmNS5cuBQUFubu7O7sjtVpJScnVq1ebNm1qLRk//r/PPvts9Z86cuSQ3oY++uij0aNHR0ZGVi4cP3681pafn19xcbFWqOWCfn5+epuonoMCoYeHx8mTJysqnHz7Yi1UUVEhpbRYWLV0GyUlJR4eHs7uRW0npSwrK+O/77fFPycbVflFVQ6KBhowYMA333yjqr/aTWj9/7mIiIjTp09rz6dPn/b19TU8c1BqYikqAACG2LFjx7Rp044cORIYGDh9+nQp5cqVK41tgkQEAFB79evXb+DAgdHR0UFBQSaT6eOPPza8CTJCAEBtd/ny5fz8/DZt2iiKYviXEwgBAC6N7RMAAJfGHKFDlZWVHTt27Ouvv/b19R01apS1XEq5fv36gwcPtmrV6r777vP09HRiJ2uDo0ePfvLJJ5cvX27Tps3EiRN9fHy08qKiouXLl58/fz4xMXHkyJHO7aTT5eXlbdmyJSsrS1GUrl27Dh061DpqtG/fvnfeecfPz+/ee++NiIhwbj9rD+0f1eTJk7W31r+7yMjIadOmufjfXWlp6dq1a61vO3fu3KVLF+352LFj69evt1gskyZNatOmjXP6V5PICB1qxYoVo0aNevXVV6vsxXnyySf//ve/R0VFvf/++/z3vaSkZMCAATk5OREREWlpaXfffbe2i0hKmZycvHv37tatW8+bN2/RokXO7qmTfffdd59//nmTJk0CAgIeeeSRRx99VCvfvXt3cnJykyZN8vLy4uPjv//+e+f2s5b4+uuvJ0yYYP0tCSHmzp2r/d299957lf/H1DVdv379gQceOHXq1OnTp0+fPn316lWt/MiRI4mJib6+vlLKhISEU6dOObefNULCgbRdgxs3buzYsaO1MD8/38fHJzMzU0pZVFTUoEGDw4cPO62LtYCqqjdu3NCeS0pKgoODP/nkEynl7t27mzZtWlpaKqX88ssvg4KCrNXw8ccfN2nSRHtOTk5evHix9pySkvK3v/3Nef2qLcrKyu6+++4lS5YEBQVpJdrf3fHjxyV/d1JKKfPy8sQvO5srmzhx4uOPP64933fffbNmzXJ412ocGaFDmUy3+IUfPHjQ39+/Xbt2Qghvb+/ExMT09HSHd60WURSl8h7ekpISX19fIUR6enqvXr3c3NyEEAkJCWVlZZmZmU7rZW2iqmpGRkaHDh20t+np6UlJSdpzUlKSi/9z0rz44ou9evW66667rCUHDhwICAiIjo4Wv/zdffbZZ87rYG2xdOnSpUuXHj161Fry2Wef9e/fX3uur/+cCITOl5ubGxwcbH3buHHjS5cuObE/tcq8efNiYmK0I3cr/6IURQkODuYXVVJSEhkZGRAQ8Pbbb69Zs0YIkZeXd+PGDesvKiQkJCcnx6l9dL6srKy33nprwYIFlQv5u6tCUZT+/ftfvnz52LFjPXr0WLp0qRBCSln5F1Vf/zmxWMb5LBZL5cPnysrKOPxJs3Tp0rS0tM8//1xbA3LzL4qzxNzd3Q8cOPDjjz8uXLhw7Nixn3/+uZYxl5eXaxXKy8td/Lekqur06dNfeeUV65IrDX93VTRs2PCTTz7RnocPHz5ixIj777/f3d3dYrHU+39OZITOFxoampOTI3/Z0Hnx4sUaOtCvblmxYsXixYvT09Otix7DwsIuXryoPZeVlX3//fehoaHO62CtoChKQEBAVFTUkiVLMjIyzp8/7+vr26BBA+svin9Op0+f3rdv39y5c+Pi4u6///78/Py4uLgzZ86EhoZeunSJv7tb6tatW3FxcW5urhAiNDS08j+nevlHRyB0vi5dupjN5t27dwshLl68uG/fvkGDBjm7U06Wmpr67LPP7tq1q0WLFtbCIUOG7N69+8qVK0KIf//736GhodrEqsuqfCvbgQMHPDw8GjduLIQYOnToO++8I4RQVfXdd98dOnSo07pYCzRr1uw///nPG2+88cYbb8yZM8fX1/eNN95o0qRJQkKC2Wzes2eP+OXvbuDAgc7urDNZb3gQQnzwwQcBAQFhYWFCiCFDhqSlpWnlaWlpQ4YMcU7/apSTF+u4mIMHD8bGxrZq1crLyys2NvbBBx/UylevXh0SEjJlypRWrVrNmTPHuZ10uitXrphMpvDw8NhfvPvuu9qP7rvvvrZt2957773BwcHvvPOOc/vpdPPmzYuLi5s0adLAgQMbNmy4YsUKrfzbb78NCQkZO3Zsz5494+LiioqKnNvP2iM9Pd26alRKuWrVKuvf3WOPPebEjtUGL7/8cseOHSdOnNi/f/8GDRqkpaVp5efPn2/WrFlKSkpycnLbtm1//PFH5/azJnDEmkMVFhaeOHHC+rZBgwZRUVHac1ZW1qFDhyIjI62bWF1WeXl55UVrQojmzZsHBQVpzxkZGefOnUtISGjVqpUzeleLlJeXHzp06NSpUw0aNIiPj698DeHVq1d37drl5+fXp0+fejmpY5+ffvrp1KlTnTp1spbwd2dVWlp64MCBs2fP+vv7x8fHW//ihBCFhYU7d+60WCz9+vXz8vJyYidrCIEQAODSmCMEALg0AiEAwKURCAEALo1ACABwaQRCAIBLIxACAFwagRAA4NIIhAAAl0YgBAC4NAIhAMClEQgBAC7t/wMyasTKVXioSgAAAABJRU5ErkJggg==",
"image/svg+xml": [
"\n",
"\n"
],
"text/html": [
"\n",
"\n"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"using Plots, DelimitedFiles\n",
"final_grid = readdlm(\"SOR_result.txt\")\n",
"M = size(final_grid, 1)\n",
"N = size(final_grid, 2)\n",
"plt = heatmap(final_grid[:, 2:N-1])\n",
"display(plt)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "22b071c2",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"\u001b[32m\u001b[1mTest Passed\u001b[22m\u001b[39m"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Test if solution is identical with analytical solution\n",
"sol = analytical_solution(M)\n",
"# Bring solution into correct form \n",
"sol = reverse(transpose(sol), dims = 1)\n",
"@test maximum(abs.(sol - final_grid[:,2:N-1])) < 0.01 * M"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3130c9fb",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.9.1",
"language": "julia",
"name": "julia-1.9"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.9.1"
}
},
"nbformat": 4,
"nbformat_minor": 5
}