Issue 1220 optimizations

CHANGELOG.md

@@ -9,6 +9,10 @@
 
 ## Optimizations
 
+-   Avoid unnecessary repeated computations in the solvers ([#1222](https://github.com/pybamm-team/PyBaMM/pull/1222))
+-   Rewrite `Symbol.is_constant` to be more efficient ([#1222](https://github.com/pybamm-team/PyBaMM/pull/1222))
+-   Cache shape and size calculations ([#1222](https://github.com/pybamm-team/PyBaMM/pull/1222))
+-   Only instantiate the geometric, electrical and thermal parameter classes once ([#1222](https://github.com/pybamm-team/PyBaMM/pull/1222))
 
 ## Bug fixes
 

pybamm/__init__.py

@@ -164,9 +164,12 @@ def version(formatted=False):
 #
 from .parameters.parameter_values import ParameterValues
 from .parameters import constants
-from .parameters.geometric_parameters import GeometricParameters
-from .parameters.electrical_parameters import ElectricalParameters
-from .parameters.thermal_parameters import ThermalParameters
+from .parameters.geometric_parameters import geometric_parameters, GeometricParameters
+from .parameters.electrical_parameters import (
+    electrical_parameters,
+    ElectricalParameters,
+)
+from .parameters.thermal_parameters import thermal_parameters, ThermalParameters
 from .parameters.lithium_ion_parameters import LithiumIonParameters
 from .parameters.lead_acid_parameters import LeadAcidParameters
 from .parameters import parameter_sets

pybamm/expression_tree/array.py

@@ -105,6 +105,10 @@ def _base_evaluate(self, t=None, y=None, y_dot=None, inputs=None):
         """ See :meth:`pybamm.Symbol._base_evaluate()`. """
         return self._entries
 
+    def is_constant(self):
+        """ See :meth:`pybamm.Symbol.is_constant()`. """
+        return True
+
 
 def linspace(start, stop, num=50, **kwargs):
     """

pybamm/expression_tree/binary_operators.py

@@ -218,6 +218,10 @@ def evaluates_on_edges(self, dimension):
             dimension
         )
 
+    def is_constant(self):
+        """ See :meth:`pybamm.Symbol.is_constant()`. """
+        return self.left.is_constant() and self.right.is_constant()
+
 
 class Power(BinaryOperator):
     """A node in the expression tree representing a `**` power operator

pybamm/expression_tree/concatenations.py

@@ -103,6 +103,10 @@ def _evaluate_for_shape(self):
                 [child.evaluate_for_shape() for child in self.children]
             )
 
+    def is_constant(self):
+        """ See :meth:`pybamm.Symbol.is_constant()`. """
+        return all(child.is_constant() for child in self.children)
+
 
 class NumpyConcatenation(Concatenation):
     """A node in the expression tree representing a concatenation of equations, when we

pybamm/expression_tree/functions.py

@@ -174,6 +174,10 @@ def evaluates_on_edges(self, dimension):
         """ See :meth:`pybamm.Symbol.evaluates_on_edges()`. """
         return any(child.evaluates_on_edges(dimension) for child in self.children)
 
+    def is_constant(self):
+        """ See :meth:`pybamm.Symbol.is_constant()`. """
+        return all(child.is_constant() for child in self.children)
+
     def _evaluate_for_shape(self):
         """
         Default behaviour: has same shape as all child

pybamm/expression_tree/input_parameter.py

@@ -36,6 +36,10 @@ def set_expected_size(self, size):
         "Specify the size that the input parameter should be"
         self._expected_size = size
 
+        # We also need to update the saved size and shape
+        self._saved_size = size
+        self._saved_shape = (size, 1)
+
     def _evaluate_for_shape(self):
         """
         Returns the scalar 'NaN' to represent the shape of a parameter.

pybamm/expression_tree/parameter.py

@@ -35,6 +35,10 @@ def _evaluate_for_shape(self):
         """
         return np.nan
 
+    def is_constant(self):
+        """ See :meth:`pybamm.Symbol.is_constant()`. """
+        return True
+
 
 class FunctionParameter(pybamm.Symbol):
     """A node in the expression tree representing a function parameter
@@ -60,7 +64,10 @@ class FunctionParameter(pybamm.Symbol):
     """
 
     def __init__(
-        self, name, inputs, diff_variable=None,
+        self,
+        name,
+        inputs,
+        diff_variable=None,
     ):
         # assign diff variable
         self.diff_variable = diff_variable

pybamm/expression_tree/scalar.py

@@ -24,9 +24,6 @@ class Scalar(pybamm.Symbol):
     """
 
     def __init__(self, value, name=None, domain=[]):
-        """
-
-        """
         # set default name if not provided
         self.value = value
         if name is None:
@@ -62,3 +59,7 @@ def _jac(self, variable):
     def new_copy(self):
         """ See :meth:`pybamm.Symbol.new_copy()`. """
         return Scalar(self.value, self.name, self.domain)
+
+    def is_constant(self):
+        """ See :meth:`pybamm.Symbol.is_constant()`. """
+        return True

pybamm/expression_tree/symbol.py

@@ -174,8 +174,8 @@ def domain(self, domain):
 
     @property
     def auxiliary_domains(self):
-        "Returns domains that are not the primary domain"
-        return {k: v for k, v in self._domains.items() if k != "primary"}
+        "Returns auxiliary domains"
+        return self._auxiliary_domains
 
     @auxiliary_domains.setter
     def auxiliary_domains(self, auxiliary_domains):
@@ -193,6 +193,7 @@ def auxiliary_domains(self, auxiliary_domains):
         if len(set(values)) != len(values):
             raise pybamm.DomainError("All auxiliary domains must be different")
 
+        self._auxiliary_domains = auxiliary_domains
         self._domains.update(auxiliary_domains)
 
     @property
@@ -203,10 +204,16 @@ def secondary_domain(self):
     def copy_domains(self, symbol):
         "Copy the domains from a given symbol, bypassing checks"
         self._domains = symbol.domains
+        self._domain = self._domains["primary"]
+        self._auxiliary_domains = {
+            k: v for k, v in self._domains.items() if k != "primary"
+        }
 
     def clear_domains(self):
         "Clear domains, bypassing checks"
         self._domains = {"primary": []}
+        self._domain = []
+        self._auxiliary_domains = {}
 
     def get_children_auxiliary_domains(self, children):
         "Combine auxiliary domains from children, at all levels"
@@ -292,13 +299,13 @@ def visualise(self, filename):
             DotExporter(
                 new_node, nodeattrfunc=lambda node: 'label="{}"'.format(node.label)
             ).to_picture(filename)
-        except FileNotFoundError:
+        except FileNotFoundError:  # pragma: no cover
             # raise error but only through logger so that test passes
             pybamm.logger.error("Please install graphviz>=2.42.2 to use dot exporter")
 
     def relabel_tree(self, symbol, counter):
-        """ Finds all children of a symbol and assigns them a new id so that they can be
-                visualised properly using the graphviz output
+        """Finds all children of a symbol and assigns them a new id so that they can be
+        visualised properly using the graphviz output
         """
         name = symbol.name
         if name == "div":
@@ -602,24 +609,15 @@ def _evaluate_for_shape(self):
 
     def is_constant(self):
         """returns true if evaluating the expression is not dependent on `t` or `y`
-        or `u`
+        or `inputs`
 
         See Also
         --------
         evaluate : evaluate the expression
 
         """
-        # if any of the nodes are instances of any of these types, then the whole
-        # expression depends on either t or y or u
-        search_types = (
-            pybamm.Variable,
-            pybamm.StateVector,
-            pybamm.Time,
-            pybamm.InputParameter,
-        )
-
-        # do the search, return true if no relevent nodes are found
-        return not self.has_symbol_of_classes(search_types)
+        # Default behaviour is False
+        return False
 
     def evaluate_ignoring_errors(self, t=0):
         """
@@ -736,37 +734,43 @@ def size(self):
         """
         Size of an object, found by evaluating it with appropriate t and y
         """
-        return np.prod(self.shape)
+        try:
+            return self._saved_size
+        except AttributeError:
+            self._saved_size = np.prod(self.shape)
+            return self._saved_size
 
     @property
     def shape(self):
         """
         Shape of an object, found by evaluating it with appropriate t and y.
         """
-        # Default behaviour is to try to evaluate the object directly
-        # Try with some large y, to avoid having to unpack (slow)
         try:
-            y = np.nan * np.ones((1000, 1))
-            evaluated_self = self.evaluate(0, y, y, inputs="shape test")
-        # If that fails, fall back to calculating how big y should really be
-        except ValueError:
-            unpacker = pybamm.SymbolUnpacker(pybamm.StateVector)
-            state_vectors_in_node = unpacker.unpack_symbol(self).values()
-            if state_vectors_in_node == []:
-                y = None
-            else:
+            return self._saved_shape
+        except AttributeError:
+            # Default behaviour is to try to evaluate the object directly
+            # Try with some large y, to avoid having to unpack (slow)
+            try:
+                y = np.nan * np.ones((1000, 1))
+                evaluated_self = self.evaluate(0, y, y, inputs="shape test")
+            # If that fails, fall back to calculating how big y should really be
+            except ValueError:
+                unpacker = pybamm.SymbolUnpacker(pybamm.StateVector)
+                state_vectors_in_node = unpacker.unpack_symbol(self).values()
                 min_y_size = max(
-                    len(x._evaluation_array) for x in state_vectors_in_node
+                    max(len(x._evaluation_array) for x in state_vectors_in_node), 1
                 )
                 # Pick a y that won't cause RuntimeWarnings
                 y = np.nan * np.ones((min_y_size, 1))
-            evaluated_self = self.evaluate(0, y, y, inputs="shape test")
+                evaluated_self = self.evaluate(0, y, y, inputs="shape test")
 
-        # Return shape of evaluated object
-        if isinstance(evaluated_self, numbers.Number):
-            return ()
-        else:
-            return evaluated_self.shape
+            # Return shape of evaluated object
+            if isinstance(evaluated_self, numbers.Number):
+                self._saved_shape = ()
+            else:
+                self._saved_shape = evaluated_self.shape
+
+        return self._saved_shape
 
     @property
     def size_for_testing(self):

pybamm/expression_tree/unary_operators.py

@@ -87,6 +87,10 @@ def evaluates_on_edges(self, dimension):
         """ See :meth:`pybamm.Symbol.evaluates_on_edges()`. """
         return self.child.evaluates_on_edges(dimension)
 
+    def is_constant(self):
+        """ See :meth:`pybamm.Symbol.is_constant()`. """
+        return self.child.is_constant()
+
 
 class Negate(UnaryOperator):
     """A node in the expression tree representing a `-` negation operator
@@ -1147,14 +1151,14 @@ def x_average(symbol):
         if a.id == b.id == c.id:
             return a
         else:
-            geo = pybamm.GeometricParameters()
+            geo = pybamm.geometric_parameters
             l_n = geo.l_n
             l_s = geo.l_s
             l_p = geo.l_p
             return (l_n * a + l_s * b + l_p * c) / (l_n + l_s + l_p)
     # Otherwise, use Integral to calculate average value
     else:
-        geo = pybamm.GeometricParameters()
+        geo = pybamm.geometric_parameters
         if symbol.domain == ["negative electrode"]:
             x = pybamm.standard_spatial_vars.x_n
             l = geo.l_n
@@ -1214,7 +1218,7 @@ def z_average(symbol):
         return symbol.orphans[0]
     # Otherwise, use Integral to calculate average value
     else:
-        geo = pybamm.GeometricParameters()
+        geo = pybamm.geometric_parameters
         z = pybamm.standard_spatial_vars.z
         l_z = geo.l_z
         return Integral(symbol, z) / l_z
@@ -1251,7 +1255,7 @@ def yz_average(symbol):
         return symbol.orphans[0]
     # Otherwise, use Integral to calculate average value
     else:
-        geo = pybamm.GeometricParameters()
+        geo = pybamm.geometric_parameters
         y = pybamm.standard_spatial_vars.y
         z = pybamm.standard_spatial_vars.z
         l_y = geo.l_y

pybamm/geometry/battery_geometry.py

@@ -22,7 +22,7 @@ def battery_geometry(include_particles=True, current_collector_dimension=0):
 
     """
     var = pybamm.standard_spatial_vars
-    geo = pybamm.GeometricParameters()
+    geo = pybamm.geometric_parameters
     l_n = geo.l_n
     l_s = geo.l_s
 

pybamm/meshes/meshes.py

@@ -94,7 +94,18 @@ def __init__(self, geometry, submesh_types, var_pts):
                 else:
                     for lim, sym in spatial_limits.items():
                         if isinstance(sym, pybamm.Symbol):
-                            sym_eval = sym.evaluate()
+                            try:
+                                sym_eval = sym.evaluate()
+                            except NotImplementedError as error:
+                                if sym.has_symbol_of_classes(pybamm.Parameter):
+                                    raise pybamm.DiscretisationError(
+                                        "Parameter values have not yet been set for "
+                                        "geometry. Make sure that something like "
+                                        "`param.process_geometry(geometry)` has been "
+                                        "run."
+                                    )
+                                else:
+                                    raise error
                         elif isinstance(sym, numbers.Number):
                             sym_eval = sym
                         geometry[domain][spatial_variable][lim] = sym_eval

pybamm/parameters/electrical_parameters.py

@@ -17,7 +17,7 @@ class ElectricalParameters:
     def __init__(self):
 
         # Get geometric parameters
-        self.geo = pybamm.GeometricParameters()
+        self.geo = pybamm.geometric_parameters
 
         # Set parameters
         self._set_dimensional_parameters()
@@ -64,3 +64,6 @@ def _set_dimensionless_parameters(self):
             / self.I_typ
             * pybamm.Function(np.sign, self.I_typ)
         )
+
+
+electrical_parameters = ElectricalParameters()

pybamm/parameters/geometric_parameters.py

@@ -101,3 +101,6 @@ def _set_dimensionless_parameters(self):
         self.l_tab_p = self.L_tab_p / self.L_z
         self.centre_y_tab_p = self.Centre_y_tab_p / self.L_z
         self.centre_z_tab_p = self.Centre_z_tab_p / self.L_z
+
+
+geometric_parameters = GeometricParameters()

pybamm/parameters/lead_acid_parameters.py

@@ -22,9 +22,9 @@ class LeadAcidParameters:
     def __init__(self):
 
         # Get geometric, electrical and thermal parameters
-        self.geo = pybamm.GeometricParameters()
-        self.elec = pybamm.ElectricalParameters()
-        self.therm = pybamm.ThermalParameters()
+        self.geo = pybamm.geometric_parameters
+        self.elec = pybamm.electrical_parameters
+        self.therm = pybamm.thermal_parameters
 
         # Set parameters and scales
         self._set_dimensional_parameters()

pybamm/parameters/lithium_ion_parameters.py

@@ -36,9 +36,9 @@ def __init__(self, options=None):
         self.options = options
 
         # Get geometric, electrical and thermal parameters
-        self.geo = pybamm.GeometricParameters()
-        self.elec = pybamm.ElectricalParameters()
-        self.therm = pybamm.ThermalParameters()
+        self.geo = pybamm.geometric_parameters
+        self.elec = pybamm.electrical_parameters
+        self.therm = pybamm.thermal_parameters
 
         # Set parameters and scales
         self._set_dimensional_parameters()