Issue 632 average second dimension

examples/scripts/compare_lead_acid.py

@@ -7,10 +7,10 @@
 
 # load models
 models = [
-    pybamm.lead_acid.LOQS(),
-    pybamm.lead_acid.FOQS(),
-    pybamm.lead_acid.CompositeExtended(),
-    pybamm.lead_acid.Full(),
+    # pybamm.lead_acid.LOQS(),
+    # pybamm.lead_acid.FOQS(),
+    pybamm.lead_acid.Composite(),
+    # pybamm.lead_acid.Full(),
 ]
 
 # create and run simulations

examples/scripts/compare_lithium_ion_2D.py

@@ -18,10 +18,10 @@
 # load models
 models = [
     pybamm.lithium_ion.SPM(
-        {"current collector": "potential pair", "dimensionality": 1}, name="2+1D SPM"
+        {"current collector": "potential pair", "dimensionality": 1}, name="1+1D SPM"
     ),
     pybamm.lithium_ion.SPMe(
-        {"current collector": "potential pair", "dimensionality": 1}, name="2+1D SPMe"
+        {"current collector": "potential pair", "dimensionality": 1}, name="1+1D SPMe"
     ),
 ]
 

pybamm/discretisations/discretisation.py

@@ -47,8 +47,18 @@ def __init__(self, mesh=None, spatial_methods=None):
                 spatial_methods["positive electrode"] = method
 
             self._spatial_methods = spatial_methods
-            for method in self._spatial_methods.values():
+            for domain, method in self._spatial_methods.items():
                 method.build(mesh)
+                # Check zero-dimensional methods are only applied to zero-dimensional
+                # meshes
+                if isinstance(method, pybamm.ZeroDimensionalSpatialMethod):
+                    if not isinstance(mesh[domain], pybamm.SubMesh0D):
+                        raise pybamm.DiscretisationError(
+                            "Zero-dimensional spatial method for the "
+                            "{} domain requires a zero-dimensional submesh".format(
+                                domain
+                            )
+                        )
 
         self.bcs = {}
         self.y_slices = {}
@@ -852,13 +862,18 @@ def _process_symbol(self, symbol):
                 )
 
             elif isinstance(symbol, pybamm.Integral):
-                out = child_spatial_method.integral(child, disc_child)
+                integral_spatial_method = self.spatial_methods[
+                    symbol.integration_variable[0].domain[0]
+                ]
+                out = integral_spatial_method.integral(
+                    child, disc_child, symbol._integration_dimension
+                )
                 out.copy_domains(symbol)
                 return out
 
             elif isinstance(symbol, pybamm.DefiniteIntegralVector):
                 return child_spatial_method.definite_integral_matrix(
-                    child.domains, vector_type=symbol.vector_type
+                    child, vector_type=symbol.vector_type
                 )
 
             elif isinstance(symbol, pybamm.BoundaryIntegral):
@@ -934,7 +949,7 @@ def _process_symbol(self, symbol):
                         domain=parent.domain,
                         auxiliary_domains=parent.auxiliary_domains,
                     )
-                    out = ext[slice(start, end)]
+                    out = pybamm.Index(ext, slice(start, end))
                     out.domain = symbol.domain
                     return out
 
