#704 added tests for linear and quadratic and on nonuniform grids

Author: Scottmar93

pybamm/discretisations/discretisation.py

@@ -707,7 +707,9 @@ 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)
+                return child_spatial_method.boundary_value_or_flux(
+                    symbol, disc_child, symbol.extrapolation
+                )
 
             else:
                 return symbol._unary_new_copy(disc_child)

pybamm/expression_tree/unary_operators.py

@@ -738,12 +738,15 @@ 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):
+    def __init__(self, child, side, extrapolation="quadratic"):
         super().__init__("boundary value", child, side)
+        self.extrapolation = extrapolation
 
 
 class BoundaryGradient(BoundaryOperator):
@@ -755,12 +758,15 @@ 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):
+    def __init__(self, child, side, extrapolation="linear"):
         super().__init__("boundary flux", child, side)
+        self.extrapolation = extrapolation
 
 
 #
@@ -850,7 +856,7 @@ def grad_squared(expression):
 #
 
 
-def surf(symbol, set_domain=False):
+def surf(symbol, set_domain=False, extrapolation="quadratic"):
     """convenience function for creating a right :class:`BoundaryValue`, usually in the
     spherical geometry
 
@@ -868,10 +874,10 @@ def surf(symbol, set_domain=False):
     if symbol.domain in [["negative electrode"], ["positive electrode"]] and isinstance(
         symbol, pybamm.PrimaryBroadcast
     ):
-        child_surf = boundary_value(symbol.orphans[0], "right")
+        child_surf = boundary_value(symbol.orphans[0], "right", extrapolation)
         out = pybamm.PrimaryBroadcast(child_surf, symbol.domain)
     else:
-        out = boundary_value(symbol, "right")
+        out = boundary_value(symbol, "right", extrapolation)
         if set_domain:
             if symbol.domain == ["negative particle"]:
                 out.domain = ["negative electrode"]
@@ -1015,7 +1021,7 @@ def yz_average(symbol):
         return Integral(symbol, [y, z]) / (l_y * l_z)
 
 
-def boundary_value(symbol, side):
+def boundary_value(symbol, side, extrapolation="quadratic"):
     """convenience function for creating a :class:`pybamm.BoundaryValue`
 
     Parameters
@@ -1040,7 +1046,7 @@ def boundary_value(symbol, side):
         return symbol.orphans[0]
     # Otherwise, calculate boundary value
     else:
-        return BoundaryValue(symbol, side)
+        return BoundaryValue(symbol, side, extrapolation)
 
 
 def r_average(symbol):

pybamm/spatial_methods/finite_volume.py

