Issue 704 surface extrapolation

examples/notebooks/create-model.ipynb

@@ -600,7 +600,7 @@
     "var_pts = {x: 100}\n",
     "mesh = pybamm.Mesh(geometry, submesh_types, var_pts)\n",
     "                   \n",
-    "spatial_methods = {\"SEI layer\": pybamm.FiniteVolume}\n",
+    "spatial_methods = {\"SEI layer\": pybamm.FiniteVolume()}\n",
     "disc = pybamm.Discretisation(mesh, spatial_methods)\n",
     "disc.process_model(model)"
    ]
@@ -710,7 +710,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.8"
+   "version": "3.7.3"
   }
  },
  "nbformat": 4,

examples/notebooks/spatial_methods/finite-volumes.ipynb

@@ -124,9 +124,9 @@
    "outputs": [],
    "source": [
     "spatial_methods = {\n",
-    "    \"macroscale\": pybamm.FiniteVolume,\n",
-    "    \"negative particle\": pybamm.FiniteVolume,\n",
-    "    \"positive particle\": pybamm.FiniteVolume,\n",
+    "    \"macroscale\": pybamm.FiniteVolume(),\n",
+    "    \"negative particle\": pybamm.FiniteVolume(),\n",
+    "    \"positive particle\": pybamm.FiniteVolume(),\n",
     "}\n",
     "disc = pybamm.Discretisation(mesh, spatial_methods)"
    ]
@@ -1282,7 +1282,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.7"
+   "version": "3.7.3"
   }
  },
  "nbformat": 4,

examples/notebooks/unsteady_heat_equation.ipynb

@@ -214,7 +214,7 @@
     "submesh_types = {\"rod\": pybamm.Uniform1DSubMesh}\n",
     "var_pts = {x: 30}\n",
     "mesh = pybamm.Mesh(geometry, submesh_types, var_pts)\n",
-    "spatial_methods = {\"rod\": pybamm.FiniteVolume}\n",
+    "spatial_methods = {\"rod\": pybamm.FiniteVolume()}\n",
     "disc = pybamm.Discretisation(mesh, spatial_methods)"
    ]
   },
@@ -421,7 +421,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.8"
+   "version": "3.7.3"
   }
  },
  "nbformat": 4,

examples/notebooks/using-model-options_thermal-example.ipynb

@@ -204,7 +204,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.8"
+   "version": "3.7.3"
   }
  },
  "nbformat": 4,

examples/scripts/compare_extrapolations.py

@@ -0,0 +1,32 @@
+import pybamm
+
+
+x_n = pybamm.standard_spatial_vars.x_n
+x_s = pybamm.standard_spatial_vars.x_s
+x_p = pybamm.standard_spatial_vars.x_p
+
+var_pts = {x_n: 10, x_s: 3, x_p: 10}
+model_lin = pybamm.lead_acid.Full()
+sim_lin = pybamm.Simulation(model_lin, var_pts=var_pts)
+sim_lin.solve()
+
+model_quad = pybamm.lead_acid.Full()
+method_options = {"extrapolation": {"order": "quadratic", "use bcs": False}}
+spatial_methods = {
+    "negative particle": pybamm.FiniteVolume(method_options),
+    "positive particle": pybamm.FiniteVolume(method_options),
+    "macroscale": pybamm.FiniteVolume(method_options),
+    "current collector": pybamm.ZeroDimensionalMethod(),
+}
+sim_quad = pybamm.Simulation(
+    model_quad, spatial_methods=spatial_methods, var_pts=var_pts
+)
+sim_quad.solve()
+
+
+# plot the two sols
+models = [sim_lin.built_model, sim_quad.built_model]
+solutions = [sim_lin.solution, sim_quad.solution]
+plot = pybamm.QuickPlot(models, sim_lin.mesh, solutions)
+plot.dynamic_plot()
+

examples/scripts/create-model.py

@@ -105,7 +105,7 @@ def Diffusivity(cc):
 var_pts = {x: 50}
 mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
 
-spatial_methods = {"SEI layer": pybamm.FiniteVolume}
+spatial_methods = {"SEI layer": pybamm.FiniteVolume()}
 disc = pybamm.Discretisation(mesh, spatial_methods)
 disc.process_model(model)
 

examples/scripts/heat_equation.py