@@ -1104,10 +1119,9 @@ def check_initial_conditions_rhs(self, model):
             assert (
                 model.rhs[var].shape == model.initial_conditions[var].shape
             ), pybamm.ModelError(
-                """
-                rhs and initial_conditions must have the same shape after discretisation
-                but rhs.shape = {} and initial_conditions.shape = {} for variable '{}'.
-                """.format(
+                "rhs and initial_conditions must have the same shape after "
+                "discretisation but rhs.shape = "
+                "{} and initial_conditions.shape = {} for variable '{}'.".format(
                     model.rhs[var].shape, model.initial_conditions[var].shape, var
                 )
             )
@@ -1145,17 +1159,20 @@ def check_variables(self, model):
                 not_concatenation = not isinstance(var, pybamm.Concatenation)
 
                 not_mult_by_one_vec = not (
-                    isinstance(var, pybamm.Multiplication)
-                    and isinstance(var.right, pybamm.Vector)
-                    and np.all(var.right.entries == 1)
+                    isinstance(
+                        var, (pybamm.Multiplication, pybamm.MatrixMultiplication)
+                    )
+                    and (
+                        pybamm.is_matrix_one(var.left)
+                        or pybamm.is_matrix_one(var.right)
+                    )
                 )
 
                 if different_shapes and not_concatenation and not_mult_by_one_vec:
                     raise pybamm.ModelError(
-                        """
-                    variable and its eqn must have the same shape after discretisation
-                    but variable.shape = {} and rhs.shape = {} for variable '{}'.
-                    """.format(
+                        "variable and its eqn must have the same shape after "
+                        "discretisation but variable.shape = "
+                        "{} and rhs.shape = {} for variable '{}'. ".format(
                             var.shape, model.rhs[rhs_var].shape, var
                         )
                     )

pybamm/expression_tree/binary_operators.py

@@ -43,6 +43,19 @@ def is_scalar_one(expr):
         return False
 
 
+def is_matrix_one(expr):
+    """
+    Utility function to test if an expression evaluates to a constant matrix one
+    """
+    if expr.is_constant():
+        result = expr.evaluate_ignoring_errors(t=None)
+        return (issparse(result) and np.all(result.toarray() == 1)) or (
+            isinstance(result, np.ndarray) and np.all(result == 1)
+        )
+    else:
+        return False
+
+
 def zeros_of_shape(shape):
     """
     Utility function to create a scalar zero, or a vector or matrix of zeros of
@@ -199,9 +212,11 @@ def _binary_evaluate(self, left, right):
         """ Perform binary operation on nodes 'left' and 'right'. """
         raise NotImplementedError
 
-    def evaluates_on_edges(self):
+    def evaluates_on_edges(self, dimension):
         """ See :meth:`pybamm.Symbol.evaluates_on_edges()`. """
-        return self.left.evaluates_on_edges() or self.right.evaluates_on_edges()
+        return self.left.evaluates_on_edges(dimension) or self.right.evaluates_on_edges(
+            dimension
+        )
 
 
 class Power(BinaryOperator):
@@ -635,7 +650,7 @@ def _binary_simplify(self, left, right):
 
         return pybamm.simplify_multiplication_division(self.__class__, left, right)
 
-    def evaluates_on_edges(self):
+    def evaluates_on_edges(self, dimension):
         """ See :meth:`pybamm.Symbol.evaluates_on_edges()`. """
         return False
 

pybamm/expression_tree/broadcasts.py

@@ -150,7 +150,7 @@ def __init__(self, child, broadcast_domain, name=None):
         super().__init__(child, broadcast_domain, name)
         self.broadcast_type = "primary to edges"
 
-    def evaluates_on_edges(self):
+    def evaluates_on_edges(self, dimension):
         return True
 
 
@@ -244,7 +244,7 @@ def __init__(self, child, broadcast_domain, name=None):
         super().__init__(child, broadcast_domain, name)
         self.broadcast_type = "secondary to edges"
 
-    def evaluates_on_edges(self):
+    def evaluates_on_edges(self, dimension):
         return True
 
 
@@ -305,7 +305,7 @@ def __init__(self, child, broadcast_domain, auxiliary_domains, name=None):
         super().__init__(child, broadcast_domain, auxiliary_domains, name)
         self.broadcast_type = "full to edges"
 
-    def evaluates_on_edges(self):
+    def evaluates_on_edges(self, dimension):
         return True
 
 

pybamm/expression_tree/concatenations.py

@@ -294,7 +294,7 @@ def _concatenation_jac(self, children_jacs):
                         a single domain"""
                     )
                 child_slice = next(iter(slices.values()))
-                jacs.append(child_jac[child_slice[i]])
+                jacs.append(pybamm.Index(child_jac, child_slice[i]))
         return SparseStack(*jacs)
 
     def _concatenation_new_copy(self, children):

pybamm/expression_tree/functions.py

@@ -169,9 +169,9 @@ def evaluate(self, t=None, y=None, y_dot=None, inputs=None, known_evals=None):
             ]
             return self._function_evaluate(evaluated_children)
 
-    def evaluates_on_edges(self):
+    def evaluates_on_edges(self, dimension):
         """ See :meth:`pybamm.Symbol.evaluates_on_edges()`. """
