diff --git a/docs/conf.py b/docs/conf.py
index 8e86dcc48d..55692309dc 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -76,6 +76,7 @@
doctest_global_setup = """
from docs import *
+import pybamm
"""
# Add any paths that contain templates here, relative to this directory.
diff --git a/docs/source/api/experiment/experiment_steps.rst b/docs/source/api/experiment/experiment_steps.rst
index 55de9d17bf..6a2e2abc31 100644
--- a/docs/source/api/experiment/experiment_steps.rst
+++ b/docs/source/api/experiment/experiment_steps.rst
@@ -18,3 +18,28 @@ directly:
.. autoclass:: pybamm.step._Step
:members:
+
+Step terminations
+-----------------
+
+Standard step termination events are implemented by the following classes, which are
+called when the termination is specified by a specific string. These classes can be
+either be called directly or via the string format specified in the class docstring
+
+.. autoclass:: pybamm.step.CrateTermination
+ :members:
+
+.. autoclass:: pybamm.step.CurrentTermination
+ :members:
+
+.. autoclass:: pybamm.step.VoltageTermination
+ :members:
+
+The following classes can be used to define custom terminations for an experiment
+step:
+
+.. autoclass:: pybamm.step.BaseTermination
+ :members:
+
+.. autoclass:: pybamm.step.CustomTermination
+ :members:
diff --git a/docs/source/api/expression_tree/unary_operator.rst b/docs/source/api/expression_tree/unary_operator.rst
index ad5bb0a48f..e6a3cbe554 100644
--- a/docs/source/api/expression_tree/unary_operator.rst
+++ b/docs/source/api/expression_tree/unary_operator.rst
@@ -34,6 +34,9 @@ Unary Operators
.. autoclass:: pybamm.Mass
:members:
+.. autoclass:: pybamm.BoundaryMass
+ :members:
+
.. autoclass:: pybamm.Integral
:members:
@@ -58,6 +61,9 @@ Unary Operators
.. autoclass:: pybamm.BoundaryGradient
:members:
+.. autoclass:: pybamm.EvaluateAt
+ :members:
+
.. autoclass:: pybamm.UpwindDownwind
:members:
diff --git a/docs/source/api/plotting/quick_plot.rst b/docs/source/api/plotting/quick_plot.rst
index ff7576a00d..870a569e9d 100644
--- a/docs/source/api/plotting/quick_plot.rst
+++ b/docs/source/api/plotting/quick_plot.rst
@@ -7,3 +7,6 @@ Quick Plot
:members:
.. autofunction:: pybamm.dynamic_plot
+
+.. autoclass:: pybamm.QuickPlotAxes
+ :members:
diff --git a/docs/source/examples/index.rst b/docs/source/examples/index.rst
index 4afaa6eeeb..e0f2bd5832 100644
--- a/docs/source/examples/index.rst
+++ b/docs/source/examples/index.rst
@@ -65,6 +65,7 @@ The notebooks are organised into subfolders, and can be viewed in the galleries
notebooks/models/rate-capability.ipynb
notebooks/models/saving_models.ipynb
notebooks/models/SEI-on-cracks.ipynb
+ notebooks/models/simulate-3E-cell.ipynb
notebooks/models/simulating-ORegan-2022-parameter-set.ipynb
notebooks/models/SPM.ipynb
notebooks/models/SPMe.ipynb
@@ -84,6 +85,17 @@ The notebooks are organised into subfolders, and can be viewed in the galleries
notebooks/parameterization/parameter-values.ipynb
notebooks/parameterization/parameterization.ipynb
+.. nbgallery::
+ :caption: Simulations and Experiments
+ :glob:
+
+ notebooks/simulations_and_experiments/callbacks.ipynb
+ notebooks/simulations_and_experiments/custom-experiments.ipynb
+ notebooks/simulations_and_experiments/experiments-start-time.ipynb
+ notebooks/simulations_and_experiments/rpt-experiment.ipynb
+ notebooks/simulations_and_experiments/simulating-long-experiments.ipynb
+ notebooks/simulations_and_experiments/simulation-class.ipynb
+
.. nbgallery::
:caption: Plotting
:glob:
@@ -110,11 +122,6 @@ The notebooks are organised into subfolders, and can be viewed in the galleries
:glob:
notebooks/batch_study.ipynb
- notebooks/callbacks.ipynb
notebooks/change-settings.ipynb
notebooks/initialize-model-with-solution.ipynb
- notebooks/rpt-experiment.ipynb
- notebooks/simulating-long-experiments.ipynb
- notebooks/simulation-class.ipynb
notebooks/solution-data-and-processed-variables.ipynb
- notebooks/experiments-start-time.ipynb
diff --git a/docs/source/examples/notebooks/batch_study.ipynb b/docs/source/examples/notebooks/batch_study.ipynb
index 807a368fcc..f02d1154ad 100644
--- a/docs/source/examples/notebooks/batch_study.ipynb
+++ b/docs/source/examples/notebooks/batch_study.ipynb
@@ -523,7 +523,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "The difference in the individual plots is not very well visible in the above slider plot, but we can access all the simulations created by `BatchStudy` (`batch_study.sims`) and pass it to `pybamm.plot_summary_variables` to plot the summary variables (more details on \"summary variables\" are available in the [`simulating-long-experiments`](https://github.com/pybamm-team/PyBaMM/blob/develop/docs/source/examples/notebooks/simulating-long-experiments.ipynb) notebook).\n",
+ "The difference in the individual plots is not very well visible in the above slider plot, but we can access all the simulations created by `BatchStudy` (`batch_study.sims`) and pass it to `pybamm.plot_summary_variables` to plot the summary variables (more details on \"summary variables\" are available in the [`simulating-long-experiments`](./simulations_and_experiments/simulating-long-experiments.ipynb) notebook).\n",
"\n",
"## Comparing summary variables"
]
diff --git a/docs/source/examples/notebooks/change-settings.ipynb b/docs/source/examples/notebooks/change-settings.ipynb
index 5b21f4dd6b..c54da8754c 100644
--- a/docs/source/examples/notebooks/change-settings.ipynb
+++ b/docs/source/examples/notebooks/change-settings.ipynb
@@ -7,7 +7,7 @@
"source": [
"# Changing settings when solving PyBaMM models\n",
"\n",
- "[This](https://github.com/pybamm-team/PyBaMM/blob/develop/docs/source/examples/notebooks/models/SPM.ipynb) example notebook showed how to run an SPM battery model, using the default parameters, discretisation and solvers that were defined for that particular model. Naturally we would like the ability to alter these options on a case by case basis, and this notebook gives an example of how to do this, again using the SPM model. In this notebook we explicitly handle all the stages of setting up, processing and solving the model in order to explain them in detail. However, it is often simpler in practice to use the `Simulation` class, which handles many of the stages automatically, as shown [here](https://github.com/pybamm-team/PyBaMM/blob/develop/docs/source/examples/notebooks/simulation-class.ipynb).\n",
+ "[This](./models/SPM.ipynb) example notebook showed how to run an SPM battery model, using the default parameters, discretisation and solvers that were defined for that particular model. Naturally we would like the ability to alter these options on a case by case basis, and this notebook gives an example of how to do this, again using the SPM model. In this notebook we explicitly handle all the stages of setting up, processing and solving the model in order to explain them in detail. However, it is often simpler in practice to use the `Simulation` class, which handles many of the stages automatically, as shown [here](./simulations_and_experiments/simulation-class.ipynb).\n",
"\n",
"\n",
"### Table of Contents\n",
diff --git a/docs/source/examples/notebooks/getting_started/tutorial-1-how-to-run-a-model.ipynb b/docs/source/examples/notebooks/getting_started/tutorial-1-how-to-run-a-model.ipynb
index aa50147343..226e016300 100644
--- a/docs/source/examples/notebooks/getting_started/tutorial-1-how-to-run-a-model.ipynb
+++ b/docs/source/examples/notebooks/getting_started/tutorial-1-how-to-run-a-model.ipynb
@@ -205,7 +205,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.9.0"
+ "version": "3.11.3"
}
},
"nbformat": 4,
diff --git a/docs/source/examples/notebooks/models/DFN.ipynb b/docs/source/examples/notebooks/models/DFN.ipynb
index 682adc8c21..d77a0856e3 100644
--- a/docs/source/examples/notebooks/models/DFN.ipynb
+++ b/docs/source/examples/notebooks/models/DFN.ipynb
@@ -107,7 +107,7 @@
"source": [
"Below we show how to solve the DFN model, using the default geometry, mesh, parameters, discretisation and solver provided with PyBaMM. For a more detailed example, see the notebook on the [SPM](https://github.com/pybamm-team/PyBaMM/blob/develop/docs/source/examples/notebooks/models/SPM.ipynb).\n",
"\n",
- "In order to show off all the different points at which the process of setting up and solving a model in PyBaMM can be customised we explicitly handle the stages of choosing a geometry, setting parameters, discretising the model and solving the model. However, it is often simpler in practice to use the `Simulation` class, which handles many of the stages automatically, as shown [here](../simulation-class.ipynb).\n",
+ "In order to show off all the different points at which the process of setting up and solving a model in PyBaMM can be customised we explicitly handle the stages of choosing a geometry, setting parameters, discretising the model and solving the model. However, it is often simpler in practice to use the `Simulation` class, which handles many of the stages automatically, as shown [here](../simulations_and_experiments/simulation-class.ipynb).\n",
"\n",
"First we need to import pybamm, along with numpy which we will use in this notebook."
]
diff --git a/docs/source/examples/notebooks/models/SPM.ipynb b/docs/source/examples/notebooks/models/SPM.ipynb
index 91a09a11b6..e373bdafb5 100644
--- a/docs/source/examples/notebooks/models/SPM.ipynb
+++ b/docs/source/examples/notebooks/models/SPM.ipynb
@@ -54,7 +54,7 @@
"source": [
"## Example solving SPM using PyBaMM\n",
"\n",
- "Below we show how to solve the Single Particle Model, using the default geometry, mesh, parameters, discretisation and solver provided with PyBaMM. In this notebook we explicitly handle all the stages of setting up, processing and solving the model in order to explain them in detail. However, it is often simpler in practice to use the `Simulation` class, which handles many of the stages automatically, as shown [here](../simulation-class.ipynb).\n",
+ "Below we show how to solve the Single Particle Model, using the default geometry, mesh, parameters, discretisation and solver provided with PyBaMM. In this notebook we explicitly handle all the stages of setting up, processing and solving the model in order to explain them in detail. However, it is often simpler in practice to use the `Simulation` class, which handles many of the stages automatically, as shown [here](../simulations_and_experiments/simulation-class.ipynb).\n",
"\n",
"First we need to import `pybamm`, and then change our working directory to the root of the pybamm folder. "
]
diff --git a/docs/source/examples/notebooks/models/simulate-3E-cell.ipynb b/docs/source/examples/notebooks/models/simulate-3E-cell.ipynb
new file mode 100644
index 0000000000..c92cb53465
--- /dev/null
+++ b/docs/source/examples/notebooks/models/simulate-3E-cell.ipynb
@@ -0,0 +1,210 @@
+{
+ "cells": [
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Simulating a 3E cell\n",
+ "\n",
+ "In this notebook we show how to insert a reference electrode to mimic a 3E cell."
+ ]
+ },
+ {
+ "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"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We first load a model"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "model = pybamm.lithium_ion.DFN()"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Next we use the helper function `insert_reference_electrode` to insert a reference electrode into the model. This function takes the position of the reference electrode as an optional argument. If no position is given, the reference electrode is inserted at the midpoint of the separator. The helper function adds the new variables \"Reference electrode potential [V]\", \"Negative electrode 3E potential [V]\" and \"Positive electrode 3E potential [V]\" to the model.\n",
+ "\n",
+ "In this example we will explicitly pass a position to show how it is done"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "L_n = model.param.n.L # Negative electrode thickness [m]\n",
+ "L_s = model.param.s.L # Separator thickness [m]\n",
+ "L_ref = L_n + L_s / 2 # Reference electrode position [m]\n",
+ "\n",
+ "model.insert_reference_electrode(L_ref)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Next we can set up a simulation and solve the model as usual"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ ""
+ ]
+ },
+ "execution_count": 4,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "sim = pybamm.Simulation(model)\n",
+ "sim.solve([0, 3600])"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Let's plot a comparison of the 3E potentials and the potential difference between the solid and electrolyte phases at the electrode/separator interfaces"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "d1f4e1ed03764660b87ee56b135b24a7",
+ "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": 5,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "sim.plot(\n",
+ " [\n",
+ " [\n",
+ " \"Negative electrode surface potential difference at separator interface [V]\",\n",
+ " \"Negative electrode 3E potential [V]\",\n",
+ " ],\n",
+ " [\n",
+ " \"Positive electrode surface potential difference at separator interface [V]\",\n",
+ " \"Positive electrode 3E potential [V]\",\n",
+ " ],\n",
+ " \"Voltage [V]\",\n",
+ " ]\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "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] Marc Doyle, Thomas F. Fuller, and John Newman. Modeling of galvanostatic charge and discharge of the lithium/polymer/insertion cell. Journal of the Electrochemical society, 140(6):1526–1533, 1993. doi:10.1149/1.2221597.\n",
+ "[3] 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",
+ "[4] Scott G. Marquis, Valentin Sulzer, Robert Timms, Colin P. Please, and S. Jon Chapman. An asymptotic derivation of a single particle model with electrolyte. Journal of The Electrochemical Society, 166(15):A3693–A3706, 2019. doi:10.1149/2.0341915jes.\n",
+ "[5] 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()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "venv",
+ "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.6"
+ },
+ "orig_nbformat": 4,
+ "vscode": {
+ "interpreter": {
+ "hash": "9ff3d0c7e37de5f5aa47f4f719e4c84fc6cba7b39c571a05173422444e82fa58"
+ }
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/docs/source/examples/notebooks/models/simulating-ORegan-2022-parameter-set.ipynb b/docs/source/examples/notebooks/models/simulating-ORegan-2022-parameter-set.ipynb
index 7eb647fc97..f20f385601 100644
--- a/docs/source/examples/notebooks/models/simulating-ORegan-2022-parameter-set.ipynb
+++ b/docs/source/examples/notebooks/models/simulating-ORegan-2022-parameter-set.ipynb
@@ -163,7 +163,7 @@
],
"metadata": {
"kernelspec": {
- "display_name": "Python 3.7.4 ('dev': venv)",
+ "display_name": "dev",
"language": "python",
"name": "python3"
},
@@ -177,7 +177,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.7.4"
+ "version": "3.9.16"
},
"toc": {
"base_numbering": 1,
@@ -194,7 +194,7 @@
},
"vscode": {
"interpreter": {
- "hash": "0f0e5a277ebcf03e91e138edc3d4774b5dee64e7d6640c0d876f99a9f6b2a4dc"
+ "hash": "bca2b99bfac80e18288b793d52fa0653ab9b5fe5d22e7b211c44eb982a41c00c"
}
}
},
diff --git a/docs/source/examples/notebooks/callbacks.ipynb b/docs/source/examples/notebooks/simulations_and_experiments/callbacks.ipynb
similarity index 100%
rename from docs/source/examples/notebooks/callbacks.ipynb
rename to docs/source/examples/notebooks/simulations_and_experiments/callbacks.ipynb
diff --git a/docs/source/examples/notebooks/simulations_and_experiments/custom-experiments.ipynb b/docs/source/examples/notebooks/simulations_and_experiments/custom-experiments.ipynb
new file mode 100644
index 0000000000..888c002c31
--- /dev/null
+++ b/docs/source/examples/notebooks/simulations_and_experiments/custom-experiments.ipynb
@@ -0,0 +1,235 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Custom experiments"
+ ]
+ },
+ {
+ "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"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Custom steps\n",
+ "\n",
+ "This feature is in development"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Custom termination\n",
+ "\n",
+ "Termination of a step can be specified using a few standard strings (e.g. \"4.2V\" for voltage, \"1 A\" for current, \"C/2\" for C-rate), or via a custom termination step. The custom termination step can be specified based on any variable in the model.\n",
+ "Below, we show an example where we specify a custom termination step based on keeping the anode potential above 0V, which is a common limit used to avoid lithium plating,"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ ""
+ ]
+ },
+ "execution_count": 2,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "# Set up model and parameters\n",
+ "model = pybamm.lithium_ion.DFN()\n",
+ "# add anode potential as a variable\n",
+ "# we use the potential at the separator interface since that is the minimum potential\n",
+ "# during charging (plating is most likely to occur first at the separator interface)\n",
+ "model.variables[\"Anode potential [V]\"] = model.variables[\n",
+ " \"Negative electrode surface potential difference at separator interface [V]\"\n",
+ "]\n",
+ "parameter_values = pybamm.ParameterValues(\"Chen2020\")\n",
+ "\n",
+ "\n",
+ "# Create a custom termination event for the anode potential cut-off at 0.02V\n",
+ "# We use 0.02V instead of 0V to give a safety factor\n",
+ "def anode_potential_cutoff(variables):\n",
+ " return variables[\"Anode potential [V]\"] - 0.02\n",
+ "\n",
+ "# The CustomTermination class takes a name and function\n",
+ "anode_potential_termination = pybamm.step.CustomTermination(\n",
+ " name=\"Anode potential cut-off [V]\", event_function=anode_potential_cutoff\n",
+ ")\n",
+ "\n",
+ "# Provide a list of termination events, each step will stop whenever the first\n",
+ "# termination event is reached\n",
+ "terminations = [anode_potential_termination, \"4.2V\"]\n",
+ "\n",
+ "# Set up multi-step CC experiment with the customer terminations followed\n",
+ "# by a voltage hold\n",
+ "experiment = pybamm.Experiment(\n",
+ " [\n",
+ " (\n",
+ " pybamm.step.c_rate(-1, termination=terminations),\n",
+ " pybamm.step.c_rate(-0.5, termination=terminations),\n",
+ " pybamm.step.c_rate(-0.25, termination=terminations),\n",
+ " \"Hold at 4.2V until C/50\",\n",
+ " )\n",
+ " ]\n",
+ ")\n",
+ "\n",
+ "# Set up simulation\n",
+ "sim = pybamm.Simulation(model, parameter_values=parameter_values, experiment=experiment)\n",
+ "\n",
+ "# for a charge we start as SOC 0\n",
+ "sim.solve(initial_soc=0)\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ ""
+ ]
+ },
+ "execution_count": 3,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAABKUAAAGGCAYAAACqvTJ0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB7f0lEQVR4nO3dd3xT9f7H8Xfa0gFpyy6rWJANshGLC7WCoigOVFAKiAsBRdQLKMJFxYKCF+4FUZB1vfIDB3ivC8RKXeyloGxBECgFgRYKdOX8/jg2NHSXJidpX8/H4zySfPM953xCQr49n3yHzTAMQwAAAAAAAIAH+VkdAAAAAAAAAMofklIAAAAAAADwOJJSAAAAAAAA8DiSUgAAAAAAAPA4klIAAAAAAADwOJJSAAAAAAAA8DiSUgAAAAAAAPA4klIAAAAAAADwOJJSAAAAAAAA8DiSUoCX279/v2w2m2w2m9q2bVvs/bP3rVy5cqnHBgDlTVRUlKZOnWp1GHmivQCA0mWz2fTJJ59YHYbHDRgwQL169Spy/ez2Z8uWLfnWiYqKcrYzp06dKlY8Xbt2de5b0Dngm0hKoUxITEzUsGHD1LBhQwUFBSkyMlI9e/ZUfHy81aEVqDgN3ddff53n6/njjz8UGBioVq1a5bnfkSNHvPYCCgA8pWfPnrrlllvyfO7777+XzWbTzz//XOzjeuMFS872YtiwYWrevHme9Q4cOCB/f3/973//k0R7AcB7rV69Wv7+/rrtttusDsUS7vpBJL9k0rRp0zR//vxSP9/LL7+sI0eOKDw8XB9//LH8/f116NChPOs2btxYI0aMkCQtWbJE69atK/V44B1ISsHn7d+/Xx06dNA333yjN954Q1u3btWyZct0ww03aMiQISU+rmEYyszMzFWenp5+KeGWWLVq1VStWrVc5fPnz9d9992nlJQUrV27NtfztWrVUnh4uCdCBACvNWjQIK1YsUJ//PFHrufmzZunjh07qnXr1hZEVvpytheDBg3Sjh07tGrVqlz15s+fr5o1a6pHjx6SaC8AeK85c+Zo2LBh+u6773T48GGrwynzwsPD3dJrNjQ0VLVq1ZLNZtMdd9yhatWqacGCBbnqfffdd9qzZ48GDRokSapatapq1KhR6vHAO5CUgs978sknZbPZtG7dOt1zzz1q0qSJWrZsqREjRmjNmjWS8v4V4NSpU7LZbEpISJAkJSQkyGaz6csvv1SHDh0UFBSkH374QV27dtXQoUM1fPhwVa9eXd27d5ckbdu2TbfeeqvsdrsiIiLUr18/HT9+3Hn8rl276qmnntLf/vY3Va1aVbVq1dLf//535/NRUVGSpLvuuks2m835uDgMw9C8efPUr18/9e3bV3PmzCn2MQCgPLj99ttVo0aNXL/8njlzRh9++KHzD9+PP/5YLVu2VFBQkKKiojRlypR8j5nf9/jevXt15513KiIiQna7XZ06ddLXX3/tsu+RI0d02223KSQkRA0aNNDChQtz/RJ+6tQpPfLII6pRo4bCwsJ044036qeffirW627btq3at2+vuXPnupQbhqH58+erf//+CggIKNYxAcCTzpw5o8WLF2vw4MG67bbbcn2PZ/8NHx8fr44dO6pixYrq0qWLdu7c6VJv5syZuvzyyxUYGKimTZvqvffec3l+9+7duu666xQcHKwWLVpoxYoVuWI5ePCg7rvvPlWuXFlVq1bVnXfeqf379+cbe3Zsn3/+uVq3bq3g4GBdddVV2rZtm0u9gtqerl276vfff9czzzzjHMKW7YcfftC1116rkJAQRUZG6qmnnlJqaqrz+aioKL322mt6+OGHFRoaqvr162vWrFnO5xs0aCBJateunWw2m7p27Sop9/C9ZcuW6ZprrlHlypVVrVo13X777dq7d2++r7soKlSooH79+uXZI2vu3Lnq3LmzWrZseUnngG8gKQWfduLECS1btkxDhgxRpUqVcj1fkgz/qFGjNHHiRG3fvt35q/mCBQsUGBioH3/8UW+//bZOnTqlG2+8Ue3atdOGDRu0bNkyHT16VPfdd5/LsRYsWKBKlSpp7dq1ev311/Xyyy87G7j169dLMn+hP3LkiPNxcaxcuVJnz55VTEyMHnroIS1atMilIQIAmAICAhQbG6v58+fLMAxn+YcffqisrCz16dNHGzdu1H333acHHnhAW7du1d///ne99NJL+Q5hyO97/MyZM+rRo4fi4+O1efNm3XLLLerZs6cOHDjg3Dc2NlaHDx9WQkKCPv74Y82aNUtJSUkux+/du7eSkpL05ZdfauPGjWrfvr1uuukmnThxolivfdCgQfrggw9c2oeEhATt27dPDz/8cLGOBQCe9sEHH6hZs2Zq2rSpHnroIc2dO9flezzbiy++qClTpmjDhg0KCAhw+X5bunSpnn76aT377LPatm2bHn/8cQ0cOFArV66UJDkcDt19990KDAzU2rVr9fbbb2vkyJEux8/IyFD37t0VGhqq77//Xj/++KPsdrtuueWWQkdSPP/885oyZYrWr1+vGjVqqGfPnsrIyJCkQtueJUuWqF69es6hb0eOHJFk/gByyy236J577tHPP/+sxYsX64cfftDQoUNdzj1lyhR17NhRmzdv1pNPPqnBgwc7E3bZQ+K+/vprHTlyREuWLMkz/tTUVI0YMUIbNmxQfHy8/Pz8dNddd8nhcBT4ugszaNAg7d69W999952z7MyZM/roo4+cPxahHDAAH7Z27VpDkrFkyZIC6+3bt8+QZGzevNlZdvLkSUOSsXLlSsMwDGPlypWGJOOTTz5x2ff666832rVr51L2yiuvGN26dXMpO3jwoCHJ2Llzp3O/a665xqVOp06djJEjRzofSzKWLl1a7Niz9e3b1xg+fLjzcZs2bYx58+blqjdv3jwjPDy8wPMAQFm3fft2l+99wzCMa6+91njooYcMwzC/U2+++WaXfZ5//nmjRYsWzseXXXaZ8Y9//MP5uCjf44ZhGC1btjT+9a9/ucSxfv165/O7d+82JDmP/f333xthYWHG+fPnXY5z+eWXG++8806e58ivvTh58qQRHBzs0j7069cvVxtlGLQXALxPly5djKlTpxqGYRgZGRlG9erVXb7Hs/+G//rrr51ln3/+uSHJOHfunPMYjz76qMtxe/fubfTo0cMwDMNYvny5ERAQYBw6dMj5/JdffunyHf/ee+8ZTZs2NRwOh7NOWlqaERISYixfvjzP2LNjW7RokbPszz//NEJCQozFixcbhlGytscwDGPQoEHGY4895lL2/fffG35+fs7XfdlllznbOMMwDIfDYdSsWdOYOXOmYRj5txv9+/c37rzzzjxfk2EYxrFjxwxJxtatWws8Tk55vQbDMIyrrrrK6N+/v/PxnDlzjIoVKxopKSku9YpyDvgmekrBpxl5/EpyqTp27JirrEOHDi6Pf/rpJ61cuVJ2u925NWvWTJJcurJePD9J7dq1c/0SXlKnTp3SkiVL9NBDDznLHnroIYbwAUA+mjVrpi5dujiHsu3Zs0fff/+989fY7du36+qrr3bZ5+qrr9bu3buVlZVV5POcOXNGzz33nJo3b67KlSvLbrdr+/btzp5SO3fuVEBAgNq3b+/cp1GjRqpSpYrz8U8//aQzZ86oWrVqLm3Nvn37ij1konLlyrr77rudrzslJUUff/wxv0ID8Ho7d+7UunXr1KdPH0lmr9f7778/z793c/7dXbt2bUly/t2d3/f79u3bnc9HRkaqTp06zuejo6Nd6v/000/as2ePQkNDnd/JVatW1fnz5wv9Xs55rKpVq6pp06Yu5y5J2/PTTz9p/vz5Lm1E9+7d5XA4tG/fvjz/XWw2m2rVqlXs65Hdu3erT58+atiwocLCwpzD1XP2AC6phx9+WB999JFOnz4tyRy617t3b4WGhl7yseEbmEQAPq1x48ay2WzasWNHgfX8/Mz8a84kVnaX2YvlNQzw4rIzZ86oZ8+emjRpUq662Y2gZI6Vzslms11yN9dsCxcu1Pnz59W5c2dnmWEYcjgc2rVrl5o0aVIq5wGAsmTQoEEaNmyYZsyYoXnz5unyyy/X9ddfX6rneO6557RixQpNnjxZjRo1UkhIiO69995iLZRx5swZ1a5d2znvYU4lGZo+aNAg3XTTTdqzZ49Wrlwpf39/9e7du9jHAQBPmjNnjjIzM12SRYZhKCgoSNOnT3dZnCHn393Z8y6V1t/dkvm93KFDB73//vu5nrNiEu4zZ87o8ccf11NPPZXrufr16zvvl8b1SM+ePXXZZZdp9uzZqlOnjhwOh1q1alUqC0A98MADeuaZZ/TBBx/ouuuu048//qi4uLhLPi58B0kp+LSqVauqe/fumjFjhp566qlcyaNTp06pcuXKzobiyJEjateunSTlWvq0ONq3b6+PP/5YUVFRlzRBbIUKFYr163tOc+bM0bPPPqsBAwa4lD/55JOaO3euJk6cWOK4AKCsuu+++/T0009r4cKF+ve//63Bgwc7L16aN2+uH3/80aX+jz/+qCZNmsjf3z/P4+X1Pf7jjz9qwIABuuuuuySZFw45J8Jt2rSpMjMztXnzZmdP3D179ujkyZPOOu3bt1diYqICAgJKtBDGxW644QY1aNBA8+bN08qVK/XAAw/k+SMMAHiLzMxM/fvf/9aUKVPUrVs3l+d69eql//u//9MTTzxRpGNlf7/379/fWfbjjz+qRYsWzucPHjyoI0eOOH9gzl4wKVv79u21ePFi1axZU2FhYcV6LWvWrHEmik6ePKldu3apefPmLrHldHHbExgYmKutad++vX799Vc1atSoWLHkFBgYKEkFXo/8+eef2rlzp2bPnq1rr71WkjnBemkJDQ1V7969NXfuXO3du1dNmjRxngflA8P34PNmzJihrKwsXXnllfr444+1e/dubd++Xf/85z+dXWVDQkJ01VVXOScw//bbbzVmzJgSn3PIkCE6ceKE+vTpo/Xr12vv3r1avny5Bg4cWKwkU1RUlOLj45WYmOhyMVKYLVu2aNOmTXrkkUfUqlUrl61Pnz5asGCBMjMzS/LSAKBMs9vtuv/++zV69GgdOXLEJbH/7LPPKj4+Xq+88op27dqlBQsWaPr06XruuefyPV5e3+ONGzfWkiVLtGXLFv3000/q27evy6/SzZo1U0xMjB577DGtW7dOmzdv1mOPPaaQkBBngiwmJkbR0dHq1auXvvrqK+3fv1+rVq3Siy++qA0bNhT7ddtsNj388MOaOXOmVq9ezdA9AF7vs88+08mTJzVo0KBcf+/ec889xZqy4vnnn9f8+fM1c+ZM7d69W2+++aaWLFni/H6PiYlRkyZN1L9/f/3000/6/vvv9eKLL7oc48EHH1T16tV155136vvvv9e+ffuUkJCgp556Sn/88UeB53/55ZcVHx+vbdu2acCAAapevbpzdbuitD1RUVH67rvvdOjQIedq3yNHjtSqVas0dOhQbdmyRbt379Z///vfXBOdF6RmzZoKCQlxLtqUnJycq06VKlVUrVo1zZo1S3v27NE333yjESNGFPkcRTFo0CCtWrVKb7/9NgtwlEMkpeDzGjZsqE2bNumGG27Qs88+q1atWunmm29WfHy8Zs6c6aw3d+5cZWZmqkOHDho+fLheffXVEp+zTp06+vHHH5WVlaVu3brpiiuu0PDhw1W5cmXnUMGimDJlilasWKHIyEhnD66imDNnjlq0aOGcxyqnu+66S0lJSfriiy+KfDwAKE8GDRqkkydPqnv37i5DQtq3b68PPvhAixYtUqtWrTR27Fi9/PLLuXqk5pTX9/ibb76pKlWqqEuXLurZs6e6d+/uMn+UJP373/9WRESErrvuOt1111169NFHFRoaquDgYElmEumLL77Qddddp4EDB6pJkyZ64IEH9PvvvysiIqJEr3vAgAFKTk5Wy5YtXYZ+A4A3mjNnjmJiYlyG6GW75557tGHDBv38889FOlavXr00bdo0TZ48WS1bttQ777yjefPmqWvXrpLMqT6WLl2qc+fO6corr9QjjzyiCRMmuByjYsWK+u6771S/fn3dfffdat68uQYNGqTz588X2nNq4sSJevrpp9WhQwclJibq008/dfZSKkrb8/LLL2v//v26/PLLnSNAWrdurW+//Va7du3Stddeq3bt2mns2LEu7VphAgIC9M9//lPvvPOO6tSpozvvvDNXHT8/Py1atEgbN25Uq1at9Mwzz+iNN94o8jmK4pprrlHTpk2VkpKi2NjYUj02vJ/NcMdM0QBKzf79+9WgQQNt3rxZbdu2LdEx5s+fr+HDh+vUqVOlGhsAoHT88ccfioyM1Ndff62bbrqpRMegvQAA75KQkKAbbrhBJ0+eLNF8gGVJVFSUhg8fruHDh5do/9Jo4+Cd6CkF+IguXbqoS5cuxd7PbrcXebw9AMAzvvnmG/3vf//Tvn37tGrVKj3wwAOKiorSddddd8nHpr0AAHijkSNHym635zlMsCC33nqrWrZs6aaoYDUmOge8XL169bR7925JUlBQULH3z57QPb9JegEAnpeRkaEXXnhBv/32m0JDQ9WlSxe9//77uVZJKg7aCwCAt/r222+dq5+HhoYWa993331X586dk+S6siDKBobvAQAAAAAAwOMYvgcAAAAAAACPIykFAAAAAAAAjyMpBQAAAMCjvvvuO/Xs2VN16tSRzWbTJ598Uug+CQkJat++vYKCgtSoUSPNnz/f7XECANzLqyc6dzgcOnz4sEJDQ2Wz2awOBwDKNMMwdPr0adWpU0d+fr71mwXtBQB4Tmm0F6mpqWrTpo0efvhh3X333YXW37dvn2677TY98cQTev/99xUfH69HHnlEtWvXVvfu3Yt0TtoKAPCcorYVXj3R+R9//KHIyEirwwCAcuXgwYOqV6+e1WEUC+0FAHheabUXNptNS5cuVa9evfKtM3LkSH3++efatm2bs+yBBx7QqVOntGzZsiKdh7YCADyvsLbCq3tKZS8VefDgQYWFhVkcDQCUbSkpKYqMjCz2Mr3egPYCADzHivZi9erViomJcSnr3r27hg8fnu8+aWlpSktLcz7O/i2etgIA3K+obYVXJ6Wyu9UGBATQcACAh/jikAbaCwDwPE+2F4mJiYqIiHApi4iIUEpKis6dO6eQkJBc+8TFxWn8+PG5ymkrAMBzCmsrfGLSEC8eYQgA8CK0FwCAbKNHj1ZycrJzO3jwoCTaCgDwJh5JSs2YMUNRUVEKDg5W586dtW7dOk+cFgAAAEAZUKtWLR09etSl7OjRowoLC8uzl5QkBQUFKSwszGUDAHgXtyelFi9erBEjRmjcuHHatGmT2rRpo+7duyspKcndpwYAAABQBkRHRys+Pt6lbMWKFYqOjrYoIgBAaXB7UurNN9/Uo48+qoEDB6pFixZ6++23VbFiRc2dO9fdpwYAAADghc6cOaMtW7Zoy5YtkqR9+/Zpy5YtOnDggCRz6F1sbKyz/hNPPKHffvtNf/vb37Rjxw699dZb+uCDD/TMM89YET4AoJS4NSmVnp6ujRs3uqyU4efnp5iYGK1evdqdpwYAAADgpTZs2KB27dqpXbt2kqQRI0aoXbt2Gjt2rCTpyJEjzgSVJDVo0ECff/65VqxYoTZt2mjKlCl699131b17d0viBwCUDreuvnf8+HFlZWXluVLGjh07ctW/eNnWlJQUd4YHAAAAwAJdu3YtcMLx+fPn57nP5s2b3RgVAMDTvGr1vbi4OIWHhzu3yMhISb65PDkAwPNoLwAAhaGtAADv4dakVPXq1eXv75/nShm1atXKVT+/ZVsrVqzozjABAG4yceJE2Ww2DR8+PN86s2fP1rXXXqsqVaqoSpUqiomJKfEqrbQXAIDC0FYAgPdwa1IqMDBQHTp0cFkpw+FwKD4+Ps+VMli2FQDKjvXr1+udd95R69atC6yXkJCgPn36aOXKlVq9erUiIyPVrVs3HTp0yEORAgAAALCCW+eUksxJC/v376+OHTvqyiuv1NSpU5WamqqBAwe6+9QAAIucOXNGDz74oGbPnq1XX321wLrvv/++y+N3331XH3/8seLj411WXiqK1NRUhYaGOodmpKenKyMjQwEBAQoKCnKpJ0khISHy8zN/n8nIyFB6err8/f0VHBxcorpnz56VYRgKDg6Wv7+/JCkzM1NpaWny8/NTSEhIieqeO3dODodDQUFBCggwm+6srCydP3++WHVtNptLD4Hz588rKytLgYGBqlChQrHrOhwOnTt3TpJUqVIlZ920tDRlZmaqQoUKCgwMLHZdwzB09uxZSWaPhovfz+LULcp7Xxqfk7zez9L4nGS/n5f6Obn4/bzUz0l+7+elfk5yvp+X+jnJ7/0s6eeE74gL3xHZrxkAgEvl9jml7r//fk2ePFljx45V27ZttWXLFi1btizX5OcFyf5DBgBQSgxDSkqSfvxRmjdPGj1aeuihUjv8kCFDdNttt7msvlpUZ8+eVUZGhqpWrZpvnbS0NKWkpLhsklSnTh0dP37cWe+NN96Q3W7X0KFDXfavWbOm7Ha7y8pOM2bMkN1u16BBg1zqRkVFyW63a/v27c6y+fPny26364EHHnCp26JFC9ntdm3atMlZtnjxYtntdt1xxx0udTt16iS73a7vv//eWfbZZ5/Jbrfn+ne77rrrZLfbtXz5cmfZN998I7vdnqvn8a233iq73a6lS5c6y9asWSO73a42bdq41L3nnntkt9tdEoNbt26V3W5X48aNXer269dPdrtds2bNcpbt3btXdrtddevWdan7+OOPy263a9q0ac6yI0eOyG63q3Llyi51R4wYIbvdrtdee81ZlpycLLvdLrvdrszMTGf5iy++KLvdrhdffNFZlpmZ6aybnJzsLH/ttddkt9s1YsQIl/NVrlxZdrtdR44ccZZNmzZNdrtdjz/+uEvdunXrym63a+/evc6yWbNmyW63q1+/fi51GzduLLvdrq1btzrL3n//fdntdt1zzz0uddu0aSO73a41a9Y4y5YuXSq73a5bb73VpW50dLTsdru++eYbZ9ny5ctlt9t13XXXudSNiYmR3W7XZ5995iz7/vvvZbfb1alTJ5e6d9xxh+x2uxYvXuws27Rpk+x2u1q0aOFS94EHHpDdbneZeHr79u2y2+2KiopyqTto0CDZ7XbNmDHDWXbgwAHZ7XbVrFnTpe7QoUNlt9v1xhtvOMuOHz/ufD8lSX/+Kb36qkY2aSK73a7xTZtKMTFSTIzO3nijs+7ZG290lo9v2lR2u10jGzd2lumvfxu73a7j11/vLHujeXPzO6JhQ5e6NcPCzO+Ia691ls1o2dL8joiKcqkbVaWK+R3RpYuzbH7r1uZ3RGSkS90W1aqZ3xHR0c6yxW3bmt8R9eq51O301/fU91de6Sz7rEMH8zuidm2XutdFRJjfER07Osu+6djR/I6oWdOl7q116pjfEe3aOcvWdO5sfkdUr+5S95569czviDZtnGVbu3RRnTp15MvOlXCIOACg9Lm9p5Rk/tFx8QVBcTgcjlKMBgDKkcxMac8e6eefpe3bpd27pV27zC3HBXxpWrRokTZt2qT169eXaP+RI0eqTp06BSa04uLiNH78+JKGCMAXDBsmzZ0r/dUbSpK0f7+5XSwhIXfZwYPmdrEciWCnw4fN7WKrVuUuO3rU3C62dm3usuPHpRzTWDjl9f144kTedXMkuZ2Sk/Ou+9NPucvOnMm77rZt5pbTuXN5192+3dzKCMeff1odAgDgLzajoLVYLZaSkqLw8HAdPnxYtWvXtjocAPBux49LmzebCaitW83bX3+V0tLyrm+zSfXrS02aSE2aKCUyUuGjRik5ObnEc/odPHhQHTt21IoVK5xzSXXt2lVt27bV1KlTC91/4sSJev3115WQkFDgXFRpaWlKy/G6UlJSFBkZqT179qhhw4YMzWH4HsP3fHH4XmamjC++0NmZM6Xly1Up+0/Utm2VNnCgMkNDVSEgQIF/xWsYhs7+9T1QMSjownufmamMzEwF+Psr6K8YJCn1/Pli1w0JDLzw3mdmKj0zU/5+fgr+6/NX3Lpn09LM9z4wUP5/1c3MylJaRob5fpaw7rn0dPP9rFBBAX99TrIcDp1PTy9WXZvNpoo5/g+cT09XlsOhwIAAVcj+nDgcSkpOVp0hQy6pvbCC89ri/fdVu29fq8MBgDIt+zu3sLaCpBQA+KIzZ8xfz9evl9atM2/37cu7bsWK0hVXSC1bSk2bSo0bm4moyy+XclwoF7XhKMgnn3yiu+66y3nxLJkXujabTX5+fkpLS3N5LqfJkyfr1Vdf1ddff62OHTsW67y0F4AP27/f7BE1Z45rb6VbbpGee0668UYziQ6vURrthRWcbcWCBapdzDkLAQDFU9S2wiPD9wAAl+jIEenbb83txx+lX36R8hra3Lix1KaNmYRq3dq8bdBA8nP7FIKSpJtuusllTh1JGjhwoJo1a6aRI0fmm5B6/fXXNWHCBC1fvrzYCSkAPuj336WPPpI+/NB12Fv16tKAAdIjj5hJdMAdcsxVBwCwFkkpAPBGBw9eSEJ9+605F9TF6taVrrxS6tTJvO3QQbpoEmlPCw0NVatWrVzKKlWqpGrVqjnLY2NjVbduXcXFxUmSJk2apLFjx2rhwoWKiopSYmKiJLlOdgzAt2VlSRs2SCtWSJ9+avbwzGazmb2hHntMuvNOKcfwMcAtsrKsjgAA8BeSUgDgDTIyzB5Qn38uffGFORdUTjab1LatdP310nXXSVddJfnoMLUDBw44512RpJkzZyo9PV333nuvS71x48bp73//u4ejA4opJUV6913JbjeTKjBlZppz261ZY06c/c030smTF5632czvsvvuk+6+W6pVy7pYUf7QUwoAvAZJKQCwypEj0pdfmkmoFSvMi9tsfn5mz6frrze3a66xvBdUSSVctCrWxY/357WSFuDt0tOld96Rxo+X/vzTTLL06SOFhlodmeedPSvt2GEOK96yxewFtXGjuZJbTpUrmz2ibr5Z6tWLRBSsQ1IKALyGTySlcq7UAgA+7dChC/Oo/Pij63M1aki33irddpt50ValijUx+jDaC7idYZj/f194Qdq717U8I8O6uNzJMMzE28GD5oIK+/ebt7/9Jm3fbj7Oa92c8HBzaPG115rfaR07SgE+8acnyrhKfA4BwGvwjQwA7vbnn9L//Z+0aFHuRFTHjmYS6rbbzJ5RHpqQHEAJbNkiPfHEhYm5IyKkMWOkYcMsDavEzp+XEhPN7ciR3Pezb48eLTzhVqOG1KKF1KqVmYi68kpzlU++0+CN6CkFAF6DpBQAuENGhjk0b8ECc1LfnBd0V18t9e4t3XOPVK+edTECKJrUVHOY3ptvmhMk2+3S889LI0ZIISHel5RKSZEOHMidXLr4/qlTxTtujRrmap7ZW1SU1KyZmYyqUcMdrwRwD5JSAOA1fCIpde7cOYWFhVkdBgAUbv9+c56ZuXOlpKQL5e3bS/36mcmounUtC6+so71AqVu2TBo82Py/LZn/h6dNu7DQgBWreGVmSnv2mEPn9u2Tfv/djO/3382tOMmmoCBzbqfatc3bnPdzlkVESIGB7npFgEedO3dOtBQA4B18IinlcDisDgEA8udwSF99Jb31lvTZZxfmVomIkB56SOrfX7riCmtjLCdoL1BqUlKkp54yeztKUmSk+X/89ts9F0NGhpl8+uUXc0XO7NudOwsfTle1qlSnjmti6eJEU61a5uTjNptHXg7gLRxldf43APBBPpGUAgCvlJYm/fvf0uTJ0q5dF8pjYqQnn5R69mRSX8AXff+9FBtr9j7y85Oeflp6+WVz2F5ebDYzGf3TT9INNxT/fOnp0u7dromnX381v1fyu3iuWFFq3lxq1Ei67DJzKN1ll13Y8osVAMP3AMCLcLUEAMWVkmIO0fvHP8z5WSRzlakBA8xhPk2bWhoegBJKT5fGjZMmTTKTTA0aSO+9Z84Dlx9/f+mBB8zFDO69V1qzRmrcOO+6Z86Yyaddu1wTULt353+RbLebcza1aCG1bHnhfv36TCIOlBRJKQDwGiSlAKCoUlLMRNQ//iElJ5tldeuakx0/+qgUGmptfABKbs8e6f77pU2bzMcDBphzRxVljrJ33zX3X79euvlms5dktWpmEurYMbPH1a5d5uTi+QkNdU08Zd9GRjK8DihtVswFBwDIE0kpACjMuXPmXDJxcdKff5plzZtLf/ub1Lcvk/8Cvu7jj6WBA6XTp825mGbPlu6+u+j7V6xorrJ51VVmAmr69PzrVq9u9qS6uOdTvXoknwBPoacUAHgNklIAkJ+sLHMVvfHjpUOHzLKmTc25Ze69l6EzgK9LT5dGjpSmTjUfX3ONtGhRyVbIjIiQ1q6VFi82e0T9+afZ+6laNbO3U+PG5lalSqm+BAAlQFIKALwGSSkAyMsPP0jDhklbtpiP69eX/v53qV8/Ji8HyoKDB6X77jPngJKk55+XJkyQKlQo+TFr1jS/NwB4N5JSAOA1fOLKqlKlSlaHAKC8OHzYHJb3/vvm48qVzYmPBw+WgoIsDQ2Fo71AkSxbJj30kNmbqXJlacEC6Y47rI4KgIdUoqczAHgNvpEBQJIcDnNS46ZNzYSUzSY98og5OfHw4SSkgLIgK0saO1bq0cNMSHXoYE5sTkIKKF/oKQUAXsMnekoBgFvt3i09/LA5ZE+SOnc2Jyru2NHauACUnuPHzYUJVqwwHw8eLL35phQcbG1cADyP1fcAwGv4RE+p8+fPWx0CgLLI4TAnOG7TxkxI2e3SzJnSqlUkpHwU7QXytG6d1L69mZCqWNHsDfnWWySkgHKKtgIAvIdP9JTK4tcMAKXt99/NOWWye0fddJP07rtSVJSlYeHS0F7AhWFI77wjPf20udJekybSxx9LrVpZHRkAC2VlZFgdAgDgLz7RUwoAStVnn0nt2l3oHfX222YPChJSQNlx9qw0YIA5TC89XbrrLmn9ehJSABi+BwBexCd6SgFAqcjMlMaMkSZNMh936iQtXiw1aGBtXABK19690t13Sz//LPn5SRMnSs89Zy5gAAAkpQDAa5CUAlA+HDok9ekjff+9+XjYMOmNN1hVDyhrFi6UnnhCOn1aqlnTTDx37Wp1VAC8CcP3AMBrkJQCUPZt3CjdfruUmCiFhkpz5ki9e1sdFYDSlJIiDR0qvfee+fjqq82EVN261sYFwPtkZlodAQDgL8wpBaBs+/RT6brrzIRUq1ZmgoqEFFC2rFtnzhP33nvmcL2//11KSCAhBSBvJKUAwGuQlAJQdv3rX1KvXuaExzffbE5s3rix1VEBKC3p6dL48WavqN9+k+rXl779Vho3TgqgMziAfDCnFAB4DZ/4i61SpUpWhwDAl2RlSSNGSP/8p/n4kUekt96SKlSwNi64He1FObJli7m63k8/mY/vu0965x2pcmULgwLgC2gpAMB70FMKQNmSkSE99NCFhFRcnDRrFgkpoKxITzeH53XqZCakqlWT/u//pEWLSEgBKBqG7wGA1/CJnlIAUCRpadIDD0iffGIO3XnvPfMxgLJhzRrp8celn382H999t9kLMiLC2rgA+BaSUgDgNXyip9T58+etDgGAt0tPl+6910xIBQVJS5eSkCqHaC/KqJMnpSeekLp0MRNS1aqZK+t99BEJKQDFdj493eoQAAB/8YmkVBaTEQIoSGam1KeP9NlnUnCwueLe7bdbHRUkTZw4UTabTcOHDy+w3ocffqhmzZopODhYV1xxhb744osSnY/2oowxDOk//5GaNTPnizIMcx6pHTvMOaRsNqsjBOCDsugpBQBewyeSUgCQL4dD6t9fWrJECgyU/vtfc6U9WG79+vV655131Lp16wLrrVq1Sn369NGgQYO0efNm9erVS7169dK2bds8FCm80ubNUteuUr9+UlKS1Ly5lJAgzZsnVa9udXQAfBlJKQDwGiSlAPguw5CefVZauNCcQ+rjj6Vu3ayOCpLOnDmjBx98ULNnz1aVKlUKrDtt2jTdcsstev7559W8eXO98sorat++vaZPn+6haOFVkpKkRx+VOnSQvvtOCgmRXnvNXG3v+uutjg5AWUCvWgDwGiSlAPiuKVOkqVPN+//+N0P2vMiQIUN02223KSYmptC6q1evzlWve/fuWr16db77pKWlKSUlxWWDj0tPlyZPlho3lt5910w69+kj7dwpjR5t9oQEgNJATykA8BqsvgfAN334ofT88+b9yZPNi1d4hUWLFmnTpk1av359keonJiYq4qLJqiMiIpSYmJjvPnFxcRo/fvwlxQkvYRjmfHDPPivt3m2WdeggTZsmXX21tbEBKJvoKQUAXoOeUgB8z6ZN5jxSkvTUU+bFLLzCwYMH9fTTT+v9999XcHCw284zevRoJScnO7eDBw+67Vxwo7VrpRtvlO64w0xI1aolzZ0rrVtHQgqA+2RkWB0BAOAv9JQC4FsSE6U775TOnZNuucUcwgevsXHjRiUlJal9+/bOsqysLH333XeaPn260tLS5O/v77JPrVq1dPToUZeyo0ePqlatWvmeJygoSEFBQaUbPDxnxw7pxRfNBQokKShIGj7cLAsNtTQ0AOUAw/cAwGv4RE+pihUrWh0CAG+QmSk98ID0xx/mEvGLFpkTnMNr3HTTTdq6dau2bNni3Dp27KgHH3xQW7ZsyZWQkqTo6GjFx8e7lK1YsULR0dHFPj/thZc7dMicxLxVKzMh5ecnDRwo7dolTZxIQgqAR1Rk+B4AeA2fuJqz2WxWhwDAG4wbJ337rWS3S//9rxQebnVEuEhoaKhatWrlUlapUiVVq1bNWR4bG6u6desqLi5OkvT000/r+uuv15QpU3Tbbbdp0aJF2rBhg2bNmlXs89NeeKmTJ82k0z//KZ0/b5bdcYe5ql7LltbGBqDcsZGUAgCv4RM9pQBAX35pXsBK5spcTZpYGw9K7MCBAzpy5IjzcZcuXbRw4ULNmjVLbdq00UcffaRPPvkkV3ILPujcOWnSJKlhQ+n1182E1DXXSD/8YCaWSUgBsAJJKQDwGj6RlEpLS7M6BABWOnpUio017w8ZIt1/v7XxoFgSEhI0depUl8fz5893qdO7d2/t3LlTaWlp2rZtm3r06FGic9FeeInMTGn2bKlRI2nUKOnUKXPI3qefSt99xyTmACRJM2bMUFRUlIKDg9W5c2etW7euwPpTp05V06ZNFRISosjISD3zzDM6n937shjSmOgcALyGTySlMpmMECi/DMOcg+b4cal1ayY2R4FoLyxmGOZcUa1aSY89Jh0+LNWvLy1YIG3ZIt1+u8QQSwCSFi9erBEjRmjcuHHatGmT2rRpo+7duyspKSnP+gsXLtSoUaM0btw4bd++XXPmzNHixYv1wgsvFPvctBUA4D18IikFoBybN8/sXREYKL33nrlKFwDvk5AgRUdL99wj7dwpVasm/eMf5iTmsbFSHpPcAyi/3nzzTT366KMaOHCgWrRoobffflsVK1bU3Llz86y/atUqXX311erbt6+ioqLUrVs39enTp9DeVXkiKQUAXoOkFADvdfCg9PTT5v1XXjF7SgHlmcMhPfec9OyzVkdywU8/SbfeKt1wg7R2rVSxojRmjLR3rzR8OIlkALmkp6dr48aNiomJcZb5+fkpJiZGq1evznOfLl26aOPGjc4k1G+//aYvvviiZMO9HQ6zZycAwHJuW31vwoQJ+vzzz7VlyxYFBgbq1KlT7joVgLJq2DDpzBmpSxfvuggHrPLqqxeGsI4ZI1WpYl0s+/ZJL70kLVxoXtwFBJhD9l56SapVy7q4AHi948ePKysrSxERES7lERER2rFjR5779O3bV8ePH9c111wjwzCUmZmpJ554osDhe2lpaS5zDaakpFx40uGgBycAeAG39ZRKT09X7969NXjwYHedAkBZ9skn5upcAQHSrFn84QgsWyb9/e8XHjsc1sSRlCQ99ZTUtKn0/vtmQuqBB6Tt26UZM0hIAXCLhIQEvfbaa3rrrbe0adMmLVmyRJ9//rleeeWVfPeJi4tTeHi4c4uMjLzwJEP4AMAruK2n1Pjx4yUp1wpLAFCo06fNXlKS9PzzLBsPHDggPfigtcNNTp+W3nxTmjzZ7MEoSd26SXFxUvv21sUFwOdUr15d/v7+Onr0qEv50aNHVSufxPZLL72kfv366ZFHHpEkXXHFFUpNTdVjjz2mF198UX5+uX9rHz16tEaMGOF8nJKSciExlZnJ8GIA8ALMKQXA+8TFSX/8ITVoYA5RAsqzjAzp/vulEyekDh08f/70dOlf/5Iuv9zsqXXmjNSxo/T119Ly5SSkABRbYGCgOnTooPj4eGeZw+FQfHy8oqOj89zn7NmzuRJP/n/1ojbySdgHBQUpLCzMZXOipxQAeAW39ZQqifzGfVesWNGqkAB42oED5opdknnL/38UQ5lsL158UVqzRgoPlxYvlho18sx5HQ5p0SIzMbxvn1nWuLE0YYJ0772SzeaZOACUSSNGjFD//v3VsWNHXXnllZo6dapSU1M1cOBASVJsbKzq1q2ruLg4SVLPnj315ptvql27durcubP27Nmjl156ST179nQmp4qqomQm/AEAlitWUmrUqFGaNGlSgXW2b9+uZs2alSiYuLg457C/nGz84QuUHy++KJ0/L11/vXTHHVZHAx9T5tqLzz+X3njDvD93rtl70N0Mw+wBNXq0tGWLWVa7tjRunPTww1KFCu6PAUCZd//99+vYsWMaO3asEhMT1bZtWy1btsw5+fmBAwdcekaNGTNGNptNY8aM0aFDh1SjRg317NlTEyZMKPa5bRI9pQDAS9iM/Pq75uHYsWP6888/C6zTsGFDBQYGOh/Pnz9fw4cPL9Lqe3n1lIqMjFRycrJrd1sAZdPGjeawIElav/7CfXhESkqKwsPDffI715djz9cff0ht20p//ikNHWoOocu5WtTx41K1aqV7zrVrpVGjpIQE83FYmDRypPT001KlSqV7LgA+y1e/c51xSwo7eFCqV8/qkACgzCpqW1GsnlI1atRQjRo1Ljm4/AQFBSkojwkHcyaqAJRh2cs6P/QQCSmUSJlpLzIzpT59zIRU+/bm5OLutHOn+f9vyRLzcVCQmQgbPbr0E18AYLE0iZ5SAOAl3Dan1IEDB3TixAkdOHBAWVlZ2vLXEIBGjRrJbrcX61iZNBpA2bd+vfTVV2YvkJdftjoa+Kgy016MHSv98IMUGmrOI+WuFaIOHDD/v82fL2VlSX5+Uv/+5oTm9eu755wAYLFMiaQUAHgJtyWlxo4dqwULFjgft2vXTpK0cuVKde3a1V2nBeCrsueEePBBz8ybA3ir5cvNFSgl6d13XSc2t9nMYXSpqdKuXVI+q1QVKjFReu016Z13zNX1JHMOt9dek1q2vLT4AcAXkJQCAK/gV3iVkpk/f74Mw8i1kZACkMvWrdJ//2tecI8ebXU0gHUOH5b69TPvP/GEdN99rs/bbNJdd5n333uv+Mc/ccKcM+ryy805qtLTpRtvlFavNv8PkpACUF6QlAIAr+C2pBQAFFl2r5B77pFKuHon4POyssyegseOSa1bS2++mXe9/v3N20WLpKLOoXXihPTKK2YvxEmTpLNnpauukuLjze2qq0rnNQCAryApBQBegaQUAGslJUkffGDez57oHCiPXn3VXPWuUiXz/0RISN71brhBqltXOnlS+uyzgo+5e7c5YXlkpDlPVUqKmfD69FNp1SqzlxQAlEckpQDAK5CUAmCtDz80e4h07Cj9NfccUO58++2FCf7ffltq2jT/uv7+F4b45Zi70Skz01w0oFcv8zgzZpg9o9q2NXtXbd4s3X67ORQQAMorklIA4BXcNtE5ABTJwoXmbd++1sYBWOX4cfPz73CYQ/MeeqjwfWJjpYkTpS+/lObONSdD371bWrNG+t//zB6I2W67TXr2WalrVxJRAJCNpBQAeAWfSEqF5DeEAYBv27fPHEJks0kPPGB1NCgDfK69MAxpwABzgvOmTaXp04u2X/Pm0pVXSuvWSYMG5X6+WjXp/vulYcOYpw0ALhIikZQCAC/hE0kpPz9GGQJl0v/9n3l7441S7drWxoIywefai2nTpM8/l4KCpMWLJbu96PsuWCC99Zb0yy/S3r3mJOadOpn/n266SapQwX1xA4AP85NISgGAl/CJpBRgmUOHpNOnS+dYpTFsprSG3nhLLO+/b94ydA/l0YYN0t/+Zt5/802pTZvi7d+smfTPf5Z+XABQHmRkWB0BAEA+kpRKT0+3OgSURx9+KN13n9VRlH2BgdLdd1sdBcoIn2kvUlLMIasZGdJdd0mDB1sdEQCUG+kSPaUAwEv4RFIqg18yYIWffzZvAwPNJdpLm2GU/jHddVx3xWqzSUOGSJUru+f4KHd8or0wDOmJJ8whd/XrS3PmMAE5AHhQhkRSCgC8hE8kpQBLPf44Q2QAlJ5588z51Pz9zdsqVayOCADKH5JSAOAVfGxGWAAAfNj27dLQoeb9V16RunSxNh4AKK9ISgGAVyApBQAoNTNnzlTr1q0VFhamsLAwRUdH68svvyxwn6lTp6pp06YKCQlRZGSknnnmGZ0/f95DEXvQuXPS/febtzEx0siRVkcEAOUXSSkA8AoM3wMAlJp69epp4sSJaty4sQzD0IIFC3TnnXdq8+bNatmyZa76Cxcu1KhRozR37lx16dJFu3bt0oABA2Sz2fTmm29a8ArcaMQIaetWqWZN6b33JD9+FwIAy5CUAgCvQFIKAFBqevbs6fJ4woQJmjlzptasWZNnUmrVqlW6+uqr1bdvX0lSVFSU+vTpo7Vr13okXo/56CPp7bfN+//5j1SrlrXxAEB5R1IKALwCP9MCANwiKytLixYtUmpqqqKjo/Os06VLF23cuFHr1q2TJP3222/64osv1KNHD0+G6l7790uPPGLeHzVKuvlmS8MBAIikFAB4CZ/oKRUSEmJ1CACAItq6dauio6N1/vx52e12LV26VC1atMizbt++fXX8+HFdc801MgxDmZmZeuKJJ/TCCy8UeI60tDSlpaU5H6ekpEjywvYiI0Pq00dKTpaio6WXX7Y6IgAo90IkklIA4CV8oqeUH/NuAIDPaNq0qbZs2aK1a9dq8ODB6t+/v3799dc86yYkJOi1117TW2+9pU2bNmnJkiX6/PPP9corrxR4jri4OIWHhzu3yMhISV7YXrz0krRmjVS5srRwoVShgtURAUC55yeRlAIAL+ETPaUAAL4jMDBQjRo1kiR16NBB69ev17Rp0/TOO+/kqvvSSy+pX79+euSv4W1XXHGFUlNT9dhjj+nFF1/MN8k0evRojRgxwvk4JSXFmZjyGsuXS5MmmffnzJGioiwNBwCQA0kpAPAKPpGUSk9PtzoEAEAJORwOl6F2OZ09ezZX4snf31+SZBhGvscMCgpSUFBQrnKvaS+OHJH69TPvP/mkdPfd1sYDAHBKl0hKAYCX8ImkVEZGhtUhAACKYPTo0br11ltVv359nT59WgsXLlRCQoKWL18uSYqNjVXdunUVFxcnyVyt780331S7du3UuXNn7dmzRy+99JJ69uzpTE4Vh1e0Fw6HmZA6dkxq3VqaMsXqiAAAOWRI5px/AADL+URSCgDgG5KSkhQbG6sjR44oPDxcrVu31vLly3XzXyvOHThwwKVn1JgxY2Sz2TRmzBgdOnRINWrUUM+ePTVhwgSrXsKlmzhRio+XKlaUFi+WgoOtjggAcDF6SgGAVyApBQAoNXPmzCnw+YSEBJfHAQEBGjdunMaNG+fGqDzoxx+lsWPN+zNmSM2aWRsPACBvJKUAwCt42TJFAAD4qD//lPr2lbKypAcflPr3tzoiAEB+SEoBgFcgKQUAwKXKnkfqwAGpcWNp5kzJZrM6KgBAfkhKAYBXICkFAMCleu016csvzfmjPvpICg21OiIAQEFISgGAVyApBQDApfj66wvzSM2caa64BwDwbiSlAMAr+ERSKpiViwAAReDx9uLQIXMeKcOQBg2SBgzw7PkBAMUWLJGUAgAv4RNJKX9/f6tDAAD4AI+2FxkZ0v33S8eOSW3aSP/6l+fODQAoMX+JpBQAeAmfSEoBAOB1nn1W+vFHKSzMnEcqJMTqiAAARUVSCgC8gk8kpdLT060OAQDgAzzWXsyceaFn1Pz5UqNGnjkvAOCSpUskpQDAS/hEUiojI8PqEAAAPsAj7cVXX0nDhpn3J0yQ7rrL/ecEAJSaDImkFAB4CZ9ISgEA4BV+/VXq3VvKypL69ZNGj7Y6IgBASZCUAgCvQFIKAICiOHZMuv12KSVFuuYaafZsyWazOioAQEkwEgMAvAJJKQAACpOWZg7T27dPathQWrpUCgqyOioAQEnRUwoAvAJJKQAACmIY0qOPmivthYdLn30mVa9udVQAgEtBUgoAvAJJKQAAChIXJ733nuTvL334odS8udURAQAuFUkpAPAKJKUAAMjPkiXSiy+a96dPl26+2dp4AAClg6QUAHgFn0hKBQcHWx0CAMAHlGp78euvUv/+5v2nnpKeeKL0jg0AsEywRFIKALyETySl/P39rQ4BAOADSq29SE42JzY/c0bq2lWaMqV0jgsAsJy/RFIKALyETySlAADwqMcfl3btkurVkxYvlgICrI4IAFCaSEoBgFfwiaRURkaG1SEAAHxAqbQXy5ebiSh/f+njj6WaNS/9mAAAr5EhkZQCAC/hE0mp9PR0q0MAAPiAS24vzp+Xhg417w8bJl155aUHBQDwKukSSSkA8BI+kZQCAMAjXn9d2rNHql1bGj/e6mgAAO5CUgoAvAJJKQAAJOn4cSkuzrz/j39IYWHWxgMAcB+SUgDgFUhKAQAgSfPnm8P32reX7rvP6mgAAO5EUgoAvAJJKQAAHA7pnXfM+088Idls1sYDAOXAjBkzFBUVpeDgYHXu3Fnr1q0rsP6pU6c0ZMgQ1a5dW0FBQWrSpIm++OKLkp2chZQAwCu4LSm1f/9+DRo0SA0aNFBISIguv/xyjRs3jknLAQDeZ+VKcy6p0FCpTx+rowGAMm/x4sUaMWKExo0bp02bNqlNmzbq3r27kpKS8qyfnp6um2++Wfv379dHH32knTt3avbs2apbt27JAqCnFAB4hQB3HXjHjh1yOBx655131KhRI23btk2PPvqoUlNTNXnyZHedFgCA4svuJfXQQ5Ldbm0sAFAOvPnmm3r00Uc1cOBASdLbb7+tzz//XHPnztWoUaNy1Z87d65OnDihVatWqUKFCpKkqKiokgdAUgoAvILbekrdcsstmjdvnrp166aGDRvqjjvu0HPPPaclS5YU+1hBQUFuiBAAUNpmzpyp1q1bKywsTGFhYYqOjtaXX35Z4D6lORyjRO3F0aPS0qXm/ccfL9F5AQBFl56ero0bNyomJsZZ5ufnp5iYGK1evTrPff73v/8pOjpaQ4YMUUREhFq1aqXXXntNWVlZ+Z4nLS1NKSkpLpskBUkkpQDAS7itp1RekpOTVbVq1WLvFxDg0TABACVUr149TZw4UY0bN5ZhGFqwYIHuvPNObd68WS1btsxVP3s4Rs2aNfXRRx+pbt26+v3331W5cuUSnb9E7cXHH5sXJ1deKbVpU6LzAgCK7vjx48rKylJERIRLeUREhHbs2JHnPr/99pu++eYbPfjgg/riiy+0Z88ePfnkk8rIyNC4cePy3CcuLk7jx4/PVR4gkZQCAC/hsWzPnj179K9//avAoXtpaWlKS0tzPs7+NQMA4Bt69uzp8njChAmaOXOm1qxZk2dSqtSHY5TE8uXmba9enj0vAKDIHA6HatasqVmzZsnf318dOnTQoUOH9MYbb+SblBo9erRGjBjhfJySkqLIyEjzAUkpAPAKxR6+N2rUKNlstgK3i3/hOHTokG655Rb17t1bjz76aL7HjouLU3h4uHPLbjQyWB0DAHxOVlaWFi1apNTUVEVHR+dZpyTDMaT8h2QUu71IT5e++ca837178fYFAJRI9erV5e/vr6NHj7qUHz16VLVq1cpzn9q1a6tJkyby9/d3ljVv3lyJiYn5LqQUFBTkHE6evUlShiRlZUmGUSqvBwBQcsVOSj377LPavn17gVvDhg2d9Q8fPqwbbrhBXbp00axZswo89ujRo5WcnOzcDh48KEms2AcAPmTr1q2y2+0KCgrSE088oaVLl6pFixZ51v3tt9/00UcfKSsrS1988YVeeuklTZkyRa+++mqB58jvR4xitxerV0tnzkg1akht2xZvXwBAiQQGBqpDhw6Kj493ljkcDsXHx+f7I8bVV1+tPXv2yOFwOMt27dql2rVrKzAwsFjnd7YUhfwAAgBwv2IP36tRo4Zq1KhRpLqHDh3SDTfcoA4dOmjevHny8ys4BxYUFMSk5gDg45o2baotW7YoOTlZH330kfr3769vv/02z8RUSYZjSIUMySiO7KF7N98sFdJGAQBKz4gRI9S/f3917NhRV155paZOnarU1FTnanyxsbGqW7eu4uLiJEmDBw/W9OnT9fTTT2vYsGHavXu3XnvtNT311FMlDyIzU2LuWgCwlNu+hQ8dOqSuXbvqsssu0+TJk3Xs2DHnc/l1ywUA+L7AwEA1atRIktShQwetX79e06ZN0zvvvJOrbu3atVWhQoV8h2Pk9+t3qf2IkZ2UYugeAHjU/fffr2PHjmns2LFKTExU27ZttWzZMufk5wcOHHD5QTsyMlLLly/XM888o9atW6tu3bp6+umnNXLkyJIHwbxSAGA5tyWlVqxYoT179mjPnj2qV6+ey3MG47cBoNxwOBwui1jkdPXVV2vhwoVyOBzOi4+SDscotqQkadMm8363bu49FwAgl6FDh2ro0KF5PpeQkJCrLDo6WmvWrCm9AEhKAYDl3DZWYcCAATIMI88NAFA2jR49Wt99953279+vrVu3avTo0UpISNCDDz4oyRyOMXr0aGf9wYMH68SJE3r66ae1a9cuff7553rttdc0ZMgQ9wf79dfmbZs2Ej14AaD8ISkFAJZjEDUAoNQkJSUpNjZWR44cUXh4uFq3bq3ly5fr5ptvluSh4RhFlZ2UYugeAJRPJKUAwHIkpQAApWbOnDkFPu+R4RhF9cMP5m3Xrp4/NwDAOgEBZkIqI8PqSACg3POJpYZYkQ8AUBRFbi+OHZN27zbvX3WV+wICAHidoOwV9+gpBQCW84mkVABLtQIAiqDI7cWqVeZtixZSlSruCwgA4HUCKlQw75CUAgDL+URSCgCAUpWdlLr6amvjAAB4Hj2lAMBr+ERSKpMGAwBQBEVuL3780bzt0sV9wQAAvFJm9oIbXGMAgOV8IimVlpZmdQgAAB9QpPYiLU3asMG8T1IKAMqdNHpKAYDX8ImkFAAApWbzZjMxVb261Lix1dEAADyNpBQAeA2SUgCA8iXn0D2bzdpYAACeR1IKALwGSSkAQPmSPck5Q/cAoHzy9zdvSUoBgOVISgEAyg/DICkFAOUdPaUAwGuQlAIAlB8HDkiJieYFSceOVkcDALACPaUAwGuQlAIAlB9r1pi3bdtKISGWhgIAsAg9pQDAa/hEUiowMNDqEAAAPqDQ9iI7KXXVVe4PBgDglQIrVDDvkJQCAMv5RFKqQnbDAQBAAQptL0hKAUC552wrMjKsDQQA4BtJKQAALllamrRpk3mfpBQAlF8M3wMAr+ETSalMGgwAQBEU2F5s2SKlp0vVq0sNG3osJgCAd8n0++sSiGsMALCcTySl0tLSrA4BAOADCmwvcg7ds9k8ExAAwOukkZQCAK/hE0kpAAAuGfNJAQAkhu8BgBchKQUAKB9ISgEAJMnf37wlKQUAliMpBQAo+xITpf37zWF7nTpZHQ0AwErZq++RlAIAywVYHQD+snmzNHCglJJidSTIduKE1REAKC1r15q3LVtKYWHWxgIAsBbD9wDAa5CU8haffCL99JPVUSAvTZpYHQGAS8XQPQBANpJSAOA1SEp5C8Mwb++9V3ruOWtjwQV2u9SihdVRALhUJKUAANmYUwoAvIZPJKUCAwOtDsFzatWSOne2OgoA8El5theZmdL69eZ9klIAUO452wqSUgBgOZ+Y6LxC9mSEAAAUIM/24pdfpNRUcy6p5s09HxQAwKtUICkFAF7DJ5JSAACUWPbQvSuvlPxo9gCg3GNOKQDwGj7x13lWVpbVIQAAfECe7QXzSQEAcsiy2cw7GRnWBgIA8I2k1Pnz560OAQBQBDNnzlTr1q0VFhamsLAwRUdH68svvyzSvosWLZLNZlOvXr1KfP482wuSUgCAHJwtBUkpALCcTySlAAC+oV69epo4caI2btyoDRs26MYbb9Sdd96pX375pcD99u/fr+eee07XXntt6QZ08qS0Y4d5n0UkAACSFBpq3p4+bW0cAACSUgCA0tOzZ0/16NFDjRs3VpMmTTRhwgTZ7Xatye6tlIesrCw9+OCDGj9+vBo2bFi6Aa1bZ942aiRVr166xwYA+KbKlc3bEycsDQMAQFIKAOAmWVlZWrRokVJTUxUdHZ1vvZdfflk1a9bUoEGDinzstLQ0paSkuGx5yk5K0UsKAJCNpBQAeI0AqwMAAJQtW7duVXR0tM6fPy+73a6lS5eqRYsWedb94YcfNGfOHG3ZsqVY54iLi9P48eMLr7htm3nbpk2xjg8AKMOyk1InT1oaBgCAnlIAgFLWtGlTbdmyRWvXrtXgwYPVv39//frrr7nqnT59Wv369dPs2bNVvZhD60aPHq3k5GTndvDgwbwrZs9l1bJlcV8GAKCsoqcUAHgNekoBAEpVYGCgGjVqJEnq0KGD1q9fr2nTpumdd95xqbd3717t379fPXv2dJY5HA5JUkBAgHbu3KnLL788z3MEBQUpKCio4EDS06WdO837rVqV8NUAAMocklIA4DV8IilVoUIFq0MAAJSQw+FQWlparvJmzZpp69atLmVjxozR6dOnNW3aNEVGRhb7XC7txe7dUmamucpSCY4FACibKtSoYd45fVrKyJC41gAAy/hEUiowMNDqEAAARTB69Gjdeuutql+/vk6fPq2FCxcqISFBy5cvlyTFxsaqbt26iouLU3BwsFpd1IOp8l+/Xl9cXlQu7UX20L0WLSSbrUTHAwCUPYHZSSlJOnVKyvkYAOBRPpGUAgD4hqSkJMXGxurIkSMKDw9X69attXz5ct18882SpAMHDsjPz0PTGWYnpRi6BwDIyd9fCg+XkpPNIXwkpQDAMj6RlMrKyrI6BABAEcyZM6fA5xMSEgp8fv78+Zd0fpf2InvlPSY5BwDkkJWVJVWtaialWIEPACzlE6vvnT9/3uoQAAA+wKW9YOU9AEAezp8/L1WpYj5gsnMAsJRPJKUAACiW8+fNic4lhu8BAHKrWtW8JSkFAJYiKQUAKHt27pQcDnPZ79q1rY4GAOBtspNSDN8DAEuRlAIAlD05Jzln5T0AwMUYvgcAXoGkFACg7GGScwBAQRi+BwBegaQUAKDsYZJzAEBBSEoBgFcgKQUAKHt+/928bdTI2jgAAN4pe/gec0oBgKXcmpS64447VL9+fQUHB6t27drq16+fDh8+XOzjVKhQwQ3RAQDKGmd7ceSIecsk5wCAi1SoUIGeUgDgJdyalLrhhhv0wQcfaOfOnfr444+1d+9e3XvvvcU+TmBgoBuiAwCUNYGBgVJmpnTsmFlAUgoAvNqMGTMUFRWl4OBgde7cWevWrSvSfosWLZLNZlOvXr2Kfc7AwECSUgDgJdyalHrmmWd01VVX6bLLLlOXLl00atQorVmzRhkZGe48LQCgPDt6VDIMyd9fqlHD6mgAAPlYvHixRowYoXHjxmnTpk1q06aNunfvrqSkpAL3279/v5577jlde+21JT85w/cAwCt4bE6pEydO6P3331eXLl2KPRzP4XC4KSoAQFnicDikxETzQUSE5MfUiQDgrd588009+uijGjhwoFq0aKG3335bFStW1Ny5c/PdJysrSw8++KDGjx+vhg0blui8DofDtaeUYZToOACAS+f2v9ZHjhypSpUqqVq1ajpw4ID++9//5ls3LS1NKSkpLpsknTt3zt1hAgDKgHPnzl2YT6pWLWuDAQDkKz09XRs3blRMTIyzzM/PTzExMVq9enW++7388suqWbOmBg0aVOJznzt37kJSKjNTOnOmxMcCAFyaYielRo0aJZvNVuC2Y8cOZ/3nn39emzdv1ldffSV/f3/FxsbKyOfXiLi4OIWHhzu3yMjIkr8yAED5xCTnAOD1jh8/rqysLEVERLiUR0REKDG7x+tFfvjhB82ZM0ezZ88u0jny+8FbkhQSImXPW8sQPgCwTEBxd3j22Wc1YMCAAuvk7EpbvXp1Va9eXU2aNFHz5s0VGRmpNWvWKDo6Otd+o0eP1ogRI5yPU1JSSEwBAIqHpBQAlDmnT59Wv379NHv2bFWvXr1I+8TFxWn8+PF5P2mzmb2lEhPNIXz165ditACAoip2UqpGjRqqUcKJY7PnhkpLS8vz+aCgIAUFBZXo2AAASLowpxRJKQDwWtWrV5e/v7+OHj3qUn706FHVymP49d69e7V//3717NnTWZZ9bREQEKCdO3fq8ssvd9mn0B+8cyalAACWKHZSqqjWrl2r9evX65prrlGVKlW0d+9evfTSS7r88svz7CUFAECpYE4pAPB6gYGB6tChg+Lj49WrVy9JZpIpPj5eQ4cOzVW/WbNm2rp1q0vZmDFjdPr0aU2bNi3P0RWF/uCdPa8Uw/cAwDJuS0pVrFhRS5Ys0bhx45SamqratWvrlltu0ZgxY+gNBQBwH4bvAYBPGDFihPr376+OHTvqyiuv1NSpU5WamqqBAwdKkmJjY1W3bl3FxcUpODhYrVq1ctm/cuXKkpSrvMiqVDFv6SkFAJZxW1Lqiiuu0DfffOOuwwMAkDeSUgDgE+6//34dO3ZMY8eOVWJiotq2batly5Y5Jz8/cOCA/PzcuFh4dk8pklIAYBm3JaVKU4UKFawOAQDgAyoEBDCnFAD4kKFDh+Y5XE+SEhISCtx3/vz5JTqn89qC4XsAYDk3/vRQegKzl2sFAKAAgampUnq6+YA5pQAAeXBeWzB8DwAs5xNJKQAAiiR7FacqVSTmLwQAFIThewBgOZ9ISmUv9woAQEEchw+bdxi6BwDIh/PagqQUAFjOJ5JS586dszoEAIAPOPfHH+YdklIAgHw4ry0YvgcAlvOJpBQAAEWSlGTekpQCABSmbl3z9uBBa+MAgHKMpBQAoOzITkoxyTkAoDBRUebtiRNSSoqloQBAeUVSCgBQambOnKnWrVsrLCxMYWFhio6O1pdffplv/dmzZ+vaa69VlSpVVKVKFcXExGjdunUlDyB7onN6SgEAChMaKlWrZt7fv9/SUACgvCIpBQAoNfXq1dPEiRO1ceNGbdiwQTfeeKPuvPNO/fLLL3nWT0hIUJ8+fbRy5UqtXr1akZGR6tatmw4dOlSyAEhKAQCKI7u31L59loYBAOUVSSkAQKnp2bOnevToocaNG6tJkyaaMGGC7Ha71qxZk2f9999/X08++aTatm2rZs2a6d1335XD4VB8fHzJAmBOKQBAcTRoYN7SUwoALBFgdQAAgLIpKytLH374oVJTUxUdHV2kfc6ePauMjAxVzV6mu7iYUwoAUBzZSSl6SgGAJXwiKRUQ4BNhAgAkbd26VdHR0Tp//rzsdruWLl2qFi1aFGnfkSNHqk6dOoqJiSmwXlpamtLS0pyPU/6aoDbgzBmzICKiZMEDAMo8l2uL7OF79JQCAEv4xPC9oKAgq0MAABRR06ZNtWXLFq1du1aDBw9W//799euvvxa638SJE7Vo0SItXbpUwcHBBdaNi4tTeHi4c4uMjJQkOVuLSpUu8VUAAMoql2sLekoBgKV8IikFAPAdgYGBatSokTp06KC4uDi1adNG06ZNK3CfyZMna+LEifrqq6/UunXrQs8xevRoJScnO7eDBw+6VqhQ4VJeAgCgvMg50blhWBoKAJRHPjEuzqCBAACf5XA4XIbaXez111/XhAkTtHz5cnXs2LFIxwwKCsqzF60hmQkpm62E0QIAyjqXa4vspNTp09LJk1JJ5zQEAJSITySlzp49q/DwcKvDAAAUYvTo0br11ltVv359nT59WgsXLlRCQoKWL18uSYqNjVXdunUVFxcnSZo0aZLGjh2rhQsXKioqSomJiZIku90uu91e7POflRQeGFhqrwcAUPa4XFuEhJjzEB49avaWIikFAB7F8D0AQKlJSkpSbGysmjZtqptuuknr16/X8uXLdfPNN0uSDhw4oCNHjjjrz5w5U+np6br33ntVu3Zt5zZ58uSSB0FSCgBQHNnzSjHZOQB4nE/0lAIA+IY5c+YU+HxCQoLL4/3uuAAgKQUAKI4GDaQ1a5jsHAAsQE8pAEDZQlIKAFAc2fNK0VMKADyOpBQAoGwhKQUAKI7s4Xv0lAIAjyMpBQAoW0hKAQCKI7unFEkpAPA4klIAgLIlKMjqCAAAviTnROeGYWkoAFDe+ERSKiCA+dgBAIULkOgpBQAoUK5ri/r1JZtNOndOOnrUmqAAoJzyiaRUEL96AwCKIEgiKQUAKFCua4vAwAu9pbZv93xAAFCO+URSCgCAIiMpBQAorpYtzdtffrE2DgAoZ3wiKWUwthsAUASGRFIKAFCgPK8tspNSv/7q2WAAoJzziaTU2bNnrQ4BAOADzkokpQAABcrz2qJFC/OWnlIA4FE+kZQCAKDISEoBAIor5/A9RmkAgMeQlAIAlC0kpQAAxdWsmbkC359/SklJVkcDAOUGSSkAQNlCUgoAUFwVK0oNG5r3GcIHAB5DUgoAULaQlAIAlASTnQOAx5GUAgCULSSlAAAlwWTnAOBxJKUAAGULSSkAQEnknOwcAOARPpGU8vf3tzoEAIAP8JdISgEACpTvtQUr8AGAx/lEUio4ONjqEAAAPiBYkoKCrA4DAODF8r22aNZM8vOTTpyQjh71bFAAUE75RFIKAIAio6cUAKAkQkIurMDHZOcA4BEkpQAAZQtJKQBASTGvFAB4lE8kpVJTU60OAQDgA1IlklIAgAIVeG3BCnwA4FE+kZQCAKDISEoBAEqKnlIA4FEkpQAAZQtJKQBASbECHwB4FEkpAEDZQlIKAFBS2SvwnTwpJSZaHQ0AlHkkpQAAZQtJKQBASQUHS5dfbt5nBT4AcDuSUgCAsoWkFADgUjDZOQB4DEkpAEDZQlIKAHApmOwcADzGJ5JS/v7+VocAACiCmTNnqnXr1goLC1NYWJiio6P15ZdfFrjPhx9+qGbNmik4OFhXXHGFvvjiixKf318iKQUAKFCh1xYkpQDAY3wiKRUcHGx1CACAIqhXr54mTpyojRs3asOGDbrxxht155136pd8/rBftWqV+vTpo0GDBmnz5s3q1auXevXqpW3btpXo/MESSSkAQIEKvbZgBT4A8BiPJKXS0tLUtm1b2Ww2bdmyxROnBABYoGfPnurRo4caN26sJk2aaMKECbLb7VqzZk2e9adNm6ZbbrlFzz//vJo3b65XXnlF7du31/Tp00seRFBQyfcFAKBpU3MFvlOnWIEPANzMI0mpv/3tb6pTp44nTgUA8BJZWVlatGiRUlNTFR0dnWed1atXKyYmxqWse/fuWr16dclPTE8pAMClyLkCH0P4AMCt3J6U+vLLL/XVV19p8uTJJT5GampqKUYEAHCnrVu3ym63KygoSE888YSWLl2qFtkrGV0kMTFRERERLmURERFKLOSX6bS0NKWkpLhskpQqkZQCABSoSNcWzCsFAB7h1qTU0aNH9eijj+q9995TxYoV3XkqAICXaNq0qbZs2aK1a9dq8ODB6t+/v3799ddSPUdcXJzCw8OdW2Rk5IUnSUoBAC4VSSkA8Ai3JaUMw9CAAQP0xBNPqGPHjkXaJ79fvgEAviMwMFCNGjVShw4dFBcXpzZt2mjatGl51q1Vq5aOHj3qUnb06FHVqlWrwHOMHj1aycnJzu3gwYM5A7jk1wAA8IwZM2YoKipKwcHB6ty5s9atW5dv3dmzZ+vaa69VlSpVVKVKFcXExBRY/5K0amXekpQCALcqdlJq1KhRstlsBW47duzQv/71L50+fVqjR48u8rEL/OUbAOCTHA6H0tLS8nwuOjpa8fHxLmUrVqzIdw6qbEFBQQoLC3PZnEhKAYBPWLx4sUaMGKFx48Zp06ZNatOmjbp3766kpKQ86yckJKhPnz5auXKlVq9ercjISHXr1k2HDh0q/eCyk1LbtrECHwC4kc0wivcte+zYMf35558F1mnYsKHuu+8+ffrpp7LZbM7yrKws+fv768EHH9SCBQty7ZeWluZy4ZKSkqLIyEgdPnxYtWvXLk6YvmfsWOmVV6ShQ6V//cvqaACUQykpKQoPD1dycrJrkqcYRo8erVtvvVX169fX6dOntXDhQk2aNEnLly/XzTffrNjYWNWtW1dxcXGSpFWrVun666/XxIkTddttt2nRokV67bXXtGnTJrXKviAoRuyHJdU+flyqVq1E8QMAClca7YUkde7cWZ06dXKuuOpwOBQZGalhw4Zp1KhRhe6flZWlKlWqaPr06YqNjS1y3EW6tkhPlypVkjIzpd9/l+rXL9JrAgCYitpWBBT3wDVq1FCNGjUKrffPf/5Tr776qvPx4cOH1b17dy1evFidO3fOc5+goCAFsZQ3APispKQkxcbG6siRIwoPD1fr1q2dCSlJOnDggPz8LnTS7dKlixYuXKgxY8bohRdeUOPGjfXJJ58UKyGVCz2lAMDrpaena+PGjS6jKvz8/BQTE1PkFVjPnj2rjIwMVa1aNc/n8/rBu8gCA6WmTc3he9u2kZQCADcpdlKqqOpf9MVtt9slSZdffrnq1avnrtMCACw0Z86cAp9PSEjIVda7d2/17t279IIgKQUAXu/48ePKysrKcwXWHTt2FOkYI0eOVJ06dRQTE5Pn83FxcRo/fnzJg2zV6kJSqkePkh8HAJAvt66+V1py/qoOAEB+/CSpQgWrwwAAuNnEiRO1aNEiLV26VMHBwXnWyW9RjCJfW+ScVwoA4BZu6yl1saioKBVz+iqnkJCQUo4GAFAWhfj7S/yQAQBer3r16vL39y/RCqyTJ0/WxIkT9fXXX6t169b51stvapAiX1tccYV5S1IKANyGv9wBAGUHQ/cAwCcEBgaqQ4cOLiuwOhwOxcfHF7gC6+uvv65XXnlFy5YtU8eOHd0bZHZPqV9/lbKy3HsuACinSEoBAMoOFssAAJ8xYsQIzZ49WwsWLND27ds1ePBgpaamauDAgZKk2NhYl4nQJ02apJdeeklz585VVFSUEhMTlZiYqDNnzrgnwAYNpJAQKS1N2rPHPecAgHLOJ5JSqampVocAAPABqQEeG5UOALhE999/vyZPnqyxY8eqbdu22rJli5YtW+ac/PzAgQM6cuSIs/7MmTOVnp6ue++9V7Vr13ZukydPLtZ5i3xt4ecntWxp3mcIHwC4BX+9AwDKDiY5BwCfMnToUA0dOjTP5y5esXX//v3uD+hirVpJGzaYSal77vH8+QGgjPOJnlIAABQJc0oBAEoTK/ABgFuRlAIAlB30lAIAlCZW4AMAtyIpBQAoO0hKAQBKU3ZSatcuyV0TqgNAOUZSCgBQdjB8DwBQmmrXlurVkxwOc24pAECpIikFACg76CkFAChtV11l3q5da20cAFAG+URSys/PJ8IEAFjMj55SAIBCFPvaIjsptWZN6QcDAOWcT2R7QkJCrA4BAOADaC8AAIUpdlvRubN5u2aNZBilHxAAlGM+kZQCAKBIGL4HACht7dtLAQFSYqJ08KDV0QBAmUJSCgBQdjB8DwBQ2ipWlFq3Nu8zrxQAlCqfSEqdPXvW6hAAAD7gLHMQAgAKUaJrC+aVAgC38Im/3g3GbgMAisCgpxQAoBAlurZgBT4AcAufSEoBAFAkzCkFAHCH7MnON26UMjKsjQUAyhCSUgCAsoOeUgAAd2jcWKpSRTp/3kxMAQBKBUkpAEDZQVIKAOAONpvUrZt5/4MPrI0FAMoQklIAgLIjIMDqCAAAZVXfvubtokVSVpa1sQBAGUFSCgBQdtBTCgDgLrfcIlWtKh05Iq1caXU0AFAm+ERSymazWR0CAMAH2EhKAQAKUeJri8BAqXdv8/7775deQABQjvlEUqpixYpWhwAA8AEVK1WyOgQAgJe7pGuLBx80bz/+WDp3rnQCAoByzCeSUgAAFAk9pQAA7nT11VL9+tLp09L//md1NADg80hKAQDKjgoVrI4AAFCW+flJsbHm/VGjpDNnrI0HAHycTySlzp49a3UIAIAiiIuLU6dOnRQaGqqaNWuqV69e2rlzZ6H7TZ06VU2bNlVISIgiIyP1zDPP6Pz588U+/1nDKEnYAIBy5JKvLf72N+myy6T9+837AIAS84mklMFFBgD4hG+//VZDhgzRmjVrtGLFCmVkZKhbt25KTU3Nd5+FCxdq1KhRGjdunLZv3645c+Zo8eLFeuGFF4p9foPhewCAQlzytUVoqDRnjnl/5kwpPv7SgwKAcirA6gAAAGXHsmXLXB7Pnz9fNWvW1MaNG3Xdddfluc+qVat09dVXq2/fvpKkqKgo9enTR2vXri1+ACSlAACecNNN0hNPSG+/ba7I99lnUpcuVkcFAD7HJ3pKAQB8U3JysiSpatWq+dbp0qWLNm7cqHXr1kmSfvvtN33xxRfq0aNHvvukpaUpJSXFZZPEnFIAAM95/XWpc2fp5EkpJkb69FOrIwIAn0NSCgDgFg6HQ8OHD9fVV1+tVq1a5Vuvb9++evnll3XNNdeoQoUKuvzyy9W1a9cCh+/FxcUpPDzcuUVGRppP0FMKAOApoaHm0L0ePaRz56Q775ReflnKyrI6MgDwGSSlAABuMWTIEG3btk2LFi0qsF5CQoJee+01vfXWW9q0aZOWLFmizz//XK+88kq++4wePVrJycnO7eDBg+YT9JQCAHhSpUrSJ59Ijz4qGYY0bpyZpDp2zOrIAMAzsrKk06elo0fNBSB+/VXasEH68cci7c6cUgCAUjd06FB99tln+u6771SvXr0C67700kvq16+fHnnkEUnSFVdcodTUVD322GN68cUX5eeX+/eToKAgBQUF5T4YSSkAgKdVqCDNmiVdfbU0eLD01VdSu3bSokXSNddYHR0AXJCRIaWkmNvp065bQWVnzkhnz7pu586Zt+nplxSSTySlbDab1SEAAIrAMAwNGzZMS5cuVUJCgho0aFDoPmfPns2VePL393cerzhseSWqAADIwW3XFv37Sx06mBOf79ghde0qTZkiPfWUxPUMgNKUkSEdP272Tjp2TDpxwpzfLvs25/2cZWfOuDeuihUvbIGB0p49he7iE0mpihUrWh0CAKAIhgwZooULF+q///2vQkNDlZiYKEkKDw9XSEiIJCk2NlZ169ZVXFycJKlnz55688031a5dO3Xu3Fl79uzRSy+9pJ49ezqTU0VVMTS0dF8QAKDMceu1RatW0vr10uOPSwsXSsOHm0NZpk+nNy+AgmVmSomJ0h9/SEeOSElJ5nb06IX72Y9PnLi0c4WEmPPihYWZt9nbxY+zy+x2c7hyzqRTSIjr46Ag1wR8SooUHl5oKD6RlAIA+IaZM2dKkrp27epSPm/ePA0YMECSdODAAZeeUWPGjJHNZtOYMWN06NAh1ahRQz179tSECROKHwATnQMArGa3S//5j9lr6rnnzKF9v/0mLV1qPgeg/HE4zITTb79JBw+aiaeLt8REs15R+flJNWqYW9Wq5lalSuG3YWFSgPekgrwnEgCAzyvKcLuEhASXxwEBARo3bpzGjRt36QGQlAIAeAObTRoxQmrSRHrgAenrr6Wbb5a++MK8KARQ9pw9K+3bZyaecm5795rl588XfoyAAKlOHal2bSkiQqpZ88LtxferVpWKOarAG/lEUurcuXMKCwuzOgwAgJc753CI1gIAUBCPXlvcfrv0zTfSLbdIa9aY80x99ZV5YekpW7eaQwqTkqS0NLOXRM4tPNy8rVpVql7dq3pQAF7p2DFz3rjt280t+/7vvxe8n5+fVL++uUVGSvXquW5165rJpjKQaCoOn/jGcRSnCxsAoNxylLNGHABQfB6/trjySum778yeUj//LF17rdlzqn5995/7f/+TevWSirNwSNWqF3pi1Khx4X52L426dc0L6Fq1SGChbMvKknbtkjZvlrZsubAdO5b/PlWqSA0b5t4aNDD/zzO3XC58iwAAyg4aegCAN2rVSvr+eykmRtq9W7rmGik+Xmrc2H3nTEyUBg0yE1KdOkktWkjBwReWeU9OvrA0fHKydOqUOZ/NiRPmtmNHwcf38zOTVNk9PLJvs+/Xq2f2BgkOdt9rBEqLYZhD7FatMreNG81ehufO5a5rs0mXXSY1b25uzZpduK1e3fOx+ziSUgCAsoM5pQAA3qpRIzMxdfPN0s6d0nXXmYmpFi1K/1wOhzRggLlkfNu25nmDggreJyvLTEYlJZk9QXKu9nXsmLniV2KidOiQdPiwuVLYkSPmtn59/setXVuKisp7u+yywuMC3CE9Xdqw4UISatUq8zN+sUqVpDZtzP9H7dqZty1amKvNoVSQlAIAlB0kpQAA3iwy0nUo3/XXSytWmBe6pemDD6Tly81eSu+/X7TEj7//hZW8CuNwmMmqQ4fMVcNy3mbf/+MPKTX1QuJq9eq8j1WnjpmgatjQ7DnWuLE5QXzjxuZcV2WRw2H2WDt9WsrIMHt6BwaaW/b9gACzRw5Kh2GYvf+++srcvv3W/HzmVKGCuWpmly5S585mEuryy81egXAbklIAgLKD4XsAAG9Xs6a0cqXUrZs5ROiGG8yL5E6dSu8cv/5q3vbr556eWH5+5pxStWqZF/F5MQzpzz+l/fvz31JTzV5Xhw+bPVUuFhFxIVGVM1nVqJFv9lQxDDNJ+OyzZlKvMHklq3JuQUEXbgvaCqtT3GP4yhyeDoeZDP3gA2npUungQdfna9SQrr7aTEJ16WJ+lhlu6nEkpQAAZQe/ZAEAfEHVqubQvVtvNS+ab7pJ+vJL8wK5NFk5NM5mM+fXqV5d6tgx9/M5k1b79kl795rzbe3aZd4ePXph++GH3Mdu2FBq2dKcryt7a9LEe4cDnjhhDqn89NMLZQEBZrIpI8McDnmxjAxzu7hHj9X8/aWQEDMxWNLbvMpCQyW7/cJtSf+u27hReu896aOPzJ572YKCzGGz3bqZ2xVX0BvNC5CUAgAAAABPCw83h9j17GkOJbr5ZmnhQnO1vPKgsKRVSoqZnMqZqMq+f/KkmcTau9dcYTCbv7852XTHjubWqZM5H5A39H55+20zIVWhgjRunPTMM2YyJjspYhgXklDp6eaW3/3sLS3NdcurrDjPF1Qnp6ws6cwZc3OnSpUuJKlyJqzyemy3m+//f/5jzqGWLSxMuuMOqXdv8/9YSIh7Y0ax+URSqlKlSlaHAADwAbQXAIDCeFVbERoqffGFdO+9Zk+pu++W3nhDGjGCHhxhYeZwqryGByYlSb/8Im3bduF22zZzFcFffjG3BQvMugEBZo+YTp2kq64yh0tGRXn0pUi6sIrbY49JL76Y+3mb7cKwPG/6jEoXEmY5k1Tnz0tnz5qvq6DbotQ5d87sDZY9z5bDYZ43NdXc8pqAvCAVKpj/px54wOwR5Q1JSeTLJ5JSAAAAAFAmVaxo9vZ56ilp5kzpueektWuld96RqlSxOjrvVLOmud1ww4UywzCHam3ZYq4GmL0dPy5t3mxus2aZdaOipK5dpRtvlGJizBUCPcUXpxrImTALDXXvuQzDTHidOXMhSVXU+6mpZg+5J580J9CHTyApBQAAAABWCgiQZsyQmjY1k1IffiitWSPNni117251dL7BZpPq1TO32283ywxDOnDgQoLq++/N2/37pfnzzU0yh/jdcou5denCar5WstnMIXYhIUVbCRI+z61p2qioKNlsNpdt4sSJxT7OueyujgAAFID2AgBQGK9tK2w26emnzVXoGjUyVwq75RbpttvMoWkoPptNuuwycyjXpEnmv+3Jk+ZQyZEjzV41Npv000/m8zfcIFWrZs7r9fbb5gTsANzK7T2lXn75ZT366KPOx6El6O7nyB5TCgBAAWgvAACF8fq2olMnadMmaexYafp0c86pL74wV+gbPNhcsa9iRauj9F12+4VeUZJ07Ji0YoW0bJk58XxSkvTf/5qbZK7od9115vvSqZO5yl+FCtbFD5Qxbk9KhYaGqlatWu4+DQAAAACUDaGh0j/+YSahXnxRWrJEio83t5AQcx6k664zJ+5u2ZK5py5FjRpS377m5nCYc1ItX24mqX780Vztb9cu6d13zfrBwVLr1uZQy0aNLmz165vH8ve39OUAvsbtSamJEyfqlVdeUf369dW3b18988wzCggo5mm//lqqWtU9AXqL3butjgAAAAAo886ePSvDMGT7a3W79PR0ZWRkKCAgQEFBQc56qampkqSQkBD5/TU5dUZGhtLT0+Xv76/gHCt6Fadu9vmDg4Pl/1cCIzMzU2lpafLz81NIjiXrz9arJ2P+fAVPmiT/2bOl//s/Zf7+u9I+/VR+n34qZ80qVXSuQQM5GjRQUKNGCti0SZKU5XDofGpqruOeO3dODodDQUFBzmuzrKwsnT9/XjabTRVz9MQ6f/68srKyFBgYqAp/9RAqTl2Hw+EcMplz5cO0tDRlZmaqQoUKCvxrDqfi1DUMQ2fPnpUkVaxYMdf7WZy6zvfez09q316pTZtKTz2lkPR0+X33nbR2rdLXrlXGxo0KSE5W0Lp10rp15nv/V4whkvne16ihjIgIpdeoIf/atRVct66ZrNq8WWclGRkZCs7Kcr73pfE5yX4/i1O3KO/9pX5O8ns/L/VzkvP9vNTPSX7/74tT19LviBJ8TvJ6P93xHZH9mgtluNGUKVOMlStXGj/99JMxc+ZMo3LlysYzzzyTb/3z588bycnJzu3gwYOGJOOwOUVd+dieftqdbwkA5Cs5OdmQZCQnJ1sdSrFlx3748GGrQwGAMs9X24vsuCUZSUlJzvJXX33VkGQ88sgjLvUrVqxoSDL27dvnLPvHP/5hSDL69u3rUrd69eqGJGPbtm3OslmzZhmSjDvvvNOl7mWXXWZIMtatW+cs+89//mNIMmJiYlzqtmjRwpBkrFy50ixwOIylb75pSDK6VK1qGHXqOK8jOv712j7LcW3xVa9ehiSjTZs2Lse9/vrrDUnGBx984Cz74YcfDElGo0aNXOr26NHDkGTMmzfPWbZ582ZDklGnTh2Xuvfee68hyZg+fbqzbNeuXYYkIzw83KVu//79DUnG66+/7iz7448/DElGQECAS90nn3zSkGSMGzfOWXby5Enn+5menu4sf+655wxJxnPPPecsS09Pd9Y9efKks3zcuHGGJOPJJ590OV9AQIAhyfjjjz+cZa+//rohyeh/112G8eGHhhEXZxiDBhnh/v6GJGNXjn/36X+d696LrvXq/FW+efNm53HnzZtnSDJ69OjhEkOjRo0MScYPP/zgLPvggw8MScb111/vUrdNmzaGJOOrr75yln322WeGJKNjx44udbt06WJIMpYuXeosW7lypSHJaNGihUvdmJgYQ5Lxn//8x1m2bt06Q5Jx2WWXudS98847DUnGrFmznGXbtm0zJBnVq1d3qdu3b19DkvGPf/zDWbZv3z5DklGxYkWXuo888oghyXj11VedZUlJSc73M6enn37akGS88MILzrIzZ8446545c8ZZ/sILLxiSjKcvuv72+e8IwzCWLl1qfkd06eJSt2PHjuZ3xGefOcu++uort35HFKWtKHZPqVGjRmnSpEkF1tm+fbuaNWumESNGOMtat26twMBAPf7444qLi3PJMGaLi4vT+PHjcx+wdevysQJCpUpSbKzVUQAAAADwRjab1KCBeb9ZM3N42Zkz5oTcd98t7dljrjwXECCdPSt16yZ98omlIZc5YWHmxOnZPvpISk6Wfv3VHEaZmGiumvjWW+b8UzfdJB0/Lv35p7RypZSWZl3sgBeyGYZhFGeHY8eO6c8//yywTsOGDZ3d5XL65Zdf1KpVK+3YsUNNmzbN9XxaWprScvwnTUlJUWRkpA4fPqzatWsXJ0wAQDGlpKQoPDxcycnJCgsLszqcYsmOnfYCANyvNNuLGTNm6I033lBiYqLatGmjf/3rX7ryyivzrf/hhx/qpZde0v79+9W4cWNNmjRJPXr0KFbce/bsUcOGDRmaw/C9Yr33pfE5yev9ZPgew/fK6ndEUlKS6tSpU2hbUeyk1KV4//33FRsbq+PHj6tKESbj4yIDADyHpBQAoChKq71YvHixYmNj9fbbb6tz586aOnWqPvzwQ+3cuVM1a9bMVX/VqlW67rrrFBcXp9tvv10LFy7UpEmTtGnTJrVq1arIcdNWAID7FbWtcFtSavXq1Vq7dq1uuOEGhYaGavXq1XrmmWd06623asGCBUU6hi9fIAGAr/Hl71xfjh0AfE1pfed27txZnTp10vTp0yWZPSUiIyM1bNgwjRo1Klf9+++/X6mpqfrss8+cZVdddZXatm2rt99+22NxAwAKV9TvXD93BRAUFKRFixbp+uuvV8uWLTVhwgQ988wzmjVrlrtOCQAAAMAHpKena+PGjYqJiXGW+fn5KSYmRqtXr85zn9WrV7vUl6Tu3bvnWx8A4P2KPdF5UbVv315r1qxx1+EBAAAA+Kjjx48rKytLERERLuURERHasWNHnvskJibmWT8xMTHP+nnNVwsA8C5u6ylVms6fP291CAAAH0B7AQDIFhcXp/DwcOcWGRkpibYCALyJTySlsrKyrA4BAOADaC8AwDdUr15d/v7+Onr0qEv50aNHVatWrTz3qVWrVrHqjx49WsnJyc7t4MGDkmgrAMCb+ERSCgAAAEDZERgYqA4dOig+Pt5Z5nA4FB8fr+jo6Dz3iY6OdqkvSStWrMi3flBQkMLCwlw2AIB3cducUgAAAACQnxEjRqh///7q2LGjrrzySk2dOlWpqakaOHCgJCk2NlZ169ZVXFycJOnpp5/W9ddfrylTpui2227TokWLtGHDBhZSAgAfRlIKAAAAgMfdf//9OnbsmMaOHavExES1bdtWy5Ytc05mfuDAAfn5XRjY0aVLFy1cuFBjxozRCy+8oMaNG+uTTz5Rq1atrHoJAIBLZDMMw7A6iPykpKQoPDxchw8fVu3ata0OBwDKtOzv3OTkZJ8b4kB7AQCe46vtBW0FAHhOUdsK5pQCAAAAAACAx3n18L3sTlynT59WpUqVLI4GAMq2lJQUSRe+e30J7QUAeI6vthe0FQDgOUVtK7w6KfXnn39Kkpo2bWpxJABQfpw+fVrh4eFWh1EstBcA4Hm+1l7QVgCA5xXWVnh1Uqpq1aqSzEkOfanBg+ekpKQoMjJSBw8e9Kk5DeA5fEaKzjAMnT59WnXq1LE6lGKjvUBh+C5AYfiMFJ2vthe0FSgM3wMoDJ+RoitqW+HVSans1TbCw8N5w1GgsLAwPiMoEJ+RovHVP9JpL1BUfBegMHxGisYX2wvaChQV3wMoDJ+RoilKW8FE5wAAAAAAAPA4klIAAAAAAADwOK9OSgUFBWncuHEKCgqyOhR4KT4jKAyfkfKB9xmF4TOCwvAZKft4j1EYPiMoDJ+R0mczfG0tVwAAAAAAAPg8r+4pBQAAAAAAgLKJpBQAAAAAAAA8jqQUAAAAAAAAPM6rk1IzZsxQVFSUgoOD1blzZ61bt87qkOAlvvvuO/Xs2VN16tSRzWbTJ598YnVI8DJxcXHq1KmTQkNDVbNmTfXq1Us7d+60Oiy4AW0FCkJ7gYLQVpQvtBfID20FCkN74T5em5RavHixRowYoXHjxmnTpk1q06aNunfvrqSkJKtDgxdITU1VmzZtNGPGDKtDgZf69ttvNWTIEK1Zs0YrVqxQRkaGunXrptTUVKtDQymirUBhaC9QENqK8oP2AgWhrUBhaC/cx2tX3+vcubM6deqk6dOnS5IcDociIyM1bNgwjRo1yuLo4E1sNpuWLl2qXr16WR0KvNixY8dUs2ZNffvtt7ruuuusDgelhLYCxUF7gcLQVpRdtBcoKtoKFAXtRenxyp5S6enp2rhxo2JiYpxlfn5+iomJ0erVqy2MDICvSk5OliRVrVrV4khQWmgrAJQ22oqyifYCQGmjvSg9XpmUOn78uLKyshQREeFSHhERocTERIuiAuCrHA6Hhg8frquvvlqtWrWyOhyUEtoKAKWJtqLsor0AUJpoL0pXgNUBAIC7DRkyRNu2bdMPP/xgdSgAAC9FWwEAKArai9LllUmp6tWry9/fX0ePHnUpP3r0qGrVqmVRVAB80dChQ/XZZ5/pu+++U7169awOB6WItgJAaaGtKNtoLwCUFtqL0ueVw/cCAwPVoUMHxcfHO8scDofi4+MVHR1tYWQAfIVhGBo6dKiWLl2qb775Rg0aNLA6JJQy2goAl4q2onygvQBwqWgv3Mcre0pJ0ogRI9S/f3917NhRV155paZOnarU1FQNHDjQ6tDgBc6cOaM9e/Y4H+/bt09btmxR1apVVb9+fQsjg7cYMmSIFi5cqP/+978KDQ11zhkRHh6ukJAQi6NDaaGtQGFoL1AQ2oryg/YCBaGtQGFoL9zHZhiGYXUQ+Zk+fbreeOMNJSYmqm3btvrnP/+pzp07Wx0WvEBCQoJuuOGGXOX9+/fX/PnzPR8QvI7NZsuzfN68eRowYIBng4Fb0VagILQXKAhtRflCe4H80FagMLQX7uPVSSkAAAAAAACUTV45pxQAAAAAAADKNpJSAAAAAAAA8DiSUgAAAAAAAPA4klIAAAAAAADwOJJSAAAAAAAA8DiSUgAAAAAAAPA4klIAAAAAAADwOJJSAAAAAAAA8DiSUsBfBgwYoF69enn8vPPnz5fNZpPNZtPw4cOd5VFRUZo6dWqB+2bvV7lyZbfGCAAw0VYAAIqC9gIomgCrAwA8wWazFfj8uHHjNG3aNBmG4aGIXIWFhWnnzp2qVKlSsfY7cuSIFi9erHHjxrkpMgAoP2grAABFQXsBlB6SUigXjhw54ry/ePFijR07Vjt37nSW2e122e12K0KTZDZstWrVKvZ+tWrVUnh4uBsiAoDyh7YCAFAUtBdA6WH4HsqFWrVqObfw8HDnF3X2Zrfbc3Wx7dq1q4YNG6bhw4erSpUqioiI0OzZs5WamqqBAwcqNDRUjRo10pdffulyrm3btunWW2+V3W5XRESE+vXrp+PHj5co7rNnz+rhhx9WaGio6tevr1mzZl3KPwMAoAC0FQCAoqC9AEoPSSmgAAsWLFD16tW1bt06DRs2TIMHD1bv3r3VpUsXbdq0Sd26dVO/fv109uxZSdKpU6d04403ql27dtqwYYOWLVumo0eP6r777ivR+adMmaKOHTtq8+bNevLJJzV48GCXX2EAANajrQAAFAXtBZAbSSmgAG3atNGYMWPUuHFjjR49WsHBwapevboeffRRNW7cWGPHjtWff/6pn3/+WZI0ffp0tWvXTq+99pqaNWumdu3aae7cuVq5cqV27dpV7PP36NFDTz75pBo1aqSRI0eqevXqWrlyZWm/TADAJaCtAAAUBe0FkBtzSgEFaN26tfO+v7+/qlWrpiuuuMJZFhERIUlKSkqSJP30009auXJlnmPI9+7dqyZNmpT4/NndgrPPBQDwDrQVAICioL0AciMpBRSgQoUKLo9tNptLWfbKGw6HQ5J05swZ9ezZU5MmTcp1rNq1a5fK+bPPBQDwDrQVAICioL0AciMpBZSi9u3b6+OPP1ZUVJQCAvjvBQDIjbYCAFAUtBcoD5hTCihFQ4YM0YkTJ9SnTx+tX79ee/fu1fLlyzVw4EBlZWVZHR4AwAvQVgAAioL2AuUBSSmgFNWpU0c//vijsrKy1K1bN11xxRUaPny4KleuLD8//rsBAGgrAABFQ3uB8sBmGIZhdRBAeTZ//nwNHz5cp06dsmR/AID3o60AABQF7QV8DelVwAskJyfLbrdr5MiRxdrPbrfriSeecFNUAABvQlsBACgK2gv4EnpKARY7ffq0jh49KkmqXLmyqlevXuR99+zZI8lcUrZBgwZuiQ8AYD3aCgBAUdBewNeQlAIAAAAAAIDHMXwPAAAAAAAAHkdSCgAAAAAAAB5HUgoAAAAAAAAeR1IKAAAAAAAAHkdSCgAAAAAAAB5HUgoAAAAAAAAeR1IKAAAAAAAAHkdSCgAAAAAAAB5HUgoAAAAAAAAe9/+e5fvfCrx71QAAAABJRU5ErkJggg==",
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "\n",
+ "# Plot\n",
+ "plot = pybamm.QuickPlot(\n",
+ " sim.solution,\n",
+ " [\n",
+ " \"Current [A]\",\n",
+ " \"Voltage [V]\",\n",
+ " \"Anode potential [V]\",\n",
+ " ]\n",
+ ")\n",
+ "plot.plot(0)\n",
+ "\n",
+ "# Plot the limits used in the termination events to check they are not surpassed\n",
+ "plot.axes.by_variable(\"Voltage [V]\").axhline(4.2, color=\"k\", linestyle=\":\")\n",
+ "plot.axes.by_variable(\"Anode potential [V]\").axhline(0.02, color=\"k\", linestyle=\":\")\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We can check which events were reached by each step"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Step 0: event: Anode potential cut-off [V] [experiment]\n",
+ "Step 1: event: Anode potential cut-off [V] [experiment]\n",
+ "Step 2: event: Charge voltage cut-off [V] [experiment]\n",
+ "Step 3: event: C-rate cut-off [experiment]\n"
+ ]
+ }
+ ],
+ "source": [
+ "for i, step in enumerate(sim.solution.cycles[0].steps):\n",
+ " print(f\"Step {i}: {step.termination}\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "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] Chang-Hui Chen, Ferran Brosa Planella, Kieran O'Regan, Dominika Gastol, W. Dhammika Widanage, and Emma Kendrick. Development of Experimental Techniques for Parameterization of Multi-scale Lithium-ion Battery Models. Journal of The Electrochemical Society, 167(8):080534, 2020. doi:10.1149/1945-7111/ab9050.\n",
+ "[3] Marc Doyle, Thomas F. Fuller, and John Newman. Modeling of galvanostatic charge and discharge of the lithium/polymer/insertion cell. Journal of the Electrochemical society, 140(6):1526–1533, 1993. doi:10.1149/1.2221597.\n",
+ "[4] 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",
+ "[5] Peyman Mohtat, Suhak Lee, Jason B Siegel, and Anna G Stefanopoulou. Towards better estimability of electrode-specific state of health: decoding the cell expansion. Journal of Power Sources, 427:101–111, 2019.\n",
+ "[6] 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",
+ "[7] 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",
+ "[8] Andrew Weng, Jason B Siegel, and Anna Stefanopoulou. Differential voltage analysis for battery manufacturing process control. arXiv preprint arXiv:2303.07088, 2023.\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "pybamm.print_citations()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "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.11.3"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/docs/source/examples/notebooks/experiments-start-time.ipynb b/docs/source/examples/notebooks/simulations_and_experiments/experiments-start-time.ipynb
similarity index 100%
rename from docs/source/examples/notebooks/experiments-start-time.ipynb
rename to docs/source/examples/notebooks/simulations_and_experiments/experiments-start-time.ipynb
diff --git a/docs/source/examples/notebooks/rpt-experiment.ipynb b/docs/source/examples/notebooks/simulations_and_experiments/rpt-experiment.ipynb
similarity index 100%
rename from docs/source/examples/notebooks/rpt-experiment.ipynb
rename to docs/source/examples/notebooks/simulations_and_experiments/rpt-experiment.ipynb
diff --git a/docs/source/examples/notebooks/simulating-long-experiments.ipynb b/docs/source/examples/notebooks/simulations_and_experiments/simulating-long-experiments.ipynb
similarity index 100%
rename from docs/source/examples/notebooks/simulating-long-experiments.ipynb
rename to docs/source/examples/notebooks/simulations_and_experiments/simulating-long-experiments.ipynb
diff --git a/docs/source/examples/notebooks/simulation-class.ipynb b/docs/source/examples/notebooks/simulations_and_experiments/simulation-class.ipynb
similarity index 96%
rename from docs/source/examples/notebooks/simulation-class.ipynb
rename to docs/source/examples/notebooks/simulations_and_experiments/simulation-class.ipynb
index bb93ec207a..df82fa8175 100644
--- a/docs/source/examples/notebooks/simulation-class.ipynb
+++ b/docs/source/examples/notebooks/simulations_and_experiments/simulation-class.ipynb
@@ -6,7 +6,7 @@
"metadata": {},
"source": [
"# A step-by-step look at the Simulation class\n",
- "The simplest way to solve a model is to use the `Simulation` class. This automatically processes the model (setting of parameters, setting up the mesh and discretisation, etc.) for you, and provides built-in functionality for solving and plotting. Changing things such as parameters in handled by passing options to the `Simulation`, as shown in the [Getting Started](getting_started/tutorial-1-how-to-run-a-model.ipynb) guides, [example notebooks](../index.rst) and [documentation](https://docs.pybamm.org/en/latest/source/api/simulation.html).\n",
+ "The simplest way to solve a model is to use the `Simulation` class. This automatically processes the model (setting of parameters, setting up the mesh and discretisation, etc.) for you, and provides built-in functionality for solving and plotting. Changing things such as parameters in handled by passing options to the `Simulation`, as shown in the [Getting Started](../getting_started/tutorial-1-how-to-run-a-model.ipynb) guides, [example notebooks](../../index.rst) and [documentation](https://docs.pybamm.org/en/latest/source/api/simulation.html).\n",
"\n",
"In this notebook we show how to solve a model using a `Simulation` and compare this to manually handling the different stages of the process, such as setting parameters, ourselves step-by-step."
]
@@ -152,7 +152,7 @@
"metadata": {},
"source": [
"## Processing the model step-by-step\n",
- "One way of gaining more control over the simulation processing is by passing options, as outlined in the [documentation](https://docs.pybamm.org/en/latest/source/api/simulation.html). However, you can also process the model step-by-step yourself. A detailed example of this can be found in the [SPM notebook](./models/SPM.ipynb), but here we outline the basic steps.\n",
+ "One way of gaining more control over the simulation processing is by passing options, as outlined in the [documentation](https://docs.pybamm.org/en/latest/source/api/simulation.html). However, you can also process the model step-by-step yourself. A detailed example of this can be found in the [SPM notebook](../models/SPM.ipynb), but here we outline the basic steps.\n",
"\n",
"First we pick a model"
]
@@ -171,7 +171,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "Next we must set up the geometry. We'll use the default geometry for the SPM. In all of the following steps we will also use the default settings provided by the model. For a look at changing these options, see the [change settings](./change-settings.ipynb) notebook."
+ "Next we must set up the geometry. We'll use the default geometry for the SPM. In all of the following steps we will also use the default settings provided by the model. For a look at changing these options, see the [change settings](../change-settings.ipynb) notebook."
]
},
{
diff --git a/docs/source/user_guide/installation/GNU-linux.rst b/docs/source/user_guide/installation/GNU-linux.rst
index ca95bbe1b5..cf027db587 100644
--- a/docs/source/user_guide/installation/GNU-linux.rst
+++ b/docs/source/user_guide/installation/GNU-linux.rst
@@ -133,7 +133,10 @@ Optional - JaxSolver
~~~~~~~~~~~~~~~~~~~~
Users can install ``jax`` and ``jaxlib`` to use the Jax solver.
-Currently, only GNU/Linux and macOS are supported.
+
+.. note::
+
+ The Jax solver is not supported on Python 3.8. It is supported on Python 3.9, 3.10, and 3.11.
.. code:: bash
diff --git a/docs/source/user_guide/installation/index.rst b/docs/source/user_guide/installation/index.rst
index 2b8b7fe304..65cbad33fb 100644
--- a/docs/source/user_guide/installation/index.rst
+++ b/docs/source/user_guide/installation/index.rst
@@ -217,13 +217,13 @@ Dependency Minimum Version p
Jax dependencies
^^^^^^^^^^^^^^^^^
-Installable with ``pip install "pybamm[jax]"``
+Installable with ``pip install "pybamm[jax]"``, currently supported on Python 3.9-3.11.
========================================================================= ================== ================== =======================
Dependency Minimum Version pip extra Notes
========================================================================= ================== ================== =======================
-`JAX `__ 0.4.8 jax For JAX solvers
-`jaxlib `__ 0.4.7 jax Support library for JAX
+`JAX `__ 0.4.20 jax For the JAX solver
+`jaxlib `__ 0.4.20 jax Support library for JAX
========================================================================= ================== ================== =======================
.. _install.odes_dependencies:
diff --git a/docs/source/user_guide/installation/windows.rst b/docs/source/user_guide/installation/windows.rst
index 5b104e91bd..5ad77b6f7f 100644
--- a/docs/source/user_guide/installation/windows.rst
+++ b/docs/source/user_guide/installation/windows.rst
@@ -66,6 +66,21 @@ installed automatically when you install PyBaMM using ``pip``.
For an introduction to virtual environments, see
(https://realpython.com/python-virtual-environments-a-primer/).
+Optional - JaxSolver
+~~~~~~~~~~~~~~~~~~~~
+
+Users can install ``jax`` and ``jaxlib`` to use the Jax solver.
+
+.. note::
+
+ The Jax solver is not supported on Python 3.8. It is supported on Python 3.9, 3.10, and 3.11.
+
+.. code:: bash
+
+ pip install "pybamm[jax]"
+
+The ``pip install "pybamm[jax]"`` command automatically downloads and installs ``pybamm`` and the compatible versions of ``jax`` and ``jaxlib`` on your system. (``pybamm_install_jax`` is deprecated.)
+
Uninstall PyBaMM
----------------
diff --git a/examples/scripts/minimal_example_of_lookup_tables.py b/examples/scripts/minimal_example_of_lookup_tables.py
new file mode 100644
index 0000000000..1c93e311c0
--- /dev/null
+++ b/examples/scripts/minimal_example_of_lookup_tables.py
@@ -0,0 +1,51 @@
+import pybamm
+import pandas as pd
+import numpy as np
+
+
+def process_2D(name, data):
+ data = data.to_numpy()
+ x1 = np.unique(data[:, 0])
+ x2 = np.unique(data[:, 1])
+
+ value = data[:, 2]
+
+ x = (x1, x2)
+
+ value_data = value.reshape(len(x1), len(x2), order="C")
+
+ formatted_data = (name, (x, value_data))
+
+ return formatted_data
+
+
+parameter_values = pybamm.ParameterValues(pybamm.parameter_sets.Chen2020)
+
+# overwrite the diffusion coefficient with a 2D lookup table
+D_s_n = parameter_values["Negative electrode diffusivity [m2.s-1]"]
+df = pd.DataFrame(
+ {
+ "T": [0, 0, 25, 25, 45, 45],
+ "sto": [0, 1, 0, 1, 0, 1],
+ "D_s_n": [D_s_n, D_s_n, D_s_n, D_s_n, D_s_n, D_s_n],
+ }
+)
+df["T"] = df["T"] + 273.15
+D_s_n_data = process_2D("Negative electrode diffusivity [m2.s-1]", df)
+
+
+def D_s_n(sto, T):
+ name, (x, y) = D_s_n_data
+ return pybamm.Interpolant(x, y, [T, sto], name)
+
+
+parameter_values["Negative electrode diffusivity [m2.s-1]"] = D_s_n
+
+k_n = parameter_values["Negative electrode exchange-current density [A.m-2]"]
+
+model = pybamm.lithium_ion.DFN()
+sim = pybamm.Simulation(model, parameter_values=parameter_values)
+
+sim.solve([0, 30])
+
+sim.plot(["Negative particle surface concentration [mol.m-3]"])
diff --git a/noxfile.py b/noxfile.py
index fd033a4573..297fc5b3d7 100644
--- a/noxfile.py
+++ b/noxfile.py
@@ -61,9 +61,12 @@ def run_coverage(session):
set_environment_variables(PYBAMM_ENV, session=session)
session.install("coverage", silent=False)
if sys.platform != "win32":
- session.install("-e", ".[all,odes,jax]", silent=False)
+ session.install("-e", ".[all,jax,odes]", silent=False)
else:
- session.install("-e", ".[all]", silent=False)
+ if sys.version_info < (3, 9):
+ session.install("-e", ".[all]", silent=False)
+ else:
+ session.install("-e", ".[all,jax]", silent=False)
session.run("coverage", "run", "run-tests.py", "--nosub")
session.run("coverage", "combine")
session.run("coverage", "xml")
@@ -74,9 +77,12 @@ def run_integration(session):
"""Run the integration tests."""
set_environment_variables(PYBAMM_ENV, session=session)
if sys.platform != "win32":
- session.install("-e", ".[all,odes,jax]", silent=False)
+ session.install("-e", ".[all,jax,odes]", silent=False)
else:
- session.install("-e", ".[all]", silent=False)
+ if sys.version_info < (3, 9):
+ session.install("-e", ".[all]", silent=False)
+ else:
+ session.install("-e", ".[all,jax]", silent=False)
session.run("python", "run-tests.py", "--integration")
@@ -92,9 +98,12 @@ def run_unit(session):
"""Run the unit tests."""
set_environment_variables(PYBAMM_ENV, session=session)
if sys.platform != "win32":
- session.install("-e", ".[all,odes,jax]", silent=False)
+ session.install("-e", ".[all,jax,odes]", silent=False)
else:
- session.install("-e", ".[all]", silent=False)
+ if sys.version_info < (3, 9):
+ session.install("-e", ".[all]", silent=False)
+ else:
+ session.install("-e", ".[all,jax]", silent=False)
session.run("python", "run-tests.py", "--unit")
@@ -144,7 +153,24 @@ def set_dev(session):
external=True,
)
else:
- session.run(python, "-m", "pip", "install", "-e", ".[all,dev]", external=True)
+ if sys.version_info < (3, 9):
+ session.run(
+ python,
+ "-m",
+ "pip",
+ "install",
+ ".[all,dev]",
+ external=True,
+ )
+ else:
+ session.run(
+ python,
+ "-m",
+ "pip",
+ "install",
+ ".[all,dev,jax]",
+ external=True,
+ )
@nox.session(name="tests")
@@ -152,9 +178,12 @@ def run_tests(session):
"""Run the unit tests and integration tests sequentially."""
set_environment_variables(PYBAMM_ENV, session=session)
if sys.platform != "win32":
- session.install("-e", ".[all,odes,jax]", silent=False)
+ session.install("-e", ".[all,jax,odes]", silent=False)
else:
- session.install("-e", ".[all]", silent=False)
+ if sys.version_info < (3, 9):
+ session.install("-e", ".[all]", silent=False)
+ else:
+ session.install("-e", ".[all,jax]", silent=False)
session.run("python", "run-tests.py", "--all")
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index cab7914cd9..00c596314a 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -54,6 +54,7 @@
from .logger import logger, set_logging_level, get_new_logger
from .settings import settings
from .citations import Citations, citations, print_citations
+
#
# Classes for the Expression Tree
#
@@ -227,12 +228,13 @@
#
from .experiment.experiment import Experiment
from . import experiment
+from .experiment import step
#
# Plotting
#
-from .plotting.quick_plot import QuickPlot, close_plots
+from .plotting.quick_plot import QuickPlot, close_plots, QuickPlotAxes
from .plotting.plot import plot
from .plotting.plot2D import plot2D
from .plotting.plot_voltage_components import plot_voltage_components
diff --git a/pybamm/citations.py b/pybamm/citations.py
index e73351a4c6..70bf4ba9d3 100644
--- a/pybamm/citations.py
+++ b/pybamm/citations.py
@@ -21,7 +21,6 @@ class Citations:
Examples
--------
- >>> import pybamm
>>> pybamm.citations.register("Sulzer2021")
>>> pybamm.citations.register("@misc{Newton1687, title={Mathematical...}}")
>>> pybamm.print_citations("citations.txt")
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index bb6e678f4c..62110b1676 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -865,6 +865,10 @@ def _process_symbol(self, symbol):
return child_spatial_method.boundary_value_or_flux(
symbol, disc_child, self.bcs
)
+ elif isinstance(symbol, pybamm.EvaluateAt):
+ return child_spatial_method.evaluate_at(
+ symbol, disc_child, symbol.position
+ )
elif isinstance(symbol, pybamm.UpwindDownwind):
direction = symbol.name # upwind or downwind
return spatial_method.upwind_or_downwind(
diff --git a/pybamm/experiment/experiment.py b/pybamm/experiment/experiment.py
index 9b02e3a20f..898d9b0f79 100644
--- a/pybamm/experiment/experiment.py
+++ b/pybamm/experiment/experiment.py
@@ -3,7 +3,7 @@
#
import pybamm
-from pybamm.step._steps_util import (
+from .step._steps_util import (
_convert_time_to_seconds,
_convert_temperature_to_kelvin,
)
diff --git a/pybamm/step/__init__.py b/pybamm/experiment/step/__init__.py
similarity index 58%
rename from pybamm/step/__init__.py
rename to pybamm/experiment/step/__init__.py
index eea47a54a7..e3b9ff8bd0 100644
--- a/pybamm/step/__init__.py
+++ b/pybamm/experiment/step/__init__.py
@@ -1,2 +1,3 @@
from .steps import *
from .steps import _Step
+from .step_termination import *
diff --git a/pybamm/step/_steps_util.py b/pybamm/experiment/step/_steps_util.py
similarity index 96%
rename from pybamm/step/_steps_util.py
rename to pybamm/experiment/step/_steps_util.py
index e524bc6064..f44cf52113 100644
--- a/pybamm/step/_steps_util.py
+++ b/pybamm/experiment/step/_steps_util.py
@@ -4,6 +4,7 @@
import pybamm
import numpy as np
from datetime import datetime
+from .step_termination import _read_termination
_examples = """
@@ -136,8 +137,10 @@ def __init__(
termination = [termination]
self.termination = []
for term in termination:
- typ, value = _convert_electric(term)
- self.termination.append({"type": typ, "value": value})
+ if isinstance(term, str):
+ term = _convert_electric(term)
+ term = _read_termination(term)
+ self.termination.append(term)
self.temperature = _convert_temperature_to_kelvin(temperature)
@@ -193,10 +196,7 @@ def to_dict(self):
}
def __eq__(self, other):
- return (
- isinstance(other, _Step)
- and self.hash_args == other.hash_args
- )
+ return isinstance(other, _Step) and self.hash_args == other.hash_args
def __hash__(self):
return hash(self.basic_repr())
diff --git a/pybamm/experiment/step/step_termination.py b/pybamm/experiment/step/step_termination.py
new file mode 100644
index 0000000000..082711a305
--- /dev/null
+++ b/pybamm/experiment/step/step_termination.py
@@ -0,0 +1,164 @@
+import pybamm
+import numpy as np
+
+
+class BaseTermination:
+ """
+ Base class for a termination event for an experiment step. To create a custom
+ termination, a class must implement `get_event` to return a :class:`pybamm.Event`
+ corresponding to the desired termination. In most cases the class
+ :class:`pybamm.step.CustomTermination` can be used to assist with this.
+
+ Parameters
+ ----------
+ value : float
+ The value at which the event is triggered
+ """
+
+ def __init__(self, value):
+ self.value = value
+
+ def get_event(self, variables, step_value):
+ """
+ Return a :class:`pybamm.Event` object corresponding to the termination event
+
+ Parameters
+ ----------
+ variables : dict
+ Dictionary of model variables, to be used for selecting the variable(s) that
+ determine the event
+ step_value : float or :class:`pybamm.Symbol`
+ Value of the step for which this is a termination event, to be used in some
+ cases to determine the sign of the event.
+ """
+ raise NotImplementedError
+
+ def __eq__(self, other):
+ # objects are equal if they have the same type and value
+ if isinstance(other, self.__class__):
+ return self.value == other.value
+ else:
+ return False
+
+
+class CrateTermination(BaseTermination):
+ """
+ Termination based on C-rate, created when a string termination of the C-rate type
+ (e.g. "C/10") is provided
+ """
+
+ def get_event(self, variables, step_value):
+ """
+ See :meth:`BaseTermination.get_event`
+ """
+ event = pybamm.Event(
+ "C-rate cut-off [experiment]",
+ abs(variables["C-rate"]) - self.value,
+ )
+ return event
+
+
+class CurrentTermination(BaseTermination):
+ """
+ Termination based on current, created when a string termination of the current type
+ (e.g. "1A") is provided
+ """
+
+ def get_event(self, variables, step_value):
+ """
+ See :meth:`BaseTermination.get_event`
+ """
+ event = pybamm.Event(
+ "Current cut-off [A] [experiment]",
+ abs(variables["Current [A]"]) - self.value,
+ )
+ return event
+
+
+class VoltageTermination(BaseTermination):
+ """
+ Termination based on voltage, created when a string termination of the voltage type
+ (e.g. "4.2V") is provided
+ """
+
+ def get_event(self, variables, step_value):
+ """
+ See :meth:`BaseTermination.get_event`
+ """
+ # The voltage event should be positive at the start of charge/
+ # discharge. We use the sign of the current or power input to
+ # figure out whether the voltage event is greater than the starting
+ # voltage (charge) or less (discharge) and set the sign of the
+ # event accordingly
+ if isinstance(step_value, pybamm.Symbol):
+ inpt = {"start time": 0}
+ init_curr = step_value.evaluate(t=0, inputs=inpt).flatten()[0]
+ else:
+ init_curr = step_value
+ sign = np.sign(init_curr)
+ if sign > 0:
+ name = "Discharge"
+ else:
+ name = "Charge"
+ if sign != 0:
+ # Event should be positive at initial conditions for both
+ # charge and discharge
+ event = pybamm.Event(
+ f"{name} voltage cut-off [V] [experiment]",
+ sign * (variables["Battery voltage [V]"] - self.value),
+ )
+ return event
+
+
+class CustomTermination(BaseTermination):
+ """
+ Define a custom termination event using a function. This can be used to create an
+ event based on any variable in the model.
+
+ Parameters
+ ----------
+ name : str
+ Name of the event
+ event_function : callable
+ A function that takes in a dictionary of variables and evaluates the event
+ value. Must be positive before the event is triggered and zero when the
+ event is triggered.
+
+ Example
+ -------
+ Add a cut-off based on negative electrode stoichiometry. The event will trigger
+ when the negative electrode stoichiometry reaches 10%.
+
+ >>> def neg_stoich_cutoff(variables):
+ ... return variables["Negative electrode stoichiometry"] - 0.1
+
+ >>> neg_stoich_termination = pybamm.step.CustomTermination(
+ ... name="Negative stoichiometry cut-off", event_function=neg_stoich_cutoff
+ ... )
+ """
+
+ def __init__(self, name, event_function):
+ if not name.endswith(" [experiment]"):
+ name += " [experiment]"
+ self.name = name
+ self.event_function = event_function
+
+ def get_event(self, variables, step_value):
+ """
+ See :meth:`BaseTermination.get_event`
+ """
+ return pybamm.Event(self.name, self.event_function(variables))
+
+
+def _read_termination(termination):
+ if isinstance(termination, tuple):
+ typ, value = termination
+ else:
+ return termination
+
+ termination_class = {
+ "current": CurrentTermination,
+ "voltage": VoltageTermination,
+ "C-rate": CrateTermination,
+ }[typ]
+ return termination_class(value)
diff --git a/pybamm/step/steps.py b/pybamm/experiment/step/steps.py
similarity index 93%
rename from pybamm/step/steps.py
rename to pybamm/experiment/step/steps.py
index 0b8123ddb0..2f2e9e31a4 100644
--- a/pybamm/step/steps.py
+++ b/pybamm/experiment/step/steps.py
@@ -120,7 +120,7 @@ def current(value, **kwargs):
Parameters
----------
value : float
- The current value in A.
+ The current value in A. It can be a number or a 2-column array (for drive cycles).
**kwargs
Any other keyword arguments are passed to the :class:`pybamm.step._Step`
class.
@@ -141,7 +141,7 @@ def c_rate(value, **kwargs):
Parameters
----------
value : float
- The C-rate value.
+ The C-rate value. It can be a number or a 2-column array (for drive cycles).
**kwargs
Any other keyword arguments are passed to the :class:`pybamm.step._Step`
class.
@@ -162,7 +162,7 @@ def voltage(value, **kwargs):
Parameters
----------
value : float
- The voltage value in V.
+ The voltage value in V. It can be a number or a 2-column array (for drive cycles).
**kwargs
Any other keyword arguments are passed to the :class:`pybamm.step._Step`
class.
@@ -183,7 +183,7 @@ def power(value, **kwargs):
Parameters
----------
value : float
- The power value in W.
+ The power value in W. It can be a number or a 2-column array (for drive cycles).
**kwargs
Any other keyword arguments are passed to the :class:`pybamm.step._Step`
class.
@@ -204,7 +204,7 @@ def resistance(value, **kwargs):
Parameters
----------
value : float
- The resistance value in Ohm.
+ The resistance value in Ohm. It can be a number or a 2-column array (for drive cycles).
**kwargs
Any other keyword arguments are passed to the :class:`pybamm.step._Step`
class.
diff --git a/pybamm/expression_tree/interpolant.py b/pybamm/expression_tree/interpolant.py
index 8296625da0..1cb5e70d05 100644
--- a/pybamm/expression_tree/interpolant.py
+++ b/pybamm/expression_tree/interpolant.py
@@ -302,10 +302,11 @@ def _function_evaluate(self, evaluated_children):
new_evaluated_children, self.function.grid
):
nan_children.append(np.ones_like(child) * interp_range.mean())
- return self.function(np.transpose(nan_children)) * np.nan
+ nan_eval = self.function(np.transpose(nan_children))
+ return np.reshape(nan_eval, shape)
else:
res = self.function(np.transpose(new_evaluated_children))
- return res[:, np.newaxis]
+ return np.reshape(res, shape)
else: # pragma: no cover
raise ValueError("Invalid dimension: {0}".format(self.dimension))
diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py
index e67ea6a89f..2c3166582e 100644
--- a/pybamm/expression_tree/symbol.py
+++ b/pybamm/expression_tree/symbol.py
@@ -423,7 +423,12 @@ def set_id(self):
need to hash once.
"""
self._id = hash(
- (self.__class__, self.name, *tuple([child.id for child in self.children]), *tuple([(k, tuple(v)) for k, v in self.domains.items() if v != []]))
+ (
+ self.__class__,
+ self.name,
+ *tuple([child.id for child in self.children]),
+ *tuple([(k, tuple(v)) for k, v in self.domains.items() if v != []]),
+ )
)
@property
@@ -532,7 +537,6 @@ def pre_order(self):
Examples
--------
- >>> import pybamm
>>> a = pybamm.Symbol('a')
>>> b = pybamm.Symbol('b')
>>> for node in (a*b).pre_order():
diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index 2705685814..319429183c 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -325,7 +325,14 @@ def _unary_jac(self, child_jac):
def set_id(self):
"""See :meth:`pybamm.Symbol.set_id()`"""
self._id = hash(
- (self.__class__, self.name, self.slice.start, self.slice.stop, self.children[0].id, *tuple(self.domain))
+ (
+ self.__class__,
+ self.name,
+ self.slice.start,
+ self.slice.stop,
+ self.children[0].id,
+ *tuple(self.domain),
+ )
)
def _unary_evaluate(self, child):
@@ -465,7 +472,9 @@ def _unary_new_copy(self, child):
def _sympy_operator(self, child):
"""Override :meth:`pybamm.UnaryOperator._sympy_operator`"""
- sympy_Divergence = have_optional_dependency("sympy.vector.operators", "Divergence")
+ sympy_Divergence = have_optional_dependency(
+ "sympy.vector.operators", "Divergence"
+ )
return sympy_Divergence(child)
@@ -616,7 +625,18 @@ def integration_variable(self):
def set_id(self):
"""See :meth:`pybamm.Symbol.set_id()`"""
self._id = hash(
- (self.__class__, self.name, *tuple([integration_variable.id for integration_variable in self.integration_variable]), self.children[0].id, *tuple(self.domain))
+ (
+ self.__class__,
+ self.name,
+ *tuple(
+ [
+ integration_variable.id
+ for integration_variable in self.integration_variable
+ ]
+ ),
+ self.children[0].id,
+ *tuple(self.domain),
+ )
)
def _unary_new_copy(self, child):
@@ -757,7 +777,13 @@ def __init__(self, child, vector_type="row"):
def set_id(self):
"""See :meth:`pybamm.Symbol.set_id()`"""
self._id = hash(
- (self.__class__, self.name, self.vector_type, self.children[0].id, *tuple(self.domain))
+ (
+ self.__class__,
+ self.name,
+ self.vector_type,
+ self.children[0].id,
+ *tuple(self.domain),
+ )
)
def _unary_new_copy(self, child):
@@ -851,7 +877,13 @@ def __init__(self, child, side, domain):
def set_id(self):
"""See :meth:`pybamm.Symbol.set_id()`"""
self._id = hash(
- (self.__class__, self.name, self.side, self.children[0].id, *tuple([(k, tuple(v)) for k, v in self.domains.items()]))
+ (
+ self.__class__,
+ self.name,
+ self.side,
+ self.children[0].id,
+ *tuple([(k, tuple(v)) for k, v in self.domains.items()]),
+ )
)
def _evaluates_on_edges(self, dimension):
@@ -872,7 +904,8 @@ def evaluate_for_shape(self):
class BoundaryOperator(SpatialOperator):
"""
- A node in the expression tree which gets the boundary value of a variable.
+ A node in the expression tree which gets the boundary value of a variable on its
+ primary domain.
Parameters
----------
@@ -909,7 +942,13 @@ def __init__(self, name, child, side):
def set_id(self):
"""See :meth:`pybamm.Symbol.set_id()`"""
self._id = hash(
- (self.__class__, self.name, self.side, self.children[0].id, *tuple([(k, tuple(v)) for k, v in self.domains.items()]))
+ (
+ self.__class__,
+ self.name,
+ self.side,
+ self.children[0].id,
+ *tuple([(k, tuple(v)) for k, v in self.domains.items()]),
+ )
)
def _unary_new_copy(self, child):
@@ -923,7 +962,8 @@ def _evaluate_for_shape(self):
class BoundaryValue(BoundaryOperator):
"""
- A node in the expression tree which gets the boundary value of a variable.
+ A node in the expression tree which gets the boundary value of a variable on its
+ primary domain.
Parameters
----------
@@ -998,7 +1038,8 @@ def to_json(self):
class BoundaryGradient(BoundaryOperator):
"""
- A node in the expression tree which gets the boundary flux of a variable.
+ A node in the expression tree which gets the boundary flux of a variable on its
+ primary domain.
Parameters
----------
@@ -1012,6 +1053,59 @@ def __init__(self, child, side):
super().__init__("boundary flux", child, side)
+class EvaluateAt(SpatialOperator):
+ """
+ A node in the expression tree which evaluates a symbol at a given position in space
+ in its primary domain. Currently this is only implemented for 1D primary domains.
+
+ Parameters
+ ----------
+ child : :class:`pybamm.Symbol`
+ The variable to evaluate
+ position : :class:`pybamm.Symbol`
+ The position in space on the symbol's primary domain at which to evaluate
+ the symbol.
+ """
+
+ def __init__(self, child, position):
+ self.position = position
+
+ # "evaluate at" of a child takes the primary domain from secondary domain
+ # of the child
+ # tertiary auxiliary domain shift down to secondary, quarternary to tertiary
+ domains = {
+ "primary": child.domains["secondary"],
+ "secondary": child.domains["tertiary"],
+ "tertiary": child.domains["quaternary"],
+ }
+
+ super().__init__("evaluate", child, domains)
+
+ def set_id(self):
+ """See :meth:`pybamm.Symbol.set_id()`"""
+ self._id = hash(
+ (
+ self.__class__,
+ self.name,
+ self.position,
+ self.children[0].id,
+ *tuple([(k, tuple(v)) for k, v in self.domains.items()]),
+ )
+ )
+
+ def _unary_jac(self, child_jac):
+ """See :meth:`pybamm.UnaryOperator._unary_jac()`."""
+ return pybamm.Scalar(0)
+
+ def _unary_new_copy(self, child):
+ """See :meth:`UnaryOperator._unary_new_copy()`."""
+ return self.__class__(child, self.position)
+
+ def _evaluate_for_shape(self):
+ """See :meth:`pybamm.Symbol.evaluate_for_shape_using_domain()`"""
+ return pybamm.evaluate_for_shape_using_domain(self.domains)
+
+
class UpwindDownwind(SpatialOperator):
"""
A node in the expression tree representing an upwinding or downwinding operator.
diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index c4ae414e7c..257bc30ef8 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -422,6 +422,7 @@ def input_parameters(self):
return self._input_parameters
def print_parameter_info(self):
+ """Returns parameters used in the model"""
self._parameter_info = ""
parameters = self._find_symbols(pybamm.Parameter)
for param in parameters:
diff --git a/pybamm/models/full_battery_models/equivalent_circuit/thevenin.py b/pybamm/models/full_battery_models/equivalent_circuit/thevenin.py
index cc456c8765..407039b6f6 100644
--- a/pybamm/models/full_battery_models/equivalent_circuit/thevenin.py
+++ b/pybamm/models/full_battery_models/equivalent_circuit/thevenin.py
@@ -52,7 +52,6 @@ class Thevenin(pybamm.BaseModel):
Examples
--------
- >>> import pybamm
>>> model = pybamm.equivalent_circuit.Thevenin()
>>> model.name
'Thevenin Equivalent Circuit Model'
diff --git a/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py b/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py
index cc736d6d04..fbe19b0d42 100644
--- a/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py
+++ b/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py
@@ -462,3 +462,51 @@ def set_convection_submodel(self):
self.submodels[
"through-cell convection"
] = pybamm.convection.through_cell.NoConvection(self.param, self.options)
+
+ def insert_reference_electrode(self, position=None):
+ """
+ Insert a reference electrode to measure the electrolyte potential at a given
+ position in space. Adds model variables for the electrolyte potential at the
+ reference electrode and for the potential difference between the electrode
+ potentials measured at the electrode/current collector interface and the
+ reference electrode. Only implemented for 1D models (i.e. where the
+ 'dimensionality' option is 0).
+
+ Parameters
+ ----------
+ position : :class:`pybamm.Symbol`, optional
+ The position in space at which to measure the electrolyte potential. If
+ None, defaults to the mid-point of the separator.
+ """
+ if self.options["dimensionality"] != 0:
+ raise NotImplementedError(
+ "Reference electrode can only be inserted for models where "
+ "'dimensionality' is 0. For other models, please add a reference "
+ "electrode manually."
+ )
+
+ param = self.param
+ if position is None:
+ position = param.n.L + param.s.L / 2
+
+ phi_e_ref = pybamm.EvaluateAt(
+ self.variables["Electrolyte potential [V]"], position
+ )
+ phi_p = pybamm.boundary_value(
+ self.variables["Positive electrode potential [V]"], "right"
+ )
+ variables = {
+ "Positive electrode 3E potential [V]": phi_p - phi_e_ref,
+ "Reference electrode potential [V]": phi_e_ref,
+ }
+
+ if self.options["working electrode"] == "both":
+ phi_n = pybamm.boundary_value(
+ self.variables["Negative electrode potential [V]"], "left"
+ )
+ variables.update(
+ {
+ "Negative electrode 3E potential [V]": phi_n - phi_e_ref,
+ }
+ )
+ self.variables.update(variables)
diff --git a/pybamm/models/full_battery_models/lithium_ion/dfn.py b/pybamm/models/full_battery_models/lithium_ion/dfn.py
index 5f0a2cfb3e..db4e0282d8 100644
--- a/pybamm/models/full_battery_models/lithium_ion/dfn.py
+++ b/pybamm/models/full_battery_models/lithium_ion/dfn.py
@@ -13,7 +13,6 @@ class DFN(BaseModel):
Examples
--------
- >>> import pybamm
>>> model = pybamm.lithium_ion.DFN()
>>> model.name
'Doyle-Fuller-Newman model'
diff --git a/pybamm/models/full_battery_models/lithium_ion/mpm.py b/pybamm/models/full_battery_models/lithium_ion/mpm.py
index fa2062a95c..40a533536f 100644
--- a/pybamm/models/full_battery_models/lithium_ion/mpm.py
+++ b/pybamm/models/full_battery_models/lithium_ion/mpm.py
@@ -13,7 +13,6 @@ class MPM(SPM):
Examples
--------
- >>> import pybamm
>>> model = pybamm.lithium_ion.MPM()
>>> model.name
'Many-Particle Model'
diff --git a/pybamm/models/full_battery_models/lithium_ion/spm.py b/pybamm/models/full_battery_models/lithium_ion/spm.py
index bdebf12aef..386c55ded9 100644
--- a/pybamm/models/full_battery_models/lithium_ion/spm.py
+++ b/pybamm/models/full_battery_models/lithium_ion/spm.py
@@ -13,7 +13,6 @@ class SPM(BaseModel):
Examples
--------
- >>> import pybamm
>>> model = pybamm.lithium_ion.SPM()
>>> model.name
'Single Particle Model'
diff --git a/pybamm/models/full_battery_models/lithium_ion/spme.py b/pybamm/models/full_battery_models/lithium_ion/spme.py
index e293a5fabd..103f13415a 100644
--- a/pybamm/models/full_battery_models/lithium_ion/spme.py
+++ b/pybamm/models/full_battery_models/lithium_ion/spme.py
@@ -14,7 +14,6 @@ class SPMe(SPM):
Examples
--------
- >>> import pybamm
>>> model = pybamm.lithium_ion.SPMe()
>>> model.name
'Single Particle Model with electrolyte'
diff --git a/pybamm/models/submodels/interface/open_circuit_potential/base_ocp.py b/pybamm/models/submodels/interface/open_circuit_potential/base_ocp.py
index e93dcefa3e..35f3894dfe 100644
--- a/pybamm/models/submodels/interface/open_circuit_potential/base_ocp.py
+++ b/pybamm/models/submodels/interface/open_circuit_potential/base_ocp.py
@@ -90,8 +90,8 @@ def _get_standard_ocp_variables(self, ocp_surf, ocp_bulk, dUdT):
if self.reaction in ["lithium-ion main", "lead-acid main"]:
variables.update(
{
- f"{Domain} electrode entropic change [V.K-1]": dUdT,
- f"X-averaged {domain} electrode entropic change [V.K-1]": dUdT_av,
+ f"{Domain} electrode {reaction_name}entropic change [V.K-1]": dUdT,
+ f"X-averaged {domain} electrode {reaction_name}entropic change [V.K-1]": dUdT_av,
}
)
diff --git a/pybamm/models/submodels/particle/base_particle.py b/pybamm/models/submodels/particle/base_particle.py
index dd5a94afc6..0f46615724 100644
--- a/pybamm/models/submodels/particle/base_particle.py
+++ b/pybamm/models/submodels/particle/base_particle.py
@@ -56,8 +56,9 @@ def _get_effective_diffusivity(self, c, T, current):
if stress_option == "true":
# Ai2019 eq [12]
- Omega = domain_param.Omega
- E = domain_param.E
+ sto = c / phase_param.c_max
+ Omega = pybamm.r_average(domain_param.Omega(sto, T))
+ E = pybamm.r_average(domain_param.E(sto, T))
nu = domain_param.nu
theta_M = Omega / (param.R * T) * (2 * Omega * E) / (9 * (1 - nu))
stress_factor = 1 + theta_M * (c - domain_param.c_0)
diff --git a/pybamm/models/submodels/particle_mechanics/base_mechanics.py b/pybamm/models/submodels/particle_mechanics/base_mechanics.py
index feffbdd380..35adadf47d 100644
--- a/pybamm/models/submodels/particle_mechanics/base_mechanics.py
+++ b/pybamm/models/submodels/particle_mechanics/base_mechanics.py
@@ -43,12 +43,19 @@ def _get_mechanical_results(self, variables):
sto_rav = variables[f"R-averaged {domain} particle concentration"]
c_s_surf = variables[f"{Domain} particle surface concentration [mol.m-3]"]
T_xav = variables["X-averaged cell temperature [K]"]
+ phase_name = self.phase_name
+ T = pybamm.PrimaryBroadcast(
+ variables[f"{Domain} electrode temperature [K]"],
+ [f"{domain} {phase_name}particle"],
+ )
eps_s = variables[f"{Domain} electrode active material volume fraction"]
- Omega = domain_param.Omega
+ #use a tangential approximation for omega
+ sto = variables[f"{Domain} particle concentration"]
+ Omega = pybamm.r_average(domain_param.Omega(sto, T))
R0 = domain_param.prim.R
c_0 = domain_param.c_0
- E0 = domain_param.E
+ E0 = pybamm.r_average(domain_param.E(sto, T))
nu = domain_param.nu
L0 = domain_param.L
sto_init = pybamm.r_average(domain_param.prim.c_init / domain_param.prim.c_max)
diff --git a/pybamm/models/submodels/thermal/base_thermal.py b/pybamm/models/submodels/thermal/base_thermal.py
index 4c476ee897..27e52d6638 100644
--- a/pybamm/models/submodels/thermal/base_thermal.py
+++ b/pybamm/models/submodels/thermal/base_thermal.py
@@ -117,38 +117,59 @@ def _get_standard_coupled_variables(self, variables):
# Total Ohmic heating
Q_ohm = Q_ohm_s + Q_ohm_e
- # Irreversible electrochemical heating
- a_j_p = variables[
- "Positive electrode volumetric interfacial current density [A.m-3]"
- ]
- eta_r_p = variables["Positive electrode reaction overpotential [V]"]
+ num_phases = int(getattr(self.options, 'positive')["particle phases"])
+ phase_names = [""]
+ if num_phases > 1:
+ phase_names = ["primary ", "secondary "]
+
+ Q_rxn_p, Q_rev_p = 0, 0
+ T_p = variables["Positive electrode temperature [K]"]
+ for phase in phase_names:
+ a_j_p = variables[
+ f"Positive electrode {phase}volumetric interfacial current density [A.m-3]"
+ ]
+ eta_r_p = variables[f"Positive electrode {phase}reaction overpotential [V]"]
+ # Irreversible electrochemical heating
+ Q_rxn_p += a_j_p * eta_r_p
+ # Reversible electrochemical heating
+ dUdT_p = variables[f"Positive electrode {phase}entropic change [V.K-1]"]
+ Q_rev_p += a_j_p * T_p * dUdT_p
+
+
+ num_phases = int(getattr(self.options, 'negative')["particle phases"])
+ phase_names = [""]
+ if num_phases > 1:
+ phase_names = ["primary", "secondary"]
+
if self.options.electrode_types["negative"] == "planar":
Q_rxn_n = pybamm.FullBroadcast(
0, ["negative electrode"], "current collector"
)
+ Q_rev_n = pybamm.FullBroadcast(
+ 0, ["negative electrode"], "current collector"
+ )
else:
- a_j_n = variables[
- "Negative electrode volumetric interfacial current density [A.m-3]"
- ]
- eta_r_n = variables["Negative electrode reaction overpotential [V]"]
- Q_rxn_n = a_j_n * eta_r_n
- Q_rxn_p = a_j_p * eta_r_p
+ T_n = variables["Negative electrode temperature [K]"]
+ Q_rxn_n = 0
+ Q_rev_n = 0
+ for phase in phase_names:
+ a_j_n = variables[
+ f"Negative electrode {phase}volumetric interfacial current density [A.m-3]"
+ ]
+ eta_r_n = variables[f"Negative electrode {phase}reaction overpotential [V]"]
+ # Irreversible electrochemical heating
+ Q_rxn_n += a_j_n * eta_r_n
+
+ # Reversible electrochemical heating
+ dUdT_n = variables[f"Negative electrode {phase}entropic change [V.K-1]"]
+ Q_rev_n += a_j_n * T_n * dUdT_n
+
+ # Irreversible electrochemical heating
Q_rxn = pybamm.concatenation(
Q_rxn_n, pybamm.FullBroadcast(0, "separator", "current collector"), Q_rxn_p
)
# Reversible electrochemical heating
- T_p = variables["Positive electrode temperature [K]"]
- dUdT_p = variables["Positive electrode entropic change [V.K-1]"]
- if self.options.electrode_types["negative"] == "planar":
- Q_rev_n = pybamm.FullBroadcast(
- 0, ["negative electrode"], "current collector"
- )
- else:
- T_n = variables["Negative electrode temperature [K]"]
- dUdT_n = variables["Negative electrode entropic change [V.K-1]"]
- Q_rev_n = a_j_n * T_n * dUdT_n
- Q_rev_p = a_j_p * T_p * dUdT_p
Q_rev = pybamm.concatenation(
Q_rev_n, pybamm.FullBroadcast(0, "separator", "current collector"), Q_rev_p
)
diff --git a/pybamm/parameters/lithium_ion_parameters.py b/pybamm/parameters/lithium_ion_parameters.py
index c459a4ef1e..12196c4044 100644
--- a/pybamm/parameters/lithium_ion_parameters.py
+++ b/pybamm/parameters/lithium_ion_parameters.py
@@ -303,14 +303,11 @@ def _set_parameters(self):
# Mechanical parameters
self.nu = pybamm.Parameter(f"{Domain} electrode Poisson's ratio")
- self.E = pybamm.Parameter(f"{Domain} electrode Young's modulus [Pa]")
self.c_0 = pybamm.Parameter(
f"{Domain} electrode reference concentration for free of deformation "
"[mol.m-3]"
)
- self.Omega = pybamm.Parameter(
- f"{Domain} electrode partial molar volume [m3.mol-1]"
- )
+
self.l_cr_0 = pybamm.Parameter(f"{Domain} electrode initial crack length [m]")
self.w_cr = pybamm.Parameter(f"{Domain} electrode initial crack width [m]")
self.rho_cr = pybamm.Parameter(
@@ -349,6 +346,22 @@ def C_dl(self, T):
f"{Domain} electrode double-layer capacity [F.m-2]", inputs
)
+ def Omega(self, sto, T):
+ """Dimensional partial molar volume of Li in solid solution [m3.mol-1]"""
+ Domain = self.domain.capitalize()
+ inputs = {f"{Domain} particle stoichiometry": sto, "Temperature [K]": T}
+ return pybamm.FunctionParameter(
+ f"{Domain} electrode partial molar volume [m3.mol-1]", inputs
+ )
+
+ def E(self, sto, T):
+ """Dimensional Young's modulus"""
+ Domain = self.domain.capitalize()
+ inputs = {f"{Domain} particle stoichiometry": sto, "Temperature [K]": T}
+ return pybamm.FunctionParameter(
+ f"{Domain} electrode Young's modulus [Pa]", inputs
+ )
+
def sigma(self, T):
"""Dimensional electrical conductivity in electrode"""
inputs = {"Temperature [K]": T}
diff --git a/pybamm/parameters/parameter_sets.py b/pybamm/parameters/parameter_sets.py
index 6c6201d9af..ea45f2df5c 100644
--- a/pybamm/parameters/parameter_sets.py
+++ b/pybamm/parameters/parameter_sets.py
@@ -16,7 +16,6 @@ class ParameterSets(Mapping):
.. doctest::
- >>> import pybamm
>>> list(pybamm.parameter_sets)
['Ai2020', 'Chen2020', ...]
@@ -24,7 +23,6 @@ class ParameterSets(Mapping):
.. doctest::
- >>> import pybamm
>>> print(pybamm.parameter_sets.get_docstring("Ai2020"))
Parameters for the Enertech cell (Ai2020), from the papers :footcite:t:`Ai2019`,
@@ -44,7 +42,7 @@ def __init__(self):
@staticmethod
def get_entries(group_name):
# Wrapper for the importlib version logic
- if sys.version_info < (3, 10): # pragma: no cover
+ if sys.version_info < (3, 10): # pragma: no cover
return importlib.metadata.entry_points()[group_name]
else:
return importlib.metadata.entry_points(group=group_name)
diff --git a/pybamm/parameters/parameter_values.py b/pybamm/parameters/parameter_values.py
index d5f12f362f..a5ed9b66fb 100644
--- a/pybamm/parameters/parameter_values.py
+++ b/pybamm/parameters/parameter_values.py
@@ -24,7 +24,6 @@ class ParameterValues:
Examples
--------
- >>> import pybamm
>>> values = {"some parameter": 1, "another parameter": 2}
>>> param = pybamm.ParameterValues(values)
>>> param["some parameter"]
@@ -721,6 +720,15 @@ def _process_symbol(self, symbol):
# f_a_dist in the size average needs to be processed
if isinstance(new_symbol, pybamm.SizeAverage):
new_symbol.f_a_dist = self.process_symbol(new_symbol.f_a_dist)
+ # position in evaluate at needs to be processed, and should be a Scalar
+ if isinstance(new_symbol, pybamm.EvaluateAt):
+ new_symbol_position = self.process_symbol(new_symbol.position)
+ if not isinstance(new_symbol_position, pybamm.Scalar):
+ raise ValueError(
+ "'position' in 'EvaluateAt' must evaluate to a scalar"
+ )
+ else:
+ new_symbol.position = new_symbol_position
return new_symbol
# Functions
diff --git a/pybamm/plotting/quick_plot.py b/pybamm/plotting/quick_plot.py
index f731a58e0e..686c58f3c5 100644
--- a/pybamm/plotting/quick_plot.py
+++ b/pybamm/plotting/quick_plot.py
@@ -486,14 +486,14 @@ def plot(self, t, dynamic=False):
self.plots = {}
self.time_lines = {}
self.colorbars = {}
- self.axes = []
+ self.axes = QuickPlotAxes()
# initialize empty handles, to be created only if the appropriate plots are made
solution_handles = []
for k, (key, variable_lists) in enumerate(self.variables.items()):
ax = self.fig.add_subplot(self.gridspec[k])
- self.axes.append(ax)
+ self.axes.add(key, ax)
x_min, x_max, y_min, y_max = self.axis_limits[key]
ax.set_xlim(x_min, x_max)
if y_min is not None and y_max is not None:
@@ -803,3 +803,40 @@ def create_gif(self, number_of_images=80, duration=0.1, output_filename="plot.gi
# remove the generated images
for image in images:
os.remove(image)
+
+
+class QuickPlotAxes:
+ """
+ Class to store axes for the QuickPlot
+ """
+
+ def __init__(self):
+ self._by_variable = {}
+ self._axes = []
+
+ def add(self, keys, axis):
+ """
+ Add axis
+
+ Parameters
+ ----------
+ keys : iter
+ Iterable of keys of variables being plotted on the axis
+ axis : matplotlib Axis object
+ The axis object
+ """
+ self._axes.append(axis)
+ for k in keys:
+ self._by_variable[k] = axis
+
+ def __getitem__(self, index):
+ """
+ Get axis by index
+ """
+ return self._axes[index]
+
+ def by_variable(self, key):
+ """
+ Get axis by variable name
+ """
+ return self._by_variable[key]
diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index bf418f068d..5c2cf0bff1 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -192,16 +192,6 @@ def set_up_and_parameterise_experiment(self):
op_conds.type = "current"
op_conds.value = op_conds.value * capacity
- # Update terminations
- termination = op_conds.termination
- for term in termination:
- term_type = term["type"]
- if term_type == "C-rate":
- # Change type to current
- term["type"] = "current"
- # Scale C-rate with capacity to obtain current
- term["value"] = term["value"] * capacity
-
# Add time to the experiment times
dt = op_conds.duration
if dt is None:
@@ -294,46 +284,9 @@ def set_up_and_parameterise_model_for_experiment(self):
def update_new_model_events(self, new_model, op):
for term in op.termination:
- if term["type"] == "current":
- new_model.events.append(
- pybamm.Event(
- "Current cut-off [A] [experiment]",
- abs(new_model.variables["Current [A]"]) - term["value"],
- )
- )
-
- # add voltage events to the model
- if term["type"] == "voltage":
- # The voltage event should be positive at the start of charge/
- # discharge. We use the sign of the current or power input to
- # figure out whether the voltage event is greater than the starting
- # voltage (charge) or less (discharge) and set the sign of the
- # event accordingly
- if isinstance(op.value, pybamm.Interpolant) or isinstance(
- op.value, pybamm.Multiplication
- ):
- inpt = {"start time": 0}
- init_curr = op.value.evaluate(t=0, inputs=inpt).flatten()[0]
- sign = np.sign(init_curr)
- else:
- sign = np.sign(op.value)
- if sign > 0:
- name = "Discharge"
- else:
- name = "Charge"
- if sign != 0:
- # Event should be positive at initial conditions for both
- # charge and discharge
- new_model.events.append(
- pybamm.Event(
- f"{name} voltage cut-off [V] [experiment]",
- sign
- * (
- new_model.variables["Battery voltage [V]"]
- - term["value"]
- ),
- )
- )
+ event = term.get_event(new_model.variables, op.value)
+ if event is not None:
+ new_model.events.append(event)
# Keep the min and max voltages as safeguards but add some tolerances
# so that they are not triggered before the voltage limits in the
@@ -1166,7 +1119,13 @@ def solution(self):
return self._solution
def save(self, filename):
- """Save simulation using pickle"""
+ """Save simulation using pickle module.
+
+ Parameters
+ ----------
+ filename : str
+ The file extension can be arbitrary, but it is common to use ".pkl" or ".pickle"
+ """
if self._model.convert_to_format == "python":
# We currently cannot save models in the 'python' format
raise NotImplementedError(
diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py
index 5c32e5a2c0..636243f829 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -1023,6 +1023,46 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
return boundary_value
+ def evaluate_at(self, symbol, discretised_child, position):
+ """
+ Returns the symbol evaluated at a given position in space.
+
+ Parameters
+ ----------
+ symbol: :class:`pybamm.Symbol`
+ The boundary value or flux symbol
+ discretised_child : :class:`pybamm.StateVector`
+ The discretised variable from which to calculate the boundary value
+ position : :class:`pybamm.Scalar`
+ The point in one-dimensional space at which to evaluate the symbol.
+
+ Returns
+ -------
+ :class:`pybamm.MatrixMultiplication`
+ The variable representing the value at the given point.
+ """
+ # Get mesh nodes
+ domain = discretised_child.domain
+ mesh = self.mesh[domain]
+ nodes = mesh.nodes
+ repeats = self._get_auxiliary_domain_repeats(discretised_child.domains)
+
+ # Find the index of the node closest to the value
+ index = np.argmin(np.abs(nodes - position.value))
+
+ # Create a sparse matrix with a 1 at the index
+ sub_matrix = csr_matrix(([1], ([0], [index])), shape=(1, mesh.npts))
+ # repeat across auxiliary domains
+ matrix = csr_matrix(kron(eye(repeats), sub_matrix))
+
+ # Index into the discretised child
+ out = pybamm.Matrix(matrix) @ discretised_child
+
+ # `EvaluateAt` removes domain
+ out.clear_domains()
+
+ return out
+
def process_binary_operators(self, bin_op, left, right, disc_left, disc_right):
"""Discretise binary operators in model equations. Performs appropriate
averaging of diffusivities if one of the children is a gradient operator, so
diff --git a/pybamm/spatial_methods/spatial_method.py b/pybamm/spatial_methods/spatial_method.py
index acb0227bc2..a461d6c150 100644
--- a/pybamm/spatial_methods/spatial_method.py
+++ b/pybamm/spatial_methods/spatial_method.py
@@ -331,7 +331,7 @@ def internal_neumann_condition(
def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
"""
- Returns the boundary value or flux using the approriate expression for the
+ Returns the boundary value or flux using the appropriate expression for the
spatial method. To do this, we create a sparse vector 'bv_vector' that extracts
either the first (for side="left") or last (for side="right") point from
'discretised_child'.
@@ -377,6 +377,26 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
out.clear_domains()
return out
+ def evaluate_at(self, symbol, discretised_child, position):
+ """
+ Returns the symbol evaluated at a given position in space.
+
+ Parameters
+ ----------
+ symbol: :class:`pybamm.Symbol`
+ The boundary value or flux symbol
+ discretised_child : :class:`pybamm.StateVector`
+ The discretised variable from which to calculate the boundary value
+ position : :class:`pybamm.Scalar`
+ The point in one-dimensional space at which to evaluate the symbol.
+
+ Returns
+ -------
+ :class:`pybamm.MatrixMultiplication`
+ The variable representing the value at the given point.
+ """
+ raise NotImplementedError
+
def mass_matrix(self, symbol, boundary_conditions):
"""
Calculates the mass matrix for a spatial method.
diff --git a/pyproject.toml b/pyproject.toml
index 19c8800a63..7d25c8e140 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -110,13 +110,13 @@ dev = [
"nbmake",
]
# Reading CSV files
-pandas = [
+pandas = [
"pandas>=1.5.0",
]
# For the Jax solver. Note: these must be kept in sync with the versions defined in pybamm/util.py.
jax = [
- "jax>=0.4,<=0.5",
- "jaxlib>=0.4,<=0.5",
+ "jax==0.4.20; python_version >= '3.9'",
+ "jaxlib==0.4.20; python_version >= '3.9'",
]
# For the scikits.odes solver
odes = [
diff --git a/tests/unit/test_experiments/test_experiment_step_termination.py b/tests/unit/test_experiments/test_experiment_step_termination.py
new file mode 100644
index 0000000000..ee45bcc9f8
--- /dev/null
+++ b/tests/unit/test_experiments/test_experiment_step_termination.py
@@ -0,0 +1,18 @@
+#
+# Test the experiment step termination classes
+#
+import pybamm
+import unittest
+
+
+class TestExperimentStepTermination(unittest.TestCase):
+ def test_base_termination(self):
+ term = pybamm.step.BaseTermination(1)
+ self.assertEqual(term.value, 1)
+
+ with self.assertRaises(NotImplementedError):
+ term.get_event(None, None)
+
+ self.assertEqual(term, pybamm.step.BaseTermination(1))
+ self.assertNotEqual(term, pybamm.step.BaseTermination(2))
+ self.assertNotEqual(term, pybamm.step.CurrentTermination(1))
diff --git a/tests/unit/test_experiments/test_experiment_steps.py b/tests/unit/test_experiments/test_experiment_steps.py
index 53b61d637f..b99ae22395 100644
--- a/tests/unit/test_experiments/test_experiment_steps.py
+++ b/tests/unit/test_experiments/test_experiment_steps.py
@@ -34,7 +34,7 @@ def test_step(self):
self.assertEqual(step.type, "voltage")
self.assertEqual(step.value, 1)
self.assertEqual(step.duration, 3600)
- self.assertEqual(step.termination, [{"type": "voltage", "value": 2.5}])
+ self.assertEqual(step.termination, [pybamm.step.VoltageTermination(2.5)])
self.assertEqual(step.period, 60)
self.assertEqual(step.temperature, 298.15)
self.assertEqual(step.tags, ["test"])
@@ -155,25 +155,25 @@ def test_step_string(self):
"type": "C-rate",
"value": -1,
"duration": None,
- "termination": [{"type": "voltage", "value": 4.1}],
+ "termination": [pybamm.step.VoltageTermination(4.1)],
},
{
"value": 4.1,
"type": "voltage",
"duration": None,
- "termination": [{"type": "current", "value": 0.05}],
+ "termination": [pybamm.step.CurrentTermination(0.05)],
},
{
"value": 3,
"type": "voltage",
"duration": None,
- "termination": [{"type": "C-rate", "value": 0.02}],
+ "termination": [pybamm.step.CrateTermination(0.02)],
},
{
"type": "C-rate",
"value": 1 / 3,
"duration": 7200.0,
- "termination": [{"type": "voltage", "value": 2.5}],
+ "termination": [pybamm.step.VoltageTermination(2.5)],
},
]
@@ -264,6 +264,18 @@ def test_start_times(self):
with self.assertRaisesRegex(TypeError, "`start_time` should be"):
pybamm.step._Step("current", 1, duration=3600, start_time="bad start_time")
+ def test_custom_termination(self):
+ def neg_stoich_cutoff(variables):
+ return variables["Negative electrode stoichiometry"] - 1
+
+ neg_stoich_termination = pybamm.step.CustomTermination(
+ name="Negative stoichiometry cut-off", event_function=neg_stoich_cutoff
+ )
+ variables = {"Negative electrode stoichiometry": 3}
+ event = neg_stoich_termination.get_event(variables, None)
+ self.assertEqual(event.name, "Negative stoichiometry cut-off [experiment]")
+ self.assertEqual(event.expression, 2)
+
if __name__ == "__main__":
print("Add -v for more debug output")
diff --git a/tests/unit/test_experiments/test_simulation_with_experiment.py b/tests/unit/test_experiments/test_simulation_with_experiment.py
index 6688fae5b1..cc04177ba2 100644
--- a/tests/unit/test_experiments/test_simulation_with_experiment.py
+++ b/tests/unit/test_experiments/test_simulation_with_experiment.py
@@ -765,6 +765,28 @@ def test_experiment_start_time_identical_steps(self):
# Check that there are only 3 built models (unique steps + padding rest)
self.assertEqual(len(sim.op_conds_to_built_models), 3)
+ def test_experiment_custom_termination(self):
+ def neg_stoich_cutoff(variables):
+ return variables["Negative electrode stoichiometry"] - 0.5
+
+ neg_stoich_termination = pybamm.step.CustomTermination(
+ name="Negative stoichiometry cut-off", event_function=neg_stoich_cutoff
+ )
+
+ model = pybamm.lithium_ion.SPM()
+ experiment = pybamm.Experiment(
+ [pybamm.step.c_rate(1, termination=neg_stoich_termination)]
+ )
+ sim = pybamm.Simulation(model, experiment=experiment)
+ sol = sim.solve(calc_esoh=False)
+ self.assertEqual(
+ sol.cycles[0].steps[0].termination,
+ "event: Negative stoichiometry cut-off [experiment]",
+ )
+
+ neg_stoich = sol["Negative electrode stoichiometry"].data
+ self.assertAlmostEqual(neg_stoich[-1], 0.5, places=4)
+
if __name__ == "__main__":
print("Add -v for more debug output")
diff --git a/tests/unit/test_expression_tree/test_interpolant.py b/tests/unit/test_expression_tree/test_interpolant.py
index 5fa078cffc..1d1a55e85c 100644
--- a/tests/unit/test_expression_tree/test_interpolant.py
+++ b/tests/unit/test_expression_tree/test_interpolant.py
@@ -131,7 +131,7 @@ def f(x, y):
value = interp.evaluate(y=np.array([[1, 1, x[1]], [5, 4, y[1]]]))
np.testing.assert_array_equal(
- value, np.array([[f(1, 5)], [f(1, 4)], [f(x[1], y[1])]])
+ value, np.array([[f(1, 5), f(1, 4), f(x[1], y[1])]])
)
# check also works for cubic
@@ -192,6 +192,17 @@ def f(x, y):
evaluated_children = [1, 4]
value = interp._function_evaluate(evaluated_children)
+ # Test that the interpolant shape is the same as the input data shape
+ interp = pybamm.Interpolant(x_in, data, (var1, var2), interpolator="linear")
+
+ evaluated_children = [np.array([[1, 1]]), np.array([[7, 7]])]
+ value = interp._function_evaluate(evaluated_children)
+ self.assertEqual(value.shape, evaluated_children[0].shape)
+
+ evaluated_children = [np.array([[1, 1], [1, 1]]), np.array([[7, 7], [7, 7]])]
+ value = interp._function_evaluate(evaluated_children)
+ self.assertEqual(value.shape, evaluated_children[0].shape)
+
def test_interpolation_3_x(self):
def f(x, y, z):
return 2 * x**3 + 3 * y**2 - z
@@ -216,7 +227,7 @@ def f(x, y, z):
value = interp.evaluate(y=np.array([[1, 1, 1], [5, 4, 4], [8, 7, 7]]))
np.testing.assert_array_equal(
- value, np.array([[f(1, 5, 8)], [f(1, 4, 7)], [f(1, 4, 7)]])
+ value, np.array([[f(1, 5, 8), f(1, 4, 7), f(1, 4, 7)]])
)
# check also works for cubic
diff --git a/tests/unit/test_expression_tree/test_operations/test_copy.py b/tests/unit/test_expression_tree/test_operations/test_copy.py
index 0340e56bb1..6800f9092f 100644
--- a/tests/unit/test_expression_tree/test_operations/test_copy.py
+++ b/tests/unit/test_expression_tree/test_operations/test_copy.py
@@ -60,6 +60,7 @@ def test_symbol_new_copy(self):
pybamm.maximum(a, b),
pybamm.SparseStack(mat, mat),
pybamm.Equality(a, b),
+ pybamm.EvaluateAt(a, 0),
]:
self.assertEqual(symbol, symbol.new_copy())
diff --git a/tests/unit/test_expression_tree/test_operations/test_evaluate_python.py b/tests/unit/test_expression_tree/test_operations/test_evaluate_python.py
index ca36804ba0..df33e0fe27 100644
--- a/tests/unit/test_expression_tree/test_operations/test_evaluate_python.py
+++ b/tests/unit/test_expression_tree/test_operations/test_evaluate_python.py
@@ -503,7 +503,7 @@ def test_evaluator_jax(self):
expr = pybamm.exp(a * b)
evaluator = pybamm.EvaluatorJax(expr)
result = evaluator(t=None, y=np.array([[2], [3]]))
- self.assertEqual(result, np.exp(6))
+ np.testing.assert_array_almost_equal(result, np.exp(6), decimal=15)
# test a constant expression
expr = pybamm.Scalar(2) * pybamm.Scalar(3)
diff --git a/tests/unit/test_expression_tree/test_operations/test_jac.py b/tests/unit/test_expression_tree/test_operations/test_jac.py
index c6e04d331f..503e7321ea 100644
--- a/tests/unit/test_expression_tree/test_operations/test_jac.py
+++ b/tests/unit/test_expression_tree/test_operations/test_jac.py
@@ -236,6 +236,12 @@ def test_index(self):
jac = ind.jac(vec).evaluate(y=np.linspace(0, 2, 5)).toarray()
np.testing.assert_array_equal(jac, np.array([[0, 0, 0, 0, 0]]))
+ def test_evaluate_at(self):
+ y = pybamm.StateVector(slice(0, 4))
+ expr = pybamm.EvaluateAt(y, 2)
+ jac = expr.jac(y).evaluate(y=np.linspace(0, 2, 4))
+ np.testing.assert_array_equal(jac, 0)
+
def test_jac_of_number(self):
"""Jacobian of a number should be zero"""
a = pybamm.Scalar(1)
diff --git a/tests/unit/test_expression_tree/test_unary_operators.py b/tests/unit/test_expression_tree/test_unary_operators.py
index 3c34db7dcd..6ae6b62d05 100644
--- a/tests/unit/test_expression_tree/test_unary_operators.py
+++ b/tests/unit/test_expression_tree/test_unary_operators.py
@@ -454,6 +454,11 @@ def test_index(self):
pybamm.Index(vec, 5)
pybamm.settings.debug_mode = debug_mode
+ def test_evaluate_at(self):
+ a = pybamm.Symbol("a", domain=["negative electrode"])
+ f = pybamm.EvaluateAt(a, 1)
+ self.assertEqual(f.position, 1)
+
def test_upwind_downwind(self):
# upwind of scalar symbol should fail
a = pybamm.Symbol("a")
@@ -673,9 +678,10 @@ def test_not_constant(self):
self.assertFalse((2 * a).is_constant())
def test_to_equation(self):
-
sympy = have_optional_dependency("sympy")
- sympy_Divergence = have_optional_dependency("sympy.vector.operators", "Divergence")
+ sympy_Divergence = have_optional_dependency(
+ "sympy.vector.operators", "Divergence"
+ )
sympy_Gradient = have_optional_dependency("sympy.vector.operators", "Gradient")
a = pybamm.Symbol("a", domain="negative particle")
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_base_lithium_ion_model.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_base_lithium_ion_model.py
index 315896b29f..fbc916d4a5 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_base_lithium_ion_model.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_base_lithium_ion_model.py
@@ -29,6 +29,25 @@ def test_default_parameters(self):
)
os.chdir(cwd)
+ def test_insert_reference_electrode(self):
+ model = pybamm.lithium_ion.SPM()
+ model.insert_reference_electrode()
+ self.assertIn("Negative electrode 3E potential [V]", model.variables)
+ self.assertIn("Positive electrode 3E potential [V]", model.variables)
+ self.assertIn("Reference electrode potential [V]", model.variables)
+
+ model = pybamm.lithium_ion.SPM({"working electrode": "positive"})
+ model.insert_reference_electrode()
+ self.assertNotIn("Negative electrode potential [V]", model.variables)
+ self.assertIn("Positive electrode 3E potential [V]", model.variables)
+ self.assertIn("Reference electrode potential [V]", model.variables)
+
+ model = pybamm.lithium_ion.SPM({"dimensionality": 2})
+ with self.assertRaisesRegex(
+ NotImplementedError, "Reference electrode can only be"
+ ):
+ model.insert_reference_electrode()
+
if __name__ == "__main__":
print("Add -v for more debug output")
diff --git a/tests/unit/test_parameters/test_parameter_values.py b/tests/unit/test_parameters/test_parameter_values.py
index 37ec89068f..7b5fd38bf0 100644
--- a/tests/unit/test_parameters/test_parameter_values.py
+++ b/tests/unit/test_parameters/test_parameter_values.py
@@ -267,6 +267,15 @@ def test_process_symbol(self):
self.assertEqual(processed_a.value, 4)
self.assertEqual(processed_x, x)
+ # process EvaluateAt
+ evaluate_at = pybamm.EvaluateAt(x, aa)
+ processed_evaluate_at = parameter_values.process_symbol(evaluate_at)
+ self.assertIsInstance(processed_evaluate_at, pybamm.EvaluateAt)
+ self.assertEqual(processed_evaluate_at.children[0], x)
+ self.assertEqual(processed_evaluate_at.position, 4)
+ with self.assertRaisesRegex(ValueError, "'position' in 'EvaluateAt'"):
+ parameter_values.process_symbol(pybamm.EvaluateAt(x, x))
+
# process broadcast
whole_cell = ["negative electrode", "separator", "positive electrode"]
broad = pybamm.PrimaryBroadcast(a, whole_cell)
@@ -559,7 +568,7 @@ def test_process_interpolant_2d(self):
processed_func = parameter_values.process_symbol(func)
self.assertIsInstance(processed_func, pybamm.Interpolant)
self.assertAlmostEqual(
- processed_func.evaluate(inputs={"a": 3.01, "b": 4.4})[0][0], 14.82
+ processed_func.evaluate(inputs={"a": 3.01, "b": 4.4}), 14.82
)
# process differentiated function parameter
diff --git a/tests/unit/test_plotting/test_quick_plot.py b/tests/unit/test_plotting/test_quick_plot.py
index f569f00152..7e2a088de6 100644
--- a/tests/unit/test_plotting/test_quick_plot.py
+++ b/tests/unit/test_plotting/test_quick_plot.py
@@ -294,8 +294,9 @@ def test_spm_simulation(self):
with TemporaryDirectory() as dir_name:
test_stub = os.path.join(dir_name, "spm_sim_test")
test_file = f"{test_stub}.gif"
- quick_plot.create_gif(number_of_images=3, duration=3,
- output_filename=test_file)
+ quick_plot.create_gif(
+ number_of_images=3, duration=3, output_filename=test_file
+ )
assert not os.path.exists(f"{test_stub}*.png")
assert os.path.exists(test_file)
pybamm.close_plots()
@@ -508,6 +509,15 @@ def test_model_with_inputs(self):
pybamm.close_plots()
+class TestQuickPlotAxes(unittest.TestCase):
+ def test_quick_plot_axes(self):
+ axes = pybamm.QuickPlotAxes()
+ axes.add(("test 1", "test 2"), 1)
+ self.assertEqual(axes[0], 1)
+ self.assertEqual(axes.by_variable("test 1"), 1)
+ self.assertEqual(axes.by_variable("test 2"), 1)
+
+
if __name__ == "__main__":
print("Add -v for more debug output")
import sys
diff --git a/tests/unit/test_spatial_methods/test_base_spatial_method.py b/tests/unit/test_spatial_methods/test_base_spatial_method.py
index 37b4eb6a0b..d48ea69a7b 100644
--- a/tests/unit/test_spatial_methods/test_base_spatial_method.py
+++ b/tests/unit/test_spatial_methods/test_base_spatial_method.py
@@ -36,6 +36,8 @@ def test_basics(self):
spatial_method.delta_function(None, None)
with self.assertRaises(NotImplementedError):
spatial_method.internal_neumann_condition(None, None, None, None)
+ with self.assertRaises(NotImplementedError):
+ spatial_method.evaluate_at(None, None, None)
def test_get_auxiliary_domain_repeats(self):
# Test the method to read number of repeats from auxiliary domains
diff --git a/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py b/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py
index 91a5b70044..16a3bbde2c 100644
--- a/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py
+++ b/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py
@@ -551,6 +551,30 @@ def test_full_broadcast_domains(self):
disc = pybamm.Discretisation(mesh, spatial_methods)
disc.process_model(model)
+ def test_evaluate_at(self):
+ mesh = get_p2d_mesh_for_testing()
+ spatial_methods = {
+ "macroscale": pybamm.FiniteVolume(),
+ "negative particle": pybamm.FiniteVolume(),
+ "positive particle": pybamm.FiniteVolume(),
+ }
+ disc = pybamm.Discretisation(mesh, spatial_methods)
+
+ n = mesh["negative electrode"].npts
+ var = pybamm.StateVector(slice(0, n), domain="negative electrode")
+
+ idx = 3
+ position = pybamm.Scalar(mesh["negative electrode"].nodes[idx])
+ evaluate_at = pybamm.EvaluateAt(var, position)
+ evaluate_at_disc = disc.process_symbol(evaluate_at)
+
+ self.assertIsInstance(evaluate_at_disc, pybamm.MatrixMultiplication)
+ self.assertIsInstance(evaluate_at_disc.left, pybamm.Matrix)
+ self.assertIsInstance(evaluate_at_disc.right, pybamm.StateVector)
+
+ y = np.arange(n)[:, np.newaxis]
+ self.assertEqual(evaluate_at_disc.evaluate(y=y), y[idx])
+
if __name__ == "__main__":
print("Add -v for more debug output")