@@ -51,7 +51,7 @@
 submesh_types = {"rod": pybamm.Uniform1DSubMesh}
 var_pts = {x: 30}
 mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
-spatial_methods = {"rod": pybamm.FiniteVolume}
+spatial_methods = {"rod": pybamm.FiniteVolume()}
 disc = pybamm.Discretisation(mesh, spatial_methods)
 disc.process_model(model)
 

pybamm/__init__.py

@@ -229,6 +229,7 @@ def version(formatted=False):
 # Mesh and Discretisation classes
 #
 from .discretisations.discretisation import Discretisation
+from .discretisations.discretisation import has_bc_of_form
 from .meshes.meshes import Mesh, SubMesh, MeshGenerator
 from .meshes.zero_dimensional_submesh import SubMesh0D
 from .meshes.one_dimensional_submeshes import (

pybamm/discretisations/discretisation.py

@@ -7,6 +7,18 @@
 from scipy.sparse import block_diag, csr_matrix
 
 
+def has_bc_of_form(symbol, side, bcs, form):
+
+    if symbol.id in bcs:
+        if bcs[symbol.id][side][1] == form:
+            return True
+        else:
+            return False
+
+    else:
+        return False
+
+
 class Discretisation(object):
     """The discretisation class, with methods to process a model and replace
     Spatial Operators with Matrices and Variables with StateVectors
@@ -16,8 +28,9 @@ class Discretisation(object):
     mesh : pybamm.Mesh
             contains all submeshes to be used on each domain
     spatial_methods : dict
-            a dictionary of the spatial method to be used on each
-            domain. The keys correspond to the keys in a pybamm.Model
+            a dictionary of the spatial methods to be used on each
+            domain. The keys correspond to the model domains and the
+            values to the spatial method.
     """
 
     def __init__(self, mesh=None, spatial_methods=None):
@@ -31,9 +44,11 @@ def __init__(self, mesh=None, spatial_methods=None):
                 spatial_methods["negative electrode"] = method
                 spatial_methods["separator"] = method
                 spatial_methods["positive electrode"] = method
-            self._spatial_methods = {
-                dom: method(mesh) for dom, method in spatial_methods.items()
-            }
+
+            self._spatial_methods = spatial_methods
+            for method in self._spatial_methods.values():
+                method.build(mesh)
+
         self.bcs = {}
         self.y_slices = {}
         self._discretised_symbols = {}
@@ -714,7 +729,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, self.bcs
+                )
 
             else:
                 return symbol._unary_new_copy(disc_child)

pybamm/models/full_battery_models/base_battery_model.py

@@ -121,17 +121,17 @@ def default_submesh_types(self):
     @property
     def default_spatial_methods(self):
         base_spatial_methods = {
-            "macroscale": pybamm.FiniteVolume,
-            "negative particle": pybamm.FiniteVolume,
-            "positive particle": pybamm.FiniteVolume,
+            "macroscale": pybamm.FiniteVolume(),
+            "negative particle": pybamm.FiniteVolume(),
+            "positive particle": pybamm.FiniteVolume(),
         }
         if self.options["dimensionality"] == 0:
             # 0D submesh - use base spatial method
-            base_spatial_methods["current collector"] = pybamm.ZeroDimensionalMethod
+            base_spatial_methods["current collector"] = pybamm.ZeroDimensionalMethod()
         elif self.options["dimensionality"] == 1:
-            base_spatial_methods["current collector"] = pybamm.FiniteVolume
+            base_spatial_methods["current collector"] = pybamm.FiniteVolume()
         elif self.options["dimensionality"] == 2:
-            base_spatial_methods["current collector"] = pybamm.ScikitFiniteElement
+            base_spatial_methods["current collector"] = pybamm.ScikitFiniteElement()
         return base_spatial_methods
 
     @property

pybamm/models/submodels/current_collector/effective_resistance_current_collector.py

@@ -223,7 +223,7 @@ def default_submesh_types(self):
 
     @property
     def default_spatial_methods(self):
-        return {"current collector": pybamm.ScikitFiniteElement}
+        return {"current collector": pybamm.ScikitFiniteElement()}
 
     @property
     def default_solver(self):

pybamm/simulation.py

@@ -183,6 +183,10 @@ def parameter_values(self):
     def submesh_types(self):
         return self._submesh_types
 
+    @property
+    def mesh(self):
+        return self._mesh
+
     @property
     def var_pts(self):
         return self._var_pts