2D FVM

src/pybamm/__init__.py

@@ -43,6 +43,7 @@
 from .expression_tree.independent_variable import *
 from .expression_tree.independent_variable import t
 from .expression_tree.vector import Vector
+from .expression_tree.vector_field import VectorField
 from .expression_tree.state_vector import StateVectorBase, StateVector, StateVectorDot
 
 from .expression_tree.exceptions import *
@@ -139,6 +140,10 @@
     SpectralVolume1DSubMesh,
     SymbolicUniform1DSubMesh,
 )
+from .meshes.two_dimensional_submeshes import (
+    SubMesh2D,
+    Uniform2DSubMesh,
+)
 from .meshes.scikit_fem_submeshes import (
     ScikitSubMesh2D,
     ScikitUniform2DSubMesh,
@@ -154,6 +159,7 @@
 from .spatial_methods.spatial_method import SpatialMethod
 from .spatial_methods.zero_dimensional_method import ZeroDimensionalSpatialMethod
 from .spatial_methods.finite_volume import FiniteVolume
+from .spatial_methods.finite_volume_2d import FiniteVolume2D
 from .spatial_methods.spectral_volume import SpectralVolume
 from .spatial_methods.scikit_finite_element import ScikitFiniteElement
 

src/pybamm/discretisations/discretisation.py

@@ -270,27 +270,65 @@ def set_variable_slices(self, variables):
         end = 0
         lower_bounds = []
         upper_bounds = []
+
+        # Make sure we hit concatenation variables first.
+        # concat_variables = [var for var in variables if isinstance(var, pybamm.ConcatenationVariable)]
+        # not_concat_variables = [var for var in variables if not isinstance(var, pybamm.ConcatenationVariable)]
+        # variables = concat_variables + not_concat_variables
+
         # Iterate through unpacked variables, adding appropriate slices to y_slices
         for variable in variables:
+            if variable in y_slices:
+                continue
             # Add up the size of all the domains in variable.domain
             if isinstance(variable, pybamm.ConcatenationVariable):
-                start_ = start
                 spatial_method = self.spatial_methods[variable.domain[0]]
+                dimension = spatial_method.mesh[variable.domain[0]].dimension
+                start_ = start
                 children = variable.children
                 meshes = OrderedDict()
+                lr_points = OrderedDict()
+                tb_points = OrderedDict()
                 for child in children:
                     meshes[child] = [spatial_method.mesh[dom] for dom in child.domain]
+                    if dimension == 2:
+                        lr_points[child] = sum(
+                            spatial_method.mesh[dom].npts_lr for dom in child.domain
+                        )
+                        tb_points[child] = sum(
+                            spatial_method.mesh[dom].npts_tb for dom in child.domain
+                        )
                 sec_points = spatial_method._get_auxiliary_domain_repeats(
                     variable.domains
                 )
                 for _ in range(sec_points):
+                    start_this_child = start_
                     for child, mesh in meshes.items():
                         for domain_mesh in mesh:
                             end += domain_mesh.npts_for_broadcast_to_nodes
                         # Add to slices
-                        y_slices[child].append(slice(start_, end))
-                        y_slices_explicit[child].append(slice(start_, end))
-                        # Increment start_
+                        if dimension == 2:
+                            other_children = set(meshes.keys()) - {child}
+                            num_pts_to_skip = sum(
+                                lr_points[other_child] for other_child in other_children
+                            )
+                            for row in range(tb_points[child]):
+                                start_this_row = (
+                                    start_this_child
+                                    + (lr_points[child] + num_pts_to_skip) * row
+                                )
+                                end_this_row = start_this_row + lr_points[child]
+                                y_slices[child].append(
+                                    slice(start_this_row, end_this_row)
+                                )
+                                y_slices_explicit[child].append(
+                                    slice(start_this_row, end_this_row)
+                                )
+                            start_this_child += lr_points[child]
+                        else:
+                            y_slices[child].append(slice(start_, end))
+                            y_slices_explicit[child].append(slice(start_, end))
+                            # Increment start_
                         start_ = end
             else:
                 end += self._get_variable_size(variable)
@@ -808,6 +846,22 @@ def _process_symbol(self, symbol):
         if isinstance(symbol, pybamm.BinaryOperator):
             # Pre-process children
             left, right = symbol.children
+            # Catch case where diffusion is a scalar and turn it into an identity matrix vector field.
+            if len(symbol.domain) != 0:
+                spatial_method = self.spatial_methods[symbol.domain[0]]
+            else:
+                spatial_method = None
+            if isinstance(spatial_method, pybamm.FiniteVolume2D):
+                if isinstance(left, pybamm.Scalar) and (
+                    isinstance(right, pybamm.VectorField)
+                    or isinstance(right, pybamm.Gradient)
+                ):
+                    left = pybamm.VectorField(left, left)
+                elif isinstance(right, pybamm.Scalar) and (
+                    isinstance(left, pybamm.VectorField)
+                    or isinstance(left, pybamm.Gradient)
+                ):
+                    right = pybamm.VectorField(right, right)
             disc_left = self.process_symbol(left)
             disc_right = self.process_symbol(right)
             if symbol.domain == []:
@@ -875,7 +929,10 @@ def _process_symbol(self, symbol):
                     symbol.integration_variable[0].domain[0]
                 ]
                 out = integral_spatial_method.integral(
-                    child, disc_child, symbol._integration_dimension
+                    child,
+                    disc_child,
+                    symbol._integration_dimension,
+                    symbol.integration_variable,
                 )
                 out.copy_domains(symbol)
                 return out
