Merge pull request #1058 from pybamm-team/issue-1030-length-scales-mo…

…dels

Issue 1030 length scales models

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
 

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

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):

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

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,

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):

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

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,

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):

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)