@@ -596,26 +596,75 @@ def boundary_value_or_flux(self, symbol, discretised_child, extrapolation="linea
             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":
+                    # to find value at x* use formula:
+                    # f(x*) = f_1 - (dx0 / dx1) (f_2 - f_1)
+                    # where
+                    # dx0 is distance between edge and first node
+                    # dx1 is distance between first and second nodes
+
                     sub_matrix = csr_matrix(
                         ([1 + (dx0 / dx1), -(dx0 / dx1)], ([0, 0], [0, 1])),
                         shape=(1, prim_pts),
                     )
+                elif extrapolation == "quadratic":
+                    # to find value at x* use formula:
+                    # see mathematica notebook at:
+                    # https://github.com/Scottmar93/extrapolation-coefficents/tree/master
+                    # f(x*) = a f_1 + b f_2 + c f_3
+
+                    a = (dx0 + dx1) * (dx0 + dx1 + dx2) / (dx1 * (dx1 + dx2))
+                    b = -dx0 * (dx0 + dx1 + dx2) / (dx1 * dx2)
+                    c = dx0 * (dx0 + dx1) / (dx2 * (dx1 + dx2))
+
+                    sub_matrix = csr_matrix(
+                        ([a, b, c], ([0, 0, 0], [0, 1, 2])), shape=(1, prim_pts),
+                    )
                 else:
                     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":
+                    # to find value at x* use formula:
+                    # f(x*) = f_N - (dxN / dxNm1) (f_N - f_Nm1)
+                    # where
+                    # dxN is distance between edge and last node
+                    # dxNm1 is distance between second lase and last nodes
+
                     sub_matrix = csr_matrix(
                         (
                             [-(dxN / dxNm1), 1 + (dxN / dxNm1)],
                             ([0, 0], [prim_pts - 2, prim_pts - 1]),
                         ),
                         shape=(1, prim_pts),
                     )
+                elif extrapolation == "quadratic":
+                    # to find value at x* use formula:
+                    # see mathematica notebook at:
+                    # https://github.com/Scottmar93/extrapolation-coefficents/tree/master
+                    # f(x*) = a f_N + b f_Nm1 + c f_Nm2
+
+                    a = (
+                        (dxN + dxNm1)
+                        * (dxN + dxNm1 + dxNm2)
+                        / (dxNm1 * (dxNm1 + dxNm2))
+                    )
+                    b = -dxN * (dxN + dxNm1 + dxNm2) / (dxNm1 * dxNm2)
+                    c = dxN * (dxN + dxNm1) / (dxNm2 * (dxNm1 + dxNm2))
+
+                    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
 

tests/unit/test_spatial_methods/test_finite_volume.py

@@ -56,7 +56,7 @@ def test_concatenation(self):
         with self.assertRaisesRegex(pybamm.ShapeError, "child must have size n_nodes"):
             fin_vol.concatenation(edges)
 
-    def test_extrapolate_left_right(self):
+    def test_linear_extrapolate_left_right(self):
         # create discretisation
         mesh = get_mesh_for_testing()
         spatial_methods = {
@@ -74,8 +74,8 @@ def test_extrapolate_left_right(self):
         # create variable
         var = pybamm.Variable("var", domain=whole_cell)
         # boundary value should work with something more complicated than a variable
-        extrap_left = pybamm.BoundaryValue(2 * var, "left")
-        extrap_right = pybamm.BoundaryValue(4 - var, "right")
+        extrap_left = pybamm.BoundaryValue(2 * var, "left", "linear")
+        extrap_right = pybamm.BoundaryValue(4 - var, "right", "linear")
         disc.set_variable_slices([var])
         extrap_left_disc = disc.process_symbol(extrap_left)
         extrap_right_disc = disc.process_symbol(extrap_right)
@@ -104,6 +104,80 @@ def test_extrapolate_left_right(self):
         self.assertEqual(extrap_flux_left_disc.evaluate(None, linear_y), 2)
         self.assertEqual(extrap_flux_right_disc.evaluate(None, linear_y), -1)
 
+        # Microscale
+        # create variable
+        var = pybamm.Variable("var", domain="negative particle")
+        surf_eqn = pybamm.surf(var, extrapolation="linear")
+        disc.set_variable_slices([var])
+        surf_eqn_disc = disc.process_symbol(surf_eqn)
+
+        # check constant extrapolates to constant
+        constant_y = np.ones_like(micro_submesh[0].nodes[:, np.newaxis])
+        self.assertEqual(surf_eqn_disc.evaluate(None, constant_y), 1.0)
+
+        # check linear variable extrapolates correctly
+        linear_y = micro_submesh[0].nodes
+        y_surf = micro_submesh[0].edges[-1]
+        np.testing.assert_array_almost_equal(
+            surf_eqn_disc.evaluate(None, linear_y), y_surf
+        )
+
+    def test_quadratic_extrapolate_left_right(self):
+        # create discretisation
+        mesh = get_mesh_for_testing()
+        spatial_methods = {
+            "macroscale": pybamm.FiniteVolume,
+            "negative particle": pybamm.FiniteVolume,
+            "current collector": pybamm.ZeroDimensionalMethod,
+        }
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        whole_cell = ["negative electrode", "separator", "positive electrode"]
+        macro_submesh = mesh.combine_submeshes(*whole_cell)
+        micro_submesh = mesh["negative particle"]
+
+        # Macroscale
+        # create variable
+        var = pybamm.Variable("var", domain=whole_cell)
+        # boundary value should work with something more complicated than a variable
+        extrap_left = pybamm.BoundaryValue(2 * var, "left", "quadratic")
+        extrap_right = pybamm.BoundaryValue(4 - var, "right", "quadratic")
+        disc.set_variable_slices([var])
+        extrap_left_disc = disc.process_symbol(extrap_left)
+        extrap_right_disc = disc.process_symbol(extrap_right)
+
+        # check constant extrapolates to constant
+        constant_y = np.ones_like(macro_submesh[0].nodes[:, np.newaxis])
+        np.testing.assert_array_almost_equal(
+            extrap_left_disc.evaluate(None, constant_y), 2.0
+        )
+        np.testing.assert_array_almost_equal(
+            extrap_right_disc.evaluate(None, constant_y), 3.0
+        )
+
+        # check linear variable extrapolates correctly
+        linear_y = macro_submesh[0].nodes
+        np.testing.assert_array_almost_equal(
+            extrap_left_disc.evaluate(None, linear_y), 0
+        )
+        np.testing.assert_array_almost_equal(
+            extrap_right_disc.evaluate(None, linear_y), 3
+        )
+
+        # Fluxes
+        extrap_flux_left = pybamm.BoundaryGradient(2 * var, "left")
+        extrap_flux_right = pybamm.BoundaryGradient(1 - var, "right")
+        extrap_flux_left_disc = disc.process_symbol(extrap_flux_left)
+        extrap_flux_right_disc = disc.process_symbol(extrap_flux_right)
+
+        # check constant extrapolates to constant
+        self.assertEqual(extrap_flux_left_disc.evaluate(None, constant_y), 0)
+        self.assertEqual(extrap_flux_right_disc.evaluate(None, constant_y), 0)
+
+        # check linear variable extrapolates correctly
+        self.assertEqual(extrap_flux_left_disc.evaluate(None, linear_y), 2)
+        self.assertEqual(extrap_flux_right_disc.evaluate(None, linear_y), -1)
+
         # Microscale
         # create variable
         var = pybamm.Variable("var", domain="negative particle")
@@ -113,11 +187,53 @@ def test_extrapolate_left_right(self):
 
         # check constant extrapolates to constant
         constant_y = np.ones_like(micro_submesh[0].nodes[:, np.newaxis])
-        self.assertEqual(surf_eqn_disc.evaluate(None, constant_y), 1)
+        np.testing.assert_array_almost_equal(
+            surf_eqn_disc.evaluate(None, constant_y), 1
+        )
 
         # check linear variable extrapolates correctly
         linear_y = micro_submesh[0].nodes
-        y_surf = micro_submesh[0].nodes[-1] + micro_submesh[0].d_nodes[-1] / 2
+        y_surf = micro_submesh[0].edges[-1]
+        np.testing.assert_array_almost_equal(
+            surf_eqn_disc.evaluate(None, linear_y), y_surf
+        )
+
+    def test_extrapolate_on_nonuniform_grid(self):
+        geometry = pybamm.Geometry("1D micro")
+
+        submesh_types = {
+            "negative particle": pybamm.MeshGenerator(pybamm.Exponential1DSubMesh),
+            "positive particle": pybamm.MeshGenerator(pybamm.Exponential1DSubMesh),
+        }
+
+        var = pybamm.standard_spatial_vars
+        rpts = 10
+        var_pts = {
+            var.r_n: rpts,
+            var.r_p: rpts,
+        }
+        mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
+        spatial_methods = {
+            "negative particle": pybamm.FiniteVolume,
+        }
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        var = pybamm.Variable("var", domain="negative particle")
+        surf_eqn = pybamm.surf(var)
+        disc.set_variable_slices([var])
+        surf_eqn_disc = disc.process_symbol(surf_eqn)
+
+        micro_submesh = mesh["negative particle"]
+
+        # check constant extrapolates to constant
+        constant_y = np.ones_like(micro_submesh[0].nodes[:, np.newaxis])
+        np.testing.assert_array_almost_equal(
+            surf_eqn_disc.evaluate(None, constant_y), 1
+        )
+
+        # check linear variable extrapolates correctly
+        linear_y = micro_submesh[0].nodes
+        y_surf = micro_submesh[0].edges[-1]
         np.testing.assert_array_almost_equal(
             surf_eqn_disc.evaluate(None, linear_y), y_surf
         )
@@ -1033,7 +1149,7 @@ def test_indefinite_integral(self):
             }
         }
         int_grad_phi_disc = disc.process_symbol(int_grad_phi)