@@ -915,6 +972,14 @@ def _process_symbol(self, symbol):
                 return child_spatial_method.evaluate_at(
                     symbol, disc_child, symbol.position
                 )
+            elif isinstance(symbol, pybamm.UpwindDownwind2D):
+                return spatial_method.upwind_or_downwind(
+                    child,
+                    disc_child,
+                    self.bcs,
+                    symbol.lr_direction,
+                    symbol.tb_direction,
+                )
             elif isinstance(symbol, pybamm.UpwindDownwind):
                 direction = symbol.name  # upwind or downwind
                 return spatial_method.upwind_or_downwind(
@@ -923,6 +988,16 @@ def _process_symbol(self, symbol):
             elif isinstance(symbol, pybamm.NotConstant):
                 # After discretisation, we can make the symbol constant
                 return disc_child
+            elif isinstance(symbol, pybamm.Magnitude):
+                if not isinstance(disc_child, pybamm.VectorField):
+                    raise ValueError("Magnitude can only be applied to a vector field")
+                direction = symbol.direction
+                if direction == "lr":
+                    return disc_child.lr_field
+                elif direction == "tb":
+                    return disc_child.tb_field
+                else:
+                    raise ValueError("Invalid direction")
             else:
                 return symbol.create_copy(new_children=[disc_child])
 
@@ -994,6 +1069,11 @@ def _process_symbol(self, symbol):
             new_symbol = self.process_symbol(symbol.children[0])
             return new_symbol
 
+        elif isinstance(symbol, pybamm.VectorField):
+            left_symbol = self.process_symbol(symbol.lr_field)
+            right_symbol = self.process_symbol(symbol.tb_field)
+            return symbol.create_copy(new_children=[left_symbol, right_symbol])
+
         else:
             # Backup option: return the object
             return symbol

src/pybamm/expression_tree/__init__.py

@@ -2,4 +2,4 @@
            'concatenations', 'exceptions', 'functions', 'independent_variable',
            'input_parameter', 'interpolant', 'matrix', 'operations',
            'parameter', 'printing', 'scalar', 'state_vector', 'symbol',
-           'unary_operators', 'variable', 'vector', 'discrete_time_sum' ]
+           'unary_operators', 'variable', 'vector', 'discrete_time_sum', 'vector_field', 'magnitude' ]

src/pybamm/expression_tree/concatenations.py

