From a27e8648ba84bf1b460def7e3b29c26ebd2be1c0 Mon Sep 17 00:00:00 2001 From: Robert Timms Date: Thu, 25 Nov 2021 16:37:46 +0000 Subject: [PATCH] #1311 comments, add resistor notebook --- .../notebooks/models/jelly-roll-model.ipynb | 385 ++++++++++++++++++ pybamm/geometry/standard_spatial_vars.py | 8 +- pybamm/spatial_methods/finite_volume.py | 45 +- 3 files changed, 409 insertions(+), 29 deletions(-) create mode 100644 examples/notebooks/models/jelly-roll-model.ipynb diff --git a/examples/notebooks/models/jelly-roll-model.ipynb b/examples/notebooks/models/jelly-roll-model.ipynb new file mode 100644 index 0000000000..b5a932a607 --- /dev/null +++ b/examples/notebooks/models/jelly-roll-model.ipynb @@ -0,0 +1,385 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "cathedral-trance", + "metadata": {}, + "source": [ + "# Jelly roll model\n", + "\n", + "In this notebook we show how to set up and solve the \"two-potential\" model from \"Homogenisation of spirally-wound high-contrast layered materials\", S. Psaltis, R. Timms, C.P. Please, S.J. Chapman, SIAM Journal on Applied Mathematics, 2020.\n", + "\n", + "We consider a spirally-wound cell, such as the common 18650 lithium-ion cell. In practice these cells are constructed by rolling a sandwich of layers containing the active cathode, positive current collector, active cathode, separator, active anode, negative current collector, active anode, and separator. The \"two-potential\" model consists of an equation for the potential $\\phi^\\pm$ in each current collector. The potential difference drives a current $I$ through the electrode/separator/electrode sandwich (which we refer to as the \"active material\" in the original paper). Thus, in non-dimensional form, the model is \n", + "\n", + "$$ \\frac{\\delta^+\\sigma^+}{2\\pi^2}\\frac{1}{r}\\frac{\\mathrm{d}}{\\mathrm{d}r}\\left(\\frac{1}{r}\\frac{\\mathrm{d}\\phi^+}{\\mathrm{d}r}\\right) + 2I(\\phi^+-\\phi^-) = 0,$$\n", + "$$ \\frac{\\delta^-\\sigma^-}{2\\pi^2}\\frac{1}{r}\\frac{\\mathrm{d}}{\\mathrm{d}r}\\left(\\frac{1}{r}\\frac{\\mathrm{d}\\phi^-}{\\mathrm{d}r}\\right) - 2I(\\phi^+-\\phi^-) = 0,$$\n", + "with boundary conditions \n", + "$$ \\frac{\\mathrm{d}\\phi^+}{\\mathrm{d}r}(r=r_0) = 0, \\quad \\phi^+(r=1) = 1, \\quad \\phi^-(r=0) = 0, \\quad \\frac{\\mathrm{d}\\phi^-}{\\mathrm{d}r}(r=1) = 0.$$\n", + "\n", + "For a complete description of the model and parameters, please refer to the original paper.\n", + "\n", + "It can be shown that the active material can be modelled using any 1D battery model we like to describe the electrochemical/thermal behaviour in the electrode/separator/electrode sandwich. Such functionality will be added to PyBaMM in a future release and will enable efficient simulations of jelly roll cells. \n" + ] + }, + { + "cell_type": "markdown", + "id": "whole-diabetes", + "metadata": {}, + "source": [ + "## Two-potential resistor model\n", + "In this section we consider a simplified model in which we ignore the details of the anode, cathode and separator, and treat them as a single region of active material, modelled as an Ohmic conductor, with two such regions per winding. In this case the model becomes \n", + "\n", + "$$ \\frac{\\delta^+\\sigma^+}{2\\pi^2}\\frac{1}{r}\\frac{\\mathrm{d}}{\\mathrm{d}r}\\left(\\frac{1}{r}\\frac{\\mathrm{d}\\phi^+}{\\mathrm{d}r}\\right) + \\frac{2\\sigma^{a}(\\phi^--\\phi^+)}{l\\epsilon^4} = 0,$$\n", + "$$ \\frac{\\delta^-\\sigma^-}{2\\pi^2}\\frac{1}{r}\\frac{\\mathrm{d}}{\\mathrm{d}r}\\left(\\frac{1}{r}\\frac{\\mathrm{d}\\phi^-}{\\mathrm{d}r}\\right) - \\frac{2\\sigma^{a}(\\phi^--\\phi^+)}{l\\epsilon^4} = 0,$$\n", + "along with the same boundary conditions.\n", + "\n", + "We begin by importing PyBaMM along with some other useful packages" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "heard-cartridge", + "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 -q # install PyBaMM if it is not installed\n", + "import pybamm\n", + "import numpy as np \n", + "from numpy import pi\n", + "import matplotlib.pyplot as plt " + ] + }, + { + "cell_type": "markdown", + "id": "steady-chest", + "metadata": {}, + "source": [ + "First we will define the parameters in the model. Note the model is posed in non-dimensional form." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "rising-executive", + "metadata": {}, + "outputs": [], + "source": [ + "N = pybamm.Parameter(\"Number of winds\")\n", + "r0 = pybamm.Parameter(\"Inner radius\")\n", + "eps = (1 - r0) / N # ratio of sandwich thickness to cell radius\n", + "delta = pybamm.Parameter(\"Current collector thickness\")\n", + "delta_p = delta # assume same thickness\n", + "delta_n = delta # assume same thickness\n", + "l = 1/2 - delta_p - delta_n # active material thickness\n", + "sigma_p = pybamm.Parameter(\"Positive current collector conductivity\")\n", + "sigma_n = pybamm.Parameter(\"Negative current collector conductivity\")\n", + "sigma_a = pybamm.Parameter(\"Active material conductivity\")" + ] + }, + { + "cell_type": "markdown", + "id": "worth-easter", + "metadata": {}, + "source": [ + "Next we define our geometry and model" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "virgin-wrestling", + "metadata": {}, + "outputs": [], + "source": [ + "# geometry\n", + "r = pybamm.SpatialVariable(\"radius\", domain=\"cell\", coord_sys=\"cylindrical polar\")\n", + "geometry = {\"cell\": {r: {\"min\": r0, \"max\": 1}}}\n", + "\n", + "# model\n", + "model = pybamm.BaseModel()\n", + "phi_p = pybamm.Variable(\"Positive potential\", domain=\"cell\")\n", + "phi_n = pybamm.Variable(\"Negative potential\", domain=\"cell\")\n", + "\n", + "A_p = (2 * sigma_a / eps ** 4 / l) / (delta_p * sigma_p / 2 / pi ** 2)\n", + "A_n = (2 * sigma_a / eps ** 4 / l) / (delta_n * sigma_n / 2 / pi ** 2)\n", + "model.algebraic = {\n", + " phi_p: pybamm.div((1 / r ** 2) * pybamm.grad(phi_p)) + A_p * (phi_n - phi_p),\n", + " phi_n: pybamm.div((1 / r ** 2) * pybamm.grad(phi_n)) - A_n * (phi_n - phi_p),\n", + "}\n", + "\n", + "model.boundary_conditions = {\n", + " phi_p: {\n", + " \"left\": (0, \"Neumann\"),\n", + " \"right\": (1, \"Dirichlet\"),\n", + " },\n", + " phi_n: {\n", + " \"left\": (0, \"Dirichlet\"),\n", + " \"right\": (0, \"Neumann\"),\n", + " } \n", + "}\n", + "\n", + "model.initial_conditions = {phi_p: 1, phi_n: 0} # initial guess for solver\n", + "\n", + "model.variables = {\"Negative potential\": phi_n, \"Positive potential\": phi_p}" + ] + }, + { + "cell_type": "markdown", + "id": "based-slope", + "metadata": {}, + "source": [ + "Next we provide values for our parameters, and process our geometry and model, thus replacing the `Parameter` symbols with numerical values" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "technological-electric", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "params = pybamm.ParameterValues(\n", + " {\n", + " \"Number of winds\":20,\n", + " \"Inner radius\": 0.25,\n", + " \"Current collector thickness\": 0.05,\n", + " \"Positive current collector conductivity\": 5e6,\n", + " \"Negative current collector conductivity\": 5e6,\n", + " \"Active material conductivity\": 1,\n", + " }\n", + ")\n", + "params.process_geometry(geometry)\n", + "params.process_model(model)" + ] + }, + { + "cell_type": "markdown", + "id": "polyphonic-opinion", + "metadata": {}, + "source": [ + "We choose to discretise in space using the Finite Volume method on a uniform grid" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "buried-blind", + "metadata": {}, + "outputs": [], + "source": [ + "# mesh\n", + "submesh_types = {\"cell\": pybamm.Uniform1DSubMesh}\n", + "var_pts = {r: 100}\n", + "mesh = pybamm.Mesh(geometry, submesh_types, var_pts)\n", + "# method\n", + "spatial_methods = {\"cell\": pybamm.FiniteVolume()}\n", + "# discretise\n", + "disc = pybamm.Discretisation(mesh, spatial_methods)\n", + "disc.process_model(model);" + ] + }, + { + "cell_type": "markdown", + "id": "gothic-deadline", + "metadata": {}, + "source": [ + "We can now solve the model" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "straight-anime", + "metadata": {}, + "outputs": [], + "source": [ + "# solver \n", + "solver = pybamm.CasadiAlgebraicSolver()\n", + "solution = solver.solve(model)" + ] + }, + { + "cell_type": "markdown", + "id": "excessive-universal", + "metadata": {}, + "source": [ + "The model gives the homogenised potentials in the negative a positive current collectors. Interestingly, the solid potential has microscale structure, varying linearly in the active material. In order to see this we need to post-process the solution and plot the potential as a function of radial position, being careful to capture the spiral geometry. " + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "invisible-laser", + "metadata": {}, + "outputs": [], + "source": [ + "# extract numerical parameter values\n", + "# Note: this overrides the definition of the `pybamm.Parameter` objects\n", + "N = params.evaluate(N)\n", + "r0 = params.evaluate(r0)\n", + "eps = params.evaluate(eps)\n", + "delta = params.evaluate(delta)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "affecting-albuquerque", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "2021-11-25 16:36:14,089 - [WARNING] processed_variable.get_spatial_scale(520): No length scale set for cell. Using default of 1 [m].\n", + "2021-11-25 16:36:14,091 - [WARNING] processed_variable.get_spatial_scale(520): No length scale set for cell. Using default of 1 [m].\n" + ] + } + ], + "source": [ + "# post-process homogenised potential \n", + "phi_n = solution[\"Negative potential\"]\n", + "phi_p = solution[\"Positive potential\"]\n", + "\n", + "def alpha(r):\n", + " return 2 * (phi_n(x=r) - phi_p(x=r))\n", + "\n", + "def phi_am1(r, theta):\n", + " # careful here - phi always returns a column vector so we need to add a new axis to r to get the right shape \n", + " return alpha(r) * (r[:,np.newaxis]/eps - r0/eps - delta - theta / 2 / pi) / (1 - 4*delta) + phi_p(x=r)\n", + "\n", + "def phi_am2(r, theta):\n", + " # careful here - phi always returns a column vector so we need to add a new axis to r to get the right shape \n", + " return alpha(r) * (r0/eps + 1 - delta + theta / 2 / pi - r[:,np.newaxis]/eps) / (1 - 4*delta) + phi_p(x=r)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "taken-hearing", + "metadata": {}, + "outputs": [], + "source": [ + "# define spiral \n", + "spiral_pos_inner = lambda t : r0 - eps * delta + eps * t / (2 * pi)\n", + "spiral_pos_outer = lambda t : r0 + eps * delta + eps * t / (2 * pi)\n", + "\n", + "spiral_neg_inner = lambda t : r0 - eps * delta + eps/2 + eps * t / (2 * pi)\n", + "spiral_neg_outer = lambda t : r0 + eps * delta + eps/2 + eps * t / (2 * pi)\n", + "\n", + "spiral_am1_inner = lambda t : r0 + eps * delta + eps * t / (2 * pi)\n", + "spiral_am1_outer = lambda t : r0 - eps * delta + eps/2 + eps * t / (2 * pi)\n", + "\n", + "spiral_am2_inner = lambda t : r0 + eps * delta + eps/2 + eps * t / (2 * pi)\n", + "spiral_am2_outer = lambda t : r0 - eps * delta + eps + eps * t / (2 * pi)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "handled-jacksonville", + "metadata": {}, + "outputs": [], + "source": [ + "# Setup fine mesh with nr points per layer\n", + "nr = 10\n", + "rr = np.linspace(r0, 1, nr)\n", + "tt = np.arange(0, (N+1)*2*pi, 2*pi)\n", + "# N+1 winds of pos c.c.\n", + "r_mesh_pos = np.zeros((len(tt),len(rr)))\n", + "for i in range(len(tt)):\n", + " r_mesh_pos[i,:] = np.linspace(spiral_pos_inner(tt[i]), spiral_pos_outer(tt[i]), nr)\n", + "# N winds of neg, am1, am2\n", + "r_mesh_neg = np.zeros((len(tt)-1, len(rr)))\n", + "r_mesh_am1 = np.zeros((len(tt)-1, len(rr)))\n", + "r_mesh_am2 = np.zeros((len(tt)-1, len(rr)))\n", + "for i in range(len(tt)-1):\n", + " r_mesh_am2[i,:] = np.linspace(spiral_am2_inner(tt[i]), spiral_am2_outer(tt[i]), nr)\n", + " r_mesh_neg[i,:] = np.linspace(spiral_neg_inner(tt[i]), spiral_neg_outer(tt[i]), nr)\n", + " r_mesh_am1[i,:] = np.linspace(spiral_am1_inner(tt[i]), spiral_am1_outer(tt[i]), nr)\n", + "# Combine and sort \n", + "r_total_mesh = np.vstack((r_mesh_pos,r_mesh_neg,r_mesh_am1, r_mesh_am2))\n", + "r_total_mesh = np.sort(r_total_mesh,axis=None)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "monetary-belarus", + "metadata": { + "scrolled": false + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfUAAAF0CAYAAAA6pKBsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABKsUlEQVR4nO3deZxPZf/H8dc1q52xjp0soeyDSu5yW7JlvbMXJVuENrfSjdzdJS0/UsoSIZE9WVIixJCdLNmXGfs2ljH79fvjDI0MZpj5fme+3s/H4/uYmXOu7/l+jqF313Wucy5jrUVERETSPy93FyAiIiIpQ6EuIiLiIRTqIiIiHkKhLiIi4iEU6iIiIh5CoS4iIuIhfNxdwL3KnTu3LVasmLvLEBERcYmNGzeesdbmSWxfug/1YsWKsWHDBneXISIi4hLGmMO32qfhdxEREQ+hUBcREfEQCnUREREPke6vqScmOjqakJAQIiIi3F3KPcuQIQOFChXC19fX3aWIiEga55GhHhISQtasWSlWrBjGGHeXc9estZw9e5aQkBCKFy/u7nJERCSN88jh94iICHLlypWuAx3AGEOuXLk8YsRBRERSn0eGOpDuA/0aTzkPERFJfR4b6iIiIvcbhbqIiIiH8MiJcmnN3LlzWbZsGaNGjbq+bfv27bz55ps3tJswYQJ58+Z1dXkiIuIhFOousGnTJqpUqXLDtvLly7NgwQI3VSQiIp7IZaFujJkANAFOWWsfTmS/AUYCjYBwoLO1dpOr6ksNe/bsoVevXqxdu5ZcuXIRFhZGv3793F2WiIi4yM7587l46hQ1XngB45X6V7xd2VP/GvgMmHyL/Q2BUvGvGsAX8V/vSb9+sGXLvR7lRpUqwYgRt28TGRlJ69atmTJlCs2aNWPNmjWUK1eOHj16kCFDhpQtSERE0qT3+/blh0OHON6yJRlz5kz1z3PZRDlr7Urg3G2aNAMmW8daIIcxJr9rqkt5P//8MxUrVqRAgQJky5aNwMBAMmTIQGxsrLtLExERFzh/8CCzDh2i/UMPuSTQIW1dUy8IHE3wc0j8tuP3ctA79ahTy9atWylfvjzbtm2jQoUKnDp1iqxZs5I5c2b3FCQiIi41tX9/IoAX33rLZZ+ZlkI9yYwx3YBuAEWKFHFzNYnLmjUr27Ztw8fHhwoVKjB48GB69erl7rJERMQFbFwc4xcsoHLGjFRp395ln5uW7lMPBQon+LlQ/LabWGvHWmuDrLVBefLkcUlxydWxY0f27t3L0KFD+eKLL8iZMycvv/yyu8sSEREXWDJqFFsjIujarJlLPzct9dTnA72NMdNxJsiFWWvvaejdnXLmzMmKFSuoVKkSS5cuJXfu3O4uSUREXOS1IUMAaP/BBy79XFfe0jYNeBLIbYwJAQYDvgDW2i+BRTi3s+3DuaXteVfVlloiIyMJCwtToIuI3E8uX6bgpUuEZcxIdhdfInZZqFtr291hvwU86qKzv78/Bw8edHcZIiLiSjNncjU2llJlyrj8o9PSNXUREZH0b9w4Tvr6kq90aZd/tEJdREQkpezYAcHBnPDyIjAw0OUfr1AXERFJKV99RbiPD5ciI8mXL5/LP16hLiIikhIiI2HyZE7WqwegUBcREUm35s2Ds2c5+dRTABp+FxERSbfGj4eiRTlR2HmOmnrq96Ht27dz5MgRd5chIiL34sABWLoUXniBk6dPAwr1+9LGjRs5cOCAu8sQEZF7MWYMeHvDCy9w4sQJAPLmzevyMhTqLjB37tybnvu+c+dOevTowaRJk/jkk0/o0aMHJ0+edFOFIiJy1yIi4KuvoGlTKFSIEydOkCtXLvz8/FxeSlp69rvH2rRpE1WqVLlhW7ly5fjyyy/5+uuvKVasGE8++aR7ihMRkXszaxacPQsvvQTA8ePHyZ8/v1tKUainoj179tCrVy/Wrl1Lrly5CAsLo1+/fu4uS0REUtLo0VC6NPzznwCcOHHCbaGu4fdUEhkZSevWrfnkk0/IkycPa9euZejQoURERNzQrnPnzuqli4ikV5s3Q3Aw9OwJXk6kHj9+3C23s8H90FPv1w+2bEnZY1aqBCNG3LbJzz//TMWKFSlQoADZsmUjMDCQDBkyEBsbm7K1iIiI+3zxBWTMCJ06AWCtVU/dE23dupXy5cuzbds2KlSowKlTp8iaNSuZM2d2d2kiIpISwsJg6lRo1w4CAgA4f/48UVFR6qmnmjv0qFNL1qxZ2bZtGz4+PlSoUIHBgwfTq5dHrSwrInJ/mzwZwsOvT5ADrt/OpolyHqZjx460aNGCOXPmEBAQQNu2bW+6rU1ERNIpa50JctWrQ9Wq1zcfP34ccM8jYkGhnmpy5szJihUrqFSpEkuXLiV37tzuLklERFLKr7/C7t3w9dc3bHZ3T13X1FNRZGQkYWFhCnQREU/z+efOdfTWrW/YfOzYMQAKFCjgjqoU6qnJ39+fgwcPursMERFJSQcOwNy50L27M/M9gePHj5M5c2ayZs3qltIU6iIiIskxcqRzT3rv3jftOnbsmNt66aBQFxERSboLF5znvLdrBwUL3rRboS4iIpJejB0LV67Aq68muluhLiIikh5ER8OnnzrPeK9U6abd1lq3h7puaRMREUmKGTMgNNRZOz0RYWFhXL16VT11ERGRNM1a+PhjKFMGGjZMtIm7b2cD9dRFRETubMUKZ0W2MWOur8b2d9dC3V0PngH11EVERO7sgw8gTx549tlbNgkNDQWgYCKz4l1FoS4iInI7GzfCjz86M97/9rCZhBTq94m5c+dqMRcRkfTq/fche3bo2fO2zUJCQsiZMycZbxP8qU2h7gKbNm2iSpUq7i5DRESSa9cumDPHeXpc9uy3bRoaGurWXjpoolyq2rNnD7169WLt2rXkypWLsLAw+vXr5+6yREQkqYYNc4bc+/a9Y9PQ0FAKFSrkgqJuzeNDvV+/fmzZsiVFj1mpUiVGjBhx2zaRkZG0bt2aKVOm0KxZM9asWUO5cuXo0aMHGTJkSNF6REQkFRw6BFOnwssvO5Pk7iA0NJTKlSunfl23oeH3VPLzzz9TsWJFChQoQLZs2QgMDCRDhgzExsa6uzQREUmKDz90bl977bU7No2OjubkyZMafk9td+pRp5atW7dSvnx5tm3bRoUKFTh16hRZs2Ylc+bMbqlHRESS4dgxZ+GWTp0gCUPqx48fx1qr4XdPlTVrVrZt24aPjw8VKlRg8ODB9OrVy91liYhIUgwbBjEx8OabSWoeEhICuPd2NlCop5qOHTvSokUL5syZQ0BAAG3bttVtbSIi6UFIiPPkuOefhwceSNJbjh49CkDhwoVTs7I7Uqinkpw5c7JixQoqVarE0qVLyZ07t7tLEhGRpHjvPYiLg4EDk/yWtBLqmiiXiiIjIwkLC1Ogi4ikF0eOwPjx0KULFCuW5LcdPXqUrFmzkv0O97KnNoV6KvL39+fgwYPuLkNERJKoZvXq1IqOxibxWvo1R48edXsvHTT8LiIiAsDR1atZe/IkPl5eRObLR3KeKJJWQl09dREREWDoCy/gBUTFxTFv3rxkvVehLiIikkbs/fFHJu7ZQ48KFShSpAiTJk1K8nsjIyM5efKkQj01WWvdXUKK8JTzEBFJy+q3aoUBBk6dSvv27fn55585ffp0kt577R71IkWKpGKFSeORoZ4hQwbOnj2b7gPRWsvZs2f1rHgRkVRk16zhXHg4D+XLR+DDD9O+fXtiY2OZOXNmkt5/+PBhAIoWLZqaZSaJR06UK1SoECEhIUn+v6y0LEOGDG5/7KCIiMeylouvv85FoG38eunly5enTJkyzJ49m5deeumOh7gW6mmhp+6Roe7r60vx4sXdXYaIiKR1CxeyMzgYgIerVr2+uUWLFgwfPpyzZ8+SK1eu2x7iyJEjGGPSRAfMI4ffRURE7ig2FgYM4I/4ZVUfeuih67tatGhBbGwsCxYsuONhDh8+TGBgIP7+/qlWalIp1EVE5P40ZQrs2MEfQUFkzpz5hmviQUFB5M+fn0WLFt3xMIcPH04T19NBoS4iIvejq1fhP/+BatX4IzKScuXK4eX1VyQaY2jQoAE//fQTMTExtz3U4cOH08T1dFCoi4jI/ejTTyEkBPv++2zdupWKFSve1KRhw4ZcuHCBtWvX3vIwcXFxHD169P7sqRtjGhhj/jTG7DPGDEhkfxFjzHJjzGZjzDZjTCNX1iciIveBEyfgf/+Dp5/meNmynD17lgoVKtzUrG7dunh5efHzzz/f8lDHjx8nKiqKYslY/CU1uSzUjTHewOdAQ6Ac0M4YU+5vzd4GZlhrKwNtgdGuqk9ERO4Tb78NERHw0Uds27YNINGeekBAAFWrVuWXX3655aEOHToEkGbuuHJlT706sM9ae8BaGwVMB5r9rY0FssV/nx045sL6RETE023eDBMmwMsvQ+nSbN26FXDuTU9MnTp1WLduHZcvX050/7WVOO+7njpQEDia4OeQ+G0JDQE6GmNCgEXAy4kdyBjTzRizwRizwRMeMCMiIi5gLbzyCuTM6UySA7Zu3UrhwoUJCAhI9C116tQhJiaG3377LdH913rq92OoJ0U74GtrbSGgETDFGHNTjdbasdbaIGttUJ74+wtFRERua84cWLEC/vtfyJEDgC1btlC5cuVbvuXRRx/Fx8eHVatWJbr/0KFD5MuXj4wZM6ZGxcnmylAPBRIuYVMofltCXYAZANbaYCADkNsl1YmIiOe6ehXeeAMefhi6dgUgPDycP//8k0qVKt3ybZkzZ6Zq1aqsXLky0f0HDx5MM9fTwbWhvh4oZYwpbozxw5kIN/9vbY4AdQCMMWVxQl3j6yIicm+GD4eDB2HkSPBxnpC+fft24uLibttTB6hVqxa///47V69evWnfoUOH0szQO7gw1K21MUBvYAmwC2eW+w5jzFBjTNP4Zq8BXY0xW4FpQGeb3pdaExER9zpwAIYNgzZt4J//vL558+bNALftqQM8/vjjREVFsWnTphu2R0dHc/jwYR544IEUL/luuXRBF2vtIpwJcAm3DUrw/U6gpitrEhERD9evH3h7w8cf37B506ZNBAQE3PHBMY8++igAwcHB1Kz5V0QdPXqU2NhYSpQokeIl3y2PXKVNREQEgAUL4IcfnOH3gjfecLVp0yaqVKmCMea2h8ibNy8PPPDATU+W279/P0Ca6qmntdnvIiIiKSMiAvr2hbJlna8JREVFsX37dqomWG71dh555JGbQv3AgQMAaaqnrlAXERHPNHy4cz39s8/Az++GXTt27CAqKipZoR4aGkpo6F83be3fvx8/Pz8KFCiQomXfC4W6iIh4nj174L33bpocd82GDRsAqFKlSpIOFxQUdMP7wOmpFy9eHG9v7xQoOGUo1EVExLNYCz16QIYMMGJEok3Wr19PQEBAkofOK1asiLe39w2hvn///jR1PR0U6iIi4mmmTIHly+GDDyAwMNEm69evp1q1anecJHdNpkyZeOihh66HurWWvXv3UqpUqRQrOyUo1EVExHOcOQOvvgqPPXb9yXF/Fx4ezvbt26lWrVqyDh0UFMSGDRuw1nLy5EmuXLmiUBcREUk1r78OYWEwZgx4JR5xW7ZsITY2NtmhXqVKFc6cOUNoaCh79+4FoGTJkom2PXMGfvkFPvkEBg5M3incC92nLiIinmHZMpg0Cd5803nG+y1cuzWtevXqyTr8tSfPbdmyhWsrhJYqVYpLl2DpUtiwAbZsga1bIcEkeYoWhaFDneffpDaFuoiIpH+XL8OLL0LJkteXVb2VtWvXUrRoUfLnz5+sj6hQoQLghPrRo2fw8vLm5ZeLsmwZREY6j5QvWxZq14ZKlaBiReflysVEFeoiIpL+vfUWHDrkLK16h2VQ165dy2OPPZbsj8iaNSv58hXivfdGc/XqccCHnTu9eeklaN4catQAf/+7qj7FKNRFRCR9W7UKRo2Cl1+GWrVu2/TYsWMcPXqURx55JMmHtxZ+/NG5Pn7yZBxwkjJlXqRNm0cZNOiWl+7dQqEuIiLpV3g4vPACFC8O779/x+Zr1qwB/lqk5U7WrnUm0wcHQ4EC0KvXAvr2LUypUrnvqezUolAXEZF0q+ejj9J53z5qLFsGmTPfsf3q1avJkCHDHddQP3TImW83fbpzq/u4cfDcc+Dnd/v3uVsaGjQQERFJul8//5wvt23jH15ezL90KUnvWb16NdWrV8fvb8+CvyY62ll6vWxZ+P57GDQI9u515uDd4i1pikJdRETSnYgLF+j96qvk9fKi3EMP0axZM0bc4pGw14SHh7N58+ZbTpJbvx6CgpweeqNG8Oef8M47kCVLKpxAKlGoi4hIutOqShV2REUxaehQ1qxbR6tWrXjllVcYNWrULd+zdu1aYmJiqPW3yXTR0c5dcI88AmfPwty5MHs2FC6c2meR8hTqIiKSrhycNo1FBw9SMWdOGgwcSMaMGZk+fTrNmzenb9++LFiwINH3rVq1CmMMNWvWvL5t504nzN99F559FnbscG5PS68U6iIikn6EhbGhTx8APvjqq+ubfXx8+Pbbb6lcuTIdOnTg4MGDN7115cqVVKpUiezZs2MtfP45VKkCR47AnDnw9deQPburTiR1KNRFRCT96NOHhWfPEpA1K3WaNLlhV8aMGZk9ezYAzz77LHFxcdf3RUVFERwcTK1atQgLg9atoXdvZ6n17duhRQuXnkWqUaiLiEj6MGsWsZMnsyhjRho0aYKPz813ZRcrVozPPvuM1atXM3r06Ovbf//9d65evUrhwrWpUsW5bj58OCxYcMvVWdMlhbqIiKR9x45B9+4ElynD6fBwmjVrdsumHTt2pH79+gwcOJBTp04BsGzZMowxvPXWE0RHw8qV8MYbaetpcCnBw05HREQ8jrXOU+OuXuX7xx7D19eXhg0b3rK5MYaRI0cSHh7OoEGDCAuL4IMPvsTaAtSuHcDmzc5y655IoS4iImnb6NGwZAn2ww+Zt3IltWvXJlu2bLd9S5kyZejRowfjx48nMDCI8PDjFC5ckkWLIFcuF9XtBgp1ERFJu3bvdsbJGzTgj1q12LdvHy1btkzSW+vXf5XY2DgiIv6kf/95HDnyq0vWNHcnhbqIiKRNERHQpo3zTPcJE5g9Zw7GGJon4Uby77+Htm2Lkzlzd4YP/5YPPrj1NXhPogVdREQkbXrtNdi2DRYuhPz5mTVrFo8//jj58uW75VushREjnLcGBcH8+V941Oz2O1FPXURE0p45c5xr6a+9Bo0asWvXLnbs2MEzzzxzy7dER8NLLzlLpbZsCb/+6lm3qyWFQl1ERNKWw4ehSxeoVg3eew+AmTNnYoyhVatWib4lLAyaNIEvv4QBA2DGDMiUyZVFpw0afhcRkbQjOhratYO4OGcxcz8/rLVMmzaNWrVqUaBAgZvecugQNG4Me/bAV185d7/drxTqIiKSdgweDMHBTqA/8AAA27ZtY/fu3fTt2/em5hs3OsukRkXBTz9B7dquLjht0fC7iIikDQsXwvvvQ9euzqz3eN9++y3e3t7861//uqH5L7/Ak09ChgzO/wfc74EOCnUREUkLDh6Ejh2dZdM+/fT65ri4OKZOnUqDBg3InTv39e0zZjg99GLFYM0aKFPGDTWnQQp1ERFxr4gIuDYBbtYsp+sdb8WKFYSGhtKxY8fr2z7/HNq2herVnWe4Fyzo6oLTLoW6iIi418svw+bNMGUKFC9+w65JkyaRLVs2mjZtirXwn/84S6Y+/bRzDT0gwE01p1GaKCciIu4zYQKMHw8DBzr3pCVw+fJlZs2aRfv27fHzy0T37jBunHO325dfQiIrr9731FMXERH32LwZevWCOnXgnXdu2j1z5kyuXLlC27adeOYZJ9AHDnS+KtATZ6y17q7hngQFBdkNGza4uwwREUmO8+ed57hGRcGmTZAnz01NatasyenTZwkM3MVvvxlGjnRG6u93xpiN1tqgxPbp/3VERMS14uKgUyc4etSZ6ZZIoO/atYs1a9ZQoMBw1q41fPutMzlObk+hLiIirvXBB/DDD86ta488kmiTYcNGAb6cPduJ77+Hhg1dW2J6pVAXERGXuTJ/Pt4DB3KiSRMi69Ujc0gIgYGB+CS4SL5o0VYmT/4SL6+qLFmSlyeecGPB6YxCXUREXOKXr7+m/vPPEwewYIHzAnx9falQoQL169cnf/4g+vbtBmRg4MBeCvRkUqiLiEiqu3T0KL27dcMHeOn55yn/+ONkzJiRixcvsn//ftatW8cHHwwnLi4WY7Izf/4Wnn66tLvLTncU6iIikqpiIiN5tFw5/oyO5qePPqLua6/d1GbtWqhf/wRRUR34/PP+CvS7pFAXEZFUNeKpp9hx+TIda9RINNCXLYOmTSEwMJBffvmFokXdUKSH0MNnREQk9UyaxA8rVpA/Uya+Wrnypt0LF/61MMuqVSjQ75FCXUREUsfatQS/+CIrgf5Dh+Ln53fD7hkzoHlzKF8eVqyA/PndUqVHUaiLiEjKCwmB5s1538+PXDlz8mL37jfsnjAB2rWDRx911kXPlctNdXoYhbqIiKSs8HBo3pwtly7xQ3g4ffv1I0uWLNd3f/qpsyhL3brw44+QLZsba/UwCnUREUk5cXHw/POwaRPvVqpEtmzZeDnBA9vfew/69oUWLWD+fMiUyY21eiCXhroxpoEx5k9jzD5jzIBbtGltjNlpjNlhjPnWlfWJiMg9evNNmDGDbf36MXvNGvr06UOOHDmw1tk1cCB07OhcT/f3d3exnsdlt7QZY7yBz4F6QAiw3hgz31q7M0GbUsCbQE1r7XljTF5X1SciIvfoyy9h+HDo2ZN3Dh8mW7ZsvPLKK8TFQZ8+8Pnn0KOH89VL48SpwpV/rNWBfdbaA9baKGA60OxvbboCn1trzwNYa0+5sD4REblbixY5a6M3bsymzp2ZM2cO/fr1I1u2nLzwghPkr78Oo0cr0FOTKx8+UxA4muDnEKDG39qUBjDGrAa8gSHW2h//fiBjTDegG0CRIkVSpVgREUmizZuhdWuoVAmmT+c/bdoQEBBA796v0rYtzJ4NQ4fC22+DMe4u1rOltf9f8gFKAU8C7YBxxpgcf29krR1rrQ2y1gblSWQdXhERcZEjR6BxY8iZExYsYOWmTSxatIhXXx3Ac89lZ/Zs+L//g//8R4HuCq7sqYcChRP8XCh+W0IhwDprbTRw0BizByfk17umRBERSbKwMCfQr1yB1auxgYH0b9GC/PkL8OOPvVmzBsaNgxdfdHeh9w9X9tTXA6WMMcWNMX5AW2D+39rMw+mlY4zJjTMcf8CFNYqISFJER8O//gW7d8OcOfDww8yaNYt169bh7/9f1q7NxNSpCnRXc1lP3VobY4zpDSzBuV4+wVq7wxgzFNhgrZ0fv6++MWYnEAu8Ya0966oaRUQkCayF7t1h6VL4+muoU4fIyEhef/3f+PuXJzS0E3PmOIu0iGu5dJU2a+0iYNHftg1K8L0FXo1/iYhIWvTuuzBxIgweDJ06AfDOOyM5cuQg/v4/sWiRN3XrurnG+5SWXhURkaQbNw4GDYLnnnNCHZg5cx3vv/8W3t7/ZNmyejz2mJtrvI+ltdnvIiKSVs2Z4zw9pmFDGD8ejGHevF20adMEsAwf3luB7mYKdRERubNly5xl1WrUIHb6dH5dvZqGDZ+lRYvygGXu3E28+moLd1d539Pwu4iI3N7GjdCsGQeLFmVs9epMLluWY8eOAX4Yk5kvvhhH8+YV3V2loFAXEZHb2bOHzXXr8r+4OObs24fXZ59RpUpDTp/+P0qUaMTSpZkpWFBPlUkrFOoiIpKoVXPm8Fzr1hyKjSV71qz8+9//pkiRXvTpU4iKFZ210HPndneVkpCuqYuIyE0uHDxI57ZtORQbS7cWLTh05AilSr1P796FeOQR+OUXBXpapFAXEZEbnDhwgCply3I4Opo5Q4YwZs4cJk3KQZcuULcuLFkC2bO7u0pJjEJdRET+EhFBh6AgDkZG8r+OHWk+aDDvvgv9+kHLljB/PmTK5O4i5VYU6iIi4oiKIrhOHZafP0+b6tXpP3kK/fs7K6w99xx89x34+7u7SLkdTZQTERGIiSGiTRteWLOGwjlzMuanpfTsCWPGwEsvwahR4KVuYJqnUBcRud/FxkLnzgyaN4/dwMIp0+jVKytTp8KAAfDee1oLPb1QqIuI3M/i4qBHD1ZPncrHxtDlha6MHVuf7793wvzNN91doCSHQl1E5H5lLfTrx+Xx4+kUEEDhbDnYs+cjVq1yhtt793Z3gZJcukIiInI/stYZWx81ilfLl+fAhQv4+X1NcHBWpkxRoKdXCnURkfvR0KEwfDjzGjRg3Pbt5MjxBiEh/+D776FjR3cXJ3dLw+8iIveb99+HIUMIad2azkuW4utbhbi4//Lzz1CzpruLk3uhnrqIyP3k3XfhrbeIadeOJn8eIywsihw5prNqlZ8C3QMo1EVE7hfvvHP9STIdY4qydetv5Mv3BevWlaJ8eXcXJylBw+8iIp7OWhgyxLmO3rkzLS8UYu68d8mRox1bt3YkXz53FygpJVmhbozxA7yttVdTqR4REUlJ1jq98//9D7p0odW58syd9wre3tlYvfoDBbqHSfLwuzGmL3Ac2GeM2WWM0Q0PIiJpmbXw1lvwv/8R0q49hX7cwZy5/ciduzYnToRSrlxhd1coKeyOoW6MGWmM6QT0BcpaawsC/wDKGWP+m9oFiojIXYi/Dz1i2DDerVaN4jPnERq6kZw5K7F9+3fkzp3F3RVKKkhKT305UALIDawxxmwCPgT2A22NMQGpWJ+IiCSXtfDGGywePpyHs2XjP+vXExPzFN267ebMmc0EBuZ2d4WSSu54Td1aOw+YZ4x5BHgFZwi+AlARyAksM8Zks9aWSM1CRUQkCazlRLdu9B0/nhlApuj8wCw+/LAer7/u7uIktSVnolwvYAawBdgOlAW2W2ufjJ9AJyIibhQbFcWIWrX43++/E+7tTb5cgzl7rj/ffONPhw7urk5cIcmhbq3da4ypAdTD6aVvA/rH74tKnfJERCQpYsLDqZg/PzsvXqRyngIc91rKlfCyLF4Mdeu6uzpxlWTd0hYf3gvjXyIikgacP3aMNhUrsvPiRWoXKcPGC3+QKZM3K1dCpUrurk5cSU+UExFJx+zFi1QvUYKfz5yhZ9WmrD6xiwIFvAkOVqDfjxTqIiLp1blz/K9cOfZFRPCPoo/w5abvqVYNVq+GYsXcXZy4g0JdRCQ9OnmSGZUr85/QUB4t9AQrD6+heXP4+WfImdPdxYm7KNRFRNKbQ4cIDgriuSNHKJr9IYJDltCrl2HmTMiY0d3FiTsp1EVE0pNt29hfvTpNQ0Px9yvI4bDlDBvmz6hR4O3t7uLE3bRKm4hIerFyJWeaNKFR+FXCTDawy/j22zy0a+fuwiStUKiLiKQHc+cS3rYtTYw3e2O9yZJ1AQsXlqZWLXcXJmmJht9FRNK6MWOIadWKFr6ZWBcZQd58U9mw4XEFutxEoS4iklZZC0OHEtujB62y5uenKxd44IHR/PFHK0qXdndxkhZp+F1EJC2KjYXevYn58kvK+QWw9+IxypYdyKZNPciQwd3FSVqlUBcRSWsiIqBjR47Pnk1t/4LsjQwlMDCIrVvfwdfX3cVJWqbhdxGRtOTsWahXj59nz6YovvwZGUrbtqM4fnw9vr66Z01uT6EuIpJW7NtHWPXqdFmzhvpAjClIjx6fMG1ab3dXJumEht9FRNKC4GAWNWhA10uXOGYNAQH/ZsWKIZQvrwvoknTqqYuIuFnYpEk8//jjNL54kTO2JNWqrWX//mEKdEk2hbqIiLtYy5IXX6R8585MiYsjI/3o1HU7q1dXIyDA3cVJeqThdxERN4iJiCCoUCG2nj1LIa+s+MQtZPioWvTqBca4uzpJrxTqIiIudnj7drrUru0EuslHWMY9fD87G0895e7KJL1TqIuIuNC5jRt5uFo1LltLLdoTWnwq6xZA2bLurkw8ga6pi4i4yLn586lfowYR1lKM5/GvO5X16xXoknIU6iIiLnD2//6POs2asT02jtyMpdVrE1i8GHLmdHdl4kk0/C4ikppiYjjVsyd1x49nN15k8pnBhxNb0bGjuwsTT+TSnroxpoEx5k9jzD5jzIDbtGtljLHGmCBX1icikqIuXOBY3bo8OX48u/EhINcilgYr0CX1uKynbozxBj4H6gEhwHpjzHxr7c6/tcsK9AXWuao2EZEUt28fB596inoHDnAEfx58eAlLlz5BvnzuLkw8mSt76tWBfdbaA9baKGA60CyRdv8FPgAiXFibiEjK+fFHdlapyuMHDnGYLDRuvoKNGxXokvpcGeoFgaMJfg6J33adMaYKUNhau/B2BzLGdDPGbDDGbDh9+nTKVyoicjeshffe47uGDXn00hVOkpMBb69m7twa+Pm5uzi5H6SZiXLGGC/gE6Dzndpaa8cCYwGCgoJs6lYmIpIEly5B585MnTOHZwGLDzNmrOaZZ0q7uzK5j7gy1EOBwgl+LhS/7ZqswMPAr8Z5RmIgMN8Y09Rau8FlVYqIJNeff3Lu6abU2buHLUCWLJWYM2c89eop0MW1XBnq64FSxpjiOGHeFmh/bae1NgzIfe1nY8yvwOsKdBFJ0374gZ/atKXr1QiOAHnzVubgwbVkyqTxdnE9l11Tt9bGAL2BJcAuYIa1docxZqgxpqmr6hARSRFxcVwYMIAXmzblqavhnDDFGDp0OSdPblKgi9u49Jq6tXYRsOhv2wbdou2TrqhJRCTZzp9nTv369N6wgZMY8uV6hZ+WvUuFChndXZnc5/SYWBGRZDj500/8q1AhWm3YQDSFqFtvPfsOfaxAlzRBoS4ikgQ2Lo6uNWtS+qmn+CE8nCKmOwP/7wA/LqlKlizurk7EkWZuaRMRSatOHTjAq7VrM/XIEbLjR4n8vzJ5/qME6UHWksYo1EVEbmP3/PlUbt6cCGupyiMUbvorEyf5kyOHuysTuZmG30VEbmHL0KHUbdaMWGspYbrw7Ihg5sxToEvapZ66iMjfhYezoGkz2v6ylEz4UyrvIib+8E+qV3d3YSK3p566iEhCf/zBiOIlaPbLUrITSNV6e/httwJd0geFuogIgLXEfPYZvStW4pVTJ8hjHqHPsH0sWlKEgAB3FyeSNBp+FxE5d45zzz5Hu0UL+QkokL0H3y/9jKAgb3dXJpIsCnURub+tXMny5q3ocf4M+/Ci1uNfsPjHbmTO7O7CRJJPw+8icn+KiYEhQxj3xBP88/wZ9uHPoCHLWLlKgS7pl3rqInL/OXiQA61a02jzBv4E/P0K8vWkKbRt+4S7KxO5Jwp1Ebl/WEvchAmMf6kXA6KiOA8UL16HHTsWkzGjr7urE7lnGn4XkfvD6dP8/uQ/efTFF+keFUmMXxBffrmWAweWKtDFY6inLiIe7/S33zLghReZGHmVjGTh8ZqjWLCwE9mzG3eXJpKi1FMXEY8VfuoUI2vWpHSHDkyOvEoBv+f4cnIoq37rrEAXj6Seuoh4pHnDhtHmzTeJAkpRlEqPfc+3syqSP7+7KxNJPQp1EfEoV86epWGlSqwKCcEfKOHVhb6fjuOllwxGnXPxcAp1EfEYJxct4umWrVgfGUFx8lCwwgbGzyjCgw+6uzIR19A1dRFJ/yIiWN2+PZUbN2ZbZCRlTVe6vX+K5RsV6HJ/UU9dRNI1GxzMJ02bM+DMKQLIRrnSi5g8uyYPP+zuykRcTz11EUmfrl7lUt++tH3sMV4/c4qcPErXtw6z7g8Futy/1FMXkfRn1Sq2dXiW1kcPswdD0XyDmLt4EJUrq58i9zf9CxCR9OP8eWzXbvz7H//gkaNHOEgAzz73C3uODFGgi6CeuoikB9bCzJns7/kSXc6dZQWQwb8kCxesoG5d3Xguco3+11ZE0rYjRwhv1Jg+bdpQ6txZVuBFxYptOHV6qwJd5G8U6iKSNsXGEvXxx3xRshQlf/yRUYCvbz7GjJnHli3TyZo1o7srFElzNPwuImlO9O+/83Xr1rx3+DCHgExeVXn91Y8YPvxJPRVO5DYU6iKSZkSeOME37dvzv+XLOQgE8ACPPTqKmbMaUqCA0lzkTjT8LiJuZ2NjGdSsGVny5+fF5cuJJS+FAqYxbtY+Vq9ppEAXSSKFuoi41S9jxpA/Y0b+O38+vnhRmLdo9coJdh5uS6tWCnOR5NDwu4i4x4ULfNe2LV2WLCEcKEcDslWfz5djfalY0d3FiaRP6qmLiGvFxXFl7Fiez5+ftkuWkJd85M/yG33HLGZ1sAJd5F6opy4irrNmDVu6daPtjh3sAfLRhZodv+Tjj33Im9fdxYmkfwp1EUl9R49i+/dn1PTpvAF4E0DhwjOZNLkOTz7p7uJEPIeG30Uk9Vy5AoMHc7p0aZp8N4O+AF5P8cpbe9i7T4EuktLUUxeRlBcXB9OmEde/Px8eO8bbeBGDN1WqfMK8eX0oXFiz2kVSg3rqIpKyli+HGjX4rmNHCh8/wwAg1vjwv//NYOPGvgp0kVSknrqIpIzt27n6xhvMWLKEkV5+bAawUdSr15+5c98hc+YM7q5QxOMp1EXk3oSEsL13b8Z9/z3fYDgPmLji1KrVm5Ejm1C5cjF3Vyhy31Coi8hdObh1K+936sSWbdtYby0+eGNoSsWKvZk0qTYVK2qYXcTVFOoikiwXjh1j4L/+xcTgYK4CRbyzERDblxzF+/Dxx7lp3hytpCbiJgp1EUmayEj+eOcdag8bxhlrKe+bjYvR/+ZitrcYPBh69gQ/P3cXKXJ/0+x3Ebm96Gjixo1jZGAgQe+/T6zxpqR3V3YTRqtX32L/fujbV4Eukhaopy4iiYuNhenTOfH22zx/6BA/AgV8a3As+nv+2SofP34AJUq4u0gRSUg9dRG5UVwczJoFFSuyuGNHKhwJ4Wf8gM8pVDmYVavyMWuWAl0kLVKoi4gjNha+/RbKlyfymWfodfgYjYDTcWUpUnwTc+a8xNq1hscfd3ehInIrCnWR+110NEycCGXLQocOzD13iXy+ORl9+TyZM7/M6NG/s2fPQ7RooVntImmdQl3kfhUZCWPGQOnS2BdeYE5EDCWylaXliaOERZ+nffvvOH36U3r2zICPZt+IpAv6pypyv7lyBb76ivPDhrH8+HEW5srDbP9chB09CEClSr0YNao7jz9e3s2FikhyuTTUjTENgJGANzDeWjvsb/tfBV4EYoDTwAvW2sOurFHEU+1YtYr5gwdzas0aVkZGshmwAGcjMKYejz/+OEOG1KFOnQpurlRE7pbLQt0Y4w18DtQDQoD1xpj51tqdCZptBoKsteHGmJ7AcKCNq2oU8UTBs2bxyb//zdwDB4gFMnh5UTBPJbzONQf+yQsvPMKgQb4UKuTmQkXknrmyp14d2GetPQBgjJkONAOuh7q1dnmC9muBji6sT8SzbNjArD59aBccTAxQL7AguUv3Ye7vfTl0zp9OneA//4FixdxdqIikFFeGekHgaIKfQ4Aat2nfBVic2A5jTDegG0CRIkVSqj6R9C8mBr7/niuffEK/NWsYDzyUKzc1nxjL1CUtuPobdOgAgwZByZLuLlZEUlqanChnjOkIBAFPJLbfWjsWGAsQFBRkXViaSNp0/jx89RV89hnbDx+mja8vuzE8VvUVtuwaxo45vrRpA4MHO3euiYhncmWohwKFE/xcKH7bDYwxdYGBwBPW2kgX1SaSPu3eDZ9+CpMmYcPDmVC6NL19/fDyCcCHqazdXId27eDNN+Ghh9xdrIikNleG+nqglDGmOE6YtwXaJ2xgjKkMjAEaWGtPubA2kfQjNhaWLHHCfMkS8PPj8jPP0Cn0AnN+XYgxdfGK+4ZOnfIxYACUKuXugkXEVVwW6tbaGGNMb2AJzi1tE6y1O4wxQ4EN1tr5wIdAFmCmcR5ddcRa29RVNYqkaSdPOkPsY8fC4cMQGAhDhzI2Qwlefectrlw5irf3O3TtOpABA7wpWtTdBYuIqxlr0/cl6aCgILthwwZ3lyGSOqyFX3+FL7+EOXMgJoaYJ5/kSKtWfHEghomT53L27ErAh1atfuTTT+tQoIC7ixaR1GSM2WitDUpsX5qcKCdyP7PWcm7fPrq0bEnUwYPkv3KFs76+nMyTh2MYQlauIu7XX+Nbe1OhwouMGNGb2rUrurNsEUkDFOoiaUVsLBfnzqVip04cDQ8nFvAxhnwBAQTkL0iszUfYsfzExRUhIKA4Tz3lT//+j1K5su5NExGHQl3E3fbuhYkT2fPVV7Q6dYpDwCP58vHFyJEUb9CGMWNg5Eg4dgzKl4c33oC2bcHX192Fi0hao1AXcYfLl2HmTJgwAX77je+N4Vlvb/yyZuWnadMoW7ExI0bA2K5w6RLUqeM0rV9fy5+KyK1p6VURV4mOhkWLnEe65csHL7xA3MmTDKlTh+bW8mClSkyd+QdTvmtM8eIwYgQ0aQIbN8LSpfDUUwp0Ebk99dRFUpO1sGYNfPstzJgBZ85AQAB07MjlZ57hudGjmTt3LrVrdyYm5gsaNMhA5szQqxf066fnsotI8ijURVLDjh0wdSpMmwaHDkHGjNC0qdNLf+opDh8/TuPGTdm58w9y5BjB8uV9KFrUMHw4dOkCOXO6+wREJD1SqIukBGuZ+fHHVDp0iFLLl8POneDtDfXqwdCh0Lw5ZM0KwEcfzeLtt3sSGRkNLKZy5fr06QNPP+28RUTkbinURe6WtZxetozmnTpx8uRJ9sfEAJDfz48SJUpQtFIl8hcvTp7jx8k8+RuWLNnHr7+u5tKldUBWWrf+nbffLkP58u49DRHxHAp1keSIjYXgYJg/n0uzZtHx4EHWAMX9/en26KPkqlyZ4xcvcuDAAVZv3MjxBQuIjPxrXSJjMlOpUns+++xNatYs477zEBGPpFAXuZNLl+Cnn2D+fFi4EM6e5YSPD40zZWKrlxcj/vtf+r711vXm1sLq1TBmDMyYYYErPPbYRZ54Yj3du1elaNFC7jsXEfFoCnWRxBw9CgsWOEG+bBlERTmz1hs3Zl/16tT/+GNOnj7N/B9+oFGjRoCzpPmUKU6Y79wJ2bJB166G7t2zUL58FqCZe89JRDyeQl0EICICfvsNfvzRee3Y4WwvWRJeftmZuf7YY2zdsYP69esTGxvL8uXLqVatOsHBTpB/951zmOrVncXU2rSBzJnde1oicn9RqMt9KS42lsVffUXdy5fx/+UXWL4crl4FPz/4xz+gc2fnyS8PPnj9iS/BwcE0bNiQrFmzMmfOz/z+exlefBG2b4csWZy3dO8OlSq588xE5H6mUJf7R0gIwWPH8uyIEVy6fJlT1vIw8HZgIM2ffx7/xo3hiScS7V4vX76cp59+moCA/Dz22FKefLIo4eFQtaqzvHm7dk6wi4i4k0JdPNfJk85a5MuWwfLlnN+7lxeA/cAjAQHUKlqU9adO0fbYMXJ99x3tvLxolyMHjzzyCF5efz1BeezYmfTq9RxeXiUICVnKwoWBtG/v9MqDEl3RWETEPRTq4hmshX37nGnnv/3mfN2929mXLRsXH3uMBpGRHDh+nAVz5tC4SRMA4uLi+Omnn5gwYQLjxo3js88+I3/+/NSv34ALF/KycmUw58+vBIpSo8ZyunbNwzPPOJPgRETSGmOtdXcN9yQoKMhu2LDB3WWIq0VFwebNN4b4qVPOvoAAqFkTHn8catfmatmyNHz6aX777Tdmz55Ns2aJz0K/cCGMzz5bwKRJ37Nv32LgMpCZAgXKM2rUp7RsWc1lpycicivGmI3W2kTHCdVTlzTv6pUrHFi+nCljxvB+sWKYDRucQL/2UJcHHoAGDf4K8jJlIH74PDo6mtYtW7Jy5UqmTp2aaKCfOuU8pn3ixOxs394Bf/8OPPNMNLlzf8Hbb7ehQIF8rjxdEZG7plCXtMVaZwGUzZth/Xp6TJnCuNBQ4uJ3h3p7M7ZGDTL27g01ajghnj//LQ5l6d69OwsWLGD06NG0a9fu+r7wcOcW9KlTnTvYYmKgWjUYPRratoWAAF+gT6qfrohISlKoi/tERjpPadmy5a/X1q0QFgbAFm9vvo6LI0/GjAxo2pTN4eFM/uEHtly8yDfPPkvFihVve/hBgwYxceJEBg8eTM+ePYmJcdYlnzoV5s6FK1egYEFnidPOneGhh1L5fEVEUplCXVJfVBTs3es80GXnzr++7tnjdJEBMmWCihWhfXuoVImzxYrRvGtXcsfGsnHjRvLlc4bA2y5ezPPPP0+1atUYPHgw/fv3x9fX96aPHDduHO+++y4vvNCFBg0G06eP83CYU6cgRw7nFrQOHaBWLa2MJiKeQxPlJGVYC2fO8MPEiWwODmbnn39y8MQJ1ubLh9m376/wNgZKlHC6xeXKOU9qqVTJ2RafrrGxsTRq1IgVK1awatUqqlW7cYLamTNn6NWrFzNmzKBChQp88cUXPPbYY9f3//TTTzRq1IjixesTFzefAwd88Pd3niXToQM0agT+/i76cxERSWGaKCcpIyYGQkPh8GHntW+f0wOPf8WEhfE8cDbBWz7Kn583+vd3Avyhh5wntGXMeNuPeeedd/jpp58YP378TYEOkDt3br777jvatWvHyy+/TM2aNWnfvj2tWr3Mp5/OYcWKT4Ay7Nv3HXXq+PD229CyJWTPnqJ/GiIiaY566vc5ay3GGKKjo7l44gS5rlz5K7T//goNhbi4v97s5QVFi0KpUlCyJMMOH+bNhQv5ZMAAer75Jp26dmXGjBlMmjSJ5557Lkn1LFmyhIYNG9K5c2cmTJhwx/a7dl3h5ZffZ9my4VgbDYAxGeja9QcGD65LgQJ39cciIpJm3a6nrlD3dDExcPo0nDiR6KvE/PlciIwkzFryAaEJ3+vtDYUKOcF9q1f8OPaff/5JxYoVady4MbNnzwYgMjKSJk2asHz5cmbOnEmLFi1uW+qxY8eoWLEigYGBrFu3jkyZMiXa7vBhmDULZsyA3393tj344FoiI//NmDFfUrfugzc8EU5ExJMo1NO7mBi4fJm4sDBmzZ5NgYwZWfP771TJk4e6BQvCuXOJv86fd74m9jvOlo31OXJQ/cgRHsyalYDMmVl74gRz+valRatWTmAXKAA+d75CExcXxxNPPMGOHTvYuXMngYGB1/ddvnyZevXqsXHjRubOnUvjxo0TPUZsbCz16tVj3bp1bNy4kTJlytyw/1qQz5wJ69Y526pWhWeecV4PPJD0P04RkfRM19ST6PLJk5w7cIC4mBh2791LnRo18PXygthYZ9g5NjZZ3586dYq3vvmGB3PnZuPBg/SvVYsqefI463NGRDi3dF37/soVuHQJLl92Xte+v3Tp+kNWTgJtEtRbAdgKzuSzHDkgZ07naWo5czoplzMn5Mrl3McdGPjXK18+bMaM9P/nP8lz9Srr9+8nQ4YMVKtWjZe++44nBg0iZ86cSf5zGzt2LL/99hsTJ068IdABsmTJwuLFi6lbty4tW7Zk9uzZNIl/RGtCH3zwAcuXL2fChAnXA333bufWs3nz/uqRV6kCw4bBv/7lzK0TEZEErLXp+lW1alWbUsZ36mSB669lTh/3rl8jExwLsJ2u7TPG2owZrQ0IsDYw0NqiRa196CFrH3nE2rp1rW3e3Npnn7W2Z09r33jD2qFDrf3kEzvk6actYL/o3du2rlPHent7273r11sbE5Psc128eLEF7KhRo65v27x5s/Xx8bEdOnRI8nGOHTtms2fPbuvUqWPj4uJu2e7cuXM2KCjI+vr62hkzZtywb926ddbb29u2bt3GBgfH2QEDrH3wwb/+KKtXt/b9963dty/Zpyki4nGADfYWmajh9wT2/Pgjv333HWevXKH/zJkMbtGCIa1bO9eWvbycr8n4vkX//qzYuJE1c+fS8ZVXuHTlCrt37MD4+l5fozuprl69StGiRalevToLFizg+PHjPPDAA7Rp04avv/46WceKi4ujatWqhIWFsXv3bvz8/K7ve+eddxgyZAizZs2iVatWdzxWu3btmDt3Ltu3b6dUqVK3bRsWFkbjxo0JDg5m9OjRdO/endOnw3j44QqEh0OWLFs5cSIHPj7w5JPQogU0bepc1hcREcftht/d3tO+11dK9tQTevTRR23FihXv+v2nT5+2vr6+9pVXXrHWWvv11187vf9ly+7qeOPGjbOAXb58+fVtr776qvXy8rK7d+9O1rGmTZtmAfvNN9/ctC8qKspWrVrV5sqVyx47duy2x/n5558tYIcMGZLkz758+bKtX7+RBWyOHBUs+FvA+vktsa1aWTtlirXnziXrdERE7ivcpqfu9lC+11dqhfqIESMsYHft2nVX7x85cqQF7NatW6211oaHh9uAgADbunXrZB8rLi7OlitXzlauXPmGIe6TJ0/azJkz27Zt2yb5WNHR0bZ06dL24YcftrGxsYm22bVrl82YMaOtX7/+LdtERETY0qVL25IlS9qrV6/e8XP37bN2xAhr69Wz1scnykIVC8YWKtTYNmr0kr1wITLJ5yAicj9TqN+F0NBQa4xJVi80ocqVK9sqVarcsO2VV16xPj4+9vjx48k61pIlSyxgJ0+efNO+t956ywJ2y5YtSTrWxIkTLWDnzp1723ajR4+2gP3oo48S3T9s2DAL2MWLFye6PyrK2uXLrX3ttRuvj5cta+3rr1u7cOFFe+hQSJJqFhGRvyjU79KTTz5pH3zwwdtOAEvM1q1bLWA//fTTG7bv2bPHAnbo0KHJOl7Dhg1tYGCgjYiIuGnfuXPnbI4cOWzjxo3veJyoqChbvHhxW7Vq1TueU1xcnG3evLn19fW169atu2FfSEiIzZw5s23atOkN20+dsnbSJGtbt7Y2e3bnb5efn7X161s7cqQmuomIpASF+l0aM2aMBezmzZuT9b5XX33V+vr62tOnT9+0r379+rZgwYI2KioqScfavXu3Bew777xzyzbXes0rVqy47bGuXZdfsGBBkj773LlztmjRorZYsWL27Nmz17d36NDB+vv72z179tvgYGvfeceZuG+M8zcqMNDaLl2snTPH2osXk/RRIiKSRAr1u3TmzBnr4+Nj33jjjSS/Jzo62ubLl8+2aNEi0f3z58+3wE23dd1K7969rZ+fnz1x4sQt21y5csUWLFjQ1qhR45Y98KioKFusWDFbvXr1ZI08rFu3zvr6+tpGjRrZ2NhYO2nSQgvYMmXesgEB9vodekFB1g4ZYu369dbe4jK8iIikAIX6PWjSpIktXLjwLSeM/d3ChQtve806JibGFitWzNaqVeuOxwoLC7NZsmSxzz777B3bfvXVV7f9n4Xx48dbwC5cuPCOx0ooPNza3r2d6+ve3vnj77n3sfnyXbKdO1s7bZq1iQxIiIhIKlGo34Nvv/02SUPb17Rt29bmzJnTRkbeejb3Rx99lKRh/U8//dQC9vfff7/j58bExNjy5cvbBx544KZr79HR0faBBx6wQUFBSbiWbu327dZ+/LFzLTxDBmshzkIJ6+OT09as2ce+/fYXNpnTDEREJIUo1O/B5cuXbebMmW23bt3u2DYsLMxmyJDBvvTSS7dtd+7cOZspUyb7/PPP37JNXFycLVOmjK1evXqSa702S/7DDz+8YfvkyZMtYOfNm5fo+/bvt3bcOGvbtbM2Xz57w0z1fv2sXbzY2nPnIpI8WiEiIqlHoX6POnToYAMCAhKdfZ7QhAkTLGCDg4PveMyePXtaf39/e/LkyUT3//LLLxawkyZNSlatjRo1stmyZbt+3NjYWFu2bFlboUKF66EcGmrtN99Y+8IL1hYr9leI589vbYcO1o4fb+3hw8n6WBERcZHbhbrWp0yCjh07cv78eRYvXnzbdt988w0lS5akRo0adzxm3759iYyM5Isvvkh0/+jRo8mVKxetW7dOVq2ffPIJ4eHhDBw4EIB58+axa9cu6tR5kz59vChbFgoWhI4dncVSqlSBzz6DnTud5dK/+Qa6dIEiRZL1sSIikhbcKu3Ty8sVPfXo6GibN29e27Jly1u2CQkJscYYO3jw4CQft1GjRjZv3rw3PZEtNDTUent729dff/2u6u3W7RULxpYo0cYak9FCgIUYmzmztQ0bWvvhh9Zu3HhX68CIiIiboZ76vfHx8aFdu3YsWLCA8+fPJ9pm+vTpWGvp0KFDko/76quvcurUKaZOnXrD9vHjxxMbG0v37t3veAxr4eBBmDTJ6WGXKgVjxw4GDPv3z8DfPw9BQc+werU358/DokXw+utOD93bO8mliohIOqBV2pJo48aNBAUFMWbMGLp163bT/qpVq+Ll5cX69euTfExrLZUrVyYqKoo//vgDLy8vYmJiKF68OOXKlWPJkiU3vScuDnbtgt9+g5UrnVdIiLMvZ06oVQv+8Q/w8lrE888/Qfbsme/6nEVEJO253SptPq4uJr2qUqUK5cqVY/LkyTeF+u7du9m0aRP/93//l6xjGmN4/fXXefbZZ1m0aBFNmjRh8eLFhISEMHLkSAAuXoR162DNGggOhrVrISzMeX/+/E6AX3uVK+es/upodI9nLCIi6Y166snwwQcfMGDAAPbt20eJEiWubx88eDD//e9/CQkJoUCBAsk6ZnR0NCVKlKBYsWL8+usKatasy86df9C6dQjr1vnyxx/OELsx8PDD8Nhj8OijULMmlCiR7GXZRUQknbtdT13X1JOhY8eOeHl5MWnSpOvbrLVMmzaN2rVrJzvQL1yAFSt8KVWqC6tWrcLbO4C1a5dx8WJlZs70pUABGDIEfvrJabttG3z5JXTqBCVLKtBFRORGGn5PhoIFC1KvXj0mT57MkCFD8PLyYtOmTezdu5f+/fvf9r1XrsDmzbB+PWzY4Hzdu/fa3h7Ae+TN+zAlSpSie/cXePbZhEPpIiIid6ZQT6bOnTvTrl07li9fTp06dZg+fTq+vr60bNnyepuoKNi+3Qnua68dO5xJbgCFCkG1avD8887XqlXzkSNHBEZdbxERuQcuDXVjTANgJOANjLfWDvvbfn9gMlAVOAu0sdYecmWNd9KsWTOyZ8/OxIkTqV27NtOmTScoqD7TpuVkyxbYutUZJo+MdNrnyuUEd/Pmztdq1SAwMLEjK9BFROTeuCzUjTHewOdAPSAEWG+MmW+t3ZmgWRfgvLW2pDGmLfAB0MZVNd6OtXD0KGzZkpFCherz7bfTmT59NbGxIYSGdic42AnwSpWgT5+/ArxoUV37FhER13BlT706sM9aewDAGDMdaAYkDPVmwJD472cBnxljjHXxFP2oKOexqdd63te+/vXcmbrAbLJkyUGWLNV5882mNG8OBQoowEVExH1cGeoFgaMJfg4B/v6Q9OttrLUxxpgwIBdwxhUFfv89DBrkPNwlOtrZlikTlC8PrVs7vfCKFeGhh14kY8bn8fX1dUVZIiIiSZIuJ8oZY7oB3QCKpODKI1myOIudNG78V4CXLJnY41S90N2AIiKS1rgy1EOBwgl+LhS/LbE2IcYYHyA7zoS5G1hrxwJjwXn4TEoVWKeO8xIREUmPXNndXA+UMsYUN8b4AW2B+X9rMx/oFP/9v4Blrr6eLiIikl65rKcef428N7AE55a2CdbaHcaYoTjLyM0HvgKmGGP2Aedwgl9ERESSwKXX1K21i4BFf9s2KMH3EcAzrqxJRETEU2i2l4iIiIdQqIuIiHgIhbqIiIiHUKiLiIh4CIW6iIiIh1Coi4iIeAiFuoiIiIdQqIuIiHgIhbqIiIiHUKiLiIh4CJPe10sxxpwGDru7jnuQGxetF+9iOq/0ReeVfnjiOYHOKzmKWmvzJLYj3Yd6emeM2WCtDXJ3HSlN55W+6LzSD088J9B5pRQNv4uIiHgIhbqIiIiHUKi731h3F5BKdF7pi84r/fDEcwKdV4rQNXUREREPoZ66iIiIh1Cou4gxpoEx5k9jzD5jzIBE9vcwxmw3xmwxxvxmjCnnjjqT607nlaBdK2OMNcaki9mtSfh9dTbGnI7/fW0xxrzojjqTKym/L2NMa2PMTmPMDmPMt66uMbmS8Lv6vwS/pz3GmAtuKDPZknBeRYwxy40xm40x24wxjdxRZ3Il4byKGmN+iT+nX40xhdxRZ3IYYyYYY04ZY/64xX5jjPk0/py3GWOqpFox1lq9UvkFeAP7gQcAP2ArUO5vbbIl+L4p8KO7606J84pvlxVYCawFgtxddwr9vjoDn7m71lQ4r1LAZiAg/ue87q77Xs/pb+1fBia4u+4U+l2NBXrGf18OOOTuulPovGYCneK//ycwxd11J+G8/gFUAf64xf5GwGLAAI8A61KrFvXUXaM6sM9ae8BaGwVMB5olbGCtvZjgx8xAepjscMfzivdf4AMgwpXF3YOknld6k5Tz6gp8bq09D2CtPeXiGpMrub+rdsA0l1R2b5JyXhbIFv99duCYC+u7W0k5r3LAsvjvlyeyP82x1q4Ezt2mSTNgsnWsBXIYY/KnRi0KddcoCBxN8HNI/LYbGGN6GWP2A8OBPi6q7V7c8bzih5kKW2sXurKwe5Sk3xfQKn4obZYxprBrSrsnSTmv0kBpY8xqY8xaY0wDl1V3d5L6u8IYUxQozl+BkZYl5byGAB2NMSHAIpxRiLQuKee1FWgZ/30LIKsxJpcLaktNSf57eq8U6mmItfZza20J4N/A2+6u514ZY7yAT4DX3F1LKvgBKGatrQD8DExycz0pxQdnCP5JnF7tOGNMDncWlILaArOstbHuLiSFtAO+ttYWwhnenRL/by69ex14whizGXgCCAU85XeW6jzhL0B6EAok7MkVit92K9OB5qlZUAq503llBR4GfjXGHMK5ljQ/HUyWu+Pvy1p71lobGf/jeKCqi2q7F0n5exgCzLfWRltrDwJ7cEI+rUrOv622pI+hd0jaeXUBZgBYa4OBDDjPGU/LkvJv65i1tqW1tjIwMH7bBZdVmDqSmwF3TaHuGuuBUsaY4sYYP5z/uMxP2MAYk/A/nI2BvS6s727d9rystWHW2tzW2mLW2mI4E+WaWms3uKfcJEvK7yvh9bCmwC4X1ne37nhewDycXjrGmNw4w/EHXFhjciXlnDDGlAECgGAX13e3knJeR4A6AMaYsjihftqlVSZfUv5t5U4w4vAmMMHFNaaG+cBz8bPgHwHCrLXHU+ODfFLjoHIja22MMaY3sARn9ucEa+0OY8xQYIO1dj7Q2xhTF4gGzgOd3Fdx0iTxvNKdJJ5XH2NMUyAGZ4JMZ7cVnERJPK8lQH1jzE6cIc83rLVn3Vf17SXj72BbYLqNn4qc1iXxvF7DuTzyCs6kuc5p/fySeF5PAu8bYyzOXTO93FZwEhljpuHUnTt+jsNgwBfAWvslzpyHRsA+IBx4PtVqSeN/B0RERCSJNPwuIiLiIRTqIiIiHkKhLiIi4iEU6iIiIh5CoS4iIuIhFOoiIiIeQqEuIiLiIRTqInJPjDHe7q5BRBx6opyIJJsxZibOk/QqAguAd91bkYiAQl1E7k55YIa19hF3FyIif9FjYkUkWYwxGXAWEylgrY1xdz0i8hddUxeR5HoIWKdAF0l7FOoiklzlgW3uLkJEbqZQF5HkUqiLpFG6pi4iIuIh1FMXERHxEAp1ERERD6FQFxER8RAKdREREQ+hUBcREfEQCnUREREPoVAXERHxEAp1ERERD/H/OIT0EsMt4UQAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# plot homogenised potential \n", + "fig, ax = plt.subplots(1, 1, figsize=(8,6))\n", + "\n", + "ax.plot(r_total_mesh, phi_n(x=r_total_mesh), 'b', label=r\"$\\phi^-$\")\n", + "ax.plot(r_total_mesh, phi_p(x=r_total_mesh), 'r', label=r\"$\\phi^+$\")\n", + "for i in range(len(tt)):\n", + " ax.plot(r_mesh_pos[i,:], phi_p(x=r_mesh_pos[i,:]), 'k', label=r\"$\\phi$\" if i ==0 else \"\")\n", + "for i in range(len(tt)-1):\n", + " ax.plot(r_mesh_neg[i,:], phi_n(x=r_mesh_neg[i,:]), 'k')\n", + " ax.plot(r_mesh_am1[i,:], phi_am1(r_mesh_am1[i,:], tt[i]), 'k')\n", + " ax.plot(r_mesh_am2[i,:], phi_am2(r_mesh_am2[i,:], tt[i]), 'k')\n", + "ax.set_xlabel(r\"$r$\")\n", + "ax.set_ylabel(r\"$\\phi$\")\n", + "ax.legend();" + ] + } + ], + "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.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/pybamm/geometry/standard_spatial_vars.py b/pybamm/geometry/standard_spatial_vars.py index 9c2659f437..e70d773a16 100644 --- a/pybamm/geometry/standard_spatial_vars.py +++ b/pybamm/geometry/standard_spatial_vars.py @@ -30,8 +30,8 @@ y = pybamm.SpatialVariable("y", domain="current collector", coord_sys="cartesian") z = pybamm.SpatialVariable("z", domain="current collector", coord_sys="cartesian") -r = pybamm.SpatialVariable( - "r", domain="current collector", coord_sys="cylindrical polar" +r_macro = pybamm.SpatialVariable( + "r_macro", domain="current collector", coord_sys="cylindrical polar" ) r_n = pybamm.SpatialVariable( @@ -104,8 +104,8 @@ z_edge = pybamm.SpatialVariableEdge( "z", domain="current collector", coord_sys="cartesian" ) -r_edge = pybamm.SpatialVariableEdge( - "r", domain="current collector", coord_sys="cylindrical polar" +r_macro_edge = pybamm.SpatialVariableEdge( + "r_macro", domain="current collector", coord_sys="cylindrical polar" ) diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py index 9411e924ed..5b937a75cc 100644 --- a/pybamm/spatial_methods/finite_volume.py +++ b/pybamm/spatial_methods/finite_volume.py @@ -168,23 +168,16 @@ def divergence(self, symbol, discretised_symbol, boundary_conditions): divergence_matrix = self.divergence_matrix(symbol.domains) - # check for particle domain - if submesh.coord_sys == "spherical polar": + # check coordinate system + if submesh.coord_sys in ["cylindrical polar", "spherical polar"]: second_dim_repeats = self._get_auxiliary_domain_repeats(symbol.domains) - - # create np.array of repeated submesh.edges - r_edges_numpy = np.kron(np.ones(second_dim_repeats), submesh.edges) - r_edges = pybamm.Vector(r_edges_numpy) - - out = divergence_matrix @ ((r_edges ** 2) * discretised_symbol) - elif submesh.coord_sys == "cylindrical polar": - second_dim_repeats = self._get_auxiliary_domain_repeats(symbol.domains) - # create np.array of repeated submesh.edges r_edges_numpy = np.kron(np.ones(second_dim_repeats), submesh.edges) r_edges = pybamm.Vector(r_edges_numpy) - - out = divergence_matrix @ (r_edges * discretised_symbol) + if submesh.coord_sys == "spherical polar": + out = divergence_matrix @ ((r_edges ** 2) * discretised_symbol) + elif submesh.coord_sys == "cylindrical polar": + out = divergence_matrix @ (r_edges * discretised_symbol) else: out = divergence_matrix @ discretised_symbol @@ -207,14 +200,15 @@ def divergence_matrix(self, domains): """ # Create appropriate submesh by combining submeshes in domain submesh = self.mesh.combine_submeshes(*domains["primary"]) - if submesh.coord_sys == "spherical polar": - r_edges_left = submesh.edges[:-1] - r_edges_right = submesh.edges[1:] - d_edges = (r_edges_right ** 3 - r_edges_left ** 3) / 3 - elif submesh.coord_sys == "cylindrical polar": + + # check coordinate system + if submesh.coord_sys in ["cylindrical polar", "spherical polar"]: r_edges_left = submesh.edges[:-1] r_edges_right = submesh.edges[1:] - d_edges = (r_edges_right ** 2 - r_edges_left ** 2) / 2 + if submesh.coord_sys == "spherical polar": + d_edges = (r_edges_right ** 3 - r_edges_left ** 3) / 3 + elif submesh.coord_sys == "cylindrical polar": + d_edges = (r_edges_right ** 2 - r_edges_left ** 2) / 2 else: d_edges = submesh.d_edges @@ -287,14 +281,15 @@ def definite_integral_matrix( domain = child.domains[integration_dimension] submesh = self.mesh.combine_submeshes(*domain) - if submesh.coord_sys == "spherical polar": - r_edges_left = submesh.edges[:-1] - r_edges_right = submesh.edges[1:] - d_edges = 4 * np.pi * (r_edges_right ** 3 - r_edges_left ** 3) / 3 - elif submesh.coord_sys == "cylindrical polar": + + # check coordinate system + if submesh.coord_sys in ["cylindrical polar", "spherical polar"]: r_edges_left = submesh.edges[:-1] r_edges_right = submesh.edges[1:] - d_edges = 2 * np.pi * (r_edges_right ** 2 - r_edges_left ** 2) / 2 + if submesh.coord_sys == "spherical polar": + d_edges = 4 * np.pi * (r_edges_right ** 3 - r_edges_left ** 3) / 3 + elif submesh.coord_sys == "cylindrical polar": + d_edges = 2 * np.pi * (r_edges_right ** 2 - r_edges_left ** 2) / 2 else: d_edges = submesh.d_edges