diff --git a/pySDC/implementations/problem_classes/AdvectionEquation_1D_FD.py b/pySDC/implementations/problem_classes/AdvectionEquation_1D_FD.py index 7b34aa79f9..ba53e165b0 100644 --- a/pySDC/implementations/problem_classes/AdvectionEquation_1D_FD.py +++ b/pySDC/implementations/problem_classes/AdvectionEquation_1D_FD.py @@ -167,12 +167,13 @@ def solve_system(self, rhs, factor, u0, t): me[:] = L.solve(rhs) return me - def u_exact(self, t): + def u_exact(self, t, u_init=None, t_init=None): """ Routine to compute the exact solution at time t Args: t (float): current time + u_init, t_init: unused parameters for common interface reasons Returns: dtype_u: exact solution diff --git a/pySDC/implementations/sweeper_classes/Runge_Kutta.py b/pySDC/implementations/sweeper_classes/Runge_Kutta.py new file mode 100644 index 0000000000..f3ee9d9315 --- /dev/null +++ b/pySDC/implementations/sweeper_classes/Runge_Kutta.py @@ -0,0 +1,217 @@ +import numpy as np +import logging + +from pySDC.core.Sweeper import _Pars +from pySDC.core.Errors import ParameterError +from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit + + +class ButcherTableau(object): + def __init__(self, weights, nodes, matrix): + """ + Initialization routine for an collocation object + + Args: + weights (numpy.ndarray): Butcher tableau weights + nodes (numpy.ndarray): Butcher tableau nodes + matrix (numpy.ndarray): Butcher tableau entries + """ + # check if the arguments have the correct form + if type(matrix) != np.ndarray: + raise ParameterError('Runge-Kutta matrix needs to be supplied as a numpy array!') + elif len(np.unique(matrix.shape)) != 1 or len(matrix.shape) != 2: + raise ParameterError('Runge-Kutta matrix needs to be a square 2D numpy array!') + + if type(weights) != np.ndarray: + raise ParameterError('Weights need to be supplied as a numpy array!') + elif len(weights.shape) != 1: + raise ParameterError(f'Incompatible dimension of weights! Need 1, got {len(weights.shape)}') + elif len(weights) != matrix.shape[0]: + raise ParameterError(f'Incompatible number of weights! Need {matrix.shape[0]}, got {len(weights)}') + + if type(nodes) != np.ndarray: + raise ParameterError('Nodes need to be supplied as a numpy array!') + elif len(nodes.shape) != 1: + raise ParameterError(f'Incompatible dimension of nodes! Need 1, got {len(nodes.shape)}') + elif len(nodes) != matrix.shape[0]: + raise ParameterError(f'Incompatible number of nodes! Need {matrix.shape[0]}, got {len(nodes)}') + + # Set number of nodes, left and right interval boundaries + self.num_nodes = matrix.shape[0] + 1 + self.tleft = 0. + self.tright = 1. + + self.nodes = np.append(np.append([0], nodes), [1]) + self.weights = weights + self.Qmat = np.zeros([self.num_nodes + 1, self.num_nodes + 1]) + self.Qmat[1:-1, 1:-1] = matrix + self.Qmat[-1, 1:-1] = weights # this is for computing the solution to the step from the previous stages + + self.left_is_node = True + self.right_is_node = self.nodes[-1] == self.tright + + # compute distances between the nodes + if self.num_nodes > 1: + self.delta_m = self.nodes[1:] - self.nodes[:-1] + else: + self.delta_m = np.zeros(1) + self.delta_m[0] = self.nodes[0] - self.tleft + + # check if the RK scheme is implicit + self.implicit = any([matrix[i, i] != 0 for i in range(self.num_nodes - 1)]) + + +class RungeKutta(generic_implicit): + """ + Runge-Kutta scheme that fits the interface of a sweeper. + Actually, the sweeper idea fits the Runge-Kutta idea when using only lower triangular rules, where solutions + at the nodes are succesively computed from earlier nodes. However, we only perform a single iteration of this. + + We have two choices to realise a Runge-Kutta sweeper: We can choose Q = Q_Delta = , but in this + implementation, that would lead to a lot of wasted FLOPS from integrating with Q and then with Q_Delta and + subtracting the two. For that reason, we built this new sweeper, which does not have a preconditioner. + + This class only supports lower triangular Butcher tableaus such that the system can be solved with forward + subsitution. In this way, we don't get the maximum order that we could for the number of stages, but computing the + stages is much cheaper. In particular, if the Butcher tableaus is strictly lower trianglar, we get an explicit + method, which does not require us to solve a system of equations to compute the stages. + + Attribues: + butcher_tableau (ButcherTableau): Butcher tableau for the Runge-Kutta scheme that you want + """ + + def __init__(self, params): + """ + Initialization routine for the custom sweeper + + Args: + params: parameters for the sweeper + """ + # set up logger + self.logger = logging.getLogger('sweeper') + + essential_keys = ['butcher_tableau'] + for key in essential_keys: + if key not in params: + msg = 'need %s to instantiate step, only got %s' % (key, str(params.keys())) + self.logger.error(msg) + raise ParameterError(msg) + + self.params = _Pars(params) + + if 'collocation_class' in params or 'num_nodes' in params: + self.logger.warning('You supplied parameters to setup a collocation problem to the Runge-Kutta sweeper. \ +Please be aware that they are ignored since the quadrature matrix is entirely determined by the Butcher tableau.') + self.coll = params['butcher_tableau'] + + if not self.coll.right_is_node and not self.params.do_coll_update: + self.logger.warning('we need to do a collocation update here, since the right end point is not a node. ' + 'Changing this!') + self.params.do_coll_update = True + + # This will be set as soon as the sweeper is instantiated at the level + self.__level = None + + self.parallelizable = False + self.QI = self.coll.Qmat + + def update_nodes(self): + """ + Update the u- and f-values at the collocation nodes + + Returns: + None + """ + + # get current level and problem description + L = self.level + P = L.prob + + # only if the level has been touched before + assert L.status.unlocked + assert L.status.sweep <= 1, "RK schemes are direct solvers. Please perform only 1 iteration!" + + # get number of collocation nodes for easier access + M = self.coll.num_nodes + + for m in range(0, M): + # build rhs, consisting of the known values from above and new values from previous nodes (at k+1) + rhs = L.u[0] + for j in range(1, m + 1): + rhs += L.dt * self.QI[m + 1, j] * L.f[j] + + # implicit solve with prefactor stemming from the diagonal of Qd + if self.coll.implicit: + L.u[m + 1] = P.solve_system(rhs, L.dt * self.QI[m + 1, m + 1], L.u[m + 1], + L.time + L.dt * self.coll.nodes[m]) + else: + L.u[m + 1] = rhs + # update function values + L.f[m + 1] = P.eval_f(L.u[m + 1], L.time + L.dt * self.coll.nodes[m]) + + # indicate presence of new values at this level + L.status.updated = True + + return None + + +class RK1(RungeKutta): + def __init__(self, params): + implicit = params.get('implicit', False) + nodes = np.array([0.]) + weights = np.array([1.]) + if implicit: + matrix = np.array([[1.], ]) + else: + matrix = np.array([[0.], ]) + params['butcher_tableau'] = ButcherTableau(weights, nodes, matrix) + super(RK1, self).__init__(params) + + +class CrankNicholson(RungeKutta): + ''' + Implicit Runge-Kutta method of second order + ''' + def __init__(self, params): + nodes = np.array([0, 1]) + weights = np.array([0.5, 0.5]) + matrix = np.zeros((2, 2)) + matrix[1, 0] = 0.5 + matrix[1, 1] = 0.5 + params['butcher_tableau'] = ButcherTableau(weights, nodes, matrix) + super(CrankNicholson, self).__init__(params) + + +class MidpointMethod(RungeKutta): + ''' + Runge-Kutta method of second order + ''' + def __init__(self, params): + implicit = params.get('implicit', False) + if implicit: + nodes = np.array([0.5]) + weights = np.array([1]) + matrix = np.zeros((1, 1)) + matrix[0, 0] = 1. / 2. + else: + nodes = np.array([0, 0.5]) + weights = np.array([0, 1]) + matrix = np.zeros((2, 2)) + matrix[1, 0] = 0.5 + params['butcher_tableau'] = ButcherTableau(weights, nodes, matrix) + super(MidpointMethod, self).__init__(params) + + +class RK4(RungeKutta): + ''' + Explicit Runge-Kutta of fourth order: Everybodies darling. + ''' + def __init__(self, params): + nodes = np.array([0, 0.5, 0.5, 1]) + weights = np.array([1., 2., 2., 1.]) / 6. + matrix = np.zeros((4, 4)) + matrix[1, 0] = 0.5 + matrix[2, 1] = 0.5 + matrix[3, 2] = 1. + params['butcher_tableau'] = ButcherTableau(weights, nodes, matrix) + super(RK4, self).__init__(params) diff --git a/pySDC/playgrounds/Runge_Kutta/Runge-Kutta-Methods.ipynb b/pySDC/playgrounds/Runge_Kutta/Runge-Kutta-Methods.ipynb new file mode 100644 index 0000000000..f382157115 --- /dev/null +++ b/pySDC/playgrounds/Runge_Kutta/Runge-Kutta-Methods.ipynb @@ -0,0 +1,144 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "4473b4af", + "metadata": {}, + "source": [ + "# Runge-Kutta Methods\n", + "Runge-Kutta (RK) methods are single step time marching schemes that use a quadrature rule to compute cheap intermediate solutions and from a combination of these, they compute a final, high order solution to the step.\n", + "Sounds familiar?\n", + "Indeed, there are a lot of similarities between RK schemes and spectral deferred corrections (SDC)!\n", + "\n", + "So what are the differences?\n", + "The most pronounced difference is that RK methods are direct instead of iterative, but SDC is a direct method as well, if you always perform the same number of iterations.\n", + "Actually, in the special case of a single iteration and without preconditioning, SDC is an RK method!\n", + "\n", + "Significant differences arise in practice from how people deal with the fact that the quadrature matrix is dense, making a direct solution computationally expensive.\n", + "In SDC, we take a preconditioner such that we never have to solve the full dense system but iterate with cheaper problems that get us closer to the solution of the dense system.\n", + "In RK schemes, however, people usually use different quadrature rules that rely on more intermediate solutions (stages in the RK jargon) to reach the same order but which are lower diagonal like the preconditioners in SDC.\n", + "Now we can solve the systems easily again with forward substitution.\n", + "In fact, some of the most popular RK schemes are explicit, meaning the quadrature rules are striclty lower triangular and we don't have to solve anything at all!\n", + "\n", + "## Implementation as a sweeper in pySDC\n", + "The fact that there are so many similarities between the RK schemes that people use in practice and SDC means we can easily misuse the sweeper module and change some parameters to get RK behaviour within the same framework we know and love.\n", + "\n", + "The sweeper does the same thing as other sweepers: It sweeps through the collocation nodes and updates the solutions at the nodes, but it completely forgoes the preconditioner.\n", + "\n", + "A general lower triangular RK scheme is usually written like this:\n", + "$$u_{n+1} = u_n + \\Delta t \\sum_{i=1}^s b_i k_i,$$\n", + "with the stages\n", + "$$k_i = f\\left(u_n+\\Delta t\\sum_{j=1}^i a_{ij}k_j, t_n+c_i\\Delta t\\right),$$\n", + "and the nodes $c_i$.\n", + "$f$, $u$, and $\\Delta t$ follow the usual definitions of right hand side, solution and step size.\n", + "$s$ is the number of stages of the scheme.\n", + "Defining the method are the Runge-Kutta matrix $a_{ij}$ and the weights $b_i$.\n", + "These are written in a Butcher tableau.\n", + "\n", + "We are familiar with nodes.\n", + "We have them in SDC the same way as here.\n", + "What is the connection between $a_{ij}$ and $b_i$ and the quadrature matrix $q_{ij}$?\n", + "In SDC, we had a collocation matrix one larger in each direction than the number of nodes to store the initial conditions.\n", + "Here, we generate a quadrature matrix two larger in both directions to store the initial conditions as well as the final solution.\n", + "That means the last row of the quadrature matrix will carry the weights: $q_{s+1,j+1}=b_j$ and we just put the rule to compute intermediate stages in the rest of the quadrature matrix: $q_{ij}=a_{ij}$ for $1 < i,j\\leq s$.\n", + "\n", + "Don't worry about the shifted indices.\n", + "There is no interesting maths here, it's just needed to conform to the decision to add a row and column for the initial conditinons in the regular sweepers.\n", + "We could now actually set the preconditioner to the same as the quadrature matrix and use existing sweepers, but to avoid overhead by adding a lot of stuff and subtracting the same stuff again, we build a new sweeper with no preconditioner.\n", + "\n", + "## Practical tests\n", + "We confirm that this works by using a few popular RK schemes.\n", + "Namely, we try implicit Euler (RK1), implicit midpoint method, Crank-Nicholson and explicit RK4 on the van der Pol problem:" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "d710dbc2", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEMCAYAAAA8vjqRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAABF0UlEQVR4nO3deVxVdfrA8c9XUBTct1IRRUVkB1cw93JJzKWsXFqtzCYzK1P7WWnTlGbLlKNNZaO2TJo2LaamZmlOTS4ouHABN1BxAxRRBGT7/v447IuCAudeeN6vFy+9Z33AK8895/s9z6O01gghhBDlUcvsAIQQQtgeSR5CCCHKTZKHEEKIcpPkIYQQotwkeQghhCg3SR5CCCHKzd7sAKpC8+bNdfv27c0OQwghbMqePXsStNYtSlpXI5JH+/btCQkJMTsMIYSwKUqp46Wtk9tWQgghyk2ShxBCiHKT5CGEEKLcJHkIIYQoN0keQgghyk2ShxBCiHKzyam6SqnRQDDQEliitd5sbkRCCGFdLl+GsDCws4PevSv++FWePJRSy4ARQJzW2rvA8mHAB4Ad8KnWekFpx9Bafw98r5RqArwDSPIQQtRY8fEQGpr/tXcvHDkCWsPQobBxY8Wf04wrjxXAYuDz3AVKKTtgCTAYiAV2K6XWYiSS+UX2n6S1jsv5+8s5+wkhRLWnNcTGGskhN0mEhhrLcrVvDwEB8OCD0LWr8VUZqjx5aK23K6XaF1ncEziitT4GoJRaBYzSWs/HuEopRCmlgAXAT1rrvSWdRyk1GZgM4OLiUnHfgBBCVIHsbOPqoWiiOH/eWK8UdOkC/foZCSIgAPz9oWnTqonPWsY82gAnC7yOBXpdY/tngDuARkqpTlrrj4puoLX+BPgEoHv37tJrVwhhtTIywGLJTxChocZ4RXKysb5OHfD2htGj8xOFry84OZkXs7UkD1XCslJ/4WutFwGLKi8cIYSoHCkpsH9/4SuKgwchPd1Y7+RkXEE88kh+ovD0NBKINbGW5BELtC3w2hk4bVIsQghRIRITjSuIgokiKsq4JQXQrJmRHKZPN/4MCAA3N6hlAw9RWEvy2A24KaVcgVPAOGCCuSEJIUTZnTlT+LbT3r0QE5O/3tnZSA733ZefKNq2NcYubJEZU3VXAgOA5kqpWGCu1vpfSqmpwCaMGVbLtNbhVR2bEEJcj9YQHV08UZw7l7+Nmxv07AlTpuQnihYldsWwXWbMthpfyvINwIYqDkcIIUqlNZw8Cbt353/t2QNJScZ6e3tjPOLOO/OThJ8fNGxobtxVwVpuWwkhhOni4wsnit27IS7nqbLatY0ZTuPH5w9ke3tD3brmxmwWSR5CiBopKcm43VQwURzP6ZunFHh4GFcUPXoYX35+4OBgbszWRJKHEKLaS001Zj0VTBRRUfnrO3SAwEB45hkjUQQEQIMGpoV707TWxG3ZT/jbG3Bq4USvf0+r8HNI8hBCVCsZGRAeXjhRHDwImZnG+latjATxwAPGn927G1NmbZ3WmvhfDxK+cD2W/54nIbU+imz82p+7/s43QJKHEMJmZWfD4cOFE0VoKKSlGeubNDGSw8yZ+bef2rQxN+aKFr8tnPC31hH+W0JewmhXP5VeIxvi8dLdOAX6VMp5JXkIIWyC1nDiRPGZT5cuGesdHY2B7Keeyk8UHTva7nMU1xK/PYLwt9Zh2RZHfEp9QNPeKYWedzXAY9bd1L/Nr9JjkOQhhLBKcXHFZz7Fxxvratc2BrAnTMhPFB4extTZ6irhjyjCF/yIZes54q4YCaOdYwp3Dq+P56yR1O9XSeVzS1GNf9RCCFuRlGRcRRRMFCdOGOuUMp6lCA7OTxS+vjVj5tP5HYcJn7+W8F/PEZfsBGhc6qUwbJgjnjPvosHA7qbFJslDCFGl0tONwoA7d8KuXcafRWc+BQXBtGlGoujaFerXNy/eqnZ+11Es89cSvuU055KNb7xt3SsMG1IXjxdH0PD2HlZxL06ShxCi0mgNR4/mJ4ldu4wB7atXjfUtW0KvXtVv5lN5XQg5Rvj8tVh+Ps3Zy0ad9bZ1Uxh6R108Xwym4eBeVpEwCpLkIYSoMAkJRoIomCwuXDDWOTpCt24wdaqRMHr2BBcXq/udWGUSQ2MIf/MHLJtPceaSkTCcHVIYMqgOnjOG02hYkFX/cCR5CCFuSO6DdwVvPx07ZqxTCry8YMwYI0n06mW8rs4D2mWRuO84ljd/IHxTLGeSjITRpk4KQwbWxvOFO2k0/DarThgF1fB/SiFEWWRnG+MSBa8o9u3Lf/DO2dlIEk8+afzZrZttP6FdkS4eOEn4m99j2XiS0xeNhNG6TgqD+9vj+dxQGo/sZzMJoyBJHkKIYs6eLXxFsXt3/vMUDRoY4xMzZuTffmrd2tx4rU2S5RThb3yPZcNxTuUmjNop3NHXDs/nhtBk9ACbTBgFSfIQooZLTjYKBBZMFidPGuvs7IxpsRMm5N9+cnc3lovCkiJPY3njeyzrY4hNNBJGq9op3H5bLbymD6bJ3QNto0VgGUnyEKIGycoy6j4VvP108GB+W1RXV+jdO/+KIiDAGOgWJbt06CyWN74nfN0xYi8YCeNW+xRu763wfPYOmo69vVoljIIkeQhRTWkNsbGFryj27IErV4z1TZoYCWLUKCNZ9OhhTJ0V13bpyDnjCuPHo5w8bySMW+xTGRQIns/cTrNxg6ttwihIkocQ1cTly8bYxM6d+V9nzxrr6tQBf3+YNCn/9lOnTjZ/273KJEWdIWL+D1jWR3MywbgUu8U+lYE9NV7PDKLZ+CHWey8vOwtqVXxskjyEsEG5t59yk8SOHWCxGFcbYPTQvuOO/NtP0sio/BIPxBLx1losG2I4lZh7hZHCgB5ZeD09kOYPDLPehAFw5QTsfR4cXaDbexV+eEkeQtiA06cLJ4qQkPzbT02bGkni3nvzk0XTpubGa6suhB7HsmAtlk2xnEkyrjBa2acYt6SeHmjdVxi5stIg4l0If8N47TO3Uk4jyUMIK5OSYoxN5CaKnTuNsQswqsn6+8OjjxqJQm4/3byEXcewvPUjEVtOc/aSkTDa1E7mjj7g+cwdNLGlQe9T62HPs5B8FNqOha7vgpNLpZzKZpOHUsoJ2A7M1VqvMzseIW5E7sN3BRPFgQPGbSkwZj/16ZOfKAICoG5dc2OuDuL/dxjLwnVYfjlLXLKRMJwdkhnSvxYezw6h8aj+tpMwAJKPwZ7pcOpHaNgFBm6GVoMr9ZRVnjyUUsuAEUCc1tq7wPJhwAeAHfCp1nrBdQ41C1hdaYEKUQni4goPaO/ebZQjB2jY0LjlNHt2frKQ2U8VQ2tN3LZILO9uwLItjoQrjoDGxeEyw263x2P6UBoG97W9S7jMFLC8ZXzVqg0Bb0PnaWBXp9JPbcaVxwpgMfB57gKllB2wBBgMxAK7lVJrMRLJ/CL7TwJ8AQsgn8GE1UpLMyrIFkwW0dHGulq1jIfvxo3LTxRdutjWh11rp7Xm7M8Hsbz7E5b/nudCqqPRorXeJXoOqU2X54fTYIh1Fx8sldYQ+wPsnQ5XjkO78UbicKy6HrtVnjy01tuVUu2LLO4JHNFaHwNQSq0CRmmt52NcpRSilBoIOAGeQKpSaoPWOrtyIxeidLmlx3NvPe3caRQNzMgw1js7GwniqaeMP7t1AycnU0OulrTWnN6wD8vfNxHxRyKJafVQZOPqeImg4XXpMiOY+gOsox/GDbt0CPZMgzOboJE33L4Nbulf5WFYy5hHG+BkgdexQK/SNtZazwFQSj0CJJSUOJRSk4HJAC4ulTNgJGquCxfyH7zLfQjv/HljnZOT0ZfiuefyryraVN0HwhpHZ2ti1+7F8sFmIv5MIulqPWqRRYf6SfQZWY8uL47CsU/VtmitFBnJxgyqyHfBrh50fR86/8W4XWUCa0keJX0M0NfbSWu94hrrPgE+Aejevft1jyVEaTIyjM53Ba8qDh0y1uW2SM19Sjsw0Hhd00uPVzadrTn5zU7C//ELETuTuZxRl1pk0bFBEgOGO+I+ayz1evmaHWbF0BpOrIHQFyAlFlwfBv+3oN4tpoZlLW/xWKBtgdfOwGmTYhE1mNZGUcCCs5/27DHGLwBuucVIEg8/nF/So2FDc2OuKbKzsjmx6k8sS7YSEXKF5Iy62JFJp0aJeA5tS+dZY6nb1cvsMCvWxXDY8wyc2wpNAuC2r6FFb7OjAqwneewG3JRSrsApYBwwwdyQRE2QnGw8cFcwWZw5Y6xzcDD6Z+eOU/TqBe3a2fbtcluTnZlNzJe/Y/nwNyJDU7iSWRd7MnBrkojnne1wmzUOB193s8OseBmX4MBrELUIajeAHh9Cx8mVUmbkRpkxVXclMABorpSKxXhO419KqanAJowZVsu01uFVHZuo3rKzITIyP0ns2FG4omynTjBokHHrqVcvo6RHncqf8SiKyMrIInr5b1g+/i+R+66SmuVAbdLp3CwRj2BX3GZPpI5HJ7PDrBxaQ8y/IfRFSDsHHR8HvzehbnOzIyvGjNlW40tZvgHYUMXhiGosPr7wFcWuXfkNjRo1MhLEqFFGsujZE5pb3//PGiMzLYNj/9pKxNL/EXkgnbRsB+pwFfcWiXjc1ZFOsx+mtlt7s8OsXIn7IGQqxP8OzXpC/7XQrIfZUZXKWm5bCXFTrl7N76edmyxy+2nb2YGPj9HQKHdQu3NneabCbBkp6Rz9eAsRy3YQZcnianYdHLhKl1su4jHajY4zH8W+Qw2YKZmeCPtfhcMfQp2m0OtT6PAoKOt+g0ryEDZHa4iJKXz7KTQU0tON9W3aGEliyhR5psLapCdf5cg/NmL5PITDUZp0XZu6pONxaxKe93Shw6x7sWtbQ3ra6mw4tgLCZkP6eej0FPj+FRxso6qlJA9h9S5dMsp4FEwW8fHGunr1jGcqnn02f1Db2dnceEVhV5NSOfzBT1i+2Mvho4pMbY8jV/F2TsbzXi/avzgWu1bmTjutcudDIORpOL8LWtwG3RdDE3+zoyoXSR7CqmRlGX0pCiaKgn0qunSB4cPzB7W9vY1Ks8K6pJ2/QtR764hYuY8j0XZkYU99ruLfLgXPcT60e+FearVoZnaYVS8tAfbPgSNLoW5LCPwMXB+0ySl8kjyEqc6eLTxOsXu3MX0WjJ4UgYFw333Gnz16GK1ThXVKOXeJqHfWEbH6AEdP1CYbOxpyle4dr+I5MQDn6fdQq0ljs8M0R3YWHF0K++ZARhK4Pws+86BOI7Mju2GSPESVyS0UuGNHfrI4ftxYZ29v9Kl45JH8Qe2OHW3yA1mNciU2kciFP2L5TwTRp+ugqUVjlUavzsl4PtidNtPuRjVsYHaY5or/05hFlbgXWvY3blE19r7+flZOkoeoFFobs50KJoqChQJdXIwkkTtWERBgjF8I63f5WDwRC38k4odDHD/rgKYWTVUqvT1T8Hy0J62eGo2SGQqQFmcMhh9bDvVaQ++V0O7+avOJSJKHqBBJSfmFAnOTRUKCsc7Jybjl9Pzz+WMVrVqZG68on6TD54iY/wMR645yIr4eoGheK4W+vml4PhZEyydGoiT7G7IzjWm3+1+FrBTwmAner0Dt+mZHVqEkeYhyy8yE8PDCiSIiIn9Q29MT7rrLSBRSKNB2JR48RcRba7FsiOHUBaPb3i12KQzomoHn5D60eHSEPIJfVNx24xbVxQNw62DotggadTE7qkoh/6XFdZ05U3xQ+8oVY12zZkaCGD8+v1Bg48amhituwvnQ41gW/EjE5pOcuWgkjFb2yQzqmYXnXwbQbOIw+SRQkpTTEDbTKC3i6AJ9/wPOY6rNLaqSyLtAFJKWBnv35ieLHTvgxAljnb29MTbx6KP5VxUdOlTr/x81QvyfR7C8vZ6ILac5d9lIGG1qX2Zwb43HM4Nocu9g4zF9UVx2BkR9YBQxzE4Hr5fB6yWwdzQ7skonyaMGy+1+VzBR7NtXeFA7MBCmT5dB7epEa03c9igs72zAsvVcfj/vOpcY2l/h8ewQGo0aIPVbrufsLxDyDFyKgNbB0O19aFBNCzaWQJJHDZI7qF3wAbyC3e9kULv60lpz5udwLO/9RMT281xINdqztqubRI877PGYPpQGw/vKZWRZXDlpNGY6sQbqd4D+P0KbYt2yqz1JHtVU7qB2wUQREZG/3sMDRo7MTxReXnIru7rR2ZpT68KM9qz/S+RiXj/vi/S+szZdXgjGaVAvSRhllZkCEe+CZQGQDT5/Bc8Xwa6u2ZGZQn5dVBMFB7V37DAaHJU0qJ37pLYMaldPOltz8rsQLB9sIWLnJS6l57RnrX+RfiPr4v7iSBxv6yoJozy0huOrIGwWpJyEtvdAwDtQv73ZkZlKkocNKs+gdq9e8qR2dZedmc3xVX9i+XArkQXbszZMZNBdrXF/8R7q9vIzO0zblLAT9kyH8zuMNrC9v4SW/cyOyipI8rByBZ/Uzk0WRZ/UDgw0ntQODJRB7ZoiKyOLmM//i+Wj7USGppGS5ZDfnnWYC24z78PB38PsMG3XlZOw7yVj6m3dW6HXMnB9yKrawJpNkoeVSUrKLz9e9EltR0cZ1K7JMtMyOLZsGxFL/yDyQAZpWXWow1U6N7uAxwhXOs2cQB1PN7PDtG2ZV8DyNkQsNPpteP0feM42+oiLQiR5mKhg+fHcr4JPant4wIgR+c9UyKB2zZORks7Rj7ZgWb6TQ5bMvG577i0T8RzZiY6zHsa+U3uzw7R9OhtivjJqUaWeApf7wP+tGj+ucS3yq6gKnTtX+PZTSeXHx40zrih69pRB7ZoqPfkqh/+xkYjPQjh0SJOha1OPNDxaXcLzbnc6zHwcO5c2ZodZfcT/CXunG42ZmnaD21ZByz5mR2X1JHlUkqtX88uP5yaLmBhjnb09+PnBww/nX1XIoHbNdvViKof+vgHLV6EcKdBtz7dtMp5jPWk340nsWt9qdpjVy5UTxpXG8ZVQrxUErshpzCQPR5aFJI8KULCndm6yKNhTu21bI0FMnWr82bWrDGoLSE1IJurd9USs2s/RGDuysKM+aQS0T8VznC8uL4ylVvMa2G2vsmUkG2MaEW8br71eBs9Z1a7qbWWzyeShlKoFvA40BEK01p9V5fkvXy4+qB0XZ6yrV88Y1J4+PX9Qu3XrqoxOWLOUM0lEvv0jEd+Ec+xkbre9VHp0TMdjYgBtnxuLamy73eWsms6G6C+NWVSpp6HdePBfAE4uZkdmk6o8eSillgEjgDittXeB5cOADwA74FOt9YJrHGYU0Aa4AMRWYrhkZRmD2AWfqQgPL9xT+847828/eXvLoLYoLPnEBSIWriXi2yhizhjd9pqoNALdk/F8uAetp45BNZDZPJUq/g/jeY0LIdCsJ/RZAy16mx2VTTPj19wKYDHwee4CpZQdsAQYjJEMdiul1mIkkvlF9p8EuAN/aq0/Vkp9A/xSGYGGhUG/fsaVBhj9swMDYexY48+ePaWntijZpSNxRLy1loi1RzgeVxdQNFMp9PFKwWNSL26dMhrlWP0rr5ruynEInQUnvoZ6bSDoC2g/QcY1KkCVJw+t9XalVPsii3sCR7TWxwCUUquAUVrr+RhXKYUopWKBnBEFsko6j1JqMjAZwMXlxi5LO3aEBx6AoCDj9pObmwxqi9JdjDxDxIK1WNZHE5tgDGq1rHWF/v5peD5+Gy0euwtVt2bWQapyGclGDaqId4xE4T3XqENlL+1xK4q13GBpA5ws8DoW6HWN7b8F/qGU6gtsL2kDrfUnwCcA3bt31zcSVIMG8OGHN7KnqCkuhJ3AsvBHIjae4HSicSVxq91lBnZPx/PJfjR/aLh026tKOhuiP4ewlyDtLLSfCH7zwamt2ZFVO9aSPEr6PF/qL3ytdQrwWOWFI0Tp4nceI+LtdVh+PsW5S/nNk+4IysbjLwNoOm6oDHyZIW477HkOEvdCs0Do9z00v9ZnUHEzrOUdHgsU/GjgDJw2KRYhCtFaE/f7YSzvrCfi17PEJxsJo22dSwzpq/F45g4a3327dNszS3I0hM6Ek9+AY1vo/RW0Gyf3mCuZtSSP3YCbUsoVOAWMAyaYG5KoybTWnNliIeK9n7BsT+BCitELw8UhiTsH1aLLs0NpOKKfdNszU8YlCJ8Pke+Bsjf6a3i8UCNawFoDM6bqrgQGAM1zBr7naq3/pZSaCmzCmGG1TGsdXtWxiZpNZ2tOrQ/D8sHPRPxxIb95Ur2LBA21p8tzd1J/SG/5RGu27CyIXgH75kDaOaPard+b4CglW6qSGbOtxpeyfAOwoYrDETVcXvOkRVuI2JHfPKmDUyL97nLAfcYIHPt2l4RhLc5tg73PQWIYNO9ttIBt1sPsqGoka7ltJUSVyc7M5vjXO7As2UpkSHKB5kkXGDSiDe4v3k3dQH+zwxQFXT4KYTPh5Lfg6GIUL3S5T5K6iSR5iBohKyOLmC9+x/LP34gMSyMlM6d5UuMLOc2T7sUhwNPsMEVRGZfg4BsQ9T7Uqg2+f4Muz4O9FIczmyQPUW0ZzZN+y2melF64eVKwK51mjqeOV2ezwxQlyc6CY8tg/8uQFg8dHjESh6MUirMWkjxEtZKRks7Rj3/BsmxHgeZJadI8yZac/dUY17i4H1r0hQEbjD4bwqpI8hA2z2ietImIz3YXb540xh3XmY9j305m4li9S4ch7EWI/QGc2hvFC9veI+MaVkqSh7BJVy+mcuj9DVj+XaR5kvNlPO/1kuZJtiQ9CQ6+DocWQS0Ho5xIl+lgJ3XArJkkD2EzSm+elILnOD9cnr+HWi2amx2mKKvsDDj8MRx8Da6eh46TjHGNepL0bYEkD2HVUs5eInLhWiK+sXDspH2R5kn+tJ0+FtWksdlhivLQGk79aEy9vRQFtwyEgHehaYDZkYlykOQhrE5e86Tvoog5nds8KYVA9yw8H+pO62fuluZJturCHtg7A+K2QcMuxkN+rYNlXMMGSfIQVuHSkTgiFv5IxNrDHD9XpHnSo724dcoolJP0YrBZV04a5URivgCHFtDjQ+j4uPHshrBJkjyEaUptnuSXhucT0jypWsi4bDRlinzPuF3l+RJ4zoI60qfd1knyEFXKaJ60joiNxws3T+qWjucUaZ5UbWRnwtFP4cBcSIvLacr0JjjdWFdPYX0keYhKF78r2ui2V7B5kv1l7gjMwuPpgdI8qTrRGk7/ZDyvkWQxHvLrv06KF1ZD8j9WVDitNXH/PYTl3Q3SPKkmSdwHe1+Ac79AAzfo+x04j5LB8GpKkoeoEKU1T2rnkER3aZ5UvaWcNmpQHVsBDk2h2yLo9CTYye3H6kySh7hhpTdPSjSaJz0/nPqDg+STZ3WVkQwRb0PEO6AzjS5+XnOgTmOzIxNVQJKHKJfsrOyc5km/ELkzv3lSR2meVHNkZ8Gx5bD/FUg7Cy73g/98qO9qdmSiCknyENeVnZnN8VV/YvlwmzRPqunObIbQGXDxgNHJr9930DzQ7KiECSR5iBJlZWQR/fnvRHxUuHlS58YX8JDmSTXPxYNG0jizCep3kIq3QpKHyJeZlsGxf20j4tP/FWmedB6P4A7SPKkmSj0L+1+FY/8C+4ZGDarOT4Odg9mRCZNJ8qjhMlLSOfLRFiKW7yzUPKlLy0Q8Rnai48yHsHeTe9k1TmYKRLwLEW9Bdjp0ngberxizqYTARpOHUsoFWAwkAIe01gtMDsmmpF/OaZ70eeHmSZ6tkvC4uwsdZj6OnYs0T6qRdDZEf27UoUo9bdya8l8ADTqZHZmwMlWePJRSy4ARQJzW2rvA8mHAB4Ad8Ol1EkJnYL3W+mOl1OeVGnA1kZaYwqH3fyLiq/zmSU6k4ds2Gc+x0jxJYLR/DX0BEsOgWU+47Wto2cfsqISVMuPKYwXGVUPeL32llB2wBBgMxAK7lVJrMRLJ/CL7TwJCgTlKqfuBL6ogZpuUGp9M1LvriPj6QF7zpAak0bV9Ch7SPEnkSoqA0Jlweh04tYPeX0G7+0HJA52idFWePLTW25VS7Yss7gkc0VofA1BKrQJGaa3nY1ylFKKUmgHMzTnWN8DySg7bZlw5k0Tk2z8S8Y2F6JzmSY1ymid5TvTHWZoniVxpcXBgHhz5BOydwP8tcJ8m7V9FmVjLmEcb4GSB17FAr2tsvxGYp5SaAMSUtIFSajIwGcDFpXpX8rx8/ILRbe+7KGLOSPMkcR2ZqRD1PoTPh6wUcHsKvF+Fui3MjkzYkOsmD6XUn1rroAKvGwCdtNahFRhHSZPFdWkba60PAmOvdUCt9SfAJwDdu3cv9Vi2KulwfvOkE3E5zZNqpdDH6wqekwK55UlpniSK0NkQsxL2/R+knIA2IyFgITR0NzsyYYPKcuXhAKCUek9r/bzW+rJS6kMg6Dr7lUcs0LbAa2fgdAUev1q4GHEGy1triVgXTez53OZJydI8SVxf3Haj4u2FEGjSFYI+g1sGmB2VsGFlSR5KKdUSeEAp9YLWWgP1KjiO3YCbUsoVOAWMAyZU8Dls0vnQE0S8vQ7LxhOcSTR+7K3sLjOoezoeT0rzJHEdlw5B2CyI/R4cnSHoc6MxkwyGi5tUluTxEvA78BXwd6XUIeCG33lKqZXAAKC5UioWY+D7X0qpqcAmjBlWy7TW4Td6DlsXv/MYlrfXEfHzac5dMhJGm9qXuCMwE8+nB9JEmieJ60lLgIN/hcP/NAbA/d4A9+fAvqI/94ma6rq/gbTWGzGeq0ApFQTcCzx2oyfUWo8vZfkGYMONHteWaa2J+/0wlrfXY9l6joTkeoDGpU4SQ/tm4zFtMI3GDJLmSeL6stIgahGEvwmZl6HjE+DzGtS7xezIRDVTro+vWus/gT8rKZYaJbd5kuW9n4j4LYELqbnNky7SY5DCY/pQGgRL8yRRRjobYr4yngxPOQGth4P/QmjsZXZkopqSex9VKK950vs/E/G/ws2Teg8zmic53SHNk0Q5nf0VQl+ExL3GYHjgcrh1kNlRiWpOkkclM5on7cGyaAsROy9xWZoniYpyMdwYDD+9HhxdIOhLaD9eBsNFlZDkUQmu1TzJ8642dJ4hzZPETUg9A/vn5pRJbyBPhgtTSPKoILnNkywf/UZUTvOk2qTj1vgCHne64PaiNE8SNykj2egXHvlOTpn0Z3LKpDczOzJRA0nyuAm5zZMsS/9H1MHCzZM8R3Sg08wJ1PZ0MztMYeuyM+HYMuNqI+0suNwLfvOhQUezIxM1mCSPcspISefox79gWbajePOkUZ3o+KI0TxIVRGtjPCNsFiRZoMVt0jNcWA1JHmWQfvkqhxfnNE+KKtI8aYw7HWY+hl07Z7PDFNXJhT2wdwbEbYMGbtD3W3AeLRMrhNWQ5HEN5/8bzpZJK4s3T7rHi3YvSvMkUQmSY2D/yxDzb3BoDt0XQ6fJUKu22ZEJUYgkj2twqGfH6SMpdG2fhsf9vri88LQ0TxKVI/2i8VR41CLj6sLzJfCcBXUamR2ZECWS5HEN9bt3Yfql16QXhqg8Welw+EM4+DqkJ4LrQ+D7Oji1vf6+QphIksd1SOIQlUJrOLEG9r0Eycfg1sFGb40m/mZHJkSZSPIQoqrF/Q6hM+D8TmjsAwM2QuuhZkclRLlI8hCiqhTsrVGvNfRaZtymqiXVkoXtkeQhRGVLi4MDr8GRj8GuHvj+Dbo8B/aOZkcmxA2T5CFEZclMgaj3IXwBZKVApyfBZy7UbWl2ZELcNEkeQlS07CyI+QL2vQypp8B5FPgtgEZdzI5MiAojyUOIinRms9Fb4+J+aNYTblsJLfuaHZUQFU6ShxAVIXG/kTTObgYnV7htFbjcJ+VERLUlyUOIm5ESC/tfgWOfQZ3G0PU9cPsL2DmYHZkQlUqShxA3IuMSWN6CyL+DzgKPF8Dr/6BOE7MjE6JKSPIQojyyM+DIUjgwD67GQ7sJ4PcG1G9vdmRCVCmrTx5KqQ7AHKCR1npszjIn4EMgHdimtf63iSGKmkBriP3BeMjv8iFo2R8C3oFm3c2OTAhT1KrMgyulliml4pRSB4ssH6aUilJKHVFKzb7WMbTWx7TWjxVZfDfwjdb6CWBkBYctRGEJO2BLP/jvGFB20P9HuH2rJA5Ro1X2lccKYDHwee4CpZQdsAQYDMQCu5VSawE7YH6R/SdpreNKOK4zcCDn71kVHLMQhstHIOwlOPkN1L0FenwEHR+DWlZ/wS5EpavU/wVa6+1KqfZFFvcEjmitjwEopVYBo7TW84ERZTx0LEYCCaOUqyel1GRgMoCLi0u5Yxc1WFq8USL98D+NWVM+86DLC1C7vtmRCWE1KvW2VSnaACcLvI7NWVYipVQzpdRHQIBS6qWcxd8C9yil/gn8WNJ+WutPtNbdtdbdW7RoUUGhi2otMwXC58OPnYweGx0fg7uOGCVFJHEIUYgZ198lPTWlS9tYa30emFJk2RXg0QqOS9RU2VkQ/bnxvEZeOZH50MjD7MiEsFpmJI9YoGCbNGfgtAlxiJpOazizCcJmwsUDOeVEvoKW/cyOTAirZ0by2A24KaVcgVPAOGCCCXGImuzCXgidCed+gfodoc9qaDtWyokIUUaVmjyUUiuBAUBzpVQsMFdr/S+l1FRgE8YMq2Va6/DKjEOIPFeOG9VuY74Eh2bQbZFRKt2ujtmRCWFTKnu21fhSlm8ANlTmuYUoJD0Rwt+EqEWgaoHnS+A5C+o0MjsyIWySTFgX1VvWVTi0BML/BukXocPD4Ps6ODqbHZkQNk2Sh6iedDYcXwX75sCVGGg1DPzfgia+ZkcmRLUgyUNUP+e2Gr01LuyBJv7Q62e49Q6zoxKiWpHkIaqPiweNwoWnN4CjCwR9Ae0nGGMcQogKJclD2L6UU3BgLhxbDvYNwH8huD8DdnXNjkyIakuSh7BdGZfAshAi3zMaMrlPNxoyOTQzOzIhqj1JHsL2ZGfA4Y/h4GtwNQHajc9pyORqdmRC1BiSPITt0BpOfgv7XoLLh6HlAAh4W/pqCGECSR7CNsT/YcygSvgTGnlC/3XQeriUExHCJJI8hHW7FGU0ZIr9Duq1gl6fguvD0pBJCJPJHEZhnVLPwe6nYb0XnP3ZeCr8rsPSya8S2NnZ4e/vj7e3N3fddRcXL14EICYmBm9v77ztli5dSteuXUlMTGTNmjV4eXlRq1YtQkJCqjTexYsX06lTJ5RSJCQklLrdrFmz8Pb2xtvbm6+//rrY+meeeYb69aVPy42S5CGsS+YVOPC60ZDpyCfQaQqMPAreL4O9k9nRVUv16tUjLCyMgwcP0rRpU5YsWVJsmy+++IJ//OMfbN68mSZNmuDt7c23335Lv35VX77+tttuY8uWLbRr167UbdavX8/evXsJCwtj586dvP3221y6dClvfUhISF6SFDdGkoewDtmZcORT+NENDrwKrYZAcDj0WAx1W5odXY0RFBTEqVOnCi1bvXo1CxYsYPPmzTRv3hwADw8P3N3db/g88+bN48EHH2TQoEG4ubmxdOnSMu8bEBBA+/btr7mNxWKhf//+2Nvb4+TkhJ+fHxs3bgQgKyuLF198kYULF95w/ELGPITZtIbT640nw5Ms0DwI+nwDLXqbHVmVmz4dwsIq9pj+/vD++2XbNisri19++YXHHnssb9nx48eZOnUqoaGh3HrrrRUa2/79+9mxYwdXrlwhICCA4OBgGjRoQN++fUvc/quvvsLT07NMx/bz8+O1117j+eefJyUlha1bt+btu3jxYkaOHEmrVq0q7HupiSR5CPOcDzFmUMVtgwZu0Pc/4DxGZlBVsdTUVPz9/YmJiaFbt24MHjw4b12LFi1o2rQpq1ev5rnnnqvQ844aNYp69epRr149Bg4cyK5duxg9ejRhFZBBhwwZwu7du+nduzctWrQgKCgIe3t7Tp8+zZo1a9i2bdtNn6Omk+Qhql7yMaPa7fFV4NACui+BTk9ArdpmR2aqsl4hVLTcMY+kpCRGjBjBkiVLmDZtGgCOjo789NNP9OnTh5YtWzJx4sQyH3fOnDmsX78eoMSEoIp8SFBKcfny5Qq58sg9/5w5cwCYMGECbm5uhIaGcuTIETp16gRASkoKnTp14siRI2U+rsihta72X926ddPCCqTGax3yrNYra2u9ylHrsJe1Tk8yO6oaz8nJKe/ve/fu1W3bttXp6ek6Ojpae3l5aa21PnbsmHZxcdEbN24stG///v317t27y33OuXPnaj8/P52amqoTEhJ027Zt9alTp8p1jHbt2un4+PgS12VmZuqEhASttdb79u3TXl5eOiMjo9h2Bb93URwQokv5vSoD5qLyZaZC+AL4sSMc+ge4PmJMu/V7HWo3NDs6UUBAQAB+fn6sWrWq0HJXV1fWrl3LpEmT2LlzJ9999x3Ozs78+eefBAcHM3To0HKfq2fPngQHBxMYGMgrr7xC69aty7TfokWLcHZ2JjY2Fl9fXx5//HHAmEGV+/eMjAz69u2Lp6cnkydP5ssvv8TeXm60VCRlJJfqrXv37rqq56ILIDsLYr6A/a9ASiy0uQv8FxhPiIsabd68edSvX58ZM2aYHYq4BqXUHq11ifV/JBWLiqc1nNkEYTPh4gFo2gOCvoRb+psdmRCigkjyEBXrwl4InQnnfoH6HaHPamg7VmZQiULmzZtndgjiJtlE8lBKdQDmAI201mNzlo0GgoGWwBKt9WbzIhQkR8O+l+H4V+DQHLotgk5Pgl0dsyMTQlSCSh8wV0otU0rFKaUOFlk+TCkVpZQ6opSafa1jaK2Paa0fK7Lse631E8AjwP0VHrgom6sXYO8LsK4LxH5rNGO660hOJz9JHEJUV1Vx5bECWAx8nrtAKWUHLAEGA7HAbqXUWsAOmF9k/0la67hrHP/lnGOJqpSVBlH/gPA3IfOSMYPK9zVwdDY7MiFEFaj05KG13q6Ual9kcU/giNb6GIBSahUwSms9HxhRluMq4wmjBcBPWuu9FRiyuBadDdFfwv6XIeUktA42ZlA19r7+vkKIasOs5zzaACcLvI7NWVYipVQzpdRHQIBS6qWcxc8AdwBjlVJTSthnslIqRCkVEh8fX4Gh12BnNsNPXWHHw0axwtt/hQHrJHHYOKUUDz74YN7rzMxMWrRowYgRxue4tWvXsmDBghL3vZmS5sOHD79uZdsVK1Zw+vTpvNcDBgzAxcWFgo8YjB49+rpxXLx4kQ8//DDv9bZt2/K+vxtxs/uX17Bhw/Dz88PLy4spU6aQlZVV4nbz58+nU6dOuLu7s2nTJsB4ij44OJguXbrg5eXF7NnXHCUoM7OSR0lTb0p94ERrfV5rPUVr3THn6gSt9SKtdbec5R+VsM8nWuvuWuvuLVq0qMDQa6ALofDrENg6FDIuQe+VMHQX3DLQ7MhEBXBycuLgwYOkpqYC8PPPP9OmTf5nuZEjR1bYL5yCNmzYQOPGja+5TdHkAdC4cWP++OMPwEgKZ86cue65iiYPW7N69Wr27dvHwYMHiY+PZ82aNcW2sVgsrFq1ivDwcDZu3Mhf/vKXvCQzY8YMIiMjCQ0N5Y8//uCnn3666ZjMSh6xQNsCr52B06VsK8xy5Tj87yHY2A0u7IGuf4cREdB+HCgpTlCd3HnnnXl1qFauXMn48ePz1q1YsYKpU6cCEB0dTVBQED169OCVV17J22bbtm3069ePMWPG4OnpyZQpU8jOzs47no+PD97e3syaNStvn/bt25OQkEBMTAweHh488cQTeHl5MWTIEFJTU/nmm28ICQlh4sSJ+Pv75yW3cePG5T0B/+2333L33XcX+l7efvttevToga+vL3PnzgVg9uzZHD16FH9/f1588UUAkpOTGTt2LF26dGHixIl5VzO//PILAQEB+Pj4MGnSJK5evQrAxo0b6dKlC3369OHbb78t98/4kUceYcqUKfTt25fOnTuzbt26Mu/bsKFRiSEzM5P09PRidcEAfvjhB8aNG4eDgwOurq506tSJXbt24ejoyMCBxge9OnXq0LVrV2JjY8sdfzGl1S2pyC+gPXCwwGt74BjgCtQB9gFelXV+qW1VTlcvaL13htYrHbReVVfr0FlaX000O6rq79lnte7fv2K/nn32uqd1cnLS+/bt0/fcc49OTU3Vfn5+euvWrTo4OFhrrfXy5cv1008/rbXW+q677tKfffaZ1lrrxYsX59WG2rp1q3ZwcNBHjx7VmZmZ+o477tBr1qzRp06d0m3bttVxcXE6IyNDDxw4UH/33Xda6/zaVNHR0drOzk6HhoZqrbW+99579RdffKG1Ll47q3///nrHjh3ax8dHZ2Zm6sGDB+vo6Oi8ODZt2qSfeOIJnZ2drbOysnRwcLD+7bffCtXpyo23YcOG+uTJkzorK0sHBgbq//73vzo1NVU7OzvrqKgorbXWDz74oP773/+et/zQoUM6Oztb33vvvXk/n7J6+OGH9dChQ3VWVpY+dOiQbtOmjU5NTdWRkZHaz8+vxK/ExMS8/YcMGaIbN26sx48frzMzM4sd/+mnn877uWmt9aRJk/SaNWsKbZOYmKhdXV310aNHyxQzZta2UkqtBP4E3JVSsUqpx7TWmcBUYBMQAazWWodXdiziOrLSIOJdWNvR+LP9eBhxyBgQr9PY7OhEJfL19SUmJoaVK1cyfPjwUrf7448/8q5KCo6TgFGrqkOHDtjZ2TF+/Hh+//13du/ezYABA2jRogX29vZMnDiR7du3Fzuuq6sr/v7+AHTr1o2YmJhSY7Czs6NPnz58/fXXpKamFmoMtXnzZjZv3kxAQABdu3YlMjKSw4cPl3icnj174uzsTK1atfJK0kdFReHq6krnzp0BePjhh9m+fTuRkZG4urri5uaGUooHHnig1Piu5b777qNWrVq4ubnRoUMHIiMjcXd3JywsrMSvgrf1Nm3axJkzZ7h69Sq//vprsWPrEkpNFbxCyczMZPz48UybNo0OHTrcUPwFVcVsq/GlLN8AbKjs84sy0NkQsxL2zzFuVbUaBv5vQRNfsyOrWcyqyZ5j5MiRzJgxg23btnH+/PlStyvplklJy5VSJf5CK4mDg0Pe3+3s7PJuUZVm3LhxjBkzptiT6lprXnrpJZ588slCy0tKRkXPmZmZec14S/u+C3r00UcJDQ2ldevWbNhQ/NdbST+jqKgo7r+/5EfVtm3bViiB1K1bl5EjR/LDDz8U6rsC4OzszMmT+fOQYmNjCxWbnDx5Mm5ubkyfPv2630dZyI3rmu7sFtjYHf58AOo0hUE/w8CfJHHUQJMmTeLVV1/Fx8en1G1uu+22vPGGf//734XW7dq1i+joaLKzs/n666/p06cPvXr14rfffiMhIYGsrCxWrlxJ//5lr3HWoEEDLl++XGx53759eemllwqNzQAMHTqUZcuWkZycDMCpU6eIi4sr9ThFdenShZiYmLz+Hl988QX9+/enS5cuREdHc/ToUcAYxynJ8uXLCQsLKzFxAKxZs4bs7GyOHj3KsWPHcHd3v+6VR3Jyct6kgMzMTDZs2ECXLl2KHXvkyJGsWrWKq1evEh0dzeHDh+nZsycAL7/8MklJSbxfgR9QJHnUVIn7YOsw+HUwpF+A3v+GYSFw6x1mRyZM4uzszLPPPnvNbT744AOWLFlCjx49SEpKKrQuKCiI2bNn4+3tjaurK2PGjKFVq1bMnz+fgQMH4ufnR9euXRk1alSZY8odZC44YA7GJ/YZM2bk9VTPNWTIECZMmEBQUBA+Pj6MHTuWy5cv06xZM2677Ta8vb3zBsxLUrduXZYvX869996Lj48PtWrVYsqUKdStW5dPPvmE4OBg+vTpQ7t27cr8PRTk7u5O//79ufPOO/noo4+oW7fudfe5cuUKI0eOxNfXFz8/P1q2bMmUKcbTCWvXruXVV18FwMvLi/vuuw9PT0+GDRvGkiVLsLOzIzY2ljfeeAOLxULXrl3x9/fn008/vaH4C5KS7DXNlZNGifToz41xDK+XofNfwO76b2IhSrNt2zbeeeedcs0gqmkeeeQRRowYwdixY80OpcykJLuA9ItgWQCR7xuvPWaA10tQp4mZUQkhbJRceVR3WVfh8Idw8G+QngiuD4Lv6+DkYnZkQggrJ1ceNZHOhuNfw77/gysxcOsQCHgLmvibHZkQohqQ5FEdndsKoS8aT4U39oOBm6DVELOjEkJUI5I8qpOLByFsFpzeAI5tIehzaD9RSokIISqcJI/qICUW9r8Kx1ZA7UYQ8DZ0niozqIQQlUY+ktqy9CQI+z/40Q1i/g1dnoeRR42ZVJI4RBmdPXuWcePG0bFjRzw9PRk+fDiHDh2qkGM/8sgjfPPNN8WWDxgwgO7d88dhQ0JCGDBgQN7fp02bds3jlrcU/Lx583jnnXfKtY/ZPvroI3x8fPD396dPnz5YLJYSt9uzZw8+Pj506tSJadOm5T0lv337drp27Yq9vX2J/wY3S5KHLcpKh8gP4MeOYJkPbcfCiCjo+g44NDU7OmFDtNaMGTOGAQMGcPToUSwWC2+++Sbnzp0rtF1p/SNuRlxcXImlwbt3786iRYsq/Hy2ZsKECRw4cICwsDBmzpzJ888/X+J2Tz31FJ988gmHDx/m8OHDbNy4EQAXFxdWrFjBhAkTKiU+SR62JHcG1XoP2DsdGvvDsD3Q+wuo397k4IQt2rp1K7Vr1857YhnA39+fvn37sm3bNgYOHMiECRPySpaMHj2abt264eXlxSeffJK3T/369ZkzZw5+fn4EBgYWSz4Ar7zyCo888kheqfYXX3yRv/3tb8W2K9hoKTk5mUcffRQfHx98fX35z3/+k7ddSec7fvw4t99+O76+vtx+++2cOHGi2PEXLVqEp6cnvr6+jBs3DoALFy4wevRofH19CQwMZP/+/YBxxTJp0iQGDBhAhw4dyp3UrlWq/npyy7CD8ZR5SbW1zpw5w6VLlwgKCkIpxUMPPcT3338PGCXvfX19qVWrcn7Ny5iHrTi3LWcGVQg09oUBG40ZVGUo1iZsw8bpGzkbdrZCj3mr/60Me39YqesPHjxIt27dSl2/a9cuDh48iKurKwDLli2jadOmpKam0qNHD+655x6aNWvGlStXCAwM5I033mDmzJksXbqUl19+Oe84M2fOJCkpieXLl+f9EgwKCuK7775j69atNGjQoMTzv/766zRq1IgDBw4AkJiYCFDq+aZOncpDDz3Eww8/zLJly5g2bVreL9NcCxYsIDo6GgcHh7xOhnPnziUgIIDvv/+eX3/9lYceeoiwsDAAIiMj2bp1K5cvX8bd3Z2nnnqK2rVrl/5DL+FnaLFYaNeuHcOGDePbb79l7Nix3H///URFRRXb/vnnn+ehhx4CYMmSJbz33nukp6eXWEn31KlTODs75712dnbm1KlTZY7tZsiVh7W7eBC2BcMvAyHtHAR+BsP2QuuhkjhEpevZs2de4gDjU3vup/2TJ0/mlTuvU6dO3tVC0ZLqr7/+OhcvXuTjjz8u9un55ZdfLvHqI9eWLVt4+umn8143adLkmuf7888/827TPPjgg/z+++/Fjunr68vEiRP58ssvsbc3Pj///vvveSXmBw0axPnz5/NqdwUHB+Pg4EDz5s1p2bJliVdV11JSqXqAr7/+usRiiLmJA+Dpp5/m6NGjvPXWWyX+nK5Xhr0yyZWHtUqJhf1zIXoF2DcwSqR3fgbs65kdmagk17pCqCxeXl7XHEx1cnLK+/u2bdvYsmULf/75J46OjgwYMIC0tDQAateunfdLK7e8ea4ePXqwZ88eLly4QNOmhcfkBg0axCuvvMKOHTtKPL/WusRfhtc6X0El7bt+/Xq2b9/O2rVref311wkPD7/mL+GSSrcX9N133/Haa68B8OmnnxaaCFBSDLmvy3LlkWvcuHE89dRTxbZ1dnYu1BWwaBn2yiRXHtam0AyqL8F9ujGDynOmJA5R4QYNGsTVq1dZunRp3rLdu3fz22+/Fds2KSmJJk2a4OjoSGRkZKm/8IsaNmwYs2fPJjg4uMSy6HPmzGHhwoUl7jtkyBAWL16c9zr3tlVpevfuXahkfJ8+fQqtz87O5uTJkwwcOJCFCxdy8eJFkpOT6devX16J+W3bttG8efNCYw7XMmbMmLyrhqKJA0ouVQ/Xv/Io2MRq/fr1uLm5FTt2q1ataNCgATt27EBrzeeff16uqsU3Q5KHtchKh6hFBWZQ3ZMzg+pdcGhmdnSimlJK8d133/Hzzz/TsWNHvLy8mDdvXomfXocNG0ZmZia+vr688sorBAYGlvk89957L0888QQjR44s1uhp+PDhtGjRosT9Xn75ZRITE/H29sbPz4+tW7de8zyLFi1i+fLl+Pr68sUXX/DBBx8UWp+VlcUDDzyAj48PAQEBPPfcczRu3Jh58+YREhKCr68vs2fP5rPPPivz93Y9JZWqL4vFixfj5eWFv78/7733XqGYcrsuAvzzn//k8ccfp1OnTnTs2JE777wTMD4EODs7s2bNGp588km8vLwq7HsCKYxoPq3hxGqjBlXyMbjldghYCE27mh2ZEOIm2XqpeimMaK3ObYPQmXBht8ygEkLYFEkeZrgYnlODar1RgypwBbR/AGrZmR2ZEKICDRgwIO/J+epGkkdVSjkFB+bCseUyg0oIYdMkeVSF9CSIWAiRfwedZcyg8vo/GQgXQtgsq59tpZTqoJT6l1LqmyLLnZRSe5RSI8yK7boKzqAKfxPa3i0zqIQQ1UKlJg+l1DKlVJxS6mCR5cOUUlFKqSNKqdnXOobW+pjW+rESVs0CVldkvBVGazi+2qhBtedZoyHTsBDo/aXUoBJCVAuVfeWxAij02KxSyg5YAtwJeALjlVKeSikfpdS6Il8tSzqoUuoOwAKUr05AVTj3G2zqBX/cD/ZOMOAnGLQFmpZeP0gIM9nZ2eHv74+3tzd33XVXXr2nmJgYvL2987ZbunQpXbt2LfSg3jvvvINSioSEhCqL95FHHsHV1RV/f3/8/f3zalAVFBYWRlBQEF5eXvj6+vL111/nrZs4cSLu7u54e3szadIkMjIyqiz26qRSk4fWejtwocjinsCRnCuKdGAVMEprfUBrPaLIV1wphx4IBAITgCeUKt4qTyk1WSkVopQKiY+Pr8DvqhQXw2HbXfDLAEg7Y8ygGhYKrYfJ1Fth1erVq0dYWBgHDx6kadOmLFmypNg2X3zxBf/4xz/YvHlzXn2pkydP8vPPP+Pi4lLVIfP222/nPZFd8IG5XI6Ojnz++eeEh4ezceNGpk+fnpcUJ06cSGRkJAcOHCA1NZVPP/20aoOvJswY82gDnCzwOjZnWYmUUs2UUh8BAUqplwC01nO01tOBr4ClWutiNY611p9orbtrrbuX9vRqhUg5BTsfh598If6/4L8ARhyCDg/L1Fthc4KCgopVZV29ejULFixg8+bNNG/ePG/5c889x8KFC2+oEN+8efN48MEHGTRoEG5uboXKo1SEzp0755XzaN26NS1btiT3Q+Tw4cNRSqGUomfPnoVqQ4myM2O2VUnvtFIfc9danwemlLJuRQXFVH4Zl8CyECLfA50JnZ8F7zkyEC5u3J7pkBhWscds4g/d3i/TpllZWfzyyy889lj+EOPx48eZOnUqoaGh3HrrrXnL165dS5s2bfDz87vh0Pbv38+OHTu4cuUKAQEBBAcH06BBA/r27Vvi9l999RWenp6AUQ/rr3/9K7fffjsLFiwoVLywqF27dpGenk7Hjh0LLc/IyCixhIkoGzOSRyzQtsBrZ+C0CXHcmKx0OPIxHPwrXE2AdhPA729Q3/X6+wphhVJTU/H39ycmJoZu3boxePDgvHUtWrSgadOmrF69mueeew6AlJQU3njjDTZv3nxT5x01ahT16tWjXr16DBw4kF27djF69OgSxzAKmj9/Prfeeivp6elMnjyZt956i1dffbXEbc+cOcODDz7IZ599Vqwp0l/+8hf69etXarIS12ZG8tgNuCmlXIFTwDiMsQvrpjWc/AbCXoLko3DLoJwaVDIQLipIGa8QKlrumEdSUhIjRoxgyZIleT3EHR0d+emnn+jTpw8tW7Zk4sSJHD16lOjo6LyrjtjYWLp27cquXbsKXZ3MmTOH9evXA5SYEEoqVX758uXrXnm0atUKMEqlP/roo6X2Jr906RLBwcH87W9/K1bE8bXXXiM+Pp6PP/64DD8hUSKtdaV9ASuBM0AGxhXHYznLhwOHgKPAnMqMQWtNt27d9E05u03rjT21/jdar/fR+tRPWmdn39wxhbASTk5OeX/fu3evbtu2rU5PT9fR0dHay8tLa631sWPHtIuLi964cWOx/du1a6fj4+PLdc65c+dqPz8/nZqaqhMSEnTbtm31qVOnyrTv6dOntdZaZ2dn62effVbPmjWr2DZXr17VgwYN0n//+9+LrVu6dKkOCgrSKSkp5Yq5JgJCdCm/Vyt7ttV4rXUrrXVtrbWz1vpfOcs3aK07a607aq3fqMwYbkpyDPw20phBlXoaApfLDCpRrQUEBODn55fXEyOXq6sra9euZdKkSezcubNCztWzZ0+Cg4MJDAzklVdeKXMTo4kTJ+Lj44OPjw8JCQl57W5DQkJ4/PHHAWOQf/v27axYsaLYlN4pU6Zw7tw5goKC8Pf3569//WuFfD81jZRkv5aUU7CpB7g/C52nSQ0qISrIvHnzqF+/PjNmzDA7FHENUpL9Rjm2gZExYFfH7EiEEMKqSPK4HkkcQlS4efPmmR2CuElWXxhRCCGE9ZHkIYQQotwkeQghhCg3SR5CCCHKTZKHEEKIcpPkIYQQotwkeQghhCi3GvGEuVIqHjh+k4dpBCRVwX5l2f5mt7nWuuZA1bWFuzk3+m9i1jmq4j1U1m2vt528f6zvPGa8f9pprUtuiFRa0Sv5Klbk8ZOq2K8s29/sNtdZV2ohNGv7utF/E7POURXvobJue73t5P1jfeexpveP1pVcGLGa+bGK9ivL9je7zY1+L9amKr6PijxHVbyHyrrt9baT94/1ncea3j8147aVKDulVIgupRCaENcj75+aQ648RFGfmB2AsGny/qkh5MpDCCFEucmVhxBCiHKT5CGEEKLcJHkIIYQoN0keokyUUh5KqY+UUt8opZ4yOx5he5RSo5VSS5VSPyilhpgdj7g5kjxqAKXUMqVUnFLqYJHlw5RSUUqpI0qp2dc6htY6Qms9BbgPkKmYNUwFvYe+11o/ATwC3F+J4YoqILOtagClVD8gGfhca+2ds8wOOAQMBmKB3cB4wA6YX+QQk7TWcUqpkcBsYLHW+quqil+Yr6LeQzn7vQv8W2u9t4rCF5VAkkcNoZRqD6wr8B8/CJintR6a8/olAK110f/0JR1rvdY6uBLDFVboZt9DSikFLAB+1lpvqZKgRaWxNzsAYZo2wMkCr2OBXqVtrJQaANwNOAAbKjMwYTPK9R4CngHuABoppTpprT+qzOBE5ZLkUXOpEpaVehmqtd4GbKusYIRNKu97aBGwqPLCEVVJBsxrrligbYHXzsBpk2IRtkneQzWYJI+aazfgppRyVUrVAcYBa02OSdgWeQ/VYJI8agCl1ErgT8BdKRWrlHpMa50JTAU2ARHAaq11uJlxCusl7yFRlMy2EkIIUW5y5SGEEKLcJHkIIYQoN0keQgghyk2ShxBCiHKT5CGEEKLcJHkIIYQoN0keQgghyk2ShxBCiHKT5CFEFVNKjVFKaaVUlwLLnJVS0iBJ2AxJHkJUvfFACEYtqFy3A13NCUeI8pPyJEJUIaVUfeAoRve9NVprd6VUH+AH4CJwGRijtY42L0ohrk+uPISoWqOBLVrr/cAVpVRXrfXvGBVqR2mt/SVxCFsgyUOIqjUeWJ3z99U5rwHcgShTIhLiBkjyEKKKKKWaAT2BjTmLvgbuz1mepLXOMC04IcpJkocQVWcssEFrfRUg5/bUWcAT6cAnbIwMmAtRRZRS2wBf4FKBxc2AbwAvwBGYrLX+X9VHJ0T5SPIQQghRbnLbSgghRLlJ8hBCCFFukjyEEEKUmyQPIYQQ5SbJQwghRLlJ8hBCCFFukjyEEEKUmyQPIYQQ5fb/LfhFlYd+kKIAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "from pySDC.projects.Resilience.test_Runge_Kutta_sweeper import test_vdp, plot_all_stability\n", + "test_vdp()" + ] + }, + { + "cell_type": "markdown", + "id": "7fa01a8e", + "metadata": {}, + "source": [ + "In the above plot, you can see the local error versus the step size and p in the legend refers to the numerically computed order of the scheme.\n", + "Since we have a method of order 1, 2 and 4, we are happy to see the local error is of order 2, 3 and 5, just as we expect.\n", + "\n", + "Let's look at stability next." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "8a2ad4e9", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxAAAAFgCAYAAAArRJ8VAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAB1s0lEQVR4nO3dd3hUZf7+8fczk55AEiCU0DtIr4IgYEekiBXFtnZd1/Jbd9V1XVfXdf2u7q69gAUbWMEKCIKI2CiCdKSF3kN6nzm/P56AIC1lkpOZuV/XlWvI5Mw59yQhz3zmacZxHERERERERMrC43YAEREREREJHiogRERERESkzFRAiIiIiIhImamAEBERERGRMlMBISIiIiIiZaYCQkREREREykwFhIQ8Y8wKY8yQAJ+zhTHGMcZElH4+zRhztRtZRESk5ihtG9qU/vtFY8wDZXhMmdoQkZrCaB8IkfIzxrQANgKRjuOUVPAcfwfaOI5zRQCjiYiIi4wxDtDWcZx1FXz8NcD1juMMDGgwkQBSD4SIiIiIiJSZCggJecaYNGPMmcaYvxtj3jfGvGWMyTbGLDPGtDPG3GeM2W2M2WKMOfuQx80xxvzLGDPfGJNpjPnYGFPnGNeYY4y5/pDPbzDGrCq9zkpjTM/fZBkK/AW41BiTY4z5uaq/DyIicjhjTKox5kNjzB5jzEZjzO3GmDrGmK3GmBGlxyQYY9YZY64q/XxC6dCkmaV/4782xjQ/xvknGGMeOeTzUcaYJcaYLGPM+tK24GAbYozpCLwI9C9tGzKq/JsgUgEqICTcjADeBJKBxcAX2P8HjYGHgZd+c/xVwLVAKlACPH2iCxhjLgb+XvrY2sBIYN+hxziOMx14FHjXcZwEx3G6VfgZiYhIuRljPMCnwM/YNuAM4E6gD/bv/nhjTH3gf8ASx3HeOOThY4F/APWAJcDbZbheX+AN4E9AEjAISDv0GMdxVgE3A9+Xtg1JFXx6IlUqwu0AItXsG8dxvgAwxrwPXAA85jiOzxjzDjDOGJPkOE5G6fFvOo6zvPT4B4AlZZjodj3wb8dxFpR+XqFxsCIiUqX6ACmO4zxc+vkGY8x4YIzjOL8rbSNmAXWBLr957OeO48wFMMbcD2QaY5o6jrPlONe7DnjVcZyZpZ9vC9gzEalm6oGQcLPrkH/nA3sdx/Ed8jlAwiHHHNoYbAIise84HU9TYH1lQoqISJVrDqQaYzIOfGCHljYo/fo4oDPwmuM4+37z2INtg+M4OUA6tqf6eNQ2SMhQASFyfE0P+XczoBjYe4LHbAFal+HcWgJNRMQ9W4CNjuMkHfJRy3GcYcYYL3ZI6xvALQeWZT3EwbbBGJMA1AG2l+F6ahskJKiAEDm+K4wxJxlj4rBzJD44pMfiWF4G7jbG9DJWm2NMsNsFtCgdhysiItVrPpBljLnHGBNrjPEaYzobY/pgeyLAzoV4AnijtKg4YJgxZqAxJgo7F+LHEwxfAngF+J0x5gxjjMcY09gY0+Eox+0CmpSeW6RG0gsXkeN7E5gA7ARigNtP9ADHcd4H/glMBLKBj7DvTv3W+6W3+4wxPwUgq4iIlFHpm0EjgO7YfX32Yt8AOh34f8BVpcf8H7ZX4N5DHj4ReBA7dKkXdlL1ia43H/gddlJ2JvA1dhjVb80GVgA7jTEn6vEWcYU2khM5BmPMHOAtx3FedjuLiIjUDMaYCcBWx3H+6nYWEbeoB0JERERERMqs0gWEMSamdKOtn40xK4wxDwUimIiIBC+1DSIioavSQ5iMMQaIdxwnxxgTCcwD7nAc54dABBQRkeCjtkFEJHRVeiM5x1YgOaWfRpZ+aGKFiEgYU9sgIhK6ArITdenSZouANsBzjuP8eJRjbgRuBIg2kb0aRZ5oLy6RmslXVARASaMW7gYRqaQdOxbtdRwnparOX962IT4+vleHDkdb1VKk5luzZg0A7du3dzmJSMUtWlS2diGgqzAZY5KAKcAfHMdZfqzjWkanOg82ujFg1xWpTtmbNgOw78FXXU4iUjkPPWQWOY7Tu6qvU9a2oXfv3s7ChQurOo5IlRgyZAgAc+bMcTWHSGUYU7Z2ISA9EAc4jpNRuvTlUOCYjYSIiIQPtQ1So5TkQ3EG+PLB7wOPF7xxEJkIEbFupxMJCpUuIIwxKUBxaQMRC5yJ3XRFRETClNoGcZXjh8xVkL4A9i+FrNWwbz6UZIG/+NiP80TZQiLlFKh9EtTpCfVOgbjU6ssuEgQC0QPRCHi9dKyrB3jPcZzPAnBeEREJXmobpHrlboHtn8HqpyB3E/gL7P0mAqLrQUx9iGwDEQngiQZPJBgPOA74i8BXAL4cKMqCPd/B1k8Bvz1HVB1oMRaajoaUQbbXQiSMBWIVpqVAjwBkERGREKG2QapFwW5ImwirnoD8bfa+yCRI7AhxzSC2MUTXtYVCeflLoGAX5G2GnI2w9gX45RlbgLT7A7S5ERJaBPLZiASNgM6BEBEREalSjgN7voVfnoUt79vhSjENocEZULsDRNUFYyp/HU8ExDW2H/X6216K7LWQsRRWPmY/Ek+CAe9AUufKX08kiKiAEBERkZrPcWDbp7DgVtvb4I2BOn0huYcdnlTVPFGQ2Ml+FGVC+nxIXwhTu0DLq6gTW0h6fnTV5xCpAVRAiIiISM3lOLBjBvzwOyjYAZHJ0GgYJHe38xjcEJUIDc+CegNh7zxIe5s3L/Hz8sLWtkekIkOmRIKICggRERGpmTJXwtcjIWe9ndvQeBQkda05L9AjYm0hkdyLFTPHc8eAdfBJGzj7W4ht5HY6kSpTQ/4Hiogc7uGHvbz4Yneef74zkyaNoKAgA4CMjDSef/7X8caLFo3npZd6kp+/nxUr3uf55zvx0EMetm/XhmQiQaskFxb/CT7vYocrNTwH2t5mex1qSvFwqOg6/HlKfR6fWQfytsKnbWH3PLdThSSv10v37t3p3LkzI0aMICMjA4C0tDQ6d/61bRg/fjw9e/Zk//79vP/++3Tq1AmPx4M2qwyMGvi/UEQEIiJiufnmJdx663JiY+swf/5zRxzz889vMn/+M1x55QxiY5OpX78zl1wymebNB7mQWEQCYudsmNLErqyU3A3a/gHq9QuCpVMNny+vBa2vt8vEzhoMaZPcDhVyYmNjWbJkCcuXL6dOnTo899yRbcObb77JM888w4wZM0hOTqZz585MnjyZQYPUNgSKhjCJSI3XpEl/du1aeth9K1a8x7ffPsZVV80iLq4eACkpHd2IJyKBUJIPS+6xS6VG1YGW10B8c7dTlV9MfVtEbHoXvrscSrLtkq8ScP3792fp0sPbhvfee4/HHnuMWbNmUa+ebRs6dlTbEGgqIETkuKZPh507A3vOhg1h6NCyHev3+9i4cRY9elx38L7MzE1MnXobN920mISEhoENJyLVL2M5zDodCvdA3ZPtkqxuTZAOBG8stLgCNr8H82+yPRKtrnY7VUDdeScsWRLYc3bvDk8+WbZjfT4fs2bN4rrrfm0bNm3axG233cbixYtp2FBtQ1XSECYRqZFKSvJ58cXu/PvfdcnPT6dVq7MOfi0uLoXExGasWPGeiwlFJCDWvwbTeoAvD5pfAY2GBnfxcIAnAppdAvGt7ApS279wO1FIyM/Pp3v37tStW5f09HTOOuvXtiElJYVmzZrx3ntqG6qaeiBE5LjK2lMQaAfmQBQUZDJp0nAWLHiOk0++HYDIyDjGjp3Ga68NJD6+Pl27jnUnpIhUnK8AFv4B1r8M8S2hyQUQmeB2qsA6UERsfA3mjoJhy6B2W7dTBURZewoC7cAciMzMTIYPH85zzz3H7bfbtiEuLo5p06YxcOBA6tevz9ixahuqinogRKRGi4lJZOjQp/nuuyfw+YoP3h8fn8LYsdOZPfsvrFund/ZEgkr+Drvc6fqX7V4KLa4IveLhAG80NBtjV4/6chD4Ct1OFBISExN5+umneeKJJygu/rVtSElJYfr06fzlL3/hiy/UNlQVFRAiUuM1atSDhg27sXz5O4fdn5zckjFjPuGTT65l69YfWbVqCv/9bxO2bv2eiRPP4623znEpsYgcU/oi+LQ9FO6GphdDwzNq5tKsgRSVZPewKNgJyx92O03I6NGjB926deOddw5vG1q2bMknn3zCtddey48//siUKVNo0qQJ33//Peeddx7nnKO2obKM4zjVftGW0anOg420IoEEp+xNmwHY9+CrLicRqZyHHjKLHMfp7XaOA3r37u1ojfYQt+UjmHcxRMRD88sgpoHbiQJmyHUTAJjzyjXHPmjrx5DxMwz7GZK6VEsukfIwpmztQoiX/CIiIlIjrH4Svhlti4bW14dU8VBmDc8Cbwx8PQpceANXJFBUQIiIiEjV8ftg0Z3w011QuyO0vBoiQnS+w4lExEH9wZC7EXZ+6XYakQpTASEiIiJVw1cAU7vAmqfs/g5NLw6NJVorI7kXRNSGH29wO4lIhamAEBERkcAr2m8nS2etgoZn2/0djHE7lfs8EVCvP+RtshPKRYKQCggREREJrNwttnjI3wZNL7QvmOVXyd3BRMK6cW4nEakQFRAiIiISOBnL4PNOUJwFzcdCYme3E9U83hhI7AgbXte+EBKUVECISI300EOGKVOuPPi531/C44+nMHHicADWrPmEefMeO+pjH3204hM03357GAUFGcc9ZsmSCWRnbz/4+YQJQ/jf/5px6LLY77xz/glzFBRksGDB8wc/T0ubc/D5VURlHy9SaTtnw/Q+9t+tfgcJLd3NU5MldgJ/Iez6yu0kQcUYw5VX/to2lJSUkJKSwvDh9m/fJ598wmOPHb1tSEioeNswbNgwMjIyjnvMhAkT2L7917ZhyJAhNGt2eNtw/vnnnzBHRkYGzz//a9swZ86cg8+vIir7+KNRASEiNVJkZDy7dy+nuDgfgPXrZ1KrVuODX2/ffiQDB94b8OuOHTuVmJik4x7z2wICICYmiS1bvgVsYZCTs+OE1/ptASES1Da+BV+dBZG1odV14blMa3nEt7LDmLZ/7naSoBIfH8/y5cvJz7dtw8yZM2nc+Ne2YeTIkdx7b+DbhqlTp5KUlHTcY35bQAAkJSXx7be2bcjIyGDHjhO3Db8tIGoiFRAiUmO1aXMua9faxnX58kl07nzZwa8tWTKBqVNvA2D//o288kp/xo/vw+zZDxw8Ji1tDq+9Noh33x3Nc8+dxGef3Yzj+AFYtmwSL7zQheef78zMmfccfMyTT7YgL28vGRlpPPdcRz755Aaef74Tb755NsXF+axc+QHbty9k8uSxvPhi94MFTufOYw7ulL1q1WQ6dLjgsOfy7bePM358H154oStfffUgAF9+eS/796/nxRe7M2PGnwAoKsrhvfcu4tlnOzB58tiD71xt2DCLl17qwQsvdOHjj6+lpMQOe1i3bjrPPtuBV18dyKpVkwP0nRcpB8eBZQ/D91dCXDNodS1EJbqdqubzREB8M9j0rttJgs65557L55/btmHSpElcdtmvbcOECRO47TbbNmzcuJH+/fvTp08fHnjg17Zhzpw5DBo0iNGjR3PSSSdx88034/f7D56vS5cudO7cmXvu+bVtaNGiBXv37iUtLY2OHTtyww030KlTJ84++2zy8/P54IMPWLhwIWPHjqV79+4HC5wxY8Yc3Cl78uTJXHDB4W3D448/Tp8+fejatSsPPmjbhnvvvZf169fTvXt3/vQn2zbk5ORw0UUX0aFDB8aO/bVtmDVrFj169KBLly5ce+21FBbatmH69Ol06NCBgQMHMnly4NuGiICfUURCytDpd9Jw55KAnnNnw+5MH/rkCY/r3HkMX3/9MO3aDWfXrqX06HEtmzd/c8Rx06ffQe/et9Ct21XMn//cYV/btm0+v//9ShITm/P220NZtWoyTZuewpdf3sONNy4iNjaZN988m9WrP6JDh/MPe+y+fWu58MJJjBw5nvffv4RVqz6ka9crmD//Wc4++wlSU3/drLNlyzP49NMb8Pt9LF/+DiNGjGPu3H8AsH79DNLT13L99fMBh0mTRrJp01zOPPMxdu9ezs03LwFswbNz52JuvXUFtWql8uqrA9iy5VtSU3vz8cfXcNVVs6hbtx1TplzFwoUv0Lv3zXz66Q1cddVs6tRpwwcfXFqun4NIpfkKYHpfyFwGSd0gdQR4vG6nCh5xTWH3HDtfJLK222nK5847YcmSwJ6ze3d48skTHjZmzBgefvhhhg8fztKlS7n22mv55psj24Y77riDW265hauuuornnju8bZg/fz4rV66kefPmDB06lMmTJ3PKKadwzz33sGjRIpKTkzn77LP56KOPOP/88w977Nq1a5k0aRLjx4/nkksu4cMPP+SKK67g2Wef5YknnqB371/bhjPOOIMbbrgBn8/HO++8w7hx4/jHP2zbMGPGDNauXcv8+fNxHIeRI0cyd+5cHnvsMZYvX86S0u/vnDlzWLx4MStWrCA1NZUBAwbw7bff0rt3b6655hpmzZpFu3btuOqqq3jhhRe4+eabueGGG5g9ezZt2rTh0ksD3zaoB0JEaqwGDbqSkZHGsmWTaNt22DGP27Ll24O9E926XXnY1xo37ktycis8Hi+dO1/G5s3z2LZtAS1aDCE+PgWPJ4IuXcayadPcI86bnNyShg27A9CoUS8yMtKOmcHj8dKs2UBWrHiXkpJ8kpJaHPza+vUzWL9+Bi+91IOXXurJ3r2r2bdv7VHP07hxX2rXboIxHho06E5GRhp7964hKakldeu2K32OV7Np01z27l1den9bjDF07XrFMfOJBFz+TvikrS0e6p8GjUepeCiv2FR7u/9nd3MEma5du5KWlsakSZMYNuzYbcO33357sHfi0HkTAH379qVVq1Z4vV4uu+wy5s2bx4IFCxgyZAgpKSlEREQwduxY5s49sm1o2bIl3bt3B6BXr16kpaUdM4PX62XgwIG8++675Ofn06JFi4NfmzFjBjNmzKBHjx707NmT1atXs3bt0duGvn370qRJEzweD927dyctLY01a9bQsmVL2rWzbcPVV1/N3LlzWb16NS1btqRtW9s2XHFF4NsG9UCIyHGVpaegKrVvP5KZM+/m6qvnkJ+/75jHmWOsL3/k/QZwjnboEbze6IP/9ni8lJTkH/f4zp3H8O67oxk8+O+H3e84DgMH3kfv3jcddv/RCpLfXtPvLzlu3mM9b5EqtfdHmH2W7YFoejEknuR2ouAUnWJvs1ZB/VPdzVJeZegpqEojR47k7rvvZs6cOezbV/m2wRhz2GTn44mO/vXvtNfrPThc6VjGjBnD6NGj+fvf/37Y/Y7jcN9993HTTYe3DUcrSH57zZKSkuPmreq2QT0QIlKj9ehxLYMG/Y0GDboc85imTQccnH+wdOnbh31t27b57N+/Ecfxs2LFuzRrNpDGjU9m06avycvbWzrkaBLNmw8uc6bo6FoUFmYfcX+zZqcycOB9dOly2WH3t2lzDkuWvEpRUQ4AWVnbyM3dTVRULYqKjjzPb9Wr14GMjDTS09eVPsc3ad58MPXqdWD//o2kp68H7DwRkSrlOLD2JZg5AIzXzndQ8VBxkbXBeCBno9tJgs61117L3/72N7p0OXbbMGDAgIPzD95++/C2Yf78+WzcuBG/38+7777LwIEDOfnkk/n666/Zu3cvPp+PSZMmMXhw2duGWrVqkZ195N/0U089lfvuu++wuRoA55xzDq+++io5ObZt2LZtG7t37z7meX6rQ4cOpKWlsW6dbRvefPNNBg8eTIcOHdi4cSPr19u2YdKkwLcN6oEQkRqtdu0m9Ot3x3GPGTr0KSZPvpwff3yKjh0vPOxrTZr0Z9ase9m1axnNmw+iY8fRGOPhjDP+xeuvn4bjOLRtO4wOHUaVOVO3btfw+ec3ExERy3XXfX/wfmMMp5xy9xHHt259Nnv2rOKVV+xmWlFRCYwe/RZ16rSmadMBPP98Z9q0OZd27c476vUiImIYNeo13n//Yvz+ElJT+9C7981EREQzYsQ4Jk48j7i4ejRrNpDdu5eX+XmIlEtxDiy4BdLegoTW0ORCiIh1O1VwMx6ISID87Sc+Vg7TpEkT7rjj+G3DU089xeWXX85TTz3FhRce3jb079+fe++9l2XLlh2cUO3xePjXv/7FaafZtmHYsGGMGlX2tuGaa67h5ptvJjY2lu+/P7xtuPvuI9uGs88+m1WrVtG/v20bEhISeOutt2jdujUDBgygc+fOnHvuuZx33tHbhpiYGF577TUuvvhiSkpK6NOnDzfffDPR0dGMGzeO8847j3r16jFw4ECWLw9s22DK2l0TSC2jU50HG91Y7dcVCYTsTZsB2Pfgqy4nkRNJS5vDd989weWXf+Z2lBrpoYfMIsdxep/4yOrRu3dvZ+HChW7HkKNJX2yHLBWlQ/3BkDIINHzuMEOumwDAnFeuKd8D170EdXrDkE8DnkmObs6cOTzxxBN89pnaht8ypmztgoYwiYiIyNH5fbDycfiiN/iLoMVVtoBQ8RA4nigoyXE7hUi5aAiTiISsFi2G0KLFELdjiASn7PUw63TI2wy1O0LqcIiIcztV6DFeyFnvdoqwMmTIEIYMGeJ2jKCmAkJERER+5ffBmqdgyT32xW3jUXaPB/U6iEgpFRAiIiJipS+COSOgYAfUagep5wXfBmfBximBWm3dTiFSLiogREREwl1hOix9ANa+YIcpNb0QandSr0N18BVCZC23U4iUiwoIERGRcOUvgfXj4ac/2k3h6vSBBqeBN8btZOHDlw9Rdd1OIVIuWoVJRGqsnJydfPDBGJ5+ujXPPXcSb789jH37fqn0eT/66BpWrvzgiPsnTBjCuHG/rl63fftCJkwYcvDf06bdftzzPvpoQrlyzJnzd7777olyPUYkIBwHtn0OU1Jhwa0Q0wDa3ASp56p4qE6O367AFNvQ7SRBw+v10r17dzp37syIESPIyMgA7O7NnTt3Pnjc+PHj6dmzJ/v37z943xNPPIExhr1791Z37JCjAkJEaiTHcXj33dG0aDGE229fz+9/v5IzzniUnJxdB4/x+30Bv25u7m7Wrp12xP2pqb0599ynA349kWq390f4pBV8PRwcHzS7xC7PGtPA7WThpzgbcCC+udtJgkZsbCxLlixh+fLl1KlTh+eee+6IY958802eeeYZZsyYQXJyMgBbtmxh5syZNGvWrLojhyQNYRKRGikt7Ss8nkh697754H0NG3YnLW0Or79+GgkJjdi5cwm///1K3nnnfLKytlBSUsDJJ99Br152o8pHH03g5JPvYO3az4iIiGXMmI9JSDj8RdLs2Q+QlbWFUaPsxoCnnPInvvnmEdq2Pfc3eX7dlK6oKIdp0/7A9u0LAcPgwQ9y0kl2l9NZs+4/4noZGZv45JNryc3dQ3x8CqNGvUZi4uGN2I8/Ps3ChS/i8USQknISF130Dvn56Xz88bXs37+ByMg4RowYR4MGXZkz5+9kZm5m//4NZGZupl+/Ozn55OP3joiwfynMuwSy14A3HhqdC8m9wON1O1n4Kkq3twmt3c0RpPr378/SpUsPu++9997jscceY9asWdSrV+/g/XfddRf//ve/y7WztBybCggROb7p02HnzsCes2FDGDr0uIfs3r2cRo16HfVr27bN55ZblpOc3BKAUaNeJTa2DsXF+Ywf34eOHS8kLq4uxcW5NGnSjzPO+CczZ/6Zn34az6BBfz14npkz/0xBQSajRr2GKZ0s2rRpf1avnsLGjV8RHX30iY1ff/0PoqMTueWWZQDk59su8mNdb9q02+ja9Sq6d7+axYtfZdq02xkz5qPDzjlv3mPcccdGIiKiKSjIAOCrrx6kYcMejBnzERs3zmbKlKu4+eYlAOzdu5qrr/6KoqJsnn22Pb1734LXG3nc76mEqYwVtnDIWgmeaKh/GtTtB94ot5NJ4R57W7uDuzkqYtGdsH9JYM+Z3B16PVmmQ30+H7NmzeK66647eN+mTZu47bbbWLx4MQ0b/jos7JNPPqFx48Z069YtsHnDmIYwiUjQady478HiAey79y++2I1XXulHVtYW0tPXAuD1RtGu3XAAGjXqRUZG2sHHzJ37DwoKMhgx4qWDxcMBgwb9lW++eeSY19+48Uv69Pn9wc9jY5OPe70tW76nS5fLAeja9Uo2b553xDkbNOjK5MljWbr0LTyeiNLHzaNbtysBaNnydPLz91FQkAlA27bnERERTVxcPeLj65Obu+uIc0qYy1gGn3WCqZ0hZx2knArt74D6g1Q81BQFu+yck9hUt5MEjfz8fLp3707dunVJT0/nrLPOOvi1lJQUmjVrxnvvvXfwvry8PP75z3/y8MMPuxE3ZKkHQkSO7wQ9BVUlJaXTUSc6A0RGxh/8d1raHDZs+JLrrvueyMg4JkwYQklJAQAeT+TB4sDj8eL3lxx8XGpqH3bsWER+fjqxsXUOO3/Llqfz1VcPsHXrD0e9vuM4RxQdJ7reoY722Msv/5xNm+ayZs0nzJ37D269dQWO4xzzsRER0Yfcd+xrSRjatxC+uxKyV4MnClIGQt3+2kW6JsrfATGNgnO53DL2FATagTkQmZmZDB8+nOeee47bb7dDOOPi4pg2bRoDBw6kfv36jB07lvXr17Nx48aDvQ9bt26lZ8+ezJ8//7BeCikf9UCISI3UsuXp+HyFLFo0/uB927YtIC3t68OOKyjIJDY2mcjIOPbuXX3MF/2/1abNUAYMuJeJE8+jsDD7iK+feur9fPvtv4/62Natz2b+/GcPfn5gCNOxNG16CsuXvwPAsmVv06zZwMO+7jh+srK20LLlaZx11r8pKMigqCiH5s0HsXTp24AtlOLi6hEdrU295Bh2fwOftIEv+kBuGqQMgnZ3QoMzVDzURP5i2wPR/FK3kwSlxMREnn76aZ544gmKi4sP3p+SksL06dP5y1/+whdffEGXLl3YvXs3aWlppKWl0aRJE3766ScVD5WkHggRqZGMMVx66RSmT7+Tb799jIiIGJKSWtC+/fmHHdemzVAWLXqRF17oSr167WnSpF+Zr9Gp08UUFWXzzjsjufzyqYd9rW3bYcTHpxz1cYMG/ZWpU3/P8893xuPxMnjwg3TseMExr3PuuU/z8cfX8t13jx+cRH0ov9/H5MlXUFiYieM49Ot3FzExSQwZ8nc+/vh3vPBCVyIj4zj//NfL/NwkTDgO7JgO82+GvM3gjbMFQ50+4I0+8ePFPfnbAD/UO8XtJEGrR48edOvWjXfeeYdTTz314P0tW7bkk08+YdiwYUyePJmTTz7ZxZShyRyti7xcJzCmKfAG0BDwA+Mcx3nqeI9pGZ3qPNjoxkpdV8Qt2Zs2A7DvwVddTiJSOQ89ZBY5jtP7xEeWX0Xaht69ezsLFy6sijihx++DrZNh4e1QsBMia9sXosk9waPJ9G4Yct0EAOa8ck3ZHrD7a9g9By7cB9F1Tni4SHUwpmztQiB6IEqAPzqO85MxphawyBgz03GclQE4t4iIBCe1DVXBXwwb34LFf4KifXYH48YjIbGrlmMNNjkbIaahigcJSpUuIBzH2QHsKP13tjFmFdAYUCMhIhKm1DYEmK8A1r8KP/8FijPtC8+mF0HtjmA0nTHo+Aohbwt0vNvtJCIVEtA5EMaYFkAP4MejfO1G4EaAut7EQF5WRERqsLK2Ddoh9ihK8mDdS7D0b1CSA3FNIfU8SGgTnCv3iJW7EfBDqjur3IlUVsAKCGNMAvAhcKfjOFm//brjOOOAcWDnQATquiIiUnOVp23o3bu32oYDSnJh7Yuw9EHw5UJ8C2hygb1V4RD8sn+xm/qlDDzxsSI1UEAKCGNMJLaBeNtxnMmBOKeIiAQ3tQ0VUJIP616Enx8oLRxa2Y3f4pu7nUwCxfHbAqLJKE14l6BV6QLC2F2NXgFWOY7z38pHEhGRYKe2oZx8RbDhFVh8D5RkQ3xLqD8E4jWsK+TkbbU9TE3OdzuJSIUFogdiAHAlsMwYs6T0vr84jjP12A8REZEQp7ahLBw/bHoXFvweivfbOQ5NS4cqSWjKWgXGC43PczuJSIUFYhWmeYAGZIqIyEFqG8pg11fw7Vgo2AExDSD1ck2ODnWOYwuIhNZ27w6RIKWdqEVERKpT1i/w9Qg7Dj4yERqfD0ldVTiEg/xtdhneXk+7nUSkUlRAiIiIVIfiLFj+D1j9XzAR0OAMqHuyJtKGk8zldvhSk5FuJxGpFBUQIiIiVclxYNMkmH+T3cshqbstHiIT3E4m1clxIHOlHaYWleR2GpFKUQEhIiJSVbLWwlfn2I3DYlOh2RiIa+x2KnFD3ma7wlaXB91OIlJpKiBEREQCzV8Mq56ApQ/Y4UqNhkGdXmA8bicTt2Qut78LjUe4nUSk0lRAiIiIBFLGMph9jl1dqXZHaHQuRNZyO5W4yfHb4UtNL9DQNQkJKiBEREQCwe+DVY/D0vvBEwtNL4bEk9xOJTVB7kbw5UHzS91OIhIQKiBEREQqK3czzDzVjnOvfRKkngcRcW6nqlrpubBpH2zdD7uyYE82ZORDVj7s3APFfvA7gAGvgUgPNG4AibFQJw5SakGjRGhaB5rVhdgQXo0qYzl4oiB1mNtJRAJCBYSIiEhlbJlsN4TDH7p7OqTnwpItsHQrfLsCtudCTvHhx8R6ISEK4iIg2gsxEeAp/T74/FDkhw3bIbfYPrbY/+tjDVA3Bnq0gk6p0K0JdGoM0SHwMsXxQfYaaHYxeGPcTiMSECHwP1NERMQFviJY8mdY85RdYanJhRBdx+1UgVHih582wTdrYeZS2JFr7/caaBQPneva2/qxUC8WkqMhylv28zsO5JVAegHsyYddefYaP22AmSvtMREGWiXC8D4wqB20rBf451kdcjaCL98WECIhQgWEiIhIeeVthxn9IG8L1O0LDc4GTzleQNdEfgcWpMH0ZfYjr8QWDK0ToVdLaJMITWvZoUiVZQzER9qPpr+ZYJ5TDBsyYV0GrNkP/51pPxrEwQV9YFhXaBZEhVrWKjt8qdE5bicRCRgVECIiIuWx90eYfSb4C6HphZDY2e1ElbMnGyb/BO98D+mFdvhRl7rQLQU6JNvPq1NCJHStZz/A9lIs3wdL9sCLX8MLX9ui5trT4eyTIKoGv5Rx/Hb4UpPzNXxJQkoN/l8nIiJSw6RNgu+vhIha0OoKiGngdqKK+2UXTPgOpi2zvQ/tk2BEK1s8lGc4UlWrEwODGtuP/YWwcBf8sBPunwKPfQZXDoAxfe3k7JomfxuU5NoCQiSEqIAQERE5EceBFY/C0r9CXHNodknwrrK0agc88qF9Vz/aC6em2o/6QfB8kqPhrGZwRlM7vOnrbfD8HHh5LlzRH64+BZJq0PPIWgN4tPqShBwVECIiIsfj98HC22Ddi5DY1e4k7AnC5nNbBvx1Ivy0B2Ij4NzmMLgxxAXh8qkeAx3r2I/tuTBjE7z2LUz8AW4cDGP7QUwNeF7ZayG+OUQlup1EJKCC8C+giIhINfEVwbQekLUS6g2ABmcE3xKtBcXwyjx49Rub/ezSd/BjQ+QlQGo8XHMSnJMLn26Ep2fDW9/BA6PgtPbu/byKMqFwN5z0uDvXF6lCIfLXQ0REJMB8BfB5V8hZCw3PgnqnuJ2o/OZvhPveg70F0Ks+jGxlhwGFokbxcGNn+GU/fLge7nrXzuf479VQv9aJHx9oOetLcw2t/muLVLEArMUmIiISYnwF8HlnWzykDg++4iG/GB6dCje8Yd+Bv60rXN0xdIuHQ7VLhj/3hFGtYPV+GPk0TF1W/TlyN0BEAiR2qv5ri1Qx9UCIiIgcylcIn3ex7yA3HgnJPdxOVD5rd8Hv37Cbsw1pDMNb1qxVlaqD12OHaXWpC2+vgfsm203x/noexFdDEeU4kJMGTc8PviFvImWgHggREZED/CV2zkPOOkgdEXzFw6c/w2XjIL8Eft8VLmgTfsXDoerHwR3d7YTxactg9NOwYU/VX7dwL/hyocFpVX8tEReogBAREQH7rvGCW+3OwQ3PgTo93U5Udj4/PP4F/PUjaF4L/twL2ie7napm8Bg4twXc1g3yimHMSzD3l6q9Zt5me5tyatVeR8QlKiBEREQAlv8D1o+3L/rq9XM7TdnlFcFVz8NbP9hlWX/fDWpHuZ2q5mmbBHf3gpRYuH0SvDO/6q6VtxW8cVCrbdVdQ8RFKiBEREQ2vgXLHoSkblA/iIadZObDZc/Cin1wURu4sA14Neb+mJKj7ZCmTnXhX9Pg2dm25ynQ8rdCg9M1/0FClgoIEREJb3t/gB+ugfgWdt5DsLzoS8+Fy56DLTlwbScY1NjtRMEh2gvXd4JTGsH4b+zQr0AWEb5COweibp/AnVOkhtEqTCIiEr7yd8LssyCiNjS9GDxBMuE4Iw/GvgB78uGmztChjtuJgovHwKVtIcoDb/9oC4g/Dw1M8Viw097W6VX5c4nUUCogREQkPPlL4It+9h3j1ldARJzbicomtxCufBF258FNXTRZuqKMgdGt7b8nzrfLu952euXPW7DL3iZ3r/y5RGooFRAiIhKelj0IeZugyWiIaeB2mrIp9sHvXoKt2XB9ZxUPlXWgiCjy2+FMdRPgsr6VO2fBLvDGQGxqYDKK1ECaAyEiIuFn52xY8ajd5yGpq9tpysZx4J+fw5r9MKY9dK7rdqLQYAxc0tZ+P/9vWuWXeC3cC9H1g2cujUgFqIAQEZHwUrQf5p4PUXWh0VC305TdxPkwZTGc3Qz6NXQ7TWjxGLi6IzROgLvfq9xmc4V7oNE5gcsmUgOpgBARkfCy6E4oyYGmF4AnSPZL+GkzPDEdutSFYS3cThOaor1wQyc7sfqm1+xck3JKiPaBLx9qt6uCgCI1hwoIEREJH9umwsY3IGVg8IxRz8iDO9+GurFwRQf7brlUjeQYuOYku7rV7a+W++GNk0rsPxLaBDiYSM2iAkJERMJDSS58dzlEp0DKILfTlI3jwB9eg5xiuKYjxGrtkyrXNgmGNoeFu2HqsnI9tGHtAwVEy8DnEqlBVECIiEh4WPYwFGdC6nDwBMkL8c+WwtK9MLwlNK3ldprwcU5zaFkbHv4EdmWV+WENavnsP+KbV1EwkZpBBYSIiIS+rF9g9ROQ1B3im7mdpmz25cCjn0Gr2nBaE7fThBePgbHtocQP/++NMu9UXb92iZ1XE5lYxQFF3KUCQkREQt/XI8BEQMMz3E5Sdk/MgEIfXNZe8x7cUD/OTlhfvg9mry7TQ+rF+yCilpZwlZCnAkJERELbrjmQ/QuknAoRCW6nKZtFm+z4+zObQoMg2SE7FA1pDKnx8MjHUFB8wsPrJvggUkPNJPSpgBARkdDlOPDdFRBRG+qe7HaasvE78LcPIDkazgqS4VahyuuBC9tAeiG8+f0JD68T54N6/aohmIi7VECIiEjo2v455G+D+oPBE+l2mrKZtgy25tiJ01Fet9NI2yToWg/Gz4X03OMemhjrs7tQi4Q4FRAiIhKaHAfm3wSRSZDcze00ZVPih/9Og8bx0EsvRGuMES2hyAevfXvMQyI8DgnRDkTXq8ZgIu5QASEiIqFp55eQv91uGmeC5J38qUthbwGc20ITp2uSBnHQuwFM+hH2Hb0XolaM3/4juk41BhNxhwoIEREJTfNvtpOmk4Kk98HvwLMzbe9Dl7pup5HfOruZ7SGa+MNRv3ywgIhKrsZQIu5QASEiIqEnYxnkbrATp4Nl07hv1sKuPDijqZYBrYkaxNm5EBN/gLyiI76cEFVaQEQmVW8uERcEpIAwxrxqjNltjFkeiPOJiEhwc71dWPOM3fchuZcrl6+Q56ZBYhT0SHE7iRzLkCaQVwKfLz3iS3EHC4ja1RxKpPoFqgdiAjA0QOcSEZHgNwG32oWiTNgwARI7Q0SsKxHKbdM+WJMBA1Pt0qFSM7WqDY0T4NU5R+xOrQJCwklA+nUdx5lrjGkRiHOJSPVxHMjLg8xMyMmx/y4ogMJCKCkBn88eYwx4vRARAdHREBMDcXGQkAC1a0N8vEZcyOFcbRc2TQKnGOr0duXyFTJlsX1Lr19Dt5OUS3ZBLFv212dHVl125ySRnlubrII48opjKCyJxOf3YAxEeHzERBSREJ1PYmwOdeOzaFBrP40T99I4aQ+RXp/bT6VsjIEBjeC9tbByB3RKPfil2KjSgiIySDYrFKmEahsYaoy5EbgRoK43sbouKyLYIiArC7Zvh127YM8e2LsX0tNtoXA0B4oGY+zj/X77cTQREZCcDHXrQkoKNGwIjRpBUpIKCzm+Q9uGZs0CtGna8kfsWvyxqSc+tibw+WHyAuhYBxKj3U5zVJn58fy0pS2Lt7Zl2Y5W/Ly5GZsyGpOen3TU4w1+Ir0lRHh8OA74HC9FvqijHusxPhrV2kPL5C30bLGJrqkb6NX0Fzo12lgzC4ue9WHKevh4yWEFRExEaQHh1c7hEvqqrYBwHGccMA6gZXSqc4LDRaQSHAf27YP162HTJtiyxfYwgH1Bn5wMAwZA27bQogU0bWpf9Kek2K/VqgWRR9lzq6QEsrNh/35bgOzYAVu3wsaNsG4dzJsHa9b82rMfH2/P3bw5tGplz6+CQg51aNvQu3fvyrcNWWvsxnENzwqeX7aFmyCzCC6oOb0Pe7ITmb22J1+t7cFXq7uwdl9znNJRz/Xi0mlTZxPntJlL08SdpNbaRYOEfdSN209STBa1onKJjig64tvvdwy5RbFkF8WzPz+Rvbl12JlTj21ZDdicmcqG/c0Y/+155JfYYWcxEQV0bbCaszsv5fS2P3FKyxVERxZX97fiSHER0LkufLYE/jwUIuz3JSay9B2WiHj3solUkyBZmkJETsTvh7Q0WL0afvnFDksCSEyEUaOgXz/o0we6doXYCg4LP9DTkJxsC4KjKSiAZctg4UL44Qf45BObCWxh0q4ddOgALVvaHg6RgEqbaG8Tu7ibozy+WA5RHujk7v4BK3c0Z8rSU3l/0Sks3dkBBw+1onLombqC89rPplvD1XSqv5Z6cRkVOr/HONSKzqNWdB6ptfYc9Ri/Y9iUkcqyXe1ZsuMkFm3vxL9mjOWRL64iNiKfAc0XcfnJ8xjZ5TvqxmdV4tlWUs/6sHgPLNgI/VsDEH2wByJI5t2IVIIKCJEg5jiwbRv8/DOsXGnnMERG2hf3//d/cNZZx36hX1ViYmyh0qcP3HKLvW/TJpg5E6ZPh08/hUWL7HEnnQTdutleimB5s1hquLXPQ3wLiKzldpKy8flhxjLoVBeiqr+i3ro/hbcWnsWE785kzT77x6Jbw5Xc2X8Cp7ZYQJcGvxDhqb5hRB7j0DJ5Gy2TtzGyw2wAsgvj+GFrd+am9WX2hn5cO3EgEZ4SBjZbyA2DZjK66zfERh25rGqV6phsi75Zqw4WEFERDj4/eINl2WCRSgjIb7kxZhIwBKhnjNkKPOg4ziuBOLeIHKmw0BYNCxfa+QwREXDhhXDppTB0aMV7GKpK8+Zw/fX2o6DAFhPvvgvvvQc//WTnTvTqBd2717zsUjGutAuZq6BwL9TpU6WXCagV2yG7uFo3jvP7DdNX9eU/M85nTlpf/I6XXqnL+PtpTzG07VwaJOyrtixlUSs6j7Naf8dZrb/j4dNh+e52TP1lCJ+sPp2xbzxAregcrjn5C2499WM6NNhcPaGivHbOyoxlcP95YAyRXoeiEoP+hEk4CNQqTJcF4jwicnxZWfD99/ZFd1ERpKbCuHG2cKgdJCsHxsTAiBH248UX4YMP4IEHYMYM+OorW0SccoqdgC3By5V2YetH9rZ2h2q/dIXNWwsG+2K0ihUUR/H6/HN47ItLSMtoSkr8Pm7uM5FLOk+jedL2Kr9+IBgDXRr8QpcGv/CngeP5YUt33l12Hi99O4Jn5l7IkBY/8OCIdxjcZknV92qeVAd+3gtrd0O7BkR5ocinAkLCg/rZRIJAZiZ88w0sXmyHLY0ZA3fdZYcJBbOEBLjmGvuxZAk8+SS8+abtWenWDQYNsvMtRMpk7QsQ0yi41uH/YjE0rwXxR1m1IEAKiqN4cd5IHp1xGXty69K1wSqeGvYw57b9umauclRGHuNwSrPFnNJsMXvzkpi4dARvLBnNac88Sc9Gy3ls9Guc2X5R1RUSB4q+79dDuwZEeB1K/BqLKeFBBYRIDZafD3PnwoIF9vMbb4R77rErJ4Wa7t1hwgR45BF44gl47jlYutQObRo82K7oJHJMRfshbyuknOp2krLLKYTN2XBmgJav/Q2/3/DGgnP4y8fXsiOnPv2b/sST5z5C/6aLQ27OUb24DG7v9yY39X6X91cM5YX5Yzn7+f/Qr8linhnzIr2brQn8RZOioUEcTF8EV59CpMehOHjrMZFyUQEhUgP5/Xai8Vdf2TkDV18NDz0EgVomvyZr0sT2RPz5z/Dww3aI1tKltojo21crN8kx7JwNOJDQ2u0kZffzFvADbQK/N9J3Gzpx48Q7WLG7Hd0aruSJof/ilGaLA36dmiY6oogrun3CxZ2mMWnZcJ794Sr6PPESV/edzv+NfIkGtfcH9oJtEmHRbvD58XrApx4ICRMetwOIyOF27oRXXoGpU6FBAzu057XXwqN4OFRqqp0jsWKFXaVpxgwYP96uOiVyhF2zwRMJcY3dTlJ2S7bY+Q8tA1dAZOQlcMOkuxnw5HOk5yXx1LCHmXLZrWFRPBwqOqKYa3pM4atrx3JT74lMXHgGbf/xJuO+HY4/kC/yWyVCgQ/W78HrsaswiYQDFRAiNYTPZ3scxo+3cx4mToQNG+y+DeGsY0e7r8WHH9plal95xa7idKwdtCVMbX4f4pqBCaIuqrnLIDUeogOTeeqKk2n/jwm89sO5XN/rXWZeczUjO8wOueFK5VErOo97B41j2lXX0illLTe9ezcD//MfNqU3CMwFWpTOt1m21fZAOGH8zZawogJCpAZIT4dXX7XzHS6/HDZvhssu094IBxgDF1xgd9S+/nr47jt4+WW7hK0IhelQuMcWEMHCcWBLNjSt/H4V+UVR/P69Ozjvpf8jOSaLKZffwv2DXyA+Kj8AQUND6zpbmHjxXTx65uMs29WeLo++wqRFp1f+xPVi7M7UK7bjNeqBkPChAkLEZStXwksv2SLi/fftKkR13N2QtsZKTLRzIj77DLKzbW/Nzz+7nUpct/d7exsfRAXEzizILYGmCZU6zdrdjen52PM8P2801/V6l4/H3kSXBr8EKGRoMQYu6/o5U6+8jvb1NnL563/jhkl3U1AcVbmTNk6ABWvxeAjs8CiRGkwFhIhL/H47FOf996F+fVizBi66yO1UweG882DdOmjcGD76yM4X8Wn1k/C170fAQGyq20nKbu0ue5ta8QJi6oqT6fX4S+zMSeG10X/mr4NfIDqiOEABQ1fTxJ1MuuRObunzNi9/P5yT//00W/enVPyEqfGwMxcPDk7gYorUaCogRFxQVGR3Yv7uO7j5Zli/PvwmSVdWo0awdi38v/9nl7l9+227YpWEobRJEJ0Cnkq8k1zdNpSOv2sUV6GHP/nVRQx/6V80S9zOp2NvYkjL+QEMF/oiPD7+fOp4Xhp5Pxv2N6Xnv19k0eZ2FTtZo3go8hOLhjBJ+FABIVLNcnPh9dfti99nn4UXXoCoIHrdU5NERMB//mPnj2zaZFerys52O5VUu4IdENvI7RTls3Ev1IqEuPJtIOf3G/445RbumnIbZ7eZx3uX3k6TxJ1VFDL0nd3mWz4YcxtR3mIGPfUUX6yqwO6cDWwRGOv4cdQFIWFCBYRINcrKsi9yd++Gjz+G3//e7USh4Xe/gy++gIwMW0zsD/BS71KD5e+EklyIaeh2kvJZthFSYsv1EJ/fw3WT/sx/v7qUq7t/yHPD/05cpLrdKqt9vY1MuexWWiRtZfhL/+KDxYPLd4L69ucY43fwozkQEh5UQIhUk6wsu9NydjbMng0jRridKLSceSbMm2eHMU2YoCIibGQut7cx9d3NUV5786Fe2QsIn9/DhS/ex4Qfz+XO/q/x4GnP4PVovEygpMSnM+mSO+necBWXTvgb7/50WtkfnBAJ0V5iHAdNgpBwoQJCpBrk5MAbb9jhS19/Daee6nai0NSnD3z/PRQX22FiWVluJ5Iql7HC3kYHUQFR7IOsIqgTU6bD/X7Dje/8kY9Xn8WfBo7jjv6va4nnKlA7OpcJF/yZ3qnLGPv6X/l46YCyPdAYqBNNFA5+FRASJlRAiFSxwkI7wTczE2bNgn793E4U2rp3h2++sT0Rb70F+VoKP7RlrQZPDETEu52k7HZn23eqk6PLdPi9n97Iqz+cx+39JnBr34lVmy3MxUfl88ro++jSYA2XvPYgX6/tVrYHJkUT5ffjeMpWFIoEOxUQIlXI74cPPoBdu+CTT2DgQLcThYdevWD6dLu3xrvvatfqkLbjC4iuG1y7Lu4u7RpLPHEB8fw3o3h81mVc0W0Kd/afULW5BICEqHxeHX0vzRK3M3LcI6zeVYYl8pKiKd90eJHgpgJCpArNmGH3K3jpJTj3XLfThJchQ+ymfJs22X0itDpKiCpKh6gg23lxb469TTz+8muz1vTk9g9u5/RW3/H3054Jqhop2CXHZvHa6HuJ9BZz7nOPsj/vBPt11I4iAv2dkfChAkKkiixbBj/+CLffDjfc4Haa8HTZZXD//bB4MSxa5HYaCThfERRnQlSy20nKZ19pAVHr2AXE5vT6XPzKg7Sus5mnhv1DE6Zd0CRxJy+O+Bvbshpw4Uv3H3+X6VpRGCBCFYSECRUQIlVg71749FO7Odx//uN2mvD28MPQpo0d0rRTy+WHlrwt9jYqydUY5ZaRZ2/jjz7opdjnZfRLf6PYH8GLI/9KQpQm8rild+PlPDDkWb7a2J8nZl967AMT7M8yQrOoJUyogBAJMJ8PJk+2m5x99529Ffd4PPbnEBsLH35oV2iSEJG32d5GJrqbo7wy8iHWC96jv6P9yBdX8tOOzvzrrCdombytmsPJb13R7WPObTuH+z+7/ti7VcdHggGveiAkTKiAEAmwb76BHTvsykuNG7udRgBSUmzxsHcvfPWV22kkYPK22ttgKyCyCiD26L0Piza3459fXMkFJ01neHv9stYExsC/znqCenH7uey1+ygsPsrPLs6+U+RVD4SECRUQIgG0Z48tIC6/HEaPdjuNHOrss+HGG+GHH2D7drfTSEDklb47H1nL3RzllVMAsUd2TRb7vFz5+p+pG7efB4c860IwOZbEmBz+ddbjrN3Xkn/NHHvkATH25+lRD4SECRUQIgHiOPD55xAdDU8+6XYaOZp//xvi4+Gzz+wSuxLkCnaCJ8p+BJNtuyDGe8Tdz84dzao9bXjo9KeoHZPjQjA5niEt5zOqw0z+NfNy1u7+Tfdy6c9TQ5gkXKiAEAmQlSvtkqH/+58dMiM1T2KiXVJ3xw74+We300il5e+EiBMsr1kTFfog+vACYk92Ig9+fg2DW/zIOW2+cSmYnMhfBr1ApKeEWybeevgXor1gwIMKCAkPKiBEAsDns7tM168P11/vdho5nssugyZN7FwITagOcvvmB9cO1AcU+SHq8Ob3H19cRV5xLH8d/Jz2e6jB6iek8/uT32LWhgGH71IdaUsHrbYr4UIFhEgALF4M+/fDa6+B98iRCVKDGGMnuGdnw4IFbqeRSinJBW+c2ynKr9gHkb/+odicXp8X543kok7TaFN3s4vBpCx+1+MDGsTv4e7J1/26cZwx+FEPhIQPFRAileTzwbx5dsUl7TYdHAYNglat7PKu6oUIYr588Ma6naL8iv0Q+Wvz++9ZlwHwh35vuJVIyiEmsohbT36bhdu68vW67gfvd9AkagkfKiBEKmnlSsjMhBdeQEMPgsgrr0BuruZCBDVfPkQEYQFR4kCEbX735iTy8nfDGN1xBo1r73Y5mJTVpZ2nUi8unb9/etnB+xxATYCECxUQIpX0449Qty6cd57bSaQ8Bg+GRo1g/nzQm4ZByFcETgl4YtxOUn4+P0TYl5rjvzuPQl801/d+z+VQUh7REUVc2X0KX6edzC+7mwDgGDD6WyJhQgWESCXs2AHbtsGDD9odjyV4GAOPPmr37tisYefBpzjT3nqDsYBwwGPw+w3PzR1J/6Y/0bbuJrdTSTld1uUzIj3FvPTtCOBAD4QqCAkPeskjUgmLF9tJ01dc4XYSqYhLLoGoKPtzlCBTnGVvPdHu5igvx7EFhNfw1doebMtqyJgun7mdSiogJX4/Z7T+jtd/OJtin1dDmCSsqIAQqSC/H1asgAsvhORkt9NIRcTF2eJv1SooKXE7jZRLSelGa8G2idyBN6g9hrcXnklCVC5nt57naiSpuAtO+oJ9+cl8uaaXLSDUASFhQgWESAVt2gR5efZdbAlel14KRUWwfr3bSaRcDhYQke7mKC+f3SjAh4fJSwZyZutviYkscjmUVNTgFvOpFZ3D+4uHAOqBkPChAkKkgtasgYgIGDrU7SRSGaedBtHRsHq120mkXEry7G2w9UD47VvUmzIakllYW7tOB7kobwmntfiBT5b2L71HXRASHlRAiFTQ+vXQogXEB+FGuPKryEgYMcL+PLUaUxDxHSgggqwHorSAWLOnGRGeEgY0W+RyIKmsIS1/ZF9+Mo5eUkkY0W+7SAU4eNi7F266ye0kEghnnGF3pk5PdzuJlFlJvr01Ee7mKK/SKnXNzsZ0a7iKWtF5LgeSyjpQBPrxagiThA0VECIV4gXsjsYS/A78HLdscTeHlIO/wN56gqyAKLU9twF9Gy91O4YEQP2EdFokbQXH63YUkWqjAkKkAhy8eL3QvbvbSSQQOnSw8yC2bXM7iZSZv3TisQmyF22lPRB+x0P3RqtcDiOB0q3hKhyC7HdRpBJUQIhUiIeUFLuHgAQ/jwcaNIBdu9xOImXmC9ICopSD4aT6a92OIQHSqf5aHIzdjlokDKiAEKkABy+nn+52Cgmkc8+FvXvdTiFl5hTb22ArIEon6kd5imlcSxVrqGhXNw0wOI5eVkl40G+6SDk5GMDQrp3bSSSQ2raF/Hz7IUHAf2Dnv+BsxpJiszB6szpkNE86MP4xOH8fRcpLv+ki5WZb/WbNXI4hAXXg55mZ6W4OKSOntIAwwdmM1Y7OdjuCBFBqbdubpB4ICRf6TRcpN/vfplEjl2NIQB34eebkuJtDysjx2dsgfRu/VpSWbw0lUd4S7Pi04Px9FCmvgBQQxpihxpg1xph1xph7A3FOkZrLNhB167ocQwKqTh17qyFMgVOlbYPjP3CVgJ62qh3YrDA2ssDdIFIFHE2ilrBR6QLCGOMFngPOBU4CLjPGnFTZ84rUdAkJbieQQKpVy94WFbmbI1RUedsQpAVEUbFtdqO8xS4nkcBzSufIiYS+QOzA0xdY5zjOBgBjzDvAKGDlsR6wo3gfj+2cEIBLi1S/mzkHgGuuuZzY2O0up5FAKSpKAj7iu++eZNmyj1xOExLK3TasWbOGIUOGlOnk1/RK45peMOT61wMQtfpEFcUwA1iycz//fG+J23EkgGZcB47jlPl3WCSYBWIIU2Pg0P1bt5bedxhjzI3GmIXGmIW+A2NXRYKa43YACSj783Q0BCFQyt02FBeHw7vy+v0SkeAXiB6Io/01POKVleM444BxAC2jU517G14TgEuLVL+sTTvwA6+/PokOHdxOI4GybRs0aQIDBtxBr153uB2nWjz0UJW+mC1329C7d29nzpw5ZTv70r/D8oeY8/LVQTWROn9vCZzxT3o0SuLK0d3djiMBtQBjDGX+HRapgUwZ/54GogdiK9D0kM+bABrXISEvT4uohJTcXHsbGelujhBStW3DweVbg6snMCba9sAX+fSLFnoMwfb7KFJRgSggFgBtjTEtjTFRwBjgkwCcV6SGsg1EerrLMSSgDvw8Y2LczRFCqrhtCJ5eh0MdeHMvvyTa3SBSBQzGqICQ8FDpIUyO45QYY24DvgC8wKuO46yodDKRGss2EDt3uhxDAmqX3QdKq2sFSJW3DR5v6YX8QbmZXG5hnNsRJIB8fg+2qPWf6FCRkBCIORA4jjMVmBqIc4nUfLaB2LLlBIdJUNm82d7Wru1ujlBSpW2DOdB8BecLtqxCVaqhZFduXVISwZjg/H0UKa/ge9tGxGUGB3BYt87tJBJI69fb+Q/x8W4nkTI5UEA4wfmCLaOgltsRJIA2ZxxYYExDmCQ8qIAQqQCDnxkz3E4hgfTZZ1CvXlAt6BPePKWTkIOtgDg4ByKWvXlJrkaRwFm7rzkAxmiZegkPKiBEKsTH7t3gD7LXLnJ0jmPntNSv73YSKTNPlL0N0n2FDA4rd7dxO4YEyKo9rTH4UQ+EhAsVECIVYPBTVASrVrmdRAJhyxa7jGtqqttJpMyCtYA4pIvr550dXQwigWR/lv5gXRxMpNxUQIhUSAkA8+a5HEMC4ptv7G3Tpsc/TmoQb+l6u06JuzkqKCVuHwu3dXE7hgRAVmE8q/e2whBkxaxIJaiAEKkQP7VqwaxZbueQQJg92+7/0KCB20mkzLyx9tYfZAVEaQ9E+/pbWbCtC0W+gCyGKC76cWs3/I5X8x8krKiAEKkAA7RuDZ9+CiVB9vpFDuc48P770LIlePQXMXgcKCCcYndzlJfHFhDtUraQXxLLgm1dXQ4klTU3rS9xkfl48GkGhIQNNZciFdSuHRQUwNy5bieRyliwALKz7c9TgkhE6UZs/iJ3c5RXaQHRpu42or2FzFw3wOVAUhl+xzBz3QDO6big9B5NgpDwoAJCpILatLH7Brz3nttJpDLef9/2PLRv73YSKZeI0g07/MHZAxFlihncYj7T1g4u3cVYgtFP2zuxKzeFC7vZd5LUAyHhQn+1RCooMtK+6Hz9dSgsdDuNVERJCbz0ki0GY2PdTiPlElG6EVuQ9kDgwHWnzmJ3bj2+29LD3UxSYZNXnk1sRD6jun5r+x7UASFhQgWESCV0726HMU2Z4nYSqYjp0+3wpe7d3U4i5RZZWkD4gqx6N8YWET4/I7t8R2J0Fu8uO8/tVFIBecUxfLbmdC7u+TUJ0fkYQFsDSbhQASFSCa1aQXIy3HOP20mkIu66CxISNP8hKEXWtrf+ICsgALwGfA4xkUX8rv8XfLFuELtz6ridSsrpo1VnkV2UwE2nfArYzgdHXRASJlRAiFSCMdCnD2zebCfjSvBYuRLWrYPevcHrdTuNlJs3FowXfAVuJym/CFtAANx26hR8fg+vLxntcigpD79jePWni+hcfw39W64AwDjgqH6QMKECQqSSeva0ewhcdZXbSaQ8HnvMzmPp08ftJFIhxoAnBnz5bicpP68HSuxgl9Yp2zmn7Te8+fNosgrjXQ4mZTVj3UDWpzfnL+e+e3BzcdsDIRIeVECIVFJ0NJx8MqxeDYsXu51GyuKXX+Ctt6BXL4iLczuNVFhEbHAWEJEeKP51tPw/z3+L7MIEXl10sYuhpKz8juGp76+mRdIWLu7+9cH77RwIdUFIeFABIRIA/frZXoiL1f4HhfPPh4gIGDjQ7SRSKd54KMlzO0X5/aaA6Nl0LUPbfM3Liy5hb16Se7mkTD5edSar97bh0VGvEeEt3X3acfAAjlEBIeFBBYRIAMTEwKBBsH49TJvmdho5nm++gVWr4JRTIF4jRoJb3d7gy3U7RflFeaHId9hdT176MoW+KP777bUuhZKyyCuO4fF5N9C5/hou7fHVr18osdOn/aofJEyogBAJkL59oW5duOIKu7Sr1DwlJbaXqHZtW0BIkItpACVBWEBEe6Hw8AKifYMt3DZoCu8sG87PO7WrYU313I9j2ZFTnxcuexaP55AZD0U+cMCvHggJEyogRALE64Vzz4X0dHjkEbfTyNH85z+waxcMHQpRUW6nkUqLaWjnQPhL3E5SPqn1ocB3xN1/P3cCKfHp/GXm3RT7tDRYTbN6T0vGLbyMq/tOZ2DrZYd/sfTn6VMBIWFCBYRIALVuDd26waOPwqJFbqeRQ61cCfffDx072g8JAXGN7W1Jtrs5yish+ogeCIDE2FxeuuxJVu5pywvzx7oQTI6l2OflT1/cS2J0Nk+c/8KRBxTYIlYFhIQLFRAiAXbOOXZzsnPOgdwgHF0RigoL4Ywz7IpZw4a5nUYCJra0gCjOcjdHedWOgbzio37p/K7zGNVhJk//cDWLt59UzcHkWJ78/ncs392e8Zf/h3oJmUcekKcCQsKLCgiRAIuNhdGjYd8+uO02cLQwuOvuvht27oSRI21xJyEivpm9LT7KC7qarHasfcHpP/ofhwnXPEXDWnu4feoDZBboF9Ztc9N688L8y7mu3+eM7jbv6AflFYMDPo8KCAkPKiBEqkDLlnZVpgkTYNw4t9OEt7ffhmeftUvtttfc1NAS39zeFgVZAZEcZ3ccyz/63I2kuBw+vP4hduXU465p9+Pzq6l2y9bMhtw59QHa1U3j6YuePvaBObZHST0QEi70V0mkigweDG3awK23wpw5bqcJT/PnwzXXQPPmcOaZbqeRgIuIg4h4KN7vdpLyqVO6fnB20TEPObnFKp6+6Bm+2tiff8+7oZqCyaFyimK54eNH8DkePr3lAeKiCo99cLYtIErUAyFhQgWESBXxeODCC6FOHbs606pVbicKLxs2wOmn2yFLF19sV8mSEBRVBwrT3U5RPvVq2dusYxcQADcP/IRbB05h3MLLePvnkdUQTA4o9nm57bO/s3ZfS96/9iHa1t92/AdkFVECOKofJEyogBCpQjExcPnl9sVrv36waZPbicLDjh3Quzf4/TB2rDaMC2kNz4SifW6nKJ+U0gIi4/gFBMBTFz7L6S2/52+z72DqL4OrOJgA+B3Dn2fcw9dpJ/PCpf/l7I4LT/ygjEKKAdUPEi5UQIhUseRku7lcYSH06AHbTvBGllTO7t12Kd2cHFu81avndiKpUrXbQ0kO+IJo98YGte1txokzR3h9fHLrQ/RotJI7p/6VL9f3r+Jw4c3vGP765V18tOpsHjnvZW445fOyPTCjkCLjwfiD6PdQpBJUQIhUg4YNbRGRmwtdusDmzW4nCk07d0LnzrB/vy0emjRxO5FUudqlS50W7nE3R3nERkJCJKQfZ0z9IeKjC5h1x710TFnHrZ8+zIx1A6s4YHjy+T38ZeYfmbRsJPed9RZ/Ofutsj84vYAiY9AcagkXKiBEqkmTJnDllZCXZ4sIzYkIrA0b4KSTICPDFg8tWridSKpFUmd7W7Db3RzlVTcG9uaX+fDE2Fy+vutPdG7wC7d++hAfrjinCsOFn8KSSO6c9lfeXT6c+89+k38Of7nsxUB+CeSVUKgCQsKICgiRatSkiV0VyOeDnj1h7ly3E4WG+fOha1coKICrrrLL6EqYiG8Bnmgo2Ol2kvLp1AL2lL2AALu869z/dzf9mi7m7i/u4+kfrtQ+MwGQWZDANZP/zWdrTuffo17gkeGvlK8QKP05FniM5kBI2FABIVLNGjaE666zE3tPOw1eftntRMFt4kQYMAAiI+HaazVsKewYAzENIX+H20nKp0VdyCiEIl+5HpYQnc/sO+/lyj5f8L/vruOuafeTXxxdRSFD3/r0Zpw/8QV+2tGJt658hD+d8W75T7IzD4B8Y/AYVXQSHlRAiLggORmuv96+U37DDXDjjfbdcym7oiK48067ylLjxvb7qQnTYarZhbYHwinfi3FXtUqxm8mVvvgsj6iIEl6/4l/8c/h4Pll9Bhe98yybMlIDnzHETf1lMOdPfIHswgRm3fZHxvb5smIn2pkLHkMeHrQNhIQLFRAiLjmwxOuAATB+vB2zv3Kl26mCw7p1tvh66ino29fOLdFSrWGs7sm2eCjY5XaSsmtb395uz63Qw42Bv5z9Np/ddB/bshow4q1xfLL69AAGDF0FxVE8MOsOfv/ZQ7Sps4nF997IwNbLKn7CHbnQIA4fRgWEhA0VECIu8njsDsmXXw7Z2Xb50SeftPsXyJEcB156yU6WTk+HSy6xm/Rpk7gwV690adPcIFrerGkdiPLA1pxKnWZYpx/5+b4baFs3jTum/o07Pv8r+/NrByhk6FmyowPD3x7PWz+P5o+nv8uPf76dpsmVXMFraw70bI3fQUOYJGyogBCpAdq2hVtugVat4K67bG/EihVup6pZfvnFfn9uvhmaNrW3HTu6nUpqhPimEJkIeUFUQHg90CQBtmRX+lTN6+zixz/fwcPDXmHq2iGc9foEPl19uiZYHyK3KJZ/fn0LF77zHHlFscy49Y88cf4LREWUVO7EmYWQWQQnNcLnN3j1qkrChH7VRWqIhAQYMwbOPx/27rWrCt19N2RluZ3MXTk58Ne/2l6HHTtg+HC7p0ZiotvJpEZpPBJyNxFUr5oHnGQLiJLKdzlGeH08MPRNFv3pRhrX2s3tU//GlR/+h7X7mgcgaPByHPhszRDOmvA6Ly+6lOv7f86qB37HWR0WBeYCaaUFYOfG6oGQsKICQqQGMcYOY7rtNnv7n/9Ao0bw7LN20nA4KSmxc0MaNYJ//hM6dbLfl1690FrrcqSGp4MvDwqDaD+I7k2hxIEtlRvGdKiujTfw03238sxFT7FsVzvOfeNV/vrlXezOqROwawSLRds7cfG7z/CHz/9Ocmwm8+68jZfG/JfE2IrNOzmqjZkQYeCkRpT4IELDKSVMqIAQqYHi4mDkSLtCU0oK/OEP0KCBXfI11AuJkhJ44w37fG+8EZKS7PKso0fbXhqRo2p4pr3NWe9ujvLo0czerssI6Gm9Hj+3DZrC+gev4JZTP+bd5ecx+NWJPPr1zezJDf1CYvGOjlw75V9c9M5zbM5MZdyYx1l6/00MaLU88BdblwnNa0NUBD7H4FUPhIQJFRAiNVhqKlx9tV2qNC7u14LiP/+BzEy30wVWdjY884x9fldfbfd1uPRSWzw0bep2Oqnx4ppAdH3IXud2krKrmwAN42BtRpWcvl5CJs9c9DSr77+Ki3t8zSs/XcypL0/ir1/exfr00PpP5XcMX204mbHv/4cLJr3A4h0n8a8R49j497HccMrneD1VsDJFXokdgnZmNwBKfEY9EBI2ItwOICLHZwy0aQOtW8P69TBvnp0bcd99dkO6G26wu1oHq2XLYNw4O1ypsNBuBDd0KLRrp6FKUk6tfwcrHwdfAXhj3E5TNmd0gfcW2A3loqrm1WfrlO28ceW/eOCcN/j3rMt4Y/65vL10FKc2n8/lXT/ljFbfEekNoj00DrEvL5HJK89h4tKRpGU0oWHCHh4f9QI3D/yEhOjy7fRdbmv22708+rUCoNgHkV71QEh4UAEhEiQOFBJt2tjJxPPn2yFNL75oh/v8v/9n37FvHgRzJrduhfffh8cft8/F67WTpPv21U7SUgmNR8HK/4PsXyCpq9tpymZAG3j7R9sL0alulV6qbf1tjL/sCR457xXGfTec5+eO5JZP/0Gd2AxGtJ/FyA6z6N5oVY2fCJxfHM3sjf34eNVZfLWxHyX+CHqlLuOfI1/jou5fV35lpbJamQ6xXuhqe3NK/IYIT83+3okEigoIkSDUqBGMGgXnnGPfwf/5Z7jnHvuRmgo33QTDhtmeCU8NGKjo98PSpTB1Krzwgi0gwD6Pc86xK07FxbmbUUJAvZMhIgGyVgdPAdG7BUR7Ydm+Ki8gDmhQez8PDH2T+86ayBer+zDhx6FMWjaC15dcSMOEPZzZ+ltOa/kD/ZouIS6yoFoyncjO7Hp8ndaX2Rv6M3dTHwpKYqgfv5c7hnzINX2n0zl1Y/UG8juwYh8M7ggR9o9sUYkhKkIFhISHShUQxpiLgb8DHYG+juMsDEQoESmbmBjo08d+pKfbnaxXr4YHH7QfsbF2o7WBA6FfP+je3d5X1QoLbVHzww92yNXnn0Nenv1ao0Zw2mm2x6FevarPItXPtbbBeKDV72DtC+ArBG90tVy2UqIjYHAH+P4XuKQt1bmVcYTXx3mdfuC8Tj+QlR/Hx8sGMHnpID5ccQ5v/Xw+kZ5iujdaSd8mS+nVaDldG66mblzVT77yO4a0jMb8vKMjC7d34cet3VifbrtWGyXs5rr+07io+9ec2npp1cxtKIv1mZBTDKd3OHhXkc8Q6QUcv/1dFAlhle2BWA5cALwUgCwiUgl16thCYeBAyM218yU2boTZs2HyZHuMMfZF+xlnQIcOdjhUy5Z22FDDhhAVVfbrFRfDrl2wZQukpcG6dbZ4mTUL9uz5dTft2rXtdVq1svM4tJJSWHCvbWh+GfzyDGStguTu1X75Cjn7JJixwg5jap/sSoTasXlc2XcmV/adSUFxFN+s78LMNb35Ynl3Xpx/OT7Hzs9olLCb9vU20KbuJlomb6Vp4nZSa+2mfvw+EqLyyjxvye8Y0vMT2ZVTj21ZDdmc2YgN6U1Zu68lq/e2IqcoHoBaUTn0Sl3OzadO5ewOC+iSuqFmzI36aTdEeuDUtgfvKiwpDebLh4h4l4KJVI9KFRCO46wCMDXif7OIHBAfb4cFdS0dxZGdbYcN7dhhX/RPnw7vvHPk46KjbQ9FVJRdBcnrtbtip6WBz2eLhqIiKCiwH7+VmAj169udtRs1gsaNteFbOHK1bajXD6LqQMbPwVNAnNrWDmNasMu1AuJQMZFFnNVhEWd1WMS/R0FuYQwLN7dnweYOLN7ahp82teS7LT0p8h3+jkO0t5Dk2CwSonKJiywgyluMx/hwMPj8Xgp9UeQWxZJdmEBGQe2DRckBSTGZtKm7iav6zqBn07X0abaaTo3S3OtlOJZiPyzeA2ecBHG/fg8Ki0t/30tyVUBIyKu2ORDGmBuBGwHqevWKQqQ61aoFHTvajwNKSmD/fsjIsLtd5+TYYUYFBbZIKC62RcOmTXYeRWSk7T2IjrYfcXH288RE+5GcbI8RKY9D24ZmzZoF4oTQ4S5Y+gAUpkN0EOx7EBMJw7rC5z/DRW0gpmZNT4yPLmBw258Z3Pbng/f5/B62ZdRj475GbMmoz/bMuuzJSSI9rxZZBfHkFsVQWBKJz+/F4CfCW0xMRBYJ0fkkxuZSLz6T+rX20zhxL82Sd9Gq3g7qxme5+CzLYfk+u4TryG6H3Z1fXDpsqSTPhVAi1euEf6WMMV8CDY/ypfsdx/m4rBdyHGccMA6gZXSqZhmJuCwiwu65kJLidhIJRlXRNvTu3TswbUOr38HSv8H+n37dYK6mu6AnTFkMi3bDgFS305yQ1+OnWZ3dNKsTRDt/B8p3OyA5+uDyrQfkH+yBCNzO4iI11QkLCMdxguSvr4iIVJca3TbENYZa7WH/Yqg/BDw16x39o+rSGBrHwzfb4ZRG2gSlptqVZ/d/+P1p4D18ovTBHojibBeCiVQvLRMgIiKhp8+z4MuDzGVuJykbY+CGM2B7LvyS4XYaOZavt4HXwIVH7t6ZV1Ra9BUHyVAskUqoVAFhjBltjNkK9Ac+N8Z8EZhYIiISrGpE29DgdIhpAHu/BydIRs0O6wK1ImHWFreTyNFkF8GPO2FEN6h75HJyuYUHeiAyqjeXiAsqVUA4jjPFcZwmjuNEO47TwHGccwIVTEREglONaBuMgR6PQ+EeyF5b7ZevkOgI+N0gWL0fNuld7Brnq61Q4offDTjql7MLSl9SFWVUXyYRl2gIk4iIhKbmYyAyCfbMDZ5eiDF9IC4Cpqa5nUQOlV0Ec7fB0M7Q4ug7YGYdLCDSqzGYiDtUQIiISGjyREKPf0P+NshZ53aasomPhhuHwKr9dmM5qRm+2GR7H24ecsxDinweuxJT4d7qyyXiEhUQIiISulpdY3shds0Onl6Iy/pCUjR8tB78QZI5lO3Kg3k74IJe0KLucQ/NzPNAQRgubSthRwWEiIiELk8k9H4GCnYGz4pMMZHw52GwJQd+2Ol2mvDmODB5HUR54NbTTnj4/nwvFOyqhmAi7lIBISIioa3F5RDTCHbNAn+x22nKZlgXaFUbPt0AOUGSORQt2WuHk/3hTKgbf8LD9+V4Yf/PJzxOJNipgBARkdBmPDDwHbs+/55v3U5TNsbA42OhwAcfBsn8jVCTW2y/900SYEzfMj1kb64XSrSRnIQ+FRAiIhL66g+CxM6wdx4UBskqOW3qw42DYNFu+FkTc6vdh+ts788Tl0NE2V4u7c6OAF8+lORWcTgRd6mAEBGR8HDaF2C8sP3z4JlQfd2p9h3wd9ZAZqHbacLHot2wcDfcNBg6Nirzw3Zne+0/cjdXUTCRmkEFhIiIhIe4VOj5P8jdABlL3E5TNpFeePZqKPbD66vAFySFTzDbnQfv/AIta8P1p5broTsyI+w/ctMCn0ukBlEBISIi4aPtzRDXDHZ8AUWZbqcpm5b14G8jYV0mfLbR7TShrdAHr64Er4EXri3z0KUDdmRG2n9ka96KhDYVECIiEj6MB86YDfhh28fBM5RpRDe4uBfM2gILtUxolfA78PZq2JELT1wKjRLLfYr0PA94oiB7bRUEFKk5VECIiEh4qdUaej8PuRthb5CsygRwz7nQOhEmroH1QdJ7Ekw+22iXbb3rLBjQpoInMRBVF7Z9HtBoIjWNCggREQk/ra+D2ifZHarztridpmwivfDqjVAnBsYvt++US2DM3QZfboGLesFV/St3rpgUKNwTmFwiNZQKCBERCT/GwNnfQWQSbH4/eJbdTIqD126wY/OfWwp78t1OFPzm74QP1kGXunDfMPu7URnR9e1eEIX7ApNPpAZSASEiIuEpKhHOnG3X7d/8Hvh9bicqmybJ8Np1dkWmZ362qwZJxczfCW+vgfZJ8MrN5Z40fVQxDe1txtLKn0ukhlIBISIi4Su5O/R/HfI2w46pwTOpuk19mHCdXd716Z9he47biYLPN9tt8dA2CV6/FaIjAnPe2NICIv2nwJxPpAZSASEiIuGtxeVw0n2w/yfY973bacquXQN48wYwwFM/w7oMtxMFB8eBqWnw/lo4qQ689XuIjQzc+SPiITIR9i0I3DlFahgVECIiIt0egdqdYOdMyFjudpqya50Ck26B2lF2TsSPO91OVLMV++HN1TB9E4zqHtieh0PFNoYd0wN/XpEaQgWEiIiI8cC5CyGuOWybElwbgTVOgg9us0u8vr0GPlwHPr/bqWqe9AJ4agks3A23nQYPjbQrW1WFuCZQnAl526vm/CIuUwEhIiIC4I2BYUsgOgU2vwu5m9xOVHa1Y2HibTD2ZPh6mx3SlF7gdqqaY/k++Pci2JUH/7sUbhhU+dWWjieumb3d803VXUPERSogREREDohKgmFL7Rj2TRODZ48IsO+m/3ko/N+Fdo+I/1sIC3YFz8TwqlDog/fWwrjlkBwD798Cp3eo+uvGNrI7Uu+aU/XXEnGBCggREZFDxdS3RUREAqS9FVw9EQBDO8Pk30OjeDvef/wK2F/odqrqt2Y/PLYQvt0OV/aDj+6A5nWr59rGY4fDbX6/eq4nUs1UQIiIiPxWXCoMWwYRtWwRkb3e7UTl0yQZPrgD/t9Z9oX0owtg9hYoCYO5EZmF8PoqO6ncY+CVa+DucyCqCiZLH09CSyjaB7mbq/e6ItVABYSIiMjRxKXCeSsgui5sngiZQbQ6E4DXA1efAh+VTrD+aIN9R/7nvaE5rKnQB9PS4B/zYckeuOFU+Owu6NXcnTwJre3tji/cub5IFVIBISIiciyxDWD4aohtAls+hL3fB9+L7ybJMOkP8MxlduLwKyvgP4thxb7gey5HU+iDWVvgoR9h2iboWAc++QPcdjrEBHB/h/KKToHI2rB9mnsZRKpINffniYiIBJmoJFtETO8FO2dA0X5oNNSOcw8mg9rB1D/Cpz/D0zPgpeWQGg+nN4UeKRAZZM8nsxDmbbc7SueVQPtk+OsF0LWJ28ksYyChLWz/HHyF4I12O5FIwKiAEBEROZGIWDhvOSy5B1Y9AYX7oOlF9v5gEuGB0T1geFeYugyenwlvrYaP1kO/hvajfpzbKY/N78Av++H7nb8OxepcF/40Cro1dTvdkWq3h/2LYNdXkDrU7TQiAaMCQkREpCyMB3o8DrU7wvwbYP14aH4pxDRwO1n5RXrtTswju8H3G+D5aXaS9ZdboHkt6FkfutWDOjFuJ7VFw6YsWLIXftoNmUUQFwGXnwyX9oFmddxOeGzxLe1yrlunqICQkKICQkREpDxaX2uLiK/OhvUvQ+pwSO7mdqqKMQZOaQ2n3AZ7sm2vxHvfw5T19iM13s4paJcErRIhuop2bv6t/YWwLsOuILUqHbKLwWtslvtPgyHtIToIXsJ4IiChDaS9Db2fB081ff9EqlgQ/O8TERGpYVL6w4i18EU/2PYR5KZB6rn23eZglVLLrtp09SmwOR2+Wg2fzIc5W+0kZQ+QmgDNakHjeLvPRIM4SIis+K7OPsfumL07D7blwtYc29twYN+KuAgY1AEGtbVzOGrVgB6R8qrdEbJWwt5vof4gt9OIBIQKCBERkYqIbQgj18Gyv8OKf0L+Vmhygd2FONg1q/NrMZFXBIs324+5y+wSqd/t+PXYKI/d5TkxyhYTcZEQ7bHzLTylhYXPgSIfFPjshOecIsgogoxCO0TpgLox0KcNdGtil19t18AuRxvMarUDE2E3lVMBISFCBYSIiEhFeSKg2yPQ4HSYez5seBnqnwb1Tgm+VZqOJS4KBrSxH7edbicu78yCDXsgbR9s228/X78NtuTYAqHIB8WHbFpngCivHQIVF2ELjX5toVEiNK0DLepC6/pQOwh7GE7EGwW12sCGCdDzSQ1jkpCgAkJERKSyGp4OozbCzEGwaxZkrYEmoyC6ntvJAs8Y+8K/UaItKo7FcWzvgjG/9kSEq8ROkLUa9nwDDYa4nUak0kLk7RERERGXRde1S72e8jYU7YN1L9oXjI7P7WTuMMYOPwr34gFKhzFFwqZ33U4iEhAqIERERALFGGhxOYzcYF807poN68ZB3la3k4mbPFFQux1sfAP8xW6nEak0FRAiIiKBFtsQhq+EQR+BrwA2vALbPoWSPLeTiVsSO4MvD3bOcjuJSKWpgBAREakqTUbB6K3Q4Y+wfzGsfRb2LQDHf+LHSmhJaAOeaNisYUwS/FRAiIiIVKXIWtDzCRi21O5avWMqrHsJcta7nUyqkycCaneAtEngK3Q7jUilqIAQERGpDkmd7dyIgR+AvwjS3rI7FBfscjuZVJfETuAvhB0z3E4iUikqIERERKqLMdDsQrhgJ/R4wk6uXvcibP0IijLcTidVLb4VeGPspnIiQUwFhIiISHXzRkPHP8LobXZ+ROZyOz9ixzQoyXE7nVQVjxdqdYDN74GvyO00IhWmAkJERMQt0XXs/IiRG6HV7+wE6zVPwc6ZUJLrdjqpCrU72mFMu2a7nUSkwipVQBhjHjfGrDbGLDXGTDHGJAUol4iIBCm1DRUQ3xROHg/D10DzMbD3O/jlKdj5pQqJUJPQyu4LsXWK20lEKqyyPRAzgc6O43QFfgHuq3wkEREJcmobKqp2WzjlTThvJTS9CPZ+W9ojMQOKNbQpJHgi7JKuaRPBcdxOI1IhlSogHMeZ4ThOSemnPwBNKh9JRESCmdqGAEjsCAMmwnmroPklsPcH2yOxfSoUZbqdTiqrVls712X/YreTiFRIIOdAXAtMO9YXjTE3GmMWGmMWZvu0E6eISJgoc9uwZ8+eaowVJBI7wClvwYhfoOVVkL4Ifnkatn4MhfvcTicVVauNvd1+zP8aIjXaCQsIY8yXxpjlR/kYdcgx9wMlwNvHOo/jOOMcx+ntOE7vWt64wKQXERFXVEXbkJKSUh3Rg1OtNtDvFRi1EdrdWrpq03Ow5QPI3+l2OimviASIaQjrxrudRKRCIk50gOM4Zx7v68aYq4HhwBmOo8F8IiLhQG2DS+KbQe9noNNfYc3/YPX/IHMFJLSFlIH26xIcElrCvvlQkgcRemNVgktlV2EaCtwDjHQcR+OSREREbUN1iG0A3R+zG9J1fQTyt8HG12DDBMhZr8m5wSC+BTg+2Pej20lEyq2ycyCeBWoBM40xS4wxLwYgk4iIBDe1DdUlKhk63w8X7oae/4OidEh7Cza8DFlrVEjUZHFN7e2eb93NIVIBJxzCdDyO47QJVBAREQkNahtcEBEPHe6EtrfAxtdh8T2w+R2IaQApg6F2BzDG7ZRyKG8sRNezw5hEgox2ohYREQkV3mhocyNcuAf6vQ7+YtjyHqx/CbJWqUeipolJhV1z3E4hUm4qIEREREKNJwJaXQUX7IL+b4G/BDa/B+vHQ/ZaFRI1RWwDKMnWkrwSdFRAiIiIhCpPBLQcaydb95sAvnzYNBE2ToC8LW6nk+jSpYuzVrubQ6ScVECIiIiEOk8EtLra9kj0fs5Ott7wqu2VKEx3O134iqpjb7PXuptDpJxUQIiIiIQLb5TdiO6CndDlYchZB+uegx1fgK/A7XThJzLR3uZucjeHSDmpgBAREQk3EfHQ5QEYtRla/Q72/QC/PAv7l2h+RHXyRIA3HvJ3uJ1EpFxUQIiIiISr2IZw8sswdKHdU2Lbx3Z+RMEet5OFj4g4KNzrdgqRclEBISIiEu7q9ILzN9lionAPrH/RLi/q97mdLPR5oqE40+0UIuWiAkJERETAeKD1dTBqEzS/DPZ8DevHaXhNVfNEQNYvbqcQKRcVECIiIvKrmBQ45S0Y/Cn48mD9y7B7Ljh+t5OFKANo3okEFxUQIiIicqTGw21vRGJH2P0VbHwdijTUJuAcP3o5JsFGv7EiIiJydNF1YNgy6P8GFOyE9S9puE2gOcW2SBMJIiogRERE5NiMgZZXwrDldt+CzZNg55ca0hQoJXm/bignEiRUQIiIiMiJ1W4LozZC6xtg77eQ9rZ98SsV5zhQkmOX0xUJIiogREREpGy8MXDyODj5FcjbBOvHQ8Fut1MFr5JccEogvoXbSUTKRQWEiIiIlE/ra+Gsb+2L3w2vQLbmRVRIYemGfbXbu5tDpJxUQIiIiEj51TsZhq+CqLqw6R3Y96PbiYJPwU57m9TN3Rwi5aQCQkRERComrgmMXAu12sOO6bBjmiZXl0feFohMgtgGbicRKRcVECIiIlJxEfFw3nJofyfsmw+b3wd/sdupaj7HD7lp0GSk20lEyk0FhIiIiFSOxwu9/ge9noLs1bDxDa3QdCJ5W8GXD43OdTuJSLmpgBAREZHAaH87DPwACnbAhlehKMPtRDVX1kowXkhVASHBRwWEiIiIBE6zC+GMr+wSpRtegfydbieqefw+yFgGtdpBVKLbaUTKTQWEiIiIBFb9U2HoQjAe2Pga5GxwO1HNkrUSfHnQ879uJxGpEBUQIiIiEnhJneC8FXaVobS3Yf/PbieqGRwH9n5nl79tdLbbaUQqRAWEiIiIVI24JjBiDcQ3g20fwe659gV0OMtabfd/6Pkf20MjEoT0mysiIiJVJyoJhq+BFlfA7q9g2yd2DkA48pfAri8huh60GOt2GpEKi3A7gIiIiIQ4bxT0fwMSWsHyh6E4A5peDBFxbierXnu/haJ0GDIdPHoJJsFLPRAiIiJS9YyBrg9B/7fsDszrX4aC3W6nqj75O2HPXEjsDKnnuJ1GpFJUQIiIiEj1aTkWzpoHTrFd5jVzpduJqp6vELZ8AN44u8StSJBTASEiIiLVq14/GL4aouvDlvdhxwxwQnRehOPA1o/s0KXBn0JMPbcTiVSaCggRERGpfnGNYeQ6aHsr7PseNkyAov1upwosx4Ed0yF7NfR6EhoMcTuRSECogBARERF3eKOhz3Mw4F0o3APrXrL7RYTCUq+OA7tmQfp86PD/oP3tbicSCRgVECIiIuKu5pfYIU0xDex+EZvfheJst1NVnOOHHVPtqkttb4EeT7idSCSgVECIiIiI+xJawMgN0ONxyFkHa5+DfQvsi/Fg4iuATe9A+kLo+Cfo/ZxdgUokhKiAEBERkZrB44WOd8N5qyA21b6Lv/5lyN3kdrIy6dCgENaNg5z10OcF6PFvFQ8SklRAiIiISM1Suy2MXA+nTIKSXNg4ATZNgoJdbic7On8x1w/Yz3NjdgJ+OOsbaHuz26lEqowKCBEREal5jIEWY+CCHdDtn7YXYt2Ldn5E3la301mOA5nLYe3zXNE3iy/WNoRRaZByitvJRKqU9lEXERGRmisiDjr9xU5GXv0/WPk4ZK2G2CZQtw/UPgk81fxyxvFD1krYM8/2isQ04M7J3ViyI4lhDyVVbxYRF6gHQkRERGq+qGTo+jBcuBt6Pgm+PNg6BVY/Ads+gZwNVb8ZXWE67JoDa56CLR+CvwT6vwnnb2PJjqSqvbZIDaIeCBEREQkekbWgwx3Q/g/2xfyGCXZY0/7F4ImGhNYQ3wLim0N0PTCVeK/UXwz52yBnI2T/AgU77f0JraHnG9B4ROXOLxKkVECIiIhI8DEeaHi6/ej7EuycAds+hU3v2eFFAJ4ou7dEdD3bgxGZCN548MaAJxIwgGMLBV8+lORAcSYUpUPBbvuB3x4X18Tu59DsEohv6t7zFqkBVECIiIhIcIuIhSaj7Eff8XYZ1T3fQfoC2PopZP0CvtxynK8WxNSHk/4M9fpD/UEQlVRl8UWCTaUKCGPMP4BR2PJ8N3CN4zjbAxFMRESCk9oGcZUxUKuN/Wh1FfR+xt5fkmtXbyrYbXsZfPng99m9J7yxtocipgHENrYFiYgcU2V7IB53HOcBAGPM7cDfAC18LCIS3tQ2SM0TEQ+129sPEamUSs38cRwn65BP4wGncnFERCTYqW0QEQltxnEq93fdGPNP4CogEzjNcZw9xzjuRuDG0k87A8srdeHgUQ/Y63aIaqLnGrrC6fmG03Nt7zhOrao4sdqGEwqn3zM919Ck5xqaytQunLCAMMZ8CTQ8ypfudxzn40OOuw+IcRznwRNe1JiFjuP0PtFxoUDPNTSF03OF8Hq+eq5lfqzahkrQcw1Neq6hSc/1SCecA+E4zpllvOZE4HPghI2EiIgEN7UNIiLhq1JzIIwxbQ/5dCSwunJxREQk2KltEBEJbZVdhekxY0x77FJ9myj7KhvjKnndYKLnGprC6blCeD1fPdfKU9twYnquoUnPNTTpuf5GpSdRi4iIiIhI+KjUECYREREREQkvKiBERERERKTMXCsgjDH/MMYsNcYsMcbMMMakupWlqhljHjfGrC59vlOMMUluZ6oqxpiLjTErjDF+Y0xILnlmjBlqjFljjFlnjLnX7TxVxRjzqjFmtzEm5NflN8Y0NcZ8ZYxZVfr7e4fbmaqKMSbGGDPfGPNz6XN9yO1MB6hdCF2h3jaES7sAahvczlRVyts2uDYHwhhT+8BupcaY24GTHMcp60S7oGKMORuY7ThOiTHm/wAcx7nH5VhVwhjTETtx8iXgbsdxFrocKaCMMV7gF+AsYCuwALjMcZyVrgarAsaYQUAO8IbjOJ3dzlOVjDGNgEaO4/xkjKkFLALOD9GfqwHiHcfJMcZEAvOAOxzH+cHlaGoXQrRdgNBuG8KpXQC1DahtAFzsgTjQSJSKB0J2NrfjODMcxykp/fQHoImbeaqS4zirHMdZ43aOKtQXWOc4zgbHcYqAd4BRLmeqEo7jzAXS3c5RHRzH2eE4zk+l/84GVgGN3U1VNRwrp/TTyNKPGvH3V+1C6ArxtiFs2gVQ24DaBsDlORDGmH8aY7YAY4G/uZmlGl0LTHM7hFRYY2DLIZ9vJUT/mIQrY0wLoAfwo8tRqowxxmuMWQLsBmY6jlNjnqvaBQlCahfCgNqGw1VpAWGM+dIYs/woH6MAHMe533GcpsDbwG1VmaWqnei5lh5zP1CCfb5BqyzPNYSZo9wXsu+ShhtjTALwIXDnb94NDymO4/gcx+mOfde7rzGm2oYhqF0IzXYBwrptULsQ4tQ2HKmyG8mdKMiZZTx0IvA58GAVxqlSJ3quxpirgeHAGU6Qb75Rjp9rKNoKND3k8ybAdpeySACVjvn8EHjbcZzJbuepDo7jZBhj5gBDgWqZEKl24Veh1C5AWLcNahdCmNqGo7cNbq7C1PaQT0cCq93KUtWMMUOBe4CRjuPkuZ1HKmUB0NYY09IYEwWMAT5xOZNUUunksVeAVY7j/NftPFXJGJNyYMUfY0wscCY15O+v2gUJUmoXQpTahmP/DXZzFaYPgfbYVRk2ATc7jrPNlTBVzBizDogG9pXe9UMIrywyGngGSAEygCWO45zjaqgAM8YMA54EvMCrjuP8091EVcMYMwkYAtQDdgEPOo7ziquhqogxZiDwDbAM+zcJ4C+O40x1L1XVMMZ0BV7H/v56gPccx3nY3VSW2oXQbBcg9NuGcGkXQG0Dahvs8SHQayoiIiIiItVEO1GLiIiIiEiZqYAQEREREZEyUwEhIiIiIiJlpgJCRERERETKTAWEiIiIiIiUmQoIEREREREpMxUQIiIiIiJSZv8fnC2eUYgKWoEAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plot_all_stability()" + ] + }, + { + "cell_type": "markdown", + "id": "ad7d65e7", + "metadata": {}, + "source": [ + "in the above plot you can see the regions of stability for implicit and explicit methods.\n", + "Please be aware that implicit Euler is also stable in the left half plane.\n", + "The plot clearly shows why implicit" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.9.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/pySDC/projects/Resilience/accuracy_check.py b/pySDC/projects/Resilience/accuracy_check.py index 8900c8cf45..dbd8a22c0e 100644 --- a/pySDC/projects/Resilience/accuracy_check.py +++ b/pySDC/projects/Resilience/accuracy_check.py @@ -12,6 +12,10 @@ from pySDC.projects.Resilience.piline import run_piline +class do_nothing(hooks): + pass + + class log_errors(hooks): def post_step(self, step, level_number): @@ -44,28 +48,41 @@ def setup_mpl(font_size=8): mpl.rcParams.update(style_options) -def get_results_from_stats(stats, var, val): - e_extrapolated = np.array(get_sorted(stats, type='e_extrapolated'))[:, 1] - - e_loc = np.array(get_sorted(stats, type='e_loc'))[:, 1] - +def get_results_from_stats(stats, var, val, hook_class=log_errors): results = { - 'e_embedded': get_sorted(stats, type='e_embedded')[-1][1], - 'e_extrapolated': e_extrapolated[e_extrapolated != [None]][-1], - 'e': max([e_loc[-1], np.finfo(float).eps]), + 'e_embedded': 0., + 'e_extrapolated': 0., + 'e': 0., var: val, } + if hook_class == log_errors: + e_extrapolated = np.array(get_sorted(stats, type='e_extrapolated'))[:, 1] + e_embedded = np.array(get_sorted(stats, type='e_embedded'))[:, 1] + e_loc = np.array(get_sorted(stats, type='e_loc'))[:, 1] + + if len(e_extrapolated[e_extrapolated != [None]]) > 0: + results['e_extrapolated'] = e_extrapolated[e_extrapolated != [None]][-1] + + if len(e_loc[e_loc != [None]]) > 0: + results['e'] = max([e_loc[e_loc != [None]][-1], np.finfo(float).eps]) + + if len(e_embedded[e_embedded != [None]]) > 0: + results['e_embedded'] = e_embedded[e_embedded != [None]][-1] + return results -def multiple_runs(ax, k=5, serial=True, Tend_fixed=None): +def multiple_runs(k=5, serial=True, Tend_fixed=None, custom_description=None, prob=run_piline, dt_list=None, + hook_class=log_errors, custom_controller_params=None): """ A simple test program to compute the order of accuracy in time """ # assemble list of dt - if Tend_fixed: + if dt_list is not None: + pass + elif Tend_fixed: dt_list = 0.1 * 10.**-(np.arange(5) / 2) else: dt_list = 0.01 * 10.**-(np.arange(20) / 10.) @@ -82,11 +99,19 @@ def multiple_runs(ax, k=5, serial=True, Tend_fixed=None): EstimateExtrapolationErrorNonMPI: {'no_storage': not serial}, } } + if custom_description is not None: + desc = {**desc, **custom_description} Tend = Tend_fixed if Tend_fixed else 30 * dt_list[i] - stats, _, _ = run_piline(custom_description=desc, num_procs=num_procs, Tend=Tend, - hook_class=log_errors) + stats, controller, _ = prob(custom_description=desc, num_procs=num_procs, Tend=Tend, + hook_class=hook_class, custom_controller_params=custom_controller_params) + + level = controller.MS[-1].levels[-1] + e_glob = abs(level.prob.u_exact(t=level.time + level.dt) - level.u[-1]) + e_loc = abs(level.prob.u_exact(t=level.time + level.dt, u_init=level.u[0], t_init=level.time) - level.u[-1]) - res_ = get_results_from_stats(stats, 'dt', dt_list[i]) + res_ = get_results_from_stats(stats, 'dt', dt_list[i], hook_class) + res_['e_glob'] = e_glob + res_['e_loc'] = e_loc if i == 0: res = res_.copy() @@ -95,9 +120,19 @@ def multiple_runs(ax, k=5, serial=True, Tend_fixed=None): else: for key in res_.keys(): res[key].append(res_[key]) + return res - # visualize results - plot(res, ax, k) + +def plot_order(res, ax, k): + color = plt.rcParams['axes.prop_cycle'].by_key()['color'][k - 2] + + key = 'e_loc' + order = get_accuracy_order(res, key=key, thresh=1e-11) + label = f'k={k}, p={np.mean(order):.2f}' + ax.loglog(res['dt'], res[key], color=color, ls='-', label=label) + ax.set_xlabel(r'$\Delta t$') + ax.set_ylabel(r'$\epsilon$') + ax.legend(frameon=False, loc='lower right') def plot(res, ax, k): @@ -106,7 +141,7 @@ def plot(res, ax, k): color = plt.rcParams['axes.prop_cycle'].by_key()['color'][k - 2] for i in range(len(keys)): - order = get_accuracy_order(res, key=keys[i], order=k) + order = get_accuracy_order(res, key=keys[i]) if keys[i] == 'e_embedded': label = rf'$k={{{np.mean(order):.2f}}}$' assert np.isclose(np.mean(order), k, atol=3e-1), f'Expected embedded error estimate to have order {k} \ @@ -123,12 +158,13 @@ def plot(res, ax, k): ax.legend(frameon=False, loc='lower right') -def get_accuracy_order(results, key='e_embedded', order=5): +def get_accuracy_order(results, key='e_embedded', thresh=1e-14): """ Routine to compute the order of accuracy in time Args: - results: the dictionary containing the errors + results (dict): the dictionary containing the errors + key (str): The key in the dictionary correspdoning to a specific error Returns: the list of orders @@ -144,7 +180,7 @@ def get_accuracy_order(results, key='e_embedded', order=5): # compute order as log(prev_error/this_error)/log(this_dt/old_dt) <-- depends on the sorting of the list! try: tmp = np.log(results[key][i] / results[key][i - 1]) / np.log(dt_list[i] / dt_list[i - 1]) - if results[key][i] > 1e-14 and results[key][i - 1] > 1e-14: + if results[key][i] > thresh and results[key][i - 1] > thresh: order.append(tmp) except TypeError: print('Type Warning', results[key]) @@ -152,10 +188,25 @@ def get_accuracy_order(results, key='e_embedded', order=5): return order -def plot_all_errors(ax, ks, serial, Tend_fixed=None): +def plot_orders(ax, ks, serial, Tend_fixed=None, custom_description=None, prob=run_piline, dt_list=None, + custom_controller_params=None): + for i in range(len(ks)): + k = ks[i] + res = multiple_runs(k=k, serial=serial, Tend_fixed=Tend_fixed, custom_description=custom_description, + prob=prob, dt_list=dt_list, hook_class=do_nothing, + custom_controller_params=custom_controller_params) + plot_order(res, ax, k) + + +def plot_all_errors(ax, ks, serial, Tend_fixed=None, custom_description=None, prob=run_piline): for i in range(len(ks)): k = ks[i] - multiple_runs(k=k, ax=ax, serial=serial, Tend_fixed=Tend_fixed) + res = multiple_runs(k=k, serial=serial, Tend_fixed=Tend_fixed, custom_description=custom_description, + prob=prob) + + # visualize results + plot(res, ax, k) + ax.plot([None, None], color='black', label=r'$\epsilon_\mathrm{embedded}$', ls='-') ax.plot([None, None], color='black', label=r'$\epsilon_\mathrm{extrapolated}$', ls=':') ax.plot([None, None], color='black', label=r'$e$', ls='-.') diff --git a/pySDC/projects/Resilience/dahlquist.py b/pySDC/projects/Resilience/dahlquist.py new file mode 100644 index 0000000000..a73fb417ee --- /dev/null +++ b/pySDC/projects/Resilience/dahlquist.py @@ -0,0 +1,143 @@ +# script to run a simple advection problem +from pySDC.implementations.collocation_classes.gauss_radau_right import CollGaussRadau_Right +from pySDC.implementations.problem_classes.TestEquation_0D import testequation0d +from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit +from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI +from pySDC.core.Hooks import hooks +from pySDC.helpers.stats_helper import get_sorted +import numpy as np +import matplotlib.pyplot as plt + + +class log_data(hooks): + + def post_iteration(self, step, level_number): + + super(log_data, self).post_iteration(step, level_number) + + # some abbreviations + L = step.levels[level_number] + + L.sweep.compute_end_point() + + self.add_to_stats(process=step.status.slot, time=L.time + L.dt, level=L.level_index, iter=step.status.iter, + sweep=L.status.sweep, type='u', value=L.uend) + self.add_to_stats(process=step.status.slot, time=L.time, level=L.level_index, iter=0, + sweep=L.status.sweep, type='dt', value=L.dt) + + def pre_run(self, step, level_number): + super(log_data, self).pre_run(step, level_number) + L = step.levels[level_number] + self.add_to_stats(process=0, time=0, level=0, iter=0, sweep=0, type='lambdas', value=L.prob.params.lambdas) + + +def run_dahlquist(custom_description=None, num_procs=1, Tend=1., hook_class=log_data, fault_stuff=None, + custom_controller_params=None, custom_problem_params=None): + + # initialize level parameters + level_params = dict() + level_params['dt'] = 1. + + # initialize sweeper parameters + sweeper_params = dict() + sweeper_params['collocation_class'] = CollGaussRadau_Right + sweeper_params['num_nodes'] = 3 + sweeper_params['QI'] = 'LMMpar' + + # build lambdas + re = np.linspace(-30, 30, 400) + im = np.linspace(-50, 50, 400) + lambdas = np.array([[complex(re[i], im[j]) for i in range(len(re))] for j in range(len(im))]).\ + reshape((len(re) * len(im))) + + problem_params = { + 'lambdas': lambdas, + 'u0': 1., + } + + if custom_problem_params is not None: + problem_params = {**problem_params, **custom_problem_params} + + # initialize step parameters + step_params = dict() + step_params['maxiter'] = 5 + + # initialize controller parameters + controller_params = dict() + controller_params['logger_level'] = 30 + controller_params['hook_class'] = hook_class + controller_params['mssdc_jac'] = False + + if custom_controller_params is not None: + controller_params = {**controller_params, **custom_controller_params} + + # fill description dictionary for easy step instantiation + description = dict() + description['problem_class'] = testequation0d # pass problem class + description['problem_params'] = problem_params # pass problem parameters + description['sweeper_class'] = generic_implicit # pass sweeper + description['sweeper_params'] = sweeper_params # pass sweeper parameters + description['level_params'] = level_params # pass level parameters + description['step_params'] = step_params + + if custom_description is not None: + for k in custom_description.keys(): + if k == 'sweeper_class': + description[k] = custom_description[k] + continue + description[k] = {**description.get(k, {}), **custom_description.get(k, {})} + + # set time parameters + t0 = 0.0 + + # instantiate controller + controller = controller_nonMPI(num_procs=num_procs, controller_params=controller_params, + description=description) + + # insert faults + if fault_stuff is not None: + raise NotImplementedError('No fault stuff here...') + + # get initial values on finest level + P = controller.MS[0].levels[0].prob + uinit = P.u_exact(t0) + + # call main function to get things done... + uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) + return stats, controller, Tend + + +def plot_stability(stats, ax=None, iter=None, colors=None, crosshair=True, fill=False): + lambdas = get_sorted(stats, type='lambdas')[0][1] + u = get_sorted(stats, type='u', sortby='iter') + + # decorate + if crosshair: + ax.axhline(0, color='black', alpha=1.) + ax.axvline(0, color='black', alpha=1.) + + if ax is None: + fig, ax = plt.subplots(1, 1) + + iter = [1] if iter is None else iter + colors = ['blue', 'red', 'violet', 'green'] if colors is None else colors + + for i in iter: + # isolate the solutions from the iteration you want + U = np.reshape([me[1] for me in u if me[0] == i], (len(np.unique(lambdas.real)), len(np.unique(lambdas.imag)))) + + # get a grid for plotting + X, Y = np.meshgrid(np.unique(lambdas.real), np.unique(lambdas.imag)) + if fill: + ax.contourf(X, Y, abs(U), levels=[-np.inf, 1 - np.finfo(float).eps], colors=colors[i - 1], alpha=0.5) + ax.contour(X, Y, abs(U), levels=[1], colors=colors[i - 1]) + ax.plot([None], [None], color=colors[i - 1], label=f'k={i}') + + ax.legend(frameon=False) + + +if __name__ == '__main__': + custom_description = None + stats, controller, Tend = run_dahlquist(custom_description=custom_description) + plot_stability(stats, iter=[1, 2, 3]) + plt.show() diff --git a/pySDC/projects/Resilience/test_Runge_Kutta_sweeper.py b/pySDC/projects/Resilience/test_Runge_Kutta_sweeper.py new file mode 100644 index 0000000000..75ed35f41b --- /dev/null +++ b/pySDC/projects/Resilience/test_Runge_Kutta_sweeper.py @@ -0,0 +1,117 @@ +import matplotlib.pyplot as plt +import numpy as np + +from pySDC.projects.Resilience.accuracy_check import plot_orders +from pySDC.projects.Resilience.dahlquist import run_dahlquist, plot_stability + +from pySDC.projects.Resilience.advection import run_advection +from pySDC.projects.Resilience.vdp import run_vdp + +from pySDC.implementations.sweeper_classes.Runge_Kutta import RK1, RK4, MidpointMethod, CrankNicholson + + +colors = { + RK1: 'blue', + MidpointMethod: 'red', + RK4: 'orange', + CrankNicholson: 'purple', +} + + +def plot_order(sweeper, prob, dt_list, description=None, ax=None, Tend_fixed=None): + if ax is None: + fig, ax = plt.subplots(1, 1) + + description = dict() if description is None else description + description['sweeper_class'] = sweeper + description['sweeper_params'] = {'implicit': True} + + custom_controller_params = {'logger_level': 40} + + # determine the order + plot_orders(ax, [1], True, Tend_fixed=Tend_fixed, custom_description=description, dt_list=dt_list, prob=prob, + custom_controller_params=custom_controller_params) + + # check if we got the expected order for the local error + orders = { + RK1: 2, + MidpointMethod: 3, + RK4: 5, + CrankNicholson: 3, + } + numerical_order = float(ax.get_lines()[-1].get_label()[7:]) + expected_order = orders.get(sweeper, numerical_order) + assert np.isclose(numerical_order, expected_order, atol=2.5e-1),\ + f"Expected order {expected_order}, got {numerical_order}!" + + # decorate + ax.get_lines()[-1].set_color(colors.get(sweeper, 'black')) + + label = f'{sweeper.__name__} - {ax.get_lines()[-1].get_label()[5:]}' + ax.get_lines()[-1].set_label(label) + ax.legend(frameon=False) + + +def plot_stability_single(sweeper, ax=None, description=None, implicit=True, re=None, im=None, crosshair=True): + if ax is None: + fig, ax = plt.subplots(1, 1) + + description = dict() if description is None else description + description['sweeper_class'] = sweeper + description['sweeper_params'] = {'implicit': implicit} + + custom_controller_params = {'logger_level': 40} + + re = np.linspace(-30, 30, 400) if re is None else re + im = np.linspace(-50, 50, 400) if im is None else im + lambdas = np.array([[complex(re[i], im[j]) for i in range(len(re))] for j in range(len(im))]).\ + reshape((len(re) * len(im))) + custom_problem_params = {'lambdas': lambdas} + + stats, _, _ = run_dahlquist(custom_description=description, custom_problem_params=custom_problem_params, + custom_controller_params=custom_controller_params) + plot_stability(stats, ax=ax, iter=[1], colors=[colors.get(sweeper, 'black')], crosshair=crosshair, fill=True) + + ax.get_lines()[-1].set_label(sweeper.__name__) + ax.legend(frameon=False) + + +def plot_all_stability(): + fig, axs = plt.subplots(1, 2, figsize=(11, 5)) + + impl = [True, False] + sweepers = [[RK1, MidpointMethod, CrankNicholson], [RK1, MidpointMethod, RK4]] + titles = ['implicit', 'explicit'] + re = np.linspace(-3, 3, 400) + im = np.linspace(-3, 3, 400) + crosshair = [True, False, False] + + for j in range(len(impl)): + for i in range(len(sweepers[j])): + plot_stability_single(sweepers[j][i], implicit=impl[j], ax=axs[j], re=re, im=im, + crosshair=crosshair[i]) + axs[j].set_title(titles[j]) + + fig.tight_layout() + + +def plot_all_orders(prob, dt_list, Tend, sweepers): + fig, ax = plt.subplots(1, 1) + for i in range(len(sweepers)): + plot_order(sweepers[i], prob, dt_list, Tend_fixed=Tend, ax=ax) + + +def test_vdp(): + Tend = 5e-2 + plot_all_orders(run_vdp, Tend * 2.**(-np.arange(8)), Tend, [RK1, MidpointMethod, CrankNicholson, RK4]) + + +def test_advection(): + plot_all_orders(run_advection, 1.e-3 * 2.**(-np.arange(8)), None, [RK1, MidpointMethod, CrankNicholson]) + + +if __name__ == '__main__': + test_vdp() + test_advection() + plot_all_stability() + plt.show() diff --git a/pySDC/tests/test_projects/test_resilience/test_rk.py b/pySDC/tests/test_projects/test_resilience/test_rk.py new file mode 100644 index 0000000000..4ad7bb09cc --- /dev/null +++ b/pySDC/tests/test_projects/test_resilience/test_rk.py @@ -0,0 +1,7 @@ +import pytest + +from pySDC.projects.Resilience.test_Runge_Kutta_sweeper import test_vdp, test_advection + +def test_main(): + test_vdp() + test_advection()