![](\"https://raw.githack.com/spectralDNS/spectralutilities/master/movies/KleinGordon.gif\")
Figure 1
\n", "\n", "\n", "Movie showing the evolution of the solution $u$ from the Klein-Gordon equation, in a slice through the center of the domain, computed with the code described in this demo." ] }, { "cell_type": "markdown", "id": "8ae995a4", "metadata": { "editable": true }, "source": [ "## The nonlinear Klein-Gordon equation\n", "\n", "The cubic nonlinear Klein-Gordon equation is a wave equation important for many\n", "scientific applications such as solid state physics, nonlinear optics and\n", "quantum field theory [[abdul08]](#abdul08). The equation is given as" ] }, { "cell_type": "markdown", "id": "83b245d3", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{\\partial^2 u}{\\partial t^2} = \\nabla^2 u - \\gamma(u - u|u|^2) \\quad\n", "\\text{for} \\, u \\in\n", "\\Omega, \\label{eq:kg} \\tag{1}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "10fb3c95", "metadata": { "editable": true }, "source": [ "with initial conditions" ] }, { "cell_type": "markdown", "id": "e4bb65b4", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "u(\\boldsymbol{x}, t=0) = u^0 \\quad \\text{and} \\quad \\frac{\\partial u(\\boldsymbol{x},\n", "t=0)}{\\partial t} = u_t^0. \\label{eq:init} \\tag{2}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "10adb615", "metadata": { "editable": true }, "source": [ "The spatial coordinates are here denoted as $\\boldsymbol{x} = (x, y, z)$, and\n", "$t$ is time. The parameter $\\gamma=\\pm 1$ determines whether the equations are focusing\n", "($+1$) or defocusing ($-1$) (in the movie we have used $\\gamma=1$).\n", "The domain $\\Omega=[-2\\pi, 2\\pi)^3$ is triply\n", "periodic and initial conditions will here be set as" ] }, { "cell_type": "markdown", "id": "40131bb3", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "u^0 = 0.1 \\exp \\left( -\\boldsymbol{x} \\cdot \\boldsymbol{x} \\right), \n", "\\label{_auto1} \\tag{3}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "53595408", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", "u_t^0 = 0.\n", "\\label{_auto2} \\tag{4}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "4c365136", "metadata": { "editable": true }, "source": [ "We will solve these equations using a mixed formulation and a spectral Galerkin\n", "method. The mixed formulation reads" ] }, { "cell_type": "markdown", "id": "66f1256e", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{\\partial f}{\\partial t} = \\nabla^2 u - \\gamma (u - u|u|^2), \\label{eq:df} \\tag{5}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "8eb65116", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", "\\frac{\\partial u}{\\partial t} = f. \\label{eq:du} \\tag{6}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "ea17eddc", "metadata": { "editable": true }, "source": [ "The energy of the solution can be computed as" ] }, { "cell_type": "markdown", "id": "9e81800c", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "E(u) = \\int_{\\Omega} \\left( \\frac{1}{2} f^2 + \\frac{1}{2}|\\nabla u|^2 + \\gamma(\\frac{1}{2}u^2 - \\frac{1}{4}u^4) \\right) dx\n", "\\label{_auto3} \\tag{7}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "b53f8a5b", "metadata": { "editable": true }, "source": [ "and it is crucial that this energy remains constant in time.\n", "\n", "The movie above is showing the solution $u$, computed with the\n", "code shown below." ] }, { "cell_type": "markdown", "id": "c6e00bc7", "metadata": { "editable": true }, "source": [ "## Spectral Galerkin formulation\n", "\n", "The PDEs in ([5](#eq:df)) and ([6](#eq:du)) can be solved with many different\n", "numerical methods. We will here use the [shenfun](https://github.com/spectralDNS/shenfun) software and this software makes use of\n", "the spectral Galerkin method. Being a Galerkin method, we need to reshape the\n", "governing equations into proper variational forms, and this is done by\n", "multiplying ([5](#eq:df)) and ([6](#eq:du)) with the complex conjugate of proper\n", "test functions and then integrating\n", "over the domain. To this end we make use of the triply periodic tensor product\n", "function space $W^{\\boldsymbol{N}}(\\Omega)$ (defined in Eq. ([14](#eq:kg:Wn)))\n", "and use testfunctions $g \\in W^{\\boldsymbol{N}}$\n", "with Eq. ([5](#eq:df)) and $v \\in W^{\\boldsymbol{N}}$ with Eq. ([6](#eq:du)),\n", "and obtain" ] }, { "cell_type": "markdown", "id": "fedc6701", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{\\partial}{\\partial t} \\int_{\\Omega} f\\, \\overline{g}\\, w \\,d\\Omega = \\int_{\\Omega}\n", "\\left(\\nabla^2 u - \\gamma( u\\, - u|u|^2) \\right) \\overline{g} \\, w \\,d\\Omega, \\label{eq:df_var} \\tag{8} \n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "08ad1175", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", "\\frac{\\partial }{\\partial t} \\int_{\\Omega} u\\, \\overline{v}\\, w \\, dx =\n", "\\int_{\\Omega} f\\, \\overline{v} \\, w \\, d\\Omega. \\label{eq:kg:du_var} \\tag{9}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "fbf5fa50", "metadata": { "editable": true }, "source": [ "Note that the overline is used to indicate a complex conjugate, and\n", "$w$ is a weight function associated with the test functions. The functions\n", "$f$ and $u$ are now\n", "to be considered as trial functions, and the integrals over the\n", "domain are referred to as inner products. With inner product notation" ] }, { "cell_type": "markdown", "id": "4979f581", "metadata": { "editable": true }, "source": [ "$$\n", "\\left(u, v\\right) = \\int_{\\Omega} u \\, \\overline{v} \\, w\\, dx.\n", "$$" ] }, { "cell_type": "markdown", "id": "3e695b09", "metadata": { "editable": true }, "source": [ "and an integration by parts on the Laplacian, the variational problem can be\n", "formulated as:" ] }, { "cell_type": "markdown", "id": "36f5f08f", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{\\partial}{\\partial t} (f, g) = -(\\nabla u, \\nabla g)\n", "-\\gamma \\left( u - u|u|^2, g \\right), \\label{eq:df_var2} \\tag{10} \n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "781f3602", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", "\\frac{\\partial }{\\partial t} (u, v) = (f, v). \\label{eq:kg:du_var2} \\tag{11}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "24226db7", "metadata": { "editable": true }, "source": [ "The time and space discretizations are\n", "still left open. There are numerous different approaches that one could take for\n", "discretizing in time, and the first two terms on the right hand side of\n", "([10](#eq:df_var2)) can easily be treated implicitly as well as explicitly. However,\n", "the approach we will follow in Sec. ([Runge-Kutta integrator](#sec:rk)) is a fully explicit 4th order [Runge-Kutta](https://en.wikipedia.org/wiki/Runge-Kutta_methods) method. Also note that\n", "the inner product in the demo will be computed numerically with quadrature\n", "through fast Fourier transforms, and the integrals are thus not computed exactly\n", "for all terms." ] }, { "cell_type": "markdown", "id": "1465f0a3", "metadata": { "editable": true }, "source": [ "## Discretization\n", "To find a numerical solution we need to discretize the continuous problem\n", "([10](#eq:df_var2)) and ([11](#eq:kg:du_var2)) in space as well as time. Since the\n", "problem is triply periodic, Fourier exponentials are normally the best choice\n", "for trial and test functions, and as such we use basis functions" ] }, { "cell_type": "markdown", "id": "4d5445cd", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\phi_l(x) = e^{\\imath \\underline{l} x}, \\quad -\\infty < l < \\infty,\n", "\\label{_auto4} \\tag{12}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "5f2a57dd", "metadata": { "editable": true }, "source": [ "where $l$ is the wavenumber, and\n", "$\\underline{l}=\\frac{2\\pi}{L}l$ is the scaled wavenumber, scaled with domain\n", "length $L$ (here $4\\pi$). Since we want to solve these equations on a computer, we need to choose\n", "a finite number of test functions. A function space $V^N$ can be defined as" ] }, { "cell_type": "markdown", "id": "b90d937b", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "V^N(x) = \\text{span} \\{\\phi_l(x)\\}_{l\\in \\boldsymbol{l}}, \\label{eq:kg:Vn} \\tag{13}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "01bd92d5", "metadata": { "editable": true }, "source": [ "where $N$ is chosen as an even positive integer and $\\boldsymbol{l} = -N/2,\n", "-N/2+1, \\ldots, N/2-1$. And now, since $\\Omega$ is a\n", "three-dimensional domain, we can create tensor products of such bases to get,\n", "e.g., for three dimensions" ] }, { "cell_type": "markdown", "id": "59d5f2f3", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "W^{\\boldsymbol{N}}(x, y, z) = V^N(x) \\otimes V^N(y) \\otimes V^N(z), \\label{eq:kg:Wn} \\tag{14}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "68961d5f", "metadata": { "editable": true }, "source": [ "where $\\boldsymbol{N} = (N, N, N)$. Obviously, it is not necessary to use the\n", "same number ($N$) of basis functions for each direction, but it is done here\n", "for simplicity. A 3D tensor product basis function is now defined as" ] }, { "cell_type": "markdown", "id": "04821397", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\Phi_{lmn}(x,y,z) = e^{\\imath \\underline{l} x} e^{\\imath \\underline{m} y}\n", "e^{\\imath \\underline{n} z} = e^{\\imath\n", "(\\underline{l}x + \\underline{m}y + \\underline{n}z)},\n", "\\label{_auto5} \\tag{15}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "e9a3ffab", "metadata": { "editable": true }, "source": [ "where the indices for $y$- and $z$-direction are $\\underline{m}=\\frac{2\\pi}{L}m,\n", "\\underline{n}=\\frac{2\\pi}{L}n$, and $\\boldsymbol{m}$ and $\\boldsymbol{n}$ are the same as\n", "$\\boldsymbol{l}$ due to using the same number of basis functions for each direction. One\n", "distinction, though, is that for the $z$-direction expansion coefficients are only stored for\n", "$n=0, 1, \\ldots, N/2$ due to Hermitian symmetry (real input data). However, for simplicity,\n", "we still write the sum in Eq. ([16](#eq:usg)) over the entire range of basis functions.\n", "\n", "We now look for solutions of the form" ] }, { "cell_type": "markdown", "id": "382c00cb", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "u(x, y, z, t) = \\sum_{l=-N/2}^{N/2-1}\\sum_{m=-N/2}^{N/2-1}\\sum_{n=-N/2}^{N/2-1}\n", "\\hat{u}_{lmn} (t)\\Phi_{lmn}(x,y,z). \\label{eq:usg} \\tag{16}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "016c1249", "metadata": { "editable": true }, "source": [ "The expansion coefficients $\\hat{\\boldsymbol{u}} = \\{\\hat{u}_{lmn}(t)\\}_{(l,m,n) \\in \\boldsymbol{l} \\times \\boldsymbol{m} \\times \\boldsymbol{n}}$\n", "can be related directly to the solution $u(x, y, z, t)$ using Fast\n", "Fourier Transforms (FFTs) if we are satisfied with obtaining\n", "the solution in quadrature points corresponding to" ] }, { "cell_type": "markdown", "id": "84463f92", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", " x_i = \\frac{4 \\pi i}{N}-2\\pi \\quad \\forall \\, i \\in \\boldsymbol{i},\n", "\\text{where}\\, \\boldsymbol{i}=(0,1,\\ldots,N-1), \n", "\\label{_auto6} \\tag{17}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "1f16ed28", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", " y_j = \\frac{4 \\pi j}{N}-2\\pi \\quad \\forall \\, j \\in \\boldsymbol{j},\n", "\\text{where}\\, \\boldsymbol{j}=(0,1,\\ldots,N-1), \n", "\\label{_auto7} \\tag{18}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "28257dfb", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", " z_k = \\frac{4 \\pi k}{N}-2\\pi \\quad \\forall \\, k \\in \\boldsymbol{k},\n", "\\text{where}\\, \\boldsymbol{k}=(0,1,\\ldots,N-1).\n", "\\label{_auto8} \\tag{19}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "df5dd8f4", "metadata": { "editable": true }, "source": [ "Note that these points are different from the standard (like $2\\pi j/N$) since\n", "the domain\n", "is set to $[-2\\pi, 2\\pi]^3$ and not the more common $[0, 2\\pi]^3$. We have" ] }, { "cell_type": "markdown", "id": "e5ad43ee", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\boldsymbol{u} = \\mathcal{F}_x^{-1}\\left(\\mathcal{F}_y^{-1}\\left(\\mathcal{F}_z^{-1}\\left(\\hat{\\boldsymbol{u}}\\right)\\right)\\right) \\label{eq:uxyz} \\tag{20}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "332d71ce", "metadata": { "editable": true }, "source": [ "with $\\boldsymbol{u} = \\{u(x_i, y_j, z_k)\\}_{(i,j,k)\\in \\boldsymbol{i} \\times \\boldsymbol{j} \\times \\boldsymbol{k}}$\n", "and where $\\mathcal{F}_x^{-1}$ is the inverse Fourier transform along the direction $x$, for\n", "all indices in the other direction. Note that the three\n", "inverse FFTs are performed sequentially, one direction at the time, and that there is no\n", "scaling factor due to\n", "the definition used for the inverse [Fourier transform](https://mpi4py-fft.readthedocs.io/en/latest/dft.html)" ] }, { "cell_type": "markdown", "id": "e18610fd", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "u(x_j) = \\sum_{l=-N/2}^{N/2-1} \\hat{u}_l e^{\\imath \\underline{l}\n", "x_j}, \\quad \\,\\, \\forall \\, j \\in \\, \\boldsymbol{j}.\n", "\\label{_auto9} \\tag{21}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "f384ac43", "metadata": { "editable": true }, "source": [ "Note that this differs from the definition used by, e.g.,\n", "[Numpy](https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.fft.html).\n", "\n", "The inner products used in Eqs. ([10](#eq:df_var2)), ([11](#eq:kg:du_var2)) may be\n", "computed using forward FFTs. However, there is a tiny detail that deserves\n", "a comment. The regular Fourier inner product is given as" ] }, { "cell_type": "markdown", "id": "928dbba0", "metadata": { "editable": true }, "source": [ "$$\n", "\\int_{0}^{L} e^{\\imath \\underline{k}x} e^{- \\imath \\underline{l}x} dx = L\\, \\delta_{kl}\n", "$$" ] }, { "cell_type": "markdown", "id": "157f8783", "metadata": { "editable": true }, "source": [ "where a weight function is chosen as $w(x) = 1$ and $\\delta_{kl}$ equals unity\n", "for $k=l$ and zero otherwise. In Shenfun we choose instead to use a weight\n", "function $w(x)=1/L$, such that the weighted inner product integrates to\n", "unity:" ] }, { "cell_type": "markdown", "id": "6b9f5897", "metadata": { "editable": true }, "source": [ "$$\n", "\\int_{0}^{L} e^{\\imath \\underline{k}x} e^{- \\imath \\underline{l}x} \\frac{1}{L} dx = \\delta_{kl}.\n", "$$" ] }, { "cell_type": "markdown", "id": "84026af7", "metadata": { "editable": true }, "source": [ "With this weight function the (discrete) scalar product and the forward transform\n", "are the same and we obtain:" ] }, { "cell_type": "markdown", "id": "98b58094", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\left(u, v \\right) = \\boldsymbol{\\hat{u}} =\n", "\\left(\\frac{1}{N}\\right)^3\n", "\\mathcal{F}_z\\left(\\mathcal{F}_y\\left(\\mathcal{F}_x\\left(\\boldsymbol{u}\\right)\\right)\\right).\n", "\\label{_auto10} \\tag{22}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "736c335d", "metadata": { "editable": true }, "source": [ "From this we see that the variational forms ([10](#eq:df_var2)) and ([11](#eq:kg:du_var2))\n", "may be written in terms of the Fourier transformed quantities $\\hat{\\boldsymbol{u}}$ and\n", "$\\hat{\\boldsymbol{f}}$. Expanding the exact derivatives of the nabla operator, we have" ] }, { "cell_type": "markdown", "id": "e0edb86f", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "(\\nabla u, \\nabla v) =\n", "\\left((\\underline{l}^2+\\underline{m}^2+\\underline{n}^2)\\hat{u}_{lmn}\\right), \n", "\\label{_auto11} \\tag{23}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "3d2ba41c", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", "(u, v) = \\left(\\hat{u}_{lmn}\\right), \n", "\\label{_auto12} \\tag{24}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "6ecb4c2c", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", "(u|u|^2, v) = \\left(\\widehat{u|u|^2}_{lmn}\\right)\n", "\\label{_auto13} \\tag{25}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "fdc7f2d9", "metadata": { "editable": true }, "source": [ "where the indices on the right hand side run over $(l, m, n) \\in \\boldsymbol{l} \\times \\boldsymbol{m} \\times \\boldsymbol{n}$.\n", "The equations to be solved for each wavenumber can now be found directly as" ] }, { "cell_type": "markdown", "id": "19dc93b3", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{\\partial \\hat{f}_{lmn}}{\\partial t} =\n", "\\left(-(\\underline{l}^2+\\underline{m}^2+\\underline{n}^2+\\gamma)\\hat{u}_{lnm} + \\gamma \\widehat{u|u|^2}_{lnm}\\right), \\label{eq:df_var3} \\tag{26} \n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "95c68d2a", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", "\\frac{\\partial \\hat{u}_{lnm}}{\\partial t} = \\hat{f}_{lnm}. \\label{eq:kg:du_var3} \\tag{27}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "04ae7faf", "metadata": { "editable": true }, "source": [ "There is more than one way to arrive at these equations. Taking the 3D Fourier\n", "transform of both equations ([5](#eq:df)) and ([6](#eq:du)) is one obvious way.\n", "With the Python module [shenfun](https://github.com/spectralDNS/shenfun), one can work with the\n", "inner products as seen in ([10](#eq:df_var2)) and ([11](#eq:kg:du_var2)), or the Fourier\n", "transforms directly. See for example Sec. [Runge-Kutta integrator](#sec:rk) for how $(\\nabla u, \\nabla\n", "v)$ can be\n", "implemented. In short, [shenfun](https://shenfun.readthedocs.io/en/latest/shenfun.html#module-shenfun) contains all the tools required to work with\n", "the spectral Galerkin method, and we will now see how [shenfun](https://shenfun.readthedocs.io/en/latest/shenfun.html#module-shenfun) can be used to solve\n", "the Klein-Gordon equation.\n", "\n", "For completion, we note that the discretized problem to solve can be formulated\n", "with the Galerkin method as:\n", "for all $t>0$, find $(f, u) \\in W^N \\times W^N$ such that" ] }, { "cell_type": "markdown", "id": "4df82afa", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{\\partial}{\\partial t} (f, g) = -(\\nabla u, \\nabla g)\n", "-\\gamma \\left( u - u|u|^2, g \\right) \\quad \\forall \\, g \\in W^{N}, \\label{eq:dff} \\tag{28} \n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "b08ed918", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", "\\frac{\\partial }{\\partial t} (u, v) = (f, v) \\quad \\forall \\, v \\in W^N. \\label{eq:kg:duu} \\tag{29}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "a339ed0f", "metadata": { "editable": true }, "source": [ "where $u(x, y, z, 0)$ and $f(x, y, z, 0)$ are given as the initial conditions\n", "according to Eq. ([2](#eq:init))." ] }, { "cell_type": "markdown", "id": "847c4d2e", "metadata": { "editable": true }, "source": [ "## Implementation\n", "\n", "To solve the Klein-Gordon equations we need to make use of the Fourier function\n", "spaces defined in\n", "[shenfun](https://shenfun.readthedocs.io/en/latest/shenfun.html#module-shenfun), and these are found in submodule\n", "[shenfun.fourier.bases](https://shenfun.readthedocs.io/en/latest/shenfun.fourier.html#module-shenfun.fourier.bases).\n", "The triply periodic domain allows for Fourier in all three directions, and we\n", "can as such create one instance of this space using [FunctionSpace()](https://shenfun.readthedocs.io/en/latest/shenfun.forms.html#shenfun.forms.arguments.FunctionSpace) with\n", "family ``Fourier``\n", "for each direction. However, since the initial data are real, we\n", "can take advantage of Hermitian symmetries and thus make use of a\n", "real to complex class for one (but only one) of the directions, by specifying\n", "``dtype='d'``. We can only make use of the\n", "real-to-complex class for the direction that we choose to transform first with the forward\n", "FFT, and the reason is obviously that the output from a forward transform of\n", "real data is now complex. We may start implementing the solver as follows" ] }, { "cell_type": "code", "execution_count": 1, "id": "82b448d1", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:49:45.743003Z", "iopub.status.busy": "2026-06-03T08:49:45.742893Z", "iopub.status.idle": "2026-06-03T08:49:46.163285Z", "shell.execute_reply": "2026-06-03T08:49:46.162935Z" }, "tags": [ "thebe-init" ] }, "outputs": [], "source": [ "from shenfun import *\n", "import numpy as np\n", "import sympy as sp\n", "\n", "# Set size of discretization\n", "N = (32, 32, 32)\n", "\n", "# Defocusing or focusing\n", "gamma = 1\n", "\n", "rank = comm.Get_rank()\n", "\n", "# Create function spaces\n", "K0 = FunctionSpace(N[0], 'F', domain=(-2*np.pi, 2*np.pi), dtype='D')\n", "K1 = FunctionSpace(N[1], 'F', domain=(-2*np.pi, 2*np.pi), dtype='D')\n", "K2 = FunctionSpace(N[2], 'F', domain=(-2*np.pi, 2*np.pi), dtype='d')" ] }, { "cell_type": "markdown", "id": "6be0674c", "metadata": { "editable": true }, "source": [ "We now have three instances `K0`, `K1` and `K2`, corresponding to the space\n", "([13](#eq:kg:Vn)), that each can be used to solve\n", "one-dimensional problems. However, we want to solve a 3D problem, and for this\n", "we need a tensor product space, like ([14](#eq:kg:Wn)), created as a tensor\n", "product of these three spaces" ] }, { "cell_type": "code", "execution_count": 2, "id": "31bae268", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:49:46.165056Z", "iopub.status.busy": "2026-06-03T08:49:46.164880Z", "iopub.status.idle": "2026-06-03T08:49:46.184194Z", "shell.execute_reply": "2026-06-03T08:49:46.183917Z" }, "tags": [ "thebe-init" ] }, "outputs": [], "source": [ "T = TensorProductSpace(comm, (K0, K1, K2), **{'planner_effort':\n", " 'FFTW_MEASURE'})" ] }, { "cell_type": "markdown", "id": "bf3976b2", "metadata": { "editable": true }, "source": [ "Here the `planner_effort`, which is a flag used by [FFTW](http://www.fftw.org), is optional. Possibel choices are from the list\n", "(`FFTW_ESTIMATE`, `FFTW_MEASURE`, `FFTW_PATIENT`, `FFTW_EXHAUSTIVE`), and the\n", "flag determines how much effort FFTW puts in looking for an optimal algorithm\n", "for the current platform. Note that it is also possible to use FFTW [wisdom](http://www.fftw.org/fftw3_doc/Wisdom.html#Wisdom) with\n", "`shenfun`, and as such, for production, one may perform exhaustive planning once\n", "and then simply import the result of that planning later, as wisdom.\n", "\n", "The [TensorProductSpace](https://shenfun.readthedocs.io/en/latest/shenfun.html#shenfun.tensorproductspace.TensorProductSpace) instance `T` contains pretty much all we need for\n", "computing inner products or fast transforms between real and wavenumber space.\n", "However, since we are going to solve for a mixed system, it is convenient to also use the\n", "[CompositeSpace](https://shenfun.readthedocs.io/en/latest/shenfun.html#shenfun.tensorproductspace.CompositeSpace) class" ] }, { "cell_type": "code", "execution_count": 3, "id": "7eb29c91", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:49:46.185636Z", "iopub.status.busy": "2026-06-03T08:49:46.185556Z", "iopub.status.idle": "2026-06-03T08:49:46.187209Z", "shell.execute_reply": "2026-06-03T08:49:46.186994Z" }, "tags": [ "thebe-init" ] }, "outputs": [], "source": [ "TT = CompositeSpace([T, T])\n", "TV = VectorSpace(T)" ] }, { "cell_type": "markdown", "id": "5e029c38", "metadata": { "editable": true }, "source": [ "Here the space `TV` will be used to compute gradients, which\n", "explains why it is a vector.\n", "\n", "We need containers for the solution as well as intermediate work arrays for,\n", "e.g., the Runge-Kutta method. Arrays are created using\n", "[Sympy](http://www.sympy.org/en/index.html) for\n", "initialization. Below `f` is initialized to 0,\n", "whereas `u = 0.1*sp.exp(-(x**2 + y**2 + z**2))`." ] }, { "cell_type": "code", "execution_count": 4, "id": "ce384d37", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:49:46.188465Z", "iopub.status.busy": "2026-06-03T08:49:46.188396Z", "iopub.status.idle": "2026-06-03T08:49:46.340406Z", "shell.execute_reply": "2026-06-03T08:49:46.340137Z" }, "tags": [ "thebe-init" ] }, "outputs": [], "source": [ "# Use sympy to set up initial condition\n", "x, y, z = sp.symbols(\"x,y,z\", real=True)\n", "ue = 0.1*sp.exp(-(x**2 + y**2 + z**2))\n", "\n", "fu = Array(TT, buffer=(0, ue)) # Solution array in physical space\n", "f, u = fu # Split solution array by creating two views u and f\n", "dfu = Function(TT) # Array for right hand sides\n", "df, du = dfu # Split into views\n", "Tp = T.get_dealiased((1.5, 1.5, 1.5))\n", "up = Array(Tp) # Work array\n", "\n", "fu_hat = Function(TT) # Solution in spectral space\n", "fu_hat = fu.forward()\n", "f_hat, u_hat = fu_hat\n", "\n", "gradu = Array(TV) # Solution array for gradient" ] }, { "cell_type": "markdown", "id": "0f7e0745", "metadata": { "editable": true }, "source": [ "The [Array](https://shenfun.readthedocs.io/en/latest/shenfun.forms.html#shenfun.forms.arguments.Array) class is a subclass of Numpy's [ndarray](https://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.html),\n", "without much more functionality than constructors that return arrays of the\n", "correct shape according to the basis used in the construction. The\n", "[Array](https://shenfun.readthedocs.io/en/latest/shenfun.forms.html#shenfun.forms.arguments.Array) represents the left hand side of ([16](#eq:usg)),\n", "evaluated on the quadrature mesh. A different type\n", "of array is returned by the [Function](https://shenfun.readthedocs.io/en/latest/shenfun.forms.html#shenfun.forms.arguments.Function)\n", "class, that subclasses both Nympy's ndarray as well as an internal\n", "[BasisFunction](https://shenfun.readthedocs.io/en/latest/shenfun.forms.html#shenfun.forms.arguments.BasisFunction)\n", "class. An instance of the [Function](https://shenfun.readthedocs.io/en/latest/shenfun.forms.html#shenfun.forms.arguments.Function) represents the entire\n", "spectral Galerkin function ([16](#eq:usg))." ] }, { "cell_type": "markdown", "id": "cf89575c", "metadata": { "editable": true }, "source": [ "### Runge-Kutta integrator\n", "\n", "\n", "\n", "We use an explicit fourth order Runge-Kutta integrator,\n", "imported from [shenfun.utilities.integrators](https://shenfun.readthedocs.io/en/latest/shenfun.utilities.html#module-shenfun.utilities.integrators). The solver\n", "requires one function to compute nonlinear terms,\n", "and one to compute linear. But here we will make\n", "just one function that computes both, and call it\n", "`NonlinearRHS`:" ] }, { "cell_type": "code", "execution_count": 5, "id": "e982a934", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:49:46.342260Z", "iopub.status.busy": "2026-06-03T08:49:46.342163Z", "iopub.status.idle": "2026-06-03T08:49:46.354326Z", "shell.execute_reply": "2026-06-03T08:49:46.354059Z" }, "tags": [ "thebe-init" ] }, "outputs": [], "source": [ "uh = TrialFunction(T)\n", "vh = TestFunction(T)\n", "L = inner(grad(vh), -grad(uh)) + [inner(vh, -gamma*uh)]\n", "L = la.SolverDiagonal(L).mat.scale\n", "\n", "def NonlinearRHS(self, fu, fu_hat, dfu_hat, **par):\n", " global count, up\n", " dfu_hat.fill(0)\n", " f_hat, u_hat = fu_hat\n", " df_hat, du_hat = dfu_hat\n", " up = Tp.backward(u_hat, up)\n", " df_hat = Tp.forward(gamma*up**3, df_hat)\n", " df_hat += L*u_hat\n", " du_hat[:] = f_hat\n", " return dfu_hat" ] }, { "cell_type": "markdown", "id": "29c06ebf", "metadata": { "editable": true }, "source": [ "Note that `L` now is an array that represents the linear\n", "coefficients in ([27](#eq:kg:du_var3)).\n", "\n", "All that is left is to write a function that is called\n", "on each time step, which will allow us to store intermediate\n", "solutions, compute intermediate energies, and plot\n", "intermediate solutions. Since we will plot the same plot\n", "many times, we create the figure first, and then simply update\n", "the plotted arrays in the `update` function." ] }, { "cell_type": "code", "execution_count": 6, "id": "1f4ce1d4", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:49:46.355723Z", "iopub.status.busy": "2026-06-03T08:49:46.355638Z", "iopub.status.idle": "2026-06-03T08:49:46.563182Z", "shell.execute_reply": "2026-06-03T08:49:46.562932Z" }, "tags": [ "thebe-init" ] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGdCAYAAAAvwBgXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA+m0lEQVR4nO3df8glZ33//9c1c85970bv3U/iYjTNJlnT0lSCKBu7CLEmVaLBL2iUoFQCCTGwuIaE/KFuIySKdKEJbWlEa1pIpWINVlKtSYsLEqOVgEm0TQMGooQsWTWuwn2vMXvf58xc3z/mx7lmzsw5M+fHfZ1z7/MBt3vfc+acM/dZ2Xnl/X7PNcZaawUAAOBB4PsAAADA2YsgAgAAvCGIAAAAbwgiAADAG4IIAADwhiACAAC8IYgAAABvCCIAAMCbju8DGCWOY508eVJra2syxvg+HAAA0IC1VqdPn9YFF1ygIBhd81joIHLy5Ent37/f92EAAIAJnDhxQhdeeOHIfRY6iKytrUmSrtT/p47pej4aAADQRN/29AN9Oz+Pj7LQQSRrx3RMlyACAMAysWo0VsGwKgAA8IYgAgAAvCGIAAAAbwgiAADAG4IIAADwhiACAAC8IYgAAABvCCIAAMAbgggAAPCGIAIAALwhiAAAAG8IIgAAwBuCCAAA8IYgAgAAvCGIAAAAbwgiAADAG4IIAADwhiACAAC8IYgAAABvCCIAAMAbgggAAPCGIAIAALwhiAAAAG8IIgAAwBuCCAAA8IYgAgAAvCGIAAAAbwgiAADAG4IIAADwhiACAAC8IYgAAABvCCIAAMAbgggAAPCGIAIAALwhiAAAAG8IIgAAwBuCCAAA8IYgAgAAvCGIAAAAbwgiAADAG4IIAADwhiACAAC8IYgAAABvCCIAAMAbgggAAPCGIAIAALwhiAAAAG8IIgAAwBuCCAAA8IYgAgAAvNm2IHLs2DEZY3T77bdv11sCAIAFty1B5Ec/+pHuv/9+velNb9qOtwMAAEti7kHkd7/7nT7ykY/oH//xH3XuuefO++0AAMASmXsQOXLkiN773vfqXe9619h9Nzc3tbGxUfgCAAA7V2eeL/61r31NTz31lH70ox812v/YsWP6zGc+M89DAgAAC2RuFZETJ07otttu01e+8hXt2rWr0XOOHj2q9fX1/OvEiRPzOjwAALAA5lYRefLJJ/XSSy/p4MGD+bYoivTYY4/p85//vDY3NxWGYeE5q6urWl1dndchAQCABTO3IPLOd75TTz/9dGHbTTfdpMsuu0yf/OQnh0IIAAA4+8wtiKytrenyyy8vbHvVq16l17zmNUPbAQDA2YmVVQEAgDdzvWqm7NFHH93OtwMAAAuOiggAAPCGIAIAALwhiAAAAG8IIgAAwBuCCAAA8IYgAgAAvCGIAAAAbwgiAADAG4IIAADwhiACAAC8IYgAAABvCCIAAMAbgggAAPCGIAIAALwhiAAAAG8IIgAAwBuCCAAA8IYgAgAAvCGIAAAAbwgiAADAG4IIAADwhiACAAC8IYgAAABvCCIAAMAbgggAAPCGIAIAALwhiAAAAG8IIgAAwBuCCAAA8IYgAgAAvCGIAAAAbwgiAADAG4IIAADwhiACAAC8IYgAAABvCCIAAMAbgggAAPCGIAIAALwhiAAAAG8IIgAAwBuCCAAA8IYgAgAAvCGIAAAAbwgiAADAG4IIAADwhiACAAC8IYgAAABvOr4PAMDZxYThXF7XRtFcXhfAfBFEAMzcvMLGNO9JUAEWE0EEwNR8BI+23GMklACLY64zIseOHdNb3/pWra2t6bWvfa3e//7369lnn53nWwLYJiYM869l4x77Mh4/sJPMNYh873vf05EjR/T444/r+PHj6vf7uuaaa/Tyyy/P820BzMlOPXnv1N8LWAZzbc3813/9V+HnBx54QK997Wv15JNP6s/+7M/m+dYAZuBsPDHTwgG217bOiKyvr0uSzjvvvMrHNzc3tbm5mf+8sbGxLccFoOhsDCBVss+BQALMz7YFEWut7rjjDl155ZW6/PLLK/c5duyYPvOZz2zXIQEo2e4AYrrT/xNke/0ZHMloVEmA+THWWrsdb3TkyBE9/PDD+sEPfqALL7ywcp+qisj+/ft1lblOHdPdjsMEziqzDh6zCBbzNsvgQigBqvVtT4/ah7S+vq49e/aM3Hdb/tW49dZb9a1vfUuPPfZYbQiRpNXVVa2urm7HIQFnrWnDx6zDxjTHM0kQqDv+SQJK+dgJJkB7cw0i1lrdeuuteuihh/Too4/qwIED83w7ABWmOdFPGjq2q8XT5n3GhYRZBBRaOEB7cw0iR44c0Ve/+lV985vf1Nramn75y19Kkvbu3avdu3fP862Bs952BpCpgkdnin+G+pOFBKl5UMg+i7YVEwZdgWbmOiNijKnc/sADD+jGG28c+/yNjQ3t3buXGRGghe0III3eY5qAMWsNAkvTwDDpjAmBBGeThZkR2aY5WACpSUPIqADS6jUnCB/TzJw0DgXucdWEkqYVEyokwGwt0H+yAJiU1wAyInzM+yqaJq8/FBgahBJp/LzHNIGEMAIMEESAJTdJCJk6gMwyfAQT3mkijhvt5h7PtKFkVoGE6ggwQBABltS2BpAxLZfa1500ZDTR5LVLYWVkaGjRvpllICGM4GxHEAGWUNsQMo8AMrPwMelQa5MrZtxjcULJyCpJ+Zgq3mdUgGgbSKiO4GxHEAGWyNIGkEnCRpi+XlTTghn3muUAMW0oKb3euAAxSSAhjOBsRBABlkSbEDLrANI6fDQJHmHDysm4/ZoGFTdITBJKpggkVEeAegQRYMHNqgpS+zpNw0eb0FEXHtpURrL3GzeU6r7kyHZNTTAp/17p+9WGkgkCCdURoB5BBFhgswghEweQSYNH1eNN50Y6Vcda2tYfcYJeWRl8PxRgyiGg4jj7/cpgUhkkxgQSaTiUtK2OEEZwNiCIAAtq2lZMmwBSW/2oq2CUKx7ufu6JvDJYFNlO2+HWZH/TH3f5bvreWXAZ+l1KgSCKq9s5TmWmTSCRqsMErRqgiCACLKC5hJBpAsio4OE+T8rDx9iAMaJKUn5uVeiwK9VzHmWDG02ExWpK+cqY8u/oBpMskJTaNpWBJNs/e/+KMEGrBhggiAALZpoQMvTcSeY/sudMED4kJ0TUBI1RAcXWPccJHaYydNRXSaxz2Hko6UfVv3P+eDmYlKoeowKJ+3qlQEKrBhhGEAEWxMyrIE0qINJwFSSsqYqUQ0JV+HD2KweOupBRte8odlToWAkqgkqQ75uFksLtOMuhJI6rqyVZhaTUrpGaBxLCCDCMIAIsgHmHkEYBRBoOIROGj3LoqA0aQfUduuPO8PagX7yJpl0JpXj4xpqVQSULIP3seAYPNQ8lWVgotWuy/TQmkDhhRJq8VUMYwU5DEAE8m2crZmQAcfcfFUBKA6dNAsigPdM8aIzaX5LiFVMRPJL9yyFFKgYVm1ZFBhWT4SrJ4NVSbiiJnWpI3rKpnh+RRsyQzKhVQxjBTkIQATxqGkJmUgVpEkDc/cZUP9ztdQGkMnBUhI24omJiw2Q/Ew2HjKDUlqkOKZJk8pCSH1upsJFVSRTH+T6mHxcDWP4cJ3yMaddIFYFkhq0arqjBTkEQATyZWQiZdRum7qqXhtWPPHy0DByFbYH7fRpInOwRpb9/o5CS/Xp5W2ZMIFGpSpKFkqxC0rJdI1WEC1o1QI4gAmyzbW3FTNiGqbrypW0AKYeOcuCww5lEtq5lIymLHMZpw9jAFAKKlISUckAJ+rFTNUmqJFnrZrhtU1MlyV7MeazQrhlxua80OoxItGpw9iKIANvISyum6nLccVWQCasfWfjIQkeTsGGrKicVH1OQnmftSlYhyeY/BrKQUq6ixOkAa1YtaV0lSQPJUHVEGq6OtLm6ZsatGsIIlhFBBFgwcw0hLaogthNMFECy8GE7ZihkVAUMdx874l+kWMX2iWTyMJIJSu+XPW76SShxqyVtqiSmr5rqSKldM2l1xD1mwgjOMgQRYJu0vW9M7fNmFUJqFiHLAsgk1Q83fMThcLWjHDTcx8shxTo/m+zc2h1URtJHCmGkLqwEQfKn6du8ndOqStJpUB2pmh0pV0fqVoCdYZsGWDYEEWCB1N05t9H+5XmQCVsxeQgJTOPWSxZAsjBhAyPbaR403BZOXDG8mgSQJHDEXZXmQswgqJQ+BinZ16YVkSyQWLWoknQCJ4zEeaCobdXUVUcqlomfx8wIVREsG4IIsA2aVEPatmQahZAJWjGFEDKm+pF8bwrVjyyAxGHDsFHYPrx/PufRNcVqiBNGgsCpipSCiomMgqAUSKJk1qRcJakLJUFgZONAJoina9VsUxgBlglBBFgAMwkhs2jFOAFkXPVDGrRfspaLG0Bs4ASOMWEjeR/n+1JRxKQ5I+6qUP1oGlJsWBVIrBSaUpUkqYQkLzFo3cSdpGVTrI5k7zWmVeMhjFAVwTIhiABzNq4aMul9YyS1DiFNqyBxJxgaPJXUKoDE3cFhjgschW2loJKFjTyAdEphZExIUWCkuDqQmL4UKw0kUfLaWSBxZ0nqw0iDVs0kYaSEYIGdjCACeDTVFTLbEELK7Zfs+3EBxAaDQOFWN0YFjrptUhoq3KDRrQgokozzXllAycPMmEASBcNtm2yWxGYtmyCerFXjDrE2CSOlqkgVwgt2CoIIMEejqiGNQkjdc5qEkKoVUse0YqYJIAqkqJtuM04QaRg43G0qtWZs9j9OGMkrHmMCijFJKBkfSKzUN2kgsQqyg+hISls2o6sj6TGMvapmsjDCvAh2KoII4EHjEFI1FzIqhMyoClI1gCrVB5A4TFojWQDJHpecQFLVnqmoluQzJNkAajzYKYiSYzGRE04kmXL1RMoDignTP6cJJI3DyJhWjRtGkr+ouYYRqiJYBgQRYE7qqhvLFELcS3DrAogNkjU43ABSFURUMQ9SDh1VoUWyg6HUbHtXeUBxw4kkWWeNM5O1XoJpA8ng0t+2rZqxQ6wVYWTo/weEEexgBBFgDhYphLQJIJKGQsioAGJDJ4QExT9l2oUOG2bljeLPJjKl1kxxifcoa82Uwkn2Wta2DySSFMjIrgwCSay21ZH09+jkhzI6jEjJz3HMZb04qxBEgBlrdVO7BQohbQOI24YpDKim2+NQrUOHyo9Lsl2bhw/Fku3YYjiRZIJiOHErJ0HPtAskMpKy6ojS8DFpGCm1aqTRl/dmf9czDiNURbDICCLANmm0amrdZbqu8mCqY+RS7Q2MCyFVcyBVVZC4a53txdAh1QcPE1opsMmu7oyIU+2wkRkfTpzKidzLc40G4xmuSFI2o9JVZXVkqjDiKrdpCo81WA7e0SZgEEawqAgiwAxNXQ1xH6+rhrgqFisbWuNcajWYKjkhpGtqqyBVISRakRTYJLB0smrIiNAhFYKHyU782T5BLBsHitNQYWMj01FlOHEXTXPDSRRKQU+KZfJdqqojwZYGYSQ93DitjiSHP00YUfXMSPb3Wq6KpMZVRYCdgCACLIpJqyEVl+lO05KJVsaHkLpWjO3a9LVt8v2I0CGpEDxMEOffJ79qdlKOFBXCSPJ7FsJJUBFMpCxNDIJFkAQKExXnVyQpXhndqklebhBGFBXv/DtKEkBUDCPS6KpIqu3sB4OrWDYEEWBGJr27buVrlds4o6ohajAX4mgTQtq2YuKuTSohgaSVOA8abvBwqx1SEjyy0JG1QjphpDD9PkrDRT8KFcdGWTAJu0kYqQwmTtXE9tLHjWT6aaho26pJ7+KbBJokjMSxGldFTFwRPKqqIuW/3wmvogGWCUEE2AblYNF6Cfe6fUstmdp5EKcakplFCCm3YmyYVEJMN85bLW2DRyeI88e7YaReFCqKA612IvXjoBhMguFgIqnQzpEkhUZ2K/lsJmnVJIu1GYUahBEbaOJ5kdqqSEV7RprtFTEEFywagggwA1NVQ+qWcC8LnZBRXjnVrY40bMkMXSFTE0Libk0VpKIVk1VBTCdW2I2G2ixNgkeYLhqyEiYny92dnraiUFF6yc34YCJFsR1UTEKrKKuKhEamZ1q1agKlBZJeMYyYjobmRcxWdbumdl4kOfDB/w+ySseUq65yKS+WCUEEmFLbENJ42feqIVX3cl01nwuRKloyTgiJu+l9YpwQEq+0b8WoEyvsxjJBrG43qg0eUhJOumnYCE2cB4+VIFIniLWSLgayFYc6p9NTPw60FYd5UIlsoF6UfJZRHFS3cgKjKL+k10r9QFbtWjVZaCmHkUjJn1kYCWOllycPV0UUO5cjO/Mi6YeTvOeYqsiQKQZXqYpgkRBEgDkbedlu05ZMWPq5oiUz6lJdG5pCCIlXzHAIyQZQa0KI7QzWBqlrxYTdpBLS7UR5+GgaPAbfR1oJkhNs18Tq2UBbcUf9ONQ56mkrDtVP50K24nCoWiIVg8lmr6MwiNXrhcpPvS1aNem8ayGM5OuhRIPhVRuo9byIdS/vzW6O17QqUtI2XBBGsCgIIsAU5lYNyYyrhlSFj7qWjDucWhNC3FVSs5ZM21ZMtxNptdvPWy7ldktW6XCrHln46Kb7rQZ9dYLBSbIfh9qMO+rZQOdIeTDpxHGhWlLVxukEsTb7g8+9rlVjAkm94VZN3C2GkTi9kiboWUUrUriVhJFIRmFPlS2auGNGz4tIgxaNNFwVmeIuvcCiI4gAvlStoCo1qoZMMheShZB4ZXQIyWdCJmjFrHT7CgOr1U5fuzq9xlUPN3isBn3tCnrqpjeP6dlQZ+KuXqVNbaYBZNPEUqjKaonb0glNnIeS5KO16qWhpNCqSUNIvm5IqVWThREbZVfRSFVhJLuSxgbFqki5NVOeFyn83c8oXDSZE6EqgkVAEAEm1KQa0mg1VZcbQqqqIZ3ie467VDfZpzicmq2aOiqExF01uiqm3IpZ7UTa1e0pNLFe3d1qVPXIgockdU2kVSeE7DI9nbFd9WyozbirtXAQTDbj5PNxqyUdZ1tWLfl9P0kOYRBrsz8Yns1bNaXZkapWjfrp57SSXE1TF0biFUlbVloJFGzFajovYjuqrorICSUNqiIECywjggiwTQrBpa4a0oDbkrG1q6hWtGTS4VQbJgOXCpz5D/cr0ExCyDmdXuPKx6oTRHYZJ4gEPXXjSD0bFiokXRNpV5B8nmfirjpxpH4cSuFWEkqCJJT047DyX7nNXkfdbvJ6US+U7aRtEgVJq6Zr85VVg+SBpDXTSytFURpMJCk2CmQVh4PLehVZ2Y4Z26IxQ4uZpYOrVbMic0J4gW8EEWCRldsyFbJqyMh9nDvrDj2W3XguG8DM2jL594Ml2k2+JkiyEmqYD6Ima4NkQ6idIJ64ApKEjJ662e1rswENl7ut6iPKyw+JfhBoy4RSmAyyRmGkrbijILSysc2vtLFhkjhMOnxqAitFpnSfHOUzJMadYnXfPlTlVTRuNaSssLZInYaDq8AyIYgAi6jclmkgmw1xWadVk4WRLHDEpcDhPiZJCmzhbrk2HCzPHoSDhco6afDoOvMgK04IWU0DgRtCsuCRhZBdxqmIpCEk25Yci7RLPZ0Z9EQGASRWbRjphlt6WSvqBKFWwkiv9IMkMAWBgsAqipOF1kyQznVk96zJ7wqcnPfDNHjYOAkfNkzmO6wTQvJKSUcyWxXHU2Poct5JMLSKJUYQARZFVZvFfaxTUc6o4LZlpOoqiHviLocRmdJzjLN/ulqqcRYoy6ohkgrVELcVk4WRRq2YNIRkQSSbEzlju6OrIzUfX88G6gdhZVUkjo3iKFQQWkWxTSo+6SJospLJKyWD4GGDtCLiFEri0CiMbV4JkZLgl94PuFJte2ZKtFqwbAgiwARmdV+ZSedD2rLp27irh7rVj1FtmcKN65wl28MgLlRDQhMXqiGTtGLc8LFiak6mLVs13TgJRnVVkSSEJG0YG6SVnzhJGePaM1IS2oJ4EFCSz7PmpniBqRxaHYuKB3Ywggiw3Rq2WkYZt4BZvp/TqsmqHNkKqpLyZdyTfTU42TptGSkJJCathkjJ/WLGVUMmacVkAcR9vFdee11q1arpB+FQVaSnZH2RvCoSmPRymUFVRFEaxGJT256xdnhOxK2E2I4ZWuBsanU3w6vAHAmWAUEEWDRtg8qItUOkYhixFS2Z8vfu7EjWlsmGVCXlQ6rjqiGNr4qpCCHJftUVkbatms2gM1QViWyQr76aVUXqhlaT17SyoVE2P1usfmhoTiT7HMNpckenI6nfbLl3YIkRRIBtMLaVM0G7pQ03jJTnQ6raMsnPg7ZM8rzBkGqTaogbQpq2YtwQ0jXZeyfPzaojZ6xTBcl+nxGtmtWgX1kVCYNY/ThQx6mKVA2tuqHD/dmGyQJnRmm1KU7mRKT2A6v5Z56vsqpGVQ8qHtgJ5vuvX+oLX/iCDhw4oF27dungwYP6/ve/vx1vCyy3BpfujhPXVD2ye8q4Pw+1ZfIrR2xhSNWthkiqrYa4ISSreFS1YvYEZ7QnODMUQnalX2tBrLWgn26PCqEleU5SeVkLzxQqLqtBT6tBT7uCXiEgdYLkEuOu8zsE7pVAWXtGg98/W+I++xzrqkk2KH3mQXaPH8mG6Q0InQpWcUG6MX/fM2jp1ZnVzBMwibkHkQcffFC333677rzzTv34xz/W29/+dl177bV64YUX5v3WgFfjBlFbr7o6Ieu8jTsfIg2fVLO2zNDaIRXVEEnqhlFtNcQNId0sRIxoxawFW1oL+k4ICZwvm24fDiP/L/i91oJX8kBSfs8klCThaCXo56EpNHE6tDoYvk0u5bVpK2ZQEbKB0m3OZxcUP7vyv6Z2zNouMzfHoALM09yDyN/8zd/o5ptv1kc/+lH9yZ/8if7u7/5O+/fv1xe/+MV5vzWwfMonE2dp9/KKquPuLVMeVC3Ph5T/q96md9YdbLOV1ZBMvoBZTTXEDQRJxeIVrQWv1M6DZAHk/wVGa6ajrgK92qyoq+owspZWUqqqI9nju0yzqkgnLM2jZFURt1KUfo5Zeyb/nMo/h9sQQubcygO201wj9NbWlp588kl96lOfKmy/5ppr9MMf/nCebw3MzTKVsYdOiMHwZbu2dLLNfy4NqWbcIdVR1ZC6Vkz5qphiKyZQV4FWTUdd4yalfnqMcRJIrFXPjjjZO1fQuFWRulkR91Jed2hVUr6miI2NTLZGSHr1jI2Ll/G6cyLZtqkGVqXhO/G2xBwJFt1cg8ipU6cURZHOP//8wvbzzz9fv/zlL4f239zc1ObmZv7zxsbGPA8P2H4eyudV/3VedQ6PS1WSXM2Q6rhqSFUIcdWFkLI8kKRh5IyNh16nSi/tSe0yPfWCMLlx3ograKKKodV8TRHns3Cvnsk/z2y5d3dbOrBqA5P8S+texrvFgmNAZlvqe8aUlp22dmibJB07dkx79+7Nv/bv378dhwdsHw+LUpmKhbNM1VpbkXN1iCs2iqP0KzaKYpOcuG2grThMvzraTL96Nsy/zsRd9WxncLmtI7kCJlTPGp2xRmdsrN7Q5S9Sz0batH31FA+FkMHrdLWV/ll+nzO2q824qzNxV/04rYjEofpxoK0oVC9Kvu9HoaI4kI0D2dhIsZGJSv9OxabyM6rclv5Vm9jK9K1MLJnIzmYtEWAHmWsQ2bdvn8IwHKp+vPTSS0NVEkk6evSo1tfX868TJ07M8/CAiSzT8tlDISQu/pe7yX62xW2yyb62fCKW1E9P3L0o1FZ2EndO8Gfi9MRvu5VhZMuG2oh3DR6vCCObtq/fxZtDIeSMNTodBzodd9LnDYeQnu3odLQ7fd9BKMqC0lbcyQNUXg2JTbrc+/DvayKTfrnbBl+ZIJKC0mqqwQL8X4W2DBbdXOvEKysrOnjwoI4fP67rrrsu3378+HG9733vG9p/dXVVq6ur8zwkYLH1+8X2TfZzOiJhO+l/aadrXdlOoKCf3Ww2VpxNV/Ztsk+YnFiDdN0LpV0OE0nGWVXVRMnTFBopW9I9W0sjNLKxkY0DRbFVkJ5d3apIJ47ViTvaNLE6caRdaSskly4w1rMdnTHdoTZNInndM9Zql+lrlwnUs1vptiSEJIGlGECSxwchJAsg2fam1ZD8UKOkGmIjMxTSgkjpPWgG20z556i6CjWxqvmQGd+fBvBp7g3rO+64QzfccIOuuOIKve1tb9P999+vF154QYcPH573WwNe2V5/5CW64x6fFdOX7EryfRBZKTCK0lqoiSTjDLAGkRQHSmYkOjavipggqRaYIJmhiGKjTiD1ojCfFenHoXrBoCrSNVGh5tqLwqFtrp4NSyuqDk62dSHEDSCS8hCSVVuy1x1VDclUtWWy1oyJNdSWMXGpJVPKBjMNI01wLxosqbn/K/ihD31Iv/nNb/TZz35Wv/jFL3T55ZfrkUce0cUXXzzvtwaWWxRPvahZEKkQOrLAYeLkv+SzU6WJk/umSEpOuLHNfzaRGaqK9KMwH1odVRWpvCtuadsZ29WWCUs3uRt8nwWQ5PvqEJLdc8YNIZvptlHVEEmFtkyhGpJ9NmlbJmuzlFsy5XASFH6umA+JrYK+lWIr049l4limydzIHIPGMrUbsfNsywj/xz72MX3sYx/bjrcCFpKNotGX/cbxXNeGMLHN2zSKVWzRhM7dZI37WHYH2jSRxMNVkSgO1A0jbUVhZVWkfN+Xwe8rnVFXvaCjbvkSlFzWqmnWismGVLMQklVERlVD6oZU88+p1HZxfzbRYOi3PB9S+yuNkgWSftS49cL8B3YCluIDFk15TmSc2NbOidiVtLohp8XghJBs1VW3WpK1Z2ygwtBqdufdclWkF4VSqMqqSPKCGrQtxlRHXFmrpm0rJgshWZWkVTVE2edRumKmoi1T+X1peHUqfW54h7MDQQTYbm2DRgXTjwsDq7X7OZWQIErudJ9/H6SDlmkIqWzPuNWRtD2j0CYn7gZVkZyzwFibVk1W5ahqxbgDqW4rJgshm3HyGVdVQyQNVUOy3zsfUtXotkzdfEh5NmTml+22GFSlYoJlQBABJjC21TInph8nVYyWbZxsYNVEyQ3e4iA52eZBIx5uzySDrIP2jDu0miz2GRSqImEnHqqKDJmwVdOmFZOFkCwIlashkiov2R1aOyRrw8SDn/PP0wkjeUipCCOmX758ejAf0gqDqNjBCCKAR4UrZ0bNicRxcmv4zvjw47ZnrIoDq4PXU7FFE1a3Z5KDdPaPB0OrSlciLSxwFqlQFVFc8U9My1ZN21ZMFkKyiki5GlJewEzS0JBquTVjnLCRV0sKl/Xa/HPLnzMmbGSDqrPG4CmWDUEEWERu+6ZhKyfo2+S87yzpbtI5EWnQpsnmRII0eNjO+PbMYE0RWxhaDQONrIpkwyubcUergfNf9S1aNW1bMVkI6eUtmPpqiKTaIVVJQ22Z/HOtqoJk21oWL0w/nn5dEComWGIEEWCRNbiEN2nXBIUAMrRPOrDqzonkj6WX8Y5qzxTXFLFDC5zVVUUkJWFEGl0dqWnV7Ap6rVsxWQjZyisi1dUQSSOHVKvaMpUzIlUdKLeC4ly2O0qjS3idwML8B3YKggiwTQpzJU6Vo+3CZtmcSLbCqi23c2JbaM+Y2CqQke1L6liZyCiUVdR1hjLT861JX8paKdySohWjoGcVO7fntQpkunF6Ah9uFe3q9vS7XrKCWicI1Q/C9PtIPRuoG8f5ts0gqZRkoaJrIvWCML+RXdW6IOMqIFtpKNlKg8nv+1290u+qF4Xa7He02esojo16vVBRL0yqIf1A2gpkekayUtAzMj2T33/HRFLQGyxilq2mml22m31vYls/HyINrR8yJLt0161wzLnaQSsHvhFEgAk1GVhtvXqqOydS1Z4pzYmUr55JOiHFyohJl3s3saS+URCkRYggvRtsN+2gqFigiLtJGLGhkZRUQdIsIKtAio2iTnIyjSOjuJu8bz8OtNpJTp4rYaTfq6uV9G63v9dK+n1yx15FUtfEWg366qRlhNWgn961N9LvtKtR5SPbtpW2YbLvt6KkEnKm19VmP1Q/CtXrJwEkjoxsP5DtBUklpGdk+kZBT1LshJD+YDYkCyFBLw0hvSSEhFtJCAm2bLJ9K65dxKz4dxMX1w9xlS/drauG0JbBkiOIAL7UVUXKQ6vl9kw/yu874+6XtWjq1hQJAqNYVsGWpJUsrNSHERMMWjWBTHJY2WORke3a9PCMTBp+ojhQtzM4oZ7pd5NLe00SFFbCSP0gkNRVJ0iWhpeUB5OuSaolL2s1DybS+KpHPw7y77PwISXLykdxoM1+qK1eJ2kf9UJFvUDqB0mrKQ0gJkpCSNBL1wxJL202cVoNiZLPKJsbKYcQ03erIvUtmcJqquWqyJTVELe60aR1QzUEi4AgAkyh7WW8rS/7raqKpOGjtkVTs8BZ1qKJlbYSVB9GTJicfG3sXEkjKZZRGGlQGYmDPJD0Y6OwOzixRrFRGNj8ypRsOfhX1E3uT9OwWtKk6iEpnwFJ3jvI9xnViglqWjGjQkjYU34XYzeEhD2bL+leNmpJd6ohONsRRIA5G9memaIqkrdo4njojrxVshZNICNtWcUrpTASmCSEWCnYkuIV9ziTXcKoulWjKJRW4vwOMVmrJgisNiV1wkhhEGiz31EniBUGcV4tcYNJVbWkXPWoCx7ZPlE6YJoNpLZtxcgOwkfWjimHkLBn87vsJq2bQQgZd1+ZXPkeM3XVkIYDqm2rG1RDsCgIIsCUZlkVqQwjY6sig++H50WGWzRZGMmGVwMZ2TBpNSgNJnE3CSM2LF4UUteqkUnnRlRs1QShTX+V9OdCMLGVwSRMJ2hXwsGJsqrdIlUHjzi9gkdKApGNg9atmCYhJOwNQkjdVTKj5kJy5XvLNF3WfYpqCCEEi4QgAvjWdMl3tyriLHBWtdrquHmR7B406hupYxVuGUUrw2FEURI2shZN3qqxzVo1JrAKQqtIoUyQnGCD0LYKJm2DR3Yc2c+2nwzWNm3FNAkhQT7EagdzIRUtGffvozwXMtSSaVMNKYWQcrDg0l4sE4IIMAPjqiLl9szQ/k1bNKV961o0Q5f0ahBGbFCcF2kSRuqMa9XYwCrKLmcNspDQLpi4wSPbr1zxyL5PXj9bmCz9eYJWTH6JbpStJ1IKIWkrxg0h41oyQ3MhmQbVkFkGC6ohWDQEEWBGZnn/maG5kqqqiRNS8hZNxbxIVhXJFjxLLuMdtGjGhREbJOuK2KhZq8ZGyQJoVoFs2poxoU3WBAuson7ypxTIBFZRT4pCWxlMMmODRzl8ZIuSpTewa9OKKYeQ5FLdUghJc8EkcyGDv9OaakjVvq4x1RBg2RBEgEXRpEWTtWcqqiLTzItUhRG5MyBSsnZZg1aNCQYLoJmeScJIz+Srp1aGk3g4mEiSCeJim6UUOKqCh6TC3XObtmKCrVIAidLPILL5uiHucGrQ4vw/1JLJVC5qVgoaLashsxxoBbYDQQSYoTZVkcbtnKqh1Uz+WLFFU3/zPFsbRtwb45lICmWTkBEnISKS8jBSp9yqyUJK8qLZMaR/Okur29DWVk2S55ja0KFY+Y3qkh1UeLx2gbJRVZAGIaRVS6YcQqqqIRVtmaFQweW62IEIIsA2abTK6qRVkdT4Fs3o+9FkVZGkSpIIlFQ4wl5yPm3TqjFBUg1JF0dNbqzXN3lVpBBOKqom+bE1DB3J75FtS49/0lZMNrA6bQgpy0LIqGpIw5vgtalwUA3BoiKIADM2dVWkanC1XBUpX0HjVEWazIuMa9EYWZlOUn2Iw+TE61ZHskASq75VM6iOSOopeR0pqXSMCSd5S0fDocPdVg4dWdWkcMfcpgHEqYLMNIRUXSWTBY2qakj62LhqSFWw4GoZLCOCCDAHdWGkqioyVRhReuLJ1hapCSNGyayFDQIFcXqn3jhdxKwTpAHCynSSdkocqlAdUaTKQGLi5MqarDpSFUgy2bIgbcKJq1HoyD/UwbZpAoik2bdjqkJIqRIy6xBCNQSLjCACzMm2hRG3TVMTRhQEEw+wGlnZbEXWikASRFLULbZr8iDSUyFcSMmJPds2Npxk6kJHoUqS/hkXt01aAUleq3kASfafPoQMIYRghyOIAB40DiNVzxkXRqTSYmeBc2JUqzCiyCbBYCutdIwIJHF3eH5ESkJJ/jsE1dsK4aSvNImYfOYk2TH5w618DFVJVBFEpggg2fcTt2Kk6hASlQLJhAuXDe1ffowQgiVAEAHmqPFy7nX7lwZSG4eRbOVVDeZBJ2nV2I5RGEs2vXOvG0jCuDhDEkRS3DWFgdbkd3J+v+xXcyol0ujqidQscAw9XmrPTBJAkmNo14pJnjP7EFKFmRDsBAQRwKO5h5FZtGokKQseQVolyb7v2OFAEgwCSX7M2aXBTuVDfclmc6gVASUPIw0Dx9D29HlNhlCz/dzwkb3GuIXKatcImXEI4aZ22KkIIsCctV3+vfI5dZf1ThNGNL5VY9IrX8xW0oJJ5kXMoG0zIpCo6wybBsVDzn/PMQHFZUa0ZqTB/EjyJu72dgFkUH1JH59FK0ZqHkIqcIUMdjKCCLAA2oaR2rv0SkrKFzVhRO1bNepLcaccSty2jfJA4l72ayI7VN2QknU64rAipKSDqvnv7y6w5lY3MqXZziCyQ/tL9ZfhJvvZQvUj2d8JINL0rRipXQiZwWW6VEOwTAgiwDZosrZI2wHW2pvj9dMyR+iEk6yaMkGrRhqclOM8NRSrJO5lv9kcSRZIMm7ICHp2eEZEFSHF4QYNqbplM9jX3W94HRBpRPVDSkKY1L4K4jzmK4QAy4YgAiyQsauv1s2LSIMTXF4dSbePadVIyfCpVNGq2UqrIrKKO0bBVjSySpIEErdtMzj0wHkfqRhSpEFQKYcUqRg0yvd4yYJF/nN/+HH3Etxk25jw4fzcuhUjeQ0hVEOwbAgiwDaZ9O68jYdXM23mRqSxC6BJLUJJerWNpLxK4ioulWGT0JIqBJVeukdgCkGjKmS4qm5E515+K7ULIMnzp2jFuNtGLdtOJQRnMYIIsGCmupImMyqMSO1aNVJ+KXB2d98gO7dWtG7KVRJXEjbcSkgxjJTnROSsWDauEiINrnYZ7JP+OUH4cH+eqhXj7pP9Vi3X/mgaQqiGYBkRRIBt1LQqMvMwImlsq0buIGvyZ76+hpKQYrbi1lUSSYOVUyWZreSkn1+e64SNbAB2lLqwUdhWmieZNoAUttVVQaSxIWSeMyGEECwrggiwzbITRpPhVUmFcDH03CZhRKq4R41UuE9Nvk/1Zb6S0xZpWSVx9407g15NXUgpt3PyxxoEjqH3ldqHD6lZAJGmq4IQQgBJBBHAm5lVR0rtlqoAM1QdqZsdkYZWZJU0vPaINDaQZFUSV14xcfVHB5RGYSP/Pav2bV/9KGyXml+WW/E6UvsQ0mYehBCCZUcQATza9laNVN2ucasjNcOs+etmbZutdJakpm0juaGkcOTVAUUqhpS6c/GosDFi31btF2n0HIg08oqYDFUQYDyCCOBZmzAitW/VlJ8zUbsmU7jkd3yVJN+vJEjnRIaqI+m8SDmoSM3ChssNFROFj9J+Uw+jEkKASgQRYAG0ubS3basme440RbtGGgolxdZNuq003CqllRMpr5ZU/k6lsFIIKg2CxtBjdW2WNtWP/LEZXA1DKwaoRRABFsRMw4hUeX+aoUDSpF0jDYeS0sPO+qmDAFJe80NOJaWkSVjJ9x2xHkdlOJkmfGTGtGEmDSAjn1uBEIKdiCACLJCmV9RILVo1mVEVklHtGnf9Efe1Cm2LwbflS4DLstmSysdqwsqo6kfh+Ktes0n4qLjbbeWiZKXnNQ0gEiEEqEMQARbQzKsjUrOWTXkoszA/kqoLJvn7JH+Y9FLgau7wa32FZJyxAUWqvuxWqq961O3jKYDUvQawUxBEgAXVNoxIY6ojmaaBJKgIHMkTiz+Wg4kbSkZUPzLGuTqnkfJAaZ1pg0fF6zQZQnVRBQHGI4gAC6zt/WnqqiNSu0Aiafjme4W7+7pK+0VbSTip0hnxT065CuIGmFH3aSmrCwbl4FG3X5PLcEc9X/UBghACDCOIAAuuzdyIVB1GRr7OmDkSKQ0l5RN0XTCpaudkoq3B93VhZRJ17+dqETykdu0XaXRwoBUD1COIAEti2laN+zpSTbCpqJK4r1d4zfKlveXXyNRVQcaFhyyoNAkZVUZULFqHjxGvNy40UAUBRiOIAEtkkuqINGEgqQkQIxdJk4ZnQkYFAmnyoNLktV2ThI8R7zHLANLk9YCdiiACLKFJZkekloGkbgC19JpDrzsqlFRpEybamnH4yJ9PGwaYGYIIsKTaVkekGbRspMlDSZ0mYaVOmyFWTdZ6KTx/hgFk3OsBZwuCCLDk2lZHpCkCiTR2BmTUa1dqGSbamjZ8SAQQYJ4IIsAOMEl1RJoykGQaXHUzSuPAUmGSEDCL8DHxezd4XeBsQxABdpBJqiNSs0CSGfn6Y1o4o957rqZsuxT2I4AAM0UQAXaYSasjUrO2StUJdexcSZ2GYaXWFIOubYPBNIGJEALUm+GKQkXPP/+8br75Zh04cEC7d+/WpZdeqrvuuktbW1vjnwxgajaKJj4B2l6/9foX7ldj/f50X21+pwmPse1nUfV+AOrNrSLy05/+VHEc60tf+pL+8A//UP/3f/+nW265RS+//LLuvffeeb0tgJJWrZXycytOwE1mOpqefCep2kzyPo1eawYtIkIH0J6x1trterN77rlHX/ziF/Xzn/+80f4bGxvau3evrjLXqWO6cz464OwzbRAovNYUQ6fbaZYzKQQPoFrf9vSofUjr6+vas2fPyH239V+O9fV1nXfeebWPb25uanNzM/95Y2NjOw4LOGtNUy0Zeq2WJ/hZBpdtGXjN3ovwAczU3GZEyn72s5/pvvvu0+HDh2v3OXbsmPbu3Zt/7d+/f7sOD4Cmmytp/V7p7MUsvrbleJn3AOaidRC5++67ZYwZ+fXEE08UnnPy5Em95z3v0fXXX6+PfvSjta999OhRra+v518nTpxo/xsBmBon3QQDp8D8tZ4ROXXqlE6dOjVyn0suuUS7du2SlISQq6++WocOHdI///M/K2ixnDMzIsDimOU8ySIjdADTm+uMyL59+7Rv375G+7744ou6+uqrdfDgQT3wwAOtQgiAxTLLeZJFQ/gA/JnbsOrJkyd11VVX6aKLLtK9996rX//61/ljr3vd6+b1tgC2iXvyXrZQQvAAFsfcgsh3vvMdPffcc3ruued04YUXFh7bxiuGAWyDZaiWED6AxTS3XsmNN94oa23lF4CdrbyK6XYNffp6XwCTW44ViADsGIQCAC6mRwEAgDcEEQAA4A1BBAAAeEMQAQAA3hBEAACANwQRAADgDUEEAAB4QxABAADeEEQAAIA3BBEAAOANQQQAAHhDEAEAAN4QRAAAgDcEEQAA4A1BBAAAeEMQAQAA3hBEAACANwQRAADgDUEEAAB4QxABAADeEEQAAIA3BBEAAOANQQQAAHhDEAEAAN4QRAAAgDcEEQAA4A1BBAAAeEMQAQAA3hBEAACANwQRAADgDUEEAAB4QxABAADeEEQAAIA3BBEAAOANQQQAAHhDEAEAAN4QRAAAgDcEEQAA4A1BBAAAeEMQAQAA3hBEAACANwQRAADgDUEEAAB4QxABAADeEEQAAIA3BBEAAOANQQQAAHhDEAEAAN4QRAAAgDfbEkQ2Nzf15je/WcYY/eQnP9mOtwQAAEtgW4LIJz7xCV1wwQXb8VYAAGCJzD2I/Od//qe+853v6N577533WwEAgCXTmeeL/+pXv9Itt9yif//3f9c555wzdv/NzU1tbm7mP29sbMzz8AAAgGdzq4hYa3XjjTfq8OHDuuKKKxo959ixY9q7d2/+tX///nkdHgAAWACtg8jdd98tY8zIryeeeEL33XefNjY2dPTo0cavffToUa2vr+dfJ06caHt4AABgiRhrrW3zhFOnTunUqVMj97nkkkv04Q9/WP/xH/8hY0y+PYoihWGoj3zkI/ryl7889r02Nja0d+9eXWWuU8d02xwmAADwpG97etQ+pPX1de3Zs2fkvq2DSFMvvPBCYcbj5MmTeve7361/+7d/06FDh3ThhReOfQ2CCAAAy6dNEJnbsOpFF11U+PnVr361JOnSSy9tFEIAAMDOx8qqAADAm7levuu65JJLNKcuEAAAWFJURAAAgDcEEQAA4A1BBAAAeEMQAQAA3hBEAACANwQRAADgDUEEAAB4QxABAADeEEQAAIA3BBEAAOANQQQAAHhDEAEAAN4QRAAAgDcEEQAA4A1BBAAAeEMQAQAA3hBEAACANwQRAADgDUEEAAB4QxABAADeEEQAAIA3BBEAAOANQQQAAHhDEAEAAN4QRAAAgDcEEQAA4A1BBAAAeEMQAQAA3hBEAACANwQRAADgDUEEAAB4QxABAADeEEQAAIA3BBEAAOANQQQAAHhDEAEAAN4QRAAAgDcEEQAA4A1BBAAAeEMQAQAA3hBEAACANwQRAADgDUEEAAB4QxABAADeEEQAAIA3BBEAAOANQQQAAHhDEAEAAN4QRAAAgDcEEQAA4M3cg8jDDz+sQ4cOaffu3dq3b58+8IEPzPstAQDAkujM88W/8Y1v6JZbbtFf/dVf6c///M9lrdXTTz89z7cEAABLZG5BpN/v67bbbtM999yjm2++Od/+x3/8x/N6SwAAsGTm1pp56qmn9OKLLyoIAr3lLW/R61//el177bV65plnap+zubmpjY2NwhcAANi55hZEfv7zn0uS7r77bn3605/Wt7/9bZ177rl6xzveod/+9reVzzl27Jj27t2bf+3fv39ehwcAABZA6yBy9913yxgz8uuJJ55QHMeSpDvvvFMf/OAHdfDgQT3wwAMyxujrX/965WsfPXpU6+vr+deJEyem++0AAMBCaz0j8vGPf1wf/vCHR+5zySWX6PTp05KkN77xjfn21dVVveENb9ALL7xQ+bzV1VWtrq62PSQAALCkWgeRffv2ad++fWP3O3jwoFZXV/Xss8/qyiuvlCT1ej09//zzuvjii9sfKQAA2HHmdtXMnj17dPjwYd11113av3+/Lr74Yt1zzz2SpOuvv35ebwsAAJbIXNcRueeee9TpdHTDDTfolVde0aFDh/Td735X55577jzfFgAALAljrbW+D6LOxsaG9u7dq6vMdeqYru/DAQAADfRtT4/ah7S+vq49e/aM3HeuFZFpZRmpb3uejwQAADSVnbeb1DoWOohkV978QN+WFrZuAwAAqpw+fVp79+4duc9Ct2biONbJkye1trYmY8zY/Tc2NrR//36dOHFibCkI1fgMZ4PPcXp8hrPB5zg9PsP2rLU6ffq0LrjgAgXB6CXLFroiEgSBLrzwwtbP27NnD/9nmRKf4WzwOU6Pz3A2+Bynx2fYzrhKSGZuS7wDAACMQxABAADe7Kggsrq6qrvuuotl4qfAZzgbfI7T4zOcDT7H6fEZztdCD6sCAICdbUdVRAAAwHIhiAAAAG8IIgAAwBuCCAAA8GZHB5GHH35Yhw4d0u7du7Vv3z594AMf8H1IS2lzc1NvfvObZYzRT37yE9+Hs1Sef/553XzzzTpw4IB2796tSy+9VHfddZe2trZ8H9rC+8IXvqADBw5o165dOnjwoL7//e/7PqSlcezYMb31rW/V2tqaXvva1+r973+/nn32Wd+HtdSOHTsmY4xuv/1234ey4+zYIPKNb3xDN9xwg2666Sb9z//8j/77v/9bf/EXf+H7sJbSJz7xCV1wwQW+D2Mp/fSnP1Ucx/rSl76kZ555Rn/7t3+rf/iHf9Bf/uVf+j60hfbggw/q9ttv15133qkf//jHevvb365rr71WL7zwgu9DWwrf+973dOTIET3++OM6fvy4+v2+rrnmGr388su+D20p/ehHP9L999+vN73pTb4PZWeyO1Cv17N/8Ad/YP/pn/7J96EsvUceecRedtll9plnnrGS7I9//GPfh7T0/vqv/9oeOHDA92EstD/90z+1hw8fLmy77LLL7Kc+9SlPR7TcXnrpJSvJfu973/N9KEvn9OnT9o/+6I/s8ePH7Tve8Q572223+T6kHWdHVkSeeuopvfjiiwqCQG95y1v0+te/Xtdee62eeeYZ34e2VH71q1/plltu0b/8y7/onHPO8X04O8b6+rrOO+8834exsLa2tvTkk0/qmmuuKWy/5ppr9MMf/tDTUS239fV1SeL/dxM4cuSI3vve9+pd73qX70PZsXZkEPn5z38uSbr77rv16U9/Wt/+9rd17rnn6h3veId++9vfej665WCt1Y033qjDhw/riiuu8H04O8bPfvYz3XfffTp8+LDvQ1lYp06dUhRFOv/88wvbzz//fP3yl7/0dFTLy1qrO+64Q1deeaUuv/xy34ezVL72ta/pqaee0rFjx3wfyo62VEHk7rvvljFm5NcTTzyhOI4lSXfeeac++MEP6uDBg3rggQdkjNHXv/51z7+FX00/w/vuu08bGxs6evSo70NeSE0/R9fJkyf1nve8R9dff70++tGPejry5WGMKfxsrR3ahvE+/vGP63//93/1r//6r74PZamcOHFCt912m77yla9o165dvg9nR+v4PoA2Pv7xj+vDH/7wyH0uueQSnT59WpL0xje+Md++urqqN7zhDWf9sFvTz/Bzn/ucHn/88aF7K1xxxRX6yEc+oi9/+cvzPMyF1/RzzJw8eVJXX3213va2t+n++++f89Ett3379ikMw6Hqx0svvTRUJcFot956q771rW/pscce04UXXuj7cJbKk08+qZdeekkHDx7Mt0VRpMcee0yf//zntbm5qTAMPR7hzrFUQWTfvn3at2/f2P0OHjyo1dVVPfvss7ryyislSb1eT88//7wuvvjieR/mQmv6Gf793/+9Pve5z+U/nzx5Uu9+97v14IMP6tChQ/M8xKXQ9HOUpBdffFFXX311XpkLgqUqRG67lZUVHTx4UMePH9d1112Xbz9+/Lje9773eTyy5WGt1a233qqHHnpIjz76qA4cOOD7kJbOO9/5Tj399NOFbTfddJMuu+wyffKTnySEzNBSBZGm9uzZo8OHD+uuu+7S/v37dfHFF+uee+6RJF1//fWej245XHTRRYWfX/3qV0uSLr30Uv7LqoWTJ0/qqquu0kUXXaR7771Xv/71r/PHXve613k8ssV2xx136IYbbtAVV1yRV5FeeOEFZmsaOnLkiL761a/qm9/8ptbW1vLq0t69e7V7927PR7cc1tbWhmZqXvWqV+k1r3kNszYztiODiCTdc8896nQ6uuGGG/TKK6/o0KFD+u53v6tzzz3X96HhLPKd73xHzz33nJ577rmhAGe58XWtD33oQ/rNb36jz372s/rFL36hyy+/XI888shZX9Fs6otf/KIk6aqrripsf+CBB3TjjTdu/wEBIxjLv4YAAMATmtUAAMAbgggAAPCGIAIAALwhiAAAAG8IIgAAwBuCCAAA8IYgAgAAvCGIAAAAbwgiAADAG4IIAADwhiACAAC8IYgAAABv/n9xpZshkmi60AAAAABJRU5ErkJggg==", "text/plain": [ "