From 2e9c1a0da1fb2333ef0465fced7f48b6d28fc4d8 Mon Sep 17 00:00:00 2001 From: Mikael Mortensen Date: Sat, 17 Aug 2024 19:56:20 +0200 Subject: [PATCH] Use ipynb on klein gordon --- docs/environment.yml | 2 + docs/source/kleingordon.ipynb | 1442 +++++++++++++++++++++++++++++++++ docs/source/kleingordon.rst | 621 -------------- 3 files changed, 1444 insertions(+), 621 deletions(-) create mode 100644 docs/source/kleingordon.ipynb delete mode 100644 docs/source/kleingordon.rst diff --git a/docs/environment.yml b/docs/environment.yml index 0692e152..97125cc5 100644 --- a/docs/environment.yml +++ b/docs/environment.yml @@ -15,6 +15,8 @@ dependencies: - pyyaml - scipy - shenfun + - matplotlib + - pandas - mpi4py - mpi4py-fft - plotly diff --git a/docs/source/kleingordon.ipynb b/docs/source/kleingordon.ipynb new file mode 100644 index 00000000..88730833 --- /dev/null +++ b/docs/source/kleingordon.ipynb @@ -0,0 +1,1442 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "755b1597", + "metadata": { + "editable": true + }, + "source": [ + "\n", + "\n", + "# Demo - Cubic nonlinear Klein-Gordon equation\n", + "**Mikael Mortensen** (email: `mikaem@math.uio.no`), Department of Mathematics, University of Oslo.\n", + "\n", + "Date: **April 13, 2018**\n", + "\n", + "**Summary.** This is a demonstration of how the Python module [shenfun](https://github.com/spectralDNS/shenfun) can be used to solve the time-dependent,\n", + "nonlinear Klein-Gordon equation, in a triply periodic domain. The demo is implemented in\n", + "a single Python file [KleinGordon.py](https://github.com/spectralDNS/shenfun/blob/master/demo/KleinGordon.py), and it may be run\n", + "in parallel using MPI. The Klein-Gordon equation is solved using a mixed\n", + "formulation. The discretization, and some background on the spectral Galerkin\n", + "method is given first, before we turn to the actual details of the `shenfun`\n", + "implementation.\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "

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": "2024-08-17T17:13:17.882929Z", + "iopub.status.busy": "2024-08-17T17:13:17.882616Z", + "iopub.status.idle": "2024-08-17T17:13:18.254249Z", + "shell.execute_reply": "2024-08-17T17:13:18.253919Z" + }, + "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": "2024-08-17T17:13:18.255844Z", + "iopub.status.busy": "2024-08-17T17:13:18.255696Z", + "iopub.status.idle": "2024-08-17T17:13:18.273853Z", + "shell.execute_reply": "2024-08-17T17:13:18.273635Z" + }, + "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": "2024-08-17T17:13:18.275090Z", + "iopub.status.busy": "2024-08-17T17:13:18.275020Z", + "iopub.status.idle": "2024-08-17T17:13:18.276587Z", + "shell.execute_reply": "2024-08-17T17:13:18.276365Z" + }, + "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": "2024-08-17T17:13:18.277731Z", + "iopub.status.busy": "2024-08-17T17:13:18.277667Z", + "iopub.status.idle": "2024-08-17T17:13:18.427276Z", + "shell.execute_reply": "2024-08-17T17:13:18.426973Z" + }, + "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": "2024-08-17T17:13:18.428668Z", + "iopub.status.busy": "2024-08-17T17:13:18.428580Z", + "iopub.status.idle": "2024-08-17T17:13:18.440348Z", + "shell.execute_reply": "2024-08-17T17:13:18.440123Z" + }, + "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": "2024-08-17T17:13:18.441655Z", + "iopub.status.busy": "2024-08-17T17:13:18.441585Z", + "iopub.status.idle": "2024-08-17T17:13:18.626031Z", + "shell.execute_reply": "2024-08-17T17:13:18.625761Z" + }, + "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": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "%matplotlib inline\n", + "\n", + "import matplotlib.pyplot as plt\n", + "X = T.local_mesh(True)\n", + "if rank == 0:\n", + " plt.figure()\n", + " image = plt.contourf(X[1][..., 0], X[0][..., 0], u[..., N[2]//2], 100)\n", + " plt.draw()\n", + " plt.pause(1e-6)" + ] + }, + { + "cell_type": "markdown", + "id": "09b92d8a", + "metadata": { + "editable": true + }, + "source": [ + "The actual `update` function is" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "0547e688", + "metadata": { + "editable": true, + "execution": { + "iopub.execute_input": "2024-08-17T17:13:18.627397Z", + "iopub.status.busy": "2024-08-17T17:13:18.627247Z", + "iopub.status.idle": "2024-08-17T17:13:18.631148Z", + "shell.execute_reply": "2024-08-17T17:13:18.630940Z" + }, + "tags": [ + "thebe-init" + ] + }, + "outputs": [], + "source": [ + "# Get wavenumbers\n", + "K = np.array(T.local_wavenumbers(True, True, True))\n", + "\n", + "def update(self, fu, fu_hat, t, tstep, **params):\n", + " global gradu\n", + "\n", + " transformed = False\n", + " if rank == 0 and tstep % params['plot_tstep'] == 0 and params['plot_tstep'] > 0:\n", + " fu = fu_hat.backward(fu)\n", + " f, u = fu[:]\n", + " self.image = plt.contourf(X[1][..., 0], X[0][..., 0], u[..., N[2]//2], 100)\n", + " display(self.image, clear=True)\n", + " plt.pause(1e-6)\n", + " transformed = True\n", + "\n", + " if tstep % params['Compute_energy'] == 0:\n", + " if transformed is False:\n", + " fu = fu_hat.backward(fu)\n", + " f, u = fu\n", + " f_hat, u_hat = fu_hat\n", + " ekin = 0.5*energy_fourier(f_hat, T)\n", + " es = 0.5*energy_fourier(1j*(K*u_hat), T)\n", + " eg = gamma*np.sum(0.5*u**2 - 0.25*u**4)/np.prod(np.array(N))\n", + " eg = comm.allreduce(eg)\n", + " gradu = TV.backward(1j*(K[0]*u_hat[0]+K[1]*u_hat[1]+K[2]*u_hat[2]), gradu)\n", + " ep = comm.allreduce(np.sum(f*gradu)/np.prod(np.array(N)))\n", + " ea = comm.allreduce(np.sum(np.array(X)*(0.5*f**2 + 0.5*gradu**2 - (0.5*u**2 - 0.25*u**4)*f))/np.prod(np.array(N)))\n", + " if rank == 0:\n", + " params['energy'][0] += \"Time = %2.2f Total energy = %2.8e Linear momentum %2.8e Angular momentum %2.8e \\n\" %(t, ekin+es+eg, ep, ea)\n", + " print(params['energy'][0])\n", + " comm.barrier()" + ] + }, + { + "cell_type": "markdown", + "id": "2654c25b", + "metadata": { + "editable": true + }, + "source": [ + "With all functions in place, the actual integrator\n", + "can be created and called as" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "b5ab745a", + "metadata": { + "editable": true, + "execution": { + "iopub.execute_input": "2024-08-17T17:13:18.632303Z", + "iopub.status.busy": "2024-08-17T17:13:18.632236Z", + "iopub.status.idle": "2024-08-17T17:13:21.068731Z", + "shell.execute_reply": "2024-08-17T17:13:21.068485Z" + }, + "tags": [ + "thebe-init" + ] + }, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGdCAYAAAAvwBgXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABRTElEQVR4nO29f6wldX3//5yZc+65u7C7RTeglAVW2pQaYrSL3dhggWoQ4h+1GiKpIcEAycaFQEg+2i0moDHdpNAfESOVNiGmxkqssVqlDZsYwNaQwKItJZEEDWHLioh+e+/K7j33nJn5/jHnfc575syP93tm3vOemfN8JDf33Llz5j13zrnzfpzX6/V+v50wDEMQQgghhFjAtX0ChBBCCFldKCKEEEIIsQZFhBBCCCHWoIgQQgghxBoUEUIIIYRYgyJCCCGEEGtQRAghhBBiDYoIIYQQQqwxsH0CeQRBgJMnT2LXrl1wHMf26RBCCCFEgTAMcerUKZx//vlw3fyYR6tF5OTJk9i3b5/t0yCEEEJICU6cOIELLrggd59Wi8iuXbsAAJ9/8p3Ycban9JxxWO5PGjnTUs9j++Xar6s9Uo2y7xeygO9l+zR937F9n23DORS1f+bXU/y/K4/P+/E8Wn0XEumYHWd72LlLTUR2IsRWMCzR2gDr7qTE89i+TvuLNtReT2KWnYiWmir3nlld4v8rfC/bRuV9HL1m1V+rrWCIHSWeV2f7ALTPoa57r277KmUVvSxWLduh1nUz7nL7dZzDujtJPYc6RIeYQbxmfI3y4TVqN6bvPWXvj6vefhGtjohUQVx43Qsn9q/6wnW1ffGcOt646+6ktmOR5pBfL0ZKKNBdRNx7xOM66KoE2G5fhd6KiEB+Q+pQd2e8yu2T7sLXj3QV2x1wnf87Zc7Bdvs6tSi9TM0k6XKqpK72y5wDPw0TQlYZ2xJSNl3etfY7EREZh4N5MVJZupoqqTO8WCY6Und4U5eupXZMyluXrgMpD99DEbb/97sYhajzHJpsvxMiAvQrVbHK7eu2KX+3KUNtIO08utSxkGWafm/ptNeG/zeb9x1duigAbWm/MyIC1F+3UCY6YLt9+fldbV+nLRvtt0U8VKCcdIcuva8AO++trP/7pt7TqygBbWi/UyIC9CtVYTs6Ybv9vOOr/L6Oc+ha56CC/DdRSuzTp/dY8m9pqgO0fc/JousC0JZz6JyICPrSGbehfXE83fbl53cVnfPvSofS9dekb3RtOPSqvX9sS8Cqtw90WESA/nTGttsXx9Bt39QNq603wraeF+kOuu8h2wWbTWHzb9S99/UlCtGG9gWdH75b12ygQDuG2Xahfc4uSUgz8P+sGXTufXVQpd/qm4QAPRARwarJgK32eWMkhPSRog9YbRCAuuriyg5NNnX/742IAKsjA020nzwHRkEIIatA2n2uDRLSh/az6JWIAOYXbluV9sU5yN8JIWQVkO+/fZEA2+3n0TsREfQpOtHF9gkhpOv0JRViu/0iOiEiOovnyNjujE2mSppsnxBCiD62oxBNrRVTlU6ICGA/VWG7fXEONtsnhBCihu0ohG0JGofqs4N0RkQEbYhOrHL7hBDSFWzc92xHIWyngsq03zkRAeynKmxHR2y339b2CCEkSdP3vDL0JQpStv1OiojAdmdsOzphu/2841NACCE2ke9Bbb3n2U7F2G5f0GkRAex3xm1o37aQJTEx9I0QQnSQ70F9XzCvi6kgmc6LCGC/M7bdvjgH3f1NSwIlhBBikybuQWXasB2FsN1+kl6IiKAN0YkutE9BIISQ+lD9YNeGVIjt9tPolYgA3ZEBG+1zmnZCCDFHU2vVdG1UTBG9ExGg3TLQVPtpa8UQQggxC9eq0aeXIgLYr9uw3b44B/k7IYQQ83CtGj16KyKCNkQnutg+IYSQavRhrZqy7esszdJ7EQHaIQNcK4YQQogOXY2C6LbfCRHRmbM+izalSmy1TwghpBt0dYKyMu13QkQA+6mKvrRPCCGkvXR1grIqEtQZEQHakSrpSyErIYSQdrEqqZgknRIRwH6qRJxDH9onhBBin64WpNbVfudERNAXGWAhKyGErC6rGgWR6ayIAP1KlZRZK4YQQki3WeW1agSdFhFBn6Ijde5HCCGk/aziWjUyvRARYDVkhGvFEEJIf1mFtWrS6I2IAP2SEa4VQwghq0ef16rJolciAvSzboQSQgghq0Mf16rJo3ciIuhTdIQQQsjqYVNCmiwF6ISI6CyeI2NbBjjElhBCiC1sz9KqSidEBLCfKrHdPiGEEKKK7VSMzhpxnRERge3ohO32CSGEkCxsj4op037nRASwH52w3T4hhBCSxHYUpGz7nRQRge3ohO32CSGEkC5GQWQ6LSKAfRngWjGEEEJs0dUoiEznRQRoR6qEa8UQQghpmjL9X5skBOiJiAjaEB2pcz9CCCFEha6sVZNGr0QEaLeMcK0YQgghpmjzWjV59E5EgHbICNeKIYQQ0jRtXKumiF6KCNCuuhFKCCGEkKZo01o1KjQmIkePHoXjOLjzzjubahJAO6IjhBBCSNN0QUKAhkTk6aefxkMPPYR3vOMdTTS3hG0ZIYQQQrpGlflJdNaIMy4iv/71r/Gxj30Mf//3f49zzjmn1DF05qzPog2pGkIIIaQLNDk/iXEROXz4MD74wQ/i/e9/f+G+4/EYm5ubsS+B7VQJZYQQQkjfsTFLa/VQQw5f+9rX8Oyzz+Lpp59W2v/o0aP4zGc+k/l7cXGq5r3E83Uvdl3tE0IIIW3D1iytxiIiJ06cwB133IGvfOUrWF9fV3rOkSNHsLGxMf86ceJE6n6MjhBCCCH1YHutGmMRkePHj+O1117DgQMH5tt838eTTz6JL3zhCxiPx/A8L/ac0WiE0WikdPytYFjLBagSHWFkhBBCSJexvVYNYFBE3ve+9+G5556Lbfv4xz+OSy+9FJ/61KeWJKQMdaZK1t2J1gtCCSGEENJ1ynwYr7v/MyYiu3btwmWXXRbbdtZZZ+HNb37z0vaq1BkdUXkxKCGEEEL6hEr/Z6rv683MqlwrhhBCCClPE2vVpGF01EySxx9/3OjxTUZGKCCEEEL6jo3+rzcREYGJhesoIYQQQlYFuc9rov9rNCLSFHUXsRJCCCGrRJN9XyciIjpz1stwvg9CCCGk3XRCRIBqa8UQQgghpDl01ojrjIgIuHAdIYQQ0k7K9LedExGgG9ERyg8hhBCbNN0PlW2rkyIiaGt0hAJCCCGkLZju96oev9MiArRz4TqOtCGEENImTPVLdfSlvRi+24a58ps+PiGEEJJHE/2Q7jptaXQ+IiKjetEpCYQQQkg9VO1TeyUiQPEFoYQQQggh9VJlLbbeiQiQLRuUEEIIIcQcZfrZXooIsHwxKCGEEEKIeXT7296KCLC4GJQQQgghpDl0lmbptYgAlBBCCCGkzfReRAghhBDSXjohIjqL5xBCCCGkOk3NEt4JEQHszJnPqdoJIYSsMk30hZ0REUHTc+ZTRgghhKw6JvvezomIoO6LwggIIYQQEpHVH5roJzsrInXD0TWEEEJIPib6SoqIBGWEEEIISe8PTfWRFJEElBFCCCEkjsm+kSKSgnzBKSaEEEJWGdP9IEUkAwoIIYSQVaeJvpAzheVAGSGkOuvOtrW2t8I1a20T0nWa6gMpIoSQWrApHFmknRPlhJB2QREhhJSijeKhAuWEkHZBESGEFNJV6VCFckKIPSgihJAl+i4eKsjXgFJCiDkoIoSQORSQdCglhJiDIkLIikLpKAfTOITUC0WEkBWgTdJxljs2ctw3gpGR46pAOSGkPBQRQnpI0+JhSi5MnUMT0sJ0DiFqUEQI6RGmBaQNwlEHeX+HCUmhlBCSDUUkh61gyNlVSesxIR99EY4ypP3tdcqJeL0oJIREUEQI6Sh1C0hT8rHuTI0cdys0dzuTr01dUsIoCekCTXwgp4hksBUMbZ8CIUt0QT5MiUYd7dYhKyalhEJCVhGKCCEdoE4BqVM+bElHWdLOt4qc1C0lFBLSJpr6QE4RIaTF1CUgdchH16RDlbrkpE4pYdqGtAnT6RmKSApMyxDb1CEgVeXDpHisO/Xf1LbC+v5vk3+7rpjULSWUEdJnKCKEtAjbAlKHfJiQjDrarSIq4rpUiZRUERKmbIhtTEZFKCIJGA0htqgqIWUFpKp82BIPXZLnWUZM5GtVNkpCISEkDkWEEMvYEJAy8mFaOFTPqa5hull/j6qglE3f1JG2oZAQG5iKilBECLFEFQHRlQ9d8agqHWbrS/SPrSMvaX+7ipyUiZZUjZKwfoT0AYoIIQ3TlIA0IR9dGUmTdZ6qwiBfGxNSUkVIGB0hXYciIsH6EGKSsgJiMvqhIx91SsdO16/tWABwOvBKPa+MoOjWmugUulZJ21BISFehiBDSAGUkxJSAqMhHFemoWzLqbFNVWHTqP1SjJbojb8pGSSgkpGtQRAgxSBNRkCJpMBX1sCEcVUk7ZxU5UU21qEgJhYSQOBSRArgCLylDlwRERT6qSMe6E5Z+bhW2QkdpP105qUtKyghJ2RoSyghpM67Jgx89ehTvfve7sWvXLpx77rn40Ic+hBdeeMFkk4RYx7SErDvTXHlYdyaFEiKOkXWcna4f+1I7rzD1yxZVzkf17y96LRb7Zb8mqscAovdIueHa9S6WSFYXE7WURkXkiSeewOHDh/HUU0/h2LFjmE6nuOaaa/DGG2+YbJYQa5StBVHpXKoKSJF8AKgkHnWw7riFX9Xb0JeTIilRubbRfsWvjwplZYRCQtqI0dTMv//7v8d+fvjhh3Huuefi+PHj+MM//EOTTWvDETOkCiajIFU7tyJUxKOMaNQhDXUddysMCo65/PdlpXbE9cpK3xSlbsTrVbWGpErtCFM1pE00WiOysbEBAHjTm96U+vvxeIzxeHFj3tzcbOS8CKmCqRExNgVEVzxMSUddJM+vSEyi5yyuQZqUyNevSEpMCwllhHSZxkQkDEPcdddduOKKK3DZZZel7nP06FF85jOfaeqUCKmMDQmpIiBF0Q9VAalLPHY61W9Bp8My09Uvn3+enFSVEtNCUiY6wlE1pC04YRg2Uk12+PBhfPe738V//Md/4IILLkjdJy0ism/fPnzh+EHsODv9n7CuES15qRmOmiFp6EqIrShIHfJRh3jUIR2qlJETGZWISdGonLyRN3lSkTcXicoImzIjaygjJEnVPvH0KR+3/t5xbGxsYPfu3bn7NnJnuP322/Htb38bTz75ZKaEAMBoNMJoVH5lSkKaomkJsSEgVeSjSelQaV9XTOS/PUtKVKMkdUZI1p2pkdoRpmqITYzeLcIwxO23345vfvObePzxx7F//36TzRHSCE1KSN0CYir6YVs8iqgiJjpSUkZIssQiS0hM1Y5QRogtjN49Dh8+jK9+9av41re+hV27duHVV18FAOzZswc7duww2TQhRtCREFNRkKLht+nPqTf60XbxKKKsmIjrVCZKkiUkRWKx7kwqRUcoI6Ru6p7o0+jd5MEHHwQAXHXVVbHtDz/8MG666SaTTRNSO01JSFsFpOvykYeumFSJkux0/VLpGsoI6SvGUzOE9IG2SkgZAVGVjz6LRxHy364qJTpCUqZ+pEqqRrduhDJCmmR17zSEKFKnhJiOglQVkFWWjyxUpaQoSrLuhMrpmmj/bCFpIjpCGSFNwbsOITnYkpC6BITyUS9VpSQvXQOoF7RSRkif4B2IkAxUJaTuVEwdaZgiAaF8VEdHSpIRkjL1I3WnaigjpC3wbkRICnVJSB2pmLoERFc+Rk796y+NQ/OTA9o47yIpyaoj0UnX1J2qoYyQtkARUaDuoUqk3diQEJ0oiE4apg3yUXT8KnJi+nzT2qkiJWlCopuuqTNVQxkhbYAiQohEWySkjiiIqoQ01Zmrtp/X0ds+V/kcVARKvAZpQlIlXaObqqlrRA1lhJiAIqIAoyGrQZslxGQUZBxOrHXwutGQtP3bICdF7HQGStGRaHt6ukZFRqLty9ERlSJWVSgjpO4+kSJCCLonIVWjIEnkDl63Y6+z7qNozo60v69s+03/nbrRkaZlRHfiM0LqgiIyY92d5K42SPqL7tox2cfppoQkqUssqqyAu9wxu5WOWafAVCUrOqKSqmmLjDAqQuqEIkKIInnREJsS0pahuGmSkDXbqAqiA46iA+WOkyUwtq9ZlVRNW9I0lJHVxESpQjvuYIRYoo6UTBMSYjoKUoa86ESWOKStTpvbRuABOasJpxFfeE4enbK4hm0QkyqpGpMywpE0BGg2S0ARkWB6ZrWoKyWTfmyzEmKj4yxKiyQ7zzzpSJu0q8q+8nVMtrtIcQSzn5evp00x0UnVlJWRNOqSEUKqQhEhK0ldU7enyUVZCWljKkalJiNPQIokomqKQFzX3HZm1zkpJNG29EiT/Hc3cb1VUzVlZaTMPCOsFyFJTI0gpYgkSEZFOHR3temqhOhMwJWFbuolS0BUZCOtkywi6lzzj73uTBfnkhCSqN3sKIngdDgtJSO6r4FqqqZJGVGFMkKqQBEhK4epuhCTEqLTESaHpY6coZaMZAmISt1HnoAUyYZOZxh1ntnHW0zuNZjvryIk0fb01E2VyIjOa6CSqmlKRlgvQgQmP5RTRMhKsWoSoovOyJc0AZE7syxRyBMOlU7vLHeceYzFeixSVFOKnKQJSbRdvcBV5fVIex10Z2RtWkayYL0IMQ1FJAUWrfaTOopT2yohRQJS9InclIBkCUNWx1b0iXrd2U59rhDH5U53ETnJEpKdrh8bKhw/n2DpNWkqOpKWqlGRESAemVKVEdaLEFt0UkSaqNtgbcjqkhUN6aqEFKEqIWUFRFc6cj8EpJRyZMlJ9v56QqIrIyqvR5VUTZGMpEEZIWVoqh/shIiMnCnW3fS1NghRoa4p3MtgOh1TZa0YFQkRnVxW/UeagMgdVlqnlCUbp4P8Dux0sIadbvy13ELO6tizPyUvnQMUC0laUWtWqkbl9ag6q6vJNE2d69IQogLfbYQoUDYakrWKbvw52SM2VEl2bEUr2pYREJ3ox1a4FpONNMFIq1MYK6RElzvOSerxd7rbc0nZ8tcWMlp0ueXL4PqloiNpMlJFPtJqRpKoyEhVWC9CTNB7EdkKcj4tESJRNRqSNcNqfJ9mJiurWg+yFTpLAqIiH8Ai2nE6WIs9JykZaVGRMwmhGPsDjLz4+e5IRETGKf/fI3eCLX9WGxJMYlISF42lp8YJkDoPiaqMmCRt0rMkSRlpKkXD9AzRodciIm50lJHVRiUtU8comSQqKZml5xiePEtHQpICoiof4jnjYBiTjaRkAJFoxH4OBkuPR0H8nMdu/DlpoiLaXXcngAts+cPFkF45lTP709/AKPs9IC6PFB3RSdXURZmRNGkw9UKaYKzxHuvtu5GjXkiT2ErJqKKbitkKB0sCklb3kRb9EAJyJljLFY20n0UUAwC2Z5/k15JC5y0kcOROY8cYudNYmyNvih3uNtbdCcYYYisczlM5sSgJIAkH0plFR6oUslbFVIrG5Ho0ZPWI7gvqK2X3UkTSJIRREZJF1WhIm1IyaZRJxWyFQ2yFg0L5iLYvCwgQRTvGwUBJNuY/+1JBrCQU61LUQ35OnqQAUTRl7A7mQgIAYwwxSqZogIWA5KVuMlI10c/x6EhTU/GXWZemiagI0zNElV6KSBaUkdWj7oXtmk7J6M6KmqRsKkZIiNyRFAkIgHkURAiIkA5ZHmTZAOLCMZH2m0j7yduHvjQRmSQoa56/JClb/hDr3iRXSOS0DZCSukmTEaAwOlJHqibt9VddJK+IpIwwKkLqoEw2onciwpQM0UE3GlK0T9k1ZNKoMj9I2VSMnIbJGvmiKiDbgTeXjq1EiiZLOKZ+dG2SsrLmLa7ZxFMTlG3fw5rnz9rPFhIAS3Uk8vbMlI3hVE3e61/XzKtJdGddLYJRkdWibP/bKxFRuQiMiqwORdEQEwWqxeeklpKRO6E61orRTcUICUmTDwBKArLlDwqFA4hLx3T22Pfj12kqycfAWwhHkaAMfR/r3nQuJNHfMMgVEgDp6Zms6IjhVE2VqFgdKRpGRYgKVYIAvRIRVSgjpA7KRENMo5uKkQtS5VSMkJAqAjLxvULhAOLSIR4HCRHxvUB6HP3O84L5cQaePz/+mufP2514HiaBNxeS7cBbqisRQ4SVhKShVI3KwoWqKZqikTSmUzSMipAieiMiTMkQGdPRkLpSMkXREHlbHXODAGqpGCEhZSMgQkK2fU9ZOgJ/0VmGfrzj9OHC8aJrGsxEJE1OgOXoyZrnZgrJujfByF1cuzPB2rKQpMmIeCxTU6qm6pT9RfUiJiY6I6tN1f63FyJS5iIwKtJfqhSoqtSFlDuuuVEyVQtSo5+XUzGqNSCi/kOOgggB8X03Jh+50iHtJ7Y5voNwJiCYfRcCkiYnrhfM2/Nmj6eepyYkwTQmJXOS0RFLqRrVqEgauikaRkWIKnUEAXohIoTooDuDah3REFOYlpCkgADITcMkBSQpHkXSAQCYOhBdphMAcGc/DaLrGnpStCQhJ4H0u8Bz4XrpkYGJny8kon4khhQNWXcnUcc6kxF5lV8Ei/VqilI1pob4lhlFQ4gtKCKEGKbJOUPyilKjn4exT7U6EqKahpElxN+e/e0z8ViSDgDOdLZtdurOfJ/o53B2+vM/zZVEJCkn3qLj970Agecg8Fz4XgDfk6MjxRPOyYjo6VY4jMkIMJvozV0I6lY4yJSRJHXISNmoSBmqzD/CqAjJgiJCekXb0jImZ09Ndj5ZNSHR40Es1J6cIyTaJ56OUS1GlQUEiGo/ptseQt9BuB1ZxFwuMqQjeozYPnKgKXIeB+FcPKLtS3IyWHT6oecAXhhFSzT6Pzk6shtbsd+N3ElMRiIB2Z5dv+h2KiJopwNvKVqWVjNiCt1aEZX0TB4cQUPKQhEhRIM2pWXySBamLrYnZkuVRsfIEvJ/k53aEiKiIMFMQELfgTOOOt2saEdSOpyU76EHOEI8Zs8X/af404SYYCLVlLgOMAgReg58zNI2KUIy9ZdrR2TSZES+fiJlI1J+ovNeREj0UjRpK/cWkRUVaVuKhlERkgZFhKwUVetD9J6b/sm3jrRMXjRElpDFtuyUjDw6BohLyOZkXT8V47vAtgtn6sCd6IuH+L6UmpG+u+LnbWe+PU1OwlmNSQgAXogpkJqqKSIpI0I+dkojbKqkaJqizCRny8fgonmkXlb63cSRM/3CdlrGZEeTN3Q3T0KSHYY8TPd0YjXcrWCYKyGqqRjHd+BsO3B8B+5YXzzEwBVH+iAfukCQiH4kv6dFTcJB9PsAbhQdAbRTNYIsGRGPi1I0soyYKFxVrRXJg5JBbMB3HCGKlFlhN4mJaMjS7xMSkoyGyMgpGV0JyUrFOAHgbjtwpg68mRvqiIfjL7YJggHgTmZRDnexDcgXk9B3FpGSAAjhaqVqdmERQcsc2hvEH6elaFRkxBRV0jNZdSJZwqJSJ8L0DElCESG9wEQ0xERapm6S0ZA0CREkoyFLI2TkmhBFCUlLxTg+4I2j764kInniIbbPH/vh/Heht6gtmadhBovn54nJor7EwSxTopWqOYXRXEY2JjuwZ3hm+UVwgZ3IT9FkyUjdqERF6kjPEFInFBGyMujWh+hQZYXdsqSlZKLHcQmRp2/PkpDNyfpcQn49WZtLyJntYerkZHmpGHcmIEJE8sRD/CyLh/id2BZ4DjyECAZStEMSEyEgaWIiZARYpGx0UzV5MiJSu/M0V0aKJktG5KiInJ4pU7Cah05UpOroGUJ0oYgQkkJRNKRMWkaFos4nfQKz9OLU6Hv8X7yqhKikYhxJRMTPQFw8gLh8iN+5fggRoHKn0e8jCXHmdSaB58CRhvF629E2OVoCRGKyXFuin6oBsmXk/8NZOGfwRmGKRlVG6qBMVKQMnFOE1AVFhHSeJopU8/ZLRkOamMBMfLpNK05NSkhyHZkqEqKSiol+jkTAG6eLBxCPejjThXhEMiL2jYbfur4DjIWERMcNB4uZVJ1BmBktmUdKYkWtaqkaWUy2XGmIbTCQts/kUZ551VtERcSsq2kyIr+e645rdLbVPIrSMzpREc4nQnShiBDScsSIGfEpN09CBGkScjpYw4a/c2nWVDFPiFwTIkvIZHsQRUE0UjHuJPrubYdLEZFk1GMhIyGc2WN3EsLxo78znK0jEwwdOJ4zkxAHwdQpjJZEa9VEvxMpG51UzXh7OH98WrKSTaxj9zAaRfN/k534jeHp2R8U1YuI6d/PcseZMiJeQxP1Iqor8xLSBigihCSoOy2j8wlXtzYgrzg1vt/a/HdpErI9W3tFFKbGIyHVJEQ18uFOxONgvo/jBwi9SDY8P5KSYOjA8UM4noNw4MD1QwRTJ9pnGsL1Z5GSKYBRiABO9rDhqQMgjEbYTB2EngMHLvxtzAVEXj14Iq3sO4+EzK4pEE12JupF1r1tvBGMYjKSfO1Mpmiq0MQwXqZniIAiQlYCk4WqNsjL78vRkKXfSROXxdIL/uLxttTxikXr5ovUTRMTlCW/B4vRL7KE6ERA3OnisTOJDuzCRTBwIkmZLXkrf7aXu+8ADlzMCl396HE4daL16abR8xbpomj/KF0TAn6UpnG82XO9aJ4UMZJm4nsYzkR0O/Cw7sWldRwM59u2wrWltKGIipQhby4ZW3DeEVIH7dBvQnqCzU+0WdGQtGW6ZQmRxWMqTVImRsgkF6pz/Kgw1fUXwiGkI1kPIiTEnYZzCXGmobKEONMg+j6Lksi/j563kJp5m9OozahdzM559j0xcif5XYiWEC8hYrGoiFwU7C/W5ZlvS7wOqUKY03lXnZRMUBSJa8NMr4QAFBFCjFKm8FB1FlUV5ND3OEVItqVOdTLrbKexiIgD+G4kI1MHTpAfDRE/z8WkIBUjJESWDGfizyUEflxGxL7zxzMZcWbHdv3Zl9zmbHiwEKXl0Tvi57hwCQETQiYLmxxBElJ3RpqpNjlrbfRaLL8XOH8HIUzNkI5TZcSMehv1r8pbJ2kdXPJTuPwpXXx6F5/mgcWnfNHZilEyadEQID5luxwNketBhBDopmKc6cxofB/Oto9wDUgmohbpmfxUTYBZesWbicdsFE1arUjohdGQX0QSkpaeERGjoetj249W6Y1d52AYXxQvJT1ThjakZTifCDEFIyJkZWm7YBSR9mk6ubBdEvGpPa0+RE7LALNoCKTUjIh4SHUihdEQxVRMMgoiJCT2fRrEoiPy8xfHjadq5ukghahImfSMiCjJ13N+XVMiUEB6Ci2t5qcO+dCJyJmaG4eQIigihFTAVJ492QkV1Q3kfVIVaYK0jlFEQJJpmbwiVbk2BMiPhiinYmaSIcRDyAem05iMlEnVZNWKpEVFdNIz82ufqBNJS4Hp1okQsko0IiJf/OIXsX//fqyvr+PAgQP4/ve/30SzhFjFxsRUQLzTyytUFZ/m09IyQHaRavQY8+950ZBo+0JAkuKQFQXBeBxJyDQuI7HoiIiiJMRmITyL4tVkVCTaHi+2jf09U2exim9CRkTkSK4TEYiIkyyGcmSqjHzUmZZpy/BgQmSMvysfeeQR3Hnnnbj77rvxwx/+EO9973tx3XXX4eWXXzbdNCHaVFnorinS1pVZfhxPy4hP6cn6EJ20THLl3NjU7BnREFlCov0yClKTUZCpH31NJjEZ0U7VJKIi8rkn/6bUotWM9IyIIAl5S0vNAGrpmfjU/P2YcKxvw+WJWYyLyF//9V/j5ptvxi233ILf/d3fxd/+7d9i3759ePDBB003TUjtNJlHH4eTWj8NJ+tDdNIycpFqMnqwNIRXlhLFVIyQjLmETCbAZIJwMo3LiGaqZp6ikYf3pkRF5n/PNF4PE/pOZnomq05EyEdy5EyZac9tFqk2UUPVRLE5aT9GRWR7exvHjx/HNddcE9t+zTXX4Ac/+IHJpgnpNWkL3AHpaZm8+hGVtIxcpBqLIChEQ5RTMUIyZhISTqaRhEynCRnRS9WI2VqzoiLibykqWi1Kz6TNJyJTNT3TZrpe9E3sY/Q/4vXXX4fv+zjvvPNi28877zy8+uqrS/uPx2OMx4uQ3ubmpsnTI6RWTOffVcP2cqeXLFRN1ocA+WkZJ1GkKh6rRkPkKAiAeBQkanQhIMAiCgJEErK9jWgS9uUhvAJnG6lDfEPPnQnQLJrjhdFCd7Np4aMp3aPzDmfDe0NvJlweELph9FwvROg7CHwXnhdg6ntY8/z5LKtpw3hjM6wGQ6y76ZGNJtadyaOOVXgJqUojlUuOE3+jh2G4tA0Ajh49ij179sy/9u3b18TpEVILpvP7qqIjh7t3utFj0RGO3OVPrwMvOm9v9t31ZkWmXghIC8gJdwm92ZcrPZ6tdCuviBsMo/VfQs9FOPSir4ELeB7CtegLngcMBsDAi76GQzjD2eejwQDO2lr0fTgAhsPoazB7jufNjwUvOnY4mLU1Wygv9FyEszVpgoGzODfPgSjrCAaLvwEAwkGI0BN/9+JauLPrI6Z7H86+r3nL6Tp5LpEsCYleq6n0uPmZTikhpA0YFZG9e/fC87yl6Mdrr722FCUBgCNHjmBjY2P+deLECZOnR0hnkWtV5AJbuUhQdIBpBbjiE7xYN0V0pkJGHG/RKc5XrR2Ei856Jh5A9D3ewUcdf+hFEjKXESEkMxmRJSKSkUFMRpzhIC4hQlZmEiLLzFxAZhISCPEYLs4riobE5UmWqWi/xWMhYfCCuZx5cxmJvq9709n3meh56WkKWQ77lsroW6qJNI9REVlbW8OBAwdw7Nix2PZjx47hD/7gD5b2H41G2L17d+yLkDbR5JTcI2eotRJvETtm0RERFVlLfLIXn/TFJ3940cq3kCIEckctOvL544yoiCwGsozI0REhJrKMzKMjchREkpC5zEhRkCXxUYyGzP8GEQGa/V2OF86lTEjI/LoJiZt9F9dVCKCIRgnKjCSp8/XXpQnB4Oq7BGhgive77roLN954Iy6//HK85z3vwUMPPYSXX34Zhw4dMt00IdpUWR21KeS8vrz6afxxfIn1kTuJph/3phgHA6x7k9nqsVNMAg8DL8C278HzAvi+C9cL4ftO1An7zjxy4EwX0uH4UUfuBGJFWweeHyLwHDgDwBXzcCCakt3xg6g+I/H5x0FU4wFE9R5LyKkYIBYFARBLxQgJAWYSNBOQvGgIkEw9LdIyzjxdlZ+WSUt5AdlpGfk9Fo9u9WOejzIjhMjqYlxEPvrRj+KXv/wlPvvZz+JnP/sZLrvsMjz66KO46KKLTDdNiFVOh1Mrk5qd5Y7nHcG6O1kaRjpyIxlZmxVaDt2o8HLN8zGdyUjgu3C8EKE/65Bna7GEA0RiIkdG3EWhZ+jNxMMHgqkDbxoi9Jy5jADx9WEEolIhXEuRETkKAsSiIADmUZD546F4PBMhhdqQpSjPYBEJyUrLpCGiTrJoVE3LjJxhbcN4+zJPCekXjdwlP/GJT+ATn/hEE00R0iimRjokQ/I7nUHuNO95C5LtdLex5UcjN84kpGTN87HlDzCcLeg28PxZRCRA4HuAF0Qrxg1CYNuZRw5cqEdFHFlGJmEkC7Ml6rKiIzEZkVMxwFJB6lIUZCYdWdEQYDkaEu0fL1J1gMy0DJBdHzJKiYKkpWX6VitCSFn6EQckpARdL7JLm1wtPgpjOc+RrBMBFh1q7ugZxItWgeWOXd4m14qEUvFqMHTmdSMAlEbV5I2KKZIQEQ2JBMqZC4mIihQVqUbXIZ6WAbLrQ2KvhUJaZrFtWWbrqA8pWqMotm+D9U+EyFBESKdpotit7cKS9sk6+Qk8tfPzFtvSRs+4XrCQkMT35FBeuYNPjqAJPGchJN5CSFKLSzNG1WSNihGpGCEhQnoWo3eWpSMpTVlFqsByWgbIH7ablI+6Zg61WbQqyJsYj5AqUEQIMYjOJ1JBXqdTpZhRpAzkIabyRFzJ0TNzRCc8CNNTGkVRkYE8hHZRvyFHR3JH1eSMipkLiHRsISHzNlOiIckhu9E5Lc8dAsTTMslhu9F1jR7vcJfnb5FJE8YmlwwgpK1QRAipEZvFgOkh/+3UFIGcSpA/2cvpGW82f4YYwip31OEgnKczksWeQDwqImRkniaRUjUAClM1WqkYKRJSFA1J+66TlgHS5w9Jvg669SF1FTgXSTAnMyNtgSJCVoK+DSfMK5AVnVxqByhJSVqdCBAXE3lOEQCxmVaT34NB1NHLEYdgsJCRuZDopmo0UjHzdIyXXRsSmw3WWxSp6qRlACxN6x5d0/TRMott5Ue/tCE9k6TtaUvSDfguIiRB0VwipwNPK6SuM4xXt7PZ6fo4HXjzOUSyRs/sdLeBABhjGEshxOpEfDkysngceG40xsULEQKAG63+EowAjB1gDXC3gWCpXMeBtx0iwGx9lykALxISZzrb5kWr4kZSEi2Uh6GHcLZqLwZebgQEs2OL/lAeniskJGtK+nAuIYshu1gLAC+AtxbVyHhegIHnz9MyQ9fHujfFmudj3Ztg5E4xcqfY4cYjT+KxkMG8CIiQyjbNIcLJzEiTUEQIaTliHgkxhHfdcbEVBvOJzZIyAiwmNzvLHYsZxQAsFxyOgqyoiNQprk3h+y6mQDSHyOw7toFgBIQ+ACxW5o1HS+Tts0JX3wFGs8XmBoAPJ1osT0QifJEGcuH4YWb9R2wG14R8zB/PIiD+2uznNennQYhgbVb3MooKc501P5KjmYSM1iZzCdk53J5LyO7h1lxCfmN4GuvuBCN3gnVngp3uNtad7SUJEXIrft7p+sbWl0lLy3AOEdJWKCKk82yFa6VHJ8irn5bdLzmXSCQJ8U+3dU9uliUj0e8ms6jOsoysB5N4QnYIbE7WlxtYA5C8pGtA4LvwAcB3EcJFOAWcmXjIK/WmfQ8Tvw+9cL7Sr+sB/ihaFTeYClFZXNO8qIfYNn8s/U6kigJJRAIPCEazdMxa9N1Z8+F4IQZr/rw+JhYJ8RaRkDXXn0vIyJsuS0gsPZMtIcnXE6ivPkSXoqG7OiNm+pYGJeahiBCSQt3pmbpIm9hMlhEEiKVo5CgJsEjRlJaR2fbAj9I0oecAfjRZmTOLgLgJARGlKEIU0sQkK1oCxKMeyWJYWT6SBamxWhURCRlFqZhgGM5TMU4iFZMmIbuG47mE7BmemUvIDnd7SULkaEj0+qRLiKmUjMpIrToKVaukb5iWITIUEbIyvBGMSi08poLKDKsqUZFxOFGuExFRkejx8vozsozARVxANGRk4AeY+tF08FPfg++7sVSNHB0BslM1eWLiI6ohkaMl/uwYyXQLsBz1SG6b/87TS8XoSkiyLiQ5nbuOhMjvjboLU3XSMknB4PwhxDQUEdILTKRnqiyAl5aeMUFSRpBTvJqsF0nKyNhd3A62VWbZnKVq5OhIALcwVaMiJiJaIlAVj+R3MTJGNRUDQEtC5GiI4Cx3nJjhNl1CTFAmGsIZVYltKCKEKJIUljLpmTpqRYrWnUkbSQOkh9LXnUXNyG8MT+P/JjsBALuHW9hESoQkBdfDUiFrWqoGSK8fSROTYBitXyMQtR5AvnjI3xdDc/VTMfLomDVXGiGTISFZKZk8CWlqlIyJItWstIxKfQjTMiTJSotI1loQpJs0UbSafwwzC+AB+SuwpqVo0kbSAMtREXll3nV3oiQjA285VZNVyAqIVE2EiJIk60jyxERHPABpunYP8/lBdFIxaRKye7BVKCGL12Mai6SpSEgdRaplZvFNwrQMscFKiwhZPXTrREykZ0xERfJkJNoWj4qsO9tLKZpxEHU6SRnZcodY8yLx2HIHmOSF8hOFrHKqBliMsBFyoiMmqdOyJyZYix5L68YMikfFACgtITJpKZkm0jE61JGW4SRmpG74jiJEgzrSM02QNpIm2h4VrooQ+jwqKM81EgxjU5br4HlB+pwjM4Q7zetICsREkCUeMekQzPcR09PrpWIA5EqITFFKJvo5LiFF0ZAyhapZ0ZC2zR3CtAxJgyJCekXb0jMmi1bzoiLRz/GRNNH5xAtX5Y5BrhdRZej5S6NqACylagRiZlQhJyFCYOoAoiZkLYQjakwkGckUD0k6BGKKdme2enCZVIw8T0hSQnRTMjJNzp5aJCHJaEjVtAznDyFloYgQYpgmJjgTpMkIFFM0p4M1JRlZc/1FqsYv+BtmqZo5IjqSJyeJqElSOqLnxFfIdaT1YOS1YtIEBMhPxaRJiKDulEyTE5hx7hDSVigiZOWoWidSlJ4xWbSaJH2Cs3QZkSc7AzBP0RTJyA53G2eCNYz9AUbuFOMg/baRFh0Rc44IgtnjPDmZrycDwPGdJemAF0B0qckF6lxJSMRidWWiIOLvzoqERD9vxyREpGSyJMREgSrA6dxJ9+mFiKy7E2wFemFEjpjpL7bTM+nHNRcV0ZWRGFkTnCV/RrqQxApZpehIcmSNQEhJnpw4AEKprkRHOpKPdQQEwFIURBSmAlhKxwgJkSctS1s/Ji0dU+Z1Txs1pTpSpqhIVSUtkxcNKUrLMBrSX8r0v0l6ISJAPReD9IciGcmLiqTJiKmoSJqMpM2umjV0VyCOkawZic4lPuEZgHQhUZARQEFIEiNrRDpk2/fmk4apyknsFKUIh0A8llcLXpMeFwkIgNQ0TDIKAiwXpgoJyVvEroyE6MyumwajIaRpqva/vRERHRgNIXWgKyNNzLZaFB2Zn09WdCSQHiPlZwklIfEWwjHwFw1O/cV12p5JSZqcCJJRDkFSOgTD2fakgADIrQORBQSAtVRMUkZUoyFpEsJoCGk7vRIRFSujhKwOpqMi5c5JLUUjd0RF0ZAkxlI1iDrxZEdVJCQAsOUPMJx12JPAm4sCsBCUNDkRFAkHgPnxo/Ocxp5XVIiaJiAArKZiBLqvv4yuhKQfo1fdBDFElahI795hTNEQHXRlpGgflRSNamSkSoheNVWz7k2jzqogVRP7vwoWQjJyJ/OJ0ICoxiIpJNuBFxMJWTJigpIRPQGKhQOIywqA0hEQ8fdG282nYrLIe/1VoiFlJi+re8guoyGrRdn+t3cikgejIatHlcLV9OMVR0XK1ItkFa5W+TQM1BAdCRI/y0hpm5E7wXoQvwmJzl8eZbPlD+eCACBTUDBErPg1SzjkY0X7xa9XnoBkyQeApShItK35UTFVUjJFmE7JkNWkjIz0UkTSLgQlhGRRd+Fq+nGK60VMzS2iIyM74cc/OaemZ7YjwZv9n+3EdrRmjSwk7hBnZuvYiKLQ5GMA8aHAw0hUBFnCkSYbWT/rCIgsH4vr0nwqRpc66kJMwGgIUaWXIgIwRUOapczU700UrwqUR9UA6dERIFVK5G1CSKJUyHAuJIIzwdrS1PHjxIRoskQISckTDQBLx9whTTiWFJAs+QAQS8HMn99gKkYFleG6daVkGA0hZVl3JzgD9Qn0eisiwEJGGA1ZbVTSMyYKV+tM0dSFbqomuXrvnFn/exbG0aRoiSiJqCMZzURknPJ/uBUMY9IAYB5FARaSkicaQHq0cxRLsywPwRUkox9if3lb0xOU6VBUF5L+nF7f9klLGGnMx9T7dyQlhDRFm1I0efOQaMkIkB4dkUkM9V33tpfSNsl0CoC5pMTOI5D2GyD1g0Ry4bnofJe3qaZeks+XX8MqURDduWCKSL5mKnUhKikZE9EQpmWIDr0XEUJUqSMqYitFUzS6ZuQMl2QE0EjVAIVCEo3AGcz3XfeWi0AFQlJk0oRlT0pWIbnWy7z9VEHJlo/o8WRpm/x61SUhyW26UmIqJaMLUzLEBJ0QkXE4gCsVtTHKQXRRHT2juw6NCqZTNFVm4VSOjgALIQHmUiKvXROJ2exYafUkchuzVE5sW4V5gLJeW93ox2J7+VoQldcjKYa61JWS0Y2GEGKCTr7j5JuVKSkRbVB6Vo8sGSkbFQHKzbqqKiNVpwRXkZFoWzg7V2fxtxQIyVY4SL2W88X2khQEhlSEY/k55QUk2qZXkKryeuhISFFKJk1CbKVkomMwLdMXmuhrgY6KiGk42qaf1DGniGkZAeIdn46MANmfxutYq2axXV9IZLLkRD5OHirSEd9WTkCi7fWPilGVkDKTlgHlJSQPSggxCUWErBR1pGhMyki0LSglI4C56Ig4r/h2NSERiM4vq6A3V1ASFBUFJ9uuGgEBqq0VI7apUtfMqToSkhUNYV0IMQ1FhKwcXZURQK0zVFkwLY+06AgQ76DTakjShCTaPpjtl30eWdGTPFTW/UkTEJ3ox/y5DUVBgO5KCKMh/cbkVBgUkQTJtAznIVlt2iQjQH2pGhXSoiOCNClJE5LTgZcrGCqSkoeKvOQJiMpopaprxVRJxQD2JEQHSgipAkWErCQ69SJdlRGgnrVqRJtZZA37jUVIEojOUzcKknmeOUOkVUbALB2vpsnJ6qwHibaVk5AsWJxKdDD1wZwiIsEi1dWi7gXx4scuLyNAvPPMkhGgudk8k+2kdZzJKMmyVC060aK5VebtKMzDsnwezadeqtCEhHAKd9JmKCJkpTFVLxJtLycj0X7FI2oA81PCZ5HWZtGIm6K5U5LIERXd5ybPI4lN8RCopmKibe2UEEZD+kuTH8zt/zcS0hFsy0i0rblVe3VJi5qUnTE2LaJSZfbZNlwfGdUoSLS9OQnRgRKymphIz7Trv9MiTMusLnWlaGzLCFBPh6tb4JpVB6F7LmkRlTqOm4fO31q13kZQNhUDmJcQ1oUQG1BECEE9KZroOM3JCFBfqqbK6Jqyz0127FUEo8r5l22jrvViuiYhhNRNtZW2VgRGS1YD1U95RTfsqp1DeofjpHZOaZ3Y6XCqtEiaoIlOPKvdtK+6n9MGsqIgbZEQHRgNIXX3iYyIECLRpsgIsDy6xER0RP5kb6NT14ksJPdt+/lWjYIA5ucJYXEqsQ1FRAFOaLZa2JIRYHlejaxUTbSveu0IoD4rq8BEJ19XnUXe8Wyfd140qqqARPtTQki/oIgQkkKdMgIsC8Zi3ZVq0ZHoGMXDfIFqUtIlbJx3USpMdUQMoCcg0XZKCOk2FBFCMtCRESB/VdiqqRqgnnTN/FiaUkKWUa3DaVsUBKCEkHbBOxAhOdQ1FXx0rOqpGiBdSPLSNdH27Lp0Sok6VeQj2q4nINFzKCGk3/CuQ0gBdcsIoJeqSdsf0K8fEVBK9FERkCz5WPxePQ0T7V9vKgaghJB2wjsNIQrUKSPR8dSjI3n769SPLH5HKSlCZ/hz3wQkOiYlhDTHat1dCKlAG2QEUE/XRM/JFpLo98WpG2C5Y+6bmOiIB1AsH9E+1etAFr+jhJD+YmxCs5deegk333wz9u/fjx07duCSSy7BPffcg+1tM6udEtIEOjdplQ5gKxxk1gDkdT5ZHVB2J5c+Idri90HmBFup7cwmTdOdPK0tlD3/omskrnNWFCRrThBKCOkSnVlr5sc//jGCIMCXvvQl/NZv/Rb+53/+B7feeiveeOMN3H///aaaJcQ4upERIH9ETXTM7OgIsFw7Ip4T/U5tuG/0nEUHWTVKEmuz5RGTKrJUJfoBlEvDRL8rLyAAJYR0B2N3i2uvvRbXXnvt/Oe3ve1teOGFF/Dggw9SREjn0V0or0qqJvpderpGPA/QE5LoeWppm2gfveBpWsffpJxUjdJUlQ+gfgEpep6gzJoxlBBik0Y/tmxsbOBNb3pT5u/H4zHG48XNenNzs4nTIqQUZWQEKJ5vBEivA8mLjojn5tWPANWjJIv99LO6puSkrtSQalrKhoAUPVdACSF1se5OGltnrTER+clPfoIHHngAf/VXf5W5z9GjR/GZz3ymqVMipDK6MgLUEx0B9NI1gqpRksV+1cUEqE8idFGVjsX+5eUjej4FhJAstO8e9957LxzHyf165plnYs85efIkrr32Wlx//fW45ZZbMo995MgRbGxszL9OnDih/xcR0jBlbuhvBCOlVXyLOrAyBa3AonCyqLi1qMh1sX+w9NUmypxb0d9fdA2jY2S/DnmvX9FzZSghpOtoR0Ruu+023HDDDbn7XHzxxfPHJ0+exNVXX433vOc9eOihh3KfNxqNMBrp/1PVQZNhKNI/xI3dVHQkOnb5CEne84tSN9Fx4p1xUcQkeo6ejKhGVUxIjopsCapEPxb7VKsDAcoJSHR8SghpF9oisnfvXuzdu1dp31deeQVXX301Dhw4gIcffhiua2y0MCGtoGyqBlAbWQOUExL5+XnHKErdLI613HGryEn+MZuJouhIh6AO+Yj2qy4gAKMgpF8YqxE5efIkrrrqKlx44YW4//778Ytf/GL+u7e85S2mmiXEOmVkBFCLjkTHryYk8jGyjpPseIvEJDpmegdfVVCqUEY6gGLxWBy//QIStUEJIfVQ9xwigEEReeyxx/Diiy/ixRdfxAUXXBD7XRjauzHpYuKik/5TJVUDFEdHojaqC4l8nLxjpXXMKnISHb+cDDSFqnQI6pAPneMIKCGkrxjLldx0000IwzD1i5BVYStcK13Mqt6GWlFrUceociwZuVizqGizLZQ9Z3Fdiq6NyrXWucaAWmFzejvl3nuENE27pj8kpKeYrB1ZtJEfIYn2iXeQKvUki32Lh9qqdOyqkRRd6hQhHVFQETzdYwKMgJDVgSIiwZEzxCRV0zVAPSmb+L6L97tOCmfxHP15QNoUOdGVg+g5evcICggh+VBECGmYskIClKshidqqV0rS2pApIyimKSMdi+e2Uz6itiggpNtQRAixRB1CAuhLSdRmviikdbyqcpLWngo68lJFKoqPXS4qWuacqghI1CYlhHQfigghlqkiJIB+LUnUpl60JHpOegetIyiq59QUZaUjem65860qH1HbFBDSPKZGkVJEErBOhNii7PwjgjJCErWrLyXx56sVwNqminQsjlH+lkkBISQdigghLaJqdAQoLyRR+9WkJDpGfodvSlTqEI3sY1NAyGpjck4tikgKjIoQ29QpJEB1KYnOpZ4CVJPCUBdV00R1yIeAEkL6DkWEkBZTh5AA1aUkOpd6hu+2jbpqUygfpK+YnmGcIkJIB6hLSIB6pERgKmpiiroLYuuUD4ACQtpFU9kBikgGTM+QNiJ3VG2TEkCtozcpK02MvKlbPgAKCFltKCKEdJQ6oyRA/VKShY1hulWhfJBVpYmFX7t3R2gQrrxLukDdURIgveM1KSdtwoR0yFBACIlDESGkR9QdJZHJ66C7JimmZSMJ5YOQbCgihPQQE1GSPHQ6dpPS0rRg5EH5IEQNigghPSerQ2xCUNJokyzUBaWDkPJQRAhZUZKdpy0x6RqUDkLqhSJCCAHQfDqnS1A+CDEHRYQQsgSjJZQPQpqCIkIIKaTvYkLpIMQeFBFCiDZpHXdX5ITSQUi7oIgQQmqhjXJC6SCk/VBEctgKhpxdlZAKUAQI6S5N9YGu8RY6Che8I4QQsuo00RdSRFKQLzyFhBBCyCpjuh+kiCSgeBBCCCFxTPaNFBEJSgghhBCS3h+a6iMpIoQQQghRwoSMdFZE1t1J7dW8HCFDCCGEZPeHJvrJzg3fNS0L4vjC+ignhBBCVh2TfWFnRKRpIaCAEEIIWXU4j8iMkTO1fQqEEELIStHUB/JOiAghhBBC+knvRYRDcgkhhJD20msRERJCGSGEEEKaYxyql6D2VkSS8kEZIYQQQsyj29/2UkSyLgJlhBBCCDFHmX62dyJSdBEoI4QQQkj9lO1feyUiqheBMkIIIYTUw1YwrNSv9kZEdC+CaRmp+sIQQgghVWiiH6rj+J0XkbIX2uRELRQQQgghbcGkkNTRl3ZaRMpe2KbWqyGEEEJsY2KR2DqP35m1ZmTaKiC22iKEEEKS2FijrUz/3DkRaVsahhBCCCERor89A0f5OZ0RkS5EQQghhBCit1htJ2pEdKaKlaGEEEIIIe2mEyJShrokhCNgCCGErCJN9X+dSc2oYkJAtoIhoyuEEEJWgqb7v15FRExGQRgZIYQQ0nds9H+9EZEmUjGcLZUQQkhfKer/TNGIiIzHY7zzne+E4zj40Y9+VPvxm64HoYwQQgjpEyr9mqkP442IyCc/+Umcf/75tR+3ztni2rZWDSGEEGKaMnJRd/9nXET+7d/+DY899hjuv//+Wo9bp4BwkjRCCCGrSNm+rE4ZMTpq5uc//zluvfVW/Mu//At27txZuP94PMZ4PJ7/vLm5mbqf7aG5lBBCCCF9QfRpZSMjVftEYxGRMAxx00034dChQ7j88suVnnP06FHs2bNn/rVv377Y7+tKxVSJglBCCCGE9BFb0RFtEbn33nvhOE7u1zPPPIMHHngAm5ubOHLkiPKxjxw5go2NjfnXiRMn5r9jFIQQQggxS9kP3FUKWZ0wDEOdJ7z++ut4/fXXc/e5+OKLccMNN+Bf//Vf4TiLhW9834fnefjYxz6GL3/5y4VtbW5uYs+ePfj7Zw9g5y5P5zSXoIAQQggh6lTpN0+f8nHr7x3HxsYGdu/enbu/toio8vLLL8dqPE6ePIkPfOAD+Od//mccPHgQF1xwQeEx6hIRSgghhBBSjjJ9aPDGlrKIGCtWvfDCC2M/n3322QCASy65RElC6sK2hHB6eEIIITaoq/9ZdyfafanOYrW9mVk1SRsKUkX7nHOEEEJIU8j9X139j8kP1I0tenfxxRfDUBZoiTZEQdK2MTJCCCHEJCb7nzKRERV6FxFpo4TIv2N0hBBCiAmaWCvGxAfqXolImyWkzH6EEEKICk2uFVP3nFq9EJE21YOY2p8QQghJYnOtmLr6z86LSBuiIFyrhhBCiA1srxVTR1/WaRFpg4TYbJ8QQgixMRtqHe0LGhs1Uye2BcB2+4QQQkiSsqNabI+q6VxExLYE2G6fEEIIyaJKdMRW+52KiNiuxbDdPiGEEKJCmeiE2L+OfmvkTJX37UREZBwOrEoAC1IJIYR0DduFrKp0QkTK0JdUDIf5EkLIamKzkLTJCTh7JyI25wYx1T5lhBBCVou2rBXTRP/TKxHpSyomrX3KCCGE9B+T9/+2ykhvRMR2KoRr1RBCCKlCm9eKMdn39EJEVkFCyuxHCCGk/ah+yLQ9AZmpD8OdFpG66jG4Vg0hhBAblJ2ArA7aEh3prIj0KQrCocGEELKa2JYB2+0DHZvQDLA/QVmd52C7fUIIIfYR93RbE5DZbr9TEZE2RCH6kgoihBDSLmxHJ2y13xkR6UsUwnb7hBBC2ksX14oR7Zc9h06IiM6c9XnYjkLYbp8QQkg3sD2qpcnoSCdEpCq2C0Jtt08IIaR7rEqqpvciYjsV0tX2CSGE2Mf2nB9l2x+H6mNhOjdqRhXbAlD2HEy0vxUMGV0hhJCGSN7767j/rruT0nOO2GxfhV5GRGxLiO1UDNeqIYQQO6Tda7uWKqm7/SJ6JyJtkJC2ts+1agghxByrsFaMiQEUvRKRNktAm9qnjBBCSH2s4loxdcpIL0TE9gRhttsX56C7P4WEENJnmrjHrfJaMXX1X50Xka5EIUy2b1uCkohzouwQQmwhF+ubvA91VQbalKrptIjY7oDbICE22887PkfpEEJsIt+D2nrPY6omopPDd213wKveflvbI4SQJE3eh7o6xNZ2+52LiNjuhFe9fUII6Qo27nuruFZM1fY7ExFpQwfclgnKbLVPCCFEjTLRAbF/H6IjOmvEdSIiojNVrEwbCkL70D4hhBB92lBIarN9VTohImXoSyqkq+0TQkjX6UMhqe32VeidiKzy3CCm2qeMEEJWCfn+35fohO328+iViNiOAvQlFZO1Vg2FhBDSd7hWjLn2s+iNiPRJQtrcPmWEENJHij5s9UUG2piq6byI2E6FiHNYpfYZHSGE9Imm1+lqgwzYFiKZzgzfTWPVBMBU+2XPoa5hXlnHzsPWiKC2CxhHSrUbvn/SUbkutu41Wfv3YYit7fYFnRWRvkjAqrefhso51fWP0PaOQZfk30MxsU+X3mO659qH/0HbnXEb2hfH021ffn4VOicitqMAttuv8xxst593/Lxzq7P95LFs3xTrgPLRLsre6NuIifeWyvVp4p5TdA5p9EkGbApRp0SkLx3wqrev2paNT/dpbbS5A6F0dIeuvbeA5v/nAbtRvVWfDdVW+50Rkb50wqvevm6bJutQdM5DFZMdi+3rQOqnaTnpwntIFhJb9x3Rvg59iY7YaL8TIhLNWe9VOobtDth2+2XPwfaNy3b7unTtfEn74HsowvZ1YHSkufY7P3xXBdsS0Ib2uyghhBBiE9tDXNvQftlhxjp0IiJShTZIwKq3T6HpLm2vYWgCvn+7SZdTFW1qXxxDt32dxWp7GxGpEgWoa0iazfbFOZSh7vY5AVq3EK8XX7MIXovuIb9efYpOdLF9FXopIm3pgG2230YJ4s28vVA+iuE1aj9Zr4/tznjV2y+id6mZNkgA2y/+PcPd7YCdajnk68b3sn2anASxDUWk4ni67cvPb7r9PHoTEWlrFED1HFatfX6ytAc/2dcLr6VddO87dVCliNN2dMJ2+2kYF5Hvfve7OHjwIHbs2IG9e/fiwx/+cO1ttKED7mI9ijiHOijbPm/gzcEO0yy8vs1S9nrX+Rp1VQZst5/EaGrmG9/4Bm699Vb8xV/8Bf7oj/4IYRjiueeeq7WNrnbAfWm/LedAipGvNzvM+uH7uVmqpAhsjygR58BUTYQxEZlOp7jjjjtw33334eabb55v/53f+Z1ajt+Gzs/23By2r0GVNx5v2nYxkeddRfg+tk9VIemDDIhj2L4GZe8nxlIzzz77LF555RW4rot3vetdeOtb34rrrrsOzz//fOZzxuMxNjc3Y19ptKEDtikhtlNB4hzKUOc5kOqI14OviR68Zu2DdRvdbd+YiPz0pz8FANx777349Kc/je985zs455xzcOWVV+JXv/pV6nOOHj2KPXv2zL/27du3tE8bJGTV27ctQcQMfI3yobR1gzZ0xjaFqIvta4vIvffeC8dxcr+eeeYZBEEAALj77rvxkY98BAcOHMDDDz8Mx3Hw9a9/PfXYR44cwcbGxvzrxIkT89/Z7gBXvX1xDjbbJ83ADjcOr0X3sN0Zi3MoQ1/aj9aIU0O7RuS2227DDTfckLvPxRdfjFOnTgEA3v72ty9ObDTC2972Nrz88supzxuNRhiNRkvbx+EAO3RPFP2KQqxy+8QefA1Jl2lD3YQ43iq2r4q2iOzduxd79+4t3O/AgQMYjUZ44YUXcMUVVwAAJpMJXnrpJVx00UX6Z6pJXzphts9iyipQJOzD93B5bHfGdZ7DKrdfhLFRM7t378ahQ4dwzz33YN++fbjoootw3333AQCuv/56U81a7wD70n4bzoE3cNJ1+B6uRp2fzNsQHVnl9vMwOo/Ifffdh8FggBtvvBFnzpzBwYMH8b3vfQ/nnHOOkfZsd4Bs374EEUL6BWWkP+1n4YRhGBo7ekU2NzexZ88efOH4Qew4O9+Z+tIJr3r7Vc6BLMPUjF34Xq6PNtxj+nQOpts/fcrHrb93HBsbG9i9e3fuvq1e9E440plfZ1ffisrc0xXbGofiUqhX+tbZ/uIc9NoX59CH9qucA0nnDByt6nVSL+PQsX0KvaHe97Iv3fP1zgHQGxFi4hy60P6ZX/sAFv14Hq2OiPzv//5v6lwihBBCCGk/J06cwAUXXJC7T6tFJAgCnDx5Ert27YLjFH+y2NzcxL59+3DixInCUBBJh9ewHngdq8NrWA+8jtXhNdQnDEOcOnUK559/Plw3f8qyVqdmXNctNKk0du/ezTdLRXgN64HXsTq8hvXA61gdXkM99uzZo7SfsSneCSGEEEKKoIgQQgghxBq9EpHRaIR77rkndZp4ogavYT3wOlaH17AeeB2rw2tollYXqxJCCCGk3/QqIkIIIYSQbkERIYQQQog1KCKEEEIIsQZFhBBCCCHW6LWIfPe738XBgwexY8cO7N27Fx/+8Idtn1InGY/HeOc73wnHcfCjH/3I9ul0ipdeegk333wz9u/fjx07duCSSy7BPffcg+3tbdun1nq++MUvYv/+/VhfX8eBAwfw/e9/3/YpdYajR4/i3e9+N3bt2oVzzz0XH/rQh/DCCy/YPq1Oc/ToUTiOgzvvvNP2qfSO3orIN77xDdx44434+Mc/jv/6r//Cf/7nf+JP//RPbZ9WJ/nkJz+J888/3/ZpdJIf//jHCIIAX/rSl/D888/jb/7mb/B3f/d3+PM//3Pbp9ZqHnnkEdx55524++678cMf/hDvfe97cd111+Hll1+2fWqd4IknnsDhw4fx1FNP4dixY5hOp7jmmmvwxhtv2D61TvL000/joYcewjve8Q7bp9JPwh4ymUzC3/zN3wz/4R/+wfapdJ5HH300vPTSS8Pnn38+BBD+8Ic/tH1Knecv//Ivw/3799s+jVbz+7//++GhQ4di2y699NLwz/7szyydUbd57bXXQgDhE088YftUOsepU6fC3/7t3w6PHTsWXnnlleEdd9xh+5R6Ry8jIs8++yxeeeUVuK6Ld73rXXjrW9+K6667Ds8//7ztU+sUP//5z3HrrbfiH//xH7Fz507bp9MbNjY28KY3vcn2abSW7e1tHD9+HNdcc01s+zXXXIMf/OAHls6q22xsbAAA33clOHz4MD74wQ/i/e9/v+1T6S29FJGf/vSnAIB7770Xn/70p/Gd73wH55xzDq688kr86le/snx23SAMQ9x00004dOgQLr/8ctun0xt+8pOf4IEHHsChQ4dsn0pref311+H7Ps4777zY9vPOOw+vvvqqpbPqLmEY4q677sIVV1yByy67zPbpdIqvfe1rePbZZ3H06FHbp9JrOiUi9957LxzHyf165plnEAQBAODuu+/GRz7yERw4cAAPP/wwHMfB17/+dct/hV1Ur+EDDzyAzc1NHDlyxPYptxLV6yhz8uRJXHvttbj++utxyy23WDrz7uA4TuznMAyXtpFibrvtNvz3f/83/umf/sn2qXSKEydO4I477sBXvvIVrK+v2z6dXjOwfQI63Hbbbbjhhhty97n44otx6tQpAMDb3/72+fbRaIS3ve1tK1/spnoNP/e5z+Gpp55aWlvh8ssvx8c+9jF8+ctfNnmarUf1OgpOnjyJq6++Gu95z3vw0EMPGT67brN37154nrcU/XjttdeWoiQkn9tvvx3f/va38eSTT+KCCy6wfTqd4vjx43jttddw4MCB+Tbf9/Hkk0/iC1/4AsbjMTzPs3iG/aFTIrJ3717s3bu3cL8DBw5gNBrhhRdewBVXXAEAmEwmeOmll3DRRReZPs1Wo3oNP//5z+Nzn/vc/OeTJ0/iAx/4AB555BEcPHjQ5Cl2AtXrCACvvPIKrr766nlkznU7FYhsnLW1NRw4cADHjh3Dn/zJn8y3Hzt2DH/8x39s8cy6QxiGuP322/HNb34Tjz/+OPbv32/7lDrH+973Pjz33HOxbR//+Mdx6aWX4lOf+hQlpEY6JSKq7N69G4cOHcI999yDffv24aKLLsJ9990HALj++ustn103uPDCC2M/n3322QCASy65hJ+sNDh58iSuuuoqXHjhhbj//vvxi1/8Yv67t7zlLRbPrN3cdddduPHGG3H55ZfPo0gvv/wya2sUOXz4ML761a/iW9/6Fnbt2jWPLu3Zswc7duywfHbdYNeuXUs1NWeddRbe/OY3s9amZnopIgBw3333YTAY4MYbb8SZM2dw8OBBfO9738M555xj+9TICvHYY4/hxRdfxIsvvrgkcCEXvs7kox/9KH75y1/is5/9LH72s5/hsssuw6OPPrryEU1VHnzwQQDAVVddFdv+8MMP46abbmr+hAjJwQl5NySEEEKIJZisJoQQQog1KCKEEEIIsQZFhBBCCCHWoIgQQgghxBoUEUIIIYRYgyJCCCGEEGtQRAghhBBiDYoIIYQQQqxBESGEEEKINSgihBBCCLEGRYQQQggh1qCIEEIIIcQa/z99L9i5HVrQcgAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Time = 0.05 Total energy = 1.98329923e-05 Linear momentum 6.42238800e-24 Angular momentum -2.92145793e-04 \n", + "Time = 0.10 Total energy = 1.98329923e-05 Linear momentum 1.31556755e-23 Angular momentum -2.80169377e-04 \n", + "Time = 0.15 Total energy = 1.98329923e-05 Linear momentum 2.04511117e-23 Angular momentum -2.61131233e-04 \n", + "Time = 0.20 Total energy = 1.98329923e-05 Linear momentum 2.83829882e-23 Angular momentum -2.36332458e-04 \n", + "Time = 0.25 Total energy = 1.98329923e-05 Linear momentum 3.67386534e-23 Angular momentum -2.07437862e-04 \n", + "Time = 0.30 Total energy = 1.98329923e-05 Linear momentum 4.48989055e-23 Angular momentum -1.76338524e-04 \n", + "Time = 0.35 Total energy = 1.98329923e-05 Linear momentum 5.17647244e-23 Angular momentum -1.45000052e-04 \n", + "Time = 0.40 Total energy = 1.98329923e-05 Linear momentum 5.58193938e-23 Angular momentum -1.15311113e-04 \n", + "Time = 0.45 Total energy = 1.98329923e-05 Linear momentum 5.53440577e-23 Angular momentum -8.89461368e-05 \n", + "Time = 0.50 Total energy = 1.98329923e-05 Linear momentum 4.88057766e-23 Angular momentum -6.72539478e-05 \n", + "Time = 0.55 Total energy = 1.98329923e-05 Linear momentum 3.53380529e-23 Angular momentum -5.11807565e-05 \n", + "Time = 0.60 Total energy = 1.98329923e-05 Linear momentum 1.52052687e-23 Angular momentum -4.12319404e-05 \n", + "Time = 0.65 Total energy = 1.98329923e-05 Linear momentum -9.90025800e-24 Angular momentum -3.74728598e-05 \n", + "Time = 0.70 Total energy = 1.98329923e-05 Linear momentum -3.68563893e-23 Angular momentum -3.95650954e-05 \n", + "Time = 0.75 Total energy = 1.98329923e-05 Linear momentum -6.14741507e-23 Angular momentum -4.68313792e-05 \n", + "Time = 0.80 Total energy = 1.98329923e-05 Linear momentum -7.91975202e-23 Angular momentum -5.83404077e-05 \n", + "Time = 0.85 Total energy = 1.98329923e-05 Linear momentum -8.60149697e-23 Angular momentum -7.30018097e-05 \n", + "Time = 0.90 Total energy = 1.98329923e-05 Linear momentum -7.93771785e-23 Angular momentum -8.96617491e-05 \n", + "Time = 0.95 Total energy = 1.98329923e-05 Linear momentum -5.89005961e-23 Angular momentum -1.07190826e-04 \n", + "Time = 1.00 Total energy = 1.98329923e-05 Linear momentum -2.66537293e-23 Angular momentum -1.24557824e-04 \n", + "\n" + ] + } + ], + "source": [ + "par = {'Compute_energy': 10,\n", + " 'plot_tstep': 10,\n", + " 'end_time': 1}\n", + "dt = 0.005\n", + "integrator = RK4(TT, N=NonlinearRHS, update=update, energy=[\"\"], **par)\n", + "integrator.setup(dt)\n", + "fu_hat = integrator.solve(fu, fu_hat, dt, (0, par['end_time']))" + ] + }, + { + "cell_type": "markdown", + "id": "57859b90", + "metadata": { + "editable": true + }, + "source": [ + "A complete solver is found [here](https://github.com/spectralDNS/shenfun/blob/master/demo/KleinGordon.py).\n", + "\n", + "\n", + "\n", + "1. **A.-M. Wazwaz**. New Travelling Wave Solutions to the Boussinesq and the Klein-Gordon Equations, *Communications in Nonlinear Science and Numerical Simulation*, 13(5), pp. 889-901, [doi: 10.1016/j.cnsns.2006.08.005](https://dx.doi.org/10.1016/j.cnsns.2006.08.005), 2008." + ] + } + ], + "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.12.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/source/kleingordon.rst b/docs/source/kleingordon.rst deleted file mode 100644 index 23dbfbee..00000000 --- a/docs/source/kleingordon.rst +++ /dev/null @@ -1,621 +0,0 @@ -.. File automatically generated using DocOnce (https://github.com/doconce/doconce/): - -.. doconce format sphinx kleingordon.do.txt --sphinx_preserve_bib_keys - -.. Document title: - -Demo - Cubic nonlinear Klein-Gordon equation -============================================ - -:Authors: Mikael Mortensen (mikaem at math.uio.no) -:Date: April 13, 2018 - -*Summary.* This is a demonstration of how the Python module `shenfun `__ can be used to solve the time-dependent, -nonlinear Klein-Gordon equation, in a triply periodic domain. The demo is implemented in -a single Python file `KleinGordon.py `__, and it may be run -in parallel using MPI. The Klein-Gordon equation is solved using a mixed -formulation. The discretization, and some background on the spectral Galerkin -method is given first, before we turn to the actual details of the ``shenfun`` -implementation. - -.. raw:: html - - -

- -Movie showing the evolution of the solution :math:`u` from Eq. :eq:`eq:kg`, in a slice through the center of the domain, computed with the code described in this demo. - -The nonlinear Klein-Gordon equation ------------------------------------ - -The cubic nonlinear Klein-Gordon equation is a wave equation important for many -scientific applications such as solid state physics, nonlinear optics and -quantum field theory :cite:`abdul08`. The equation is given as - -.. math:: - :label: eq:kg - - - \frac{\partial^2 u}{\partial t^2} = \nabla^2 u - \gamma(u - u|u|^2) \quad - \text{for} \, u \in - \Omega, - - -with initial conditions - -.. math:: - :label: eq:init - - - u(\boldsymbol{x}, t=0) = u^0 \quad \text{and} \quad \frac{\partial u(\boldsymbol{x}, - t=0)}{\partial t} = u_t^0. - - -The spatial coordinates are here denoted as :math:`\boldsymbol{x} = (x, y, z)`, and -:math:`t` is time. The parameter :math:`\gamma=\pm 1` determines whether the equations are focusing -(:math:`+1`) or defocusing (:math:`-1`) (in the movie we have used :math:`\gamma=1`). -The domain :math:`\Omega=[-2\pi, 2\pi)^3` is triply -periodic and initial conditions will here be set as - -.. math:: - :label: _auto1 - - - u^0 = 0.1 \exp \left( -\boldsymbol{x} \cdot \boldsymbol{x} \right), - - - -.. math:: - :label: _auto2 - - - u_t^0 = 0. - - - -We will solve these equations using a mixed formulation and a spectral Galerkin -method. The mixed formulation reads - -.. math:: - :label: eq:df - - - \frac{\partial f}{\partial t} = \nabla^2 u - \gamma (u - u|u|^2), - - -.. math:: - :label: eq:du - - - \frac{\partial u}{\partial t} = f. - - -The energy of the solution can be computed as - -.. math:: - :label: _auto3 - - - 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 - - - -and it is crucial that this energy remains constant in time. - -The movie above is showing the solution :math:`u`, computed with the -code shown below. - -.. _sec:specgal: - -Spectral Galerkin formulation ------------------------------ -The PDEs in :eq:`eq:df` and :eq:`eq:du` can be solved with many different -numerical methods. We will here use the `shenfun `__ software and this software makes use of -the spectral Galerkin method. Being a Galerkin method, we need to reshape the -governing equations into proper variational forms, and this is done by -multiplying :eq:`eq:df` and :eq:`eq:du` with the complex conjugate of proper -test functions and then integrating -over the domain. To this end we make use of the triply periodic tensor product -function space :math:`W^{\boldsymbol{N}}(\Omega)` (defined in Eq. :eq:`eq:kg:Wn`) -and use testfunctions :math:`g \in W^{\boldsymbol{N}}` -with Eq. :eq:`eq:df` and :math:`v \in W^{\boldsymbol{N}}` with Eq. :eq:`eq:du`, -and obtain - -.. math:: - :label: eq:df_var - - - \frac{\partial}{\partial t} \int_{\Omega} f\, \overline{g}\, w \,d\Omega = \int_{\Omega} - \left(\nabla^2 u - \gamma( u\, - u|u|^2) \right) \overline{g} \, w \,d\Omega, - - -.. math:: - :label: eq:kg:du_var - - - \frac{\partial }{\partial t} \int_{\Omega} u\, \overline{v}\, w \, dx = - \int_{\Omega} f\, \overline{v} \, w \, d\Omega. - - -Note that the overline is used to indicate a complex conjugate, and -:math:`w` is a weight function associated with the test functions. The functions -:math:`f` and :math:`u` are now -to be considered as trial functions, and the integrals over the -domain are referred to as inner products. With inner product notation - -.. math:: - - \left(u, v\right) = \int_{\Omega} u \, \overline{v} \, w\, dx. - - -and an integration by parts on the Laplacian, the variational problem can be -formulated as: - -.. math:: - :label: eq:df_var2 - - - \frac{\partial}{\partial t} (f, g) = -(\nabla u, \nabla g) - -\gamma \left( u - u|u|^2, g \right), - - -.. math:: - :label: eq:kg:du_var2 - - - \frac{\partial }{\partial t} (u, v) = (f, v). - - -The time and space discretizations are -still left open. There are numerous different approaches that one could take for -discretizing in time, and the first two terms on the right hand side of -:eq:`eq:df_var2` can easily be treated implicitly as well as explicitly. However, -the approach we will follow in Sec. (:ref:`sec:rk`) is a fully explicit 4th order `Runge-Kutta `__ method. Also note that -the inner product in the demo will be computed numerically with quadrature -through fast Fourier transforms, and the integrals are thus not computed exactly -for all terms. - -Discretization --------------- -To find a numerical solution we need to discretize the continuous problem -:eq:`eq:df_var2` and :eq:`eq:kg:du_var2` in space as well as time. Since the -problem is triply periodic, Fourier exponentials are normally the best choice -for trial and test functions, and as such we use basis functions - -.. math:: - :label: _auto4 - - - \phi_l(x) = e^{\imath \underline{l} x}, \quad -\infty < l < \infty, - - - -where :math:`l` is the wavenumber, and -:math:`\underline{l}=\frac{2\pi}{L}l` is the scaled wavenumber, scaled with domain -length :math:`L` (here :math:`4\pi`). Since we want to solve these equations on a computer, we need to choose -a finite number of test functions. A function space :math:`V^N` can be defined as - -.. math:: - :label: eq:kg:Vn - - - V^N(x) = \text{span} \{\phi_l(x)\}_{l\in \boldsymbol{l}}, - - -where :math:`N` is chosen as an even positive integer and :math:`\boldsymbol{l} = -N/2, --N/2+1, \ldots, N/2-1`. And now, since :math:`\Omega` is a -three-dimensional domain, we can create tensor products of such bases to get, -e.g., for three dimensions - -.. math:: - :label: eq:kg:Wn - - - W^{\boldsymbol{N}}(x, y, z) = V^N(x) \otimes V^N(y) \otimes V^N(z), - - -where :math:`\boldsymbol{N} = (N, N, N)`. Obviously, it is not necessary to use the -same number (:math:`N`) of basis functions for each direction, but it is done here -for simplicity. A 3D tensor product basis function is now defined as - -.. math:: - :label: _auto5 - - - \Phi_{lmn}(x,y,z) = e^{\imath \underline{l} x} e^{\imath \underline{m} y} - e^{\imath \underline{n} z} = e^{\imath - (\underline{l}x + \underline{m}y + \underline{n}z)}, - - - -where the indices for :math:`y`- and :math:`z`-direction are :math:`\underline{m}=\frac{2\pi}{L}m, -\underline{n}=\frac{2\pi}{L}n`, and :math:`\boldsymbol{m}` and :math:`\boldsymbol{n}` are the same as -:math:`\boldsymbol{l}` due to using the same number of basis functions for each direction. One -distinction, though, is that for the :math:`z`-direction expansion coefficients are only stored for -:math:`n=0, 1, \ldots, N/2` due to Hermitian symmetry (real input data). However, for simplicity, -we still write the sum in Eq. :eq:`eq:usg` over the entire range of basis functions. - -We now look for solutions of the form - -.. math:: - :label: eq:usg - - - 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} - \hat{u}_{lmn} (t)\Phi_{lmn}(x,y,z). - - -The expansion coefficients :math:`\hat{\boldsymbol{u}} = \{\hat{u}_{lmn}(t)\}_{(l,m,n) \in \boldsymbol{l} \times \boldsymbol{m} \times \boldsymbol{n}}` -can be related directly to the solution :math:`u(x, y, z, t)` using Fast -Fourier Transforms (FFTs) if we are satisfied with obtaining -the solution in quadrature points corresponding to - -.. math:: - :label: _auto6 - - - x_i = \frac{4 \pi i}{N}-2\pi \quad \forall \, i \in \boldsymbol{i}, - \text{where}\, \boldsymbol{i}=(0,1,\ldots,N-1), - - - -.. math:: - :label: _auto7 - - - y_j = \frac{4 \pi j}{N}-2\pi \quad \forall \, j \in \boldsymbol{j}, - \text{where}\, \boldsymbol{j}=(0,1,\ldots,N-1), - - - -.. math:: - :label: _auto8 - - - z_k = \frac{4 \pi k}{N}-2\pi \quad \forall \, k \in \boldsymbol{k}, - \text{where}\, \boldsymbol{k}=(0,1,\ldots,N-1). - - - -Note that these points are different from the standard (like :math:`2\pi j/N`) since -the domain -is set to :math:`[-2\pi, 2\pi]^3` and not the more common :math:`[0, 2\pi]^3`. We have - -.. math:: - :label: eq:uxyz - - - \boldsymbol{u} = \mathcal{F}_x^{-1}\left(\mathcal{F}_y^{-1}\left(\mathcal{F}_z^{-1}\left(\hat{\boldsymbol{u}}\right)\right)\right) - - -with :math:`\boldsymbol{u} = \{u(x_i, y_j, z_k)\}_{(i,j,k)\in \boldsymbol{i} \times \boldsymbol{j} \times \boldsymbol{k}}` -and where :math:`\mathcal{F}_x^{-1}` is the inverse Fourier transform along the direction :math:`x`, for -all indices in the other direction. Note that the three -inverse FFTs are performed sequentially, one direction at the time, and that there is no -scaling factor due to -the definition used for the inverse `Fourier transform `__ - -.. math:: - :label: _auto9 - - - u(x_j) = \sum_{l=-N/2}^{N/2-1} \hat{u}_l e^{\imath \underline{l} - x_j}, \quad \,\, \forall \, j \in \, \boldsymbol{j}. - - - -Note that this differs from the definition used by, e.g., -`Numpy `__. - -The inner products used in Eqs. :eq:`eq:df_var2`, :eq:`eq:kg:du_var2` may be -computed using forward FFTs. However, there is a tiny detail that deserves -a comment. The regular Fourier inner product is given as - -.. math:: - \int_{0}^{L} e^{\imath \underline{k}x} e^{- \imath \underline{l}x} dx = L\, \delta_{kl} - -where a weight function is chosen as :math:`w(x) = 1` and :math:`\delta_{kl}` equals unity -for :math:`k=l` and zero otherwise. In Shenfun we choose instead to use a weight -function :math:`w(x)=1/L`, such that the weighted inner product integrates to -unity: - -.. math:: - \int_{0}^{L} e^{\imath \underline{k}x} e^{- \imath \underline{l}x} \frac{1}{L} dx = \delta_{kl}. - -With this weight function the (discrete) scalar product and the forward transform -are the same and we obtain: - -.. math:: - :label: _auto10 - - - \left(u, v \right) = \boldsymbol{\hat{u}} = - \left(\frac{1}{N}\right)^3 - \mathcal{F}_z\left(\mathcal{F}_y\left(\mathcal{F}_x\left(\boldsymbol{u}\right)\right)\right). - - - -From this we see that the variational forms :eq:`eq:df_var2` and :eq:`eq:kg:du_var2` -may be written in terms of the Fourier transformed quantities :math:`\hat{\boldsymbol{u}}` and -:math:`\hat{\boldsymbol{f}}`. Expanding the exact derivatives of the nabla operator, we have - -.. math:: - :label: _auto11 - - - (\nabla u, \nabla v) = - \left((\underline{l}^2+\underline{m}^2+\underline{n}^2)\hat{u}_{lmn}\right), - - - -.. math:: - :label: _auto12 - - - (u, v) = \left(\hat{u}_{lmn}\right), - - - -.. math:: - :label: _auto13 - - - (u|u|^2, v) = \left(\widehat{u|u|^2}_{lmn}\right) - - - -where the indices on the right hand side run over :math:`(l, m, n) \in \boldsymbol{l} \times \boldsymbol{m} \times \boldsymbol{n}`. -The equations to be solved for each wavenumber can now be found directly as - -.. math:: - :label: eq:df_var3 - - - \frac{\partial \hat{f}_{lmn}}{\partial t} = - \left(-(\underline{l}^2+\underline{m}^2+\underline{n}^2+\gamma)\hat{u}_{lnm} + \gamma \widehat{u|u|^2}_{lnm}\right), - - -.. math:: - :label: eq:kg:du_var3 - - - \frac{\partial \hat{u}_{lnm}}{\partial t} = \hat{f}_{lnm}. - - -There is more than one way to arrive at these equations. Taking the 3D Fourier -transform of both equations :eq:`eq:df` and :eq:`eq:du` is one obvious way. -With the Python module `shenfun `__, one can work with the -inner products as seen in :eq:`eq:df_var2` and :eq:`eq:kg:du_var2`, or the Fourier -transforms directly. See for example Sec. :ref:`sec:rk` for how :math:`(\nabla u, \nabla -v)` can be -implemented. In short, :mod:`.shenfun` contains all the tools required to work with -the spectral Galerkin method, and we will now see how :mod:`.shenfun` can be used to solve -the Klein-Gordon equation. - -For completion, we note that the discretized problem to solve can be formulated -with the Galerkin method as: -for all :math:`t>0`, find :math:`(f, u) \in W^N \times W^N` such that - -.. math:: - :label: eq:dff - - - \frac{\partial}{\partial t} (f, g) = -(\nabla u, \nabla g) - -\gamma \left( u - u|u|^2, g \right) \quad \forall \, g \in W^{N}, - - -.. math:: - :label: eq:kg:duu - - - \frac{\partial }{\partial t} (u, v) = (f, v) \quad \forall \, v \in W^N. - - -where :math:`u(x, y, z, 0)` and :math:`f(x, y, z, 0)` are given as the initial conditions -according to Eq. :eq:`eq:init`. - -Implementation --------------- - -To solve the Klein-Gordon equations we need to make use of the Fourier function -spaces defined in -:mod:`.shenfun`, and these are found in submodule -:mod:`shenfun.fourier.bases`. -The triply periodic domain allows for Fourier in all three directions, and we -can as such create one instance of this space using :func:`.FunctionSpace` with -family ``Fourier`` -for each direction. However, since the initial data are real, we -can take advantage of Hermitian symmetries and thus make use of a -real to complex class for one (but only one) of the directions, by specifying -``dtype='d'``. We can only make use of the -real-to-complex class for the direction that we choose to transform first with the forward -FFT, and the reason is obviously that the output from a forward transform of -real data is now complex. We may start implementing the solver as follows - -.. code-block:: python - - from shenfun import * - import numpy as np - import sympy as sp - - # Set size of discretization - N = (32, 32, 32) - - # Defocusing or focusing - gamma = 1 - - rank = comm.Get_rank() - - # Create function spaces - K0 = FunctionSpace(N[0], 'F', domain=(-2*np.pi, 2*np.pi), dtype='D') - K1 = FunctionSpace(N[1], 'F', domain=(-2*np.pi, 2*np.pi), dtype='D') - K2 = FunctionSpace(N[2], 'F', domain=(-2*np.pi, 2*np.pi), dtype='d') - -We now have three instances ``K0``, ``K1`` and ``K2``, corresponding to the space -:eq:`eq:kg:Vn`, that each can be used to solve -one-dimensional problems. However, we want to solve a 3D problem, and for this -we need a tensor product space, like :eq:`eq:kg:Wn`, created as a tensor -product of these three spaces - -.. code-block:: python - - T = TensorProductSpace(comm, (K0, K1, K2), **{'planner_effort': - 'FFTW_MEASURE'}) - -Here the ``planner_effort``, which is a flag used by `FFTW `__, is optional. Possibel choices are from the list -(``FFTW_ESTIMATE``, ``FFTW_MEASURE``, ``FFTW_PATIENT``, ``FFTW_EXHAUSTIVE``), and the -flag determines how much effort FFTW puts in looking for an optimal algorithm -for the current platform. Note that it is also possible to use FFTW `wisdom `__ with -``shenfun``, and as such, for production, one may perform exhaustive planning once -and then simply import the result of that planning later, as wisdom. - -The :class:`.TensorProductSpace` instance ``T`` contains pretty much all we need for -computing inner products or fast transforms between real and wavenumber space. -However, since we are going to solve for a mixed system, it is convenient to also use the -:class:`.CompositeSpace` class - -.. code-block:: python - - TT = CompositeSpace([T, T]) - TV = VectorSpace(T) - -Here the space ``TV`` will be used to compute gradients, which -explains why it is a vector. - -We need containers for the solution as well as intermediate work arrays for, -e.g., the Runge-Kutta method. Arrays are created using -`Sympy `__ for -initialization. Below ``f`` is initialized to 0, -whereas ``u = 0.1*sp.exp(-(x**2 + y**2 + z**2))``. - -.. code-block:: python - - # Use sympy to set up initial condition - x, y, z = sp.symbols("x,y,z", real=True) - ue = 0.1*sp.exp(-(x**2 + y**2 + z**2)) - - fu = Array(TT, buffer=(0, ue)) # Solution array in physical space - f, u = fu # Split solution array by creating two views u and f - dfu = Function(TT) # Array for right hand sides - df, du = dfu # Split into views - Tp = T.get_dealiased((1.5, 1.5, 1.5)) - up = Array(Tp) # Work array - - fu_hat = Function(TT) # Solution in spectral space - fu_hat = fu.forward() - f_hat, u_hat = fu_hat - - gradu = Array(TV) # Solution array for gradient - -The :class:`.Array` class is a subclass of Numpy's `ndarray `__, -without much more functionality than constructors that return arrays of the -correct shape according to the basis used in the construction. The -:class:`.Array` represents the left hand side of :eq:`eq:usg`, -evaluated on the quadrature mesh. A different type -of array is returned by the :class:`.Function` -class, that subclasses both Nympy's ndarray as well as an internal -:class:`.BasisFunction` -class. An instance of the :class:`.Function` represents the entire -spectral Galerkin function :eq:`eq:usg`. - -.. _sec:rk: - -Runge-Kutta integrator -~~~~~~~~~~~~~~~~~~~~~~ - -We use an explicit fourth order Runge-Kutta integrator, -imported from :mod:`shenfun.utilities.integrators`. The solver -requires one function to compute nonlinear terms, -and one to compute linear. But here we will make -just one function that computes both, and call it -``NonlinearRHS``: - -.. code-block:: python - - uh = TrialFunction(T) - vh = TestFunction(T) - L = inner(grad(vh), -grad(uh)) + [inner(vh, -gamma*uh)] - L = la.SolverDiagonal(L).mat.scale - - def NonlinearRHS(self, fu, fu_hat, dfu_hat, **par): - global count, up - dfu_hat.fill(0) - f_hat, u_hat = fu_hat - df_hat, du_hat = dfu_hat - up = Tp.backward(u_hat, up) - df_hat = Tp.forward(gamma*up**3, df_hat) - df_hat += L*u_hat - du_hat[:] = f_hat - return dfu_hat - -Note that ``L`` now is an array that represents the linear -coefficients in :eq:`eq:kg:du_var3`. - -All that is left is to write a function that is called -on each time step, which will allow us to store intermediate -solutions, compute intermediate energies, and plot -intermediate solutions. Since we will plot the same plot -many times, we create the figure first, and then simply update -the plotted arrays in the ``update`` function. - -.. code-block:: python - - import matplotlib.pyplot as plt - X = T.local_mesh(True) - if rank == 0: - plt.figure() - image = plt.contourf(X[1][..., 0], X[0][..., 0], u[..., N[2]//2], 100) - plt.draw() - plt.pause(1e-6) - -The actual ``update`` function is - -.. code-block:: python - - # Get wavenumbers - K = np.array(T.local_wavenumbers(True, True, True)) - - def update(self, fu, fu_hat, t, tstep, **params): - global gradu - - transformed = False - if rank == 0 and tstep % params['plot_tstep'] == 0 and params['plot_tstep'] > 0: - fu = fu_hat.backward(fu) - f, u = fu[:] - image.axes.clear() - image.axes.contourf(X[1][..., 0], X[0][..., 0], u[..., N[2]//2], 100) - plt.pause(1e-6) - transformed = True - - if tstep % params['Compute_energy'] == 0: - if transformed is False: - fu = fu_hat.backward(fu) - f, u = fu - f_hat, u_hat = fu_hat - ekin = 0.5*energy_fourier(f_hat, T) - es = 0.5*energy_fourier(1j*(K*u_hat), T) - eg = gamma*np.sum(0.5*u**2 - 0.25*u**4)/np.prod(np.array(N)) - eg = comm.allreduce(eg) - gradu = TV.backward(1j*(K[0]*u_hat[0]+K[1]*u_hat[1]+K[2]*u_hat[2]), gradu) - ep = comm.allreduce(np.sum(f*gradu)/np.prod(np.array(N))) - ea = comm.allreduce(np.sum(np.array(X)*(0.5*f**2 + 0.5*gradu**2 - (0.5*u**2 - 0.25*u**4)*f))/np.prod(np.array(N))) - if rank == 0: - print("Time = %2.2f Total energy = %2.8e Linear momentum %2.8e Angular momentum %2.8e" %(t, ekin+es+eg, ep, ea)) - comm.barrier() - -With all functions in place, the actual integrator -can be created and called as - -.. code-block:: python - - par = {'Compute_energy': 10, - 'plot_tstep': 10, - 'end_time': 1} - dt = 0.005 - integrator = RK4(TT, N=NonlinearRHS, update=update, **par) - integrator.setup(dt) - fu_hat = integrator.solve(fu, fu_hat, dt, (0, par['end_time'])) - -A complete solver is found `here `__. - -.. ======= Bibliography =======