{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Tutorial 1 : Thermomechanical Model\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# A quick Lithospheric Model with visco-plastic non-linear rheologies\n", "\n", "**Romain Beucher**\n", "romain.beucher (at) anu.edu.au\n", "\n", "The following tutorial presents a simple usage of the geodynamics module.\n", "The *geodynamics* module intents to facilitate rapid prototyping of geodynamics models. It can be seen as a set of high-level functions within the underworld ecosystem. It is a means to quickly get the user into Underworld modelling and assumes very little knowledge in coding. The module make some assumptions based on how the user defines the boundary conditions and the properties of the materials (rocks, phases). \n", "\n", "Its simplicity comes with a relatively more rigid workflow (compared to the classic Underworld functions).\n", "However, the user can easily break the high level objects and get back to core Underworld function at any step of model design. As we think the low-level interface is more flexible, and in so allows for more complex models, we strongly encourage users to explore and break the High Level functions. We hope that the user will naturally move to the low-level functionalities as he or her gets more confident, and by doing so will access the wide range of possibilities offered by Underworld.\n", "\n", "![Tutorial 1](./images/Tutorial1.gif)\n", "\n", "\n", "The module can be imported as follows:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "loaded rc file /home/romain/Projects/Project_UWGeodynamics/UWGeodynamics/UWGeodynamics/uwgeo-data/uwgeodynamicsrc\n" ] } ], "source": [ "import UWGeodynamics as GEO\n", "from UWGeodynamics import visualisation as vis" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Working with units\n", "\n", "Note that this is not an obligation and you can use values without units \n", "\n", "\n", "The geodynamics module enables usage of a broad range of units using a *UnitRegistry*. You can attach a unit to any value that you need to define. A value with a units attached becomes a *Quantity* python object. The geodynamics module take care of the conversion internally so you may use any units that you think are appropriate. You can also mix them.\n", "\n", "The module will also try to work out a scaling of the values to help the computation process. The user can chose to alter the way his or her values are scaled or can rely on the default options.\n", "\n", "To use the units system, you can link the unit registry as follow:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "u = GEO.UnitRegistry" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Scaling" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "half_rate = 1.8 * u.centimeter / u.year\n", "model_length = 360e3 * u.meter\n", "surfaceTemp = 273.15 * u.degK\n", "baseModelTemp = 1603.15 * u.degK\n", "bodyforce = 3300 * u.kilogram / u.metre**3 * 9.81 * u.meter / u.second**2\n", "\n", "KL = model_length\n", "Kt = KL / half_rate\n", "KM = bodyforce * KL**2 * Kt**2\n", "KT = (baseModelTemp - surfaceTemp)\n", "\n", "GEO.scaling_coefficients[\"[length]\"] = KL\n", "GEO.scaling_coefficients[\"[time]\"] = Kt\n", "GEO.scaling_coefficients[\"[mass]\"]= KM\n", "GEO.scaling_coefficients[\"[temperature]\"] = KT" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Define the external geometry\n", "\n", "The first step is to define the geometry of our problem, essentially a box on which we will apply some physical constraints and that will contain a set of materials. We can think of it as an \"universe\".\n", "The \"laws\" and their effects are calculated on a mesh, that mesh discretized our universe into finite elements.\n", "\n", "The geodynamics module allows you to quickly set up a model by creating a *Model* object.\n", "A series of arguments are required to define a *Model*:\n", "\n", " - The number of elements in each direction elementRes=(nx, ny);\n", " - The minimum coordinates for each axis minCoord=(minX, minY)\n", " - The maximum coordinates for each axis maxCoord=(maxX, maxY)\n", " - A vector that defines the magnitude and direction of the gravity components gravity=(gx, gy)\n", " " ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "Model = GEO.Model(elementRes=(128, 64), \n", " minCoord=(0. * u.kilometer, -110. * u.kilometer), \n", " maxCoord=(360. * u.kilometer, 10. * u.kilometer), \n", " gravity=(0.0, -9.81 * u.meter / u.second**2))" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "Model.outputDir = \"outputs_tutorial1\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Add some Materials\n", "\n", "Now that we have our \"universe\" (box, sand pit) ready, we need to fill it with some materials.\n", "The *geodynamics* module is designed around that idea of materials, which are essentially a way to define physical properties across the Model domain.\n", "\n", "A material (or a phase) is first defined by the space it takes in the box (its shape)." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "Model.diffusivity = 1e-6 * u.metre**2 / u.second \n", "Model.capacity = 1000. * u.joule / (u.kelvin * u.kilogram)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Model we are building is essentially a layered cake. The geodynamics module provide and easy way to define a layer by defining shape as *layer* and specifying its *top* and *bottom*. The order is important: when 2 shapes overlay each other, only the second is used." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "air = Model.add_material(name=\"Air\", shape=GEO.shapes.Layer(top=Model.top, bottom=2.0 * u.kilometer))\n", "stickyAir = Model.add_material(name=\"StickyAir\", shape=GEO.shapes.Layer(top=air.bottom, bottom= 0.0 * u.kilometer))\n", "crust = Model.add_material(name=\"Crust\", shape=GEO.shapes.Layer(top=stickyAir.bottom, bottom=-35.0 * u.kilometer))\n", "mantleLithosphere = Model.add_material(name=\"MantleLithosphere\", shape=GEO.shapes.Layer(top=crust.bottom, bottom=-100.0 * u.kilometer))\n", "mantle = Model.add_material(name=\"Mantle\", shape=GEO.shapes.Layer(top=mantleLithosphere.bottom, bottom=Model.bottom))" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "air.diffusivity = 1.0e-6 * u.metre**2 / u.second\n", "stickyAir.diffusivity = 1.0e-6 * u.metre**2 / u.second\n", "\n", "air.capacity = 100. * u.joule / (u.kelvin * u.kilogram)\n", "stickyAir.capacity = 100. * u.joule / (u.kelvin * u.kilogram)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "air.density = 1. * u.kilogram / u.metre**3\n", "stickyAir.density = 1. * u.kilogram / u.metre**3\n", "crust.density = GEO.LinearDensity(2620. * u.kilogram / u.metre**3, thermalExpansivity=3e-5 / u.kelvin)\n", "mantleLithosphere.density = GEO.LinearDensity(3370. * u.kilogram / u.metre**3, thermalExpansivity=3e-5 / u.kelvin)\n", "mantle.density = GEO.LinearDensity(3370. * u.kilogram / u.metre**3, thermalExpansivity=3e-5 / u.kelvin)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "crust.radiogenicHeatProd = 0.7 * u.microwatt / u.meter**3\n", "mantleLithosphere.radiogenicHeatProd = 0.02e-6 * u.microwatt / u.meter**3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Define Viscosities\n", "\n", "The rheology library contains some commonly used rheologies stored in a python dictionary structure. We can list the keys defining the rheologies as follows:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "rh = GEO.ViscousCreepRegistry()" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "Model.minViscosity = 1e18 * u.pascal * u.second\n", "Model.maxViscosity = 1e23 * u.pascal * u.second\n", "\n", "air.viscosity = 1e19 * u.pascal * u.second\n", "stickyAir.viscosity = 1e20 * u.pascal * u.second\n", "crust.viscosity = rh.Wet_Quartz_Dislocation_Gleason_and_Tullis_1995\n", "mantleLithosphere.viscosity = rh.Dry_Olivine_Dislocation_Karato_and_Wu_1993\n", "mantle.viscosity = 0.2 * rh.Dry_Olivine_Dislocation_Karato_and_Wu_1993" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Define Plasticity\n", "\n", "Plastic behavior is assigned using the same approach as for viscosities." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "pl = GEO.PlasticityRegistry()" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "crust.plasticity = pl.Huismans_et_al_2011_Crust\n", "mantleLithosphere.plasticity = pl.Huismans_et_al_2011_Crust\n", "mantle.plasticity = pl.Huismans_et_al_2011_Crust" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " \n", " \n", " \n", " \n", "
Air
DensityConstant (1.0 kilogram / meter ** 3)
Radiogenic Heat Production0.0 microwatt / meter ** 3
Diffusivity1e-06 meter ** 2 / second
Capacity100.0 joule / kelvin / kilogram
Min Viscosity Limit1e+18 pascal * second
Max Viscosity Limit1e+23 pascal * second
Rheology (Viscous)
ViscosityConstant (1e+19 pascal * second)
" ], "text/plain": [ "" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "air" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " \n", " \n", " \n", " \n", "
Crust
DensityLinear (ref: 2620.0 kilogram / meter ** 3)
Radiogenic Heat Production0.7 microwatt / meter ** 3
Diffusivity1e-06 meter ** 2 / second
Capacity1000.0 joule / kelvin / kilogram
Min Viscosity Limit1e+18 pascal * second
Max Viscosity Limit1e+23 pascal * second
Rheology (Visco-plastic)
ViscosityWet Quartz, Viscous Dislocation Creep, Gleason and Tullis, 1995
PlasticityHuismans et al. 2011, (Crust)
" ], "text/plain": [ "" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "crust" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Temperature Boundary Conditions" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Model.set_temperatureBCs(top=293.15 * u.degK, \n", " bottom=1603.15 * u.degK, \n", " materials=[(air, 293.15 * u.degK), (stickyAir, 293.15 * u.degK)])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Velocity Boundary Conditions" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Model.set_velocityBCs(left=[-2.5 * u.centimeter / u.year, None],\n", " right=[2.5 * u.centimeter / u.year, None],\n", " bottom=GEO.LecodeIsostasy(reference_mat=mantle, average=True))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Initialize plastic strain" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "def gaussian(xx, centre, width):\n", " return ( np.exp( -(xx - centre)**2 / width**2 ))\n", "\n", "maxDamage = 0.7\n", "Model.plasticStrain.data[:] = maxDamage * np.random.rand(*Model.plasticStrain.data.shape[:])\n", "Model.plasticStrain.data[:,0] *= gaussian(Model.swarm.particleCoordinates.data[:,0], (GEO.nd(Model.maxCoord[0] - Model.minCoord[0])) / 2.0, GEO.nd(5.0 * u.kilometer))\n", "Model.plasticStrain.data[:,0] *= gaussian(Model.swarm.particleCoordinates.data[:,1], GEO.nd(-35. * u.kilometer) , GEO.nd(5.0 * u.kilometer))" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "tags": [ "analysis" ] }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "Fig = vis.Figure(figsize=(1200, 400))\n", "Fig.Points(Model.swarm, Model.plasticStrain, fn_size=4.0)\n", "Fig.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Passive Tracers\n", "\n", "### Single tracer\n", "\n", "We create a single tracer that we place at the surface, in a central position. " ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "vertices = np.ndarray((1, 2))\n", "vertices[:, 0] = GEO.nd((Model.maxCoord[0] - Model.minCoord[0]) / 2.0)\n", "vertices[:, 1] = 0.\n", "\n", "Model.add_passive_tracers(name=\"Central\", vertices=vertices)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Interface tracers\n", "\n", "We add 2 sets of tracers for tracking the interface at the surface and the moho." ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "tags": [ "outputPrepend" ] }, "outputs": [], "source": [ "import numpy as np\n", "\n", "npoints = 1000 # This is the number of points used to define the surface\n", "coords = np.ndarray((npoints, 2))\n", "\n", "coords[:, 0] = np.linspace(GEO.nd(Model.minCoord[0]), GEO.nd(Model.maxCoord[0]), npoints)\n", "coords[:, 1] = 0.\n", "\n", "Model.add_passive_tracers(name=\"Surface\", vertices=coords)\n", "\n", "coords[:, 1] -= GEO.nd(35.*u.km)\n", "Model.add_passive_tracers(name=\"Moho\", vertices=coords)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Grid Tracers" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "pts = GEO.circles_grid(radius=2.0*u.kilometer, \n", " minCoord=[Model.minCoord[0], crust.bottom], \n", " maxCoord=[Model.maxCoord[0], 0.*u.kilometer])\n", "\n", "Model.add_passive_tracers(name=\"FSE_Crust\", vertices=pts)\n", "\n", "pts = GEO.circles_grid(radius=2.0*u.kilometer, \n", " minCoord=[Model.minCoord[0], mantleLithosphere.bottom], \n", " maxCoord=[Model.maxCoord[0], mantleLithosphere.top])\n", "\n", "Model.add_passive_tracers(name=\"FSE_Mantle\", vertices=pts)" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "tags": [ "analysis" ] }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "Fig = vis.Figure(figsize=(1200,400), title=\"Material Field\", quality=2)\n", "Fig.Points(Model.Central_tracers, pointSize=10, colour=\"blue\")\n", "Fig.Points(Model.Surface_tracers, pointSize=2.0)\n", "Fig.Points(Model.Moho_tracers, pointSize=2.0, colour=\"red\")\n", "Fig.Points(Model.FSE_Crust_tracers, pointSize=2.0)\n", "Fig.Points(Model.FSE_Mantle_tracers, pointSize=2.0)\n", "Fig.Points(Model.swarm, Model.materialField, fn_size=2.0)\n", "Fig.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Compute initial condition" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [], "source": [ "Model.init_model(temperature=\"steady-state\", pressure=\"lithostatic\")" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "tags": [ "analysis" ] }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "Fig = vis.Figure(figsize=(1200,400), title=\"Temperature Field (Celsius)\", quality=3)\n", "Fig.Surface(Model.mesh, GEO.dimensionalise(Model.temperature, u.degC), colours=\"coolwarm\")\n", "Fig.show()" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "tags": [ "analysis" ] }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "Fig = vis.Figure(figsize=(1200,400), title=\"Pressure Field (MPa)\", quality=3)\n", "Fig.Surface(Model.mesh, GEO.dimensionalise(Model.pressureField, u.megapascal))\n", "Fig.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Basic Analysis of the initial set-up\n", "\n", "Let's analyze the pressure and temperature field by creating a vertical profile at the center of the model." ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "tags": [ "analysis" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Average Temperature at Moho: 607 degree_Celsius\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/home/romain/Projects/Project_UWGeodynamics/virtualenv/lib/python3.9/site-packages/numpy/core/_asarray.py:102: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.\n", " return array(a, dtype, copy=False, order=order)\n", "/home/romain/Projects/Project_UWGeodynamics/virtualenv/lib/python3.9/site-packages/numpy/core/_asarray.py:102: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.\n", " return array(a, dtype, copy=False, order=order)\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3sAAAG5CAYAAAA3ci11AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABa5klEQVR4nO3dd3gVZfrG8ftJ7wkBQm9SRaRJt2HH3hV7x95XV1fddXfdn3Wt66rYOyr2rmsDQUB6kw5KJ7QkQCp5f3+cCQaEECDnzCnfz3Xlyjkzc07uYTRvnjMz72POOQEAAAAAokuc3wEAAAAAAHWPYg8AAAAAohDFHgAAAABEIYo9AAAAAIhCFHsAAAAAEIUo9gAAAAAgClHsAQgKM+toZpPNrMjMrjOzp83sLm/dQDNb4ndGAAD2FOMdwlmC3wGAumBmG6o9TZNUKmmz9/xy59zroU+1Z8xskaRLnXP/8zvLbrpV0nfOue5+BwEA7BpvDGqkwFi6UdLnkq5xzm2o6XUxivEOYYsze4gKzrmMqi9Jv0k6vtqysCv0zCzoH7QE82fU8r1bSZoRrAwAgKA73htXe0rqJenObTcIxXi2M4x3wI5R7CGqmVmcmd1mZvPNbI2ZvW1mud661mbmzOwiM1tsZuvM7Aoz621mU81svZn9p9p7XWhmo8zsP2ZWYGazzOywauuzzex5M1tuZkvN7B4zi9/mtY+Y2RpJd5tZWzP71su12sxeN7Mcb/tXJbWU9LGZbTCzW7d3KYiZLTKzw73Hd5vZcDN7zcwKJV1YU6bt/FtVvf4t71KUiWbWbZuf9Wczmyppo5klmNkJZjbD+7f63sz29rb9VtIhkv7j5e9gZi+Z2T07+NlNzexdM8s3s4Vmdt2uHmsAQHA455YqcGaviyR5Y+fVZjZX0lxv2XHepYzrzWy0mXWter03diz1xpbZVWOnmfUxs/FmVmhmK83sYW854x1QRyj2EO2ulXSSpIMlNZW0TtKT22zTV1J7SWdKelTSHZIOl7SPpDPM7OBttp0vqYGkv0l6z7ziUdJLkioktZPUQ9KRki7d5rULFLgs5l+STNK9Xq69JbWQdLckOefO09ZnKB+o5f6eKGm4pBxJr9ci0/Ze/46kXElvSPrAzBKrrT9L0rHe++8l6U1JN0hqKOkzBYrTJOfcoZJGKnDJT4Zzbs6OfqCZxUn6WNIUSc0kHSbpBjM7qpb7DAAIIjNrIekYSZOqLT5JgXGts5n1kPSCpMsl1Zf0jKSPzCzZzDpKukZSb+dcpqSjJC3y3uMxSY8557IktZX09i7EYrwDaoFiD9HuCkl3OOeWOOdKFSimTrOtL8v4p3OuxDn3lQL3JbzpnFvlfZI5UoFBo8oqSY8658qdc29Jmi3pWDNrpMBAeINzbqNzbpWkRyQNrvbaZc65J5xzFc65YufcPOfc1865UudcvqSHFShK98RPzrkPnHOVkrJqkWlbE5xzw51z5V6eFEn9qq1/3Dm32DlXrEBx/Km3D+WSHpKUKmnALmbuLamhc+4fzrky59wCSc/uJCcAIPg+MLP1kn6U9IOk/6u27l7n3FpvPBgi6Rnn3Fjn3Gbn3MsK3DvfT4F7/pIVKAoTnXOLnHPzvfcol9TOzBo45zY458bsQjbGO6AWfL/OGgiyVpLeN7PKass2K3B2rcrKao+Lt/M8o9rzpc45V+35rwqcmWslKVHScjOrWhcnaXG1bas/llcgPibpQEmZ3vbrarVXO1b9Z9Qm0w5f75yr9C6jabqD92+qwP5X336xAp9W7opWkpp6f1BUiVeg0AYA+OekGiYJ23a8ucDMrq22LElSU+fcD2Z2gwIftu5jZl9Kusk5t0zSJZL+IWmWmS2U9Hfn3Ce1zMZ4B9QCxR6i3WJJFzvnRm27wsxa78b7NTMzq1bwtZT0kfdzSiU1cM5V7OC1bpvn/+ct29c5t9bMTpL0nxq236jATKNV+eMVuJxkRz+jNpm21aLa+8dJai5p2Q7ef5mkfattb97rl9byZ1XPudA5134XXwcA8M+2482/nHP/2u6Gzr0h6Q0zy1LgEs/7JZ3nnJsr6SxvvDlF0nAzqy/GO6DOcBknot3Tkv5lZq0kycwamtmJe/B+eZKuM7NEMztdgXvtPnPOLZf0laR/m1mWBSaGabvN/X7bypS0QVKBmTWTdMs261cqcJ9AlTmSUszsWO++gjsVuDRmu3Yz035mdop3mesNCgyeO7qs5m0FLmE9zMtzs7f96Bref3vGSSryboZPNbN4M+tiZr138X0AAP54VtIVZtbXAtK9sSrTAj3oDjWzZEklClwxUylJZnaumTX0LsVc771XpRjvgDpDsYdo95gCZ96+MrMiBX6R992D9xurwGQuqxWYZOU059wab935Cly2MlOByzGHS2pSw3v9XYHprAskfSrpvW3W3yvpTm/mrz855wokXSXpOQU+TdwoaWeNWnc104cK3JuwTtJ5kk7x7k/4A+fcbEnnSnpCgX+P4xWYUKZsJ5m2fZ/Nko6T1F3SQu+9npOUvSvvAwDwh3NuvKTLFLg6ZZ2keZIu9FYnS7pPgd/tKxT40PR2b90gSTMs0Cv3MUmDvXvaGe+AOmJb334EYEfM7EIFmpwf4HeWYDCzuyW1c86d63cWAACChfEOsYQzewAAAAAQhcKq2DOzQRZotjnPzG7zOw8AAOGCMRIAsKvC5jJOb6alOZKOUOC67J8lneWcm+lrMAAAfMYYCQDYHeF0Zq+PpHnOuQXeDa/DJO3JrIkAAEQLxkgAwC4Lpz57zbR1A8sl2s6siWY2RNIQSUpPT9+vU6dOu/0D124s09L1xUpJjFfL3DQlJ4RT7QsAqG7ChAmrnXPb9tqKFSEfI+esLNoyPgIAwldN42M4FXu14pwbKmmoJPXq1cuNHz9+j97v65krdevwKSopr9Qdx3fWmb1bKNArEwAQTszsV78zhLu6HCMPf/gHdWyUqSfP6VlX8QAAQVDT+BhOp7KWSmpR7Xlzb1lQHdG5kb644SD1bJWj296bpqvfmKiCTdttswIAgF98GSMBAJEtnIq9nyW1N7M2ZpYkabACzbCDrlFWil69uK9uO7qTvpqxUkc/NkLjFq4NxY8GAKA2fBsjAQCRK2yKPedchaRrJH0p6RdJbzvnZoTq58fFma44uK3evXKAkhLiNHjoT3r4q9mq2FwZqggAAGyXX2Pk6PmrdenL43XXB9P11PfztXpDabB/JACgDoXVPXvOuc8kfeZnhm4tcvTJdQfq7o9m6PFv5+nHeav12OAeasEN6gAAH4V6jBxy4F76fPpyLVm3ST8vWquC4nKZSVcc3DZUEQAAeyisir1wkZGcoIdO76aDOjTUHe9N0zGPjdQ9J3fRid2b+R0NAICQOKN3C53RO3CbYFlFpTrc+bk2V4ZHb14AQO2EzWWc4eiEbk312fUHqn2jDF0/bLJufnuKNpRW+B0LAAAAAHaKYm8nWuSm6e3L++u6w9rr/UlLdOzjIzVl8Xq/YwEAAABAjSj2aiEhPk43HdFBw4b0V3lFpU59arSe+n6+KrmcBQAAAECYotjbBX3a5Orz6w/Skfs00v1fzNIFL45jZjIAAAAAYYkJWnZRdlqinjy7p94ct1h3fzxDxzw2Uo8N7qH+bev7HQ0AgKD6eMoyLVtfrPoZyWqQkaRGWSk6tFOeEuP57BgAwhHF3m4wM53dt6V6tMzR1W9M1DnPjdF1h7XXtYe2V3yc+R0PAIA6lRhvOrVnc01bul5fTF+htZvK5Lw7GZ4+dz8N6tLY34AAgO2i2NsDezfJ0sfXHKC7PpiuR/83V+MWrtWjZ3ZXXlaK39EAAKgzZqZ/n9Fty/PNlU5TlqzXKf8drdKKzT4mAwDUhOsu9lB6coL+fUY3PXBaV038bZ2OeXykRs7N9zsWAABBEx9nyklN9DsGAGAnKPbqgJnpjF4t9PE1Byg3PUnnvzBOD305WxWbK/2OBgAAACBGUezVofaNMvXh1QfojP1a6D/fzdPZz47VioISv2MBAAAAiEEUe3UsNSle95/WVY+e2V3TlxXouCdGatS81X7HAgAAABBjmKAlSE7q0UxdmmXpytcm6rznx+rGwzvo6kPaKY7ZOgEAUeSjycu0dH2x6qcnKTc90JKha/McZqcGgDBAsRdE7fIy9eE1++sv703Tv7+eowm/rdMjZ3RXvfQkv6MBALBH8rJS1KlxpkbNX61vZq3aat3fju+si/Zv41MyAEAVir0gS0tK0CNndlfvNrn6+0czdezjI/XkOT3Vo2U9v6MBALDbMpIT9MUNB0mSNpVVaM2GMq3eUKqT/ztaRSUVPqcDAEjcsxcSZqZz+rbSu1cOUFyc6YxnftJLoxbKVXWkBQAggqUlJahFbpq6Ns/xOwoAoBqKvRDat3m2Pr32QB3coaHu/nimrnljkopKyv2OBQAAACAKUeyFWHZaooae10u3Hd1JX8xYoRP/M0qzVxT5HQsAAABAlKHY80FcnOmKg9vqjUv7qqi0Qic9OUofT1nmdywAAAAAUYRiz0d996qvT689QPs0zdK1b07SPZ/MVMXmSr9jAQCwRyYvXq8PJi3ViDn5mr60QMsLirlPHQB8wGycPsvLStEbl/XT/332i577caGmLyvQf87uqQYZyX5HAwBgl8SZ1Lp+mr6dtUrfbtOO4aqBbXXroE4+JQOA2ESxFwaSEuJ09wn7qFuLbN3+3jQd9/iPeupc2jMAACKLmenbmweqsKRcazaWae3GMq3ZUKa/vD9NKwpK/I4HADGHyzjDyMk9muvdKwcoMcF05jNj9PrYX7nsBQAQUeLiTDlpSWrbMEO9W+dqUJfGSk+O9zsWAMQkir0ws0/TbH18zQHq37a+7nh/uv787lSVlG/2OxYAAACACEOxF4Zy0pL0woW9dd2h7fT2+CU6/emftGTdJr9jAQAAAIggFHthKj7OdNORHfXc+b20aPVGnfCfURq7YI3fsQAAAABECIq9MHd450b64Jr9lZOWqHOeG6tXx3AfHwAg8iwrKNaYBWs0d2WR1mwo1eZKxjIACDZm44wAbRtm6IOr99cNwybrrg+ma+ayAv39hC5KSqBWBwCEv/rpyRqzYK0GDx2zZZmZdF6/VvrHiV18TAYA0Y1iL0JkpSTq2fN76eGvZ+vJ7+ZrzsoNeurcnsrLTPE7GgAANXrjsr5atHpToBXDxlKt3Viml0Yv0vSlBX5HA4CoRrEXQeLjTLcc1Ul7N8nSLe9M1QlPjNLQ8/dT1+Y5fkcDAGCH0pIS1Llp1lbLvp21ShtLK3xKBACxgesAI9BxXZtq+JX9FR9nOu3pn/TexCV+RwIAAAAQZij2ItQ+TbP10TX7q2fLHN309hT969OZ3OwOAAAAYAuKvQhWPyNZr17SV+f3b6VnRy7UZa+MV1FJud+xAAAAAIQBir0Ilxgfp3+c2EX/PKmLfpiTr1OfGq3Fa2nADgAIf5vKNuu3NZu0obSCtkIAEARM0BIlzuvXSm3qp+uq1yfoxCdH6elz91OfNrl+xwIAYLvSkuI1a0WRDnrwO0lSUkKc6qcnaWDHPN17yr4+pwOA6MCZvShyQPsG+uDq/ZWTmqhznhujt8cv9jsSAADbdd8pXfXChb304GlddfvRnXTRgNZKT07Q1zNX+B0NAKIGZ/aizF4NM/T+Vfvr6jcm6tbhUzVv1Qb9eVAnxceZ39EAANiiXnqSDu3UaKtlG0or9OUMij0AqCuc2YtC2WmJevGi3jq/fysNHbFAQ14Zrw30MgIAAABiCsVelNoyccuJ++j7Ofk69b+jtWQdE7cAAAAAsYJiL8qd17+1Xr6oj5YVFOukJ0dryuL1fkcCAAAAEAIUezHggPYN9N6VA5SSGKczh/7E/RAAgLBV6aSyikq/YwBAVGCClhjRvlGm3r9qf136ynhd8doE3XHM3rrkgDYyY+IWAEB4SIyP09qNZepw5+fKTE5QbkaSctOT1KVptv5x4j6MWQCwiyj2YkjDzGQNu6yfbnp7su759Bf9tnaT/npcZyXEc4IXAOC/qw5pq3Z5GVq7sWzL1/SlBXp1zK/6yzF7KzUp3u+IABBRKPZiTGpSvJ48u6fu+2KWho5YoCXrivXEWT2Unsx/CgAAf+Vlpujcfq22Wvb0D/N13+ezfEoEAJGNUzoxKC7O9Jdj9tY9J3XRD3PydfrTP2lFQYnfsQAAAADUIYq9GHZuv1Z6/oJe+nXNRp305CjNXFbodyQAAAAAdYRiL8YN7Jind64YIEk645mfNGreap8TAQAAAKgLFHtQ56ZZeu+qAWqak6ILXxynDycv9TsSAAAAgD3ErByQJDXNSdU7VwzQkFfG6/phk7WysESXHbgX01wDAHwV741DBz7wneqnB1ox5GYkqVVumm46ogMzSgNADSj2sEV2aqJeuaSPbnp7iv7vs1latr5Edx3XWfFxFHwAAH+c2L2pikrKlb+hVGs2BNoxTFi0Tp9OXa5TejZTu7xMvyMCQNii2MNWkhPi9cTgHmqUmaIXRi3UqqISPXxGd6Uk0tsIABB6eVkpuunIjlst+2TqMl3zxiQ551MoAIgQFHv4g7g401+P76ymOSm659NftHrDOD17Xi9lpyX6HQ0AAABALXGhO3bo0gP30uNn9dCk39bptKdHa9n6Yr8jAQAAAKglij3U6IRuTfXyxX20oqBEpz41WvNWFfkdCQAAAEAtUOxhpwa0baBhl/dT+eZKnf70T5qyeL3fkQAAAADsBPfsoVb2aZqt4VcM0HkvjNVZz47R0PN66YD2DfyOBQCIYde8MUmNs1O2tGRokpOqc/q2ZFIxAPBwZg+11rpBut69YoBa5qbp4pd+1mfTlvsdCQAQg/rtVV9n9mqhxtkpWrepTGMXrtVrY3/VPz+ZqZ8XrfU7HgCEDc7sYZfkZaXorSH9dfHLP+vqNybqXyftq7P7tvQ7FgAghjTISNb9p3Xdatmk39bp5P+OVkUl/RgAoApn9rDLstMS9dolfXVwh4b6y/vT9OR38+RodgQAAACEFYo97JbUpHg9e34vndS9qR78crbu+fQXVfJpKgAAABA2uIwTuy0xPk4Pn9FdOWlJev7HhSoqKde9p3RVfJz5HQ0AAACIeRR72CNxcaa/Hd9ZWamJevybudpUtlmPnNldifGcNAYAAAD8RLGHPWZmuumIDkpPite9n89SSflm/efsnkx9DQAIuSe/naevZqxQbnqSctOTlZeZrCM6N2JMAhCTKPZQZy4/uK3SkuJ114czdOnL4zX0/P2UlsR/YgCA4Gubl6HDOuVp8bpN+nrmSq3bVK7N3r3kD53eTaft19znhAAQevwljjp1Xv/WSk1K0K3Dp+j858fphYt6Kysl0e9YAIAol5WSqOcv7L3leWWl04LVG3T4wyNUWrHZx2QA4B9urEKdO22/5nrirJ6avHi9znl2rNZtLPM7EgAgxsTFGR82Aoh5FHsIimO7NtHQ8/fT7JVFOnPoT1pVVOJ3JAAAACCmhLzYM7MWZvadmc00sxlmdr23PNfMvjazud73eqHOhrp1aKdGeunC3lq8tlhnDR2jVYUUfABQE8ZIAEBd8uPMXoWkm51znSX1k3S1mXWWdJukb5xz7SV94z1HhBvQroFevriPlheUaPDQMVpRQMEHADVgjAQA1JmQF3vOueXOuYne4yJJv0hqJulESS97m70s6aRQZ0Nw9GmTq1cu7qNVRaUaPPQnLS8o9jsSAIQlxsjg+Gzacj3+zVy9NuZXfT5tuX5etFblmyv9jgUAQefrbJxm1lpSD0ljJTVyzi33Vq2Q1GgHrxkiaYgktWzZMgQpURd6tc7Vyxf30YUvjNOZz4zRm0P6qVlOqt+xACBsMUbuuey0RPVomaOpSwo0at6ardbddnQnXXFwW5+SAUBo+DZBi5llSHpX0g3OucLq65xzTpLb3uucc0Odc72cc70aNmwYgqSoK/u1qqdXL+2rdZvKdOYzP2nx2k1+RwKAsMQYWTeSE+L1/lX7a9rdR2nOPUdr7F8O02fXHag4kzaUVPgdDwCCzpdiz8wSFRjEXnfOvectXmlmTbz1TSSt8iMbgqt7ixy9fmlfFRaXa/DQMRR8ALANxsjgSEqIU6OsFHVumiUz8zsOAISEH7NxmqTnJf3inHu42qqPJF3gPb5A0oehzobQ6No8R29c1k8bSit05jM/6dc1G/2OBABhgTESAFCX/Dizt7+k8yQdamaTva9jJN0n6QgzmyvpcO85olSXZtl647K+Ki7frLOfHcsZPgAIYIwEANSZkE/Q4pz7UdKOrp84LJRZ4K99mmbr1Uv66uxnx+js58borSH91ZRJWwDEMMZIAEBd8m2CFkAKnOF79ZK+Wr+xXOc8N1YrabwOAAiBKUvWa/iEJfp21kpNXrxei9duUmDuGwCIHr62XgAkqVuLHL10cR+d//xYnf3sGA0b0l8NM5P9jgUAiFLt8zI0cu5qjZy7eqvlVw1sq1sHdfIpFQDUPYo9hIX9WtXTCxf21oUv/qxznhujNy/rp/oZFHwAgLr3+fUHakNphdZuLNOajWVau6FMtwyfopWFpX5HA4A6xWWcCBt996qv5y/opV/XbNK5z4/T+k1lfkcCAEQhM1NmSqJa1U9Xz5b1dHjnRkpL4vNvANGHYg9hZUC7Bnr2/F6av2qDznt+nAqKy/2OBAAAAEQkij2EnYM6NNTT5/XUrBWFuuSln7WprMLvSAAAAEDEodhDWDq0UyM9NriHJv62Tpe/OkGlFZv9jgQAAABEFIo9hK1j9m2i+07pqpFzV+uGYZNVsbnS70gAgCi2dP0mjZq3WrNWFGpVUQnjDoCIx93ICGtn9G6hotIK/fOTmbrtvWl64NSuiovbUb9hAAB2T4PMZI1ZsFZjFozdavmlB7TRncd19ikVAOwZij2EvUsOaKPC4nI99s1cZaYk6K/HdZYZBR8AoO68NaSfFq/dpDUby7RmQ5nWbizV0z8s0KwVRX5HA4DdRrGHiHDD4e1VVFKhF0YtVHZqom44vIPfkQAAUSQlMV7tG2WqfbVlH0xe5lseAKgLFHuICGamO4/dW0Ul5Xr0f3OVmZKoSw5o43csAAAAIGxR7CFixMWZ7j1lX23w7uGrl5aoU3o29zsWAAAAEJaYjRMRJSE+To8O7q4Bbevr1uFT9f3sVX5HAgAAAMISxR4iTnJCvJ45bz91aJSpq16fqCmL1/sdCQAQpYpKK7Rw9UYVFJfLOed3HADYJVzGiYiUmZKoly7urVOfGq2LXvpZ7145QG0apPsdCwAQRdKTEzRiTr4Oeeh7SVJivKleWpJO7N5UdxxLOwYA4Y8ze4hYeZkpevmiPpKk818Yq1VFJT4nAgBEk0fP7K6XLuqth8/opjuP3VuXHriXUhLj9f3sfL+jAUCtcGYPEW2vhhl68cLeGjx0jC584We9dXk/ZaYk+h0LABAFctOTNLBj3lbLfl2zUXNXbvApEQDsGs7sIeJ1a5Gjp87tqTkri3TFaxNUWrHZ70gAAACA7yj2EBUGdszT/ad21ah5a3Tr8KncRA8AAICYx2WciBqn7tdcKwpL9OCXs9UqN003HdnR70gAAACAbyj2EFWuGthWv63ZpMe/nacWuWk6vVcLvyMBAKLM5kqnkvLNSkmM9zsKANSIYg9Rxcx0z8ldtHR9sW5/b5qa5aRqQLsGfscCAESJpPg4LVi9UZ3u+kLpSfHKzUhSbnqyDmhXX7cc1cnveACwFe7ZQ9RJjI/Tf8/tqb0apuvy1yZo7soivyMBAKLErYM66d5T9tUtR3XUmb1bar+W9bRmQ6le/elXv6MBwB9Q7CEqZaUk6oULeyslMV4XvfSz8otK/Y4EAIgCTXNSdVaflrr6kHb66/Gd9ejgHjp870Z+xwKA7aLYQ9RqXi9Nz1/QS2s2lOnSl39WcRktGQAAABA7KPYQ1bo2z9Fjg7tr6tIC3fT2ZFVW0pIBAAAAsYFiD1HvyH0a645j9tbn01fo0W/m+h0HAAAACAlm40RMuOSANpqzskiPfzNXHRpl6LiuTf2OBAAAAAQVxR5igpnpnyd10YL8jfrTO1PUKjdd+zbP9jsWACAKxMeZCksq1POfXys3PUm56Umqn56kfZtn66qB7fyOByCGcRknYkZyQryePm8/1U9P1mWvjNeqwhK/IwEAosD5/VvpukPbaVCXxmrXMEOS9POidXrgi9kqrWByMAD+4cweYkqDjGQ9e34vnfb0aF326gS9NaSfUhLj/Y4FAIhgreqn66YjO2617Mnv5unBL2f7lAgAAjizh5jTuWmWHj6ju6YsXq/b35sm55ihEwAAANGHYg8xaVCXxrr5iA56f9JSPTdyod9xAAAAgDpHsYeYdc2h7XR0l8a674tZGj1/td9xAAAAgDpFsYeYZWZ68PRual0/Tde8MUnL1hf7HQkAAACoM0zQgpiWkZygZ87rpZOeHKUrX5ugty7vz4QtAIA9Zhb4fsp/R6t+RrLqey0ZWtdP0zl9WykuzvwNCCAmcGYPMa9dXob+fUY3TVlSoLs/muF3HABAFDi6SxOdvl9z5WUmq2BTmX5etFavj/1Vd304QwvXbPQ7HoAYwZk9QNJR+zTWNYe003++m6duLXJ0Vp+WfkcCAESwNg3S9eDp3bZa9snUZbrmjUmqrGQWaAChwZk9wHPjER10UIeG+tuHMzTpt3V+xwEAAAD2CMUe4ImPMz0+uLsaZSfrytcmau3GMr8jAQAAALuNYg+oJictSU+ds5/WbizTTW9P5lIbAAAARCyKPWAbXZpl667j9tb3s/P19Ij5fscBAAAAdgsTtADbcW6/VhqzcK3+/dUc9WqVqz5tcv2OBACIEn/7aIaa10tVvfQk5aYFWjIc2bmxstMS/Y4GIMpQ7AHbYWa675R9NWNpga59c6I+u+5A1c9I9jsWACCCdW+Ro4EdG2plYakW5G/U2o1lKttcKUlaclixbjyig88JAUQbij1gBzJTEvXkOT118n9H68a3p+ilC3vTBBcAsNua10vTSxf12fLcOadNZZvV/R9fbSn6AKAucc8eUIN9mmbrb8d31og5+XrqB+7fAwDUHTNTenKCTHyQCCA4KPaAnTi7T0sd362p/v3VbI1buNbvOAAAAECtUOwBO2FmuveUfdUiN003DJukgk3lfkcCAAAAdopiD6iFjOQEPTa4h1YVleovH0yTc/TfAwAAQHhjghaglrq3yNFNR3bQA1/M1sHtG+qM3i38jgQAiBIj5uTLOaleWqJy0hKVnZqkvKxk9WiRIzPu6QOweyj2gF1w+UFtNXLOav3toxnar3U9tW2Y4XckAECEO2zvPI3/dZ2e/3GByjdvfeXIixf21iGd8nxKBiDSUewBuyA+zvTImd016LERun7YJL135f5KSuBqaADA7nvq3P0k/d6KYd2mMk1fWqArXpuoDaUVPqcDEMn4KxXYRY2zU3T/qV01fWmhHvpqtt9xAABRoqoVQ/N6aWqXx5UjAPYcxR6wG47ap7HO6dtSQ0cs0Mi5+X7HAQAAAP6AYg/YTXce21nt8zJ009tTtGZDqd9xAAAAgK1Q7AG7KTUpXo+f1UMFxeW6ZfhU2jEAAAAgrDBBC7AH9m6SpduP7qS/fzxTr/z0qy4Y0NrvSACAKPK/X1aqoLhc2alVLRkSlZOapOb1UhUXR0sGADWj2AP20IUDWmvEnHz967Nf1HevXHVqnOV3JABAhGuQkawGGUn6cPIyfTh52R/WX3toO918ZEcfkgGIJBR7wB4yMz14ejcNenSkrn9zsj68Zn+lJMb7HQsAEMFy0pI0/s4jVFK+WQXF5Vq/qVwFxYGvm9+erPwi7hUHsHPcswfUgQYZyXro9K6avbJI930+y+84AIAokZIYr0ZZKerYOFN92uTqiM6NlJrEB4oAaodiD6gjAzvm6aL9W+ul0Yv03axVfscBAABAjKPYA+rQnwd1UqfGmbpl+BQusQEAAICvKPaAOpSSGGjHUFRSoVuGT6EdAwAAAHxDsQfUsQ6NMnXHsXvr+9n5emn0Ir/jAACi0K9rNul/M1dq3MK1mr2iSMsLirWprIIPGQFshdk4gSA4r18rfT87X/d+Pkv929anHQMAoM7kZabopwVr9NOCNX9Y1zgrRd/fMpBZoQFIotgDgsLM9MBpXWnHAACoc+9c0V8rC0tUUFyuwuKKLS0ZRszJ1xczVqiwpJwxB4AkHy/jNLN4M5tkZp94z9uY2Vgzm2dmb5lZkl/ZgLpAOwYAu4sxEjVJSYxXq/rp6to8Rwe0b6BjuzbR2X1b6oD2DfyOBiDM+HnP3vWSfqn2/H5Jjzjn2klaJ+kSX1IBdYh2DAB2E2MkAGCP+VLsmVlzScdKes57bpIOlTTc2+RlSSf5kQ2oa7RjALArGCMBAHXFrzN7j0q6VVKl97y+pPXOuQrv+RJJzbb3QjMbYmbjzWx8fn5+0IMCe4p2DAB20aNijAQA1IGQF3tmdpykVc65CbvzeufcUOdcL+dcr4YNG9ZxOiA4aMcAoDYYI1EXZiwr1IL8DVqzoVTlmyt3/gIAUcuP2Tj3l3SCmR0jKUVSlqTHJOWYWYL3yWVzSUt9yAYEzXn9WukH2jEAqBljJHZbZkrgz7qLXvx5q+VpSfHKSknUTUd00Bm9W/gRDYBPQl7sOedul3S7JJnZQEl/cs6dY2bvSDpN0jBJF0j6MNTZgGAyM93vtWO47s1J+uiaA5gaG8BWGCOxJ47v2lQtc9O0blPZlrYMhcXlKiwp1xtjf9PYhWsp9oAYE0599v4saZiZ3SNpkqTnfc4D1LkGGcn69xnddMEL43TvZ7/o7yd28TsSgMjAGImdiosz9WhZb7vrPpu2IsRpAIQDX4s959z3kr73Hi+Q1MfPPEAoHNyhoS7ev41eGLVQAzvm6ZBOeX5HAhCGGCMBAHvKzz57QMy6dVBH2jEAAAAgqCj2AB/QjgEAAADBttNiz8weMLMsM0s0s2/MLN/Mzg1FOCCa0Y4BiHyMkYgkJeWbVVhSrspKPmAEYkVt7tk70jl3q5mdLGmRpFMkjZD0WjCDAbGAdgxAxGOMRERITozTp9OW69NpyxVnUmZKorJSE5SdmqgW9dL06ODuSk5ghmgg2tSm2Kva5lhJ7zjnCswsiJGA2EE7BiDiMUYiIvz3nJ6aurhAhSXlXluGwPd5+Rv0+fQV+tO6YrVtmOF3TAB1rDbF3idmNktSsaQrzayhpJLgxgJiB+0YgIjGGImI0Klx1navHvlw8lJdP2xy6AMBCImd3rPnnLtN0gBJvZxz5ZI2Sjox2MGAWFLVjuHln37Vt7NW+h0HQC0xRgIAwtlOz+yZWbykAyS1NrPq2z8ctFRADLp1UEeNnr9at7wzVZ/fcKDyMlP8jgRgJxgjAQDhrDatFz6WdKGk+pIyq30BqENV7Rg2lFbolnem0o4BiAyMkQCAsFWbe/aaO+e6Bj0JgC3tGP764Qy9NHqRLtq/jd+RANSMMRIAELZqU+x9bmZHOue+CnoaALRjACILYySiwtPfz1eL3DRlp/7ekiErJVH7NM1WahKzRAORqjbF3hhJ75tZnKRySSbJOef4CxQIAtoxABGFMRIRrUOjTLXITdWn05ZrU9nmP6w/fb/mevD0bj4kA1AXalPsPSypv6RpjpuIgJCgHQMQMRgjEdH2bpKlkbceKkkqq6hUUUm5CksqVFBcrmvemKiikgqfEwLYE7WZoGWxpOkMYkBo0Y4BiAiMkYgaSQlxqp+RrDYN0tW9RY7Sk2pzTgBAOKvN/8ULJH1vZp9LKq1a6JxjWmkgyGjHAIQ9xkgAQNiqzZm9hZK+kZSk36eUzghmKAABKYnxeqJaO4bKSk4eAGGGMRIAELZqc2bvFefcwuoLzKx3kPIA2Eb7Rpm689i9dZfXjuHiA2jHAIQRxkgAQNiqzZm94WbWrOqJmR0k6YXgRQKwrXP7tdJhnfJ03+ez9MvyQr/jAPgdYySi2qI1G/XmuN/06dTlGjk3X1OXrNei1RtVVFLudzQAtVCbM3tXSPrAzI6X1FPSvZKOCWoqAFup3o7h+mG0YwDCCGMkotZeDdP1+fQVuv29aX9YlxQfp+9vGaimOak+JANQWzst9pxzP5vZdZK+klQi6XDnXH7QkwHYCu0YgPDDGIlo9t9zeqq4fLMKistVWFzhfS/XTwvW6PkfF2rtxjKKPSDM7bDYM7OPJVWfDSJNUoGk581MzrkTgh0OwNaq2jG8MGqhDu7YUId2auR3JCAmMUYiFpiZ0pISlJaUoCbZvy93kp7/ceEOXwcgfNR0Zu+hkKUAUGu0YwDCAmMkACDs7bDYc879EMogAGqnqh3DcU/8qFvemaoXL+ytuDjzOxYQUxgjAQCRoDazcQIIM1XtGH6Yk6+XRi/yOw4AAADCUG1m4wQQhs7t10o/zMnXfZ/PUv+29bV3kyy/IwEAYsjYhWu1sbRCWamJyk5NVFZqotKT4mXG1SZAuKDYAyKUmen+U7tq0GMjdd2bk/TxtbRjAAAEX256kiTpn5/M/MO6+DjTFQfvpVuO6hTqWAC2Y6fFnpntL+luSa287U2Sc87tFdxoAHamfkay/n16N53/wjj932e/6B+0YwBCijESsWi/VvX00+2HanVRmQpLAu0YCorLVVhSrud/XKgZywr9jgjAU5sze89LulHSBEmbgxsHwK46qENDXXJAGz3/40Id3KGhDtubdgxACDFGIiY1yU5Vk+w/9tj7dOpyH9IA2JHaTNBS4Jz73Dm3yjm3puor6MkA1NotR3VUp8aZunX4VK0qKvE7DhBLGCMBAGGrNsXed2b2oJn1N7OeVV9BTwag1qraMWwordAt70xVZaXb+YsA1AXGSABA2KrNZZx9ve+9qi1zkg6t+zgAdldVO4a7Ppyhl0Yv0sUHtPE7EhALGCMBAGFrp8Wec+6QUAQBsOdoxwCEFmMk8EdFJRWat2qD144hQckJzBQN+GWHxZ6Zneuce83Mbtreeufcw8GLBWB30I4BCA3GSGD70pIS9NOCNTr84R+2LEtJjFNWSqKa5KTq+Qt6qUFGso8JgdhS05m9dO97ZiiCAKgbtGMAQoIxEtiOx8/qoenLClRY/HtLhoLics1dtUHfz87Xr2s2UuwBIbTDYs8594z3/e+hiwOgLtCOAQguxkhg+xpmJuuQjnl/WD5iTr6+n53vQyIgttVmNk4AEejWQR21d5Ms3UI7BgAAgJhEsQdEqeSEeD0+uLs2llboT7RjAAAAiDkUe0AUq2rHMGJOvl4avcjvOAAAAAihnbZeMLNkSadKal19e+fcP4IXC0BdoR0DEDyMkcCu2VS2WZWVTnFx5ncUICbUpqn6h5IKJE2QVBrcOADqGu0YgKBijARqISkhcDHZec+Pk5mUmZygrNTEQC++lESd0bu5Tu7R3OeUQPSpTbHX3Dk3KOhJAAQN7RiAoGGMBGqhV6t6evysHlpVWBJoy1BSsaU1w7iFaxUXJ4o9IAhqU+yNNrN9nXPTgp4GQNDQjgEICsZIoBYS4uN0Qrem21132lOjQ5wGiB07LPbMbJok521zkZktUOASFZPknHNdQxMRQF25dVBHjZ6/RrcMn6ovbjhQeZkpfkcCIhJjJAAgEtR0Zu+4kKUAEBJV7RiOe+JH/emdqXrpwt7cJA/sHsZIAEDY22HrBefcr865XyXdU/W4+rLQRQRQl9o3ytSdx3XWiDn5epF2DMBuYYwEAESC2vTZ26f6EzOLl7RfcOIACIVz+7bU4Xvn6f7PZ2nmskK/4wCRjDESABC2arpn73ZJf5GUamaFCtyHIEllkoaGIBuAIKnejuH6YbRjAHYVYyRQd+LiTD/NX6Ne93ytrJTE31sypCaqSXaKbjy8g1KTGKOA3VHTZZz3OucyJT3onMtyzmV6X/Wdc7eHMCOAIKhqxzB31Qb932e/+B0HiCiMkUDd+dORHTXkoLY6onNj7d00S5kpCVq/qUzjF63V0BELNH1Zgd8RgYhVm9YLfzGzUyQdoMDMYyOdcx8ENRWAkKAdA7DHGCOBPdSnTa76tMn9w/JR81brnOfGyjkfQgFRojb37D0p6QpJ0yRNl3SFmT0Z1FQAQubWQR21d5Ms3TJ8qlYVlfgdB4g0jJEAgLBVm2LvUElHOededM69KOkYbxmAKFDVjmFjaYX+9M5UVVbyESqwCxgjAQBhqzbF3jxJLas9b+EtAxAlaMcA7DbGSABA2KpNsZcp6Rcz+97MvpM0U1KWmX1kZh8FNx6AUAm0Y2hEOwZg1zBGAgDCVm0maPlr0FMA8F2gHcO+tGMAdg1jJBBktw6forysFGWlVLVkSFBWSqIO37uR9m2e7Xc8IKzttNhzzv1gZq0ktXfO/c/MUiUlOOeKgh8PQCjVz0jWw2d003nPj9O/Pv1F/zypi9+RgLDGGAkET9fm2TqvXyutLCxRQXG5lq4v1i/LC1VYXK6i0gpN/G2dXr2kr98xgbC202LPzC6TNERSrqS2kppLelrSYcGNBsAPB7ZvqEsPaKPnflyogR1pxwDUhDESCJ7MlMQdfuh4xtM/qWIzE4oBO1Obe/aulrS/pEJJcs7NlZQXzFAA/HUL7RiA2mKMBACErdoUe6XOubKqJ2aWoEDjWABRinYMQK0xRgIAwlZtir0fzOwvklLN7AhJ70j6OLixAPiNdgxArTBGAgDCVm2Kvdsk5UuaJulySZ9JujOYoQCEB9oxADvFGAkACFu1mY2z0sw+kPSBcy4/+JEAhAvaMQA1Y4wE/PPLikJd/frEQCuG1MQtrRmaZKfokI55ioszvyMCvtthsWdmJulvkq6RdwbQzDZLesI594/QxAPgN9oxAH/EGAn464TuTfXuxCX6ZUWhCosrVFhcrrLNlVvWf3D1/ureIse/gECYqOnM3o0KzDDW2zm3UJLMbC9JT5nZjc65R0IREID/qrdjOLhDQx3emXYMiHmMkYCPzu3XSuf2a7XluXNOpRWV+nbWKl31+kSVlm/2MR0QPmq6Z+88SWdVDWKS5JxbIOlcSecHOxiA8HLLoI7q3CRLt747VasKaceAmMcYCYQRM1NKYrxyUhP9jgKElZqKvUTn3OptF3r3JPB/EhBjkhPi9fhZgXYMN78zhXYMiHWMkQCAsFdTsVe2m+sARKl2eYF2DCPnrqYdA2IdYyQAIOzVdM9eNzPb3lzrJiklSHkAhLlz+7bUD7Pzdf/ns9R/r/rq3DTL70iAHxgjAQBhb4fFnnOO+dUB/AHtGADGSCDcvTR6kUbNWx1oyVCtLcM+zbKUlcKV1ogdO+2zFwxmliPpOUldJDlJF0uaLektSa0lLZJ0hnNunR/5ANSMdgxA8DBGAruvdYN0dWiUodHz1+iLGSvktrm9/Kh9GumZ83r5Ew7wgS/FnqTHJH3hnDvNzJIkpUn6i6RvnHP3mdltkm6T9Gef8gHYCdoxAEHDGAnspqY5qfrqxoMlSZWVThvKKlSwqVwFxeX60ztTtKG0wueEQGjVNEFLUJhZtqSDJD0vSc65MufcekknSnrZ2+xlSSeFOhuAXUM7BqBuMUYCdScuzpSVkqgWuWnq0ixbGcl+neMA/BPyYk9SG0n5kl40s0lm9pyZpUtq5Jxb7m2zQtJ2TxOY2RAzG29m4/Pz80MUGcD2VLVj2FRGOwagjjBGAgDqjB/FXoKknpKecs71kLRRgctRtnDOOQXuU/gD59xQ51wv51yvhg0bBj0sgJq1y8vUnccG2jG8MGrhzl8AoCaMkQCAOuNHsbdE0hLn3Fjv+XAFBraVZtZEkrzvq3zIBmA3nNO3pQ7fu5Ee+GK2Ziwr8DsOEMkYIwEAdSbkxZ5zboWkxWbW0Vt0mKSZkj6SdIG37AJJH4Y6G4DdU9WOITstUdcPm6ziss1+RwIiEmMkEFxL1hXrlZ8W6YNJS/XdrFWa8Os6zVtVpIJN5X5HA4LCrztVr5X0ujfL2AJJFylQeL5tZpdI+lXSGT5lA7Abqrdj+L/PaMcA7AHGSCAI2jfK0Phf1+mvH874w7r4ONP/bjpYbRqk+5AMCB5fij3n3GRJ22tycliIowCoQ7RjAPYcYyQQHPee0lV3n7CPikoqVFBcrsLiQEuGnxet1ZPfzdfajaUUe4g6ftyzByCK0Y4BABCukhPi1SAjWW0bZqhHy3oa2DFPfdvU9zsWEDQUewDqFO0YAAAAwgPFHoA6RzsGAAAA/1HsAQiKc/q21BGdaccAAADgF4o9AEERaMfQVTm0YwAARIAxC9Zq1LzVmr60QL+t2aT1m8q0mVsREOH8ar0AIAbkpifp3147hn99NlP3nLSv35EAANhKbnqSzKQHv5y93fWXH7yXbj967xCnAuoGxR6AoDqwfUNddmAbPTtyoQZ2yKMdAwAgrHRplq2f7zhc+UWlW9oxFHrtGZ4dsUCzlhf5HRHYbRR7AILuT0d11Kh5a3Tru1P1RfMDlZeV4nckAAC2aJCRrAYZyX9Y/tGUZT6kAeoO9+wBCDraMQAAAIQexR6AkKAdAwAAQGhR7AEIGdoxAAAAhA7FHoCQoR0DACDSFJaUa/aKIq0oKNGmsgo5x60IiBxM0AIgpGjHAACIFBnJ8Ro1b42OenTElmWJ8abs1EQ1zUnVKxf3UU5ako8JgZpR7AEIuertGA7ukKcjaMcAAAhDT5zVUzOWFQTaMRRXeG0ZyvXL8kJ9Pztfi9cWU+whrFHsAfBFVTuGP787Vd1oxwAACEO56Uk6sH3DPyz/38yV+n52vg+JgF3DPXsAfEE7BgAAgOCi2APgG9oxAAAABA/FHgBf0Y4BAAAgOCj2APiKdgwAgEhVVFKuis2VfscAdogJWgD4jnYMAIBIkpQQOF9y9nNjJUnpSfHKTk1UVmqislISdf6AVjqua1M/IwKSKPYAhAnaMQAAIkX/tvX133N6alVhiQqKK1RYUq6C4sDXqHmr9fGUZRR7CAsUewDCBu0YAACRIDE+Tsfs22S76wZVa8AO+I179gCEjUA7hh60YwAAAKgDFHsAwkq7vAzddRztGAAAAPYUxR6AsHN2H9oxAAAA7CmKPQBhh3YMAAAAe44JWgCEJdoxAAAiUZyZvp65Uj3+8ZWyUhOV7X1lpSSqeW6qbj2qk+LjzO+YiBEUewDCFu0YAACR5q7jOuvHefkqLK5QQXH5lrYMc1du0KfTluuMXi3UtmGG3zERIyj2AIQ12jEAACJJ/7b11b9t/T8s/2jKMl335iQ5JppGCHHPHoCwRjsGAACA3UOxByDs0Y4BAABg11HsAYgItGMAAADYNRR7ACJC9XYM1705iXYMAAAAO8EELQAiRm56kh4+o7vOfX6s7vl0pv51Mu0YAACRoarZwqUv/6wGGclb2jJkpSQoOzVRJ3RvqnZ5mb5mRPSh2AMQUQ5o30BDDtpLQ0cs0MCOtGMAAESG/m3r69x+LbW6qEyFJeVaWViiOSuLVFhcrsKSCi0vKNGDp3fzOyaiDMUegIhz85Ed9OPc1bRjAABEjAYZybrnpO1fkXLA/d9qMz0ZEATcswcg4tCOAQAAYOco9gBEJNoxAAAA1IxiD0DEoh0DAADAjlHsAYhYtGMAAADYMSZoARDRaMcAAIgGo+at1pWvTVBWSqKy035vy9CmQYYOaN/A73iIUBR7ACJe9XYMB3doqCP3aex3JAAAau2sPi313axVmrdqgwpLylVQXK6S8sot6yfddYTqpSf5mBCRimIPQFT405EdNWqe146hRY4a0Y4BABAhrj6kna4+pN1Wy0orNuuV0b/qX5/9ovLNlTt4JVAz7tkDEBWSEuL02OAeKi7frD/RjgEAEOGSE+KVlhzvdwxEOIo9AFGDdgwAAAC/o9gDEFVoxwAAABBAsQcgqtCOAQAAIIAJWgBEHdoxAACiyQNfzlbjrBRlpSZ4LRkC7Rl6tqynlETu68OOUewBiEq0YwAARLouTbO1V8N0/e+XlSosLte2c49de2g73XxkR3/CISJQ7AGIWrRjAABEsm4tcvTtzQMlSc45bSitUEFxuQqLK3Ta06O1obTC34AIe9yzByBqVW/HcPPbtGMAAEQuM1NmSqKa10tT56ZZio8zvyMhAlDsAYhq7fIy9Nfj9tGP81br+R9pxwAAAGIHxR6AqHdWnxY6ap9GeuDLWZq+lHYMAAAgNlDsAYh6Zqb7Tumq3PQkXT+MdgwAACA2UOwBiAn10pP079O7a37+Rt3z6Uy/4wAAsMemLinQy6MX6f1JS/TtrJUav2it5q4sUlFJud/RECaYjRNAzKAdAwAgWuzdJEvjFq7VhF/X/WFdvbREjb/zCCZxAcUegNhCOwYAQDR4+/L+Kt9cqcLichWWBFoyFBSX690JS/TRlGWqqKxUfBwN12Mdl3ECiCm0YwAARIvE+DjVz0hWmwbp6t4iRwd3aKiOjTP9joUwQrEHIObQjgEAAMQCij0AMYl2DAAAINpR7AGISbRjAAAA0Y5iD0DMqpeepIfP6K4Fqzfqn7RjAABEkS+mr9CPc1dr6pL1+nXNRq3fVKbN3Kcec5iNE0BM279dAw05cC89M2KBBtKOAQAQ4fIykyVJ1w+b/Id18XGmxwf30LFdm4Q4FfxCsQcg5t18ZEf9SDsGAEAUOL1XCx3UoaHWbSpTYXGgJUNhcbnWbSrTPZ/+ogX5G/yOiBDiMk4AMY92DACAaNIoK0WdGmepT5tcHdG5kU7dr7kuHNDa71jwAcUeAIh2DAAAIPpQ7AGAh3YMAAAgmvhS7JnZjWY2w8ymm9mbZpZiZm3MbKyZzTOzt8wsyY9sAGIX7RgQDhgjAQB1JeTFnpk1k3SdpF7OuS6S4iUNlnS/pEecc+0krZN0SaizAQDtGOAnxkgAwbasoESzVxRpRUGJNpVVyDnuU49mfs3GmSAp1czKJaVJWi7pUElne+tflnS3pKd8SQcgplVvx3Bwh4Y6inYMCC3GSAB1Ls5MaUnxenPcb3pz3G9blifGm7JSEtW/bX395+yePiZEMIS82HPOLTWzhyT9JqlY0leSJkha75yr8DZbIqnZ9l5vZkMkDZGkli1bBj8wgJh085EdNWr+at327lR1px0DQoQxEkCwxMWZvrrxIC1avUkFxeVbvgpLyvXdrFX6YU6+3xERBH5cxllP0omS2khqKild0qDavt45N9Q518s516thw4ZBSgkg1lVvx3DT25Npx4CQYIwEEEzN66XpgPYNdGzXJjq7b0tdObCt/jyok/q3re93NASJHxO0HC5poXMu3zlXLuk9SftLyjGzqjONzSUt9SEbAGzRtmGgHcOoeWtox4BQYYwEANQZP4q93yT1M7M0MzNJh0maKek7Sad521wg6UMfsgHAVmjHgBBjjAQA1JmQF3vOubGShkuaKGmal2GopD9LusnM5kmqL+n5UGcDgG3RjgGhxBgJAKhLvvTZc879zTnXyTnXxTl3nnOu1Dm3wDnXxznXzjl3unOu1I9sALAt2jEglBgjAfjBOWntxjKVb670OwrqkF+tFwAgotCOAQAQrZIT4rWhtEI9//m1JCk9KV7ZqYnKSk1U/Ywk/f2ELmqXl+FzSuwOij0AqCXaMQAAotEVB++ljo0zVLCpXAXFFSosCbRlWFVUqhFz8jV+0VqKvQhFsQcAtVTVjuG4x3/UTW9P1qsX91VcnPkdCwCAPZKTlqSTezT/w/LlBcXqf++3PiRCXfHlnj0AiFRtG2bor8d31qh5a/Tcjwv8jgMAALBDFHsAsIsG9w60Y3jwy9m0YwAAAGGLYg8AdlH1dgzX0Y4BAACEKYo9ANgNVe0YFtKOAQAAhCkmaAGA3UQ7BgBANIu3wCRkd3wwXQ98OVtZKQlbWjJkpSbqsE55OqXnHyd2Qfig2AOAPUA7BgBAtMrLStGDp3XVojUbVVBcrsLiChUUB9oyTP5tvWYtL6TYC3MUewCwB2jHAACIZqf3arHd5Ve/MVGzlheGOA12FffsAcAeoh0DAAAIRxR7AFAHaMcAAADCDcUeANSBqnYM9dOTdd2wSdpUVuF3JAAAEOMo9gCgjgTaMXQLtGP45Be/4wAAgBjHBC0AUIcGtGugIQftpWd+WKCBHWnHAACITibp1zWbdOQjPyg7NTHQkiEl0JKhUVaKLtq/tVIS4/2OGfMo9gCgjt18REeNmrdaf353qro1z1HjbNoxAACiywUDWis1MV6FJYGWDMvWl+iX4iKt31SmjWWb1aVZlg5s39DvmDGPYg8A6lj1dgw3v0M7BgBA9OndOle9W+f+YfnE39bplP+O1uZK50MqbIt79gAgCNo2zNDfvHYMz46kHQMAAAg9ij0ACJIze7fQoH0a66GvaMcAAABCj2IPAILEzHTfqfsG2jG8STsGAAAQWhR7ABBEOWlJevjMblq4ZqP+8fFMv+MAAIAYwgQtABBkA9o20BUHt9VT38/XwR0a6uh9m/gdCQCAoHroq9l6e/zirVoyZKcm6oB2DdS6Qbrf8WIGxR4AhMCNh3fQqHmrddt709S9ZY6aZKf6HQkAgDrXPi9Dx3ZtouXrizV7RZEKSypUUFyusopKSdLhe+fpuQt6+5wydlDsAUAIVLVjOPbxkbrxrcl6/dJ+iqcdAwAgymSmJOrJs3v+YXlJ+Wad/vRPKvWKPoQG9+wBQIi0aZCuu0/YR2MWrNUzI+b7HQcAgJBJSYxXQjwfcoYaxR4AhNDp+zXXsfs20cNfzdGUxev9jgMAAKIYxR4AhJCZ6f9O3ld5mcm6ftgkbSylHQMAAAgOij0ACLHstEQ9cmZ3/bZ2k+7+aIbfcQAAQJRighYA8EHfverrqoHt9J/v5ungjg11XNemfkcCACDo5qws0t8+nB5oyZD6e0uG5vVStU/TbL/jRR2KPQDwyfWHt9eP81br9vemqUfLemqWQzsGAED0OqRjnt6ZsFjvT1qqotIKObf1+tG3HaqmjIV1imIPAHySGB+nxwZ31zGPjdSNwybrzSG0YwAARK/rDmuv6w5rL0mqrHQqKq1QYXG5Pp22XPd9PkubyriPva5xzx4A+KhV/XT948QuGrdorZ76fp7fcQAACIm4OFN2aqJa5KZxZUsQUewBgM9O6dlMJ3Rrqkf+N1cTf1vndxwAABAlKPYAwGdmpntO7qLGWSm6YdhkFZWU+x0JAABEAYo9AAgDWSmJemxwdy1Zt0l/ox0DAACoA0zQAgBholfrXF17aHs99s1cHdyhoU7s3szvSAAAhMwbYxerTYO0rVoy5KQmqnX9dMUxgdluodgDgDBy7aHt9OO81brz/enq2bKeWuSm+R0JAICgal0/XZkpCXph1MLtrr/u0Ha66ciOIU4VHSj2ACCMJMTH6dEzvXYMb03WsCH9lBDPFfcAgOi1b/NsTbv7KJVWbFZBcbkKiyu87+W6btgkrd5Y5nfEiMVfEAAQZlrkpumek7to/K/r9OR38/2OAwBASCQnxCsvM0Xt8jK0X6t6OqRTnpIT4v2OFdEo9gAgDJ3YvZlO7tFMj30zRxN+Xet3HAAAEIEo9gAgTP3jxH3UrF6qrh82WYW0YwAAALuIYg8AwlRmSqIePbOHlheU6K8fTPc7DgAAiDAUewAQxvZrVU/XH9ZeH0xepvcnLfE7DgAAIbcwf6O+mL5co+ev1vSlBVq8dpMKS8pVWen8jhb2mI0TAMLc1Ye008i5+brrgxnar2WuWtanHQMAIDY0yU7RTwvW6KcFa/6wrmVumr7/00B68NWAYg8Awlx8nOmRM7vr6MdG6vq3Jumdy/vTjgEAEBPevXKAVhWVbGnHUNWS4YsZK/TtrFWqqHRKotjbIYo9AIgAzeul6V8n76vr3pykx7+ZS3NZAEBMSEqIU/N6aVK9rZfnbyjVt7NW+RMqgvDRMABEiBO6NdWpPZvrP9/N07iFtGMAAAA1o9gDgAjy9xP3UYvcNN341mQVFNOOAQAA7BjFHgBEkIzkBD02uIdWFpbojvenyTlmIgMAANtHsQcAEaZ7ixzdeEQHfTJ1ud6duNTvOAAA+GbSb+s0b1WRVhWVqKR8s99xwg4TtABABLri4LYaMSdff/twunq1qqfWDdL9jgQAQMhkpgTKmDOHjtlqeXJCnHLTk/TvM7ppQNsGfkQLKxR7ABCBqtoxDHp0hK4fNknDrxygRNoxAABixDl9W6lb8xyt21QWaMdQUqHC4nKt3lCqF0ct0oylhRR7otgDgIjVNCdV953aVVe9PlGP/m+Objmqk9+RAAAIifg4U7cWOX9YvqG0Qi+OWhTyPOGKj4EBIIIds28TndGruf77/XyNWbDG7zgAACCMUOwBQIT72/H7qHX99EA7hk20YwAAAAEUewAQ4dKTE/TY4O7KLyrV7e9PpR0DAACQRLEHAFGha/Mc3XxkR302bYXeGb/E7zgAAPiqsKRcBcXlqqyM7Q9AmaAFAKLE5QftpZFz83X3xzPUq3U97dUww+9IAACEVEKcKT7O9MS38/TEt/NkJmUmJygrNVHZqYk6bO9GuumIDn7HDBmKPQCIEnFxpofP6K5Bj43Q9cMm690rBygpgQs4AACxIyUxXu9c0V/zV20ItGTw2jIUFJdr7II1em/iEoo9AEBkapydovtO6aorXpugh7+eo9uOph0DACC29GxZTz1b1vvD8pvenqxxC9f6kMg/fOQLAFFmUJfGOqtPSz0zYr5Gz1vtdxwAAOATij0AiEJ3Hbe32jRI141vT9a6jWV+xwEAAD6g2AOAKJSWlKDHB/fQ2o1luu092jEAABCLKPYAIEp1aZatW4/qpC9nrNSwnxf7HQcAAN/F2mefTNACAFHskgPaaMTcfP394xnq3TpX7fJoxwAAiE0Jcaal64vV4c7PlZWSqOzUBGV7LRkaZaXotqM7KSctye+YdYozewAQxeLiTA+d3k2pifG6ftgklVZs9jsSAAC+uGpgO91yVEddNKC1juicp46NM5WWlKDF64o17OfFmvTber8j1jnO7AFAlGuUlaIHTuumy14Zr39/NUd/OWZvvyMBABByrRuk6+pD2v1h+eTF63XSk6N8SBR8QTuzZ2YvmNkqM5tebVmumX1tZnO97/W85WZmj5vZPDObamY9g5ULAGLREZ0b6dx+LTV0xAKNnJvvd5yYxxgJAAiFYF7G+ZKkQdssu03SN8659pK+8Z5L0tGS2ntfQyQ9FcRcABCT7jims9rlZejmt6doLe0Y/PaSGCMBAEEWtGLPOTdC0rYt6k+U9LL3+GVJJ1Vb/ooLGCMpx8yaBCsbAMSi1KR4PT64h9ZvKtetw2nH4CfGSABAKIT6nr1Gzrnl3uMVkhp5j5tJqj4v+BJv2XJtw8yGKPDJpiRtMLPZQcoqSQ0krQ7i+/st2vdPiv59jPb9k6J/H33Zv7mSnr8wJD+qLvevVR29T7gKlzEyWv+fY78iRzTuk8R+hb1D79/qaSTt1w7HR98maHHOOTPb5Y+VnXNDJQ0NQqQ/MLPxzrleofhZfoj2/ZOifx+jff+k6N9H9g/b4+cYGa3HjP2KHNG4TxL7FWmiZb9C3XphZdWlJ973Vd7ypZJaVNuuubcMAIBYwRgJAKhToS72PpJ0gff4AkkfVlt+vjfjWD9JBdUuZQEAIBYwRgIA6lTQLuM0szclDZTUwMyWSPqbpPskvW1ml0j6VdIZ3uafSTpG0jxJmyRdFKxcuygkl4v6KNr3T4r+fYz2/ZOifx/ZvxgU5mNktB4z9ityROM+SexXpImK/TJmYwMAAACA6BPqyzgBAAAAACFAsQcAAAAAUSjmiz0zizezSWb2ife8jZmNNbN5ZvaWmSV5y5O95/O89a19DV5LZpZjZsPNbJaZ/WJm/c0s18y+NrO53vd63rZmZo97+zjVzHr6nX9nzOxGM5thZtPN7E0zS4n0Y2hmL5jZKjObXm3ZLh8zM7vA236umV2wvZ/lhx3s34Pef6NTzex9M8uptu52b/9mm9lR1ZYP8pbNM7PbQrwbNdrePlZbd7OZOTNr4D2PimPoLb/WO44zzOyBassj7hjGqkg+Jma2yMymmdlkMxvvLYu48S5ax4Ad7NfdZrbUO2aTzeyYauvC/veGmbUws+/MbKb3e+96b3lEH68a9ivSj1eKmY0zsynefv3dW97GdvHvxh3tb1hyzsX0l6SbJL0h6RPv+duSBnuPn5Z0pff4KklPe48HS3rL7+y13L+XJV3qPU6SlCPpAUm3ectuk3S/9/gYSZ9LMkn9JI31O/9O9q2ZpIWSUqsduwsj/RhKOkhST0nTqy3bpWMmKVfSAu97Pe9xPb/3rYb9O1JSgvf4/mr711nSFEnJktpImi8p3vuaL2kv77/rKZI6+71vNe2jt7yFpC8VmHyjQZQdw0Mk/U9Ssvc8L5KPYSx+RfoxkbSo6v+rassibryL1jFgB/t1t6Q/bWfbiPi9IamJpJ7e40xJc7zsEX28ativSD9eJinDe5woaax3HHbp78Yd7a9f+7Wzr5g+s2dmzSUdK+k577lJOlTScG+TlyWd5D0+0Xsub/1h3vZhy8yyFfjl+rwkOefKnHPrtfW+bLuPr7iAMZJyzOv5FMYSJKWaWYKkNEnLFeHH0Dk3QtLabRbv6jE7StLXzrm1zrl1kr6WNCjo4Wthe/vnnPvKOVfhPR2jQB8xKbB/w5xzpc65hQrMRtjH+5rnnFvgnCuTNMzbNizs4BhK0iOSbpVUfWasqDiGkq6UdJ9zrtTbpqpHXEQewxgVjcck4sa7aB0Davi9uD0R8XvDObfcOTfRe1wk6RcFPoiO6ONVw37tSKQcL+ec2+A9TfS+nHb978Yd7W9YiuliT9KjCvzhVek9ry9pfbU/Opfo9/+4m0laLEne+gJv+3DWRlK+pBctcKnqc2aWLqmR+71H0wpJjbzHW/bRU33/w45zbqmkhyT9pkCRVyBpgqLrGFbZ1WMWUcdyGxcr8MmnFEX7Z2YnSlrqnJuyzapo2ccOkg70LnX5wcx6e8ujZf9iQaQfEyfpKzObYGZDvGVRMd4puseAa7xLGl+outxREbhf3iV+PRQ4WxQ1x2ub/ZIi/HhZ4PatyZJWKVBUz9eu/90YdvtVk5gt9szsOEmrnHMT/M4SRAkKXDLxlHOuh6SNClxOsIVzzmnrswwRw/slc6ICRW1TSekKkzMfwRTJx2xnzOwOSRWSXvc7S10yszRJf5H0V7+zBFGCApcg9ZN0iwL94sLuzDmi2gHOuZ6SjpZ0tZkdVH1ltPzujJb98Dwlqa2k7gp8aPtvX9PsJjPLkPSupBucc4XV10Xy8drOfkX88XLObXbOdVfgCqI+kjr5myj4YrbYk7S/pBPMbJECp5UPlfSYAqfUq5rNN5e01Hu8VIH7beStz5a0JpSBd8MSSUucc1WfxgxXoPhbWXW5ive96nKrLfvoqb7/4ehwSQudc/nOuXJJ7ylwXKPpGFbZ1WMWacdSZnahpOMkneMNjlL07F9bBT6UmOL9zmkuaaKZNVb07OMSSe95l8mMU+CKiQaKnv2LBRF9TLyrPaouIX5fgT/komW8i8oxwDm30vvju1LSs/r9UriI2S8zS1SgIHrdOfeetzjij9f29isajlcVF7it6TtJ/bXrfzeG7X5tT8wWe865251zzZ1zrRW46fJb59w5Chz407zNLpD0off4I++5vPXfVvuDNCw551ZIWmxmHb1Fh0maqa33Zdt9PN+bLaqfpIJqlyGEo98k9TOzNO8MQtX+Rc0xrGZXj9mXko40s3reGdAjvWVhycwGKXBJ9QnOuU3VVn0kabA3I1YbSe0ljZP0s6T23gxaSQr8P/xRqHPXlnNumnMuzznX2vuds0SBm99XKEqOoaQPFJikRWbWQYGb8VcrSo5hjIjYY2Jm6WaWWfVYgf9fpit6xruoHAO2uU/yZAWOmRQhvze8vz2el/SLc+7haqsi+njtaL+i4Hg1NG+2bzNLlXSEAvcj7urfjTva3/DkwmCWGL+/JA3U77Nx7qXAAZsn6R39PrNcivd8nrd+L79z13LfuksaL2mqAn+M1VPgeuNvJM1VYPa8XG9bk/SkAtcvT5PUy+/8tdi/v0uapcAvnFcVmBkpoo+hpDcVuDyiXIGi4JLdOWYK3Ps2z/u6yO/92sn+zVPg+vfJ3tfT1ba/w9u/2ZKOrrb8GAVmCJsv6Q6/92tn+7jN+kX6fTbOaDmGSZJe8/5fnCjp0Eg+hrH6FanHxPu9P8X7mlGVPRLHu2gdA3awX696uacq8Ad0k2rbh/3vDUkHKHCJ5tRq49cxkX68ativSD9eXSVN8vJPl/RXb/ku/924o/0Nxy/zAgMAAAAAokjMXsYJAAAAANGMYg8AAAAAohDFHgAAAABEIYo9AAAAAIhCFHsAAAAAEIUo9hC1zKy+mU32vlaY2dJqz5P8zledmQ00swFBfP/Ru7j9S2Z2mvc418wmmdlFNWx/t5n9aSfv+ZyZdd6VHACAumFmm73xb7qZvWNmaX5n2hkza2pmw/3OEW68vxk+8TsHIgPFHqKWc26Nc667c667pKclPVL13DlXFuo8ZpZQw+qBknap2NvJ+23FObdbhaSZZSvQ2HWoc+7F3XmPahkudc7N3JP3AADstmJv/OsiqUzSFdVX7sqYsqdq+7Occ8ucc6ftfEsAO0Kxh5hiZvuZ2Q9mNsHMvjSzJt7y783sETMbb2a/mFlvM3vPzOaa2T3eNq3NbJaZve5tM7zqk9GdvO+jZjZe0vVmdryZjfXOlP3PzBqZWWsFBt0bvU9dD6x+Zs17nw3e94FmNtLMPpI008zizexBM/vZzKaa2eU72O/qr//ey161L7aDf64MSZ9LesM595T3+rZm9oW3nyPNrNM2P6eTmY2r9ry1mU2r9m/RqyqPmf3LzKaY2Rgza7QrxxEAsEdGSmpX2zHFzJqY2YhqZwYP9LZ9yXs+zcxu9Lat/ru+gZkt8h5faGYfmdm3kr4xs3Qze8HMxnlj4onbhvTGkOnVXv+eNwbNNbMHtrdjZrbIzO71so43s57euDzfzK6ott0t1fbz79WW32Vms83sRzN707yrVszsMm/7KWb2brXx/yUze9r7WXPM7Lhq2Uea2UTva8CO/i295YO87aaY2Tfesj5m9pP37zPazDruyUFHbArZpzhAGDBJT0g60TmXb2ZnSvqXpIu99WXOuV5mdr2kDyXtJ2mtpPlm9oi3TUdJlzjnRpnZC5KuMrPHdvK+Sc65qoGvnqR+zjlnZpdKutU5d7OZPS1pg3PuIW+7S2rYj56SujjnFprZEEkFzrneZpYsaZSZfeWcW1jD63tI2kfSMkmjJO0v6cftbPewpOecc49UWzZU0hXOublm1lfSfyUdWrXSOTfLzJLMrI2X4UxJb23nvdMljXHO3eEN2JdJuqeGzACAOmCBs2pHS/rCW7TTMUXSKZK+dM79y8ziJaVJ6i6pmXemUGaWU4sf31NSV+fcWjP7P0nfOucu9l47zsz+55zbWMPruyswhpVKmm1mTzjnFm9nu9+cc929sfslBca5FEnTJT1tZkdKai+pjwJ/G3xkZgdJKpZ0qqRukhIlTZQ0wXvP95xzz3r7eo+kSxQY+yWptfdebSV9Z2btJK2SdIRzrsTM2kt6U1IvSWdrm39LM2so6VlJB3nHIdd731mSDnTOVZjZ4ZL+z8sH1BrFHmJJsqQukr62wMmseEnLq63/yPs+TdIM59xySTKzBZJaSFovabFzbpS33WuSrlNgwKzpfasXO80lvWWBM39JkmoqynZkXLVi7khJXe33s4DZCgxgNb3vOOfcEm/fJiswSG2v2PtW0olm9pBzbpWZZShwqek79vvJwOTtvO5tBYq8+7zvZ25nmzJJVfcbTJB0RA15AQB7LtX7nS8Fzuw9r8Dv9NqMKT9LesHMEiV94Jyb7I2Ne5nZE5I+lfRVLTJ87ZxbW+1nnWC/3++dIqmlpF9qeP03zrkCSTKzmZJaSdpesVd9PM9wzhVJKjKzUq+wPNL7muRtl+HtZ6akD51zJZJKzOzjau/ZxSvycrztv6y27m3nXKWkud6/SycFxuH/mFl3SZsldfC23d6/5UBJI6qOQ7V/o2xJL3vFolOgAAV2CcUeYokpUMT138H6Uu97ZbXHVc+r/l9x27zG1eJ9q39K+YSkh51zH3m/3O/ewWsq5F1mbWZxChSG23s/k3Stc676oLMz1fdts3b8e2CYAmf+PjOzQ7yftd67B7ImbylQEL4nyTnn5m5nm3LnXNW/ZU0ZAAB1o3jb39/eB3e1GlO8M1/HSnrJzB52zr1iZt0kHaXArQhnKHBFy5bxS4ECrrptf9apzrnZu7APtR2/djaem6R7nXPPVH+Rmd1Qw89+SdJJzrkpZnahAvfaV9ne3wY3SlqpwFnCOEklkuScG7Htv6WkdTv4mf+U9J1z7mQL3PLxfQ35gO3inj3EklJJDc2svySZWaKZ7bOL79Gy6vUKXIrxo6TZu/C+2ZKWeo8vqLa8SIFPFKssUuAyUkk6QTv+NO9LSVd6nxDKzDqYWXrtd6dm3iWc30h6T4GBaqGZne79LPMG+m1fM1+BQfgubf8STgBAeNrumGJmrSSt9C5jfE5STzNrICnOOfeupDsVuERT2nr8qmlylS8lXWtexWlmPep8b2r+2Rd7V6zIzJqZWZ4CH3Aeb2Yp3rrjqr0mU9Jy79/mnG3e73QzizOztpL2UuDvgmxJy70zfucpcNWPtvdvKWmMpIPMrI23TdVlnNX/ZriwzvYeMYViD7GkUoGB534zmyJpsnZxBkwFfoFfbWa/SKon6SlvZs/avu/dCpz1miBpdbXlH0s62bwJWhS4dv9g7/36a+tPQ6t7TtJMSRMtcBP7M6rjs2TOuT9LWiLpVQUGrEu8XDMk/eGGes9bks5V4JJOAEBk2NGYMlDSFDObpMCl+Y9Jaibpe+/S0Nck3e69x0MKFIyTJDWo4Wf9U4EPMqea2QzveUg4576S9IaknywwidhwSZnOuZ8VuAR0qgITlE2TVOC97C5JYxUoCGdt85a/SRrnveYK7zLQ/0q6wBsvO+n3cXygtvm3dM7lSxoi6T1v+6oPSh+QdK+3LVfAYLfY71dSAaiJdwnFJ1U3owMAgOhiZhnOuQ0WmG1zhKQhzrmJNWz/kgJ/G9APEGGJTwkAAACAgKFm1lmB+w1frqnQAyIBZ/YAAAAAIApxzx4AAAAARCGKPQAAAACIQhR7AAAAABCFKPYAAAAAIApR7AEAAABAFPp/iPc/QhoQVOUAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Only run this when in serial. Will fail in parallel\n", "\n", "if GEO.nProcs == 1:\n", "\n", " import matplotlib.pyplot as plt\n", "\n", " moho_average_temperature = Model.temperature.evaluate(Model.Moho_tracers).mean()\n", " moho_average_temperature = GEO.dimensionalise(moho_average_temperature, u.degC)\n", "\n", " print(\"Average Temperature at Moho: {0:5.0f}\".format(moho_average_temperature))\n", "\n", " distances, temperature = GEO.extract_profile(Model.temperature, line = [(180.* u.kilometer, 0.), (180.* u.kilometer, Model.bottom)])\n", " distances, pressure = GEO.extract_profile(Model.pressureField, line = [(180.* u.kilometer, 0.), (180.* u.kilometer, Model.bottom)])\n", "\n", " Fig, (ax1, ax2) = plt.subplots(1,2,figsize=(15,7))\n", " ax1.plot(GEO.dimensionalise(temperature, u.degK), GEO.dimensionalise(distances, u.kilometer))\n", " ax1.set_xlabel(\"Temperature in Kelvin\")\n", " ax1.set_ylabel(\"Depth in kms\")\n", " ax1.set_ylim(100, 0)\n", " ax1.set_title(\"Temperature profile\")\n", "\n", " ax2.plot(GEO.dimensionalise(pressure, u.megapascal), GEO.dimensionalise(distances, u.kilometer))\n", " ax2.set_xlabel(\"Pressure in megapascal\")\n", " ax2.set_ylabel(\"Depth in kms\")\n", " ax2.set_title(\"Pressure profile\")\n", " ax2.set_ylim(100, 0)\n", " plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In a similar fashion, one can extract a vertical profile of the viscosity field." ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "tags": [ "analysis" ] }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "Fig = vis.Figure(figsize=(1200,400), title=\"Viscosity Field (Pa.s)\", quality=3)\n", "Fig.Points(Model.swarm, \n", " GEO.dimensionalise(Model.viscosityField, u.pascal * u.second),\n", " logScale=True,\n", " fn_size=3.0)\n", "Fig.show()" ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "tags": [ "analysis" ] }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/romain/Projects/Project_UWGeodynamics/virtualenv/lib/python3.9/site-packages/numpy/core/_asarray.py:102: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.\n", " return array(a, dtype, copy=False, order=order)\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbwAAAG5CAYAAAD1fYFsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABAzElEQVR4nO3deXhdZbn+8e+TeU6aZmiTNm3TiQ60FAqUwbbQIigooIKgoCBSEEVFnOXn7DnqET0oBwRBBEQEARkEUaQt89CRzvOUoU3SNkkzj+/vj72SpiFJM6+d7PtzXb2ys/baez3Z0Nx93/UO5pxDRERkuAvzuwAREZHBoMATEZGQoMATEZGQoMATEZGQoMATEZGQoMATEZGQoMAT6QEz22hmC324bqWZ5Q72dbvDzL5gZkVejSPb1mpmfzKzn/pdowgo8ERamdmLZvbjDo5fbGYHzCzCOTfDObd8sGtzziU453Z59QRNiJhZJPBr4INejYfa1ioSTBR4Ikc9CFxlZtbu+NXAI865Rh9q8pWZRRznlEwgBtg4COWI9IkCT+Sop4GRwAdaDpjZCOAi4CHv+z1mtth7fJqZrTSzI16X3q/bvO5sM3vTzMrMLM/MrvGOJ5vZQ2ZWYmZ7zew2MwvznptkZq+YWbmZHTSzx9q8n/OeXwJ8Gvim13X4nJl9w8yebPuDmNlvzeyOjn5I72f4jpltMrNSM3vAzGK85xaaWb6ZfcvMDgAPmFm0mf2vmRV6f/7XOzYF2Oq9bZmZLW1bayfXvsjM1nqfy5tmNqub/21E+kyBJ+JxztUAjwOfaXP4cmCLc+69Dl5yB3CHcy4JmOi9FjMbB/wT+B2QDpwErPVe8zsgGcgFFnjXutZ77ifAv4ERwBjv3PY13gs8AvzS6zr8CPBn4AIzS/GuHwFcgRfSnfg0cL5X9xTgtjbPjQJSgXHAEuB7wDzv55gNnAbc5pzbBszwXpPinDu3i+thZnOAPwI3EPiHxT3As2YW3dXrRPqLAk/kWA8Cn2hp8RAIpAc7ObcBmGRmac65Sufc297xTwH/cc496pxr8O5rrTWzcAJB9B3nXIVzbg9wO4Eu05b3GwdkOedqnXOvd6dg59x+4FXgMu/QBcBB59yqLl52p3Muzzl3GPgZcGWb55qBHzjn6rx/BHwa+LFzrtg5VwL8qE3NPbEEuMc5945zrsk59yBQRyBMRQacAk+kDS9kDgKXmNlEAq2Zv3Ry+nUEWkdbzGyFmV3kHR8L7Ozg/DQgEtjb5theINt7/E3AgHe90aCf60HpDwJXeY+vAh4+zvl57WrIavN9iXOuts33WR3U3Pb87hoH3Op1Z5aZWRmBz6o37yXSYwo8kfd7iEDL7irgX865oo5Ocs5td85dCWQAvwCeMLN4AmEysYOXHORoK65FDlDgvd8B59z1zrksAt1+d3VyL6yjLU6eBmaZ2UwC9xwfOc7POLZdDYVdvH9hBzUX0nN5wM+ccylt/sQ55x7txXuJ9JgCT+T9HgIWA9fTeXcmZnaVmaU755qBMu9wM4GwWWxml5tZhDc37STnXBOB+3w/M7NE717f1wjcg8PMLjOzMd77lBIInuYOLl1E4B5gK69F9gSB1ui7zrl9x/kZv2hmY8wslcA9use6OPdR4DYzSzezNOD7LTX30B+AG83sdAuIN7MLzSyxF+8l0mMKPJF2vHtrbwLxwLNdnHoBsNHMKgkMYLnCOVfjhc2HgVuBwwQGrMz2XnMzUAXsAl4nEFB/9J47FXjHe79nga90Mp/tfmC61y34dJvjDwIncvzuTLzr/turYyfQ1by+nwIrgXXAemD1cc7vkHNuJYF/RNxJINB3ANf09H1Eesu0AazI8GBmOcAWYJRz7kgX5+0BPu+c+89g1SYSDNTCExkGvLl8XwP+2lXYiYSy462iICJBzhsoU0Rg9OQFPpcjErTUpSkiIiFBXZoiIhIShnSXZlpamhs/frzfZYiISBBZtWrVQedcevvjQzrwxo8fz8qVK/0uQ0REgoiZ7e3ouLo0RUQkJCjwREQkJCjwREQkJCjwREQkJCjwREQkJCjwREQkJCjwREQkJCjwREQkJCjwREQkJCjwREQkJCjwREQkJCjwREQkJCjwREQkJCjwREQkJCjwREQkJCjwREQkJCjwREQkJCjwREQkJCjwREQkJCjwREQkJCjwREQkJCjwREQkJCjwREQkJCjwREQkJCjwREQkJCjwREQkJCjwREQkJCjwREQkJCjwREQkJCjwREQkJARV4JnZBWa21cx2mNm3/a5HRESGj6AJPDMLB/4P+BAwHbjSzKb7W1XvOef8LkFERNqI8LuANk4DdjjndgGY2V+Bi4FNvlbVS5fe9SZhBtOzkpg+OpnpWUlMzUwkNirc79JEREJSMAVeNpDX5vt84PT2J5nZEmAJQE5OzuBU1kPOOeaOG8G6gnKeWVPIn9/eB0CYwYS0eKZnJTN9dJIXhkmkJ0b7XLGIyPAXTIHXLc65e4F7AebOnRuU/YZmxm0XBXpjnXPkl9awsfAIm/YfYVPhEVbvLeW59wpbz09PjD4mAKdnJTF+ZDzhYebXjyAiMuwEU+AVAGPbfD/GOzakmRljU+MYmxrHBTNHtR4vq65n0/4jbN5fwSYvDN94dReNzYEMj40M54TRiccE4QmjktQlKiLSSxYsgyvMLALYBiwiEHQrgE855zZ29pq5c+e6lStXDlKFA6+usYkdxZWtAdjytaK2ETjaJTp7bAoXzRrNOVMzMFMrUESkLTNb5Zyb2/540LTwnHONZvYl4F9AOPDHrsJuOIqOCGdGVjIzspJbj7V0ibYNwKVbinlqdQGnjh/BDz4yg5nZyV28q4iIQBC18HpjuLXwuquhqZknV+XzP//ayuHqeq48LYevf3AqqfFRfpcmIuK7zlp4QTMPT7ovMjyMK07LYenXF/K5sybw2Io8Fv7PMv70xm4am5r9Lk9EJCgp8Iaw5NhI/t9F03nxKx9g1pgUfvjcJi787eu8ueOg36WJiAQdBd4wMDkzkYevO417rj6FqvpGPnXfO9z0yCryS6v9Lk1EJGgo8IYJM+P8GaP4z9cWcOt5U1i6pZhFt7/Cb17aRk19k9/liYj4ToE3zMREhnPzosksvXUhH5wxijte3s7iX7/CC+v3a31PEQlpCrxhKisllt9dOYfHlswjMSaCmx5ZzZV/eJstB474XZqIiC8UeMPc6bkj+cfNZ/OTS2ay5UAFF/72dX72/CYq6xr9Lk1EZFAp8EJARHgYV88bx7JbF3L53LH84bXdLL79FZ5fp25OEQkdCrwQMiI+iv/+2Ik8ddOZpMZH8cW/rOYzf3yX3Qer/C5NRGTAKfBC0Mk5I3j2S2fxw49MZ+2+Ms7/zav8+qVt1DZoNKeIDF8KvBAVER7GNWdN4OVbF/ChE0fx25e388HfvMqyrcV+lyYiMiAUeCEuIymGO66Yw18+fzoR4ca1D6zgxodXUVhW43dpIiL9SoEnAJw5KY0XvzKfb5w/leXbApPWf//KThq0NqeIDBMKPGkVFRHGF8+ZxEu3LOCsSWn8/J9b+PAdr/H2rkN+lyYi0mcKPHmfsalx3PfZudz3mbnUNDRxxb1v87XH1lJSUed3aSIivabAk04tnp7JS7cs4EvnTOK5dYWce/tyHnxzj7YgEpEhSYEnXYqNCufr50/lxa/OZ/aYFH7w7EY+eucbrNpb6ndpIiI9osCTbpmYnsDD153GnZ+aw6GqOj5+95t884n3OFSpbk4RGRoUeNJtZsZFs7J4+daFLJmfy1OrCzj39ld45J29NDVriTIRCW4KPOmxhOgIvvvhabzwlQ8wbXQi3/v7Bi696w3eyyvzuzQRkU4p8KTXpmQm8uj187jjipPYX17LJXe9wXf/vp7Sqnq/SxMReR8FnvSJmXHxSdksvXUB1545gcdW5HHu7ct5bMU+mtXNKSJBRIEn/SIxJpLvf2Q6/7j5bCZlJPCtJ9fz8d+/yYaCcr9LExEBFHjSz6aNTuLxG87g9stmk3e4mo/e+Trff2YD5dUNfpcmIiFOgSf9zsz4+CljePnWhVw9bxx/fnsv596+nCdW5WvDWRHxjQJPBkxybCQ/ungmz37pbHJGxvH1v73H5fe8xeb9R/wuTURCkAJPBtzM7GSevPFMfvnxWewsqeKi373Oj5/bxJFadXOKyOBR4MmgCAszLj91LEtvXcAVp47lgTd3s+j2V3h6TYG6OUVkUCjwZFClxEXxs0tP5JkvnkVWcgxffWwtV9z7NtuKKvwuTUSGOQWe+GLWmBSeuuks/uvSE9lyoIIP3/Ea//XCZirrGv0uTUSGKQWe+CY8zPjU6Tks+/pCPnHKGO59dReLbl/Oc+8VqptTRPqdAk98lxofxc8/PounbjqTtIRobn50DVff/y47iiv9Lk1EhhEFngSNk3NG8OyXzubHF8/gvfwyPnTHq/zixS1U16ubU0T6ToEnQSU8zPjMGeNZ9vWFXHxSNncv38ni21/hxQ371c0pIn2iwJOglJYQza8um83fbjyDpNhIbvzzaj77wAp2H6zyuzQRGaIUeBLUTh2fyj9uPpvvXzSd1XtLOf83r3L7v7dSU9/kd2kiMsQo8CToRYSH8bmzJ7D01gVcOGs0v1u6g/N+8wovbSryuzQRGUIUeDJkZCTF8JtPnsRfl8wjLiqc6x9ayXV/WkFhWY3fpYnIEKDAkyFnXu5Inv/yB/jeh6fx5s5DfPA3r/Lou/s0qEVEuqTAkyEpMjyM6+fn8q+vzmfWmGS+89R6rrr/HfIOV/tdmogEKQWeDGk5I+N45POn81+Xnsh7eeV88Dev8qc3dtPcrNaeiBxLgSdDnllgibJ/3TKf0yak8sPnNvHJe99iV4lWahGRoxR4Mmxkp8Typ2tP5VeXzWbrgQo+dMdr3PfaLt3bExFAgSfDjJlx0azRfHnRZOoam/np85spqazzuywRCQIRfhcg0l82Fpbz+Io8nl5bSHlNA2NGxHLd2RNIT4j2uzQRCQIKPBnSyqsbeOa9Ah5bkcfGwiNERYRxwYxRfPLUsZyRO5KwMPO7RBEJEgo8GXKamx1v7zrEYyvzeHHDAeoam5k+OokffXQGF5+URUpclN8likgQUuDJkFFYVsOTq/J5fFUeeYdrSIyJ4PK5Y/nkqWOZmZ3sd3kiEuQUeBLU6hub+c/mIh5bkcdr20todnDmxJF8/YNTOX/GKGIiw/0uUUSGCAWeBJXahia2HqhgQ2E5K3Yf5pVtJZRWNzAqKYYvnjOJy04ZS87IOL/LFJEhSIEnvqmsa2Tz/iNsKChnQ8ERNhaWs724kiZvlZS0hCgWTEnn4jnZzJ+cTrgGoIhIHyjwZFCUVdezsdALt8IjbCwoZ/ehKlrmhKclRHNidhKLp2UyMzuJGVnJjBkRi5lCTkT6hwJP+l1xRS0bC1rCrZyNhUfILz26hU92SiwzspK4ZE42M7OTmJmVTEZSjI8Vi0goUOBJrznnKCirae2O3FAQCLfiiqMrm0xIi+eksSlcNW8cM7OSmZGVxIh4TRsQkcGnwJNuaW527DlU1dodubHwCBsKyymrbgAgzGByRiJnT05jZlYyM7OTmTY6kcSYSJ8rFxEJUODJ+zQ2NbOjpLK15bbR+1pV3wRAVHgYU0cl8qGZo5jhhdsJoxI1RUBEgpoCTwAor2lg+dZiXtpUxCtbS6ioawQgNjKc6VlJfOKUMczISmZGdhKTMxKJitC64yIytCjwQlje4Wr+s7mIlzYV8e7uwzQ2O9ISovjwiaOZNzGVE7OTmZCWoOkAIjIsKPBCSHOzY0NhOS9tCoTclgMVAEzOSOD6+bmcNz2Tk8akaMFlERmWFHjDXG1DE2/tOsRLm4p4eXMRRUfqCDOYOz6V2y6cxqJpmUxIi/e7TBGRATfogWdmY4GHgEzAAfc65+4ws1TgMWA8sAe43DlXOtj1DQelVfUs3VLMfzYX8eq2Eqrqm4iLCmfBlHQWT8vk3BMyNDVAREKOHy28RuBW59xqM0sEVpnZS8A1wMvOuZ+b2beBbwPf8qG+IWnPwSr+s7mIf28qYuWewzQ7yEiM5uI52Zw3PZMzckdqFKWIhLRBDzzn3H5gv/e4wsw2A9nAxcBC77QHgeUo8LpUXFHLw2/t5Z8bDrCjuBKAE0Yl8sVzJnHe9ExmZiXrfpyIiMfXe3hmNh6YA7wDZHphCHCAQJdnR69ZAiwByMnJGYQqg09BWQ33vrKTR1fk0djUzLzckXz69BwWT8tkbKp2EhAR6YhvgWdmCcCTwFedc0faLhLsnHNm5jp6nXPuXuBegLlz53Z4znC152AVdy/fyVNr8nEOPnZyNl9YOEmDTkREusGXwDOzSAJh94hz7invcJGZjXbO7Tez0UCxH7UFo+1FFfzfsh08+14hEeFhXHlaDkvm5zJmhFpzIiLd5ccoTQPuBzY7537d5qlngc8CP/e+PjPYtQWbDQXl3Ll0By9uPEBcVDif/0Aunz97gnYWEBHpBT9aeGcBVwPrzWytd+y7BILucTO7DtgLXO5DbUFh1d5S7ly6nWVbS0iMjuDmcydx7VkTSNVUAhGRXvNjlObrQGdDBxcNZi3BxDnHW7sOcefSHby58xAj4iL5xvlTufqMcSRpxwERkT7TSitB4EhtA995aj3Pr9tPRmI0t104jU+dnkNclP7ziIj0F/1G9dnavDJufnQ1+8tq+cb5U7nu7AmaIC4iMgAUeD5pbnbc//pufvHiFjKTYnj8xjM4OWeE32WJiAxbCjwfHKqs49a/vcfyrSVcMGMUv/j4LJLjdJ9ORGQgKfAG2Vs7D/HVx9ZQWt3ATy6ewVXzxtF20r2IiAwMBd4gaWp23PHydn63dDsTRsbzx2tOZUZWst9liYiEDAXeICivbmDJwyt5Z/dhPnZyNj+5eCbx0froRUQGk37rDrCa+iY+9+AK1ueXc/tls/n4KWP8LklEJCQp8AZQQ1MzX3hkFWv2lfJ/nzqZD5042u+SRERClgJvgDQ3O77hjcT874+dqLATEfFZmN8FDEfOOX7y/CaeXlvIN86fypWnhea+fSIiwUSBNwD+b9kOHnhjD587awI3LZzodzkiIoICr9898s5efvXvbVw6J5vbLpymOXYiIkFCgdePXt5cxG1Pb+Ccqen88hOzCAtT2ImIBAsFXj8pq67nW0+uZ9qoJO769ClEhuujFREJJhql2U/+64XNlFbX8+DnTiU2SrsdiIgEGzVD+sEbOw7y+Mp8lszP1XJhIiJBSoHXRzX1TXznqfVMSIvnK4sm+12OiIh0Ql2affTrl7ay73A1f10yTxu3iogEMbXw+mBdfhn3v76bK0/LYV7uSL/LERGRLijw+uCn/9hMWkI03/nwCX6XIiIix6HA66XV+0p5d89hblwwkaQY7VYuIhLsFHi99IdXd5EcG8knTx3rdykiItINCrxe2HOwihc3HuCqeTnayFVEZIhQ4PXC/a/vJjIsjM+eMd7vUkREpJsUeD10uKqev63K49I52WQkxfhdjoiIdJMCr4ceemsPtQ3NXD9/gt+liIhIDyjweqC52fHXd/M4Z2o6kzIS/S5HRER6QIHXA6v3lXLgSC2XzMn2uxQREekhBV4PPL9+P1ERYSyalul3KSIi0kMKvG5qbnb8c/0BFkxJJ0FTEUREhhwFXjetyQt0Z1544mi/SxERkV5Q4HXT8+sOEBUexrnTMvwuRUREekGB1w3NzY5/btjP/ClpWjdTRGSIUuB1w+YDR9hfXsv5M0b5XYqIiPSSAq8b3tp5CICzJ6f5XImIiPSWAq8b3thxkNy0eEYnx/pdioiI9JIC7zgampp5d/dhzpykHc1FRIYyTSg7jnX55VTVNzEiLoqtBypIjIkgISaC+KgIwsPM7/JERKSbFHjHERMZRkJ0BL9buoPfLd1xzHPxUeEkxESQEB1BQkwkSS2PowOhmOh9TYiObA3Ko8ciSIyOJD46nIhwNbRFRAaaAu84ZmQls+J7i1m1t5TymgYq6xqoqG2koraRyrpGKr2vFXWNVNY2cKC89ujx+kacO/414qLCOwjJo0GZGHM0RBOiW76PbH2cGBNBfHQEkQpOEZFOKfC6ITYqvFcjNJubHdUNTV4oNnCk9mhAVta2hGTj0RBtE6AHK6oDQVrbQGVdI83dCM5AazSyNSDbBmXb1mZCTEQHrdHI1kCNilBwisjwo8AbQGFh1hoq0PvNYp1zVNc3eQHYtmXZLkSPeT4Qovuqqo8eq2ukqRvJGRURdjQQ27U2W4O0fYi2aW22vC46IrzXP7OISH9T4A0BZkZ8dKDbMjOp9+/jnKO2oZkKr0XZ2h3bJiRbvm/b2qysbaSgrKa1JVpZ20hjd4IzPKxNYEYc27KMOTZEjw3SyGO6b6MjwjDTACER6RsFXggxM2KjwomNCqcv+9c656hrbD6mtVlR1xD4+r7WZkNrcB6pbaSwrLbN8w00NB0/OCPD7diQbNvy9FqVia2hGvn+7lwvRGMiFZwioUyBJz1mZsREhhMTGU56YnSf3quusemYoDwamA2t9znf3xptoLiill0lR1uj9Y3Nx71WeJsu5mO7Z9uF5DGt0ch2g4UiiIsKV3CKDEEKPPFVdEQ40QnhjEzoe3BW1TV12NpsOziofYgerKxnz6Hq1iCtbTh+cIYZxEdFeN3M4V5ghhMfFQjElu7nhOhw72v7497rojW6VmQwKfBkWIiOCCc6IpzU+Kg+vU9DU/P772223Lts6b71HlfVNVJV30hlXROVtQ0crKgPHK8PPNed7tpA7WEdBmV8dKCrtqMAbQnclikpLQEaG6nWp0hnFHgibUSGhzEiPooRfQxOONrqrPLCMxCOXlDWBYLy6OOjz1XWNXKosp59h6rbBGtTt67ZtvWZENMmKDtofSYc0+JU61OGPwWeyADpr1YnBOZ0BlqOTe1Cs20rs82xNser6ho5VFl9TKj2tfX5/m7aY1ufLfc91fqUYHLcwDOzXwI/BWqAF4FZwC3OuT8PcG0i4gkLMxJjIknspw2I27Y+KztsZR7b+mw5VlnX0OfWZ0JMu6Bs1/pMjIkgPiq809anluST3upOC++DzrlvmtmlwB7gY8CrgAJPZIjq79ZndUObrtsOWp/HdN3Wdtz6bHm+L63Pjrtpj22VnjgmmYzE3i8EIUNXdwKv5ZwLgb8558rVLSEiLdquKNSXhRFadNb6bL0f2r5Lt02L9FBVPXsPV7d273bU+owKD+PSOdlcPz+XSRkJfS9YhozuBN4/zGwLgS7NL5hZOlA7sGWJSKgaqNZnZV0jZdUNPL2mgMdX5vHYyjzOm57JjQtyOWVcaj9ULsHOXDeW8zezVKDcOddkZnFAknPuwIBXdxxz5851K1eu9LsMERliDlXW8dBbe3nwrT2UVTcwd9wIblgwkUUnZBCmfS6HPDNb5Zyb+77jxws8Mwsn0J05njYtQufcr/u5xh5T4IlIX1TXN/K3lfn84bVd5JfWMDE9niXzc7lkTrYWPx/C+hJ4LxDowlwPtC5D4Zz7UX8X2VMKPBHpD41Nzbyw4QD3vLKTjYVHyEiM5tqzJvDpeTkk9dPIWBk8fQm8dc65WQNWWR8o8ESkPznneGPHIe55dSevbT9IQnQEnzo9h8+dNYFRyRrZOVR0FnjdGbTyTzP7oHPu3wNQl4hI0DAzzp6cxtmT09hQUM69r+7ivtd28cAbu7n4pGyWzM9lSmYfthoRX3WnhXcpgTl3YUADYIBzzvXDAOS+UQtPRAZa3uFq7n99N39dsY/ahmYWnZDBDQsmcur4EVo5Jkj1pUtzN3AxsN51Z0hn9wsKB1YCBc65i8xsAvBXYCSwCrjaOVff1Xso8ERksByuqudhb2Tn4ap65uSkcMP8iXxweqZGdgaZzgKvO2vz5AEb+jPsPF8BNrf5/hfAb5xzk4BS4Lp+vp6ISK+lxkfxlcWTeeNb5/KTS2ZyqLKeG/+8isW/eYXHV+Z1a09G8Vd3Wnh/AnKBfwJ1Lcf7Mi3BzMYADwI/A74GfAQoAUY55xrN7Azgh86587t6H7XwRMQvTc2Of27Yz13LdrJp/xGykmO4fn4uV5yaQ2yUpjT4qS8tvN3Ay0AUkOj96et6PP8LfJOj0xxGAmXOuUbv+3wgu6MXmtkSM1tpZitLSkr6WIaISO+EhxkXzcri+S+fzQPXnsqYEXH86LlNnPWLpfzu5e2UVzf4XaK0051Rmg8553a3PWBmp/b2gmZ2EVDsnFtlZgt7+nrn3L3AvRBo4fW2DhGR/mBmnDM1g3OmZrByz2HuWr6T21/axj2v7uLT83K47uwJWqw6SHQn8J4ws4865woAzGw+8H/Aib285lnAR83sw0AMkATcAaSYWYTXyhsDFPTy/UVEfDF3fCp/vCaVTYVHuPuVnfzh1V088MYeLjtlDDfMn0jOyDi/Swxp3bmHdypwF4H7bCcD/w1c5JzL6/PFAy28r3ujNP8GPOmc+6uZ/R5Y55y7q6vX6x6eiASzPQeruOfVXTy5Kp8m5/jIrNF8YeEkpo7SXL6B1OtpCd6LzwDuIbDE2IXOuX65edYu8HIJTEtIBdYAVznn6rp4uQJPRIaEoiO13PfaLh55Zx/V9U0snpbBTedM4uScEX6XNiz1OPDM7Dmg7ZPTgf0EpgzgnPvoANTZIwo8ERlKSqvqefCtPfzpzcAuDfNyU7lp4SQ+MDlNk9j7UW8Cb0FXb+ice6Wfaus1BZ6IDEVVdY08+u4+7nttNweO1HJidjI3LZzI+TNGaRJ7P+hTl2awUuCJyFBW19jE31cX8PtXdrLnUDW56fHcuGAil5yUTVREd2aNSUcUeCIiQUqT2PuXAk9EJMg551i+rYS7l+3k3T2HSY2P4tozx/OZM8aTHKd9+bpLgSciMoSs2HOYu5btYNnWEhKiIzSJvQf6slvCWcAPgXEEJqq3bA+UOwB19ogCT0SGu5ZJ7M+vKyQiPEyT2LuhL4G3BbiFwJY9TS3HnXOH+rvInlLgiUioCExi38mTqwo0if04+hJ47zjnTh+wyvpAgScioeb9k9gzuemciZrE3kZfAu/nQDjwFMduD7S6v4vsKQWeiIQqTWLvXF8Cb1kHh51z7tz+Kq63FHgiEupaJrH/4bVdFB2p0yR2NEpTRGRY0yT2o3qztNhVzrk/m9nXOnq+Lzue9xcFnojIsZqaHS+s389dy3eyOUQnsXcWeF3thxfvfdUQIBGRISI8zPjI7CwumjW6dRL7j57bxO+W7uBzZ43n6jPGkxwbmpPY1aUpIjLMhdokdt3DExEJce0nsV8+NzCJfWzq8JrErsATERFg+E9iV+CJiMgxDpTXcv/rw28Se1/m4UUDHwfG02aQi3Pux/1cY48p8ERE+m64TWLvS+C9CJTz/rU0b+/vIntKgSci0n+GyyT2vgTeBufczAGrrA8UeCIi/W+oT2LvLPC6U/mbZnbiANQkIiJBKDoinCtOy+HlWxfyuyvnEB0RzjefWMfC/1nGA2/spqa+6fhvEoS6WmllPeAI3LebDOwisHh0y354swaryM6ohSciMvA62ok9mCex92ZpsXFdvaFzbm8/1dZrCjwRkcE1FCax9+Ue3sPOuauPd8wPCjwREX9sLCzn7uU7eWH9/qCbxN6btTRbzGj3RuHAKf1VmIiIDD0zspK581Mnt05if3xFPo++mxfUk9i76tL8DvBdIBaoJnDvDqAeuNc5951BqbALauGJiASHA+WBndj/8q7/k9j70qX538EQbh1R4ImIBJdgmMTel8Az4FLgbAKjNl9zzj09EEX2lAJPRCQ4+TmJvS+BdxcwCXjUO/RJYKdz7ov9XmUPKfBERIJb+0nss8Yk88drTiUtIXrArtmXwNsCTHPeiWYWBmx0zk0bkEp7QIEnIjI0NDU7nllbwHf/vp6s5Fj+/PnTyUqJHZBr9WWllR1ATpvvx3rHREREuiU8zPjYyWN4+LrTKamo47Lfv8Weg1WDWkN3Ai8R2Gxmy81sGbAJSDKzZ83s2YEtT0REhpNTx6fy6JJ5VNc3cvk9b7G9qGLQrt2dLs0FXT3vnHulXyvqAXVpiogMTduKKvj0fe/Q1Ox4+LrTmJGV3G/v3esuTS/Q9gCR3uN3gdXOuVf8DDsRERm6pmQm8vgNZxATEcaV977Nmn2lA37N4waemV0PPAHc4x0aAzw9gDWJiEgIyEyK5meXnkhdYzNX3fcOq/YeHtDrdWdpsS8CpwHvADjntptZxoBWJSIiQ5pzjvKaBgrKaigsq6WwrIYC709+aQ35h6s5VFXfen4dsGxLCaeMSx2wmroTeHXOufqWGfJmFkFgArqIiISohqZmio7UUlBaQ2F5INTyS2soLDv6p6rdvnlREWFkp8SSnRLLB2dkMmZEHGNT4xg7IpaxqXGMjI8a0Jq7E3ivmNl3gVgzOw+4CXhuQKsSERFfHaltaA2ugtIaCtq00grLaig6Uktzu6bPyPgoslJiyU2P5wOT08lKiSE7JZaslFiyR8QyMj5q0JYX60h3Au/bwHXAeuAG4AXgvoEsSkREBk5Ts6PoSNsAq6WgrPpo12NpDRV1jce8JjLcGJ0caJ2dOTGN7JQYskcEwiwrJZas5Fhio8J9+om657iB55xrNrOngaedcyUDX5KIiPRFVV0jhWU15Je17WKs9VpqNRw4UktTu+ZZSlwkWcmxjBkRx7zckWSlxARaZt6ftIToAV8Dc6B1GnjeotE/AL6EN5rTzJqA3znnfjw45YmISFvNzY6SyrrWrsUC775Z2y7H8pqGY14TEWaMSg4E2GkTUr2uxrhjuhzjo7vT4Te0dfUT3gKcBZzqnNsNYGa5wN1mdotz7jeDUaCISChoGdVYUlFHSWVd4GtFHQcr6ymuqG1tpe0vr6Gh6djWWWJMRGtL7JRxI7xuxhjGeF2OGYkxhA/x1ll/6CrwrgbOc84dbDngnNtlZlcB/wYUeCIiXXDOUVXf1Ca86t7/uLKOg97X9kEGgXtn6QnRZKXEctLYFD584miyR8SS7XU5ZqXEkhQT6cNPN/R0FXiRbcOuhXOuxMz06YpIyKptaDqmJdZZkJVU1FHb0Py+14cZjEyIJj0hmvTEaCZnJJKeGHiclhBFemI0GYnRpCfEkBQb4evIxuGkq8Cr7+VzIiJDTn1jM4eq6jhYUU9JZW1rd2JLkLW2xCrq3jeCsUVqfFRrYJ2SM4I0L9COhlng64i4KHUx+qCrwJttZkc6OG5AzADVIyLSb5qaHYer6jvtQmzbKiutbujwPRJjIlrDalpWEvNbQizh2CAbmRBFZHh3NqARv3QaeM654J5QISIh6ZjBHW2Dq7KldXY0yA5X1b1vcjRAbGR4a1jlpsdzem4q6Qkxx3QptgRZTKR+FQ4Xw38cqogMCbUNTewvr+1wgEdJ5dHvD3YyuCMqPKw1rLJTYjhpbPLRLsWEaNLatMpCYQi+vJ/+q4vIoKisa/QmPldTUOotIOzNI8svreFgZd37XhNmkJZw9N7XlMzEY+6FBQIsSoM7pFsUeCLSZy3djPleeBW0Bll16wr5Ze3ukUWFhwUmPo+IZdEJGa3LVGUmHe1O1OAO6U8KPBE5Lue81T1Kj27vcvRxoMXWfmX82MhwxowILBp80tgUxoyI8+aPxTJ2xPBYqkqGFgWeiLQuJpzfrsuxtaVWVkN947HzyZJiIsgeEce4kfGcOTGNMSNiAwGXEgi2EXGR6mKUoKLAEwkB9Y3N7C+v6eDeWaDL8UB5LY3thjOmJUSRnRLLCaMTWTQtI9BC87Z5yR6h1T1k6FHgiQwDNfVNFJRVt7bK2nc5FlfU4drkmRlkJgbun50ybkRrkLWGWkrwb/Ui0lMKPJEhoq6xiT0Hq9lWVMH2ogp2llSR590/O1R17OJHEWHGaG8l/LMnpbfeSxvjBdvo5FiiIjRJWkKLAk8kyNQ1NrH7YBXbiirZUVTBtqJKthdXsOdQdeseZmEGY1PjyEmNY0ZW0vtaaJlJWh1fpD0FnohP6hqb2FVSxfbiSrYXVbC9qJJtxRXsbRds40bGMzkjgQtmjmJKZiKTMhKYmJ6gFUBEekiBJzLAahtaWmwV7Ciu9LokK9l7+NhgGz8ynsmZCXx45mgmZyYwOSOR3PR4BZtIP1HgifST2oaWFlsg0Fq+7jlU1bqeY3iYMW5kHJMzE7hw1mgmZSQwJTORCWkKNpGBpsAT6aG2wdbSWtteXMneDoJtSmYiF80azaTMRKZkJjAhLZ7oCAWbiB98CTwzSwHuA2YCDvgcsBV4DBgP7AEud86V+lGfCEBjUzPbiyvZeqDCC7dKdnQQbONHxnHCqEQ+Mms0kzMTmaxgEwlKfrXw7gBedM59wsyigDjgu8DLzrmfm9m3gW8D3/KpPglBxUdqWb2vjDV5pazZV8b6/HJqGgLLZUWEGePT4pk2OpGPzM5iinePbUJavIb3iwwRgx54ZpYMzAeuAXDO1QP1ZnYxsNA77UFgOQo8GSC1DU1sLCxnzb4y1uSVsXZfGQVlNQBEhhszspL55KljmZOTwrTRSYwfqWATGer8aOFNAEqAB8xsNrAK+AqQ6Zzb751zAMjs6MVmtgRYApCTkzPw1cqQ55wj73BNa8ttzb5SNu0/0rqnWnZKLHNyUvjc2ROYk5PC9NFJGkAiMgz5EXgRwMnAzc65d8zsDgLdl62cc87MOtinGJxz9wL3AsydO7fDcyS0VdQ2sC6/nDX7SlmbV8aafWWtK5HERYUza0wyn/9ALnPGpnBSTgoZiTE+Vywig8GPwMsH8p1z73jfP0Eg8IrMbLRzbr+ZjQaKfahNhpjmZseOkkrW7GtpvZWxrbiidd3IienxnHNCBnNyUpgzdgRTMhOICFfXpEgoGvTAc84dMLM8M5vqnNsKLAI2eX8+C/zc+/rMYNcmwe9QZV1rq21NXinr8sqpqGsEIDk2kjk5KXz4xNHMyUlh9pgUkuO0or+IBPg1SvNm4BFvhOYu4FogDHjczK4D9gKX+1SbBJHymgZe217C8q0lrNhzmL2HqoHAdIBpoxO5ZE42J41NYU5OChPS4rX/moh0ypfAc86tBeZ28NSiQS5Fgoxzjq1FFSzbUsKyrcWs2ltKU7MjOTaSebmpfOq0HObkjODE7GRtXyMiPaKVVsR3VXWNvLHjIMu2lrB8azH7y2sBmJGVxBcWTOScE9KZPSZF995EpE8UeDLonHPsOljFsi3FLN9awru7D1Pf1ExCdAQfmJzGLYszWDA1ncwkjZ4Ukf6jwJNBUdvQxNu7DrFsSzHLtpaw73DgXtzkjASuOWs8C6emM3dcqiZ3i8iAUeDJgMk7XM3yrYGAe3PnQWobmomJDOOsiWlcPz+XhVPSGZsa53eZIhIiFHjSb5qaHe/sPtqK21FcCUBOahxXnJrDOSdkcPqEVK1iIiK+UOBJn+0oruCJVQU8vaaAA0dqiQoP4/TcVK48LYdzpqZruoCIBAUFnvRKaVU9z60r5MlV+byXX054mLFgSjq3XTSNc6ZmEB+t/7VEJLjot5J0W0NTM8u2FPPk6nyWbimmockxbXQSt104jYtPyiY9MdrvEkVEOqXAky4559hYeIQnVuXz7HuFHK6qJy0his+cMZ6PnzyG6VlJfpcoItItCjzpUPGRWp5eW8CTqwrYWlRBVHgYi6dn8PGTxzB/SjqRmgQuIkOMAk9a1TY08e9NRTy1Op9Xt5XQ7GBOTgo/vWQmF80aTUpclN8lioj0mgJPOFBeyz2v7uSJVflU1DaSlRzDFxZO5GMnj2FieoLf5YmI9AsFXgjLL63m7uU7+dvKfJqd46JZo7ls7ljOyB1JWJimEYjI8KLAC0F7DlZx1/IdPLW6ADO4bO5YvrBgolY9EZFhTYEXQnYUV3Dn0h08+14hkeFhXDVvHDcsyGV0cqzfpYmIDDgFXgjYvP8Idy7dwQsb9hMTEc7nP5DL5z8wgYxE7UYgIqFDgTeMrc8v57dLt/PSpiISoiO4aeFErjs7l9R4jbYUkdCjwBuG1uwr5Y6Xt7N8awlJMRF8dfFkrj1zAslxkX6XJiLiGwXeMFJWXc/P/7mFv67IIzU+im+cP5XPnDGOxBgFnYiIAm8YcM7xzNpCfvKPTZTVNHDD/Fy+vGiyFnAWEWlDvxGHuL2Hqrjt6Q28tv0gs8em8PClJ2p9SxGRDijwhqj6xmb+8NoufvvydqLCw/jJxTP41OnjCNeEcRGRDinwhqAVew7z3afWs724kgtPHM33PzKdzCRNMRAR6YoCbwgpr27g5y9u5tF388hOieX+z85l0bRMv8sSERkSFHhDxGvbS7jlsbWUVjewZH4uX108mbgo/ecTEeku/cYMcs457nttN//9z81Mzkjkwc+dxoysZL/LEhEZchR4Qay2oYnvPLWev68p4EMzR/Gry2ZrqoGISC/pt2eQKiyr4YaHV7G+oJxbz5vCF8+ZpC17RET6QIEXhFbsOcwX/ryK2oZm/vCZuZw3XQNTRET6SoEXZB55Zy8/fHYjY0bE8dclpzApI9HvkkREhgUFXpCob2zmh89t5C/v7GPBlHR+e+UckmO1BqaISH9R4AWBmvomrntwBW/uPMSNCybyjfOnasUUEZF+psDzWW1DE59/aAVv7TrEry6bzSdOGeN3SSIiw1KY3wWEstqGJpY8vIo3dx7ifz6hsBMRGUgKPJ/UNTZx0yOreXVbCb/42CyFnYjIAFPg+aC+sZkvPrKGpVuK+dmlM7n81LF+lyQiMuwp8AZZQ1MzX350Df/ZXMSPL57Bp08f53dJIiIhQYE3iJqaHbc8tpYXNx7g+xdN5zNnjPe7JBGRkKHAG0R3/Gcb/1i3n29/6AQ+d/YEv8sREQkpCrxB8vLmIn67dAeXnTKGG+bn+l2OiEjIUeANgj0Hq/jqY2uZmZ3ETy6ZiZkmlYuIDDYF3gCrqW/ixj+vIjzMuPvTpxATGe53SSIiIUkrrQwg5xzf/ft6thZV8MA1pzI2Nc7vkkREQpZaeAPo4bf38vc1BdyyeAoLp2b4XY6ISEhT4A2QjYXl/OQfm1h0QgZfOmeS3+WIiIQ8Bd4AaG52fO/vG0iKieT2y2drp3IRkSCgwBsAj67Yx9q8Mr534TRS4qL8LkdERFDg9buDlXX84p9bmJebyqVzsv0uR0REPAq8fvZfL2ympqGJn2q+nYhIUFHg9aO3dh7iqdUFLJmfy6SMRL/LERGRNhR4/aS+sZn/98wGxoyI5UvnTPa7HBERaUcTz/vJQ2/tYUdxJX+8Zi6xUVpNRUQk2KiF1w9qG5r4/Su7OCN3JOeekOl3OSIi0gEFXj/4yzv7OFhZx1cWqytTRCRYKfD6KNC628npE1KZlzvS73JERKQTCrw+emxFHsUVat2JiAQ7BV4f1DU2cffynZw2PpUz1LoTEQlqCrw+eHxlPgeO1PKVxZM1yVxEJMgp8HrJOceDb+5h9phkzpyo1p2ISLBT4PXSqr2l7Ciu5NOnj1PrTkRkCFDg9dJf3t1HQnQEF80e7XcpIiLSDb4EnpndYmYbzWyDmT1qZjFmNsHM3jGzHWb2mJkF7b465dUNPL9uPxeflEVclBarEREZCgY98MwsG/gyMNc5NxMIB64AfgH8xjk3CSgFrhvs2rrr6bUF1DU2c+VpOX6XIiIi3eRXl2YEEGtmEUAcsB84F3jCe/5B4BJ/Suuac45H393HzOwkZmYn+12OiIh006AHnnOuAPgVsI9A0JUDq4Ay51yjd1o+0OHuqWa2xMxWmtnKkpKSwSj5GO/ll7PlQIVadyIiQ4wfXZojgIuBCUAWEA9c0N3XO+fudc7Ndc7NTU9PH6AqO/f8ukIiw42PzM4a9GuLiEjv+dGluRjY7Zwrcc41AE8BZwEpXhcnwBigwIfauuSc418bizhrUhpJMZF+lyMiIj3gR+DtA+aZWZwFJrAtAjYBy4BPeOd8FnjGh9q6tHl/BfsOV3P+jFF+lyIiIj3kxz28dwgMTlkNrPdquBf4FvA1M9sBjATuH+zajudfGw9gBudN1553IiJDjS+TyJxzPwB+0O7wLuA0H8rptn9tPMCp41JJS4j2uxQREekhrbTSTXsPVbHlQAXnz1R3pojIUKTA66Z/bywC4IPqzhQRGZIUeN30+o6DTM5IYGxqnN+liIhILyjwuqGhqZkVew5zhrYBEhEZshR43bAuv4zq+ibteyciMoQp8LrhzR2HMIPTJyjwRESGKgVeN7y16xDTRiUxIj5odywSEZHjUOAdR21DEyv3lur+nYjIEKfAO441+8qob2zW/TsRkSFOgXccq/eVAjB3XKrPlYiISF8o8I5jbV4ZuWnxJMdpdwQRkaFMgXcc7+WVMXtsit9liIhIHynwunCgvJbiijpmj0n2uxQREekjBV4X1uaVATBLLTwRkSFPgdeF9/LLiAgzpo9O8rsUERHpIwVeFzYUlDN1VCIxkeF+lyIiIn2kwOvC1gMVnDBKrTsRkeFAgdeJsup6iivqmDoqwe9SRESkHyjwOrGtqBKAyZmJPlciIiL9QYHXia1FFQBMUeCJiAwLCrxObC+qICE6gqzkGL9LERGRfqDA68S2ogomZyZgZn6XIiIi/UCB14ltRZVMVXemiMiwocDrQFl1PYer6pmYrhGaIiLDhQKvA3sOVQMwbmScz5WIiEh/UeB1YO+hKgDGjYz3uRIREekvCrwO7PNaeDmpauGJiAwXCrwO7D1cTWZSNLFRWkNTRGS4UOB1YO+hKsalqjtTRGQ4UeB1YO+hag1YEREZZhR47dQ2NFFcUaf7dyIiw4wCr5395bUAjE6J9bkSERHpTwq8dvaX1QBoDU0RkWFGgddOoVp4IiLDkgKvnZYW3mi18EREhhUFXjuF5bWkxkcRE6k5eCIiw4kCr5395TVq3YmIDEMKvHYOlNcyOln370REhhsFXjtFR2rJTIr2uwwREelnCrw2GpqaKa1uID1RgSciMtwo8No4VFkPoMATERmGFHhtHKysAyAtQYEnIjLcKPDaKFHgiYgMWwq8Ng5WBAIvQ12aIiLDjgKvDbXwRESGLwVeG6VV9cRGhmuncxGRYUiB18bhqgZS46P8LkNERAaAAq+N0up6RsRH+l2GiIgMAAVeG4er6hkRpxaeiMhwpMBro7S6Xl2aIiLDlAKvDbXwRESGLwWep6GpmYraRgWeiMgwpcDzlFYH1tFM1aAVEZFhSYHnKa1qAGCE7uGJiAxLCjzP4SqvhacuTRGRYUmB5zlUFVhWLDVBgSciMhwp8DxFRwKBl5kY43MlIiIyEBR4nuIjtUSFh5ESp0ErIiLDkQLPU3SklvTEaMzM71JERGQAKPA86wrKmZKZ4HcZIiIyQAYs8Mzsj2ZWbGYb2hxLNbOXzGy793WEd9zM7LdmtsPM1pnZyQNVV0eKK2rZVVLF6bkjB/OyIiIyiAayhfcn4IJ2x74NvOycmwy87H0P8CFgsvdnCXD3ANb1Pu/uPgzA6RNSB/OyIiIyiAYs8JxzrwKH2x2+GHjQe/wgcEmb4w+5gLeBFDMbPVC1tffu7sPERYUzMzt5sC4pIiKDLGKQr5fpnNvvPT4AZHqPs4G8Nufle8f2046ZLSHQCgSoNLOt/VBXGnAw6if98E7DXxpw0O8ihhB9Xj2jz6tn9Hl1bFxHBwc78Fo555yZuV687l7g3v6sxcxWOufm9ud7Dlf6rHpGn1fP6PPqGX1ePTPYozSLWroqva/F3vECYGyb88Z4x0RERPrFYAfes8BnvcefBZ5pc/wz3mjNeUB5m65PERGRPhuwLk0zexRYCKSZWT7wA+DnwONmdh2wF7jcO/0F4MPADqAauHag6upEv3aRDnP6rHpGn1fP6PPqGX1ePWDO9fg2moiIyJCjlVZERCQkKPBERCQkhHTgmdkFZrbVW9Ls28d/Regys7FmtszMNpnZRjP7it81BTszCzezNWb2D79rCXZmlmJmT5jZFjPbbGZn+F1TMDOzW7y/hxvM7FEz075m3RCygWdm4cD/EVjWbDpwpZlN97eqoNYI3Oqcmw7MA76oz+u4vgJs9ruIIeIO4EXn3AnAbPS5dcrMsoEvA3OdczOBcOAKf6saGkI28IDTgB3OuV3OuXrgrwSWOJMOOOf2O+dWe48rCPxCyva3quBlZmOAC4H7/K4l2JlZMjAfuB/AOVfvnCvztajgFwHEmlkEEAcU+lzPkBDKgdfZcmZyHGY2HpgDvONzKcHsf4FvAs0+1zEUTABKgAe8LuD7zCze76KClXOuAPgVsI/A8ovlzrl/+1vV0BDKgSe9YGYJwJPAV51zR/yuJxiZ2UVAsXNuld+1DBERwMnA3c65OUAVR3dSkXa8bdUuJvAPhSwg3syu8reqoSGUA0/LmfWQmUUSCLtHnHNP+V1PEDsL+KiZ7SHQVX6umf3Z35KCWj6Q75xr6TF4gkAASscWA7udcyXOuQbgKeBMn2saEkI58FYAk81sgplFEbjp+6zPNQUtMzMC91g2O+d+7Xc9wcw59x3n3Bjn3HgC/18tdc7pX+CdcM4dAPLMbKp3aBGwyceSgt0+YJ6ZxXl/LxehQT7d4ttuCX5zzjWa2ZeAfxEY5fRH59xGn8sKZmcBVwPrzWytd+y7zrkX/CtJhpGbgUe8f3zuYvCXFxwynHPvmNkTwGoCo6fXoCXGukVLi4mISEgI5S5NEREJIQo8EREJCQo8EREJCQo8EREJCQo8EREZVGb2RzMrNrMN3Tj3a96i9evM7GUzG+cdP8nM3vIW0V5nZp883nsp8ES6wdsp4vx2x75qZneb2UcHercNM7vRzD7jPb7GzLJ6+Pr7erLYt3eNEjNb6/2yub6nNYt04U/ABd08dw2BhbJnEViU4Jfe8WrgM865Gd57/a+ZpXT1RpqWININZrYEOMM5d22bY28D33TOvTrItSwHvu6cWzmA17iGwC+ZL5lZBrARmOmcKxqoa0po8dbk/Ye34wNmNpHADjbpBMLseufclnavmQPc6Zw7q4P3ew/4hHNue2fXVAtPpHueAC70Jka3/GXNAl7zWkN3escv8/Yoe8/MXvWOhZvZr7zj68zsZu/4Im+x5PVeF0+0d/znbbpwfuUd+6GZfd3MPgHMJTBJe62ZXWhmT7cUaWbnmdnf2xdvZsvNbK73uNLMfubV+LaZZXb1gzvnioGdwDivRbvS60b6UUfnd/QZiHTDvcDNzrlTgK8Dd3VwznXAP9sfNLPTgCgC/592KmRXWhHpCefcYTN7l8D+ic8QWDLsceecC6zu1Or7wPnOuYI23StLgPHASd4KP6nehp1/AhY557aZ2UPAF8zsYeBS4ATvvVPavrlz7glvhaCvO+dWektL3W5m6c65EgIrlPzxOD9OPPC2c+57ZvZL4Hrgp52dbGa5QC6wA/ie91mEAy+b2Szn3Lp2L+noMxDplLco/ZnA39r8fYpud85VBP6xt6Dd8dHAw8BnnXNd7k6iFp5I9z3K0Y02r/C+b+8N4E/ePa9w79hi4B7nXCMEwhOYSmAB4G3eOQ8S2BOuHKgF7jezjxHo2umUC9yTeBi4yguXM+jgX8Dt1AMtu7CvIhDGHfmkt4zco8ANXt2Xm9lqAvdVZhDYPLm9jj4Dka6EAWXOuZPa/JnW8qSZLQa+B3zUOVfX5ngS8DyBf4i93Z2LiEj3PAMsMrOTgbiOtv9xzt0I3EZgJ45VZjayJxfwQvE0Al2oFwEvduNlDwBXAVcCf2sJ1i40uKM375vovKfnMe8Xz+nOub+b2QQCXU2LvAEEzwMxHfwMffoMJPR4W43tNrPLILBYvZnN9h7PAe4hEHbFLa/xbi/8HXjIOfdEd66jwBPpJudcJbCMQJdhR607zGyic+4d59z3CWxqOhZ4CbjBArtTY2apwFZgvJlN8l56NfCK17WT7C3KfQswu4PLVACJbeoqJLDj9W0Ewm+gJBHYq67cu+/3oY5O6uQzEGllZo8CbwFTzSzfzK4DPg1c5w0+2Uhgzz+A/wESCHR3rjWzll1tLifQK3KNd3ytmZ3U1XV1D0+kZx4l8K/KKzp5/n/MbDJgwMvAe8AGYAqwzswagD845+40s2sJ/CWOILBd1e+BVOAZ7x6fAV/r4Bp/An5vZjUERo7WAI8A6c65Adsmxjn3npmtAbYAeQS6LgEwsx8DK51zz9LxZyDSyjl3ZSdPvW+qgnNucSfv8WegR/tMalqCyDDgjRJd45y73+9aRIKVAk9kiDOzVQS6Gs9re0NfRI6lwBMRkZCgQSsiIhISFHgiIhISFHgiIhISFHgiIhISFHgiIhIS/j+j/Ny15CqAlwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Only run this when in serial. Will fail in parallel\n", "\n", "if GEO.nProcs == 1:\n", "\n", " distances, viscosities = GEO.extract_profile(Model.projViscosityField, \n", " line = [(180.* u.kilometer, 0.), (180.* u.kilometer, Model.bottom)])\n", " distances, stresses = GEO.extract_profile(Model.projStressField, \n", " line = [(180.* u.kilometer, 0.), (180.* u.kilometer, Model.bottom)])\n", "\n", "\n", " Fig, (ax1) = plt.subplots(1,1,figsize=(7,7))\n", " ax1.plot(GEO.dimensionalise(viscosities, u.pascal * u.second), GEO.dimensionalise(distances, u.kilometer))\n", " ax1.set_xlabel(\"Viscosity in Pa.s\")\n", " ax1.set_ylabel(\"Depth in kms\")\n", " ax1.set_ylim(100,-10)\n", " ax1.set_title(\"Viscosity profile\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solver options" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [], "source": [ "Model.solver.set_inner_method(\"mumps\")\n", "Model.solver.set_penalty(1e6)\n", "GEO.rcParams[\"initial.nonlinear.tolerance\"] = 1e-4" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Run the Model" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Model.run_for(nstep=2, checkpoint_interval=1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "analysis" ] }, "outputs": [], "source": [ "Fig = vis.Figure(figsize=(1200,400), title=\"Viscosity Field (Pa.s)\", quality=3)\n", "Fig.Points(Model.swarm, \n", " GEO.dimensionalise(Model.viscosityField, u.pascal * u.second),\n", " logScale=True,\n", " fn_size=3.0)\n", "Fig.VectorArrows(Model.mesh, Model.velocityField)\n", "Fig.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "analysis" ] }, "outputs": [], "source": [ "Fig = vis.Figure(figsize=(1200,400), title=\"Viscosity Field (Pa.s)\", quality=3)\n", "Fig.Surface(Model.mesh, \n", " GEO.dimensionalise(Model.strainRateField, 1.0 / u.second),\n", " logScale=True)\n", "Fig.show()" ] } ], "metadata": { "kernelspec": { "display_name": "virtualenv", "language": "python", "name": "virtualenv" }, "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.6" } }, "nbformat": 4, "nbformat_minor": 2 }