Skip to content

Commit b1395e5

Browse files
valentinsulzervalentinsulzer
authored
Merge pull request #776 from pybamm-team/issue-771-secondary-broadcast
Issue 771 secondary broadcast
2 parents 1e39620 + a8c241c commit b1395e5

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

43 files changed

+796
-309
lines changed

‎CHANGELOG.md

+2
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,7 @@
22

33
## Features
44

5+
- Added notebook to explain broadcasts ([#776](https://github.com/pybamm-team/PyBaMM/pull/776))
56
- Added a step to discretisation that automatically compute the inverse of the mass matrix of the differential part of the problem so that the underlying DAEs can be provided in semi-explicit form, as required by the CasADi solver ([#769](https://github.com/pybamm-team/PyBaMM/pull/769))
67
- Added the gradient operation for the Finite Element Method ([#767](https://github.com/pybamm-team/PyBaMM/pull/767))
78
- Added `InputParameter` node for quickly changing parameter values ([#752](https://github.com/pybamm-team/PyBaMM/pull/752))
@@ -37,6 +38,7 @@
3738

3839
## Bug fixes
3940

41+
- Improved implementation of broadcasts ([#776](https://github.com/pybamm-team/PyBaMM/pull/776))
4042
- Fixed a bug which meant that the Ohmic heating in the current collectors was incorrect if using the Finite Element Method ([#767](https://github.com/pybamm-team/PyBaMM/pull/767))
4143
- Improved automatic broadcasting ([#747](https://github.com/pybamm-team/PyBaMM/pull/747))
4244
- Fixed bug with wrong temperature in initial conditions ([#737](https://github.com/pybamm-team/PyBaMM/pull/737))

‎INSTALL-LINUX.md

+8-8
Original file line numberDiff line numberDiff line change
@@ -94,25 +94,25 @@ Before installing scikits.odes, you need to have installed:
9494
- Fortran compiler (e.g. gfortran)
9595
- BLAS/LAPACK install (OpenBLAS is recommended by the scikits.odes developers)
9696
- CMake (for building Sundials)
97-
- Sundials 4.1.0
97+
- Sundials 5.0.0
9898

9999
You can install these on Ubuntu or Debian using apt-get:
100100

101101
```bash
102102
sudo apt-get install python3-dev gfortran gcc cmake libopenblas-dev
103103
```
104104

105-
To install Sundials 4.1.0, on the command-line type:
105+
To install Sundials 5.0.0, on the command-line type:
106106

107107
```bash
108108
INSTALL_DIR=`pwd`/sundials
109-
wget https://computation.llnl.gov/projects/sundials/download/sundials-4.1.0.tar.gz
110-
tar -xvf sundials-4.1.0.tar.gz
111-
mkdir build-sundials-4.1.0
112-
cd build-sundials-4.1.0/
113-
cmake -DLAPACK_ENABLE=ON -DSUNDIALS_INDEX_TYPE=int32_t -DBUILD_ARKODE:BOOL=OFF -DEXAMPLES_ENABLE:BOOL=OFF -DCMAKE_INSTALL_PREFIX=$INSTALL_DIR ../sundials-4.1.0/
109+
wget https://computation.llnl.gov/projects/sundials/download/sundials-5.0.0.tar.gz
110+
tar -xvf sundials-5.0.0.tar.gz
111+
mkdir build-sundials-5.0.0
112+
cd build-sundials-5.0.0/
113+
cmake -DLAPACK_ENABLE=ON -DSUNDIALS_INDEX_TYPE=int32_t -DBUILD_ARKODE:BOOL=OFF -DEXAMPLES_ENABLE:BOOL=OFF -DCMAKE_INSTALL_PREFIX=$INSTALL_DIR ../sundials-5.0.0/
114114
make install
115-
rm -r ../sundials-4.1.0
115+
rm -r ../sundials-5.0.0
116116
```
117117

118118
Then install [scikits.odes](https://github.com/bmcage/odes), letting it know the sundials install location:

‎docs/source/expression_tree/broadcasts.rst

+3
Original file line numberDiff line numberDiff line change
@@ -9,3 +9,6 @@ Broadcasting Operators
99

1010
.. autoclass:: pybamm.PrimaryBroadcast
1111
:members:
12+
13+
.. autoclass:: pybamm.SecondaryBroadcast
14+
:members:

‎examples/notebooks/README.md

+3-2
Original file line numberDiff line numberDiff line change
@@ -40,8 +40,9 @@ For more advanced usage, new sets of parameters, spatial methods and solvers can
4040
## Expression tree structure
4141

4242
PyBaMM is built around an expression tree structure.
43-
[This](expression_tree/expression-tree.ipynb) notebook explains how this works, from
44-
model creation to solution.
43+
44+
- [The expression tree notebook](expression_tree/expression-tree.ipynb) explains how this works, from model creation to solution.
45+
- [The broadcast notebook](expression_tree/broadcasts.ipynb) explains the different types of broadcast.
4546

4647
The following notebooks are specific to different stages of the PyBaMM pipeline, such as choosing a model, spatial method, or solver.
4748

‎examples/notebooks/expression_tree/broadcasts.ipynb

+269
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,269 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"metadata": {},
6+
"source": [
7+
"# Broadcasts\n",
8+
"\n",
9+
"This notebook explains the different types of broadcast available in PyBaMM.\n",
10+
"Understanding of the [expression_tree](./expression-tree.ipynb) and [discretisation](../spatial_methods/finite-volumes.ipynb) notebooks is assumed."
11+
]
12+
},
13+
{
14+
"cell_type": "code",
15+
"execution_count": 1,
16+
"metadata": {},
17+
"outputs": [],
18+
"source": [
19+
"import pybamm\n",
20+
"import numpy as np"
21+
]
22+
},
23+
{
24+
"cell_type": "markdown",
25+
"metadata": {},
26+
"source": [
27+
"We also explicitly set up the discretisation that is used for this notebook. We use a small number of points in each domain, in order to easily visualise the results."
28+
]
29+
},
30+
{
31+
"cell_type": "code",
32+
"execution_count": 2,
33+
"metadata": {},
34+
"outputs": [],
35+
"source": [
36+
"var = pybamm.standard_spatial_vars\n",
37+
"geometry = {\n",
38+
" \"negative electrode\": {\"primary\": {var.x_n: {\"min\": pybamm.Scalar(0), \"max\": pybamm.Scalar(1)}}},\n",
39+
" \"negative particle\": {\"primary\": {var.r_n: {\"min\": pybamm.Scalar(0), \"max\": pybamm.Scalar(1)}}},\n",
40+
"\n",
41+
"}\n",
42+
"\n",
43+
"submesh_types = {\n",
44+
" \"negative electrode\": pybamm.Uniform1DSubMesh,\n",
45+
" \"negative particle\": pybamm.Uniform1DSubMesh,\n",
46+
"}\n",
47+
"\n",
48+
"var_pts = {var.x_n: 5, var.r_n: 3}\n",
49+
"mesh = pybamm.Mesh(geometry, submesh_types, var_pts)\n",
50+
"\n",
51+
"spatial_methods = {\n",
52+
" \"negative electrode\": pybamm.FiniteVolume(),\n",
53+
" \"negative particle\": pybamm.FiniteVolume(),\n",
54+
"}\n",
55+
"disc = pybamm.Discretisation(mesh, spatial_methods)"
56+
]
57+
},
58+
{
59+
"cell_type": "markdown",
60+
"metadata": {},
61+
"source": [
62+
"## Primary broadcasts"
63+
]
64+
},
65+
{
66+
"cell_type": "markdown",
67+
"metadata": {},
68+
"source": [
69+
"Primary broadcasts are used to broadcast from a \"larger\" scale to a \"smaller\" scale, for example broadcasting temperature T(x) from the electrode to the particles, or broadcasting current collector current i(y, z) from the current collector to the electrodes.\n",
70+
"To demonstrate this, we first create a variable `T` on the negative electrode domain, discretise it, and evaluate it with a simple linear vector"
71+
]
72+
},
73+
{
74+
"cell_type": "code",
75+
"execution_count": 3,
76+
"metadata": {},
77+
"outputs": [
78+
{
79+
"data": {
80+
"text/plain": [
81+
"array([[0. ],\n",
82+
" [0.25],\n",
83+
" [0.5 ],\n",
84+
" [0.75],\n",
85+
" [1. ]])"
86+
]
87+
},
88+
"execution_count": 3,
89+
"metadata": {},
90+
"output_type": "execute_result"
91+
}
92+
],
93+
"source": [
94+
"T = pybamm.Variable(\"T\", domain=\"negative electrode\")\n",
95+
"disc.set_variable_slices([T])\n",
96+
"disc_T = disc.process_symbol(T)\n",
97+
"disc_T.evaluate(y=np.linspace(0,1,5))"
98+
]
99+
},
100+
{
101+
"cell_type": "markdown",
102+
"metadata": {},
103+
"source": [
104+
"We then broadcast `T` onto the \"negative particle\" domain (using primary broadcast as we are going from the larger electrode scale to the smaller particle scale), and discretise and evaluate the resulting object."
105+
]
106+
},
107+
{
108+
"cell_type": "code",
109+
"execution_count": 4,
110+
"metadata": {},
111+
"outputs": [
112+
{
113+
"data": {
114+
"text/plain": [
115+
"array([[0. ],\n",
116+
" [0. ],\n",
117+
" [0. ],\n",
118+
" [0.25],\n",
119+
" [0.25],\n",
120+
" [0.25],\n",
121+
" [0.5 ],\n",
122+
" [0.5 ],\n",
123+
" [0.5 ],\n",
124+
" [0.75],\n",
125+
" [0.75],\n",
126+
" [0.75],\n",
127+
" [1. ],\n",
128+
" [1. ],\n",
129+
" [1. ]])"
130+
]
131+
},
132+
"execution_count": 4,
133+
"metadata": {},
134+
"output_type": "execute_result"
135+
}
136+
],
137+
"source": [
138+
"primary_broad_T = pybamm.PrimaryBroadcast(T, \"negative particle\")\n",
139+
"disc_T = disc.process_symbol(primary_broad_T)\n",
140+
"disc_T.evaluate(y=np.linspace(0,1,5))"
141+
]
142+
},
143+
{
144+
"cell_type": "markdown",
145+
"metadata": {},
146+
"source": [
147+
"The broadcasted object makes 3 (since the r-grid has 3 points) copies of each element of `T` and stacks them all up to give an object with size 3x5=15. In the resulting vector, the first 3 entries correspond to the 3 points in the r-domain at the first x-grid point (where T=0 uniformly in r), the next 3 entries correspond to the next 3 points in the r-domain at the second x-grid point (where T=0.25 uniformly in r), etc"
148+
]
149+
},
150+
{
151+
"cell_type": "markdown",
152+
"metadata": {},
153+
"source": [
154+
"## Secondary broadcasts"
155+
]
156+
},
157+
{
158+
"cell_type": "markdown",
159+
"metadata": {},
160+
"source": [
161+
"Secondary broadcasts are used to broadcast from a \"smaller\" scale to a \"larger\" scale, for example broadcasting SPM particle concentrations c_s(r) from the particles to the electrodes. Note that this wouldn't be used to broadcast particle concentrations in the DFN, since these already depend on both x and r.\n",
162+
"To demonstrate this, we first create a variable `c_s` on the negative particle domain, discretise it, and evaluate it with a simple linear vector"
163+
]
164+
},
165+
{
166+
"cell_type": "code",
167+
"execution_count": 5,
168+
"metadata": {},
169+
"outputs": [
170+
{
171+
"data": {
172+
"text/plain": [
173+
"array([[0. ],\n",
174+
" [0.5],\n",
175+
" [1. ]])"
176+
]
177+
},
178+
"execution_count": 5,
179+
"metadata": {},
180+
"output_type": "execute_result"
181+
}
182+
],
183+
"source": [
184+
"c_s = pybamm.Variable(\"c_s\", domain=\"negative particle\")\n",
185+
"disc.set_variable_slices([c_s])\n",
186+
"disc_c_s = disc.process_symbol(c_s)\n",
187+
"disc_c_s.evaluate(y=np.linspace(0,1,3))"
188+
]
189+
},
190+
{
191+
"cell_type": "markdown",
192+
"metadata": {},
193+
"source": [
194+
"We then broadcast `c_s` onto the \"negative electrode\" domain (using secondary broadcast as we are going from the smaller particle scale to the large electrode scale), and discretise and evaluate the resulting object."
195+
]
196+
},
197+
{
198+
"cell_type": "code",
199+
"execution_count": 6,
200+
"metadata": {},
201+
"outputs": [
202+
{
203+
"data": {
204+
"text/plain": [
205+
"array([[0. ],\n",
206+
" [0.5],\n",
207+
" [1. ],\n",
208+
" [0. ],\n",
209+
" [0.5],\n",
210+
" [1. ],\n",
211+
" [0. ],\n",
212+
" [0.5],\n",
213+
" [1. ],\n",
214+
" [0. ],\n",
215+
" [0.5],\n",
216+
" [1. ],\n",
217+
" [0. ],\n",
218+
" [0.5],\n",
219+
" [1. ]])"
220+
]
221+
},
222+
"execution_count": 6,
223+
"metadata": {},
224+
"output_type": "execute_result"
225+
}
226+
],
227+
"source": [
228+
"secondary_broad_c_s = pybamm.SecondaryBroadcast(c_s, \"negative electrode\")\n",
229+
"disc_broad_c_s = disc.process_symbol(secondary_broad_c_s)\n",
230+
"disc_broad_c_s.evaluate(y=np.linspace(0,1,3))"
231+
]
232+
},
233+
{
234+
"cell_type": "markdown",
235+
"metadata": {},
236+
"source": [
237+
"The broadcasted object makes 5 (since the x-grid has 5 points) identical copies of the whole variable `c_s` to give an object with size 5x3=15. In the resulting vector, the first 3 entries correspond to the 3 points in the r-domain at the first x-grid point (where c_s varies in r), the next 3 entries correspond to the next 3 points in the r-domain at the second x-grid point (where c_s varies in r), etc"
238+
]
239+
},
240+
{
241+
"cell_type": "code",
242+
"execution_count": null,
243+
"metadata": {},
244+
"outputs": [],
245+
"source": []
246+
}
247+
],
248+
"metadata": {
249+
"kernelspec": {
250+
"display_name": "Python 3",
251+
"language": "python",
252+
"name": "python3"
253+
},
254+
"language_info": {
255+
"codemirror_mode": {
256+
"name": "ipython",
257+
"version": 3
258+
},
259+
"file_extension": ".py",
260+
"mimetype": "text/x-python",
261+
"name": "python",
262+
"nbconvert_exporter": "python",
263+
"pygments_lexer": "ipython3",
264+
"version": "3.7.3"
265+
}
266+
},
267+
"nbformat": 4,
268+
"nbformat_minor": 2
269+
}

‎pybamm/__init__.py

+1
Original file line numberDiff line numberDiff line change
@@ -103,6 +103,7 @@ def version(formatted=False):
103103
from .expression_tree.broadcasts import (
104104
Broadcast,
105105
PrimaryBroadcast,
106+
SecondaryBroadcast,
106107
FullBroadcast,
107108
ones_like,
108109
)

‎pybamm/discretisations/discretisation.py

+6-3
Original file line numberDiff line numberDiff line change
@@ -682,7 +682,11 @@ def process_dict(self, var_eqn_dict):
682682
for eqn_key, eqn in var_eqn_dict.items():
683683
# Broadcast if the equation evaluates to a number(e.g. Scalar)
684684
if eqn.evaluates_to_number() and not isinstance(eqn_key, str):
685-
eqn = pybamm.Broadcast(eqn, eqn_key.domain)
685+
eqn = pybamm.FullBroadcast(
686+
eqn,
687+
eqn_key.domain,
688+
eqn_key.auxiliary_domains,
689+
)
686690

687691
# note we are sending in the key.id here so we don't have to
688692
# keep calling .id
@@ -771,8 +775,7 @@ def _process_symbol(self, symbol):
771775

772776
elif isinstance(symbol, pybamm.Integral):
773777
out = child_spatial_method.integral(child, disc_child)
774-
out.domain = symbol.domain
775-
out.auxiliary_domains = symbol.auxiliary_domains
778+
out.copy_domains(symbol)
776779
return out
777780

778781
elif isinstance(symbol, pybamm.DefiniteIntegralVector):

‎pybamm/expression_tree/binary_operators.py

+5-4
Original file line numberDiff line numberDiff line change
@@ -84,7 +84,9 @@ def __init__(self, name, left, right):
8484
# right child.
8585
if isinstance(self, (pybamm.Outer, pybamm.Kron)):
8686
domain = right.domain
87-
auxiliary_domains = {"secondary": left.domain}
87+
auxiliary_domains = {}
88+
if domain != []:
89+
auxiliary_domains["secondary"] = left.domain
8890
else:
8991
domain = self.get_children_domains(left.domain, right.domain)
9092
auxiliary_domains = self.get_children_auxiliary_domains([left, right])
@@ -165,8 +167,7 @@ def new_copy(self):
165167

166168
# make new symbol, ensure domain(s) remain the same
167169
out = self._binary_new_copy(new_left, new_right)
168-
out.domain = self.domain
169-
out.auxiliary_domains = self.auxiliary_domains
170+
out.copy_domains(self)
170171

171172
return out
172173

@@ -823,7 +824,7 @@ def source(left, right, boundary=False):
823824
"""
824825
# Broadcast if left is number
825826
if isinstance(left, numbers.Number):
826-
left = pybamm.Broadcast(left, "current collector")
827+
left = pybamm.PrimaryBroadcast(left, "current collector")
827828

828829
if left.domain != ["current collector"] or right.domain != ["current collector"]:
829830
raise pybamm.DomainError(

0 commit comments

Comments
 (0)