{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "In the previous post, we explored how difficult it is to solve the simplest hyperbolic conservation law, the scalar advection equation.\n", "To solve this PDE accurately, we had to understand how the partly arbitrary choice of a numerical flux can make or break the stability of the spatial discretization, how low order schemes are very diffusive, and how higher-order explicit schemes introduce spurious maxima and minima that we can only control through a nonlinear flux limiting procedure.\n", "The scalar advection equation is comparatively simple in that signals can only propagate along the predetermined velocity field.\n", "In this post, we'll look at something more realistic and much more difficult: the shallow water equations.\n", "The shallow water equations are a system of equations rather than a scalar problem, and as such they can exhibit non-trivial wave propagation in a way that the advection equation can't.\n", "They're a great model for testing numerical solvers because they're both simple enough to keep in your head all at once, and yet at the same time they exhibit many of the complexities of more \"serious\" models -- nonlinearity, non-trivial conserved quantities, mimetic properties, the list goes on.\n", "\n", "The shallow water equations can be derived from the incompressible Euler equations of fluid dynamics with a free surface under the assumption that the horizontal length scale is much longer than the vertical one.\n", "This approximation reduces the unknowns to the thickness $h$ of the fluid and the depth-averaged velocity $u$.\n", "The conservation laws are for mass and momentum:\n", "\n", "\\begin{align}\n", "\\frac{\\partial}{\\partial t}h + \\nabla\\cdot hu & = 0 \\\\\n", "\\frac{\\partial}{\\partial t}hu + \\nabla\\cdot\\left(hu\\otimes u + \\frac{1}{2}gh^2I\\right) & = -gh\\nabla b\n", "\\end{align}\n", "\n", "where $g$ is the acceleration due to gravity, $b$ is the bathymetry, and $I$ is the identity matrix.\n", "This problem is a little more complicated because of the time derivative on $h\\cdot u$, a combination of two of the state variables.\n", "To make things a little easier, we'll instead work with the momentum $q = h\\cdot u$ and rewrite the system as\n", "\n", "\\begin{align}\n", "\\frac{\\partial}{\\partial t}h + \\nabla\\cdot q & = 0 \\\\\n", "\\frac{\\partial}{\\partial t}q + \\nabla\\cdot\\left(h^{-1}q\\otimes q + \\frac{1}{2}gh^2I\\right) & = -gh\\nabla b.\n", "\\end{align}\n", "\n", "As in the previous post, we'll use a discontinuous Galerkin basis.\n", "We showed that there is more than one way to come up with a discretized problem that is consistent with the idealized one and this is manifested in which numerical flux to use.\n", "Things get much more interesting for systems of PDE, which can have more than one characteristic speed besides that of the background flow field.\n", "In the following, I'll assume you're familiar with the fact that the characteristic wave speed for the shallow water equations is\n", "\n", "$$c = |u| + \\sqrt{gh}.$$\n", "\n", "The fact that the wave speed now depends on the solution and that waves propagate in all directions instead of only along a pre-set vector field has several consequences.\n", "First, we can't pick a CFL-stable timestep from the outset because the fluid velocity and thickness could increase well beyond their initial values.\n", "The only options for timestepping are to use an adaptive procedure or a whole mess of trial and error.\n", "Second, we have to think harder about numerical fluxes.\n", "For scalar equations, we can use the numerical flux to mimic the properties of upwind finite difference schemes, but for systems this reasoning doesn't work -- there can be waves simultaneously propagating in both normal directions at a given internal facet." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Spatial discretization\n", "\n", "We'll examine several different test problems for the shallow water equations, so to avoid a huge amount of repeated code we'll first write a few Python functions to set up the weak form of the problem." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import firedrake\n", "from firedrake import Constant\n", "g = Constant(9.81)\n", "I = firedrake.Identity(2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The component of the fluxes in the interior of the cells is exactly what you get from applying the divergence theorem to the original system." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from firedrake import inner, grad, dx\n", "def cell_flux(z):\n", " Z = z.function_space()\n", " h, q = firedrake.split(z)\n", " ϕ, v = firedrake.TestFunctions(Z)\n", " \n", " f_h = -inner(q, grad(ϕ)) * dx\n", "\n", " F = outer(q, q) / h + 0.5 * g * h**2 * I\n", " f_q = -inner(F, grad(v)) * dx\n", "\n", " return f_h + f_q" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we need to add the facet terms and this is where the numerical fluxes come in.\n", "We'll start with a function to compute the central flux.\n", "This is the easy part -- there are no choices to make." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "from firedrake import avg, outer, dS\n", "def central_facet_flux(z):\n", " Z = z.function_space()\n", " h, q = firedrake.split(z)\n", " ϕ, v = firedrake.TestFunctions(Z)\n", " \n", " mesh = z.ufl_domain()\n", " n = firedrake.FacetNormal(mesh)\n", "\n", " f_h = inner(avg(q), ϕ('+') * n('+') + ϕ('-') * n('-')) * dS\n", "\n", " F = outer(q, q) / h + 0.5 * g * h**2 * I\n", " f_q = inner(avg(F), outer(v('+'), n('+')) + outer(v('-'), n('-'))) * dS\n", " \n", " return f_h + f_q" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The central flux by itself is unstable with explicit timestepping schemes and to reocver stability we need an extra diffusive flux term.\n", "The subtle choices are all in this diffusive flux.\n", "**For the remainder of this demo, we'll use the Lax-Friedrichs numerical flux.**\n", "This flux uses the maximum outgoing wave speed to set the local diffusion coefficient.\n", "An upper bound for the outgoing wavespeed is then $|c| = |q/h \\cdot n| + \\sqrt{gh}$.\n", "I say an *upper bound* and not the actual maximum value because the system could exhibit the shallow water equivalent of supersonic flow -- where the speed $|u|$ exceeds $\\sqrt{gh}$, both waves could be moving in the same direction rather than opposite each other.\n", "\n", "The vast majority of the literature you'll find on DG methods uses the Lax-Friedrichs flux.\n", "For example, [Cockburn and Shu (1998)](https://doi.org/10.1006/jcph.1998.5892) in the last of their famous series of papers on the Runge-Kutta DG method consider only the Lax-Friedrichs flux and neglect to even mention the alternatives.\n", "This can be confusing for beginners to the subject because it isn't clear at what points in the process you (the programmer) have choices to make and where you don't.\n", "Part of the reason for this singular focus is that the Lax-Friedrichs flux is the simplest to implement, but several studies have found -- at least for problems without shock waves -- that other flux functions offer negligible gains in accuracy; see for example [Qiu (2006)](https://doi.org/10.1007/s10915-006-9109-5) and [Bernard et al. (2009)](https://doi.org/10.1016/j.jcp.2009.05.046).\n", "In that case, there isn't much value in using a different flux function that may be more expensive to compute and make the code harder to understand and maintain.\n", "The choice of numerical flux is more consequential for problems with shock waves or for more complex systems (see for example [Beck et al. (2014)](https://doi.org/10.1007/978-3-319-01601-6_11), which studied turbulent flow)." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "from firedrake import sqrt, max_value\n", "def lax_friedrichs_facet_flux(z):\n", " Z = z.function_space()\n", " h, q = firedrake.split(z)\n", " ϕ, v = firedrake.TestFunctions(Z)\n", " \n", " mesh = h.ufl_domain()\n", " n = firedrake.FacetNormal(mesh)\n", " \n", " c = abs(inner(q / h, n)) + sqrt(g * h)\n", " α = avg(c)\n", " \n", " f_h = -α * (h('+') - h('-')) * (ϕ('+') - ϕ('-')) * dS\n", " f_q = -α * inner(q('+') - q('-'), v('+') - v('-')) * dS\n", "\n", " return f_h + f_q" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we'll do a few experiments with variable bottom topography as well." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def topographic_forcing(z, b):\n", " Z = z.function_space()\n", " h = firedrake.split(z)[0]\n", " v = firedrake.TestFunctions(Z)[1]\n", "\n", " return -g * h * inner(grad(b), v) * dx" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Time discretization\n", "\n", "We'll consider two different timestepping methods: the forward Euler scheme and the strong stability-preserving Runge-Kutta method of order 3, or SSPRK3 for short.\n", "Since we'll be using these quite a bit we'll factor them out into classes that store the needed internal state.\n", "\n", "In the previous demo on the conservative advection equation, we used solver parameters that specified a block ILU preconditioner for the linear system.\n", "One application of this preconditioner gives an exact solution to the linear system because the mass matrix for DG discretizations is block diagonal.\n", "We're doing something a little different here by using a mixed function space for the thickness and momentum because it makes the time discretization much easier.\n", "But as a consequence the mass matrix that Firedrake will build under the hood uses a [*nested*](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MATNEST.html) storage format -- probably compressed row storage for the thicknesses and block CRS for the momentum, but this shouldn't matter.\n", "In order to achieve the same exact linear solve effect here that we had for the conservative advection equation, we'll specify an outer-level [*fieldsplit*](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCFIELDSPLIT.html) preconditioner.\n", "Fieldsplit lets us separately specify preconditioners for each block, and those inner preconditioners will be our good old ILU + block Jacobi.\n", "You're likely to encounter fieldsplit preconditioners again if you ever have solved mixed finite element problems like the Stokes equations." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "block_parameters = {\n", " 'ksp_type': 'preonly',\n", " 'pc_type': 'ilu',\n", " 'sub_pc_type': 'bjacobi'\n", "}\n", "\n", "parameters = {\n", " 'solver_parameters': {\n", " 'ksp_type': 'preonly',\n", " 'pc_type': 'fieldsplit',\n", " 'fieldsplit_0': block_parameters,\n", " 'fieldsplit_1': block_parameters\n", " }\n", "}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "I've defined a time integrator class below that provides the bare minimum amount of functionality to do what we need.\n", "The data stored in this class consists of the current value of the state variables, an auxiliary state for the value at the next timestep, and a cached solver object for the mass matrix solve.\n", "The step method lets us advance the simulation by a single timestep which we can change on one invocation to the next.\n", "In principle we could do adaptive timestepping with this implementation." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "from firedrake import (\n", " NonlinearVariationalProblem as Problem,\n", " NonlinearVariationalSolver as Solver\n", ")\n", "\n", "class ForwardEuler:\n", " def __init__(self, state, equation):\n", " z = state.copy(deepcopy=True)\n", " F = equation(z)\n", " \n", " z_n = z.copy(deepcopy=True)\n", " Z = z.function_space()\n", " w = firedrake.TestFunction(Z)\n", " \n", " dt = firedrake.Constant(1.)\n", "\n", " problem = Problem(inner(z_n - z, w) * dx - dt * F, z_n)\n", " solver = Solver(problem, **parameters)\n", " \n", " self.state = z\n", " self.next_state = z_n\n", " self.timestep = dt\n", " self.solver = solver\n", " \n", " def step(self, timestep):\n", " self.timestep.assign(timestep)\n", " self.solver.solve()\n", " self.state.assign(self.next_state)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Demonstration\n", "\n", "Our first test problem will be a periodic domain with a flat bottom.\n", "The initial state of the system will consist of a spherical perturbation to the water thickness, and we'll look at how this disturbance evolves in time." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "nx, ny = 32, 32\n", "Lx, Ly = 20., 20.\n", "mesh = firedrake.PeriodicRectangleMesh(\n", " nx, ny, Lx, Ly, diagonal='crossed'\n", ")\n", "\n", "x = firedrake.SpatialCoordinate(mesh)\n", "lx = 5.\n", "y = Constant((lx, lx))\n", "r = Constant(2.5)\n", "\n", "h_0 = Constant(1.)\n", "δh = Constant(1/16)\n", "h_expr = h_0 + δh * max_value(0, 1 - inner(x - y, x - y) / r**2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### DG(0) basis; or, the finite volume method\n", "\n", "We'll start by using the simplest discretization possible: piecewise constant basis functions for both the thickness and momentum.\n", "This method is identical to the lowest-order finite volume method.\n", "We'll also use a *mixed* function space $Z = Q \\times V$ that includes both the thickness and momentum.\n", "This choice isn't strictly necessary but it makes it that much easier to write the time integrator you saw above.\n", "The code is equivalent to a single update for the combined state vector $z = (h, q)$ rather than two separate updates for each." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "Q0 = firedrake.FunctionSpace(mesh, family='DG', degree=0)\n", "V0 = firedrake.VectorFunctionSpace(mesh, family='DG', degree=0)\n", "Z0 = Q0 * V0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The split method of functions in mixed spaces give us the tuple of components.\n", "That way we can initialize the thickness to the expression defined just above.\n", "Note that the method split of functions in mixed spaces is different from the module-level function split, which gives us symbolic objects." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "z0 = firedrake.Function(Z0)\n", "h0, q0 = z0.split()\n", "h0.project(h_expr);" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "solver = ForwardEuler(\n", " z0,\n", " lambda z: (\n", " cell_flux(z) +\n", " central_facet_flux(z) +\n", " lax_friedrichs_facet_flux(z)\n", " )\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since we'll be running the same simulation many times, we'll wrap up the loop in a function." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "import tqdm\n", "\n", "def run_simulation(solver, final_time, num_steps, output_freq):\n", " hs, qs = [], []\n", " qs = []\n", " pbar = tqdm.trange(num_steps)\n", " for step in pbar:\n", " if step % output_freq == 0:\n", " h, q = solver.state.split()\n", " hmin, hmax = h.dat.data_ro.min(), h.dat.data_ro.max()\n", " pbar.set_description(f'{hmin:5.3f}, {hmax:5.3f}')\n", " hs.append(h.copy(deepcopy=True))\n", " qs.append(q.copy(deepcopy=True))\n", "\n", " solver.step(timestep)\n", " \n", " return hs, qs" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "1.000, 1.003: 100%|██████████| 2000/2000 [00:13<00:00, 147.03it/s]\n" ] } ], "source": [ "final_time = 10.\n", "timestep = 5e-3\n", "num_steps = int(final_time / timestep)\n", "output_freq = 10\n", "hs, qs = run_simulation(solver, final_time, num_steps, output_freq)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "from matplotlib.animation import FuncAnimation\n", "from IPython.display import HTML\n", "\n", "def make_animation(hs, b, timestep, output_freq, **kwargs):\n", " fig, axes = plt.subplots()\n", " axes.set_aspect('equal')\n", " axes.set_xlim((0.0, Lx))\n", " axes.set_ylim((0.0, Ly))\n", " axes.get_xaxis().set_visible(False)\n", " axes.get_yaxis().set_visible(False)\n", " η = firedrake.project(hs[0] + b, hs[0].function_space())\n", " colors = firedrake.tripcolor(\n", " hs[0], num_sample_points=1, axes=axes, **kwargs\n", " )\n", " \n", " def animate(h):\n", " η.project(h + b)\n", " colors.set_array(η.dat.data_ro[:])\n", "\n", " interval = 1e3 * output_freq * timestep\n", " animation = FuncAnimation(fig, animate, frames=hs, interval=interval)\n", " \n", " plt.close(fig)\n", " return HTML(animation.to_html5_video())" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "make_animation(\n", " hs, Constant(0), timestep, output_freq, vmin=0.98, vmax=1.03\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We get the expected propagation at first, but the method is so diffusive that the waves quickly get washed out.\n", "Another sanity check that we'll repeat through the following is to track the total energy of the system.\n", "The shallow water equations conserve the quantity\n", "\n", "$$E = \\frac{1}{2}\\int_\\Omega\\left(\\frac{|q|^2}{h} + g(h + b)^2\\right)dx.$$\n", "\n", "(There's a Hamiltonian formulation of the shallow water equations based on this energy functional, although the Poisson bracket is a little weird.)\n", "Approximately conserving the total energy is especially important for long-running simulations of systems where physical dissipation mechanisms are relatively weak or non-existent.\n", "The plot below shows that the explicit Euler scheme with DG(0) basis functions and the Lax-Friedrichs flux dissipates energy." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "energies = [\n", " firedrake.assemble(\n", " 0.5 * (inner(q, q) / h + g * h**2) * dx\n", " )\n", " for h, q in zip(hs, qs)\n", "]" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEDCAYAAAA7jc+ZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXhc9X3v8fdXGo32zfJ4k3dsAzYGDMJAylZoiU0SnBAI0LRAy700aehtmyYNPLk3zSX3yQ0pCZckNIQGCiELIWRzHpwYypqmYCyDsTFekI1XbEuWZcvat+/9Y47s8SBZIyTNjGc+r+fRM2d+5zcz3zkafc7R75w5x9wdERHJXDmpLkBERMaWgl5EJMMp6EVEMpyCXkQkwynoRUQynIJeRCTDnXRBb2bXmdkGM+szs5oT9Ps7M3sz6Pv3cfP+1sw2BfO+HrTlmdmjZrbezDaa2Z0J1PI7M3sjeJ4HzCx35O9QRGR0pXXQm9llZvZIXPObwDXASyd43BnAfwcWA2cBHzazOcG8PwaWAWe5+wLgnuBh1wH57r4QOBf4azObOUSJn3D3s4AzgEjwHCIiaSWtg34g7r7R3TcP0e10YJW7t7l7D/Ai0ZUDwKeBr7l7Z/B89f1PDRSbWQgoBLqAZgAz+3Mze9XM1prZ9/q33N29OXhsCAgHzyEiklZOuqBP0JvAxWZWZWZFwFXAtGDevGDeKjN70czOC9qfBFqBvcBO4B53P2hmpwPXA3/k7mcDvcAn+1/IzFYC9cCR4DlERNJKKNUFDMTMVgH5QAkwzszWBrO+4O4rh3q8u280s7uBp4mG91qiAQ3R9zwOuAA4D3jCzGYTHebpBaYAlcDvzew/gCuIDuWsNjOIbu3Xx7zWB82sAPgRcDnwzAjeuojIqEvLoHf38yE6Rg/c4u63vI/neAh4KHierwK7g1m7gV949CQ/r5pZHzAe+DPgd+7eDdSb2R+AGsCAR9190J2z7t5hZr8mOvavoBeRtJKpQzeY2YTgdjrR8fkfB7N+BfxxMG8e0bH1A0SHay4P2ouJbvFvAp4Fro15vnFmNsPMSsxsctAWAj4U9BcRSSsnXdCb2cfMbDdwIfBUMEaOmU0xsxUxXX9uZm8BvwE+4+6HgvaHgdlm9ibwOHBzsHV/P1BiZhuA1cC/u/s6d38L+J/A02a2jugW+2SgGFgetK0lOpzzwNi+exGR4TOdplhEJLOddFv0IiIyPGm3M3b8+PE+c+bMVJchInJSWbNmzQF3jww0L+2CfubMmdTW1qa6DBGRk4qZ7RhsnoZuREQynIJeRCTDKehFRDKcgl5EJMMlFPRmtsTMNptZnZndMcD8S8zsNTPrMbNrB5hfZma7zew7o1G0iIgkbsigD07Jez+wFJgP3Ghm8+O67QRu4dhpBuJ9hROcP15ERMZOIlv0i4E6d9/m7l1ETxuwLLaDu29393VAX/yDzexcYCLRM0mKiEiSJRL01cCumPu7g7YhmVkO8A3gc0P0u83Mas2stqGhIZGnfo/mjm7ufWYLa3cdGrqziEgWGeudsX8DrHD33Sfq5O4PunuNu9dEIgN+sWtIfX3Ofc++zZodTe/r8SIimSqRb8bu4djVmQCmBm2JuJDo1Zz+huhFRMJm1uLu79mhO1JlBXnkGBxq6xrtpxYROaklEvSrgblmNotowN9A9CIdQ3L32Evu3QLUjEXIA+TkGBVFYQ62KuhFRGINOXQTXFz7dmAlsBF4wt03mNldZnY1gJmdF5wj/jrge8E53ZOuoiiPQ23dqXhpEZG0ldBJzdx9BbAiru1LMdOriQ7pnOg5HgEeGXaFw1BZFKZJQzciIsfJqG/GVhbl0aQtehGR42RY0Idp0hi9iMhxMivoizV0IyISL6OCvqIoj86ePtq7elNdiohI2siooK8sCgNoq15EJEZGBr2OpRcROSbDgj4PQMfSi4jEyKygL9bQjYhIvIwK+oqjW/QKehGRfpkV9IX9Y/QauhER6ZdRQR8O5VCaH9LQjYhIjIwKeoCK4jwN3YiIxMi4oI+e2ExDNyIi/TIu6Ct0BksRkeNkXNCPK8pT0IuIxMi4oK8oCtOko25ERI7KuKCPlObT0tlDR7dObCYiApkY9CX5ADQc6UxxJSIi6SHzgr4sGvT1CnoRESATg15b9CIix8m4oJ9QGgR9i4JeRAQyMOjHFYcx0xa9iEi/jAv6UG4OVcVhBb2ISCDjgh5gfEk+DUc6Ul2GiEhayMign1BWoC16EZFAQkFvZkvMbLOZ1ZnZHQPMv8TMXjOzHjO7Nqb9bDN72cw2mNk6M7t+NIsfTKQkX0EvIhIYMujNLBe4H1gKzAduNLP5cd12ArcAP45rbwNucvcFwBLg/5lZxUiLHkqkNJ+Glk7cfaxfSkQk7SWyRb8YqHP3be7eBTwOLIvt4O7b3X0d0BfXvsXd3w6m3wXqgcioVH4CkdJ8unudw+06542ISCJBXw3sirm/O2gbFjNbDISBrQPMu83Mas2stqGhYbhP/R6RUn1pSkSkX1J2xprZZOAx4C/dvS9+vrs/6O417l4TiYx8g1/fjhUROSaRoN8DTIu5PzVoS4iZlQFPAV9091eGV977M0HnuxEROSqRoF8NzDWzWWYWBm4Alify5EH/XwI/cPcn33+Zw9M/dFOvY+lFRIYOenfvAW4HVgIbgSfcfYOZ3WVmVwOY2Xlmthu4DviemW0IHv4J4BLgFjNbG/ycPSbvJEZpfojCvFz2HdYWvYhIKJFO7r4CWBHX9qWY6dVEh3TiH/dD4IcjrHHYzIzJFQXsPdye7JcWEUk7GfnNWIAp5YW8e1hDNyIiGRv0k8sL2HtIW/QiIpkb9BWFNLR00tXznqM5RUSySsYG/ZTyAtxhf7OGb0Qku2Vs0E+uKARgr8bpRSTLZWzQTykvANCRNyKS9TI26Pu36N89pC16EcluGRv0JfkhSgtC2qIXkayXsUEPwbH02qIXkSyX0UGvb8eKiGR60JcX6qgbEcl6GR301RUFHGztoq2rJ9WliIikTEYH/bRxRQDsbtLwjYhkr4wO+hlVxQDsbGxLcSUiIqmT0UE/Pdii33lQQS8i2Sujg76yKI+S/JCCXkSyWkYHvZkxbVwRuxT0IpLFMjroAaaPK9QWvYhktSwI+iJ2HmzD3VNdiohISmRF0Hf29NFwRBcKF5HslPFB338s/Q4N34hIlsr4oD96iKWOpReRLJXxQV9dWYiZtuhFJHtlfNDnh3Kprihk+4HWVJciIpISGR/0AKdEStja0JLqMkREUiKhoDezJWa22czqzOyOAeZfYmavmVmPmV0bN+9mM3s7+Ll5tAofjtmRYrY1tNLXp0MsRST7DBn0ZpYL3A8sBeYDN5rZ/LhuO4FbgB/HPXYc8M/A+cBi4J/NrHLkZQ/PKZES2rt72desc9OLSPZJZIt+MVDn7tvcvQt4HFgW28Hdt7v7OqAv7rEfBJ5x94Pu3gQ8AywZhbqHZXYkehbLbQ0apxeR7JNI0FcDu2Lu7w7aEpHQY83sNjOrNbPahoaGBJ86cXMiJQAapxeRrJQWO2Pd/UF3r3H3mkgkMurPHynNpyQ/xDYFvYhkoUSCfg8wLeb+1KAtESN57KgxM06JFLNVQzcikoUSCfrVwFwzm2VmYeAGYHmCz78SuNLMKoOdsFcGbUk3W4dYikiWGjLo3b0HuJ1oQG8EnnD3DWZ2l5ldDWBm55nZbuA64HtmtiF47EHgK0RXFquBu4K2pDslUszewx20dupC4SKSXUKJdHL3FcCKuLYvxUyvJjosM9BjHwYeHkGNo2LuxFIAtuw/wqLpST/CU0QkZdJiZ2wynD6pDIDN+46kuBIRkeTKmqCfWllIUTiXTQp6EckyWRP0OTnGvImlbNrXnOpSRESSKmuCHuD0yaVs3ndElxUUkaySVUF/6sRSmtq6qddlBUUki2RX0Ac7ZDVOLyLZJKuC/rRJ0UMsN2ucXkSySFYFfWVxmEllBbz1roJeRLJHVgU9wBnV5azfczjVZYiIJE3WBf3C6nK2HWilRadCEJEskXVBf+bUctxhg7bqRSRLZF3Qn1FdDqDhGxHJGlkX9JHSfCaVFSjoRSRrZF3QAyycqh2yIpI9sjPoq8vZ1tDKkY7uVJciIjLmsjLoz5pWAcAbu7RVLyKZLyuDftH0CsxgzY6mVJciIjLmsjLoywryOHViKWt2KuhFJPNlZdADnDOjktd3NNHXp1MWi0hmy9qgP3d6JUc6e3i7viXVpYiIjKnsDfoZ0QuEv6bhGxHJcFkb9DOqiqgqDrN6+8FUlyIiMqayNujNjMWzxrFq20FdWlBEMlrWBj3AhadUsedQO7sOtqe6FBGRMZNQ0JvZEjPbbGZ1ZnbHAPPzzeynwfxVZjYzaM8zs0fNbL2ZbTSzO0e3/JG5cHYVAC9vO5DiSkRExs6QQW9mucD9wFJgPnCjmc2P63Yr0OTuc4B7gbuD9uuAfHdfCJwL/HX/SiAdzJlQwviSMC9vbUx1KSIiYyaRLfrFQJ27b3P3LuBxYFlcn2XAo8H0k8AVZmaAA8VmFgIKgS4gba7jZ2acP7uKVzROLyIZLJGgrwZ2xdzfHbQN2Mfde4DDQBXR0G8F9gI7gXvcPa0Oc7lwdhX7mjt450BrqksRERkTY70zdjHQC0wBZgH/aGaz4zuZ2W1mVmtmtQ0NDWNc0vEumRsB4KUtyX1dEZFkSSTo9wDTYu5PDdoG7BMM05QDjcCfAb9z9253rwf+ANTEv4C7P+juNe5eE4lEhv8uRmB6VRGzxhfzgoJeRDJUIkG/GphrZrPMLAzcACyP67McuDmYvhZ4zqOD3juBywHMrBi4ANg0GoWPpkvnRXhlWyMd3b2pLkVEZNQNGfTBmPvtwEpgI/CEu28ws7vM7Oqg20NAlZnVAZ8F+g/BvB8oMbMNRFcY/+7u60b7TYzUpadG6OjuY9U7abX7QERkVIQS6eTuK4AVcW1fipnuIHooZfzjWgZqTzcXzKoiHMrhhc31XDovuUNHIiJjLau/GduvMJzLhbOreG5TvQ6zFJGMo6APXLlgIjsa29iyX6ctFpHMoqAP/OnpEwF4esO+FFciIjK6FPSBCWUFLJpewTMb96e6FBGRUaWgj3Hl/Ems232Ydw/pbJYikjkU9DGWnDEJgBXr96a4EhGR0aOgjzFrfDFnVJfxm3UKehHJHAr6OB85cwpv7DrEzsa2VJciIjIqFPRxPnzWFAB+s+7dFFciIjI6FPRxqisKqZlRya9e36MvT4lIRlDQD+Cac6bydn0L63YfTnUpIiIjpqAfwIfPmkx+KIefrdk1dGcRkTSnoB9AWUEeS86YxPK17+rUxSJy0lPQD+Lac6fS3NHDSp0SQUROcgr6QfzRKeOZPq6IH72yM9WliIiMiIJ+EDk5xp9fMJ1Xtx9k077mVJcjIvK+KehP4LpzpxEO5fDDV3akuhQRkfdNQX8ClcVhPnLmFH6+Zg9NrV2pLkdE5H1R0A/htktm097dyw9e1la9iJycFPRDOHVSKZefNoFHX95Oe5cOtRSRk4+CPgGfvuwUDrZ28UStvkAlIicfBX0Czps5jnNnVPLgS9vo7u1LdTkiIsOioE/Qpy49hT2H2nlK56oXkZOMgj5BV5w2gbkTSvjuC1vp69NZLUXk5KGgT1BOjvG3V8xl8/4j/GrtnlSXIyKSsISC3syWmNlmM6szszsGmJ9vZj8N5q8ys5kx8840s5fNbIOZrTezgtErP7k+vHAyC6vL+cbTW3SyMxE5aQwZ9GaWC9wPLAXmAzea2fy4brcCTe4+B7gXuDt4bAj4IfApd18AXAZ0j1r1SZaTY9y59DT2HGrnBy9vT3U5IiIJSWSLfjFQ5+7b3L0LeBxYFtdnGfBoMP0kcIWZGXAlsM7d3wBw90Z3P6k3hT8wZzyXzovwnefqONSmb8uKSPpLJOirgdgDyHcHbQP2cfce4DBQBcwD3MxWmtlrZvZPA72Amd1mZrVmVtvQ0DDc95B0dyw9jSOdPfzrC1tTXYqIyJDGemdsCLgI+GRw+zEzuyK+k7s/6O417l4TiUTGuKSRO31yGR8/ZyqP/GE7dfUtqS5HROSEEgn6PcC0mPtTg7YB+wTj8uVAI9Gt/5fc/YC7twErgHNGWnQ6+MKS0yjIy+GLv1yvi4iLSFpLJOhXA3PNbJaZhYEbgOVxfZYDNwfT1wLPeTT9VgILzawoWAFcCrw1OqWnVqQ0nzuvOp1V7xzkZ2t2p7ocEZFBDRn0wZj77URDeyPwhLtvMLO7zOzqoNtDQJWZ1QGfBe4IHtsEfJPoymIt8Jq7PzX6byM1rq+ZRs2MSr66YiONLZ2pLkdEZECWbsMONTU1Xltbm+oyEvb2/iNc9a3f85Ezp/DN689OdTkikqXMbI271ww0T9+MHaG5E0v59KWn8IvX9/C7N3UeHBFJPwr6UXD75XNZWF3OF36+nr2H21NdjojIcRT0oyAcyuG+G86mu7ePf/jpWnp10jMRSSMK+lEyO1LCl69ewCvbDvLAi/oilYikDwX9KLru3Kl8aOFkvvnMFl7e2pjqckREAAX9qDIz/u/HFzKzqojP/Pg1dh1sS3VJIiIK+tFWVpDHv91UQ3dvH7c9toa2rp5UlyQiWU5BPwZmR0r49o2L2Lyvmc//bJ2uSCUiKaWgHyOXnTqBO5eezlPr9/LVFRt1PhwRSZlQqgvIZP/t4lnsOdTO9//zHapK8vn0ZaekuiQRyUIK+jFkZnzpw/M52NrF3b/bxLjiPK4/b3qqyxKRLKOgH2M5OcY9153FofZu7vzFevJDuXx0Ufx1W0RExo7G6JMgHMrhgT8/h/NnVfEPT6zlSZ3WWESSSEGfJEXhEA/fch4XzRnP5598g5+8ujPVJYlIllDQJ1FhOJd/u6mGy+ZFuPMX6/nuC1t1NI6IjDkFfZIV5OXywF+cy0fOmsLdv9vE//r1mzoJmoiMKe2MTYH8UC73XX82U8oL+N5L29h3uIP7blhEcb5+HSIy+rRFnyI5OcadV53OXcsW8Nymeq751/9i+4HWVJclIhlIQZ9iN104k0f+cjH7mju4+jv/yfOb61NdkohkGAV9GrhkXoTf3H4R1ZVF/NUjq7ln5Wa6e/tSXZaIZAgFfZqYXlXELz79Aa49Zyrfeb6OT3zvZZ3mWERGhYI+jRSGc/mX687iWzcuom5/C1fd93t+VrtLh2CKyIgo6NPQ1WdNYcXfXczpk8v4/JPruOnhV7V1LyLvm4I+TU0bV8Tjt13AVz56Bq/taOLKe1/i+7/fpmPuRWTYFPRpLCfH+IsLZvDMZy/lwlOq+D9PbeSaf/0Da3Y0pbo0ETmJJBT0ZrbEzDabWZ2Z3THA/Hwz+2kwf5WZzYybP93MWszsc6NTdnaZUlHIQzfXcN8NZ7P3cAcf/+5/8T9+8jp7DrWnujQROQkMGfRmlgvcDywF5gM3mtn8uG63Ak3uPge4F7g7bv43gd+OvNzsZWYsO7ua5z93GX97+RxWbtjH5fe8wDee3kxrp65LKyKDS2SLfjFQ5+7b3L0LeBxYFtdnGfBoMP0kcIWZGYCZfRR4B9gwOiVnt+L8EP945ak897nL+OCCSXz7uTou/vrzPPDiVl2IXEQGlEjQVwO7Yu7vDtoG7OPuPcBhoMrMSoAvAP/7RC9gZreZWa2Z1TY0NCRae1arrijkWzcu4pd/8wHOqC7na7/dxMV3P8+DL22lvas31eWJSBoZ652xXwbudfeWE3Vy9wfdvcbdayKRyBiXlFkWTa/kB3+1mJ9/+kLmTynjqys2cfHXn+Pbz75NU2tXqssTkTSQyOkS9wDTYu5PDdoG6rPbzEJAOdAInA9ca2ZfByqAPjPrcPfvjLhyOc65M8bx2K3nU7v9IN9+ro5vPLOF+1+o49pzp3LrRbOZNb441SWKSIrYUN+6DIJ7C3AF0UBfDfyZu2+I6fMZYKG7f8rMbgCucfdPxD3Pl4EWd7/nRK9XU1PjtbW17+e9SIwt+4/w/d9v41evv0t3Xx9XnDaBT54/g0vmRcjNsVSXJyKjzMzWuHvNQPOG3KJ39x4zux1YCeQCD7v7BjO7C6h19+XAQ8BjZlYHHARuGL3y5f2YN7GUr197Fp/74Kk89vIOfvLqTv5jYz3VFYXcuHgan6iZxoSyglSXKSJJMOQWfbJpi35sdPX08cxb+/nxqzv4Q10joRzj0nkRPrqomj85fSKF4dxUlygiIzCiLXrJDOFQDh86czIfOnMy7xxo5fHVO/n16+/y7KZ6SvJDLDljEh9bVM0Fs6s0tCOSYbRFn8V6+5xV7zTyq9f38Nv1+zjS2cPEsnyWnjGZKxdMZPHMcYRydZYMkZPBibboFfQCQEd3L89urOdXa/fw0pYGOnv6qCzK44rTJ/LBBZO4eO54CvI0vCOSrhT0MiytnT28tKWBlRv28eymeo509FAUzuWiOeO59NQIl8yNMG1cUarLFJEYGqOXYSnOD7F04WSWLpxMV08fr2xrZOWGfbywuYGn39oPwOxIMZfOi3DJvAgXzKrSzlyRNKYtekmYu7O1oZUXtzTw0pYGXtnWSGdPH3m5xplTK7hg9jjOn1XFuTMqKc7XNoRIMmnoRsZER3cvq945yH9tPcCqbQdZv+cwvX1OKMdYOLWc82dVcf7scZwzvZLywrxUlyuS0RT0khQtnT2s2dHEqm2NrHrnIOt2H6K7N/r5OiVSzNnTKjl7egWLplVw2qRSHdEjMoo0Ri9JUZIf4tJ5ES6dFz0xXXtXL6/vbOL1XYd4fWcTL26p5+ev7QagIC+HM6srOGtaOQumlLNgShmzxhcr/EXGgIJexkxhOJcPzBnPB+aMB6Jj/Lub2o8G/9pdh3j05R109fQBkB/K4bTJZcyfXMaCKWXMn1LG6ZPKtKNXZIQ0dCMp1d3bx7aGVja8e5i33m1mw7vNvLW3mcPt3QDkGMwcX8ypE0uZO6GEOcHtrPHFOq5fJIaGbiRt5eXmcOqkUk6dVMo150Tb3J09h9qPBv/Gvc1s3neElRv20Rdsl+QYzKgqZs6EEuZOKGHuxBJmjS9hZlURFUXh1L0hkTSkoJe0Y2ZMrSxiamURVy6YdLS9s6eXdw608vb+Ft6ub6Gu/ghv72/h+U319PQd+8+0oiiPGVXFzKoqit6OL2bm+GKtBCRrKejlpJEfyuW0SWWcNqnsuPaunj52NLayvbGN7QdaeaexlR2Nraze3sSv33iX2NHJiqI8ZowrClYkhVRXFlJdUcjUyiKqKwsp0fH/koH0qZaTXjiUw9yJpcydWPqeeR3dvew62HZ0JbC9sZUdjW28tbeZZzbuP7ojuF9FUV4Q/IVUVxQdXRFMLi9gUnkB40vydXZPOeko6CWjFeTlDroS6OtzDrR0svtQO7ub2tnT1M7upjb2HGpnW0MrL205QHv38RdazzGIlOYzqayAicHPpPL+6aC9vIDS/BBmWiFIelDQS9bKyTEmlBUwoayAc6ZXvme+u9PU1s3upjb2N3eyr7mD+uYO9h3uYF9zBzsa21j1zsGjRwjFKszLZWJZPuNLoj9VJeFgOhzcj05XleRTVqCVgowtBb3IIMyMccVhxhWfeAdue1cv9UeOrQDqg5XC/uYOGlu62NrQwqvbu2hq62Kgo5nDuTlHVwSxt+OKwlQU5VFRFKayKExlMF1RlEeevlgmw6CgFxmhwnAuM6qKmVFVfMJ+Pb19HGzr4sCRLhpbOznQ0kljSxcNwW3//c37jtDY0kVXb9+gz1WaH6K8KI/KIPhjVwSxK4SywjzKCvIoKwhRVphHfihH/z1kIQW9SJKEcnOYUFrAhNKhL8ru7rR29dLU2sXh9m6a2rpoauvmUFsXTa3R+7HtOw+20dTaRXNHzwmfN5ybQ2kQ+mUFIUoL8igrDFFWkBdtL4iuHPqn+/uW5Icozg9RnJ9LfkhfVDvZKOhF0pCZUZIfoiQ/xLRhPK6nt4/mjh6a2ro41BYN/ub2bo509NDc0U1zew9HOrpp7ghu27vZ19wRTPe8Z+fzQPJyLRr64WjwFwd1FoWPTUfn5wYrh2PzS+LuF4ZzKQjlkqMjmcaUgl4kg4RycxLarzCY7t6+6EqhvZvmju6j0y2dPbR29tDa1XtsurM3aOuhpbOH/c0d0bau6Pz+M5cmoiAvh8K8XIrCoeh0OJfCvFwKwyEKg3nR6VwKwwPdDx17TNCWH8olP5RDfl5wm8XDVgp6ETkqb4QrilidPb3HrQxaO3to6b8f/HT09NHW1UtHdy/tXb20x90ebu9m/+HgftDW1tVD3/s8RVc4CPz+lUBBXjCdF9/ev5KIWWEE/QpiVh7hUA7h3BzyQjnk50bv58Xc5h933472T/YKR0EvImMiGpC5o7LSiOXudPf6sRXC0ZVDD+1dfUfbOrt76ezpC3566ew+Nt3RHbT19AXt0elD7d10dvfS1dNHR9zjh/MfylDycu3oCiKce2ylsKC6nG/fuGjUXqefgl5ETipmRjgU3TouJ3lXLuvtc7p6jl9BdPREVwrdvX3BrdPV20tXj9N1tO3YbWfc/a6ePrp6PbjtY/q4wjGpPaGgN7MlwH1ALvB9d/9a3Px84AfAuUAjcL27bzezPwW+BoSBLuDz7v7cKNYvIpIUuTkW3Q9wEl4fYchvXZhZLnA/sBSYD9xoZvPjut0KNLn7HOBe4O6g/QDwEXdfCNwMPDZahYuISGIS+XrdYqDO3be5exfwOLAsrs8y4NFg+kngCjMzd3/d3d8N2jcAhcHWv4iIJEkiQV8N7Iq5vztoG7CPu/cAh4GquD4fB15z9874FzCz28ys1sxqGxoaEq1dREQSkJQTZpjZAqLDOX890Hx3f9Dda9y9JhKJJKMkEZGskUjQ74Hjvpw3NWgbsI+ZhYByojtlMbOpwC+Bm9x960gLFhGR4Ukk6FcDc81slpmFgRuA5XF9lhPd2QpwLfCcu7uZVQBPAXe4+x9Gq2gREUnckEEfjLnfDqwENgJPuPsGM7vLzK4Ouj0EVJlZHfBZ4I6g/XZgDvAlM1sb/EwY9XchIiKDMh/oBNkpVFNT47W1takuQ0TkpGJma9y9ZsB56Rb0ZtYA7BjBU4wnemhVvRUAAASuSURBVPx+ulFdw5OudUH61qa6hidd64L3V9sMdx/waJa0C/qRMrPawdZqqaS6hidd64L0rU11DU+61gWjX5uuRyYikuEU9CIiGS4Tg/7BVBcwCNU1POlaF6RvbapreNK1Lhjl2jJujF5ERI6XiVv0IiISQ0EvIpLhMibozWyJmW02szozu2PoR4xZHdPM7Hkze8vMNpjZ3wXtXzazPTHfEL4qRfVtN7P1QQ21Qds4M3vGzN4ObiuTXNOpMctlrZk1m9nfp2KZmdnDZlZvZm/GtA24fCzqW8Fnbp2ZnZPkuv7FzDYFr/3L4JQjmNlMM2uPWW4PjFVdJ6ht0N+dmd0ZLLPNZvbBJNf105iatpvZ2qA9acvsBBkxdp8zdz/pf4he+WorMJvo1azeAOanqJbJwDnBdCmwhegFW74MfC4NltV2YHxc29eJno8IoqevuDvFv8t9wIxULDPgEuAc4M2hlg9wFfBbwIALgFVJrutKIBRM3x1T18zYfilaZgP+7oK/hTeAfGBW8Hebm6y64uZ/A/hSspfZCTJizD5nmbJFn8jFUZLC3fe6+2vB9BGi5weKP39/uom9cMyjwEdTWMsVwFZ3H8m3o983d38JOBjXPNjyWQb8wKNeASrMbHKy6nL3pz16LiqAV4ieWTbpBllmg1kGPO7une7+DlBH9O83qXWZmQGfAH4yFq99IifIiDH7nGVK0CdycZSkM7OZwCJgVdB0e/Cv18PJHh6J4cDTZrbGzG4L2ia6+95geh8wMTWlAdGzo8b+8aXDMhts+aTT5+6viG719ZtlZq+b2YtmdnGKahrod5cuy+xiYL+7vx3TlvRlFpcRY/Y5y5SgTztmVgL8HPh7d28GvgucApwN7CX6b2MqXOTu5xC9BvBnzOyS2Jke/V8xJcfcWvQ02FcDPwua0mWZHZXK5TMYM/si0AP8KGjaC0x390VEzyb7YzMrS3JZafe7i3Mjx29QJH2ZDZARR4325yxTgj6Ri6MkjZnlEf0F/sjdfwHg7vvdvdfd+4B/Y4z+XR2Ku+8JbuuJXhBmMbC//1/B4LY+FbURXfm85u77gxrTYpkx+PJJ+efOzG4BPgx8MggHgmGRxmB6DdFx8HnJrOsEv7t0WGYh4Brgp/1tyV5mA2UEY/g5y5SgT+TiKEkRjP09BGx092/GtMeOqX0MeDP+sUmordjMSvunie7Me5PjLxxzM/DrZNcWOG4rKx2WWWCw5bMcuCk4KuIC4HDMv95jzsyWAP8EXO3ubTHtETPLDaZnA3OBbcmqK3jdwX53y4EbzCzfzGYFtb2azNqAPwE2ufvu/oZkLrPBMoKx/JwlYy9zMn6I7pneQnRN/MUU1nER0X+51gFrg5+rgMeA9UH7cmByCmqbTfSIhzeADf3LieiF3J8F3gb+AxiXgtqKiV5+sjymLenLjOiKZi/QTXQs9NbBlg/RoyDuDz5z64GaJNdVR3Tstv9z9kDQ9+PB73ct8BrwkRQss0F/d8AXg2W2GViazLqC9keAT8X1TdoyO0FGjNnnTKdAEBHJcJkydCMiIoNQ0IuIZDgFvYhIhlPQi4hkOAW9iEiGU9CLiGQ4Bb2ISIb7/92uIBcxKnKhAAAAAElFTkSuQmCC\n", "text/plain": [ "