Issue 580 processed variable

examples/scripts/compare_lead_acid_3D.py

@@ -83,6 +83,7 @@
 for i, model in enumerate(models):
     solution = model.default_solver.solve(model, t_eval)
     solutions[i] = solution
+    pybamm.post_process_variables(model.variables, solution.t, solution.y, mesh=mesh)
 
 # plot
 output_variables = [

pybamm/discretisations/discretisation.py

@@ -472,8 +472,6 @@ def process_dict(self, var_eqn_dict):
 
             new_var_eqn_dict[eqn_key] = self.process_symbol(eqn)
 
-            new_var_eqn_dict[eqn_key].test_shape()
-
         return new_var_eqn_dict
 
     def process_symbol(self, symbol):
@@ -496,6 +494,7 @@ def process_symbol(self, symbol):
         except KeyError:
             discretised_symbol = self._process_symbol(symbol)
             self._discretised_symbols[symbol.id] = discretised_symbol
+            discretised_symbol.test_shape()
             return discretised_symbol
 
     def _process_symbol(self, symbol):

pybamm/expression_tree/binary_operators.py

@@ -82,11 +82,12 @@ def __init__(self, name, left, right):
         # right child.
         if isinstance(self, (pybamm.Outer, pybamm.Kron)):
             domain = right.domain
+            auxiliary_domains = {"secondary": left.domain}
         else:
             domain = self.get_children_domains(left.domain, right.domain)
-        auxiliary_domains = self.get_children_auxiliary_domains(
-            left.auxiliary_domains, right.auxiliary_domains
-        )
+            auxiliary_domains = self.get_children_auxiliary_domains(
+                left.auxiliary_domains, right.auxiliary_domains
+            )
         super().__init__(
             name,
             children=[left, right],

pybamm/expression_tree/functions.py

@@ -31,6 +31,7 @@ def __init__(self, function, *children):
             name = "function ({})".format(function.__class__)
         children_list = list(children)
         domain = self.get_children_domains(children_list)
+        auxiliary_domains = self.get_children_auxiliary_domains(children)
 
         self.function = function
 
@@ -41,7 +42,12 @@ def __init__(self, function, *children):
         else:
             self.takes_no_params = len(signature(function).parameters) == 0
 
-        super().__init__(name, children=children_list, domain=domain)
+        super().__init__(
+            name,
+            children=children_list,
+            domain=domain,
+            auxiliary_domains=auxiliary_domains,
+        )
 
     def get_children_domains(self, children_list):
         """Obtains the unique domain of the children. If the
@@ -63,6 +69,27 @@ def get_children_domains(self, children_list):
 
         return domain
 
+    def get_children_auxiliary_domains(self, children):
+        "Combine auxiliary domains from children, at all levels"
+        aux_domains = {}
+        for child in children:
+            for level in child.auxiliary_domains.keys():
+                if (
+                    not hasattr(aux_domains, level)
+                    or aux_domains[level] == []
+                    or child.auxiliary_domains[level] == aux_domains[level]
+                ):
+                    aux_domains[level] = child.auxiliary_domains[level]
+                else:
+                    raise pybamm.DomainError(
+                        """children must have same or empty auxiliary domains,
+                        not {!s} and {!s}""".format(
+                            aux_domains[level], child.auxiliary_domains[level]
+                        )
+                    )
+
+        return aux_domains
+
     def diff(self, variable):
         """ See :meth:`pybamm.Symbol.diff()`. """
         if variable.id == self.id:

pybamm/expression_tree/parameter.py

@@ -61,7 +61,10 @@ def __init__(self, name, *children, diff_variable=None):
         self.diff_variable = diff_variable
         children_list = list(children)
         domain = self.get_children_domains(children_list)
-        super().__init__(name, children=children, domain=domain)
+        auxiliary_domains = self.get_children_auxiliary_domains(children)
+        super().__init__(
+            name, children=children, domain=domain, auxiliary_domains=auxiliary_domains
+        )
 
     def set_id(self):
         """See :meth:`pybamm.Symbol.set_id` """
@@ -90,6 +93,27 @@ def get_children_domains(self, children_list):
 
         return domain
 
+    def get_children_auxiliary_domains(self, children):
+        "Combine auxiliary domains from children, at all levels"
+        aux_domains = {}
+        for child in children:
+            for level in child.auxiliary_domains.keys():
+                if (
+                    not hasattr(aux_domains, level)
+                    or aux_domains[level] == []
+                    or child.auxiliary_domains[level] == aux_domains[level]
+                ):
+                    aux_domains[level] = child.auxiliary_domains[level]
+                else:
+                    raise pybamm.DomainError(
+                        """children must have same or empty auxiliary domains,
+                        not {!s} and {!s}""".format(
+                            aux_domains[level], child.auxiliary_domains[level]
+                        )
+                    )
+
+        return aux_domains
+
     def diff(self, variable):
         """ See :meth:`pybamm.Symbol.diff()`. """
         # return a new FunctionParameter, that knows it will need to be differentiated

pybamm/expression_tree/unary_operators.py

@@ -172,10 +172,12 @@ def __init__(self, child, index, name=None):
         else:
             raise TypeError("index must be integer or slice")
 
-        if self.slice in (slice(0, 1), slice(-1, None)):
-            pass
-        elif self.slice.stop > child.size:
-            raise ValueError("slice size exceeds child size")
+        # Perform some additional checks if debug_mode is True (child.size is slow)
+        if pybamm.settings.debug_mode is True:
+            if self.slice in (slice(0, 1), slice(-1, None)):
+                pass
+            elif self.slice.stop > child.size:
+                raise ValueError("slice size exceeds child size")
 
         super().__init__(name, child)
 

pybamm/processed_variable.py

@@ -200,37 +200,80 @@ def initialise_2D(self):
 
     def initialise_3D(self):
         """
-        Initialise a 3D object that depends on x and r.
+        Initialise a 3D object that depends on x and r, or x and z.
         Needs to be generalised to deal with other domains
         """
-        if self.domain in [["negative particle"], ["negative electrode"]]:
+        # Dealt with weird particle/electrode case
+        if self.domain in [
+            ["negative electrode"],
+            ["positive electrode"],
+        ] and self.auxiliary_domains["secondary"] in [
+            ["negative particle"],
+            ["positive particle"],
+        ]:
+            # Switch domain and auxiliary domains and set order to "F"
+            dom = self.domain
+            self.domain = self.auxiliary_domains["secondary"]
+            self.auxiliary_domains["secondary"] = dom
+            order = "F"
+        else:
+            order = "C"
+
+        # Process x-r or x-z
+        if self.domain == ["negative particle"] and self.auxiliary_domains[
+            "secondary"
+        ] == ["negative electrode"]:
             x_sol = self.mesh["negative electrode"][0].nodes
             r_nodes = self.mesh["negative particle"][0].nodes
             r_edges = self.mesh["negative particle"][0].edges
-        elif self.domain in [["positive particle"], ["positive electrode"]]:
+            set_up_r = True
+        elif self.domain == ["positive particle"] and self.auxiliary_domains[
+            "secondary"
+        ] == ["positive electrode"]:
             x_sol = self.mesh["positive electrode"][0].nodes
             r_nodes = self.mesh["positive particle"][0].nodes
             r_edges = self.mesh["positive particle"][0].edges
+            set_up_r = True
+        elif self.domain[0] in [
+            "negative electrode",
+            "separator",
+            "positive electrode",
+        ] and self.auxiliary_domains["secondary"] == ["current collector"]:
+            x_nodes = self.mesh.combine_submeshes(*self.domain)[0].nodes
+            x_edges = self.mesh.combine_submeshes(*self.domain)[0].edges
+            z_sol = self.mesh["current collector"][0].nodes
+            r_sol = None
+            self.first_dimension = "x"
+            self.second_dimension = "z"
+
+            if self.base_eval.size // len(z_sol) == len(x_nodes):
+                x_sol = x_nodes
+            elif self.base_eval.size // len(z_sol) == len(x_edges):
+                x_sol = x_edges
+            first_dim_nodes = x_sol
+            second_dim_nodes = z_sol
+            set_up_r = False
         else:
             raise pybamm.DomainError(
-                """ Can only create 3D objects on electrodes and particles. Current
-                collector not yet implemented"""
+                """ Cannot process 3D object with domain '{}'
+                and auxiliary_domains '{}'""".format(
+                    self.domain, self.auxiliary_domains
+                )
             )
-        len_x = len(x_sol)
-        len_r = self.base_eval.shape[0] // len_x
-        if self.domain in [["negative particle"], ["positive particle"]]:
+        if set_up_r:
+            z_sol = None
             self.first_dimension = "x"
             self.second_dimension = "r"
-            first_dim_size = len_x
-            second_dim_size = len_r
-            transpose = False
-        else:
-            self.first_dimension = "r"
-            self.second_dimension = "x"
-            first_dim_size = len_r
-            second_dim_size = len_x
-            transpose = True
-        entries = np.empty((len_x, len_r, len(self.t_sol)))
+            if self.base_eval.size // len(x_sol) == len(r_nodes):
+                r_sol = r_nodes
+            elif self.base_eval.size // len(x_sol) == len(r_edges):
+                r_sol = r_edges
+            first_dim_nodes = x_sol
+            second_dim_nodes = r_sol
+
+        first_dim_size = len(first_dim_nodes)
+        second_dim_size = len(second_dim_nodes)
+        entries = np.empty((first_dim_size, second_dim_size, len(self.t_sol)))
 
         # Evaluate the base_variable index-by-index
         for idx in range(len(self.t_sol)):
@@ -240,36 +283,29 @@ def initialise_3D(self):
                 eval_and_known_evals = self.base_variable.evaluate(
                     t, u, self.known_evals[t]
                 )
-                temporary = np.reshape(
-                    eval_and_known_evals[0], [first_dim_size, second_dim_size]
+                entries[:, :, idx] = np.reshape(
+                    eval_and_known_evals[0],
+                    [first_dim_size, second_dim_size],
+                    order=order,
                 )
                 self.known_evals[t] = eval_and_known_evals[1]
             else:
-                temporary = np.reshape(
-                    self.base_variable.evaluate(t, u), [first_dim_size, second_dim_size]
+                entries[:, :, idx] = np.reshape(
+                    self.base_variable.evaluate(t, u),
+                    [first_dim_size, second_dim_size],
+                    order=order,
                 )
-            if transpose is True:
-                entries[:, :, idx] = np.transpose(temporary)
-            else:
-                entries[:, :, idx] = temporary
-
-        # Assess whether on nodes or edges
-        if entries.shape[1] == len(r_nodes):
-            r_sol = r_nodes
-        elif entries.shape[1] == len(r_edges):
-            r_sol = r_edges
-        else:
-            raise ValueError("3D variable shape does not match domain shape")
 
         # assign attributes for reference
         self.entries = entries
         self.dimensions = 3
         self.x_sol = x_sol
         self.r_sol = r_sol
+        self.z_sol = z_sol
 
         # set up interpolation
         self._interpolation_function = interp.RegularGridInterpolator(
-            (x_sol, r_sol, self.t_sol),
+            (first_dim_nodes, second_dim_nodes, self.t_sol),
             entries,
             method=self.interp_kind,
             fill_value=np.nan,
@@ -324,32 +360,22 @@ def __call__(self, t, x=None, r=None, y=None, z=None):
 
     def call_2D(self, t, x, r, z):
         "Evaluate a 2D variable"
-        if self.spatial_var_name == "r":
-            if r is not None:
-                return self._interpolation_function(t, r)
-            else:
-                raise ValueError("r cannot be None for microscale variable")
-        elif self.spatial_var_name == "x":
-            if x is not None:
-                return self._interpolation_function(t, x)
-            else:
-                raise ValueError("x cannot be None for macroscale variable")
+        spatial_var = eval(self.spatial_var_name)
+        if spatial_var is not None:
+            return self._interpolation_function(t, spatial_var)
         else:
-            if z is not None:
-                return self._interpolation_function(t, z)
-            else:
-                raise ValueError("z cannot be None for macroscale variable")
+            raise ValueError("input {} cannot be None".format(self.spatial_var_name))
 
     def call_3D(self, t, x, r, y, z):
         "Evaluate a 3D variable"
-        if (self.first_dimension == "x" and self.second_dimension == "r") or (
-            self.first_dimension == "r" and self.second_dimension == "x"
-        ):
-            first_dim = x
-            second_dim = r
-        elif self.first_dimension == "y" and self.second_dimension == "z":
-            first_dim = y
-            second_dim = z
+        first_dim = eval(self.first_dimension)
+        second_dim = eval(self.second_dimension)
+        if first_dim is None or second_dim is None:
+            raise ValueError(
+                "inputs {} and {} cannot be None".format(
+                    self.first_dimension, self.second_dimension
+                )
+            )
         if isinstance(first_dim, np.ndarray):
             if isinstance(second_dim, np.ndarray) and isinstance(t, np.ndarray):
                 first_dim = first_dim[:, np.newaxis, np.newaxis]

pybamm/spatial_methods/spatial_method.py

@@ -84,11 +84,11 @@ def broadcast(self, symbol, domain, auxiliary_domains, broadcast_type):
             out = pybamm.Outer(
                 symbol, pybamm.Vector(np.ones(primary_pts_for_broadcast), domain=domain)
             )
-            out.auxiliary_domains = auxiliary_domains
 
         elif broadcast_type == "full":
             out = symbol * pybamm.Vector(np.ones(full_pts_for_broadcast), domain=domain)
 
+        out.auxiliary_domains = auxiliary_domains
         return out
 
     def gradient(self, symbol, discretised_symbol, boundary_conditions):

tests/unit/test_expression_tree/test_unary_operators.py

@@ -139,8 +139,11 @@ def test_index(self):
         # errors
         with self.assertRaisesRegex(TypeError, "index must be integer or slice"):
             pybamm.Index(vec, 0.0)
+        debug_mode = pybamm.settings.debug_mode
+        pybamm.settings.debug_mode = True
         with self.assertRaisesRegex(ValueError, "slice size exceeds child size"):
             pybamm.Index(vec, 5)
+        pybamm.settings.debug_mode = debug_mode
 
     def test_diff(self):
         a = pybamm.StateVector(slice(0, 1))