-        return any(child.evaluates_on_edges() for child in self.children)
+        return any(child.evaluates_on_edges(dimension) for child in self.children)
 
     def _evaluate_for_shape(self):
         """
@@ -450,5 +450,3 @@ def _function_diff(self, children, idx):
 def arctan(child):
     " Returns hyperbolic tan function of child. "
     return pybamm.simplify_if_constant(Arctan(child), keep_domains=True)
-
-

pybamm/expression_tree/independent_variable.py

@@ -127,7 +127,7 @@ class SpatialVariableEdge(SpatialVariable):
     def __init__(self, name, domain=None, auxiliary_domains=None, coord_sys=None):
         super().__init__(name, domain, auxiliary_domains, coord_sys)
 
-    def evaluates_on_edges(self):
+    def evaluates_on_edges(self, dimension):
         return True
 
 

pybamm/expression_tree/input_parameter.py

@@ -41,7 +41,7 @@ def _evaluate_for_shape(self):
         Returns the scalar 'NaN' to represent the shape of a parameter.
         See :meth:`pybamm.Symbol.evaluate_for_shape()`
         """
-        return np.nan * np.ones_like(self._expected_size)
+        return np.nan * np.ones((self._expected_size, 1))
 
     def _jac(self, variable):
         """ See :meth:`pybamm.Symbol._jac()`. """
@@ -55,7 +55,7 @@ def _base_evaluate(self, t=None, y=None, y_dot=None, inputs=None):
         if not isinstance(inputs, dict):
             # if the special input "shape test" is passed, just return 1
             if inputs == "shape test":
-                return np.ones_like(self._expected_size)
+                return np.ones((self._expected_size, 1))
             raise TypeError("inputs should be a dictionary")
         try:
             input_eval = inputs[self.name]
@@ -64,15 +64,20 @@ def _base_evaluate(self, t=None, y=None, y_dot=None, inputs=None):
             raise KeyError("Input parameter '{}' not found".format(self.name))
 
         if isinstance(input_eval, numbers.Number):
-            input_shape = 1
+            input_size = 1
+            input_ndim = 0
         else:
-            input_shape = input_eval.shape[0]
-        if input_shape == self._expected_size:
-            return input_eval
+            input_size = input_eval.shape[0]
+            input_ndim = len(input_eval.shape)
+        if input_size == self._expected_size:
+            if input_ndim == 1:
+                return input_eval[:, np.newaxis]
+            else:
+                return input_eval
         else:
             raise ValueError(
                 "Input parameter '{}' was given an object of size '{}'".format(
-                    self.name, input_shape
+                    self.name, input_size
                 )
                 + " but was expecting an object of size '{}'.".format(
                     self._expected_size

pybamm/expression_tree/operations/simplify.py

@@ -611,7 +611,9 @@ def _simplify(self, symbol, clear_domains=True):
 
         elif isinstance(symbol, pybamm.UnaryOperator):
             # Reassign domain for gradient and divergence
-            if isinstance(symbol, (pybamm.Gradient, pybamm.Divergence)):
+            if isinstance(
+                symbol, (pybamm.Gradient, pybamm.Divergence, pybamm.Integral)
+            ):
                 new_child = self.simplify(symbol.child, clear_domains=False)
             else:
                 new_child = self.simplify(symbol.child)

pybamm/expression_tree/symbol.py

@@ -467,10 +467,6 @@ def __abs__(self):
             pybamm.AbsoluteValue(self), keep_domains=True
         )
 
-    def __getitem__(self, key):
-        """return a :class:`Index` object"""
-        return pybamm.simplify_if_constant(pybamm.Index(self, key), keep_domains=True)
-
     def diff(self, variable):
         """
         Differentiate a symbol with respect to a variable. For any symbol that can be
@@ -670,16 +666,28 @@ def evaluates_to_number(self):
         result = self.evaluate_ignoring_errors()
 
         if isinstance(result, numbers.Number) or (
-            isinstance(result, np.ndarray) and result.shape == ()
+            isinstance(result, np.ndarray) and np.prod(result.shape) == 1
         ):
             return True
         else:
             return False
 
-    def evaluates_on_edges(self):
+    def evaluates_on_edges(self, dimension):
         """
         Returns True if a symbol evaluates on an edge, i.e. symbol contains a gradient
         operator, but not a divergence operator, and is not an IndefiniteIntegral.
+
+        Parameters
+        ----------
+        dimension : str
+            The dimension (primary, secondary, etc) in which to query evaluation on
+            edges
+
+        Returns
+        -------
+        bool
+            Whether the symbol evaluates on edges (in the finite volume discretisation
+            sense)
         """
         # Default behaviour: return False
         return False