separate boundary grad and value extrapolation, make grad default quadratic

CHANGELOG.md

@@ -4,6 +4,7 @@
 - Added two more submodels (options) for the SEI: Lars von Kolzenberg (2020) model and Tunneling Limit model ([#4394](https://github.com/pybamm-team/PyBaMM/pull/4394))
 
 ## Breaking changes
+- Separated extrapolation options for `pybamm.BoundaryValue` and `pybamm.BoundaryGradient`, and updated the default to be "linear" for the value and "quadratic" for the gradient. ([#4614](https://github.com/pybamm-team/PyBaMM/pull/4614))
 
 - Double-layer SEI models have been removed (with the corresponding parameters). All models assume now a single SEI layer. ([#4470](https://github.com/pybamm-team/PyBaMM/pull/4470))
 

examples/scripts/compare_extrapolations.py

@@ -11,7 +11,12 @@
 sim_lin.solve([0, 3600])
 
 model_quad = pybamm.lead_acid.Full()
-method_options = {"extrapolation": {"order": "quadratic", "use bcs": False}}
+method_options = {
+    "extrapolation": {
+        "order": {"gradient": "quadratic", "value": "quadratic"},
+        "use bcs": False,
+    }
+}
 spatial_methods = {
     "negative particle": pybamm.FiniteVolume(method_options),
     "positive particle": pybamm.FiniteVolume(method_options),

src/pybamm/spatial_methods/finite_volume.py

@@ -26,8 +26,11 @@ class FiniteVolume(pybamm.SpatialMethod):
 
     Parameters
     ----------
-    mesh : :class:`pybamm.Mesh`
-        Contains all the submeshes for discretisation
+    options : dict-like, optional
+        A dictionary of options to be passed to the spatial method. The only option
+        currently available is "extrapolation", which has options for "order" and "use_bcs".
+        It sets the order separately for `pybamm.BoundaryValue` and `pybamm.BoundaryGradient`.
+        Default is "linear" for the value and quadratic for the gradient.
     """
 
     def __init__(self, options=None):
@@ -824,7 +827,8 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
         if bcs is None:
             bcs = {}
 
-        extrap_order = self.options["extrapolation"]["order"]
+        extrap_order_gradient = self.options["extrapolation"]["order"]["gradient"]
+        extrap_order_value = self.options["extrapolation"]["order"]["value"]
         use_bcs = self.options["extrapolation"]["use bcs"]
 
         nodes = submesh.nodes
@@ -850,7 +854,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                 additive = bcs[child][symbol.side][0]
 
             elif symbol.side == "left":
-                if extrap_order == "linear":
+                if extrap_order_value == "linear":
                     # to find value at x* use formula:
                     # f(x*) = f_1 - (dx0 / dx1) (f_2 - f_1)
 
@@ -868,7 +872,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                         )
                         additive = pybamm.Scalar(0)
 
-                elif extrap_order == "quadratic":
+                elif extrap_order_value == "quadratic":
                     if use_bcs and pybamm.has_bc_of_form(
                         child, symbol.side, bcs, "Neumann"
                     ):
@@ -895,7 +899,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                     raise NotImplementedError
 
             elif symbol.side == "right":
-                if extrap_order == "linear":
+                if extrap_order_value == "linear":
                     if use_bcs and pybamm.has_bc_of_form(
                         child, symbol.side, bcs, "Neumann"
                     ):
@@ -917,7 +921,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                             shape=(1, prim_pts),
                         )
                         additive = pybamm.Scalar(0)
-                elif extrap_order == "quadratic":
+                elif extrap_order_value == "quadratic":
                     if use_bcs and pybamm.has_bc_of_form(
                         child, symbol.side, bcs, "Neumann"
                     ):
@@ -958,14 +962,14 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                 additive = bcs[child][symbol.side][0]
 
             elif symbol.side == "left":
-                if extrap_order == "linear":
+                if extrap_order_gradient == "linear":
                     # f'(x*) = (f_2 - f_1) / dx1
                     sub_matrix = (1 / dx1) * csr_matrix(
                         ([-1, 1], ([0, 0], [0, 1])), shape=(1, prim_pts)
                     )
                     additive = pybamm.Scalar(0)
 
-                elif extrap_order == "quadratic":
+                elif extrap_order_gradient == "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)
@@ -978,7 +982,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                     raise NotImplementedError
 
             elif symbol.side == "right":
-                if extrap_order == "linear":
+                if extrap_order_gradient == "linear":
                     # use formula:
                     # f'(x*) = (f_N - f_Nm1) / dxNm1
                     sub_matrix = (1 / dxNm1) * csr_matrix(
@@ -987,7 +991,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
                     )
                     additive = pybamm.Scalar(0)
 
-                elif extrap_order == "quadratic":
+                elif extrap_order_gradient == "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)

src/pybamm/spatial_methods/spatial_method.py

@@ -13,7 +13,12 @@ class SpatialMethod:
     """
 
     def __init__(self, options=None):
-        self.options = {"extrapolation": {"order": "linear", "use bcs": False}}
+        self.options = {
+            "extrapolation": {
+                "order": {"gradient": "quadratic", "value": "linear"},
+                "use bcs": False,
+            }
+        }
 
         # update double-layered dict
         if options:

tests/integration/test_models/test_full_battery_models/test_lead_acid/test_asymptotics_convergence.py

@@ -28,7 +28,12 @@ def test_leading_order_convergence(self):
         var_pts = {"x_n": 3, "x_s": 3, "x_p": 3}
         mesh = pybamm.Mesh(geometry, full_model.default_submesh_types, var_pts)
 
-        method_options = {"extrapolation": {"order": "linear", "use bcs": False}}
+        method_options = {
+            "extrapolation": {
+                "order": {"gradient": "linear", "value": "linear"},
+                "use bcs": False,
+            }
+        }
         spatial_methods = {
             "macroscale": pybamm.FiniteVolume(method_options),
             "current collector": pybamm.ZeroDimensionalSpatialMethod(method_options),

tests/unit/test_spatial_methods/test_finite_volume/test_extrapolation.py

@@ -59,8 +59,10 @@ def get_errors(function, method_options, pts, bcs=None):
 class TestExtrapolation:
     def test_convergence_without_bcs(self):
         # all tests are performed on x in [0, 1]
-        linear = {"extrapolation": {"order": "linear"}}
-        quad = {"extrapolation": {"order": "quadratic"}}
+        linear = {"extrapolation": {"order": {"gradient": "linear", "value": "linear"}}}
+        quad = {
+            "extrapolation": {"order": {"gradient": "quadratic", "value": "quadratic"}}
+        }
 
         def x_squared(x):
             y = x**2
@@ -147,8 +149,18 @@ def x_cubed(x):
 
         bcs = {"left": (left_val, "Dirichlet"), "right": (right_flux, "Neumann")}
 
-        linear = {"extrapolation": {"order": "linear", "use bcs": True}}
-        quad = {"extrapolation": {"order": "quadratic", "use bcs": True}}
+        linear = {
+            "extrapolation": {
+                "order": {"gradient": "linear", "value": "linear"},
+                "use bcs": True,
+            }
+        }
+        quad = {
+            "extrapolation": {
+                "order": {"gradient": "quadratic", "value": "quadratic"},
+                "use bcs": True,
+            }
+        }
         l_errors_lin_no_bc, r_errors_lin_no_bc = get_errors(x_cubed, linear, pts)
         l_errors_quad_no_bc, r_errors_quad_no_bc = get_errors(x_cubed, quad, pts)
 
@@ -202,8 +214,18 @@ def x_cubed(x):
 
         bcs = {"left": (left_flux, "Neumann"), "right": (right_val, "Dirichlet")}
 
-        linear = {"extrapolation": {"order": "linear", "use bcs": True}}
-        quad = {"extrapolation": {"order": "quadratic", "use bcs": True}}
+        linear = {
+            "extrapolation": {
+                "order": {"gradient": "linear", "value": "linear"},
+                "use bcs": True,
+            }
+        }
+        quad = {
+            "extrapolation": {
+                "order": {"gradient": "quadratic", "value": "quadratic"},
+                "use bcs": True,
+            }
+        }
         l_errors_lin_no_bc, r_errors_lin_no_bc = get_errors(x_cubed, linear, pts)
         l_errors_quad_no_bc, r_errors_quad_no_bc = get_errors(x_cubed, quad, pts)
 
@@ -237,7 +259,12 @@ def x_cubed(x):
     def test_linear_extrapolate_left_right(self):
         # create discretisation
         mesh = get_mesh_for_testing()
-        method_options = {"extrapolation": {"order": "linear", "use bcs": True}}
+        method_options = {
+            "extrapolation": {
+                "order": {"gradient": "linear", "value": "linear"},
+                "use bcs": True,
+            }
+        }
         spatial_methods = {
             "macroscale": pybamm.FiniteVolume(method_options),
             "negative particle": pybamm.FiniteVolume(method_options),
@@ -308,7 +335,12 @@ def test_linear_extrapolate_left_right(self):
     def test_quadratic_extrapolate_left_right(self):
         # create discretisation
         mesh = get_mesh_for_testing()
-        method_options = {"extrapolation": {"order": "quadratic", "use bcs": False}}
+        method_options = {
+            "extrapolation": {
+                "order": {"gradient": "quadratic", "value": "quadratic"},
+                "use bcs": False,
+            }
+        }
         spatial_methods = {
             "macroscale": pybamm.FiniteVolume(method_options),
             "negative particle": pybamm.FiniteVolume(method_options),
@@ -402,7 +434,12 @@ def test_extrapolate_on_nonuniform_grid(self):
         rpts = 10
         var_pts = {"r_n": rpts, "r_p": rpts}
         mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
-        method_options = {"extrapolation": {"order": "linear", "use bcs": False}}
+        method_options = {
+            "extrapolation": {
+                "order": {"gradient": "linear", "value": "linear"},
+                "use bcs": False,
+            }
+        }
         spatial_methods = {"negative particle": pybamm.FiniteVolume(method_options)}
         disc = pybamm.Discretisation(mesh, spatial_methods)
 
@@ -429,7 +466,12 @@ def test_extrapolate_on_nonuniform_grid(self):
     def test_extrapolate_2d_models(self):
         # create discretisation
         mesh = get_p2d_mesh_for_testing()
-        method_options = {"extrapolation": {"order": "linear", "use bcs": False}}
+        method_options = {
+            "extrapolation": {
+                "order": {"gradient": "linear", "value": "linear"},
+                "use bcs": False,
+            }
+        }
         spatial_methods = {
             "macroscale": pybamm.FiniteVolume(method_options),
             "negative particle": pybamm.FiniteVolume(method_options),