@@ -451,15 +451,15 @@ class SparseStack(Concatenation):
         The equations to concatenate
     """
 
-    def __init__(self, *children):
+    def __init__(self, *children, name="sparse_stack"):
         children = list(children)
         if not any(issparse(child.evaluate_for_shape()) for child in children):
             concatenation_function = np.vstack
         else:
             concatenation_function = vstack
         super().__init__(
             *children,
-            name="sparse_stack",
+            name=name,
             check_domain=False,
             concat_fun=concatenation_function,
         )

src/pybamm/expression_tree/independent_variable.py

@@ -142,12 +142,14 @@ def __init__(
         auxiliary_domains: AuxiliaryDomainType = None,
         domains: DomainsType = None,
         coord_sys=None,
+        direction=None,
     ) -> None:
         self.coord_sys = coord_sys
         super().__init__(
             name, domain=domain, auxiliary_domains=auxiliary_domains, domains=domains
         )
         domain = self.domain
+        self.direction = direction
 
         if domain == []:
             raise ValueError("domain must be provided")

src/pybamm/expression_tree/unary_operators.py

@@ -904,8 +904,17 @@ def __init__(self, child, region="entire"):
             name += "negative tab"
         elif region == "positive tab":
             name += "positive tab"
+        elif region == "top":
+            name += "top"
+        elif region == "bottom":
+            name += "bottom"
+        elif region == "left":
+            name += "left"
+        elif region == "right":
+            name += "right"
         self.region = region
         super().__init__(name, child, domains)
+        self.domains = {}
 
     def set_id(self):
         """See :meth:`pybamm.Symbol.set_id()`"""
@@ -1187,6 +1196,10 @@ class UpwindDownwind(SpatialOperator):
     """
 
     def __init__(self, name, child):
+        self._perform_checks(child)
+        super().__init__(name, child)
+
+    def _perform_checks(self, child):
         if child.domain == []:
             raise pybamm.DomainError(
                 f"Cannot upwind '{child}' since its domain is empty. "
@@ -1197,13 +1210,42 @@ def __init__(self, name, child):
             raise TypeError(
                 f"Cannot upwind '{child}' since it does not " + "evaluate on nodes."
             )
-        super().__init__(name, child)
 
     def _evaluates_on_edges(self, dimension: str) -> bool:
         """See :meth:`pybamm.Symbol._evaluates_on_edges()`."""
         return True
 
 
+class UpwindDownwind2D(UpwindDownwind):
+    """
+    A node in the expression tree representing an upwinding or downwinding operator.
+    Usually to be used for better stability in convection-dominated equations.
+    """
+
+    def __init__(self, child, lr_direction, tb_direction):
+        super().__init__("upwind_downwind_2d", child)
+        self.lr_direction = lr_direction
+        self.tb_direction = tb_direction
+
+    def _unary_new_copy(self, child, perform_simplifications=True):
+        """See :meth:`UnaryOperator._unary_new_copy()`."""
+        return self.__class__(child, self.lr_direction, self.tb_direction)
+
+
+class Magnitude(UnaryOperator):
+    """
+    A node in the expression tree representing the magnitude of a vector field.
+    """
+
+    def __init__(self, child, direction):
+        super().__init__("magnitude" + f"({direction})", child)
+        self.direction = direction
+
+    def _unary_new_copy(self, child, perform_simplifications=True):
+        """See :meth:`UnaryOperator._unary_new_copy()`."""
+        return self.__class__(child, self.direction)
+
+
 class Upwind(UpwindDownwind):
     """
     Upwinding operator. To be used if flow velocity is positive (left to right).

src/pybamm/expression_tree/vector_field.py

@@ -0,0 +1,32 @@
+import pybamm
+from typing import Optional
+
+
+class VectorField(pybamm.Symbol):
+    """
+    A node in the expression tree representing a vector field.
+    """
+
+    def __init__(self, lr_field, tb_field):
+        children = [lr_field, tb_field]
+        if lr_field.domain != tb_field.domain:
+            raise ValueError("lr_field and tb_field must have the same domain")
+        super().__init__(name="vector_field", children=children, domain=lr_field.domain)
+        self.lr_field = lr_field
+        self.tb_field = tb_field
+
+    def create_copy(self, new_children: Optional[list[pybamm.Symbol]] = None):
+        if new_children is None:
+            new_children = [self.lr_field, self.tb_field]
+        return VectorField(*new_children)
+
+    def _evaluate_for_shape(self):
+        return self.lr_field.evaluate_for_shape()
+
+    def evaluates_on_edges(self, dimension: str) -> bool:
+        left_evaluates_on_edges = self.lr_field.evaluates_on_edges(dimension)
+        right_evaluates_on_edges = self.tb_field.evaluates_on_edges(dimension)
+        if left_evaluates_on_edges == right_evaluates_on_edges:
+            return left_evaluates_on_edges
+        else:
+            raise ValueError("lr_field and tb_field must have the same domain")