Merge pull request #2440 from pybamm-team/dimensional-dfn

Variable scale and reference

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 
+- Added `scale` and `reference` attributes to `Variable` objects, which can be use to make the ODE/DAE solver better conditioned ([#2440](https://github.com/pybamm-team/PyBaMM/pull/2440))
 - SEI reactions can now be asymmetric ([#2425](https://github.com/pybamm-team/PyBaMM/pull/2425))
 
 ## Optimizations

examples/scripts/experimental_protocols/cccv.py

@@ -4,7 +4,7 @@
 import pybamm
 import matplotlib.pyplot as plt
 
-pybamm.set_logging_level("INFO")
+pybamm.set_logging_level("NOTICE")
 experiment = pybamm.Experiment(
     [
         (

pybamm/discretisations/discretisation.py

@@ -60,7 +60,7 @@ def __init__(self, mesh=None, spatial_methods=None):
                             )
                         )
 
-        self.bcs = {}
+        self._bcs = {}
         self.y_slices = {}
         self._discretised_symbols = {}
         self.external_variables = {}
@@ -187,7 +187,7 @@ def process_model(
         pybamm.logger.verbose(
             "Discretise boundary conditions for {}".format(model.name)
         )
-        self.bcs = self.process_boundary_conditions(model)
+        self._bcs = self.process_boundary_conditions(model)
         pybamm.logger.verbose(
             "Set internal boundary conditions for {}".format(model.name)
         )
@@ -489,7 +489,9 @@ def process_initial_conditions(self, model):
 
         """
         # Discretise initial conditions
-        processed_initial_conditions = self.process_dict(model.initial_conditions)
+        processed_initial_conditions = self.process_dict(
+            model.initial_conditions, ics=True
+        )
 
         # Concatenate initial conditions into a single vector
         # check that all initial conditions are set
@@ -728,7 +730,7 @@ def create_mass_matrix(self, model):
 
         return mass_matrix, mass_matrix_inv
 
-    def process_dict(self, var_eqn_dict):
+    def process_dict(self, var_eqn_dict, ics=False):
         """Discretise a dictionary of {variable: equation}, broadcasting if necessary
         (can be model.rhs, model.algebraic, model.initial_conditions or
         model.variables).
@@ -739,6 +741,9 @@ def process_dict(self, var_eqn_dict):
             Equations ({variable: equation} dict) to dicretise
             (can be model.rhs, model.algebraic, model.initial_conditions or
             model.variables)
+        ics : bool, optional
+            Whether the equations are initial conditions. If True, the equations are
+            scaled by the reference value of the variable, if given
 
         Returns
         -------
@@ -756,11 +761,19 @@ def process_dict(self, var_eqn_dict):
                     eqn = pybamm.FullBroadcast(eqn, broadcast_domains=eqn_key.domains)
 
             pybamm.logger.debug("Discretise {!r}".format(eqn_key))
-
             processed_eqn = self.process_symbol(eqn)
 
-            new_var_eqn_dict[eqn_key] = processed_eqn
+            # Calculate scale if the key has a scale
+            scale = getattr(eqn_key, "scale", 1)
+            if ics:
+                reference = getattr(eqn_key, "reference", 0)
+            else:
+                reference = 0
+
+            if scale != 1 or reference != 0:
+                processed_eqn = (processed_eqn - reference) / scale
 
+            new_var_eqn_dict[eqn_key] = processed_eqn
         return new_var_eqn_dict
 
     def process_symbol(self, symbol):
@@ -790,6 +803,7 @@ def process_symbol(self, symbol):
                 discretised_symbol.mesh = self.mesh[symbol.domain]
             else:
                 discretised_symbol.mesh = None
+
             # Assign secondary mesh
             if symbol.domains["secondary"] != []:
                 discretised_symbol.secondary_mesh = self.mesh[
@@ -934,7 +948,9 @@ def _process_symbol(self, symbol):
             return symbol._function_new_copy(disc_children)
 
         elif isinstance(symbol, pybamm.VariableDot):
-            return pybamm.StateVectorDot(
+            # Add symbol's reference and multiply by the symbol's scale
+            # so that the state vector is of order 1
+            return symbol.reference + symbol.scale * pybamm.StateVectorDot(
                 *self.y_slices[symbol.get_variable()],
                 domains=symbol.domains,
             )
@@ -981,11 +997,29 @@ def _process_symbol(self, symbol):
                             symbol.name
                         )
                     )
-                return pybamm.StateVector(*y_slices, domains=symbol.domains)
+                # Add symbol's reference and multiply by the symbol's scale
+                # so that the state vector is of order 1
+                return symbol.reference + symbol.scale * pybamm.StateVector(
+                    *y_slices, domains=symbol.domains
+                )
 
         elif isinstance(symbol, pybamm.SpatialVariable):
             return spatial_method.spatial_variable(symbol)
 
+        elif isinstance(symbol, pybamm.ConcatenationVariable):
+            # create new children without scale and reference
+            # the scale and reference will be applied to the concatenation instead
+            new_children = []
+            for child in symbol.children:
+                child = child.create_copy()
+                child._scale = 1
+                child._reference = 0
+                child.set_id()
+                new_children.append(self.process_symbol(child))
+            new_symbol = spatial_method.concatenation(new_children)
+            # apply scale to the whole concatenation
+            return symbol.reference + symbol.scale * new_symbol
+
         elif isinstance(symbol, pybamm.Concatenation):
             new_children = [self.process_symbol(child) for child in symbol.children]
             new_symbol = spatial_method.concatenation(new_children)

pybamm/expression_tree/concatenations.py

@@ -42,6 +42,7 @@ def __init__(self, *children, name=None, check_domain=True, concat_fun=None):
         else:
             domains = {"primary": []}
         self.concatenation_function = concat_fun
+
         super().__init__(name, children, domains=domains)
 
     def __str__(self):
@@ -364,9 +365,22 @@ def __init__(self, *children):
         name = intersect(children[0].name, children[1].name)
         for child in children[2:]:
             name = intersect(name, child.name)
-        name = name.capitalize()
-        if name == "":
+        if len(name) == 0:
             name = None
+        # name is unchanged if its length is 1
+        elif len(name) > 1:
+            name = name[0].capitalize() + name[1:]
+
+        if len(children) > 0:
+            if all(child.scale == children[0].scale for child in children):
+                self._scale = children[0].scale
+            else:
+                raise ValueError("Cannot concatenate symbols with different scales")
+            if all(child.reference == children[0].reference for child in children):
+                self._reference = children[0].reference
+            else:
+                raise ValueError("Cannot concatenate symbols with different references")
+
         super().__init__(*children, name=name)
         # Overly tight bounds, can edit later if required
         self.bounds = (

pybamm/expression_tree/symbol.py

@@ -200,7 +200,12 @@ class Symbol:
     """
 
     def __init__(
-        self, name, children=None, domain=None, auxiliary_domains=None, domains=None
+        self,
+        name,
+        children=None,
+        domain=None,
+        auxiliary_domains=None,
+        domains=None,
     ):
         super(Symbol, self).__init__()
         self.name = name
@@ -404,6 +409,14 @@ def set_id(self):
             + tuple([(k, tuple(v)) for k, v in self.domains.items() if v != []])
         )
 
+    @property
+    def scale(self):
+        return self._scale
+
+    @property
+    def reference(self):
+        return self._reference
+
     def __eq__(self, other):
         try:
             return self._id == other._id

pybamm/expression_tree/variable.py

@@ -18,7 +18,6 @@ class VariableBase(pybamm.Symbol):
 
     Parameters
     ----------
-
     name : str
         name of the node
     domain : iterable of str
@@ -39,6 +38,14 @@ class VariableBase(pybamm.Symbol):
         Physical bounds on the variable
     print_name : str, optional
         The name to use for printing. Default is None, in which case self.name is used.
+    scale : float or :class:`pybamm.Symbol`, optional
+        The scale of the variable, used for scaling the model when solving. The state
+        vector representing this variable will be multiplied by this scale.
+        Default is 1.
+    reference : float or :class:`pybamm.Symbol`, optional
+        The reference value of the variable, used for scaling the model when solving.
+        This value will be added to the state vector representing this variable.
+        Default is 0.
 
     *Extends:* :class:`Symbol`
     """
@@ -51,9 +58,16 @@ def __init__(
         domains=None,
         bounds=None,
         print_name=None,
+        scale=1,
+        reference=0,
     ):
+        self._scale = scale
+        self._reference = reference
         super().__init__(
-            name, domain=domain, auxiliary_domains=auxiliary_domains, domains=domains
+            name,
+            domain=domain,
+            auxiliary_domains=auxiliary_domains,
+            domains=domains,
         )
         if bounds is None:
             bounds = (-np.inf, np.inf)
@@ -66,13 +80,21 @@ def __init__(
         self.bounds = bounds
         self.print_name = print_name
 
+    def set_id(self):
+        self._id = hash(
+            (self.__class__, self.name, self.scale, self.reference)
+            + tuple([(k, tuple(v)) for k, v in self.domains.items() if v != []])
+        )
+
     def create_copy(self):
         """See :meth:`pybamm.Symbol.new_copy()`."""
         return self.__class__(
             self.name,
             domains=self.domains,
             bounds=self.bounds,
             print_name=self._raw_print_name,
+            scale=self.scale,
+            reference=self.reference,
         )
 
     def _evaluate_for_shape(self):
@@ -117,33 +139,26 @@ class Variable(VariableBase):
         Physical bounds on the variable
     print_name : str, optional
         The name to use for printing. Default is None, in which case self.name is used.
+    scale : float or :class:`pybamm.Symbol`, optional
+        The scale of the variable, used for scaling the model when solving. The state
+        vector representing this variable will be multiplied by this scale.
+        Default is 1.
+    reference : float or :class:`pybamm.Symbol`, optional
+        The reference value of the variable, used for scaling the model when solving.
+        This value will be added to the state vector representing this variable.
+        Default is 0.
 
     *Extends:* :class:`VariableBase`
     """
 
-    def __init__(
-        self,
-        name,
-        domain=None,
-        auxiliary_domains=None,
-        domains=None,
-        bounds=None,
-        print_name=None,
-    ):
-        super().__init__(
-            name,
-            domain=domain,
-            auxiliary_domains=auxiliary_domains,
-            domains=domains,
-            bounds=bounds,
-            print_name=print_name,
-        )
-
     def diff(self, variable):
         if variable == self:
             return pybamm.Scalar(1)
         elif variable == pybamm.t:
-            return pybamm.VariableDot(self.name + "'", domains=self.domains)
+            # reference gets differentiated out
+            return pybamm.VariableDot(
+                self.name + "'", domains=self.domains, scale=self.scale
+            )
         else:
             return pybamm.Scalar(0)
 
@@ -180,36 +195,26 @@ class VariableDot(VariableBase):
         but ignored.
     print_name : str, optional
         The name to use for printing. Default is None, in which case self.name is used.
+    scale : float or :class:`pybamm.Symbol`, optional
+        The scale of the variable, used for scaling the model when solving. The state
+        vector representing this variable will be multiplied by this scale.
+        Default is 1.
+    reference : float or :class:`pybamm.Symbol`, optional
+        The reference value of the variable, used for scaling the model when solving.
+        This value will be added to the state vector representing this variable.
+        Default is 0.
 
     *Extends:* :class:`VariableBase`
     """
 
-    def __init__(
-        self,
-        name,
-        domain=None,
-        auxiliary_domains=None,
-        domains=None,
-        bounds=None,
-        print_name=None,
-    ):
-        super().__init__(
-            name,
-            domain=domain,
-            auxiliary_domains=auxiliary_domains,
-            domains=domains,
-            bounds=bounds,
-            print_name=print_name,
-        )
-
     def get_variable(self):
         """
         return a :class:`.Variable` corresponding to this VariableDot
 
         Note: Variable._jac adds a dash to the name of the corresponding VariableDot, so
         we remove this here
         """
-        return Variable(self.name[:-1], domains=self.domains)
+        return Variable(self.name[:-1], domains=self.domains, scale=self.scale)
 
     def diff(self, variable):
         if variable == self:
@@ -246,13 +251,32 @@ class ExternalVariable(Variable):
         'domain' and 'auxiliary_domains', or just 'domains', should be provided
         (not both). In future, the 'domain' and 'auxiliary_domains' arguments may be
         deprecated.
+    scale : float or :class:`pybamm.Symbol`, optional
+        The scale of the variable, used for scaling the model when solving. The state
+        vector representing this variable will be multiplied by this scale.
+        Default is 1.
+    reference : float or :class:`pybamm.Symbol`, optional
+        The reference value of the variable, used for scaling the model when solving.
+        This value will be added to the state vector representing this variable.
+        Default is 0.
 
     *Extends:* :class:`pybamm.Variable`
     """
 
-    def __init__(self, name, size, domain=None, auxiliary_domains=None, domains=None):
+    def __init__(
+        self,
+        name,
+        size,
+        domain=None,
+        auxiliary_domains=None,
+        domains=None,
+        scale=1,
+        reference=0,
+    ):
         self._size = size
-        super().__init__(name, domain, auxiliary_domains, domains)
+        super().__init__(
+            name, domain, auxiliary_domains, domains, scale=scale, reference=reference
+        )
 
     @property
     def size(self):

pybamm/models/full_battery_models/lead_acid/basic_full.py

@@ -181,7 +181,7 @@ def __init__(self, name="Basic full model"):
         # Current in the electrolyte
         ######################
         i_e = (param.kappa_e(c_e, T) * tor * param.gamma_e / param.C_e) * (
-            param.chiT_over_c(c_e, T) * pybamm.grad(c_e) - pybamm.grad(phi_e)
+            param.chiRT_over_Fc(c_e, T) * pybamm.grad(c_e) - pybamm.grad(phi_e)
         )
         self.algebraic[phi_e] = pybamm.div(i_e) - j
         self.boundary_conditions[phi_e] = {