diff --git a/doc/notebooks/demo_central_moments.ipynb b/doc/notebooks/demo_central_moments.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..c5997429f0eb6903528875a5ef4a2f3eb5d68726 --- /dev/null +++ b/doc/notebooks/demo_central_moments.ipynb @@ -0,0 +1,1069 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "# imports\n", + "from lbmpy.session import *\n", + "from lbmpy.moments import *\n", + "\n", + "from lbmpy.forcemodels import Guo\n", + "from lbmpy.methods.abstractlbmethod import LbmCollisionRule\n", + "\n", + "from collections import OrderedDict\n", + "from lbmpy.methods.creationfunctions import create_with_discrete_maxwellian_eq_moments\n", + "\n", + "from functions import get_shift_matrix, get_central_moments" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Tutorial: Central Moments\n", + "\n", + "In this tutorial, we will learn what central moments are, what advantages they have and how we can use them with lbmpy." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1) Theoretical Background" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### What are Central Moments and why use them?\n", + "\n", + "The most common collision operators are the SRT and MRT operators. Both of these operators are based on a second-order approximation of the Maxwell-Boltzmann distribution, where the relaxation is done in either single or multiple steps. However, both approaches have the disadvantage that at relaxation rates that are close to the limit, both methods become unstable. For these kinds of problems, other operators can be used, e.g. the entropic LB operator, which enforces an H-theorem on the lattice and is therefore unconditionally stable. In order to get this feature of unconditional stability, the relaxation time is modulated in dependence of the local entropy, which removes high frequencies from the flux fields. Nevertheless, the method that is presented in this notebook is another one, the Cascaded LBM (= CLBM), which uses central moments. It does not remove any frequencies, but instead tries to deal with the high frequencies with sufficient accuracy, such that these do not threaten the overall stability." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "So CLBM has the goal to remove instabilities which can be traced back to an insufficient level of Galilean invariance (GI). GI describes the independence of physical processes from the speed of the reference system." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "From tutorial 3 we know, that _lbmpy_ is using the following equation for moment based relaxations\n", + "\n", + "$$K(f) = f - C^{-1}S\\left(Cf - c^{(eq)}\\right)$$\n", + "\n", + "A few things are happening here:\n", + "* $Cf$ is transforming the distribution function $f$ into the moment space. The solution of this linear transformation is the moment-vector $c = Cf$.\n", + "* $c^{(eq)}$ defines the equilibrium in the moment space (alternatively $Cf^{(eq)}$, but it is more natural to define it directly in the moment space)\n", + "* In the moment space, the relaxation is executed by multiplying $S$ from left, to get $\\Delta m = S \\cdot \\left(m - m^{(eq)}\\right)$\n", + "* then an inverse transformation of the relaxed moments is done by $\\Delta f = C^{-1} \\Delta m$\n", + "* and in the end, the streaming is happening by $K(f) = f - \\Delta f$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now, in order to restore the Galilean invariance, we will use central moments. For that, we have to perform an additional transformation; we have to transform the raw moments into the central moment space. First, we define the continuous central moments:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\\kappa_{\\alpha \\beta \\gamma} = \\int_{\\mathbb{R}^3} (\\xi_x - u_x)^\\alpha \\cdot (\\xi_y - u_y)^\\beta \\cdot (\\xi_z - u_z)^\\gamma \\cdot f(t, x, \\xi) \\,\\, d\\xi$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "where $u = \\left( u_x, u_y, u_z \\right)$ describes the macroscopic velocities. The discrete central moments are defined as" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\\kappa_{\\alpha \\beta \\gamma} = \\sum_{i = 0}^{n_v} [(c_i)_x - u_x]^\\alpha \\cdot [(c_i)_y - u_y]^\\beta \\cdot [(c_i)_z - u_z]^\\gamma \\cdot f_i \\\\\n", + "\\qquad\\qquad\\quad = \\sum_{i, j, k = -1}^{1} [(c_{ijk})_x - u_x]^\\alpha \\cdot [(c_{ijk})_y - u_y]^\\beta \\cdot [(c_{ijk})_z - u_z]^\\gamma \\cdot f_{ijk} $$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The main difference between the transformation to raw moments and the one to central moments is that the first is a linear transformation, while the second is a polynomial transformation. As a consequence, we cannot easily formulate this transformation as a simple matrix-vector multiplication; nevertheless, we can use the fast central moment transformation:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$ \\kappa_{ij\\,\\mid\\,\\gamma} = \\sum_{k = -1}^{1} [(c_{ijk})_z - u_z]^\\gamma \\cdot f_{ijk} $$\n", + "$$ \\kappa_{i\\,\\mid\\,\\beta \\gamma} = \\sum_{j = -1}^{1} [(c_{ijk})_y - u_y]^\\beta \\cdot \\kappa_{ij\\,\\mid\\,\\gamma} $$\n", + "$$ \\kappa_{\\alpha \\beta \\gamma} = \\sum_{i = -1}^{1} [(c_{ijk})_x - u_x]^\\alpha \\cdot \\kappa_{i\\,\\mid\\, \\beta \\gamma} $$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we can replace the transformation from the distribution function into the raw moment space with the fast central moment transformation into the central moment space, but because it is not a linear transformation, we are not able to write this as one equation as before, at least not directly. Therefore we have to do computations step by step and put the results together, such that we get the following algorithm:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "* Transform distribution function $f_{ijk}$ with the fast central moment transformation into the central moments $\\kappa_{\\alpha \\beta \\gamma}$ $\\left(f_{ijk} \\Rightarrow \\kappa_{\\alpha\\beta\\gamma}\\right)$\n", + "\n", + "* Calculate the central equilibrium moments $\\kappa_{i}^{(eq)}$ $\\left(f_{ijk}^{(eq)} \\Rightarrow \\kappa_{\\alpha\\beta\\gamma}^{(eq)}\\right)$ \n", + "\n", + "* Relax the central moments $\\kappa_i$ each with the relaxation frequency $\\omega_i$ (diagonal elements of $S$) $\\left(\\kappa^{*} = S \\left( \\kappa_{\\alpha\\beta\\gamma} - \\kappa_{\\alpha\\beta\\gamma}^{(eq)} \\right)\\right)$\n", + "\n", + "* Transform the relaxed central moments with the inverse central moments transformation into the relaxed distribution function $f_{ijk}^{*}$ $\\left(\\kappa^{*} \\Rightarrow f^{*}\\right)$\n", + "\n", + "* Do the update $\\left(f(t+\\Delta t, x+c_x \\Delta t) = f(t, x) - f^{*}(t, x)\\right)$ and streaming step " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Possible Solutions to Fix the Problem" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "However, there are two possibilities on how to solve the dilemma and write this again in one equation. One is to modify the above-described approach, and the other one uses features of lbmpy for a workaround: \n", + "1. Extract with the help of the fast central moment transformation a shift matrix, which is lower triangular, to transform the moment space into the central moment space. \n", + "2. Set up a method directly with the central moments. However, the problem with this approach is, that the transformation matrix which brings the distribution function into the moment space changes, and is not sparse anymore but full. As a consequence, the subexpression elimination is not able to reduce the operation count as much as in the case of the shift matrix (D2Q9, in D3Q27 case the subexpression elimination is not working anymore), which we will later see in an example.\n", + "\n", + "Let us first look at the first possibility. We introduce a little modification to the approach described above. With the help of the fast central moment transformation, we can extract a shift matrix $N$, which transforms the raw moments into the central moment space via matrix multiplication. The shift matrix consists of the coefficients of the moments from the polynomials, which we get from the fast central moment transformation. Furthermore, with this, we can again write the collision equation in a form with matrix-matrix/vector multiplications and do not have to execute multiple things step after step." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$K(f) = f - C^{-1}N^{-1}SN\\left(Cf - c^{(eq)}\\right)$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2) Using Central Moments with *lbmpy* by Using the Shift Matrix $N$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First let us define variables which we later need, like the velocity $u$, the distribution function $f$, the relaxation rates $\\omega$ and the viscosity $\\rho$." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "u = sp.Matrix(sp.symbols(\"u_:2\")) # velocity\n", + "f = sp.Matrix(sp.symbols(\"f_:9\")) # distribution function\n", + "ω = sp.Matrix(sp.symbols(\"omega_:9\")) # relaxation rates\n", + "ρ = sp.symbols(\"rho\") # viscosity" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We will show the usage of the central moments in an easy setup, the simulation of a force-driven channel flow. First, we define a method with raw moments." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + " <table style=\"border:none; width: 100%\">\n", + " <tr style=\"border:none\">\n", + " <th style=\"border:none\" >Moment</th>\n", + " <th style=\"border:none\" >Eq. Value </th>\n", + " <th style=\"border:none\" >Relaxation Rate</th>\n", + " </tr>\n", + " <tr style=\"border:none\">\n", + " <td style=\"border:none\">$1$</td>\n", + " <td style=\"border:none\">$\\rho$</td>\n", + " <td style=\"border:none\">$\\omega_{0}$</td>\n", + " </tr>\n", + "<tr style=\"border:none\">\n", + " <td style=\"border:none\">$x$</td>\n", + " <td style=\"border:none\">$\\rho u_{0}$</td>\n", + " <td style=\"border:none\">$\\omega_{1}$</td>\n", + " </tr>\n", + "<tr style=\"border:none\">\n", + " <td style=\"border:none\">$y$</td>\n", + " <td style=\"border:none\">$\\rho u_{1}$</td>\n", + " <td style=\"border:none\">$\\omega_{2}$</td>\n", + " </tr>\n", + "<tr style=\"border:none\">\n", + " <td style=\"border:none\">$x^{2}$</td>\n", + " <td style=\"border:none\">$\\rho u_{0}^{2} + \\frac{\\rho}{3}$</td>\n", + " <td style=\"border:none\">$\\omega_{3}$</td>\n", + " </tr>\n", + "<tr style=\"border:none\">\n", + " <td style=\"border:none\">$y^{2}$</td>\n", + " <td style=\"border:none\">$\\rho u_{1}^{2} + \\frac{\\rho}{3}$</td>\n", + " <td style=\"border:none\">$\\omega_{4}$</td>\n", + " </tr>\n", + "<tr style=\"border:none\">\n", + " <td style=\"border:none\">$x y$</td>\n", + " <td style=\"border:none\">$\\rho u_{0} u_{1}$</td>\n", + " <td style=\"border:none\">$\\omega_{5}$</td>\n", + " </tr>\n", + "<tr style=\"border:none\">\n", + " <td style=\"border:none\">$x^{2} y$</td>\n", + " <td style=\"border:none\">$\\rho u_{0}^{2} u_{1} + \\frac{\\rho u_{1}}{3}$</td>\n", + " <td style=\"border:none\">$\\omega_{6}$</td>\n", + " </tr>\n", + "<tr style=\"border:none\">\n", + " <td style=\"border:none\">$x y^{2}$</td>\n", + " <td style=\"border:none\">$\\rho u_{0} u_{1}^{2} + \\frac{\\rho u_{0}}{3}$</td>\n", + " <td style=\"border:none\">$\\omega_{7}$</td>\n", + " </tr>\n", + "<tr style=\"border:none\">\n", + " <td style=\"border:none\">$x^{2} y^{2}$</td>\n", + " <td style=\"border:none\">$\\rho u_{0}^{2} u_{1}^{2} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}$</td>\n", + " <td style=\"border:none\">$\\omega_{8}$</td>\n", + " </tr>\n", + "\n", + " </table>\n", + " " + ], + "text/plain": [ + "<lbmpy.methods.momentbased.MomentBasedLbMethod at 0x7f9b5d1abfd0>" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "stencil = get_stencil(\"D2Q9\", ordering='walberla')\n", + "\n", + "force = [5e-5, 0]\n", + "F = Guo(force)\n", + "\n", + "method = create_lb_method(stencil=stencil, method='mrt_raw', compressible=True, \n", + " equilibrium_order=4, force_model=\"Guo\", force=force)\n", + "\n", + "method" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For the sake of simplicity, we use an SRT approach and set the single relaxation rate to 1.9." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + " <table style=\"border:none; width: 100%\">\n", + " <tr style=\"border:none\">\n", + " <th style=\"border:none\" >Moment</th>\n", + " <th style=\"border:none\" >Eq. Value </th>\n", + " <th style=\"border:none\" >Relaxation Rate</th>\n", + " </tr>\n", + " <tr style=\"border:none\">\n", + " <td style=\"border:none\">$1$</td>\n", + " <td style=\"border:none\">$\\rho$</td>\n", + " <td style=\"border:none\">$1.99$</td>\n", + " </tr>\n", + "<tr style=\"border:none\">\n", + " <td style=\"border:none\">$x$</td>\n", + " <td style=\"border:none\">$\\rho u_{0}$</td>\n", + " <td style=\"border:none\">$1.99$</td>\n", + " </tr>\n", + "<tr style=\"border:none\">\n", + " <td style=\"border:none\">$y$</td>\n", + " <td style=\"border:none\">$\\rho u_{1}$</td>\n", + " <td style=\"border:none\">$1.99$</td>\n", + " </tr>\n", + "<tr style=\"border:none\">\n", + " <td style=\"border:none\">$x^{2}$</td>\n", + " <td style=\"border:none\">$\\rho u_{0}^{2} + \\frac{\\rho}{3}$</td>\n", + " <td style=\"border:none\">$1.99$</td>\n", + " </tr>\n", + "<tr style=\"border:none\">\n", + " <td style=\"border:none\">$y^{2}$</td>\n", + " <td style=\"border:none\">$\\rho u_{1}^{2} + \\frac{\\rho}{3}$</td>\n", + " <td style=\"border:none\">$1.99$</td>\n", + " </tr>\n", + "<tr style=\"border:none\">\n", + " <td style=\"border:none\">$x y$</td>\n", + " <td style=\"border:none\">$\\rho u_{0} u_{1}$</td>\n", + " <td style=\"border:none\">$1.99$</td>\n", + " </tr>\n", + "<tr style=\"border:none\">\n", + " <td style=\"border:none\">$x^{2} y$</td>\n", + " <td style=\"border:none\">$\\rho u_{0}^{2} u_{1} + \\frac{\\rho u_{1}}{3}$</td>\n", + " <td style=\"border:none\">$1.99$</td>\n", + " </tr>\n", + "<tr style=\"border:none\">\n", + " <td style=\"border:none\">$x y^{2}$</td>\n", + " <td style=\"border:none\">$\\rho u_{0} u_{1}^{2} + \\frac{\\rho u_{0}}{3}$</td>\n", + " <td style=\"border:none\">$1.99$</td>\n", + " </tr>\n", + "<tr style=\"border:none\">\n", + " <td style=\"border:none\">$x^{2} y^{2}$</td>\n", + " <td style=\"border:none\">$\\rho u_{0}^{2} u_{1}^{2} + \\frac{\\rho u_{0}^{2}}{3} + \\frac{\\rho u_{1}^{2}}{3} + \\frac{\\rho}{9}$</td>\n", + " <td style=\"border:none\">$1.99$</td>\n", + " </tr>\n", + "\n", + " </table>\n", + " " + ], + "text/plain": [ + "<lbmpy.methods.momentbased.MomentBasedLbMethod at 0x7f9b5d2683d0>" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "def modification_func(moment, eq, rr):\n", + " rr = 1.99\n", + " return moment, eq, rr\n", + "\n", + "method = create_lb_method_from_existing(method, modification_func)\n", + "\n", + "method" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "With the help of the method, we can extract specific vectors/matrices, which we will need to define our collision equation (transformation matrix to transform distribution function $f$ into the moment space, the discrete equilibrium values in the moment space). Further, we use the extracted moments and the defined stencil to get the second transformation matrix with `get_shift_matrix`, to get from the moment space to the central moment space. " + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\- u_{0} & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\- u_{1} & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\\\u_{0}^{2} & - 2 u_{0} & 0 & 1 & 0 & 0 & 0 & 0 & 0\\\\u_{1}^{2} & 0 & - 2 u_{1} & 0 & 1 & 0 & 0 & 0 & 0\\\\u_{0} u_{1} & - u_{1} & - u_{0} & 0 & 0 & 1 & 0 & 0 & 0\\\\- u_{0}^{2} u_{1} & 2 u_{0} u_{1} & u_{0}^{2} & - u_{1} & 0 & - 2 u_{0} & 1 & 0 & 0\\\\- u_{0} u_{1}^{2} & u_{1}^{2} & 2 u_{0} u_{1} & 0 & - u_{0} & - 2 u_{1} & 0 & 1 & 0\\\\u_{0}^{2} u_{1}^{2} & - 2 u_{0} u_{1}^{2} & - 2 u_{0}^{2} u_{1} & u_{1}^{2} & u_{0}^{2} & 4 u_{0} u_{1} & - 2 u_{1} & - 2 u_{0} & 1\\end{matrix}\\right]$" + ], + "text/plain": [ + "⎡ 1 0 0 0 0 0 0 0 0⎤\n", + "⎢ ⎥\n", + "⎢ -u₀ 1 0 0 0 0 0 0 0⎥\n", + "⎢ ⎥\n", + "⎢ -u₁ 0 1 0 0 0 0 0 0⎥\n", + "⎢ ⎥\n", + "⎢ 2 ⎥\n", + "⎢ u₀ -2⋅u₀ 0 1 0 0 0 0 0⎥\n", + "⎢ ⎥\n", + "⎢ 2 ⎥\n", + "⎢ u₁ 0 -2⋅u₁ 0 1 0 0 0 0⎥\n", + "⎢ ⎥\n", + "⎢ u₀⋅u₁ -u₁ -u₀ 0 0 1 0 0 0⎥\n", + "⎢ ⎥\n", + "⎢ 2 2 ⎥\n", + "⎢-u₀ ⋅u₁ 2⋅u₀⋅u₁ u₀ -u₁ 0 -2⋅u₀ 1 0 0⎥\n", + "⎢ ⎥\n", + "⎢ 2 2 ⎥\n", + "⎢-u₀⋅u₁ u₁ 2⋅u₀⋅u₁ 0 -u₀ -2⋅u₁ 0 1 0⎥\n", + "⎢ ⎥\n", + "⎢ 2 2 2 2 2 2 ⎥\n", + "⎣u₀ ⋅u₁ -2⋅u₀⋅u₁ -2⋅u₀ ⋅u₁ u₁ u₀ 4⋅u₀⋅u₁ -2⋅u₁ -2⋅u₀ 1⎦" + ] + }, + "execution_count": 28, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "raw_moments = method.moments\n", + "\n", + "C = moment_matrix(raw_moments, stencil)\n", + "S = sp.diag(*method.relaxation_rates)\n", + "c_eq = sp.Matrix(method.moment_equilibrium_values)\n", + "\n", + "N = get_shift_matrix(raw_moments, stencil)\n", + "\n", + "N" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "With the definition of all relevant variables, vectors and matrices, we can define the collision equation as described above." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "collision_equation = f + C.inv() * N.inv() * S * N * (c_eq - C * f) + sp.Matrix(F(method))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In order to use this collision equation later on in a predefined scenario, we have to convert it into a collision rule, which we will do below, while introducing three subexpressions for the calculation of $\\rho$, $u_0$ and $u_1$." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "collision_rule = [Assignment(lhs, rhs) for lhs, rhs in zip(method.post_collision_pdf_symbols, collision_equation)]\n", + "\n", + "@ps.kernel\n", + "def subexpressions():\n", + " ρ @= sum(f) \n", + " u[0] @= sum(f.T * sp.Matrix(np.array(stencil)[:, 0]))/ρ + force[0]/(2*ρ)\n", + " u[1] @= sum(f.T * sp.Matrix(np.array(stencil)[:, 1]))/ρ + force[1]/(2*ρ)\n", + "\n", + "collision_rule = LbmCollisionRule(method, main_assignments=collision_rule, subexpressions=[subexpressions[0],\n", + " subexpressions[1],\n", + " subexpressions[2]])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "After that, we can use the subexpression elimination from pystencils to reduce the number of operations needed to perform the collision." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "<table style=\"border:none\"><tr><th>Name</th><th>Runtime</th><th>Adds</th><th>Muls</th><th>Divs</th><th>Total</th></tr><tr><td>OriginalTerm</td><td>-</td> <td>899</td> <td>1204</td> <td>3</td> <td>2106</td> </tr><tr><td>sympy_cse</td><td>195.11 ms</td> <td>393</td> <td>246</td> <td>1</td> <td>640</td> </tr></table>" + ], + "text/plain": [ + "<pystencils.simp.simplificationstrategy.SimplificationStrategy.create_simplification_report.<locals>.Report at 0x7f9b78f46610>" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "generic_strategy = ps.simp.SimplificationStrategy()\n", + "generic_strategy.add(ps.simp.sympy_cse)\n", + "generic_strategy.create_simplification_report(collision_rule)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "collision_rule = ps.simp.sympy_cse(collision_rule)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'adds': 393,\n", + " 'muls': 246,\n", + " 'divs': 1,\n", + " 'sqrts': 0,\n", + " 'fast_sqrts': 0,\n", + " 'fast_inv_sqrts': 0,\n", + " 'fast_div': 0}" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "collision_rule.operation_count" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we can pass our method and customized collision rule to the predefined scenario `create_channel` to simulate an easy channel flow." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "<matplotlib.image.AxesImage at 0x7f9b78828280>" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6UAAAFlCAYAAAATVk7bAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUjUlEQVR4nO3dX4yl9X3f8c93ZlgIxAn/tisMqEsb2ohWqu2uXCpXkRuaFrtVl0qphRUlKwtpe4Fbp4nUkNy4l47UxnWk1BI11BvJtUMcR6AKJUHEUdSLUC8Oig009YqAWbTA2gTbgtp4Z769mMfpdL0DZM7M+Z2Zeb0kdM55nnPmfC9+PKu3nuecU90dAAAAGGFp9AAAAADsX6IUAACAYUQpAAAAw4hSAAAAhhGlAAAADCNKAQAAGGZl9ABJcvXVV/fhw4dHjwEAAMAOePTRR7/W3QcvtG8hovTw4cM5efLk6DEAAADYAVX1zGb7XL4LAADAMKIUAACAYUQpAAAAw4hSAAAAhhGlAAAADCNKAQAAGEaUAgAAMMwbRmlV3VtVL1bVlzdsu7KqHqqqr0y3V0zbq6p+tapOVdWfVNU7dnJ4AAAAdrc3c6b0k0luPW/bXUke7u4bkzw8PU6S9yS5cfrveJKPb8+YAAAA7EVvGKXd/YdJXjpv89EkJ6b7J5LctmH7r/e6P0pyeVVds02zAgAAsMds9TOlh7r7zHT/+SSHpvvXJnl2w/NOT9u+T1Udr6qTVXXy7NmzWxwDAACA3WzmLzrq7k7SW3jd3d19pLuPHDx4cNYxAAAA2IVWtvi6F6rqmu4+M12e++K0/bkk12943nXTttf1vx99Kj+x9C+3OAoAAAC71VbPlD6Q5Nh0/1iS+zds/5npW3hvTvKNDZf5AgAAwP/nDc+UVtWnk7w7ydVVdTrJh5N8JMl9VXVHkmeSvG96+oNJ3pvkVJJXk3zgzQxRBy7Kyluvf+MnAgAAsPs8vfmuN4zS7n7/JrtuucBzO8mdb3Ksv7B62YF848hb/7IvAwAAYDd4evNdW/1M6bZaW0le/Sszf+cSAAAAu8xCRGkvJ9+5vEaPAQAAwJwtRJSuHei8ev250WMAAAAwZwsRpVnuLL3lu6OnAAAAYM4WIkqXljo/cNlro8cAAABgzhYiSg+srObwlS+NHgMAAIAd8OTr7FuIKF2p1Vx58SujxwAAAGDOFiRK13LlgVdHjwEAAMCcLUSUHlg6l8OXfG30GAAAAMzZQkTpctZy+bIzpQAAAPvNYkRpdd6y9O3RYwAAADBnCxGlK7WagyvfHD0GAAAAc7YQUbqUzmXld0oBAAD2m4WI0krnolobPQYAAABzthBRupzOW+rc6DEAAACYs4WI0kpyUY2eAgAAgHlbjCityoFSpQAAAPvNQkTpUpKLa2n0GAAAAMzZQkRpkizHmVIAAID9ZiGitFJZijOlAAAA+81CRGmSLPtMKQAAwL6zEFFaiTOlAAAA+5ASBAAAYJiFOFOaJEu+6AgAAGDfcaYUAACAYRbiTGmlsux3SgEAAPYdJQgAAMAwohQAAIBhRCkAAADDiFIAAACGEaUAAAAMI0oBAAAYRpQCAAAwjCgFAABgGFEKAADAMKIUAACAYUQpAAAAw4hSAAAAhhGlAAAADCNKAQAAGEaUAgAAMIwoBQAAYBhRCgAAwDCiFAAAgGFEKQAAAMOIUgAAAIYRpQAAAAwjSgEAABhGlAIAADCMKAUAAGAYUQoAAMAwM0VpVf3bqnq8qr5cVZ+uqkuq6oaqeqSqTlXVb1TVge0aFgAAgL1ly1FaVdcm+TdJjnT3306ynOT2JL+c5KPd/SNJ/jzJHdsxKAAAAHvPrJfvriT5gapaSXJpkjNJfjzJZ6f9J5LcNuN7AAAAsEdtOUq7+7kk/yHJV7Meo99I8miSl7v73PS000muvdDrq+p4VZ2sqpNnv7661TEAAADYxWa5fPeKJEeT3JDkrUkuS3Lrm319d9/d3Ue6+8jBq5a3OgYAAAC72CyX7/6jJH/W3We7+7tJPpfkXUkuny7nTZLrkjw344wAAADsUbNE6VeT3FxVl1ZVJbklyRNJPp/kJ6fnHEty/2wjAgAAsFfN8pnSR7L+hUZfTPKl6W/dneQXkvxcVZ1KclWSe7ZhTgAAAPaglTd+yua6+8NJPnze5qeSvHOWvwsAAMD+MOtPwgAAAMCWiVIAAACGEaUAAAAMI0oBAAAYRpQCAAAwjCgFAABgGFEKAADAMKIUAACAYUQpAAAAw4hSAAAAhhGlAAAADCNKAQAAGEaUAgAAMIwoBQAAYBhRCgAAwDCiFAAAgGFEKQAAAMOIUgAAAIYRpQAAAAwjSgEAABhGlAIAADCMKAUAAGAYUQoAAMAwohQAAIBhRCkAAADDiFIAAACGEaUAAAAMI0oBAAAYRpQCAAAwjCgFAABgGFEKAADAMKIUAACAYUQpAAAAw4hSAAAAhhGlAAAADCNKAQAAGEaUAgAAMIwoBQAAYBhRCgAAwDCiFAAAgGFEKQAAAMOIUgAAAIYRpQAAAAwjSgEAABhGlAIAADCMKAUAAGAYUQoAAMAwohQAAIBhRCkAAADDiFIAAACGmSlKq+ryqvpsVf2vqnqyqv5+VV1ZVQ9V1Vem2yu2a1gAAAD2llnPlH4sye90948m+TtJnkxyV5KHu/vGJA9PjwEAAOD7bDlKq+qHk/xYknuSpLtf6+6XkxxNcmJ62okkt802IgAAAHvVLGdKb0hyNsl/rao/rqpPVNVlSQ5195npOc8nOXShF1fV8ao6WVUnz359dYYxAAAA2K1midKVJO9I8vHufnuSV3Lepbrd3Un6Qi/u7ru7+0h3Hzl41fIMYwAAALBbzRKlp5Oc7u5HpsefzXqkvlBV1yTJdPvibCMCAACwV205Srv7+STPVtXfnDbdkuSJJA8kOTZtO5bk/pkmBAAAYM9amfH1/zrJp6rqQJKnknwg66F7X1XdkeSZJO+b8T0AAADYo2aK0u5+LMmRC+y6ZZa/CwAAwP4w6++UAgAAwJaJUgAAAIYRpQAAAAwjSgEAABhGlAIAADCMKAUAAGAYUQoAAMAwohQAAIBhRCkAAADDiFIAAACGEaUAAAAMI0oBAAAYRpQCAAAwjCgFAABgGFEKAADAMKIUAACAYUQpAAAAw4hSAAAAhhGlAAAADCNKAQAAGEaUAgAAMIwoBQAAYBhRCgAAwDCiFAAAgGFEKQAAAMOIUgAAAIYRpQAAAAwjSgEAABhGlAIAADCMKAUAAGAYUQoAAMAwohQAAIBhRCkAAADDiFIAAACGEaUAAAAMI0oBAAAYRpQCAAAwjCgFAABgGFEKAADAMKIUAACAYUQpAAAAw4hSAAAAhhGlAAAADCNKAQAAGEaUAgAAMIwoBQAAYBhRCgAAwDCiFAAAgGFEKQAAAMOIUgAAAIYRpQAAAAwzc5RW1XJV/XFV/ffp8Q1V9UhVnaqq36iqA7OPCQAAwF60HWdKP5TkyQ2PfznJR7v7R5L8eZI7tuE9AAAA2INmitKqui7JP03yielxJfnxJJ+dnnIiyW2zvAcAAAB716xnSv9Tkn+XZG16fFWSl7v73PT4dJJrL/TCqjpeVSer6uTZr6/OOAYAAAC70ZajtKr+WZIXu/vRrby+u+/u7iPdfeTgVctbHQMAAIBdbGWG174ryT+vqvcmuSTJDyX5WJLLq2plOlt6XZLnZh8TAACAvWjLZ0q7+xe7+7ruPpzk9iS/390/leTzSX5yetqxJPfPPCUAAAB70k78TukvJPm5qjqV9c+Y3rMD7wEAAMAeMMvlu3+hu/8gyR9M959K8s7t+LsAAADsbTtxphQAAADeFFEKAADAMKIUAACAYUQpAAAAw4hSAAAAhhGlAAAADCNKAQAAGEaUAgAAMIwoBQAAYBhRCgAAwDCiFAAAgGFEKQAAAMOIUgAAAIYRpQAAAAwjSgEAABhGlAIAADCMKAUAAGAYUQoAAMAwohQAAIBhRCkAAADDiFIAAACGEaUAAAAMI0oBAAAYRpQCAAAwjCgFAABgGFEKAADAMKIUAACAYUQpAAAAw4hSAAAAhhGlAAAADCNKAQAAGEaUAgAAMIwoBQAAYBhRCgAAwDCiFAAAgGFEKQAAAMOIUgAAAIYRpQAAAAwjSgEAABhGlAIAADCMKAUAAGAYUQoAAMAwohQAAIBhRCkAAADDiFIAAACGEaUAAAAMI0oBAAAYRpQCAAAwjCgFAABgGFEKAADAMFuO0qq6vqo+X1VPVNXjVfWhafuVVfVQVX1lur1i+8YFAABgL5nlTOm5JD/f3TcluTnJnVV1U5K7kjzc3TcmeXh6DAAAAN9ny1Ha3We6+4vT/W8leTLJtUmOJjkxPe1EkttmnBEAAIA9als+U1pVh5O8PckjSQ5195lp1/NJDm3HewAAALD3zBylVfWDSX4ryc929zc37uvuTtKbvO54VZ2sqpNnv7466xgAAADsQjNFaVVdlPUg/VR3f27a/EJVXTPtvybJixd6bXff3d1HuvvIwauWZxkDAACAXWqWb9+tJPckebK7f2XDrgeSHJvuH0ty/9bHAwAAYC9bmeG170ry00m+VFWPTdt+KclHktxXVXckeSbJ+2aaEAAAgD1ry1Ha3f8jSW2y+5at/l0AAAD2j2359l0AAADYClEKAADAMKIUAACAYUQpAAAAw4hSAAAAhhGlAAAADCNKAQAAGEaUAgAAMIwoBQAAYBhRCgAAwDCiFAAAgGFEKQAAAMOIUgAAAIYRpQAAAAwjSgEAABhGlAIAADCMKAUAAGAYUQoAAMAwohQAAIBhRCkAAADDiFIAAACGEaUAAAAMI0oBAAAYRpQCAAAwjCgFAABgGFEKAADAMKIUAACAYVZGD5Aknc5qr40eAwAAgDlzphQAAIBhFuJMaZKspUePAAAAwJw5UwoAAMAwC3GmtJOsxWdKAQAA9puFiNIkWW2X7wIAAOw3CxGlnXamFAAAYB9aiChNklVfdAQAALDvLESUriX5jt8pBQAA2HcWIkq7O6/5TCkAAMC+sxhRmuS7mhQAAGDfWYgoXUvl1V4ePQYAAABztjBR+m1RCgAAsO8sTJS+0gdGjwEAAMCcLUSUnuvlnD33Q6PHAAAAYM4WIkpXu/KttUtGjwEAAMCcLUaUZikvr146egwAAADmbCGi9LW1lTz97atHjwEAAMCcLUSUnuulvPSaM6UAAAD7zYJE6XJe+s5lo8cAAABgzhYiSl87t5ynX7py9BgAAADM2UJE6dpa5f+84ndKAQAA9puFiNKsVta+ddHoKQAAAJizhYjSpdcqlz67EKMAAAAwRztSglV1a5KPJVlO8onu/sjrPn81ufjl3olRAAAAWGDbHqVVtZzk15L8RJLTSb5QVQ909xObvWbpXHLp2bXtHgUAAIAFtxNnSt+Z5FR3P5UkVfWZJEeTbBqly698Nz988swOjAIAAMAi24kovTbJsxsen07y985/UlUdT3I8SS7JpTn39Fd3YBQAAAAW2dKoN+7uu7v7SHcfuSgXjxoDAACAgXbiTOlzSa7f8Pi6adum/sbf/Wt56ORv7sAoAAAAjFZVm+7biTOlX0hyY1XdUFUHktye5IEdeB8AAAB2uW0/U9rd56rqg0l+N+s/CXNvdz++3e8DAADA7rcjv1Pa3Q8meXAn/jYAAAB7x7AvOgIAAABRCgAAwDCiFAAAgGFEKQAAAMOIUgAAAIYRpQAAAAwjSgEAABhGlAIAADCMKAUAAGAYUQoAAMAw1d2jZ0hVnU3yTJKrk3xt8DjsHdYT282aYrtZU2w3a4rtZk2xXf5qdx+80I6FiNLvqaqT3X1k9BzsDdYT282aYrtZU2w3a4rtZk0xDy7fBQAAYBhRCgAAwDCLFqV3jx6APcV6YrtZU2w3a4rtZk2x3awpdtxCfaYUAACA/WXRzpQCAACwjyxElFbVrVX1p1V1qqruGj0Pu1NVPV1VX6qqx6rq5LTtyqp6qKq+Mt1eMXpOFldV3VtVL1bVlzdsu+AaqnW/Oh23/qSq3jFuchbVJmvq31fVc9Ox6rGqeu+Gfb84rak/rap/MmZqFlVVXV9Vn6+qJ6rq8ar60LTdcYoteZ015TjFXA2P0qpaTvJrSd6T5KYk76+qm8ZOxS72D7v7bRu+uvyuJA93941JHp4ew2Y+meTW87Zttobek+TG6b/jST4+pxnZXT6Z719TSfLR6Vj1tu5+MEmmf/tuT/K3ptf85+nfSPiec0l+vrtvSnJzkjundeM4xVZttqYSxynmaHiUJnlnklPd/VR3v5bkM0mODp6JveNokhPT/RNJbhs3Couuu/8wyUvnbd5sDR1N8uu97o+SXF5V18xlUHaNTdbUZo4m+Ux3f6e7/yzJqaz/GwlJku4+091fnO5/K8mTSa6N4xRb9DprajOOU+yIRYjSa5M8u+Hx6bz+/wywmU7ye1X1aFUdn7Yd6u4z0/3nkxwaMxq72GZryLGLWXxwupzy3g0fK7CmeNOq6nCStyd5JI5TbIPz1lTiOMUcLUKUwnb5B939jqxfrnRnVf3Yxp29/lXTvm6aLbOG2CYfT/LXk7wtyZkk/3HoNOw6VfWDSX4ryc929zc37nOcYisusKYcp5irRYjS55Jcv+HxddM2+Evp7uem2xeT/HbWLyd54XuXKk23L46bkF1qszXk2MWWdPcL3b3a3WtJ/kv+36Vv1hRvqKouyno8fKq7Pzdtdpxiyy60phynmLdFiNIvJLmxqm6oqgNZ//D0A4NnYpepqsuq6i3fu5/kHyf5ctbX0rHpaceS3D9mQnaxzdbQA0l+Zvp2y5uTfGPD5XOwqfM+0/cvsn6sStbX1O1VdXFV3ZD1L6f5n/Oej8VVVZXkniRPdvevbNjlOMWWbLamHKeYt5XRA3T3uar6YJLfTbKc5N7ufnzwWOw+h5L89vqxNStJ/lt3/05VfSHJfVV1R5Jnkrxv4IwsuKr6dJJ3J7m6qk4n+XCSj+TCa+jBJO/N+pc8vJrkA3MfmIW3yZp6d1W9LeuXWD6d5F8lSXc/XlX3JXki69+IeWd3rw4Ym8X1riQ/neRLVfXYtO2X4jjF1m22pt7vOMU81fpHDwAAAGD+FuHyXQAAAPYpUQoAAMAwohQAAIBhRCkAAADDiFIAAACGEaUAAAAMI0oBAAAYRpQCAAAwzP8FTuZGRHiwWF0AAAAASUVORK5CYII=\n", + "text/plain": [ + "<Figure size 1152x432 with 1 Axes>" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "ch_with_cm = create_channel((300,100), method=method, force=force[0], force_model=\"Guo\", initial_velocity=(0.5, 0), collision_rule=collision_rule,\n", + " optimization={\"target\": \"cpu\", \"openmp\": 4})\n", + "\n", + "ch_with_cm.run(10000)\n", + "plt.vector_field_magnitude(ch_with_cm.velocity[:, :])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Just for verification, we simulate the same setup, but with a standard SRT approach, without central moments." + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "<matplotlib.image.AxesImage at 0x7f9b73f22d00>" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6UAAAFlCAYAAAATVk7bAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUjUlEQVR4nO3dX4yl9X3f8c93ZlgIxAn/tisMqEsb2ohWqu2uXCpXkRuaFrtVl0qphRUlKwtpe4Fbp4nUkNy4l47UxnWk1BI11BvJtUMcR6AKJUHEUdSLUC8Oig009YqAWbTA2gTbgtp4Z769mMfpdL0DZM7M+Z2Zeb0kdM55nnPmfC9+PKu3nuecU90dAAAAGGFp9AAAAADsX6IUAACAYUQpAAAAw4hSAAAAhhGlAAAADCNKAQAAGGZl9ABJcvXVV/fhw4dHjwEAAMAOePTRR7/W3QcvtG8hovTw4cM5efLk6DEAAADYAVX1zGb7XL4LAADAMKIUAACAYUQpAAAAw4hSAAAAhhGlAAAADCNKAQAAGEaUAgAAMMwbRmlV3VtVL1bVlzdsu7KqHqqqr0y3V0zbq6p+tapOVdWfVNU7dnJ4AAAAdrc3c6b0k0luPW/bXUke7u4bkzw8PU6S9yS5cfrveJKPb8+YAAAA7EVvGKXd/YdJXjpv89EkJ6b7J5LctmH7r/e6P0pyeVVds02zAgAAsMds9TOlh7r7zHT/+SSHpvvXJnl2w/NOT9u+T1Udr6qTVXXy7NmzWxwDAACA3WzmLzrq7k7SW3jd3d19pLuPHDx4cNYxAAAA2IVWtvi6F6rqmu4+M12e++K0/bkk12943nXTttf1vx99Kj+x9C+3OAoAAAC71VbPlD6Q5Nh0/1iS+zds/5npW3hvTvKNDZf5AgAAwP/nDc+UVtWnk7w7ydVVdTrJh5N8JMl9VXVHkmeSvG96+oNJ3pvkVJJXk3zgzQxRBy7Kyluvf+MnAgAAsPs8vfmuN4zS7n7/JrtuucBzO8mdb3Ksv7B62YF848hb/7IvAwAAYDd4evNdW/1M6bZaW0le/Sszf+cSAAAAu8xCRGkvJ9+5vEaPAQAAwJwtRJSuHei8ev250WMAAAAwZwsRpVnuLL3lu6OnAAAAYM4WIkqXljo/cNlro8cAAABgzhYiSg+srObwlS+NHgMAAIAd8OTr7FuIKF2p1Vx58SujxwAAAGDOFiRK13LlgVdHjwEAAMCcLUSUHlg6l8OXfG30GAAAAMzZQkTpctZy+bIzpQAAAPvNYkRpdd6y9O3RYwAAADBnCxGlK7WagyvfHD0GAAAAc7YQUbqUzmXld0oBAAD2m4WI0krnolobPQYAAABzthBRupzOW+rc6DEAAACYs4WI0kpyUY2eAgAAgHlbjCityoFSpQAAAPvNQkTpUpKLa2n0GAAAAMzZQkRpkizHmVIAAID9ZiGitFJZijOlAAAA+81CRGmSLPtMKQAAwL6zEFFaiTOlAAAA+5ASBAAAYJiFOFOaJEu+6AgAAGDfcaYUAACAYRbiTGmlsux3SgEAAPYdJQgAAMAwohQAAIBhRCkAAADDiFIAAACGEaUAAAAMI0oBAAAYRpQCAAAwjCgFAABgGFEKAADAMKIUAACAYUQpAAAAw4hSAAAAhhGlAAAADCNKAQAAGEaUAgAAMIwoBQAAYBhRCgAAwDCiFAAAgGFEKQAAAMOIUgAAAIYRpQAAAAwjSgEAABhGlAIAADCMKAUAAGAYUQoAAMAwM0VpVf3bqnq8qr5cVZ+uqkuq6oaqeqSqTlXVb1TVge0aFgAAgL1ly1FaVdcm+TdJjnT3306ynOT2JL+c5KPd/SNJ/jzJHdsxKAAAAHvPrJfvriT5gapaSXJpkjNJfjzJZ6f9J5LcNuN7AAAAsEdtOUq7+7kk/yHJV7Meo99I8miSl7v73PS000muvdDrq+p4VZ2sqpNnv7661TEAAADYxWa5fPeKJEeT3JDkrUkuS3Lrm319d9/d3Ue6+8jBq5a3OgYAAAC72CyX7/6jJH/W3We7+7tJPpfkXUkuny7nTZLrkjw344wAAADsUbNE6VeT3FxVl1ZVJbklyRNJPp/kJ6fnHEty/2wjAgAAsFfN8pnSR7L+hUZfTPKl6W/dneQXkvxcVZ1KclWSe7ZhTgAAAPaglTd+yua6+8NJPnze5qeSvHOWvwsAAMD+MOtPwgAAAMCWiVIAAACGEaUAAAAMI0oBAAAYRpQCAAAwjCgFAABgGFEKAADAMKIUAACAYUQpAAAAw4hSAAAAhhGlAAAADCNKAQAAGEaUAgAAMIwoBQAAYBhRCgAAwDCiFAAAgGFEKQAAAMOIUgAAAIYRpQAAAAwjSgEAABhGlAIAADCMKAUAAGAYUQoAAMAwohQAAIBhRCkAAADDiFIAAACGEaUAAAAMI0oBAAAYRpQCAAAwjCgFAABgGFEKAADAMKIUAACAYUQpAAAAw4hSAAAAhhGlAAAADCNKAQAAGEaUAgAAMIwoBQAAYBhRCgAAwDCiFAAAgGFEKQAAAMOIUgAAAIYRpQAAAAwjSgEAABhGlAIAADCMKAUAAGAYUQoAAMAwohQAAIBhRCkAAADDiFIAAACGmSlKq+ryqvpsVf2vqnqyqv5+VV1ZVQ9V1Vem2yu2a1gAAAD2llnPlH4sye90948m+TtJnkxyV5KHu/vGJA9PjwEAAOD7bDlKq+qHk/xYknuSpLtf6+6XkxxNcmJ62okkt802IgAAAHvVLGdKb0hyNsl/rao/rqpPVNVlSQ5195npOc8nOXShF1fV8ao6WVUnz359dYYxAAAA2K1midKVJO9I8vHufnuSV3Lepbrd3Un6Qi/u7ru7+0h3Hzl41fIMYwAAALBbzRKlp5Oc7u5HpsefzXqkvlBV1yTJdPvibCMCAACwV205Srv7+STPVtXfnDbdkuSJJA8kOTZtO5bk/pkmBAAAYM9amfH1/zrJp6rqQJKnknwg66F7X1XdkeSZJO+b8T0AAADYo2aK0u5+LMmRC+y6ZZa/CwAAwP4w6++UAgAAwJaJUgAAAIYRpQAAAAwjSgEAABhGlAIAADCMKAUAAGAYUQoAAMAwohQAAIBhRCkAAADDiFIAAACGEaUAAAAMI0oBAAAYRpQCAAAwjCgFAABgGFEKAADAMKIUAACAYUQpAAAAw4hSAAAAhhGlAAAADCNKAQAAGEaUAgAAMIwoBQAAYBhRCgAAwDCiFAAAgGFEKQAAAMOIUgAAAIYRpQAAAAwjSgEAABhGlAIAADCMKAUAAGAYUQoAAMAwohQAAIBhRCkAAADDiFIAAACGEaUAAAAMI0oBAAAYRpQCAAAwjCgFAABgGFEKAADAMKIUAACAYUQpAAAAw4hSAAAAhhGlAAAADCNKAQAAGEaUAgAAMIwoBQAAYBhRCgAAwDCiFAAAgGFEKQAAAMOIUgAAAIYRpQAAAAwzc5RW1XJV/XFV/ffp8Q1V9UhVnaqq36iqA7OPCQAAwF60HWdKP5TkyQ2PfznJR7v7R5L8eZI7tuE9AAAA2INmitKqui7JP03yielxJfnxJJ+dnnIiyW2zvAcAAAB716xnSv9Tkn+XZG16fFWSl7v73PT4dJJrL/TCqjpeVSer6uTZr6/OOAYAAAC70ZajtKr+WZIXu/vRrby+u+/u7iPdfeTgVctbHQMAAIBdbGWG174ryT+vqvcmuSTJDyX5WJLLq2plOlt6XZLnZh8TAACAvWjLZ0q7+xe7+7ruPpzk9iS/390/leTzSX5yetqxJPfPPCUAAAB70k78TukvJPm5qjqV9c+Y3rMD7wEAAMAeMMvlu3+hu/8gyR9M959K8s7t+LsAAADsbTtxphQAAADeFFEKAADAMKIUAACAYUQpAAAAw4hSAAAAhhGlAAAADCNKAQAAGEaUAgAAMIwoBQAAYBhRCgAAwDCiFAAAgGFEKQAAAMOIUgAAAIYRpQAAAAwjSgEAABhGlAIAADCMKAUAAGAYUQoAAMAwohQAAIBhRCkAAADDiFIAAACGEaUAAAAMI0oBAAAYRpQCAAAwjCgFAABgGFEKAADAMKIUAACAYUQpAAAAw4hSAAAAhhGlAAAADCNKAQAAGEaUAgAAMIwoBQAAYBhRCgAAwDCiFAAAgGFEKQAAAMOIUgAAAIYRpQAAAAwjSgEAABhGlAIAADCMKAUAAGAYUQoAAMAwohQAAIBhRCkAAADDiFIAAACGEaUAAAAMI0oBAAAYRpQCAAAwjCgFAABgGFEKAADAMFuO0qq6vqo+X1VPVNXjVfWhafuVVfVQVX1lur1i+8YFAABgL5nlTOm5JD/f3TcluTnJnVV1U5K7kjzc3TcmeXh6DAAAAN9ny1Ha3We6+4vT/W8leTLJtUmOJjkxPe1EkttmnBEAAIA9als+U1pVh5O8PckjSQ5195lp1/NJDm3HewAAALD3zBylVfWDSX4ryc929zc37uvuTtKbvO54VZ2sqpNnv7466xgAAADsQjNFaVVdlPUg/VR3f27a/EJVXTPtvybJixd6bXff3d1HuvvIwauWZxkDAACAXWqWb9+tJPckebK7f2XDrgeSHJvuH0ty/9bHAwAAYC9bmeG170ry00m+VFWPTdt+KclHktxXVXckeSbJ+2aaEAAAgD1ry1Ha3f8jSW2y+5at/l0AAAD2j2359l0AAADYClEKAADAMKIUAACAYUQpAAAAw4hSAAAAhhGlAAAADCNKAQAAGEaUAgAAMIwoBQAAYBhRCgAAwDCiFAAAgGFEKQAAAMOIUgAAAIYRpQAAAAwjSgEAABhGlAIAADCMKAUAAGAYUQoAAMAwohQAAIBhRCkAAADDiFIAAACGEaUAAAAMI0oBAAAYRpQCAAAwjCgFAABgGFEKAADAMKIUAACAYVZGD5Aknc5qr40eAwAAgDlzphQAAIBhFuJMaZKspUePAAAAwJw5UwoAAMAwC3GmtJOsxWdKAQAA9puFiNIkWW2X7wIAAOw3CxGlnXamFAAAYB9aiChNklVfdAQAALDvLESUriX5jt8pBQAA2HcWIkq7O6/5TCkAAMC+sxhRmuS7mhQAAGDfWYgoXUvl1V4ePQYAAABztjBR+m1RCgAAsO8sTJS+0gdGjwEAAMCcLUSUnuvlnD33Q6PHAAAAYM4WIkpXu/KttUtGjwEAAMCcLUaUZikvr146egwAAADmbCGi9LW1lTz97atHjwEAAMCcLUSUnuulvPSaM6UAAAD7zYJE6XJe+s5lo8cAAABgzhYiSl87t5ynX7py9BgAAADM2UJE6dpa5f+84ndKAQAA9puFiNKsVta+ddHoKQAAAJizhYjSpdcqlz67EKMAAAAwRztSglV1a5KPJVlO8onu/sjrPn81ufjl3olRAAAAWGDbHqVVtZzk15L8RJLTSb5QVQ909xObvWbpXHLp2bXtHgUAAIAFtxNnSt+Z5FR3P5UkVfWZJEeTbBqly698Nz988swOjAIAAMAi24kovTbJsxsen07y985/UlUdT3I8SS7JpTn39Fd3YBQAAAAW2dKoN+7uu7v7SHcfuSgXjxoDAACAgXbiTOlzSa7f8Pi6adum/sbf/Wt56ORv7sAoAAAAjFZVm+7biTOlX0hyY1XdUFUHktye5IEdeB8AAAB2uW0/U9rd56rqg0l+N+s/CXNvdz++3e8DAADA7rcjv1Pa3Q8meXAn/jYAAAB7x7AvOgIAAABRCgAAwDCiFAAAgGFEKQAAAMOIUgAAAIYRpQAAAAwjSgEAABhGlAIAADCMKAUAAGAYUQoAAMAw1d2jZ0hVnU3yTJKrk3xt8DjsHdYT282aYrtZU2w3a4rtZk2xXf5qdx+80I6FiNLvqaqT3X1k9BzsDdYT282aYrtZU2w3a4rtZk0xDy7fBQAAYBhRCgAAwDCLFqV3jx6APcV6YrtZU2w3a4rtZk2x3awpdtxCfaYUAACA/WXRzpQCAACwjyxElFbVrVX1p1V1qqruGj0Pu1NVPV1VX6qqx6rq5LTtyqp6qKq+Mt1eMXpOFldV3VtVL1bVlzdsu+AaqnW/Oh23/qSq3jFuchbVJmvq31fVc9Ox6rGqeu+Gfb84rak/rap/MmZqFlVVXV9Vn6+qJ6rq8ar60LTdcYoteZ015TjFXA2P0qpaTvJrSd6T5KYk76+qm8ZOxS72D7v7bRu+uvyuJA93941JHp4ew2Y+meTW87Zttobek+TG6b/jST4+pxnZXT6Z719TSfLR6Vj1tu5+MEmmf/tuT/K3ptf85+nfSPiec0l+vrtvSnJzkjundeM4xVZttqYSxynmaHiUJnlnklPd/VR3v5bkM0mODp6JveNokhPT/RNJbhs3Couuu/8wyUvnbd5sDR1N8uu97o+SXF5V18xlUHaNTdbUZo4m+Ux3f6e7/yzJqaz/GwlJku4+091fnO5/K8mTSa6N4xRb9DprajOOU+yIRYjSa5M8u+Hx6bz+/wywmU7ye1X1aFUdn7Yd6u4z0/3nkxwaMxq72GZryLGLWXxwupzy3g0fK7CmeNOq6nCStyd5JI5TbIPz1lTiOMUcLUKUwnb5B939jqxfrnRnVf3Yxp29/lXTvm6aLbOG2CYfT/LXk7wtyZkk/3HoNOw6VfWDSX4ryc929zc37nOcYisusKYcp5irRYjS55Jcv+HxddM2+Evp7uem2xeT/HbWLyd54XuXKk23L46bkF1qszXk2MWWdPcL3b3a3WtJ/kv+36Vv1hRvqKouyno8fKq7Pzdtdpxiyy60phynmLdFiNIvJLmxqm6oqgNZ//D0A4NnYpepqsuq6i3fu5/kHyf5ctbX0rHpaceS3D9mQnaxzdbQA0l+Zvp2y5uTfGPD5XOwqfM+0/cvsn6sStbX1O1VdXFV3ZD1L6f5n/Oej8VVVZXkniRPdvevbNjlOMWWbLamHKeYt5XRA3T3uar6YJLfTbKc5N7ufnzwWOw+h5L89vqxNStJ/lt3/05VfSHJfVV1R5Jnkrxv4IwsuKr6dJJ3J7m6qk4n+XCSj+TCa+jBJO/N+pc8vJrkA3MfmIW3yZp6d1W9LeuXWD6d5F8lSXc/XlX3JXki69+IeWd3rw4Ym8X1riQ/neRLVfXYtO2X4jjF1m22pt7vOMU81fpHDwAAAGD+FuHyXQAAAPYpUQoAAMAwohQAAIBhRCkAAADDiFIAAACGEaUAAAAMI0oBAAAYRpQCAAAwzP8FTuZGRHiwWF0AAAAASUVORK5CYII=\n", + "text/plain": [ + "<Figure size 1152x432 with 1 Axes>" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "ch_normal_srt = create_channel(domain_size=(300,100), force=force[0], initial_velocity=(0.5,0),\n", + " relaxation_rate=1.99, optimization={'target': 'cpu'})\n", + "ch_normal_srt.run(10000)\n", + "plt.vector_field_magnitude(ch_normal_srt.velocity[:, :])" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "<matplotlib.legend.Legend at 0x7f9b73ed9370>" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxIAAAJXCAYAAAAHJmKgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAACOjklEQVR4nOzdd3hUZf7+8feTnkkyoSWhSAldQUAE7CCKBQuigiCuggV/ttVdd22rrq7u6q7rrro21rUgVlRA8WtHAVGpCioIKNKbhJYJTHqe3x+TxCSkTTKTMzO5X9eVCzLnmTN3Zk4m5zPnKcZai4iIiIiIiD+inA4gIiIiIiLhR4WEiIiIiIj4TYWEiIiIiIj4TYWEiIiIiIj4TYWEiIiIiIj4TYWEiIiIiIj4LcbpAMHSpk0b26VLF8cev6SkBICoKNVq4qNjQqrSMSHV0XEhVemYkKqa8pj4+uuvd1tr06rbFrGFRJcuXVi2bJljj5+TkwNASkqKYxkktOiYkKp0TEh1dFxIVTompKqmPCaMMZtq2qbSVkRERERE/KZCQkRERERE/KZCQkRERERE/KZCQkRERERE/KZCQkRERERE/KZCQkRERERE/KZCQkRERERE/KZCQkRERERE/KZCQkRERERE/KZCQkRERERE/KZCQkRERERE/KZCQkRERERE/KZCQkRERERE/KZCQkRERERE/KZCQkRERERE/KZCQkRERERE/KZCQkRERERE/KZCQkRERERE/KZCQkRERERE/BYShYQx5kxjzFpjzDpjzO3VbJ9kjMkyxqwo/brKiZwiIiIiIuIT43QAY0w08CRwGrAVWGqMmW2t/aFK0+nW2huaPKCIiIiIiBwiFK5IDAHWWWvXW2sLgNeB8xzOJCIiIiIitXD8igTQAdhS4futwDHVtLvQGDMU+BH4vbV2SzVtypWUlJCTkxO4lH5y8rFDVmEhZvduzC+/YHbtgl1ZHNyRg+fccXgTWlGwYBkF788j+6obybUJFH08H/PlIopKoigqgqJiQ0lx6b/WYC2UWMi/9npsQiJm/gKivv6GgptuwlqI+fBDzCrfhS1bU6aoaApuugmAmI8+xGzdRuGVV/q+nzGDqM2bav2RrDu1UntTWEjh+PEAxE57EbNnz69trS+FMebX29q1r9TetmxF0bnnAhD31FOQn1fr45d071GpfcnhvSkafgoUFxP3n8dqvS9AyYABldoXn3AixUOGYLKziX3+uTrvX7V98RlnUnzEEURt307M9Nfrvn+V9kUXXEhJ585E/fQjMf/3f3Xev2r7wksvw7ZpQ/Ty5UTPm1vn/au2L7j2WkhIJPqLL4heuqTO+1dq//XXlY6lqB+qXlStIiqanGuuASBp3ttBPfaqvb+OvZA99qKWLAYqv1dU1dhjr6ne96q9v449v4+9in8/gvu+t9T3OJQ93q/3M8Z3e+H112FcicR+8Tkx33xN4c2/JyoK4j56n7gfvic6BmKiLTHRluhoiImBhERwuSAhJZroO2/G5YKkhZ+Rmr0Vc+Wlvv3v3IlNSIDU1MoPLNUKlfPMUCgk6uNd4DVrbb4x5v8BLwKnVG1kjLkauBqgY8eOTZtQDlGwYClr/jOfFd/Hsm2vi115bnaRTglRQDxwGCnkELc8noT2cSSvd5O8MJmYvpDQOgrX7mjiPdG+N6MYiIspISbe98YUFeX7MgYKOhYTlVxCzJYiYvbmU9irBIDY9YVEFRaUvyFW+8YUFUVB3xKMgZhNBZikPAqP9N0/5ocCohLza/0ZrbtyewoKKCr9PrZ3AWb3r/cvLi4GIDo6uvy2kvb5ldrb1oUVvs/H5NX++CW9KrcvObyQ4iNLoKiYuF613xeguEr74j5FFB9ZgtlfTGx97l+lfdERhZQcWYJpWVSv+1dtX9inEJtZgoktIvanuu9ftX1hnyJseglRBUXE7KjP/Su3L+hTDK4SovcUEO2p+/6V2h/Mo6DsWNhUQFRxHfePjmL/EYUAuHcG99irjo690D322JcLVH6vqKqxx17F9jr2Qv/Yq/j3IxDHXklaCVF5hcTsKKCgTwkklh5L2b8WcbbCJ3Bl/y+xhvyuJdiEEqI2FGG2F1HQroSSEoNxF2NdUFxkKSw25OVHUVRsfP8viCKvIBqvTSS7XRx5eZA/J5Wc3QnkfZgIQNTS3aTs3Ux69G4ykg/S87CDDDg2hsOGdqFk4FHYTp1UYIQgY22Nn9U2TQBjjgPutdaeUfr9HQDW2gdraB8N7LXWpta230GDBtlly5YFOm69lVWKKSkpjmVoSkVF8O1bP7HsluksO2oyWwoySNy3nf4rX2Fgby+dusXStlMcbTJTiOmQARmlX23aQHKyryqIcM3tmJC66ZiQ6ui4kKoi8pgoKACvF1q08H3/3nsc+G49v2zwsnNLIatXW5ZtyWBTSUdiKaRf4k8MOimRIS9cR/v2jiYPCU15TBhjvrbWDqpuWyhckVgK9DDGZALbgPHAhIoNjDHtrLU7Sr8dBaxu2ohSrV27yHrkZZ5ZfyqfZvXnuCPbMbhNMfdcvpkOozMwti2YP+oTBBEREaksLs73Vebss0k+G5KBbsAJwFX5+fD99+Qv/Ibv53hZVtiD226DAzmWqw8+whl/Poaok05w6AcQCIFCwlpbZIy5AfgIiAaet9auMsbcByyz1s4GbjTGjAKKgL3AJMcCCyxbxjf3zuaJD7qyv6QLV5/zBXfM6U9UVDJwz6/tTORfZRAREZEgiY+HQYOIHzSIQb+FQcA1QNbSTfzv1AIevrY7510Nky7w4E6xvvEV0qQc79oULOraFHh22de8dclMXvjxBHrEbuT6i3bT8+5x0KuX09HCQiQeE9I4OiakOjoupCodE9UoLKSwyPD2/8Uw9a51dF0/h99P3EvXJ//gK0AinLo2SVjJn/42/++SHHokpTL9H5tIueY34HY7HUtERESao9hYYmNh7FgY2/Mgy+/cxrXPncgfl/yJ0+b+CVq3djphs6C+J1I7a9n95/9w/vg4xnT/ljt/mkTKrdeqiBAREZHQ0L8/R/3f/bw9NZsXVg1hyuGPwdq1TqdqFlRISK3WrCxi7OMn8eDp8zhn+f2Qnu50JBEREZFDJE68iJfnd2TrwRbc3P9Tij+d53SkiKdCQqq3bx+fzsrmtzfH8urSnvT/4O+QmOh0KhEREZEaRZ14PH9deT79W21h/Gm7OfD0S05HimgqJORQxcU82+8//O/ab3h7lqVd96Rmsc6DiIiIRIDMTCb+cBs3DFzI6OvasfWRN51OFLF0diiHuP+BaH46ejyvvhFLUrLWgBAREZEw06IFwxb+nadvWsulM89n3TqnA0UmFRJSydI3NvDzz/CPt3sRNfREp+OIiIiINExsLD0evZ7np8Xw+xuLsT+vdzpRxFEhIeWKZr/PHeN+5p+nf+J0FBEREZGAyMyEk9Y+y0snPwclJU7HiSgqJMTn4EEem/gNE9rNI23MMKfTiIiIiATM75/szotpf2T3Xp36BpKeTQFg0+8f5bP9R3H59DMhLs7pOCIiIiIBE3vmqfztqZbcdhu6KhFAKiQEu3wFf/hfb/510WLMSRoXISIiIpHn2GMh8fslzDvhTrDW6TgRQYVEc1dczIyxr9PXtZ7eU37vdBoRERGRoPnbmOXct+g08qa94XSUiKBCopnL/tezPPnzGdz+VCdo2dLpOCIiIiJBk/qHq7iu60f8/brNsG+f03HCngqJ5mzrVu68y/DnQR+QcNlFTqcRERERCa7oaC58czzfe7ux5v894nSasKdCohlbdO2L5JbEM3z6NWC08JyIiIhEPjPwKP49eTV/ePMY7IIvnI4T1lRINFOFhXBXzq3849WO0LWr03FEREREmkznf9/EKS2WM3X8B1BQ4HScsKVCopl66SW4aEIsbS46xekoIiIiIk0rOZmbnu/PtO0jKHjlTafThC0VEs3RqlXMvHkBE/p+53QSEREREUfEnHc2Z7f4ik+eWOt0lLClQqIZ2vqjlxYxB0juluF0FBERERFnREUxfnIKr2eP1LoSDRTjdABpetPXD2b8VEB1hIiIiDRjhz10I3tXwUEvJCU5nSb86IpEc/PLL3z0XhGnn+50EBERERHnnTfKMvuZHU7HCEsqJJqZNX/4H12/eok4NEOBiIiIyIU5U3nr5q9gwwano4QdFRLNSVERr81KYMJxGyEuzuk0IiIiIo5rPW4EZsAA9tDa6ShhR4VEM2LnfMoX3qM48bdHOR1FREREJDR07MjY27vx1sdup5OEHRUSzciyxxcyKO57os4e6XQUERERkZBx7vADzH5iM/z4o9NRwooKieYiN5dXP2nDhLP2Q3y802lEREREQobL5NJ61Xy2PDbT6ShhRYVEM1H8zv/xfWFv+t0w1OkoIiIiIqElLY2Lj/6J118tgZISp9OEDRUSzcS8x79nePIyzMnDnI4iIiIiEnJGXN+LT/YPgoULnY4SNlRINAf79vHqwkwuHlsE0dFOpxEREREJObEXjqJH9AZ+eOIzp6OEDRUSzUDenC/YbtvS9boznY4iIiIiEppSUpgwbBuvveOCwkKn04QFFRLNwAcx53LWfcfC0Uc7HUVEREQkZB3320F8lTsA+8kcp6OEBRUSzcD06XDR1S3BGKejiIiIiISsqLPOZEjCdyx5fLHTUcKCCokI53loCrnzFpHhznU6ioiIiEhoi4tjwtkeXp2TDl6v02lCngqJCDdrZQ8uaL8YEhOdjiIiIiIS8o68cTg/uI+laPsup6OEvBinA0hwzcw+lZfmnep0DBEREZHwMHQop/wR5m6A07o7HSa06YpEBMvZvI/oKIvb7XQSERERkfAxZgy8/XoeFBU5HSWkqZCIYN9M+g+Dv3zU6RgiIiIiYaX7mv9j/fNzYcUKp6OENHVtimBLv41jSB+P0zFEREREwooZdDRxhx8gv0UC8U6HCWG6IhGp9uzh671dOPqUVKeTiIiIiISXdu3oP6YH3+3r6HSSkKZCIlItW8Z+WtBiaD+nk4iIiIiEncG9PCx99SenY4Q0FRIRKmveKlqzR6tZi4iIiDTA4G/+y5JHv4IDB5yOErJUSESoZZ95GJy+GVLVtUlERETEX21POYJfSIdvvnE6SshSIRGJrGXpDy4GH6Upy0REREQaZPBgkjlAzoIVTicJWSokItG2baw40J0Bp6U5nUREREQkPKWnMzB1Pd98us/pJCFLhUQEskuWkksirhOOcjqKiIiISNgafGQeS7+NczpGyFIhEYG2th1Ex+M7woABTkcRERERCVuDTk3l671dYM8ep6OEJBUSEWjpzo4MntQXEhKcjiIiIiIStloM7cd+WsCyZU5HCUkqJCKNtSyZtobBXVU5i4iIiDTK0UfTmj1kzVvldJKQpEIi0qxbx6p3fqLPunecTiIiIiIS3lJTGZy+iWVzc5xOEpJinA4ggVXSqQvFx6cTe0GB01FEREREwt7gv49hzg/tGel0kBCkKxIR5qeNsfQYlAppmvpVREREpLEGjOvFip9TnI4RklRIRJilf/2QwYkrnY4hIiIiEhFcxTnk/rABu3iJ01FCjgqJSFJUxNLXfmbwjtlOJxERERGJDLGxHPbjZ2ydqUKiKo2RiCSrVvFTcSY9TtvvdBIRERGRyJCQwOB/jWdp5yQ6Op0lxOiKRAQpXLiMaIqJOmaw01FEREREIsaQ4UksXep0itCjQiKCrPxkB33j10H37k5HEREREYkYfYq/Y+XUZbB+vdNRQooKiQiydIllcC8PGON0FBEREZGIERtjKd65i5IvFzodJaSokIgUubks3daOwUMTnU4iIiIiEln69KFH9AZ+mrPJ6SQhRYVEpFixgq22A4ed0tPpJCIiIiKRJSaGwV33sHRhodNJQooKiQjh/XI5ieRihmigtYiIiEigDT4uhqXr20BRkdNRQoYKiQixYrllgHs9dOjgdBQRERGRiNPjtC78VJwJq1Y5HSVkqJCIEEuHXM/gl25yOoaIiIhIRIo6ZjBRlFC4cJnTUUKGCokIsWQJDD4+1ukYIiIiIpGpe3f6xq9j5Sc7nE4SMlRIRIIvvmDPR8tok/2z00lEREREIpMxDOmVrYXpKlAhEQH278ilRdFuaNnS6SgiIiIiEWvwWWksZTBY63SUkBDjdABpvGUtT2PQXUArp5OIiIiIRK7DHriOrSsArf0L6IpERFi6xDJYs76KiIiIBJUxkJgI3gMlTkcJCSokwt3evay4dxYDN8xwOomIiIhIxOv3zVS+n/Cg0zFCggqJcLd+PQcK40lpEe10EhEREZGI1+2oFNZnHOd0jJCgQiLcbdjg+zcz09kcIiIiIs1A5h8uZGPXU5yOERJUSIQ5z+ptuPFAly5ORxERERGJeJmZsGFtAZRonIQKiTC3YeVBMhN2QGqq01FEREREIl67OS+x48WPYONGp6M4ToVEmNuwrpgu6V6nY4iIiIg0C1GdDqOEqF+7lzdjKiTC3IatsWR20qU1ERERkSaRmUkUJRT/vNHpJI5TIRHOSkrYsNdNZq84p5OIiIiINA+HHUYHs51t3+1xOonjVEiEs5072VzcgU79WjidRERERKR5iIkhs8V+NqzJczqJ41RIhLPiYgo7dCFuyACnk4iIiIg0G5nt89m40Tgdw3EqJMKYPawjZsAAOPZYp6OIiIiINBuZXQ0bdiY6HcNxKiTC2O5t+bRpY52OISIiItKsZPZxseFgGuTmOh3FUSokwtiG6/5Jl9mPOx1DREREpFlp3actu2nT7NeSiHE6gDTchh6nkZm83+kYIiIiIs2KGTIYegDuJKejOEqFRBjbkHYMx5/ndAoRERGRZqZnT+L6QH4biHc6i4PUtSlcFRezYcV+MtvnO51EREREpNnplJbL5iU7nY7hKBUS4WrLFrZPX0D7z152OomIiIhIs5P59iNsvOO/TsdwlLo2hasNG7AYort1cTqJiIiISLOTeeUpbChs5XQMR+mKRJgq+XkDBguZmU5HEREREWl2Mi8+lg2xPZ2O4SgVEmFq+/d7aMdO6NjR6SgiIiIizU6XlD1s+GpHs15LQoVEmNqwOo/MFvsgNtbpKCIiIiLNjvvrueR8/g2sXet0FMeokAhTGzZAZrs8p2OIiIiINE9l3cs3bHA2h4NUSISpjTsTyOxqnI4hIiIi0jx16UISBzm4ZovTSRyjQiIc5eWx4UAbMo9IdDqJiIiISPPUqhVdYrez8fscp5M4RoVEONq0iSzSSOub4XQSERERkebJGDLTD7BhXbHTSRyjdSTCUevW2L4JmBOb74ErIiIi4rTMjsX8tLH5TnyjKxJhqDC1DTFdO0PXrk5HEREREWm2uvSMY8PuFLDW6SiOUCERhjZ/+hOdkvc4HUNERESkWevSz82movaQleV0FEeokAhDGx54jcyPn3E6hoiIiEizltCzE3kkNNspYDVGIgxtPPVKMk/f53QMERERkebt+ONhSAn0jnM6iSNUSIShDQUdGDS2g9MxRERERJq31q1p2QP2lUBLp7M4QF2bwo3Xy4ZP15MZ03wXPxEREREJFZklP7Ph9cVOx3CEColws24d2YtXk7p6kdNJRERERJq9zIWvsmHKR07HcIS6NoWbDRuAGOjSxekkIiIiIs1el/uv5JufU52O4QhdkQgz3rVbSCQXMjOdjiIiIiLS7GWe0J4Nu5KcjuEIXZEIMxu/89AltgRat3Y6ioiIiEiz19G7li1zLGS1hrQ0p+M0KRUSYWbjuiIy03LAGKejiIiIiDR7MVs3UvxjEazJanaFhLo2hZkNW6LJ7FjsdAwRERERAcjMxGCx65vfonQqJMKJtWzISiGzhy4kiYiIiISEzp3JYBe/fL/L6SRNToVEONmzh42F7elyZIrTSUREREQEID6ezJTdbPgh1+kkTU6FRDjZsIFcEkns1cnpJCIiIiJSqku7fN8M/c2M+siEk3794KR8ONk6nURERERESmVmwrxF8U7HaHK6IhFG9ufGk3qYG1Kb56InIiIiIqEo84hENma3goICp6M0qZAoJIwxZxpj1hpj1hljbq+l3YXGGGuMGdSU+ULFxhfmknnwe6djiIiIiEgFbfu0ZicZsHmz01GalOOFhDEmGngSGAkcAVxsjDmimnYpwE3A4qZNGDo2vPwlmd++43QMEREREakgqlsmFtPsColQGCMxBFhnrV0PYIx5HTgP+KFKu/uBfwC3NG280LHh4j/Rv/tBp2OIiIiISEXHH0/UyCiKh8UQ7XSWJhQKhUQHYEuF77cCx1RsYIwZCHS01r5njKlXIVFSUkJOTk7gUvopGI+99qd4RpxuHP25pOH0uklVOiakOjoupCodE+EhvX08a9YcoFOn4E+KEyrHhONdm+pijIkC/g38oR5trzbGLDPGLNu9e3fwwzUhk5XF9g9W0Wnvt05HEREREZEqMrcuYNt/3nU6RpMKhSsS24COFb4/rPS2MilAX2CeMQagLTDbGDPKWrus4o6stc8AzwAMGjTIpqQ4v3BbwDJ8+y0lW7JpkRcNIfBzScOFwnEpoUXHhFRHx4VUpWMitPXctYRte9xN+jo5fUyEQiGxFOhhjMnEV0CMByaUbbTWZgNtyr43xswD/li1iIh0dsNGoKVvomIRERERCSmZz9zBu83rgoTzXZustUXADcBHwGrgDWvtKmPMfcaYUc6mCx1ZK38hnV3QubPTUURERESkisxMmt3q1qFwRQJr7fvA+1Vu+3MNbU9uikyhZsMPuWSm7IaEBKejiIiIiEgVrdYtYe8HJfBTa+jRw+k4TcLxKxJSP5s3ltApLc/pGCIiIiJSDVNcBLt3w7p1TkdpMiokwsSuvdFkpAd/OjERERERaYD0dADsriyHgzQdFRJhIis7nrQMvVwiIiIiISktjWQOcHDbfqeTNBmdmYYDa8nyJpF2WLzTSURERESkOm43aVF7yNqc63SSJqNCIhx4POyKaktap0Snk4iIiIhIdYwhLcnLrm2FTidpMiExa5PUITWV3DPPx/XHEqeTiIiIiEgN0lMLyNrVfM7XdEUiTBgDROnlEhEREQlVaa2KydptnI7RZHRmGg6+/BL7zXLYts3pJCIiIiJSg7R0Q9b+OKdjNBkVEmGgeMcuovbvKb0sISIiIiKhKK1nS3bFd3Q6RpPRGIkwsHfY+bS+CGjvdBIRERERqUn6fTeQdcDpFE1HVyTCQFYWpKU5nUJEREREatOyJezb53SKpqMrEmEg64H/kbbWDYxzOoqIiIiI1CBq6WLsF4XwQys44gin4wSdCokwsOvbHaR5dzodQ0RERERqEx+PjbJQ2DzWklAhEQay9kbTuWOx0zFEREREpDYDBsCxQH+ngzQNjZEIA1meeNIz9FKJiIiIhLqkJDh40OkUTUNnp6HOWrIOukhrH+t0EhERERGpjbWkffIqWXf/x+kkTUKFRKjbv59dtg1pnRKdTiIiIiIitTGGtMLtZP3scTpJk1AhEeqysvDiIumwlk4nEREREZE6pLnz2bWzxOkYTUKFRKjLyvL9m57ubA4RERERqVN6qyKydhunYzSJBhcSxpg4Y0w7Y0yrQAaSKrKyMFitSCciIiISBtIyDFn7m8fYVr8LCWPMZcaYpcBBYCvwcIVt5xtjXjXGZAYwY7NWEhWDSUmGjAyno4iIiIhIHdLaxpB1oHmMbfVrHQljzFTgUsAAB4DkKk3WAuOB5cA/A5Cv2dt7/Dm0PB/o4HQSEREREalLWscEsgrckJ8P8fFOxwmqel+RMMZMBC4DvgUGAalV21hrfwC2ACMDFbC5y8pSryYRERGRcNGqUzJ7aP3rONcI5k/XpslADnCutfYba62tod33QJfGBhOfrL88RfonrzgdQ0RERETqITqjDSVENYtCwp+uTUcCi6y12+potx9o2+BEUklWVAZpbXKdjiEiIiIi9dGtGyY9H+LinE4SdP4UErH4xkXUJR0obFgcqWrXSRfSsaPTKURERESkXvr3xw4G+jgdJPj86dq0GehbWwNjTDS+p+3nxoSSX2mMhIiIiEh4SUwEr9fpFMHnTyHxEdDdGPObWtr8P6Ad8F6jUomPtWQ9+CzpM6c4nURERERE6sNa0j+cRtYd/3Y6SdD5U0j8E99g6+eNMQ8YYwaW3p5gjDncGPNn4N/AHuDxAOdsnjwesvKSSUstcDqJiIiIiNSHMaT1zSCrwwCnkwRdvcdIWGu3GmPOB2YAt5V+WWBc6ZcBPMAYa+2uIGRtfrKyOEAySR1USIiIiIiEi7TfnEFWV6dTBJ9fK1tba+cCR+BbzXoVkAsU4BsT8TjQ11o7P9Ahm61dvnrMpGuQhIiIiEi4SGtjydqS53SMoPNrZWsAa+1Ofr0iIcGUlQVEa7S1iIiISBhJf+avLPsuFq6+3ekoQeXXFQlpWiW/ZGGwkJ7udBQRERERqae0jCiyPAlOxwg6FRIhbN+WA7Rkn65IiIiIiISRtI4JZBW4oSCyx7nW2LXJGPNZI/ZrrbWnNuL+AmRtySMtZr9vMmIRERERCQutOyezh9aweze0b+90nKCpbYzEyTXcbkv/NbXcbpFGy9peSFpyrtMxRERERMQP0RltKCbaN961mRYSw6u5bTRwE/A18DKwsfT2LsBvgKOBx4C3A5SvWctq35/0ISokRERERMJKejqQXT4DZ6SqsZCoOo2rMWYocANwi7X2X9Xc5TFjzO+Bh1AhERC7jh1Fhw5OpxARERERv6SlYdhfOgNn5PJnsPVdwA81FBEAWGsfAX4A7mxsMIGsXVbjrEVERETCTVoa8eSTt32v00mCyp9CYhDwfT3afQ8MblgcKWctWX/9L2nTn3A6iYiIiIj4o2VL0sxusjZHdhd1fxakiwU616NdZz/3K9UpLiar+3GkD43sS2IiIiIiEScqivQTerKrWxwdnc4SRP6c8H8HHG+MOcta+351DYwxI4HjgYWBCNesxcSQ07U/yaOdDiIiIiIi/kobdwpZ3Z1OEVz+dG36J76pXWcZY54zxpxijMks/RpujHmWXwdZPxzooM1OXh4UFvhWthYRERGRsJLWopCsnz1Oxwiqel+RsNa+bYy5HfgbMKn0qyIDlAB3WmvfDlC+5uuTT+BjA1+3hUGDnE4jIiIiIn5Ie/ZBVvwQC9ff4XSUoPHnigTW2ofwDaSeCqwHCkq/NgAvAEOstX8PcMZmye7K8l2N0LRNIiIiImEn/eJT2TXkXKdjBJXfg6KttSuAKwMfRSravzmbFsSqkBAREREJQ2mjTyBridMpgsuvKxLSdHZtzictZj+4XE5HERERERE/tY71sOfn/b5xrxFKhUSIytpRRFpyZM89LCIiIhKpYj77mKL5X8DatU5HCZp6d20yxnzmx36ttfbUBuSRUlm7LGktCp2OISIiIiINkZYG5EBW5K4J5s8YiZPr0cbim71Jc5Y2UtaeKNqmFTsdQ0REREQaIj0dFRK/Gl7D7VH4VrM+G7gQ+AfwYSNzNXu79sdxZJ/I7VMnIiIiEtHS0ojnB/J37CXe6SxB4s86EvPraDLVGHMd8G/grUalau6sJeugi7T2uiIhIiIiEpZatSKN3WRtPMhhTmcJkoAOtrbWPgVsBO4N5H6bnZISso4YRtq5xzqdREREREQaIiqKNNcBsrYVOJ0kaPxeR6IevgdOCcJ+m4/oaDyd+uIe5XQQEREREWmodHc+u3aWOB0jaIIx/WtbIDEI+20+vF7Iy8UUFzmdREREREQaKK1lEVm7jdMxgiaghYQxZjxwPLAmkPttdj77DD79FL75xukkIiIiItJAaWmQtS8YHYBCgz/rSDxfy+ZkoDfQp/T7/zQmVHNn+x4JRwLdU5yOIiIiIiINlPabM8j6INXpGEHjT4k0qR5tcoD7rLVTG5RGAMhu0ZkW/YBWTicRERERkYZKP/cYdi12OkXw+FNIXF7LtgJgG7DUWpvbuEiStWQDaSQDaU5HEREREZEGam32sufHEshNgsTIG0LszzoSLwYziPxq1yOvkLYwDrjV6SgiIiIi0kCxn39K0YIEWNcFjjzS6TgB588YiaHATmvtj3W06wG0s9Z+3thwzVXWLyWkpUbunMMiIiIizcLQoXBMFHRJcDpJUPgza9M84LZ6tLsVmNugNAJA1p4o0lpH7pzDIiIiIs1CRga2TRqkROYEOv5O/xq5E+GGkKzsONIz9FSLiIiIhLWiIuKytpH/zSqnkwRFMBakSwc04LqhrGXXARdp7WOdTiIiIiIijWEMaUveY/e0951OEhS1jpEoHRdRUdtqbqu4r8OB04HVAcjWPB08SFZxS9I6alVrERERkbAWHU1a4gGythXQweksQVDXYOt5gK3w/RmlX7UxwH8bkal5y8oim1RSO+Y7nUREREREGinNnUfWzmKnYwRFXYXE5/xaSAwDdgFramhbtpbELGvtu4GJ1wxlZQFgMtIdDiIiIiIijZXesohdWZE59rXWQsJae3LZ/40xJcAH1torgh2qWSstJEjTYnQiIiIi4S6tjWXlj5E59tWfla2HAzuDFUR8bP8BMMBCzySno4iIiIhII6W1iyHr63inYwSFPytbzw9mEPHxpHTAfTjQ0ukkIiIiItJYaYfFk5WbAsXFEB3tdJyACsb0r9IIWV/+SHrxDqdjiIiIiEgAtOnkYjetYc8ep6MEXI2FhDGm2BhTZIzpWeH7+n5p7tIG2vXEG6R99JLTMUREREQkAOLataaQWNi1y+koAVdb1yZD5ZWs/RluHplD05tA1qgrSeun9fxEREREIsIpp8AJ0dA1wekkAVdjIWGtjartewmOrJh2pB3jdAoRERERCYi0NN/YV5fTQQJPxUGIyfp4OWm7tTC4iIiISEQoLCRm+2YKln7rdJKAUyERYrJmfUH6V287HUNEREREAsEY0r75iN2vfeJ0koDzZx0JCTavl11FLUnrpLHqIiIiIhEhJoa0a8eQdVEC7Z3OEmB+XZEwxnQxxvzXGLPOGOPVrE0BlpXFflrQoqPb6SQiIiIiEiBpPVuSdSDR6RgBV+8rEsaYPsAXgJu6Z2XSrE0NUTotmMlIdziIiIiIiARK2s+LyFp5AEaMcDpKQPlzReJvQCrwAXAMkGqtjarpKyhpI11Wlu/ftDRnc4iIiIhIwKSv+Jis/1vidIyA8+eEfyiwETjfWrvUWpsTnEjNmAoJERERkYiT1jaaXZ7IW0fCn0IiHlhqrS0MVpjmLm/3ARLIgxYtnI4iIiIiIgHSok0M2fnxYK3TUQLKn0LiR3xdmyRIPFn5uPGAW4OtRURERCKFOy0eT0ky5Oc7HSWg/Ckk/gcMNcZ0CVKWZs+zpxB3tBdiY52OIiIiIiIBkpKWgAc3ZGc7HSWg6j1rk7X2KWPMEGCOMea3wEfW2pLgRWt+si+8glTUc0xEREQkksS2TKaIGPB4ICPD6TgB48/0r+tL/9sF+D+gyBizA6iumLDW2m6Nj9e8eOLTcPd2OoWIiIiIBFRZt3WPx9kcAebPytZdKvzfALFApxraRtZIkibi+XgR7pxEoL/TUUREREQkUNxuIKf5dm0CMoOWQgDInjGH1LgSVEiIiIiIRJCyQqK5XpGw1m4KZhABz1U3k96pwOkYIiIiIhJIffvCaUfA2dFOJwkof65ISJB5ily4O7mcjiEiIiIigRQXR2wyFFiIczpLAPkz/asEmef9L3Cv/MrpGCIiIiISSNbiXr8cz+x5TicJqIbM2lSXAmA3sAyYZq39piHBmh1ryf5qJakDc4DjnU4jIiIiIoFiDKkrv8IzN482Y052Ok3ANGTWJotv1qbqlG3rie9s+AZjzF+ttfc2NGCz4fXisSm4W+c5nUREREREAsx9+3V4xtR0Ch2e/OnalAk8jG/diDeA0cBR+KYYOg+YXrrt38BJwD1APnC3MWZk4CJHKI+HHFJITtcYCREREZFI4041kTZpk19XJPoDfwAusNa+U2Xb98C7xpjXgFnAF9ba+40xX+NbvO5a4INABI5Y2dmUEEV0ixSnk4iIiIhIgKV++T7Zq3Jh6IVORwkYf65I3AIsqqaIKGetnQ0sAv5Y+v37wGpgcGNCNgtlJWrZyociIiIiEjHcPy7Ds+gHp2MElD+FRD9gQz3abQCOrPD9WqCVP6GapbJCIjXV2RwiIiIiEnBuN3gONN8xEgC9GtDGAl4/H6f50RUJERERkYiVmgrZB2OdjhFQ/hQSS4GBxpiJNTUwxlwGHA0sqXBzZ+CXhsVrPuz+bN9/VEiIiIiIRBx3qxg8eZG0HJ1/g63/CpwMPG+MmQC8BmzCd8WhM3AxcBq+mZv+BmCMSQcGAFMDFThS5e05SCLJ6tokIiIiEoHcrWLw5MeDtWAio4tTvQsJa+08Y8wk4Gl8BcOIKk0MkAtca62dV3pbPPD/AC3XXAfP+KtxLy+G1Hino4iIiIhIgLnbxJJtUyA3F1yRMd2/P1cksNa+bIz5DLgKGAp0KN20HfgceM5au7VC+y3AcwHKGtGyc+NIbYv/o1ZEREREJOS50xPx4PaNi22OhQSAtXY7cF8QsjRrnjc/wv19MnCC01FEREREJMCiW6RQQpSvkGjb1uk4AaHPv0OEZ8G3uL/7wukYIiIiIhIMbjdERcHBg04nCRi/r0hIcGT/v1tJ3WOdjiEiIiIiwXDOOXCWgaOcDhI4fhcSxpixwBigJ+DGN8i6Kmut7dbIbM2KxwPu1MgYwS8iIiIiVUTITE0V1buQMMZEAW8B51F98QC+qWBN6b/iB8+Ls2jbJxbGnuN0FBEREREJtD17YPlm+GwfnHKK02kCwp8xEtcAo4FvgdOBmfgKhl7A2fjWlQB4AOgauIjNQ/bX60hdv9zpGCIiIiISDMYQv28n+Vt2OZ0kYPzp2nQpkAeMtNb+Yoy5BMBa+xPwE/CBMWYO8CwwH99idVJPnrw43C0LnY4hIiIiIsHQqhXucSPxnAVpTmcJEH+uSBwOLLTW/lL6vQUw5tcOX9baF4BVwC3+hDDGnGmMWWuMWWeMub2a7dcYY743xqwwxnxhjDnCn/2HPGvxFCTgbhNZy6aLiIiIyK/cpctIRAp/Col4YGeF7/NK/02t0u574Oj67tQYEw08CYwEjgAurqZQeNVae6S1dgDwEPBvP3KHvgMH8JBCaptYp5OIiIiISJCkznmL7MemOh0jYPwpJHYAGRW+Lysqeldp1xbw54x4CLDOWrveWlsAvI5vQHc5a23F2i2JSBvM7fGQQwpJbRKdTiIiIiIiQeL+ZR2elZudjhEw/oyRWIvvikGZhfhmaLrVGHOhtdYaY04ChgEr/NhvB2BLhe+3AsdUbWSMuR64GYgD6hzqXlJSQk5Ojh8xAsufx47avh2LIT8+liIHM0twOXk8SmjSMSHV0XEhVemYiBzJiYXs21fc6Nc0VI4Jf65IfAh0NMYMLv3+M2ANvqsH240xXwNz8BUXTwc0JWCtfbJ0bYrbgLuqa2OMudoYs8wYs2z37t2BjhA8pZ3lrNvtcBARERERCZaUpGKyc/w5/Q5t/lyReAXYDXgArLXFxpjzgBlAX3zdnkqAJ621z/mx321AxwrfH1Z6W01ep4ZCxVr7DPAMwKBBg2xKSoofMYKjXhmKiwFwtW0LIZBZgisUjksJLTompDo6LqQqHRPhr1XLKH7eFB2w19LpY6LehYS1dje+YqLibT8B/YwxvYBWwE+l7fyxFOhhjMnEV0CMByZUbGCM6VH6WOBbs+InIkl+PkTH+Ibyi4iIiEhEcqcaPN7ImVzHnysSNbLWrm3EfYuMMTcAHwHRwPPW2lXGmPuAZdba2cANxpgRQCGwD5gYiNyhwp5zLpwJHOl0EhEREREJFnerGLJzI2e6/4AUEo1lrX0feL/KbX+u8P+bmjxUE/J6ISnJ6RQiIiIiEkyprWPwFMSDtfDrUmxhq8ZCwhgztDE7ttZ+3pj7NyeeaW/j/joV7MkRcVCJiIiIyKHcrWPx2JSI+RS5tisS82j4eg22jn1LBZ6Ne3Hvz1YRISIiIhLBkrumk5PcFnJzI76QKPMzkB/sIM2ZZ8wVpIb/sSQiIiIitYi+dAIlbwBtnE4SGHUVEgbfbEyvAy9Ya78OfqTmJztbEzaJiIiISHipbUWMnsCDgBe4DlhijPnWGPN7Y0xak6RrJjyPPIf73VfqbigiIiIi4WvVKvjyS1i0yOkkAVFjIWGtXWetvRPojG9y0jeAHsC/gK3GmLeNMecZY6KbJmrk8qzdgXv3eqdjiIiIiEgwxcVBVJRv1qYIUOca3dbnY2vtxUA7fFcnVgCjgJnAdmPMv40x/YKaNIJlH4gm1R0ZB5SIiIiI1KBHDzjuOOyxxzmdJCDqLCQqstZmW2unWGuPAfoADwNFwE3AJ0HI1yx4cmNxt/DrpRARERGRMJSQAPkRMo1RY85etwMbSv81pV/SAJ7cWNyt1ENMREREJKKVlOD+6E2y/zHF6SQB4ddaD8YYA5wOTALOA+KBYuBd4NlAh2sWSkrILkwktXWh00lEREREJJiiokjN3Yln80EynM4SAPUqJIwxPfEVD5cC7fFdfVgJTAVettbuClK+yHfgAB7cuNvkOZ1ERERERILMnVCAZ2+R0zECosZCwhjjBsbjKyCOwVc87AOeRmtKBE52NgdJIqlNotNJRERERCTI3ImFZO+PjEl2arsisRNf16US4CPgeWC2tbagKYI1Gx4PAKZFqsNBRERERCTYUpOK8GRHfiGRAFjgW3wzM10GXOYbJlEna609r/HxmgFrscluSE9xOomIiIiIBJk7xeLJiYw5iuoaI2GAgaVf/oiMMqsp9O0Lw4CTnQ4iIiIiIsHmdsPO3X7NdxSyavspLm+yFCIiIiIizYA71eDxRsa0/zUWEtbaF5sySHNlZ8zEfJUOWb0gLc3pOCIiIiISRKkto/DkxjkdIyC0nLLDDhbFkxRX6FvmUEREREQimntAVzzteoEN/5EAkdFBK4x5Tjob9yhAY61FREREIp776vFkf4NvJHKY0xUJh2VnQ6pmfhURERFpFpKT4cABp1MEhgoJh3n+8gju6f9zOoaIiIiINIGoD97DfvwJrF3rdJRGU9cmh3m25eAuyHM6hoiIiIg0hfbtIT0qIsbHqpBwWLbHkJpU5HQMEREREWkKRx0F/YHOTgdpPHVtcpgnx+BOCf9R+yIiIiJSXxaKi50O0WgqJBzmORiN2+10ChERERFpElu2YP/vfezU8F+yrd6FhDHmdWPMicEM0xx5cmNwt4yM1Q1FREREpA4pKbg4SO7ug04naTR/rkhcBMw3xnxrjPl/xpikYIVqTrLzEkhtqQtDIiIiIs1CSgqpZOPJync6SaP5cwZ7E/AjcCTwFLDNGPOYMaZ3UJI1B8XFeAoTcLeOdTqJiIiIiDSF6GjcsXlk7wn/yXbqXUhYax+31h4OjADeBlzAb4FVxpg5xpjzjTH6aN0fOTl4cONOi3c6iYiIiIg0EXdCPp59zXCwtbX2M2vthUAX4H5gJ3AK8BawyRhzlzEmI6ApI1VxMQfTu+Lq29XpJCIiIiLSRFITC5tnIVHGWrvdWnsPvllwxwGfAx2Av+ArKF41xgwOTMwI1bo1DB6MOedsp5OIiIiISBNxJxXj8TidovEC0RXJAPFAQoXv44DxwCJjzGsamC0iIiIi4uNOsWTnhP+IgAb/BMaYLsaYvwPbgBeBY4BFwAQgHfg9sB3fbE8PNz5qBPrkE5jzCaxY4XQSEREREWkiqW6L50D4FxIx/t7BGHMWcB1wBhAN5AMvAY9ba7+u0PQxY8wLwArgfODaRqeNNOnpkGGgTSunk4iIiIhIE3EPHYAnJvw77NS7kDDG3Ar8P3yDrA2+KxFPA89Ya3dXdx9rrccY8zlwaeOjRp6SI/tj+gGHOZ1ERERERJqK+zejyC5wOkXj+XNF4u+l/34BPA7MtNbWZ7j5SnwDsaWKg/sLSUmOwVeXiYiIiEhzkJqKb9YmGwUmfM8D/emc9TxwlLV2qLX2zXoWEVhrH7bWDm9YvMiW/ZdHcU9/BkpKnI4iIiIiIk3EPe0JPM+9ATk5TkdplHpfkbDWXhXMIM2RZ08h7tg8iAr/wTYiIiIiUj9Jw4dwoFdLiI52Okqj1PsM1hhTbIx5rh7t/meMCf81v5uAZ18x7oQI6CAnIiIiIvVmjhmC7d4DksJ7wLU/H4Ub6t+ZP3w7ezWh7GxIdRU6HUNEREREmlJBASbPC16v00kaJRh9apIBnR3Xgyfb4k4K/+XRRURERMQPX3+N/fQzWLDA6SSNErBCwhgTZYzpA5wCbA3UfiOZ50AU7hTrdAwRERERaUput+9fj8fZHI1UayFROi6i2BhT9rH5xIq3VdleCHwHtAFmBTl3RMg+GENqqtMpRERERKRJud0YLHZ/ttNJGqWuWZsqjnWw1D72oRDfInWzgLsbmatZ8HhjcLfQuHQRERGRZsXtJomDeHd7Cefh1rUWEtba8isWxpgSYKq19oqgp2omPPnxuFtpOImIiIhIs5KSghsP2bsLI7eQqOIvwPJgBWl2rMXTaxCpw8L7kpaIiIiI+CkqitTYXDy7C2jvdJZG8GdBur8EM0izYwzZhx+H+yKng4iIiIhIU3MnFODZF96zd/pzRUICqbiY3IOQEB+Flt0QERERaV7criKy94f37J01FhLGmM/wDbCeaK3dWvp9fVlr7amNThfJfvgBPtiEmZUPF17odBoRERERaULupOJwn/211isSJ+MrJFwVvq+v8C6vmkJaGvRKhH56qkRERESam9RzTiI7P8HpGI1SWyExvPTfzVW+l0Bo2xa6Az2cDiIiIiIiTc197jA2r3Q6RePUWEhYa+fX9r000v79kBsHRXEQo6EqIiIiIs2JGw+eDYVAa6ejNFitK1tL8BS/9gZRn82BXbucjiIiIiIiTSz15SfxPP2K0zEapd4fhRtjWgJHAuustdtraNMB6AZ8Z63dH5CEEerALi8p5IDb7XQUEREREWli7vFnkb052ekYjeLPFYmbgLlAu1ratC1tc0NjQjUHnj2FuMmBpHBez1BEREREGsJ9Un88ad2cjtEo/hQSZ+O7GvF1TQ1Kt/0MnNPYYJHOs6eQ1Pg8MFpDQkRERKS5cXl34928G7xep6M0mD+FRBfgx3q0WwtkNihNM5K9rwR3QoHTMURERETEAWbeXOyiRbB+vdNRGsyfQiIFyKlHuxwgtWFxmg/P/hLcriKnY4iIiIiIE1JLT5fDeFU6fwqJnUDferTrA+xuWJzmw+MBd3KJ0zFERERExAllE+40k0LiS6CPMeasmhoYY0bim9npi8YGi3TZB6JJTVEhISIiItIsud0YLGRnO52kwfwpJB4r/fc1Y8xkY0x82QZjTLwxZjLwGmCB/wQwY0TyHIzWzK8iIiIizVVp16aS/c3gioS1dglwF76xElOAbGPMj8aYH4H9pbe5gXustV8FIWtE8Rx7Ou7zhjsdQ0RERESc4HaTzAEO7s51OkmD1XtBOgBr7YPGmDXAPUA/oHuFzd8Bf7HWzgpgvoiV3W0gqZokV0RERKR5SkoiFQ+erHxSnM7SQH4VEgClhcIsY0wG0BlfV6bN1tpfAh0uYpWU4NmcjZtofBdxRERERKRZiYrCHZdH9p4iOjidpYH8GSNRibX2F2vtEmvtUhURftq3D8+sT3G/85LTSURERETEIe6EAjz7ip2O0WB+X5EoY4xpB+UF1DZr7Y7ARGoGXC7yBx5HwjkHnE4iIiIiIg5JHXcmnr6xTsdoML8LidLZmf5I5fERGGN+Ah621j4boGyRKzER264D9HI6iIiIiIg4xX3CkWS7nE7RcH51bTLGTMU3O1OP0pu2l34B9AT+a4x5IWDpItXu3b6vMF6AREREREQax71/M55vfnI6RoPVu5AwxlwMXAZkAdcCLmttR2ttR8BVetsu4DJjzPhghI0YixbB4kWwdq3TSURERETEIe5ZL+J57i2nYzSYP12bJgMFwCnW2h8qbrDW5uO7GrEAWA5cDbwesJSRxuMB3GhFOhEREZHmK/W6S8j+NNHpGA3mT9emAcC8qkVERaXb5pa2lRoU7/MQRYkKCREREZFmzD2gK56kdk7HaDB/CgkXsLce7fYC4VtaNYGcrDzceFRIiIiIiDRj7h1r8SxbC8XhOQWsP4XENmCIMcbU1KB022B+HYAt1cjeXUiq8YArjIfpi4iIiEijpC7+GM+Cb8N2Ah5/ComPgEzgn8aY6KobjTFRwD+ArsCHgYkXmTx7CnHH50PNNZmIiIiIRLiE1kl4cYVtIeHPYOu/A+OB3wPnG2NeBTYAFl/xcDG+QmN/aVupgWdfMe6EQqdjiIiIiIiDTItU338ivZCw1m42xpwFvIGvYPhTlSYG2AJcZK3dEriIkceTbUl1qZAQERERadbcbgx5kJ3tdJIG8Wtla2vtImNMD2AsMAzoULppGzAfeLN0KlipRbbH4E4ucTqGiIiIiDjJ7caSH/lXJMqUFgovl35JA3hGjqNtizynY4iIiIiIk1JTgazmU0hI43na9qTXYKdTiIiIiIij3G6iKKFkv8evGZBCRThmDnvZy3/GvX+z0zFERERExEluNynkkLMr1+kkDVLjFQljzPpG7Ndaa7s14v4RzTP9A9xJB+G825yOIiIiIiJOSUrCTQ6e3QWkOp2lAWrr2tSlEfu1jbhvZLMWz8hxuK/2Op1ERERERJxkDO5LziX7vBg6Op2lAWorJDKbLEVzYgwHEtJI7uN0EBERERFxWnKvDhxMdjpFw9RYSFhrNzVlkGYjJ4fcH/eRsLUAenR3Oo2IiIiIOMi1fiXez/PgmEFOR/GbBls3tR074LvvMEsWO51ERERERBzmmjMb74z3nY7RIH5P/2qMiQJGAscBacBia+3zpdvSgJbAz9ba4kAGjRje0rERLpezOURERETEca4/XIvXHZ4rMviV2hgzEHgd6AYYfIOqY4HnS5uMwLdQ3Wjg3YCljCQqJERERESklKtDy/LTw3BT765NxpjOwCdAd+B94FZ8xURF7wAF+AoJqY4KCREREREp5VrxFd4PP3c6RoP4M0biTnzdlm6w1p5rrX24agNrrRf4FtC6zTVRISEiIiIipVwLP8U75yunYzSIP4XEGcBqa+1TdbTbCLRrcKJI5/VisCokRERERISkZIO3IDzHSPhTSGQAK+vRzgApDYvTDOiKhIiIiIiUciVHNYtCIgdfMVGXrsDuhsWJfIU5ecRSqEJCRERERHClROMtjHU6RoP4U0gsBwYZY2rstmSM6QUMAJY0MlfE8u4vwIVXhYSIiIiI+AqJ4jgoKXE6it/8KSSeB1zAK8aY1lU3GmPcwDOl+3wuMPEij3f0BFznnQaJiU5HERERERGHuVJj8eKC3Fyno/it3h2yrLWvG2PGAucD640x80s3HWuMmY5vDYmWwHRr7XuBjxoZvMnpuLqgNcVFREREhMTUOF8hcfAgJCU5Hccv/p7OjgMewleAnFN6W29gLJAMPApcGqhwkcj7+TJcG1Y5HUNEREREQkBsSgJFxBCOq9LVeEXCGDMDXxelD621JQDW2iLgdmPMP4Dh+AZWRwNbgDnW2l3BjxzevDM/xPV5PNDH6SgiIiIi4rSycbORVEjg68I0GthpjHkReMFa+xOAtXYfMDP48SKP9+rf4TqqyOkYIiIiIhIKzjkHRkRDz/CbAra2rk1PAfvxLS53G7DGGPO5MWaSMUZTDjWQNyoZV/sWTscQERERkVCQlAQJCRATQYWEtfYGfEXEOOBjwAIn4uvutNMY86wx5oQmSRlBvO/NxfXNF07HEBEREZFQsH07/LgW1qxxOonfah1sba0tsNa+aa0dCXQC7gR+xDew+grgc2PMGmPMrcaYtsGPG/68H3+Ba9FnTscQERERkVCQlQU//girVzudxG/1nrXJWrvdWvugtfZw4AR8VyZygJ7Ag8BmY8xsY8xoY0x0cOKGP29elNaiExERERGfI4+Es8+B8893OonfGrSagbV2obV2Mr6uTxOBefhmbzobmAFsC1TASOPNj1YhISIiIiI+UVFgDNY6HcR/jVoWzVqba619yVp7KnAmsBswQFogwkUib340rmStRiciIiIiQG4uiT98Td67nzidxG+NOqM1xiQbY640xiwAPuTXAmJLo5NFKG9hjAoJEREREfGJisK1/nu8S1Y6ncRvDZpnyhgzHLgc31oTLnxXIfKBd4DngfArqZqItzAWV0qB0zFEREREJBTExeEiF6+niNZOZ/FTvQsJY0wXYBK+MRGd8BUPACvwFQ+vlC5UJzUpLsZbHI8rJc/pJCIiIiISCozBFVuI1xN+CxbXWkiULjw3Fl8BcRK+4sEAe4FXgeettSuCGzGC5ObixYUrNfyWQBcRERGR4HDFFuI9UOJ0DL/VWEgYY54HxgBJ+IqHEnxdlp4H3rbWqn+Ov7xevLhITI1zOomIiIiIhAhXXFFkFRL4rkIArAemAlOttVuDHSiitWlD0Wkjib280OkkIiIiIhIiXPHFeA+G3/yvtU0f9DJwirW2u7X2r8EsIowxZxpj1hpj1hljbq9m+83GmB+MMd8ZYz41xnQOVpagioqC+ARISXE6iYiIiIiECFdCCd4w7PleYyFhrb3MWjsv2AFKV8F+EhgJHAFcbIw5okqz5cAga20/4C3goWDnCootW2DtWli3zukkIiIiIhIiXIkWb67TKfwXCgsaDAHWWWvXl467eB04r2IDa+1ca21ZnbYIOKyJMwbGpk3w04++f0VEREREKC0k8hu0KoOjQiFxByovYLcVOKaW9lcCH9S105KSEnJychoZreGqfez+/Sk6M4GcQV5wMJs4w8njUUKTjgmpjo4LqUrHRDNw9x/Yuy663q91qBwToVBI1Jsx5jfAIGBYDduvBq4G6NixYxMm84fxjZUQEREREQFcSYbcXFN3wxATCoXENqDiWf9hpbdVYowZAdwJDLPW5le3I2vtM8AzAIMGDbIpITCouVKGL78k5odkUvLaQ1qac6HEUaFwXEpo0TEh1dFxIVXpmIhcbZZ+QPH/QcqfR/p1P6ePiVD4aHwp0MMYk2mMiQPGA7MrNjDGHAX8Fxhlrd3lQMaAsCtXwebNUKAlOERERETEx7XtJ7zrdzodw2+OX5Gw1hYZY24APgKi8a2WvcoYcx+wzFo7G/gnkAy8aYwB2GytHeVY6AbKy84nkVxwuZyOIiIiIiIhwnX7jXjDr2eT84UEgLX2feD9Krf9ucL/RzR5qCDweopw4VUhISIiIiLlXC4iax0JCbyD2UW4yIW4OKejiIiIiEiIcM17n4MffB52s3qqkGhC3pxiXLGFYMLw2pWIiIiIBEVC1hZyf8lWISE1Ky8kRERERERKRSW7sJiw69+kQqIJeQ9aXHFFTscQERERkVBSNn5WhYTUxHvQ4oovdjqGiIiIiIQSFRJSF29+NK6EEqdjiIiIiEgoCdNCIiSmf20uvL/7Ey5rnY4hIiIiIqHE5QJywq6Q0BWJJuT1gitJMzaJiIiISAVhekVChUQT8k5/F9fc95yOISIiIiKhRIWE1MW7bjuuLWudjiEiIiIioSQ5mZiEWApNeC1arDESTch76f/DNdrpFCIiIiISUtLScJ1/BrnnQ6zTWfygKxJNyOv99cqViIiIiEgZlyvsejapkGgy1vrGSHw40+kkIiIiIhJiXJ/9H97nXnM6hl/UtampFBbi3boHV9Y+p5OIiIiISIhxFefgzQmvhYtVSDQVrxcvLlypB5xOIiIiIiIhxnXlxXjPdDqFf9S1qal4veSSSII7vEbji4iIiEjwaYyE1MzrxWKIStZoaxERERGpzPXqs3j/+aTTMfyirk1NxevFYDVtk4iIiIgcIsmzA294DZHQFYkmU3atSoWEiIiIiFThSizBm2ucjuEXFRJNRYWEiIiIiNTAlQjePBUSUh0VEiIiIiJSA1eSwZsX7XQMv2iMRFOJjoaEWEhJcTqJiIiIiIQYV3IU3nwVElKdkSNhBNDL6SAiIiIiEmpcyVF4C8Pr1Fxdm5pIYSHEhNexISIiIiJNxJUchbcg1ukYflEh0URyX3oL1+LPIC/P6SgiIiIiEmJc7hi8xXFQUuJ0lHpTIdFEvAdKcBXlQGx4VZoiIiIiEnyuHh3wpnXxdWMJE+ps00S851yE62cgvMbQiIiIiEgTcF15Md6PgXink9Sfrkg0Ea9XM7+KiIiISPViY8PqYgSgQqLJeB9+CtebLzodQ0RERERC0Ucfwdy5sG6d00nqTV2bmoh3yx5c2R6nY4iIiIhIKGrRAtwxYTXNZ/gkDXNeL7gSwmcUvoiIiIg0oWOOgaOBLk4HqT91bWoi3oNWhYSIiIiIRAwVEk3Em2twJVqnY4iIiIhIKFq5Ej78EGbPdjpJvamQaCLevChcLhUSIiIiIlKNuDgoKoQDB5xOUm8qJJqINz8Kl8s4HUNEREREQlHpOgH2oNfhIPWnQqKJePNjcCWpkBARERGRarhcJJBHviff6ST1pkKiiXgLYnAl6+kWERERkWq4XLjw4s0On1XpdGbbRLzpXXD16OB0DBEREREJRfHxuMjF6ylyOkm9aR2JJuI943xc1zqdQkRERERCkjG4YgrCqpDQFYkm4vWWj6ERERERETmEK64Q74HwWXdMhURT2LMH7zsf43p3utNJRERERCREueKKw6qQUNemJlLobk1sW9VtIiIiIlI9V/8eeLuGzyyfKiSaQuvWMLA1jHA6iIiIiIiEKtcl5+PNcDpF/ekjchERERGREOBy+cbVhgsVEk3hyy/hww/g88+dTiIiIiIiIcr15D/x/vnvTseoN3VtagoHDkBREcTGOp1EREREREKUa0hftm2OdzpGvamQaApeLxCj+V9FREREpEauC0fiXeh0ivpT16amUNbZTYWEiIiIiNTAFVeId2+e0zHqTYVEE7AHVUiIiIiISO1cTz2M95H/Oh2j3lRINIG87HwSyFMhISIiIiI1cqXEcLAozukY9aZCogl4PUW48KqQEBEREZEaudwxeEsSoLDQ6Sj1okKiCXizC3GRC3HhU2GKiIiISNNyuWPw4oLcXKej1IsKiSbgzSnGFVsIJnyWPBcRERGRppWYGucrJMJkVTpN/9oEvJl9cPVxOoWIiIiIhLKoZBcWo0JCfuU99VxcWotORERERGpTNp42TAoJdW1qAl6NsxYRERGRuqiQkKq8dz2A66mHnY4hIiIiIqEszAoJdW1qAt5+x+LKK3A6hoiIiIiEss6doasX2ofHKXp4pAxz3mNPoVUrp1OIiIiISEjr0gUOB3o6HaR+1LWpCXj35OKK1RUJEREREamFtVBUqHUk5Ffef/wH1+P/cDqGiIiIiISyXbuI/uh9ip570ekk9aKuTU3Amx+DK6nI6RgiIiIiEspSU3EN6EXuUZmkOJ2lHlRINAFvYSyuFHVtEhEREZFaJCTgGtgbb3fCopBQ16Ym4CskVLOJiIiISO1cBfvwbvjF6Rj1orPbYCssxGsTcLnznE4iIiIiIiHO9eY0vNYLx97hdJQ66YpEsHm9eHGRlKqaTURERERqlxRXgPdAidMx6kWFRLB5veSSSII7zukkIiIiIhLiXHHFKiSklNeLxRCV7HI6iYiIiIiEOFd8MV6v0ynqR4VEsJUdCS4VEiIiIiJSO1dCSdgUEuq4H2wtW0K3aOhtnU4iIiIiIiHOlWjx5obHZ/0qJILtsMPgcKCP00FEREREJNS5XLBvp3E6Rr2ER7kTzvLzobAQSsJj0IyIiIiIOMflAm9+tNMx6kWFRLDNng0ffQirVzudRERERERCnCvJhE0hoa5Nwda/PxwRA+1bOp1EREREREKca9QIvCY8JulRIRFkRV17Et0DUB0hIiIiInVwnXoc3jVOp6gfdW0KstxNu3AVZWuMhIiIiIjUyZW3F++WPWBDf8ZPFRJB5n3mZVzvvalCQkRERETq5HprGt4Z70NentNR6qSuTUHm9RThisqHGD3VIiIiIlI71wVn4v0wJSzOHUM/YZjz5hTjii10OoY0gLWWnJwcPB4PXq+X4uLiRu2vpPSqVFSULgSKj44JqY6OC6mqIcdETEwMqamptGrVipgwOCGVX8Ue2ZuCNCDW6SR105EVZN4DJbjiVEiEG2stu3bt4uDBg7Rq1Yq2bdsSHR2NMQ1fIKasEImODo8p3ST4dExIdXRcSFX+HhPWWgoKCtizZw9btmyhc+fOKkzDiMnaBVmAJwHcbqfj1EpHVZD5ConGfZItTS8nJ4eDBw/SuXNnWrRoQUxMTKOKCBERkaZijCE+Pp527doRExPDvn37nI4k/li4EBYvhnXrnE5SJxUSQeb1gitehUS48Xg8tGrVSp8IiohI2DLG0KJFCw4ePOh0FPGHq3QNCa/X2Rz1oEIiyLxecCVoxqZw4/V6SU5OdjqGiIhIo7hcLnJzc52OIf5QISFlvF5wJYb+PMBSWXFxsa5GiIhI2IuKiiofrC1hQoWElPHmRZUfDxJeNCZCRETCnf6WhaEwKiQ0a1OQeYeNxNUr9BcUEREREZEQUFpI2INeQr0MVCERZN7eA3Gd6nQKEREREQkLLhfx5FOQk0+801nqoK5NQebduAtX7h6nY4gEzKWXXooxhkmTJtWr/Q033IAxhvPPP79Bj7dx40aMMXTp0qVB92+sefPmYYzh5JNPduTxK5o6dSrGmEpfUVFRtGzZkuOPP55HH32UgoKCJsuzb98+rr32Wjp16kRsbCzGGEaPHg3ApEmTMMYwderUan+G+h4/gch47733MmjQINxuN3FxcbRr144BAwZw1VVXMXXq1EYvNum0mp5rpzTmNX7zzTc5++yzadu2LXFxcbRs2ZKePXsyatQo/vGPf7Bx48ZK7ct+P6v+TqSkpNC/f39uvfVWfvnll/L2J5988iHt6/py6r1HmjGXCxdevNlN937eULoiEWTeF6bjsgdh+O1ORxEJiCuuuIKXX36Zt956iyeeeKLW2a3y8/N59dVXy+8Xacr6HlvbtBMqZGRkcOaZZwJQVFTEzz//zKJFi1i4cCGvv/46n376KUlJSUHPMXnyZGbMmEFmZiZjxowhPj6egQMHBv1x6+uHH35gxIgR7Nixg+TkZIYMGUJGRgYHDhzg+++/57nnnuO5555jzJgxlY7jSZMm8eKLL/LCCy80WcHT3BUVFTF+/HhmzJgBwMCBAznxxBOJjo5m/fr1fPjhh7z77rskJSVxww03HHL/pKQkxowZA/hWgd6yZQuLFy/mu+++Y+rUqXzxxRf07NmTM888s9rC4MUXXwTgwgsvPOQ9rU2bNgH+aUXqkJBQWkgU0dLpLHVQIRFk3lPPwTVe065J5Dj55JPp2rUr69ev58033+Tyyy+vse0777zDvn37aNu2LSNHjmzClIEzZMgQVq9ejSuEZk3o3bv3IZ8+f/XVV5x++uksXryYBx98kL/+9a9BzVBYWMg777xDQkICK1aswF1l9dUHH3yQ22+/nXbt2gU1R20uvfRSduzYwYQJE3j66acPybhmzRqef/55zdAWAp5++mlmzJhB+/bt+eCDD+jXr1+l7dnZ2cyYMaPG46lNmzaH/E5s3LiRESNG8PPPP/O73/2O999/n9tvr/5DvbJC4uGHH9YVCHGeMbjOORXvGaE+QkJdm4LOm56Ja9ARTscQCZiKXRbq6krxwgsvAHDZZZcRExOen1u4XC569+5Np06dnI5Sq+OPP56bb74Z8HUPCbYdO3ZQVFRERkbGISfoAO3ataN3796kpqYGPUt11q1bxzfffENMTAzPPPNMtRl79+7NQw89RGJiogMJpaLp06cDcM899xxSRACkpqZyxRVX+PWBRJcuXfjLX/4CwJw5c8jPzw9MWJEm4OrfA2/77k7HqJMKiWDKy8O7fieuvVudTiISUJMmTSIqKooFCxawfv36atts27aNjz/+GKjcrclay+uvv87pp59OmzZtiI+Pp1OnTkyePPmQ/s/1sWnTJq677jq6du1KfHw8LVu2ZPjw4eVdqmqyePFiLrnkEjp37kx8fDxt2rRh0KBB3HPPPezZ8+u4purGSNx7772VplSs2qca4Morr8QYw9///vcaMzzxxBPExMRw0UUX+f1zV2fQoEGA7zkpU7H//HfffcfYsWNp27Yt0dHRPProo+Xt/HkejTF07ty5/H4Vf/ay17Ch/fa3bNnCTTfdRK9evUhMTMTtdnPCCScwdepUv7qQ7dq1C4Dk5OR6d/MqG49T9un05ZdfXulnq/izzJkzh+uvv57+/fvTunVr4uPj6dy5MxMnTmT16tXV7r/ic7Ju3TomTJhARkYG8fHx9O7dm3/84x81zvd/8OBB7rzzTrp160Z8fDwdO3bkuuuuq3SsVpWVlcVjjz3GmWeeSWZmJgkJCaSmpnLsscfy5JNPVjs2pOKYpKKiIh5++GH69+9PUlISLVq0KG9nreW5555j4MCBJCYm0qZNG0aPHs13331Xj2f6UGWvV3p6eoPuX5P+/fsDvitoe/fuDei+RYLJtWsj3uVrnY5Rp/D8iDBMmKwsCr76mdjPdkHXK52OIxIwHTt2ZMSIEXz88cdMnTqV++6775A206ZNo6SkhOOPP55evXoBvj/m48ePZ+bMmSQmJjJo0CAyMjJYuXIlzz77LDNmzODjjz8uPyGuy6JFixg5ciT79+8nMzOT888/n7179zJv3jzmzZvHhx9+yIsvvnjIPOoPPvggd955J9Za+vTpw3HHHUdOTg4//vgj9913H8OHD691cPWAAQOYOHFi+QnnxIkTD2nz29/+lueff57//ve/3HrrrURFHfq5zZQpUwC4/vrr6/Xz1sXj8QAQH3/oPB9ffvkl11xzDR06dODkk08mJyenvLuWv8/jxIkTOXDgADNmzKjUNx1o1Irwc+fO5fzzzyc7O5vu3btz5plncuDAARYtWsTll1/OZ599xrRp0+q1r7IrSPv372fq1Kn1GuuQnJzMxIkT+eKLL/j555854YQT6N79108EK/7/mmuuYevWrfTp04ehQ4cCsHLlSqZNm8Zbb73FRx99xIknnljt46xYsYKbbrqJNm3aMHz4cHbt2sWCBQu4/fbb2bp1a6UCD3xFxPDhw1m6dClut5uRI0cSHR3N66+/zscff0yfPn2qfZyPPvqI3/3udxx22GH06NGDY489lp07d7Jw4UIWL17MJ598wqxZs6pdZ8Bay4UXXsiHH37I0KFDOeKII9i8eXP59uuvv56nn36a6Ohohg0bRnp6OkuWLOGYY46ptbtjTTp16sRPP/3ElClTGDlyZLXHcENkZ2cDvkXZNNZBwonrndfwri2ESX92OkrtrLUR+XX00UdbJ3k8Hntg2TJ7DrOtffVVR7OI/3744YeA77OoqMgWFRUFfL9Oef311y1gO3XqZEtKSg7Z3qtXLwvYZ599tvy22267zQJ26NChdsuWLZXaP/744xaw3bp1s4WFheW3b9iwwQK2c+fOldrn5ubajh07WsD+7ne/q/Tcfv/99zY9Pd0CdsqUKZXuN3PmTAvY5ORkO3v27ENyL1mypFK2uXPnWsAOGzbskLaA9b2NVu/EE0+0gH333XcP2fbpp59awPbp06fG+1f1wgsv1JjFWmsvuOCC8ue3zMSJE8tz3nnnnba4uLjSfRr6PNb0ulR93BdeeKHan2HixImVbt++fbtt2bKljY6OtlOnTq10TG3evNkOGDCg2v3V5txzzy3/2QcPHmzvvPNOO2vWrEOOvfpmr2jWrFl23759lW4rKSmxU6ZMsYA9/PDDD/m9qPha3HPPPZVei/nz59uoqCgbFRVlN2zYUOl1uPnmmy1gjzzySPvLL7+U375v3z57/PHHl++zat4ffvjBLlq06JDs27dvL38+X3/99Urbyl7Xst/tn3766ZD7z5492wLW7XbbxYsXl99eVFRkf/vb35bfv+prXJsZM2aU3y8jI8NOnjzZPvfcc/abb76p9X2z7PezpuOw7D1n5MiRtT5+2WNv2LCh3pmbUmP/fgTjb5oE13/v2mxnT9lW43aPx2M9Hk+TZAGW2RrOtx0/4Q/WV0gUEp9/7isk3n7b0Sziv3q96Q4bVvfXP/9Z3rxk6FBb/Nxzvm+ysup3/7ITg7L2ZSe+a9bU7/4V2wdYXl6ebdWqlQXsnDlzKm378ssvLWCTkpJsTk6OtdbaPXv22MTERJucnFzpZKiis88+2wKVTvBrOmGdNm2aBWyXLl1sQUHBIft6+umnLWC7d+9e6fb+/ftbwD799NP1+jkbU0hMnz69xpOYspP+xx9/vF45rK2+kCgsLLRr16611113XXmemTNnlm8vO3nt3bt3tSciDX0eA11I3HrrrRawt956a7X7W7p0qQXswIEDq91enf3799vx48dbY0z5c1P21bNnT/v3v//der3eemevr7KT+5UrV1a738GDB1dbfI8cObL8ccteK6/Xa5OTky1g586de8h9vv322/Kfz5+8H3/8sQXsmDFjKt1esZB45ZVXqr3vKaecUl4MVZWXl2fbt2/vdyFhrbXPPPOMbdmy5SGvVUpKir3sssvsmmrex6orJIqLi+2mTZvsgw8+aGNjY+1hhx1WbUFUkQoJCTUvvWRtlTq/klApJDRGIohMbulsTSE024tIoMTHxzNhwgTg10HVZcq+Hzt2bHlXl7lz55Kbm1veDaI6w4YNA2DhwoV1Pv78+fMBmDBhArGxsYdsL+uPvm7dOrZt2wbAzp07+fbbb4mNja22O1KgXXDBBXTo0IGPPvqo0liSbdu2MXv2bFJSUvjNb37j937nz59f3m8/NjaWXr168dRTTxEXF8e//vWvatfsOO+886qdnaghz2MwvP/++4DvmKnO0UcfTXJyMitWrCAvL69e+0xNTeW1117jxx9/5J///CejR48u7/L0448/cvvtt3Pcccexf//+BmXeunUr//3vf/n973/PlVdeyaRJk5g0aRI7d+4sf4zqnHXWWdV2J+rduzcA27dvL7/t66+/5sCBA+Vd0qrq169ftYOTyxQVFfHxxx/zl7/8hWuvvZbLL7+cSZMmlXerqykjUO1xVFRUxJdffglQ7bEbHx9f42tYl8mTJ7N582ZeffVVJk+ezFFHHUVMTAw5OTlMmzaNo446qvw4qariWJ3o6Gg6d+7MHXfcQd++ffnuu+8qdUsTCQeun77Fu+Brp2PUSWMkgsnrBaJVSESqefP8al7y2We/ftOmjX/3r9q+Vy//7l86RiHQrrjiCp544glmzZqFx+PB7Xbj9Xp54403yreXKTuRfu+996o9iaooKyurzscuO6nNzMysdntCQgLt27dn27ZtbNu2jQ4dOpQPQu7UqVOTzNQTExPDtddey1133cWUKVN46KGHAHjmmWcoKiri0ksvJSUlxe/9VlxHIioqCrfbzRFHHMGoUaNo27ZttfcpGxxdVUOex2AoOz4GDx5cZ9s9e/b4laN79+788Y9/5I9//CMAa9eu5cknn+SJJ57g22+/5c477+TJJ5/0K+8999zDAw88QFFRUY1tysasVFXTDGBlM0tVLJS2bvVN1lHT6wO+2Ym+/fbbQ27/8ccfGT16dI2Dv2vLmJ6eXu3vyO7du8nPzycqKqrGY6ox06cmJydz8cUXc/HFFwO+MQ6zZs3iT3/6Ezt27GDixIls2rTpkOmYK47Vyc/PZ82aNaxYsYLly5dz1VVX8dZbb9X5viMSSlyfvMOOTdHwxNFOR6mVCokgMrm5WFKgCRaGEnHCUUcdxYABA1ixYgXTp08vX6DM4/HQo0cPTjrppPK2ZTPE9OrVi2OPPbbW/R5zzDH1zuDPyYETJxJXX301999/P88//zz3338/UVFR/O9//wN8A3Yborp1JOpSV+Hk9ElW2fExbtw4EhISam3b2IG4vXr14j//+Q/GGP7zn//w9ttv+1VIzJgxg/vuu4+UlBT+/e9/c8opp9CuXbvy53jChAm89tprNc4yVd3A+2AYM2YMq1evZtSoUdx6660cfvjhpKamEh0dzY8//kivXr1qzBgqU+KmpqYyadIk+vfvz8CBA9m9ezdffvklp512WqV21a0j8fbbbzN27FhmzpzJU089FbBJDUSagivJcDAv9E/TQz9hOPN6gRRdkZCIdsUVV3DjjTcydepUJk+eXP7HvOrMLR07dgTgyCOP9PskuDpln0jXNP1sXl5eeReRsrZlnwRv2bKF3NzcJjlZSktLY9y4cUybNo3p06eTkJDAjh07OPnkkzniCOfXmGnI8xgMHTt2ZN26ddx99901zkIUaKeffjr/+c9/6nUFrKKydToeeOABrrrqqkO2r1u3LiD54NfnvLapkavbtmbNGr7//nvS09OZOXPmId3aGpqxbMrm/Px8Nm/eTLdu3eqVp7GOOuoo2rRpw+7du+v9eo0ePZo77riD+++/n3vuuYff/OY3jq1rIuIvV5LBWxD6p+kaIxFMXo2RkMh3ySWXEB8fz1dffcXHH3/M3LlziY6OPmQMwogRI4iNjWXOnDkN7pNeUdl4itdee63a7iUvvvgi1lq6d+9efjLWtm1b+vXrR0FBQb2nEa1N2ZiC2rq3gG8qWICnnnqKp556CgjclK+N1ZDnMRjKFhoL1GJ6NX3SXlHZdKaHHXZYpdvj4uKAml/XsvUIyorjilavXs3y5cv9ylqbo48+mqSkJLZu3crnn39+yPaVK1dWu3ZDWcb27dtXOzbmlVdeaVCemJgYjj/++Br3UVBQwFtvveX3fut6vbKzs8u7YVV9vWpTtrr6nj17eOSRR/zOJeIUV3KUConmLt+TTwJ5KiQkorVq1YrzzjsP8A2+tNZyxhln0L59+0rtMjIyuP7669m/fz+jRo1izZo1h+zr4MGDvPrqq/zyyy91Pu7YsWPp2LEjGzZs4I477qi0kNcPP/zAPffcA1DeL75M2e233HJLtQM3ly1bVt4vvS5lJ9a19UEH30Jxxx57LIsXL2b+/Pm0b9+e0aNH1+sxgq2hz2Og3XLLLbjdbh544AGefPLJak/iV61axcyZM+u1v++++45TTjmFd999l8LCwkO2L1iwoHzV43HjxlXaVtfrWjYo+n//+x8FBQXlt+/atYuJEyfWWVj6w+VylV/1uOmmmyp9Gp+dnc11111X7Ul4jx49iIqKYuXKlYcUIC+88AKvvfZagzPdeOONAPz73/9m2bJl5beXlJRw2223NWhQ/jnnnMPDDz9c7e/+L7/8wqRJkygoKKBTp04cd9xx9d6vy+UqP4YfffTRgHyIIdIUXCnReIsOnQAj1KiQCCLP0DNxDRsMpQPoRCJV2aDqspOcioOsK3rooYe46KKLWLBgAX379mXQoEFcdNFFjBs3jmOOOYbWrVtzySWXsG/fvjofMyEhgTfeeIMWLVrw8MMP07NnTy6++GLOOOMMBgwYwC+//MKll17K1VdfXel+F1xwAX/5y1/Iycnh7LPPpl+/flx88cWcc8459OjRg8GDB9e720fZrDannnoq48eP56qrrqq2qwv8evIFvnETMTGh8UlTQ5/HQOvYsSNvv/02KSkp3HDDDXTq1InTTjuNSy65hLPPPptOnTrRt2/f8oH8dbHWMnfuXEaNGkWrVq0YNmwYF198MaNHjy5fRO6XX35h+PDh3HXXXZXue9555xEVFcWjjz7KGWecwZVXXslVV13FV199BcDvfvc7UlNTee+99+jevTtjx47lnHPOoVu3bhw4cCDgReJf//pXBg4cyIoVK+jevTvnn38+Y8aMoWvXrmzbto1Ro0Ydcp+0tDSuu+46ioqKGD58OKeccgoTJkzgyCOP5IorruD2229vcJ7Ro0dz9dVXk52dzXHHHceIESOYMGECPXv2ZMqUKVx77bV+73Pbtm3ccssttG/fniOPPJILL7yQ8ePHM3ToUDp37szbb79Ny5Ytee2116qdXaw2V155JT169CA7O5t///vffmcTcYLLHYO3JAGq+SAklKiQCCJvWmdcvTtDiJwwiATLaaedVt7No02bNtWe2ICvK9D06dOZPXs255xzDtu3b+ftt99mzpw5HDx4kIsvvphZs2ZV2++6OsceeywrVqzgmmuuobi4mJkzZ7J48WKOO+44Xn755WpXtQb485//zIIFCxg7diy7d+9mxowZLF68mJYtW3LvvffWOp1mRX/729+4+eabSU5OZubMmTz33HM899xz1bYdMWJE+XMQ7JNyfzX0eQy04cOHs2rVKv70pz+Rnp7OokWLmDFjBqtWraJr1648+OCD/O1vf6vXvvr27cu8efO46667GDhwIFu2bOGdd97hww8/xOPxcO655/LKK68wZ84ckqpMiDFgwACmT5/O4MGD+eqrr3j++ed57rnnyqdK7dq1K8uXL2f8+PFYa3n33XdZvXo1V199NQsXLgx4P/zk5GTmz5/P7bffTqtWrXj//fdZtGgRY8aMKT9uq/PYY4/xzDPP0L9/f5YsWcIHH3xARkYGH3zwQaOPwSlTpvDMM8/Qt29fvvjiCz788EMOP/xwFi5cyJAhQ/ze34wZM3jyyScZPXo0JSUlzJ07lxkzZrBy5UoGDBjAn//8Z9asWVPercofMTEx5cfNY489Vq8PKkSc5kqNxYsLypYSCFGmPv1Iw9GgQYNsxUuuTS0nJ4efP1zPtOkt+fdb1U/1J6Fr9erVHH744QHdZ9msNNX1V5bI99hjj/G73/2Oiy66iOnTpwM6JqR6Oi6kqsYeE8H4mybBVfzkFM6/oT2zdwyBaqb1zsnJAWjQFOL+MsZ8ba0dVN02XZEIosJX38b1gf+DzkQksng8Hh5++GEAbr75ZofTiIhIqItOTqSEqNIZQEOX+twEUfaFl+JqH+d0DBFxyD//+c/ywa5bt25l7Nixfq2RISIizdRZZ8GJsdC+9nV1nKZCIoi87ra4emjGJpHm6r333mP+/PmkpaUxefJk/vWvfzkdSUREwkFaGqQCoV1HqJAIpvwvl+M6mALUb+CmiESWefPmOR1BRETC0S+/wOY82BIF1axZEypUSARRwew5uGwUKiREREREpN42boTvs+B7FRLNVm6eoVULp1OIiIiISFgZMABOtXBq8KfeboyQmLXJGHOmMWatMWadMeaQVXKMMUONMd8YY4qMMWOcyNgQ3rxoLWotIiIiIv6Jj4eEBN+/IczxQsIYEw08CYwEjgAuNsYcUaXZZmAS8GrTpmuc3HwVEiIiIiLip4MHiVq3luIvFzmdpFaOFxLAEGCdtXa9tbYAeB04r2IDa+1Ga+13QIkTARvKWxBDUnJoX5ISERERkRBTVETS2m/I/Xyp00lqFQpjJDoAWyp8vxVo9ETrJSUl5av+OSEnJwdvYSwxCV5Hc0jDlJSUlK8kGiiB3p+EPx0TUh0dF1JVY48Jp8+JpAGKi3HhZV/WQWw1r12ovJ6hcEUiYIwxVxtjlhljlu3evdvpOHiLYklMadhy9iIiIiLSTMXGkmjy8OaE9gcLoXBFYhtQcV6rw0pv85u19hngGYBBgwbZlJSUxqdrqMJCcksSaNk6EUdzSINERUURHR2cIjBY+5XwpWNCqqPjQqpq6DERFRWlc5EwlBRbSHFe7a+d069rKFyRWAr0MMZkGmPigPHAbIczNV5uLl5cuNyhUKuJiIiISDhxxRXhPRDaw4MdLySstUXADcBHwGrgDWvtKmPMfcaYUQDGmMHGmK3AWOC/xphVziWuJ5eLnBNOI3HiRU4nEREREZEwo0Kinqy171tre1pru1lr/1Z625+ttbNL/7/UWnuYtTbJWtvaWtvH2cT1EBNDSUoqUe3bOp1EJKAuvfRSjDFMmjSpXu1vuOEGjDGcf/75DXq8jRs3YoyhS5cuDbp/Y82bNw9jDCeffLIjj1/R1KlTMcZU+oqKiqJly5Ycf/zxPProoxQUFDRZnn379nHttdfSqVMnYmNjMcYwevRoACZNmoQxhqlTp1b7M9T3+AmU77//nhtvvJF+/frRsmVLYmNjadOmDSeddBJ/+ctf2LRpU5PmCSX33nsvxhjuvffeet+n7Pei7GvNmjU1tvV4PCQlJZW3rXpMiEj1XAnFeL1Op6hdSBQSkcjs2UPU5s3QjP84SWS64oorAHjrrbc4cOBArW3z8/N59dVXK90vkpSdGDW1jIwMJk6cyMSJE5kwYQK9e/dm0aJF/P73v2fo0KEcPHiwSXJMnjyZKVOmEBMTw5gxY5g4cSKnnHJKkzx2fRUWFnLDDTfQv39/Hn/8cXbv3s0JJ5zA2LFjGThwICtXruTee++lZ8+ezJo1y7GcJ598MsYY5s2b51iGxqitOJg+fTreUD8bcohT7yESHlwJJSFfSKgDf5CYDRuI+uEArIqGzp2djiMSMCeffDJdu3Zl/fr1vPnmm1x++eU1tn3nnXfYt28fbdu2ZeTIkU2YMnCGDBnC6tWrcYXQ6pK9e/c+5MTtq6++4vTTT2fx4sU8+OCD/PWvfw1qhsLCQt555x0SEhJYsWIFbre70vYHH3yQ22+/nXbt2gU1R10uueQS3nzzTdq1a8eUKVMYNWpUpe1FRUXMmjWLP/3pT836qkRDdevWjT179vDSSy/xt7/9rdrBwC+88ALR0dH069eP5cuXO5BSJDy5Eize3NAuNHVFIkhK+venaPhwGD7c6SgiAVWxW0pdXRReeOEFAC677DJiYsLzcwuXy0Xv3r3p1KmT01Fqdfzxx3PzzTcD8Oabbwb98Xbs2EFRUREZGRmHFBEA7dq1o3fv3qSmpgY9S02effZZ3nzzTVq0aMGXX355SBEBEBMTw9ixY1m+fDnDhg1zIGV4S0hIYPz48Wzfvp1PPvnkkO0//vgjCxcu5PTTT6d9+/YOJBQJX647bsJ7zR+cjlErFRLBEhsL8QmQmOh0EpGAmzRpElFRUSxYsID169dX22bbtm18/PHHQOVuTdZaXn/9dU4//XTatGlDfHw8nTp1YvLkyWzcuNHvLJs2beK6666ja9euxMfH07JlS4YPH17epaomixcv5pJLLqFz587Ex8fTpk0bBg0axD333MOePXvK21U3RqKsT3mZquMWAK688kqMMfz973+vMcMTTzxBTEwMF10UmEkZBg0aBFDpk/WKYxW+++47xo4dS9u2bYmOjubRRx8tb+fP82iMoXPpldZNmzZV+tnLXsOaxkjUZcuWLdx000306tWLxMRE3G43J5xwAlOnTsVaW+/9WGt54IEHAPjzn/9MZmZmre2Tk5M56qijDrl99erVXHnllWRmZpKQkEDLli0ZMWIEs2dXP7lgly5dyp+HTz75hFNPPZXU1FRcLhfHHnvsIfcrO77mz58PwPDhw4mJiSEmJqZSV6eKx6HX6+Wuu+6id+/eJCYmMmDAgPL9zZkzh+uvv57+/fvTunVr4uPj6dy5MxMnTmT16tX1ffr8UnZVsrrXuuzDhNquXILv9XrppZc4+eSTadmyJQkJCXTr1o3rr7+eLVu2VHufir9vU6dOZdCgQSQlJdG2bVuuvPJKsrKyAMjLy+Oee+6hZ8+eJCQk0KlTJ+68804KCwtrzPPRRx8xatQoMjIyiIuLo127dlx88cV8//33h7StOI7LWstTTz3FgAEDcLlctGzZkvPOO4+VK1dWuk993kPAtxDdlClTOPHEE2nVqhVxcXFkZGQwcOBA/vCHP5T/jBKZXO1S8RI6V8OrZa2NyK+jjz7aOunAV1/Zs7r9YG1WlqM5pGF++OGHgO+zqKjIFhUVBXy/Tjn99NMtYO++++5qtz/wwAMWsMcff3z5bQUFBfaCCy6wgE1MTLQnnXSSHTNmjO3du7cFbMuWLe3SpUsr7WfDhg0WsJ07dz7kMRYuXGhbtGhhAZuZmWnHjRtnTzvtNBsbG2sBe+mll9qSkpJqsxljLGD79Oljx40bZ8866yzbvXt3C9i5c+eWt507d64F7LBhw8pvmzVrlp04caIFLGAnTpxY6ctaa5cvX24B26VLF1tcXFztc1T2c8+bN6+GZ7myF1544ZAsFb300ksWsG63u/y2spxXXXWVjY+Pt127di3/ef/73/826HmcOHGivfDCCy1gk5KSKv3sWaXveWWP+8ILL1T7M5Q9TxV99tlnNjU11QK2e/fudvTo0XbEiBE2OTm5PEd9rVixwgLWGGN3795d7/tV9Nprr9m4uLjy4+TCCy+0Q4cOLb+tumO/c+fOFrB33XWXNcbYIUOG2PHjx9ujjjqqPM+bb75Z3n716tV24sSJNiMjwwL2jDPOsJdddpm97LLL7MSJE+3q1auttb8eh8ccc4wdPHiwTU5OtmeddZa96KKL7OjRo8v3161bNxsfH28HDhxoR48ebUePHl1+XLtcLrtgwYJDMt9zzz0WsPfcc0+9n5uyPH369LHWWnvEEUfY+Ph4u2/fvvI2xcXFtkOHDrZVq1Y2Ly/Pnn322dUeEyUlJXbChAkWsLGxsfa0006z48aNs5mZmRawrVq1skuWLDkkQ9nv36233mrj4uLsaaedZi+44ALbtm1bC9h+/frZnJwce/zxx9uWLVva0aNH25EjR1qXy2UBO3ny5Gp/thtvvNECNiYmxh533HF27Nix5a9fQkKCfe+99yq1r/geddlll9n4+Hh72mmn2TFjxpQfD2632/7888/l96nPe4i1v/4eJSYm2hEjRtiLL77YnnbaabZr164WsAsXLqzX6xWMv2kSfF89tMD+7Yz51W7zeDzW4/E0SQ5gma3hfNvxE/5gfTldSHifftqew2xrK7xxSPhQIVG3119/3QK2U6dO1Z6s9+rVywL22WefLb/ttttus4AdOnSo3bJlS6X2jz/+uAVst27dbGFhYfntNRUSubm5tmPHjhawv/vd7yo9t99//71NT0+3gJ0yZUql+82cOdMCNjk52c6ePfuQ3EuWLKmUrbpCokzZSUBNTjzxRAvYd99995Btn376aaUTsfqoq5AoK9KGDh1aflvFk5U777zzkKKmoc9jbQVexcetbyGxfft227JlSxsdHW2nTp1a6ZjavHmzHTBgQLX7q8lzzz1nAdu1a9d6ta/q22+/tXFxcTY5Odm+//77lbatXLmy/Dn77LPPKm0rO3GMi4uzH3zwQaVt999/f3mRVNWwYcPKi9jq3ivKjkPADhgwwO7cubPa3LNmzap0Mm+t70R9ypQpFrCHH374Ib+vgSgkHnroIQvYp556qrzNBx98YAF7/fXXW2ttjYXEk08+aQGbkZFhV65cWX57UVGR/e1vf1t+nOXl5VW6X9nzkZGRUek9e+/eveXvP3379rUnnnii3b9/f/n25cuX25iYGGuMsRs3bqy0z6effrr85yor4srMmjXLxsTE2BYtWti9e/eW3172u1BWiK9bt658W15enj3rrLPKi/mqansP2bhxowVsx44d7bZt2w45JpYvX25/+eWXau9blQqJ8LRi7F/tncmPVLtNhUSEFxIHHn7EjuJta3fscDSHNIwKibrl5eXZVq1aWcDOmTOn0rYvv/yy/NPqnJwca621e/bssYmJiTY5ObnGP35lJxoVT/BrOmGdNm1a+Sf+BQUFh+yr7ISg6klb//79LWCffvrpev2cjSkkpk+fbgE7cuTIQ7aVnfQ//vjj9cphbfWFRGFhoV27dq297rrryvPMnDmzfHvZCX3v3r2rPf4a+jwGupC49dZbyz9drs7SpUstYAcOHFjt9qr+/ve/W8Aee+yx9Wpf1UUXXXTIiXFFb775pgXsBRdcUOn2skLiD3/4wyH3yc/PL7/ismnTpkrb/Ckkvvzyywb9TMcff7wFKp2sWxuYQmLHjh02OjraDhkypLxN2XO4bNkya23NhUTZp+vPPPPMIY+Tn59vO3XqZAH78ssvV9pW9nyUXVmr6JFHHrGAjYqKqvb9fNSoURawL774YvltRUVFtl27dhawq1atqvbnvv766y1g//Of/5TfVrGQqO5Dg8WLF5cXGVXV9h6yZMkSC9jzzjuv0X8/VEiEpx/XFNvf/776baFSSITn6McwkOspxIUXQmimFwmsq6+G7dvr395a35Akp2b6a98ennkmcPuLj49nwoQJPPHEE7zwwguceuqp5dvK+kWPHTuW5ORkAObOnUtubi5nn3026enp1e5z2LBhvPfeeyxcuJBzzz231scv61M+YcIEYmNjD9k+adIkrrvuOtatW8e2bdvo0KEDO3fu5NtvvyU2NpaJEyc26Of2xwUXXECHDh346KOPWL9+PV27dgV840dmz55NSkoKv/nNb/ze7/z586udMjIuLo4HH3yw2jU7zjvvvGpn1GnI8xgM77//PuA7Zqpz9NFHk5yczIoVK8jLyyMhISEoOQBKSkr48MMPMcYwZsyYatuUDcxeuHBhtdvPOeecQ26Li4uja9euLF++nO3btzdoAH9GRgbHH398rW22bt3Ke++9x5o1a/B4PBQXFwOwc+dOwDcAuk+fwC7H1LZtW84880zee+89Vq9eTdu2bXnnnXc48sgjOfroo2vNun79eqKiorj00ksP2R4XF8cll1zCgw8+yLx587jkkksOaXPmmWceclv37t0B6Ny5M4cffvgh23v06AHA9gpv4itWrGDHjh306dOHI444otq8w4YN48knn2ThwoX89re/rbQtJiam2iy9e/c+5LHqo3fv3qSkpPDee+/x4IMPMmHChPL3EGkeXMlRmv61ucr1FPkKCQ22jlj+npQXF/tWp6zuZC5cXXHFFTzxxBPMmjULj8eD2+3G6/XyxhtvlG8vUzYo+7333qtz3vT6DCDctm0bQI2DaBMSEmjfvj3btm0rPwEuG4TcqVMnEpvgdzMmJoZrr72Wu+66iylTpvDQQw8B8Mwzz1BUVMSll15KSkqK3/vNyMgoP2GJiorC7XZzxBFHMGrUKNq2rX4RzM41TEPdkOcxGMqOj8GDB9fZds+ePXXmSEtLA2DXrl1+Z9mzZw8ejwegxqK3TE3Hak1FQtkMV3l5eX7ngppfxzL33HMPDzzwAEVFRTW2KfvZAu3yyy/nvffe44UXXqBLly7k5+fXufBg2fHXrl27GovDigV4dQ477LBDbiv7AKO6bRW3V3wdyo7BVatWNeg9ql27dtXOTlf2mufn59e6z6pSUlJ4/vnnueKKK7j77ru5++676dChA8cddxxnn30248ePD2pBLc5zLfoM75x4yB8E8fFOx6mWCokg8eYUkxiV55u9SSRCHXXUUQwYMIAVK1Ywffp0Jk+ezIwZM/B4PPTo0YOTTjqpvG3Zp6K9evXi2GOPrXW/xxxzTL0z+LOYkxMLP1199dXcf//9PP/889x///1ERUXxv//9D4BrrrmmQfusbh2JutRVODm9KFbZ8TFu3Lg6T47i6/EHtexT8A0bNrBnzx5at27td5bo6OgGXTECX4EXDLW9jjNmzOC+++4jJSWFf//735xyyim0a9eu/D4TJkzgtdde8/VrDoJzzz2X1q1b8/LLL9O+fXtiYmLq/fw15vir7bn253Uoe907dOjAiBEjam1bdpWhoY9VX2PGjGHEiBHMmjWLBQsW8NVXX/HWW2/x1ltvce+997JgwQI6duwY8MeV0ODatBrvzxlw8KAKieYmL6cYV0zNU8uJRIorrriCG2+8kalTpzJ58uTyE9yq0z2W/bE78sgj/T4Jrk7ZJ9I1TT+bl5dX3pWgrG3Zp8RbtmwhNze3Sa5KpKWlMW7cOKZNm8b06dNJSEhgx44dnHzyyTV2n2hKDXkeg6Fjx46sW7eOu+++OyDdbvr160eXLl3YuHEj06ZN4/e//32979umTRsSExPJzc3liSeeKP/0OtSVrR/ywAMPcNVVVx2yfd26dUF9/Li4OCZMmMDjjz/Ojh07GDVqVJ1XdMqOqe3bt5Ofn19tkVh2bAbz+INf36PatWsXkPeoQGnRogWXXXYZl112GdHR0fz8889MnjyZuXPnctttt9U51bWEr7iUePKJB68XWrVyOk61tI5EkHgPWlxxKiQk8l1yySXEx8fz1Vdf8fHHHzN37lyio6MPGYMwYsQIYmNjmTNnDvv372/045b1UX/ttdeq7cbx4osvYq2le/fu5Scgbdu2pV+/fhQUFDBt2rRGZygbU1BbNxKgvC/1U089xVNPPQXA9ddf3+jHD4SGPI/BULbyeaAW0zPGcMcddwBw3333sWHDhlrbHzhwoHzV5ZiYmPJPpN96662A5KlLXFwcUPexVJu9e/cCVPsJ9erVq5tkVekrrriC1q1b07p162qLmaoOO+wwunbtSklJCS+//PIh2wsLC3nllVcAKq3lEgxDhgyhdevWLF++POhFV5n6vodU1K1bN+68804Avv3226DkktBgkkrH2YbwQAkVEkGSe9Diim34HwSRcNGqVSvOO+88AH7zm99greWMM844ZBXbjIwMrr/+evbv38+oUaNYs2bNIfs6ePAgr776Kr/88kudjzt27Fg6duzIhg0buOOOOygpKSnf9sMPP3DPPfcA8Mc//rHS/cpuv+WWW8oH+Fa0bNkytm7dWufjw6+fkNa10NegQYM49thjWbx4MfPnz6d9+/aMHj26Xo8RbA19HgPtlltuwe1288ADD/Dkk09We2K1atUqZs6cWe99Tp48mfPPP5/9+/dz4okn8u677x7Spri4mFmzZnH00UeXDzwH3yJ2sbGx3HTTTbz++uuHdAey1rJkyZLyRRcbq77HUm3Kutv873//o6CgoPz2Xbt2MXHixEYVKfU1YMAAdu/eze7du+ucMKFM2Yrsd999d6X3heLiYm699VY2b95M586daxz4HiixsbHcfffdFBcXM3r0aJYsWXJIm4KCAmbPnl3t+1dD1Pa6L1++nOnTp5Obm3vItrJjua4xMxLmXKFfSKhrU5B4vZAYX+x0DJEmccUVV/DGG2+UD0CsOMi6ooceeojt27fzxhtv0LdvXwYMGEDXrl3LVwL+9ttvyc/PZ/Xq1WRkZNT6mAkJCbzxxhuMHDmShx9+mFmzZjF48GD27t3L3LlzKSws5NJLL+Xqq6+udL8LLriAv/zlL9xzzz2cffbZHHnkkfTp04ecnBzWrl3LunXrmDt3bo2DNCs6//zzeeSRRzj11FM55ZRTyrvAPPvss4e0vfHGG1m0aBHgGzcRExNT3ifbSQ19HgOtY8eOvP3224wZM4YbbriBv/3tb/Tp04f09HT279/P999/z5YtWxg3bhwXXHBBvfZpjGH69OnceOON/Pe//2XUqFG0b9+eo48+mpSUFPbs2cOSJUvYt28f8fHxlQacDxo0iGnTpnHFFVdw8cUXc/vtt3PEEUfQqlUrsrKyWLFiBbt27eK2227j9NNPb/TPf/755zN16lRuueUWPv74Y9LT0zHGcMstt9CrV6967eN3v/sd06ZN47333qN79+4cc8wx5ObmMn/+fDp27Mjo0aN5++23G5010K677jq+/PJLXnvtNfr378/JJ59Mq1atWLJkCevXr6dly5a8+eab9Rob01g33XQTmzZt4pFHHuGYY46hX79+dOvWjbi4OLZt28by5cs5ePAgH3zwQbXjJPxV23vIpk2bGD9+PC6Xi6OOOoqOHTtSWFjI8uXLWb9+PSkpKdx3332NziAhzOUC8kO6kHB8vYdgfTm9jsSKL3bYRe9vdjSDNJzWkfBPcXFx+QJdbdq0qXY9gopmz55tzzvvPNuuXTsbGxtrW7VqZfv06WMnTZpkZ82aVen+da1XsHHjRnvNNdfYLl262Li4OJuammqHDh1qX3755WoXyiuzYMECO3bs2PIMbdq0sYMHD7b33nuv3bNnT3m72taR8Hq99uabb7aZmZnlq0BTw5zwu3btspSu3Lt9+3Zrrf/HRF0L0lWnpvUcqvL3eQz0OhJlduzYYf/0pz/Z/v372+TkZBsfH287d+5shw0bZh988MFKi335Y8WKFfaGG26wffv2tampqTYmJsa2bt3annTSSfavf/3rIQskllm3bp298cYb7eGHH25dLpd1uVy2a9eu9vTTT7ePPfaY3bZtW6X2ZetIbNiwodr9VVwvoqqnnnrK9u/f3yYmJpYfS2XtajsOK1q/fr0dP368Peyww8pXMr/55pvt/v37a3xNArGORH3UtI6Etb5F86ZNm2ZPOukkm5qaauPi4myXLl3stddeazdvrv5vaW2/b3U9X3X9zPPnz7fjx4+3HTt2LP996N27tx03bpx95ZVX7IEDB8rb1vW7UFvW2t5DduzYYR988EF75pln2i5dutjExESbmppq+/bta//whz8csphebbSORJjKzbVvT91nbTV/U0NlHQljgzR7g9MGDRpkly1b5tjj5+TkAPz/9u49WM66vuP4+2MiBQ4EEBJUsESLA0KdiqUWFSil7QiViq2IOlrvdCxM5dpCO7YMll4slzYKtaVKgVYpeItUrFopAdGKYhHLzfECgUAgQAwJFwMh3/7xPAe2657IJmfPbs6+XzNnnn3u37Pzm9/ud5/fZaOGdtTw3XLLLT3HHt8UnSPBaPwsWrSI4447jiOPPJJLLrkEsEyoN8uFum1qmRjEZ5qGaya/Zyb5VlXt22uffSQkacBWr17NmWeeCTzVHlySpM2dfSQkaUDOOOMMbrzxRq6++mqWLVvG61//+r7myJAkaZSZSEjSgFx++eVcddVVzJ8/n6OOOoqzzjpr2CFJkjRtTCQkaUCWLFky7BAkSRoY+0hIkiRJ6puJhCRJkqS+mUhIkiRJ6puJhDSF2TrHiiRpfPhZpkEykZB6mDNnzpMTAEmStLlav349z3iGX/c0GJYsqYett96ahx56aNhhSJK0SR555BG22mqrYYehWcpEQuph3rx5rFy50qcSkqTNVlWxatUqJiYmhh2KZikTCamHbbfdlomJCZYuXcqqVatYt26d7UwlSZuFqmLt2rUsX76cdevWscMOOww7JM1STkgn9ZCEBQsWsGbNGlavXs2KFSs2+enE+vXrAWyrqidZJtSL5ULdNqZMzJ07l+22244FCxZYljQwJhLSFJIwb9485s2bNy3XW7NmDdA87ZDAMqHeLBfqZpnQqDJFlSRJktQ3EwlJkiRJfTORkCRJktQ3EwlJkiRJfTORkCRJktQ3EwlJkiRJfTORkCRJktQ3EwlJkiRJfTORkCRJktS3VNWwYxiIJPcBS4ccxk7A/UOOQaPFMqFulgn1YrlQN8uEus1Umditqub32jFrE4lRkOS6qtp32HFodFgm1M0yoV4sF+pmmVC3USgTNm2SJEmS1DcTCUmSJEl9M5EYrPOGHYBGjmVC3SwT6sVyoW6WCXUbepmwj4QkSZKkvvlEQpIkSVLfTCQGJMkhSb6b5PtJThl2PJp5SZ6X5MokNye5Kcmx7fZnJfnPJN9rlzsMO1bNrCRzklyf5HPt+vOTXNvWF5ck2WLYMWrmJNk+ySeT3JrkliQvt54Yb0mObz83bkxycZItrSfGT5Lzk6xIcmPHtp51QxofbMvHd5K8dCZiNJEYgCRzgHOBQ4G9gDcl2Wu4UWkI1gEnVtVewH7AMW05OAW4oqpeCFzRrmu8HAvc0rH+AeBvq2p34EfAu4YSlYZlEfCFqtoT+AWasmE9MaaS7AK8F9i3qn4emAO8EeuJcXQBcEjXtqnqhkOBF7Z/vwd8eCYCNJEYjJcB36+qH1bVY8C/AYcPOSbNsKpaXlX/075eQ/PlYBeasnBhe9iFwGuHEqCGIsmuwKuBj7TrAQ4GPtkeYpkYI0m2Aw4EPgpQVY9V1SqsJ8bdXGCrJHOBrYHlWE+Mnaq6GljZtXmquuFw4KJqfB3YPslzBh2jicRg7ALc2bG+rN2mMZVkIbAPcC2wc1Utb3fdA+w8rLg0FH8H/BGwvl3fEVhVVevadeuL8fJ84D7gn9vmbh9JMoH1xNiqqruAM4E7aBKIB4FvYT2hxlR1w1C+e5pISAOWZBvgU8BxVbW6c181w6Y5dNqYSHIYsKKqvjXsWDQy5gIvBT5cVfsAD9PVjMl6Yry0bd4Pp0kynwtM8JPNW6SRqBtMJAbjLuB5Heu7tts0ZpI8kyaJ+FhVfbrdfO/k48Z2uWJY8WnGvRJ4TZLbaZo8HkzTPn77tgkDWF+Mm2XAsqq6tl3/JE1iYT0xvn4duK2q7quqx4FP09Qd1hOCqeuGoXz3NJEYjG8CL2xHWNiCppPUZUOOSTOsbfv+UeCWqjq7Y9dlwNva128DPjvTsWk4quqPq2rXqlpIUy/8V1W9GbgSOKI9zDIxRqrqHuDOJHu0m34NuBnriXF2B7Bfkq3bz5HJMmE9IZi6brgMeGs7etN+wIMdTaAGxgnpBiTJb9K0hZ4DnF9VfzHciDTTkuwPfAX4X55qD/8nNP0kLgV+FlgKHFlV3Z2pNMslOQg4qaoOS/ICmicUzwKuB95SVWuHGJ5mUJKX0HS+3wL4IfAOmh/6rCfGVJLTgDfQjP53PfBumvbu1hNjJMnFwEHATsC9wKnAYnrUDW3SeQ5NM7hHgHdU1XUDj9FEQpIkSVK/bNokSZIkqW8mEpIkSZL6ZiIhSZIkqW8mEpIkSZL6ZiIhSZIkqW8mEpI0yyU5KEl1/T2e5IEktya5NMnRSbbfwDXe3p53wcxFPr063oclw45FkmYDEwlJGh8PAxe2fx+nmefkUeC3gXOBu5Mc345HvtlJcnubKCwcdiySNA7m/vRDJEmzxP1V9fbujUl2BN5LM2Hi2cCuwIldh30G+Drw4IBjHKRvAC+imaxJkrSJfCIhSWOuqh6oqlOBI4ACTmhn3u485sGqurWqlg8hxGlRVY+0/8Mdw45FkmYDEwlJGiFJFiZ5IsnKJFtNccwzkyxvm/HsPV33rqrPAp9qV0/qumfPPhKd/Q6SbJ3k9LbfxaNJvt117KuSXJbk3iSPtf/DxUlePFVMSXZM8v4k1ydZneThJN9LckGSV3TGBuzWnnZbV3+Qhd2xTnGvvZNclOTOJGuT3J/k80kOneL4C9rrvT3J7kk+3v5va9v34OQkfs5KmrVs2iRJI6Sqbk/y78DhwJuA83sc9jrg2cCSqrppmkP4V5onEwclmVtV657meVsCS2iaDl0N3ABsMbkzySKa5lPrgG8Cy4DdgTcCr03yuqr6fOcFk+wDXA48B1jZXv/HNAnDm9rDvgZ8n6bfxxHABE0y9FDHpTpf95TkNcClwM8AN9H0H9kVeBVwaJLTq+pPpzj9JcAi4H7gSmABcADw1+01/uCn3V+SNkcmEpI0ej5Ek0j8Pr0TiaPb5bkDuPd17XICWEjzJf3p+GXg28DuVXVv544k76FJIm4CjqiqWzv2vRb4BPCxJC+oqh+127cBLqNJIv4BOKGqHu04bz6wB0BVXQNc0zbHmgBOqqrbn2bcJHk28C80ScSJVXV2x76DaJKZ9yW5pqq+2OMSxwKnAe+vqvXteQfSJBVHJ/mbqrrz6cYjSZsLH7lK0oipqiuAm4F9k7ysc1/bDOgA4G5g8QBuf3/H6x37PPeYHknEHODP2tUjO5MIgKpaDPwjsD3wlo5d76b5Nf+/gaM7k4j2vPvaBGI6HAXMA77amUS091lCk9hBV3OvDt8ETptMItrzrga+SPM5+6vTFKckjRQTCUkaTee0y6O7tk+un9dHs6N+dH4urJ/yqJ90b1V9rcf2l9A8Vbipqm6e4tyr2uXLO7Yd0i7Pr6rqI46N8Svt8sIp9k8+Fdq/TYy6fX6KGCeTpuduSnCSNKps2iRJo+ki4K+ANyQ5oapWJplH86v948B5AEn2BE7pcf7i9tf+fu3U8XplH+ctnWL7C9rl3m2H6A2Z3/F6suP0rb0OnGa7tMvbpth/O01StSXNU5oVXfunGgVqdbvcclOCk6RRZSIhSSOoqh5Ocj5wPPBO4EzgrcA2wCc6hmF9NvC2Hpe4nY1r+vSL7XJNe42n69Eptk/+gn8X8OWfco3OpGHQTyF62dh79vPkRpJmDRMJSRpd59J05H1PkrNpOl9PbgeebMM/nTNRT/ZTuLKqnpiG6012Ml7eazK8DbiDZgSoPYDp6gsxlbuAPWmenlzRY/9CmiZfP6a/pzSSNKvZR0KSRlRV/QD4D+DngL8E9qLpa3DVBk/cSO0QqL/Trp4xTZf9BvAAsE+S3fs4b3J0pHcmebqJ0mPtst8fySbfz7dOsf8d7fKaAfVLkaTNkomEJI22yRGDTm6Xfz/dN0jyrCSn0sy/EOAD0zUiUlU9Dvw5TROnxd2jULX33yLJa9r+HpM+QjMy1SuADyXZsuuc+Un277rUXe3yRX2G+U80Tbn2T/LervscyFPzQJzV53UlaVazaZMkjbYvAd+laeKzhma+g421U8fM1AG2pWm282Kaz4NHaIY4/eAm3OMnVNWiJLvR9Pe4Nsl3gB/QPEHYBdiHZv6HQ2n7SVTVmiSH08zhcAzwxiRf5akJ6fYBLub/N3v6DHAQzZwUXwJWtdtPrqoHNhDfPUl+F7gEWJTk3cCNNKMtHUDzo9vpVfWFTXwrJGlWMZGQpBFWVZXkyzSJxEVVtWYTLjfBUx2z19EkJitovoBfCVxcVas24fpTqqoTkiym6efxSuDVNB20lwOfo5l87itd51zXzptxPPBbwG/QdGy+G/g4zfwTnc6hmQ/izcBhNBPMAZxO07xqQ/F9Nsm+NE9+DqaZJXsNTSL3oe5ZtyVJkMEPzy1J2lhJtqDpeLwzsPcG5mKQJGlG2UdCkkbbMTRJxBdMIiRJo8QnEpI0YpLsAfwhTRv9VwFPAL9UVTcMNTBJkjrYR0KSRs9zgHcBa4EbgPeZREiSRo1PJCRJkiT1zT4SkiRJkvpmIiFJkiSpbyYSkiRJkvpmIiFJkiSpbyYSkiRJkvpmIiFJkiSpb/8HOZDzeX4L8vQAAAAASUVORK5CYII=\n", + "text/plain": [ + "<Figure size 936x720 with 1 Axes>" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "fig = plt.figure(figsize=(13, 10))\n", + "ax = plt.gca()\n", + "fontsize = 22\n", + "\n", + "vel_profile_ch_with_cm = ch_with_cm.velocity[0.5, :, 0]\n", + "vel_profile_ch_normal_srt = ch_normal_srt.velocity[0.5, :, 0]\n", + "\n", + "plt.xlabel('y-Direction', fontsize=fontsize)\n", + "plt.ylabel('Velocity Magnitude', fontsize=fontsize)\n", + "\n", + "plt.grid(color='black', linestyle='-', linewidth=0.1)\n", + "\n", + "plt.rc('font', size=fontsize)\n", + "\n", + "plt.plot(vel_profile_ch_normal_srt, linestyle='-.', color='r')\n", + "plt.plot(vel_profile_ch_with_cm, color='b', linewidth=0.7)\n", + "\n", + "plt.legend(['Velocity Profile Standard SRT', 'Velocity Profile Central Moments'], fontsize=fontsize)" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[<matplotlib.lines.Line2D at 0x7f9b73e689a0>]" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxAAAAJ0CAYAAACC68UFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABeaUlEQVR4nO3deZwkdX3/8ddnZ0/2mmWXYxHkFgQP1MVbwSuoUVFj8IxoYg41UWOMJv6i0XjEeMYzhkRFY0i8It6aeIFnZPFGMbKcwi6wy87e93x+f1QX0zvM0T3TM9Xd9Xo+HvWo6e7qns9M1/TUu77fb30jM5EkSZKkVsypugBJkiRJvcMAIUmSJKllBghJkiRJLTNASJIkSWqZAUKSJElSywwQkiRJklrW1wEiIk6JiBdHxEcj4sqIGI6IjIinVF3beCLi3hHx8oj4eERc06g3I2LNBM85JCIeHxHvi4ifRcT2iNgTEddGxEci4l6z+TNIkiSpf82tuoAZ9nzgxVUX0aZXA+e2+ZxnAP/S+Po64KvAfuCewO8Bz4iI52fmv4zzfEmSJKklfd0CAfwceAvwVOAk4JJqy2nJ94DXA08CjqEIBJPZB3wQuHdmHpeZT8zMpwB3Af4CGADeGxF3maGaJUmSVBNRp5moI+KbwFnA72bmJysupyURcS1wLHBmZq6d4mt8FXgE8OrMfF0Hy5MkSVLN9HsLxJRExOLGOITLImJrROyKiCsi4jURsaTq+qbgR4310ZVWIUmSpJ7X72Mg2hYRRwNfAU4DbqXoUrQbOBP4W+BJEXF2Zm6ursq2ndxYr6+0CkmSJPU8WyCaREQAH6cID+8Bjs3MczLzXOBE4KPAPYB3VFdleyLi7sBvAwl8uuJyJEmS1OMMEAd7NPAA4PvAizNzV/lA4+s/AW4BnhkRK6opsXWN7lYXUbQ0fSgzf1JxSZIkSepxdmE62GMb609l5vDoBzNzR0SsbWx3JvDfABFx5/LrNn06M/96qsVOJCLmAZ8A7gb8GHjRTHwfSZIk1YsB4mAnNNZviYi3TLLtYU1fzwdOmcL3Wz2F50wqIuYC/0nRovJL4JzM3DET30uSJEn1YoA42EBjfQlw7STb3j4/Q2ZeBcQM1dSWiBgA/h14MnAV8MjMvKXaqiRJktQvDBAHu6Gx/kRmvrfSSqagER4+CpwHXAM8PDNvqrYqSZIk9RMHUR/sS43171ZaxRRExBzgI8DTKFpHHpaZN0z8LEmSJKk9BoiDXQxcDpwVEe+PiENHbxARR0bEH856ZRNohIcPAc8ArqcID9dN/CxJkiSpfZGZVdcwYyLi3sD7mu46DVgK/Bq4rbwzM+/f9JyjgS8Cdwe2AT+h6Nq0ELhL4zVuycwjZ6jm3wZe1XTXvSgGaf8cKAdCr8/MJzU950XAOxs3v0nT+IxRvp2Z/9rRgiVJklQr/T4GYhlwvzHuP3mM+wDIzN9ExH2BP6AYS3D3xmtsAm4E3sbMTsh2GGPXfLemr0cHhOaWkrMneX0DhCRJkqasr1sgJEmSJHWWYyAkSZIktcwAIUmSJKllfTcGYtWqVXncccdV9v2Hh4cBmDPHbKaC+4TG4n6h0dwnNJr7hEab7X3i8ssv35iZh42+v+8CxHHHHcfatWsr+/7btm0DYOnSpZXVoO7iPqGxuF9oNPcJjeY+odFme5+IiDGv7GmklSRJktQyA4QkSZKklhkgJEmSJLXMACFJkiSpZQYISZIkSS0zQEiSJElqmQFCkiRJUssMEJIkSZJaZoCQJEmS1DIDhCRJkqSWGSAkSZIktcwAIUmSJKllBghJkiRJLTNASJIkSWqZAUKSJElSywwQkiRJklpmgJAkSZLUMgOEJEmSpJYZICRJkiS1zAAhSZIkqWUGCEmSJEktM0BIkiRJatncqguQRsuE3bth2zbYurVYlwvA3LnFMm/eyNdz5sD27cX25XPK9Y4dsHcv7Nlz8HrvXliwAJYsgaVLR5YlS+CQQ0a2LZfdu0fWu3fDrl3FUn69e3dR05IlsHjxyDJv3gIWLUqWLYP58++4DAxAxNi/i4GB4jWbl/I5Bw4cvOzfX6znzi1qGL0sXlz8nlTsY+P9zvvBgQNw662wfn2xbNsGK1bAypUjy5Ild/wdDA+P7Ne7dhWPDwwU+9TAwMjXEbBvX7HP7dt38NcDA8Xfz6JFxTJ/fud+18PDxffYu7dYz58PCxcWNWl6MmFoCK6/vvjsLD+/Djlk5OupvpeZI/tH+dlbLpljPydi7M/6efOK/XvbtuIzv/n/w/btxT5S7qvlMmdO8Xo7dxb/D0Yve/cWrzk8fMdl3rxiHxu9LF4MRx8Nd74zHHts8b+jE/bsGfk/tnUr3HzzAJs3B9u2FX/Tt94Kt9xSrDduLGqcP3/kf0P5f2Lu3Dv+vsvbc+aM/E8Y/T9ivL+lgQFYtQoOP7xYDjusWK9aVXxf1Y8fu5pVe/fCz34GP/1pcWBzyy3FcvPNI1/fdlvxz6YT5swp/gEuWDCyzJ8/st69++B/Qnv3jv9a5fMWLBg5OFq4cOTrQw8tPqC3bIGbbhr557Rz5zx27Bj/H+VsWr0ajj8eTjihWJfL0UePBLGIYl1+vW/fyO9o9HrnzuJAc+fOg7/evbv4x5ZZLOXXw8MH11MejEQUy8KFdwxgS5YU90eMvB6MfH3gwEjIaw6Iu3cX/4CHhopl8+aRr3fuHHkfy/ewXA8M3LHm8nuVB9Tl72fOnJGDlOb9o/nrclm48OCvM+cxZ07xPcufv/l3cuDAyIF587Jnz8EH+eWyc2fx97N+ffH3NPp3Pdq8eUWQGBgYee927+7Mfjb6PS7/Rsrf8+hl7tyRYN687Np18Hs63ufCwMAdX68M1eVS3l60aCREHXroyNerVhV/H0cdNbKsXDm1A+b9+4u/kdHvT/Pt5p+xXG/ZMp/9+4v9qvmkQLke6wREuZT73eivm/8+ypMge/YUn0033lgEhnLZvn3in6v8jCj/HpqXsd730mT7YpXKv8vy77r577v8/Cvfqz17xn+dFSuKIHHsscW+1Pz50Px6u3bd8URX8+07/g865KBby5YVB++HHVZ8bg8MHByqd+wY+dxo3j+WLBkJGQcOFNtt3158VmzfPnL7wIGxf75yHxzL0qWwfHmxLFs28vXy5cVjy5YdvF66tKj/HvcwfPSyyG44qumgNWvW5Nq1ayv7/tsap8mXdup0RA8bHoZf/Qouuwx+8INi/ZOfHPwhvHRpcRbjiCNG1oceOvYHztKlxQdweSarPDDYt6/4XkuW3PF5hxzS3gHA3r3Fh+jOnXc8IJzq2ftyn1i0aOkdzr5NdFBUHhyXZ3ibl/37Dz7D1nyGeP/+kX8GzcvWrXDDDXDNNcVyww2d+8fefNa5DGxly0pzKGk+UG4OAlDUUga6MnxNJUg2H0QtXw6DgyPLihXF+pBDRg6oRrckHThwx5rL9eizk+VZy/37R97P0WGm+eBtvH/A7SpDSPOB+aJFxd/Q6tUjB8Ll10uXFgFq06Y7Lpkj71vze1iEnDu2cJU/c/OZzuavDxy448FyeQBdhp/RIaE8sC8DQPl1+XOOdcBcnmFtPrgrv963b+SMdbmUfx87d479e9ix446/53nzit9fGbRGB+w5c4rv2/x3tmPHxAeak5k/PxkYiNtrbm5lHX1Gebr70+GHF2fQy+WYY4r18uXjn60fHh75Ox69lEYfVkSM/R7Omzf+52r59zX6s74MWGVr8ejW4zJ8Ne+rBw4UNZX7eHPrysBA67+v4eGRkxPbthWfodddV4Sv664bWTZvHrs1Y3i4+J7l/6hyKesvD8Cb/4/NnbuT5cuT449fzKpVxd9DFcoWqrL1ozzhV57027q1OHk2etm2rfgbH8uiRfCAB8BZZ8FDHwr3u19xnyY228eZEXF5Zq65w/0GiM4yQBQfMu97H7zrXcXZDSg+rO9zH7jvfeHMM+He94Y73akeHxbduk/s21f8A7z66qLFpPwn23zmfXi4OHhp/gddrsuWgkWLigOBmegSVAa6sjvN6LP0ZYtAGfRmqo5OKVsQ9uyBjRu3MzwMixcvuUOrChx88FsenJctHXZF67zdu2HDhuJvoXlZv34kaJV/E81fL1hwcFeQ8uvyAHV0yCvD2egWmUWLYM+ebUS0/llx4MBIUN23746tcHv2jLTOjG4FK/9u1d269f9HO8oWueYuyb/5DXz723DJJcWJxczis+2+94UnPhH++I+LvyPdkQFihhggqnPjjfCOd8A//3PxYXHOOfDUpxYfCKee2t6Znn5S531C43O/0GjuExqtDvvE0FARJi69FL7xDVi7tugG9hd/AS98YefGl/SLbgkQnsfStF15JfzBHxR96d/xDnjc4+CHP4Qvfxme+1w4/fT6hgdJkjS+wcHiuOHNby66On/ve7BmDfz1X8Nxx8Eb31i0Xqi7GCA0ZZnwspfBaafBRRfBH/4h/PrX8B//Afe6V9XVSZKkXnP/+8OXvgTf/37x9f/7f0WQeP3ri7E56g4GCE1JJrz0pfC2t8HznlcMHHvve4ur+0iSJE3H/e4HX/hCcRGWBz0IXvUqeOQji0Hbqp4BQm3LhL/6K/jHf4QXvagY83D44VVXJUmS+s2ZZ8LnPgef+hRcfnlx1aabbqq6Khkg1LZXv7roq/j85xchopuveiNJknrfk58MX/wiXHstPPjBsG5d1RXVmwFCbXnd64p+iM97HrznPYYHSZI0Ox7xCPj614tB1Q9+cDEpraphgFDL3vSmovXh/POLbktei16SJM2mM8+Eb32ruLrjQx8K3/lO1RXVk4eAasnb3lZcUu2Zz4QPfMDwIEmSqnHXuxbB4Ygj4FGPKq7apNnlYaAm9ZnPFJdrPe88uPBC53SQJEnVOvbYoiXirneFc8+FX/yi6orqxQChCe3cWVxp6e53h49+FObOrboiSZKk4gqQX/4yLF4Mf/7nxVUiNTsMEJrQ3/89XH99MWB63ryqq5EkSRpx2GHwmtfAf/83fP7zVVdTHwYIjeuqq4rLtT7zmcVAJUmSpG7zghfAqacWE9zu2VN1NfVggNCYMouuSwsWwFveUnU1kiRJY5s3r5iX6qqr4J3vrLqaejBAaEyf/WxxVYPXvAZWr666GkmSpPGdcw48/vHFfFUbNlRdTf8zQOgOdu2Cl7wETj8d/uzPqq5GkiRpcm97W9GF6ZWvrLqS/meA0B286U3FVPEOnJYkSb3i5JOLE6Af+hBcdlnV1fQ3A4QOsm4d/MM/wNOfDmefXXU1kiRJrfubvykmmHvRi7ys60wyQOggL3lJ0erw1rdWXYkkSVJ7li0rLkH//e/DRRdVXU3/MkDodp//fLH87d/CUUdVXY0kSVL7zj8f1qyBl78ctm+vupr+ZIAQAMPDxSyOd70rvPjFVVcjSZI0NXPmwLveBTfdVHTLVucZIATA975XXD/5la904LQkSeptD3gAPPGJcMEFsH9/1dX0HwOEAPjYx4pJ4849t+pKJEmSpu9Zz4JbboFLL626kv5jgBAHDsAnPwmPfSwsXVp1NZIkSdP3mMfA4sXFSVJ1lgFCfPvbsH49PPWpVVciSZLUGYccUsxO/alP2Y2p0wwQ4mMfg0WL4HGPq7oSSZKkzjnvPNi0Cb7xjaor6S8GiJrbv79I5o97XNHMJ0mS1C8e/WhYsgQ+/vGqK+kvBoiau+SSYoCR3ZckSVK/WbQInvAE+K//gn37qq6mfxggau5jHytaHh772KorkSRJ6rzzzoPbboOvf73qSvqHAaLG9u0rEvkTnlAkdEmSpH5zzjmwbJndmDrJAFFjX/96MbDI7kuSJKlfLVxYzHP16U/D3r1VV9MfDBA19vGPF4n8nHOqrkSSJGnmnHcebN4MX/ta1ZX0BwNETe3dW3RfOvfcIplLkiT1q0c9CpYvtxtTpxggaup//geGhuy+JEmS+t+CBfDEJ9qNqVMMEDX18Y/D4GCRyCVJkvrdeefBli3w3/9ddSW9zwBRQ7t3w8UXw5OeBPPnV12NJEnSzHvkI4uTp3Zjmj4DRA195SuwdavdlyRJUn3Mn1+cPP3MZ4qTqZo6A0QNffzjsHIlPPzhVVciSZI0e847rziJajem6TFA1MyuXfDZz8KTnwzz5lVdjSRJ0ux5xCNgxQq7MU2XAaJmvvhF2L7d7kuSJKl+5s0rTqJ+5jPFSVVNjQGiZj73OVi1Cs46q+pKJEmSZt955xUnU+3GNHUGiJr5wQ/gAQ+AuXOrrkSSJGn2nX12MS/Et75VdSW9ywBRI9u2wZVXwplnVl2JJElSNebPhzPOgLVrq66kdxkgauSHP4RMWLOm6kokSZKqs2YNXH45HDhQdSW9yQBRI2XSNkBIkqQ6O/PMYhzE//1f1ZX0JgNEjVx2GRx7LBx2WNWVSJIkVac8mXrZZdXW0asMEDWydq3jHyRJkk49FRYvNkBMlQGiJm67Ddats/uSJEnSwADc+94OpJ4qA0RNXH55sTZASJIkFb0yfvxj2Lev6kp6jwGiJsomuvvcp9o6JEmSusGaNbB7N1xxRdWV9B4DRE2sXQsnnwyDg1VXIkmSVL1yXKjjINpngKiJyy5zALUkSVLpxBOLE6uOg2ifAaIGNmyA3/zG8Q+SJEmliOLYyBaI9hkgaqBM1rZASJIkjTjzTPjZz4qxEGqdAaIG1q6FOXPgXvequhJJkqTusWYN7N8PP/lJ1ZX0FgNEDVx2GZx2WjFhiiRJkgoOpJ4aA0SfyyxaIBz/IEmSdLCjj4bDD3cgdbsMEH3uhhvgllsMEJIkSaNFFK0QtkC0xwDR5xxALUmSNL4zz4Rf/hK2b6+6kt5hgOhzl10Gc+fCPe5RdSWSJEndZ82aosv3D39YdSW9wwDR59auLcLDwoVVVyJJktR9ym7ejoNonQGijzmAWpIkaWJHHAHHHOM4iHZ0bYCIiDdGRDaWl1VdTy9atw6Ghhz/IEmSNBEHUrenKwNERJwJvBzIqmvpZWVTnC0QkiRJ41uzpjjxunlz1ZX0hq4LEBGxAPgwcDPwmYrL6WmXXVaMfTj99KorkSRJ6l5lbw3HQbSm6wIE8HfAXYE/AbZUXEtPW7sWzjgD5s2ruhJJkqTudZ/7FGsDRGumHCAi4qSIeEBE3KVTxUTE/YC/AC7KzM916nXr6MABuPxyuy9JkiRNZsUKOOkkx0G0qq0AERFzI+LVEXEz8Cvg28BfNT3+zIj4bkTcrd1CImIhRdel24AXt/t8HexXv4IdOxxALUmS1Iozz7QFolUtB4iImAt8EfhbYBD4JRCjNvsOcH/gd6ZQyxuAU4A/y8yNU3i+mpQJ2hYISZKkya1ZAzfcADffXHUl3W9uG9v+KfBI4KvA+Zm5PiKGmzfIzGsj4irgt4DXtvrCEfFA4CXAxZn5sTZquoPh4WG2bds2nZeYliq/d7PvfncBS5bM46ijttMlJdVWt+wT6i7uFxrNfUKjuU/MrtNPHwAO4dJLd/LoRx+oupwxdcs+0U4Xpt8DNgHnZeb6Cbb7JXBMqy8aEYuAC4GtwAvaqKf5Nf4oItZGxNqNG228APjhDwe45z0PMDBQdSWSJEnd7x73OMCcOckPf+jB02TaaYE4BfhmZg5Nst024LA2XveNwMnA708STMaVmRcAFwCsWbMmly5dOpWX6agqa9i3D372M/jTP622Dh3M90Jjcb/QaO4TGs19YnYsXQp3vSv89KcLWLp0QdXlTKjqfaKdAJHA8KRbwVHA7jZe90mN1z0/Is4f9dipjfXzI+JxwFWZ+bw2XruWrrgC9uxx/IMkSVI71qyBL30JMiFGj/TV7doJENcA94yIOZk5ZpBodEe6B0U3pnbMAc6a4PETGstgm69bS7/4RbG++92rrUOSJKmX3Oc+8OEPw4YNsHp11dV0r3bGQHwWOJpinobxvBxYQRszSGfmcZkZYy0Ul3UF+MvGfWe0UW9trVtXrE84odo6JEmSeslJJxXr8lhKY2snQLwd2AC8KSIuiognN+5fFRGPiYgPAq8Grgfe1+E61Yarr4ajjoJFi6quRJIkqXeceGKxvvrqauvodi13YcrM2yLi0RStC08DnkoxLuK3G0sANwCPz8zuuMZUTa1bN/IHIEmSpNYcdxzMmWMLxGTaGQNBZv4sIk4Dngs8hmJcwgBFcPgScEFm7uh4lWrLunXwqEdVXYUkSVJvmT8fjjnGADGZtgIEQGbuBv6pscyozHwO8JyZ/j79ZNcuuOkmWyAkSZKm4oQTDBCTaWcMhHrANdcUawOEJElS+0480TEQkzFA9JkyMRsgJEmS2nfiiXDLLbDNEb3jGrcLU0RMJ3tlZnoIWwEv4SpJkjR15THU1VfDPe9ZbS3daqIxEMeNc39SXHFposdyGjVpGtatK6ZiX7Wq6kokSZJ6T/OlXA0QY5uoC9PxYyzvoAgH/wU8CbhXY3kS8KnGY2+nuDqTKnD11cWO7/TrkiRJ7SsDhAOpxzduC0RmXtd8OyKeCLwEeFpmfmLU5j8BPhMRTwE+BnwHuA7NunXr4G53q7oKSZKk3jQ4CCtWGCAm0s4g6pcBPxgjPNwuMz8J/KCxrWbZgQPFVZgc/yBJkjR1XolpYu0EiHsArWSxqwHPgVfgpptg716vwCRJkjQdJ55oC8RE2gkQCZzawnanTLEWTZOXcJUkSZq+E0+E666D/furrqQ7tRMgfgDcKyL+cLwNIuJ5wL2B/51uYWqfl3CVJEmavhNOKMLDDTdUXUl3mugyrqP9HfAw4P0R8TTg34HGvMccBzyz8fgB4PUdrFEtWrcO5s6FO9+56kokSZJ6V/OVmI4/vtpaulHLASIzvxURvwf8M0VQOHvUJgHsAP4kMy/tWIVq2dVXw7HHFiFCkiRJU9McIB75yGpr6UZtHWpm5n9ExDeB5wEPBY5uPHQjcAnwgcy8qaMVqmXr1jn+QZIkabqOOgrmz3cg9XjaPledmeuB181ALZqmdetgzZqqq5AkSeptAwNF1yUv5Tq2dgZRq4tt3lwstkBIkiRNn5dyHZ8Bok+UCdkAIUmSNH1lgMisupLu03IXpohopxEnM9ND2VnkHBCSJEmdc8IJsG0bbNoEq1ZVXU13aWcMxHEtbJMUV2Myq82yMkB4qTFJkqTpa74SkwHiYO10YTp+nOVE4OHA24B9FHNAOJXZLFu3Dg4/HJYurboSSZKk3tccIHSwduaBuG6Ch68BvhkR3wE+CVwKTLS9Ouzqq+2+JEmS1Cllrw6vxHRHHR1EnZkXAz8DXtnJ19Xk1q0r+upJkiRp+hYtKuaDsAXijmbiKky/Bu49A6+rcezZAzfcYAuEJElSJ3kp17HNRIA4gSlMUKepu+664hJjBghJkqTOOeEEuzCNpWMBIiIGIuLlFK0PP+nU62pyXsJVkiSp8048EW68EXbtqrqS7tLOPBBfn+DhJRRXYxoEhoG/n15ZakcZIBwDIUmS1DnlydlrroHTTqu2lm7STlejs1vYZh3w15n5+amVo6lYtw4OOQSOPLLqSiRJkvpH86VcDRAj2gkQD5vgsb3AjZl5/TTr0RRcfXXR+hBRdSWSJEn9o+zd4TiIg7UzD8QlM1mIpm7dOjjppKqrkCRJ6i+rVhWT9HolpoO1PIg6Ip4dEQ9sYbv7R8Szp1eWWpU50gIhSZKkzonwUq5jaecqTBcCz2thuz8APjSlatS29euLKwN4BSZJkqTO81KudzQT80DYE38WlTu0AUKSJKnzTjyxuArT8HDVlXSPmQgQRwPbZ+B1NQYv4SpJkjRzTjwR9uwp5oNQYcJB1GOMZThpgvENc4G7Ao8ALutAbWrBunUwZw4cd1zVlUiSJPWfspfH1VfDMcdUW0u3mOwqTBcC2XT7QY1lPEExkdxbp1eWWrVuXbEzz59fdSWSJEn9p+zlsW4dnHVWtbV0i8kCxEcYCRDnU0wU951xtt0L3Ah8JjN/0pnyNJmrr3b8gyRJ0ky5851h7lyvxNRswgCRmc8pv46I84FvZ+bvz3RRat26dXDuuVVXIUmS1J/mzoVjj/VKTM3amYn6eBwc3VW2bYNbb7UFQpIkaSadcIItEM1avgpTZl6XmZtmshi1x0u4SpIkzTwnkzvYuC0QEXHnxpc3ZuaBptstyczrp1WZJuUlXCVJkmbeiSfCbbfB0BAMDlZdTfUm6sJ0LcUVlU4D/q9xOyfYvllO8trqgDJA2AIhSZI0c8qTtVdfDfe+d7W1dIOJDvKvpwgC+0bdVpdYtw4OPdQkLEmSNJPKk7Xr1hkgYIIAkZnHTXRb1fMSrpIkSTOveS4ItTGIWt1n3TrHP0iSJM20pUvh8MO9lGvJANGjMuGGG+C446quRJIkqf8df7wBojSlgc4RMQCsBBaOt41XYZpZmzbBvn2wenXVlUiSJPW/O90Jrryy6iq6Q1sBIiLuB/wd8BBgwQSbehWmGbZ+fbE2QEiSJM281avhG9+ouoru0PJBfkQ8CPgqI8FhM7B1JorS5MoAceSR1dYhSZJUB6tXw+bNsGsXLFpUdTXVaqeV4LUU4eFfgFdl5i0zU5JasWFDsbYFQpIkaeaVx1wbNhTjIeqsnUHU9wV+mZl/bHionl2YJEmSZs9RRxXr8hisztoJEAH8dKYKUXvWr4clS4pFkiRJM6s8aWuAaC9A/Aywx32XWL/e1gdJkqTZYoAY0U6AeCfwkIg4Y4ZqURsMEJIkSbNn1SqYO9cAAW0EiMz8GPAG4H8i4vkRceeZK0uTMUBIkiTNnjlz4IgjDBDQ3mVcDzTdfA/wnogYb/PMTOeBmEEbNngJV0mSpNm0erUBAtq7jOu4aWGa26pN27cXiy0QkiRJs2f1arj++qqrqF47XZjmtLPMZNF15yVcJUmSZp8tEAUP9HuQAUKSJGn2rV4Nt9wC+/ZVXUm1DBA9yAAhSZI0+8pjr5tvrraOqhkgepABQpIkafY5G3Vhqldhmsg+YCOwFrgwMy+eQl2awIYNMG8eHHpo1ZVIkiTVh5PJFdppgYgWl/nAUcATgE9FxIc6WbCKnfbII2H8q+hKkiSp0wwQhbauwgS8FdgOvBk4A1gBLAfuCfwDsA14O3Bn4HzgVuDZEfGMjlZdc04iJ0mSNPuOOKI4gVv3ANFOF6bnAi8BHpqZ3x/18M+Av46Ii4FvAb/MzA9ExK+B7wLPAS7qRMEqdtoTT6y6CkmSpHqZOxcOO8wA0U4XphcC3xojPNwuM/+XIkC8oHH7+8CPgHtNp0gdzBYISZKkajgXRHsB4lRgQwvbbQBOabp9NbCsnaI0vr17YdMmA4QkSVIVDBDtBYg9FOMeJnNGY9vSfIqxEeqA8rrDBghJkqTZt3o13HRT1VVUq50A8W3g1Ih49XgbRMTfAHel6MZUOh6oeU7rHOeAkCRJqs7q1cUJ3QOtTnDQh1oeRA28GngU8LcR8XTgY8B1QALHAudRdHPaDbwGICLuDNwNeF/nSq63MkAceWS1dUiSJNXRUUcV4WHjxuKqTHXUcoDIzJ9ExOOAj1KMcXjVqE0CuBn4vcz8ceO+XRSh48rplyqwBUKSJKlKzXNBGCBakJlfj4gTgacAZwF3ajx0E3Ap8InM3Nm0/a3A1zpUqyh21oj67rCSJElVag4QZ5xRaSmVaStAAGTmLuDfGotm2fr1xfWH57b9zkmSJGm6nI26vUHU6gLOASFJklSdchyqAUI9Y8MGA4QkSVJVFi6EFSsMEC2LiPkR8ZcR8b8RsTkiDoyz7J+pguvOFghJkqRq1X0yuZZ70kfEQuAbwH0prrg04ebTKUpjGx4urjvsJVwlSZKqU/cA0U4LxEuB+wFfBu4CfIRiDogFwOnA31PMAfGGzLRr1AzYuBH277cFQpIkqUp1n426nWv5PAXYCjw9M7dGRAJk5j7gl8D/i4hvAV+IiCsy8z87X269OQeEJElS9coWiMzi8vp1005LwcnA/2bm1sbtBIiIgXKDzPwycBnwpx2rULczQEiSJFXvqKNg717YvLnqSqrRToCYA2xqur2rsR4ctd064G7TqEnj2LChWBsgJEmSqlP3uSDaCRA3AUc13f5NY32PUdsdR6N1Qp1lC4QkSVL1DBCt+zlwStPtSymutvSaiFgKEBFPBx4A/KJjFep269fD8uWwaFHVlUiSJNWXAaJ1XwKOiIizATLzO8D3gIcAmyJiE/BRitaHt3a2TEGxk3oJV0mSpGoZIFp3EUVY+HXTfU8CvtB4nRXAEPCXmfnpThWoEU4iJ0mSVL0lS4qlrgGi5cu4ZuZ24Duj7rsFeHxEHAIsB27OzOHOlqjS+vVwv/tVXYUkSZLqPJlcO/NAjCszdwI7O/FaGltmcRUmWyAkSZKqV+fJ5Jwxukds2wY7dxogJEmSuoEtEGOIiGdP54Uz8yPTeb4O5iVcJUmSusdRR8HnP191FdWYqAvThUxvPgcDRAeVAcKrMEmSJFVv9WrYsaPoJbJ0adXVzK5WxkD8CKhpA033sAVCkiSpezRfytUAMWIXsIhipukNwIeAz2bmvtkoTAczQEiSJHWP5gBxl7tUW8tsm2gQ9ZHAHwOXAY8FPg7cFBHvjIgzZqE2NdmwARYsgMHBqiuRJElSnSeTGzdAZOa2zPyXzHwQcArwZmA38GfA5RHx44h4cUSsmqVaa62cRC6i6kokSZJkgJhEZv46M/8auDNFa8QnKULFO4AbI+K/IuK3Z65MOQu1JElS9xgcLHqHGCAmkYUvZ+ZTgdXAC4GfAk8EPtj58lQyQEiSJHWPiPrOBTGdieSOAI6lCBIAdq6ZQevXewlXSZKkblLX2ahbuYzr7SJiGfB04DnAfSlCwybg3dgCMWN274bNm22BkCRJ6iarV8Mvf1l1FbNv0gAREQE8iiI0PBFYAAwDX6S4tOvnvLTrzLr55mJtgJAkSeoeRx0FX/961VXMvnEDRETchSI0/B5wFEVrwy8oZqj+t8y8eRbqE84BIUmS1I1Wr4ahIdi1CxYtqrqa2TNRC8SVQAJDwPuBD2Xm2tkoSgczQEiSJHWf8thswwY4/vhqa5lNrYyB2AQ8BHhItD4JQWbmPadclQ5igJAkSeo+zXNBGCBGBHDSFF43p/AcjWP9epgzBw47rOpKJEmSVKrrZHITBYiHzVoVmtD69XD44TAwUHUlkiRJKhkgRsnMS2azEI1vwwa7L0mSJHWbVatg7tz6BYjpTCTXURHxZxHx8Yj4ZURsioh9EXFrRHw1Ip4VbQzA6DfOQi1JktR95syBI46o32RyXRMggFdQzDOxC/gu8CngKuDhwL8Bn46Ibqp31hggJEmSutPq1fVrgWhrJuoZ9jTgR5m5o/nOiDgd+BpwLnA+xeR1tXHgQDGRnAFCkiSp+xx1FFx7bdVVzK6uOaOfmd8eHR4a918BvLdx81GzW1X1br0VhocNEJIkSd2oji0QXRMgJrG/sd5TaRUVcA4ISZKk7rV6dXHCd9++qiuZPV0fICLieOBPGjc/W2UtVSgDxJFHVluHJEmS7qg8yXvzzdXWMZu6aQwEABHxXOAsYB5wNPBAiqDzxsz8dJW1VWHDhmJtC4QkSVL3aZ4L4uijq61ltnRdgAAeRDFYurQfeBXw9laePDw8zLZt22airpZ0+ntfe+18YAGLF2+jwh9L01Dl/qju5X6h0dwnNJr7RG9YtmwOsJirr97JqacemNHv1S37RMtdmCLi+435GObPZEGZ+bzMDOAQ4HTgH4HXAN+PiKPGqe2PImJtRKzduHHjTJY3626+ORgcTBYurLoSSZIkjXbkkQnAhg1dPzKgY9ppgbgv8GHg7RHxAeD9mXndzJQFmbkL+AXwlxGxAXgr8B7gyWNsewFwAcCaNWty6dKlM1VWyzpVw5YtcPjhnXs9Vcf3UGNxv9Bo7hMazX2iux17bLHesWMhS5fOzhnfqveJdqLSucD/ACspJn27KiI+ExHnzEhlB7uwsX58RMybhe/XNTZuLKZJlyRJUvdZuBCWLCmO2eqi5QCRmZ/LzEcDdwHeAWwBHg98MSL+LyL+PCIGZ6ZMNlOMhZgLHDpD36MrGSAkSZK628qVsGlT1VXMnrY7a2Xmusz8C4orJP0B8EPgJIouRjdGxL9GxL07WyYPpQgPQ0CN8l2xM65cWXUVkiRJGs+qVbZAtCQzd2fmhzLzTOB+wEXAIuC5wGUR8d2IeEorrxURD46Ix0XEHcZkRMSDgA80bn4gM2d2eHuXsQVCkiSpu9WtBWLal3GNiCOAcyjmbgAIYCtwf+BjEfF94MmZOdH0GicBHwKGIuKHwAZgKXAicFpjmy9QXM61NnbuhN27bYGQJEnqZqtWwbp1VVcxe6bcAhERD42I/wSuA14LHAH8J8XEb4cCTwJ+BDyAYszERC4BXgf8GDiZ4kpLvwUsBj4FPCkzH9e4MlNtlE1htkBIkiR1r7p1YWqrBSIilgDPBp5P0TIQwM3AP1Nc1nVD0+afiYjPAz+hCAPjysxrgFe3U0sdGCAkSZK638qVxaX39+2DeTW4XmjLASIi3gc8E1hCERx+ALwL+ERm7hvrOZl5ICJ+wMEzS6tFZV86uzBJkiR1r/Jk7223wRFHVFvLbGinBeJPgL3AvwPvzszLWnzepRSBQ22yBUKSJKn7lcdqGzcaIEZ7NXBBZt7SzjfIzAsZmQhObTBASJIkdb+yt0hdrsTUziDqj1BM5jahiFgREXeeekkqbdoEEbBiRdWVSJIkaTzNLRB10E6AuAZ4SwvbvRm4emrlqNnGjUV4GBiouhJJkiSNxxaI8QWtj2VwzEMHOAu1JElS9yuP12yBmLpBYM8MvG7tOAu1JElS91u0CBYvrk+AmHAQ9RhjGZZMML5hLnBXijkfrulAbbW3cSMcc0zVVUiSJGkyK1fWpwvTZFdhuhbIptu/01gmEhSXetU0bdoEZ5xRdRWSJEmaTJ1mo54sQFzPSIC4M7ATGO9Xsxe4Efg08J6OVFdzdmGSJEnqDStXGiAAyMzjyq8jYphi1unfn+miBDt3wq5dBghJkqResGoVXFOTTvztTCT3XOCqmSpEByv70HkVJkmSpO5nF6YxZOaHZ7IQHcxZqCVJknrHypUwNAT798Pcdk7R96CZuIyrOqAMELZASJIkdb/ypO9tt1Vbx2wYNx9FxNUUA6gfmZnXNG63KjPzxGlXV2NlFyZbICRJkrpf82Ryhx9ebS0zbaIGluMoAsS8ptutysk30UTswiRJktQ7ymO2OswFMVGAOL6xvnHUbc2CcudbsaLaOiRJkjS5MkDUYSD1uAEiM6+b6LZm1saNRXjo90E4kiRJ/aDswlSHFggHUXcpJ5GTJEnqHc1jIPpdywEiIo6JiGdHxCkTbHNKY5ujO1NefW3a5BWYJEmSesUhhxSLAeJgLwI+1MJ2FwIvnFI1up0tEJIkSb1l5Uq7MI32W8AVmfmr8TZoPHYFcM50C6u7jRttgZAkSeoldZmNup0AcQxwVQvbXQXceWrlqLRpky0QkiRJvWTlSgPEaAuBvS1stxdYPLVyBLBzZ7EYICRJknrHqlV2YRrtRuA+LWx3b2DD1MoRjOx4dmGSJEnqHXZhuqNvACdExHPG2yAizgdOBL4+zbpqrQwQtkBIkiT1jpUrYWgI9u+vupKZ1U6AeDuwD7ggIt4QESeUD0TE8RHxBuCCxjZv72yZ9VImVwOEJElS71i1CjJh8+aqK5lZLQeIzLwS+KPGzb8Cfh0ReyJiD8XA6b9qvN4fZ+YVHa+0RuzCJEmS1HvqMplcWzNRZ+ZHgAcCnwN2AfMay67GfQ/MzAs7XGPt2AIhSZLUe8pjt34fSD233Sdk5lrgiRExB1gFJLApM4c7XVxdlQHi0EOrrUOSJEmtKwNEv7dAtB0gSo3AcEsHa1HDpk0wOAhzp/zuSJIkabbVpQvTlA5RI2I+xSVd79S460bg8sxsZZ4ITWLjRrsvSZIk9Rq7MI0hIuYBrwFeCCwd9fD2iHg38NrM3NeZ8upp40YHUEuSJPWaQw6BRYtsgbhdRAwAnwceCQSwHri68fAJwGrgr4EzI+KxmXmgw7XWxqZNsHp11VVIkiSpXStX9n8LRDtXYfoj4FHAr4HHZOadMvMhjeVOwGOA/6MIGH/Y+VLrwy5MkiRJvakOs1G3EyCeDewAHpGZXxn9YOO+RwI7gfM7U149bdpkFyZJkqRetHKlAaLZacA3MvPG8TZoPPaNxraagl27YMcOWyAkSZJ60apVdmFqNo+idWEyOxvbagrKHc4AIUmS1HvswnSw64CHNC7hOqbGYw9pbKspKAOEXZgkSZJ6z8qVsHkzHOjjywm1EyA+S3GlpQ9HxODoByNiOfBB4EjgMx2probKxGoLhCRJUu9ZtQoyixDRr9qZB+LNwNOB84DHRMTngGuApLiM6+Mp5ob4TWNbTUEZIGyBkCRJ6j3Ns1H36wnhlgNEZm6KiIcDFwFrgGdShAco5oUAuAx4Rmbe1tEqa8QxEJIkSb2rDrNRtzUTdWZeBdw3Ih4MnAXcqfHQjcAlmfntDtdXO2ULxKGHVluHJEmS2lcGiH4eSN1WgCg1goJhYQZs3AjLl8M8r2MlSZLUc5q7MPWrdgZRaxZs2mT3JUmSpF5Vhy5MBogu088DbiRJkvrdIYfAwoX93QIxbhemiJjO1WszM6fUParuNm2CI46ougpJkiRNRUTRjamuLRAxjcWWjSmyBUKSJKm39fts1OO2EmSmIaACGzc6B4QkSVIvW7myvwOEIaGL7N4NO3bYAiFJktTLVq2qbxcmzTInkZMkSep9te3CNJ6IWA48C3gAcBjwtcx8c+OxuwDHAd/KzF0drLMWyh3NLkySJEm9a+VKuO02OHAABgaqrqbz2goQEfFo4N+BQYrB0kkxC3XpFOBi4BnAxzpSYY3YAiFJktT7Vq2CTBga6s8Twy13YYqIuwH/BSwF3gc8lSJENPsysBM4t1MF1knZAmGAkCRJ6l39Pht1Oy0QrwQWAE/KzM8CRMRBrQyZuS8ifgTcs3Ml1kfZAtGPSVWSJKku+n026nYGUZ8N/KgMDxO4EVg95YpqzDEQkiRJva8MEP3aAtFOgFgJXNXCdvOBRVMrp942boRly2DevKorkSRJ0lT1exemdgLEZuDoFrY7Ebh5auXU26ZNjn+QJEnqdXZhGvED4MyIOHm8DSLiTOAewHemW1gdbdxogJAkSep1ixfD/Pm2QAC8F5gHfDIiThn9YEScAHyQ4tKu/9SZ8upl40bHP0iSJPW6iP6ejbrlAJGZXwHeDdwd+EVE/JQiLDwyIv4XuBI4HXhHZn57Jortd3ZhkiRJ6g/9PBt1Oy0QZOaLgRdQjHG4G8U8EEcDZwJbgJdk5ss6XWRd2IVJkiSpP6xc2b8BYtx5ICLicuADwEWZOVTen5nvj4gLgDOAE4AB4AbgB5m5f0ar7WN79sD27XZhkiRJ6gerVsHPflZ1FTNjoonk7kXRZemtEXEx8KHM/B+AzBwGfthY1AFlHzlbICRJknpfXbswvRz4JbAQeBrw5Yi4NiJeExHHzUZxdeIkcpIkSf1j5Uq47TYYHq66ks4bN0Bk5lsz827A/YB/phjjcGfgVcBVEfHViHhGRCyYnVL7my0QkiRJ/WPVqiI8DA1VXUnnTTqIOjMvy8znA6uBZwJfbTz0cODfgA0R8b7GHBCaorIFwgAhSZLU+/p5Nup2LuO6JzP/IzPPAY6laIlYBywH/gT4fkT8NCJeEhEeBrfJLkySJEn9o59no27rMq6lzLwxM9+QmXcBzgI+BOyguLTr2yiuyqQ2lDuXAUKSJKn3lQGi1i0Q48nMb2XmHwDPBm6lmBti/nRft242boRly4ppzyVJktTb+rkL00SXcZ1URJwIPAf4PeAYivAA8N3plVU/Gzfa+iBJktQv+rkLU9sBIiIWA+cBzwUeVN4NbAA+AnwwM/+vYxXWxKZNDqCWJEnqF0uWwLx5NW+BiIizKVobfgc4hCI07AO+AHwQ+FJmHuh4hTWxcaMBQpIkqV9E9O9kchMGiIg4liI0PBs4jpEuSldQDJz+t8y8dQbrq42NG+HUU6uuQpIkSZ2yalXNujBFxNeBh1KEhgC2Av9J0UXpB7NTXn3YhUmSJKm/rFxZvxaIs4EELqHoovTJzNw9G0XVzd69sG2bg6glSZL6yapVcMUVVVfReRMFiNcDH8rMa2armLoqm7ZsgZAkSeoftRsDkZmvns1C6qzcsQwQkiRJ/WPlyuJE8fAwzJn27Gvdo49+lN5VBgi7MEmSJPWPlSuL8LBlS9WVdJYBogts3lysV6yotg5JkiR1TnlsVx7r9QsDRBcYGirWBghJkqT+UR7blcd6/cIA0QVsgZAkSeo/g4PF2hYIddzmzcXAmiVLqq5EkiRJnWIXJs2YoaEiofbT6HxJkqS6swuTZszmzXZfkiRJ6jf92oVpoonkxhURRwNHAQvH2yYzL51qUXWzefPIDiZJkqT+sGQJDAzUPEBExJOBvwdOmmTTbPe168wWCEmSpP4TURzj9VsXppYP8iPi8cDHKbo9bQGuBrbOUF21MjQExxxTdRWSJEnqtBUr6t0C8UoggL8B3pKZ+2ampPqxC5MkSVJ/Ghysd4C4B/CjzHzjTBVTR5l2YZIkSepX/diFqZ2rMO0DfjVThdTV7t2wd68BQpIkqR/1YxemdgLE5cAJM1VIXTkLtSRJUv/qxy5M7QSINwH3jYhHzVQxdVTuUI6BkCRJ6j9lF6bMqivpnHbGQPwKeAPw2Yh4F/AF4HpgeKyNM/P66ZfX/8o+cbZASJIk9Z8VK2DfPti5ExYvrrqazmgnQFxLMb9DAC9rLONxHogW2YVJkiSpfzXPRl3HAHE9RTBQB9mFSZIkqX+VJ4mHhuDooystpWNaDhCZedwM1lFbtkBIkiT1r/IYr58GUrcziFozoBwDYQuEJElS/2nuwtQvDBAV27wZliyBuY4YkSRJ6jvNXZj6hQGiYs5CLUmS1L/6sQvTuOe9I+IAxaDp0zLz/xq3W5WZ6Tn1FgwNGSAkSZL61fLlxboWAYLicq0x6nar2tm21myBkCRJ6l8DA7BsWU0CRGbOmeh2J0XEPOChwGOBs4C7AAuBW4HvAe/JzG/O1Pev0ubNcMIJVVchSZKkmVLORt0vuqWb0VnA/zS+3gBcCuwATgN+B/idiHhdZr66ovpmjF2YJEmS+tvgYH+1QHTLIOph4FPAQzNzdWY+LjOfmpl3B54GHABeFREPq7TKGWAXJkmSpP62YoUBouMy8+uZ+ZTM/NYYj30MuLBx81mzWtgM27cPtm93DghJkqR+1m9dmLoiQLTgR411n0wAXih3JFsgJEmS+pctENU4ubFeX2kVHWaAkCRJ6n+OgZhlEXEk8JzGzU9VWErHlTuSXZgkSZL614oVsHMn7N1bdSWd0S1XYRpTRMwFPgosB76WmZ+b7DnDw8Ns27ZtxmsbTzvf+8YbB4BDWLBgJ9u2tTNPn3pJlfujupf7hUZzn9Bo7hP9Y9GiecBCfvOb7Rx2WE75dbpln+j2Foj3A48AbmCCAdQR8UcRsTYi1m7cuHHWipuuLVuK+fYGB6e+I0mSJKm7LV9eHOv1y0Dqrm2BiIh3An9AMS/EIzJzw3jbZuYFwAUAa9asyaVLl85OkRNopYbdu4v10UcvpgtK1gzrhv1S3cf9QqO5T2g094ned9RRxXrfviUdOearep9oO0BExHKK1oAHAIdRdC16c+OxuwDHAd/KzF1TLSoi3ga8iGIm6kdk5q+n+lrdzDEQkiRJ/a+8YE6/DKRuK0BExKOBfwcGgQASuLFpk1OAi4FnAB+bSkER8WbgpcAm4JGZ+YupvE4vGBqCBQtg0aKqK5EkSdJMKU8W90uAaHkMRETcDfgvYCnwPuCpFCGi2ZeBncC5UykmIt4E/CWwGXhUZv50Kq/TK5yFWpIkqf+Vx3t1HAPxSmAB8KTM/CxARBzUypCZ+yLiR8A92y0kIl4PvAIYoggPP5r4Gb1v82a7L0mSJPW7fmuBaCdAnA38qAwPE7gRuFs7RUTEE4D/17h5FfBnEaMbNwC4MjPf1M5rd7OhIVsgJEmS+t3ChcVSxwCxEri0he3mA+326j+06es1jWUslwB9EyA2b4Yjjqi6CkmSJM20FSv6pwtTO/NAbAaObmG7E4Gb2ykiMy/MzGhhObud1+12dmGSJEmqhxUr+qcFop0A8QPgzIg4ebwNIuJM4B7Ad6ZbWB04iFqSJKkeBgfrGSDeC8wDPhkRp4x+MCJOAD5IcWnXf+pMef1reBi2bDFASJIk1UEtuzBl5leAdwN3B34RET+lCAuPjIj/Ba4ETgfekZnfnoli+8m2bUWIMEBIkiT1v7p2YSIzXwy8gGKMw90o5oE4GjgT2AK8JDNf1uki+5GzUEuSJNVHP3VhamsmaoDMfH9EXACcAZwADAA3AD/IzP2dLa9/lU1YtkBIkiT1vxUriu7rw8Mwp61T+N2n7QABkJnDwA8bi6agTKAGCEmSpP63YgVkwtatvd8DpcfzT++yC5MkSVJ9lCeN+6EbU8sBIiL+NCIORMTjJtjmcY1t/rgz5fUvuzBJkiTVR3nSuFYBAngicAvwhQm2+SJwK/DkadRUC3ZhkiRJqo/ymK8fLuXaToA4Ffh5ZuZ4GzTGRvwMuOt0C+t3mzcXA2iWLKm6EkmSJM20WnZhAg6juHzrZG4BDp9aOfWxeXPRlNXro/AlSZI0ubp2YRoC7tzCdkcD26dUTY0MDdl9SZIkqS7q2oXph8D9I+Lk8TZoPPYA4EfTLazfbd5sgJAkSaqLJUtgYKB+LRAfopg34jMRceroByPiFOBiionlPtSR6vpY2YVJkiRJ/S+if2ajbnkiucz8eEQ8E3g88LOI+B5wZePhU4AHUoSHL2TmRR2vtM8MDcExx1RdhSRJkmbLihX90YWp3ZmonwK8BfgT4MGNpbQPeB/wl50prb/ZhUmSJKleVqyoWQsEQGbuA14SEW8AHg4cCyRwPfD1zLy18yX2n0y7MEmSJNVN7bowNWsEhY91uJba2L0b9u61BUKSJKlOVqyA66+vuorpcxaCCjgLtSRJUv3UsgsTQEQsBNYARwELx9suMz8yjbr6Wrnj2IVJkiSpPsouTJnFVZl6VVsBIiL+HHg1sKyFzQ0Q4yhH39sCIUmSVB8rVsC+fbBrFxxySNXVTF3LASIifh94W+PmLyku4bp1Jorqd3ZhkiRJqp/y2G/z5poECOBFFFdc+j3neZgeA4QkSVL9lN3XN2+GO92p0lKmpZ1B1KcA3zU8TJ9jICRJkuqnPHnc65PJtRMgdlDM96BpKncaA4QkSVJ9NHdh6mXtBIjvAnebqULqZPNmWLoU5k5pFg5JkiT1ojoGiNcCp0bE+TNVTF04C7UkSVL9NI+B6GXtnANfDLwd+GBEPBb4AkWXpuGxNs7MS6dfXn8aGnIAtSRJUt2UAaLXx0C0EyC+SXEVpgCe0ljGk22+dq1s3myAkCRJqpuBAVi2rF4tEJdSBANN0+bNcMIJVVchSZKk2VbORt3LWg4QmXn2DNZRK3ZhkiRJqqcVK3q/C1M7g6jVIXZhkiRJqqcVK3q/BWLKASIi5kfE6og4tJMF9bt9+2D7dgOEJElSHfVDF6a2A0REPDsiLqOYWO43wFubHntSRFwUEcd3sMa+4iRykiRJ9VW7LkwRcSHwIeA+wC6KKzI1+xXwNCa+QlOtlTuMLRCSJEn1U6suTI0J5J4N/ARYAywfvU1m/gK4AXhMpwrsN+UOY4CQJEmqn8FB2LGj6Nbeq9ppgfhDYBvw+Mz8YWaOd0nXnwHHTbewflUGCLswSZIk1U95ErmXuzG1EyDuDnw/M2+cZLsh4MgpV9Tn7MIkSZJUX+UxYC93Y2onQMwDtrew3eFADzfKzCy7MEmSJNVX2QulLgHieuBuE20QEQPA6cC66RTVz+zCJEmSVF9168L0FeCkiHjWBNv8MbAa+MK0qupjQ0OwYAEsWlR1JZIkSZpt/dCFaW4b274FOB/4YEScBnyycf/CiLgr8LvAK4FNwLs7WmUfcRZqSZKk+qpVF6bM/A3wJIpxEK8ALgMSeCrwc+A1wG7gKZl5S8cr7RObN9t9SZIkqa7q1oWJzPwGcBrF7NNXUEwmt5dizMO7gbtl5iWdLrKf2AIhSZJUXwsXFksvt0C03IUpIu4BDGfmzylaIF4xY1X1saEhOOKIqquQJElSVXp9Nup2WiB+DLxnhuqoDVsgJEmS6m1wsD5dmIaA38xQHbXhGAhJkqR6q1sLxIkzVEctDA/Dli22QEiSJNVZnQLEu4D7RcSjZ6qYfrdtWxEiDBCSJEn1NTjY2wGinXkgfkgxBuIzEfFB4NPAdRRXYrqDzLx++uX1F2ehliRJ0ooVvT0Gop0AcU1jHcAfNZbxZJuvXQvljmILhCRJUn2VAWJ4GOa0NalCd2jnIP8GimCgKSpbIAwQkiRJ9TU4CJmwdWtv9kxpOUBk5nEzWEct2IVJkiRJzbNR9+JxYQ82mvQuuzBJkiSpPBbs1YHUBohZZBcmSZIkla0OtQkQEXFWRHwiIn4TEXsi4gNNjz0qIt4YEUd2tsz+sHlzMVBm6dKqK5EkSVJVmrsw9aK2rpQUEa8BXkVxJabb7276egh4BXAj8N5p1tZ3ylmoIybdVJIkSX2qNl2YIuLxwKuB3wC/CxwxepvMvAy4FXhcpwrsJ0NDdl+SJEmqu17vwtROC8SLgD3AozPzlwAx9qn0HwMnTbuyPrR5swFCkiSp7pYuhYGB3u3C1M4YiPsA3y/DwwRuBRwDMYayC5MkSZLqK6I4JuzVFoh2AsQiinAwmUOnWEvfswuTJEmSoD4BYj1wagvbnQZcN7Vy+pstEJIkSYLipHIdAsQ3gNMj4rfG2yAingocC/zPdAvrR70626AkSZI6a3AQtmypuoqpaSdAvAXYB3wiIv4oIlaWD0TEIRHxbOCfgZ3AuzpbZu/bs6dYDBCSJElavrwGASIzrwSeAywA/gm4GUjgWcA24EMU4yR+PzOv6XilPa7cQZYvr7YOSZIkVa8WAQIgM/8TOBP4JLCdYhK5ucBu4PPAAzPz450ush+Ul+kyQEiSJGn58t69jOu480BExDJgb2bubr4/M38GPDWKSSBWAgPAxsw8MKOV9rgyYdqFSZIkSYODsGMH7N8Pc9uZma0LTNQCsRl4b3kjIl4dEU8ob2dhY2bebHiYnF2YJEmSVCqPCbdurbaOqZgoQERjKb0GeOJMFtPPyiYqWyAkSZJUBoheHAcxUYDYiZPCdYwtEJIkSSqVJ5V7cRzERD2urgQeFRG/D1zVuO/IiHhoKy+cmZdOt7h+YoCQJElSqZdbICYKEO8HLgD+pem+cxrLZHKS166doSGIgKVLq65EkiRJVevLAJGZ/xoRG4CnAMcADwNuoWiZUJu2bIFly2BOWxfOlSRJUj8qA0S/dWEiMz9PMb8DETEMfCkzf382Cus3W7bYfUmSJEmFcgxEL7ZAjHs+PCKeEBFnNN31YeDbM15Rnxoa8gpMkiRJKixbVqz7KkAAFwMvarp9HLBqJovpZ7ZASJIkqTRvHhxySG92YZooQOSox88CTp3ZcvqXAUKSJEnNBgf7rwViC3DsbBXS7+zCJEmSpGbLl/dmgJhoEPXlwMMi4iOMzANxRkS8uoXXzcx83bSr6yO2QEiSJKlZPwaIvwPWAM9quu+MxjKeBKKxNkA0ZBogJEmSdLDly+G226quon0TzQPxrYg4Hfgt4M7Aa4AfA5+Zlcr6yI4dcOCAXZgkSZI0YnAQrrmm6iraN9k8EDcBFwJExGuAH2fma2e+rP5SNk3ZAiFJkqRSP3ZhGu25jIyFUBsMEJIkSRpt+fLevIxrywEiMz88k4X0s3LHsAuTJEmSSoODsGdPsSxYUHU1rZvoMq7qEFsgJEmSNFp5bNhr3ZjGDRARcSAi9kfEXZput7rsn70fofsZICRJkjRarwaIibowRWNpvt2qdrbte3ZhkiRJ0mhlgOi1cRATXcZ1zkS31TpbICRJkjRaeXK511ogDAWzYMsWmDsXFi2quhJJkiR1i17twmSAmAVDQ0XCDDt2SZIkqaHvujCNFhGHAw8DTgdWAsPAbcDPgG9m5sYZqbAPbNli9yVJkiQdrFe7ME0aICJiBfA24FnAwDib7YuIDwMvz8we+xXMPAOEJEmSRlu6tOih0lcBIiKOAL4J3IXiykq3AT8ENlJ0f1oF3AtYATwPeFBEnG1rxMHKLkySJElSac6cIkT0WxemC4BTgKuAl2TmF8faKCIeB7wDuCvwfuApnSyy123ZAiefXHUVkiRJ6jbLl/deC8REE8ndHXg8sA44c7zwAJCZnwfuC1wDPCkiTut0ob3MLkySJEkay+BgHwUI4OlAAi9tZVxDZm4GXkrR1enpnSmvP9iFSZIkSWPpqxYI4ExgS2Z+ro3X+xwwBNxvOkX1kwMHYNs2WyAkSZJ0R8uX994YiIkCxCnAj9p5scxMikHWp0ynqH6ybVuxNkBIkiRptH7rwjQI3DqF17y18VwxkijtwiRJkqTR+q0L02Jg5xRec3fjuWJkh7AFQpIkSaOVXZgyq66kdRMFiJjG607nuX3FACFJkqTxLF9ejJndOZXT9hWZbB6IIyPioW2+5pFTLaYf2YVJkiRJ4ymPEbdsgcU90odnsgBxTmPRFNkCIUmSpPGUx4hbtsBRR1VbS6smChDXU8wDMSsi4hTg0RSXj10D3IWiK9TvZuYnZ6uOTjNASJIkaTzlMWIvXcp13ACRmcfNYh0AzwdePMvfc8aVO4MBQpIkSaM1d2HqFRMNop5tPwfeAjwVOAm4pNpyOmPLFli0CObPr7oSSZIkdZvmLky9YrIxELMmM/+1+XZEf1zIacsWWx8kSZI0tl7swtRNLRB9aWjIKzBJkiRpbL3YAmGAmGG2QEiSJGk8ixfDwIABQk0MEJIkSRpPxMhs1L2ia8ZAdMrw8DDbtm2r7PuP/t633XYId7rTMNu27a6oIlWtyv1R3cv9QqO5T2g094n6WLZsMRs3Hpj0eLFb9om+aIGIiD+KiLURsXbjxo1Vl3OQrVuDZctmbToNSZIk9Zjly5OtW3vnAkJ90QKRmRcAFwCsWbMmly5dWnFFUNawdSusWjWfpUu9jmvddcN+qe7jfqHR3Cc0mvtE/1uxAnbsaP29rnqf6IsWiG61dy/s2uVVmCRJkjS+XhsDYYCYQeVoegdRS5IkaTzLl3sVJjUYICRJkjSZwUEDhBrKpii7MEmSJGk8y5cX42aHh6uupDVdM4g6Iu4NvK/prtMa6zdGxMvKOzPz/rNa2DTYAiFJkqTJLF8OmbBtW28cN3ZNgACWAfcb4/6TZ7uQTjFASJIkaTJlb5VemYC4awJEZn4T6J0L4LbALkySJEmaTBkaemUchGMgZpAtEJIkSZpMeazYK5dyNUDMoDJAOP+LJEmSxmMLhG43NATLlsHAQNWVSJIkqVs1j4HoBQaIGdQrA2EkSZJUHbsw6XYGCEmSJE3GLky63dCQV2CSJEnSxBYuhAULDBDCFghJkiS1ZvlyuzAJA4QkSZJas3y5LRDCLkySJElqjQFCZNoCIUmSpNYMDhogam/XLti/3wAhSZKkyTkGQrfvAHZhkiRJ0mTswqTbdwBbICRJkjQZuzDJACFJkqSWLV8OO3bAvn1VVzI5A8QMsQuTJEmSWlWedN66tdo6WmGAmCG2QEiSJKlV5TFjL3RjMkDMEAOEJEmSWlX2WjFA1JhdmCRJktSq8qRzL1zK1QAxQ7ZsgYEBOOSQqiuRJElSt7MLk26fhTqi6kokSZLU7ezCJIaG7L4kSZKk1tiFSbe3QEiSJEmTWbasWNsCUWMGCEmSJLVq3rxi7KwBosbswiRJkqR2DA7ahanWbIGQJElSO5YvtwWi1gwQkiRJaocBosaGh2HrVrswSZIkqXWDgwaI2tq2DTJtgZAkSVLrli93DERtbd1azB5ngJAkSVKr7MJUY1u2FAHCLkySJElqlQGixmyBkCRJUrsGB2HPHti9u+pKJmaAmAFlcjRASJIkqVXlsWO3t0IYIGaAXZgkSZLULgNEjdmFSZIkSe0qTz4bIGrIACFJkqR2lceO3X4pVwPEDNiyJVi4EBYsqLoSSZIk9Qq7MNXY1q22PkiSJKk9Boga27IlDBCSJElqSzkGwi5MNbRlS3gFJkmSJLVl6VKIsAWilrZutQVCkiRJ7ZkzpwgRBogacgyEJEmSpmJw0C5MtWQXJkmSJE3F8uW2QNSSXZgkSZI0FQaIGtq3D3buNEBIkiSpfQaIGipnobYLkyRJktrlGIgaKhOjLRCSJElqly0QNbRlS9ECYYCQJElSu8oAkVl1JeMzQHSYXZgkSZI0VYODcOAA7NhRdSXjM0B0WBkgbIGQJElSu8pjyG7uxmSA6DDHQEiSJGmqDBA1VI6BsAuTJEmS2mWAqKGyC9OyZRUXIkmSpJ5TnoTu5ku5GiA6bMuWYMmSZGCg6kokSZLUa2yBqKG73nWYc8/dX3UZkiRJ6kGHHgp3uQvMn191JeObW3UB/eb88/dx/vn7gHlVlyJJkqQec+SR8KtfVV3FxGyBkCRJktQyA4QkSZKklhkgJEmSJLXMACFJkiSpZQYISZIkSS0zQEiSJElqmQFCkiRJUssMEJIkSZJaZoCQJEmS1DIDhCRJkqSWGSAkSZIktcwAIUmSJKllBghJkiRJLTNASJIkSWqZAUKSJElSywwQkiRJklpmgJAkSZLUMgOEJEmSpJYZICRJkiS1zAAhSZIkqWUGCEmSJEktM0BIkiRJapkBQpIkSVLLIjOrrqGjIuJW4LqKy1gFbKy4BnUX9wmNxf1Co7lPaDT3CY02m/vEsZl52Og7+y5AdIOIWJuZa6quQ93DfUJjcb/QaO4TGs19QqN1wz5hFyZJkiRJLTNASJIkSWqZAWJmXFB1Aeo67hMai/uFRnOf0GjuExqt8n3CMRCSJEmSWmYLhCRJkqSWGSA6KCKeERHfiogtEbE9ItZGxAsjwt9zn4mIeRHxiIh4W+N93hoReyPixoj4ZEScPcnz3VdqIiLeGBHZWF42wXbuE30sIhZFxMsj4rKIGIqInRFxTUR8IiIeNMb2cxrv/9rG/rClsX88vYr61XkRcXREvDsifhURuyJid0T8OiLeHxEnTPA8Pyt6VEScEhEvjoiPRsSVETHc+N/wlBaeO6X3PSIeHRH/HRG3NT53fh4R/y8iFkzrZ7ELU2dExHuBFwC7ga8B+4BHAEuBTwNPyczh6ipUJ0XEI4H/adzcAFwO7ABOA+7WuP91mfnqMZ7rvlITEXEm8D2KkzUB/GVmvnWM7dwn+lhEHA/8N3ASsB74X2A/cCxwL+C1mfn6pu0HgP8CngBspdgnFlDsEwuAd2Xmi2fzZ1BnRcS9gK8Dg8BvKP6HAKwB7gRsB87JzO+Oep6fFT0sIv4RGOtv93cz85MTPG9K73tEvBz4B+AA8E1gM3AWcBjwfeARmblzSj9MZrpMcwF+B0iKfwwnN91/BPCLxmMvrrpOl46+5w8HPgk8ZIzHnkpxcJDAw9xX6rlQHOj9Arix8QGfwMvG2M59oo8XYDFwFTAMvAIYGPX4SuAuo+77i8b7fgVwRNP9J1OcsEjg3Kp/Npdp7RffbbyPFwDzmu6fB3yg8dhPRj3Hz4oeX4DnAW8GzgNOpDioT4oAMN5zpvS+U4TRYYqTm/drun8JcEnjee+Y8s9S9S+zHxZgbeONePYYj53V9MbPqbpWl1nbJ/618b5/YNT97is1WSjO+iTweODCCQKE+0QfL8DfN97Dd7e4/QBwc+M5Dx3j8fMbj/2g6p/NZcr7xMLGe5jA6jEeX930+CFN9/tZ0WdLiwFiSu87xUnOBF49xvNOoGiV2AMMTqV2+8tNU0QcDdwH2At8YvTjmXkJxRnII4H7z251qtCPGuujyzvcV+ojIu5HcRb5osz83ATbuU/0sYiYD/xh4+bbW3zaA4DDgd9k5qVjPP4Jiu4LZ0bEnaZfpSpwgKKVejI7gF3gZ0VdTfV9b3z2PKZx89/HeN7VFN1r5wOPnUptBojpu1djfUVm7hpnm8tGbav+d3Jjvb7pPveVGoiIhcCHgdsYu69rM/eJ/nYfii5KN2bmNRFx74h4XUT8c0T8XUQ8eIznlO/zZWM8Rhb9la9o3Dyj4xVrxmXmPop+7ACvjYh55WONr1/XuPmBbJwuxs+Kuprq+34KcAhwW2aua+N5LZs7lSfpIMc31tdNsM31o7ZVH4uII4HnNG5+qukh95V6eAPFh/fTMnPjJNu6T/S3uzfWN0bEWylapZq9KiIuBp6VmTsa97W6T5yB+0QvewHwZYoWqsdExNrG/WcCK4B/BF7etL2fFfU01ff9+FGPtfq8ltkCMX1LGusdE2yzvbFeOsO1qGIRMRf4KLAc+Nqo7ivuK30uIh4IvAS4ODM/1sJT3Cf626GN9b0owsM/UlyJaQVwLkXXgycC72t6jvtEDTS6kDwQ+BJFV9cnNpY7UQyM/VajpaLkflFPU33fZ3x/MUBInfV+ikur3QA8q+JaNIsiYhHFYOmtFGcXpfJ/7Dzgo5n555m5LjOHMvOzFAeMCfxeRJxYVZGafY2TDT+nCJTnUlxW8zCKfWIF8KmIuMNlwKVuYYCYvjLBLZ5gmzIJbpvhWlShiHgn8AcUl1l8RGZuGLWJ+0p/eyPF2JeXZub6yTZucJ/ob83v2b+MfjAz11Jc/z8orqYC7hN9LyIGgYspzvw+OjM/m5kbG8tngEdTDJ5+VUSU4+ncL+ppqu/7jO8vBojpu7axPnaCbY4Zta36TES8DXgRcCtFePj1GJtd21i7r/SnJ1Fcc/v8iPhm80JxQADw/MZ9/9q4fW1j7T7Rn64Z5+uxtjmysb62sXaf6F+/TWMir0ZXpoNk5lUUkw3OBc5u3H1tY+1+US/XNtbtvu/l13du83ktcxD19JWX6zw9IhaNM0r+zFHbqo9ExJuBlwKbgEdm5i/G2dR9pf/NYeRM8lhOaCyDjdvuE/2t+T1bSdG1cbRVjXV5xvCHjfWZY2xLRBzCyGz37hO9qTyo2zLBNkONdTmOxs+Keprq+34lRSvWoRFx4jhXYrrvGM9rmS0Q05SZN1B84M8Hfnf04xFxFsUAqQ0U19xVH4mINwF/STE9/KMy86fjbeu+0t8y87jMjLEWisu6Avxl474zGs9xn+hjmXkjxZlkKMZGHSQiVgD3btwsr8LzPYqWzKMj4qFjvOzvUoypuKzx+uo9NzXW92m+hGupcd99GjevAT8r6mqq73tm7qUYoA/wzDGedwLFnDN7gS9MpTYDRGf8fWP9DxFxUnlnRBzOyNU13pSZw7NemWZMRLweeAXFmaJHZWYrKd59RaO5T/S3NzTWr4yINeWdjflC/oniim2X0/jnn5kHgDc3Nvunxn5QPudk4E2jXle950vAToqWiHdExILygcbX76LoXrIZ+ErT8/ysqKepvu9vorhIwysi4r5Nz1sCfJAiA7wvM4emUlSMzFGi6YiI9wHPB3YDX6WYKfQRwDKKwVJPafxjUB+IiCcAn2ncXMvIxE6jXZmZb2q+w32lfiLiQuB8ihaIt47xuPtEH2uaA2If8H2K7o73BY6iuJTrw5rHTUXEAPBp4PEUV/X6GkWrwyOBhcC7M/NFs/kzqLMi4nzgA8AARYtE2XXtPsBqYA/FXDIXj3qenxU9LCLuzcGXbT6NYjD9rykmHwUgM+8/6nlTet8j4uXAP1DMfv51ihOeZ1HMdv+/wMMbk1O2/7MYIDonIp4BvJBi8qABij5oHwT+yTMC/SUingN8qIVNL8nMs8d4vvtKjUwWIBrbuE/0sYh4MvCnFHNCHEIxidNnKc4c3jrG9nMoLgf8XOBUigOAn1KcMbxoturWzGkcTL4EeAhFaIAiUH4DePt44+n8rOhdEXE2xfs7oUbX19HPndL7HhGPpjiBsYbiBMTVwEXAWzNzT9s/RPm6BghJkiRJrXIMhCRJkqSWGSAkSZIktcwAIUmSJKllBghJkiRJLTNASJIkSWqZAUKSJElSywwQkiRJklpmgJCkPhURZ0dEjlr2RcSmiLgyIj4eES+IiMEJXuM5jeddOHuVd1bT7+GbVdciSf3AACFJ/W8H8OHGchHwLWAX8CTgvcBNEfHnEXGH2U97QURc2wgIx1VdiyTVwdyqC5AkzbiNmfmc0XdGxErgRcArgbcDRwN/MWqzTwPfB7bMcI0z6QfAXYGdVRciSf3AFghJqqnM3JSZfws8BUjgpRFx9qhttmTmlZm5voISOyIzdzZ+huurrkWS+oEBQpK6QEQcFxEHIuK2iFg0zjbzImJ9o7vO6Z363pn5GeBTjZsvG/U9xxwD0TyuICIOiYjXN8ZV7IqIH4/a9pyI+GxE3BwRexs/w39ExN3HqykiVkbE30XEjyJia0TsiIhfR8SFEfHA5tqAYxtPu2bUeI/jRtc6zvc6PSI+EhE3RMSeiNgYEV+MiMeMs/2Fjdd7TkScFBEXNX62PY3fwSsiwv+vkvqWXZgkqQtk5rUR8TngXODpwAfH2Ox3gCOBb2bmFR0u4aMULRFnR8TczNzf4vMWAt+k6CJ0KfATYH75YES8k6Kb1H7gMuA3wEnA04AnRsTvZOYXm18wIu4FfAFYDdzWeP3dFEHh6Y3NvgtcRTGu4ynAYooQtL3ppZq/HlNEPAH4OLAAuIJifMjRwDnAYyLi9Zn5qnGefgbwTmAj8A3gcOAhwJsar/Fnk31/SepFBghJ6h7vpggQz2fsAPGCxvq9M/C91zbWi4HjKA7OW3E/4MfASZl5c/MDEfEnFOHhCuApmXll02NPBD4B/HtEnJCZmxv3LwE+SxEe3g+8NDN3NT3vMOAUgMz8NvDtRrerxcDLMvPaFusmIo4E/o0iPPxFZr696bGzKULM30TEtzPzK2O8xIuB1wJ/l5nDjec9lCJMvCAi3pyZN7RajyT1CptYJalLZObXgF8AayLivs2PNbr7PAS4Cbh4Br79xqavV7b53BeOER4GgFc3bp7XHB4AMvNi4J+BQeBZTQ89j+Ls/feAFzSHh8bzbm0Eh074Q2AZ8J3m8ND4Pt+kCHQwqltXk8uA15bhofG8S4GvUPx/fViH6pSkrmKAkKTu8p7G+gWj7i9vX9BG96J2NP8/GB53qzu6OTO/O8b9Z1C0IlyRmb8Y57mXNNYPaLrv0Y31BzMz26hjKs5qrD88zuNlK9CDG4FotC+OU2MZlo6aTnGS1K3swiRJ3eUjwN8DT42Il2bmbRGxjOIs/T7gAoCIOBX4qzGef3Hj7H67VjV9fVsbz7tunPtPaKxPbwx0nshhTV+XA6KvHGvDDrtTY33NOI9fSxGmFlK0ytwy6vHxruq0tbFeOJ3iJKlbGSAkqYtk5o6I+CDw58DvA28Fng0sAT7RdDnVI4Hzx3iJa5laF6f7NNbbGq/Rql3j3F+esb8R+Ookr9EcFma61WEsU/2e7bTUSFLfMEBIUvd5L8UA3T+JiLdTDKou7wdu76PfyZmjy3EI38jMAx14vXLw8PqxJrGbwPUUV3Q6BejUWIfx3AicStFa8rUxHj+OomvXbtprlZGkvuYYCEnqMpm5DvgScCLwRuA0irEEl0z4xClqXMr0yY2bb+nQy/4A2ATcKyJOauN55dWOfj8iWg1Iexvrdk+Klb/PZ4/z+HMb62/P0LgTSepJBghJ6k7lFYBe0Vi/r9PfICIOjYi/pZg/IYB/6NQVjjJzH/A6iq5MF4++qlTj+8+PiCc0xnOU/pXiSlMPBN4dEQtHPeewiHjwqJe6sbG+a5tl/gtFl60HR8SLRn2fhzIyj8Pb2nxdSeprdmGSpO7038CvKLrybKOYr2CqVjXNJB3AUoruOXen+D+wk+JSpe+axve4g8x8Z0QcSzGe438j4qfAOooWgzsB96KYv+ExNMZBZOa2iDiXYg6GFwJPi4jvMDKR3L2A/+Dg7k2fBs6mmFPiv4Ghxv2vyMxNE9S3ISJ+D/gY8M6IeB7wc4qrJz2E4iTb6zPzy9P8VUhSXzFASFIXysyMiK9SBIiPZOa2abzcYkYGXO+nCCS3UBx4fwP4j8wcmsbrjyszXxoRF1OM43gQ8NsUA6/XA5+nmDTuW6Oes7Yx78WfA48HHkUxYPkm4CKK+SOavYdiPodnAo+jmBgO4PUU3agmqu8zEbGGoqXn4RSzWm+jCHDvHj1LtiQJYuYvsy1JaldEzKcYUHwEcPoEcylIkjSrHAMhSd3phRTh4cuGB0lSN7EFQpK6REScAvwlRR/8c4ADwJmZ+ZNKC5MkqYljICSpe6wG/gDYA/wE+BvDgySp29gCIUmSJKlljoGQJEmS1DIDhCRJkqSWGSAkSZIktcwAIUmSJKllBghJkiRJLTNASJIkSWrZ/wfL8OwdidBOMgAAAABJRU5ErkJggg==\n", + "text/plain": [ + "<Figure size 936x720 with 1 Axes>" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "fig = plt.figure(figsize=(13, 10))\n", + "ax = plt.gca()\n", + "fontsize = 22\n", + "\n", + "plt.xlabel('y-Direction', fontsize=fontsize)\n", + "plt.ylabel('Difference in Velocity Magnitude', fontsize=fontsize)\n", + "\n", + "plt.grid(color='black', linestyle='-', linewidth=0.1)\n", + "\n", + "plt.rc('font', size=fontsize)\n", + "\n", + "plt.plot(vel_profile_ch_with_cm - vel_profile_ch_normal_srt, color='b')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Compare Operation Counts of Both Described Methods" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First we set up a method directly with the central moments." + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "cm = get_central_moments(u = u)\n", + "mtrr = OrderedDict(zip(cm, ω))" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + " <table style=\"border:none; width: 100%\">\n", + " <tr style=\"border:none\">\n", + " <th style=\"border:none\" >Moment</th>\n", + " <th style=\"border:none\" >Eq. Value </th>\n", + " <th style=\"border:none\" >Relaxation Rate</th>\n", + " </tr>\n", + " <tr style=\"border:none\">\n", + " <td style=\"border:none\">$1$</td>\n", + " <td style=\"border:none\">$\\rho$</td>\n", + " <td style=\"border:none\">$1.99$</td>\n", + " </tr>\n", + "<tr style=\"border:none\">\n", + " <td style=\"border:none\">$- u_{0} + x$</td>\n", + " <td style=\"border:none\">$0$</td>\n", + " <td style=\"border:none\">$1.99$</td>\n", + " </tr>\n", + "<tr style=\"border:none\">\n", + " <td style=\"border:none\">$- u_{1} + y$</td>\n", + " <td style=\"border:none\">$0$</td>\n", + " <td style=\"border:none\">$1.99$</td>\n", + " </tr>\n", + "<tr style=\"border:none\">\n", + " <td style=\"border:none\">$\\left(- u_{0} + x\\right)^{2}$</td>\n", + " <td style=\"border:none\">$\\frac{\\rho}{3}$</td>\n", + " <td style=\"border:none\">$1.99$</td>\n", + " </tr>\n", + "<tr style=\"border:none\">\n", + " <td style=\"border:none\">$\\left(- u_{1} + y\\right)^{2}$</td>\n", + " <td style=\"border:none\">$\\frac{\\rho}{3}$</td>\n", + " <td style=\"border:none\">$1.99$</td>\n", + " </tr>\n", + "<tr style=\"border:none\">\n", + " <td style=\"border:none\">$\\left(- u_{0} + x\\right) \\left(- u_{1} + y\\right)$</td>\n", + " <td style=\"border:none\">$0$</td>\n", + " <td style=\"border:none\">$1.99$</td>\n", + " </tr>\n", + "<tr style=\"border:none\">\n", + " <td style=\"border:none\">$\\left(- u_{0} + x\\right)^{2} \\left(- u_{1} + y\\right)$</td>\n", + " <td style=\"border:none\">$0$</td>\n", + " <td style=\"border:none\">$1.99$</td>\n", + " </tr>\n", + "<tr style=\"border:none\">\n", + " <td style=\"border:none\">$\\left(- u_{0} + x\\right) \\left(- u_{1} + y\\right)^{2}$</td>\n", + " <td style=\"border:none\">$0$</td>\n", + " <td style=\"border:none\">$1.99$</td>\n", + " </tr>\n", + "<tr style=\"border:none\">\n", + " <td style=\"border:none\">$\\left(- u_{0} + x\\right)^{2} \\left(- u_{1} + y\\right)^{2}$</td>\n", + " <td style=\"border:none\">$\\frac{\\rho}{9}$</td>\n", + " <td style=\"border:none\">$1.99$</td>\n", + " </tr>\n", + "\n", + " </table>\n", + " " + ], + "text/plain": [ + "<lbmpy.methods.momentbased.MomentBasedLbMethod at 0x7f9b73948f70>" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "method_dir = create_with_discrete_maxwellian_eq_moments(stencil, mtrr, compressible=True, equilibrium_order=4)\n", + "\n", + "\n", + "def modification_func(moment, eq, rr):\n", + " rr = 1.99\n", + " return moment, eq.subs(u[0], 0).subs(u[1], 0), rr\n", + "method_dir = create_lb_method_from_existing(method_dir, modification_func)\n", + "method_dir" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now let us first look at the moment matrix, to show, that it is indeed not sparse anymore." + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\left[\\begin{matrix}1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1\\\\- u_{0} & - u_{0} & - u_{0} & - u_{0} - 1 & 1 - u_{0} & - u_{0} - 1 & 1 - u_{0} & - u_{0} - 1 & 1 - u_{0}\\\\- u_{1} & 1 - u_{1} & - u_{1} - 1 & - u_{1} & - u_{1} & 1 - u_{1} & 1 - u_{1} & - u_{1} - 1 & - u_{1} - 1\\\\u_{0}^{2} & u_{0}^{2} & u_{0}^{2} & \\left(- u_{0} - 1\\right)^{2} & \\left(1 - u_{0}\\right)^{2} & \\left(- u_{0} - 1\\right)^{2} & \\left(1 - u_{0}\\right)^{2} & \\left(- u_{0} - 1\\right)^{2} & \\left(1 - u_{0}\\right)^{2}\\\\u_{1}^{2} & \\left(1 - u_{1}\\right)^{2} & \\left(- u_{1} - 1\\right)^{2} & u_{1}^{2} & u_{1}^{2} & \\left(1 - u_{1}\\right)^{2} & \\left(1 - u_{1}\\right)^{2} & \\left(- u_{1} - 1\\right)^{2} & \\left(- u_{1} - 1\\right)^{2}\\\\u_{0} u_{1} & - u_{0} \\left(1 - u_{1}\\right) & - u_{0} \\left(- u_{1} - 1\\right) & - u_{1} \\left(- u_{0} - 1\\right) & - u_{1} \\left(1 - u_{0}\\right) & \\left(1 - u_{1}\\right) \\left(- u_{0} - 1\\right) & \\left(1 - u_{0}\\right) \\left(1 - u_{1}\\right) & \\left(- u_{0} - 1\\right) \\left(- u_{1} - 1\\right) & \\left(1 - u_{0}\\right) \\left(- u_{1} - 1\\right)\\\\- u_{0}^{2} u_{1} & u_{0}^{2} \\left(1 - u_{1}\\right) & u_{0}^{2} \\left(- u_{1} - 1\\right) & - u_{1} \\left(- u_{0} - 1\\right)^{2} & - u_{1} \\left(1 - u_{0}\\right)^{2} & \\left(1 - u_{1}\\right) \\left(- u_{0} - 1\\right)^{2} & \\left(1 - u_{0}\\right)^{2} \\left(1 - u_{1}\\right) & \\left(- u_{0} - 1\\right)^{2} \\left(- u_{1} - 1\\right) & \\left(1 - u_{0}\\right)^{2} \\left(- u_{1} - 1\\right)\\\\- u_{0} u_{1}^{2} & - u_{0} \\left(1 - u_{1}\\right)^{2} & - u_{0} \\left(- u_{1} - 1\\right)^{2} & u_{1}^{2} \\left(- u_{0} - 1\\right) & u_{1}^{2} \\left(1 - u_{0}\\right) & \\left(1 - u_{1}\\right)^{2} \\left(- u_{0} - 1\\right) & \\left(1 - u_{0}\\right) \\left(1 - u_{1}\\right)^{2} & \\left(- u_{0} - 1\\right) \\left(- u_{1} - 1\\right)^{2} & \\left(1 - u_{0}\\right) \\left(- u_{1} - 1\\right)^{2}\\\\u_{0}^{2} u_{1}^{2} & u_{0}^{2} \\left(1 - u_{1}\\right)^{2} & u_{0}^{2} \\left(- u_{1} - 1\\right)^{2} & u_{1}^{2} \\left(- u_{0} - 1\\right)^{2} & u_{1}^{2} \\left(1 - u_{0}\\right)^{2} & \\left(1 - u_{1}\\right)^{2} \\left(- u_{0} - 1\\right)^{2} & \\left(1 - u_{0}\\right)^{2} \\left(1 - u_{1}\\right)^{2} & \\left(- u_{0} - 1\\right)^{2} \\left(- u_{1} - 1\\right)^{2} & \\left(1 - u_{0}\\right)^{2} \\left(- u_{1} - 1\\right)^{2}\\end{matrix}\\right]$" + ], + "text/plain": [ + "⎡ 1 1 1 1 1 \n", + "⎢ \n", + "⎢ -u₀ -u₀ -u₀ -u₀ - 1 1 - u₀ \n", + "⎢ \n", + "⎢ -u₁ 1 - u₁ -u₁ - 1 -u₁ -u₁ \n", + "⎢ \n", + "⎢ 2 2 2 2 2 \n", + "⎢ u₀ u₀ u₀ (-u₀ - 1) (1 - u₀) (\n", + "⎢ \n", + "⎢ 2 2 2 2 2 \n", + "⎢ u₁ (1 - u₁) (-u₁ - 1) u₁ u₁ (\n", + "⎢ \n", + "⎢ u₀⋅u₁ -u₀⋅(1 - u₁) -u₀⋅(-u₁ - 1) -u₁⋅(-u₀ - 1) -u₁⋅(1 - u₀) (1 - \n", + "⎢ \n", + "⎢ 2 2 2 2 2 \n", + "⎢-u₀ ⋅u₁ u₀ ⋅(1 - u₁) u₀ ⋅(-u₁ - 1) -u₁⋅(-u₀ - 1) -u₁⋅(1 - u₀) (1 - u\n", + "⎢ \n", + "⎢ 2 2 2 2 2 \n", + "⎢-u₀⋅u₁ -u₀⋅(1 - u₁) -u₀⋅(-u₁ - 1) u₁ ⋅(-u₀ - 1) u₁ ⋅(1 - u₀) (1 - u\n", + "⎢ \n", + "⎢ 2 2 2 2 2 2 2 2 2 2 \n", + "⎣u₀ ⋅u₁ u₀ ⋅(1 - u₁) u₀ ⋅(-u₁ - 1) u₁ ⋅(-u₀ - 1) u₁ ⋅(1 - u₀) (1 - u\n", + "\n", + " 1 1 1 1 \n", + " \n", + "-u₀ - 1 1 - u₀ -u₀ - 1 1 - u₀ \n", + " \n", + " 1 - u₁ 1 - u₁ -u₁ - 1 -u₁ - 1 \n", + " \n", + " 2 2 2 2 \n", + "-u₀ - 1) (1 - u₀) (-u₀ - 1) (1 - u₀) \n", + " \n", + " 2 2 2 2 \n", + "1 - u₁) (1 - u₁) (-u₁ - 1) (-u₁ - 1) \n", + " \n", + "u₁)⋅(-u₀ - 1) (1 - u₀)⋅(1 - u₁) (-u₀ - 1)⋅(-u₁ - 1) (1 - u₀)⋅(-u₁ - 1\n", + " \n", + " 2 2 2 2 \n", + "₁)⋅(-u₀ - 1) (1 - u₀) ⋅(1 - u₁) (-u₀ - 1) ⋅(-u₁ - 1) (1 - u₀) ⋅(-u₁ - 1\n", + " \n", + " 2 2 2 \n", + "₁) ⋅(-u₀ - 1) (1 - u₀)⋅(1 - u₁) (-u₀ - 1)⋅(-u₁ - 1) (1 - u₀)⋅(-u₁ - 1)\n", + " \n", + " 2 2 2 2 2 2 2 \n", + "₁) ⋅(-u₀ - 1) (1 - u₀) ⋅(1 - u₁) (-u₀ - 1) ⋅(-u₁ - 1) (1 - u₀) ⋅(-u₁ - 1\n", + "\n", + " ⎤\n", + " ⎥\n", + " ⎥\n", + " ⎥\n", + " ⎥\n", + " ⎥\n", + " ⎥\n", + " ⎥\n", + " ⎥\n", + " ⎥\n", + " ⎥\n", + " ⎥\n", + ") ⎥\n", + " ⎥\n", + " ⎥\n", + ") ⎥\n", + " ⎥\n", + "2 ⎥\n", + " ⎥\n", + " ⎥\n", + " 2⎥\n", + ") ⎦" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "method_dir.moment_matrix" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, we want to know how many operations the collision rule has and how much we can reduce them." + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [], + "source": [ + "collision_rule_dir = method_dir.get_collision_rule()" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "<table style=\"border:none\"><tr><th>Name</th><th>Runtime</th><th>Adds</th><th>Muls</th><th>Divs</th><th>Total</th></tr><tr><td>OriginalTerm</td><td>-</td> <td>1543</td> <td>1917</td> <td>2</td> <td>3462</td> </tr><tr><td>sympy_cse</td><td>155.96 ms</td> <td>248</td> <td>179</td> <td>1</td> <td>428</td> </tr></table>" + ], + "text/plain": [ + "<pystencils.simp.simplificationstrategy.SimplificationStrategy.create_simplification_report.<locals>.Report at 0x7f9b738c1100>" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "generic_strategy_dir = ps.simp.SimplificationStrategy()\n", + "generic_strategy_dir.add(ps.simp.sympy_cse)\n", + "generic_strategy_dir.create_simplification_report(collision_rule_dir)" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [], + "source": [ + "collision_rule_dir = ps.simp.sympy_cse(collision_rule_dir)" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'adds': 248,\n", + " 'muls': 179,\n", + " 'divs': 1,\n", + " 'sqrts': 0,\n", + " 'fast_sqrts': 0,\n", + " 'fast_inv_sqrts': 0,\n", + " 'fast_div': 0}" + ] + }, + "execution_count": 21, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "collision_rule_dir.operation_count" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'adds': 393,\n", + " 'muls': 246,\n", + " 'divs': 1,\n", + " 'sqrts': 0,\n", + " 'fast_sqrts': 0,\n", + " 'fast_inv_sqrts': 0,\n", + " 'fast_div': 0}" + ] + }, + "execution_count": 22, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "collision_rule.operation_count" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As we can see, not only is the initial operation count - before the subexpression elimination - much larger than the operation count with the shift matrix (116002 vs 2106 operations), but also the operation count after the common subexpression elimination is larger, and that is already the case for D2Q9, for D3Q27 the difference will get even more extensive." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Literature" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "- Martin Geier, Andreas Greiner, and Jan G. Korvink. \"Cascaded digital lattice Boltzmann automata for high Reynolds number flow\". In: _Phys. Rev. E_ 73 (6. June 2006), p. 066705. DOI: 10.1103/PhysRevE.73.066705. URL: https://link.aps.prg/doi/10.1103/PhysRevE.73.066705.\n", + "\n", + "- G. Gruszczyński et al. \"A cascaded phase-field lattice Boltzmann model for the simulation of incompressible, immiscible fluids with high density contrast\". In: _Computers Mathematics with Applications_ 79.4 (2020), pp. 1049-1071. ISSN:0898-1221. DOI: https://doi.org/10.1016/j.camwa.2019.08.018. URL: http://www.sciencedirect.com/science/article/pii/S0898122119304158\n", + "\n", + "- Markus Muhr. \"Kumulanten-basierte Lattice-Boltzmann-Methode\". Seminar Ausarbeitung in _Transportprozesse: Mathematische Modellierung und Numerische Methoden_. Technische Universität München Fakultät für Mathematik." + ] + } + ], + "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.8.6" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/lbmpy/methods/momentbased.py b/lbmpy/methods/momentbased.py index 3a6bb217c9341eacfd8bcae048a8dcc4058e18d8..53b901a25dfee07d82bfa6b75228607026714515 100644 --- a/lbmpy/methods/momentbased.py +++ b/lbmpy/methods/momentbased.py @@ -5,7 +5,7 @@ import sympy as sp from lbmpy.maxwellian_equilibrium import get_weights from lbmpy.methods.abstractlbmethod import AbstractLbMethod, LbmCollisionRule, RelaxationInfo from lbmpy.methods.conservedquantitycomputation import AbstractConservedQuantityComputation -from lbmpy.moments import MOMENT_SYMBOLS, moment_matrix +from lbmpy.moments import MOMENT_SYMBOLS, moment_matrix, shift_matrix from pystencils import Assignment from pystencils.sympyextensions import subs_additive @@ -132,6 +132,10 @@ class MomentBasedLbMethod(AbstractLbMethod): def moment_matrix(self): return moment_matrix(self.moments, self.stencil) + @property + def shift_matrix(self): + return shift_matrix(self.moments, self.stencil) + @property def relaxation_matrix(self): d = sp.zeros(len(self.relaxation_rates)) diff --git a/lbmpy/moments.py b/lbmpy/moments.py index 13a84ce62c52cff5ac155c96ec703dd20f04d8fa..b1e5af5833c808a71dec6c4e1347926b62a752dc 100644 --- a/lbmpy/moments.py +++ b/lbmpy/moments.py @@ -39,6 +39,7 @@ import math from collections import Counter, defaultdict from copy import copy from typing import Iterable, List, Optional, Sequence, Tuple, TypeVar +from numpy import asarray import sympy as sp @@ -368,6 +369,55 @@ def moment_matrix(moments, stencil, shift_velocity=None): return sp.Matrix(len(moments), len(stencil), generator) +def shift_matrix(moments, stencil): + """ + Returns shift matrix to shift raw moments to central moment space + """ + x, y, z = MOMENT_SYMBOLS + dim = len(stencil[0]) + nr_directions = len(stencil) + + directions = asarray(stencil) + + u = sp.symbols("u_:{dim}".format(dim=dim)) + f = sp.symbols("f_:{dim}".format(dim=nr_directions)) + m = sp.symbols("m_:{dim}".format(dim=nr_directions)) + + Mf = moment_matrix(moments, stencil=stencil) * sp.Matrix(f) + + shift_matrix_list = [] + shift_equations = [] + for nr_moment in range(len(moments)): + exponent_central_moment = moments[nr_moment].as_powers_dict() + + exponent = [exponent_central_moment[x], exponent_central_moment[y], exponent_central_moment[z]] + + shift_equation = 1 + for i in range(dim): + shift_equation *= (directions[:, i] - u[i]) ** exponent[i] + shift_equation = sum(shift_equation * f) + + collected_expr = sp.Poly(shift_equation, u).as_expr() + terms_of_collected_expr = collected_expr.args + + constants = set() + for i in range(len(terms_of_collected_expr)): + constants.add(sp.Abs(sp.LC(terms_of_collected_expr[i]))) + + for i in range(len(stencil)): + for c in constants: + collected_expr = collected_expr.subs(c * Mf[i], c * m[i]) + + collected_expr = sp.collect(collected_expr.expand(), m) + + inner_list = [] + for i in range(len(stencil)): + inner_list.append(collected_expr.coeff(m[i], 1)) + + shift_matrix_list.append(inner_list) + return sp.Matrix(shift_matrix_list) + + def gram_schmidt(moments, stencil, weights=None): r""" Computes orthogonal set of moments using the method by Gram-Schmidt