From 03508042c092ad1ef84816419d4c1d388cc76870 Mon Sep 17 00:00:00 2001 From: Ferran Brosa Planella Date: Thu, 27 Jun 2024 22:02:43 +0100 Subject: [PATCH] #3844 remove old notebooks --- .../old_notebooks/2-a-pde-model.ipynb | 382 ---------- .../3-negative-particle-problem.ipynb | 400 ---------- ...paring-full-and-reduced-order-models.ipynb | 474 ------------ .../old_notebooks/5-half-cell-model.ipynb | 693 ------------------ .../old_notebooks/6-a-simple-SEI-model.ipynb | 681 ----------------- .../old_notebooks/7-creating-a-submodel.ipynb | 486 ------------ 6 files changed, 3116 deletions(-) delete mode 100644 docs/source/examples/notebooks/creating_models/old_notebooks/2-a-pde-model.ipynb delete mode 100644 docs/source/examples/notebooks/creating_models/old_notebooks/3-negative-particle-problem.ipynb delete mode 100644 docs/source/examples/notebooks/creating_models/old_notebooks/4-comparing-full-and-reduced-order-models.ipynb delete mode 100644 docs/source/examples/notebooks/creating_models/old_notebooks/5-half-cell-model.ipynb delete mode 100644 docs/source/examples/notebooks/creating_models/old_notebooks/6-a-simple-SEI-model.ipynb delete mode 100644 docs/source/examples/notebooks/creating_models/old_notebooks/7-creating-a-submodel.ipynb diff --git a/docs/source/examples/notebooks/creating_models/old_notebooks/2-a-pde-model.ipynb b/docs/source/examples/notebooks/creating_models/old_notebooks/2-a-pde-model.ipynb deleted file mode 100644 index 4926d19432..0000000000 --- a/docs/source/examples/notebooks/creating_models/old_notebooks/2-a-pde-model.ipynb +++ /dev/null @@ -1,382 +0,0 @@ -{ - "cells": [ - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Creating a simple PDE model" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "In the [previous notebook](./1-an-ode-model.ipynb) we show how to create, discretise and solve an ODE model in pybamm. In this notebook we show how to create and solve a PDE problem, which will require meshing of the spatial domain.\n", - "\n", - "As an example, we consider the problem of linear diffusion on a unit sphere,\n", - "$$\n", - " \\frac{\\partial c}{\\partial t} = \\nabla \\cdot (\\nabla c),\n", - "$$\n", - "with the following boundary and initial conditions:\n", - "$$\n", - " \\left.\\frac{\\partial c}{\\partial r}\\right\\vert_{r=0} = 0, \\quad \\left.\\frac{\\partial c}{\\partial r}\\right\\vert_{r=1} = 2, \\quad \\left.c\\right\\vert_{t=0} = 1.\n", - "$$\n", - "\n", - "As before, we begin by importing the PyBaMM library into this notebook, along with any other packages we require:\n" - ] - }, - { - "cell_type": "code", - "execution_count": 1, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Note: you may need to restart the kernel to use updated packages.\n" - ] - } - ], - "source": [ - "%pip install \"pybamm[plot,cite]\" -q # install PyBaMM if it is not installed\n", - "import pybamm\n", - "import numpy as np\n", - "import matplotlib.pyplot as plt" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Setting up the model" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "As in the previous example, we start with a `pybamm.BaseModel` object and define our model variables. Since we are now solving a PDE we need to tell pybamm the domain each variable belongs to so that it can be discretised in space in the correct way. This is done by passing the keyword argument `domain`, and in this example we choose the domain \"negative particle\"." - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": {}, - "outputs": [], - "source": [ - "model = pybamm.BaseModel()\n", - "\n", - "c = pybamm.Variable(\"Concentration\", domain=\"negative particle\")" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Note that we have given our variable the (useful) name \"Concentration\", but the symbol representing this variable is simply `c`." - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We then state out governing equations. Sometime it is useful to define intermediate quantities in order to express the governing equations more easily. In this example we define the flux, then define the rhs to be minus the divergence of the flux. The equation is then added to the dictionary `model.rhs`" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "metadata": {}, - "outputs": [], - "source": [ - "N = -pybamm.grad(c) # define the flux\n", - "dcdt = -pybamm.div(N) # define the rhs equation\n", - "\n", - "model.rhs = {c: dcdt} # add the equation to rhs dictionary" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Unlike ODE models, PDE models require both initial and boundary conditions. Similar to initial conditions, boundary conditions can be added using the dictionary `model.boundary_conditions`. Boundary conditions for each variable are provided as a dictionary of the form `{side: (value, type)`, where, in 1D, side can be \"left\" or \"right\", value is the value of the boundary conditions, and type is the type of boundary condition (at present, this can be \"Dirichlet\" or \"Neumann\")." - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": {}, - "outputs": [], - "source": [ - "# initial conditions\n", - "model.initial_conditions = {c: pybamm.Scalar(1)}\n", - "\n", - "# boundary conditions\n", - "lbc = pybamm.Scalar(0)\n", - "rbc = pybamm.Scalar(2)\n", - "model.boundary_conditions = {c: {\"left\": (lbc, \"Neumann\"), \"right\": (rbc, \"Neumann\")}}" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Note that in our example the boundary conditions take constant values, but the value can be any valid pybamm expression.\n", - "\n", - "Finally, we add any variables of interest to the dictionary `model.variables`" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, - "outputs": [], - "source": [ - "model.variables = {\"Concentration\": c, \"Flux\": N}" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Using the model" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now the model is now completely defined all that remains is to discretise and solve. Since this model is a PDE we need to define the geometry on which it will be solved, and choose how to mesh the geometry and discretise in space." - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Defining a geometry and mesh\n", - "\n", - "We can define spatial variables in a similar way to how we defined model variables, providing a domain and a coordinate system. The geometry on which we wish to solve the model is defined using a nested dictionary. The first key is the domain name (here \"negative particle\") and the entry is a dictionary giving the limits of the domain." - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": {}, - "outputs": [], - "source": [ - "# define geometry\n", - "r = pybamm.SpatialVariable(\n", - " \"r\", domain=[\"negative particle\"], coord_sys=\"spherical polar\"\n", - ")\n", - "geometry = {\n", - " \"negative particle\": {r: {\"min\": pybamm.Scalar(0), \"max\": pybamm.Scalar(1)}}\n", - "}" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We then create a uniform one-dimensional mesh with 20 points. " - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": {}, - "outputs": [], - "source": [ - "# mesh and discretise\n", - "submesh_types = {\"negative particle\": pybamm.Uniform1DSubMesh}\n", - "var_pts = {r: 20}\n", - "mesh = pybamm.Mesh(geometry, submesh_types, var_pts)" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Example of meshes that do require parameters include the `pybamm.Exponential1DSubMesh` which clusters points close to one or both boundaries using an exponential rule. It takes a parameter which sets how closely the points are clustered together, and also lets the users select the side on which more points should be clustered. For example, to create a mesh with more nodes clustered to the right (i.e. the surface in the particle problem), using a stretch factor of 2, we pass an instance of the exponential submesh class and a dictionary of parameters into the `MeshGenerator` class as follows: `pybamm.MeshGenerator(pybamm.Exponential1DSubMesh, submesh_params={\"side\": \"right\", \"stretch\": 2})`" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "After defining a mesh we choose a spatial method. Here we choose the Finite Volume Method. We then set up a discretisation by passing the mesh and spatial methods to the class `pybamm.Discretisation`. The model is then processed, turning the variables into (slices of) a statevector, spatial variables into vector and spatial operators into matrix-vector multiplications." - ] - }, - { - "cell_type": "code", - "execution_count": 8, - "metadata": {}, - "outputs": [], - "source": [ - "spatial_methods = {\"negative particle\": pybamm.FiniteVolume()}\n", - "disc = pybamm.Discretisation(mesh, spatial_methods)\n", - "disc.process_model(model);" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now that the model has been discretised we are ready to solve. " - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Solving the model" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "As before, we choose a solver and times at which we want the solution returned. We then solve, extract the variables we are interested in, and plot the result." - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "2021-11-19 15:31:50,774 - [WARNING] processed_variable.get_spatial_scale(520): No length scale set for negative particle. Using default of 1 [m].\n" - ] - }, - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6AAAAEYCAYAAABCw5uAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABQiUlEQVR4nO3dd3yV9fn/8deVRQKEGXYIewjIDOCeKO7RuveurbXaVm3t/nVacVstRdzbumqtioiCoLIFBCFhhD3DCBlknuv3R45+KWWcQM65T5L38/HIg5z73Oect7eQT677/tzXx9wdERERERERkWhLCDqAiIiIiIiINAwqQEVERERERCQmVICKiIiIiIhITKgAFRERERERkZhQASoiIiIiIiIxkRR0gN1lZGR4165dg44hIiIN2Jw5c/LdvU3QOQ6WxlIREYkH+xpP46oA7dq1K7Nnzw46hoiINGBmtiroDIdCY6mIiMSDfY2nmoIrIiIiIiIiMaECVERERERERGJCBaiIiIiIiIjEhApQERERERERiQkVoCIiIiIiIhITKkBFREREREQkJlSAioiIBMTMEs3sSzN7dz/7DDezKjO7IJbZREREokEFqIiISHBuAxbv60kzSwT+CkyIWSIREZEoUgEqIiJ1VlllFU9/lsfv3lkUdJQaM7NM4Exg/H52uxV4A9gck1AiItIgbSks4+qnZrJmW0nUPysp6p8gIiJSy6pCzptz1/LQR0tZt2MXR/VoTXlliJSkOnVe9SHgLiB9b0+aWSfgfOAkYPj+3sjMbgJuAsjKyqrVkCIiUr8VlFRw5ZMzWLW1hC1FZXRu1Tiqn6cCVERE6gx354OFG7l/Yi7LNhcxMLM593z3cI7pmYGZBR0vYmZ2FrDZ3eeY2Qn72O0h4GfuXnWg/zZ3HweMA8jOzvbaSyoiIvVZcVkl1z4zkxVbinnymmyGZrWM+meqABURkbjn7kxdms+YCTl8ta6Anm2bMvaKoYzu375OFZ67ORo4x8zOAFKBZmb2grtfsds+2cAr4f++DOAMM6t097djnlZEROqd0ooqvvf8HOat2cHjlw/j2F5tYvK5KkBFRCSuzVm1nTETljB9xTY6tUjjvgsHcf6QTiQm1MnCEwB3vxu4GyB8BfSOPYpP3L3bN9+b2TPAuyo+RUSkNlRWhbjtlS+Ztiyf+y4cxGkD2sfss1WAiohIXFqycSf3Tcjho8WbyWiawu/O7selI7NolJQYdLSoMbObAdx9bNBZRESkfgqFnLveWMCERZv43dn9uGBYZkw/XwWoiIjElVVbi3lgYi7vzF9P00ZJ3Dm6D9ce3ZXGKfVzyHL3ycDk8Pd7LTzd/ZrYJRIRkfrK3fn9u1/z5tx1/PSU3lxzdLcDv6iW1c/RXERE6pxNO0t5ZNJSXp21hqRE4+bje3DzcT1o3jg56GgiIiL1woMTc3nm85XceGw3fnhSz0AyqAAVEZFAbS8uZ+yU5Tzz+UqqQs6lI7K49aSetG2WGnQ0ERGRemP81BU88vEyLs7uzC/OOCywJn4qQEVEJBDFZZU8OS2PJz5dQVF5JecP7sTto3qT1Tq664+JiIg0NK/MXM0f/7OYMw/vwJ+/c3igHeRVgIqISEyVVVbx4vTVPPbJMrYWl3Nqv3b89NQ+9GmfHnQ0ERGReufdBeu5+62vOL53Gx68eHDgXeSjVoCaWR/g1d02dQd+4+4PReszRUQkflVWhXhz7joe+iiX9QWlHNWjNXeO7sOQGCx6LSIi0hB9krOZH786j+wuLRl7xTBSkhKCjhS9AtTdc4DBAGaWCKwD3orW54mISHwKhZwPFm3kvg9zWLGlmEGZzbn3gkEc0ysj6GgiIiL11sy8bdz8/Bz6tE/nyWuGk5YSH8uYxWoK7snAcndfFaPPExGRgLk7ny7NZ8yEJSxct5NebZsy9ophjO7fLtB7T0REROq7r9YWcN0zs8hsmcaz146gWWr8dJSPVQF6CfDy3p4ws5uAmwCysrJiFEdERKJpzqpt3PtBDjPytpHZMo37LxzEeUM6BX7fiYiISH23bHMhVz89k+Zpybxww0haN20UdKT/EvUC1MxSgHOAu/f2vLuPA8YBZGdne7TziIhI9CzesJP7JuQwaclmMpo24v+d059LRnSmUVJ8TPsRERGpz9ZsK+GK8TNJMOPFG0bSoXla0JH+RyyugJ4OzHX3TTH4LBERCcCqrcU8MDGXd+avJ71REneO7sO1R3elcYqarYuIiMTC5p2lXPHkDHZVVPHq946ga0aToCPtVSx+M7iUfUy/FRGRum1jQSmPfLyU12atITkxge8f34PvHdeD5o3j514TERGR+m5HSTlXPjmTLYVlvHDDSPq2bxZ0pH2KagFqZo2BU4DvRfNzREQktrYVl/P3yct47otVhNy5fGQWt5zUk7bpqUFHExERaVCKyiq5+ulZ5OUX8/S1wxka58ubRbUAdfcSoHU0P0NERGKnqKySJ6fm8cTUFZSUV3L+kExuH9WLzq0aBx1NRESkwSmtqOKm52azcF0Bf798KEf3jP8lznRzjoiIHFBpRRUvTF/F45OXs624nNH92/HTU/vQu1160NFEREQapIqqED986Us+X76VBy8exKn92wcdKSIqQEVEZJ8qq0K8PmctD09ayoaCUo7pmcGdo/swqHOLoKOJiIg0WKGQc+c/5/PR4k38/tz+nD8kM+hIEVMBKiIi/yMUct5buIEHPsxlRX4xgzu34P4LB3FUHZjaIyIiUp+5O799ZxFvz1vPHaf25qojuwYdqUZUgIqIyLfcncm5W7hvQg6L1u+kd7umjLtyGKf0a4eZBR1PRESkwbvvwxyen76K7x3XnVtO7Bl0nBpTASoiIgDMXrmNez/IYebKbXRulcYDFw3i3MGdSExQ4SkiIhIPxk5ZzmOfLOfSEZ35+el96+TJYRWgIiIN3KL1Bdw3IYdPcrbQJr0Rfzi3PxcPzyIlKSHoaCIiIhL20ozV3PP+Es4a2IE/nnd4nSw+QQWoiEiDlZdfzAMTc/n3/PU0T0vmZ6f15eqjutA4RUODiIhIPHln/np++fZXnNS3LQ9ePLhOz07SbxkiIg3MhoJdPDJpKa/NXktKYgK3nNiDm47rQfO05KCjiYiIyB4mLd7ET16dx/CurXj88qEkJ9btGUoqQEVEGohtxeU8/skynpu+CnfnyiO6cMuJPWmT3ijoaCIiIrIXXyzfyg9enMthHZrx5NXZpCYnBh3pkKkAFRGp54rKKhk/dQXjp+ZRUl7J+UMyuX1ULzq3ahx0NBEREdmHeWt2cONzs+ncqjHPXjeC9NT6MVNJBaiISD1VWlHFC9NX8fjk5WwrLuf0Ae35ySm96dUuPehoIiIish8L1xVw1ZMzaNUkhReuH0mrJilBR6o1KkBFROqZyqoQ/5yzlkcmLWVDQSnH9srgztF9GJjZIuhoIiIicgA5Gwu58skZpKcm89KNI2nfPDXoSLVKBaiISD0RCjnvfrWBByfmkpdfzJCsFtx/0SCO6pERdDQRERGJwPItRVw+fgbJiQm8eMNIMlvWv9tlVICKiNRx7s4nOZsZMyGXxRt20rd9OuOvyubkw9rW2TXCREREGpqV+cVc/sQMwHnpxiPpmtEk6EhRoQJURKQOm5m3jTETljBr5XayWjXmoYsHc/agjnV6fTAREZGGZtXWYi59YjpllVW8dOMR9GzbNOhIUaMCVESkDlq4roAxE3KYkruFtumN+ON5A7h4eOc6vzaYiIhIQ7N6awmXjpvOrooqXrrhCA7r0CzoSFGlAlREpA5ZvqWIBybm8p8FG2ielszPT+/L1Ud2JS2l7q8LJiIi0tCs2VbCpU9Mp7i8ihdvGEm/jvW7+AQVoCIidcL6Hbt4+KOlvD53LY2SErj1pJ7ceFx3mtWTNcFEREQamrXbq4vPwtIKXrrxCAZ0ah50pJhQASoiEse2FpXx2CfLeWH6KgCuOrILPzihJ23SGwWcTERERA7W+h27uPSJ6RTsquDFG0Y2mOITVICKiMSlwtIKnpiax5NTV7CroooLhmXyo5N71ct27CIiIg3JhoJdXDJuOjuKK3jhhpENbp1uFaAiInGktKKK579YxeOTl7G9pIIzDm/PT07pTc+26UFHExERkUO0saCUS8dNZ1txOc9fP4JBnVsEHSnmVICKiMSBiqoQ/5y9lkcmLWXjzlKO692GO0/tw+GZDWdKTkNkZonAbGCdu5+1x3OXAz8LPywCvu/u82McUUREasmmnaVc9sR08ovKefa6EQzJahl0pECoABURCVAo5Px7wXoenJjLyq0lDM1qwYMXD+bIHq2DjiaxcRuwGNhb28M84Hh3325mpwPjgJGxDCciIrVj885SLn1iOpt2lvLc9SMY1qVhFp+gAlREJBDuzsdLNjNmQg5LNhbSt306T16dzUl922JmQceTGDCzTOBM4E/AT/Z83t0/3+3hdCAzRtFERKQWbS6sLj43FpTy7HUjGNalVdCRAqUCVEQkxqav2MqYCTnMWbWdLq0b8/Algzl7YEcSElR4NjAPAXcBkdzgez3w/r6eNLObgJsAsrKyaiObiIjUgi2FZVz+xAzW7yjlmWuHM7xrwy4+IcoFqJm1AMYDAwAHrnP3L6L5mSIi8eqrtQWM+TCHT3O30K5ZI/58/uFcmJ1JcmJC0NEkxszsLGCzu88xsxMOsO+JVBegx+xrH3cfR/UUXbKzs732koqIyMHKLyrj8vHTWbO9hKevGcHI7rq9BqJ/BfRh4AN3v8DMUgCtHyAiDc6yzUU8MDGH977aSMvGyfzijL5cdWRXUpMTg44mwTkaOMfMzgBSgWZm9oK7X7H7TmY2kOoTuae7+9YAcoqIyEHYVlzOFeNnsHpbCU9dPVy9HXYTtQLUzJoBxwHXALh7OVAerc8TEYk363bs4qGJubwxdy1pyYn86ORe3HhsN9JTk4OOJgFz97uBuwHCV0Dv2EvxmQW8CVzp7rmxzigiIgdne3E5lz0xnbz8Yp66ZjhH9cwIOlJcieYV0O7AFuBpMxsEzAFuc/fi3XfSfSsiUt/kF5Xx2CfLeHH6ajC45qhu3HJiD1o3bRR0NIlzZnYzgLuPBX4DtAYeDzemqnT37ADjiYjIAewoKefy8TPIyy/myauHc7SKz/9h7tG5VcTMsqnu2ne0u88ws4eBne7+6329Jjs722fPnh2VPCIi0baztIInPl3Bk9PyKK2o4sJhnbltVC86tkgLOprUgJnNqcuFnsZSEZFgbC0q44onZ7J8SxHjr8rmuN5tgo4UqH2Np9G8AroWWOvuM8KPXwd+HsXPExEJxK7yKp77YiV/n7KcHSUVnHl4B35yam96tGkadDQRERGJgc2FpVwxfgartpbw5NXZHNurYRef+xO1AtTdN5rZGjPr4+45wMnA19H6PBGRWKuoCvHqrDU8MmkpmwvLOL53G+4c3YcBnZoHHU1ERERiZGNBKZeNn86GHaU8fe1wjuqhabf7E+0uuLcCL4Y74K4Aro3y54mIRF0o5Lwzfz0PTMxl9bYSsru05NFLh6i9umBmrdx9W9A5REQkNtbt2MVlT0wnv7CM564foXU+IxDVAtTd5wF19j4aEZHduTuTFm/mvg9zWLKxkMM6NOOpa7I5sU9bwk1ipAExs1+5+x/D3/cD3gaSrfovw8W73YIiIiL1UF5+MVeMn8HO0gqev2EkQ7NaBh2pToj2FVARkXrhi+VbGTNhCXNX76Br68Y8fMlgzh7YkYQEFZ4N2HeAP4a/H0N1p/f3zWwE8BBwVFDBREQkunI3FXL5+BlUVoV4+cYjdPtNDRywADWzRsB3ga677+/uv49eLBGR+PDV2gLunbCEqUvzad8slb9853AuGJZJcmJC0NEkvnR09/cB3H2mman1sYhIPfXV2gKuemoGyYkJvPa9I+nVLj3oSHVKJFdA/wUUUL2OZ1l044iIxIdlmwu5/8Nc3l+4kZaNk/nlGYdx5ZFdSE1ODDqaxI/uZvYOYECmmTV295Lwc8kB5hIRkSiZs2ob1zw1i2Zpybx4w0i6ZjQJOlKdE0kBmunup0U9iYhIHFi7vYSHP1rKG3PXkpacyG0n9+KGY7uRnqp6Qv7HuXs8TgAws3bA32MfR0REoumzZfnc8Oxs2jdP5YUbRtJJ63wflEgK0M/N7HB3/yrqaUREArKlsIzHPlnGSzNWg8F1R3fjByf2pFWTlKCjSZxy9yn72L4JeCzGcUREJIomLd7E91+cS7fWTXj+hhG0TU8NOlKdFUkBegxwjZnlUT0F1wB394FRTSYiEgMFuyp44tMVPPVZHmWVIS7KzuRHJ/eiQ3Od1ZSDZ2Y3ufu4oHOIiMih+8+CDdz2ypcc1qEZz103gpY6OX1IIilAT496ChGRGNtVXsUzn69k7JTlFOyq4KyBHfjJKb3p3qZp0NGkflB7ZBGReuD1OWu56/X5DOvSkievGU4z3ZJzyA5YgLr7KjMbBBwb3jTV3edHN5aISHSUV4Z4ddZqHvl4GVsKyzixTxvuGN2H/h3VPl1qj7v/I+gMIiJyaJ7/YiW//tcijumZwbirhtE4RStY1oZIlmG5DbgReDO86QUzG+fuj0Y1mYhILaoKOf+at44HP8plzbZdDO/akscvH8rwrq2CjiZ1mJn1pboZUSfAgfXAO+6+ONBgIiJySP4xZTl/eX8Jow5ry98uG6ou+LUokjL+emCkuxcDmNlfgS8AFaAiEvfcnYlfb+L+D3PJ2VRIvw7NePraAZzQuw1mmiUpB8/MfgZcCrwCzAxvzgReNrNX3P2ewMKJiMhBcXce+mgpD09aylkDO/DgxYO19ncti6QANaBqt8dV6N4WEakDPl+Wz70Tcpi3ZgfdMprw6KVDOPPwDiQk6EeY1Irrgf7uXrH7RjN7AFgEqAAVEalD3J0//Wcx46flceGwTO757kAS9TtDrYukAH0amGFmb4Ufnwc8GbVEIiKHaP6aHYyZkMO0Zfl0aJ7KPd85nAuGZZKkM5hSu0JAR2DVHts7hJ8TEZE6IhRyfvWvhbw0YzVXH9mF357dXyesoySSJkQPmNlkqpdjMeBad/8y2sFERGpq6aZC7vswhwmLNtGqSQq/OvMwrjiii+7bkGi5HZhkZkuBNeFtWUBP4IdBhRIRkZqprApx1+sLePPLdXz/hB7cNbqPbtOJon0WoGbWzN13mlkrYGX465vnWrn7tujHExE5sDXbSnjoo6W89eVaGqck8eNRvbn+2G40baRudRI97v6BmfUGRlDdhMiAtcAsd6/a74tFRCQulFZUcdsrXzJh0SbuOLU3PzypV9CR6r39/Xb2EnAWMIfqzn7fsPDj7lHMJSJyQFsKy/jbx0t5aeZqzIzrj+nG90/oSSstEC0x4u4hYHrQOUREpOYKSyu46bk5fLFiK785qx/XHdMt6EgNwj4LUHc/K/yn/k+ISFwp2FXBuE+X89S0lZRXhbgouzM/OrknHZqnBR1NGigzm+bux3zzZ9B5RERk/7YUlnHN0zPJ2VjIQxcP5rwhnYKO1GBEsg7oJHc/+UDbRESiraS8kmc+X8nYycvZWVrJ2YM68pNTetMto0nQ0UQah//UX0YRkTi3ZlsJVz45g407S3ni6mxO7NM26EgNyv7uAU2lekDNMLOW/N/SK82o7vonIhIT5ZUhXpm1mkc/XsaWwjJO7tuWn57ah34dmwUdTUREROqQResLuObpWZRXhnjxhiMY1qVl0JEanP1dAf0e1R3+OlJ9H+g3BehO4LHoxhIRgaqQ8/aX63jwo1zWbt/FiG6t+PvlQ8nu2iroaCIiIlLHfL4sn5uen0Oz1CReuvlIerVLDzpSg7S/e0AfBh42s1vd/dEYZhKRBs7dmbBoE/d/mMPSzUUM6NSMP51/OMf1ylBbdBEREamxf89fz09em0f3jKY8c91w9Y0IUCTrgD5qZgOAfkDqbtufi2YwEWmYpi3NZ8yEJcxfW0D3Nk147LKhnD6gvRaDlninv6AiInHqyWl5/PE/XzO8ayueuCqb5mnJQUdq0CJpQvRb4ASqC9D3gNOBaYAKUBGpNV+u3s6YCTl8vnwrHZuncu93B/KdoZ1ISkwIOppIJH68x58iIhKwUMj56wdL+MenKzitf3seumQwqcmJQcdq8CJZpf0CYBDwpbtfa2btgPHRjSUiDUXOxkLu+zCHiV9vonWTFH5zVj8uG5mlAULqDDO70N3/CeDuk/fcJiIisVdeGeKu1+fz9rz1XHVkF357dn8SNZsqLkRSgO5y95CZVZpZM2Az0D3KuUSknluzrYQHJ+by1rx1NE1J4qen9ObaY7rRtFEkP5ZE4srdwJ7F5t62iYhIDBSVVfL9F+YwdWk+d47uww9O6KEeEnEkkt/0ZptZC+AJqrvhFgEzI3lzM1sJFAJVQKW7Zx9cTBGpLzbvLOVvnyzj5ZmrSTDjpmO7c/PxPWjZJCXoaCI1YmanA2cAnczskd2eagZUBpNKRKRh21JYxrXPzGTxhkLuvWAgF2V3DjqS7GG/BahVnyr4i7vvAMaa2QdAM3dfUIPPONHd8w8ho4jUAwUlFYz9dDlPf5ZHZZVz8fDO3HpSL9o3Tz3wi0Xi03pgNnAO1Sdov1GI7gUVEYm5lfnFXPXUTDYXlvLEVcM4qW+7oCPJXuy3AHV3N7O3gWHhxytjkElE6pGS8kqe/mwlY6csp6isknMGdeTHo3rTNaNJ0NFEDom7zwfmm9lL7l4RdB4RkYZswdodXPv0LELuvHzjEQzJahl0JNmHSKbgTjez4e4+6yDe34EPzcyBf7j7uD13MLObgJsAsrKyDuIjRCQelVVW8crMNTz68TLyi8oYdVhbfnpqHw7r0CzoaCK1SsWniEiwpuRu4fsvzKFVkxSevW4EPdo0DTqS7EckBeiJwPfMbBVQTPVaZ+7uAyN47dHuvt7M2gITzWyJu3+6+w7honQcQHZ2ttcsvojEm6qQ89aX63hwYi7rduxiZLdW/OPKoQzr0iroaCIiIlLPvDl3LXe9voBe7dJ59trhtG2mW3viXSQF6OkH++buvj7852YzewsYAXy6/1eJSF3k7kxYtIn7P8xh6eYiDu/UnL9853CO7ZWhznNSr+1tyZVIl2Exs0Sq7yNd5+5n7fGcAQ9T3eioBLjG3efWXnIRkbrL3fnHpyu45/0lHNWjNWOvHEaz1OSgY0kEIilA/+juV+6+wcyeB67cx/7f7NMESHD3wvD3pwK/P+ikIhKX3J1py/IZMyGHBWsL6N6mCY9fPpTTB7RX4SkNxaEsw3IbsJjqzrl7Oh3oFf4aCfw9/KeISINWFXL+8O7XPPP5Ss4e1JH7LhxIoyStH15XRFKA9t/9Qfhs7bAIXtcOeCv8C2gS8JK7f1DjhCISt+au3s6YD3L4YsVWOrVI494LBvKdIZ1ISkwIOppI1B3qMixmlgmcCfwJ+MledjkXeM7dnep+DC3MrIO7bzj09CIidVNJeSU/enkeHy3exPXHdOOXZxxGQoJOeNcl+yxAzexu4BdAmpnt/GYzUE74ns39cfcVwKDaCCki8WXJxp3cNyGXjxZvIqNpCr89ux+XjczS2UdpaA51GZaHgLuA9H083wlYs9vjteFtKkBFpEHatLOU65+dxdfrd/K7s/txzdHdgo4kB2GfBai7/wX4i5n9xd3vjmEmEYlTq7YW8+DEXP41fz1NGyVxx6m9ufbobjRpFMlkCpH65VCWYTGzs4DN7j7HzE7Y1257+9h9vJ86yotIvbZ4w06uf2YWO3ZVMP7qbK3xWYcd8LdGd7/bzDoBXXbff89utiJSf23aWcojk5by6qw1JCUa3zuuBzcf350WjVOCjiYSD7qa2V+AfsC37Rfdvft+XnM0cI6ZnRF+TTMze8Hdr9htn7VA590eZ1J91fV/qKO8iNRnHy/ZxK0vfUl6ajL/vPlI+ndsHnQkOQQHLEDN7B7gEuBroCq82VE3W5F6b0dJOX+fspxnP19JZZVzyYjO3HpSL9qpxbnI7p4Gfgs8SPXSZdey96uX3wrPLLobIHwF9I49ik+Ad4AfmtkrVDcfKtD9nyLSkLg7T07L48/vLaZ/x+Y8cVU27Zvrd5C6LpJ5c+cDfdy9LNphRCQ+FJdV8tS0PMZ9uoKi8krOG9yJH4/qTVbrxkFHE4lHae4+yczM3VcBvzOzqVQXpTViZjcDuPtY4D2qmxwto3oZlmtrMbOISFyrqArxm38t5OWZazh9QHseuGgwaSnqNVEfRFKArgCSARWgIvVcWWUVL81YzWOfLCO/qJxT+rXjp6f2pm/7va0QISJhpWaWACw1sx8C64C2kb7Y3ScDk8Pfj91tuwO31GpSEZE6YEdJOd9/YS5frNjKLSf24Ken9FGn23okkgK0BJhnZpPYrQh19x9FLZWIxFRlVYg3v1zHwx8tZd2OXRzZvTXjrurD0KyWQUcTqQtuBxoDPwL+QPU03KuDDCQiUlet2FLE9c/OZt32Xdx/4SC+Oywz6EhSyyIpQN8Jf4lIPePufLBwI/d9mMPyLcUMymzOX787kKN7tia8hq+IHIC7zwp/W4SmyYqIHLTPl+Xz/RfnkphgvHjjSIZ3bRV0JImCSLrgPmtmaUCWu+fEIJOIRJm7M3VpPmMm5PDVugJ6tm3K2CuGMrp/exWeIiIiEnMvz1zNr99eSLeMJjx1zXA6t1Lfifoqki64ZwP3ASlANzMbDPze3c+JcjYRiYI5q7Zz7wdLmJG3jU4t0rjvwkGcP6QTibq3QkRERGKsKuT8+b3FPDktj+N7t+HRy4bQLDU56FgSRZFMwf0dMIL/a5Awz8y6RTGTiETB4g07uf/DHD5avJmMpo34f+f055IRnWmUpI5yIiIiEntFZZX86OUv+XjJZq45qiu/OvMwkhITgo4lURZJAVrp7gV7TMvTItcidcSqrcU8MDGXd+avp2mjJO4c3YdrjupKk0aR/PMXkQMxszbAjUBXdhtX3f26oDKJiMS7tdtLuOHZ2SzdXMQfzu3PlUd2DTqSxEgkv4EuNLPLgEQz60V1l7/PoxtLRA7Vpp2lPDxpKa/NWkNSonHz8T343nHdadE4JehoIvXNv4CpwEdAVcBZRETi3tzV27npudmUVYZ45trhHNurTdCRJIYiKUBvBX5J9RIsLwETgD9GM5SIHLztxeWMnbKcZz5fSVXIuXREFree1JO2zVKDjiZSXzV2958FHUJEpC7417x13Pn6Ato3S+WVm7Lp2TY96EgSY5F0wS2hugD9ZfTjiMjBKi6r5MlpeTzx6QqKyis5f3Anbh/Vm6zW6iInEmXvmtkZ7v5e0EFEROKVu/PQR0t5eNJSRnRrxdgrhtGqiWZlNUSRdMGdCFzo7jvCj1sCr7j76ChnE5EIlFVW8eL01Tz2yTK2Fpdzar92/PTUPvRprzOKIjFyG/ALMysHKsLb3N2bBZhJRCRuFJVVcsdr8/lg0UYuHJbJn84/nJQkNRtqqCKZgpvxTfEJ4O7bzaxt9CKJSCQqq0K8OXcdD32Uy/qCUo7u2Zo7R/dlcOcWQUcTaVDcXWd7RET2YcWWIm56fg55+cX86szDuP6YblpzvIGLpAANmVmWu68GMLMuqAuuSGBCIef9hRu5f2IOK7YUM6hzC8ZcOIije2YEHU2kwTKzc4Djwg8nu/u7QeYREYkHkxZv4vZX5pGclMDz14/gqB76XUUiK0B/CUwzsynhx8cBN0UvkojsjbszJXcL932Yw8J1O+ndrin/uHIYp/ZrpzOJIgEys3uA4cCL4U23mdkx7v7zAGOJiAQmFHIenlR9v+eATs34x5XZdGqRFnQsiRORNCH6wMyGAkcABvzY3fOjnkxEvjV75TbunZDDzLxtZLZM4/4LB3HekE4kJqjwFIkDZwCD3T0EYGbPAl8CKkBFpMEpKKng9le/5JOcLXx3aCZ/On8AqcmJQceSOBLpSvSNgG3h/fuZGe7+afRiiQjA1+t3ct+HOXy8ZDNt0hvxh3P7c/HwLN24LxJ/WlA9TgI0DzCHiEhgvl6/k5tfmMOGgl388bwBXD4yS7O05H9E0gX3r8DFwCIgFN7sgApQkShZmV/MAxNzeWf+epqlJnHXaX245qiuNE6J9JyRiMTQX4AvzewTqmcKHQfcHWwkEZHYeuvLtdz95le0SEvh1e8dydCslkFHkjgVyW+z5wF93L0syllEGrwNBbt4ZNIyXpu9hpTEBH5wQg++d1wPmjdODjqaiOyDu79sZpOpvg/UgJ+5+8ZgU4mIxEZ5ZYg//edrnv1iFSO7teJvlw2lTXqjoGNJHIukAF0BJAMqQEWiZFtxOX+fvIxnv1iFu3PFyCxuOaknbdNTg44mIvtgZn3dfUm4TwLA2vCfHc2so7vPDSqbiEgsbNpZyg9enMucVdu54Zhu/Pz0viQl6jYh2b9ICtASYJ6ZTWK3ItTdfxS1VCINRFFZJU9OzeOJqSsoKa/k/CGZ3D6qF51bNQ46mogc2E+o7gp//16ec+Ck2MYREYmdmXnbuOWluRSXVfLopUM4e1DHoCNJHRFJAfpO+EtEaklpRRUvTF/F45OXs624nNH92/HTU/vQu53WsxepK9z9myXJTnf30t2fMzNNXxCResndefqzlfz5vcV0btWYF64fSZ/2+v1FIhfJMizPmlkK0Du8KcfdKyL9ADNLBGYD69z9rIOLKVI/VFaFeH3OWh6etJQNBaUc0zODO0f3YVDnFkFHE5GD9zkwNIJtIiJ1Wkl5JXe/+RX/mreeU/q14/6LBtEsVX0qpGYi6YJ7AvAssJLq5gqdzezqGizDchuwGGh2cBFF6r5QyHlv4QYe+DCXFfnFDOrcgvsvHMRRPTOCjiYiB8nM2gOdgDQzG0L1GAnV453m0YtIvbIyv5ibX5hDzqZC7hzdh+8f34MErUcuByGSKbj3A6e6ew6AmfUGXgaGHeiFZpYJnAn8iep7ZUQaFHdncu4W7puQw6L1O+ndrinjrhzGKf3aaV0skbpvNHANkAk8sNv2QuAXQQQSEYmGSYs3cfur80hMMJ65dgTH924TdCSpwyIpQJO/KT4B3D3XzCK91v4QcBewz4nhZnYT1U0cyMrKivBtReLf7JXbuPeDHGau3EbnVmk8cNEgzh3ciUSdLRSpF9z9WeBZM/uuu78RdB4RkdoWCjkPTVrKI5OW0r9jM8ZeMUyNEuWQRVKAzjazJ4Hnw48vB+Yc6EVmdhaw2d3nhKfx7pW7jwPGAWRnZ3sEeUTi2qL1Bdw3IYdPcrbQJr0Rfzi3PxcPzyIlSW3JReojd3/DzM4E+gOpu23/fXCpREQOzZbCMn7y2jymLs3ngmGZ/PG8AaQmJwYdS+qBSArQ7wO3AD+i+v6WT4HHI3jd0cA5ZnYG1QNyMzN7wd2vONiwIvEsL7+YBybm8u/562melsxdp/XhmqO60jglkn9mIlJXmdlYqu/5PBEYD1wAzAw0lIjIIfhsWT63vzqPnbsq+PP5h3PpiM66dUhqTSS/GScBD7v7A/BtV9tGB3qRu98N3B1+zQnAHSo+pT7aULCLRyYt5bXZa0lJTOCHJ/bkxuO60zxNXeFEGoij3H2gmS1w9/9nZvcDbwYdSkSkpiqrQjz00VIem7yMHm2a8vz1I+jbXn1EpXZFUoBOAkYBReHHacCHwFHRCiVSF2wrLufxT5bx3PRVuDtXHtGFW07sSZv0A56fEZH65Zs1QEvMrCOwFegWYB4RkRpbu72E21+Zx+xV27koO5PfndNfs7gkKiL5W5Xq7t8Un7h7kZnV6O5jd58MTK5ZNJH4VFRWyfipKxg/NY+S8krOH5LJ7aN66aZ8kYbr32bWAhgDzAUceCLQRCIiNfD+Vxv42RsLCDk8fMlgzh3cKehIUo9FUoAWm9lQd58LYGbDgF3RjSUSf0orqnhh+ioen7ycbcXlnNa/PT89tTe92u2zybOI1HNmlgBMcvcdwBtm9i7VJ24Lgk0mInJgpRVV/P7dr3lpxmoGZTbnkUuH0KV1k6BjST0XSQF6O/BPM1sfftwBuDhqiUTiTGVViNfnrOXhSUvZUFDKsb0yuOPUPgzq3CLoaCISMHcPhe/5PDL8uAwoCzaViMiB5Wws5NaX55K7qYjvHdedn57aRx37JSYOWIC6+ywz6wv0oboL7hJ3r4h6MpGAhULOf77awAMTc8nLL2ZIVgvuv2gQR/XICDqaiMSXD83su8Cb7q7lxEQkrrk7L85YzR/e/Zr01CSeu24Ex/VuE3QsaUAiurM4XHAujHIWkbjg7kzO2cKYCTl8vWEnfdunM/6qbE4+rK1akIvI3vwEaAJUmlkp1Sdr3d3VOlJE4kpBSQU/f3MB7y/cyLG9MnjgosFqnigxp9ZWIruZtXIb936whFkrt5PVqjEPXTyYcwZ1JCFBhaeI7J2760ZwEYl7s1Zu47aXv2RzYRl3n96XG4/trt9vJBAqQEWAResLuG9CDp/kbKFteiP+eN4ALh7emeRE3QshIvtnZpPc/eQDbdvL61KBT6leWzsJeN3df7vHPs2BF4Cs8D73ufvTtZlfROq3qpDz2CfLeOijXDJbNub17x/FYPWxkAAdsAC16jmHlwPd3f33ZpYFtHf3mVFPJxJlK7YU8cDEXN5dsIHmacn8/PS+XH1kV9JSEoOOJiJxLlxANgYyzKwl1VNvAZoBHSN4izLgpPDyZsnANDN7392n77bPLcDX7n62mbUBcszsRXcvr8X/FBGppzYWlHL7q18yfcU2zh3ckT+eN4D01OSgY0kDF8kV0MeBEHAS8HugEHgDGB7FXCJRtX7HLh79eCmvzV5Lo6QEbj2pJzce151m+qEsIpH7HtWd4jsCc/i/AnQn8NiBXhxuWPTNOtvJ4a89mxg5kB4+GdwU2AZUHmpwEan/Jn69iTtfn095ZYj7LhzEd4d2Ui8LiQuRFKAj3X2omX0J4O7bzSwlyrlEomJrURmPT17O89NXgcOVR3ThlhN76gZ8Eakxd38YeNjMbnX3Rw/mPcwskeritSfwmLvP2GOXvwHvAOuBdOBidw/t5X1uAm4CyMrKOpgoIlJPlFZUcc/7S3jm85X079iMRy4dQo82TYOOJfKtSArQivAA6QDhKUD/M/iJxLPC0grGT83jyWl5lJRX8t2hmdw2qheZLRsHHU1E6jh3f9TMjgK6stu46u7PRfDaKmCwmbUA3jKzAe6+e9f50cA8qmch9QAmmtlUd9+5x/uMA8YBZGdnaykYkQbq6/U7+clr81iysZBrj+7Kz0/vS6Mk3VYk8SWSAvQR4C2grZn9CbgA+FVUU4nUktKKKp7/YhWPT17G9pIKTuvfnjtG96ZnWzWtFJHaYWbPU10czgOqwpsdOGAB+g1332Fmk4HT+O9lz64F7glP111mZnlAX0B9GETkW5VVIcZOWc7Dk5bSonEKT12TzUl92wUdS2SvDliAuvuLZjYHOJnq+1vOc/fFUU8mcggqqkK8Pmctj0xayoaCUo7tlcGdo/swMLNF0NFEpP7JBvqFi8SIhWcUVYSLzzRgFPDXPXZbTfX4O9XM2gF9gBW1kFlE6ollm4v46WvzmL+2gLMGduAP5w6gZRPdLSfxK5IuuEcAi9z9sfDjdDMbuZf7VEQCFwo57361gQcn5pKXX8yQrBbcf9EgjuqREXQ0Eam/FgLtgQ01fF0H4NnwbS4JwGvu/q6Z3Qzg7mOBPwDPmNlXVJ8E/pm759dedBGpq6pCzlPT8hjzYQ6NUxJ59NIhnD0okgbcIsGKZAru34Ghuz0u3ss2kUC5O5NztjBmQg5fb9hJn3bpPHFVNqMOa6uObyISbRnA12Y2k+qlVQBw93P29yJ3XwAM2cv2sbt9vx44tfaiikh9sDK/mDv+OZ/Zq7Yz6rB2/Pk7A2ibnhp0LJGIRFKA2u7Titw9ZGaRvE4kJmbmbWPMhCXMWrmdrFaNeejiwZw9qCOJCSo8RSQmfhd0ABFpGEIh57kvVnLPB0tISUzggYsGcf4QLa8idUskheQKM/sR1Vc9AX6A7j+ROLBofQH3Tcjhk5wttE1vxB/PG8DFwzuTnJgQdDQRaUDcfYqZdQF6uftHZtYYUNtJEalVq7eWcOfr85mRt40T+7ThL98ZSPvmuuopdU8kBejNVHfC/RXVXf0mEV5rTCQIK7YU8cDEXN5dsIHmacn8/PS+XH1kV9JS9PueiMSemd1I9bjYiupuuJ2AsVQ3DxIROSTuzoszVvPn9xaTaMa93x3IhdmZuuopdVYkXXA3A5fEIIvIfm0o2MXDHy3ln3PW0igpgR+e2JMbj+tO87TkoKOJSMN2CzACmAHg7kvNrG2wkUSkPli3Yxc/e30B05blc2yvDO757kA6tUgLOpbIIYmkC24qcD3QH/j2Or+7XxfFXCLf2lZczuOfLOO56avA4cojunDLiT1pk94o6GgiIgBl7l7+zdWIcJ+EGi3JIiKyO3fntdlr+MO7i3F3/nT+AC4bkaWrnlIvRDIF93lgCTAa+D1wOaB1QCXqisoqGT91BeOn5lFSXsl3hmZy+6heZLZsHHQ0EZHdTTGzXwBpZnYK1b0S/h1wJhGpozYWlPLzNxcwOWcLR3Zvzb0XDKRzK/3uI/VHJAVoT3e/0MzOdfdnzewlYEK0g0nDVVpRxQvTV/H45OVsKy7ntP7tuWN0b3q2TQ86mojI3vyc6plCXwHfA94DxgeaSETqHHfnzbnr+N2/F1FZ5fz+3P5cMbILCerqL/VMJAVoRfjPHWY2ANgIdI1aImmwKqtCvD5nLQ9PWsqGglKO7ZXBnaP7MDCzRdDRRET2Jw14yt2fADCzxPC2kkBTiUidsXlnKb946ys+WryZEV1bMebCgXRp3SToWCJREUkBOs7MWlLdBfcdoCnw66imkgYlFHLeW7iBBz7MZUV+MYM7t+D+CwdxVM+MoKOJiERiEjAKKAo/TgM+BI4KLJGI1AmhUPW9nn95fwmlFVX8+qx+XHtUV131lHptnwWomd3m7g8Di919O/Ap0D1myaTec3em5G5hzIQcFq3fSe92TRl35TBO6ddON9mLSF2S6u7fFJ+4e1F4LVARkX3K3VTIL978itmrtjOiWyv+8p3D6dGmadCxRKJuf1dArwUeBh4Fhtb0jcPdcz8FGoU/53V3/+3BhJT6Z86qbfz1gxxm5m0js2UaD1w0iHMHdyJRZ/xEpO4pNrOh7j4XwMyGAbsCziQicWpXeRWPfLyUJz5dQXpqEmMuGMgFw7SupzQc+ytAF5vZSqCtmS3YbbsB7u4DD/DeZcBJ4TPBycA0M3vf3acfWmSpy5Zs3Ml9E3L4aPFmMpo24vfn9ueS4VmkJCUEHU1E5GDdDvzTzNaHH3cALg4ujojEq4+XbOI3/1rE2u27uGBYJr844zBaNUkJOpZITO2zAHX3S82sPdUdb8+p6Ru7u/N/98Mkh7+0LloDtWprMQ9MzOWd+etJb5TEnaP7cO3RXWmcEsltyCIi8cvdZ5lZX6AP1Sdpl7h7xQFeJiINyIaCXfy/d77mg0Ub6dW2Ka/edAQju7cOOpZIIA702/8W4Ct3X3Uwbx7uBDgH6Ak85u4z9rLPTcBNAFlZWQfzMRLHNu0s5ZFJS3l11hqSEo2bj+/Bzcf1oHnj5KCjiYjUpuFUd4hPAoaYGe7+XLCRRCRolVUhnvl8JQ9OzKXKnbtO68MNx3TXzC9p0PZbgLp7lZllmFmKu5fX9M3dvQoYbGYtgLfMbIC7L9xjn3HAOIDs7GxdIa0nCkoq+PuU5TzzeR6VVc4lIzrzo5N60bZZatDRRERqlZk9D/QA5gFV4c0OqAAVacDmrt7OL99ayOINOzmxTxt+f+4AOrdSfzKRSOY/rgI+M7N3gOJvNrr7A5F+iLvvMLPJwGnAwgPsLnVYSXklT3+2krFTllNUVsm5gzry41N6ay0rEanPsoF+4VtPRKSBKyip4K8TlvDyzNW0S09l7BVDGd2/vZoMiYRFUoCuD38lAOmRvrGZtQEqwsVnGtVrpP31oFJK3CuvDPHyzNU8+vEy8ovKGHVYW+4Y3Ye+7ZsFHU1EJNoWAu2BDUEHEZHguDtvz1vHn/6zmG3F5Vx7VDd+cmpvmjZSvwuR3R3wX4S7/7+DfO8OwLPh+0ATgNfc/d2DfC+JU1Uh51/z1vHgR7ms2baLEd1a8Y8rhzKsS6ugo4mIxEoG8LWZzaS6AzwA7l7jBn4iUjct31LEr99eyOfLtzKocwueuXYEAzo1DzqWSFw6YAFqZp+wl+617n7S/l7n7guAIQcfTeKZu/PR4s3cNyGHnE2F9OvQjGeuHcDxvdtoiomINDS/CzqAiASjtKKKxycvZ+zk5TRKTuCP5w3g0hFZWtdcZD8imRNwx27fpwLfBSqjE0fqgukrtnLvB0uYu3oH3TKa8OilQzjz8A4k6IetiDRA7j7FzNpR3QkXYKa7bw4yk4hE36e5W/j1vxayamsJ5w7uyK/O7Eeb9EZBxxKJe5FMwZ2zx6bPzGxKlPJIHFu4roAxE3KYkruF9s1S+ct3DueCYZkkJ6qVuIg0XGZ2ETAGmEz1OqCPmtmd7v56oMFEJCry8ov5y3uL+fDrTXTLaMKLN4zk6J4ZQccSqTMimYK7+818CcAwqpstSAORl1/MAxNz+ff89bRonMwvzujLVUd2JTU5MehoIiLx4JfA8G+ueoab8H0EqAAVqUd2lJTzyKRlPPfFSholJXDn6D5cf0w3/T4kUkORTMGdQ/U9oEb11Ns84PpohpL4sGlnKQ9PWsqrs9aQkpjArSf15MbjutMsNTnoaCIi8SRhjym3W6k+YSsi9UB5ZYjnp6/ikUlLKSyt4OLhnfnxKb1pm661zUUORiRTcLvFIojEj4KSCv4+ZTlPf5ZHyJ0rRmbxw5N66b4GEZG9+8DMJgAvhx9fDLwfYB4RqQXuzoRFm7jn/cWs3FrCMT0z+OWZh3FYBy0xJ3Io9lmAmtlwYI27bww/vorqBkSrgN+5+7bYRJRYKSmv5OnPVjJ2ynKKyio5f3AnfnxKbzq3ahx0NBGRuOXud5rZd4BjqJ4tNM7d3wo4logcgq/WFvCH/3zNzLxt9GzblKevHc4J6vQvUiv2dwX0H8AoADM7DrgHuBUYDIwDLoh2OImNiqoQr8xawyOTlrKlsIxRh7XljtF96NteZ/hERPbFzHoC7dz9M3d/E3gzvP04M+vh7suDTSgiNbWhYBdjPsjhzS/X0bpJCn84bwCXDu9MkhouitSa/RWgibtd5byY6jO6bwBvmNm8qCeTqAuFnHe/2sD9H+awamsJw7u25O+XDyW7a6sDv1hERB4CfrGX7SXh586OZRgROXjFZZX8Y8pyxk1dQcjh+yf04Psn9FDfC5Eo2G8BamZJ7l4JnAzcFOHrJM65O58uzefeD5awaP1O+rZP56lrsjmxT1tNLRERiVxXd1+w50Z3n21mXQPIIyI1VBVyXp+zhvs+zGVLYRnnDOrInaP76PYjkSjaXyH5MjDFzPKBXcBU+HbKUUEMskkUfLl6O/d+kMMXK7aS2TKNBy8exLmDOpGQoMJTRKSG9tcCMy1mKUTkoExbms8f//M1SzYWMjSrBf+4chhDs1oGHUuk3ttnAerufzKzSUAH4EN39/BTCVTfCyp1yLLNRdw3IYcPFm2kdZMUfnd2Py4b2YWUJN3TICJykGaZ2Y3u/sTuG83seqqXMBOROLR0UyF/fm8xn+RsIbNlGn+7bAhnHt5Bs8BEYmS/U2ndffpetuVGL47Uto0FpTz0US6vzV5DWnIit4/qxQ3HdqdpI82iFhE5RLcDb5nZ5fxfwZkNpADnBxVKRPZuS2EZD0/K5eWZa2icksgvzujL1Ud1pVFSYtDRRBoUVSH11J5reV51ZFd+eFJPMppqLU8Rkdrg7puAo8zsRGBAePN/3P3jAGOJyB62FpUx7tMVPPfFKsqrQlwxMovbRvWmVZOUoKOJNEgqQOuZ0ooqnv18JY99sozCskrOG9yJn2gtTxGRqHH3T4BPgs4hIv9tW3F5uPBcSWlFFecM6siPTu5F9zZNg44m0qCpAK0nqkLOG3PX8uDEXDYUlHJCnzbcNbov/TpqLU8RERFpOHaUlPPE1BU889lKSiqqOHtgdeHZs60KT5F4oAK0jnN3Pl6ymb9+sITcTUUMymzO/RcN4qgeGUFHExGR/TCzVOBToBHV4/Hr7v7bvex3AtXriiYD+e5+fOxSitQdBSUVPDltBU99tpLi8krOOLwDt5/ci17t0oOOJiK7UQFah81dvZ173l/CzLxtdMtowmOXDeWMw9uri5uISN1QBpzk7kVmlgxMM7P3d28AaGYtgMeB09x9tZm1DSirSNwq2FXBU9PyeGpaHoVllZxxeHtuO7k3fdqr8BSJRypA66AVW4oYMyGH9xduJKNpI/5w3gAuGd6Z5EQtqSIiUleElzcrCj9MDn/5HrtdBrzp7qvDr9kcu4Qi8W1naQVPT1vJk9NWsLO0ktH923Hbyb11+5FInFMBWodsKSzjkUlLeXnmalKSErh9VC9uPLY7TbSkiohInWRmiVQv4dITeMzdZ+yxS28g2cwmA+nAw+7+3F7e5ybgJoCsrKyoZhYJWlFZJc98lscTU/Mo2FXBKf3acfuoXvTv2DzoaCISAVUudUBJeSXjp+bxjynLKasMcemILH50ci/apGtJFRGRuszdq4DB4am2b5nZAHdfuNsuScAw4GQgDfjCzKbvuSa3u48DxgFkZ2fveRVVpF4oKqvk2c9X8sTUFewoqeDkvm25fVRvDs9U4SlSl6gAjWOVVSFen7OWBybmsrmwjNP6t+eu0/qofbiISD3j7jvCVzlPA3YvQNdS3XioGCg2s0+BQUDu/76LSP1UVFbJc1+s5IlPV7C9pIIT+rThx6N6M6hzi6CjichBUAEah9ydyTlb+PN7i1m6uYihWS34+xVDGdalVdDRRESklphZG6AiXHymAaOAv+6x27+Av5lZEpACjAQejG1SkWDsKCnnxRmrGT/1/wrP207uxZCslkFHE5FDoAI0zixaX8Cf31vMZ8u20rV1Y/5++VBOG6DOtiIi9VAH4NnwfaAJwGvu/q6Z3Qzg7mPdfbGZfQAsAELA+D2m6IrUO6u2FvPUtDxem72WXRVVKjxF6hkVoHFiY0EpYybk8OaXa2mRlszvzu7HZSO7kJKkzrYiIvWRuy8Ahuxl+9g9Ho8BxsQql0hQ5qzaxhOf5jHh640kJRjnDu7E9cd047AO6morUp9ErQA1s87Ac0B7qs/ajnP3h6P1eXVVcVkl/5iynHFTVxAKwU3HducHJ/akeVpy0NFEREREoqqyKsSERZsYP20FX67eQfO0ZH5wQg+uPrIrbZulBh1PRKIgmldAK4GfuvtcM0sH5pjZRHf/OoqfWWdUhZw35qxlzIc5bCks46yBHfjZaX3p3Kpx0NFEREREoqqorJLXZq3hqc/yWLt9F11aN+b35/bngmGZNE7RBD2R+ixq/8LdfQOwIfx9oZktBjoBDb4Anb5iK39492sWrd/J4M4tGHvFMIZ10X0NIiIiUr9tKNjFM5+v5KUZqyksrSS7S0t+deZhnNKvPYkJ6nch0hDE5BSTmXWl+j6XPRfYblBWby3hz+8t5oNFG+nYPJWHLxnMOYM6qsGQiIiI1GsL1xUwfuoK3l2wgZA7pw/owA3HdlNjIZEGKOoFqJk1Bd4Abnf3nXt5/ibgJoCsrKxoxwlEcVklj32yjPFT80hKNH56Sm9uPK47qcmJQUcTERERiYpQyPkkZzNPTsvj8+VbaZKSyFVHduXao7vqliORBiyqBaiZJVNdfL7o7m/ubR93HweMA8jOzvZo5om1UMh5e9467nl/CZsLy/jOkE7cdVpf2jfXTfUiIiJSP5WUV/LG3HU8PS2PFfnFdGieyt2n9+WSEVlqsigiUe2Ca8CTwGJ3fyBanxOvvlpbwG/eWciXq3cwKLM5Y68cxlBNMxEREZF6asnGnbwycw1vfbmOgl0VDMpszsOXDOaMwzuQnKhl5USkWjSvgB4NXAl8ZWbzwtt+4e7vRfEzA7etuJwxE3J4ZdZqWjdJ4d4LBnLB0EwSdGO9iIiI1DNFZZX8e/56Xpm1hvlrdpCSmMDoAe25+sguDOvSUn0uROR/RLML7jSgwfzUqQo5r8xazb0f5FBUVsl1R3fjtlG9aJaqqSYiIiJSf7g789bs4JWZa/j3gvWUlFfRq21Tfn1WP74zpBMtm6QEHVFE4pgWWqoFX60t4Fdvf8X8tQUc0b0Vvz93AL3bpQcdS0RERKTW7Cgp5+0v1/HKrDUs2VhIWnIiZw/qwMXDsxia1UJXO0UkIipAD8HO0grum5DD89NX0bpJIy2rIiIiIvWKuzN9xTZenbWa9xZupLwyxMDM5vzp/AGcM6gj6ZrpJSI1pAL0ILg77y/cyO/eWUR+URlXHdGFn47uo+m2IiIiUi9sLizljTnreHXWalZuLSE9NYlLhnfm4uGd6d+xedDxRKQOUwFaQ2u3l/Cbfy3i4yWb6d+xGeOvzmZgZougY4mIiIgcksqqEJNztvDa7DV8vGQzlSFnRNdW/OjkXpw+oANpKVq/XEQOnQrQCIVCzgszVnHP+0twh1+deRjXHNWVJLUVFxERkTrK3VmysZC3563jzbnr2FJYRkbTFK47phsXD+9MjzZNg44oIvWMCtAIrMwv5q43FjAzbxvH9srgL985nMyWjYOOJSIiIlJj7k7upiL+s2A97361gRVbiklMME7s05aLsjM5sW9brdspIlGjAnQ/QiHnqc/yuO/DHJITE7j3goFcOCxTTYZERESkzlm2uZB3F2zg3QUbWLa5iASDI7q35rqju3HagPZkNG0UdEQRaQBUgO7Duh27uOO1+XyxYisn923Ln84/nPbNU4OOJSIiIhKx5VuK+M+CDfxnwQZyNhViBiO6tuLqc/tz2oAOtElX0SkisaUCdA/uztvz1vGbtxcRcufe7w7kwmxd9RQREZG6YWV+Mf/5qvpK5+INOzGD7C4t+X/n9Of0Ae1p20wn1EUkOCpAd1NUVskv3vyKd+avJ7tLSx64aDBZrXWvp4iIiMS31VtLwkXnehat3wnA0KwW/OasfpxxeAfN4hKRuKECNOzr9Tu55aW5rNpazB2n9ub7J/QkMUFXPUVERCQ+Ld1UyPsLN/LBwo18vaG66ByS1YJfnXkYZxzegY4t0gJOKCLyvxp8AeruvDJrDb97ZxHN05J5+cYjGNm9ddCxRERERP5LSXkl01dsZUrOFqbkbmHl1hLMYFhWS3515mGcNqC9uvSLSNxr0AVoeWWIX7+9kFdnr+HYXhk8ePFgdYATERGRuPDNcilTcjczJXcLs/K2U14VIi05kSN7tOb6Y7oxur/u6RSRuqXBFqA7Ssr5/gtz+WLFVm49qSe3j+qtKbciIiISqIKSCj5bns+UnC18unQLGwpKAejdrinXHN2V43u3IbtrSxolJQacVETk4DTIAjQvv5jrn5nF2u27eOjiwZw3pFPQkURERKQBCoWcr9YVMCW3elrtl6u3E3JIT03imJ4Z3D6qDcf1bkOH5rqfU0TqhwZXgH65ejvXPjOLBDNevHEkw7u2CjqSiIiINCDrduzis6X5TFtW/bWtuBwzGNipObec2JPje7dhcOcWJCUmBB1VRKTWNagCdOG6Aq56aiatmqTw3HUj6NK6SdCRREREpB5zd1ZuLWHuqu3MXrWd6Su2kpdfDEBG00ac0KcNx/duwzE9M2itPhQi0gA0mAJ06aZCrnpqJumNknjxhpHqEiciIiK1bld5FfPX7mDu6u3MXbWduat3sK24HID0RkkM79aKK4/owjG9MujVtilm6j8hIg1LgyhAV+YXc/n4GSQlGC/deISKTxERETlk7s7a7bv+q9j8esNOqkIOQPc2TTipb1uGdWnJsC4t6dmmKQlqeCgiDVy9L0Dzi8q4fPwMKkPOqzcdQdcMTbsVERGRmiutqGLhugLmrt7OnHDBuaWwDIDGKYkMymzBzcd3Z1iXlgzp3JKWTVICTiwiEn/qfQE67tMVbCjYxb9uOYZe7dKDjiMiIiJ1QMGuCpZuKmTJxkIWritg/toCcjcVfnt1M6tVY47u0ZqhXVoyNKslfdunq2mQiEgE6nUBur24nBemr+LsQR05PLN50HFEREQkzhSWVrB0cxFLNxWSu6mI3E2FLN1UxMadpd/u0zwtmYGZzTm5bw8Oz2zO0KyWtElXwyARkYNRrwvQZz5fSUl5FT84oWfQUURERCRAxWWVLNtcXWDmhovNpZsKWV/wf4VmanICPds25agerenVLp3e7ZrSu106mS3T1CxIRKSW1NsCtKiskmc+X8kp/drRp72m3oqIiNR3u8qrWF+wi3Xbd7FyazErthSTl1/M8i1FrN2+69v9UpIS6NGmKcO7taJ3u3R6tW1Kn/bpZLZsTKKaBImIRFW9LUBfmL6Kgl0V/PBEXf0UERGp60orqthYUMr6gl1s2FHKhoJdbCgoZUNBKet37GLjzlJ2lFT812uapCTSNaMJQ7JacnF252+vama1aqz7NUVEAhK1AtTMngLOAja7+4Bofc7elFZUMX5qHsf2ymBQ5xax/GgREZGImFkq8CnQiOrx+HV3/+0+9h0OTAcudvfXY5cyusorQxSWVlBUVsn2kgq2FpWRX1RGflE5G8PF5TeF5jdrae6uZeNkOjRPo1OLNLK7tqRD8zQ6NE+lY4s0umc0oU16I02dFRGJM9G8AvoM8DfguSh+xl69NnsN+UVl3HLikFh/tIiISKTKgJPcvcjMkoFpZva+u0/ffSczSwT+CkwIIuTeVFaFKCqrpLC0+qv6+4pvt337OPx84Z7Ph7eVV4b2+RnN05Lp0DyVDs1TGdS5BR2apdKhRRodm1f/2b5ZKmkpiTH8rxYRkdoQtQLU3T81s67Rev99qagK8Y8pKxjWpSUju7WK9ceLiIhExN0dKAo/TA5/+V52vRV4Axgeo2gAPPHpCuas2k5hWcW3BeM3xeOuiqoDvj4xwUhPTaJpo+qv9NQk2qan0j0jiaapSaSHt6WnJtO0URLN05LJSG9ERtMUMpo2IjVZxaWISH0U+D2gZnYTcBNAVlbWIb/fuu27SElK4Icn9tS0GxERiWvhq5tzgJ7AY+4+Y4/nOwHnAyexnwK0tsdSgFXbqhv4NE1NokXjFDq3avxtIZmemkyT8PfNUpNo2ii5uqj85qtRMqnJCRqHRUTkfwRegLr7OGAcQHZ29t7O/NZI14wmfPST41ETOxERiXfuXgUMNrMWwFtmNsDdF+62y0PAz9y9an/FXG2PpQB/PO/w2ngbERGR/xJ4ARoNaqEuIiJ1ibvvMLPJwGnA7gVoNvBKuPjMAM4ws0p3fzvmIUVERGpBvSxARURE4p2ZtQEqwsVnGjCK6mZD33L3brvt/wzwropPERGpy6K2CJaZvQx8AfQxs7Vmdn20PktERKQO6gB8YmYLgFnARHd/18xuNrObA84mIiISFdHsgntptN5bRESkrnP3BcD/rBfm7mP3sf810c4kIiISbVG7AioiIiIiIiKyOxWgIiIiIiIiEhMqQEVERERERCQmVICKiIiIiIhITKgAFRERERERkZgwdw86w7fMbAuwqpbeLgPIr6X3agh0vGpGx6vmdMxqRserZmrzeHVx9za19F4xp7E0UDpeNadjVjM6XjWj41VzUR9P46oArU1mNtvds4POUVfoeNWMjlfN6ZjVjI5Xzeh4RYeOa83oeNWcjlnN6HjVjI5XzcXimGkKroiIiIiIiMSEClARERERERGJifpcgI4LOkAdo+NVMzpeNadjVjM6XjWj4xUdOq41o+NVczpmNaPjVTM6XjUX9WNWb+8BFRERERERkfhSn6+AioiIiIiISBxRASoiIiIiIiIxUacLUDM7zcxyzGyZmf18L8+bmT0Sfn6BmQ0NImc8ieCYXR4+VgvM7HMzGxREznhxoOO1237DzazKzC6IZb54E8nxMrMTzGyemS0ysymxzhhPIvj32NzM/m1m88PH69ogcsYLM3vKzDab2cJ9PK+f+QdJ42nNaCytGY2lNafxtGY0ntZM4OOpu9fJLyARWA50B1KA+UC/PfY5A3gfMOAIYEbQuevAMTsKaBn+/vSGfMwiOV677fcx8B5wQdC54/l4AS2Ar4Gs8OO2QeeO8+P1C+Cv4e/bANuAlKCzB3jMjgOGAgv38bx+5h/ccdV4WvvHS2NpDY7Xbvs1+LE00mOm8bTGx0vj6X8fj0DH07p8BXQEsMzdV7h7OfAKcO4e+5wLPOfVpgMtzKxDrIPGkQMeM3f/3N23hx9OBzJjnDGeRPJ3DOBW4A1gcyzDxaFIjtdlwJvuvhrA3RvyMYvkeDmQbmYGNKV6wKyMbcz44e6fUn0M9kU/8w+OxtOa0VhaMxpLa07jac1oPK2hoMfTulyAdgLW7PZ4bXhbTfdpSGp6PK6n+uxHQ3XA42VmnYDzgbExzBWvIvn71RtoaWaTzWyOmV0Vs3TxJ5Lj9TfgMGA98BVwm7uHYhOvTtLP/IOj8bRmNJbWjMbSmtN4WjMaT2tfVH/mJ9XWGwXA9rJtzzVlItmnIYn4eJjZiVQPmsdENVF8i+R4PQT8zN2rqk+qNWiRHK8kYBhwMpAGfGFm0909N9rh4lAkx2s0MA84CegBTDSzqe6+M8rZ6ir9zD84Gk9rRmNpzWgsrTmNpzWj8bT2RfVnfl0uQNcCnXd7nEn1WY2a7tOQRHQ8zGwgMB443d23xihbPIrkeGUDr4QHzAzgDDOrdPe3Y5IwvkT6bzLf3YuBYjP7FBgENMQBM5LjdS1wj1ffkLHMzPKAvsDM2ESsc/Qz/+BoPK0ZjaU1o7G05jSe1ozG09oX1Z/5dXkK7iygl5l1M7MU4BLgnT32eQe4KtzJ6QigwN03xDpoHDngMTOzLOBN4MoGehZtdwc8Xu7ezd27untX4HXgBw14wIzk3+S/gGPNLMnMGgMjgcUxzhkvIjleq6k+u42ZtQP6ACtimrJu0c/8g6PxtGY0ltaMxtKa03haMxpPa19Uf+bX2Sug7l5pZj8EJlDd/eopd19kZjeHnx9LdSe1M4BlQAnVZz8arAiP2W+A1sDj4TORle6eHVTmIEV4vCQskuPl7ovN7ANgARACxrv7XluA13cR/v36A/CMmX1F9XSYn7l7fmChA2ZmLwMnABlmthb4LZAM+pl/KDSe1ozG0prRWFpzGk9rRuNpzQU9nlr1lWgRERERERGR6KrLU3BFRERERESkDlEBKiIiIiIiIjGhAlRERERERERiQgWoiIiIiIiIxIQKUBEREREREYkJFaAi9YSZtTCzHwSdQ0RERERkX1SAitQfLQAVoCIiIrXEqun3ZZFapH9QIvXHPUAPM5tnZmOCDiMiIlIXmVlXM1tsZo8Dc4HOQWcSqU/M3YPOICK1wMy6Au+6+4Cgs4iIiNRV4fF0BXCUu08POI5IvaMroCIiIiIi/22Vik+R6FABKiIiIiLy34qDDiBSX6kAFak/CoH0oEOIiIiIiOyLClCResLdtwKfmdlCNSESERERkXikJkQiIiIiIiISE7oCKiIiIiIiIjGhAlRERERERERiQgWoiIiIiIiIxIQKUBEREREREYkJFaAiIiIiIiISEypARUREREREJCZUgIqIiIiIiEhM/H8my20HSaJRvwAAAABJRU5ErkJggg==", - "text/plain": [ - "
" - ] - }, - "metadata": { - "needs_background": "light" - }, - "output_type": "display_data" - } - ], - "source": [ - "# solve\n", - "solver = pybamm.ScipySolver()\n", - "t = np.linspace(0, 1, 100)\n", - "solution = solver.solve(model, t)\n", - "\n", - "# post-process, so that the solution can be called at any time t or space r\n", - "# (using interpolation)\n", - "c = solution[\"Concentration\"]\n", - "\n", - "# plot\n", - "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4))\n", - "\n", - "ax1.plot(solution.t, c(solution.t, r=1))\n", - "ax1.set_xlabel(\"t\")\n", - "ax1.set_ylabel(\"Surface concentration\")\n", - "r = np.linspace(0, 1, 100)\n", - "ax2.plot(r, c(t=0.5, r=r))\n", - "ax2.set_xlabel(\"r\")\n", - "ax2.set_ylabel(\"Concentration at t=0.5\")\n", - "plt.tight_layout()\n", - "plt.show()" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "In the [next notebook](./3-negative-particle-problem.ipynb) we build on the example here to to solve the problem of diffusion in the negative electrode particle within the single particle model. In doing so we will also cover how to include parameters in a model. " - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## References\n", - "\n", - "The relevant papers for this notebook are:" - ] - }, - { - "cell_type": "code", - "execution_count": 10, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "[1] Joel A. E. Andersson, Joris Gillis, Greg Horn, James B. Rawlings, and Moritz Diehl. CasADi – A software framework for nonlinear optimization and optimal control. Mathematical Programming Computation, 11(1):1–36, 2019. doi:10.1007/s12532-018-0139-4.\n", - "[2] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, and others. Array programming with NumPy. Nature, 585(7825):357–362, 2020. doi:10.1038/s41586-020-2649-2.\n", - "[3] Valentin Sulzer, Scott G. Marquis, Robert Timms, Martin Robinson, and S. Jon Chapman. Python Battery Mathematical Modelling (PyBaMM). Journal of Open Research Software, 9(1):14, 2021. doi:10.5334/jors.309.\n", - "[4] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, and others. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature Methods, 17(3):261–272, 2020. doi:10.1038/s41592-019-0686-2.\n", - "\n" - ] - } - ], - "source": [ - "pybamm.print_citations()" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.8.13" - }, - "toc": { - "base_numbering": 1, - "nav_menu": {}, - "number_sections": true, - "sideBar": true, - "skip_h1_title": false, - "title_cell": "Table of Contents", - "title_sidebar": "Contents", - "toc_cell": false, - "toc_position": {}, - "toc_section_display": true, - "toc_window_display": true - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/docs/source/examples/notebooks/creating_models/old_notebooks/3-negative-particle-problem.ipynb b/docs/source/examples/notebooks/creating_models/old_notebooks/3-negative-particle-problem.ipynb deleted file mode 100644 index c14c1279e6..0000000000 --- a/docs/source/examples/notebooks/creating_models/old_notebooks/3-negative-particle-problem.ipynb +++ /dev/null @@ -1,400 +0,0 @@ -{ - "cells": [ - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# A step towards the Single Particle Model" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "In the [previous notebook](./2-a-pde-model.ipynb) we saw how to solve a PDE model in pybamm. Now it is time to solve a real-life battery problem! We consider the problem of spherical diffusion in the negative electrode particle within the single particle model. That is,\n", - "$$\n", - " \\frac{\\partial c}{\\partial t} = \\nabla \\cdot (D \\nabla c),\n", - "$$\n", - "with the following boundary and initial conditions:\n", - "$$\n", - " \\left.\\frac{\\partial c}{\\partial r}\\right\\vert_{r=0} = 0, \\quad \\left.\\frac{\\partial c}{\\partial r}\\right\\vert_{r=R} = -\\frac{j}{FD}, \\quad \\left.c\\right\\vert_{t=0} = c_0,\n", - "$$\n", - "where $c$ is the concentration, $r$ the radial coordinate, $t$ time, $R$ the particle radius, $D$ the diffusion coefficient, $j$ the interfacial current density, $F$ Faraday's constant, and $c_0$ the initial concentration. \n", - "\n", - "In this example we use the following parameters:\n", - "\n", - "| Symbol | Units | Value |\n", - "|:-------|:-------------------|:-----------------------------------------------|\n", - "| $R$ | m | $10 \\times 10^{-6}$ |\n", - "| $D$ | m${^2}$ s$^{-1}$ | $3.9 \\times 10^{-14}$ |\n", - "| $j$ | A m$^{-2}$ | $1.4$ |\n", - "| $F$ | C mol$^{-1}$ | $96485$ |\n", - "| $c_0$ | mol m$^{-3}$ | $2.5 \\times 10^{4}$ |\n", - "\n" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Setting up the model\n", - "As before, we begin by importing the PyBaMM library into this notebook, along with any other packages we require, and start with an empty `pybamm.BaseModel`\n" - ] - }, - { - "cell_type": "code", - "execution_count": 1, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Note: you may need to restart the kernel to use updated packages.\n" - ] - } - ], - "source": [ - "%pip install \"pybamm[plot,cite]\" -q # install PyBaMM if it is not installed\n", - "import pybamm\n", - "import numpy as np\n", - "import matplotlib.pyplot as plt\n", - "\n", - "model = pybamm.BaseModel()" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We then define all of the model variables and parameters. Parameters are created using the `pybamm.Parameter` class and are given informative names (with units). Later, we will provide parameter values and the `Parameter` objects will be turned into numerical values. For more information please see the [parameter values notebook](../parameterization/parameter-values.ipynb)." - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": {}, - "outputs": [], - "source": [ - "R = pybamm.Parameter(\"Particle radius [m]\")\n", - "D = pybamm.Parameter(\"Diffusion coefficient [m2.s-1]\")\n", - "j = pybamm.Parameter(\"Interfacial current density [A.m-2]\")\n", - "F = pybamm.Parameter(\"Faraday constant [C.mol-1]\")\n", - "c0 = pybamm.Parameter(\"Initial concentration [mol.m-3]\")\n", - "\n", - "c = pybamm.Variable(\"Concentration [mol.m-3]\", domain=\"negative particle\")" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now we define our model equations, boundary and initial conditions, as in the previous example. " - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "metadata": {}, - "outputs": [], - "source": [ - "# governing equations\n", - "N = -D * pybamm.grad(c) # flux\n", - "dcdt = -pybamm.div(N)\n", - "model.rhs = {c: dcdt}\n", - "\n", - "# boundary conditions\n", - "lbc = pybamm.Scalar(0)\n", - "rbc = -j / F / D\n", - "model.boundary_conditions = {c: {\"left\": (lbc, \"Neumann\"), \"right\": (rbc, \"Neumann\")}}\n", - "\n", - "# initial conditions\n", - "model.initial_conditions = {c: c0}" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Finally, we add any variables of interest to the dictionary `model.variables`" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": {}, - "outputs": [], - "source": [ - "model.variables = {\n", - " \"Concentration [mol.m-3]\": c,\n", - " \"Surface concentration [mol.m-3]\": pybamm.surf(c),\n", - " \"Flux [mol.m-2.s-1]\": N,\n", - "}" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Using the model" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "In order to discretise and solve the model we need to provide values for all of the parameters. This is done via the `pybamm.ParameterValues` class, which accepts a dictionary of parameter names and values" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, - "outputs": [], - "source": [ - "param = pybamm.ParameterValues(\n", - " {\n", - " \"Particle radius [m]\": 10e-6,\n", - " \"Diffusion coefficient [m2.s-1]\": 3.9e-14,\n", - " \"Interfacial current density [A.m-2]\": 1.4,\n", - " \"Faraday constant [C.mol-1]\": 96485,\n", - " \"Initial concentration [mol.m-3]\": 2.5e4,\n", - " }\n", - ")" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Here all of the parameters are simply scalars, but they can also be functions or read in from data (see [parameter values notebook](../parameterization/parameter-values.ipynb))." - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "As in the previous example, we define the particle geometry. Note that in this example the definition of the geometry contains a parameter, the particle radius $R$" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": {}, - "outputs": [], - "source": [ - "r = pybamm.SpatialVariable(\n", - " \"r\", domain=[\"negative particle\"], coord_sys=\"spherical polar\"\n", - ")\n", - "geometry = {\"negative particle\": {r: {\"min\": pybamm.Scalar(0), \"max\": R}}}" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Both the model and geometry can now be processed by the parameter class. This replaces the parameters with the values" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": {}, - "outputs": [], - "source": [ - "param.process_model(model)\n", - "param.process_geometry(geometry)" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We can now set up our mesh, choose a spatial method, and discretise our model" - ] - }, - { - "cell_type": "code", - "execution_count": 8, - "metadata": {}, - "outputs": [], - "source": [ - "submesh_types = {\"negative particle\": pybamm.Uniform1DSubMesh}\n", - "var_pts = {r: 20}\n", - "mesh = pybamm.Mesh(geometry, submesh_types, var_pts)\n", - "\n", - "spatial_methods = {\"negative particle\": pybamm.FiniteVolume()}\n", - "disc = pybamm.Discretisation(mesh, spatial_methods)\n", - "disc.process_model(model);" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The model is now discretised and ready to be solved." - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Solving the model" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "As is the previous example, we choose a solver and times at which we want the solution returned." - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "2021-11-19 15:29:22,931 - [WARNING] processed_variable.get_spatial_scale(520): No length scale set for negative particle. Using default of 1 [m].\n" - ] - }, - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6AAAAEYCAYAAABCw5uAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABw/UlEQVR4nO3dd3hVVdbH8e9KIaH3XqT3KkUQ6RbEXrGMXbGgYp1BneI7o2MvY28oWLEPdkSkiDSDIL0X6b2EnrLeP+6JRiYkN5Dk5ia/z/Oc5567T7nr5EJ21tn77G3ujoiIiIiIiEh+i4l0ACIiIiIiIlI8KAEVERERERGRAqEEVERERERERAqEElAREREREREpEEpARUREREREpEDERTqAglalShWvX79+pMMQEZEoMmPGjC3uXjXScUSK6k4REcmtw9WdxS4BrV+/PklJSZEOQ0REooiZrYp0DJGkulNERHLrcHWnuuCKiIiIiIhIgVACKiIiIiIiIgVCCaiIiIiIiIgUiGL3DKiIiIiIiBQ/KSkprFmzhv3790c6lCIlMTGROnXqEB8fH9b+SkBFRERERKTIW7NmDWXLlqV+/fqYWaTDKRLcna1bt7JmzRoaNGgQ1jHqgisiIiIiIkXe/v37qVy5spLPPGRmVK5cOVetyvmWgJpZXTMbZ2YLzGyemQ0Jyu83s7VmNitYBmQ65h4zW2pmi8zslEzlHc1sTrDtGQv+1ZhZgpm9H5RPM7P6+XU9IiIiIiIS3ZR85r3c/kzzswU0FbjT3VsAXYHBZtYy2PaUu7cPlq8Agm0XAa2A/sALZhYb7P8iMAhoEiz9g/JrgO3u3hh4CngkH69HREQkX2Vz8/aC4H26mXXKtP+lmW7ozgq2tw+2jQ9u6GZsqxaU6+atiIhETL49A+ru64H1wXqymS0AamdzyFnASHc/AKwws6VAFzNbCZRz9ykAZvYmcDbwdXDM/cHxHwHPmZm5u+f9Ff1u1Ky1TFy8hScubJefHyMiIsVPxs3bn82sLDDDzMYAc4FzgZcz7+zu7wDvAJhZG2CUu8/KtMul7p50yGf8dvPWzC4idPN2YL5cTSZ//ugXAMomxlM2Me6313LBernfykPvS8TpKSERkaKoQAYhCu6udgCmAd2Bm83sciCJUEW7nVByOjXTYWuCspRg/dBygtfVAO6eamY7gcrAlkM+fxChFlTq1at31Nezbsd+Pv55DUNPbU7VsglHfT4RERE4/M1bdx8DOXZzuhh4L4yPicjN23nrdrF190GS96ew52BajvsnxMUEiWncHxLWjPWqZROoVaEktSskUqtCSaqVTSQ2Rl3rRKTw2rFjB++++y433XRTWPs/99xzPP300yxbtozNmzdTpUoVIDTwz5AhQ/jqq68oVaoUw4cP59hjjwXgm2++YciQIaSlpXHttdcydOhQALZt28bAgQNZuXIl9evX54MPPqBixYoMHz6cu+++m7POOovXXnstyzj27dtHt27dmD9/PuvWrfstjiOV7wmomZUBPgZuc/ddZvYi8C/Ag9cngKuBrGoNz6acHLb9XuD+CvAKQKdOnY66gj2+UWUAJi/bwlnts2vUFREROTKH3LwNx0BCyWVmb5hZGqF6+IEgyYzIzdsvb+3x23pqWjq7D6SSvD+VXftTSN6fGiwpf3jd9YftKWzYtZ/k/Sns2pfKvpQ/JrFxMUaN8olBUlqSWhUSqV2hVPBakloVSlI6QYP/i0jk7NixgxdeeCHsBLR79+6cfvrp9O7d+w/lX3/9NUuWLGHJkiVMmzaNG2+8kWnTppGWlsbgwYMZM2YMderUoXPnzpx55pm0bNmShx9+mH79+jF06FAefvhhHn74YR55JPT04sCBA3nuuecOG0fJkiWZNWsW9evXP9JL/4N8/U1sZvGEKr133P0TAHffmGn7q8AXwds1QN1Mh9cB1gXldbIoz3zMGjOLA8oD2/L+Sv6ode3yVCgVzzdzNygBFRGRPHfozdsw9j8O2OvuczMVX+rua4OuvB8DlwFvEqGbt5nFxcZQoVQJKpQqccTnSN6fwvqd+1m7Yx/rduxj7fbQ67od+5m+Yhsbdu0nLf2PYVcoFU+t8iX/0HJaq0JJGlQpTeNqZUiMjz3Mp4lIUfN/n89j/rocf73mSsta5fjHGa0Ou33o0KEsW7aM9u3bc9JJJ/HYY49le74OHTpkWT5q1Cguv/xyzIyuXbuyY8cO1q9fz8qVK2ncuDENGzYE4KKLLmLUqFG0bNmSUaNGMX78eACuuOIKevfu/VsCmtm8efO46qqrOHjwIOnp6Xz88cc0adIkzJ9AePItAQ1Gqh0GLHD3JzOV1wy6GAGcQ+i5FoDPgHfN7EmgFqHBhqa7e5qZJZtZV0J3gS8Hns10zBXAFOB84Pv87kIEEBtjXNS5Hq9MXMbijck0rV42vz9SRESKiaxu3obhIg7pfuvua4PXZDN7F+hCKAGNyM3bvBbqkht/2Do4Ld3ZuGt/KDkNEtOM9TXb9zJtxVaS96f+tn+MQf0qpWleoyzNa5SjWY2yNK9RlroVSxGjrr0ikgcefvhh5s6dy6xZs0hOTqZ9+/ZZ7vfuu+/SsmXLLLcBrF27lrp1f2+3q1OnDmvXrs2yfNq0UCeajRs3UrNmTQBq1qzJpk2bsjz3Sy+9xJAhQ7j00ks5ePAgaWk5PzKRW/nZAtqd0N3WOWY2Kyi7F7g4GKHPgZXA9QDuPs/MPgDmExqEYbC7Z1zxjcBwoCShwYe+DsqHAW8FAxZtI1QBF4jrejRg5E+/cu8nc/jg+m6qnERE5Kgd7uZtDsfEABcAPTOVxQEV3H1LkNCeDnwXbI7IzduCFhtjv7VwdjrMPrv2p7Buxz6Wb97Dwg3JLFy/i3nrdvH13A1k/ERKlYilafVQMtosSE6b1yhLxdJH3norIpGXXUtlQShbtiyzZs06omOz+pVtZoctz41u3brx4IMPsmbNGs4999w8b/2E/B0FdxJZd/P5KptjHgQezKI8CWidRfl+QpVugatcJoH7BrTg7o9m89LEZdzUu3EkwhARkaLlcDdvEwj1/qkKfGlms9w9Y77snsAad1+e6TwJwOgg+YwllHy+GmyL2M3bwqZcYjzlasTTvEY5BrSp+Vv5ngOpLN6YzKINySzcEHodPW8DI39a/ds+1com0KxGWVrULEez6qHkVN14RSRcycnJ9OjRI8ttObWA1qlTh9Wrf/99tGbNGmrVqsXBgwezLAeoXr0669evp2bNmqxfv55q1aplee5LLrmE4447ji+//JJTTjmF1157jb59+x7JJR6WnsY/Cud3rMOExZt5bPQiGlUtwymtakQ6JBERiWLZ3LwF+PQwx4wnNN925rI9QMfD7B+xm7fRonRCHB3qVaRDvYq/lbk7m5MP/JaQLtiwi0Ubkhk+eSUHU9OBUKtrk2pl6FS/Ip3rV6Jz/UrUqlAyUpchIoVM2bJlSU5O/m39SFtAzzzzTJ577jkuuugipk2bRvny5alZsyZVq1ZlyZIlrFixgtq1azNy5Ejefffd344ZMWIEQ4cOZcSIEZx11qFj1oUsX76chg0bcuutt7J8+XJmz56tBLQwMTMev6Adq7fv47aRs3j72i50PKZSpMMSERGRPGZmVCuXSLVyifRsWvW38tS0dFZu3fNbYjpr9Q4+/Xktb0/9FYDaFUrSuX5FOtWvRJcGlWhctYwe2xEppipXrkz37t1p3bo1p556ao6DED3zzDM8+uijbNiwgbZt2zJgwABee+01BgwYwFdffUXjxo0pVaoUb7zxBgBxcXE899xznHLKKaSlpXH11VfTqlWoq/HQoUO58MILGTZsGPXq1ePDDz/M8jPff/993n77beLj46lRowZ///vf8/aHAFgRfOwjW506dfKkpEPn5D46m5L3M/DlqWzatZ8RV3ehU30loSIiRYmZzXD3wz1KWOTlR91ZlKWmpbNwQzI/rdzGTyu3MX3FdrbsPgCERuLtdEwoIe1cvyJtalegRFxMhCMWKR4WLFhAixYtIh1GoTJ8+HCSkpKynYYlQ/369UlKSspyHtCsfraHqzvVApoHqpVNZOSgrlz8ylQuf306z19yLH2aZ92vWkRERIq2uNgYWtcuT+va5bmqewPcnVVb9/LTym0krdzOTyu38d2C0AiUCXExtKtbgS71K9GpfkWOPaYi5RLjI3wFIlJclCxZkq+//pprr72W1157Lct99u3bR7du3UhJSSEm5uhvmKkFNA9t2rWfq4b/xMINyfzrrNZcctzRT9wtIiKRpxZQtYDmtS27D5C0chs/BQnpvHW7SEt3Ygya1yhH5/oV6dqwMic0qUJZJaQieWLBggU0b9481yPDSvbcnYULF6oFNBKqlUvkg+u7Mfjdn7n30zks27yboac2Jz5WXWtERETkd1XKJNC/dU36tw6NvrvnQCqzVu9g+optJK3axgdJaxgxZRVxMUaXBpXo27wafZtXo2HVMhGOXCR6JSYmsnXrVipXrqwkNI+4O1u3biUxMTHsY9QCmg9S09J54MsFDJ+8ki71K/HcJR2oVi78L0VERAoXtYCqBbSgpaSlM/PXHYxduJFxCzexeONuAOpXLkXf5tXp27waXRpU0vOjIrmQkpLCmjVr2L9/f6RDKVISExOpU6cO8fF/7K1xuLpTCWg+GjVrLUM/nkOZxDiev+RYujTQ4EQiItFICagS0EhbvW0v4xZt4vuFm5i8bCsHU9MpkxDHCY2r0LdFNfo0q0bVsgmRDlNE5DdKQAMFXYku2pDMjW/PYNW2vdzStzE392lMnLrkiohEFSWgSkALk70HU5m8dCtjF25i3MJNbNgVas1pV6c8fZpXo1/z6rSqVU7TvYhIRCkBDUSiEk3en8I/Rs3jk5lr6VCvAk8PbM8xlUsXaAwiInLklIAqAS2s3J3563cxbuEmxi7cxKzVO3CHqmUT6NusGn2aV+OEJlUok6BhP0SkYCkBDUSyEv38l3Xc9+kc0tKdf5zZigs61tED0CIiUUAJqBLQaLF19wEmLN7M2IWbmLh4M8n7UykRG0PXRpU5s10tTmlVXaPqikiBUAIaiHQlum7HPu74YBZTl2/jxBbVeODsNtQorwGKREQKMyWgSkCjUUpaOkkrtzNu0Sa+nrue1dv2kRAXw4ktqnNW+1r0alaVhLjYSIcpIkWUEtBAYahE09KdN35cwePfLiI+Noa/ntaCCzvVVWuoiEghpQQ08nWnHB135+dfd/DZrLV8MXs9W/ccpHzJeAa0qcFZ7WvTpX4lPTMqInlKCWigMFWiK7fsYegns5m6fBsnNK7CQ+e2oW6lUpEOS0REDqEEtPDUnXL0UtLSmbR0C5/NWsfoeRvYezCNmuUTObNdLc5sX4uWNcvppriIHDUloIHCVommpzvvTv+Vh75agAN3n9KMy7vVJ1Z3IUVECg0loIWr7pS8s/dgKt8t2MSomWuZsHgzqelOk2plOKt9Lc5qX1s3xkXkiCkBDRTWSnTtjn3c88kcJi7eTJva5fn3OW1oU6d8pMMSERGUgBbWulPy1rY9B/lqznpGzVrLTyu3A3BsvQqc3aE2p7WpSeUymmdURMKnBDRQmCtRd+eL2ev55xfz2br7AJd3q88dJzelnEarExGJKCWghbfulPyxZvtePvtlHaNmrmPRxmRiY4weTapwdvvanNSyOqU1rYuI5EAJaCAaKtFd+1N4YvQi3py6iqplEvjHGa0Y0KaGnscQEYkQJaCFv+6U/LNwwy7+O3Mdn81ay7qd+ylVIpazO9TmT8cdQ8ta5SIdnogUUkpAA9FUif6yegf3fjqHeet20atpVf55ViuOqVw60mGJiBQ7SkCjp+6U/JOe7iSt2s4HSav5/Jd1HEhNp+MxFbms6zGc2qaGpnQRkT9QAhqItko0NS2dt6au4olvF3MwLZ0bejXixl6NKFlCv+RFRAqKEtDoqjsl/+3Ye5CPZqzh7amrWLl1L5VKl+DCTnW59Lh6GrhIRAAloL+J1kp0w879/PurBXz2yzpqVyjJ305vwSmt1C1XRKQgKAGNzrpT8l96uvPjsi28NWUV3y3YiAN9mlXjT13r0atpNY3qL1KMKQENRHslOnX5Vu7/bB4LNyRzQuMq3H9mSxpXKxvpsEREijQloNFdd0rBWLdjHyOn/8p7P61mc/IB6lQsyaXHHcOFnepoBF2RYuhwdWdMDgftymFJNrPF+Re2HKprw8p8ccsJ/N+ZrZi9Zgf9n/6BB76YT/L+lEiHJiIiR8nM6prZODNbYGbzzGxIUH5B8D7dzDpl2r++me0zs1nB8lKmbR3NbI6ZLTWzZyzoMmNmCWb2flA+zczqF/iFSpFUq0JJ7ji5GT/+pS/PXdKBOhVL8sg3C+n20Pfc/v4sZqzaRnFr+BCR/5VtAgosc/dy2SxlgT1ZHZhNJfqYmS00s9lm9qmZVQjKVYmGKS42hiuOr8+4u3pzfsc6DPtxBX0en8DHM9aQnq5f7CIiUSwVuNPdWwBdgcFm1hKYC5wLTMzimGXu3j5YbshU/iIwCGgSLP2D8muA7e7eGHgKeCR/LkWKqxJxMZzethYjB3VjzO09ubhLXb6bv5HzXpzCgGcm8e60X9lzIDXSYYpIhOSUgJ4XxjkOt8/hKtExQGt3bwssBu7JdIwq0VyoXCaBh89ry39v6k6diiW588NfOP+lycxduzPSoYmIyBFw9/Xu/nOwngwsAGq7+wJ3XxTuecysJlDO3ad4qMnpTeDsYPNZwIhg/SOgX8aNXZG81qR6Wf7vrNZMvbcf/z6nDQD3fjqHrv8eyz9GzWXJxuQIRygiBS3bBNTdl+d0gsPtk00l+q27Z9z2mgrUye78qkRz1q5uBT658XgePb8tv27byxnPTWLox7PZnHwg0qGJiMgRCnr1dACm5bBrAzObaWYTzKxHUFYbWJNpnzVBWca21QBBfbwTqJxXcYtkpXRCHJccV4+vbj2Bj288nhNbVue96as56amJXD38J5JWbot0iCJSQHJ6BrS5mX1tZl+aWSMzG25mO8xsupm1CPdDsqlErwa+zvRelegRiokxLuxUl+/v6s3V3Rvw0Yw19Hl8PC9NWMaB1LRIhyciIrlgZmWAj4Hb3H1XNruuB+q5ewfgDuBdMysHZHUzNuMZjey2ZY5hkJklmVnS5s2bc3cBIodhZnQ8piJPDWzPlHv6csdJTZm1egfnvzSFC1+awrhFm/ScqEgRl1MX3FeAF4C3ge+Bb4CKwL+A58L5gMNVomZ2H6Fuuu8ERapE80C5xHj+dnpLvr29J8c1qMTDXy/k5KcmMnreBv1CFxGJAmYWT6jefMfdP8luX3c/4O5bg/UZwDKgKaGbtZl7GNUB1gXra4C6wWfFAeWB/2l+cvdX3L2Tu3eqWrXq0V2USBYql0ng1n5NmPSXPvzjjJas2b6Xq974iQHPTOKzX9aRmpYe6RBFJB/klICWdffP3f09IMXdR3rI54QS0WwdrhI1syuA04FLg261qkTzWMOqZRh2ZWfevLoLJWJjuP6tGVz62jQWrM/uRrqIiERS8BjJMGCBuz8Zxv5VzSw2WG9IaJyE5e6+Hkg2s67BOS8HRgWHfQZcEayfD3zvukMpEVSqRBxXdW/A+Lv78Nj5bTmYmsat782k35MTeGfaKvanqCeXSFGSUwIam2n90IqwRHYHHq4SNbP+wF+AM919b6ZyVaL5oGfTqnw9pAf/PKsV89fv4rRnfuDeT+ewdbeeDxURKYS6A5cBfTONCj/AzM4xszVAN+BLMxsd7N8TmG1mvxAaC+EGd8+4EXsj8BqwlNBN3YxHXoYBlc1sKaEeR0ML5MpEclAiLoYLOtVlzO29eOlPHalQMp77Pp1Lz0fH8fKEZezWyLkiRYJll6+Z2fWEWi93H1LeGLjZ3W/L5tgTgB+AOUBGH4p7gWeABGBrUDbV3W8ws/OAfxLqlpsG/CNoaSWY82w4UJJQBXqLu7uZJQJvEXq+dBtwUU4DJxXnybR37D3I098t4a2pqygVH8ut/ZpwxfH1KRGX030IEZHizQ4zmXZxUZzrTokcd2fysq28OH4Zk5ZuoVxiHFccX58rj69P5TIJkQ5PRHJwuLoz2wS0KFIlCks3JfOvLxYwYfFmGlQpzX0DWtCvRTWK4QDCIiJhUQKqulMi65fVO3hx/DJGz99AQlwMF3Wux3U9G1K7QslIhyYih3G4ujPXTV9m9nPehCSR0rhaWUZc3YU3rupMjMG1byZx+evTWay5uERERKQQale3Ai9d1pExt/fk9La1eHvqKno9Oo47PpiluURFokyuW0DNbGYwUm1U0l3cP0pJS+etKat4+rvF7DmYxiVd6nH7SU2pVDrbR3xFRIoVtYCq7pTCZe2Ofbz2w3Lem/4r+1PSOblldW7q05j2dStEOjQRCeRZCyjwZR7EI4VEfGwMV58QGnnu0uPq8e70X+n92Dhen7SCFA1/LiIiIoVQ7Qol+ccZrZg8tB+39m3M1OVbOfv5H7lm+E8a8V+kkMtVC2gwL2dcxvtMI+1FDd3Fzd7ijcn864v5/LBkCw2rlOZePR8qIqIWUNWdUsjtPpDKiMkreSkYLffs9rW5/cSm1KtcKtKhiRRbR9UCambXm9lGYDaQBMwIXqWIaVq9LG9e3YVhV3SC4PnQPw2bxvx1upsoIiIihVOZhDgG92nMD3/uw/U9G/HVnPX0e3I8/xg1l83JmnpOpDAJqwXUzJYA3dx9S/6HlL90Fzd8KWnpvDN1FU+PXcLOfSkM7FSXO05uSrWyiZEOTUSkQKkFVHWnRJcNO/fzzPdLeP+n1STExXDNCQ24rmdDyiXGRzo0kWLjaJ8BXQbszduQpLCLj43hyu4NmHBXH67u3oCPf15Dn8fG8/y4pexPSYt0eCIiIiJZqlE+kX+f04bv7uhF3+bVePb7pfR8dByvTlyuv2FEIizcFtAOwBvANOC3fgzufmv+hZY/dBf3yK3YsoeHvlrAt/M3UrtCSf7cvxlntqul50NFpMhTC6jqToluc9fu5NHRi5i4eDM1yydy24lNOO/YOsTFHsl4nCISjqNtAX0Z+B6YSuj5z4xFipEGVUrzyuWdePe64yhfMp4hI2dx7ouTmbFqe6RDExERETms1rXL8+bVXXj3uuOoXi6Rv3w8h5OfnsjXc9aT2ykJReTohNsCOtndjy+AePKd7uLmjbR05+Of1/DY6EVsTj7AGe1q8Zf+zahTUaPNiUjRk5sWUDOrFMZu6e6+4+iiKjiqO6UocXe+nb+Rx0YvYumm3bStU54/n9KcE5pUiXRoIkXK4erOcBPQB4FVwOf8sQuupmEp5vYcSOXlCct4eeJyHLj2hAbc1KcxZRLicjxWRCRa5DIB3Q+sA7J7PiHW3evlSXAFQHWnFEVp6c4nP6/h6e+WsHbHPro3rsyfT2lOu7oVIh2aSJFwtAnoiiyK3d0b5kVwBUmVaP5Yt2Mfj36zkP/OWkeVMgncdXJTLuhUl9gYPR8qItEvlwnoTHfvcLT7FCaqO6UoO5CaxjtTf+W5cUvZtucgp7auwZ0nN6NxtTKRDk0kqh1VAlqUqBLNX7NW7+CBL+aTtGo7zWuU5a+ntVSXFhGJerlMQBPdff/R7lOYqO6U4mD3gVRe+2E5r05czr6UNC497hjuPLkpFUqViHRoIlHpaAchyuqENY4uJCmK2tetwIc3dOP5S45l94FU/jRsGtcM/4llm3dHOjQRkQKRVWJ56HOh0ZR8ihQXZRLiuO3Epkz8cx/+1PUY3pm2ij6Pj+fdab+Sll68GmxE8tPRjD09LM+ikCLFzDitbU2+u6MXf+nfnGkrtnHKUxO5/7N5bN9zMNLhiYjkKzPrbmYLzGyemR1nZmOAJDNbbWbdIh2fiGSvcpkE/nlWa768tQdNqpfl3k/ncPbzP2rUf5E8oi64ku+27D7Ak2MWM3L6r5RNjOfWfk24rOsxlIjT3FsiEh1y2QV3OnANUIbQ4H1nu/skMzsWeNbdu+djqPlCdacUV+7O57PX8+8vF7Bh137OPbY2Q09tTrWyiZEOTaTQO6IuuGZWKbsl/8KVoqRKmQT+fU4bvhrSg7Z1yvOvL+ZzytMT+XbeBs29JSJFUby7z3H3KcBmd58E4O4/AyUjG5qI5IaZcWa7Woy9sxc39W7EF7+sp+/jE3h14nJS0tIjHZ5IVMqpCWoGkBS8HrroVqjkSvMa5Xjz6i68cWVnYgwGvTWDS16dxrx1OyMdmohIXspct95zyDaNZiIShUonxPHn/s0ZfXtPujSoxINfLaD/0xP5YcnmSIcmEnXUBVciIiUtnfem/8pTYxazY18KF3Ssw10nN6NaOXVpEZHCJ5ddcM8EvnP3vYeUNwLOc/dH8yPG/KS6U+SPxi7YyD+/mM+qrXs5pVV1/npaS+pWKhXpsEQKlaOehiWoUHsGb8e7+xd5GF+BUSVauOzcl8Jz3y9h+OSVxMfGcGOvRlzXsyGJ8bGRDk1E5De5SUAPc3wNd9+QlzEVJNWdIv9rf0oawyat4Lnvl5Luzg29GnFj70b6G0YkcFTTsJjZw8AQYH6wDDGzh/I2RCmOypeM577TWjLm9l70bFKVJ8Yspu/j4/nvzLWka8hzESk6vop0ACKStxLjYxncpzFj7+zFSS2r85+xS+j3xAS+mbteY1yIZCPcYUgHACe5++vu/jrQHzgt/8KS4qZ+ldK8dFlHRg7qSqUyJbjt/Vmc8+JkZqzaFunQRETygoW1k1ldMxuXaRqXIUH5BcH7dDPrlGn/k8xshpnNCV77Zto23swWmdmsYKkWlCeY2ftmttTMpplZ/Ty+VpFipVaFkjx3ybG8d11XyibGccPbP3PZsOks3ZQc6dBECqXczINRIdN6+TyOQwSArg0r89ngE3j8gnZs2LmP816cwuB3f2b1tr05HywiUni9GuZ+qcCd7t4C6AoMNrOWwFzgXGDiIftvAc5w9zbAFcBbh2y/1N3bB8umoOwaYLu7NwaeAh7J/eWIyKG6NarMF7ecwP1ntGT2mh30f/oHHvxyPsn7UyIdmkihEm4C+hAw08yGm9kIQqPg/jv/wpLiLCbGOL9jHcbd1Ztb+zVh7IKN9HtyAo98s1C/xEUkqphZRTNrC0w1s2ODuUAPy93XB9O14O7JwAKgtrsvcPdFWew/093XBW/nAYlmlpBDWGcBI4L1j4B+ZhZWC62IZC8uNoYruzdg3F29Ob9jHV6btII+j0/gs1/WqVuuSCCsBNTd3yN0J/aTYOnm7iOzOyabbkSVzGyMmS0JXitmOuaeoEvQIjM7JVN5x6B70VIzeyajolQ3oqKtVIk47jipKePu6s3pbWry4vhl9Hl8PO9O+5VUzb0lIoWcmf0LmA08AzwRLI/n4vj6QAdgWpiHnAfMdPcDmcreCLrf/i1TklkbWA3g7qnATqByFp8/yMySzCxp82ZNNSGSG5XLJPDweW0ZNbg7tSskcut7M7nuzRls2Lk/0qGJRFxuuuBWDV5jgePN7Nwc9j9cN6KhwFh3bwKMDd4TbLsIaEXoGdMXzCxjGLEXgUFAk2DpH5SrG1ExULN8SZ4c2J7Pbu5OgyqluffTOZz2zCTNvSUihd2FQCN37+3ufYKlb45HAWZWBvgYuM3dd4WxfytCdeD1mYovDbrm9giWyzJ2z+IU/9M04+6vuHsnd+9UtWrVLA4RkZy0rVOBT27qzl9Pa8GkpZs56ckJvDf9V7WGSrEW7ii4rwOvE7q7ekawnJ7dMYfrRsQfu/6MAM4O1s8CRrr7AXdfASwFuphZTaCcu0/x0P/WNw85Rt2Iiom2dSrwwfXdePHSY9mbksplw6Zz1Rt6yF9ECq25/HH8hLCYWTyh5PMdd/8kjP3rAJ8Cl7v7soxyd18bvCYD7wJdgk1rgLrBsXGExnXQiG8i+SQ2xri2R0NG39aT1rXLc88nc7jk1Wms3LIn0qGJRERcmPt1dfeWR/ohh3Qjqu7u6yGUpGaMykcoOZ2a6bA1QVlKsH5oecYxv3UjMrOMbkRbDvn8QYRaUKlXr96RXoYUAmbGqW1q0rdFNUZMXsmzY5dyytM/cOlx9bjtxKZUKl0i0iGKiGTIGD9hLvBbt1h3P/NwBwQ3UYcBC9z9yZw+wMwqAF8C97j7j5nK44AK7r4lSGhPB74LNn9GaMCiKcD5wPeu5hiRfHdM5dK8e91xvP/Tah78cgH9/zORO09qxtUnNCA2Ru0nUnyEm4BOMbOW7j4/tx9waDeibBooD9clKLuuQmF3IwJegdBk2jnFLIVfQlwsg3o24rxj6/D0d0t4Z9qvfDpzLbf0bcwVx9cnIU6TQItIxI0g1C12DhDug+vdCXWVnWNms4Kye4EE4FlCj8N8aWaz3P0U4GagMfA3M/tbsP/JwB5gdJB8xhJKPjNG4h0GvGVmSwm1fF50xFcoIrliZlzUpR69m1Xjr/+dy4NfLeCL2et45Py2NK9RLtLhiRSIcBPQEYSS0A2E7uIa4O7eNruDDtONaKOZ1QxaP2sCGcPC/9YlKFAHWBeU18miPPMxa9SNqHiqXCaBf53dmsu7HcO/v1rAv79ayFtTV3HPqS04tXUN1CNbRCJoi7s/k5sD3H0Sh58z9NMs9n8AeOAw+3c8zGfsBy7ITVwikrdqlE/k1cs78sXs9dz/2TxOf2YSN/VpzOA+jXQTXYq8cAchep3QHdn+/P785xnZHZBNN6KMrj8Er6MylV8UjGzbgNBgQ9OD7rrJZtY1OOflhxyTcS51IyrGmlQvyxtXdeHNq7tQKj6Om975mQtfnsIvq3dEOjQRKb5mmNlDZtYtYwqWnKZhEZHiw8w4o10txtzRizPa1eKZsUs4/ZlJ/Pzr9kiHJpKvLJx8zcy+D3fkvkzHnAD8wB+7Ht1L6DnQD4B6wK/ABe6+LTjmPuBqQiPo3ubuXwflnYDhQEnga+AWd3czSyQ06XYHgm5E7r48u7g6derkSUlJubkUiTKpael8kLSGJ8csYsvug5zToTZ3n9KMWhVKRjo0EYlSZjbD3Tvl8phxWRR7buvTwkB1p0j+G7dwE/d+OocNu/ZzdfcG3HlyU0qVCLezokjhc7i6M9wE9AVCI/l9zh8HUshxdL7CRpVo8ZG8P4UXxy/jtUkriDEY1KMh1/dqROkE/TIXkdw5kgS0KFHdKVIwkven8Og3i3hr6irqVirJw+e2pXvjKpEOS+SIHK7uDLcLbklCiefJhDkNi0iklU2M58/9m/P9nb04qWUNnvl+Kb0fH88HP60mLV09tUWk4KkLrohkp2xiPP86uzXvD+pKXEwMl742jb98NJud+1IiHZpIngmrBbQo0V3c4mvGqu088OV8Zv66gxY1y/G301pwvO4qikgY8qoF1Mxedffr8iKmgqS6U6Tg7U9J4+nvlvDqD8upXLoE/zq7Nae0qhHpsETCdkQtoMH8mTmdOMd9RAqDjsdU5JMbj+fZizuwa18Kl7w2jWtHJLF88+5IhyYixUQ0Jp8iEhmJ8bEMPbU5/72pO5XLJHD9WzMY/M7PbNl9IOeDRQqxbFtAzWw5cFd2xwP/dPdWeR1YftFdXIHQXcXXf1zBC+OWsT8ljT91PYYh/ZpQsXSJSIcmIoVQblpAc+pm6+4/501UBUd1p0hkpaSl88rE5fznuyWUKxnHo+e3pW/z6pEOSyRbRzQIkZm9Eca5d7r7bUcRW4FSJSqZbU4+wFPfLWbk9F8pmxjPrf2acFnXYygRF+7j0SJSHOQyAc1q9NsMGgVXRI7Yog3JDBk5k4Ubkrms6zHcO6AFJUto3lApnI5qFNyiRJWoZGXRhmQe+HI+PyzZQv3KpbhnQAtOblmd0NSzIlLcaRRc1Z0ihcX+lDQeH72I1yatoHG1Mjw9sD2ta5ePdFgi/+NoR8EVKdKa1SjLm1d34Y2rOhMXG8P1b83g4lenMnftzkiHJiJRyszizexWM/soWG42s/hIxyUi0S0xPpa/nt6St67pwq59KZzzwo+8PGEZ6RrhX6KEElCRgJnRp1k1vhnSg3+d3ZrFG3dzxnOTuOvDX9i4a3+kwxOR6PMi0BF4IVg6BmUiIketR5OqjL6tJ/2aV+ehrxdy6WvTWLdjX6TDEsmRuuCKHMau/Sk8//1S3vhxJbExxvW9GjKoZ0NKlYiLdGgiUsCOpAuumf3i7u1yKosGqjtFCi9358OkNdz/+TziYoyHzm3LaW1rRjoskaPrgmtmCWZ2iZnda2Z/z1jyPkyRwqNcYjz3DGjBd3f0om/zajz93RL6Pj6Bj2esUTcXEQlHmpk1ynhjZg2BtAjGIyJFkJlxYee6fHVrDxpULcPgd3/mzg9+IXl/SqRDE8lSuF1wRwFnAanAnkyLSJFXr3Ipnr/0WD66oRvVyyVw54e/cObzk5i2fGukQxORwu1uYJyZjTezCcD3wJ0RjklEiqj6VUrz0Q3duLVvYz6duYYBz/zAjFXbIx2WyP8Iqwuumc1199YFEE++UzciORrp6c5nv6zjkW8Wsn7nfvq3qsHQU5tTv0rpSIcmIvnoSEfBNbMEoBmhebMXuntUziCvulMkuiSt3MZt789i/c793NynMbf0bUxcrIZ+kYJ1tKPgTjazNnkck0jUiYkxzu5Qm+/v7M1dJzdl4pLNnPTUBB74Yj4796qri4j8zsxigVOA3kA/YLCZ3RHRoESkWOhUvxJfDenBWe1q8Z+xS7jg5Sms2qrOi1I4hJuAngDMMLNFZjbbzOaY2ez8DEykMCtZIpab+zZh/F29ObdDHYb9uILej49jxOSVpKSlRzo8ESkcPgeuBCoDZTMtIiL5rlxiPE8ObM+zF3dg6abdDPjPD3yYtJriNgCpFD7hdsE9Jqtyd1+V5xHlM3Ujkvwwf90uHvhyPpOXbaVh1dLcN6AFfZtXw8wiHZqI5IEjHAV3tru3za+YCpLqTpHotnbHPu54fxbTVmxjQJsa/PucNlQoVSLSYUkRd1RdcINEswJwRrBUiMbkUyS/tKxVjneuPY5hV4T+j10zIok/DZvG/HW7IhyZiETQ12Z2cqSDEBGpXaEk717Xlb/0b86Y+Rvp//QPTF66JdJhSTEV7jQsQ4B3gGrB8raZ3ZKfgYlEGzOjX4vqjL6tJ/ef0ZJ563Zx2rM/MPTj2WxK3h/p8ESk4E0FPjWzfWa2y8ySzUx3pUQkImJjjBt7N+LTm7pTKiGWS16bxsNfLyRVjw5JAQu3C+5soJu77wnelwamRGPXInUjkoKyc28Kz36/hBFTVhIfG8NNvRtxbY+GJMbHRjo0EcmlI+yCuxw4G5jjUf7QlepOkaJl38E0/vnFfN6b/ivHNajEs5d0oFrZxEiHJUXM0Y6Ca/xx8uy0oExEDqN8qXj+enpLvr29Fz2aVOHxbxfT9/HxjJq1lvT0qP5bVETCswSYm5vk08zqmtk4M1tgZvOCHkiY2QXB+3Qz63TIMfeY2dJgoMBTMpV3DAYNXGpmz1jwULqZJZjZ+0H5NDOrnzeXKyLRomSJWB46tw1PXtiOX9bs4LRnJjFV85tLAQk3AX0DmGZm95vZ/YS6FQ3Lt6hEipAGVUrz8mWdGDmoK5XKlGDIyFmc8+JkklZui3RoIpK/1gPjgwTxjowlh2NSgTvdvQXQldDULS2BucC5wMTMOwfbLgJaAf2BF4LpXwBeBAYBTYKlf1B+DbDd3RsDTwGPHOV1ikiUOvfYOowafAJlE+K45NWpvDh+mW6SS74LdxCiJ4GrgG3AduAqd386H+MSKXK6NqzMZ4NP4PEL2rFh5z7Of2kKg9/5mdXb9kY6NBHJHyuAsUAJwpyGxd3Xu/vPwXoysACo7e4L3H1RFoecBYx09wPuvgJYCnQxs5pAOXefErTAvkmoO3DGMSOC9Y+AfhmtoyJS/DSrUZbPbjmBU9vU5JFvFjLorRma21zyVVx2G82snLvvMrNKwMpgydhWyd3VhCOSCzExxvkd6zCgTQ1embiclycsZ8z8jVx1Qn0G92lMucT4SIcoInnE3f/vaI4PusZ2AKZls1ttQr2SMqwJylKC9UPLM45ZHcSYamY7Cc1V+ochMc1sEKEWVOrVq3eklyEiUaBMQhzPXdyBzsdU5MGvFnD6cz/w4qUdaV27fKRDkyIopxbQd4PXGUBSpiXjvYgcgVIl4rjtxKaMu6s3Z7avxSsTl9P7sfG8NXWVRqMTiXLBoypHtY+ZlQE+Bm5z9+xGzs2q5dKzKc/umD8WuL/i7p3cvVPVqlWzC1dEigAz48ruDXj/+m6kpjnnvjiZ96b/SpSPoSaFULYtoO5+evDaoGDCESleapRP5PEL2nHl8fX51xfz+dt/5/Lm5JXcd1oLejerFunwROTIXJvDdCtG6LnN+7PcaBZPKPl8x90/yeGz1gB1M72vA6wLyutkUZ75mDVmFgeUJ/SIjYgIx9aryJe39mDIyJnc88kcflq5jQfPbkPJEhrFX/JGuPOAjg2n7JDtr5vZJjObm6nsfTObFSwrzWxWUF4/mCctY9tLmY7RKH5S5LWuXZ6Rg7ry8mUdSUlL58o3fuLy16ezaENypEMTkdx7lT8+83noUibY538EddwwYEEw/kJOPgMuCurEBoQGG5ru7uuBZDPrGpzzcmBUpmOuCNbPB76P9mliRCRvVSpdguFXdeG2E5vw6cy1nPPCjyzfvDvSYUkRkdMzoIlAKaCKmVXk92475YBaOZx7OPAcoYEPAHD3gZnO/QSwM9P+y9y9fRbnyRjFbyrwFaFR/L4m0yh+ZnYRoVH8BmZxvEhUMDNOaVWDPs2q8eaUlTwzdgmn/mciF3Wpxx0nNaVKmYRIhygiYTjKZz+7A5cBczJu0gL3AgnAs0BV4Eszm+Xup7j7PDP7AJhPaATdwe6eMW3ajYTq4pKE6s2vg/JhwFtmtpRQy+dFRxGviBRRsTHGbSc2pUO9itw2ciZnPvcjj57flgFtakY6NIlylt1Nz2D+sdsIJZtr+T0B3QW86u7PZXvyUKvkF+7e+pByA34F+rr7kmz2qwmMc/fmwfuLgd7ufr2ZjQbud/cpQReiDUDVnO7iajJtiRbb9xzkP2OX8PbUVSTGxzK4T2Ou6l6fxHh1gREpaIebTLu4UN0pUryt27GPwe/+zMxfd3B19wYMPbU5JeLCnc1RiqvD1Z3Z/stx9/8Ez3/e5e4N3b1BsLTLKfnMQQ9go7svyVTWwMxmmtkEM+sRlNUmzFH8CLWmVs7qw8xskJklmVnS5s2bjyJskYJTsXQJ7j+zFaNv70nXhpV55JuF9HtiAp//sk4DAoiIiEiBqVWhJO8P6sZV3evz+o8ruOiVKazfuS/SYUmUCnce0GfNrLWZXWhml2csR/G5FwPvZXq/Hqjn7h2AO4B3zawceTCKXxC/RvKTqNWoahleu6IT71x7HOVKxnPLezM578XJzPx1e6RDExERkWKiRFwM/zijFc9d0oFFG5I57ZlJ/LBEDTuSe+EOQvQPQs+ePAv0AR4FzjySDwy6y54LvJ9RFkygvTVYnwEsA5oS3ih+GefUKH5SpHVvXIUvbjmBR89ry+rt+zjnhcnc+t5M1mzfG+nQRCQLZlbVzO41s1eCgfleN7PXIx2XiMjROL1tLT675QSqlCnB5a9P5z/fLSE9XT2zJHzhdt4+H+gHbHD3q4B2hAZEOBInAgvd/beutUElHRusNyQ0it9yjeIn8kexMcaFnesy7q7e3NK3MaPnbaDfExN4bPRCdh9IjXR4IvJHowjdHP0O+DLTIiIS1RpVLcN/B3fn7Pa1eeq7xVw5/Ce27TkY6bAkSoSbgO5z93QgNegauwlomN0BZvYeMAVoZmZrzOyaYNNF/LH7LUBPYLaZ/QJ8BNzg7hmtmTcCrwFLCbWMZh7Fr3Iwit8dwNAwr0Uk6pVJiOPOk5sx7q7eDGhTk+fHLaP3Y+N4b/qvpOkupEhhUcrd/+LuH7j7xxlLpIMSEckLpUrE8eSF7fj3OW2YumwrZz43iYUbspsCWSQk21Fwf9vJ7AVCw8BfBNwJ7AZmBa2hUUUj+UlRNGv1Dh74Yj5Jq7bTvEZZ7jutBT2a6HlnkbxyJKPgmtkDwGR3/yqfwiowqjtFJDu/rN7BdW8msedAKk8NbM/JrWpEOiQpBA5Xd+aYgAZdX+u4++rgfX2gnLvPzo9A85sqUSmq3J2v527goa8XsHrbPvo0q8p9p7WgcbWykQ5NJOodYQKaDJQGDgIpQbG7e7m8ji+/qe4UkZxs3LWfQW8mMXvtTu46uRk39W5EKI2Q4uqIpmGBUE0J/DfT+5XRmnyKFGVmxoA2Nfnujl7cO6A5SSu3c8rTP/D3UXP1XIZIBLh7WXePcffEYL1sNCafIiLhqF4ukfev78YZbWvx2OhF3Pb+LPanpEU6LCmEwn0GdKqZdc7XSEQkTyTExTKoZyPG392bS7rU451pv9LrsXG8OnE5B1JVEYgUJDM708weD5bTIx2PiEh+SoyP5T8XtefuU5oxatY6Br48hY279kc6LClkwk1A+wBTzGyZmc02szlmplZQkUKscpkE/nV2a74Z0oNOx1Tkwa8WcNKTE/l6zno0YLRI/jOzh4EhwPxgGRKUiYgUWWbG4D6NeeWyjizZtJszn5vE7DU7Ih2WFCLhDkJ0TFbl7r4qzyPKZ3qORYqriYs388CX81m8cTdd6lfir6e3oG2dCpEOSyQqHOEzoLOB9sEo8gTTjc1097b5EWN+Ut0pIkdiwfpdXDsiiS27D/DYBe04s12tSIckBeiInwENPODuqzIvwAN5G6KI5KeeTavy1a09+Pc5bVi+ZTdnPvcjd7w/i/U790U6NJGirEKm9fKRCkJEJBJa1CzHZzd3p12dCtz63kweH72IdE0XV+yFm4C2yvwmuIvbMe/DEZH8FBcbwyXH1WPcXb25qXcjvpiznj6Pj+fJbxex50BqpMMTKWoeAmaa2XAzGwHMAP4d4ZhERApU5TIJvH3tcVzUuS7PjVvKDW/P0N8cxVy2CaiZ3RMMI9/WzHYFSzKwCRhVIBGKSJ4rmxjPn/s3Z+wdvTipZQ2e+X4pfR4fzwdJq0nTnUmRPOHu7wFdgU+CpZu7j4xsVCIiBa9EXAwPnduGf5zRku8WbOS8FyezetveSIclEZJtAuruD7l7WeAxdy8XLGXdvbK731NAMYpIPqlbqRTPXtyBj288ntoVS/Lnj2ZzxrOTmLxsS6RDE4laZtY8eD0WqAmsAVYDtYIyEZFix8y4qnsDhl/VhbU79nHW8z8yfcW2SIclERDWIEQAZlYbOAaIyyhz94n5FFe+0UAKIllzdz6fvZ5Hvl7I2h37OLFFde4d0JyGVctEOjSRiMvNIERm9oq7DzKzcVlsdnfvm8fh5TvVnSKSl5Zt3s11I5JYvX0vD5zdmoGd60U6JMkHh6s7wx0F92HgIkLDyGdMJOjufmaeRlkAVImKZG9/Shqv/7iCF8YtY39KGpd1O4Yh/ZpQoVSJSIcmEjFHOApuorvvz6ksGqjuFJG8tnNvCje/9zM/LNnCVd3rc9+AFsTFhjs8jUSDox0F9xygmbsPcPczgiXqkk8RyVlifCw39W7MuLt6c2HnuoyYvJJej41n2KQVHExNj3R4ItFkcphlIiLFTvlS8bxxZWeu6l6fN35cyVXDf2LnvpRIhyUFINwEdDkQn5+BiEjhUrVsAv8+pw1fDelB2zrl+dcX8znl6Yl8O28D4XbdFymOzKyGmXUESppZBzM7Nlh6A6UiG52ISOERFxvDP85oxcPntmHq8q2c8/yPLN+8O9JhST6Ly3kXAPYCs8xsLHAgo9Ddb82XqESk0GheoxxvXt2F8Ys38+CXCxj01gy6NqzEX09rSevamtZQJAunAFcCdYAnM5UnA/dGIiARkcLsoi71aFClNDe+8zNnP/8jz196LD2aVI10WJJPwn0G9Iqsyt19RJ5HlM/0HIvIkUtNS+e96b/y1HdL2L73IOcdW4e7T2lG9XKJkQ5NJF8d4TOg57n7x/kVU0FS3SkiBWH1tr1c92YSSzft5qFz23BBp7qRDkmOwlENQhScoCRQz90X5XVwBUmVqMjR27kvhRfGLeWNH1cSG2Pc0KsR1/VsQKkS4XaqEIkuR5KABsedBrQCfrtL4+7/zGb/usCbQA0gHXjF3f9jZpWA94H6wErgQnffbmaXAndnOkVb4Fh3n2Vm4wlNA7Mv2Hayu28ys4TgMzoCW4GB7r4yu+tQ3SkiBSV5fwo3vv0zk5Zu4Y6TmnJL38aYWaTDkiNwVIMQmdkZwCzgm+B9ezP7LE8jFJGoUb5kPPcMaMF3d/SiT/OqPPXdYvo+PoGPZ6whPV3Ph4oAmNlLwEDgFsCACwhNZ5adVOBOd28BdAUGm1lLYCgw1t2bAGOD97j7O+7e3t3bA5cBK919VqbzXZqx3d03BWXXANvdvTHwFPDI0V+tiEjeKJsYz+tXdubcDrV5csxi7v10DqlpGgSxKAl3EKL7gS7ADoCgcmuQLxGJSNSoV7kUL1zakQ9v6Ea1cgnc+eEvnPX8j0xbvjXSoYkUBse7++WEkr3/A7oB2fYnc/f17v5zsJ4MLABqA2cBGY+9jADOzuLwi4H3wogr87k+AvqZmhdEpBApERfDExe246bejXhv+moGvTWDvQdTIx2W5JFwE9BUd995SJmaOUQEgM71K/Hfm7rz9MD2bNl9gIGvTOWGt2awcsueSIcmEkkZ833uNbNaQAq5uHlrZvWBDsA0oLq7r4dQkgpUy+KQgfxvAvqGmc0ys79lSjJrA6uDc6UCO4HKWXz+IDNLMrOkzZs3hxu2iEieMDP+3L85/zq7NeMXbeLiV6ayZfeBnA+UQi/cBHSumV0CxJpZEzN7Fs1lJiKZxMQYZ3eozfd39ubOk5oycclmTnpqAg98MV/zeklx9bmZVQAeA34m9OxmOC2UmFkZ4GPgNnffFcb+xwF73X1upuJL3b0N0CNYLsvYPYtT/M9NZXd/xd07uXunqlU1GqWIRMZlXY/hpT91ZNHGZM57cbJubhcB4SagtxAaROEA8C6hu6W35VNMIhLFSpaI5ZZ+TRh/V2/O7VCHYT+uoPdj4xgxeSUpeoZDigkziyH0zOaOYCTcY4Dm7v73MI6NJ5R8vuPunwTFG82sZrC9JrDpkMMu4pDk1t3XBq/JhOruLsGmNQRdgc0sDigPbMv1RYqIFJCTW9Xg3eu6krw/lXNfnMzMX7dHOiQ5CmEloO6+193vc/fOwfJXd9+f85EiUlxVK5fII+e35ctbetCiZjn+8dk8+j89ke8XbiTc0bdFopW7pwNPZHp/IItHWf5H0E12GLDA3TPPIfoZkDEl2hXAqEzHxBAa4GhkprI4M6sSrMcDpwNzszjX+cD3rv+UIlLIHVuvIh/feDxlEuK4+NWpjJm/MdIhyREKdxTcMUE3ooz3Fc1sdL5FJSJFRsta5Xjn2uN47fJOuMPVw5O4bNh0FqzPsVehSLT71szOy+UAP90JdZXtGzy7OcvMBgAPAyeZ2RLgpOB9hp7AGndfnqksARhtZrMJjWK/Fng12DYMqGxmS4E7CEbUFREp7BpUKc3HNx5P0+pluf6tJN6euirSIckRCGseUDOb6e4dcio7ZPvrhO64bnL31kHZ/cB1QMZoBve6+1fBtnsIDQ2fBtzq7qOD8o7AcKAk8BUwxN39SOYxA81lJhJJKWnpvDN1FU+PXcKufSlc2Kkud5zclGplE3M+WCSCjmQeUDNLBkoTmlplP6FnL93dy+VDiPlKdaeIFCZ7D6Yy+J2fGbdoM4P7NOKuk5tprtBC6KjmAQXSzaxeppMdQ86j4A4H+mdR/lSmOckyks+WhJ5faRUc84KZxQb7vwgMApoES8Y5NY+ZSJSJj43hyu4NmHBXH67q3oCPf15Dn8fG8/y4pexPSYt0eCJ5yt3LunuMu5dw93LB+6hLPkVECptSJeJ49fJOXNS5Ls+PW8adH/7CwVSNMxEtwk1A7wMmmdlbZvYWMBG4J7sD3H0i4Q9qcBYwMnhGZgWwFOgSDLRQzt2nBM+nvMnvc59pHjORKFW+VDx/O70l397eixOaVOGx0Yvo+/h4Rs1aq+dDpcgws7HhlImISO7Fxcbw0LltuOOkpnzy81quGfETyfs16n40CHcQom+AY4H3gQ+AjhldZI/AzWY228xeN7OKQdlvc5IF1gRltYP1Q8v/cEx285iB5jITKawaVCnNy5d1YuSgrlQqU4IhI2dxzguTmbFKA3JK9DKzRDOrBFQJxkyoFCz1gVoRDk9EpMgwM27t14RHz2/L5GVbGfjyVDbu0jiphV24LaAQGtBgG6FEr6WZ9TyCz3sRaAS0B9bz+wiBh5uTLLu5ysKaxww0l5lIYde1YWU+G3wCj1/QjvU793Hei1MY/O7PrN62N9KhiRyJ64EZQPPgNWMZBTwfwbhERIqkCzvV5fUrO7Ny6x7OfWEySzYmRzokyUa4o+A+AvxIqCvu3cFyV24/zN03untaMDz9q2QxJ1mgDrAuKK+TRfkfjtE8ZiLRLybGOL9jHcbd1ZvbTmzC9ws20e+JCTz09QJ2qUuNRBF3/4+7NwDucveG7t4gWNq5+3ORjk9EpCjq1bQqH1zfjQOp6Zz34mSmr1BaUFiF2wJ6NtDM3U9z9zOC5czcfljGJNqBc/jjnGQXmVmCmTUgNNjQdHdfDySbWdfg+c7L+X3uM81jJlIElSoRx20nNmXcXb05o10tXp6wnD6PjeftqatITdMAAxI93P1ZMzvezC4xs8szlkjHJSJSVLWuXZ5PbzqeKmUT+NOwaXw1Z32kQ5IshJuALgfic3NiM3sPmAI0M7M1ZnYN8KiZzQnmJesD3A7g7vMIPVs6H/gGGOzuGUNi3gi8RmhgomXA10G55jETKcJqlE/kiQvb8fnNJ9CoWhn++t+5nPqfHxi/aFOkQxMJSzBo3+PACUDnYMnVVC4iIpI7dSuV4uMbjqdN7fIMfvdnXp+0ItIhySHCnQf0Y6AdMBY4kFHu7rfmX2j5Q3OZiUQfd2f0vI089PUCVm3dS6+mVbnvtBY0rV420qFJMXGE84AuAFoWhd45qjtFJNrsT0ljyMiZjJ63kdtObMKQfk00V2gBO1zdGRfm8Z8Fi4hIgTMz+reuQd/m1XhzykqeGbuE/k9P5OIu9bj9pKZUKZMQ6RBFsjIXqEFo0D0RESlAifGxPH/JsQz9ZA5Pf7eE3ftTue+0FkpCC4GwElB3H2FmJYCmQdEid9eoICJSoErExXBtj4acd2wd/jN2CW9PXcWoWesY3KcxV3WvT2J8bKRDFMmsCjDfzKbzx95DuR5DQUREci8uNoZHz2tLmYQ4Xpu0gt0HUnnwnDbExigJjaSwElAz6w2MAFYSmv6krpld4e4T8y0yEZHDqFi6BPef2YrLuh3DQ18t4JFvFvLOtFUMPbU5p7WpqbubUljcH+kARESKu5gY4x9ntKRsYhzPfr+U3QdSeWpge+JjczMbpeSlcH/yTwAnu3svd+8JnAI8lX9hiYjkrFHVMrx2RWfeufY4yiTEcfO7MznvxcnM/HV7pEMTwd0nELpxGx+s/wT8HNGgRESKITPjzpObMfTU5nwxez03vDWD/SlpOR8o+SLcBDTe3RdlvHH3xeRyVFwRkfzSvXEVvry1B4+c14bV2/dxzguTGTJyJmt37It0aFKMmdl1wEfAy0FRbeC/EQtIRKSYu6FXIx44uzXfL9rEVW/8xO4DqZEOqVgKNwFNMrNhZtY7WF4FZuRnYCIiuREbYwzsXI9xd/Xmlr6N+WbuBvo+Pp7HRi9UBSORMhjoDuwCcPclQLWIRiQiUsz9qesxPHVhe6av3MafXpvGjr0HIx1SsRNuAnojMA+4FRhCaL7OG/IrKBGRI1UmIY47T27GuLt6M6BNTZ4ft4zej41n5PRfSUuP+tkwJLoccPff/rIxszhA/whFRCLs7A61efHSY5m/bhcXvTKVTcn7Ix1SsRJuAhoH/Mfdz3X3c4BnAA03KSKFVq0KJXlqYHv+O7g79SuXYugnczjtmR+YtGRLpEOT4mOCmd0LlDSzk4APgc8jHJOIiAAnt6rB61d2ZtXWvQx8eaoe2ylA4SagY4GSmd6XBL7L+3BERPJW+7oV+PCGbrxw6bHsOZjKn4ZN4+rhP7F0U3KkQ5OibyiwGZgDXA98Bfw1ohGJiMhvTmhShbev7cKW3Qe44MXJLN+8O9IhFQvhJqCJ7v7bNxKsl8qfkERE8paZMaBNTcbc3ot7Tm3OTyu2ccrTP/CPUXPZtkfPfki+KQm87u4XuPv5wOv88WauiIhEWMdjKjFyUFcOpKZz4ctTWLB+V6RDKvLCTUD3mNmxGW/MrCOgdmoRiSqJ8bFc36sR4+/uzSVd6vH2tF/p9dg4Xp24nAOpGo5d8px6D4mIRIFWtcrz/vXdiI+NYeDLUzSdWz4LNwG9DfjQzH4wsx+A94Gb8y0qEZF8VLlMAv86uzXfDOlBp2Mq8uBXCzj5qYl8M3c97hojRvKMeg+JiESJxtXK8MH13ahYugSXvjaNycs0ZkR+CSsBdfefgOaERsO9CWjh7pqGRUSiWpPqZXnjqi68eXUXEuJiuOHtnxn48lRmr9kR6dCkaFDvIRGRKFK3Uik+vL4bdSqW5Mo3fmLsgo2RDqlICrcFFHdPcfe57j7H3VPyMygRkYLUs2lVvrq1B/8+pw3Lt+zmzOd+5I73Z7F+p3IFOSq3kcveQ2ZW18zGmdkCM5tnZkOC8kpmNsbMlgSvFYPy+ma2z8xmBctLmc7V0czmmNlSM3vGzCwoTzCz94PyaWZWP5+uX0Qk6lQrl8j7g7rRvEZZrn9rBp//si7SIRU5YSegIiJFWVxsDJccV49xd/Xmxt6N+GLOevo8Pp4nxyxmz4HUSIcnUegIew+lAne6ewugKzDYzFoSGlF3rLs3IfRs6dBMxyxz9/bBknmO7heBQUCTYOkflF8DbHf3xsBTwCNHc50iIkVNxdIleOfa4zj2mIrcOnImI6f/GumQihQloCIimZRNjOcv/Zsz9o5enNSyBs+MXUKfx8fzQdJq0tL1fKjkWmegLdABuNjMLs9uZ3df7+4/B+vJwAKgNnAWMCLYbQRwdnbnMbOaQDl3n+KhB5vfzHRM5nN9BPTLaB0VEZGQsonxjLiqCz2bVGXoJ3N47YflkQ6pyAgrAbWQP5nZ34P39cysS/6GJiISOXUrleLZizvw8Y3HU6tCSf780WzOeHaSBiWQsJnZW8DjwAmEEtHOQKdcHF+fUOI6Daju7ushlKQC1TLt2sDMZprZBDPrEZTVBtZk2mdNUJaxbXVwrlRgJ1A5VxcnIlIMlCwRy6uXd+K0NjV54MsFPP3dYg1WmAfiwtzvBSAd6Av8E0gGPiZUmYqIFFkdj6nIpzcdz+ez1/PI1wu55NVpnNSyOvec2pyGVctEOjwp3DoBLf0I/loxszKE6tnb3H1XNg2U64F67r41GOTov2bWCsjqgIw4stuWOYZBhLrwUq9evVxegYhI0VAiLoZnLu5AqRKxPP3dEtId7jipaaTDimrhdsE9zt0HA/sB3H07UCLfohIRKUTMjDPb1WLsnb34c/9mTFm2lZOfmsg/P5/Pjr0HIx2eFF5zgRq5PcjM4gkln++4+ydB8cagW21G99pNAO5+wN23BuszgGVAU0ItnnUynbYOkDGSxhqgbnCuOKA8sO3QONz9FXfv5O6dqlatmtvLEBEpMmJjjEfOa8vATnV5ZuwSnhqzONIhRbVwE9AUM4sluENqZlUJtYiKiBQbifGx3NS7MePu6s0FneoyfPIKej02ntcnreBgqn4lyv+oAsw3s9Fm9lnGkt0BwbOYw4AF7v5kpk2fAVcE61cAo4L9qwb1M2bWkNBgQ8uDbrrJZtY1OOflGccccq7zge+PpJVWRKQ4iYkxHjq3DRd2qsN/lIQelXC74D4DfApUM7MHCVVYf823qERECrGqZRN46Nw2XHH8MTz45QL++cV83pq6intObc5JLauj8VwkcP8RHNMduAyYY2azgrJ7gYeBD8zsGuBX4IJgW0/gn2aWCqQBN7h7RmvmjcBwoCTwdbBAKMF9y8yWEmr5vOgI4hQRKXZiYoyHz22LO/xn7BLM4LYT1R03tyzcm55m1hzoR+jZkbHuviA/A8svnTp18qSkpEiHISJFhLszftFmHvhyPss276Fbw8rcd1oLWtcuH+nQJA+Z2Qx3D3sAoUzHVef38RKmu/umvI2sYKjuFBH5XXq685ePZ/PhjDXcfmJThpzYJNIhFUqHqzvDHQW3K7DW3Z939+eANWZ2XF4HKSISbcyMPs2r8c1tPfnXWa1YtDGZM56bxN0f/sLGXfsjHZ5EkJldCEwn1Fp5ITDNzM6PbFQiInK0YoJnQs/vWIenvlvMM2OXRDqkqBJuF9wXgWMzvd+TRZmISLEVHxvDZd3qc2b72rwwbilv/LiSL+es54ZejbiuR0NKloiNdIhS8O4DOme0egbjJ3xHaO5NERGJYhlJqDs8OWYxBtzSTy2h4Qh3ECLLPECBu6eTQ/JqZq+b2SYzm5up7DEzW2hms83sUzOrEJTXN7N9ZjYrWF7KdExHM5tjZkvN7JmMybLNLMHM3g/KpwXzpYmIRFT5kvHcM6AFY+7oSe9mVXlyzGL6PD6eT35eQ3q6xnkpZmIO6XK7lfDrXRERKeRiY4xHz2/LucfW5okxi3lWLaFhCbciXG5mt5pZfLAMAZbncMxwoP8hZWOA1u7eFlgM3JNp2zJ3bx8sN2Qqf5HQPGRNgiXjnNcA2929MfAU8EiY1yIiku+OqVyaFy7tyIc3dKNauQTu+OAXzn7hR6av+J/ZLqTo+iYYAfdKM7sS+JLfBwISEZEiIDbGeOz8dpzbIZSEPve9ktCchJuA3gAcD6wlNH/YcQSTUx+Ou0/kkHnF3P1bd08N3k7lj3OU/Y9grrNy7j4laIF9Ezg72HwWMCJY/wjol9E6KiJSWHSuX4n/3tSdpwa2Y3PyAS58eQo3vj2DVVv3RDo0yWfufjfwMtAWaAe84u5/jmxUIiKS12JjjMcuCCWhj3+7mOfHLY10SIVaWM+ABl2I8nqY9quB9zO9b2BmM4FdwF/d/QegNqGEN8OaoIzgdXUQX6qZ7QQqA1sO/SAzG0SQMNerVy+PL0NEJHsxMcY5HerQv1VNXvthOS9OWMZ3CzZy5fH1ublvE8qXjI90iJKHzKwxUN3df3T3T4BPgvKeZtbI3ZdFNkIREclrGUmoA4+NXgTA4D6NIxtUIRVWAmpmiYS6vLYCEjPK3f3qI/lQM7sPSAXeCYrWA/XcfauZdQT+a2atCE35cqiMh6iy2/bHQvdXgFcgNJT8kcQsInK0SpaI5ZZ+TRjYuS5PfLuY1yat4KMZa7j9pKZc3KUe8bF6PLCIeJrQ3J2H2htsO6MggxERkYIRG2M8fkE73J3HRi/CDG7qrST0UOH+tfMWUAM4BZhAqOts8pF8oJldAZwOXJoxsJG7H3D3rcH6DGAZ0JRQi2fmbrp1gHXB+hqgbnDOOKA8h3T5FREpjKqVS+SR89vyxS0n0LxGOf4+ah79n57I9ws3Eu7czFKo1Xf32YcWunsSUL/gwxERkYISG2M8cWF7zmpfi0e/WcSL49Xp5VDhJqCN3f1vwB53HwGcBrTJ7YeZWX/gL8CZ7r43U3lVM4sN1hsSGmxoubuvB5LNrGvwfOflwKjgsM+AK4L184HvXX+5iUgUaVWrPO9edxyvXd4Jd7h6eBKXDZvOwg27Ih2aHJ3EbLaVLLAoREQkImJjjCcuaMeZ7WrxyDcLlYQeItx5QFOC1x1m1hrYQA53cc3sPaA3UMXM1gD/IDTqbQIwJhgvaGow4m1P4J9mlgqkATe4e0Zr5o2ERtQtSWj0wIwRBIcBb5nZUkItn3n9jKqISL4zM05sWZ1ezary9tRVPP3dEgb85wcGdq7L7Sc1pVrZ7HIZKaR+MrPr3P3VzIVmdg0wI0IxiYhIAYqLjeHJC0PPhD7yzULM4IZejSIdVqFg4TQamtm1wMeEWj2HA2WAv7n7y/kaXT7o1KmTJyUlRToMEZEs7dh7kGe/X8qbU1ZSIjaGm/o05poTGpAYHxvp0Io1M5vh7p3C3Lc68ClwkN8Tzk5ACeAcd9+QP1HmH9WdIiJHJjUtnds/+IXPf1nHPac25/pilIQeru7MtgXUzIa4+3+ABe6+HZgINMynGEVEir0KpUrwt9Nb8qeux/Dw1wt4bPQi3p32K3/u34wz29VCs00Vfu6+ETjezPoArYPiL939+wiGJSIiERAXG8NTF4YGJnro61BL6KCexScJzUpOz4BeFbw+m9+BiIjI7xpUKc3Ll3Xiveu6UqFUPENGzuKcFyYzY5XGWosW7j7O3Z8NFiWfIiLFVFxsDE8PbM9pbWvy768W8urE5ZEOKaJyegZ0gZmtBKqZWeYR/Qxwd2+bb5GJiAjdGlXm85tP4JOZa3ls9ELOe3EKp7WtydD+zalbqVSkwxMREZEwxMXG8J+B7cHhwa8WAHBdz+LZsTTbBNTdLzazGsBo4MyCCUlERDKLiTHO71iHAW1q8PKE5bw8cRlj5m/k6u4NuKlPI8olxkc6RBEREclBXGwMT1/UHggloYnxMVzWrX5EY4qEcEbB3QzMcfdV+R2MiIgcXqkScdx+UlMu7lKPx0Yv4qUJy/gwaTW3n9SUizrXJS423Jm1REREJBLigyR0f0oafxs1jzKJcZzToU6kwypQOf614u5phKZSKVEA8YiISA5qlE/kiQvb8fnNJ9CoWhn++t+5DHjmByYs3hzp0ERERCQH8bExPH/psXRrWJm7PpzN6HlRNzj6UQn3dvkq4Ecz+5uZ3ZGx5GdgIiKSvTZ1yvP+oK689KeOHEhN54rXp3PF69NZvDE50qGJiIhINhLjY3n1ik60rl2eW96dyY9Lt0Q6pAITbgK6Dvgi2L9spkVERCLIzOjfugZjbu/FX09rwcxft9P/6Ync9+kctuw+EOnwRERE5DDKJMQx4qrONKhSmuveTGLGqu2RDqlAmLtHOoYCpcm0RaQo277nIP8Zu4S3pq6iVHwsg/s25srj65MYHxvp0KLa4SbTLi5Ud4qI5J9Nu/ZzwctT2L7nICMHdaNlrXKRDilPHK7uDKsF1MzGmdn3hy55H6aIiByNiqVLcP+ZrRh9W0+Oa1iJh79eyElPTeDL2espbjccRUREokG1com8fc1xlE6I4/LXp7F88+5Ih5Svwu2Cexdwd7D8DZgF6FaoiEgh1bhaGV67ojPvXHscpUvEMfjdnzn/pSnMWr0j0qGJiIjIIepWKsVb1xyHO/zptWms3bEv0iHlm7ASUHefkWn50d3vAI7L59hEROQodW9chS9v7cEj57Vh1da9nP38jwwZObNIV2zRzMzqBr2OFpjZPDMbEpRXMrMxZrYkeK0YlJ9kZjPMbE7w2jfTucab2SIzmxUs1YLyBDN738yWmtk0M6sfkYsVEZE/aFytDCOu7kLygVT+9No0NicXzbEcwu2CWynTUsXMTgFq5HNsIiKSB2JjjIGd6zH+7t7c0rcx38zdQN/Hx/PY6IXsPpAa6fDkj1KBO929BdAVGGxmLYGhwFh3bwKMDd4DbAHOcPc2wBXAW4ec71J3bx8sm4Kya4Dt7t4YeAp4JH8vSUREwtW6dnneuLIzG3bu57Jh09i5NyXSIeW5cLvgziDU5XYGMAW4k1AFJiIiUaJMQhx3ntyMcXf15tTWNXh+3DJ6PzaekdN/JS1dz4cWBu6+3t1/DtaTgQVAbeAsYESw2wjg7GCfme6+LiifBySaWUIOH5P5XB8B/czM8uwiRETkqHSqX4mXL+vI8s17uHL4dPYUsZvF4XbBbeDuDYPXJu5+srtPyu/gREQk79WqUJKnL+rAfwd3p37lUgz9ZA6nPfMDk5YUnznIokHQNbYDMA2o7u7rIZSkAtWyOOQ8YKa7Z+6z9UbQ/fZvmZLM2sDq4FypwE6gchafP8jMkswsafPmzXl1WSIiEoaeTavyzMXt+WX1Dga9lcT+lLRIh5Rnsk1AzayzmdXI9P5yMxtlZs+YWaX8D09ERPJL+7oV+PCGbrxw6bHsOZjKn4ZN45rhP7F0U9EefS8amFkZ4GPgNnffFcb+rQh1pb0+U/GlQdfcHsFyWcbuWZzif5rA3f0Vd+/k7p2qVq2a20sQEZGj1L91TR49vx0/Lt3KLe/NJCUtPdIh5YmcWkBfBg4CmFlP4GHgTUJ3S1/J39BERCS/mRkD2tRkzO29uOfU5kxfsY1Tnp7IP0bNZfueg5EOr1gys3hCyec77v5JULzRzGoG22sCmzLtXwf4FLjc3ZdllLv72uA1GXgX6BJsWgPUDY6NA8oD2/LzmkRE5Mic37EO/3dmK8bM38ifP5pNehF4ZCanBDTW3TMqpYHAK+7+sbv/DWicv6GJiEhBSYyP5fpejRh/d28u6VKPt6f9Sq/HxvHaD8s5kFp0uv0UdkE32WHAAnd/MtOmzwgNMkTwOirYvwLwJXCPu/+Y6TxxZlYlWI8HTgfmZnGu84HvXZPEiogUWlccX5+7Tm7KpzPX8vfP5kb9vN45JqDB3VGAfsD3mbbFZbG/iIhEscplEvjX2a35ZkgPjj2mIg98uYCTn5rIN3PXR32FFyW6E+oq2zfT9CkDCPVAOsnMlgAnBe8BbiZ0Q/hvh0y3kgCMNrPZhObuXgu8GhwzDKhsZkuBO/h9RF0RESmkBvdpzPU9G/L21F95dPSiSIdzVHJKIt8DJpjZFmAf8AOAmTUm1A1XRESKoCbVyzL8qi5MWLyZB7+czw1v/0yXBpX422ktaVOnfKTDK7KCAf4ONyJtvyz2fwB44DD7dzzMZ+wHLjiiAEVEJCLMjKGnNif5QCovjl9G2cQ4buodnR1Ss01A3f1BMxsL1AS+zdRFJwa4Jb+DExGRyOrVtCrdG/Xgg6Q1PDlmEWc8N4lzj63N3ac0o2b5kpEOT0REpNgwM/51Vmv2HEjl0W8WUTYxnsu6HhPpsHItx2607j41i7LF+ROOiIgUNnGxMVxyXD3OaFeTF8YvY9ikFXw1Zz2Dejbihl4NKVVCT2SIiIgUhNgY4/EL2rHnQCp/HzWXMgmxnNOhTqTDypWw5gEVEREpmxjPX/o3Z+wdvTipZQ2eGbuE3o+N58Ok1UViVD4REZFoEB8bw3OXHEvXBpW568PZfDtvQ6RDypV8S0DN7HUz22RmczOVVTKzMWa2JHitmGnbPWa21MwWmdkpmco7mtmcYNszGRNpm1mCmb0flE8LJuwWEZF8VrdSKZ69uAMf33g8tSqU5O6PZnPGc5OYsmxrpEMTEREpFhLjY3n1ik60rl2em9+dydTl0VMH52cL6HCg/yFlQ4Gx7t4EGBu8x8xaAhcBrYJjXjCz2OCYF4FBQJNgyTjnNcB2d28MPEVoAm4RESkgHY+pyKc3Hc8zF3dgx94ULn51KoPeTGLFlj2RDk1ERKTIK5MQx4irOlOvcimuezOJRRuSIx1SWPItAXX3ifzvxNZnASOC9RHA2ZnKR7r7AXdfASwFugSTbZdz9ynBAEhvHnJMxrk+AvpltI6KiEjBMDPObFeLsXf24s/9mzF52VZOenIC//x8Pjv2Hox0eCIiIkVahVIlGH5VZ0rGx3LF69NZt2NfpEPKUUE/A1rd3dcDBK/VgvLawOpM+60JymoH64eW/+EYd08lNC1M5XyLXEREDisxPpabejdm3F29uaBTXYZPXkGvx8bz+qQVpKSlRzo8ERGRIqtOxVIMv6oLew6kcuUb09m5NyXSIWWrsAxClFXLpWdTnt0x/3tys0FmlmRmSZs3bz7CEEVEJCdVyybw0Llt+GpID9rWKc8/v5jPKU9NZMz8jfw+k5eIiIjkpZa1yvHyZR1ZsWUP172VxP6UtEiHdFgFnYBuDLrVErxuCsrXAHUz7VcHWBeU18mi/A/HmFkcUJ7/7fILgLu/4u6d3L1T1apV8+hSRETkcJrXKMebV3fhjSs7YwbXvZnEJa9OY966nZEOTUREpEg6vnEVHr+gHdNXbOOOD2YV2hHqCzoB/Qy4Ili/AhiVqfyiYGTbBoQGG5oedNNNNrOuwfOdlx9yTMa5zge+d91eFxEpNMyMPs2r8c1tPfnXWa1YuGEXpz87iT9/9Aubdu2PdHgiIiJFzlnta3PfgBZ8NWcD//xifqHsfZRvs4eb2XtAb6CKma0B/gE8DHxgZtcAvwIXALj7PDP7AJgPpAKD3T2j3fhGQiPqlgS+DhaAYcBbZraUUMvnRfl1LSIicuTiY2O4rFt9zmxfmxfGLeWNH1fyxez13NCrEdf1aEjJErE5n0RERETCcm2PBqzfuZ/Xf1xBzfKJXN+rUaRD+gMrjFlxfurUqZMnJSVFOgwRkWJr1dY9PPLNQr6as4Ea5RL5c/9mnN2+NjExhXcgczOb4e6dIh1HpKjuFBGJLunpzi0jZ/Ll7PU8PbA9Z3eonfNBeexwdWdhGYRIRESKiWMql+aFSzvy4Q3dqFYugTs++IWzX/iR6SuyfIxfREREcikmxnjywnZ0bViJuz/6hUlLtkQ6pN8oARURkYjoXL8S/72pO08NbMfm5ANc+PIUbnx7Bqu27ol0aCIiIlEvIS6Wly/rRMMqZbjh7RmFZiBAJaAiIhIxMTHGOR3q8P2dvbnzpKZMWLyZk56cyL+/WsDOfYV7HjMREZHCrnzJeIZf3ZmyiXFc+cZPrN62N9IhKQEVEZHIK1killv6NWH8Xb05p0NtXv1hOb0fG8ebU1aSkpYe6fBERESiVs3yJRlxdRcOpKRxxRvT2b7nYETjUQIqIiKFRrVyiTxyflu+uOUEmtcox99HzaP/0xP5fuHGQjmUvIiISDRoWr0sr13RmTXb93HNiJ/Yn5KW80H5RAmoiIgUOq1qlefd647jtcs74Q5XD0/i8tens3DDrkiHJiIiEpW6NKjE0wPbM3P1Dm55byZp6ZG5sasEVERECiUz48SW1fnmtp7844yWzF6zkwH/+YF7PpnNpuT9kQ5PREQk6gxoU5N/nN6SMfM38o/P5kakd1FcgX+iiIhILpSIi+Gq7g04p0Ntnv1+KSMmr+SzWeu4qU9jrjmhAYnxsZEOUUREJGpc2b0B63ft5+UJy6lRLpGb+zYp0M9XC6iIiESFCqVK8LfTWzLmjl6c0KQKj41eRL8nJjBq1lo9HyoiIpILfzmlOed0qM3j3y7mw6TVBfrZSkBFRCSqNKhSmpcv68R713WlQql4hoycxbkvTmbGqu2RDu2omVldMxtnZgvMbJ6ZDQnKK5nZGDNbErxWzHTMPWa21MwWmdkpmco7mtmcYNszZmZBeYKZvR+UTzOz+gV+oSIiElExMcYj57XlhMZVGPrJHMYt2lRwn11gnyQiIpKHujWqzOc3n8Bj57dl7fZ9nPfiZG5+9+dCMcfZUUgF7nT3FkBXYLCZtQSGAmPdvQkwNnhPsO0ioBXQH3jBzDL6JL8IDAKaBEv/oPwaYLu7NwaeAh4piAsTEZHCpURcDC/+6ViaVS/L4Hd+ZvaaHQXyuUpARUQkasXEGBd0qsv4u3szpF8TvluwkX5PTuDhrxeSvD8l0uHlmruvd/efg/VkYAFQGzgLGBHsNgI4O1g/Cxjp7gfcfQWwFOhiZjWBcu4+xUP9k9885JiMc30E9MtoHRURkeKlbGI8w6/qTKXSJbh6+E+s2ron3z9TCaiIiES9UiXiuP2kpoy/qw9ntK3FSxOW0fux8QX+XEteCrrGdgCmAdXdfT2EklSgWrBbbSDzRa4JymoH64eW/+EYd08FdgKVs/j8QWaWZGZJmzdvzqOrEhGRwqZauURGXN2F1HTn8tens2X3gXz9PCWgIiJSZNQon8gTF7bj85tPoFG1MqzZvi/SIR0RMysDfAzc5u7ZTX6aVculZ1Oe3TF/LHB/xd07uXunqlWr5hSyiIhEsUZVyzDsis6kpKazfkf+TnWmaVhERKTIaVOnPO8P6kpqhCbZPhpmFk8o+XzH3T8JijeaWU13Xx90r80YLWINUDfT4XWAdUF5nSzKMx+zxszigPLAtny5GBERiRodj6nIuLt7kxCXv9ObqQVURESKJDMjPja6qrngWcxhwAJ3fzLTps+AK4L1K4BRmcovCka2bUBosKHpQTfdZDPrGpzz8kOOyTjX+cD3rnlsREQE8j35BLWAioiIFCbdgcuAOWY2Kyi7F3gY+MDMrgF+BS4AcPd5ZvYBMJ/QCLqD3T0tOO5GYDhQEvg6WCCU4L5lZksJtXxelM/XJCIi8hsloCIiIoWEu08i62c0Afod5pgHgQezKE8CWmdRvp8ggRURESlo0dU3SURERERERKKWElAREREREREpEEpARUREREREpEAoARUREREREZECoQRURERERERECoQSUBERERERESkQVtzmnjazzcCqPDhVFWBLHpynMCgq16LrKHyKyrXoOgqfgr6WY9y9agF+XqESRt1ZlP5tFRT9zHJHP6/c0c8r9/Qzy51wfl5Z1p3FLgHNK2aW5O6dIh1HXigq16LrKHyKyrXoOgqfonQtRYG+j9zTzyx39PPKHf28ck8/s9w5mp+XuuCKiIiIiIhIgVACKiIiIiIiIgVCCeiReyXSAeShonItuo7Cp6hci66j8ClK11IU6PvIPf3Mckc/r9zRzyv39DPLnSP+eekZUBERERERESkQagEVERERERGRAqEEVERERERERAqEEtAjYGb9zWyRmS01s6GRjicnZrbSzOaY2SwzSwrKKpnZGDNbErxWzLT/PcG1LTKzUyIY9+tmtsnM5mYqy3XcZtYxuP6lZvaMmVkhuZb7zWxt8L3MMrMBhf1azKyumY0zswVmNs/MhgTlUfW9ZHMdUfWdmFmimU03s1+C6/i/oDyqvo8criWqvpPiyKKsToykw/3ukeyZWayZzTSzLyIdSzQwswpm9pGZLQz+rXWLdEyFmZndHvx/nGtm75lZYqRjKmwO83fsYf/WyJG7a8nFAsQCy4CGQAngF6BlpOPKIeaVQJVDyh4FhgbrQ4FHgvWWwTUlAA2Ca42NUNw9gWOBuUcTNzAd6AYY8DVwaiG5lvuBu7LYt9BeC1ATODZYLwssDuKNqu8lm+uIqu8k+MwywXo8MA3oGm3fRw7XElXfSXFbiMI6McI/ryx/90Q6rsK+AHcA7wJfRDqWaFiAEcC1wXoJoEKkYyqsC1AbWAGUDN5/AFwZ6bgK20Iu/iYPZ1ELaO51AZa6+3J3PwiMBM6KcExH4ixCv6AIXs/OVD7S3Q+4+wpgKaFrLnDuPhHYdkhxruI2s5pAOXef4qH/IW9mOqbAHOZaDqfQXou7r3f3n4P1ZGABoV/eUfW9ZHMdh1NYr8PdfXfwNj5YnCj7PiDbazmcQnstxUxRqRMLxBH87in2zKwOcBrwWqRjiQZmVo5QsjAMwN0PuvuOiAZV+MUBJc0sDigFrItwPIVOLv8mz5ES0NyrDazO9H4Nhb/ycOBbM5thZoOCsuruvh5CFSJQLSgv7NeX27hrB+uHlhcWN5vZ7KBrQ0bXhai4FjOrD3Qg1FIVtd/LIdcBUfadBF3TZgGbgDHuHrXfx2GuBaLsOylmCnudUWhl8btHsvY08GcgPcJxRIuGwGbgjaDb8mtmVjrSQRVW7r4WeBz4FVgP7HT3byMbVdQ43N8aOVICmntZPUtU2Oey6e7uxwKnAoPNrGc2+0bj9cHh4y7M1/Mi0AhoT+iX3hNBeaG/FjMrA3wM3Obuu7LbNYuyQnMtWVxH1H0n7p7m7u2BOoRaAFtns3uhvQ447LVE3XdSzOjnfQRy8Tu0WDOz04FN7j4j0rFEkThCXSVfdPcOwB5C3SMlC8FNzbMIPcpRCyhtZn+KbFRFnxLQ3FsD1M30vg6FvKne3dcFr5uATwl1mdoYdFUjeN0U7F7Yry+3ca8J1g8tjzh33xj8wZ0OvMrvXZ0L9bWYWTyhP5zecfdPguKo+16yuo5o/U4Agi5W44H+ROH3kVnma4nm76SYKOx1RqFzmN+hkrXuwJlmtpJQ9+6+ZvZ2ZEMq9NYAazL1IPmIUEIqWTsRWOHum909BfgEOD7CMUWLw/2tkSMloLn3E9DEzBqYWQngIuCzCMd0WGZW2szKZqwDJwNzCcV8RbDbFcCoYP0z4CIzSzCzBkATQgN6FBa5ijvoEpBsZl2DkTAvz3RMRGX8pw2cQ+h7gUJ8LcHnDgMWuPuTmTZF1fdyuOuItu/EzKqaWYVgvSShinQhUfZ9BPFneS3R9p0UQ1FVJ0ZaNr9DJQvufo+713H3+oT+bX3v7mqdyoa7bwBWm1mzoKgfMD+CIRV2vwJdzaxU8P+zH6FnsyVnh/tbI2fhjlak5Q8jQQ0gNHLdMuC+SMeTQ6wNCY1K+AswLyNeoDIwFlgSvFbKdMx9wbUtIoKjRwLvEepyl0Lojt41RxI30InQH63LgOcAKyTX8hYwB5gd/CeuWdivBTiBUPe62cCsYBkQbd9LNtcRVd8J0BaYGcQ7F/h7UB5V30cO1xJV30lxXIiiOjHSy+F+90Q6rmhYgN5oFNxwf1btgaTg39l/gYqRjqkwL8D/Ebp5OzeocxIiHVNhW8jl3+Q5LRacVERERERERCRfqQuuiIiIiIiIFAgloCIiIiIiIlIglICKiIiIiIhIgVACKiIiIiIiIgVCCaiIiIiIiIgUCCWgIiIiIiIiUiCUgIpEKTOrbGazgmWDma0N1neb2Qv58HnDzWyFmd2QzT49zGy+mc3N688XEZHCz8zSgrporpl9aGalcnFsezMbkOn9mWY2NIdjdh9NvGHGVT+jXjOzTmb2TB6cc7yZLTKzM3N53OSj/ezcMrNGGX9fFPRnS9GkeUBFigAzux/Y7e6P5+NnDCc0CfhHOexXP9ivdX7FIiIihZOZ7Xb3MsH6O8AMd38yjOPigD8Bndz95iP5vCOINdbd08LYrz55XK+Z2XjgLndPyqtzZvEZce6emofnO+KftUhmagEVKWLMrLeZfRGs329mI8zsWzNbaWbnmtmjZjbHzL4xs/hgv45mNsHMZpjZaDOrGcbnXBDc4f7FzCbm93WJiEjU+QFobGZnmNk0M5tpZt+ZWXX4rY56xcy+Bd4E/gkMDFrbBprZlWb2XLBvdTP7NKhzfjGz4w/9MDO728x+MrPZZvZ/WQUU9BL6p5lNA7qZ2d+DY+YGsViwX8fgc6YAgzMdf2gde1embXOD1tLSZvZlcPxcMxuY0w8qaBF9yswmmtkCM+tsZp+Y2RIzeyBz/JnW/xzU57+Y2cOZzvNvM5sADDGzfsHPfY6ZvW5mCcF+K83s/8zs52Bb86C8l/3eu2qmmZXNKXaR3FICKlL0NQJOA84C3gbGuXsbYB9wWpCEPguc7+4dgdeBB8M479+BU9y9HZCrLkQiIlK0BS2apwJzgElAV3fvAIwE/pxp147AWe5+CaF65X13b+/u7x9yymeACUGdcyww75DPOxloAnQB2gMdzaxnFqGVBua6+3HuPgl4zt07B62bJYHTg/3eAG51925HcPn9gXXu3i447zdhHnfQ3XsCLwGjCCW+rYErzaxy5h3N7FTgbOC44GfyaKbNFdy9F/A8MBwYGNT7ccCNmfbb4u7HAi8CGYn0XcBgd28P9CD0t4JInlICKlL0fe3uKYT+CIjl94pwDlAfaEaoghtjZrOAvwJ1wjjvj8BwM7suOK+IiEjJoC5JAn4FhhGqU0ab2RzgbqBVpv0/c/dwkpy+hBIl3D3N3Xcesv3kYJkJ/Aw0J5SQHioN+DjT+z5B6+yc4DNamVl5QknchGCft8KIL7M5wIlm9oiZ9cgi1sP5LNPx89x9vbsfAJYDdQ/Z90TgDXffC+Du2zJty0jemwEr3H1x8H4EkDkp/yR4nUHo7wEI1e1PmtmthH4GedaFVyRDXKQDEJF8dwDA3dPNLMV/f/A7ndDvACNU0eXqLq+732BmxxFqXZ1lZu3dfWteBi4iIlFnX9B69hszexZ40t0/M7PewP2ZNu/Jo8814CF3fzmH/fZnPPdpZonAC4SeO11tofEUEoNzhTNISip/bMxJBHD3xWbWERgAPGRm37r7P8M434HgNT3Tesb7Q/9mzy7GPZn2Cefz0jLO7+4Pm9mXQexTzexEd18YRuwiYVMLqIgsAqqaWTcAM4s3s1Y5HIOZNXL3ae7+d2AL/3t3VkREBKA8sDZYvyKb/ZKBwz1zOJag+6iZxZpZuUO2jwauNrOMAZBqm1m1HOJKDF63BMedD+DuO4CdZnZCsP3Swxy/klB3YMzsWKBBsF4L2OvubwOPZ+yTx74ldL2lgs+slMU+C4H6ZtY4eH8ZMCGL/X4T1O1z3P0RQq3YzfMwZhFACahIsefuBwlVuo+Y2S/ALOB/BnfIwmPBwAVzgYnAL/kXpYiIRLH7gQ/N7AdCNywPZxzQMmMQokO2DSHUXXYOoS6jf7hR6u7fAu8CU4J9PuLwyWzGMTuAVwl1ef0v8FOmzVcBzweDEB2ui/DHQKWgy/GNQEZX1zbA9KD8PuCBLI8+Cu7+DaEuu0nB59yVxT77CV3Hh8HPJJ3Q86XZuS0YOOkXQtf9dZ4GLoKmYRGRMJmmYREREckTVgDTsOQ10zQskkfUAioi4doJ/MvMbjjcDmbWA/ic7O9wi4iIFHfbCA3kV+hHkTezRkEr68ZIxyJFg1pARUREREREpECoBVREREREREQKhBJQERERERERKRBKQEVERERERKRAKAEVERERERGRAvH/yO2f1MOhcwUAAAAASUVORK5CYII=", - "text/plain": [ - "
" - ] - }, - "metadata": { - "needs_background": "light" - }, - "output_type": "display_data" - } - ], - "source": [ - "# solve\n", - "solver = pybamm.ScipySolver()\n", - "t = np.linspace(0, 3600, 600)\n", - "solution = solver.solve(model, t)\n", - "\n", - "# post-process, so that the solution can be called at any time t or space r\n", - "# (using interpolation)\n", - "c = solution[\"Concentration [mol.m-3]\"]\n", - "c_surf = solution[\"Surface concentration [mol.m-3]\"]\n", - "\n", - "# plot\n", - "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4))\n", - "\n", - "ax1.plot(solution.t, c_surf(solution.t))\n", - "ax1.set_xlabel(\"Time [s]\")\n", - "ax1.set_ylabel(\"Surface concentration [mol.m-3]\")\n", - "\n", - "r = mesh[\"negative particle\"].nodes # radial position\n", - "time = 1000 # time in seconds\n", - "ax2.plot(r * 1e6, c(t=time, r=r), label=f\"t={time}[s]\")\n", - "ax2.set_xlabel(\"Particle radius [microns]\")\n", - "ax2.set_ylabel(\"Concentration [mol.m-3]\")\n", - "ax2.legend()\n", - "\n", - "plt.tight_layout()\n", - "plt.show()" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "In the [next notebook](./4-comparing-full-and-reduced-order-models.ipynb) we consider the limit of fast diffusion in the particle. This leads to a reduced-order model for the particle behaviour, which we compare with the full (Fickian diffusion) model. " - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## References\n", - "\n", - "The relevant papers for this notebook are:" - ] - }, - { - "cell_type": "code", - "execution_count": 10, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "[1] Joel A. E. Andersson, Joris Gillis, Greg Horn, James B. Rawlings, and Moritz Diehl. CasADi – A software framework for nonlinear optimization and optimal control. Mathematical Programming Computation, 11(1):1–36, 2019. doi:10.1007/s12532-018-0139-4.\n", - "[2] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, and others. Array programming with NumPy. Nature, 585(7825):357–362, 2020. doi:10.1038/s41586-020-2649-2.\n", - "[3] Valentin Sulzer, Scott G. Marquis, Robert Timms, Martin Robinson, and S. Jon Chapman. Python Battery Mathematical Modelling (PyBaMM). Journal of Open Research Software, 9(1):14, 2021. doi:10.5334/jors.309.\n", - "[4] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, and others. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature Methods, 17(3):261–272, 2020. doi:10.1038/s41592-019-0686-2.\n", - "\n" - ] - } - ], - "source": [ - "pybamm.print_citations()" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "pybamm", - "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.16" - }, - "toc": { - "base_numbering": 1, - "nav_menu": {}, - "number_sections": true, - "sideBar": true, - "skip_h1_title": false, - "title_cell": "Table of Contents", - "title_sidebar": "Contents", - "toc_cell": false, - "toc_position": {}, - "toc_section_display": true, - "toc_window_display": true - }, - "vscode": { - "interpreter": { - "hash": "187972e187ab8dfbecfab9e8e194ae6d08262b2d51a54fa40644e3ddb6b5f74c" - } - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/docs/source/examples/notebooks/creating_models/old_notebooks/4-comparing-full-and-reduced-order-models.ipynb b/docs/source/examples/notebooks/creating_models/old_notebooks/4-comparing-full-and-reduced-order-models.ipynb deleted file mode 100644 index f180d16f0d..0000000000 --- a/docs/source/examples/notebooks/creating_models/old_notebooks/4-comparing-full-and-reduced-order-models.ipynb +++ /dev/null @@ -1,474 +0,0 @@ -{ - "cells": [ - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Comparing full and reduced-order models" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "In the [previous notebook](./3-negative-particle-problem.ipynb) we saw how to solve the problem of diffusion on a sphere, motivated by the problem in the negative particle in battery modelling. In this notebook we consider a reduced-order ODE model for the particle behaviour, suitable in the limit of fast diffusion. We also show how to compare the results of the full and reduced-order models. \n", - "\n", - "In the limit of fast diffusion in the particles the concentration is uniform in $r$. This result in the following ODE model for the (uniform) concentration in the particle \n", - "\n", - "$$\n", - " \\frac{\\textrm{d} c}{\\textrm{d} t} = -3\\frac{j}{RF}\n", - "$$\n", - "with the initial condition:\n", - "$$\n", - "\\left.c\\right\\vert_{t=0} = c_0,\n", - "$$\n", - "where $c$$ is the concentration, $r$ the radial coordinate, $t$ time, $R$ the particle radius, $D$ the diffusion coefficient, $j$ the interfacial current density, $F$ Faraday's constant, and $c_0$ the initial concentration. \n", - "\n", - "As in the previous example we use the following parameters:\n", - "\n", - "| Symbol | Units | Value |\n", - "|:-------|:-------------------|:-----------------------------------------------|\n", - "| $R$ | m | $10 \\times 10^{-6}$ |\n", - "| $D$ | m${^2}$ s$^{-1}$ | $3.9 \\times 10^{-14}$ |\n", - "| $j$ | A m$^{-2}$ | $1.4$ |\n", - "| $F$ | C mol$^{-1}$ | $96485$ |\n", - "| $c_0$ | mol m$^{-3}$ | $2.5 \\times 10^{4}$ |" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Setting up the models\n", - "As in the single particle diffusion example, we begin by importing the pybamm library into this notebook, along with any other packages we require. In this notebook we want to compare the results of the full and reduced-order models, so we create two empty `pybamm.BaseModel` objects. We can pass in a name when we create the model, for easy reference. " - ] - }, - { - "cell_type": "code", - "execution_count": 11, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "\u001b[33mWARNING: You are using pip version 22.0.4; however, version 22.3.1 is available.\n", - "You should consider upgrading via the '/home/mrobins/git/PyBaMM/env/bin/python -m pip install --upgrade pip' command.\u001b[0m\u001b[33m\n", - "\u001b[0mNote: you may need to restart the kernel to use updated packages.\n" - ] - } - ], - "source": [ - "%pip install \"pybamm[plot,cite]\" -q # install PyBaMM if it is not installed\n", - "import pybamm\n", - "import numpy as np\n", - "import matplotlib.pyplot as plt\n", - "\n", - "full_model = pybamm.BaseModel(name=\"full model\")\n", - "reduced_model = pybamm.BaseModel(name=\"reduced model\")" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "It can be useful to add the models to a list so that we can perform the same operations on each model easily" - ] - }, - { - "cell_type": "code", - "execution_count": 12, - "metadata": {}, - "outputs": [], - "source": [ - "models = [full_model, reduced_model]" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We then define our parameters, as seen previously, " - ] - }, - { - "cell_type": "code", - "execution_count": 13, - "metadata": {}, - "outputs": [], - "source": [ - "R = pybamm.Parameter(\"Particle radius [m]\")\n", - "D = pybamm.Parameter(\"Diffusion coefficient [m2.s-1]\")\n", - "j = pybamm.Parameter(\"Interfacial current density [A.m-2]\")\n", - "F = pybamm.Parameter(\"Faraday constant [C.mol-1]\")\n", - "c0 = pybamm.Parameter(\"Initial concentration [mol.m-3]\")" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The reduced order model solves and ODE for the (uniform) concentration in the particle. In the parameter regime where this is valid, we expect that the solution of the ODE model should agree with the average concentration in the PDE mode. In anticipation of this, we create two variables: the concentration (which we will use in the PDE model), and the average concentration (which we will use in the ODE model). This will make it straightforward to compare the results in a consistent way. Note that the average concentration doesn't have a domain since it is a scalar quantity." - ] - }, - { - "cell_type": "code", - "execution_count": 14, - "metadata": {}, - "outputs": [], - "source": [ - "c = pybamm.Variable(\"Concentration [mol.m-3]\", domain=\"negative particle\")\n", - "c_av = pybamm.Variable(\"Average concentration [mol.m-3]\")" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now we define our model equations, initial and boundary conditions (where appropriate)" - ] - }, - { - "cell_type": "code", - "execution_count": 15, - "metadata": {}, - "outputs": [], - "source": [ - "# governing equations for full model\n", - "N = -D * pybamm.grad(c) # flux\n", - "dcdt = -pybamm.div(N)\n", - "full_model.rhs = {c: dcdt}\n", - "\n", - "# governing equations for reduced model\n", - "dc_avdt = -3 * j / R / F\n", - "reduced_model.rhs = {c_av: dc_avdt}\n", - "\n", - "# initial conditions (these are the same for both models)\n", - "full_model.initial_conditions = {c: c0}\n", - "reduced_model.initial_conditions = {c_av: c0}\n", - "\n", - "# boundary conditions (only required for full model)\n", - "lbc = pybamm.Scalar(0)\n", - "rbc = -j / F / D\n", - "full_model.boundary_conditions = {\n", - " c: {\"left\": (lbc, \"Neumann\"), \"right\": (rbc, \"Neumann\")}\n", - "}" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We can now populate the variables dictionary of both models with any variables of interest. We can compute the average concentration in the full model using the operator `pybamm.r_average`. We may also wish to compare the concentration profile predicted by the full model with the uniform concentration profile predicted by the reduced model. We can use the operator `pybamm.PrimaryBroadcast` to broadcast the scalar valued uniform concentration across the particle domain so that it can be visualised as a function of $r$. \n", - "\n", - "Note: the \"Primary\" refers to the fact the we are broadcasting in only one dimension. For some models, such as the DFN, variables may depend on a \"pseudo-dimension\" (e.g. the position in $x$ across the cell), but spatial operators only act in the \"primary dimension\" (e.g. the position in $r$ within the particle). If you are unfamiliar with battery models, do not worry, the details of this are not important for this example. For more information see the [broadcasts notebook](../expression_tree/broadcasts.ipynb)." - ] - }, - { - "cell_type": "code", - "execution_count": 16, - "metadata": {}, - "outputs": [], - "source": [ - "# full model\n", - "full_model.variables = {\n", - " \"Concentration [mol.m-3]\": c,\n", - " \"Surface concentration [mol.m-3]\": pybamm.surf(c),\n", - " \"Average concentration [mol.m-3]\": pybamm.r_average(c),\n", - "}\n", - "\n", - "# reduced model\n", - "reduced_model.variables = {\n", - " \"Concentration [mol.m-3]\": pybamm.PrimaryBroadcast(c_av, \"negative particle\"),\n", - " \"Surface concentration [mol.m-3]\": c_av, # in this model the surface concentration is just equal to the scalar average concentration\n", - " \"Average concentration [mol.m-3]\": c_av,\n", - "}" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Using the model" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "As before, we provide our parameter values" - ] - }, - { - "cell_type": "code", - "execution_count": 17, - "metadata": {}, - "outputs": [], - "source": [ - "param = pybamm.ParameterValues(\n", - " {\n", - " \"Particle radius [m]\": 10e-6,\n", - " \"Diffusion coefficient [m2.s-1]\": 3.9e-14,\n", - " \"Interfacial current density [A.m-2]\": 1.4,\n", - " \"Faraday constant [C.mol-1]\": 96485,\n", - " \"Initial concentration [mol.m-3]\": 2.5e4,\n", - " }\n", - ")" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We then define and process our geometry, and process both of the models" - ] - }, - { - "cell_type": "code", - "execution_count": 18, - "metadata": {}, - "outputs": [], - "source": [ - "# geometry\n", - "r = pybamm.SpatialVariable(\n", - " \"r\", domain=[\"negative particle\"], coord_sys=\"spherical polar\"\n", - ")\n", - "geometry = {\"negative particle\": {r: {\"min\": pybamm.Scalar(0), \"max\": R}}}\n", - "param.process_geometry(geometry)\n", - "\n", - "# models\n", - "for model in models:\n", - " param.process_model(model)" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We can now set up our mesh, choose a spatial method, and discretise our models. Note that, even though the reduced-order model is an ODE model, we discretise using the mesh for the particle so that our `PrimaryBroadcast` operator is discretised in the correct way." - ] - }, - { - "cell_type": "code", - "execution_count": 19, - "metadata": {}, - "outputs": [], - "source": [ - "# mesh\n", - "submesh_types = {\"negative particle\": pybamm.Uniform1DSubMesh}\n", - "var_pts = {r: 20}\n", - "mesh = pybamm.Mesh(geometry, submesh_types, var_pts)\n", - "\n", - "# discretisation\n", - "spatial_methods = {\"negative particle\": pybamm.FiniteVolume()}\n", - "disc = pybamm.Discretisation(mesh, spatial_methods)\n", - "\n", - "# process models\n", - "for model in models:\n", - " disc.process_model(model)" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Both models are now discretised and ready to be solved." - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Solving the model" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "As before, we choose a solver and times at which we want the solution returned. We then solve both models, post-process the results, and create a slider plot to compare the results." - ] - }, - { - "cell_type": "code", - "execution_count": 20, - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "2022-12-12 12:41:59.589 - [WARNING] processed_variable.get_spatial_scale(518): No length scale set for negative particle. Using default of 1 [m].\n", - "2022-12-12 12:41:59.609 - [WARNING] processed_variable.get_spatial_scale(518): No length scale set for negative particle. Using default of 1 [m].\n" - ] - }, - { - "data": { - "application/vnd.jupyter.widget-view+json": { - "model_id": "6fa9ebee57924ba5b79a2c51313fba25", - "version_major": 2, - "version_minor": 0 - }, - "text/plain": [ - "interactive(children=(FloatSlider(value=0.0, description='t', max=3600.0, step=1.0), Output()), _dom_classes=(…" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "# loop over models to solve\n", - "t = np.linspace(0, 3600, 600)\n", - "solutions = [None] * len(models) # create list to hold solutions\n", - "for i, model in enumerate(models):\n", - " solver = pybamm.ScipySolver()\n", - " solutions[i] = solver.solve(model, t)\n", - "\n", - "# post-process the solution of the full model\n", - "c_full = solutions[0][\"Concentration [mol.m-3]\"]\n", - "c_av_full = solutions[0][\"Average concentration [mol.m-3]\"]\n", - "\n", - "\n", - "# post-process the solution of the reduced model\n", - "c_reduced = solutions[1][\"Concentration [mol.m-3]\"]\n", - "c_av_reduced = solutions[1][\"Average concentration [mol.m-3]\"]\n", - "\n", - "# plot\n", - "r = mesh[\"negative particle\"].nodes # radial position\n", - "\n", - "\n", - "def plot(t):\n", - " fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4))\n", - "\n", - " # Plot concetration as a function of r\n", - " ax1.plot(r * 1e6, c_full(t=t, r=r), label=\"Full Model\")\n", - " ax1.plot(r * 1e6, c_reduced(t=t, r=r), label=\"Reduced Model\")\n", - " ax1.set_xlabel(\"Particle radius [microns]\")\n", - " ax1.set_ylabel(\"Concentration [mol.m-3]\")\n", - " ax1.legend()\n", - "\n", - " # Plot average concentration over time\n", - " t_hour = np.linspace(0, 3600, 600) # plot over full hour\n", - " c_min = c_av_reduced(t=3600) * 0.98 # minimum axes limit\n", - " c_max = param[\"Initial concentration [mol.m-3]\"] * 1.02 # maximum axes limit\n", - "\n", - " ax2.plot(t_hour, c_av_full(t=t_hour), label=\"Full Model\")\n", - " ax2.plot(t_hour, c_av_reduced(t=t_hour), label=\"Reduced Model\")\n", - " ax2.plot([t, t], [c_min, c_max], \"k--\") # plot line to track time\n", - " ax2.set_xlabel(\"Time [s]\")\n", - " ax2.set_ylabel(\"Average concentration [mol.m-3]\")\n", - " ax2.legend()\n", - "\n", - " plt.tight_layout()\n", - " plt.show()\n", - "\n", - "\n", - "import ipywidgets as widgets\n", - "\n", - "widgets.interact(plot, t=widgets.FloatSlider(min=0, max=3600, step=1, value=0));" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "From the results we observe that the reduced-order model does a good job of predicting the average concentration, but, since it is only an ODE model, cannot predicted the spatial distribution." - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "In the [next notebook](./5-half-cell-model.ipynb) we will show how to set up and solve a model which contains multiple domains." - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## References\n", - "\n", - "The relevant papers for this notebook are:" - ] - }, - { - "cell_type": "code", - "execution_count": 21, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "[1] Joel A. E. Andersson, Joris Gillis, Greg Horn, James B. Rawlings, and Moritz Diehl. CasADi – A software framework for nonlinear optimization and optimal control. Mathematical Programming Computation, 11(1):1–36, 2019. doi:10.1007/s12532-018-0139-4.\n", - "[2] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, and others. Array programming with NumPy. Nature, 585(7825):357–362, 2020. doi:10.1038/s41586-020-2649-2.\n", - "[3] Valentin Sulzer, Scott G. Marquis, Robert Timms, Martin Robinson, and S. Jon Chapman. Python Battery Mathematical Modelling (PyBaMM). Journal of Open Research Software, 9(1):14, 2021. doi:10.5334/jors.309.\n", - "[4] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, and others. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature Methods, 17(3):261–272, 2020. doi:10.1038/s41592-019-0686-2.\n", - "\n" - ] - } - ], - "source": [ - "pybamm.print_citations()" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "env", - "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.15" - }, - "toc": { - "base_numbering": 1, - "nav_menu": {}, - "number_sections": true, - "sideBar": true, - "skip_h1_title": false, - "title_cell": "Table of Contents", - "title_sidebar": "Contents", - "toc_cell": false, - "toc_position": {}, - "toc_section_display": true, - "toc_window_display": true - }, - "vscode": { - "interpreter": { - "hash": "19e5ebaa8d5a3277b4deed2928f02ad0cad6c3ab0b2beced644d557f155bce64" - } - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/docs/source/examples/notebooks/creating_models/old_notebooks/5-half-cell-model.ipynb b/docs/source/examples/notebooks/creating_models/old_notebooks/5-half-cell-model.ipynb deleted file mode 100644 index 685ac0b8d1..0000000000 --- a/docs/source/examples/notebooks/creating_models/old_notebooks/5-half-cell-model.ipynb +++ /dev/null @@ -1,693 +0,0 @@ -{ - "cells": [ - { - "attachments": {}, - "cell_type": "markdown", - "id": "professional-composer", - "metadata": {}, - "source": [ - "# A half cell model" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "id": "naval-management", - "metadata": {}, - "source": [ - "In the [previous notebook](./4-comparing-full-and-reduced-order-models.ipynb) we saw how to compare full and reduced-order models. Both of these models were posed on a single domain: the negative electrode particle. Here we will see how to create a model which contains multiple domains.\n", - "\n", - "We consider a problem posed on a half-cell geometry, which consists of a separator ($-L_s" - ] - }, - "execution_count": 35, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "# plot\n", - "pybamm.dynamic_plot(\n", - " solution,\n", - " [\n", - " \"Positive electrode potential [V]\",\n", - " \"Electrolyte potential [V]\",\n", - " \"Positive electrode interfacial current density [A.m-2]\",\n", - " \"Positive particle surface concentration [mol.m-3]\",\n", - " \"Average positive particle surface concentration [mol.m-3]\",\n", - " [\"Positive electrode OCP [V]\", \"Voltage [V]\"],\n", - " ],\n", - ")" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "id": "abandoned-shirt", - "metadata": {}, - "source": [ - "In the [next notebook](./6-a-simple-SEI-model.ipynb) we consider a simple model for SEI growth, and show how to correctly pose the model in non-dimensional form and then create and solve it using pybamm." - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "id": "independent-development", - "metadata": {}, - "source": [ - "## References\n", - "\n", - "The relevant papers for this notebook are:" - ] - }, - { - "cell_type": "code", - "execution_count": 36, - "id": "laden-replica", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "[1] Joel A. E. Andersson, Joris Gillis, Greg Horn, James B. Rawlings, and Moritz Diehl. CasADi – A software framework for nonlinear optimization and optimal control. Mathematical Programming Computation, 11(1):1–36, 2019. doi:10.1007/s12532-018-0139-4.\n", - "[2] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, and others. Array programming with NumPy. Nature, 585(7825):357–362, 2020. doi:10.1038/s41586-020-2649-2.\n", - "[3] Valentin Sulzer, Scott G. Marquis, Robert Timms, Martin Robinson, and S. Jon Chapman. Python Battery Mathematical Modelling (PyBaMM). Journal of Open Research Software, 9(1):14, 2021. doi:10.5334/jors.309.\n", - "\n" - ] - } - ], - "source": [ - "pybamm.print_citations()" - ] - } - ], - "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.0" - }, - "toc": { - "base_numbering": 1, - "nav_menu": {}, - "number_sections": true, - "sideBar": true, - "skip_h1_title": false, - "title_cell": "Table of Contents", - "title_sidebar": "Contents", - "toc_cell": false, - "toc_position": {}, - "toc_section_display": true, - "toc_window_display": true - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/docs/source/examples/notebooks/creating_models/old_notebooks/6-a-simple-SEI-model.ipynb b/docs/source/examples/notebooks/creating_models/old_notebooks/6-a-simple-SEI-model.ipynb deleted file mode 100644 index fbbe0cac5e..0000000000 --- a/docs/source/examples/notebooks/creating_models/old_notebooks/6-a-simple-SEI-model.ipynb +++ /dev/null @@ -1,681 +0,0 @@ -{ - "cells": [ - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Creating a Simple Model for SEI Growth\n", - "Before adding a new model, please read the [contribution guidelines](https://github.com/pybamm-team/PyBaMM/blob/develop/CONTRIBUTING.md)" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "In this notebook, we will run through the steps involved in creating a new model within pybamm. We will then solve and plot the outputs of the model. We have chosen to implement a very simple model of SEI growth. We first give a brief derivation of the model and discuss how to nondimensionalise the model so that we can show the full process of model conception to solution within a single notebook. " - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Note: if you run the entire notebook and then try to evaluate the earlier cells, you will likely receive an error. This is because the state of objects is mutated as it is passed through various stages of processing. In this case, we recommend that you restart the Kernel and then evaluate cells in turn through the notebook. " - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## A Simple Model of Solid Electrolyte Interphase (SEI) Growth" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The SEI is a porous layer that forms on the surfaces of negative electrode particles from the products of electrochemical reactions which consume lithium and electrolyte solvents. In the first few cycles of use, a lithium-ion battery loses a large amount of capacity; this is generally attributed to lithium being consumed to produce SEI. However, after a few cycles, the rate of capacity loss slows at a rate often (but not always) reported to scale with the square root of time. SEI growth is therefore often considered to be limited in some way by a diffusion process." - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Dimensional Model\n", - "\n", - "In our simple SEI model, we consider a one-dimensional SEI which extends from the surface of a planar negative electrode at $x=0$ until $x=L$, where $L$ is the thickness of the SEI. Since the SEI is porous, there is some electrolyte within the region $x\\in[0, L]$ and therefore some concentration of solvent, $c$. Within the porous SEI, the solvent is transported via a diffusion process according to:\n", - "$$\n", - "\\frac{\\partial c}{\\partial t} = - \\frac{\\partial N}{\\partial x}, \\quad N = - D(c) \\frac{\\partial c}{\\partial x} \\label{1}\\\\\n", - "$$\n", - "where $t$ is the time, $N$ is the solvent flux, and $D(c)$ is the effective solvent diffusivity (a function of the solvent concentration).\n", - "\n", - "On the electrode-SEI surface ($x=0$) the solvent is consumed by the SEI growth reaction, $R$. We assume that diffusion of solvent in the bulk electrolyte ($x>L$) is fast so that on the SEI-electrolyte surface ($x=L$) the concentration of solvent is fixed at the value $c_{\\infty}$. Therefore, the boundary conditions are\n", - "$$\n", - " N|_{x=0} = - R, \\quad c|_{x=L} = c_{\\infty},\n", - "$$\n", - "We also assume that the concentration of solvent within the SEI is initially uniform and equal to the bulk electrolyte solvent concentration, so that the initial condition is\n", - "$$\n", - " c|_{t=0} = c_{\\infty}\n", - "$$\n", - "\n", - "Since the SEI is growing, we require an additional equation for the SEI thickness. The thickness of the SEI grows at a rate proportional to the SEI growth reaction $R$, where the constant of proportionality is the partial molar volume of the reaction products, $\\hat{V}$. We also assume that the SEI is initially of thickness $L_0$. Therefore, we have\n", - "$$\n", - " \\frac{d L}{d t} = \\hat{V} R, \\quad L|_{t=0} = L_0\n", - "$$\n", - "\n", - "Finally, we assume for the sake of simplicity that the SEI growth reaction is irreversible and that the potential difference across the SEI is constant. The reaction is also assumed to be proportional to the concentration of solvent at the electrode-SEI surface ($x=0$). Therefore, the reaction flux is given by\n", - "$$\n", - " R = k c|_{x=0}\n", - "$$\n", - "where $k$ is the reaction rate constant (which is in general dependent upon the potential difference across the SEI).\n", - "\n", - "### Fixing the moving boundary\n", - "The model above is a moving boundary problem as it is defined in $x\\in[0, L]$. In order to solve it, we need to fix the boundary by introducing the scaling\n", - "$$\n", - " x = L \\xi.\n", - "$$\n", - "Then, applying the chain rule we have\n", - "$$\n", - " \\frac{\\partial}{\\partial x} \\rightarrow \\frac{1}{L} \\frac{\\partial}{\\partial \\xi}, \\quad \\text{ and } \\quad \\frac{\\partial}{\\partial t} \\rightarrow \\frac{\\partial}{\\partial t} - \\frac{L'(t)}{L(t)} \\xi \\frac{\\partial}{\\partial \\xi}.\n", - "$$\n", - "\n", - "Then, (1) can be rewritten as\n", - "$$\n", - " \\frac{\\partial c}{\\partial t} = \\frac{\\hat{V} R}{L} \\xi \\frac{\\partial c}{\\partial \\xi} + \\frac{1}{L^2} \\frac{\\partial}{\\partial \\xi} \\left( D(c) \\frac{\\partial c}{\\partial \\xi}\\right)\n", - "$$" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Entering the Model into PyBaMM" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "As always, we begin by importing PyBaMM and changing our working directory to the root of the `pybamm/` folder." - ] - }, - { - "cell_type": "code", - "execution_count": 1, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "\n", - "\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m A new release of pip is available: \u001b[0m\u001b[31;49m23.0.1\u001b[0m\u001b[39;49m -> \u001b[0m\u001b[32;49m23.2\u001b[0m\n", - "\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m To update, run: \u001b[0m\u001b[32;49mpip install --upgrade pip\u001b[0m\n", - "Note: you may need to restart the kernel to use updated packages.\n" - ] - } - ], - "source": [ - "%pip install \"pybamm[plot,cite]\" -q # install PyBaMM if it is not installed\n", - "import pybamm\n", - "import numpy as np\n", - "import os\n", - "\n", - "os.chdir(pybamm.__path__[0] + \"/..\")" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "A model is defined in six steps:\n", - "1. Initialise model\n", - "2. Define parameters and variables\n", - "3. State governing equations\n", - "4. State boundary conditions\n", - "5. State initial conditions\n", - "6. State output variables\n", - "\n", - "We shall proceed through each step to enter our simple SEI growth model." - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Initialise model" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We first initialise the model using the `BaseModel` class. This sets up the required structure for our model. " - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": {}, - "outputs": [], - "source": [ - "model = pybamm.BaseModel()" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Define parameters and variables" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "In our SEI model, we have two dimensionless parameters, $k$ and $\\hat{V}$, and one dimensionless function $D(c)$, which are all given in terms of the dimensional parameters, see (5). In pybamm, inputs are dimensional, so we first state all the dimensional parameters. We then define the dimensionless parameters, which are expressed an non-dimensional groupings of dimensional parameters. To define the dimensional parameters, we use the `Parameter` object to create parameter symbols. Parameters which are functions are defined using `FunctionParameter` object and should be defined within a python function as shown. " - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "metadata": {}, - "outputs": [], - "source": [ - "# dimensional parameters\n", - "k = pybamm.Parameter(\"Reaction rate constant [m.s-1]\")\n", - "L_0 = pybamm.Parameter(\"Initial thickness [m]\")\n", - "V_hat = pybamm.Parameter(\"Partial molar volume [m3.mol-1]\")\n", - "c_inf = pybamm.Parameter(\"Bulk electrolyte solvent concentration [mol.m-3]\")\n", - "\n", - "\n", - "def D(cc):\n", - " return pybamm.FunctionParameter(\n", - " \"Diffusivity [m2.s-1]\", {\"Solvent concentration [mol.m-3]\": cc}\n", - " )" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We now define the dimensionless variables in our model. Since these are the variables we solve for directly, we do not need to write them in terms of the dimensional variables. We simply use `SpatialVariable` and `Variable` to create the required symbols: " - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": {}, - "outputs": [], - "source": [ - "xi = pybamm.SpatialVariable(\"xi\", domain=\"SEI layer\", coord_sys=\"cartesian\")\n", - "c = pybamm.Variable(\"Solvent concentration [mol.m-3]\", domain=\"SEI layer\")\n", - "L = pybamm.Variable(\"SEI thickness [m]\")" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### State governing equations" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We can now use the symbols we have created for our parameters and variables to write out our governing equations. Note that before we use the reaction flux and solvent flux, we must derive new symbols for them from the defined parameter and variable symbols. Each governing equation must also be stated in the explicit form `d/dt = rhs` since pybamm only stores the right hand side (rhs) and assumes that the left hand side is the time derivative. The governing equations are then simply" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, - "outputs": [], - "source": [ - "# SEI reaction flux\n", - "R = k * pybamm.BoundaryValue(c, \"left\")\n", - "\n", - "# solvent concentration equation\n", - "N = -1 / L * D(c) * pybamm.grad(c)\n", - "dcdt = (V_hat * R) / L * pybamm.inner(xi, pybamm.grad(c)) - 1 / L * pybamm.div(N)\n", - "\n", - "# SEI thickness equation\n", - "dLdt = V_hat * R" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Once we have stated the equations, we can add them to the `model.rhs` dictionary. This is a dictionary whose keys are the variables being solved for, and whose values correspond right hand sides of the governing equations for each variable." - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": {}, - "outputs": [], - "source": [ - "model.rhs = {c: dcdt, L: dLdt}" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### State boundary conditions" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We only have boundary conditions on the solvent concentration equation. We must state where a condition is Neumann (on the gradient) or Dirichlet (on the variable itself). \n", - "\n", - "The boundary condition on the electrode-SEI (x=0) boundary is: \n", - "$$\n", - " N|_{\\xi=0} = - R, \\quad N|_{\\xi=0} = - D(c|_{\\xi=0} )\\left.\\frac{\\partial c}{\\partial \\xi}\\right|_{\\xi=0}\n", - "$$\n", - "which is a Neumann condition. To implement this boundary condition in pybamm, we must first rearrange the equation so that the gradient of the concentration, $\\nabla c|_{x=0}$, is the subject. Therefore we have\n", - "$$\n", - " \\left.\\frac{\\partial c}{\\partial \\xi}\\right|_{\\xi=0} = \\frac{R L}{D(c|_{\\xi=0} )}\n", - "$$\n", - "which we enter into pybamm as " - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": {}, - "outputs": [], - "source": [ - "D_left = pybamm.BoundaryValue(\n", - " D(c), \"left\"\n", - ") # pybamm requires BoundaryValue(D(c)) and not D(BoundaryValue(c))\n", - "grad_c_left = R * L / D_left" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "On the SEI-electrolyte boundary ($\\xi=1$), we have the boundary condition\n", - "$$\n", - " c|_{\\xi=1} = c_∞\n", - "$$" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "which is a Dirichlet condition and is just entered as" - ] - }, - { - "cell_type": "code", - "execution_count": 8, - "metadata": {}, - "outputs": [], - "source": [ - "c_right = c_inf" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We now load these boundary conditions into the `model.boundary_conditions` dictionary in the following way, being careful to state the type of boundary condition: " - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "metadata": {}, - "outputs": [], - "source": [ - "model.boundary_conditions = {\n", - " c: {\"left\": (grad_c_left, \"Neumann\"), \"right\": (c_right, \"Dirichlet\")}\n", - "}" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### State initial conditions" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "There are two initial conditions in our model:\n", - "$$\n", - " c|_{t=0} = c_∞, \\quad L|_{t=0} = L_0\n", - "$$" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "which are simply written in pybamm as" - ] - }, - { - "cell_type": "code", - "execution_count": 10, - "metadata": {}, - "outputs": [], - "source": [ - "c_init = c_inf\n", - "L_init = L_0" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "and then included into the `model.initial_conditions` dictionary:" - ] - }, - { - "cell_type": "code", - "execution_count": 11, - "metadata": {}, - "outputs": [], - "source": [ - "model.initial_conditions = {c: c_init, L: L_init}" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### State output variables" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We already have everything required in model for the model to be used and solved, but we have not yet stated what we actually want to output from the model. PyBaMM allows users to output any combination of symbols as an output variable therefore allowing the user the flexibility to output important quantities without further tedious postprocessing steps. \n", - "\n", - "Some useful outputs for this simple model are:\n", - "- the SEI thickness\n", - "- the SEI growth rate\n", - "- the solvent concentration\n", - "\n", - "These are added to the model by adding entries to the `model.variables` dictionary" - ] - }, - { - "cell_type": "code", - "execution_count": 12, - "metadata": {}, - "outputs": [], - "source": [ - "model.variables = {\n", - " \"SEI thickness [m]\": L,\n", - " \"SEI growth rate [m]\": dLdt,\n", - " \"Solvent concentration [mol.m-3]\": c,\n", - "}" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The model is now fully defined and ready to be used. If you plan on reusing the model several times, you can additionally set model defaults which may include: a default geometry to run the model on, a default set of parameter values, a default solver, etc." - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Using the Model" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The model will now behave in the same way as any of the inbuilt PyBaMM models. However, to demonstrate that the model works we display the steps involved in solving the model but we will not go into details within this notebook." - ] - }, - { - "cell_type": "code", - "execution_count": 13, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "" - ] - }, - "execution_count": 13, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "# define geometry\n", - "geometry = pybamm.Geometry(\n", - " {\"SEI layer\": {xi: {\"min\": pybamm.Scalar(0), \"max\": pybamm.Scalar(1)}}}\n", - ")\n", - "\n", - "\n", - "def Diffusivity(cc):\n", - " return cc * 10 ** (-12)\n", - "\n", - "\n", - "# parameter values (not physically based, for example only!)\n", - "param = pybamm.ParameterValues(\n", - " {\n", - " \"Reaction rate constant [m.s-1]\": 1e-6,\n", - " \"Initial thickness [m]\": 1e-6,\n", - " \"Partial molar volume [m3.mol-1]\": 10,\n", - " \"Bulk electrolyte solvent concentration [mol.m-3]\": 1,\n", - " \"Diffusivity [m2.s-1]\": Diffusivity,\n", - " }\n", - ")\n", - "\n", - "# process model and geometry\n", - "param.process_model(model)\n", - "param.process_geometry(geometry)\n", - "\n", - "# mesh and discretise\n", - "submesh_types = {\"SEI layer\": pybamm.Uniform1DSubMesh}\n", - "var_pts = {xi: 100}\n", - "mesh = pybamm.Mesh(geometry, submesh_types, var_pts)\n", - "\n", - "spatial_methods = {\"SEI layer\": pybamm.FiniteVolume()}\n", - "disc = pybamm.Discretisation(mesh, spatial_methods)\n", - "disc.process_model(model)" - ] - }, - { - "cell_type": "code", - "execution_count": 14, - "metadata": {}, - "outputs": [], - "source": [ - "# solve\n", - "solver = pybamm.ScipySolver()\n", - "t = [0, 100] # solve for 100s\n", - "solution = solver.solve(model, t)\n", - "\n", - "# post-process output variables\n", - "L_out = solution[\"SEI thickness [m]\"]\n", - "c_out = solution[\"Solvent concentration [mol.m-3]\"]" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Using these outputs, we can now plot the SEI thickness as a function of time and also the solvent concentration profile within the SEI. We use a slider to plot the concentration profile at different times. Note that, even though our model is written in nondimensional form, the processed variables are functions of dimensional space and time (in SI units). " - ] - }, - { - "cell_type": "code", - "execution_count": 15, - "metadata": {}, - "outputs": [ - { - "data": { - "application/vnd.jupyter.widget-view+json": { - "model_id": "0f4941d4712049e494267074dca70b4b", - "version_major": 2, - "version_minor": 0 - }, - "text/plain": [ - "interactive(children=(FloatSlider(value=0.0, description='t'), Output()), _dom_classes=('widget-interact',))" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "import matplotlib.pyplot as plt\n", - "\n", - "# plot SEI thickness in microns as a function of t in microseconds\n", - "# and concentration in mol/m3 as a function of x in microns\n", - "L_0_eval = param.evaluate(L_0)\n", - "xi = np.linspace(0, 1, 100) # dimensionless space\n", - "\n", - "\n", - "def plot(t):\n", - " _, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5))\n", - " ax1.plot(solution.t, L_out(solution.t) * 1e6)\n", - " ax1.plot(t, L_out(t) * 1e6, \"r.\")\n", - " ax1.set_ylabel(r\"SEI thickness [$\\mu$m]\")\n", - " ax1.set_xlabel(r\"t [s]\")\n", - "\n", - " ax2.plot(xi * L_out(t) * 1e6, c_out(t, xi))\n", - " ax2.set_ylim(0, 1.1)\n", - " ax2.set_xlim(0, L_out(solution.t[-1]) * 1e6)\n", - " ax2.set_ylabel(\"Solvent concentration [mol.m-3]\")\n", - " ax2.set_xlabel(r\"x [$\\mu$m]\")\n", - "\n", - " plt.tight_layout()\n", - " plt.show()\n", - "\n", - "\n", - "import ipywidgets as widgets\n", - "\n", - "widgets.interact(\n", - " plot, t=widgets.FloatSlider(min=0, max=solution.t[-1], step=0.1, value=0)\n", - ");" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Formally adding your model" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The purpose of this notebook has been to go through the steps involved in getting a simple model working within PyBaMM. However, if you plan on reusing your model and want greater flexibility then we recommend that you create a new class for your model. We have set out instructions on how to do this in the \"Adding a Model\" tutorial in the documentation. " - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## References\n", - "\n", - "The relevant papers for this notebook are:" - ] - }, - { - "cell_type": "code", - "execution_count": 16, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "[1] Joel A. E. Andersson, Joris Gillis, Greg Horn, James B. Rawlings, and Moritz Diehl. CasADi – A software framework for nonlinear optimization and optimal control. Mathematical Programming Computation, 11(1):1–36, 2019. doi:10.1007/s12532-018-0139-4.\n", - "[2] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, and others. Array programming with NumPy. Nature, 585(7825):357–362, 2020. doi:10.1038/s41586-020-2649-2.\n", - "[3] Valentin Sulzer, Scott G. Marquis, Robert Timms, Martin Robinson, and S. Jon Chapman. Python Battery Mathematical Modelling (PyBaMM). Journal of Open Research Software, 9(1):14, 2021. doi:10.5334/jors.309.\n", - "[4] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, and others. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature Methods, 17(3):261–272, 2020. doi:10.1038/s41592-019-0686-2.\n", - "\n" - ] - } - ], - "source": [ - "pybamm.print_citations()" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3.9.13 ('conda_jl')", - "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.11.4" - }, - "vscode": { - "interpreter": { - "hash": "612adcc456652826e82b485a1edaef831aa6d5abc680d008e93d513dd8724f14" - } - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/docs/source/examples/notebooks/creating_models/old_notebooks/7-creating-a-submodel.ipynb b/docs/source/examples/notebooks/creating_models/old_notebooks/7-creating-a-submodel.ipynb deleted file mode 100644 index 6303d6bbae..0000000000 --- a/docs/source/examples/notebooks/creating_models/old_notebooks/7-creating-a-submodel.ipynb +++ /dev/null @@ -1,486 +0,0 @@ -{ - "cells": [ - { - "attachments": {}, - "cell_type": "markdown", - "id": "86a13ed9", - "metadata": {}, - "source": [ - "Creating a submodel\n", - "\n", - "In [Tutorial 3](./3-negative-particle-problem.ipynb) we showed how to create a simple model from scratch in PyBaMM. In this tutorial we will solve a very similar problem, but using separate submodels for the linear diffusion problem and the model for the surface flux. In this simple example the surface flux is just some known function of the concentration, so we could just explicitly define it in the model for diffusion. However, we write it as a separate model to show how submodels interact. \n", - "\n", - "We solved the problem of linear diffusion on a unit sphere with a flux at the boundary that depends on the concentration\n", - "$$\n", - " \\frac{\\partial c}{\\partial t} = \\nabla \\cdot (\\nabla c),\n", - "$$\n", - "with the following boundary and initial conditions:\n", - "$$\n", - " \\left.\\frac{\\partial c}{\\partial r}\\right\\vert_{r=0} = 0, \\quad \\left.\\frac{\\partial c}{\\partial r}\\right\\vert_{r=1} = -j, \\quad \\left.c\\right\\vert_{t=0} = c_0,\n", - "$$\n", - "where\n", - "$$\n", - "j = \\left.j_0(1-c)^{1/2}c^{1/2}\\right\\vert_{r=1}\n", - "$$\n", - "Here $c_0$ and $j_0$ are parameters we can control. Again we will assume that everything is non-dimensional and focus on how to set up and solve the model rather than any specific physical interpretation." - ] - }, - { - "cell_type": "code", - "execution_count": 1, - "id": "9e8cf744", - "metadata": {}, - "outputs": [], - "source": [ - "%pip install \"pybamm[plot,cite]\" -q # install PyBaMM if it is not installed\n", - "import pybamm" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "id": "2d074d41", - "metadata": {}, - "source": [ - "## Setting up the model\n", - "\n", - "Again we start with an empty `BaseModel` class" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "id": "1dff6cd1", - "metadata": {}, - "outputs": [], - "source": [ - "model = pybamm.BaseModel()" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "id": "c44d3484", - "metadata": {}, - "source": [ - "Next we set up the submodel for diffusion in the particle using a `pybamm.BaseSubModel`. Each submodel has methods that define the variables, equations, initial and boundary conditions corresponding to the physics the model describes.\n", - "\n", - "First `get_fundamental_variables` defines any variables that can be defined independently of any other submodels that may be included in our model. Here we can define the concentration, surface concentration and flux, and add them to the dictionary of variables. \n", - "\n", - "Next we can use `get_coupled_variables` to define any variables that _do_ depend on variables defined in another submodel. In this simple example we don't have any variables to define here. However, if we had included a temperature dependent diffusivity, for example, then we would have needed to define the flux in `get_coupled_variables` since it would now depend on the temperature which would be defined in a separate submodel.\n", - "\n", - "Once we have defined all the variables we need we can write down our equations. Any equations that include time derivatives will turn into ordinary differential equations after discretisation. We set the right hand sides of such equations in the `set_rhs` method. In this example we add the right hand side of the diffusion equation to the `rhs` dictionary. Equations that don't contain time derivatives give algebraic constraints in our model. These equations are set in the `set_algebraic` method. In this example we don't have any algebraic equations, so we can skip this method.\n", - "\n", - "Finally we set the boundary and initial conditions using the methods `set_boundary_conditions` and `set_initial_conditions`, respectively. " - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "id": "337c1475", - "metadata": {}, - "outputs": [], - "source": [ - "class Particle(pybamm.BaseSubModel):\n", - " def __init__(self, param, domain, options=None):\n", - " super().__init__(param, domain, options=options)\n", - "\n", - " def get_fundamental_variables(self):\n", - " # create concentration variable\n", - " c = pybamm.Variable(\"Concentration\", domain=\"negative particle\")\n", - "\n", - " # define concentration at the surface of the sphere\n", - " c_surf = pybamm.surf(c)\n", - "\n", - " # define flux\n", - " N = -pybamm.grad(c)\n", - "\n", - " # create dictionary of model variables\n", - " variables = {\n", - " \"Concentration\": c,\n", - " \"Surface concentration\": c_surf,\n", - " \"Flux\": N,\n", - " }\n", - "\n", - " return variables\n", - "\n", - " def get_coupled_variables(self, variables):\n", - " return variables\n", - "\n", - " def set_rhs(self, variables):\n", - " # extract the variables we need\n", - " c = variables[\"Concentration\"]\n", - " N = variables[\"Flux\"]\n", - "\n", - " # define the rhs of the PDE\n", - " dcdt = -pybamm.div(N)\n", - "\n", - " # add it to the submodel dictionary\n", - " self.rhs = {c: dcdt}\n", - "\n", - " def set_algebraic(self, variables):\n", - " pass\n", - "\n", - " def set_boundary_conditions(self, variables):\n", - " # extract the variables we need\n", - " c = variables[\"Concentration\"]\n", - " j = variables[\"Boundary flux\"]\n", - "\n", - " # add the boundary conditions to the submodel dictionary\n", - " self.boundary_conditions = {\n", - " c: {\"left\": (0, \"Neumann\"), \"right\": (-j, \"Neumann\")}\n", - " }\n", - "\n", - " def set_initial_conditions(self, variables):\n", - " # extract the variable we need\n", - " c = variables[\"Concentration\"]\n", - "\n", - " # define the initial concentration parameter\n", - " c0 = pybamm.Parameter(\"Initial concentration\")\n", - "\n", - " # add the initial conditions to the submodel dictionary\n", - " self.initial_conditions = {c: c0}" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "id": "642bbd46", - "metadata": {}, - "outputs": [], - "source": [ - "class BoundaryFlux(pybamm.BaseSubModel):\n", - " def __init__(self, param, domain, options=None):\n", - " super().__init__(param, domain, options=options)\n", - "\n", - " def get_coupled_variables(self, variables):\n", - " # extract the variable we need\n", - " c_surf = variables[\"Surface concentration\"]\n", - "\n", - " # define the flux parameter\n", - " j0 = pybamm.Parameter(\"Flux parameter\")\n", - " j = j0 * (1 - c_surf) ** (1 / 2) * c_surf ** (1 / 2) # prescribed boundary flux\n", - "\n", - " # update dictionary of model variables\n", - " variables.update({\"Boundary flux\": j})\n", - "\n", - " return variables" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "id": "62907bc5", - "metadata": {}, - "source": [ - "We can now set the submodels in a model by assigning a dictionary to `model.submodels`. The dictionary key is the name we want to give to the submodel and the value is an instance of the submodel class we want to use.\n", - "\n", - "When we instantiate a submodel we are required to pass in `param`, a class of parameter symbols we are going to call, and `domain`, the domain on which the submodel lives. In this example we will simply set `param` to `None` and hard-code the definition of our parameters into the submodel. When writing lots of submodels it is more efficient to define _all_ the parameters in a shared class, and pass this to each submodel. For the domain we will choose \"Negative\". " - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "id": "7fed4b8a", - "metadata": {}, - "outputs": [], - "source": [ - "model.submodels = {\n", - " \"Particle\": Particle(None, \"Negative\"),\n", - " \"Boundary flux\": BoundaryFlux(None, \"Negative\"),\n", - "}" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "id": "9e910625", - "metadata": {}, - "source": [ - "At this stage we have just told the model which submodels it is constructed from, but the variables and equations have not yet been created. For example if we look at the `rhs` dictionary it is empty." - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "id": "93541aa0", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "{}" - ] - }, - "execution_count": 6, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "model.rhs" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "id": "c6e07559", - "metadata": {}, - "source": [ - "To populate the model variables, equations, boundary and initial conditions we need to \"build\" the model. To do this we call `build_model`" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "id": "1a8b9548", - "metadata": {}, - "outputs": [], - "source": [ - "model.build_model()" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "id": "649c0260", - "metadata": {}, - "source": [ - "This loops through all of the submodels, first creating the \"fundamental variables\", followed by the \"coupled variables\" and finally the equations (`rhs` and `algebraic`) and the boundary and initial conditions. Now we see that `model.rhs` contains our diffusion equation." - ] - }, - { - "cell_type": "code", - "execution_count": 8, - "id": "374fae5c", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "{Variable(0x48a91fa41dea72d3, Concentration, children=[], domains={'primary': ['negative particle']}): Divergence(0x4befa5e7cf20d0c5, div, children=['grad(Concentration)'], domains={'primary': ['negative particle']})}" - ] - }, - "execution_count": 8, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "model.rhs" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "id": "a7cab3ac", - "metadata": {}, - "source": [ - "## Using the model\n", - "\n", - "We can now use our model as in the previous tutorial.\n", - "\n", - "We first set up our geometry and mesh" - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "id": "2f2da3bd", - "metadata": {}, - "outputs": [], - "source": [ - "r = pybamm.SpatialVariable(\n", - " \"r\", domain=[\"negative particle\"], coord_sys=\"spherical polar\"\n", - ")\n", - "geometry = {\"negative particle\": {r: {\"min\": 0, \"max\": 1}}}\n", - "spatial_methods = {\"negative particle\": pybamm.FiniteVolume()}\n", - "submesh_types = {\"negative particle\": pybamm.Uniform1DSubMesh}\n", - "var_pts = {r: 20}\n", - "mesh = pybamm.Mesh(geometry, submesh_types, var_pts)" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "id": "2a89e132", - "metadata": {}, - "source": [ - "and then set up our simulation, remembering to set values for our parameters" - ] - }, - { - "cell_type": "code", - "execution_count": 10, - "id": "7589648c", - "metadata": {}, - "outputs": [], - "source": [ - "parameter_values = pybamm.ParameterValues(\n", - " {\n", - " \"Initial concentration\": 0.9,\n", - " \"Flux parameter\": 0.8,\n", - " }\n", - ")\n", - "\n", - "solver = pybamm.ScipySolver()\n", - "\n", - "sim = pybamm.Simulation(\n", - " model,\n", - " geometry=geometry,\n", - " parameter_values=parameter_values,\n", - " submesh_types=submesh_types,\n", - " var_pts=var_pts,\n", - " spatial_methods=spatial_methods,\n", - " solver=solver,\n", - ")" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "id": "1cec6e66", - "metadata": {}, - "source": [ - "Finally we can solve the model" - ] - }, - { - "cell_type": "code", - "execution_count": 11, - "id": "cc9bdc87", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "" - ] - }, - "execution_count": 11, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "sim.solve([0, 1])" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "id": "e6bafd61", - "metadata": {}, - "source": [ - "and plot the results" - ] - }, - { - "cell_type": "code", - "execution_count": 12, - "id": "f549d023", - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "2022-07-11 13:20:04.194 - [WARNING] processed_variable.get_spatial_scale(521): No length scale set for negative particle. Using default of 1 [m].\n", - "2022-07-11 13:20:04.209 - [WARNING] processed_variable.get_spatial_scale(521): No length scale set for negative particle. Using default of 1 [m].\n" - ] - }, - { - "data": { - "application/vnd.jupyter.widget-view+json": { - "model_id": "d9780b7718bd4a9f932f5dd92c065756", - "version_major": 2, - "version_minor": 0 - }, - "text/plain": [ - "interactive(children=(FloatSlider(value=0.0, description='t', max=1.0, step=0.01), Output()), _dom_classes=('w…" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/plain": [ - "" - ] - }, - "execution_count": 12, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "# pass in a list of the variables we want to plot\n", - "sim.plot([\"Concentration\", \"Surface concentration\", \"Flux\", \"Boundary flux\"])" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "id": "3d88c346", - "metadata": {}, - "source": [ - "In this notebook we saw how to split a model up into submodels. Although this was a simple example it let us understand how to construct submodels and see how they interact via coupled variables." - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "id": "f20b7969", - "metadata": {}, - "source": [ - "## References\n", - "\n", - "The relevant papers for this notebook are:" - ] - }, - { - "cell_type": "code", - "execution_count": 13, - "id": "6301f16b", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "[1] Joel A. E. Andersson, Joris Gillis, Greg Horn, James B. Rawlings, and Moritz Diehl. CasADi – A software framework for nonlinear optimization and optimal control. Mathematical Programming Computation, 11(1):1–36, 2019. doi:10.1007/s12532-018-0139-4.\n", - "[2] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, and others. Array programming with NumPy. Nature, 585(7825):357–362, 2020. doi:10.1038/s41586-020-2649-2.\n", - "[3] Valentin Sulzer, Scott G. Marquis, Robert Timms, Martin Robinson, and S. Jon Chapman. Python Battery Mathematical Modelling (PyBaMM). Journal of Open Research Software, 9(1):14, 2021. doi:10.5334/jors.309.\n", - "[4] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, and others. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature Methods, 17(3):261–272, 2020. doi:10.1038/s41592-019-0686-2.\n", - "\n" - ] - } - ], - "source": [ - "pybamm.print_citations()" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.8.13" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -}