Merge pull request #1893 from pybamm-team/issue-1889-scalar-timescale

#1889 fix some parameter bugs and add timescale option

pybamm/expression_tree/operations/replace_symbols.py

@@ -110,13 +110,13 @@ def process_model(self, unprocessed_model, inplace=True):
         ]
 
         # Process timescale
-        model.timescale = self.process_symbol(unprocessed_model.timescale)
+        model._timescale = self.process_symbol(unprocessed_model.timescale)
 
         # Process length scales
         new_length_scales = {}
         for domain, scale in unprocessed_model.length_scales.items():
             new_length_scales[domain] = self.process_symbol(scale)
-        model.length_scales = new_length_scales
+        model._length_scales = new_length_scales
 
         pybamm.logger.info("Finish replacing symbols in {}".format(model.name))
 

pybamm/models/base_model.py

@@ -119,8 +119,8 @@ def __init__(self, name="Unnamed model"):
         self.y_slices = None
 
         # Default timescale is 1 second
-        self.timescale = pybamm.Scalar(1)
-        self.length_scales = {}
+        self._timescale = pybamm.Scalar(1)
+        self._length_scales = {}
 
     @property
     def name(self):
@@ -290,15 +290,15 @@ def timescale(self, value):
     @property
     def length_scales(self):
         "Length scales of model"
-        return self._length_scale
+        return self._length_scales
 
     @length_scales.setter
     def length_scales(self, values):
         "Set the length scale, converting any numbers to pybamm.Scalar"
         for domain, scale in values.items():
             if isinstance(scale, numbers.Number):
                 values[domain] = pybamm.Scalar(scale)
-        self._length_scale = values
+        self._length_scales = values
 
     @property
     def parameters(self):
@@ -376,8 +376,8 @@ def new_empty_copy(self):
         new_model = self.__class__(name=self.name)
         new_model.use_jacobian = self.use_jacobian
         new_model.convert_to_format = self.convert_to_format
-        new_model.timescale = self.timescale
-        new_model.length_scales = self.length_scales
+        new_model._timescale = self.timescale
+        new_model._length_scales = self.length_scales
 
         # Variables from discretisation
         new_model.is_discretised = self.is_discretised

pybamm/models/full_battery_models/base_battery_model.py

@@ -3,6 +3,7 @@
 #
 
 import pybamm
+import numbers
 
 
 class BatteryModelOptions(pybamm.FuzzyDict):
@@ -133,6 +134,9 @@ class BatteryModelOptions(pybamm.FuzzyDict):
             * "thermal" : str
                 Sets the thermal model to use. Can be "isothermal" (default), "lumped",
                 "x-lumped", or "x-full".
+            * "timescale" : str or number
+                Sets the timescale of the model. If "default", the discharge timescale,
+                as defined by other parameters, is used. Otherwise, the number is used.
             * "total interfacial current density as a state" : str
                 Whether to make a state for the total interfacial current density and
                 solve an algebraic equation for it. Default is "false", unless "SEI film
@@ -213,6 +217,7 @@ def __init__(self, extra_options):
             name: options[0] for name, options in self.possible_options.items()
         }
         default_options["external submodels"] = []
+        default_options["timescale"] = "default"
 
         # Change the default for cell geometry based on which thermal option is provided
         extra_options = extra_options or {}
@@ -371,14 +376,16 @@ def __init__(self, extra_options):
                     "mechanics model"
                 )
 
+        # Check options are valid
         for option, value in options.items():
             if option == "external submodels" or option == "working electrode":
                 pass
             else:
                 if isinstance(value, str) or option in [
                     "dimensionality",
                     "operating mode",
-                ]:  # some options don't take strings
+                    "timescale",
+                ]:  # some options accept non-strings
                     value = (value,)
                 else:
                     if not (
@@ -403,7 +410,13 @@ def __init__(self, extra_options):
                             "2-tuples of strings"
                         )
                 for val in value:
-                    if val not in self.possible_options[option]:
+                    if option == "timescale":
+                        if not (val == "default" or isinstance(val, numbers.Number)):
+                            raise pybamm.OptionError(
+                                "'timescale' option must be either 'default' "
+                                "or a number"
+                            )
+                    elif val not in self.possible_options[option]:
                         if not (option == "operating mode" and callable(val)):
                             raise pybamm.OptionError(
                                 f"\n'{val}' is not recognized in option '{option}'. "
@@ -468,6 +481,21 @@ def __init__(self, options=None, name="Unnamed battery model"):
         self._built = False
         self._built_fundamental_and_external = False
 
+    @pybamm.BaseModel.timescale.setter
+    def timescale(self, value):
+        """Set the timescale"""
+        raise NotImplementedError(
+            "Timescale cannot be directly overwritten for this model. "
+            "Pass a timescale to the 'timescale' option instead."
+        )
+
+    @pybamm.BaseModel.length_scales.setter
+    def length_scales(self, value):
+        """Set the length scales"""
+        raise NotImplementedError(
+            "Length scales cannot be directly overwritten for this model. "
+        )
+
     @property
     def default_geometry(self):
         return pybamm.battery_geometry(
@@ -796,8 +824,8 @@ def new_empty_copy(self):
         new_model = self.__class__(name=self.name, options=self.options)
         new_model.use_jacobian = self.use_jacobian
         new_model.convert_to_format = self.convert_to_format
-        new_model.timescale = self.timescale
-        new_model.length_scales = self.length_scales
+        new_model._timescale = self.timescale
+        new_model._length_scales = self.length_scales
         return new_model
 
     @property

pybamm/models/full_battery_models/lead_acid/base_lead_acid_model.py

@@ -23,11 +23,11 @@ def __init__(self, options=None, name="Unnamed lead-acid model", build=False):
         super().__init__(options, name)
         self.param = pybamm.LeadAcidParameters()
 
-        # Default timescale is discharge timescale
-        self.timescale = self.param.tau_discharge
+        # Default timescale
+        self._timescale = self.param.timescale
 
         # Set default length scales
-        self.length_scales = {
+        self._length_scales = {
             "negative electrode": self.param.L_x,
             "separator": self.param.L_x,
             "positive electrode": self.param.L_x,

pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py

@@ -20,11 +20,11 @@ def __init__(self, options=None, name="Unnamed lithium-ion model", build=False):
         # Assess whether the submodel is a half-cell model
         self.half_cell = self.options["working electrode"] != "both"
 
-        # Default timescale is discharge timescale
-        self.timescale = self.param.tau_discharge
+        # Default timescale
+        self._timescale = self.param.timescale
 
         # Set default length scales
-        self.length_scales = {
+        self._length_scales = {
             "negative electrode": self.param.L_x,
             "separator": self.param.L_x,
             "positive electrode": self.param.L_x,

pybamm/models/full_battery_models/lithium_ion/basic_dfn.py

@@ -123,7 +123,7 @@ def __init__(self, name="Doyle-Fuller-Newman model"):
         # right side. This is also accessible via `boundary_value(x, "right")`, with
         # "left" providing the boundary value of the left side
         c_s_surf_n = pybamm.surf(c_s_n)
-        j0_n = param.j0_n(c_e_n, c_s_surf_n, T) / param.C_r_n
+        j0_n = param.gamma_n * param.j0_n(c_e_n, c_s_surf_n, T) / param.C_r_n
         j_n = (
             2
             * j0_n
@@ -167,7 +167,11 @@ def __init__(self, name="Doyle-Fuller-Newman model"):
         self.boundary_conditions[c_s_n] = {
             "left": (pybamm.Scalar(0), "Neumann"),
             "right": (
-                -param.C_n * j_n / param.a_R_n / param.D_n(c_s_surf_n, T),
+                -param.C_n
+                * j_n
+                / param.a_R_n
+                / param.gamma_n
+                / param.D_n(c_s_surf_n, T),
                 "Neumann",
             ),
         }

pybamm/models/full_battery_models/lithium_ion/basic_dfn_half_cell.py

@@ -60,7 +60,7 @@ def __init__(self, options=None, name="Doyle-Fuller-Newman half cell model"):
             R_w_typ = param.R_p_typ
 
         # Set default length scales
-        self.length_scales = {
+        self._length_scales = {
             "separator": param.L_x,
             "positive electrode": param.L_x,
             "positive particle": R_w_typ,
@@ -425,6 +425,6 @@ def new_copy(self, build=False):
         new_model = self.__class__(name=self.name, options=self.options)
         new_model.use_jacobian = self.use_jacobian
         new_model.convert_to_format = self.convert_to_format
-        new_model.timescale = self.timescale
-        new_model.length_scales = self.length_scales
+        new_model._timescale = self.timescale
+        new_model._length_scales = self.length_scales
         return new_model

pybamm/models/full_battery_models/lithium_ion/basic_spm.py

@@ -89,7 +89,11 @@ def __init__(self, name="Single Particle Model"):
         self.boundary_conditions[c_s_n] = {
             "left": (pybamm.Scalar(0), "Neumann"),
             "right": (
-                -param.C_n * j_n / param.a_R_n / param.D_n(c_s_surf_n, T),
+                -param.C_n
+                * j_n
+                / param.a_R_n
+                / param.gamma_n
+                / param.D_n(c_s_surf_n, T),
                 "Neumann",
             ),
         }
@@ -135,7 +139,7 @@ def __init__(self, name="Single Particle Model"):
         # (Some) variables
         ######################
         # Interfacial reactions
-        j0_n = param.j0_n(1, c_s_surf_n, T) / param.C_r_n
+        j0_n = param.gamma_n * param.j0_n(1, c_s_surf_n, T) / param.C_r_n
         j0_p = param.gamma_p * param.j0_p(1, c_s_surf_p, T) / param.C_r_p
         eta_n = (2 / param.ne_n) * pybamm.arcsinh(j_n / (2 * j0_n))
         eta_p = (2 / param.ne_p) * pybamm.arcsinh(j_p / (2 * j0_p))
@@ -144,12 +148,17 @@ def __init__(self, name="Single Particle Model"):
         phi_s_p = eta_p + phi_e + param.U_p(c_s_surf_p, T)
         V = phi_s_p
 
+        pot_scale = self.param.potential_scale
+        U_ref = self.param.U_p_ref - self.param.U_n_ref
+        V_dim = U_ref + pot_scale * V
+
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         # The `variables` dictionary contains all variables that might be useful for
         # visualising the solution of the model
         # Primary broadcasts are used to broadcast scalar quantities across a domain
         # into a vector of the right shape, for multiplying with other vectors
         self.variables = {
+            "Discharge capacity [A.h]": Q,
             "Negative particle surface concentration": pybamm.PrimaryBroadcast(
                 c_s_surf_n, "negative electrode"
             ),
@@ -166,6 +175,7 @@ def __init__(self, name="Single Particle Model"):
                 phi_s_p, "positive electrode"
             ),
             "Terminal voltage": V,
+            "Terminal voltage [V]": V_dim,
         }
         self.events += [
             pybamm.Event("Minimum voltage", V - param.voltage_low_cut),

pybamm/models/submodels/interface/base_interface.py

@@ -110,7 +110,7 @@ def _get_exchange_current_density(self, variables):
             c_s_surf = pybamm.maximum(tol, pybamm.minimum(c_s_surf, 1 - tol))
 
             if self.domain == "Negative":
-                j0 = param.j0_n(c_e, c_s_surf, T) / param.C_r_n
+                j0 = param.gamma_n * param.j0_n(c_e, c_s_surf, T) / param.C_r_n
             elif self.domain == "Positive":
                 j0 = param.gamma_p * param.j0_p(c_e, c_s_surf, T) / param.C_r_p
 

pybamm/models/submodels/particle/no_distribution/fickian_diffusion.py

@@ -73,7 +73,14 @@ def set_boundary_conditions(self, variables):
         R = variables[self.domain + " particle radius"]
 
         if self.domain == "Negative":
-            rbc = -self.param.C_n * j * R / self.param.a_R_n / pybamm.surf(D_eff)
+            rbc = (
+                -self.param.C_n
+                * j
+                * R
+                / self.param.a_R_n
+                / self.param.gamma_n
+                / pybamm.surf(D_eff)
+            )
 
         elif self.domain == "Positive":
             rbc = (

pybamm/models/submodels/particle/no_distribution/polynomial_profile.py

@@ -194,7 +194,7 @@ def set_rhs(self, variables):
         R = variables[self.domain + " particle radius"]
 
         if self.domain == "Negative":
-            self.rhs = {c_s_rav: -3 * j / self.param.a_R_n / R}
+            self.rhs = {c_s_rav: -3 * j / self.param.a_R_n / self.param.gamma_n / R}
 
         elif self.domain == "Positive":
             self.rhs = {c_s_rav: -3 * j / self.param.a_R_p / self.param.gamma_p / R}
@@ -216,7 +216,7 @@ def set_rhs(self, variables):
                         * pybamm.r_average(D_eff)
                         * q_s_rav
                         / self.param.C_n
-                        - 45 * j / self.param.a_R_n / 2
+                        - 45 * j / self.param.a_R_n / self.param.gamma_n / 2
                     }
                 )
             elif self.domain == "Positive":
@@ -247,7 +247,8 @@ def set_algebraic(self, variables):
             if self.domain == "Negative":
                 self.algebraic = {
                     c_s_surf: pybamm.surf(D_eff) * (c_s_surf - c_s_rav)
-                    + self.param.C_n * (j * R / self.param.a_R_n / 5)
+                    + self.param.C_n
+                    * (j * R / self.param.a_R_n / self.param.gamma_n / 5)
                 }
 
             elif self.domain == "Positive":
@@ -266,7 +267,7 @@ def set_algebraic(self, variables):
                 self.algebraic = {
                     c_s_surf: pybamm.surf(D_eff)
                     * (35 * (c_s_surf - c_s_rav) - 8 * q_s_rav)
-                    + self.param.C_n * (j * R / self.param.a_R_n)
+                    + self.param.C_n * (j * R / self.param.a_R_n / self.param.gamma_n)
                 }
 
             elif self.domain == "Positive":

pybamm/models/submodels/particle/no_distribution/x_averaged_fickian_diffusion.py

@@ -90,7 +90,13 @@ def set_boundary_conditions(self, variables):
         ]
 
         if self.domain == "Negative":
-            rbc = -self.param.C_n * j_xav / self.param.a_R_n / pybamm.surf(D_eff_xav)
+            rbc = (
+                -self.param.C_n
+                * j_xav
+                / self.param.a_R_n
+                / self.param.gamma_n
+                / pybamm.surf(D_eff_xav)
+            )
 
         elif self.domain == "Positive":
             rbc = (

pybamm/models/submodels/particle/no_distribution/x_averaged_polynomial_profile.py

@@ -107,7 +107,7 @@ def get_coupled_variables(self, variables):
             if self.domain == "Negative":
                 j_xav = i_boundary_cc / (a_av * self.param.l_n)
                 c_s_surf_xav = c_s_av - self.param.C_n * (
-                    j_xav / 5 / self.param.a_R_n / D_eff_av
+                    j_xav / 5 / self.param.a_R_n / self.param.gamma_n / D_eff_av
                 )
 
             if self.domain == "Positive":
@@ -127,7 +127,8 @@ def get_coupled_variables(self, variables):
                 c_s_surf_xav = (
                     c_s_av
                     + 8 * q_s_av / 35
-                    - self.param.C_n * (j_xav / 35 / self.param.a_R_n / D_eff_av)
+                    - self.param.C_n
+                    * (j_xav / 35 / self.param.a_R_n / self.param.gamma_n / D_eff_av)
                 )
 
             if self.domain == "Positive":
@@ -283,7 +284,11 @@ def set_rhs(self, variables):
         ]
 
         if self.domain == "Negative":
-            self.rhs = {c_s_av: pybamm.source(-3 * j_xav / self.param.a_R_n, c_s_av)}
+            self.rhs = {
+                c_s_av: pybamm.source(
+                    -3 * j_xav / self.param.a_R_n / self.param.gamma_n, c_s_av
+                )
+            }
 
         elif self.domain == "Positive":
             self.rhs = {
@@ -306,7 +311,7 @@ def set_rhs(self, variables):
                     {
                         q_s_av: pybamm.source(
                             -30 * pybamm.surf(D_eff_xav) * q_s_av / self.param.C_n
-                            - 45 * j_xav / self.param.a_R_n / 2,
+                            - 45 * j_xav / self.param.a_R_n / self.param.gamma_n / 2,
                             q_s_av,
                         )
                     }