{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.linalg import solve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution strategy for the channel flow\n",
    "\n",
    "Rearranging the discretized equation so that all unknowns are on the left hand side, on the right hand side there are only known values.\n",
    "\n",
    "$$v_{i-1} - 2v_i + v_{i+1} = \\frac{dP/dx}{\\eta}\\Delta y^2$$\n",
    "\n",
    "There is one such equation for each internal grid point. At the boundaries, boundary conditions have to be used, e.g. $v_0 = v_\\mathrm{surf}$ and $v_{\\mathrm{ny}-1} = v_\\mathrm{bott}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For five grid points we thus have five equations in total:\n",
    "\n",
    "$$\\begin{eqnarray}\n",
    "i=0:\\qquad\\qquad & v_0 &    &      &    &     &    &     &   &     & = & v_\\mathrm{surf} \\\\\n",
    "i=1:\\qquad\\qquad & v_0 & -2 & v_1  & +  & v_2 &    &     &   &     & = & \\frac{dP/dx}{\\eta}\\Delta y^2 \\\\\n",
    "i=2:\\qquad\\qquad &     &    & v_1  & -2 & v_2 & +  & v_3 &   &     & = & \\frac{dP/dx}{\\eta}\\Delta y^2 \\\\\n",
    "i=3:\\qquad\\qquad &     &    &      &    & v_2 & -2 & v_3 & + & v_4 & = & \\frac{dP/dx}{\\eta}\\Delta y^2 \\\\\n",
    "i=4:\\qquad\\qquad &     &    &      &    &     &    &     &   & v_4 & = & v_\\mathrm{bott} \\\\\n",
    "\\end{eqnarray}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This system of equations can be written in matrix format as\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "1 & & & & \\\\ \n",
    "1 & -2 &  1 & & \\\\ \n",
    " & 1 & -2 & 1 & \\\\ \n",
    " & & 1 & -2 & 1 \\\\ \n",
    " & & & & 1 \\\\ \n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "v_0 \\\\ \n",
    "v_1 \\\\ \n",
    "v_2 \\\\ \n",
    "v_3 \\\\ \n",
    "v_4 \\\\ \n",
    "\\end{pmatrix}\n",
    "=\n",
    "\\begin{pmatrix}\n",
    "v_\\mathrm{surf} \\\\ \n",
    "\\frac{dP/dx}{\\eta}\\Delta y^2 \\\\ \n",
    "\\frac{dP/dx}{\\eta}\\Delta y^2 \\\\ \n",
    "\\frac{dP/dx}{\\eta}\\Delta y^2 \\\\ \n",
    "v_\\mathrm{bott} \\\\ \n",
    "\\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and is typically written with symbols\n",
    "\n",
    "$$ \\mathbf{Ax}=\\mathbf{b}$$\n",
    "\n",
    "Python comes with a module, `scipy.linalg` that includes function `solve`, which can be used to solve this type of systems of linear equations:\n",
    "\n",
    "    from scipy.linalg import solve\n",
    "    x = solve(A, b)\n",
    "  \n",
    "where `b` has to be a 1-D NumPy array of size `ny`, and includes the numerical values of the right hand side expressions; and `A` a 2-D NumPy array of size (`ny`,`ny`), and includes the numerical coefficients, as above. `x` will then contain the wanted velocity values. The following script utilizes this function to solve the channel flow problem.\n",
    "\n",
    "### Do\n",
    "\n",
    "- Complete the following script\n",
    "- Once complete:\n",
    "    - Which direction is the flow mostly going?\n",
    "    - How can you make it go more to the left/ro the right?\n",
    "    - What happens if you decrease the viscosity? Why?\n",
    "    - What happens if you remove (or comment out) the lines where B.C. are set?\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A is:\n",
      " [[ 1.  0.  0.  0.  0.]\n",
      " [ 1. -2.  1.  0.  0.]\n",
      " [ 0.  1. -2.  1.  0.]\n",
      " [ 0.  0.  1. -2.  1.]\n",
      " [ 0.  0.  0.  0.  1.]] \n",
      "b is:\n",
      " [ 0.00000000e+00 -1.37340000e-10 -1.37340000e-10 -1.37340000e-10\n",
      " -3.16880878e-10]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAD8CAYAAACPWyg8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VfWd//HXJwlJ2EluWLOYBVBBFCFAolgWN1BHumhF22qt87NabWut7ejYeUzrTKfTSketC4rWqm0ddGpbaW1LpYArAYILiqiEJEDYsxGSkP37++OepCkGSEhuzr3J+/l43Afnfs/3nvshuck753zPOV9zziEiItIVUX4XICIikUfhISIiXabwEBGRLlN4iIhIlyk8RESkyxQeIiLSZQoPERHpMoWHiIh0mcJDRES6LMbvAkIlKSnJpaen+12GiEhE2bRpU6lzbuSJ+vXZ8EhPTyc/P9/vMkREIoqZ7ehMPx22EhGRLlN4iIhIlyk8RESkyxQeIiLSZQoPERHpsogJDzNbYGYfmVmBmd3pdz0iIv1ZRISHmUUDDwMLgUnA1WY2yd+qRET6r4gID2AmUOCcK3TONQDLgUU+1yQiJ6G+qZkf/GELuyuP+F2KdEOkhEcysKvd8xKv7R+Y2Y1mlm9m+QcPHuy14kSk894rOcSz63cyf8lafvrXj6ipb/K7JDkJkRIe1kGb+0SDc8ucc9nOueyRI094db2I+CA7PZE1d8xl4RljeHB1AfOWrOX/8nfR0vKJH2kJY5ESHiVAarvnKcAen2oRkW4aN2Ig9y8+m9997RySEwbynd9s5vKHX2d9YZnfpUknRUp4bAQmmFmGmcUCi4EVPtckIt10dloCv735HB5YPJXy6gauWpbHzb/axM6yWr9LkxOIiBsjOueazOxWYCUQDTzpnNvic1ki0gPMjEVTk7l48hgef7WQpa9s529bD3D97HRunTeeofED/C5ROmDO9c3jjNnZ2U531RWJPPur6rh35Ue88FYJiYNi+fZFp3LVjFSiozoa+pSeZmabnHPZJ+oXKYetRKSfGD0sniVXnsWKW2aTNXII//q797j0Z6/xRkGp36VJOwoPEQlLU1KG89xXc1j6hWnUNDTxhSfW889Pb6TwYLXfpQkKDxEJY2bGwiljeflbc7hz4WnkFZZz0X2vcs8fPuBQbaPf5fVrCg8RCXvxA6K5aU4Wa+6Yy5XZqTz1ZhFzlqzh6TeLaWxu8bu8fknhISIRY+TQOH702Sm89I3zmDxuGP++YgsL7n+VNR8d8Lu0fkfhISIR5/Sxw/jVDbN44tpsWhxc/4uNXPvkBj7ef9jv0voNhYeIRCQz44JJo1l526f4t8sm8c7OChY+8Br/9vv3Ka9p8Lu8Pk/hISIRLTYmihtmZ/DKd+bxxVlpPLthJ3PuXcMTrxXS0KTxkFBReIhIn5AwOJYfLDqDv3zzPKafksB/vrSVi+57hb9u2UdfvRjaTwoPEelTJoweylPXz+Sp62cQEx3Fjb/cxDWPr+eDPVV+l9anKDxEpE+ae+oo/vLN87hn0WQ+3FfFpQ++xp0vbObg4Xq/S+sTFB4i0mfFREdxbW46a++Yxw3nZvCbTSXMW7KWR9YWUNfY7Hd5EU3hISJ93vBBA/jeZZN4+fY55GYF+MlfPuKC/3mFlzbv1XjISVJ4iEi/kZE0mMevzebZf57FkLgYbnn2LT7/2Do2l1T6XVrEUXiISL9zzvgkXvrGefzos1MoKq3h8ofe4Pbn32HfoTq/S4sYCg8R6Zeio4yrZ6ax5o653DQniz++u5d5S9bywKptHGnQeMiJKDxEpF8bGj+AOxeexqrb5zDvtJHct+pj5v90Lb9/ezctLRoPORaFh4gIkBYYxCNfmM7zX80lMCSW2557h88sfZNNOyr8Li0sKTxERNqZmZHIiltms+TKs9hbeYTPLX2Tr//v25RU1PpdWlhReIiIHCUqyrhiegpr7pjLN+aP569b9nH+T19hycqPqKlv8ru8sKDwEBE5hsFxMdx+0amsvmMuC84Yw0NrCpi7ZC3P5+/q9+MhCg8RkRNIHjGQBxafzW+/dg4pCQP57m82c/nDr7O+sMzv0nyj8BAR6aRpaQn89uZzeGDxVMqrG7hqWR43/XITO8pq/C6t1yk8RES6wMxYNDWZv317Lt++cCKvbjvIhf/zKj/601aq6hr9Lq/XKDxERE7CwNhovn7+BNbcMZfLp47jsVcLmXfvWn69fgdNzX1/EiqFh4hIN4weFs+SK8/iD7fOJmvkEO7+3ftc+rPXeW3bQb9LCymFh4hID5iSMpznvprD0i9Mo7axiS/9fAM3PLWR7Qer/S4tJBQeIiI9xMxYOGUsL39rDncuPI31ReVcfN+r/OAPW6isbfC7vB6l8BAR6WHxA6K5aU4Wa+6Yy5XZqTz9ZjFzl6zlqTeKaOwj4yEKDxGREBk5NI4ffXYKL33jPCaPG8b3//ABC+5/lTUfHoj4SahCFh5mdq+ZfWhmm83sd2Y2ot26u8yswMw+MrOL27Uv8NoKzOzOdu0ZZrbezLaZ2XNmFhuqukVEetrpY4fxqxtm8fi12bQ4uP6pjVz75AY+3n/Y79JOWij3PF4GznDOnQl8DNwFYGaTgMXAZGAB8IiZRZtZNPAwsBCYBFzt9QX4MXCfc24CUAHcEMK6RUR6nJlx4aTRrLztU3zv0tN5d1clC+5/le/9/j3Kquv9Lq/LQhYezrm/Ouda7yCWB6R4y4uA5c65eudcEVAAzPQeBc65QudcA7AcWGRmBswHfuO9/mng06GqW0QklGJjovjn8zJZ+515fDHnFP53wy7mLlnL468W0tAUOeMhvTXm8RXgz95yMrCr3boSr+1Y7QGgsl0QtbaLiESsxMGx3LPoDP7yzfOYlpbAD/+0lYvue4WVW/ZFxHhIt8LDzFaZ2fsdPBa163M30AT8urWpg025k2jvqJ4bzSzfzPIPHuzbF+iISN8wYfRQnv7KTJ66fgYx0VF89ZebuObx9WzZc8jv0o4rpjsvds5dcLz1ZnYdcBlwvvt7lJYAqe26pQB7vOWO2kuBEWYW4+19tO9/dD3LgGUA2dnZ4R/dIiKeuaeOYvb4JJ7dsJP7Xv6Yyx58nc9PT+XbF09k1NB4v8v7hFCebbUA+Bfgcudc+ym4VgCLzSzOzDKACcAGYCMwwTuzKpbgoPoKL3TWAFd4r78OeDFUdYuI+CUmOoprc9NZe8c8vnJuBi+8VcK8e9fy8JoC6hqb/S7vH1iojq2ZWQEQB7Te8D7POXeTt+5uguMgTcBtzrk/e+2XAPcD0cCTzrkfeu2ZBAfQE4G3gS865457ekJ2drbLz8/v8f+XiEhvKTxYzX/96UNWbd1PSsJA7lx4GpdOGUvwPKLQMLNNzrnsE/aLhIGZk6HwEJG+4o2CUv7jjx/w4b7DZJ+SwL9dNomzUkec+IUnobPhoSvMRUTC3Lnjk3jpG+fxo89OobishkUPv8Htz73DvkN1vtWk8BARiQDRUcbVM9NYc8dcbpqTxR8372XekrXcv+pjjjT0/niIwkNEJIIMjR/AnQtPY9Xtc5h32kjuX7WN+T9dy+/eLqGlpfeGIRQeIiIRKC0wiEe+MJ3nbswhMCSWbz33Lp9Z+iabdpT3yvsrPEREItiszAArbpnNkivPYm/lET63dF2vXGCo8BARiXBRUcYV01P410tOByA2OvS/2hUeIiJ9RFFpDWaQmjgo5O+l8BAR6SOKy2pIHjGQ+AHRIX8vhYeISB9RVFpDRtLgXnkvhYeISB/gnFN4iIhI15TVNHC4ron0gMJDREQ6qbi0BoCMkQoPERHppMLW8NCeh4iIdFZxaQ0xUUZKwsBeeT+Fh4hIH1BUWkNa4iBieuECQVB4iIj0Cb15phUoPEREIl5Li6O4rIZ0hYeIiHTW/sN11DW2aM9DREQ6r+igd6aVwkNERDqrqEzhISIiXVR0sIa4mCjGDIvvtfdUeIiIRLjisuCZVlFR1mvvqfAQEYlwhaU1vXZPq1YKDxGRCNbU3MKu8tpeu6dVK4WHiEgE2115hMZm12v3tGql8BARiWBFvXw33VYKDxGRCNYaHhrzEBGRTisurWFoXAxJQ2J79X0VHiIiEayorJb0pMGY9d5puqDwEBGJaEWl1b16ZXkrhYeISISqb2pmd8WRXr2bbquQh4eZ3WFmzsySvOdmZj8zswIz22xm09r1vc7MtnmP69q1Tzez97zX/Mx6e/9MRCQM7SqvpcVBZl8LDzNLBS4EdrZrXghM8B43Aku9vonAvwOzgJnAv5tZgveapV7f1tctCGXdIiKRoKi0FqBP7nncB3wXcO3aFgHPuKA8YISZjQUuBl52zpU75yqAl4EF3rphzrl1zjkHPAN8OsR1i4iEvaLSaoBev0AQQhgeZnY5sNs59+5Rq5KBXe2el3htx2sv6aC9o/e80czyzSz/4MGD3fwfiIiEt6LSWhIHxzJ80IBef++Y7rzYzFYBYzpYdTfwr8BFHb2sgzZ3Eu2fbHRuGbAMIDs7u8M+IiJ9hV9nWkE3w8M5d0FH7WY2BcgA3vXGtlOAt8xsJsE9h9R23VOAPV773KPa13rtKR30FxHp14pLazl3fJIv7x2Sw1bOufecc6Occ+nOuXSCATDNObcPWAFc6511lQMccs7tBVYCF5lZgjdQfhGw0lt32MxyvLOsrgVeDEXdIiKRorahiX1VdWT28j2tWnVrz+Mk/Qm4BCgAaoHrAZxz5Wb2H8BGr989zrlyb/lm4ClgIPBn7yEi0m8Vt55p5cNgOfRSeHh7H63LDrjlGP2eBJ7soD0fOCNU9YmIRJq2u+n6NOahK8xFRCJQcZl3N92kQb68v8JDRCQCFR6sYcyweAbF+jH6oPAQEYlIxWU1vu11gMJDRCQiFZXWkJE0xLf3V3iIiESYQ7WNlNc0kKE9DxER6ayiMn+mnm1P4SEiEmGKvdN0/bpAEBQeIiIRp7C0hiiD1EQdthIRkU4qLq0hOWEgcTHRvtWg8BARiTBFpTW+jneAwkNEJKI45ygurfFl6tn2FB4iIhGktLqBw/VNvkw9257CQ0QkgrTe08qvGyK2UniIiESIxuYWVm3dD/gfHv7cUUtERDqtvqmZFzbtZukrBewqP0JOZiIpCf6dpgsKDxGRsFXX2MzyDTt57NVC9h6q46zUEXz/nyYz/7RReFN8+0bhISISZmrqm3h2/U6WvVbIwcP1zEhP4MefO5PzJiT5HhqtFB4iImHicF0jz6zbwROvFVJR28i54wM8ePXZ5GQG/C7tExQeIiI+q6xt4BdvFPOLN4qoqmti3qkjuXX+BKafkuB3acek8BAR8UlpdT0/f72IX67bQXV9ExdNGs3X509gSspwv0s7IYWHiEgvO1BVx2OvFvLr9Tuob2rh0iljuXX+eE4bM8zv0jpN4SEi0kt2Vx7hsVe2s3zjLppbHIumjuNrc8czfpR/MwKeLIWHiEiI7Syr5ZG1BbzwVgkAn5uWws1zszjF55sbdofCQ0QkRAoOVPPI2gJefGcP0VHG1TPT+OqcLJJHDPS7tG5TeIiI9LAP91Xx0OoCXnpvL3ExUXz5nHRu/FQmo4fF+11aj1F4iIj0kPdKDvHg6m389YP9DI6N5qY5WdwwO4OkIXF+l9bjFB4iIt20aUcFD63expqPDjIsPoZvnj+B689NZ8SgWL9LCxmFh4jIScorLOPB1dt4o6CMhEED+M7Fp/Kl3FMYFj/A79JCTuEhItIFzjle21bKQ6sL2FBcTtKQOO6+5HSumZXG4Lj+8yu1//xPRUS6wTnH6g8P8LPVBby7q5Ixw+L5/j9NYvHMNOIHRPtdXq8LaXiY2deBW4Em4CXn3He99ruAG4Bm4BvOuZVe+wLgASAaeMI5999eewawHEgE3gK+5JxrCGXtIiIALS2OlVv28eDqAj7YW0VKwkD+6zNT+Nz0ZOJi+l9otApZeJjZPGARcKZzrt7MRnntk4DFwGRgHLDKzCZ6L3sYuBAoATaa2Qrn3AfAj4H7nHPLzexRgsGzNFS1i4g0tzj+uHkPD60uYNuBajKSBrPkyrNYNHUcA6I1CWso9zxuBv7bOVcP4Jw74LUvApZ77UVmVgDM9NYVOOcKAcxsObDIzLYC84FrvD5PA99H4SEiIdDY3MLv397NI2u3U1Raw8TRQ3hg8VQuO3Mc0VHhMZdGOAhleEwEzjOzHwJ1wB3OuY1AMpDXrl+J1waw66j2WUAAqHTONXXQX0SkR9Q3NfObTSUsXbudkoojTBo7jEe/OI2LJo0hSqHxCd0KDzNbBYzpYNXd3rYTgBxgBvC8mWUCHX0XHNDRfqA7Tv+O6rkRuBEgLS3tROWLiHQ41esPLg+PqV7DWbfCwzl3wbHWmdnNwG+dcw7YYGYtQBLBPYfUdl1TgD3eckftpcAIM4vx9j7a9z+6nmXAMoDs7OwOA0ZEBIJTvf56/Q6WvVpEaXU9M9MT+ckVZzJ7fPhM9RrOQnnY6vcExyrWegPisQSDYAXwrJn9D8EB8wnABoJ7GBO8M6t2ExxUv8Y558xsDXAFwTOurgNeDGHdItKHVdU18st2U73OHp/E1+efzawwnOo1nIUyPJ4EnjSz94EG4DpvL2SLmT0PfEDwFN5bnHPNAGZ2K7CS4Km6Tzrntnjb+hdguZn9J/A28PMQ1i0ifVBlbQNPvlHMU95Ur/NPG8Ut88aH9VSv4cyCv8/7nuzsbJefn+93GSLis9Lqep54rYhfriumpqGZiycHp3o9Izn8p3r1g5ltcs5ln6ifrjAXkT5pf1Udy9pN9XrZmeO4ZV5WRE31Gs4UHiLSp+yuPMKja7fzXP7fp3q9Zd54skZG3lSv4UzhISJ9wo6yGh5Zs50X3irBDK6YnsLNc8aTFhjkd2l9ksJDRCJawYFqHllTwIvvBqd6vWZW35nqNZwpPEQkIn24r4oHVxfwp/f2Eh8TzfXeVK+j+tBUr+FM4SEiEaX9VK9D4mK42ZvqNdAHp3oNZwoPEYkIm3ZU8ODqbaztR1O9hjOFh4iELecceYXlPLh6G29uLyNxcCzfufhUrs09haH9YKrXcKbwEJGw0zrV64Ort7GxuKJtqtcv5KQxKFa/tsKBvgsiEjacc/xt6wEeXBOc6nXs8Hh+cPlkrpqR2i+neg1nCg8R8V1Li+Mv3lSvW/dWkZo4kB99dgqfnda/p3oNZwoPEfFNU3MLL723t22q10xN9RoxFB4i0usam1v43du7eWRNAcVltUwcPYSfXX02l04Zq6leI4TCQ0R6zdFTvU4ep6leI5XCQ0RCrq6xmf/dsJPHXilkX1UdU1NHcM+iycw7VVO9RiqFh4iETE19E7/K28Hjr3lTvWYksuTKszh3fEChEeEUHiLS46rqGnnmzWJ+/noRFbWNnDchiVvnaarXvkThISI9prK2gSdfL+IXbxZz2Jvq9db545mWpqle+xqFh4h0W2l1PY+/Vsiv1u2gpqGZBZPHcOv88ZrqtQ9TeIjISdtfVcdjrxTy7Ia/T/V667zxnDpmqN+lSYgpPESky0oqann0le08v7GEZuf49NRkvjYvS1O99iMKDxHptOLSGpau1VSvovAQkU4oOHCYh9ds58V3dhMTHcUXvKlex2mq135L4SEiHWpucWwuqeSJ14r40/vBqV5vmJ3B/ztPU72KwkNEPC0tjg/3HWZdYRl5hWWsLyyjqq5JU71KhxQeIv1US4tj24Fq1m0vZV1hGeuLyqmsbQQgLXEQC88YS05WIvNPHc3wQZq1T/6RwkOkn3DOUXCgmrzCMm/vopzymgYAkkcM5ILTR5ObGSAnK0CyxjLkBBQeIn2Uc47C0hrWbQ8ehsorLKe0uh6AccPjmXvqSHIyA+RmBkhN1NlS0jUKD5E+wjlHcVltcM/CC4wDh4NhMXpYHLPHB8jNCpCTGSAtcZBuTCjdovAQiVDOOXaVH2FdYSl5heWs217Gvqo6AEYOjWvbq8jNCpAeUFhIz1J4iESQkopa1m0PjlmsLyxnd+URAJKGxDIrM9AWGFkjByssJKRCFh5mNhV4FIgHmoCvOec2WPAT/QBwCVALfNk595b3muuA73mb+E/n3NNe+3TgKWAg8Cfgm845F6raRcLFnsojbYeh1hWWUVIRDIuEQQPIyQzw1TmZ5GQGmDBqiMJCelUo9zx+AvzAOfdnM7vEez4XWAhM8B6zgKXALDNLBP4dyAYcsMnMVjjnKrw+NwJ5BMNjAfDnENYu4ov9VXXBoNheRl5RGTvKagEYPnAAszISuWF2BrlZASaOGqppW8VXoQwPBwzzlocDe7zlRcAz3p5DnpmNMLOxBIPlZedcOYCZvQwsMLO1wDDn3Dqv/Rng0yg8pA84cLiubbxifWEZhaU1AAyNj2FWRoBrc9PJyUzk9DHDFBYSVkIZHrcBK81sCRAFnOO1JwO72vUr8dqO117SQfsnmNmNBPdQSEtL6/7/QKSHlVbXe6fNBvcuth8MhsWQuBhmZiRy9cw0crMCnD52GNEKCwlj3QoPM1sFjOlg1d3A+cC3nHMvmNnngZ8DFwAd/US4k2j/ZKNzy4BlANnZ2RoTEd+V1zSwvu2ivDI+3l8NwODYaGZkJHJldiq5mQEmjxtGTHSUz9WKdF63wsM5d8Gx1nmHl77pPf0/4AlvuQRIbdc1heAhrRKCh67at6/12lM66C8SdiprG1hfVN52ncWH+w4DMHBANNnpCXz67GRyMgNMSR7OAIWFRLBQHrbaA8whGADzgW1e+wrgVjNbTnDA/JBzbq+ZrQT+y8xaJzu+CLjLOVduZofNLAdYD1wLPBjCukU67dCRRjYUlbcdhtq6rwrnIC4miuz0BO64aCK5WQGmJI8gNkZhIX1HKMPj/wEPmFkMUIc3FkHwbKlLgAKCp+peD+CFxH8AG71+97QOngM38/dTdf+MBsvFJ4frGtlY3LpnUc6WPYdocRAbE8X0tARuOz8YFmelDicuJtrvckVCxvrq5RLZ2dkuPz/f7zIkwlXXN7GxOLhnkbe9jPd2e2ERHcXUtBHBGwlmBjg7bQTxAxQWEvnMbJNzLvtE/XSFuUg7tQ1N5BdXtA1wby45RHOLY0C0MTV1BLfMG09uZoBppyQoLKRfU3hIv3akoZm3dla0DXC/W1JJY7MjJso4M2U4N3lXcE8/JYFBsfpxEWmlnwbpV+oag2GRV1hO3vYy3tlVSUNzC9FRxhnJw7lhdia5WQGyT0lgcJx+PESORT8d0qfVNzXzzs7KtsNQb+2spKGphSiDM5KH8+Vz08nNDJCdnsDQeM2WJ9JZCg/pUxqaWthcUtl2I8FNOyqob2rBDCaNHca1OaeQkxlgRkYiwwcqLEROlsJDIlpjcwubSw613fIjv7iCI43NAJw2ZijXzEojNzPAzIxERgyK9blakb5D4SERpam5hff3VLXtWeQXl1PbEAyLU0cP5aoZqeRkJjIrI0DCYIWFSKgoPCSsNbc4tuw51HYF98biCqrrmwAYP2oIn5uWQm5WcM8iaUicz9WK9B8KDwkrLS2OD/ZWtR2GWl9UzuG6YFhkjhzM5VPHtV2YN3KowkLELwoP8VVLi+Oj/YfbDkNtKCrn0JFGANIDg7jszLHkeGExeli8z9WKSCuFh/Qq5xwf769uOwy1vqiMitpgWKQmDuTiyaPJzQqGxdjhA32uVkSOReEhIeWcY/vB6rYbCeYVllFW0wBA8oiBzD+tNSwSSUkY5HO1ItJZCg/pUc45Cktr2vYs8grLKa2uB2Ds8HjmTBxJTmaA3KwAqYkKC5FIpfCQbnHOsaOstu0K7rzCMvZXBcNi1NA4zh0fINcLi7TEQZhpalWRvkDhIV3inKOk4kjbAHdeYRl7D9UBkDQkru0QVG5mgIykwQoLkT5K4SEnVFJRS17h36dW3V15BIDA4NjgmVBZAXIzE8kaOURhIdJPKDzkE/YeOtIWFOsKy9hVHgyLhEEDmJUR4MZPBe88O2GUwkKkv1J4CPur6toNcJdRXFYLwPCBA5iVkcj152SQmxXg1NFDiYpSWIiIwqNfOnC4ru202bztZRSW1gAwND6GWRmJfNG78+zpY4cRrbAQkQ4oPPqBsur64JhFYSl5heUUHKgGYEhcDDPSE1g8M5XczCQmjVNYiEjnKDz6oIqaBtYXlbWdEfXx/mBYDIqNZkZ6IldMTyEnM8AZ44YREx3lc7UiEokUHn3AodpG8orK2sYtPtx3GICBA6LJTk9g0dRkcrMCTEkezgCFhYj0AIVHBKqqa2RDYXnbdRYf7K3COYiLiSI7PYFvXziR3KwAZ6aMIDZGYSEiPU/hEQEO1zWSX1zBOm/PYsueQ7Q4iI2JYlraCG47fyI5mYlMTRtBXEy03+WKSD+g8AhDNfVNbCwu9wa5y3h/9yGaWxwDoo2zUxO4df4EcjMDnJ02gvgBCgsR6X0KjzBQ29DEph0VbddZbC45RFOLIybKmJo6gpvnZJGbFWBaWgIDYxUWIuI/hYcP6hqb2bSjom2A+92SShqbHdFRxpkpw9uu4J5+SgKDYvUtEpHwo99MvaCusZm3d1a2DXC/s7OShuYWogymJA/nK7MzyM0MkJ2eyJA4fUtEJPzpN1UI1Dc18+6uQ951FqW8tbOShqZgWEweN5wvn5tOTmYiM9ITGRo/wO9yRUS6TOHRAxqaWthcUtl2I8FNOyqoa2zBDE4fM4wv5ZxCbmaAGRmJDB+osBCRyNet8DCzK4HvA6cDM51z+e3W3QXcADQD33DOrfTaFwAPANHAE865//baM4DlQCLwFvAl51yDmcUBzwDTgTLgKudccXfq7q7G5hbe232obYA7v7iCI43NAJw2ZiiLZ6SRmxVgVkYiIwbF+lmqiEhIdHfP433gs8Bj7RvNbBKwGJgMjANWmdlEb/XDwIVACbDRzFY45z4Afgzc55xbbmaPEgyepd6/Fc658Wa22Ot3VTfr7pKm5hbe31PVNsCdX1xOTUMwLCaOHsLns4O3+5iVGSBxsMJCRPq+boWHc24r0NGcDouA5c65eqDIzAqAmd66Audcofe65cAiM9sKzAeu8fo8TXCPZqm3re977b8BHjIzc8657tTSWu/HAAAHGklEQVR+PM0tjg/2VLXdSHBjUTmH65sAyBo5mM9MSyY3M4lZmYkkDYkLVRkiImErVGMeyUBeu+clXhvArqPaZwEBoNI519RB/+TW1zjnmszskNe/NBSFP7BqG0+8XsjhumApmUmDueyscW3Tq44aGh+KtxURiSgnDA8zWwWM6WDV3c65F4/1sg7aHNDRjZbccfofb1uffFOzG4EbAdLS0o5R2vGNHR7PpVPGBqdXzQwwZrjCQkTkaCcMD+fcBSex3RIgtd3zFGCPt9xReykwwsxivL2P9v1bt1ViZjHAcKD8GLUuA5YBZGdnn9Rhrc/PSOXzM1JP3FFEpB8L1S1XVwCLzSzOO4tqArAB2AhMMLMMM4slOKi+whu/WANc4b3+OuDFdtu6zlu+AlgdyvEOERE5sW6Fh5l9xsxKgFzgJTNbCeCc2wI8D3wA/AW4xTnX7O1V3AqsBLYCz3t9Af4FuN0bXA8AP/fafw4EvPbbgTu7U7OIiHSf9dU/4rOzs11+fv6JO4qISBsz2+Scyz5RP80UJCIiXabwEBGRLlN4iIhIlyk8RESkyxQeIiLSZX32bCszOwjsCOFbJBGiW6T0sEioMxJqBNXZ0yKhzkioEXq2zlOccyNP1KnPhkeomVl+Z05n81sk1BkJNYLq7GmRUGck1Aj+1KnDViIi0mUKDxER6TKFx8lb5ncBnRQJdUZCjaA6e1ok1BkJNYIPdWrMQ0REukx7HiIi0mUKj6OYWaKZvWxm27x/E47R7zqvzzYzu65d+w/NbJeZVR/VP87MnjOzAjNbb2bpPtY43cze82r5mXnzCJvZVDPLM7N3zCzfzGZ2tF2/6/TWfd3MPjKzLWb2k3Ct01t/h5k5M0sKxzrN7F4z+9DMNpvZ78xsxEnUtsD7fhSY2SfufH28z7+Z3eW1f2RmF3d2myejp+s0s1QzW2NmW73P4jfDsc5266LN7G0z+2O3i3TO6dHuAfwEuNNbvhP4cQd9EoFC798EbznBW5cDjAWqj3rN14BHveXFwHM+1riB4G30DfgzsNBr/2u75UuAtT5/LY9V5zxgFRDnPR8VjnV661IJTkGwA0gKxzqBi4AYb/nHHW33BHVFA9uBTCAWeBeY1JnPPzDJ6x8HZHjbie7MNk/i6xeKOscC07w+Q4GPw7HOdq+7HXgW+GN3anTOKTw6+MZ9BIz1lscCH3XQ52rgsXbPHwOuPqrP0eGxEsj1lmMIXtBjvV2j1//Djvp5NV7Vrv1Zv76WJ6jzeeCCcPieH69O7/lvgLOAYrofHiGrs137Z4Bfd7GuXGBlu+d3AXd15vN/dN/Wfp3Z5kl8/Xq8zg7e40XgwnCsk+AMrX8D5tMD4aHDVp802jm3F8D7d1QHfZKBXe2el3htx9P2GhecFOsQwUmvervGZG/56HaA24B7zWwXsITgB7E7QlXnROA8b3f9FTObEY51mtnlwG7n3LvdrC+kdR7lKwT3SrqiMz8Px/r8H6/erv6M+VFnG+/Q0dnA+jCt837gu0BLN+sDOjGHeV9kZquAMR2suruzm+ig7USnrXXpNSGs8Xh13Ax8yzn3gpl9nuAsjsedw96nOmMIHpLJAWYAz5tZpvP+vAqHOs1skLftizq5/eCb+PP1bH3vu4Em4NedfK8TvWd36uroD9vunhoaijqDLzIbArwA3OacqzrpCo9fQ2f6HOvzeBlwwDm3yczmdrM+oJ+Gh3PumL8QzWy/mY11zu01s7HAgQ66lQBz2z1PAdae4G1LCB7/LjGzGGA4UO5DjSXecvv2Pd7ydUDrgN//AU8c5//jZ50lwG+9sNhgZi0E7+1zMIzqzCJ4zPldb1w6BXjLzGY65/aFUZ2t274OuAw4/3ghfAytn+0Ot31Un6M//8d77Ym22VUhqdPMBhAMjl87537bzRpDVeflwOVmdgkQDwwzs18557540lV297hXX3sA9/KPg5I/6aBPIlBE8K/fBG858ag+R4953MI/DnA971eNwEaCf7W3Dpxe4rVvBeZ6y+cDm/z8Wh6nzpuAe7zliQR3009q/CiUdR71+mK6P+YRqq/nAuADYORJ1hVDcGA+g78P8E7uzOcfmMw/DvAWEhwwPuE2w6ROA54B7u9ObaGu86jXzkUD5j3/IHjc8G/ANu/f1h+8bOCJdv2+AhR4j+vbtf+EYPq3eP9+32uPJ/jXfAHBs14yfawxG3if4JkYD/H3i0VnA5u8D996YLrPX8tj1RkL/Mpb9xYwPxzrPOo9iul+eITq61lAMIDf8R6PnkRtlxA802g7cLfXdg9w+Yk+/wQPyW0neELAwuNtswd+vnu0Tu9nxgGb2339PvHHg991HrXtufRAeOgKcxER6TKdbSUiIl2m8BARkS5TeIiISJcpPEREpMsUHiIi0mUKDxER6TKFh4iIdJnCQ0REuuz/A6s8UlZxciYwAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Define physical parameters\n",
    "visc = 1e19                    # Viscosity of the rock in the channel, Pa s\n",
    "dpdx = -2800*9.81*4000/500e3   # Pressure gradient within the channel, Pa m^-1\n",
    "\n",
    "# Define problem geometry\n",
    "ny = 5      # number of grid points\n",
    "L = 10e3    # size of the problem (width of the channel)\n",
    "\n",
    "# Define values for boundary conditions, m s^-1\n",
    "vx_surf = 0\n",
    "vx_bott = -0.01/(60*60*24*365.25)\n",
    "\n",
    "# Calculate grid spacing and y coordinates\n",
    "dy = L / (ny-1)\n",
    "y = np.linspace(0, L, ny)\n",
    "\n",
    "# Create the empty (zero) coefficient and right hand side arrays\n",
    "A = np.zeros((ny,ny))  # 2-dimensional array, ny rows, ny columns\n",
    "b = np.zeros(ny)\n",
    "\n",
    "# Set B.C. values in the coefficient array and in the r.h.s. array\n",
    "A[0, 0] = 1\n",
    "b[0] = vx_surf\n",
    "A[ny-1, ny-1] = 1\n",
    "b[ny-1] = vx_bott\n",
    "\n",
    "# Set remaining (internal grid point) coefficients in the array `A`.\n",
    "# TODO: Complete the for loop so that it will write the coefficient \n",
    "# values in the array `A`. The for loop loops over all the rows\n",
    "# of the matrix. On each row, you need to fill in three coefficients.\n",
    "# In python, elements of 2D arrays are referenced to like 'A[row, col]'\n",
    "for iy in range(1, ny-1):\n",
    "    .\n",
    "    .\n",
    "    .\n",
    "    \n",
    "\n",
    "# Set remaining (internal grid point) coefficients in the r.h.s. array b\n",
    "# TODO: Create a for loop that loops over the internal grid points (i.e.\n",
    "# from 1 to ny-1), and fills in the `b` array with r.h.s. values\n",
    ".\n",
    ".\n",
    ".\n",
    "\n",
    "# For debugging purposes, plot A (if it is small enough) and b\n",
    "if ny < 15:\n",
    "    print(\"A is:\\n\", A, \"\\nb is:\\n\", b.T)\n",
    "    \n",
    "# Solve it!\n",
    "vx = solve(A, b)\n",
    "\n",
    "# Plot it!\n",
    "plt.plot(vx * (60*60*25*365.25), -y)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating strain rates\n",
    "\n",
    "Strain rate in 2D is defined as \n",
    "\n",
    "$$\\underline{\\underline{\\dot\\epsilon}} = \\begin{pmatrix}\n",
    "\\frac{\\partial v_x}{\\partial x} & \\frac{\\partial v_x}{\\partial y} \\\\\n",
    "\\frac{\\partial v_y}{\\partial x} & \\frac{\\partial v_y}{\\partial y} \\\\\n",
    "\\end{pmatrix}$$\n",
    "\n",
    "In our case the only relevant component is the shear strain rate $\\dot\\epsilon_{xy} = \\frac{\\partial v_x}{\\partial y}$ (why?).\n",
    "\n",
    "### Do\n",
    " \n",
    " - Discretize the shear strain rate expression\n",
    " - Complete the script below so that it will loop over all the grid points, calculate the shear strain, and store it in the array `e`\n",
    " - How does the resulting plot look like? Note that the _shear stress_ is defined as $\\sigma_{xy} = \\eta\\dot\\epsilon_{xy}$. Where is the shear stress smallest?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEJCAYAAABhbdtlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XeYVPXZ//H3Te+9l6W3pamMYIk1KFhR0diiWCLG6C/tiQKiBruYRKOJJdh9EkWlyIoFQSEaFWVR2YWlLUVYeu8sW+7fH3PIMyEzLDC7O7O7n9d1zcWZ7/me2XsPs/vZc87MPebuiIiIRFMp0QWIiEjyUkiIiEhMCgkREYlJISEiIjEpJEREJCaFhIiIxKSQEBGRmBQSIiISk0JCRERiqhLPxmZ2BTAG6AH0d/f0iHWjgJuBAuCX7j4tGB8MPAVUBl5098eC8Q7AeKAR8C1wnbsfKKqGJk2aePv27eP5NkREKpy5c+dudvemRc2LKySA+cBlwN8iB80sFbgK6Am0AmaYWddg9TPAOUAOMMfM0tw9CxgLPOnu483secIB81xRBbRv35709PSipomISAQz++FI5sV1usndF7r74iirhgDj3T3X3VcA2UD/4Jbt7suDo4TxwBAzM+BsYEKw/WvAJfHUJiIi8SupaxKtgdUR93OCsVjjjYHt7p5/yLiIiCRQkaebzGwG0CLKqtHuPiXWZlHGnOih5IeZH6um4cBwgJSUlFjTREQkTkWGhLsPPIbHzQHaRtxvA6wNlqONbwYamFmV4Ggicn60msYB4wBCoZB6nYuIlJCSOt2UBlxlZtWDVy11Ab4B5gBdzKyDmVUjfHE7zcMfajETuDzYfhgQ6yhFRERKSVwhYWaXmlkOcDLwvplNA3D3BcDbQBbwEXC7uxcERwl3ANOAhcDbwVyAEcBvzSyb8DWKl+KpTURE4mdl/ZPpQqGQ6yWwIiJHx8zmunuoqHl6x7WISBmzeP0uHv9oEaXxR368b6YTEZFSciC/kGdnZfPMzGzq1qjKT09qR6sGNUv0ayokRETKgHmrt3PXhAwWb9jFkONacd+FqTSuU73Ev65CQkQkie07UMAT0xfz0r9W0KxuDV4aFuLHPZqX2tdXSIiIJKkvl21m5MRMVm3dyzUDUhh5Xnfq1ahaqjUoJEREkszO/Xk8+sEi3vxmFe0a1+LNW07i5E6NE1KLQkJEJInMyNrA6Hcz2bQrl+Gnd+Q3A7tSs1rlhNWjkBARSQJbdudy/3tZpM1bS/cWdRl3XYi+bRskuiyFhIhIIrk7afPWMiZtAbtz8/ntOV35+RmdqFYlOd7GppAQEUmQtdv3cc+78/l00UaOa9uAxy/vQ9fmdRNd1n9QSIiIlLLCQufNOat49INFFBQ6916Yyg2ntKdypWifmpBYCgkRkVK0YvMeRk7M4OsVWzm1c2MevbQPKY1rJbqsmBQSIiKlIL+gkJe/WMGfPl5CtSqVGDu0Nz8JtSX86c3JSyEhIlLCFq7byYiJGWTk7OCc1OY8dEkvmterkeiyjohCQkSkhOTmF/DMp9k8O2sZDWpV5ZlrTuD83i2S/ughkkJCRKQEfLtqGyMmZLB0424uO741916YSsPa1RJd1lFTSIiIFKO9B/L547QlvPLlClrWq8ErN57IWd2aJbqsYxbvx5f+wcwWmVmGmU02swYR60aZWbaZLTazQRHjg4OxbDMbGTHewcy+NrOlZvZW8BnYIiJlxhfZmxn05894+YsV/HRAO6b95vQyHRAQ/yfTTQd6uXsfYAkwCsDMUoGrgJ7AYOBZM6tsZpWBZ4DzgFTg6mAuwFjgSXfvAmwDbo6zNhGRUrFjXx4jJmRw7YtfU6VSJd6+9WQevKQXdUu5Y2tJiOt0k7t/HHF3NnB5sDwEGO/uucAKM8sG+gfrst19OYCZjQeGmNlC4GzgmmDOa8AY4Ll46hMRKWnTFqzn3nfns2XPAW47sxO/+nEXalRNXEO+4lac1yRuAt4KllsTDo2DcoIxgNWHjA8AGgPb3T0/ynwRkaSzaVcuY9IW8H7mOnq0rMdLw06kd5v6iS6r2BUZEmY2A2gRZdVod58SzBkN5AP/OLhZlPlO9NNbfpj5sWoaDgwHSElJiVm7iEhxc3cmf7eGB6ZmsTe3gDsHdWP46R2pWjk5GvIVtyJDwt0HHm69mQ0DLgR+7O4Hf7HnAG0jprUB1gbL0cY3Aw3MrEpwNBE5P1pN44BxAKFQKGaYiIgUpzXb9zF6ciazFm+iX7uGjB3am87NkqshX3GL63STmQ0GRgBnuPveiFVpwBtm9gTQCugCfEP4iKGLmXUA1hC+uH2Nu7uZzSR8TWM8MAyYEk9tIiLFpbDQ+cfXP/DYh4twYMxFqVx/cnsqJWFDvuIW7zWJvwLVgenBOwhnu/vP3X2Bmb0NZBE+DXW7uxcAmNkdwDSgMvCyuy8IHmsEMN7MHgK+A16KszYRkbgt27SbkRMzmLNyG6d1acIjl/ambaPkbchX3Oz/zhCVTaFQyNPT0xNdhoiUM/kFhYz7fDl/nrGUGlUqce+FqVzer02ZaqlxOGY2191DRc3TO65FRA6xYO0ORkzMYP6anQzu2YIHLulJs7ployFfcVNIiIgE9ucV8JdPl/L8P5fTsFY1nrv2BM7r3TLRZSWUQkJEBEhfuZUREzNYtmkPl/drwz0X9KBBLXUHUkiISIW2JzefP0xbzGtfraRV/Zq8flN/Tu/aNNFlJQ2FhIhUWJ8t2cSoSZms3bGPYSe3585B3ahdXb8WI2lviEiFs33vAR56fyET5ubQsWlt3rn1ZELtGyW6rKSkkBCRCuXDzHXcO2UB2/Ye4I6zOnPH2Z3LVUO+4qaQEJEKYeOu/fx+ygI+nL+enq3q8dpNJ9KzVflryFfcFBIiUq65OxPm5vDg1Cz25xcyYnB3bjmtA1XKaUO+4qaQEJFya/XWvdw9OZPPl27mxPYNeWxoHzo1rZPossoUhYSIlDuFhc7rX63k8WmLMeDBIT25dkC7CtGQr7gpJESkXMneuIsREzOZ+8M2zujalEcu603rBjUTXVaZpZAQkXIhr6CQcZ8t56kZS6lVvTJP/KQvlx7futw05EsUhYSIlHnz1+zgzgkZLFy3kwv6tGTMRT1pWrd6ossqFxQSIlJm7c8r4M8zlvLC58tpVLsaf7uuH4N6Rvu0ZTlWCgkRKZO+WbGVkRMzWL55D1eG2nL3+T2oX6tqossqdxQSIlKm7M7NZ+yHi/jf2T/QtlFN/n7zAH7UpUmiyyq3FBIiUmbMXLyR0ZMyWbdzPzed2oHfDepKrWr6NVaS4nrLoZk9aGYZZva9mX1sZq2CcTOzp80sO1h/QsQ2w8xsaXAbFjHez8wyg22eNr0kQUQC2/Yc4Ldvfc+Nr8yhdvUqTLztFO67KFUBUQrifV/6H9y9j7sfB0wF7gvGzwO6BLfhwHMAZtYI+D0wAOgP/N7MGgbbPBfMPbjd4DhrE5Eyzt2ZmrGWgU/8k7R5a/nl2Z2Z+ssfcUJKw6I3lmIRVwy7+86Iu7UBD5aHAK+7uwOzzayBmbUEzgSmu/tWADObDgw2s1lAPXf/Khh/HbgE+DCe+kSk7Nqwcz/3vjufj7M20Lt1ff7+swH0aFkv0WVVOHEfq5nZw8D1wA7grGC4NbA6YlpOMHa48Zwo47G+5nDCRx2kpKTE9w2ISFJxd95OX81D7y/kQH4ho87rzs0/UkO+RClyr5vZDDObH+U2BMDdR7t7W+AfwB0HN4vyUH4M41G5+zh3D7l7qGlTfcygSHmxastefvrS14yYmElqy3p89OvTufWMTgqIBCrySMLdBx7hY70BvE/4mkMO0DZiXRtgbTB+5iHjs4LxNlHmi0gFUFDovPrlSv44bTGVKxkPX9qLq09MUUO+JBDvq5u6RNy9GFgULKcB1wevcjoJ2OHu64BpwLlm1jC4YH0uMC1Yt8vMTgpe1XQ9MCWe2kSkbFiyYRdDn/uSB6dmcXKnxkz/7enq2JpE4r0m8ZiZdQMKgR+AnwfjHwDnA9nAXuBGAHffamYPAnOCeQ8cvIgN3Aa8CtQkfMFaF61FyrED+YU8/89l/OXTpdSpXoWnrjqOi/u2UkO+JGPhFyCVXaFQyNPT0xNdhogchXmrtzNiYgaL1u/ior6tGHNRKo3rqCFfaTKzue4eKmqe3okiIqVm34ECnpyxhBc/X07TutV54foQ56Q2T3RZchgKCREpFV8t28KoSRms3LKXq/unMOr87tSroYZ8yU4hISIlauf+PB77cBFvfL2Kdo1r8cYtAzilkxrylRUKCREpMZ8u2sDdk+azcdd+bjmtA789pxs1q1VOdFlyFBQSIlLstuzO5YGpWUz5fi3dmtfl+ev6cVzbBokuS46BQkJEio27817GOsakLWDX/jx+PbALvzizM9Wq6B3TZZVCQkSKxbod+7j33fnMWLiRvm0b8PjQPnRrUTfRZUmcFBIiEpfCQmf8nNU8+sFC8goLueeCHtx4agcq6x3T5YJCQkSO2crNexg5KYPZy7dycsfGPDa0N+0a1050WVKMFBIictQKCp2X/7WCP01fTNVKlXjsst5ceWJbtdQohxQSInJUFq/fxV0T5jEvZwcDezTnoUt60aJ+jUSXJSVEISEiRyQ3v4BnZy7j2VnZ1KtRlb9cfTwX9mmpo4dyTiEhIkX6btU2RkzMYMmG3Vx6fGvuvTCVRrWrJbosKQUKCRGJae+BfP708RJe/mIFLerV4OUbQpzdXQ35KhKFhIhE9WX2ZkZOymTV1r1cOyCFked1p64a8lU4CgkR+Q879uXx6AcLGT9nNR2a1Gb88JM4qWPjRJclCaKQEJF/m561gXvezWTTrlxuPaMjvxnYlRpV1ZCvIiuWhipm9jszczNrEtw3M3vazLLNLMPMToiYO8zMlga3YRHj/cwsM9jmadNLJkRKzebdudzxxrfc8no6DWtV493bT2XUeT0UEBL/kYSZtQXOAVZFDJ8HdAluA4DngAFm1gj4PRACHJhrZmnuvi2YMxyYTfgzsgejz7kWKVHuzrvfr+H+97LYm1vA/5zTlZ+f2YmqldWQT8KK43TTk8BdwJSIsSHA6x7+AO3ZZtbAzFoCZwLT3X0rgJlNBwab2Sygnrt/FYy/DlyCQkKkxKzdvo/RkzOZuXgTx6eEG/J1aa6GfPKf4goJM7sYWOPu8w45O9QaWB1xPycYO9x4TpRxESlmhYXOP75ZxdgPF1FQ6Nx3YSrDTmmvhnwSVZEhYWYzgBZRVo0G7gbOjbZZlDE/hvFYNQ0nfGqKlJSUWNNE5BArNu9hxMQMvlmxlR91bsKjl/WmbaNaiS5LkliRIeHuA6ONm1lvoANw8CiiDfCtmfUnfCTQNmJ6G2BtMH7mIeOzgvE2UebHqmkcMA4gFArFDBMRCcsvKOTFf63gyelLqF6lEo9f3ocr+rVRSw0p0jGfbnL3TKDZwftmthIIuftmM0sD7jCz8YQvXO9w93VmNg14xMwaBpudC4xy961mtsvMTgK+Bq4H/nKstYnI/8lau5O7Js5j/pqdDOrZnAeH9KJZPTXkkyNTUu+T+AA4H8gG9gI3AgRh8CAwJ5j3wMGL2MBtwKtATcIXrHXRWiQOufkF/PXTbJ6btYwGtary7LUncF6vFjp6kKNi4RcglV2hUMjT09MTXYZIUpn7Q7ghX/bG3Vx2QmvuvSCVhmrIJxHMbK67h4qap3dci5Qje3Lz+ePHi3n1y5W0ql+TV288kTO7NSt6Q5EYFBIi5cTnSzcxalImOdv2cf3J7bhrcHfqVNePuMRHzyCRMm7H3jwe/iCLt9Nz6NikNm/fejL9OzRKdFlSTigkRMqwj+av594p89m65wC/OLMTv/xxF/VbkmKlkBApgzbu2s+YtAV8kLme1Jb1eOWGE+nVun6iy5JySCEhUoa4O5O+XcMDU7PYl1fAnYO6Mfz0jmrIJyVGISFSRuRs28vdk+fz2ZJN9GvXkLFD+9C5WZ1ElyXlnEJCJMkVFjp///oHxn64CAfuv7gn153UjkpqyCelQCEhksSWbdrNyIkZzFm5jdO6NOGRS9WQT0qXQkIkCeUVFDLus+U89clSalatzB+v6MvQE1qrpYaUOoWESJKZv2YHIyZmsGDtTs7v3YIxF/ekWV015JPEUEiIJIn9eQU8/clS/vbZchrWqsbzPz2Bwb1aJrosqeAUEiJJIH3lVu6amMHyTXu4ol8b7rkglfq1qia6LBGFhEgi7c7N5w8fLeL12T/Qqn5NXr+pP6d3bZroskT+TSEhkiD/XLKJuydlsnbHPoad3J47B3WjthrySZLRM1KklG3fe4AHpmYx6ds1dGpam3duPZlQezXkk+SkkBApRR9kruO+KfPZtjePO87qzB1nd1ZDPklqcTV8MbMxZrbGzL4PbudHrBtlZtlmttjMBkWMDw7Gss1sZMR4BzP72syWmtlbZqaP0ZJyY+PO/fz8f+fyi398S4v6NUi741R+N6ibAkKSXnEcSTzp7n+MHDCzVOAqoCfQCphhZl2D1c8A5wA5wBwzS3P3LGBs8Fjjzex54GbguWKoTyRh3J135ubw0NQs9ucXMmJwd245rQNV1JBPyoiSOt00BBjv7rnACjPLBvoH67LdfTmAmY0HhpjZQuBs4JpgzmvAGBQSUoat3rqXuydn8vnSzfRv34jHhvamY1M15JOypThC4g4zux5IB/7H3bcBrYHZEXNygjGA1YeMDwAaA9vdPT/KfJEypaDQef2rlfxh2mIMeHBIT64doIZ8UjYVGRJmNgNoEWXVaMJ/6T8IePDvn4CbgGg/DU70ayB+mPmxahoODAdISUk5TPUipSt74y7umpDBt6u2c0bXpjxyWW9aN6iZ6LJEjlmRIeHuA4/kgczsBWBqcDcHaBuxug2wNliONr4ZaGBmVYKjicj50WoaB4wDCIVCMcNEpLTkFRTyt38u4+lPsqlVvTJP/KQvlx6vhnxS9sV1usnMWrr7uuDupcD8YDkNeMPMniB84boL8A3hI4YuZtYBWEP44vY17u5mNhO4HBgPDAOmxFObSGnJzNnBnRPmsWj9Li7o05IxF/Wkad3qiS5LpFjEe03icTM7jvCpoZXArQDuvsDM3gaygHzgdncvADCzO4BpQGXgZXdfEDzWCGC8mT0EfAe8FGdtIiVqf14Bf56xlBc+X07j2tX423X9GNQz2plZkbLL3Mv22ZpQKOTp6emJLkMqmK+Xb2HkpExWbN7DlaG23H1BD+rXVEM+KTvMbK67h4qap3dcixyFXfvzGPvRIv4+exVtG9XkHz8bwKmdmyS6LJESo5AQOUIzF21k9ORM1u3cz02nduB3g7pSq5p+hKR80zNcpAhb9xzgwalZTP5uDV2a1WHibadwQkrDRJclUioUEiIxuDvvZ67j91MWsGNfHr88uzO3n92Z6lXUb0kqDoWESBQbdu7nnnfnMz1rA33a1OfvPxtAj5b1El2WSKlTSIhEcHfeTl/NQ+8v5EB+IXef352bTlVDPqm4FBIigVVb9jJyUgZfLtvCgA6NGDu0D+2b1E50WSIJpZCQCq+g0HnlixX88ePFVKlUiYcv7cXVJ6aoIZ8ICgmp4JZsCDfk+371ds7u3oyHL+1Fy/pqyCdykEJCKqQD+YU8N2sZf525lDrVq/DUVcdxcd9WasgncgiFhFQ481ZvZ8TEDBat38VFfVsx5qJUGtdRQz6RaBQSUmHsO1DAkzOW8OLny2latzovXB/inNTmiS5LJKkpJKRC+GrZFkZOyuCHLXu5un8Ko87vTr0aasgnUhSFhJRrO/fn8egHi3jzm1W0a1yLN24ZwCmd1JBP5EgpJKTc+mThBkZPns/GXfu55bQO/PacbtSsppYaIkdDISHlzpbdudz/XhZp89bSrXldnr+uH8e1bZDoskTKJIWElBvuTtq8tdz/Xha79ufx64Fd+MWZnalWRS01RI6VQkLKhXU79nHP5Pl8smgjfds24PGhfejWom6iyxIp8+L+E8vM/p+ZLTazBWb2eMT4KDPLDtYNihgfHIxlm9nIiPEOZva1mS01s7fMrFq8tUn5V1jovPH1Ks594jO+WLaZey7owaTbTlFAiBSTuI4kzOwsYAjQx91zzaxZMJ4KXAX0BFoBM8ysa7DZM8A5QA4wx8zS3D0LGAs86e7jzex54GbguXjqk/Jt5eY9jJyUwezlWzm5Y2MeG9qbdo3VkE+kOMV7uuk24DF3zwVw943B+BBgfDC+wsyygf7Bumx3Xw5gZuOBIWa2EDgbuCaY8xowBoWERJFfUMjLX6zgTx8voVrlSjx2WW+uPLGtWmqIlIB4Q6IrcJqZPQzsB37n7nOA1sDsiHk5wRjA6kPGBwCNge3unh9l/n8xs+HAcICUlJQ4vwUpSxat38mICRnMy9nBwB7NeOiS3rSoXyPRZYmUW0WGhJnNAFpEWTU62L4hcBJwIvC2mXUEov1J50S/BuKHmR+Vu48DxgGEQqGY86T8yM0v4JmZy3h2Zjb1a1blL1cfz4V9WuroQaSEFRkS7j4w1jozuw2Y5O4OfGNmhUATwkcCbSOmtgHWBsvRxjcDDcysSnA0ETlfKrjvVm1jxMQMlmzYzSXHteK+i3rSqLZe1yBSGuJ9ddO7hK8lEFyYrkb4F34acJWZVTezDkAX4BtgDtAleCVTNcIXt9OCkJkJXB487jBgSpy1SRm390A+D07N4rLnvmTX/nxeviHEn686XgEhUorivSbxMvCymc0HDgDDgl/4C8zsbSALyAdud/cCADO7A5gGVAZedvcFwWONAMab2UPAd8BLcdYmZdgX2ZsZOSmD1Vv3ce2AFEae1526asgnUuos/Du97AqFQp6enp7oMqSY7NiXx6MfLGT8nNW0b1yLx4b24aSOjRNdlki5Y2Zz3T1U1Dy941qSxvSsDdzzbiabduVy6xkd+c3ArtSoqoZ8IomkkJCE27w7lzFpC5iasY7uLerywvUh+rRRQz6RZKCQkIRxd979fg33v5fF3twC/uecrtx6Ric15BNJIgoJSYg12/cxenImsxZv4viUcEO+Ls3Vb0kk2SgkpFQVFjr/+GYVj32wkEKH+y5MZdgp7alcSW+KE0lGCgkpNcs37WbkxEy+WbmVH3VuwqOX9aZto1qJLktEDkMhISUuv6CQF/+1gienL6FalUo8PrQPV4TaqKWGSBmgkJASlbV2J3dNnMf8NTs5N7U5D17Si+b11JBPpKxQSEiJyM0v4K+fZvPcrGU0qFWVZ645gfN7t9DRg0gZo5CQYjf3h63cNSGDZZv2cNkJrbn3glQaqt+SSJmkkJBisyc3nz9MW8xrX62kVf2avHrjiZzZrVmiyxKROCgkpFh8vnQToyZlkrNtH9ef3I67BnenTnU9vUTKOv0US1x27M3jofezeGduDh2b1ObtW0+mf4dGiS5LRIqJQkKO2Ufz13PvlPls3XOA287sxK9+3EUN+UTKGYWEHLWNu/YzJm0BH2SuJ7VlPV654UR6ta6f6LJEpAQoJOSIuTsTv13Dg1Oz2JdXwJ2DujH89I5UrayGfCLllUJCjkjOtr3cPXk+ny3ZRL92DRk7tA+dm9VJdFkiUsLi+hPQzN4ys++D20oz+z5i3SgzyzazxWY2KGJ8cDCWbWYjI8Y7mNnXZrY0eFy9sD4JFBY6r325knOf/Iz0lVu5/+KevHPryQoIkQoiriMJd7/y4LKZ/QnYESynAlcBPYFWwAwz6xpMfQY4B8gB5phZmrtnAWOBJ919vJk9D9wMPBdPfRKfZZt2M2JCBuk/bOO0Lk145FI15BOpaIrldJOFey38BDg7GBoCjHf3XGCFmWUD/YN12e6+PNhuPDDEzBYG214TzHkNGINCIiHyCgoZ99lynvpkKTWrVuaPV/Rl6Amt1VJDpAIqrmsSpwEb3H1pcL81MDtifU4wBrD6kPEBQGNgu7vnR5kvpWj+mh3cNSGDrHU7Oa9XC+4f0pNmddWQT6SiKjIkzGwG0CLKqtHuPiVYvhp4M3KzKPOd6NdA/DDzY9U0HBgOkJKSEmuaHIX9eQU89clSxn22nIa1qvHctSdwXu+WiS5LRBKsyJBw94GHW29mVYDLgH4RwzlA24j7bYC1wXK08c1AAzOrEhxNRM6PVtM4YBxAKBSKGSZyZOas3MqICRks37yHK/q14Z4LUqlfq2qiyxKRJFAcp5sGAovcPSdiLA14w8yeIHzhugvwDeEjhi5m1gFYQ/ji9jXu7mY2E7gcGA8MA6YgJWp3bj6Pf7SI17/6gdYNavL6Tf05vWvTRJclIkmkOELiKv7zVBPuvsDM3gaygHzgdncvADCzO4BpQGXgZXdfEGw2AhhvZg8B3wEvFUNtEsM/l2zi7kmZrN2xjxtOac+dg7pRWw35ROQQ5l62z9aEQiFPT09PdBllxva9B3hgahaTvl1Dp6a1GTu0D6H2asgnUtGY2Vx3DxU1T386VhDuzofz13PflPls35vHHWd15o6zO6shn4gclkKiAti4cz/3TpnPtAUb6NW6Hq/d1J+erdSQT0SKppAox9ydd+bm8NDULPbnFzJicHduOa0DVdSQT0SOkEKinFq9dS+jJmXyr+zN9G/fiMeG9qZjU/VbEpGjo5AoZwoKnde/WsnjHy2mksGDQ3py7YB2VKqklhoicvQUEuXI0g27GDExg29XbeeMrk155LLetG5QM9FliUgZppAoB/IKCnl+1jL+8mk2tapX5skr+3LJcWrIJyLxU0iUcZk5O7hzwjwWrd/FBX1acv/FPWlSp3qiyxKRckIhUUbtzyvgyRlLeOGz5TSpU52/XdePQT2j9WEUETl2Coky6OvlWxg5KZMVm/dwZagtd1/Qg/o11ZBPRIqfQqIM2bU/j7EfLeLvs1fRtlFN/vGzAZzauUmiyxKRckwhUUbMXLSRuydnsn7nfm7+UQf+59yu1Kqm/z4RKVn6LZPktu45wAPvLeDd79fSpVkdJt52CiekNEx0WSJSQSgkkpS7MzVjHWPSFrBjXx6//HEXbj+rE9WrqCGfiJQehUQS2rBzP6Mnz2fGwg30aVOfv/9sAD1a1kt0WSJSASkkkoi789ac1TzRWihGAAAMlUlEQVT8wUIO5Bdy9/nduelUNeQTkcRRSCSJVVv2MnJSBl8u28KADo0YO7QP7ZvUTnRZIlLBxfUnqpkdZ2azzex7M0s3s/7BuJnZ02aWbWYZZnZCxDbDzGxpcBsWMd7PzDKDbZ62CtJToqDQefHz5Zz753+SkbODhy/txZu3nKSAEJGkEO+RxOPA/e7+oZmdH9w/EzgP6BLcBgDPAQPMrBHweyAEODDXzNLcfVswZzgwG/gAGAx8GGd9SW3x+l3cNTGDeau3c3b3Zjx8aS9a1ldDPhFJHvGGhAMHr6jWB9YGy0OA1z38AdqzzayBmbUkHCDT3X0rgJlNBwab2Sygnrt/FYy/DlxCOQ2JA/mFPDsrm2dmZlO3RlWeuuo4Lu7bSg35RCTpxBsSvwammdkfCZ+6OiUYbw2sjpiXE4wdbjwnyni5M2/1du6akMHiDbu4uG8rfn9RKo3VkE9EklSRIWFmM4BoneNGAz8GfuPuE83sJ8BLwEAg2p/EfgzjsWoaTvjUFCkpKYetP1nsO1DAE9MX89K/VtCsbg1evD7EwNTmiS5LROSwigwJdx8Ya11wWuhXwd13gBeD5RygbcTUNoRPReUQPuUUOT4rGG8TZX6smsYB4wBCoVDMMEkWXy3bwshJGfywZS/XDEhh5HndqVdDDflEJPnF+wL8tcAZwfLZwNJgOQ24PniV00nADndfB0wDzjWzhmbWEDgXmBas22VmJwWvaroemBJnbQm3c38eoyZlcvULswF445YBPHJpbwWEiJQZ8V6TuAV4ysyqAPsJTgERfnXS+UA2sBe4EcDdt5rZg8CcYN4DBy9iA7cBrwI1CV+wLtMXrWdkbWD0u5ls2pXL8NM78puBXalZTS01RKRssfALkMquUCjk6enpiS7j37bszuX+97JIm7eWbs3rMvbyPhzXtkGiyxIR+Q9mNtfdQ0XN0zuui4m7kzZvLWPSFrA7N5/fDOzKbWd2oloVtdQQkbJLIVEM1u3Yxz2T5/PJoo30bduAx4f2oVuLuokuS0QkbgqJOBQWOm/OWcWjHywiv7CQey7owY2ndqByJb0pTkTKB4XEMVqxeQ8jJ2bw9YqtnNKpMY9d1oeUxrUSXZaISLFSSByl/IJCXv5iBX/6eAnVKlfisct6c+WJbdVSQ0TKJYXEUVi4bicjJmaQkbODgT2a89AlvWhRv0aiyxIRKTEKiSOQm1/AMzOX8ezMbOrXrMpfrj6eC/u01NGDiJR7CokifLtqGyMmZLB0424uPb41916YSqPa1RJdlohIqVBIxLD3QD5/+ngJL3+xghb1avDKDSdyVvdmiS5LRKRUKSSi+CJ7MyMnZbB66z5+elIKIwZ3p676LYlIBaSQiLBjXx6PvL+Qt9JX06FJbd4afhIDOjZOdFkiIgmjkAh8vGA997w7n827c7n1jHBDvhpV1ZBPRCq2Ch8Sm3blMua9BbyfsY7uLery4rAQfdqoIZ+ICFTgkHB33v1+Dfe/l8Xe3AJ+d25Xbj2jE1UrqyGfiMhBFTIk8goKGf56OjMXb+KElAY8fnkfOjdTQz4RkUNVyJCoWrkSHZvW4fSuTbn+5PZqyCciEkOFDAmAey9MTXQJIiJJTyfgRUQkprhCwsz6mtlXZpZpZu+ZWb2IdaPMLNvMFpvZoIjxwcFYtpmNjBjvYGZfm9lSM3vLzNT7QkQkweI9kngRGOnuvYHJwJ0AZpYKXAX0BAYDz5pZZTOrDDwDnAekAlcHcwHGAk+6exdgG3BznLWJiEic4g2JbsBnwfJ0YGiwPAQY7+657r4CyAb6B7dsd1/u7geA8cAQC7dTPRuYEGz/GnBJnLWJiEic4g2J+cDFwfIVQNtguTWwOmJeTjAWa7wxsN3d8w8ZFxGRBCry1U1mNgNoEWXVaOAm4Gkzuw9IAw4c3CzKfCd6KPlh5seqaTgwHCAlJSVm7SIiEp8iQ8LdBxYx5VwAM+sKXBCM5fB/RxUAbYC1wXK08c1AAzOrEhxNRM6PVtM4YBxAKBSKGSYiIhKfeF/d1Cz4txJwD/B8sCoNuMrMqptZB6AL8A0wB+gSvJKpGuGL22nu7sBM4PJg+2HAlHhqExGR+Fn49/Mxbmz2K+D24O4kYFTwCx8zO3g6Kh/4tbt/GIyfD/wZqAy87O4PB+MdCV/IbgR8B/zU3XOPoIZNwA/H/E2UjiaEj5bKGtVduspq3VB2a6/Idbdz96ZFTYorJOTImFm6u4cSXcfRUt2lq6zWDWW3dtVdNL3jWkREYlJIiIhITAqJ0jEu0QUcI9Vduspq3VB2a1fdRdA1CRERiUlHEiIiEpNCopiY2RVmtsDMCs0s5qsOkq0Lrpk1MrPpwdedbmYNo8w5y8y+j7jtN7NLgnWvmtmKiHXHJUvdwbyCiNrSIsaTeX8fF3RXXmBmGWZ2ZcS6Ut3fsZ6vEeurB/svO9if7SPWRe0EXRqOoO7fmllWsH8/MbN2EeuiPmdKyxHUfoOZbYqo8WcR64YFz62lZjasWApyd92K4Qb0INzwcBYQijGnMrAM6AhUA+YBqcG6t4GrguXngdtKqe7HCXfyBRgJjC1ifiNgK1AruP8qcHkC9vcR1Q3sjjGetPsb6Ap0CZZbAeuABqW9vw/3fI2Y8wvg+WD5KuCtYDk1mF8d6BA8TuUkqvusiOfwbQfrPtxzJolqvwH4a5RtGwHLg38bBssN461JRxLFxN0XuvviIqYlYxfcIcHXO9KveznwobvvLdGqina0df9bsu9vd1/i7kuD5bXARqDINz2VgKjP10PmRH4/E4AfB/s3VifopKjb3WdGPIdnE24FlAyOZJ/HMgiY7u5b3X0b4c7cg+MtSCFRupKxC25zd18HEPzbrIj5VwFvHjL2cHDY/qSZVS+JIqM40rprmFm6mc0+eIqMMrS/zaw/4b8ol0UMl9b+jvV8jTon2J87CO/fI9m2pBzt174Z+DDifrTnTGk50tqHBs+BCWZWVPftuFTYz7g+FnaYjrjufiS9pmJ1uz2qLrhH63B1H+XjtAR6A9MihkcB6wn/IhsHjAAeOLZK/+vrFUfdKe6+1sJtXz41s0xgZ5R5ybq//xcY5u6FwXCJ7e9oJUQZO3Q/JeQ5XYQj/tpm9lMgBJwRMfxfzxl3XxZt+xJwJLW/B7zp7rlm9nPCR3JnH+G2R00hcRS86I64RYnVHfeouuAercPVbWYbzKylu68LfiltPMxD/QSY7O55EY+9LljMNbNXgN8VS9EUT93B6RrcfbmZzQKOByaS5Pvbwh8F/D5wj7vPjnjsEtvfURyum/Ohc3LMrApQn/A1qyPZtqQc0dc2s4GEg/sMj+gTF+M5U1ohUWTt7r4l4u4LhD/V8+C2Zx6y7ax4C9LpptKVjF1w04KvdyRf92oOOdUU/KI7eJ7/EsIfRFUaiqzbzBoePB1jZk2AU4GsZN/fwXNjMvC6u79zyLrS3N9Rn6+HzIn8fi4HPg32b6xO0KWhyLrN7Hjgb8DF7r4xYjzqc6aU6oYjq71lxN2LgYXB8jTg3OB7aEj4Yxwij/qPTaKu4pe3G3Ap4STPBTYA04LxVsAHEfPOB5YQ/stkdMR4R8I/RNnAO0D1Uqq7MfAJsDT4t1EwHgJejJjXHlgDVDpk+0+BTMK/rP4O1EmWuoFTgtrmBf/eXBb2N/BTIA/4PuJ2XCL2d7TnK+HTWxcHyzWC/Zcd7M+OEduODrZbDJxXGvv3KOqeEfycHty/aUU9Z5Ko9keBBUGNM4HuEdveFPxfZAM3Fkc9ese1iIjEpNNNIiISk0JCRERiUkiIiEhMCgkREYlJISEiUsrM7GUz22hmxfISZjP7yMy2m9nUQ8ZfMrN5Ee/OrnO0j62QEBEpfa9SDH2VIvwBuC7K+G/cva+79wFWAXcc7QMrJERESpm7f0b4nen/ZmadgiOCuWb2uZl1P4rH+wTYFWV8Z/DYBtTkGNp0KCRERJLDOOD/uXs/wu1Wni2OBw3at6wHugN/Odrt1btJRCTBgmsFpwDvhP/oB8KfxYGZXUb0Jo5r3L3ID3Ny9xvNrDLhgLgSeOVoalNIiIgkXiXC7ev/65MG3X0SMCmeB3f3AjN7C7iTowwJnW4SEUmw4NrBCjO7AsLXEMysbzyPGTxG54PLwEXAoqN+HPVuEhEpXWb2JuG23k0INxr8PeHmjc8BLYGqhD/Z74g+K8TMPid8zaEOsIXwBylNBz4H6hH+rIl5hD+mN9rnqcR+bIWEiIjEotNNIiISk0JCRERiUkiIiEhMCgkREYlJISEiIjEpJEREJCaFhIiIxKSQEBGRmP4/XAnv3EJeW58AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "e = np.zeros(ny-1)\n",
    "\n",
    "# TODO: Write here a for-loop that loops over the grid points\n",
    "# and calculate the shear strain rate on every step\n",
    ".\n",
    ".\n",
    ".\n",
    " \n",
    "plt.plot(e, -(y[1:] + y[:-1]) / 2.0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Velocity boundary conditions\n",
    "\n",
    "Like in the case of heat equation, the FD formulation for velocity field may have different boundary conditions:\n",
    "\n",
    " 1. No-slip:\n",
    "     - Velocity parallel to the surface is fixed at zero\n",
    "     - $v_x = 0$\n",
    " 2. Free-slip: \n",
    "     - Shear stress at the surface is zero\n",
    "     - $0 = \\sigma_{xy} = \\eta\\dot\\epsilon_{xy} = \\eta\\frac{\\partial v_x}{\\partial y} \\Rightarrow \\frac{\\partial v_x}{\\partial y} = 0$\n",
    " 3. Other:\n",
    "     - Velocity is fixed at non-zero value\n",
    "     - e.g. $v_x = a \\neq 0$, also incoming/outgoing flow ($v_y \\neq 0$) possible (in 2D models)\n",
    " \n",
    " \n",
    " We used cases 1 and 3 previously.\n",
    " \n",
    " ### Do\n",
    " \n",
    " - Discretize the free-slip boundary condition $\\frac{\\partial v_x}{\\partial y} = 0$\n",
    " - If you want to apply this at the upper surface (i.e. at $x=x_0$), how would your system of equations change? How would the coefficients in the matrix $A$ change?\n",
    " - Modify your script for the channel flow model, and change the upper surface velocity boundary condition from $v_x=0$ to a free-slip B.C.\n",
    "     - How does the resulting velocity plot change?\n",
    "     - How does the plot for shear strain rate change?"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}