-        left_boundary_value = pybamm.BoundaryValue(int_grad_phi, "left")
+        left_boundary_value = pybamm.BoundaryValue(int_grad_phi, "left", "linear")
         left_boundary_value_disc = disc.process_symbol(left_boundary_value)
 
         combined_submesh = mesh.combine_submeshes("negative electrode", "separator")
@@ -1087,15 +1203,19 @@ def test_indefinite_integral(self):
         )
         phi_approx = int_grad_phi_disc.evaluate(None, phi_exact)
         np.testing.assert_array_almost_equal(phi_exact, phi_approx)
-        self.assertEqual(left_boundary_value_disc.evaluate(y=phi_exact), 0)
+        np.testing.assert_array_almost_equal(
+            left_boundary_value_disc.evaluate(y=phi_exact), 0
+        )
 
         # sine case
         phi_exact = np.sin(
             combined_submesh[0].nodes[:, np.newaxis] - combined_submesh[0].edges[0]
         )
         phi_approx = int_grad_phi_disc.evaluate(None, phi_exact)
         np.testing.assert_array_almost_equal(phi_exact, phi_approx)
-        self.assertEqual(left_boundary_value_disc.evaluate(y=phi_exact), 0)
+        np.testing.assert_array_almost_equal(
+            left_boundary_value_disc.evaluate(y=phi_exact), 0
+        )
 
         # --------------------------------------------------------------------
         # micrsoscale case
