Issue 1126 upwind

CHANGELOG.md

@@ -2,7 +2,7 @@
 
 ## Features
 
-
+-   Added `Upwind` and `Downwind` operators for convection ([#1134](https://github.com/pybamm-team/PyBaMM/pull/1134))
 -   Added Getting Started notebook on solver options and changing the mesh. Also added a notebook detailing the different thermal options, and a notebook explaining the steps that occur behind the scenes in the `Simulation` class ([#1131](https://github.com/pybamm-team/PyBaMM/pull/1131))
 -   Added particle submodel that use a polynomial approximation to the concentration within the electrode particles ([#1130](https://github.com/pybamm-team/PyBaMM/pull/1130))
 -   Added Modulo, Floor and Ceiling operators ([#1121](https://github.com/pybamm-team/PyBaMM/pull/1121))

docs/source/expression_tree/binary_operator.rst

@@ -34,6 +34,9 @@ Binary Operators
 .. autoclass:: pybamm.NotEqualHeaviside
   :members:
 
+.. autoclass:: pybamm.Modulo
+  :members:
+
 .. autoclass:: pybamm.Minimum
   :members:
 

docs/source/expression_tree/unary_operator.rst

@@ -58,6 +58,15 @@ Unary Operators
 .. autoclass:: pybamm.BoundaryGradient
   :members:
 
+.. autoclass:: pybamm.UpwindDownwind
+  :members:
+
+.. autoclass:: pybamm.Upwind
+  :members:
+
+.. autoclass:: pybamm.Downwind
+  :members:
+
 .. autofunction:: pybamm.grad
 
 .. autofunction:: pybamm.div
@@ -79,3 +88,7 @@ Unary Operators
 .. autofunction:: pybamm.boundary_value
 
 .. autofunction:: pybamm.sign
+
+.. autofunction:: pybamm.upwind
+
+.. autofunction:: pybamm.downwind

pybamm/discretisations/discretisation.py

@@ -927,7 +927,11 @@ def _process_symbol(self, symbol):
                 return child_spatial_method.boundary_value_or_flux(
                     symbol, disc_child, self.bcs
                 )
-
+            elif isinstance(symbol, pybamm.UpwindDownwind):
+                direction = symbol.name  # upwind or downwind
+                return spatial_method.upwind_or_downwind(
+                    child, disc_child, self.bcs, direction
+                )
             else:
                 return symbol._unary_new_copy(disc_child)
 

pybamm/expression_tree/binary_operators.py

@@ -745,7 +745,7 @@ def _diff(self, variable):
         # derivative if variable is in the right term (rare, check separately to avoid
         # unecessarily big tree)
         if any(variable.id == x.id for x in right.pre_order()):
-            diff += - pybamm.Floor(left / right) * right.diff(variable)
+            diff += -pybamm.Floor(left / right) * right.diff(variable)
         return diff
 
     def _binary_jac(self, left_jac, right_jac):
@@ -757,7 +757,7 @@ def _binary_jac(self, left_jac, right_jac):
         elif right.evaluates_to_number():
             return left_jac
         elif left.evaluates_to_number():
-            return - right_jac * pybamm.Floor(left / right)
+            return -right_jac * pybamm.Floor(left / right)
         else:
             return left_jac - right_jac * pybamm.Floor(left / right)
 

pybamm/expression_tree/unary_operators.py

@@ -391,7 +391,7 @@ def __init__(self, child):
         if child.evaluates_on_edges("primary") is False:
             raise TypeError(
                 "Cannot take divergence of '{}' since it does not ".format(child)
-                + "evaluates on nodes. Usually, a gradient should be taken before the "
+                + "evaluate on edges. Usually, a gradient should be taken before the "
                 "divergence."
             )
         super().__init__("div", child)
@@ -942,94 +942,152 @@ def __init__(self, child, side):
         super().__init__("boundary flux", child, side)
 
 
+class UpwindDownwind(SpatialOperator):
+    """A node in the expression tree representing an upwinding or downwinding operator.
+    Usually to be used for better stability in convection-dominated equations.
+
+    **Extends:** :class:`SpatialOperator`
+    """
+
+    def __init__(self, name, child):
+        if child.domain == []:
+            raise pybamm.DomainError(
+                "Cannot upwind '{}' since its domain is empty. ".format(child)
+                + "Try broadcasting the object first, e.g.\n\n"
+                "\tpybamm.div(pybamm.PrimaryBroadcast(symbol, 'domain'))"
+            )
+        if child.evaluates_on_edges("primary") is True:
+            raise TypeError(
+                "Cannot upwind '{}' since it does not ".format(child)
+                + "evaluate on nodes."
+            )
+        super().__init__(name, child)
+
+    def evaluates_on_edges(self, dimension):
+        """ See :meth:`pybamm.Symbol.evaluates_on_edges()`. """
+        return True
+
+
+class Upwind(UpwindDownwind):
+    """
+    Upwinding operator. To be used if flow velocity is positive (left to right).
+
+    **Extends:** :class:`UpwindDownwind`
+    """
+
+    def __init__(self, child):
+        super().__init__("upwind", child)
+
+
+class Downwind(UpwindDownwind):
+    """
+    Downwinding operator. To be used if flow velocity is negative (right to left).
+
+    **Extends:** :class:`UpwindDownwind`
+    """
+
+    def __init__(self, child):
+        super().__init__("downwind", child)
+
+
 #
 # Methods to call Gradient, Divergence, Laplacian and Gradient_Squared
 #
 
 
-def grad(expression):
+def grad(symbol):
     """convenience function for creating a :class:`Gradient`
 
     Parameters
     ----------
 
-    expression : :class:`Symbol`
-        the gradient will be performed on this sub-expression
+    symbol : :class:`Symbol`
+        the gradient will be performed on this sub-symbol
 
     Returns
     -------
 
     :class:`Gradient`
-        the gradient of ``expression``
+        the gradient of ``symbol``
     """
     # Gradient of a broadcast is zero
-    if isinstance(expression, pybamm.PrimaryBroadcast):
-        new_child = pybamm.PrimaryBroadcast(0, expression.child.domain)
-        return pybamm.PrimaryBroadcastToEdges(new_child, expression.domain)
+    if isinstance(symbol, pybamm.PrimaryBroadcast):
+        new_child = pybamm.PrimaryBroadcast(0, symbol.child.domain)
+        return pybamm.PrimaryBroadcastToEdges(new_child, symbol.domain)
     else:
-        return Gradient(expression)
+        return Gradient(symbol)
 
 
-def div(expression):
+def div(symbol):
     """convenience function for creating a :class:`Divergence`
 
     Parameters
     ----------
 
-    expression : :class:`Symbol`
-        the divergence will be performed on this sub-expression
+    symbol : :class:`Symbol`
+        the divergence will be performed on this sub-symbol
 
     Returns
     -------
 
     :class:`Divergence`
-        the divergence of ``expression``
+        the divergence of ``symbol``
     """
     # Divergence of a broadcast is zero
-    if isinstance(expression, pybamm.PrimaryBroadcastToEdges):
-        new_child = pybamm.PrimaryBroadcast(0, expression.child.domain)
-        return pybamm.PrimaryBroadcast(new_child, expression.domain)
+    if isinstance(symbol, pybamm.PrimaryBroadcastToEdges):
+        new_child = pybamm.PrimaryBroadcast(0, symbol.child.domain)
+        return pybamm.PrimaryBroadcast(new_child, symbol.domain)
     else:
-        return Divergence(expression)
+        return Divergence(symbol)
 
 
-def laplacian(expression):
+def laplacian(symbol):
     """convenience function for creating a :class:`Laplacian`
 
     Parameters
     ----------
 
-    expression : :class:`Symbol`
-        the laplacian will be performed on this sub-expression
+    symbol : :class:`Symbol`
+        the laplacian will be performed on this sub-symbol
 
     Returns
     -------
 
     :class:`Laplacian`
-        the laplacian of ``expression``
+        the laplacian of ``symbol``
     """
 
-    return Laplacian(expression)
+    return Laplacian(symbol)
 
 
-def grad_squared(expression):
+def grad_squared(symbol):
     """convenience function for creating a :class:`Gradient_Squared`
 
     Parameters
     ----------
 
-    expression : :class:`Symbol`
+    symbol : :class:`Symbol`
         the inner product of the gradient with itself will be performed on this
-        sub-expression
+        sub-symbol
 
     Returns
     -------
 
     :class:`Gradient_Squared`
-        inner product of the gradient of ``expression`` with itself
+        inner product of the gradient of ``symbol`` with itself
     """
 
-    return Gradient_Squared(expression)
+    return Gradient_Squared(symbol)
+
+
+def upwind(symbol):
+    "convenience function for creating a :class:`Upwind`"
+    return Upwind(symbol)
+
+
+def downwind(symbol):
+    "convenience function for creating a :class:`Downwind`"
+    return Downwind(symbol)
 
 
 #

pybamm/plotting/quick_plot.py

@@ -742,7 +742,7 @@ def slider_update(self, t):
                     vmin = ax_min(var)
                     vmax = ax_max(var)
                     cb = self.colorbars[key]
-                    cb.update_bruteforce(
+                    cb.update_normal(
                         cm.ScalarMappable(
                             colors.Normalize(vmin=vmin, vmax=vmax), cmap="coolwarm"
                         )

pybamm/spatial_methods/finite_volume.py

@@ -1384,3 +1384,66 @@ def harmonic_mean(array):
         else:
             raise ValueError("method '{}' not recognised".format(method))
         return out
+
+    def upwind_or_downwind(self, symbol, discretised_symbol, bcs, direction):
+        """
+        Implement an upwinding operator. Currently, this requires the symbol to have
+        a Dirichlet boundary condition on the left side (for upwinding) or right side
+        (for downwinding).
+
+        Parameters
+        ----------
+        symbol : :class:`pybamm.SpatialVariable`
+            The variable to be discretised
+        discretised_gradient : :class:`pybamm.Vector`
+            Contains the discretised gradient of symbol
+        bcs : dict of tuples (:class:`pybamm.Scalar`, str)
+            Dictionary (with keys "left" and "right") of boundary conditions. Each
+            boundary condition consists of a value and a flag indicating its type
+            (e.g. "Dirichlet")
+        direction : str
+            Direction in which to apply the operator (upwind or downwind)
+        """
+        submesh = self.mesh.combine_submeshes(*symbol.domain)
+        n = submesh.npts
+
+        if symbol.id not in bcs:
+            raise pybamm.ModelError(
+                "Boundary conditions must be provided for "
+                "{}ing '{}'".format(direction, symbol)
+            )
+
+        if direction == "upwind":
+            bc, typ = bcs[symbol.id]["left"]
+            if typ != "Dirichlet":
+                raise pybamm.ModelError(
+                    "Dirichlet boundary conditions must be provided for "
+                    "upwinding '{}'".format(symbol)
+                )
+
+            concat_bc = pybamm.NumpyConcatenation(bc, discretised_symbol)
+
+            upwind_mat = vstack(
+                [
+                    csr_matrix(([1], ([0], [0])), shape=(1, n + 1)),
+                    diags([-0.5, 1.5], [0, 1], shape=(n, n + 1)),
+                ]
+            )
+            symbol_out = pybamm.Matrix(upwind_mat) @ concat_bc
+        elif direction == "downwind":
+            bc, typ = bcs[symbol.id]["right"]
+            if typ != "Dirichlet":
+                raise pybamm.ModelError(
+                    "Dirichlet boundary conditions must be provided for "
+                    "downwinding '{}'".format(symbol)
+                )
+
+            concat_bc = pybamm.NumpyConcatenation(discretised_symbol, bc)
+            downwind_mat = vstack(
+                [
+                    diags([1.5, -0.5], [0, 1], shape=(n, n + 1)),
+                    csr_matrix(([1], ([0], [n])), shape=(1, n + 1)),
+                ]
+            )
+            symbol_out = pybamm.Matrix(downwind_mat) @ concat_bc
+        return symbol_out

tests/unit/test_expression_tree/test_operations/test_copy.py

@@ -29,6 +29,7 @@ def test_symbol_new_copy(self):
             pybamm.FunctionParameter("function", {"a": a}),
             pybamm.grad(v_n),
             pybamm.div(pybamm.grad(v_n)),
+            pybamm.upwind(v_n),
             pybamm.IndefiniteIntegral(v_n, x_n),
             pybamm.BackwardIndefiniteIntegral(v_n, x_n),
             pybamm.BoundaryValue(v_n, "right"),