#704 removed the option in other functions

Author: Scottmar93

pybamm/discretisations/discretisation.py

@@ -707,9 +707,7 @@ def _process_symbol(self, symbol):
                     mesh = self.mesh[symbol.children[0].domain[0]][0]
                     if isinstance(mesh, pybamm.SubMesh1D):
                         symbol.side = mesh.tabs[symbol.side]
-                return child_spatial_method.boundary_value_or_flux(
-                    symbol, disc_child, symbol.extrapolation
-                )
+                return child_spatial_method.boundary_value_or_flux(symbol, disc_child)
 
             else:
                 return symbol._unary_new_copy(disc_child)

pybamm/expression_tree/unary_operators.py

@@ -738,15 +738,12 @@ class BoundaryValue(BoundaryOperator):
         The variable whose boundary value to take
     side : str
         Which side to take the boundary value on ("left" or "right")
-    extrapolation : str
-        What extrapolation method to use for this variable ("linear" or "quadratic")
 
     **Extends:** :class:`BoundaryOperator`
     """
 
-    def __init__(self, child, side, extrapolation="quadratic"):
+    def __init__(self, child, side):
         super().__init__("boundary value", child, side)
-        self.extrapolation = extrapolation
 
 
 class BoundaryGradient(BoundaryOperator):
@@ -758,15 +755,12 @@ class BoundaryGradient(BoundaryOperator):
         The variable whose boundary flux to take
     side : str
         Which side to take the boundary flux on ("left" or "right")
-    extrapolation : str
-        What extrapolation method to use on this variable ("linear" or "quadratic")
 
     **Extends:** :class:`BoundaryOperator`
     """
 
-    def __init__(self, child, side, extrapolation="linear"):
+    def __init__(self, child, side):
         super().__init__("boundary flux", child, side)
-        self.extrapolation = extrapolation
 
 
 #
@@ -856,7 +850,7 @@ def grad_squared(expression):
 #
 
 
-def surf(symbol, set_domain=False, extrapolation="quadratic"):
+def surf(symbol, set_domain=False):
     """convenience function for creating a right :class:`BoundaryValue`, usually in the
     spherical geometry
 
@@ -865,8 +859,6 @@ def surf(symbol, set_domain=False, extrapolation="quadratic"):
 
     symbol : :class:`pybamm.Symbol`
         the surface value of this symbol will be returned
-    extrapolation : str
-        The type of extrapolation method to use ("linear" or "quadratic")
 
     Returns
     -------
@@ -876,10 +868,10 @@ def surf(symbol, set_domain=False, extrapolation="quadratic"):
     if symbol.domain in [["negative electrode"], ["positive electrode"]] and isinstance(
         symbol, pybamm.PrimaryBroadcast
     ):
-        child_surf = boundary_value(symbol.orphans[0], "right", extrapolation)
+        child_surf = boundary_value(symbol.orphans[0], "right")
         out = pybamm.PrimaryBroadcast(child_surf, symbol.domain)
     else:
-        out = boundary_value(symbol, "right", extrapolation)
+        out = boundary_value(symbol, "right")
         if set_domain:
             if symbol.domain == ["negative particle"]:
                 out.domain = ["negative electrode"]
@@ -1023,7 +1015,7 @@ def yz_average(symbol):
         return Integral(symbol, [y, z]) / (l_y * l_z)
 
 
-def boundary_value(symbol, side, extrapolation="quadratic"):
+def boundary_value(symbol, side):
     """convenience function for creating a :class:`pybamm.BoundaryValue`
 
     Parameters
@@ -1032,8 +1024,6 @@ def boundary_value(symbol, side, extrapolation="quadratic"):
         The symbol whose boundary value to take
     side : str
         Which side to take the boundary value on ("left" or "right")
-    extrapolation : str
-        The type of extrapolation method to use ("linear" or "quadratic")
 
     Returns
     -------
@@ -1050,7 +1040,7 @@ def boundary_value(symbol, side, extrapolation="quadratic"):
         return symbol.orphans[0]
     # Otherwise, calculate boundary value
     else:
-        return BoundaryValue(symbol, side, extrapolation)
+        return BoundaryValue(symbol, side)
 
 
 def r_average(symbol):

pybamm/spatial_methods/finite_volume.py

@@ -34,6 +34,10 @@ class FiniteVolume(pybamm.SpatialMethod):
 
     def __init__(self, mesh):
         super().__init__(mesh)
+
+        # there is no way to set this at the moment
+        self.extrapolation = "linear"
+
         # add npts_for_broadcast to mesh domains for this particular discretisation
         for dom in mesh.keys():
             for i in range(len(mesh[dom])):
@@ -572,7 +576,7 @@ def add_ghost_nodes(self, symbol, discretised_symbol, bcs):
 
         return pybamm.Matrix(matrix) @ discretised_symbol + bcs_vector
 