@@ -1112,7 +1232,7 @@ def test_indefinite_integral(self):
         }
 
         c_integral_disc = disc.process_symbol(c_integral)
-        left_boundary_value = pybamm.BoundaryValue(c_integral, "left")
+        left_boundary_value = pybamm.BoundaryValue(c_integral, "left", "linear")
         left_boundary_value_disc = disc.process_symbol(left_boundary_value)
         combined_submesh = mesh["negative particle"]
 
@@ -1127,13 +1247,17 @@ def test_indefinite_integral(self):
         c_exact = combined_submesh[0].nodes[:, np.newaxis]
         c_approx = c_integral_disc.evaluate(None, c_exact)
         np.testing.assert_array_almost_equal(c_exact, c_approx)
-        self.assertEqual(left_boundary_value_disc.evaluate(y=c_exact), 0)
+        np.testing.assert_array_almost_equal(
+            left_boundary_value_disc.evaluate(y=c_exact), 0
+        )
 
         # sine case
         c_exact = np.sin(combined_submesh[0].nodes[:, np.newaxis])
         c_approx = c_integral_disc.evaluate(None, c_exact)
         np.testing.assert_array_almost_equal(c_exact, c_approx, decimal=3)
-        self.assertEqual(left_boundary_value_disc.evaluate(y=c_exact), 0)
+        np.testing.assert_array_almost_equal(
+            left_boundary_value_disc.evaluate(y=c_exact), 0
+        )
 
     def test_indefinite_integral_on_nodes(self):
         mesh = get_mesh_for_testing()