diff --git a/CHANGELOG.md b/CHANGELOG.md index 852cadcc0c..62242ee1d7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,7 +2,8 @@ ## Features -- Add averaging in secondary dimensions ([#1057](https://github.com/pybamm-team/PyBaMM/pull/1057)) +- Added `length_scales` attribute to models ([#1058](https://github.com/pybamm-team/PyBaMM/pull/1058)) +- Added averaging in secondary dimensions ([#1057](https://github.com/pybamm-team/PyBaMM/pull/1057)) ## Optimizations diff --git a/examples/scripts/compare_comsol/compare_comsol_DFN.py b/examples/scripts/compare_comsol/compare_comsol_DFN.py index 876c59f3ac..a9cafb0f1b 100644 --- a/examples/scripts/compare_comsol/compare_comsol_DFN.py +++ b/examples/scripts/compare_comsol/compare_comsol_DFN.py @@ -119,7 +119,7 @@ def myinterp(t): comsol_voltage.secondary_mesh = None # Create comsol model with dictionary of Matrix variables -comsol_model = pybamm.BaseModel() +comsol_model = pybamm.lithium_ion.BaseModel() comsol_model.variables = { "Negative particle surface concentration [mol.m-3]": comsol_c_n_surf, "Electrolyte concentration [mol.m-3]": comsol_c_e, @@ -129,15 +129,13 @@ def myinterp(t): "Electrolyte potential [V]": comsol_phi_e, "Positive electrode potential [V]": comsol_phi_p, "Terminal voltage [V]": comsol_voltage, - # Add spatial variables for use in QuickPlot - "x": pybamm_model.variables["x"], - "x [m]": pybamm_model.variables["x [m]"], } # Make new solution with same t and y comsol_solution = pybamm.Solution(pybamm_solution.t, pybamm_solution.y) -# Update model timescale to match the pybamm model +# Update model scales to match the pybamm model comsol_model.timescale = pybamm_model.timescale +comsol_model.length_scales = pybamm_model.length_scales comsol_solution.model = comsol_model # plot diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py index 0776e40ba8..bb158382b1 100644 --- a/pybamm/models/base_model.py +++ b/pybamm/models/base_model.py @@ -128,6 +128,7 @@ def __init__(self, name="Unnamed model"): # Default timescale is 1 second self.timescale = pybamm.Scalar(1) + self.length_scales = {} @property def name(self): @@ -335,6 +336,7 @@ def new_copy(self, options=None): new_model.use_simplify = self.use_simplify new_model.convert_to_format = self.convert_to_format new_model.timescale = self.timescale + new_model.length_scales = self.length_scales return new_model def update(self, *submodels): diff --git a/pybamm/models/full_battery_models/lead_acid/base_lead_acid_model.py b/pybamm/models/full_battery_models/lead_acid/base_lead_acid_model.py index e961af9232..463a2c8b67 100644 --- a/pybamm/models/full_battery_models/lead_acid/base_lead_acid_model.py +++ b/pybamm/models/full_battery_models/lead_acid/base_lead_acid_model.py @@ -21,6 +21,16 @@ def __init__(self, options=None, name="Unnamed lead-acid model"): # Default timescale is discharge timescale self.timescale = self.param.tau_discharge + + # Set default length scales + self.length_scales = { + "negative electrode": self.param.L_x, + "separator": self.param.L_x, + "positive electrode": self.param.L_x, + "current collector y": self.param.L_y, + "current collector z": self.param.L_z, + } + self.set_standard_output_variables() @property diff --git a/pybamm/models/full_battery_models/lead_acid/basic_full.py b/pybamm/models/full_battery_models/lead_acid/basic_full.py index 9b6408752e..60addb80fe 100644 --- a/pybamm/models/full_battery_models/lead_acid/basic_full.py +++ b/pybamm/models/full_battery_models/lead_acid/basic_full.py @@ -288,8 +288,6 @@ def __init__(self, name="Basic full model"): - param.U_n_ref + pot * phi_s_p, "Terminal voltage [V]": param.U_p_ref - param.U_n_ref + pot * voltage, - "x [m]": pybamm.standard_spatial_vars.x * param.L_x, - "x": pybamm.standard_spatial_vars.x, "Porosity": eps, "Volume-averaged velocity": v, "X-averaged separator transverse volume-averaged velocity": div_V_s, diff --git a/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py b/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py index 5176da2174..bfa6d963c8 100644 --- a/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py +++ b/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py @@ -19,6 +19,17 @@ def __init__(self, options=None, name="Unnamed lithium-ion model"): # Default timescale is discharge timescale self.timescale = self.param.tau_discharge + + # Set default length scales + self.length_scales = { + "negative electrode": self.param.L_x, + "separator": self.param.L_x, + "positive electrode": self.param.L_x, + "negative particle": self.param.R_n, + "positive particle": self.param.R_p, + "current collector y": self.param.L_y, + "current collector z": self.param.L_z, + } self.set_standard_output_variables() def set_standard_output_variables(self): diff --git a/pybamm/parameters/parameter_values.py b/pybamm/parameters/parameter_values.py index d9afc30840..48361f8bfd 100644 --- a/pybamm/parameters/parameter_values.py +++ b/pybamm/parameters/parameter_values.py @@ -386,6 +386,10 @@ def process_model(self, unprocessed_model, inplace=True): # Process timescale model.timescale = self.process_symbol(model.timescale) + # Process length scales + for domain, scale in model.length_scales.items(): + model.length_scales[domain] = self.process_symbol(scale) + pybamm.logger.info("Finish setting parameters for {}".format(model.name)) return model diff --git a/pybamm/plotting/quick_plot.py b/pybamm/plotting/quick_plot.py index ed5649dd42..a965a4a0ba 100644 --- a/pybamm/plotting/quick_plot.py +++ b/pybamm/plotting/quick_plot.py @@ -146,30 +146,11 @@ def __init__( else: raise ValueError("spatial unit '{}' not recognized".format(spatial_unit)) - variables = models[0].variables - # empty spatial scales, will raise error later if can't find a particular one - self.spatial_scales = {} - if "x [m]" in variables and "x" in variables: - x_scale = (variables["x [m]"] / variables["x"]).evaluate()[ - -1 - ] * self.spatial_factor - self.spatial_scales.update({dom: x_scale for dom in variables["x"].domain}) - if "y [m]" in variables and "y" in variables: - self.spatial_scales["current collector y"] = ( - variables["y [m]"] / variables["y"] - ).evaluate()[-1] * self.spatial_factor - if "z [m]" in variables and "z" in variables: - self.spatial_scales["current collector z"] = ( - variables["z [m]"] / variables["z"] - ).evaluate()[-1] * self.spatial_factor - if "r_n [m]" in variables and "r_n" in variables: - self.spatial_scales["negative particle"] = ( - variables["r_n [m]"] / variables["r_n"] - ).evaluate()[-1] * self.spatial_factor - if "r_p [m]" in variables and "r_p" in variables: - self.spatial_scales["positive particle"] = ( - variables["r_p [m]"] / variables["r_p"] - ).evaluate()[-1] * self.spatial_factor + # Set length scales + self.length_scales = { + domain: scale.evaluate() * self.spatial_factor + for domain, scale in models[0].length_scales.items() + } # Time parameters model_timescale_in_seconds = models[0].timescale_eval @@ -269,8 +250,6 @@ def set_output_variables(self, output_variables, solutions): self.spatial_variable_dict = {} self.first_dimensional_spatial_variable = {} self.second_dimensional_spatial_variable = {} - self.first_spatial_scale = {} - self.second_spatial_scale = {} self.is_x_r = {} self.is_y_z = {} @@ -317,18 +296,15 @@ def set_output_variables(self, output_variables, solutions): # Set the x variable (i.e. "x" or "r" for any one-dimensional variables) if first_variable.dimensions == 1: - ( - spatial_var_name, - spatial_var_value, - spatial_scale, - ) = self.get_spatial_var(variable_tuple, first_variable, "first") + (spatial_var_name, spatial_var_value,) = self.get_spatial_var( + variable_tuple, first_variable, "first" + ) self.spatial_variable_dict[variable_tuple] = { spatial_var_name: spatial_var_value } self.first_dimensional_spatial_variable[variable_tuple] = ( spatial_var_value * self.spatial_factor ) - self.first_spatial_scale[variable_tuple] = spatial_scale elif first_variable.dimensions == 2: # Don't allow 2D variables if there are multiple solutions @@ -343,12 +319,10 @@ def set_output_variables(self, output_variables, solutions): ( first_spatial_var_name, first_spatial_var_value, - first_spatial_scale, ) = self.get_spatial_var(variable_tuple, first_variable, "first") ( second_spatial_var_name, second_spatial_var_value, - second_spatial_scale, ) = self.get_spatial_var(variable_tuple, first_variable, "second") self.spatial_variable_dict[variable_tuple] = { first_spatial_var_name: first_spatial_var_value, @@ -397,19 +371,7 @@ def get_spatial_var(self, key, variable, dimension): if domain == "current collector": domain += " {}".format(spatial_var_name) - # Get scale to go from dimensionless to dimensional in the units - # specified by spatial_unit - try: - spatial_scale = self.spatial_scales[domain] - except KeyError: - raise KeyError( - ( - "Can't find spatial scale for '{}', make sure both '{} [m]' " - + "and '{}' are defined in the model variables" - ).format(domain, *[spatial_var_name] * 2) - ) - - return spatial_var_name, spatial_var_value, spatial_scale + return spatial_var_name, spatial_var_value def reset_axis(self): """ @@ -578,10 +540,14 @@ def plot(self, t): # add dashed lines for boundaries between subdomains y_min, y_max = ax.get_ylim() ax.set_ylim(y_min, y_max) - for bnd in variable_lists[0][0].internal_boundaries: - bnd_dim = bnd * self.first_spatial_scale[key] + for boundary in variable_lists[0][0].internal_boundaries: + boundary_scaled = boundary * self.spatial_factor ax.plot( - [bnd_dim, bnd_dim], [y_min, y_max], color="0.5", lw=1, zorder=0 + [boundary_scaled, boundary_scaled], + [y_min, y_max], + color="0.5", + lw=1, + zorder=0, ) elif variable_lists[0][0].dimensions == 2: # Read dictionary of spatial variables @@ -716,10 +682,10 @@ def slider_update(self, t): if y_min is None and y_max is None: y_min, y_max = ax_min(var_min), ax_max(var_max) ax.set_ylim(y_min, y_max) - for bnd in self.variables[key][0][0].internal_boundaries: - bnd_dim = bnd * self.first_spatial_scale[key] + for boundary in self.variables[key][0][0].internal_boundaries: + boundary_scaled = boundary * self.spatial_factor ax.plot( - [bnd_dim, bnd_dim], + [boundary_scaled, boundary_scaled], [y_min, y_max], color="0.5", lw=1, diff --git a/pybamm/solvers/processed_variable.py b/pybamm/solvers/processed_variable.py index bec8946883..8e0096f412 100644 --- a/pybamm/solvers/processed_variable.py +++ b/pybamm/solvers/processed_variable.py @@ -63,18 +63,12 @@ def __init__(self, base_variable, solution, known_evals=None, warn=True): self.timescale = solution.model.timescale.evaluate() self.t_pts = self.t_sol * self.timescale - # Store spatial variables to get scales - self.spatial_vars = {} + # Store length scales if solution.model: - for var in ["x", "y", "z", "r_n", "r_p"]: - if ( - var in solution.model.variables - and var + " [m]" in solution.model.variables - ): - self.spatial_vars[var] = solution.model.variables[var] - self.spatial_vars[var + " [m]"] = solution.model.variables[ - var + " [m]" - ] + self.length_scales = { + domain: scale.evaluate() + for domain, scale in solution.model.length_scales.items() + } # Evaluate base variable at initial time if self.known_evals: @@ -222,16 +216,20 @@ def initialise_1D(self, fixed_t=False): self.x_sol = space # assign attributes for reference - self.first_dim_pts = space * self.get_spatial_scale( - self.first_dimension, self.domain[0] - ) - self.internal_boundaries = self.mesh.internal_boundaries + length_scale = self.get_spatial_scale(self.first_dimension, self.domain[0]) + pts_for_interp = space * length_scale + self.internal_boundaries = [ + bnd * length_scale for bnd in self.mesh.internal_boundaries + ] + + # Set first_dim_pts to edges for nicer plotting + self.first_dim_pts = edges * length_scale # set up interpolation if len(self.t_sol) == 1: # function of space only interpolant = interp.interp1d( - self.first_dim_pts, + pts_for_interp, entries_for_interp[:, 0], kind="linear", fill_value=np.nan, @@ -250,7 +248,7 @@ def interp_fun(t, z): # is the reverse of what you'd expect self._interpolation_function = interp.interp2d( self.t_pts, - self.first_dim_pts, + pts_for_interp, entries_for_interp, kind="linear", fill_value=np.nan, @@ -325,12 +323,15 @@ def initialise_2D(self): # assign attributes for reference self.entries = entries self.dimensions = 2 - self.first_dim_pts = first_dim_pts * self.get_spatial_scale( + first_length_scale = self.get_spatial_scale( self.first_dimension, self.domain[0] ) - self.second_dim_pts = second_dim_pts * self.get_spatial_scale( - self.second_dimension + self.first_dim_pts = first_dim_pts * first_length_scale + + second_length_scale = self.get_spatial_scale( + self.second_dimension, self.auxiliary_domains["secondary"][0] ) + self.second_dim_pts = second_dim_pts * second_length_scale # set up interpolation if len(self.t_sol) == 1: @@ -394,8 +395,8 @@ def initialise_2D_scikit_fem(self): self.z_sol = z_sol self.first_dimension = "y" self.second_dimension = "z" - self.first_dim_pts = y_sol * self.get_spatial_scale("y") - self.second_dim_pts = z_sol * self.get_spatial_scale("z") + self.first_dim_pts = y_sol * self.get_spatial_scale("y", "current collector") + self.second_dim_pts = z_sol * self.get_spatial_scale("z", "current collector") # set up interpolation if len(self.t_sol) == 1: @@ -477,27 +478,22 @@ def call_2D(self, t, x, r, y, z): second_dim = second_dim[:, np.newaxis] return self._interpolation_function((first_dim, second_dim, t)) - def get_spatial_scale(self, name, domain=None): + def get_spatial_scale(self, name, domain): "Returns the spatial scale for a named spatial variable" - # Different scale in negative and positive particles - if domain == "negative particle": - name = "r_n" - elif domain == "positive particle": - name = "r_p" - - # Try to get length scale - if name + " [m]" in self.spatial_vars and name in self.spatial_vars: - scale = ( - self.spatial_vars[name + " [m]"] / self.spatial_vars[name] - ).evaluate()[-1] - else: + try: + if name == "y" and domain == "current collector": + return self.length_scales["current collector y"] + elif name == "z" and domain == "current collector": + return self.length_scales["current collector z"] + else: + return self.length_scales[domain] + except KeyError: if self.warn: pybamm.logger.warning( - "No scale set for spatial variable {}. " - "Using default of 1 [m].".format(name) + "No length scale set for {}. " + "Using default of 1 [m].".format(domain) ) - scale = 1 - return scale + return 1 @property def data(self): diff --git a/tests/integration/test_models/standard_output_tests.py b/tests/integration/test_models/standard_output_tests.py index a20d5f4e5c..4a48d27915 100644 --- a/tests/integration/test_models/standard_output_tests.py +++ b/tests/integration/test_models/standard_output_tests.py @@ -629,13 +629,13 @@ def __init__(self, model, param, disc, solution, operating_condition): def test_velocity_boundaries(self): """Test the boundary values of the current densities""" - L_x = self.v_box.first_dim_pts[-1] / self.v_box.x_sol[-1] + L_x = self.x_edge[-1] np.testing.assert_array_almost_equal(self.v_box(self.t, 0), 0, decimal=4) np.testing.assert_array_almost_equal(self.v_box(self.t, L_x), 0, decimal=4) def test_vertical_velocity(self): """Test the boundary values of the current densities""" - L_x = self.v_box.first_dim_pts[-1] / self.v_box.x_sol[-1] + L_x = self.x_edge[-1] np.testing.assert_array_equal(self.dVbox_dz(self.t, 0), 0) np.testing.assert_array_less(self.dVbox_dz(self.t, 0.5 * L_x), 0) np.testing.assert_array_equal(self.dVbox_dz(self.t, L_x), 0) diff --git a/tests/unit/test_parameters/test_parameter_values.py b/tests/unit/test_parameters/test_parameter_values.py index 94204b6a37..bd85c292db 100644 --- a/tests/unit/test_parameters/test_parameter_values.py +++ b/tests/unit/test_parameters/test_parameter_values.py @@ -511,6 +511,9 @@ def test_process_model(self): "grad_var1": pybamm.grad(var1), "d_var1": d * var1, } + model.timescale = b + model.length_scales = {"test": c} + parameter_values = pybamm.ParameterValues({"a": 1, "b": 2, "c": 3, "d": 42}) parameter_values.process_model(model) # rhs @@ -547,6 +550,9 @@ def test_process_model(self): self.assertTrue( isinstance(model.variables["d_var1"].children[1], pybamm.Variable) ) + # timescale and length scales + self.assertEqual(model.timescale.evaluate(), 2) + self.assertEqual(model.length_scales["test"].evaluate(), 3) # bad boundary conditions model = pybamm.BaseModel() diff --git a/tests/unit/test_quick_plot.py b/tests/unit/test_quick_plot.py index dfda895664..ed69408f1b 100644 --- a/tests/unit/test_quick_plot.py +++ b/tests/unit/test_quick_plot.py @@ -5,7 +5,7 @@ class TestQuickPlot(unittest.TestCase): def test_simple_ode_model(self): - model = pybamm.BaseBatteryModel(name="Simple ODE Model") + model = pybamm.lithium_ion.BaseModel(name="Simple ODE Model") whole_cell = ["negative electrode", "separator", "positive electrode"] # Create variables: domain is explicitly empty since these variables are only @@ -37,13 +37,8 @@ def test_simple_ode_model(self): "c broadcasted positive electrode": pybamm.PrimaryBroadcast( c, "positive particle" ), - "x [m]": pybamm.standard_spatial_vars.x, - "x": pybamm.standard_spatial_vars.x, - "r_n [m]": pybamm.standard_spatial_vars.r_n, - "r_n": pybamm.standard_spatial_vars.r_n, - "r_p [m]": pybamm.standard_spatial_vars.r_p, - "r_p": pybamm.standard_spatial_vars.r_p, } + model.timescale = pybamm.Scalar(1) # ODEs only (don't use jacobian) model.use_jacobian = False @@ -255,13 +250,6 @@ def test_simple_ode_model(self): [solution, solution], ["a"], labels=["sol 1", "sol 2", "sol 3"] ) - # Remove 'x [m]' from the variables and make sure a key error is raise - del solution.model.variables["x [m]"] - with self.assertRaisesRegex( - KeyError, "Can't find spatial scale for 'negative electrode'" - ): - pybamm.QuickPlot(solution, ["b broadcasted"]) - # No variable can be NaN model.variables["NaN variable"] = disc.process_symbol(pybamm.Scalar(np.nan)) with self.assertRaisesRegex( @@ -304,26 +292,30 @@ def test_loqs_spme(self): pybamm.QuickPlot(solution) # check 1D (space) variables update properly for different time units - c_e = solution["Electrolyte concentration [mol.m-3]"].entries + t = solution["Time [s]"].entries + c_e_var = solution["Electrolyte concentration [mol.m-3]"] + # 1D variables should be evaluated on edges + L_x = param.evaluate(pybamm.geometric_parameters.L_x) + c_e = c_e_var(t=t, x=mesh.combine_submeshes(*c_e_var.domain).edges * L_x) for unit, scale in zip(["seconds", "minutes", "hours"], [1, 60, 3600]): quick_plot = pybamm.QuickPlot( solution, ["Electrolyte concentration [mol.m-3]"], time_unit=unit ) quick_plot.plot(0) - # take off extrapolated points + qp_data = ( quick_plot.plots[("Electrolyte concentration [mol.m-3]",)][0][ 0 - ].get_ydata()[1:-1], + ].get_ydata(), )[0] np.testing.assert_array_almost_equal(qp_data, c_e[:, 0]) quick_plot.slider_update(t_eval[-1] / scale) - # take off extrapolated points + qp_data = ( quick_plot.plots[("Electrolyte concentration [mol.m-3]",)][0][ 0 - ].get_ydata()[1:-1], + ].get_ydata(), )[0][:, 0] np.testing.assert_array_almost_equal(qp_data, c_e[:, 1])