{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# A time-dependent problem, Burgers' equation\n", "\n", "We will solve the viscous Burgers equation, a nonlinear equation for the advection and diffusion on momentum in one dimension.\n", "\n", "$$\n", "\\frac{\\partial u}{\\partial t} + u \\frac{\\partial u}{\\partial x} - \\nu \\frac{\\partial^2 u}{\\partial x^2} = 0\n", "$$\n", "\n", "We will solve on a periodic interval mesh, and therefore do not impose any boundary conditions. As usual, we need to derive a variational form.\n", "\n", "## Spatial discretisation\n", "\n", "We first discretise in space, mulitplying by a test function $v \\in V$ and integrating the viscosity term by parts to obtain the semi-discrete problem. Find $u(x, t) \\in V$ such that\n", "\n", "$$\n", "\\int_\\Omega \\frac{\\partial u}{\\partial t} v + u \\frac{\\partial u}{\\partial x} v + \\nu \\frac{\\partial u}{\\partial x}\\frac{\\partial v}{\\partial x} \\, \\mathrm{d}x = 0 \\quad \\forall v \\in V.\n", "$$\n", "\n", "## Time discretisation\n", "We now need to discretise in time. For simplicity, and stability we'll use backward Euler, replacing all instances of $u$ with $u^{n+1}$ and the time derivative by $\\frac{u^{n+1} - u^n}{\\Delta t}$. We end up with the discrete problem, find $u^{n+1} \\in V$ such that\n", "\n", "$$\n", "\\int_\\Omega \\frac{u^{n+1} - u^n}{\\Delta t} v + u^{n+1} \\frac{\\partial u^{n+1}}{\\partial x} v + \\nu \\frac{\\partial u^{n+1}}{\\partial x}\\frac{\\partial v}{\\partial x} \\, \\mathrm{d}x = 0 \\quad \\forall v \\in V.\n", "$$\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Implementation\n", "\n", "To solve the problem in a concrete setting, we need two things. A domain, and an initial condition for $u$. For the former, we'll choose a periodic interval of length 2, for the latter, we'll start with $u = \\sin(2 \\pi x)$.\n", "\n", "In addition we need to choose the viscosity, which we will set to a small constant value $\\nu = 10^{-2}$\n", "\n", "As ever, we begin by importing Firedrake" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "from firedrake import *\n", "\n", "n = 100\n", "mesh = PeriodicIntervalMesh(n, length=2)\n", "\n", "x = SpatialCoordinate(mesh)[0]\n", "\n", "u_init = sin(2*pi*x)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "nu = Constant(1e-2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We choose degree 2 piecewise continuous Lagrange polynomials for our solution and test space:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "V = FunctionSpace(mesh, \"Lagrange\", 2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We also need solution functions for $u^{n+1}$ and $u^n$, along with a test function $v$." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "u_n1 = Function(V, name=\"u^{n+1}\")\n", "u_n = Function(V, name=\"u^{n}\")\n", "v = TestFunction(V)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We provide the initial condition for $u_n$, and choose a $\\Delta t$ such that the advective Courant number is around 1. This is more restrictive than required for stability of the time integration, but gives us enough accuracy to see the temporal evolution of the system." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "u_n.interpolate(u_init)\n", "dt = 1.0 / n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we're ready to define the variational form. Since this problem is nonlinear, note that we do not have a trial function anywhere. We just write down the residual, Firedrake will automatically compute the Jacobian by differentiating the residual inside the nonlinear solver." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "F = (((u_n1 - u_n)/dt) * v +\n", " u_n1 * u_n1.dx(0) * v + \n", " nu*u_n1.dx(0)*v.dx(0))*dx" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For visualisation purposes, we will save a copy of the state $u_n$ at each timestep, we can plot and animate these in the notebook if the `ipywidgets` package is installed." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# If passed an existing Function object, the Function \n", "# constructor makes a copy.\n", "results = [Function(u_n)]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we loop over the timesteps, solving the equation and advancing in time.\n" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "t = 0.0\n", "t_end = 0.5\n", "while t <= t_end:\n", " solve(F == 0, u_n1)\n", " u_n.assign(u_n1)\n", " results.append(Function(u_n))\n", " t += dt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To visualize the results, we'll create a movie using matplotlib's animation tools.\n", "First, we'll create a figure and axes to draw on and plot the initial values." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD8CAYAAAB6paOMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXRV9bn4//dzEpKQQEyYB5FBhiAYDAQIYkVFEVFBK1qxWlunX+319tfae3t11eX9Xqvr2q+rV1t/drDVVmtvRakDTkVkFg0YxjA0ECaZiYAQIAlJzvP74+wdjjEhwxn2GZ7XWmflnD0+Z+ezz7M/+7P3Z4uqYowxJnn5vA7AGGOMtywRGGNMkrNEYIwxSc4SgTHGJDlLBMYYk+QsERhjTJILSyIQkRdF5JCIbGhmvIjIr0WkXETWi8jooHF3ishW53VnOOIxxhjTeuGqEfwZmHqW8dcAQ5zXfcBvAUSkC/CfwHhgHPCfIpIbppiMMca0QlgSgaouBY6cZZIZwMsaUAzkiEhv4GpgvqoeUdWjwHzOnlCMMcaEWWqU1tMX2B30eY8zrLnhXyMi9xGoTZCVlTUmLy8vMpFGQX19PXv27OHw4cM0dWf3Oeecw3nnnUdaWpoH0Zl4p6ocOHCAAwcO4Pf7vzY+MzOT/v37k5mZ6UF0xkurVq36QlW7Nx4erUQQMlV9HngeoLCwUEtKSjyOqH02b97MzTffzBdffIHP56OgoIBZs2aRlpbG22+/zSeffMKxY8fYt28fb7zxBpdeeqnXIZs4cuDAAe666y5Wr14NwLBhw/jWt75Fz549mT9/PkuXLuXIkSNs3bqVF198kTvuuMPjiE00iciuJkeoalhewABgQzPjfg/MCvpcBvQGZgG/b2665l5jxozReFRSUqK5ubkK6IgRI3Tjxo1fm2b//v167bXXKqCpqak6Z84cDyI18Wjv3r06bNgwBTQ7O1s/+uijr01z8uRJ/dd//VcFFNDHHnvMg0iNV4ASbeo3uqmB7Xm1kAiuBT4ABCgCVjrDuwA7gFzntQPo0tK64jERbNu2rSEJzJgxQ0+ePNnstH6/X3/yk58ooCkpKbpw4cIoRmri0fHjx3XEiBEKaH5+vu7fv/+s07/44osqIgroc889F6UojdcimgiAvwH7gVoC5/nvBr4PfN8ZL8BzwDagFCgMmvcuoNx5fa8164u3RHDixAm98MILFdBp06ZpTU1Ni/P4/X598MEHFdBu3brpzp07oxCpiUd+v1+/+c1vKqDDhw/XL774olXz/fnPf26oeS5dujTCUZpYEPEaQTRf8ZYI3Kr40KFD9csvv2z1fHV1dXr11VcroFdddZXW1dVFMEoTr373u99pamqqdu7cWcvKyto0r3uwMWTIED169GiEIjSxwhKBR+bPn99w1FVSUtLm+Q8dOqTdunVTQH/1q19FIEITz3bs2KFZWVkK6CuvvNLm+Wtra7WgoEABveuuuyIQoYkllgg8cOrUKR0wYIAC+vjjj7d7OW+//bYCmpmZaaeITAO/399QY7zlllvavZzNmzdrenq6Ajp//vwwRmhiTXOJwPoaiqCnn36anTt3kp+fz09/+tN2L2f69OnccsstnDp1iocffjiMEZp49v777zNv3jxycnJ49tln272cvLw8Hn30UQAefPBB6uvrwxWiiRdNZYdYf8VDjeDAgQPaqVMnBZq8jK+tdu7c2XDUVlxcHIYITTw7ffq05uXlKaD/8z//E/LyqqqqtH///groH/7whzBEaGIRViOIrieffJITJ05w3XXXMXny5JCX179/fx588EEAHnnkkZCXZ+LbX/7yF8rKyhg0aBD/8i//EvLyMjIyePLJJwF4/PHHqaqqCnmZJo40lR1i/RXrNYK9e/c2HL2vXbs2bMs9cuSIZmdnK6CLFy8O23JNfDl9+nRD29PLL7/c8gyt5Pf7NT8/XwH99a9/HbblmtiB1Qii56mnnqKmpoaZM2cyatSosC03Nze3oVbw85//PGzLNfHllVdeYefOnQwbNozbbrstbMsVEf7rv/4LgF/84hfU1NSEbdkmtlkiCLOKigp+//vfA5E5hfPDH/6Q7OxsFixYQHFxcdiXb2Kb3+/nv//7v4FA+UpJSQnr8qdPn05+fj579+7lpZdeCuuyTeyyRBBmzz//PFVVVVx77bVhrQ24cnNzuf/++wF45plnwr58E9vee+89tm7dSv/+/bn11lvDvnyfz8dDDz0EBMpX4GyCSXSWCMLo9OnTPPfccwD8+Mc/jth6HnjgAVJSUpgzZw67d+9ueQaTMJ5++mkgUDNMTY1M58EzZ86kb9++bN68mQ8//DAi6zCxxRJBGM2ZM4f9+/czcuRIrrjiioit59xzz+Xmm2+mvr6e3/72txFbj4ktGzZsYNGiRXTq1Im77747Yuvp0KEDDzzwAAC//vWvI7YeEzssEYTRb37zGyBwxC4iEV3XD3/4QwBeeOEFa9RLEm75uv322znnnHMiuq57772X1NRU5s2bx/bt2yO6LuM9SwRhsmHDBpYvX052dja33357xNdXVFTEqFGjOHToEG+99VbE12e8dfLkSV555RWAhqP1SOratSu33XYb9fX1vPDCCxFfn/GWJYIw+cMf/gDAbbfdRlZWVsTXJyLce++9QKCB2iS22bNnU1lZyYQJExgxYkRU1umWrxdffJHa2tqorNN4wxJBGFRXV/OXv/wFgPvuuy9q6/32t79Nx44dWbhwIeXl5VFbr4k+90AjmuVr4sSJDB8+nAMHDvDuu+9Gbb0m+iwRhMHcuXM5evQoo0ePpqCgIGrrzcnJ4eabbwbg5Zdfjtp6TXSVlZVRXFxM586dG/7f0SAi3HPPPQB2T0GCs0QQBm5t4Dvf+U7U1+2u85VXXrFrvhOUW75mzpwZldOOwW677TZ8Ph/vv/8+hw8fjuq6TfRYIghRRUUF//jHP0hJSWHWrFlRX/9ll11G37592bFjB8uXL4/6+k1k+f3+hkbiO+64I+rr79WrF1OmTKG2tpbZs2dHff0mOsKSCERkqoiUiUi5iDzUxPinRWSt89oiIl8GjasPGjc3HPFE0+zZs6mrq+Pqq6+mR48eUV9/SkpKw1VK7pGjSRwff/wxu3btol+/fkyaNMmTGNwEZOUrcYWcCEQkhcCD6a8BLgBmicgFwdOo6o9V9SJVvQh4FngjaHSVO05Vp4caT7S9+uqrQKDh1ituInjjjTfs6o4E45avWbNm4fN5U4G/4YYb6NixIytXrrR7ChJUOErWOKBcVber6mngVWDGWaafBfwtDOv13Oeff87y5cvp2LEjM2ac7StH1siRIxkxYgRffPEFH330kWdxmPCqq6tjzpw5AGHtZbStMjMzufHGG/H7/bz++uuexWEiJxyJoC8Q3OHNHmfY14hIf2AgsDBocIaIlIhIsYjcEIZ4osbdSa+77rqoN+I19q1vfQuA1157zdM4TPgsXryYiooKhg0bRn5+vqex3HLLLYCVr0QV7brmrcAcVQ1+KGp/VS0EbgOeEZHzm5pRRO5zEkZJRUVFNGJtkdt45u4kXnIvK3zzzTc5ffq0x9GYcHB/dG+55ZaId1nSkquvvprOnTuzevVqu2clAYUjEewF+gV9PtcZ1pRbaXRaSFX3On+3A4uBJi/EV9XnVbVQVQu7d+8easwh27lzJytXriQrK4tp06Z5HQ55eXnk5+dz7Ngx5s+f73U4JkR1dXW88UagKS0WDjQyMjIaTn9arSDxhCMRfAYMEZGBIpJG4Mf+a1f/iEgekAt8GjQsV0TSnffdgInApjDEFHFu/z7Tpk0jMzPT42gCbrrpJiBQKzDxbdmyZRw+fJihQ4dGrUuJlsycOROw8pWIQk4EqloHPADMAzYDr6nqRhF5TESCrwK6FXhVv3rX03CgRETWAYuAJ1U1LhKBuzPceOONHkdyhhvL3Llzqa+vb2FqE8uCy5fXp4VcU6ZMITMzk5KSEnsORqJp6kHGsf7y+uH1hw4dUp/Ppx06dNAvv/zS01iC+f1+Pf/88xXQpUuXeh2OaSe/36/9+vVTQIuLi70O5yu++c1vKqDPPvus16GYdsAeXh8+77zzDn6/n8mTJ0e8X/i2EJGGWoFV3+PXqlWr2L17N3369GHs2LFeh/MVVr4SkyWCdpg7N9AE4uW9A8254YbAFbhz5861vofilFu+rr/+es9uImvOtddei8/nY9myZRw9etTrcEyYxFYpiwNVVVUNV+VMnx57N0IXFRXRvXt3tm3bxubNm70Ox7TDO++8A8RW+5MrNzeXSZMmUVtby7x587wOx4SJJYI2Wrx4MadOnWLMmDH06dPH63C+JiUlpeFyVutDPv7s2bOHtWvXkpWV5VnfQi257rrrACtficQSQRu5hd/dGWKR7ajx67333gPgqquuIiMjw+NomuaWrw8++IC6ujqPozHhYImgDVQ1LhLBlClTSE1NZfny5Rw5csTrcEwbxEP5Gjp0KEOHDuXIkSMUFxd7HY4JA0sEbbBhwwY+//xzevXqxejRo70Op1nZ2dlMmjQJv9/PBx984HU4ppWqqqpYsGABQEzcrX42bqJy2zNMfLNE0Abuj+rUqVNj7mqOxtwfkn/84x8eR2Jaa8mSJVRVVVFQUEDv3r29DuesrHwlltj+NYsxbqG/5pprPI6kZVOnTgVg3rx5+P1+j6MxrRFP5euSSy4hKyuL9evXs2/fPq/DMSGyRNBKlZWVfPzxx/h8Pq666iqvw2nR8OHDOe+886ioqGDNmjVeh2NaIZ4SQXp6OldccQWAXUaaACwRtNKiRYuora2lqKiI3Nxcr8NpkYg01Aqs+h77duzYQVlZGeeccw5FRUVeh9MqVr4ShyWCVnILu1v444HtqPHDPaq+8sorSU1N9Tia1nHL1/z58+0y0jhniaCV3B11ypQpHkfSeldccQU+n49PP/2UY8eOeR2OOYt4LF+DBg1i0KBBHD16lJKSEq/DMSGwRNAK27dvZ/v27eTm5lJYWOh1OK3mnmaor69n8eLFXodjmlFXV8fChYGnt8ZTIoDAk8sAexhSnLNE0ApuIb/iiitISUnxOJq2cX9YbEeNXZ999hnHjx9nyJAhDBgwwOtw2sS9cMLKV3yzRNAKbiGPh6uFGrMdNfbFc/m6/PLLG04/VlZWeh2OaSdLBC2or69vqLbH4446btw4srOz2bJlC59//rnX4ZgmxHMiyMnJYdy4cdTV1bFkyRKvwzHtZImgBatWreLo0aMNDWPxJjU1lcsvvxywWkEsqqyspLi4mJSUlIb/U7yxWmf8s0TQArc2MHnyZI8jab8rr7wSOPNdTOxYtmwZdXV1FBYWxtTT7trCylf8C0siEJGpIlImIuUi8lAT478rIhUistZ53RM07k4R2eq87gxHPOHkdgIWz4nAvQN04cKF9tSyGJMI5Wv8+PF07NiRDRs2cPDgQa/DMe0QciIQkRTgOeAa4AJglohc0MSks1X1Iuf1R2feLsB/AuOBccB/ikjM3LZbU1PDxx9/DBC31XYIdDfRq1cvDhw4wD//+U+vwzFBEqHGmZ6ezje+8Q0gcAe+iT/hqBGMA8pVdbuqngZeBVr7MN+rgfmqekRVjwLzgZi5dbe4uJjq6mouvPBCevTo4XU47SYiX6kVmNhw+PBh1q5dS3p6OhMmTPA6nJBY+Ypv4UgEfYHdQZ/3OMMau0lE1ovIHBHp18Z5EZH7RKREREoqKirCEHbL3ELtFvJ45n4H91SE8Z57k9/FF19Mx44dvQ0mRJYI4lu0GovfAQaoaj6Bo/6X2roAVX1eVQtVtbB79+5hD7ApbjU3nk8LudwddcmSJdTX13scjYHEKl+jR48mOzubbdu2sWvXLq/DMW0UjkSwF+gX9PlcZ1gDVT2sqjXOxz8CY1o7r1dOnjxJcXExPp8vZh8i3hYDBw6kf//+HDlyhNLSUq/DMZxJBIlQ40xJSeGyyy4DsO5M4lA4EsFnwBARGSgiacCtwNzgCUQk+HFL04HNzvt5wBQRyXUaiac4wzz36aefUltbS0FBATk5OV6HExa2o8aOQ4cOsWnTJjp27MjYsWO9DicsrHzFr5ATgarWAQ8Q+AHfDLymqhtF5DERme5M9kMR2Sgi64AfAt915j0C/JxAMvkMeMwZ5jm3MLuFOxHYjho73LtwJ06cSFpamsfRhIeVr/gVlo7PVfV94P1Gwx4Nev8w8HAz874IvBiOOMIpkRPB0qVL8fv9Mf/c5USWiOUrPz+fnJwcdu7cya5du+jfv7/XIZlWsl+CJpw6dYqVK1fi8/m45JJLvA4nbAYMGMCAAQM4evQo69ev9zqcpJaIiSAlJYVLL70UwPodijOWCJqQiO0DLveHx2788Y7bPpCZmZkw7QMuK1/xyRJBE9yjGffoJpG4V0AtXbrU40iSl7vtJ0yYkDDtAy4rX/HJEkET3EKcCJeNNuYmt2XLluH3+z2OJjklcvkaNWoUnTt3Zvv27ezZs8frcEwrWSJopKamhuLiYoCEah9wDRw4kL59+3L48GE2b97c8gwm7NxEkIg1zpSUlIb9ZtmyZR5HY1rLEkEjn332GTU1NVx44YV07drV63DCTkQafoCs+h59bkN9Wloa48aN8zqciLDTQ/HHEkEjiXy05rJE4J3ly5ejqowbNy7u+xdqjpWv+GOJoJFkSwT2fILoSobyNWbMGDp27MimTZuIVgeRJjSWCILU1dWxfPlyIDHbB1zDhw+nW7du7Nu3j+3bt3sdTlJxE4Hbf38iSktLo6ioCKDheR4mtlkiCLJu3TpOnDjB+eefT58+fbwOJ2JEpCHR2Y4aPadOnWLVqlX4fD4uvvhir8OJKDfRWfmKD5YIgrhXOSRybcBlV3ZE34oVK6irq2PUqFFkZ2d7HU5EWfmKL5YIgrhHL4lcbXe539F21Ohxt3UylK8JEyaQkpLC6tWrOXHihNfhmBZYInCoalLtqAUFBWRmZrJlyxYOHTrkdThJwT3QSIYaZ6dOnSgoKKC+vp4VK1Z4HY5pgSUCR3l5OYcOHaJHjx4MGTLE63AirkOHDtagF0V1dXV8+umnQHIkArBaZzyxROAIPloTEY+jiQ47jxs9wRci9O7du+UZEoCVr/hhicCRTNV2l/td3UtmTeQkY/maOHEicKaR3MQuSwQO98fQLbzJoKioCJ/Px5o1azh16pTX4SS0ZLg/pbGePXsyePBgTp48ybp167wOx5yFJQKgoqKCsrIyOnbsSEFBgdfhRE3nzp0ZNWoUdXV1rFy50utwEpaqJuWBBmD3q8SJsCQCEZkqImUiUi4iDzUx/kER2SQi60VkgYj0DxpXLyJrndfcxvNGwyeffALA+PHj6dChgxcheMZ21MjbuXMn+/bto2vXruTl5XkdTlS5ic9OP8a2kBOBiKQAzwHXABcAs0TkgkaTrQEKVTUfmAP836BxVap6kfOajgfcH8FkO1oD21Gjwd22F198cdJciOByy9fHH39s/VrFsHDUCMYB5aq6XVVPA68CM4InUNVFquqehC4Gzg3DesPGrREk0/lbl7ujfvrpp9TX13scTWJK1tNCAHl5eXTt2pX9+/ezc+dOr8MxzQhHIugL7A76vMcZ1py7gQ+CPmeISImIFIvIDc3NJCL3OdOVhLNHw+rqakpKSgAarqtPJueeey79+vXj2LFjbNq0yetwElJwjSDZiAgTJkwAzhxwmdgT1cZiEbkdKASeChrcX1ULgduAZ0Tk/KbmVdXnVbVQVQu7d+8etphWrVrF6dOnGTlyZMI9qL617PRQ5Hz55Zds2LCBDh06UFhY6HU4nrDyFfvCkQj2Av2CPp/rDPsKEbkS+BkwXVVr3OGqutf5ux1YDET1sh33KCUZq+0u97vbEVv4rVixAlVt6KM/GVn5in3hSASfAUNEZKCIpAG3Al+5+kdECoDfE0gCh4KG54pIuvO+GzARiOr5iWQ+f+uyHTVy7EADCgsLSU1NpbS0lOPHj3sdjmlCyIlAVeuAB4B5wGbgNVXdKCKPiYh7FdBTQCfg9UaXiQ4HSkRkHbAIeFJVo5YIVLVhR03G87euCy+8kKysLLZt28bBgwe9DiehJHP7gKtjx46MHj0av99vHdDFqLC0Eajq+6o6VFXPV9UnnGGPqupc5/2Vqtqz8WWiqvqJql6oqqOcvy+EI57WKi8vp6Kigp49ezJo0KBorjqmpKamNjSUW60gfOrq6hp++JI5EYC1E8S6pL6zOLg2kGzXdzfm/lDZjho+GzZs4MSJEwwaNIhevXp5HY6nrHzFNksE2NEanNkGblfJJnRWvs5wt8GKFSvsfpUYZIkA21HhzD0UJSUl1NTUtDC1aQ0rX2f06dOH/v37U1lZycaNG70OxzSStIng2LFjbNy4kbS0NEaPHu11OJ7LyclhxIgRnD59mtWrV3sdTkKwRPBVVuuMXUmbCNzru0ePHk1GRobX4cQEuwM0fA4cOMCOHTvo1KkTI0eO9DqcmOAmAitfsSdpE4EdrX2dHbGFj7sNi4qKSElJ8Tia2GCJIHYlbSJwd1RLBGcEX9lhPUWGxv2xc2tZBvLz88nMzGy4bNvEjqRMBPX19RQXFwO2owYbOnQoXbp04cCBA+zatcvrcOKaJYKvS01NZezYsYDVCmJNUiaCTZs2cfz4cfr370+fPn28DidmBPcUaaeH2q+mpoZVq1YBydmj7dnY6cfYlJSJwC2EdrT2ddZgHLo1a9ZQU1PD8OHDyc3N9TqcmGKJIDYlZSKwanvzrEYQOjvQaJ5bQ/rss8+ora31OBrjSspEYA3FzRs3bhw+n4+1a9dy8uRJr8OJS1a+mtetWzeGDBlCVVUV69at8zoc40i6RHD48GG2bNlCx44dGTVqlNfhxJxOnTqRn59PfX19w5PbTNtYjfPsrNYZe5IuEbhXCxUWFtKhQwePo4lNtqO23+7du9m7dy85OTnk5eV5HU5MsvIVe5IuEdj525ZZg177udts/Pjx+HxJt3u1it1YFnuSrqTaHcUtCz5isxvL2sbaB1o2YsQIOnfuzK5du9i/f7/X4RiSLBHU1dWxcuVKwGoEZzNo0CC6d+9ORUUF27dv9zqcuGI1zpalpKQwbtw4wGqdsSKpEsGGDRs4efIkgwYNokePHl6HE7PsxrL2qa6uZvXq1YgI48eP9zqcmGblK7aEJRGIyFQRKRORchF5qInx6SIy2xm/QkQGBI172BleJiJXhyOe5tjRWuvZjtp2q1evpra2lhEjRpCdne11ODHNyldsCTkRiEgK8BxwDXABMEtELmg02d3AUVUdDDwN/MKZ9wLgVmAEMBX4jbO8iLBE0Hq2o7adla/WC34Q0unTpz2OxoSjRjAOKFfV7ap6GngVmNFomhnAS877OcBkCTwkeAbwqqrWqOoOoNxZXkSsX78esB21NQoLC0lJSWH9+vWcOHHC63DigjUUt16XLl3Iy8ujpqaGtWvXeh1OXDhx4kTELt4IRyLoC+wO+rzHGdbkNKpaBxwDurZyXgBE5D4RKRGRkvZ2Ybtq1SrWrVtHfn5+u+ZPJllZWYwaNcpuLGslVbUaQRtZrbNtbr75Zvr06cPy5cvDvuy4aSxW1edVtVBVC7t3796uZaSkpJCfn09qamqYo0tMtqO23u7du9m3bx9dunRh6NChXocTF6x8tZ7f76e4uJgDBw5w3nnnhX354UgEe4F+QZ/PdYY1OY2IpALnAIdbOa/xiPVE2nruNioqKiJw1tO0xBJB65WVlfHll1/St29f+vXr1/IMbRSORPAZMEREBopIGoHG37mNppkL3Om8nwks1MDJrrnArc5VRQOBIcDKMMRkwsDdUYuLi+3GshYEP5rStM4FF1xAdnY2n3/+Ofv27fM6nJgW6dOOIScC55z/A8A8YDPwmqpuFJHHRGS6M9kLQFcRKQceBB5y5t0IvAZsAv4B/Iuq1ocakwmPgQMH0qNHD7744gvKy8u9DiemWftA2/l8vob7LaxWcHYxnwgAVPV9VR2qquer6hPOsEdVda7zvlpVb1bVwao6TlW3B837hDPfMFX9IBzxmPCwG8tap6qqijVr1iAiDXfMmtax8tU6cZEITOKyHbVlq1atoq6ujpEjR9qNZG1k5atlx44dY9OmTaSlpTF69OiIrMMSgTkr64m0ZXZaqP3cNpVVq1ZRU1PjcTSxacWKFagqo0ePJj09PSLrsERgzqqwsJDU1FRKS0uprKz0OpyYZD3atl9OTg4XXHABNTU1rFmzxutwYlI0ypclAnNWHTt25KKLLsLv9zf03GrOsBvJQmenh84uGuXLEoFpkZ0eat7OnTs5ePAgXbt2ZciQIV6HE5csETTPvZEMLBEYj9mO2rzgozW7kax97MbF5m3atInjx4/Tr18/+vZtsvedsLBEYFoUfGOZ3+/3OJrYYqeFQpeXl0dOTg579+5l9+7dLc+QRKLVkaElAtOi8847j969e3PkyBG2bNnidTgxxT2KtUTQfj6fr+HqIasVfFW0ypclAtMiEbF2giacPHmSdevWfeXRi6Z9rHw1zWoEJqa4BdGO2M747LPPqK+vZ9SoUWRlZXkdTlyzdoKv++KLLygrKyMjI4NRo0ZFdF2WCEyrWCL4Ort/IHzGjRuHz+djzZo1VFVVeR1OTHCvFho7dixpaWkRXZclAtMqBQUFpKWlsWnTJo4ePep1ODHBEkH4ZGdnc+GFF1JXV2cPQnJE84l3lghMq6Snp1NYWAicOVJJZsE3klkiCA+rdX5VNA80LBGYVrMGvTO2bNnCkSNH6NOnT0SeGJWMLBGcUVdX13AnfzSuSLNEYFrN3VEj8czUeBN8WZ/dSBYewYkg2R+EtG7dOk6dOsXgwYNp76N528ISgWk198hkxYoV1NXVeRyNt9xkOHHiRI8jSRwDBw6kZ8+e9iAkzhxoRKt8WSIwrdarVy8GDRrEyZMnKS0t9TocT1lDcfgF36+S7LVO9/tHq3xZIjBt4h6hJPOOeuTIETZv3kxGRgYFBQVeh5NQ3PKV7O0EcVUjEJEuIjJfRLY6f3ObmOYiEflURDaKyHoR+VbQuD+LyA4RWeu8LgolHhN51qB3prE8Gtd3J/BSHqgAABukSURBVBurEcDu3bvZvXs3OTk5DB8+PCrrDLVG8BCwQFWHAAucz42dAr6jqiOAqcAzIpITNP7fVfUi57U2xHhMhNkRm50WiiT3KVzJfL9K8IUIPl90TtqEupYZwEvO+5eAGxpPoKpbVHWr834fcAiIfDO4iYgLLriA7Oxsdu3axd69e70OxxPRPn+bTILvV0nWgw0vyleoiaCnqu533h8Aep5tYhEZB6QB24IGP+GcMnpaRJp9IKeI3CciJSJSUlFREWLYpr1SUlKSuvpeW1sb1eu7k1Gy1zpjMhGIyEcisqGJ14zg6TRw4W+zF/+KSG/gL8D3VNXt1P5hIA8YC3QB/qO5+VX1eVUtVNXCaFxXa5qXzInA7Qtn2LBhUbm+Oxkl8wUJlZWVrF27lpSUFMaPHx+19aa2NIGqXtncOBE5KCK9VXW/80N/qJnpsoH3gJ+pakP/BEG1iRoR+RPwb22K3njC3VE//vhjjyOJPvc72/0DkeMeaKxYsYLTp08nVYP8ihUr8Pv9FBYWRrVH21BPDc0F7nTe3wm83XgCEUkD3gReVtU5jcb1dv4KgfaFDSHGY6Jg/PjxpKSksG7dOk6cOOF1OFHlHqVecsklHkeSuLp160ZeXh7V1dWsWbPG63CiyqvyFWoieBK4SkS2Alc6nxGRQhH5ozPNLcClwHebuEz0ryJSCpQC3YDHQ4zHREFWVhYFBQXU19ezYsUKr8OJGlW1GkGUJOvpIa/uWA8pEajqYVWdrKpDVPVKVT3iDC9R1Xuc96+oaoegS0QbLhNV1StU9UJVHamqt6tqch1exjH3iCWZTg9t27aNQ4cO0b17d4YMGeJ1OAktGctXXV1dwz0qcZUITPJKxnaC4NqAdTQXWcE1gmTpgK60tJQTJ04waNAgevfuHdV1WyIw7eIesRUXFydNB3RuIrD2gcgbPHgwPXr04NChQ0nTAZ2X5csSgWmXXr16MXjwYE6cOMG6deu8Dicqli1bBsA3vvENjyNJfCLSsJ3d7Z7ovCxflghMuyXTjnrw4EG2bNlCZmamdTQXJe6RcTKUL1Vt+J5WIzBxJZl2VPdqjqKiIjp06OBxNMnBPdBIhnao7du3c+DAAbp168awYcOivn5LBKbdgmsEid6gZ6eFom/UqFF06tSJ8vJyDhw44HU4ERVcG/DiQgRLBKbdBg8eTM+ePamoqGDLli1ehxNRlgiiLzU1taE/p0SvdXpdviwRmHYTES699FIgsXfUyspK1qxZQ2pqalT7fzEkRfkCSwQmzrkFd+nSpR5HEjmffPIJfr+fMWPG0KlTJ6/DSSrJUL7279/P1q1b6dSpk2cXIlgiMCFxj9gSeUd1v5v7XU30jB8/nrS0NNavX5+wD6pxawMTJ04kNbXFfkAjwhKBCcnIkSPJyclh165d7Nq1y+twIsISgXcyMjIYP378V/p5SjRLliwBvC1flghMSFJSUhL6MtKqqipWrlyJiFhHcx5J9FpnLBxoWCIwIUvkHXXlypWcPn2a/Px8cnNzvQ4nKSVy+Tp8+DAbNmwgPT2dsWPHehaHJQITskmTJgGwePFibwOJgFiotie7CRMmkJKSwqpVq6isrPQ6nLBya9Hjx48nPb3ZJ/VGnCUCE7LRo0fTqVMntm7dyr59+7wOJ6zc5Hb55Zd7G0gS69y5M4WFhdTX1yfc8wlipXxZIjAhS01NbbjMzz2CTgTV1dUN/cNbjcBbl112GZB4tU73+7jfzyuWCExYJOKOunLlSqqrq8nPz6dr165eh5PUErF8HTlyhPXr15Oenk5RUZGnsVgiMGGRiDtqrBytmcA19ikpKZSUlCRMO8HSpUtRVYqKisjIyPA0lpASgYh0EZH5IrLV+dvkZRUiUh/0vOK5QcMHisgKESkXkdnOg+5NHHLbCbZs2ZIw7QSWCGJHIrYTxFL5CrVG8BCwQFWHAAucz02pCnpe8fSg4b8AnlbVwcBR4O4Q4zEeSU1NbTiPvmjRIo+jCV11dTWffPLJVx6QYrzlNqguXLjQ40jCw/0eiZAIZgAvOe9fAm5o7YwS6Gv1CmBOe+Y3scct0ImQCD799FNqamrIz8+nW7duXodjOJMIEqF8VVRUUFpaSkZGhuftAxB6Iuipqvud9weAns1MlyEiJSJSLCLuj31X4EtVdR94uwfo29yKROQ+ZxklFRUVIYZtIuGKK64AYMGCBR5HEjr3O7jfyXhv4sSJdOjQgdWrV8d9v0PuaaGJEyd63j4ArUgEIvKRiGxo4jUjeDoNPJmkuaeT9FfVQuA24BkROb+tgarq86paqKqF3bt3b+vsJgouuugicnNz2blzJzt27PA6nJC41fbJkyd7HIlxZWVlMWHCBPx+f9zfZRxr5avFRKCqV6rqyCZebwMHRaQ3gPP3UDPL2Ov83Q4sBgqAw0COiLjd7Z0L7A35GxnPpKSkNJweiufzuJWVlaxcuZKUlBRrH4gxiVLrjLUaZ6inhuYCdzrv7wTebjyBiOSKSLrzvhswEdjk1CAWATPPNr+JL27BjudEsGzZMurr6xk7dizZ2dleh2OCJEL52r17N1u3bqVz586MGTPG63CA0BPBk8BVIrIVuNL5jIgUisgfnWmGAyUiso7AD/+TqrrJGfcfwIMiUk6gzeCFEOMxHgs+YovX5xh/9NFHQOwcrZkzxo8fT2ZmJhs3bmT//v0tzxCD3NrApEmTPHv+QGMhJQJVPayqk1V1iHMK6YgzvERV73Hef6KqF6rqKOfvC0Hzb1fVcao6WFVvVtWa0L6O8drw4cPp06cPBw8epLS01Otw2mX+/PkAXHXVVR5HYhpLS0tr6OTQTdjxJhbLl91ZbMJKRBoKuFvg48n+/fvZsGEDmZmZDQ9ON7ElnsuX3+9vSGCWCExCi+cd1d1JJ02a5Gm3wKZ5bvn66KOP4u70Y2lpKYcOHaJv377k5eV5HU4DSwQm7K688kog0JdKdXW1x9G0TSxW281XjRgxgt69e7N//342btzodThtEly+AvfUxgZLBCbsevbsSX5+PlVVVXHVL4yq8uGHHwKWCGKZiDQcbLj/r3gRq+XLEoGJiKuvvhqAefPmeRxJ661fv56DBw/St29fRowY4XU45iymTJkCxFciqKqqYtmyZV9JZLHCEoGJCDcR/OMf//A4ktZzY50yZUpMVdvN17lH1EuWLKGqqsrjaFpnyZIlVFdXU1BQQI8ePbwO5yssEZiIuOSSS8jKyqK0tJS9e+PjhnE3EUydOtXjSExLevbsyejRo6muro6bp+LFcvmyRGAiIj09veGGrHg4PVRZWcnHH3+Mz+eLufO3pmnXXHMNED+1TjdON+5YYonARIx75BMPO+rChQupq6ujqKiI3Nwmn69kYkw8la8dO3ZQVlbGOeecExPdTjdmicBEjLujfvjhh9TW1noczdm9//77wJm2DRP7ioqKOOeccygrK2Pbtm1eh3NWH3zwARDobTRWupUIZonARMygQYPIy8vj2LFjMX0Zqary7rvvAnDdddd5HI1prdTU1IbE/d5773kczdnFevmyRGAiyi347o4Qi9auXcu+ffvo06cPBQUFXodj2iAeytfJkycbekudNm2ax9E0zRKBiah42FHd2K699lq7bDTOXHPNNYgIixcvprKy0utwmrRgwQJqamoYN24cPXs29xBHb1kiMBF18cUXk5OTQ1lZGVu3bvU6nCbFerXdNK9bt25MmDCB2tramO3bKh7KlyUCE1EdOnRoqA6//XbsPXdo//79rFy5ko4dO8bMYwNN20yfPh2At956y+NIvs7v9zN37lxExBKBSW7XX389AHPnzvU4kq9zY5o8eTJZWVkeR2Paw00E7733HnV1dR5H81UrVqzg4MGDnHfeeVx00UVeh9MsSwQm4qZNm0ZaWhoff/wxhw41+Vhrz7z55psA3HjjjR5HYtpr+PDhDBs2jCNHjrBs2TKvw/kKt5Zyww03xHT7kyUCE3HZ2dlMnjwZVeWdd97xOpwGx44dY+HChfh8voZai4lPN9xwA3AmsccCVW2Ix40vVoWUCESki4jMF5Gtzt+v3ZIpIpeLyNqgV7WI3OCM+7OI7AgaF7t1JxMS94g7lnbU999/n9raWr7xjW/QvXt3r8MxIXDL11tvvRUzD6vZtGkTW7dupWvXrlxyySVeh3NWodYIHgIWqOoQYIHz+StUdZGqXqSqFwFXAKeA4L5j/90dr6prQ4zHxKjp06cjIsyfP59jx455HQ4Af//734HYP1ozLRs7dix9+vRh9+7dfPbZZ16HA5wpX9OnT4/Ju4mDhZoIZgAvOe9fAlrao2YCH6jqqRDXa+JMz549mTRpEqdPn46JRuMTJ0403I06c+ZMj6MxofL5fA3/x9dee83jaALcOG6++WaPI2lZqImgp6rud94fAFq6W+JW4G+Nhj0hIutF5GkRsYfEJrBbbrkFiI0d9d1336W6upqJEydy7rnneh2OCYPg8uX16aGNGzeyceNGcnNz4+Ky5BYTgYh8JCIbmnjNCJ5OA1u+2a0vIr2BC4HgPokfBvKAsUAX4D/OMv99IlIiIiUVFRUthW1i0E033YTP52PevHkcPXrU01hmz54NwLe+9S1P4zDhM2HCBM4991x2795NcXGxp7G4Bzvf/OY3SUtL8zSW1mgxEajqlao6sonX28BB5wfe/aE/27WBtwBvqmpDN5Squl8DaoA/AePOEsfzqlqoqoXWsBefevToweWXX05tbS1vvPGGZ3EcO3aMd999F5/Px0033eRZHCa8fD5fw2mY//3f//UsDlVtWL9bS4l1oZ4amgvc6by/EzjbraOzaHRaKCiJCIH2hQ0hxmNi3Le//W3A2x3173//O3V1dVx22WX06dPHszhM+N12220AvP76657dXLZq1SrKy8vp2bNnw8OZYl2oieBJ4CoR2Qpc6XxGRApF5I/uRCIyAOgHNH6m3F9FpBQoBboBj4cYj4lxN954I+np6SxatIg9e/Z4EsPLL78MnElKJnGMGTOGoUOHcvDgQc/6HvrLX/4CBE47xvrVQq6QEoGqHlbVyao6xDmFdMQZXqKq9wRNt1NV+6qqv9H8V6jqhc6ppttV9UQo8ZjYl5OTw/Tp01FV/vrXv0Z9/bt27WLJkiVkZGTY1UIJSES44447gDM/yNFUW1vL3/4WOPHhxhEP7M5iE3XuDvLyyy9H/eoON/nMmDGD7OzsqK7bRMftt98OBG4uO378eFTX/eGHH1JRUUFeXh5jxoyJ6rpDYYnARN3UqVPp0aMHmzZtYsWKFVFbr9/v54UXXgDgO9/5TtTWa6JrwIABTJo0iaqqqoaj82j54x8DZ8S/853vxHTfQo1ZIjBR16FDB7773e8C8Pzzz0dtvYsWLWL79u3069fPnk2c4O69914A/vCHP0Rtnfv37+edd94hNTWV733ve1FbbzhYIjCeuOeeQBPS7Nmzo9blhJt07r77blJSUqKyTuONm266idzcXFatWsXq1aujss4//elP1NfXc/3119OrV6+orDNcLBEYTwwZMoTLL7+cU6dO8dJLL7U8Q4gOHjzI66+/js/n46677or4+oy3MjIyGk7/PffccxFfX319Pb/97W8BuO+++yK+vnCzRGA888ADDwDwm9/8Br/f38LUofn973+PqnLDDTfQr1+/iK7LxIYf/OAH+Hw+/va3v3HkyJGIruudd95hz549DB06lKuuuiqi64oESwTGM9OnT6dfv36UlZUxb968lmdop5qaGn7zm98AZ5KPSXxDhw5lypQpVFVVRbyt4Fe/+hUQSD7xeNrREoHxTGpqasMP89NPPx2x9cyePZuDBw+Sn5/PZZddFrH1mNjzox/9CIBnn32W2traFqZun7Vr17J48WI6d+4cd43ELksExlP33HMPmZmZzJ8/PyKNen6/n1/84hdA4Echni7pM6GbMmUKw4cPZ+/evRHr1sQtX3fffXfc3ptiicB4qkuXLtx///0A/PznPw/78t944w02bdpEv379rEuJJCQiPPzwwwA88cQT1NfXh3X5//znP5k9ezZpaWn85Cc/Ceuyo8kSgfHcv/3bv5GRkcFbb73F+vXrw7Zcv9/fkFweeuihuOgO2ITfrFmzOP/889m6dSuvvvpqWJf9xBNPoKp873vfi+vnWlgiMJ7r1asXP/jBDwD46U9/Grblvvbaa6xfv54BAwbYJaNJLDU1lUceeQSfz8fDDz/M6dOnw7Lc0tJSXnnlFbKysnjooa89pTeuWCIwMeGnP/0pnTt3Zt68eSxYsCDk5VVXVzfsnI888ggZGRkhL9PEr9tvv51hw4axe/fusNxXoKoNp4LuuusuBgwYEPIyvWSJwMSEnj17NpzLffDBB0PuS/7pp59m165d5OfnN3RnYZJXamoqTz31FACPPfYYoT7l8P3332f+/Pmcc845PProo+EI0VOWCEzM+NGPfkT//v1Zv349v/zlL9u9nC1btvDYY48B8Mtf/jIur+s24Tdt2jSuuuoqvvzyy4bLStujsrKy4VTmo48+Srdu3cIVondUNe5eY8aMUZOY5s2bp4Cmp6draWlpm+evq6vTCRMmKKB33nln+AM0ca28vFwzMzMV0L///e/tWsa9996rgI4ZM0Zra2vDHGFkASXaxG+q5z/q7XlZIkhs9957r/p8Ps3Ly9PKyso2zfvII48ooAMGDNDDhw9HKEITz5555hlNT0/XTp066fbt29s07+zZsxXQ7OxsXb9+fYQijBxLBCZuVFZWal5engI6c+ZMraura9V8c+bMUUB9Pp8uWLAgwlGaeOX3+/W6665TQC+66CI9fvx4q+Zbu3ZtQ23i2WefjXCUkWGJwMSVzZs3a+fOnRXQO+64o8Vk8Oabb2qHDh0U0KeeeipKUZp4dfToUR08eLACevHFF+uxY8fOOv2aNWu0a9euDeXR7/dHKdLwikgiAG4GNgJ+oPAs000FyoBy4KGg4QOBFc7w2UBaa9ZriSA5LF68WLOyshTQoqIi3bt379emqaur08cff1wBBfTHP/5x3O6kJrrKy8u1X79+CmivXr107dq1TU7317/+taF8TZs2TauqqqIcafhEKhEMB4YBi5tLBEAKsA0YBKQB64ALnHGvAbc6738H3N+a9VoiSB7Lly/XIUOGKKCdO3fWW2+9Vd9880396KOP9KGHHmo4hZSSkqK//OUvLQmYNtmxY4dOnDhRAU1NTdWZM2fq7373O124cKE+8cQTDeMAvf/++7W6utrrkEMS0VNDLSSCCcC8oM8POy8BvgBSm5rubC9LBMll3759eu211zbskI1f/fr103nz5nkdpolTJ06c0O9///vq8/maLF/Z2dn6/PPPJ8RBRnOJQALjQiMii4F/U9WSJsbNBKaq6j3O5zuA8cD/AYpVdbAzvB/wgaqObGYd9wHuo3+GETjV1B7dCCSgWGNxtY3F1TYWV9skalz9VbV744GpLc0lIh8BTT2A82eq+nYIAbWJqj4PhPykcxEpUdXCMIQUVhZX21hcbWNxtU2yxdViIlDVK0Ncx14g+NmA5zrDDgM5IpKqqnVBw40xxkRRNLqY+AwYIiIDRSQNuBWY65yvWgTMdKa7E4haDcMYY0xASIlARG4UkT0EGnrfE5F5zvA+IvI+gHO0/wAwD9gMvKaqG51F/AfwoIiUA12BF0KJp5VCPr0UIRZX21hcbWNxtU1SxRWWxmJjjDHxy3ofNcaYJGeJwBhjklxCJQIRmSoiZSJSLiJfe3aciKSLyGxn/AoRGRA07mFneJmIXB3luB4UkU0isl5EFohI/6Bx9SKy1nnNjXJc3xWRiqD13xM07k4R2eq87oxyXE8HxbRFRL4MGheR7SUiL4rIIRHZ0Mx4EZFfOzGvF5HRQeMiua1aiuvbTjylIvKJiIwKGrfTGb5WRL52D1CE47pMRI4F/a8eDRp31v9/hOP696CYNjjlqYszLpLbq5+ILHJ+BzaKyP/bxDSRK2NN3WUWjy/O0pVF0DQ/AH7nvL8VmO28v8CZPp1A/0fbgJQoxnU5kOm8v9+Ny/l8wsPt9V3g/2ti3i7AdudvrvM+N1pxNZr+X4EXo7C9LgVGAxuaGT8N+IDAHfNFwIpIb6tWxnWxuz7gGjcu5/NOoJtH2+sy4N1Q///hjqvRtNcDC6O0vXoDo533nYEtTeyPEStjiVQjGAeUq+p2VT0NvArMaDTNDOAl5/0cYLKIiDP8VVWtUdUdBDrBGxetuFR1kaqecj4WE7inItJas72aczUwX1WPqOpRYD6BjgW9iGsW8LcwrbtZqroUOHKWSWYAL2tAMYF7ZHoT2W3VYlyq+omzXohe2WrN9mpOKOUy3HFFpWwBqOp+VV3tvK8kcIVl30aTRayMJVIi6AvsDvq8h69vyIZpNHBZ6zECl622Zt5IxhXsbgJZ35UhIiUiUiwiN4QpprbEdZNTDZ0jgW5A2jJvJOPCOYU2EFgYNDhS26slzcUdyW3VVo3LlgIfisgqCXThEm0TRGSdiHwgIiOcYTGxvUQkk8CP6d+DBkdle0nglHUBgZ6Zg0WsjLV4Z7GJHhG5HSgEJgUN7q+qe0VkELBQREpVdVuUQnoH+Juq1ojI/0OgNnVFlNbdGrcCc1S1PmiYl9srZonI5QQSwSVBgy9xtlUPYL6I/NM5Yo6G1QT+VydEZBrwFjAkSutujeuB5aoaXHuI+PYSkU4Eks+PVPV4OJd9NolUI2iuK4smpxGRVOAcAl1dtGbeSMaFiFwJ/AyYrqo17nBV3ev83U6gl9eCaMWlqoeDYvkjMKa180YyriC30qjqHsHt1ZLm4o7ktmoVEckn8P+boaqH3eFB2+oQ8CbhOx3aIlU9rqonnPfvAx1EpBsxsL0cZytbEdleItKBQBL4q6q+0cQkkStjkWj48OJFoHazncCpAreRaUSjaf6FrzYWv+a8H8FXG4u3E77G4tbEVUCggWxIo+G5QLrzvhuwlTA1nLUyrt5B728k0FssBBqldjjx5Trvu0QrLme6PAKNdxKN7eUscwDNN35ey1cb8lZGelu1Mq7zCLR5XdxoeBbQOej9JwR6CY5WXL3c/x2BH9TPnW3Xqv9/pOJyxp9DoB0hK1rby/nuLwPPnGWaiJWxsG3cWHgRaFXfQuBH9WfOsMcIHGUDZACvOzvGSmBQ0Lw/c+YrA66JclwfAQeBtc5rrjP8YqDU2RlKgbujHNd/E3gC3ToC/ULlBc17l7Mdy4HvRTMu5/P/AZ5sNF/EtheBo8P9QC2Bc7B3A98Hvu+MF+A5J+ZSgp7PEeFt1VJcfwSOBpWtEmf4IGc7rXP+xz+LclwPBJWtYoISVVP//2jF5UzzXQIXjwTPF+ntdQmBNoj1Qf+radEqY9bFhDHGJLlEaiMwxhjTDpYIjDEmyVkiMMaYJGeJwBhjkpwlAmOMSXKWCIwxJslZIjDGmCT3/wMkNO3Em+ZKOgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# NBVAL_IGNORE_OUTPUT\n", "import matplotlib.pyplot as plt\n", "fig, axes = plt.subplots()\n", "axes.set_ylim((-1., 1.))\n", "plot(results[0], axes=axes)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we'll create a function that tells matplotlib how to draw each frame of the animation, which in our case will just be plotting the value at that timestep.\n", "The `FuncAnimation` function will call this on the list of results that we pass in, together with a given interval in milliseconds between each frame.\n", "Finally, we'll use the IPython API to render the animation in the notebook." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "
\n", " \n", "
\n", " \n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "
\n", "
\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# NBVAL_IGNORE_OUTPUT\n", "from matplotlib.animation import FuncAnimation\n", "\n", "def animate(u):\n", " axes.clear()\n", " firedrake.plot(u, axes=axes)\n", " axes.set_ylim((-1., 1.))\n", "\n", "interval = 4e3 * float(dt)\n", "animation = FuncAnimation(fig, animate, frames=results, interval=interval)\n", "\n", "from IPython.display import HTML\n", "HTML(animation.to_jshtml())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## A faster implementation\n", "\n", "Although the code we wrote above works fine, it can be quite slow. In particular, each call to `solve` necessitates rederiving the symbolic Jacobian, building new matrices and vectors and solver objects, using them once, and then destroying them. To avoid this, we can create a solver object and reuse it.\n", "\n", "This is what the `solve` call does internally, only it then immediately discards all of this work.\n", "\n", "We start by creating a `NonlinearVariationalProblem` which gathers the information about the problem. The residual, the solution variable, any boundary conditions, and so forth." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "problem = NonlinearVariationalProblem(F, u_n1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we create a `NonlinearVariationalSolver`. Here we provide the problem to be solved, and any options to the solver.\n", "\n", "In this case, we will modify the solver options used, noting that in one dimension, an LU factorisation produces no fill and is, obviously, an exact solve." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "solver = NonlinearVariationalSolver(problem, solver_parameters={\"ksp_type\": \"preonly\",\n", " \"pc_type\": \"lu\"})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we just write the time loop as before, but instead of writing `solve(F == 0, u_n1)`, we just call the `solve` method on our `solver` object." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "t = 0\n", "t_end = 0.5\n", "while t <= t_end:\n", " solver.solve()\n", " u_n.assign(u_n1)\n", " t += dt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercise 1\n", "\n", "Compare the speed of the two implementation choices on a mesh with 1000 elements.\n", "\n", "- Hint: You can use the \"notebook magic\" `%%timeit` to time the execution of a notebook cell." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercise 2\n", "\n", "Implement Crank-Nicolson timestepping instead of backward Euler.\n", "\n", "- Hint 1: The Crank-Nicolson scheme writes:\n", "\n", " $$\\frac{\\partial u}{\\partial t} + G(u) = 0$$\n", "\n", " as\n", "\n", " $$ \\frac{u^{n+1} - u^n}{\\Delta t} + \\frac{1}{2}\\left[G(u^{n+1}) + G(u^n)\\right] = 0$$\n", "\n", "\n", "- Hint 2: It might be convenient to write a python function that returns $G(u)$ given a $u$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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.7.4" }, "widgets": { "application/vnd.jupyter.widget-state+json": { "state": { "483f8c5b55a94cb5b9b3c6223be10e49": { "model_module": "@jupyter-widgets/controls", "model_module_version": "1.5.0", "model_name": "VBoxModel", "state": { "_dom_classes": [ "widget-interact" ], "_model_module": "@jupyter-widgets/controls", "_model_module_version": "1.5.0", "_model_name": "VBoxModel", "_view_count": null, "_view_module": "@jupyter-widgets/controls", "_view_module_version": "1.5.0", "_view_name": "VBoxView", "box_style": "", "children": [ "IPY_MODEL_8f43ada5d62046f7a6a50fbcc6587acf", "IPY_MODEL_ff45eb9314f84c0b99d6cbc4dafaeea7" ], "layout": "IPY_MODEL_c987f4ac87ac4cae9b51332042a41dd2" } }, "8f43ada5d62046f7a6a50fbcc6587acf": { "model_module": "@jupyter-widgets/controls", "model_module_version": "1.5.0", "model_name": "IntSliderModel", "state": { "_dom_classes": [], "_model_module": "@jupyter-widgets/controls", "_model_module_version": "1.5.0", "_model_name": "IntSliderModel", "_view_count": null, "_view_module": "@jupyter-widgets/controls", "_view_module_version": "1.5.0", "_view_name": "IntSliderView", "continuous_update": true, "description": "index", "description_tooltip": null, "disabled": false, "layout": "IPY_MODEL_e34369af8020497289732bb687d2e8e2", "max": 50, "min": 0, "orientation": "horizontal", "readout": true, "readout_format": "d", "step": 1, "style": "IPY_MODEL_eb5e810f3cf7471c8f32d7e6ef99b66b", "value": 0 } }, "c987f4ac87ac4cae9b51332042a41dd2": { "model_module": "@jupyter-widgets/base", "model_module_version": "1.2.0", "model_name": "LayoutModel", "state": { "_model_module": "@jupyter-widgets/base", "_model_module_version": "1.2.0", "_model_name": "LayoutModel", "_view_count": null, "_view_module": "@jupyter-widgets/base", "_view_module_version": "1.2.0", "_view_name": "LayoutView", "align_content": null, "align_items": null, "align_self": null, "border": null, "bottom": null, "display": null, "flex": null, "flex_flow": null, "grid_area": null, "grid_auto_columns": null, "grid_auto_flow": null, "grid_auto_rows": null, "grid_column": null, "grid_gap": null, "grid_row": null, "grid_template_areas": null, "grid_template_columns": null, "grid_template_rows": null, "height": null, "justify_content": null, "justify_items": null, "left": null, "margin": null, "max_height": null, "max_width": null, "min_height": null, "min_width": null, "object_fit": null, "object_position": null, "order": null, "overflow": null, "overflow_x": null, "overflow_y": null, "padding": null, "right": null, "top": null, "visibility": null, "width": null } }, "e34369af8020497289732bb687d2e8e2": { "model_module": "@jupyter-widgets/base", "model_module_version": "1.2.0", "model_name": "LayoutModel", "state": { "_model_module": "@jupyter-widgets/base", "_model_module_version": "1.2.0", "_model_name": "LayoutModel", "_view_count": null, "_view_module": "@jupyter-widgets/base", "_view_module_version": "1.2.0", "_view_name": "LayoutView", "align_content": null, "align_items": null, "align_self": null, "border": null, "bottom": null, "display": null, "flex": null, "flex_flow": null, "grid_area": null, "grid_auto_columns": null, "grid_auto_flow": null, "grid_auto_rows": null, "grid_column": null, "grid_gap": null, "grid_row": null, "grid_template_areas": null, "grid_template_columns": null, "grid_template_rows": null, "height": null, "justify_content": null, "justify_items": null, "left": null, "margin": null, "max_height": null, "max_width": null, "min_height": null, "min_width": null, "object_fit": null, "object_position": null, "order": null, "overflow": null, "overflow_x": null, "overflow_y": null, "padding": null, "right": null, "top": null, "visibility": null, "width": null } }, "eaddcef7c5f24f0bae9a44e7536afc39": { "model_module": "@jupyter-widgets/base", "model_module_version": "1.2.0", "model_name": "LayoutModel", "state": { "_model_module": "@jupyter-widgets/base", "_model_module_version": "1.2.0", "_model_name": "LayoutModel", "_view_count": null, "_view_module": "@jupyter-widgets/base", "_view_module_version": "1.2.0", "_view_name": "LayoutView", "align_content": null, "align_items": null, "align_self": null, "border": null, "bottom": null, "display": null, "flex": null, "flex_flow": null, "grid_area": null, "grid_auto_columns": null, "grid_auto_flow": null, "grid_auto_rows": null, "grid_column": null, "grid_gap": null, "grid_row": null, "grid_template_areas": null, "grid_template_columns": null, "grid_template_rows": null, "height": null, "justify_content": null, "justify_items": null, "left": null, "margin": null, "max_height": null, "max_width": null, "min_height": null, "min_width": null, "object_fit": null, "object_position": null, "order": null, "overflow": null, "overflow_x": null, "overflow_y": null, "padding": null, "right": null, "top": null, "visibility": null, "width": null } }, "eb5e810f3cf7471c8f32d7e6ef99b66b": { "model_module": "@jupyter-widgets/controls", "model_module_version": "1.5.0", "model_name": "SliderStyleModel", "state": { "_model_module": "@jupyter-widgets/controls", "_model_module_version": "1.5.0", "_model_name": "SliderStyleModel", "_view_count": null, "_view_module": "@jupyter-widgets/base", "_view_module_version": "1.2.0", "_view_name": "StyleView", "description_width": "", "handle_color": null } }, "ff45eb9314f84c0b99d6cbc4dafaeea7": { "model_module": "@jupyter-widgets/output", "model_module_version": "1.0.0", "model_name": "OutputModel", "state": { "_dom_classes": [], "_model_module": "@jupyter-widgets/output", "_model_module_version": "1.0.0", "_model_name": "OutputModel", "_view_count": null, "_view_module": "@jupyter-widgets/output", "_view_module_version": "1.0.0", "_view_name": "OutputView", "layout": "IPY_MODEL_eaddcef7c5f24f0bae9a44e7536afc39", "msg_id": "", "outputs": [] } } }, "version_major": 2, "version_minor": 0 } } }, "nbformat": 4, "nbformat_minor": 1 }