From 71834c0d4ce6674578e9ce57d8fd6b6abb901fb3 Mon Sep 17 00:00:00 2001 From: Ferran Brosa Planella Date: Mon, 29 Jun 2020 17:30:49 +0200 Subject: [PATCH 1/4] #1087 added extrapolation to endpoints for 2D processed variables --- pybamm/solvers/processed_variable.py | 67 +++++++++++++++++++++++++--- 1 file changed, 61 insertions(+), 6 deletions(-) diff --git a/pybamm/solvers/processed_variable.py b/pybamm/solvers/processed_variable.py index 8e0096f412..6b4e248f60 100644 --- a/pybamm/solvers/processed_variable.py +++ b/pybamm/solvers/processed_variable.py @@ -320,27 +320,82 @@ def initialise_2D(self): order="F", ) + # Get node and edge values + nodes = self.mesh.nodes + edges = self.mesh.edges + if entries.shape[0] == len(nodes): + space = nodes + elif entries.shape[0] == len(edges): + space = edges + + # add points outside first dimension domain for extrapolation to + # boundaries + extrap_space_first_dim_left = np.array([2 * space[0] - space[1]]) + extrap_space_first_dim_right = np.array([2 * space[-1] - space[-2]]) + space_first_dim = np.concatenate( + [extrap_space_first_dim_left, space, extrap_space_first_dim_right] + ) + extrap_entries_left = np.expand_dims(2 * entries[0] - entries[1], axis=0) + extrap_entries_right = np.expand_dims(2 * entries[-1] - entries[-2], axis=0) + entries_for_interp = np.concatenate( + [extrap_entries_left, entries, extrap_entries_right], axis=0 + ) + + # add points outside second dimension domain for extrapolation to + # boundaries + extrap_space_second_dim_left = np.array( + [2 * second_dim_pts[0] - second_dim_pts[1]] + ) + extrap_space_second_dim_right = np.array( + [2 * second_dim_pts[-1] - second_dim_pts[-2]] + ) + space_second_dim = np.concatenate( + [ + extrap_space_second_dim_left, + second_dim_pts, + extrap_space_second_dim_right, + ] + ) + extrap_entries_second_dim_left = np.expand_dims( + 2 * entries_for_interp[:, 0, :] - entries_for_interp[:, 1, :], axis=1 + ) + extrap_entries_second_dim_right = np.expand_dims( + 2 * entries_for_interp[:, -1, :] - entries_for_interp[:, -2, :], axis=1 + ) + entries_for_interp = np.concatenate( + [ + extrap_entries_second_dim_left, + entries_for_interp, + extrap_entries_second_dim_right, + ], + axis=1, + ) + # assign attributes for reference self.entries = entries self.dimensions = 2 first_length_scale = self.get_spatial_scale( self.first_dimension, self.domain[0] ) - self.first_dim_pts = first_dim_pts * first_length_scale + first_dim_pts_for_interp = space_first_dim * first_length_scale second_length_scale = self.get_spatial_scale( self.second_dimension, self.auxiliary_domains["secondary"][0] ) + second_dim_pts_for_interp = space_second_dim * second_length_scale self.second_dim_pts = second_dim_pts * second_length_scale + # Set first_dim_pts to edges for nicer plotting + self.first_dim_pts = edges * first_length_scale + # set up interpolation if len(self.t_sol) == 1: # function of space only. Note the order of the points is the reverse # of what you'd expect interpolant = interp.interp2d( - self.second_dim_pts, - self.first_dim_pts, - entries[:, :, 0], + second_dim_pts_for_interp, + first_dim_pts_for_interp, + entries_for_interp[:, :, 0], kind="linear", fill_value=np.nan, bounds_error=False, @@ -353,8 +408,8 @@ def interp_fun(input): else: # function of space and time. self._interpolation_function = interp.RegularGridInterpolator( - (self.first_dim_pts, self.second_dim_pts, self.t_pts), - entries, + (first_dim_pts_for_interp, second_dim_pts_for_interp, self.t_pts), + entries_for_interp, method="linear", fill_value=np.nan, bounds_error=False, From 5d52fe3c826956b5a17c3f43e0f5827dd6815c39 Mon Sep 17 00:00:00 2001 From: Ferran Brosa Planella Date: Mon, 29 Jun 2020 19:36:38 +0200 Subject: [PATCH 2/4] #1087 fixed tests --- pybamm/solvers/processed_variable.py | 92 +++++++++++++--------------- tests/unit/test_quick_plot.py | 16 +++-- 2 files changed, 54 insertions(+), 54 deletions(-) diff --git a/pybamm/solvers/processed_variable.py b/pybamm/solvers/processed_variable.py index 6b4e248f60..b7d79ec3d7 100644 --- a/pybamm/solvers/processed_variable.py +++ b/pybamm/solvers/processed_variable.py @@ -261,39 +261,14 @@ def initialise_2D(self): """ first_dim_nodes = self.mesh.nodes first_dim_edges = self.mesh.edges - second_dim_pts = self.base_variable.secondary_mesh.nodes - if self.base_eval.size // len(second_dim_pts) == len(first_dim_nodes): + second_dim_nodes = self.base_variable.secondary_mesh.nodes + second_dim_edges = self.base_variable.secondary_mesh.edges + if self.base_eval.size // len(second_dim_nodes) == len(first_dim_nodes): first_dim_pts = first_dim_nodes - elif self.base_eval.size // len(second_dim_pts) == len(first_dim_edges): + elif self.base_eval.size // len(second_dim_nodes) == len(first_dim_edges): first_dim_pts = first_dim_edges - # Process r-x or x-z - if self.domain[0] in [ - "negative particle", - "positive particle", - ] and self.auxiliary_domains["secondary"][0] in [ - "negative electrode", - "positive electrode", - ]: - self.first_dimension = "r" - self.second_dimension = "x" - self.r_sol = first_dim_pts - self.x_sol = second_dim_pts - elif self.domain[0] in [ - "negative electrode", - "separator", - "positive electrode", - ] and self.auxiliary_domains["secondary"] == ["current collector"]: - self.first_dimension = "x" - self.second_dimension = "z" - self.x_sol = first_dim_pts - self.z_sol = second_dim_pts - else: - raise pybamm.DomainError( - "Cannot process 3D object with domain '{}' " - "and auxiliary_domains '{}'".format(self.domain, self.auxiliary_domains) - ) - + second_dim_pts = second_dim_nodes first_dim_size = len(first_dim_pts) second_dim_size = len(second_dim_pts) entries = np.empty((first_dim_size, second_dim_size, len(self.t_sol))) @@ -320,20 +295,12 @@ def initialise_2D(self): order="F", ) - # Get node and edge values - nodes = self.mesh.nodes - edges = self.mesh.edges - if entries.shape[0] == len(nodes): - space = nodes - elif entries.shape[0] == len(edges): - space = edges - # add points outside first dimension domain for extrapolation to # boundaries - extrap_space_first_dim_left = np.array([2 * space[0] - space[1]]) - extrap_space_first_dim_right = np.array([2 * space[-1] - space[-2]]) - space_first_dim = np.concatenate( - [extrap_space_first_dim_left, space, extrap_space_first_dim_right] + extrap_space_first_dim_left = np.array([2 * first_dim_pts[0] - first_dim_pts[1]]) + extrap_space_first_dim_right = np.array([2 * first_dim_pts[-1] - first_dim_pts[-2]]) + first_dim_pts = np.concatenate( + [extrap_space_first_dim_left, first_dim_pts, extrap_space_first_dim_right] ) extrap_entries_left = np.expand_dims(2 * entries[0] - entries[1], axis=0) extrap_entries_right = np.expand_dims(2 * entries[-1] - entries[-2], axis=0) @@ -349,7 +316,7 @@ def initialise_2D(self): extrap_space_second_dim_right = np.array( [2 * second_dim_pts[-1] - second_dim_pts[-2]] ) - space_second_dim = np.concatenate( + second_dim_pts = np.concatenate( [ extrap_space_second_dim_left, second_dim_pts, @@ -371,22 +338,51 @@ def initialise_2D(self): axis=1, ) + # Process r-x or x-z + if self.domain[0] in [ + "negative particle", + "positive particle", + ] and self.auxiliary_domains["secondary"][0] in [ + "negative electrode", + "positive electrode", + ]: + self.first_dimension = "r" + self.second_dimension = "x" + self.r_sol = first_dim_pts + self.x_sol = second_dim_pts + elif self.domain[0] in [ + "negative electrode", + "separator", + "positive electrode", + ] and self.auxiliary_domains["secondary"] == ["current collector"]: + self.first_dimension = "x" + self.second_dimension = "z" + self.x_sol = first_dim_pts + self.z_sol = second_dim_pts + else: + raise pybamm.DomainError( + "Cannot process 3D object with domain '{}' " + "and auxiliary_domains '{}'".format(self.domain, self.auxiliary_domains) + ) + + + # assign attributes for reference self.entries = entries self.dimensions = 2 first_length_scale = self.get_spatial_scale( self.first_dimension, self.domain[0] ) - first_dim_pts_for_interp = space_first_dim * first_length_scale + first_dim_pts_for_interp = first_dim_pts * first_length_scale second_length_scale = self.get_spatial_scale( self.second_dimension, self.auxiliary_domains["secondary"][0] ) - second_dim_pts_for_interp = space_second_dim * second_length_scale - self.second_dim_pts = second_dim_pts * second_length_scale + second_dim_pts_for_interp = second_dim_pts * second_length_scale - # Set first_dim_pts to edges for nicer plotting - self.first_dim_pts = edges * first_length_scale + # Set pts to edges for nicer plotting + self.first_dim_pts = first_dim_edges * first_length_scale + self.second_dim_pts = second_dim_edges * second_length_scale # set up interpolation if len(self.t_sol) == 1: diff --git a/tests/unit/test_quick_plot.py b/tests/unit/test_quick_plot.py index ed69408f1b..2e08582790 100644 --- a/tests/unit/test_quick_plot.py +++ b/tests/unit/test_quick_plot.py @@ -330,7 +330,7 @@ def test_loqs_spme(self): pybamm.QuickPlot(solution, output_variables) # check 2D (space) variables update properly for different time units - c_n = solution["Negative particle concentration [mol.m-3]"].entries + c_n = solution["Negative particle concentration [mol.m-3]"] for unit, scale in zip(["seconds", "minutes", "hours"], [1, 60, 3600]): quick_plot = pybamm.QuickPlot( @@ -342,12 +342,14 @@ def test_loqs_spme(self): qp_data = quick_plot.plots[ ("Negative particle concentration [mol.m-3]",) ][0][1] - np.testing.assert_array_almost_equal(qp_data, c_n[:, :, 0]) + c_n_eval = c_n(t_eval[0], r=c_n.first_dim_pts, x=c_n.second_dim_pts) + np.testing.assert_array_almost_equal(qp_data, c_n_eval) quick_plot.slider_update(t_eval[-1] / scale) qp_data = quick_plot.plots[ ("Negative particle concentration [mol.m-3]",) ][0][1] - np.testing.assert_array_almost_equal(qp_data, c_n[:, :, 1]) + c_n_eval = c_n(t_eval[-1], r=c_n.first_dim_pts, x=c_n.second_dim_pts) + np.testing.assert_array_almost_equal(qp_data, c_n_eval) pybamm.close_plots() @@ -369,7 +371,7 @@ def test_plot_1plus1D_spme(self): # check 2D (x,z space) variables update properly for different time units # Note: these should be the transpose of the entries in the processed variable - c_e = solution["Electrolyte concentration [mol.m-3]"].entries + c_e = solution["Electrolyte concentration [mol.m-3]"] for unit, scale in zip(["seconds", "minutes", "hours"], [1, 60, 3600]): quick_plot = pybamm.QuickPlot( @@ -377,10 +379,12 @@ def test_plot_1plus1D_spme(self): ) quick_plot.plot(0) qp_data = quick_plot.plots[("Electrolyte concentration [mol.m-3]",)][0][1] - np.testing.assert_array_almost_equal(qp_data.T, c_e[:, :, 0]) + c_e_eval = c_e(t_eval[0], x=c_e.first_dim_pts, z=c_e.second_dim_pts) + np.testing.assert_array_almost_equal(qp_data.T, c_e_eval) quick_plot.slider_update(t_eval[-1] / scale) qp_data = quick_plot.plots[("Electrolyte concentration [mol.m-3]",)][0][1] - np.testing.assert_array_almost_equal(qp_data.T, c_e[:, :, -1]) + c_e_eval = c_e(t_eval[-1], x=c_e.first_dim_pts, z=c_e.second_dim_pts) + np.testing.assert_array_almost_equal(qp_data.T, c_e_eval) pybamm.close_plots() From bfd52bf1cbce2f490e2a95700b9bd2999c94e3c4 Mon Sep 17 00:00:00 2001 From: Ferran Brosa Planella Date: Mon, 29 Jun 2020 19:39:42 +0200 Subject: [PATCH 3/4] #1087 flake8 --- pybamm/solvers/processed_variable.py | 10 ++++++---- tests/unit/test_quick_plot.py | 4 +++- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/pybamm/solvers/processed_variable.py b/pybamm/solvers/processed_variable.py index b7d79ec3d7..5dce767dd9 100644 --- a/pybamm/solvers/processed_variable.py +++ b/pybamm/solvers/processed_variable.py @@ -297,8 +297,12 @@ def initialise_2D(self): # add points outside first dimension domain for extrapolation to # boundaries - extrap_space_first_dim_left = np.array([2 * first_dim_pts[0] - first_dim_pts[1]]) - extrap_space_first_dim_right = np.array([2 * first_dim_pts[-1] - first_dim_pts[-2]]) + extrap_space_first_dim_left = np.array( + [2 * first_dim_pts[0] - first_dim_pts[1]] + ) + extrap_space_first_dim_right = np.array( + [2 * first_dim_pts[-1] - first_dim_pts[-2]] + ) first_dim_pts = np.concatenate( [extrap_space_first_dim_left, first_dim_pts, extrap_space_first_dim_right] ) @@ -365,8 +369,6 @@ def initialise_2D(self): "and auxiliary_domains '{}'".format(self.domain, self.auxiliary_domains) ) - - # assign attributes for reference self.entries = entries self.dimensions = 2 diff --git a/tests/unit/test_quick_plot.py b/tests/unit/test_quick_plot.py index 2e08582790..fde665c3b7 100644 --- a/tests/unit/test_quick_plot.py +++ b/tests/unit/test_quick_plot.py @@ -348,7 +348,9 @@ def test_loqs_spme(self): qp_data = quick_plot.plots[ ("Negative particle concentration [mol.m-3]",) ][0][1] - c_n_eval = c_n(t_eval[-1], r=c_n.first_dim_pts, x=c_n.second_dim_pts) + c_n_eval = c_n( + t_eval[-1], r=c_n.first_dim_pts, x=c_n.second_dim_pts + ) np.testing.assert_array_almost_equal(qp_data, c_n_eval) pybamm.close_plots() From 2d5f326309b5f4a052ddceeff7e67d768b09f6c8 Mon Sep 17 00:00:00 2001 From: Ferran Brosa Planella Date: Mon, 29 Jun 2020 20:00:22 +0200 Subject: [PATCH 4/4] #1087 CHANGELOG --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 60f90d399b..67251fd77e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,6 +14,7 @@ ## Bug fixes +- 2D processed variables can now be evaluated at the domain boundaries ([#1088](https://github.com/pybamm-team/PyBaMM/pull/1088)) - Update the default variable points to better capture behaviour in the solid particles in li-ion models ([#1081](https://github.com/pybamm-team/PyBaMM/pull/1081)) - Fix `QuickPlot` to display variables discretised by FEM (in y-z) properly ([#1078](https://github.com/pybamm-team/PyBaMM/pull/1078)) - Add length scales to `EffectiveResistance` models ([#1071](https://github.com/pybamm-team/PyBaMM/pull/1071))