-    def boundary_value_or_flux(self, symbol, discretised_child, extrapolation="linear"):
+    def boundary_value_or_flux(self, symbol, discretised_child):
         """
         Uses linear extrapolation to get the boundary value or flux of a variable in the
         Finite Volume Method.
@@ -592,13 +596,17 @@ def boundary_value_or_flux(self, symbol, discretised_child, extrapolation="linea
             nodes = submesh_list[0].nodes
             edges = submesh_list[0].edges
 
-            # assumes uniform!
+            dx0 = nodes[0] - edges[0]
+            dx1 = submesh_list[0].d_nodes[0]
+            dx2 = submesh_list[0].d_nodes[1]
+
+            dxN = edges[-1] - nodes[-1]
+            dxNm1 = submesh_list[0].d_nodes[-1]
+            dxNm2 = submesh_list[0].d_nodes[-2]
+
             if symbol.side == "left":
-                dx0 = nodes[0] - edges[0]
-                dx1 = submesh_list[0].d_nodes[0]
-                dx2 = submesh_list[0].d_nodes[1]
 
-                if extrapolation == "linear":
+                if self.extrapolation == "linear":
                     # to find value at x* use formula:
                     # f(x*) = f_1 - (dx0 / dx1) (f_2 - f_1)
                     # where
@@ -609,7 +617,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, extrapolation="linea
                         ([1 + (dx0 / dx1), -(dx0 / dx1)], ([0, 0], [0, 1])),
                         shape=(1, prim_pts),
                     )
-                elif extrapolation == "quadratic":
+                elif self.extrapolation == "quadratic":
                     # to find value at x* use formula:
                     # see mathematica notebook at:
                     # https://github.com/Scottmar93/extrapolation-coefficents/tree/master
@@ -626,11 +634,8 @@ def boundary_value_or_flux(self, symbol, discretised_child, extrapolation="linea
                     raise NotImplementedError
 
             elif symbol.side == "right":
-                dxN = edges[-1] - nodes[-1]
-                dxNm1 = submesh_list[0].d_nodes[-1]
-                dxNm2 = submesh_list[0].d_nodes[-2]
 
-                if extrapolation == "linear":
+                if self.extrapolation == "linear":
                     # to find value at x* use formula:
                     # f(x*) = f_N - (dxN / dxNm1) (f_N - f_Nm1)
                     # where
@@ -644,7 +649,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, extrapolation="linea
                         ),
                         shape=(1, prim_pts),
                     )
-                elif extrapolation == "quadratic":
+                elif self.extrapolation == "quadratic":
                     # to find value at x* use formula:
                     # see mathematica notebook at:
                     # https://github.com/Scottmar93/extrapolation-coefficents/tree/master
@@ -670,23 +675,46 @@ def boundary_value_or_flux(self, symbol, discretised_child, extrapolation="linea
 
         elif isinstance(symbol, pybamm.BoundaryGradient):
             if symbol.side == "left":
-                dx = submesh_list[0].d_nodes[0]
+                # use formula:
+                # f'(x*) = (f_2 - f_1) / dx1
+                # where dx1 = x_2 - x_1
 
-                if extrapolation == "linear":
-                    sub_matrix = (1 / dx) * csr_matrix(
+                if self.extrapolation == "linear":
+                    sub_matrix = (1 / dx1) * csr_matrix(
                         ([-1, 1], ([0, 0], [0, 1])), shape=(1, prim_pts)
                     )
+                elif self.extrapolation == "quadratic":
+
+                    a = -(2 * dx0 + 2 * dx1 + dx2) / (dx1 ** 2 + dx1 * dx2)
+                    b = (2 * dx0 + dx1 + dx2) / (dx1 * dx2)
+                    c = -(2 * dx0 + dx1) / (dx1 * dx2 + dx2 ** 2)
+
+                    sub_matrix = csr_matrix(
+                        ([a, b, c], ([0, 0, 0], [0, 1, 2])), shape=(1, prim_pts)
+                    )
                 else:
                     raise NotImplementedError
 
             elif symbol.side == "right":
                 dx = submesh_list[0].d_nodes[-1]
 
-                if extrapolation == "linear":
+                if self.extrapolation == "linear":
                     sub_matrix = (1 / dx) * csr_matrix(
                         ([-1, 1], ([0, 0], [prim_pts - 2, prim_pts - 1])),
                         shape=(1, prim_pts),
                     )
+                elif self.extrapolation == "quadratic":
+                    a = (2 * dxN + 2 * dxNm1 + dxNm2) / (dxNm1 ** 2 + dxNm1 * dxNm2)
+                    b = -(2 * dxN + dxNm1 + dxNm2) / (dxNm1 * dxNm2)
+                    c = (2 * dxN + dxNm1) / (dxNm1 * dxNm2 + dxNm2 ** 2)
+
+                    sub_matrix = csr_matrix(
+                        (
+                            [c, b, a],
+                            ([0, 0, 0], [prim_pts - 3, prim_pts - 2, prim_pts - 1]),
+                        ),
+                        shape=(1, prim_pts),
+                    )
                 else:
                     raise NotImplementedError
 

pybamm/spatial_methods/scikit_finite_element.py

@@ -322,7 +322,7 @@ def integral_form(v, dv, w):
 
         return pybamm.Matrix(integration_vector[np.newaxis, :])
 
-    def boundary_value_or_flux(self, symbol, discretised_child, extrapolation="linear"):
+    def boundary_value_or_flux(self, symbol, discretised_child):
         """
         Returns the average value of the symbol over the negative tab ("negative tab")
         or the positive tab ("positive tab") in the Finite Element Method.

pybamm/spatial_methods/spatial_method.py

@@ -278,7 +278,7 @@ def internal_neumann_condition(
 
         raise NotImplementedError
 
-    def boundary_value_or_flux(self, symbol, discretised_child, extrapolation="linear"):
+    def boundary_value_or_flux(self, symbol, discretised_child):
         """
         Returns the boundary value or flux using the approriate expression for the
         spatial method. To do this, we create a sparse vector 'bv_vector' that extracts