diff --git a/.all-contributorsrc b/.all-contributorsrc index ab735fcd78..a069737f2d 100644 --- a/.all-contributorsrc +++ b/.all-contributorsrc @@ -1,6 +1,7 @@ { + "skipCi": false, "files": [ - "all_contributors.md", + "all_contributors.md" ], "imageSize": 100, "commit": false, @@ -999,6 +1000,15 @@ "contributions": [ "doc" ] + }, + { + "login": "AndEsterson", + "name": "Andrew Esterson", + "avatar_url": "https://avatars.githubusercontent.com/u/55912083?v=4", + "profile": "https://github.com/AndEsterson", + "contributions": [ + "infra" + ] } ], "contributorsPerLine": 7, @@ -1006,7 +1016,6 @@ "projectOwner": "pybamm-team", "repoType": "github", "repoHost": "https://github.com", - "skipCi": true, "commitConvention": "angular", "commitType": "docs" } diff --git a/.github/workflows/docker.yml b/.github/workflows/docker.yml index 719267912c..f914bc87b3 100644 --- a/.github/workflows/docker.yml +++ b/.github/workflows/docker.yml @@ -21,7 +21,7 @@ jobs: uses: docker/setup-qemu-action@49b3bc8e6bdd4a60e6116a5414239cba5943d3cf # v3.2.0 - name: Set up Docker Buildx - uses: docker/setup-buildx-action@c47758b77c9736f4b2ef4073d4d51994fabfe349 # v3.7.1 + uses: docker/setup-buildx-action@6524bf65af31da8d45b59e8c27de4bd072b392f5 # v3.8.0 - name: Login to Docker Hub uses: docker/login-action@9780b0c442fbb1117ed29e0efdff1e18412f7567 # v3.3.0 diff --git a/.github/workflows/scorecard.yml b/.github/workflows/scorecard.yml index 03822b3bd6..975be1bc1e 100644 --- a/.github/workflows/scorecard.yml +++ b/.github/workflows/scorecard.yml @@ -68,6 +68,6 @@ jobs: # Upload the results to GitHub's code scanning dashboard (optional). # Commenting out will disable upload of results to your repo's Code Scanning dashboard - name: "Upload to code-scanning" - uses: github/codeql-action/upload-sarif@aa578102511db1f4524ed59b8cc2bae4f6e88195 # v3.27.6 + uses: github/codeql-action/upload-sarif@df409f7d9260372bd5f19e5b04e83cb3c43714ae # v3.27.9 with: sarif_file: results.sarif diff --git a/.github/workflows/test_on_push.yml b/.github/workflows/test_on_push.yml index 1e759b0f2a..3443279414 100644 --- a/.github/workflows/test_on_push.yml +++ b/.github/workflows/test_on_push.yml @@ -3,7 +3,6 @@ name: PyBaMM on: workflow_dispatch: pull_request: - paths-ignore: [".all-contributorsrc", "all_contributors.md"] env: PYBAMM_DISABLE_TELEMETRY: "true" diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 4c3ef44e6e..a2e599cc00 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -4,7 +4,7 @@ ci: repos: - repo: https://github.com/astral-sh/ruff-pre-commit - rev: "v0.8.2" + rev: "v0.8.3" hooks: - id: ruff args: [--fix, --show-fixes] diff --git a/CHANGELOG.md b/CHANGELOG.md index 920a292683..8322884eda 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,10 +1,15 @@ # [Unreleased](https://github.com/pybamm-team/PyBaMM/) ## Features + +- Made composite electrode model compatible with particle size distribution ([#4687](https://github.com/pybamm-team/PyBaMM/pull/4687)) +- Added `Symbol.post_order()` method to return an iterable that steps through the tree in post-order fashion. ([#4684](https://github.com/pybamm-team/PyBaMM/pull/4684)) - Added two more submodels (options) for the SEI: Lars von Kolzenberg (2020) model and Tunneling Limit model ([#4394](https://github.com/pybamm-team/PyBaMM/pull/4394)) + ## Breaking changes +- Summary variables now calculated only when called, accessed via a class in the same manner as other variables rather than a dictionary. ([#4621](https://github.com/pybamm-team/PyBaMM/pull/4621)) - The conda distribution (`pybamm`) now installs all optional dependencies available on conda-forge. Use the new `pybamm-base` conda package to install PyBaMM with only the required dependencies. ([conda-forge/pybamm-feedstock#70](https://github.com/conda-forge/pybamm-feedstock/pull/70)) - Separated extrapolation options for `pybamm.BoundaryValue` and `pybamm.BoundaryGradient`, and updated the default to be "linear" for the value and "quadratic" for the gradient. ([#4614](https://github.com/pybamm-team/PyBaMM/pull/4614)) diff --git a/all_contributors.md b/all_contributors.md index 61434b6750..2d6cdb10b4 100644 --- a/all_contributors.md +++ b/all_contributors.md @@ -1,6 +1,6 @@ -[![All Contributors](https://img.shields.io/badge/all_contributors-94-orange.svg)](#-contributors) +[![All Contributors](https://img.shields.io/badge/all_contributors-95-orange.svg)](#-contributors) Thanks goes to these wonderful people ([emoji key](https://allcontributors.org/docs/en/emoji-key)): @@ -68,7 +68,7 @@ Thanks goes to these wonderful people ([emoji key](https://allcontributors.org/d bardsleypt
bardsleypt
πŸ› πŸ’» ndrewwang
ndrewwang
πŸ› πŸ’» MichaPhilipp
MichaPhilipp
πŸ› - Alec Bills
Alec Bills
πŸ’» + Alec Bills
Alec Bills
πŸ’» Agriya Khetarpal
Agriya Khetarpal
πŸš‡ πŸ’» πŸ“– πŸ‘€ Alex Wadell
Alex Wadell
πŸ’» ⚠️ πŸ“– iatzak
iatzak
πŸ“– πŸ› πŸ’» @@ -131,6 +131,7 @@ Thanks goes to these wonderful people ([emoji key](https://allcontributors.org/d Marc Berliner
Marc Berliner
πŸ’» πŸ“– πŸš‡ 🚧 Aswinr24
Aswinr24
⚠️ isaacbasil
isaacbasil
πŸ“– + Andrew Esterson
Andrew Esterson
πŸš‡ diff --git a/docs/source/api/solvers/index.rst b/docs/source/api/solvers/index.rst index a9aa8ac1dd..2bb15503bc 100644 --- a/docs/source/api/solvers/index.rst +++ b/docs/source/api/solvers/index.rst @@ -13,3 +13,4 @@ Solvers algebraic_solvers solution processed_variable + summary_variables diff --git a/docs/source/api/solvers/summary_variables.rst b/docs/source/api/solvers/summary_variables.rst new file mode 100644 index 0000000000..60344660c4 --- /dev/null +++ b/docs/source/api/solvers/summary_variables.rst @@ -0,0 +1,5 @@ +Summary Variables +====================== + +.. autoclass:: pybamm.SummaryVariables + :members: diff --git a/docs/source/examples/notebooks/models/loss_of_active_materials.ipynb b/docs/source/examples/notebooks/models/loss_of_active_materials.ipynb index 54ace50fc6..249ae509a6 100644 --- a/docs/source/examples/notebooks/models/loss_of_active_materials.ipynb +++ b/docs/source/examples/notebooks/models/loss_of_active_materials.ipynb @@ -443,6 +443,8 @@ "second = 0.1\n", "parameter_values.update(\n", " {\n", + " \"Primary: Negative electrode reference concentration for free of deformation [mol.m-3]\": 0.0,\n", + " \"Secondary: Negative electrode reference concentration for free of deformation [mol.m-3]\": 0.0,\n", " \"Primary: Negative electrode LAM constant proportional term [s-1]\": 1e-4 / 3600,\n", " \"Secondary: Negative electrode LAM constant proportional term [s-1]\": 1e-4\n", " / 3600\n", diff --git a/docs/source/examples/notebooks/models/using-submodels.ipynb b/docs/source/examples/notebooks/models/using-submodels.ipynb index cce32fcefd..db889972f8 100644 --- a/docs/source/examples/notebooks/models/using-submodels.ipynb +++ b/docs/source/examples/notebooks/models/using-submodels.ipynb @@ -527,10 +527,10 @@ " model.param, \"positive\", model.options, cracks=True\n", ")\n", "model.submodels[\"Negative lithium plating\"] = pybamm.lithium_plating.NoPlating(\n", - " model.param, \"Negative\"\n", + " model.param, \"negative\"\n", ")\n", "model.submodels[\"Positive lithium plating\"] = pybamm.lithium_plating.NoPlating(\n", - " model.param, \"Positive\"\n", + " model.param, \"positive\"\n", ")" ] }, diff --git a/docs/source/examples/notebooks/simulations_and_experiments/simulating-long-experiments.ipynb b/docs/source/examples/notebooks/simulations_and_experiments/simulating-long-experiments.ipynb index c7f1f0e634..3dcb46b1b5 100644 --- a/docs/source/examples/notebooks/simulations_and_experiments/simulating-long-experiments.ipynb +++ b/docs/source/examples/notebooks/simulations_and_experiments/simulating-long-experiments.ipynb @@ -541,7 +541,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": null, "id": "right-skiing", "metadata": {}, "outputs": [ @@ -638,7 +638,7 @@ } ], "source": [ - "sorted(sol.summary_variables.keys())" + "sorted(sol.summary_variables.all_variables)" ] }, { @@ -1936,7 +1936,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "venv", "language": "python", "name": "python3" }, @@ -1950,7 +1950,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.10" + "version": "3.11.10" }, "toc": { "base_numbering": 1, @@ -1964,11 +1964,6 @@ "toc_position": {}, "toc_section_display": true, "toc_window_display": true - }, - "vscode": { - "interpreter": { - "hash": "612adcc456652826e82b485a1edaef831aa6d5abc680d008e93d513dd8724f14" - } } }, "nbformat": 4, diff --git a/docs/source/user_guide/installation/index.rst b/docs/source/user_guide/installation/index.rst index 9fb81c8532..4c493d6ee4 100644 --- a/docs/source/user_guide/installation/index.rst +++ b/docs/source/user_guide/installation/index.rst @@ -63,6 +63,8 @@ Package Minimum supp `pandas `__ 1.5.0 `pooch `__ 1.8.1 `posthog `__ 3.6.5 +`pyyaml `__ +`platformdirs `__ =================================================================== ========================== .. _install.optional_dependencies: @@ -102,6 +104,7 @@ Installable with ``pip install "pybamm[docs]"`` Dependency Minimum Version pip extra Notes ================================================================================================= ================== ================== ======================================================================= `sphinx `__ \- docs Sphinx makes it easy to create intelligent and beautiful documentation. +`sphinx_rtd_theme `__ \- docs This Sphinx theme provides a great reader experience for documentation. `pydata-sphinx-theme `__ \- docs A clean, Bootstrap-based Sphinx theme. `sphinx_design `__ \- docs A sphinx extension for designing. `sphinx-copybutton `__ \- docs To copy codeblocks. @@ -109,6 +112,13 @@ Dependency `sphinx-inline-tabs `__ \- docs Add inline tabbed content to your Sphinx documentation. `sphinxcontrib-bibtex `__ \- docs For BibTeX citations. `sphinx-autobuild `__ \- docs For re-building docs once triggered. +`sphinx-last-updated-by-git `__ \- docs To get the "last updated" time for each Sphinx page from Git. +`nbsphinx `__ \- docs Sphinx extension that provides a source parser for .ipynb files +`ipykernel `__ \- docs Provides the IPython kernel for Jupyter. +`ipywidgets `__ \- docs Interactive HTML widgets for Jupyter notebooks and the IPython kernel. +`sphinx-gallery `__ \- docs Builds an HTML gallery of examples from any set of Python scripts. +`sphinx-hoverxref `__ \- docs Sphinx extension to show a floating window. +`sphinx-docsearch `__ \- docs To replaces Sphinx’s built-in search with Algolia DocSearch. ================================================================================================= ================== ================== ======================================================================= .. _install.examples_dependencies: @@ -142,7 +152,9 @@ Dependency `pytest `__ 6.0.0 dev For running the test suites. `pytest-doctestplus `__ \- dev For running doctests. `pytest-xdist `__ \- dev For running tests in parallel across distributed workers. +`pytest-mock `__ \- dev Provides a mocker fixture. `nbmake `__ \- dev A ``pytest`` plugin for executing Jupyter notebooks. +`importlib-metadata `__ \- dev Used to read metadata from Python packages. ================================================================================ ================== ================== ============================================================= .. _install.cite_dependencies: diff --git a/pyproject.toml b/pyproject.toml index 923f72438d..9d5d788129 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -91,7 +91,7 @@ cite = [ ] # Battery Parameter eXchange format bpx = [ - "bpx>=0.4.0", + "bpx==0.4.0", ] # Low-overhead progress bars tqdm = [ diff --git a/src/pybamm/__init__.py b/src/pybamm/__init__.py index e4ab341c1e..ef9b10e384 100644 --- a/src/pybamm/__init__.py +++ b/src/pybamm/__init__.py @@ -163,6 +163,7 @@ from .solvers.processed_variable_time_integral import ProcessedVariableTimeIntegral from .solvers.processed_variable import ProcessedVariable, process_variable from .solvers.processed_variable_computed import ProcessedVariableComputed +from .solvers.summary_variable import SummaryVariables from .solvers.base_solver import BaseSolver from .solvers.dummy_solver import DummySolver from .solvers.algebraic_solver import AlgebraicSolver diff --git a/src/pybamm/callbacks.py b/src/pybamm/callbacks.py index 4e8c67c8be..29286ab5a3 100644 --- a/src/pybamm/callbacks.py +++ b/src/pybamm/callbacks.py @@ -212,7 +212,7 @@ def on_cycle_end(self, logs): voltage_stop = logs["stopping conditions"]["voltage"] if voltage_stop is not None: - min_voltage = logs["summary variables"]["Minimum voltage [V]"] + min_voltage = logs["Minimum voltage [V]"] if min_voltage > voltage_stop[0]: self.logger.notice( f"Minimum voltage is now {min_voltage:.3f} V " diff --git a/src/pybamm/expression_tree/averages.py b/src/pybamm/expression_tree/averages.py index 11538ea153..c95ccb4e5e 100644 --- a/src/pybamm/expression_tree/averages.py +++ b/src/pybamm/expression_tree/averages.py @@ -306,6 +306,10 @@ def xyz_average(symbol: pybamm.Symbol) -> pybamm.Symbol: return yz_average(x_average(symbol)) +def xyzs_average(symbol: pybamm.Symbol) -> pybamm.Symbol: + return xyz_average(size_average(symbol)) + + def r_average(symbol: pybamm.Symbol) -> pybamm.Symbol: """ Convenience function for creating an average in the r-direction. @@ -373,7 +377,15 @@ def size_average( # If symbol doesn't have a domain, or doesn't have "negative particle size" # or "positive particle size" as a domain, it's average value is itself if symbol.domain == [] or not any( - domain in [["negative particle size"], ["positive particle size"]] + domain + in [ + ["negative particle size"], + ["positive particle size"], + ["negative primary particle size"], + ["positive primary particle size"], + ["negative secondary particle size"], + ["positive secondary particle size"], + ] for domain in list(symbol.domains.values()) ): return symbol @@ -394,13 +406,30 @@ def size_average( else: if f_a_dist is None: geo = pybamm.geometric_parameters + name = "R" + if "negative" in symbol.domain[0]: + name += "_n" + elif "positive" in symbol.domain[0]: + name += "_p" + if "primary" in symbol.domain[0]: + name += "_prim" + elif "secondary" in symbol.domain[0]: + name += "_sec" R = pybamm.SpatialVariable( - "R", domains=symbol.domains, coord_sys="cartesian" + name, domains=symbol.domains, coord_sys="cartesian" ) - if ["negative particle size"] in symbol.domains.values(): + if ["negative particle size"] in symbol.domains.values() or [ + "negative primary particle size" + ] in symbol.domains.values(): f_a_dist = geo.n.prim.f_a_dist(R) - elif ["positive particle size"] in symbol.domains.values(): + elif ["negative secondary particle size"] in symbol.domains.values(): + f_a_dist = geo.n.sec.f_a_dist(R) + elif ["positive particle size"] in symbol.domains.values() or [ + "positive primary particle size" + ] in symbol.domains.values(): f_a_dist = geo.p.prim.f_a_dist(R) + elif ["positive secondary particle size"] in symbol.domains.values(): + f_a_dist = geo.p.sec.f_a_dist(R) return SizeAverage(symbol, f_a_dist) diff --git a/src/pybamm/expression_tree/symbol.py b/src/pybamm/expression_tree/symbol.py index c4e0043f17..087acf8986 100644 --- a/src/pybamm/expression_tree/symbol.py +++ b/src/pybamm/expression_tree/symbol.py @@ -570,6 +570,13 @@ def pre_order(self): anytree = import_optional_dependency("anytree") return anytree.PreOrderIter(self) + def post_order(self): + """ + returns an iterable that steps through the tree in post-order fashion. + """ + anytree = import_optional_dependency("anytree") + return anytree.PostOrderIter(self) + def __str__(self): """return a string representation of the node and its children.""" return self._name diff --git a/src/pybamm/geometry/battery_geometry.py b/src/pybamm/geometry/battery_geometry.py index 9be08ff619..3e1986507a 100644 --- a/src/pybamm/geometry/battery_geometry.py +++ b/src/pybamm/geometry/battery_geometry.py @@ -83,6 +83,18 @@ def battery_geometry( "negative particle size": {"R_n": {"min": R_min_n, "max": R_max_n}}, } ) + phases = int(options.negative["particle phases"]) + if phases >= 2: + geometry.update( + { + "negative primary particle size": { + "R_n_prim": {"min": R_min_n, "max": R_max_n} + }, + "negative secondary particle size": { + "R_n_sec": {"min": R_min_n, "max": R_max_n} + }, + } + ) if ( options is not None and options.positive["particle size"] == "distribution" @@ -95,6 +107,18 @@ def battery_geometry( "positive particle size": {"R_p": {"min": R_min_p, "max": R_max_p}}, } ) + phases = int(options.positive["particle phases"]) + if phases >= 2: + geometry.update( + { + "positive primary particle size": { + "R_p_prim": {"min": R_min_p, "max": R_max_p} + }, + "positive secondary particle size": { + "R_p_sec": {"min": R_min_p, "max": R_max_p} + }, + } + ) # Add current collector domains current_collector_dimension = options["dimensionality"] diff --git a/src/pybamm/geometry/standard_spatial_vars.py b/src/pybamm/geometry/standard_spatial_vars.py index 565908deb9..b6638694f7 100644 --- a/src/pybamm/geometry/standard_spatial_vars.py +++ b/src/pybamm/geometry/standard_spatial_vars.py @@ -98,6 +98,7 @@ }, coord_sys="cartesian", ) + R_p = pybamm.SpatialVariable( "R_p", domain=["positive particle size"], @@ -108,6 +109,44 @@ coord_sys="cartesian", ) +R_n_prim = pybamm.SpatialVariable( + "R_n_prim", + domain=["negative primary particle size"], + auxiliary_domains={ + "secondary": "negative electrode", + "tertiary": "current collector", + }, + coord_sys="cartesian", +) +R_p_prim = pybamm.SpatialVariable( + "R_p_prim", + domain=["positive primary particle size"], + auxiliary_domains={ + "secondary": "positive electrode", + "tertiary": "current collector", + }, + coord_sys="cartesian", +) + +R_n_sec = pybamm.SpatialVariable( + "R_n_sec", + domain=["negative secondary particle size"], + auxiliary_domains={ + "secondary": "negative electrode", + "tertiary": "current collector", + }, + coord_sys="cartesian", +) +R_p_sec = pybamm.SpatialVariable( + "R_p_sec", + domain=["positive secondary particle size"], + auxiliary_domains={ + "secondary": "positive electrode", + "tertiary": "current collector", + }, + coord_sys="cartesian", +) + # Domains at cell edges x_n_edge = pybamm.SpatialVariableEdge( "x_n", diff --git a/src/pybamm/models/full_battery_models/base_battery_model.py b/src/pybamm/models/full_battery_models/base_battery_model.py index fd404cec75..83d62df278 100644 --- a/src/pybamm/models/full_battery_models/base_battery_model.py +++ b/src/pybamm/models/full_battery_models/base_battery_model.py @@ -525,14 +525,14 @@ def __init__(self, extra_options): # Options not yet compatible with particle-size distributions if options["particle size"] == "distribution": - if options["heat of mixing"] != "false": + if options["lithium plating porosity change"] != "false": raise NotImplementedError( - "Heat of mixing submodels do not yet support particle-size " - "distributions." + "Lithium plating porosity change not yet supported for particle-size" + " distributions." ) - if options["lithium plating"] != "none": + if options["heat of mixing"] != "false": raise NotImplementedError( - "Lithium plating submodels do not yet support particle-size " + "Heat of mixing submodels do not yet support particle-size " "distributions." ) if options["particle"] in ["quadratic profile", "quartic profile"]: @@ -623,7 +623,6 @@ def __init__(self, extra_options): if options["particle phases"] not in ["1", ("1", "1")]: if not ( options["surface form"] != "false" - and options["particle size"] == "single" and options["particle"] == "Fickian diffusion" ): raise pybamm.OptionError( @@ -868,6 +867,10 @@ def default_var_pts(self): "z": 10, "R_n": 30, "R_p": 30, + "R_n_prim": 30, + "R_p_prim": 30, + "R_n_sec": 30, + "R_p_sec": 30, } # Reduce the default points for 2D current collectors if self.options["dimensionality"] == 2: @@ -888,6 +891,10 @@ def default_submesh_types(self): "positive secondary particle": pybamm.Uniform1DSubMesh, "negative particle size": pybamm.Uniform1DSubMesh, "positive particle size": pybamm.Uniform1DSubMesh, + "negative primary particle size": pybamm.Uniform1DSubMesh, + "positive primary particle size": pybamm.Uniform1DSubMesh, + "negative secondary particle size": pybamm.Uniform1DSubMesh, + "positive secondary particle size": pybamm.Uniform1DSubMesh, } if self.options["dimensionality"] == 0: base_submeshes["current collector"] = pybamm.SubMesh0D @@ -910,6 +917,10 @@ def default_spatial_methods(self): "positive secondary particle": pybamm.FiniteVolume(), "negative particle size": pybamm.FiniteVolume(), "positive particle size": pybamm.FiniteVolume(), + "negative primary particle size": pybamm.FiniteVolume(), + "positive primary particle size": pybamm.FiniteVolume(), + "negative secondary particle size": pybamm.FiniteVolume(), + "positive secondary particle size": pybamm.FiniteVolume(), } if self.options["dimensionality"] == 0: # 0D submesh - use base spatial method diff --git a/src/pybamm/models/full_battery_models/lithium_ion/electrode_soh.py b/src/pybamm/models/full_battery_models/lithium_ion/electrode_soh.py index a743910905..cdcb9a55fa 100644 --- a/src/pybamm/models/full_battery_models/lithium_ion/electrode_soh.py +++ b/src/pybamm/models/full_battery_models/lithium_ion/electrode_soh.py @@ -312,6 +312,27 @@ def __init__( self.__get_electrode_soh_sims_split ) + def __getstate__(self): + """ + Return dictionary of picklable items + """ + result = self.__dict__.copy() + result["_get_electrode_soh_sims_full"] = None # Exclude LRU cache + result["_get_electrode_soh_sims_split"] = None # Exclude LRU cache + return result + + def __setstate__(self, state): + """ + Unpickle, restoring unpicklable relationships + """ + self.__dict__ = state + self._get_electrode_soh_sims_full = lru_cache()( + self.__get_electrode_soh_sims_full + ) + self._get_electrode_soh_sims_split = lru_cache()( + self.__get_electrode_soh_sims_split + ) + def _get_lims_ocp(self): parameter_values = self.parameter_values diff --git a/src/pybamm/models/submodels/active_material/base_active_material.py b/src/pybamm/models/submodels/active_material/base_active_material.py index ea0e826e09..4827c1f6d4 100644 --- a/src/pybamm/models/submodels/active_material/base_active_material.py +++ b/src/pybamm/models/submodels/active_material/base_active_material.py @@ -81,9 +81,19 @@ def _get_standard_active_material_variables(self, eps_solid): R = self.phase_param.R elif domain_options["particle size"] == "distribution": if self.domain == "negative": - R_ = pybamm.standard_spatial_vars.R_n + if self.phase_name == "primary ": + R_ = pybamm.standard_spatial_vars.R_n_prim + elif self.phase_name == "secondary ": + R_ = pybamm.standard_spatial_vars.R_n_sec + else: + R_ = pybamm.standard_spatial_vars.R_n elif self.domain == "positive": - R_ = pybamm.standard_spatial_vars.R_p + if self.phase_name == "primary ": + R_ = pybamm.standard_spatial_vars.R_p_prim + elif self.phase_name == "secondary ": + R_ = pybamm.standard_spatial_vars.R_p_sec + else: + R_ = pybamm.standard_spatial_vars.R_p R = pybamm.size_average(R_) R_av = pybamm.x_average(R) diff --git a/src/pybamm/models/submodels/interface/base_interface.py b/src/pybamm/models/submodels/interface/base_interface.py index 0ad08d5454..229aa5fdf2 100644 --- a/src/pybamm/models/submodels/interface/base_interface.py +++ b/src/pybamm/models/submodels/interface/base_interface.py @@ -46,6 +46,8 @@ def __init__(self, param, domain, reaction, options, phase="primary"): self.reaction_name = self.phase_name + self.reaction_name self.reaction = reaction + domain_options = getattr(self.options, domain) + self.size_distribution = domain_options["particle size"] == "distribution" def _get_exchange_current_density(self, variables): """ @@ -93,8 +95,10 @@ def _get_exchange_current_density(self, variables): # "current collector" c_e = pybamm.PrimaryBroadcast(c_e, ["current collector"]) # broadcast c_e, T onto "particle size" - c_e = pybamm.PrimaryBroadcast(c_e, [f"{domain} particle size"]) - T = pybamm.PrimaryBroadcast(T, [f"{domain} particle size"]) + c_e = pybamm.PrimaryBroadcast( + c_e, [f"{domain} {phase_name}particle size"] + ) + T = pybamm.PrimaryBroadcast(T, [f"{domain} {phase_name}particle size"]) else: c_s_surf = variables[ @@ -215,13 +219,23 @@ def _get_standard_interfacial_current_variables(self, j): # Size average. For j variables that depend on particle size, see # "_get_standard_size_distribution_interfacial_current_variables" - if j.domain in [["negative particle size"], ["positive particle size"]]: + if j.domain in [ + ["negative particle size"], + ["positive particle size"], + ["negative primary particle size"], + ["positive primary particle size"], + ["negative secondary particle size"], + ["positive secondary particle size"], + ]: j = pybamm.size_average(j) # Average, and broadcast if necessary j_av = pybamm.x_average(j) if j.domain == []: j = pybamm.FullBroadcast(j, f"{domain} electrode", "current collector") - elif j.domain == ["current collector"]: + elif ( + j.domain == ["current collector"] + and getattr(self, "reaction_loc", None) != "interface" + ): j = pybamm.PrimaryBroadcast(j, f"{domain} electrode") variables = { @@ -230,6 +244,13 @@ def _get_standard_interfacial_current_variables(self, j): f"X-averaged {domain} electrode {reaction_name}" "interfacial current density [A.m-2]": j_av, } + if self.phase_name == reaction_name: + variables.update( + { + f"{Domain} electrode {reaction_name}interfacial current density [A.m-2]": j, + f"X-averaged {domain} electrode {reaction_name}interfacial current density [A.m-2]": j_av, + } + ) return variables @@ -298,7 +319,6 @@ def _get_standard_volumetric_current_density_variables(self, variables): a = variables[ f"{Domain} electrode {phase_name}surface area to volume ratio [m-1]" ] - j = variables[ f"{Domain} electrode {reaction_name}interfacial current density [A.m-2]" ] @@ -306,9 +326,10 @@ def _get_standard_volumetric_current_density_variables(self, variables): a_j_av = pybamm.x_average(a_j) if reaction_name == "SEI on cracks ": - roughness = variables[f"{Domain} electrode roughness ratio"] - 1 + roughness = variables[f"{Domain} {phase_name}electrode roughness ratio"] - 1 roughness_av = ( - variables[f"X-averaged {domain} electrode roughness ratio"] - 1 + variables[f"X-averaged {domain} {phase_name}electrode roughness ratio"] + - 1 ) else: roughness = 1 @@ -335,7 +356,14 @@ def _get_standard_overpotential_variables(self, eta_r): # Size average. For eta_r variables that depend on particle size, see # "_get_standard_size_distribution_overpotential_variables" - if eta_r.domain in [["negative particle size"], ["positive particle size"]]: + if eta_r.domain in [ + ["negative particle size"], + ["positive particle size"], + ["negative primary particle size"], + ["positive primary particle size"], + ["negative secondary particle size"], + ["positive secondary particle size"], + ]: eta_r = pybamm.size_average(eta_r) # X-average, and broadcast if necessary @@ -410,9 +438,11 @@ def _get_standard_size_distribution_interfacial_current_variables(self, j): if j.domains["secondary"] == [f"{domain} electrode"]: # x-average j_xav = pybamm.x_average(j) - else: + elif getattr(self, "reaction_loc", None) != "interface": j_xav = j j = pybamm.SecondaryBroadcast(j_xav, [f"{domain} electrode"]) + else: + j_xav = j variables = { f"{Domain} electrode {reaction_name}" diff --git a/src/pybamm/models/submodels/interface/kinetics/base_kinetics.py b/src/pybamm/models/submodels/interface/kinetics/base_kinetics.py index 7cfa83631b..c3b6ed5329 100644 --- a/src/pybamm/models/submodels/interface/kinetics/base_kinetics.py +++ b/src/pybamm/models/submodels/interface/kinetics/base_kinetics.py @@ -76,7 +76,9 @@ def get_coupled_variables(self, variables): self.reaction == "lithium-ion main" and domain_options["particle size"] == "distribution" ): - delta_phi = pybamm.PrimaryBroadcast(delta_phi, [f"{domain} particle size"]) + delta_phi = pybamm.PrimaryBroadcast( + delta_phi, [f"{domain} {phase_name}particle size"] + ) # Get exchange-current density. For MSMR models we calculate the exchange # current density for each reaction, then sum these to give a total exchange @@ -169,7 +171,7 @@ def get_coupled_variables(self, variables): elif j0.domain in ["current collector", ["current collector"]]: T = variables["X-averaged cell temperature [K]"] u = variables[f"X-averaged {domain} electrode interface utilisation"] - elif j0.domain == [f"{domain} particle size"]: + elif j0.domain == [f"{domain} {phase_name}particle size"]: if j0.domains["secondary"] != [f"{domain} electrode"]: T = variables["X-averaged cell temperature [K]"] u = variables[f"X-averaged {domain} electrode interface utilisation"] @@ -178,7 +180,7 @@ def get_coupled_variables(self, variables): u = variables[f"{Domain} electrode interface utilisation"] # Broadcast T onto "particle size" domain - T = pybamm.PrimaryBroadcast(T, [f"{domain} particle size"]) + T = pybamm.PrimaryBroadcast(T, [f"{domain} {phase_name}particle size"]) else: T = variables[f"{Domain} electrode temperature [K]"] u = variables[f"{Domain} electrode interface utilisation"] @@ -198,7 +200,9 @@ def get_coupled_variables(self, variables): else: j = self._get_kinetics(j0, ne, eta_r, T, u) - if j.domain == [f"{domain} particle size"]: + if j.domain == [f"{domain} particle size"] or j.domain == [ + f"{domain} {phase_name}particle size" + ]: # If j depends on particle size, get size-dependent "distribution" # variables first variables.update( diff --git a/src/pybamm/models/submodels/interface/lithium_plating/base_plating.py b/src/pybamm/models/submodels/interface/lithium_plating/base_plating.py index 0530476f9d..c34eca9e05 100644 --- a/src/pybamm/models/submodels/interface/lithium_plating/base_plating.py +++ b/src/pybamm/models/submodels/interface/lithium_plating/base_plating.py @@ -64,7 +64,9 @@ def _get_standard_concentration_variables(self, c_plated_Li, c_dead_Li): phase_name = self.phase_name phase_param = self.phase_param domain, Domain = self.domain_Domain - + if self.size_distribution is True: + c_plated_Li = pybamm.size_average(c_plated_Li) + c_dead_Li = pybamm.size_average(c_dead_Li) # Set scales to one for the "no plating" model so that they are not required # by parameter values in general if isinstance(self, pybamm.lithium_plating.NoPlating): @@ -137,3 +139,24 @@ def _get_standard_reaction_variables(self, j_stripping): } return variables + + def _get_standard_size_distribution_reaction_variables(self, j_stripping): + """ + A private function to obtain the standard variables which + can be derived from the lithum stripping interfacial reaction current + """ + domain, Domain = self.domain_Domain + phase_name = self.phase_name + j_stripping_sav = pybamm.size_average(j_stripping) + j_stripping_av = pybamm.x_average(j_stripping_sav) + + variables = { + f"{Domain} electrode {phase_name}lithium plating " + "interfacial current density distribution [A.m-2]": j_stripping, + f"{Domain} electrode {phase_name}lithium plating " + "interfacial current density [A.m-2]": j_stripping_sav, + f"X-averaged {domain} electrode {phase_name}lithium plating " + "interfacial current density [A.m-2]": j_stripping_av, + } + + return variables diff --git a/src/pybamm/models/submodels/interface/lithium_plating/no_plating.py b/src/pybamm/models/submodels/interface/lithium_plating/no_plating.py index ddc13cfb97..8efa664a92 100644 --- a/src/pybamm/models/submodels/interface/lithium_plating/no_plating.py +++ b/src/pybamm/models/submodels/interface/lithium_plating/no_plating.py @@ -20,12 +20,32 @@ def __init__(self, param, domain, options=None, phase="primary"): super().__init__(param, domain, options=options, phase=phase) def get_fundamental_variables(self): - zero = pybamm.FullBroadcast( - pybamm.Scalar(0), f"{self.domain} electrode", "current collector" - ) + phase_name = self.phase_name + if self.size_distribution is False: + zero = pybamm.FullBroadcast( + pybamm.Scalar(0), f"{self.domain} electrode", "current collector" + ) + else: + zero = pybamm.FullBroadcast( + pybamm.Scalar(0), + f"{self.domain} {phase_name}particle size", + { + "secondary": f"{self.domain} electrode", + "tertiary": "current collector", + }, + ) variables = self._get_standard_concentration_variables(zero, zero) + if self.size_distribution: + variables.update( + self._get_standard_size_distribution_overpotential_variables(zero) + ) + variables.update( + self._get_standard_size_distribution_reaction_variables(zero) + ) + else: + variables.update(self._get_standard_reaction_variables(zero)) variables.update(self._get_standard_overpotential_variables(zero)) - variables.update(self._get_standard_reaction_variables(zero)) + return variables def get_coupled_variables(self, variables): diff --git a/src/pybamm/models/submodels/interface/open_circuit_potential/base_ocp.py b/src/pybamm/models/submodels/interface/open_circuit_potential/base_ocp.py index 309eb343c0..e2c4a43551 100644 --- a/src/pybamm/models/submodels/interface/open_circuit_potential/base_ocp.py +++ b/src/pybamm/models/submodels/interface/open_circuit_potential/base_ocp.py @@ -49,7 +49,14 @@ def _get_standard_ocp_variables(self, ocp_surf, ocp_bulk, dUdT): reaction_name = self.reaction_name # Update size variables then size average. - if ocp_surf.domain in [["negative particle size"], ["positive particle size"]]: + if ocp_surf.domain in [ + ["negative particle size"], + ["positive particle size"], + ["negative primary particle size"], + ["positive primary particle size"], + ["negative secondary particle size"], + ["positive secondary particle size"], + ]: variables = self._get_standard_size_distribution_ocp_variables( ocp_surf, dUdT ) diff --git a/src/pybamm/models/submodels/interface/open_circuit_potential/current_sigmoid_ocp.py b/src/pybamm/models/submodels/interface/open_circuit_potential/current_sigmoid_ocp.py index fa55ba1079..8fc41d1711 100644 --- a/src/pybamm/models/submodels/interface/open_circuit_potential/current_sigmoid_ocp.py +++ b/src/pybamm/models/submodels/interface/open_circuit_potential/current_sigmoid_ocp.py @@ -36,7 +36,7 @@ def get_coupled_variables(self, variables): ): sto_surf = sto_surf.orphans[0] T = T.orphans[0] - T = pybamm.PrimaryBroadcast(T, [f"{domain} particle size"]) + T = pybamm.PrimaryBroadcast(T, [f"{domain} {phase_name}particle size"]) else: sto_surf = variables[ f"{Domain} {phase_name}particle surface stoichiometry" diff --git a/src/pybamm/models/submodels/interface/open_circuit_potential/single_ocp.py b/src/pybamm/models/submodels/interface/open_circuit_potential/single_ocp.py index 932e2e3eff..42a4b53ebb 100644 --- a/src/pybamm/models/submodels/interface/open_circuit_potential/single_ocp.py +++ b/src/pybamm/models/submodels/interface/open_circuit_potential/single_ocp.py @@ -25,7 +25,7 @@ def get_coupled_variables(self, variables): ): sto_surf = sto_surf.orphans[0] T = T.orphans[0] - T = pybamm.PrimaryBroadcast(T, [f"{domain} particle size"]) + T = pybamm.PrimaryBroadcast(T, [f"{domain} {phase_name}particle size"]) else: sto_surf = variables[ f"{Domain} {phase_name}particle surface stoichiometry" @@ -42,7 +42,7 @@ def get_coupled_variables(self, variables): # Bulk OCP is from the average SOC and temperature sto_bulk = variables[f"{Domain} electrode {phase_name}stoichiometry"] - T_bulk = pybamm.xyz_average(pybamm.size_average(T)) + T_bulk = pybamm.xyzs_average(T) ocp_bulk = self.phase_param.U(sto_bulk, T_bulk) elif self.reaction == "lithium metal plating": T = variables[f"{Domain} electrode temperature [K]"] diff --git a/src/pybamm/models/submodels/interface/sei/base_sei.py b/src/pybamm/models/submodels/interface/sei/base_sei.py index 8a4176d416..d801e25e61 100644 --- a/src/pybamm/models/submodels/interface/sei/base_sei.py +++ b/src/pybamm/models/submodels/interface/sei/base_sei.py @@ -120,6 +120,7 @@ def _get_standard_concentration_variables(self, variables): """Update variables related to the SEI concentration.""" domain, Domain = self.domain_Domain phase_param = self.phase_param + phase_name = phase_param.phase_name reaction_name = self.reaction_name # Set scales to one for the "no SEI" model so that they are not required @@ -152,7 +153,11 @@ def _get_standard_concentration_variables(self, variables): L = variables[f"{Domain} {reaction_name}thickness [m]"] n_SEI = L * L_to_n - n_SEI_xav = pybamm.x_average(n_SEI) + if self.size_distribution: + n_SEI_sav = pybamm.size_average(n_SEI) + n_SEI_xav = pybamm.x_average(n_SEI_sav) + else: + n_SEI_xav = pybamm.x_average(n_SEI) n_SEI_av = pybamm.yz_average(n_SEI_xav) # Calculate change in SEI concentration with respect to initial state @@ -185,11 +190,16 @@ def _get_standard_concentration_variables(self, variables): # Concentration variables are handled slightly differently for SEI on cracks elif self.reaction == "SEI on cracks": L_cr = variables[f"{Domain} {reaction_name}thickness [m]"] - roughness = variables[f"{Domain} electrode roughness ratio"] + roughness = variables[f"{Domain} {phase_name}electrode roughness ratio"] n_SEI_cr = L_cr * L_to_n * (roughness - 1) # SEI on cracks concentration - n_SEI_cr_xav = pybamm.x_average(n_SEI_cr) - n_SEI_cr_av = pybamm.yz_average(n_SEI_cr_xav) + if self.size_distribution: + n_SEI_cr_sav = pybamm.size_average(n_SEI_cr) + n_SEI_cr_xav = pybamm.x_average(n_SEI_cr_sav) + n_SEI_cr_av = pybamm.yz_average(n_SEI_cr_xav) + else: + n_SEI_cr_xav = pybamm.x_average(n_SEI_cr) + n_SEI_cr_av = pybamm.yz_average(n_SEI_cr_xav) # Calculate change in SEI cracks concentration # Initial state depends on roughness (to avoid division by zero) @@ -253,3 +263,36 @@ def _get_standard_reaction_variables(self, j_sei): ) return variables + + def _get_standard_reaction_distribution_variables(self, j_sei): + """ + A private function to obtain the standard variables which + can be derived from the SEI interfacial reaction current + + Parameters + ---------- + j_sei : :class:`pybamm.Symbol` + The SEI interfacial reaction current. + + Returns + ------- + variables : dict + The variables which can be derived from the SEI currents. + """ + domain, Domain = self.domain_Domain + j_sei = pybamm.size_average(j_sei) + variables = { + f"{Domain} electrode {self.reaction_name}" + "interfacial current density [A.m-2]": j_sei, + } + + if self.reaction_loc != "interface": + j_sei_av = pybamm.x_average(j_sei) + variables.update( + { + f"X-averaged {domain} electrode {self.reaction_name}" + "interfacial current density [A.m-2]": j_sei_av, + } + ) + + return variables diff --git a/src/pybamm/models/submodels/interface/sei/no_sei.py b/src/pybamm/models/submodels/interface/sei/no_sei.py index 27cea3b69b..10192a485f 100644 --- a/src/pybamm/models/submodels/interface/sei/no_sei.py +++ b/src/pybamm/models/submodels/interface/sei/no_sei.py @@ -32,12 +32,26 @@ def get_fundamental_variables(self): domain = self.domain.lower() if self.reaction_loc == "interface": zero = pybamm.PrimaryBroadcast(pybamm.Scalar(0), "current collector") + elif self.size_distribution: + zero = pybamm.FullBroadcast( + pybamm.Scalar(0), + f"{domain} {self.phase_name}particle size", + {"secondary": f"{domain} electrode", "tertiary": "current collector"}, + ) else: zero = pybamm.FullBroadcast( pybamm.Scalar(0), f"{domain} electrode", "current collector" ) variables = self._get_standard_thickness_variables(zero) - variables.update(self._get_standard_reaction_variables(zero)) + + if self.size_distribution: + variables.update(self._get_standard_reaction_distribution_variables(zero)) + variables.update( + self._get_standard_size_distribution_interfacial_current_variables(zero) + ) + else: + variables.update(self._get_standard_reaction_variables(zero)) + variables.update(self._get_standard_interfacial_current_variables(zero)) return variables def get_coupled_variables(self, variables): diff --git a/src/pybamm/models/submodels/interface/total_interfacial_current.py b/src/pybamm/models/submodels/interface/total_interfacial_current.py index 79e13e37f6..3df4ff757a 100644 --- a/src/pybamm/models/submodels/interface/total_interfacial_current.py +++ b/src/pybamm/models/submodels/interface/total_interfacial_current.py @@ -165,7 +165,8 @@ def _get_whole_cell_coupled_variables(self, variables): if domain == "separator": var_dict[domain] = zero_s else: - var_dict[domain] = variables[variable_template.format(domain + " ")] + var = variables[variable_template.format(domain + " ")] + var_dict[domain] = var var = pybamm.concatenation(*var_dict.values()) variables.update({variable_template.format(""): var}) diff --git a/src/pybamm/models/submodels/particle/base_particle.py b/src/pybamm/models/submodels/particle/base_particle.py index b774e58a0c..e47039b973 100644 --- a/src/pybamm/models/submodels/particle/base_particle.py +++ b/src/pybamm/models/submodels/particle/base_particle.py @@ -29,7 +29,6 @@ def __init__(self, param, domain, options, phase="primary"): def _get_effective_diffusivity(self, c, T, current): domain, Domain = self.domain_Domain - domain_param = self.domain_param phase_param = self.phase_param domain_options = getattr(self.options, domain) @@ -60,7 +59,7 @@ def _get_effective_diffusivity(self, c, T, current): E = pybamm.r_average(phase_param.E(sto, T)) nu = phase_param.nu theta_M = Omega / (self.param.R * T) * (2 * Omega * E) / (9 * (1 - nu)) - stress_factor = 1 + theta_M * (c - domain_param.c_0) + stress_factor = 1 + theta_M * (c - phase_param.c_0) else: stress_factor = 1 @@ -185,6 +184,7 @@ def _get_distribution_variables(self, R): """ domain, Domain = self.domain_Domain phase_name = self.phase_name + Phase_prefactor = self.phase_param.phase_prefactor R_typ = self.phase_param.R_typ # [m] # Particle-size distribution (area-weighted) @@ -237,21 +237,21 @@ def _get_distribution_variables(self, R): variables = { f"{Domain} {phase_name}particle sizes": R / R_typ, - f"{Domain} {phase_name}particle sizes [m]": R, - f"{Domain} area-weighted {phase_name}particle-size" + f"{Phase_prefactor}{Domain} {phase_name}particle sizes [m]": R, + f"{Phase_prefactor}{Domain} area-weighted {phase_name}particle-size" " distribution [m-1]": f_a_dist, - f"{Domain} volume-weighted {phase_name}particle-size" + f"{Phase_prefactor}{Domain} volume-weighted {phase_name}particle-size" " distribution [m-1]": f_v_dist, - f"{Domain} number-based {phase_name}particle-size" + f"{Phase_prefactor}{Domain} number-based {phase_name}particle-size" " distribution [m-1]": f_num_dist, - f"{Domain} area-weighted mean particle radius [m]": R_a_mean, - f"{Domain} volume-weighted mean particle radius [m]": R_v_mean, - f"{Domain} number-based mean particle radius [m]": R_num_mean, - f"{Domain} area-weighted {phase_name}particle-size" + f"{Phase_prefactor}{Domain} area-weighted mean particle radius [m]": R_a_mean, + f"{Phase_prefactor}{Domain} volume-weighted mean particle radius [m]": R_v_mean, + f"{Phase_prefactor}{Domain} number-based mean particle radius [m]": R_num_mean, + f"{Phase_prefactor}{Domain} area-weighted {phase_name}particle-size" " standard deviation [m]": sd_a, - f"{Domain} volume-weighted {phase_name}particle-size" + f"{Phase_prefactor}{Domain} volume-weighted {phase_name}particle-size" " standard deviation [m]": sd_v, - f"{Domain} number-based {phase_name}particle-size" + f"{Phase_prefactor}{Domain} number-based {phase_name}particle-size" " standard deviation [m]": sd_num, # X-averaged sizes and distributions f"X-averaged {domain} {phase_name}particle sizes [m]": pybamm.x_average(R), diff --git a/src/pybamm/models/submodels/particle/fickian_diffusion.py b/src/pybamm/models/submodels/particle/fickian_diffusion.py index 634e2ce730..5ae5d3a4da 100644 --- a/src/pybamm/models/submodels/particle/fickian_diffusion.py +++ b/src/pybamm/models/submodels/particle/fickian_diffusion.py @@ -31,7 +31,7 @@ def __init__(self, param, domain, options, phase="primary", x_average=False): def get_fundamental_variables(self): domain, Domain = self.domain_Domain phase_name = self.phase_name - + Phase_prefactor = self.phase_param.phase_prefactor variables = {} if self.size_distribution is False: if self.x_average is False: @@ -81,7 +81,7 @@ def get_fundamental_variables(self): ) variables = self._get_distribution_variables(R) f_v_dist = variables[ - f"{Domain} volume-weighted {phase_name}" + f"{Phase_prefactor}{Domain} volume-weighted {phase_name}" "particle-size distribution [m-1]" ] else: @@ -130,6 +130,7 @@ def get_fundamental_variables(self): def get_coupled_variables(self, variables): domain, Domain = self.domain_Domain phase_name = self.phase_name + Phase_prefactor = self.phase_param.phase_prefactor if self.size_distribution is False: if self.x_average is False: @@ -219,7 +220,7 @@ def get_coupled_variables(self, variables): ) variables.update(self._get_standard_flux_distribution_variables(N_s)) # Size-averaged flux variables - R = variables[f"{Domain} {phase_name}particle sizes [m]"] + R = variables[f"{Phase_prefactor}{Domain} {phase_name}particle sizes [m]"] f_a_dist = self.phase_param.f_a_dist(R) D_eff = pybamm.Integral(f_a_dist * D_eff, R) N_s = pybamm.Integral(f_a_dist * N_s, R) diff --git a/src/pybamm/models/submodels/particle_mechanics/base_mechanics.py b/src/pybamm/models/submodels/particle_mechanics/base_mechanics.py index 1301722da0..c41c12fe1b 100644 --- a/src/pybamm/models/submodels/particle_mechanics/base_mechanics.py +++ b/src/pybamm/models/submodels/particle_mechanics/base_mechanics.py @@ -30,8 +30,8 @@ def _get_standard_variables(self, l_cr): domain, Domain = self.domain_Domain l_cr_av = pybamm.x_average(l_cr) variables = { - f"{Domain} particle crack length [m]": l_cr, - f"X-averaged {domain} particle crack length [m]": l_cr_av, + f"{Domain} {self.phase_param.phase_name}particle crack length [m]": l_cr, + f"X-averaged {domain} {self.phase_param.phase_name}particle crack length [m]": l_cr_av, } return variables @@ -61,7 +61,7 @@ def _get_mechanical_results(self, variables): sto = variables[f"{Domain} {phase_name}particle concentration"] Omega = pybamm.r_average(phase_param.Omega(sto, T)) R0 = phase_param.R - c_0 = domain_param.c_0 + c_0 = phase_param.c_0 E0 = pybamm.r_average(phase_param.E(sto, T)) nu = phase_param.nu L0 = domain_param.L @@ -127,21 +127,21 @@ def _get_mechanical_results(self, variables): def _get_standard_surface_variables(self, variables): domain, Domain = self.domain_Domain - phase_name = self.phase_name + phase_name = self.phase_param.phase_name - l_cr = variables[f"{Domain} particle crack length [m]"] + l_cr = variables[f"{Domain} {phase_name}particle crack length [m]"] a = variables[ f"{Domain} electrode {phase_name}surface area to volume ratio [m-1]" ] - rho_cr = self.domain_param.rho_cr - w_cr = self.domain_param.w_cr + rho_cr = self.phase_param.rho_cr + w_cr = self.phase_param.w_cr roughness = 1 + 2 * l_cr * rho_cr * w_cr # ratio of cracks to normal surface a_cr = (roughness - 1) * a # crack surface area to volume ratio roughness_xavg = pybamm.x_average(roughness) variables = { - f"{Domain} crack surface to volume ratio [m-1]": a_cr, - f"{Domain} electrode roughness ratio": roughness, - f"X-averaged {domain} electrode roughness ratio": roughness_xavg, + f"{Domain} {phase_name}crack surface to volume ratio [m-1]": a_cr, + f"{Domain} {phase_name}electrode roughness ratio": roughness, + f"X-averaged {domain} {phase_name}electrode roughness ratio": roughness_xavg, } return variables diff --git a/src/pybamm/models/submodels/particle_mechanics/crack_propagation.py b/src/pybamm/models/submodels/particle_mechanics/crack_propagation.py index b51f1d1ebd..5a84f870d8 100644 --- a/src/pybamm/models/submodels/particle_mechanics/crack_propagation.py +++ b/src/pybamm/models/submodels/particle_mechanics/crack_propagation.py @@ -36,20 +36,20 @@ def __init__(self, param, domain, x_average, options, phase="primary"): def get_fundamental_variables(self): domain, Domain = self.domain_Domain - + phase_name = self.phase_param.phase_name if self.x_average is True: l_cr_av = pybamm.Variable( - f"X-averaged {domain} particle crack length [m]", + f"X-averaged {domain} {phase_name}particle crack length [m]", domain="current collector", - scale=self.domain_param.l_cr_0, + scale=self.phase_param.l_cr_0, ) l_cr = pybamm.PrimaryBroadcast(l_cr_av, f"{domain} electrode") else: l_cr = pybamm.Variable( - f"{Domain} particle crack length [m]", + f"{Domain} {phase_name}particle crack length [m]", domain=f"{domain} electrode", auxiliary_domains={"secondary": "current collector"}, - scale=self.domain_param.l_cr_0, + scale=self.phase_param.l_cr_0, ) variables = self._get_standard_variables(l_cr) @@ -58,22 +58,24 @@ def get_fundamental_variables(self): def get_coupled_variables(self, variables): domain, Domain = self.domain_Domain - + phase_name = self.phase_param.phase_name variables.update(self._get_standard_surface_variables(variables)) variables.update(self._get_mechanical_results(variables)) T = variables[f"{Domain} electrode temperature [K]"] - k_cr = self.domain_param.k_cr(T) - m_cr = self.domain_param.m_cr - b_cr = self.domain_param.b_cr - stress_t_surf = variables[f"{Domain} particle surface tangential stress [Pa]"] - l_cr = variables[f"{Domain} particle crack length [m]"] + k_cr = self.phase_param.k_cr(T) + m_cr = self.phase_param.m_cr + b_cr = self.phase_param.b_cr + stress_t_surf = variables[ + f"{Domain} {phase_name}particle surface tangential stress [Pa]" + ] + l_cr = variables[f"{Domain} {phase_name}particle crack length [m]"] # # compressive stress will not lead to crack propagation dK_SIF = stress_t_surf * b_cr * pybamm.sqrt(np.pi * l_cr) * (stress_t_surf >= 0) dl_cr = k_cr * (dK_SIF**m_cr) / 3600 # divide by 3600 to replace t0_cr variables.update( { - f"{Domain} particle cracking rate [m.s-1]": dl_cr, - f"X-averaged {domain} particle cracking rate [m.s-1]": pybamm.x_average( + f"{Domain} {phase_name}particle cracking rate [m.s-1]": dl_cr, + f"X-averaged {domain} {phase_name}particle cracking rate [m.s-1]": pybamm.x_average( dl_cr ), } @@ -82,36 +84,44 @@ def get_coupled_variables(self, variables): def set_rhs(self, variables): domain, Domain = self.domain_Domain - + phase_name = self.phase_param.phase_name if self.x_average is True: - l_cr = variables[f"X-averaged {domain} particle crack length [m]"] - dl_cr = variables[f"X-averaged {domain} particle cracking rate [m.s-1]"] + l_cr = variables[ + f"X-averaged {domain} {phase_name}particle crack length [m]" + ] + dl_cr = variables[ + f"X-averaged {domain} {phase_name}particle cracking rate [m.s-1]" + ] else: - l_cr = variables[f"{Domain} particle crack length [m]"] - dl_cr = variables[f"{Domain} particle cracking rate [m.s-1]"] + l_cr = variables[f"{Domain} {phase_name}particle crack length [m]"] + dl_cr = variables[f"{Domain} {phase_name}particle cracking rate [m.s-1]"] self.rhs = {l_cr: dl_cr} def set_initial_conditions(self, variables): domain, Domain = self.domain_Domain - - l_cr_0 = self.domain_param.l_cr_0 + phase_name = self.phase_param.phase_name + l_cr_0 = self.phase_param.l_cr_0 if self.x_average is True: - l_cr = variables[f"X-averaged {domain} particle crack length [m]"] + l_cr = variables[ + f"X-averaged {domain} {phase_name}particle crack length [m]" + ] else: - l_cr = variables[f"{Domain} particle crack length [m]"] + l_cr = variables[f"{Domain} {phase_name}particle crack length [m]"] l_cr_0 = pybamm.PrimaryBroadcast(l_cr_0, f"{domain} electrode") self.initial_conditions = {l_cr: l_cr_0} def add_events_from(self, variables): domain, Domain = self.domain_Domain - + phase_name = self.phase_param.phase_name if self.x_average is True: - l_cr = variables[f"X-averaged {domain} particle crack length [m]"] + l_cr = variables[ + f"X-averaged {domain} {phase_name}particle crack length [m]" + ] else: - l_cr = variables[f"{Domain} particle crack length [m]"] + l_cr = variables[f"{Domain} {phase_name}particle crack length [m]"] self.events.append( pybamm.Event( - f"{domain} particle crack length larger than particle radius", - 1 - pybamm.max(l_cr) / self.domain_param.prim.R_typ, + f"{domain} {phase_name} particle crack length larger than particle radius", + 1 - pybamm.max(l_cr) / self.phase_param.R_typ, ) ) diff --git a/src/pybamm/parameters/geometric_parameters.py b/src/pybamm/parameters/geometric_parameters.py index 0155bbe03c..2998e97318 100644 --- a/src/pybamm/parameters/geometric_parameters.py +++ b/src/pybamm/parameters/geometric_parameters.py @@ -165,6 +165,14 @@ def f_a_dist(self, R): Dimensional electrode area-weighted particle-size distribution """ Domain = self.domain.capitalize() + if len(R.domain) == 0: + self.phase_prefactor = "" + elif "primary" in R.domains["primary"][0]: + self.phase_prefactor = "Primary: " + elif "secondary" in R.domains["primary"][0]: + self.phase_prefactor = "Secondary: " + else: + self.phase_prefactor = "" inputs = {f"{self.phase_prefactor}{Domain} particle-size variable [m]": R} return pybamm.FunctionParameter( f"{self.phase_prefactor}{Domain} " diff --git a/src/pybamm/parameters/lead_acid_parameters.py b/src/pybamm/parameters/lead_acid_parameters.py index 7da6a19410..7beef48142 100644 --- a/src/pybamm/parameters/lead_acid_parameters.py +++ b/src/pybamm/parameters/lead_acid_parameters.py @@ -344,6 +344,7 @@ def __init__(self, phase, domain_param): self.domain = domain_param.domain self.main_param = domain_param.main_param self.geo = domain_param.geo.prim + self.phase_name = phase def _set_parameters(self): main = self.main_param diff --git a/src/pybamm/parameters/lithium_ion_parameters.py b/src/pybamm/parameters/lithium_ion_parameters.py index 1ea20aeb3e..1ccd21975a 100644 --- a/src/pybamm/parameters/lithium_ion_parameters.py +++ b/src/pybamm/parameters/lithium_ion_parameters.py @@ -269,20 +269,6 @@ def _set_parameters(self): self.b_s = self.geo.b_s self.tau_s = self.geo.tau_s - # Mechanical parameters - self.c_0 = pybamm.Parameter( - f"{Domain} electrode reference concentration for free of deformation " - "[mol.m-3]" - ) - - 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( - f"{Domain} electrode number of cracks per unit area [m-2]" - ) - self.b_cr = pybamm.Parameter(f"{Domain} electrode Paris' law constant b") - self.m_cr = pybamm.Parameter(f"{Domain} electrode Paris' law constant m") - # Utilisation parameters self.u_init = pybamm.Parameter( f"Initial {domain} electrode interface utilisation" @@ -307,15 +293,6 @@ def sigma(self, T): f"{Domain} electrode conductivity [S.m-1]", inputs ) - def k_cr(self, T): - """ - Cracking rate for the electrode; - """ - Domain = self.domain.capitalize() - return pybamm.FunctionParameter( - f"{Domain} electrode cracking rate", {"Temperature [K]": T} - ) - def LAM_rate_current(self, i, T): """ Dimensional rate of loss of active material as a function of applied current @@ -516,6 +493,34 @@ def _set_parameters(self): f"{pref}{Domain} electrode reaction-driven LAM factor [m3.mol-1]" ) + # Mechanical parameters + self.c_0 = pybamm.Parameter( + f"{pref}{Domain} electrode reference concentration for free of deformation " + "[mol.m-3]" + ) + + self.l_cr_0 = pybamm.Parameter( + f"{pref}{Domain} electrode initial crack length [m]" + ) + self.w_cr = pybamm.Parameter( + f"{pref}{Domain} electrode initial crack width [m]" + ) + self.rho_cr = pybamm.Parameter( + f"{pref}{Domain} electrode number of cracks per unit area [m-2]" + ) + self.b_cr = pybamm.Parameter(f"{pref}{Domain} electrode Paris' law constant b") + self.m_cr = pybamm.Parameter(f"{pref}{Domain} electrode Paris' law constant m") + + def k_cr(self, T): + """ + Cracking rate for the electrode; + """ + phase_prefactor = self.phase_prefactor + Domain = self.domain.capitalize() + return pybamm.FunctionParameter( + f"{phase_prefactor}{Domain} electrode cracking rate", {"Temperature [K]": T} + ) + def D(self, c_s, T, lithiation=None): """ Dimensional diffusivity in particle. In the parameter sets this is defined as diff --git a/src/pybamm/parameters/size_distribution_parameters.py b/src/pybamm/parameters/size_distribution_parameters.py index 60383fff50..6b7087b8c8 100644 --- a/src/pybamm/parameters/size_distribution_parameters.py +++ b/src/pybamm/parameters/size_distribution_parameters.py @@ -9,14 +9,31 @@ def get_size_distribution_parameters( param, R_n_av=None, + R_n_av_prim=None, + R_n_av_sec=None, R_p_av=None, + R_p_av_prim=None, + R_p_av_sec=None, sd_n=0.3, + sd_n_prim=0.3, + sd_n_sec=0.3, sd_p=0.3, + sd_p_prim=0.3, + sd_p_sec=0.3, R_min_n=None, + R_min_n_prim=None, + R_min_n_sec=None, R_min_p=None, + R_min_p_prim=None, + R_min_p_sec=None, R_max_n=None, + R_max_n_prim=None, + R_max_n_sec=None, R_max_p=None, + R_max_p_prim=None, + R_max_p_sec=None, working_electrode="both", + composite="neither", ): """ A convenience method to add standard area-weighted particle-size distribution @@ -62,54 +79,132 @@ def get_size_distribution_parameters( """ if working_electrode == "both": # Radii from given parameter set - R_n_typ = param["Negative particle radius [m]"] - + if composite == "negative" or composite == "both": + R_n_typ_prim = param["Primary: Negative particle radius [m]"] + R_n_typ_sec = param["Secondary: Negative particle radius [m]"] + # Set the mean particle radii + R_n_av_prim = R_n_av_prim or R_n_typ_prim + R_n_av_sec = R_n_av_sec or R_n_typ_sec + # Minimum radii + R_min_n_prim = R_min_n_prim or np.max([0, 1 - sd_n_prim * 5]) + R_min_n_sec = R_min_n_sec or np.max([0, 1 - sd_n_sec * 5]) + + # Max radii + R_max_n_prim = R_max_n_prim or (1 + sd_n_prim * 5) + R_max_n_sec = R_max_n_sec or (1 + sd_n_sec * 5) + + # Area-weighted particle-size distribution + def f_a_dist_n_prim(R): + return lognormal(R, R_n_av_prim, sd_n_prim * R_n_av_prim) + + def f_a_dist_n_sec(R): + return lognormal(R, R_n_av_sec, sd_n_sec * R_n_av_sec) + + param.update( + { + "Primary: Negative minimum particle radius [m]": R_min_n_prim + * R_n_av_prim, + "Primary: Negative maximum particle radius [m]": R_max_n_prim + * R_n_av_prim, + "Primary: Negative area-weighted " + + "particle-size distribution [m-1]": f_a_dist_n_prim, + "Secondary: Negative minimum particle radius [m]": R_min_n_sec + * R_n_av_sec, + "Secondary: Negative maximum particle radius [m]": R_max_n_sec + * R_n_av_sec, + "Secondary: Negative area-weighted " + + "particle-size distribution [m-1]": f_a_dist_n_sec, + }, + check_already_exists=False, + ) + else: + R_n_typ = param["Negative particle radius [m]"] + # Set the mean particle radii + R_n_av = R_n_av or R_n_typ + # Minimum radii + R_min_n = R_min_n or np.max([0, 1 - sd_n * 5]) + + # Max radii + R_max_n = R_max_n or (1 + sd_n * 5) + + # Area-weighted particle-size distribution + def f_a_dist_n(R): + return lognormal(R, R_n_av, sd_n * R_n_av) + + param.update( + { + "Negative minimum particle radius [m]": R_min_n * R_n_av, + "Negative maximum particle radius [m]": R_max_n * R_n_av, + "Negative area-weighted " + + "particle-size distribution [m-1]": f_a_dist_n, + }, + check_already_exists=False, + ) + + if composite == "positive" or composite == "both": + R_p_typ_prim = param["Primary: Positive particle radius [m]"] + R_p_typ_sec = param["Secondary: Positive particle radius [m]"] # Set the mean particle radii - R_n_av = R_n_av or R_n_typ - + R_p_av_prim = R_p_av_prim or R_p_typ_prim + R_p_av_sec = R_p_av_sec or R_p_typ_sec # Minimum radii - R_min_n = R_min_n or np.max([0, 1 - sd_n * 5]) + R_min_p_prim = R_min_p_prim or np.max([0, 1 - sd_p_prim * 5]) + R_min_p_sec = R_min_p_sec or np.max([0, 1 - sd_p_sec * 5]) # Max radii - R_max_n = R_max_n or (1 + sd_n * 5) + R_max_p_prim = R_max_p_prim or (1 + sd_p_prim * 5) + R_max_p_sec = R_max_p_sec or (1 + sd_p_sec * 5) # Area-weighted particle-size distribution - def f_a_dist_n(R): - return lognormal(R, R_n_av, sd_n * R_n_av) + def f_a_dist_p_prim(R): + return lognormal(R, R_p_av_prim, sd_p_prim * R_p_av_prim) + + def f_a_dist_p_sec(R): + return lognormal(R, R_p_av_sec, sd_p_sec * R_p_av_sec) param.update( { - "Negative minimum particle radius [m]": R_min_n * R_n_av, - "Negative maximum particle radius [m]": R_max_n * R_n_av, - "Negative area-weighted " - + "particle-size distribution [m-1]": f_a_dist_n, + "Primary: Positive minimum particle radius [m]": R_min_p_prim + * R_p_av_prim, + "Primary: Positive maximum particle radius [m]": R_max_p_prim + * R_p_av_prim, + "Primary: Positive area-weighted " + + "particle-size distribution [m-1]": f_a_dist_p_prim, + "Secondary: Positive minimum particle radius [m]": R_min_p_sec + * R_p_av_sec, + "Secondary: Positive maximum particle radius [m]": R_max_p_sec + * R_p_av_sec, + "Secondary: Positive area-weighted " + + "particle-size distribution [m-1]": f_a_dist_p_sec, }, check_already_exists=False, ) - # Radii from given parameter set - R_p_typ = param["Positive particle radius [m]"] + else: + # Radii from given parameter set + R_p_typ = param["Positive particle radius [m]"] - # Set the mean particle radii - R_p_av = R_p_av or R_p_typ + # Set the mean particle radii + R_p_av = R_p_av or R_p_typ - # Minimum radii - R_min_p = R_min_p or np.max([0, 1 - sd_p * 5]) + # Minimum radii + R_min_p = R_min_p or np.max([0, 1 - sd_p * 5]) - # Max radii - R_max_p = R_max_p or (1 + sd_p * 5) + # Max radii + R_max_p = R_max_p or (1 + sd_p * 5) - # Area-weighted particle-size distribution - def f_a_dist_p(R): - return lognormal(R, R_p_av, sd_p * R_p_av) + # Area-weighted particle-size distribution + def f_a_dist_p(R): + return lognormal(R, R_p_av, sd_p * R_p_av) - param.update( - { - "Positive minimum particle radius [m]": R_min_p * R_p_av, - "Positive maximum particle radius [m]": R_max_p * R_p_av, - "Positive area-weighted " + "particle-size distribution [m-1]": f_a_dist_p, - }, - check_already_exists=False, - ) + param.update( + { + "Positive minimum particle radius [m]": R_min_p * R_p_av, + "Positive maximum particle radius [m]": R_max_p * R_p_av, + "Positive area-weighted " + + "particle-size distribution [m-1]": f_a_dist_p, + }, + check_already_exists=False, + ) return param diff --git a/src/pybamm/plotting/plot_summary_variables.py b/src/pybamm/plotting/plot_summary_variables.py index bd4db0ee6c..7181563fc3 100644 --- a/src/pybamm/plotting/plot_summary_variables.py +++ b/src/pybamm/plotting/plot_summary_variables.py @@ -62,13 +62,13 @@ def plot_summary_variables( for solution in solutions: # plot summary variable v/s cycle number ax.plot( - solution.summary_variables["Cycle number"], + solution.summary_variables.cycle_number, solution.summary_variables[var], ) # label the axes ax.set_xlabel("Cycle number") ax.set_ylabel(var) - ax.set_xlim([1, solution.summary_variables["Cycle number"][-1]]) + ax.set_xlim([1, solution.summary_variables.cycle_number[-1]]) fig.tight_layout() diff --git a/src/pybamm/simulation.py b/src/pybamm/simulation.py index 9dafb23ce2..d499f77403 100644 --- a/src/pybamm/simulation.py +++ b/src/pybamm/simulation.py @@ -920,7 +920,7 @@ def solve( # See PR #3995 if voltage_stop is not None: min_voltage = np.min(cycle_solution["Battery voltage [V]"].data) - logs["summary variables"]["Minimum voltage [V]"] = min_voltage + logs["Minimum voltage [V]"] = min_voltage callbacks.on_cycle_end(logs) @@ -941,7 +941,7 @@ def solve( if self._solution is not None and len(all_cycle_solutions) > 0: self._solution.cycles = all_cycle_solutions - self._solution.set_summary_variables(all_summary_variables) + self._solution.update_summary_variables(all_summary_variables) self._solution.all_first_states = all_first_states callbacks.on_experiment_end(logs) diff --git a/src/pybamm/solvers/solution.py b/src/pybamm/solvers/solution.py index 256d596fd4..49e8209830 100644 --- a/src/pybamm/solvers/solution.py +++ b/src/pybamm/solvers/solution.py @@ -563,16 +563,10 @@ def initial_start_time(self, value): """Updates the initial start time of the experiment""" self._initial_start_time = value - def set_summary_variables(self, all_summary_variables): - summary_variables = {var: [] for var in all_summary_variables[0]} - for sum_vars in all_summary_variables: - for name, value in sum_vars.items(): - summary_variables[name].append(value) - - summary_variables["Cycle number"] = range(1, len(all_summary_variables) + 1) + def update_summary_variables(self, all_summary_variables): self.all_summary_variables = all_summary_variables - self._summary_variables = pybamm.FuzzyDict( - {name: np.array(value) for name, value in summary_variables.items()} + self._summary_variables = pybamm.SummaryVariables( + self, cycle_summary_variables=all_summary_variables ) def update(self, variables): @@ -1142,8 +1136,8 @@ def make_cycle_solution( cycle_solution.steps = step_solutions - cycle_summary_variables = _get_cycle_summary_variables( - cycle_solution, esoh_solver, user_inputs=inputs + cycle_summary_variables = pybamm.SummaryVariables( + cycle_solution, esoh_solver=esoh_solver, user_inputs=inputs ) cycle_first_state = cycle_solution.first_state @@ -1154,46 +1148,3 @@ def make_cycle_solution( cycle_solution = None return cycle_solution, cycle_summary_variables, cycle_first_state - - -def _get_cycle_summary_variables(cycle_solution, esoh_solver, user_inputs=None): - user_inputs = user_inputs or {} - model = cycle_solution.all_models[0] - cycle_summary_variables = pybamm.FuzzyDict({}) - - # Summary variables - summary_variables = model.summary_variables - first_state = cycle_solution.first_state - last_state = cycle_solution.last_state - for var in summary_variables: - data_first = first_state[var].data - data_last = last_state[var].data - cycle_summary_variables[var] = data_last[0] - var_lowercase = var[0].lower() + var[1:] - cycle_summary_variables["Change in " + var_lowercase] = ( - data_last[0] - data_first[0] - ) - - # eSOH variables (full-cell lithium-ion model only, for now) - if ( - esoh_solver is not None - and isinstance(model, pybamm.lithium_ion.BaseModel) - and model.options.electrode_types["negative"] == "porous" - and "Negative electrode capacity [A.h]" in model.variables - and "Positive electrode capacity [A.h]" in model.variables - ): - Q_n = last_state["Negative electrode capacity [A.h]"].data[0] - Q_p = last_state["Positive electrode capacity [A.h]"].data[0] - Q_Li = last_state["Total lithium capacity in particles [A.h]"].data[0] - all_inputs = {**user_inputs, "Q_n": Q_n, "Q_p": Q_p, "Q_Li": Q_Li} - try: - esoh_sol = esoh_solver.solve(inputs=all_inputs) - except pybamm.SolverError as error: # pragma: no cover - raise pybamm.SolverError( - "Could not solve for summary variables, run " - "`sim.solve(calc_esoh=False)` to skip this step" - ) from error - - cycle_summary_variables.update(esoh_sol) - - return cycle_summary_variables diff --git a/src/pybamm/solvers/summary_variable.py b/src/pybamm/solvers/summary_variable.py new file mode 100644 index 0000000000..812dea72dd --- /dev/null +++ b/src/pybamm/solvers/summary_variable.py @@ -0,0 +1,188 @@ +# +# Summary Variable class +# +from __future__ import annotations +import pybamm +import numpy as np +from typing import Any + + +class SummaryVariables: + """ + Class for managing and calculating summary variables from a PyBaMM solution. + Summary variables are only calculated when simulations are run with PyBaMM + Experiments. + + Parameters + ---------- + solution : :class:`pybamm.Solution` + The solution object to be used for creating the processed variables. + cycle_summary_variables : list[pybamm.SummaryVariables], optional + A list of cycle summary variables. + esoh_solver : :class:`pybamm.lithium_ion.ElectrodeSOHSolver`, optional + Solver for electrode state-of-health (eSOH) calculations. + user_inputs : dict, optional + Additional user inputs for calculations. + + Attributes + ---------- + cycle_number : array[int] + Stores the cycle number for each saved cycle, for use when plotting. + Length is equal to the number of cycles in a solution. + """ + + def __init__( + self, + solution: pybamm.Solution, + cycle_summary_variables: list[SummaryVariables] | None = None, + esoh_solver: pybamm.lithium_ion.ElectrodeSOHSolver | None = None, + user_inputs: dict[str, Any] | None = None, + ): + self.user_inputs = user_inputs or {} + self.esoh_solver = esoh_solver + self._variables = {} # Store computed variables + self.cycle_number = None + + model = solution.all_models[0] + self._possible_variables = model.summary_variables # minus esoh variables + self._esoh_variables = None # Store eSOH variable names + + # Flag if eSOH calculations are needed + self.calc_esoh = ( + self.esoh_solver is not None + and isinstance(model, pybamm.lithium_ion.BaseModel) + and model.options.electrode_types["negative"] == "porous" + and "Negative electrode capacity [A.h]" in model.variables + and "Positive electrode capacity [A.h]" in model.variables + ) + + # Initialize based on cycle information + if cycle_summary_variables: + self._initialize_for_cycles(cycle_summary_variables) + else: + self.first_state = solution.first_state + self.last_state = solution.last_state + self.cycles = None + + def _initialize_for_cycles(self, cycle_summary_variables: list[SummaryVariables]): + """Initialize attributes for when multiple cycles are provided.""" + self.first_state = None + self.last_state = None + self.cycles = cycle_summary_variables + self.cycle_number = np.arange(1, len(self.cycles) + 1) + first_cycle = self.cycles[0] + self.calc_esoh = first_cycle.calc_esoh + self.esoh_solver = first_cycle.esoh_solver + self.user_inputs = first_cycle.user_inputs + + @property + def all_variables(self) -> list[str]: + """ + Return names of all possible summary variables, including eSOH variables + if appropriate. + """ + try: + return self._all_variables + except AttributeError: + base_vars = self._possible_variables.copy() + base_vars.extend( + f"Change in {var[0].lower() + var[1:]}" + for var in self._possible_variables + ) + + if self.calc_esoh: + base_vars.extend(self.esoh_variables) + + self._all_variables = base_vars + return self._all_variables + + @property + def esoh_variables(self) -> list[str] | None: + """Return names of all eSOH variables.""" + if self.calc_esoh and self._esoh_variables is None: + esoh_model = self.esoh_solver._get_electrode_soh_sims_full().model + esoh_vars = list(esoh_model.variables.keys()) + self._esoh_variables = esoh_vars + return self._esoh_variables + + def __getitem__(self, key: str) -> float | list[float]: + """ + Access or compute a summary variable by its name. + + Parameters + ---------- + key : str + The name of the variable + + Returns + ------- + float or list[float] + """ + + if key in self._variables: + # return it if it exists + return self._variables[key] + elif key not in self.all_variables: + # check it's listed as a summary variable + raise KeyError(f"Variable '{key}' is not a summary variable.") + else: + # otherwise create it, save it and then return it + if self.calc_esoh and key in self._esoh_variables: + self.update_esoh() + else: + base_key = key.removeprefix("Change in ") + base_key = base_key[0].upper() + base_key[1:] + # this will create 'X' and 'Change in x' at the same time + self.update(base_key) + return self._variables[key] + + def update(self, var: str): + """Compute and store a variable and its change.""" + var_lowercase = var[0].lower() + var[1:] + if self.cycles: + self._update_multiple_cycles(var, var_lowercase) + else: + self._update(var, var_lowercase) + + def _update_multiple_cycles(self, var: str, var_lowercase: str): + """Creates aggregated summary variables for where more than one cycle exists.""" + var_cycle = [cycle[var] for cycle in self.cycles] + change_var_cycle = [ + cycle[f"Change in {var_lowercase}"] for cycle in self.cycles + ] + self._variables[var] = var_cycle + self._variables[f"Change in {var_lowercase}"] = change_var_cycle + + def _update(self, var: str, var_lowercase: str): + """Create variable `var` for a single cycle.""" + data_first = self.first_state[var].data + data_last = self.last_state[var].data + self._variables[var] = data_last[0] + self._variables[f"Change in {var_lowercase}"] = data_last[0] - data_first[0] + + def update_esoh(self): + """Create all aggregated eSOH variables""" + if self.cycles is not None: + var_cycle = [cycle._get_esoh_variables() for cycle in self.cycles] + aggregated_vars = {k: [] for k in var_cycle[0].keys()} + for cycle in var_cycle: + for k, v in cycle.items(): + aggregated_vars[k].append(v) + self._variables.update(aggregated_vars) + else: + self._variables.update(self._get_esoh_variables()) + + def _get_esoh_variables(self) -> dict[str, float]: + """Compute eSOH variables for a single solution.""" + Q_n = self.last_state["Negative electrode capacity [A.h]"].data[0] + Q_p = self.last_state["Positive electrode capacity [A.h]"].data[0] + Q_Li = self.last_state["Total lithium capacity in particles [A.h]"].data[0] + all_inputs = {**self.user_inputs, "Q_n": Q_n, "Q_p": Q_p, "Q_Li": Q_Li} + try: + esoh_sol = self.esoh_solver.solve(inputs=all_inputs) + except pybamm.SolverError as error: # pragma: no cover + raise pybamm.SolverError( + "Could not solve for eSOH summary variables" + ) from error + + return esoh_sol diff --git a/src/pybamm/util.py b/src/pybamm/util.py index 5b10b23fcb..39fcdde6f1 100644 --- a/src/pybamm/util.py +++ b/src/pybamm/util.py @@ -94,10 +94,16 @@ def __getitem__(self, key): ) from error best_matches = self.get_best_matches(key) for k in best_matches: - if key in k and k.endswith("]"): + if key in k and k.endswith("]") and not key.endswith("]"): raise KeyError( f"'{key}' not found. Use the dimensional version '{k}' instead." ) from error + elif key in k and ( + k.startswith("Primary") or k.startswith("Secondary") + ): + raise KeyError( + f"'{key}' not found. If you are using a composite model, you may need to use {k} instead. Otherwise, best matches are {best_matches}" + ) from error raise KeyError( f"'{key}' not found. Best matches are {best_matches}" ) from error diff --git a/tests/integration/test_models/standard_output_tests.py b/tests/integration/test_models/standard_output_tests.py index ca5e27607d..4e6788724f 100644 --- a/tests/integration/test_models/standard_output_tests.py +++ b/tests/integration/test_models/standard_output_tests.py @@ -322,17 +322,17 @@ def __init__(self, model, param, disc, solution, operating_condition): # Take only the x-averaged of these for now, since variables cannot have # 4 domains yet self.c_s_n_dist = solution[ - "X-averaged negative particle concentration distribution" + f"X-averaged negative {self.phase_name_n}particle concentration distribution" ] self.c_s_p_dist = solution[ - "X-averaged positive particle concentration distribution" + f"X-averaged positive {self.phase_name_p}particle concentration distribution" ] self.c_s_n_surf_dist = solution[ - "Negative particle surface concentration distribution" + f"Negative {self.phase_name_n}particle surface concentration distribution" ] self.c_s_p_surf_dist = solution[ - "Positive particle surface concentration distribution" + f"Positive {self.phase_name_p}particle surface concentration distribution" ] def test_concentration_increase_decrease(self): diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py index 768872e454..2456a89b60 100644 --- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py +++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py @@ -462,7 +462,8 @@ def lico2_volume_change_Ai2020(sto): "Secondary: Negative electrode partial molar volume [m3.mol-1]": 3.1e-06, "Secondary: Negative electrode Young's modulus [Pa]": 15000000000.0, "Secondary: Negative electrode Poisson's ratio": 0.3, - "Negative electrode reference concentration for free of deformation [mol.m-3]": 0.0, + "Primary: Negative electrode reference concentration for free of deformation [mol.m-3]": 0.0, + "Secondary: Negative electrode reference concentration for free of deformation [mol.m-3]": 0.0, "Primary: Negative electrode volume change": graphite_volume_change_Ai2020, "Secondary: Negative electrode volume change": graphite_volume_change_Ai2020, "Positive electrode partial molar volume [m3.mol-1]": -7.28e-07, diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py index 5ede55ff6d..7a728a744b 100644 --- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py +++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py @@ -52,6 +52,10 @@ def setup(self): "R_p": 3, "y": 5, "z": 5, + "R_n_prim": 3, + "R_n_sec": 3, + "R_p_prim": 3, + "R_p_sec": 3, } def test_basic_processing(self): @@ -62,6 +66,31 @@ def test_basic_processing(self): ) modeltest.test_all() + def test_composite(self): + options = { + "particle phases": ("2", "1"), + "open-circuit potential": (("single", "current sigmoid"), "single"), + "particle size": "distribution", + } + parameter_values = pybamm.ParameterValues("Chen2020_composite") + name = "Negative electrode active material volume fraction" + x = 0.1 + parameter_values.update( + {f"Primary: {name}": (1 - x) * 0.75, f"Secondary: {name}": x * 0.75} + ) + parameter_values = pybamm.get_size_distribution_parameters( + parameter_values, + composite="negative", + R_min_n_prim=0.9, + R_min_n_sec=0.9, + R_max_n_prim=1.1, + R_max_n_sec=1.1, + ) + # self.run_basic_processing_test(options, parameter_values=parameter_values) + model = pybamm.lithium_ion.DFN(options) + modeltest = tests.StandardModelTest(model, parameter_values=parameter_values) + modeltest.test_all() + def test_basic_processing_tuple(self): options = {"particle size": ("single", "distribution")} model = pybamm.lithium_ion.DFN(options) diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_newman_tobias.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_newman_tobias.py index fb091ef56d..c24a063f70 100644 --- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_newman_tobias.py +++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_newman_tobias.py @@ -31,6 +31,9 @@ def test_composite_graphite_silicon_sei(self): def test_composite_reaction_driven_LAM(self): pass # skip this test + def test_particle_size_composite(self): + pass # skip this test + def test_composite_stress_driven_LAM(self): pass # skip this test diff --git a/tests/integration/test_solvers/test_idaklu.py b/tests/integration/test_solvers/test_idaklu.py index 3ee96a9ccb..1766c1f2aa 100644 --- a/tests/integration/test_solvers/test_idaklu.py +++ b/tests/integration/test_solvers/test_idaklu.py @@ -202,5 +202,5 @@ def test_with_experiments(self): ) # check summary variables are the same if output variables are specified - for var in summary_vars[0].keys(): + for var in model.summary_variables: assert summary_vars[0][var] == summary_vars[1][var] diff --git a/tests/unit/test_citations.py b/tests/unit/test_citations.py index 9b1c9fa2d3..dccef9b51c 100644 --- a/tests/unit/test_citations.py +++ b/tests/unit/test_citations.py @@ -357,13 +357,13 @@ def test_sripad_2020(self): citations._reset() assert "Sripad2020" not in citations._papers_to_cite - pybamm.kinetics.Marcus(None, None, None, None, None) + pybamm.kinetics.Marcus(None, "negative", None, None, None) assert "Sripad2020" in citations._papers_to_cite assert "Sripad2020" in citations._citation_tags.keys() citations._reset() assert "Sripad2020" not in citations._papers_to_cite - pybamm.kinetics.MarcusHushChidsey(None, None, None, None, None) + pybamm.kinetics.MarcusHushChidsey(None, "negative", None, None, None) assert "Sripad2020" in citations._papers_to_cite assert "Sripad2020" in citations._citation_tags.keys() diff --git a/tests/unit/test_expression_tree/test_averages.py b/tests/unit/test_expression_tree/test_averages.py index 5f9736693e..9065b85ca9 100644 --- a/tests/unit/test_expression_tree/test_averages.py +++ b/tests/unit/test_expression_tree/test_averages.py @@ -201,7 +201,14 @@ def test_size_average(self): ) assert average_a == a - for domain in [["negative particle size"], ["positive particle size"]]: + for domain in [ + ["negative particle size"], + ["positive particle size"], + ["negative primary particle size"], + ["positive primary particle size"], + ["negative secondary particle size"], + ["positive secondary particle size"], + ]: a = pybamm.Symbol("a", domain=domain) R = pybamm.SpatialVariable("R", domain) av_a = pybamm.size_average(a) diff --git a/tests/unit/test_expression_tree/test_symbol.py b/tests/unit/test_expression_tree/test_symbol.py index d7374e8da7..735724ef11 100644 --- a/tests/unit/test_expression_tree/test_symbol.py +++ b/tests/unit/test_expression_tree/test_symbol.py @@ -227,8 +227,27 @@ def test_multiple_symbols(self): "c", "a", ] + expected_postorder = [ + "a", + "c", + "*", + "a", + "b", + "*", + "c", + "*", + "a", + "+", + "c", + "a", + "*", + "-", + "*", + ] for node, expect in zip(exp.pre_order(), expected_preorder): assert node.name == expect + for node, expect in zip(exp.post_order(), expected_postorder): + assert node.name == expect def test_symbol_diff(self): a = pybamm.Symbol("a") diff --git a/tests/unit/test_geometry/test_battery_geometry.py b/tests/unit/test_geometry/test_battery_geometry.py index c73ef89af1..30d417fcbd 100644 --- a/tests/unit/test_geometry/test_battery_geometry.py +++ b/tests/unit/test_geometry/test_battery_geometry.py @@ -77,6 +77,16 @@ def test_geometry(self): geometry["positive secondary particle"]["r_p_sec"]["max"] == geo.p.sec.R_typ ) + options = { + "particle phases": "2", + "particle size": "distribution", + } + geometry = pybamm.battery_geometry(options=options) + assert "negative primary particle size" in geometry + assert "positive primary particle size" in geometry + assert "negative secondary particle size" in geometry + assert "positive secondary particle size" in geometry + def test_geometry_error(self): with pytest.raises(pybamm.GeometryError, match="Invalid current"): pybamm.battery_geometry( diff --git a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py index be7ee861c4..554b65d880 100644 --- a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py +++ b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py @@ -140,6 +140,10 @@ def test_default_var_pts(self): "z": 10, "R_n": 30, "R_p": 30, + "R_n_prim": 30, + "R_p_prim": 30, + "R_n_sec": 30, + "R_p_sec": 30, } model = pybamm.BaseBatteryModel({"dimensionality": 0}) assert var_pts == model.default_var_pts diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py index 973a0f348b..56f4c76f44 100644 --- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py +++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py @@ -27,6 +27,10 @@ def test_well_posed_size_distribution_tuple(self): options = {"particle size": ("single", "distribution")} self.check_well_posedness(options) + def test_well_posed_size_distribution_composite(self): + options = {"particle size": "distribution", "particle phases": "2"} + self.check_well_posedness(options) + def test_well_posed_current_sigmoid_ocp_with_psd(self): options = { "open-circuit potential": "current sigmoid", diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py index ad02212840..5f693264d3 100644 --- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py +++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py @@ -123,6 +123,18 @@ def test_wycisk_ocp(self): model = pybamm.lithium_ion.MPM(options) model.check_well_posedness() + def test_mpm_with_lithium_plating(self): + options = { + "lithium plating": "irreversible", + } + model = pybamm.lithium_ion.MPM(options) + model.check_well_posedness() + options = { + "lithium plating": "reversible", + } + model = pybamm.lithium_ion.MPM(options) + model.check_well_posedness() + class TestMPMExternalCircuits: def test_well_posed_voltage(self): @@ -196,15 +208,3 @@ def test_well_posed_both_swelling_only_not_implemented(self): options = {"particle mechanics": "swelling only"} with pytest.raises(NotImplementedError): pybamm.lithium_ion.MPM(options) - - -class TestMPMWithPlating: - def test_well_posed_reversible_plating_not_implemented(self): - options = {"lithium plating": "reversible"} - with pytest.raises(NotImplementedError): - pybamm.lithium_ion.MPM(options) - - def test_well_posed_irreversible_plating_not_implemented(self): - options = {"lithium plating": "irreversible"} - with pytest.raises(NotImplementedError): - pybamm.lithium_ion.MPM(options) diff --git a/tests/unit/test_parameters/test_parameter_values.py b/tests/unit/test_parameters/test_parameter_values.py index cdcfb30ede..0e32f60476 100644 --- a/tests/unit/test_parameters/test_parameter_values.py +++ b/tests/unit/test_parameters/test_parameter_values.py @@ -893,7 +893,7 @@ def dist(R): var_av_proc = param.process_symbol(var_av) assert isinstance(var_av_proc, pybamm.SizeAverage) - R = pybamm.SpatialVariable("R", "negative particle size") + R = pybamm.SpatialVariable("R_n", "negative particle size") assert var_av_proc.f_a_dist == R**2 def test_process_not_constant(self): diff --git a/tests/unit/test_parameters/test_size_distribution_parameters.py b/tests/unit/test_parameters/test_size_distribution_parameters.py index 414b422055..a0f7847a32 100644 --- a/tests/unit/test_parameters/test_size_distribution_parameters.py +++ b/tests/unit/test_parameters/test_size_distribution_parameters.py @@ -40,3 +40,76 @@ def test_parameter_values(self): R_test = pybamm.Scalar(1.0) values.evaluate(param.n.prim.f_a_dist(R_test)) values.evaluate(param.p.prim.f_a_dist(R_test)) + + # check that the composite parameters works as well + params_composite = pybamm.ParameterValues("Chen2020_composite") + params_composite.update( + { + "Primary: Positive particle radius [m]": 1e-6, + "Secondary: Positive particle radius [m]": 1e-6, + "Negative particle radius [m]": 1e-6, + }, + check_already_exists=False, + ) + values_composite = pybamm.get_size_distribution_parameters( + params_composite, + composite="negative", + ) + assert "Primary: Negative maximum particle radius [m]" in values_composite + param = pybamm.LithiumIonParameters({"particle phases": ("2", "1")}) + R_test_n_prim = pybamm.Scalar(1e-6) + R_test_n_prim.domains = {"primary": ["negative primary particle size"]} + R_test_n_sec = pybamm.Scalar(1e-6) + R_test_n_sec.domains = {"primary": ["negative secondary particle size"]} + values_composite.evaluate(param.n.prim.f_a_dist(R_test_n_prim)) + values_composite.evaluate(param.n.sec.f_a_dist(R_test_n_sec)) + params_composite = pybamm.ParameterValues("Chen2020_composite") + params_composite.update( + { + "Primary: Positive particle radius [m]": 1e-6, + "Secondary: Positive particle radius [m]": 1e-6, + "Negative particle radius [m]": 1e-6, + }, + check_already_exists=False, + ) + values_composite = pybamm.get_size_distribution_parameters( + params_composite, + composite="positive", + ) + assert "Primary: Positive maximum particle radius [m]" in values_composite + param = pybamm.LithiumIonParameters({"particle phases": ("1", "2")}) + R_test_p_prim = pybamm.Scalar(1e-6) + R_test_p_prim.domains = {"primary": ["positive primary particle size"]} + R_test_p_sec = pybamm.Scalar(1e-6) + R_test_p_sec.domains = {"primary": ["positive secondary particle size"]} + values_composite.evaluate(param.p.prim.f_a_dist(R_test_p_prim)) + values_composite.evaluate(param.p.sec.f_a_dist(R_test_p_sec)) + params_composite = pybamm.ParameterValues("Chen2020_composite") + params_composite.update( + { + "Primary: Positive particle radius [m]": 1e-6, + "Secondary: Positive particle radius [m]": 1e-6, + "Negative particle radius [m]": 1e-6, + }, + check_already_exists=False, + ) + values_composite = pybamm.get_size_distribution_parameters( + params_composite, + composite="both", + ) + assert "Primary: Negative maximum particle radius [m]" in values_composite + assert "Primary: Positive maximum particle radius [m]" in values_composite + params_composite = pybamm.ParameterValues("Chen2020_composite") + params_composite.update( + { + "Primary: Positive particle radius [m]": 1e-6, + "Secondary: Positive particle radius [m]": 1e-6, + "Negative particle radius [m]": 1e-6, + }, + check_already_exists=False, + ) + values_composite = pybamm.get_size_distribution_parameters( + params_composite, + ) + assert "Primary: Negative maximum particle radius [m]" not in values_composite + assert "Primary: Positive maximum particle radius [m]" not in values_composite diff --git a/tests/unit/test_plotting/test_plot_summary_variables.py b/tests/unit/test_plotting/test_plot_summary_variables.py index 5f1a650ced..418eeaf71a 100644 --- a/tests/unit/test_plotting/test_plot_summary_variables.py +++ b/tests/unit/test_plotting/test_plot_summary_variables.py @@ -45,9 +45,9 @@ def test_plot(self): cycle_number, var = ax.get_lines()[0].get_data() np.testing.assert_array_equal( - cycle_number, sol.summary_variables["Cycle number"] + cycle_number, sol.summary_variables.cycle_number ) - np.testing.assert_array_equal(var, sol.summary_variables[output_var]) + np.testing.assert_allclose(sol.summary_variables[output_var], var) axes = pybamm.plot_summary_variables( [sol, sol], labels=["SPM", "SPM"], show_plot=False @@ -62,12 +62,12 @@ def test_plot(self): cycle_number, var = ax.get_lines()[0].get_data() np.testing.assert_array_equal( - cycle_number, sol.summary_variables["Cycle number"] + cycle_number, sol.summary_variables.cycle_number ) - np.testing.assert_array_equal(var, sol.summary_variables[output_var]) + np.testing.assert_allclose(sol.summary_variables[output_var], var) cycle_number, var = ax.get_lines()[1].get_data() np.testing.assert_array_equal( - cycle_number, sol.summary_variables["Cycle number"] + cycle_number, sol.summary_variables.cycle_number ) - np.testing.assert_array_equal(var, sol.summary_variables[output_var]) + np.testing.assert_allclose(sol.summary_variables[output_var], var) diff --git a/tests/unit/test_solvers/test_summary_variables.py b/tests/unit/test_solvers/test_summary_variables.py new file mode 100644 index 0000000000..0132c42ea5 --- /dev/null +++ b/tests/unit/test_solvers/test_summary_variables.py @@ -0,0 +1,171 @@ +# +# Tests for the Summary Variables class +# + +import pybamm +import numpy as np +import pytest + + +class TestSummaryVariables: + @staticmethod + def create_sum_vars(): + model = pybamm.BaseModel() + c = pybamm.Variable("c") + model.rhs = {c: -c} + model.initial_conditions = {c: 1} + model.variables["c"] = c + model.variables["2c"] = 2 * c + model.summary_variables = ["2c"] + + solution = pybamm.ScipySolver().solve(model, np.linspace(0, 1)) + + sum_vars = pybamm.SummaryVariables(solution) + + return sum_vars, solution + + def test_init(self): + model = pybamm.BaseModel() + c = pybamm.Variable("c") + model.rhs = {c: -c} + model.initial_conditions = {c: 1} + model.variables["c"] = c + model.variables["2c"] = 2 * c + model.summary_variables = ["2c"] + + solution = pybamm.ScipySolver().solve(model, np.linspace(0, 1)) + + sum_vars = pybamm.SummaryVariables(solution) + + # no variables should have been calculated until called + assert sum_vars._variables == {} + + assert sum_vars.first_state == solution.first_state + assert sum_vars.last_state == solution.last_state + assert sum_vars.cycles is None + + def test_init_with_cycle_summary_variables(self): + model = pybamm.BaseModel() + c = pybamm.Variable("c") + model.rhs = {c: -c} + model.initial_conditions = {c: 1} + model.variables["c"] = c + model.variables["2c"] = 2 * c + model.summary_variables = ["2c"] + + sol1 = pybamm.ScipySolver().solve(model, np.linspace(0, 1)) + sol2 = pybamm.ScipySolver().solve(model, np.linspace(1, 2)) + sol3 = pybamm.ScipySolver().solve(model, np.linspace(2, 3)) + + cycle_sol = sol1 + sol2 + sol3 + + all_sum_vars = [ + pybamm.SummaryVariables(sol1), + pybamm.SummaryVariables(sol2), + pybamm.SummaryVariables(sol3), + ] + + cycle_sum_vars = pybamm.SummaryVariables( + cycle_sol, + cycle_summary_variables=all_sum_vars, + ) + + assert cycle_sum_vars.first_state is None + assert cycle_sum_vars.last_state is None + assert cycle_sum_vars._variables == {} + assert cycle_sum_vars.cycles == all_sum_vars + np.testing.assert_array_equal(cycle_sum_vars.cycle_number, np.array([1, 2, 3])) + + def test_get_variable(self): + sum_vars, solution = self.create_sum_vars() + + summary_c = sum_vars["2c"] + + assert summary_c == solution["2c"].data[-1] + assert list(sum_vars._variables.keys()) == ["2c", "Change in 2c"] + + def test_get_esoh_variable(self): + model = pybamm.lithium_ion.SPM() + sim = pybamm.Simulation(model) + sol = sim.solve(np.linspace(0, 1)) + esoh_solver = sim.get_esoh_solver(True) + sum_vars_esoh = pybamm.SummaryVariables(sol, esoh_solver=esoh_solver) + + assert np.isclose(sum_vars_esoh["x_100"], 0.9493) + + # all esoh vars should be calculated at the same time to reduce solver calls + assert "Practical NPR" in sum_vars_esoh._variables + + def test_get_variable_error_not_summary_variable(self): + sum_vars, _ = self.create_sum_vars() + + with pytest.raises(KeyError, match="Variable 'c' is not a summary variable"): + sum_vars["c"] + + def test_summary_vars_all_variables(self): + # no esoh + sum_vars, _ = self.create_sum_vars() + + assert sum_vars.all_variables == ["2c", "Change in 2c"] + + # with esoh + model = pybamm.lithium_ion.SPM() + sim = pybamm.Simulation(model) + sol = sim.solve(np.linspace(0, 1)) + esoh_solver = sim.get_esoh_solver(True) + sum_vars_esoh = pybamm.SummaryVariables(sol, esoh_solver=esoh_solver) + + assert sum_vars_esoh.calc_esoh is True + assert "Total lithium lost [mol]" in sum_vars_esoh.all_variables + assert "x_100" in sum_vars_esoh.all_variables + + assert "x_100" in sum_vars_esoh._esoh_variables + + def test_get_with_cycle_summary_variables(self): + model = pybamm.BaseModel() + c = pybamm.Variable("c") + model.rhs = {c: -c} + model.initial_conditions = {c: 1} + model.variables["c"] = c + model.variables["2c"] = 2 * c + model.summary_variables = ["2c"] + + sol1 = pybamm.ScipySolver().solve(model, np.linspace(0, 1)) + sol2 = pybamm.ScipySolver().solve(model, np.linspace(1, 2)) + sol3 = pybamm.ScipySolver().solve(model, np.linspace(2, 3)) + + cycle_sol = sol1 + sol2 + sol3 + + all_sum_vars = [ + pybamm.SummaryVariables(sol1), + pybamm.SummaryVariables(sol2), + pybamm.SummaryVariables(sol3), + ] + + cycle_sum_vars = pybamm.SummaryVariables( + cycle_sol, + cycle_summary_variables=all_sum_vars, + ) + + np.testing.assert_array_equal(cycle_sum_vars.cycle_number, np.array([1, 2, 3])) + np.testing.assert_allclose( + cycle_sum_vars["2c"], np.array([0.735758, 0.735758, 0.735758]) + ) + + def test_get_esoh_cycle_summary_vars(self): + experiment = pybamm.Experiment( + [ + ( + "Discharge at 1C for 1 sec", + "Charge at 1C for 1 sec", + ), + ] + * 10, + ) + model = pybamm.lithium_ion.SPM() + sim = pybamm.Simulation(model, experiment=experiment) + sol = sim.solve() + + assert len(sol.summary_variables.cycles) == 10 + assert len(sol.summary_variables["x_100"]) == 10 + assert np.isclose(sol.summary_variables["x_100"][0], 0.9493) diff --git a/tests/unit/test_util.py b/tests/unit/test_util.py index 058b7d4a14..1bf9d6db82 100644 --- a/tests/unit/test_util.py +++ b/tests/unit/test_util.py @@ -35,6 +35,7 @@ def test_fuzzy_dict(self): "Lithium plating current": 4, "A dimensional variable [m]": 5, "Positive particle diffusivity [m2.s-1]": 6, + "Primary: Open circuit voltage [V]": 7, } ) d2 = pybamm.FuzzyDict( @@ -52,6 +53,9 @@ def test_fuzzy_dict(self): with pytest.raises(KeyError, match="dimensional version"): d.__getitem__("A dimensional variable") + with pytest.raises(KeyError, match="composite model"): + d.__getitem__("Open circuit voltage [V]") + with pytest.raises(KeyError, match="open circuit voltage"): d.__getitem__("Measured open circuit voltage [V]")