Repository: pybamm-team/PyBaMM

Issue:
Refactor `add_ghost_nodes` to work with different symbols
`add_ghost_nodes` currently needs to access `y_slices` and so only works with `StateVectors` or concatenations of.
Refactor (e.g. using matrix multiplication instead) so that it can work with any kind of symbol. Needed for some of the different models.

Discussion:

Patch:
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 587a739d5..3f69f29ba 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -506,7 +506,12 @@ class Discretisation(object):
 
     def check_model(self, model):
         """ Perform some basic checks to make sure the discretised model makes sense."""
-        # Check initial conditions are a numpy array
+        self.check_initial_conditions(model)
+        self.check_initial_conditions_rhs(model)
+        self.check_variables(model)
+
+    def check_initial_conditions(self, model):
+        """Check initial conditions are a numpy array"""
         # Individual
         for var, eqn in model.initial_conditions.items():
             assert type(eqn.evaluate(0, None)) is np.ndarray, pybamm.ModelError(
@@ -529,7 +534,8 @@ class Discretisation(object):
             )
         )
 
-        # Check initial conditions and rhs have the same shape
+    def check_initial_conditions_rhs(self, model):
+        """Check initial conditions and rhs have the same shape"""
         y0 = model.concatenated_initial_conditions
         # Individual
         for var in model.rhs.keys():
@@ -563,9 +569,13 @@ class Discretisation(object):
             )
         )
 
-        # Check variables in variable list against rhs
-        # Be lenient with size check if the variable in model.variables is broadcasted
-        # If broadcasted, variable is a multiplication with a vector of ones
+    def check_variables(self, model):
+        """
+        Check variables in variable list against rhs
+        Be lenient with size check if the variable in model.variables is broadcasted
+        (if broadcasted, variable is a multiplication with a vector of ones)
+        """
+        y0 = model.concatenated_initial_conditions
         for rhs_var in model.rhs.keys():
             if rhs_var.name in model.variables.keys():
                 var = model.variables[rhs_var.name]
diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index d486fa6a1..17ad835ee 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -202,23 +202,59 @@ class Function(UnaryOperator):
 class Index(UnaryOperator):
     """A node in the expression tree, which stores the index that should be
     extracted from its child after the child has been evaluated.
+
+    Parameters
+    ----------
+    child : :class:`pybamm.Symbol`
+        The symbol of which to take the index
+    index : int or slice
+        The index (if int) or indices (if slice) to extract from the symbol
+    name : str, optional
+        The name of the symbol
     """
 
     def __init__(self, child, index, name=None):
-        if name is None:
-            name = child.name + "[" + str(index) + "]"
         self.index = index
+        if isinstance(index, int):
+            self.slice = slice(index, index + 1)
+            if name is None:
+                name = "Index[" + str(index) + "]"
+        elif isinstance(index, slice):
+            self.slice = index
+            if name is None:
+                if index.start is None:
+                    name = "Index[:{:d}]".format(index.stop)
+                else:
+                    name = "Index[{:d}:{:d}]".format(index.start, index.stop)
+        else:
+            raise TypeError("index must be integer or slice")
+
+        if self.slice.stop > child.size:
+            raise ValueError("slice size exceeds child size")
+
         super().__init__(name, child)
 
+    def jac(self, variable):
+        """ See :meth:`pybamm.Symbol.jac()`. """
+        child_jac = self.child.jac(variable)
+        return Index(child_jac, self.index)
+
     def set_id(self):
         """ See :meth:`pybamm.Symbol.set_id()` """
         self._id = hash(
-            (self.__class__, self.name, self.index, self.child.id) + tuple(self.domain)
+            (
+                self.__class__,
+                self.name,
+                self.slice.start,
+                self.slice.stop,
+                self.child.id,
+            )
+            + tuple(self.domain)
         )
 
     def _unary_evaluate(self, child):
         """ See :meth:`UnaryOperator._unary_evaluate()`. """
-        return child[self.index]
+        return child[self.slice]
 
     def _unary_new_copy(self, child):
         """ See :meth:`UnaryOperator._unary_new_copy()`. """
diff --git a/pybamm/models/base_models.py b/pybamm/models/base_models.py
index 859ec3fb7..9401567ec 100644
--- a/pybamm/models/base_models.py
+++ b/pybamm/models/base_models.py
@@ -244,6 +244,13 @@ class BaseModel(object):
         post_discretisation : boolean
             A flag indicating tests to be skipped after discretisation
         """
+        self.check_well_determined(post_discretisation)
+        self.check_algebraic_equations(post_discretisation)
+        self.check_ics_bcs()
+        self.check_variables()
+
+    def check_well_determined(self, post_discretisation):
+        """ Check that the model is not under- or over-determined. """
         # Equations (differential and algebraic)
         # Get all the variables from differential and algebraic equations
         vars_in_rhs_keys = set()
@@ -285,9 +292,14 @@ class BaseModel(object):
         extra_variables = vars_in_eqns.difference(vars_in_keys)
         if extra_variables:
             raise pybamm.ModelError("model is underdetermined (too many variables)")
-        # Before discretisation, each algebraic equation key must appear in the equation
-        # After discretisation, there must be at least one StateVector in each algebraic
-        # equation
+
+    def check_algebraic_equations(self, post_discretisation):
+        """
+        Check that the algebraic equations are well-posed.
+        Before discretisation, each algebraic equation key must appear in the equation
+        After discretisation, there must be at least one StateVector in each algebraic
+        equation
+        """
         if not post_discretisation:
             # After the model has been defined, each algebraic equation key should
             # appear in that algebraic equation
@@ -306,6 +318,8 @@ class BaseModel(object):
                         "each algebraic equation must contain at least one StateVector"
                     )
 
+    def check_ics_bcs(self):
+        """ Check that the initial and boundary conditions are well-posed. """
         # Initial conditions
         for var in self.rhs.keys():
             if var not in self.initial_conditions.keys():
@@ -328,7 +342,8 @@ class BaseModel(object):
                         )
                     )
 
-        # Standard Output Variables
+    def check_variables(self):
+        """ Chec that the right variables are provided. """
         missing_vars = []
         for output, expression in self._variables.items():
             if expression is None:
@@ -355,11 +370,14 @@ class StandardBatteryBaseModel(BaseModel):
 
     def __init__(self):
         super().__init__()
+        self.set_standard_output_variables()
 
+    @property
+    def default_parameter_values(self):
         # Default parameter values, geometry, submesh, spatial methods and solver
         # Lion parameters left as default parameter set for tests
         input_path = os.path.join(os.getcwd(), "input", "parameters", "lithium-ion")
-        self.default_parameter_values = pybamm.ParameterValues(
+        return pybamm.ParameterValues(
             os.path.join(
                 input_path, "mcmb2528_lif6-in-ecdmc_lico2_parameters_Dualfoil.csv"
             ),
@@ -392,32 +410,47 @@ class StandardBatteryBaseModel(BaseModel):
                 ),
             },
         )
-        self.default_geometry = pybamm.Geometry("1D macro", "1+1D micro")
+
+    @property
+    def default_geometry(self):
+        return pybamm.Geometry("1D macro", "1+1D micro")
+
+    @property
+    def default_var_pts(self):
         var = pybamm.standard_spatial_vars
-        self.default_var_pts = {
-            var.x_n: 40,
-            var.x_s: 25,
-            var.x_p: 35,
-            var.r_n: 10,
-            var.r_p: 10,
-        }
-        self.default_submesh_types = {
+        return {var.x_n: 40, var.x_s: 25, var.x_p: 35, var.r_n: 10, var.r_p: 10}
+
+    @property
+    def default_submesh_types(self):
+        return {
             "negative electrode": pybamm.Uniform1DSubMesh,
             "separator": pybamm.Uniform1DSubMesh,
             "positive electrode": pybamm.Uniform1DSubMesh,
             "negative particle": pybamm.Uniform1DSubMesh,
             "positive particle": pybamm.Uniform1DSubMesh,
         }
-        self.default_spatial_methods = {
+
+    @property
+    def default_spatial_methods(self):
+        return {
             "macroscale": pybamm.FiniteVolume,
             "negative particle": pybamm.FiniteVolume,
             "positive particle": pybamm.FiniteVolume,
         }
+
+    @property
+    def default_solver(self):
+        """
+        Create and return the default solver for this model
+        """
         try:
-            self.default_solver = pybamm.ScikitsOdeSolver()
+            default_solver = pybamm.ScikitsOdeSolver()
         except ImportError:
-            self.default_solver = pybamm.ScipySolver()
+            default_solver = pybamm.ScipySolver()
+
+        return default_solver
 
+    def set_standard_output_variables(self):
         # Standard output variables
         # Interfacial current
         self.variables.update(
@@ -555,9 +588,10 @@ class LeadAcidBaseModel(StandardBatteryBaseModel):
     def __init__(self):
         super().__init__()
 
-        # Overwrite default parameter values
+    @property
+    def default_parameter_values(self):
         input_path = os.path.join(os.getcwd(), "input", "parameters", "lead-acid")
-        self.default_parameter_values = pybamm.ParameterValues(
+        return pybamm.ParameterValues(
             "input/parameters/lead-acid/default.csv",
             {
                 "Typical current [A]": 1,
@@ -589,9 +623,12 @@ class LeadAcidBaseModel(StandardBatteryBaseModel):
             },
         )
 
-        # Overwrite geometry
-        self.default_geometry = pybamm.Geometry("1D macro")
+    @property
+    def default_geometry(self):
+        return pybamm.Geometry("1D macro")
 
+    def set_standard_output_variables(self):
+        super().set_standard_output_variables()
         # Standard time variable
         time_scale = pybamm.standard_parameters_lead_acid.tau_discharge
         I = pybamm.electrical_parameters.dimensional_current_with_time
@@ -617,6 +654,8 @@ class LithiumIonBaseModel(StandardBatteryBaseModel):
     def __init__(self):
         super().__init__()
 
+    def set_standard_output_variables(self):
+        super().set_standard_output_variables()
         # Additional standard output variables
         # Time
         time_scale = pybamm.standard_parameters_lithium_ion.tau_discharge
diff --git a/pybamm/models/lead_acid/composite_capacitance.py b/pybamm/models/lead_acid/composite_capacitance.py
index 8151b31a9..17d2ffba0 100644
--- a/pybamm/models/lead_acid/composite_capacitance.py
+++ b/pybamm/models/lead_acid/composite_capacitance.py
@@ -162,8 +162,15 @@ class CompositeCapacitance(pybamm.LeadAcidBaseModel):
         # Don't use jacobian for now (simplifications failing)
         self.use_jacobian = False
 
+    @property
+    def default_solver(self):
+        """
+        Create and return the default solver for this model
+        """
         # Different solver depending on whether we solve ODEs or DAEs
-        if use_capacitance:
-            self.default_solver = pybamm.ScikitsOdeSolver()
+        if self._use_capacitance:
+            default_solver = pybamm.ScikitsOdeSolver()
         else:
-            self.default_solver = pybamm.ScikitsDaeSolver()
+            default_solver = pybamm.ScikitsDaeSolver()
+
+        return default_solver
diff --git a/pybamm/models/lead_acid/loqs.py b/pybamm/models/lead_acid/loqs.py
index e6e6d1989..a0c96bbd4 100644
--- a/pybamm/models/lead_acid/loqs.py
+++ b/pybamm/models/lead_acid/loqs.py
@@ -99,7 +99,10 @@ class LOQS(pybamm.LeadAcidBaseModel):
         # ODEs only (don't use jacobian, use base spatial method)
         self.use_jacobian = False
 
-        self.default_spatial_methods = {
+    @property
+    def default_spatial_methods(self):
+        # ODEs only, so use base spatial method
+        return {
             "macroscale": pybamm.SpatialMethod,
             "negative particle": pybamm.SpatialMethod,
             "positive particle": pybamm.SpatialMethod,
diff --git a/pybamm/models/lead_acid/loqs_capacitance.py b/pybamm/models/lead_acid/loqs_capacitance.py
index 886b52c7c..9b4156c5b 100644
--- a/pybamm/models/lead_acid/loqs_capacitance.py
+++ b/pybamm/models/lead_acid/loqs_capacitance.py
@@ -130,8 +130,15 @@ class LOQSCapacitance(pybamm.LeadAcidBaseModel):
         # ODE model, don't use Jacobian
         self.use_jacobian = False
 
+    @property
+    def default_solver(self):
+        """
+        Create and return the default solver for this model
+        """
         # Different solver depending on whether we solve ODEs or DAEs
-        if use_capacitance:
-            self.default_solver = pybamm.ScikitsOdeSolver()
+        if self._use_capacitance:
+            default_solver = pybamm.ScikitsOdeSolver()
         else:
-            self.default_solver = pybamm.ScikitsDaeSolver()
+            default_solver = pybamm.ScikitsDaeSolver()
+
+        return default_solver
diff --git a/pybamm/models/lead_acid/newman_tiedemann.py b/pybamm/models/lead_acid/newman_tiedemann.py
index fed9b38b7..786f13d14 100644
--- a/pybamm/models/lead_acid/newman_tiedemann.py
+++ b/pybamm/models/lead_acid/newman_tiedemann.py
@@ -115,7 +115,10 @@ class NewmanTiedemann(pybamm.LeadAcidBaseModel):
         voltage = self.variables["Terminal voltage"]
         self.events.append(voltage - param.voltage_low_cut)
 
-        "-----------------------------------------------------------------------------"
-        "Defaults and Solver Conditions"
+    @property
+    def default_solver(self):
+        """
+        Create and return the default solver for this model
+        """
         # Default solver to DAE
-        self.default_solver = pybamm.ScikitsDaeSolver()
+        return pybamm.ScikitsDaeSolver()
diff --git a/pybamm/models/lead_acid/newman_tiedemann_capacitance.py b/pybamm/models/lead_acid/newman_tiedemann_capacitance.py
index 01bac033e..4d8108a03 100644
--- a/pybamm/models/lead_acid/newman_tiedemann_capacitance.py
+++ b/pybamm/models/lead_acid/newman_tiedemann_capacitance.py
@@ -133,11 +133,15 @@ class NewmanTiedemannCapacitance(pybamm.LeadAcidBaseModel):
         )
         self.events.append(voltage - param.voltage_low_cut)
 
-        "-----------------------------------------------------------------------------"
-        "Extra settings"
-
+    @property
+    def default_solver(self):
+        """
+        Create and return the default solver for this model
+        """
         # Different solver depending on whether we solve ODEs or DAEs
-        if use_capacitance:
-            self.default_solver = pybamm.ScikitsOdeSolver()
+        if self._use_capacitance:
+            default_solver = pybamm.ScikitsOdeSolver()
         else:
-            self.default_solver = pybamm.ScikitsDaeSolver()
+            default_solver = pybamm.ScikitsDaeSolver()
+
+        return default_solver
diff --git a/pybamm/models/lithium_ion/dfn.py b/pybamm/models/lithium_ion/dfn.py
index 3ecf49ec0..3bf8a3e0b 100644
--- a/pybamm/models/lithium_ion/dfn.py
+++ b/pybamm/models/lithium_ion/dfn.py
@@ -109,9 +109,15 @@ class DFN(pybamm.LithiumIonBaseModel):
         voltage = self.variables["Terminal voltage"]
         self.events.append(voltage - param.voltage_low_cut)
 
-        "-----------------------------------------------------------------------------"
-        "Defaults and Solver Conditions"
-        # default geometry
-        self.default_geometry = pybamm.Geometry("1D macro", "1+1D micro")
+    @property
+    def default_geometry(self):
+        return pybamm.Geometry("1D macro", "1+1D micro")
+
+    @property
+    def default_solver(self):
+        """
+        Create and return the default solver for this model
+        """
+
         # Default solver to DAE
-        self.default_solver = pybamm.ScikitsDaeSolver()
+        return pybamm.ScikitsDaeSolver()
diff --git a/pybamm/models/lithium_ion/spm.py b/pybamm/models/lithium_ion/spm.py
index dbd15a841..f0f23a539 100644
--- a/pybamm/models/lithium_ion/spm.py
+++ b/pybamm/models/lithium_ion/spm.py
@@ -85,7 +85,6 @@ class SPM(pybamm.LithiumIonBaseModel):
         voltage = self.variables["Terminal voltage"]
         self.events.append(voltage - param.voltage_low_cut)
 
-        "-----------------------------------------------------------------------------"
-        "Defaults and Solver Conditions"
-        # default geometry
-        self.default_geometry = pybamm.Geometry("1D macro", "1D micro")
+    @property
+    def default_geometry(self):
+        return pybamm.Geometry("1D macro", "1D micro")
diff --git a/pybamm/models/lithium_ion/spme.py b/pybamm/models/lithium_ion/spme.py
index ae637a384..51a0599e9 100644
--- a/pybamm/models/lithium_ion/spme.py
+++ b/pybamm/models/lithium_ion/spme.py
@@ -105,7 +105,6 @@ class SPMe(pybamm.LithiumIonBaseModel):
         voltage = self.variables["Terminal voltage"]
         self.events.append(voltage - param.voltage_low_cut)
 
-        "-----------------------------------------------------------------------------"
-        "Defaults and Solver Conditions"
-        # default geometry
-        self.default_geometry = pybamm.Geometry("1D macro", "1D micro")
+    @property
+    def default_geometry(self):
+        return pybamm.Geometry("1D macro", "1D micro")
diff --git a/pybamm/models/reaction_diffusion.py b/pybamm/models/reaction_diffusion.py
index 8b4a34557..d0293e797 100644
--- a/pybamm/models/reaction_diffusion.py
+++ b/pybamm/models/reaction_diffusion.py
@@ -61,9 +61,6 @@ class ReactionDiffusionModel(pybamm.StandardBatteryBaseModel):
         j_vars = int_curr_model.get_derived_interfacial_currents(j_n, j_p, j0_n, j0_p)
         self.variables.update(j_vars)
 
-        "-----------------------------------------------------------------------------"
-        "Defaults and Solver Conditions"
-        # default geometry
-        self.default_parameter_values = (
-            pybamm.LeadAcidBaseModel().default_parameter_values
-        )
+    @property
+    def default_parameter_values(self):
+        return pybamm.LeadAcidBaseModel().default_parameter_values
diff --git a/pybamm/models/submodels/electrode.py b/pybamm/models/submodels/electrode.py
index 26001447c..635f8614c 100644
--- a/pybamm/models/submodels/electrode.py
+++ b/pybamm/models/submodels/electrode.py
@@ -78,8 +78,12 @@ class Ohm(pybamm.SubModel):
         else:
             raise pybamm.DomainError("domain '{}' not recognised".format(phi_s.domain))
 
-        # Set default solver to DAE
-        self.default_solver = pybamm.ScikitsDaeSolver()
+    @property
+    def default_solver(self):
+        """
+        Create and return the default solver for this model
+        """
+        return pybamm.ScikitsDaeSolver()
 
     def get_explicit_leading_order(self, ocp_p, eta_r_p, phi_e):
         """
diff --git a/pybamm/models/submodels/electrolyte_current.py b/pybamm/models/submodels/electrolyte_current.py
index a68b71bbc..ba39032ec 100644
--- a/pybamm/models/submodels/electrolyte_current.py
+++ b/pybamm/models/submodels/electrolyte_current.py
@@ -332,8 +332,12 @@ class MacInnesStefanMaxwell(ElectrolyteCurrentBaseModel):
 
         self.variables = self.get_variables(phi_e, i_e, eta_e_av)
 
-        # Set default solver to DAE
-        self.default_solver = pybamm.ScikitsDaeSolver()
+    @property
+    def default_solver(self):
+        """
+        Create and return the default solver for this model
+        """
+        return pybamm.ScikitsDaeSolver()
 
 
 class MacInnesCapacitance(ElectrolyteCurrentBaseModel):
@@ -355,11 +359,19 @@ class MacInnesCapacitance(ElectrolyteCurrentBaseModel):
     def __init__(self, set_of_parameters, use_capacitance=True):
         super().__init__(set_of_parameters)
         self._use_capacitance = use_capacitance
+
+    @property
+    def default_solver(self):
+        """
+        Create and return the default solver for this model
+        """
         # Different solver depending on whether we solve ODEs or DAEs
-        if use_capacitance:
-            self.default_solver = pybamm.ScikitsOdeSolver()
+        if self._use_capacitance:
+            default_solver = pybamm.ScikitsOdeSolver()
         else:
-            self.default_solver = pybamm.ScikitsDaeSolver()
+            default_solver = pybamm.ScikitsDaeSolver()
+
+        return default_solver
 
     def set_full_system(self, delta_phi, c_e, reactions, eps=None):
         param = self.set_of_parameters
diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py
index 2c60ca360..1690ffb09 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -3,7 +3,7 @@
 #
 import pybamm
 
-from scipy.sparse import diags, eye, kron, csr_matrix, vstack
+from scipy.sparse import diags, eye, kron, csr_matrix, vstack, coo_matrix
 import autograd.numpy as np
 from autograd.builtins import isinstance
 
@@ -284,7 +284,7 @@ class FiniteVolume(pybamm.SpatialMethod):
 
     def add_ghost_nodes(self, symbol, discretised_symbol, bcs):
         """
-        Add boundary conditions via ghost nodes.
+        Add ghost nodes to a symbol.
 
         For Dirichlet bcs, for a boundary condition "y = a at the left-hand boundary",
         we concatenate a ghost node to the start of the vector y with value "2*a - y1"
@@ -301,19 +301,10 @@ class FiniteVolume(pybamm.SpatialMethod):
         where y1 is the value of the first node and h is the mesh size.
         Similarly for the right-hand boundary condition.
 
-        Currently, boundary conditions can only be applied on state
-        variables (e.g. concentration, temperature) or concatenations, and not on
-        arbitrary expressions.
-        To access the value of the first node (y1), we create a "first_node" object
-        which is a StateVector whose y_slice is the start of the y_slice of
-        discretised_symbol.
-        Similarly, the last node is a StateVector whose y_slice is the end of the
-        y_slice of discretised_symbol
-
         Parameters
         ----------
-        discretised_symbol : :class:`pybamm.StateVector` (size n)
-            The discretised variable (a state vector) to which to add ghost nodes
+        domain : list of strings
+            The domain of the symbol for which to add ghost nodes
         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
@@ -321,54 +312,26 @@ class FiniteVolume(pybamm.SpatialMethod):
 
         Returns
         -------
-        :class:`pybamm.Concatenation` (size n+1 or n+2)
-            Concatenation of the variable (a state vector) and ghost nodes
+        :class:`pybamm.Symbol` (shape (n+2, n))
+            `Matrix @ discretised_symbol + bcs_vector`. When evaluated, this gives the
+            discretised_symbol, with appropriate ghost nodes concatenated at each end.
 
         """
-        if isinstance(discretised_symbol, pybamm.StateVector):
-            y_slice_start = discretised_symbol.y_slice.start
-            y_slice_stop = discretised_symbol.y_slice.stop
-        elif isinstance(discretised_symbol, pybamm.Concatenation):
-            y_slice_start = discretised_symbol.children[0].y_slice.start
-            y_slice_stop = discretised_symbol.children[-1].y_slice.stop
-        else:
-            raise TypeError(
-                """
-                discretised_symbol must be a StateVector or Concatenation, not '{}'
-                """.format(
-                    type(discretised_symbol)
-                )
-            )
-        y = np.arange(y_slice_start, y_slice_stop)
-
-        # reshape y_slices into more helpful form
+        # get relevant grid points
         submesh_list = self.mesh.combine_submeshes(*symbol.domain)
         if isinstance(submesh_list[0].npts, list):
             NotImplementedError("Can only take in 1D primary directions")
 
-        size = [len(submesh_list), submesh_list[0].npts]
-        y = np.reshape(y, size)
-        y_left = y[:, 0]
-        y_right = y[:, -1]
+        # Prepare sizes and empty bcs_vector
+        n = submesh_list[0].npts
+        sec_pts = len(submesh_list)
 
-        new_discretised_symbol = pybamm.Vector(np.array([]))  # starts empty
+        bcs_vector = pybamm.Vector(np.array([]))  # starts empty
 
         lbc_value, lbc_type = bcs["left"]
         rbc_value, rbc_type = bcs["right"]
 
-        for i in range(len(submesh_list)):
-            y_slice_start = y_left[i]
-            y_slice_stop = y_right[i]
-
-            # left ghost cell
-            first_node = pybamm.StateVector(slice(y_slice_start, y_slice_start + 1))
-
-            # middle symbol
-            sub_disc_symbol = pybamm.StateVector(slice(y_slice_start, y_slice_stop + 1))
-
-            # right ghost cell
-            last_node = pybamm.StateVector(slice(y_slice_stop, y_slice_stop + 1))
-
+        for i in range(sec_pts):
             if lbc_value.evaluates_to_number():
                 lbc_i = lbc_value
             else:
@@ -379,10 +342,10 @@ class FiniteVolume(pybamm.SpatialMethod):
                 rbc_i = pybamm.Index(rbc_value, i)
 
             if lbc_type == "Dirichlet":
-                left_ghost_cell = 2 * lbc_i - first_node
+                left_ghost_constant = 2 * lbc_i
             elif lbc_type == "Neumann":
                 dx = 2 * (submesh_list[0].nodes[0] - submesh_list[0].edges[0])
-                left_ghost_cell = first_node - dx * lbc_i
+                left_ghost_constant = -dx * lbc_i
             else:
                 raise ValueError(
                     "boundary condition must be Dirichlet or Neumann, not '{}'".format(
@@ -390,10 +353,10 @@ class FiniteVolume(pybamm.SpatialMethod):
                     )
                 )
             if rbc_type == "Dirichlet":
-                right_ghost_cell = 2 * rbc_i - last_node
+                right_ghost_constant = 2 * rbc_i
             elif rbc_type == "Neumann":
                 dx = 2 * (submesh_list[0].edges[-1] - submesh_list[0].nodes[-1])
-                right_ghost_cell = last_node + dx * rbc_i
+                right_ghost_constant = dx * rbc_i
             else:
                 raise ValueError(
                     "boundary condition must be Dirichlet or Neumann, not '{}'".format(
@@ -401,17 +364,32 @@ class FiniteVolume(pybamm.SpatialMethod):
                     )
                 )
             # concatenate
-            concatenated_sub_disc_symbol = pybamm.NumpyConcatenation(
-                left_ghost_cell, sub_disc_symbol, right_ghost_cell
-            )
-            new_discretised_symbol = pybamm.NumpyConcatenation(
-                new_discretised_symbol, concatenated_sub_disc_symbol
+            bcs_vector = pybamm.NumpyConcatenation(
+                bcs_vector,
+                left_ghost_constant,
+                pybamm.Vector(np.zeros(n)),
+                right_ghost_constant,
             )
 
-        # Keep same domain
-        new_discretised_symbol.domain = discretised_symbol.domain
+        # Make matrix to calculate ghost nodes
+        if lbc_type == "Dirichlet":
+            left_factor = -1
+        else:
+            left_factor = 1
+        if rbc_type == "Dirichlet":
+            right_factor = -1
+        else:
+            right_factor = 1
+        # coo_matrix takes inputs (data, (row, col)) and puts data[i] at the point
+        # (row[i], col[i]) for each index of data.
+        left_ghost_vector = coo_matrix(([left_factor], ([0], [0])), shape=(1, n))
+        right_ghost_vector = coo_matrix(([right_factor], ([0], [n - 1])), shape=(1, n))
+        sub_matrix = vstack([left_ghost_vector, eye(n), right_ghost_vector])
+
+        # repeat matrix for secondary dimensions
+        matrix = kron(eye(sec_pts), sub_matrix)
 
-        return new_discretised_symbol
+        return pybamm.Matrix(matrix) @ discretised_symbol + bcs_vector
 
     def boundary_value_or_flux(self, symbol, discretised_child):
         """
@@ -454,7 +432,11 @@ class FiniteVolume(pybamm.SpatialMethod):
                 )
 
         # Generate full matrix from the submatrix
-        matrix = kron(eye(sec_pts), sub_matrix)
+        # Convert to csr_matrix so that we can take the index (row-slicing), which is
+        # not supported by the default kron format
+        # Note that this makes column-slicing inefficient, but this should not be an
+        # issue
+        matrix = csr_matrix(kron(eye(sec_pts), sub_matrix))
 
         # Return boundary value with domain removed
         boundary_value = pybamm.Matrix(matrix) @ discretised_child


Test patch:
diff --git a/tests/integration/test_models/standard_model_tests.py b/tests/integration/test_models/standard_model_tests.py
index eb968124b..a2291dcf2 100644
--- a/tests/integration/test_models/standard_model_tests.py
+++ b/tests/integration/test_models/standard_model_tests.py
@@ -131,8 +131,8 @@ class OptimisationsTest(object):
             parameter_values = model.default_parameter_values
         # Process model and geometry
         parameter_values.process_model(model)
-        parameter_values.process_geometry(model.default_geometry)
         geometry = model.default_geometry
+        parameter_values.process_geometry(geometry)
         # Set discretisation
         if disc is None:
             mesh = pybamm.Mesh(
diff --git a/tests/integration/test_models/standard_output_tests.py b/tests/integration/test_models/standard_output_tests.py
index 000830612..8222e8b22 100644
--- a/tests/integration/test_models/standard_output_tests.py
+++ b/tests/integration/test_models/standard_output_tests.py
@@ -51,10 +51,10 @@ class StandardOutputTests(object):
 class BaseOutputTest(object):
     def __init__(self, model, param, disc, solver, operating_condition):
         self.model = model
+        self.param = param
         self.disc = disc
         self.solver = solver
         self.operating_condition = operating_condition
-
         self.t = solver.t
 
         self.x_n = disc.mesh["negative electrode"][0].nodes
@@ -476,8 +476,7 @@ class CurrentTests(BaseOutputTest):
         t, x_n, x_p = self.t, self.x_n_edge, self.x_p_edge
 
         current_param = pybamm.electrical_parameters.current_with_time
-        parameter_values = self.model.default_parameter_values
-        i_cell = parameter_values.process_symbol(current_param).evaluate(t=t)
+        i_cell = self.param.process_symbol(current_param).evaluate(t=t)
         np.testing.assert_array_almost_equal(self.i_s_n(t, x_n[0]), i_cell, decimal=4)
         np.testing.assert_array_almost_equal(self.i_s_n(t, x_n[-1]), 0, decimal=4)
         np.testing.assert_array_almost_equal(self.i_s_p(t, x_p[-1]), i_cell, decimal=4)
diff --git a/tests/unit/test_expression_tree/test_jac.py b/tests/unit/test_expression_tree/test_jac.py
index 7015300d4..8d7ee590f 100644
--- a/tests/unit/test_expression_tree/test_jac.py
+++ b/tests/unit/test_expression_tree/test_jac.py
@@ -134,6 +134,12 @@ class TestJacobian(unittest.TestCase):
         dfunc_dy = func.jac(y).evaluate(y=y0)
         np.testing.assert_array_equal(jacobian, dfunc_dy.toarray())
 
+    def test_index(self):
+        vec = pybamm.StateVector(slice(0, 5))
+        ind = pybamm.Index(vec, 3)
+        jac = ind.jac(vec).evaluate(y=np.linspace(0, 2, 5)).toarray()
+        np.testing.assert_array_equal(jac, np.array([[0, 0, 0, 1, 0]]))
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_expression_tree/test_unary_operators.py b/tests/unit/test_expression_tree/test_unary_operators.py
index 27c923645..2ea52b63d 100644
--- a/tests/unit/test_expression_tree/test_unary_operators.py
+++ b/tests/unit/test_expression_tree/test_unary_operators.py
@@ -110,8 +110,17 @@ class TestUnaryOperators(unittest.TestCase):
 
     def test_index(self):
         vec = pybamm.Vector(np.array([1, 2, 3, 4, 5]))
+        # with integer
         ind = pybamm.Index(vec, 3)
         self.assertEqual(ind.evaluate(), 4)
+        # with slice
+        ind = pybamm.Index(vec, slice(1, 3))
+        np.testing.assert_array_equal(ind.evaluate(), np.array([[2], [3]]))
+        # errors
+        with self.assertRaisesRegex(TypeError, "index must be integer or slice"):
+            pybamm.Index(vec, 0.0)
+        with self.assertRaisesRegex(ValueError, "slice size exceeds child size"):
+            pybamm.Index(vec, 5)
 
     def test_diff(self):
         a = pybamm.StateVector(slice(0, 1))
diff --git a/tests/unit/test_models/test_base_model.py b/tests/unit/test_models/test_base_model.py
index bdee45f3f..c94116b0a 100644
--- a/tests/unit/test_models/test_base_model.py
+++ b/tests/unit/test_models/test_base_model.py
@@ -331,6 +331,22 @@ class TestBaseModel(unittest.TestCase):
             self.assertIsNotNone(item)
 
 
+class TestStandardBatteryBaseModel(unittest.TestCase):
+    def test_default_solver(self):
+        model = pybamm.StandardBatteryBaseModel()
+        self.assertIsInstance(
+            model.default_solver, (pybamm.ScipySolver, pybamm.ScikitsOdeSolver)
+        )
+
+        # check that default_solver gives you a new solver, not an internal object
+        solver = model.default_solver
+        solver = pybamm.BaseModel()
+        self.assertIsInstance(
+            model.default_solver, (pybamm.ScipySolver, pybamm.ScikitsOdeSolver)
+        )
+        self.assertIsInstance(solver, pybamm.BaseModel)
+
+
 if __name__ == "__main__":
     print("Add -v for more debug output")
     import sys
diff --git a/tests/unit/test_models/test_lead_acid/test_composite_capacitance.py b/tests/unit/test_models/test_lead_acid/test_composite_capacitance.py
index ef2f5d5ed..e1becfbbe 100644
--- a/tests/unit/test_models/test_lead_acid/test_composite_capacitance.py
+++ b/tests/unit/test_models/test_lead_acid/test_composite_capacitance.py
@@ -14,6 +14,12 @@ class TestLeadAcidCompositeCapacitance(unittest.TestCase):
         model = pybamm.lead_acid.CompositeCapacitance(use_capacitance=False)
         model.check_well_posedness()
 
+    def test_default_solver(self):
+        model = pybamm.lead_acid.CompositeCapacitance(use_capacitance=True)
+        self.assertIsInstance(model.default_solver, pybamm.ScikitsOdeSolver)
+        model = pybamm.lead_acid.CompositeCapacitance(use_capacitance=False)
+        self.assertIsInstance(model.default_solver, pybamm.ScikitsDaeSolver)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_models/test_lead_acid/test_loqs.py b/tests/unit/test_models/test_lead_acid/test_loqs.py
index 07b6f1abc..1292687cc 100644
--- a/tests/unit/test_models/test_lead_acid/test_loqs.py
+++ b/tests/unit/test_models/test_lead_acid/test_loqs.py
@@ -10,6 +10,16 @@ class TestLeadAcidLOQS(unittest.TestCase):
         model = pybamm.lead_acid.LOQS()
         model.check_well_posedness()
 
+    def test_default_geometry(self):
+        model = pybamm.lead_acid.LOQS()
+        self.assertIsInstance(model.default_geometry, pybamm.Geometry)
+        self.assertTrue("negative particle" not in model.default_geometry)
+
+    def test_default_spatial_methods(self):
+        model = pybamm.lead_acid.LOQS()
+        self.assertIsInstance(model.default_spatial_methods, dict)
+        self.assertTrue("negative particle" not in model.default_geometry)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_models/test_lead_acid/test_loqs_capacitance.py b/tests/unit/test_models/test_lead_acid/test_loqs_capacitance.py
index 73e915eed..1c9713c14 100644
--- a/tests/unit/test_models/test_lead_acid/test_loqs_capacitance.py
+++ b/tests/unit/test_models/test_lead_acid/test_loqs_capacitance.py
@@ -14,6 +14,12 @@ class TestLeadAcidLOQSCapacitance(unittest.TestCase):
         model = pybamm.lead_acid.LOQSCapacitance(use_capacitance=False)
         model.check_well_posedness()
 
+    def test_default_solver(self):
+        model = pybamm.lead_acid.LOQSCapacitance(use_capacitance=True)
+        self.assertIsInstance(model.default_solver, pybamm.ScikitsOdeSolver)
+        model = pybamm.lead_acid.LOQSCapacitance(use_capacitance=False)
+        self.assertIsInstance(model.default_solver, pybamm.ScikitsDaeSolver)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_models/test_lead_acid/test_newman_tiedemann.py b/tests/unit/test_models/test_lead_acid/test_newman_tiedemann.py
index a30d77bd4..c95e5b02e 100644
--- a/tests/unit/test_models/test_lead_acid/test_newman_tiedemann.py
+++ b/tests/unit/test_models/test_lead_acid/test_newman_tiedemann.py
@@ -10,6 +10,10 @@ class TestLeadAcidNewmanTiedemann(unittest.TestCase):
         model = pybamm.lead_acid.NewmanTiedemann()
         model.check_well_posedness()
 
+    def test_default_solver(self):
+        model = pybamm.lead_acid.NewmanTiedemann()
+        self.assertIsInstance(model.default_solver, pybamm.ScikitsDaeSolver)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_models/test_lead_acid/test_newman_tiedemann_capacitance.py b/tests/unit/test_models/test_lead_acid/test_newman_tiedemann_capacitance.py
index d725b32d5..a60457092 100644
--- a/tests/unit/test_models/test_lead_acid/test_newman_tiedemann_capacitance.py
+++ b/tests/unit/test_models/test_lead_acid/test_newman_tiedemann_capacitance.py
@@ -14,6 +14,12 @@ class TestLeadAcidNewmanTiedemannCapacitance(unittest.TestCase):
         model = pybamm.lead_acid.NewmanTiedemannCapacitance(use_capacitance=False)
         model.check_well_posedness()
 
+    def test_default_solver(self):
+        model = pybamm.lead_acid.NewmanTiedemannCapacitance(use_capacitance=True)
+        self.assertIsInstance(model.default_solver, pybamm.ScikitsOdeSolver)
+        model = pybamm.lead_acid.NewmanTiedemannCapacitance(use_capacitance=False)
+        self.assertIsInstance(model.default_solver, pybamm.ScikitsDaeSolver)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_models/test_lithium_ion/test_lithium_ion_dfn.py b/tests/unit/test_models/test_lithium_ion/test_lithium_ion_dfn.py
index 862bf9245..f63605187 100644
--- a/tests/unit/test_models/test_lithium_ion/test_lithium_ion_dfn.py
+++ b/tests/unit/test_models/test_lithium_ion/test_lithium_ion_dfn.py
@@ -10,6 +10,15 @@ class TestDFN(unittest.TestCase):
         model = pybamm.lithium_ion.DFN()
         model.check_well_posedness()
 
+    def test_default_geometry(self):
+        model = pybamm.lithium_ion.DFN()
+        self.assertIsInstance(model.default_geometry, pybamm.Geometry)
+        self.assertTrue("secondary" in model.default_geometry["negative particle"])
+
+    def test_default_solver(self):
+        model = pybamm.lithium_ion.DFN()
+        self.assertIsInstance(model.default_solver, pybamm.ScikitsDaeSolver)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_models/test_lithium_ion/test_lithium_ion_spm.py b/tests/unit/test_models/test_lithium_ion/test_lithium_ion_spm.py
index 951e93c93..b529c28fd 100644
--- a/tests/unit/test_models/test_lithium_ion/test_lithium_ion_spm.py
+++ b/tests/unit/test_models/test_lithium_ion/test_lithium_ion_spm.py
@@ -10,6 +10,11 @@ class TestSPM(unittest.TestCase):
         model = pybamm.lithium_ion.SPM()
         model.check_well_posedness()
 
+    def test_default_geometry(self):
+        model = pybamm.lithium_ion.SPM()
+        self.assertIsInstance(model.default_geometry, pybamm.Geometry)
+        self.assertTrue("negative particle" in model.default_geometry)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_models/test_lithium_ion/test_lithium_ion_spme.py b/tests/unit/test_models/test_lithium_ion/test_lithium_ion_spme.py
index 9b2424ee2..c050e7d9e 100644
--- a/tests/unit/test_models/test_lithium_ion/test_lithium_ion_spme.py
+++ b/tests/unit/test_models/test_lithium_ion/test_lithium_ion_spme.py
@@ -10,6 +10,11 @@ class TestSPMe(unittest.TestCase):
         model = pybamm.lithium_ion.SPMe()
         model.check_well_posedness()
 
+    def test_default_geometry(self):
+        model = pybamm.lithium_ion.SPMe()
+        self.assertIsInstance(model.default_geometry, pybamm.Geometry)
+        self.assertTrue("negative particle" in model.default_geometry)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_models/test_reaction_diffusion.py b/tests/unit/test_models/test_reaction_diffusion.py
index a0922e200..1ecbb3bc8 100644
--- a/tests/unit/test_models/test_reaction_diffusion.py
+++ b/tests/unit/test_models/test_reaction_diffusion.py
@@ -10,6 +10,16 @@ class TestReactionDiffusionModel(unittest.TestCase):
         model = pybamm.ReactionDiffusionModel()
         model.check_well_posedness()
 
+    def test_default_parameter_values(self):
+        model = pybamm.ReactionDiffusionModel()
+        self.assertTrue(
+            isinstance(model.default_parameter_values, pybamm.ParameterValues)
+        )
+        self.assertTrue(
+            "Partial molar volume of cations [m3.mol-1]"
+            in model.default_parameter_values
+        )
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_models/test_submodels/__init__.py b/tests/unit/test_models/test_submodels/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/tests/unit/test_models/test_submodels/test_electrode.py b/tests/unit/test_models/test_submodels/test_electrode.py
new file mode 100644
index 000000000..00746fe94
--- /dev/null
+++ b/tests/unit/test_models/test_submodels/test_electrode.py
@@ -0,0 +1,23 @@
+#
+# Tests for the electrode submodels
+#
+import pybamm
+from pybamm.solvers.scikits_ode_solver import scikits_odes_spec
+import unittest
+
+
+@unittest.skipIf(scikits_odes_spec is None, "scikits.odes not installed")
+class TestOhm(unittest.TestCase):
+    def test_default_solver(self):
+        param = pybamm.standard_parameters_lithium_ion
+        model = pybamm.electrode.Ohm(param)
+        self.assertIsInstance(model.default_solver, pybamm.ScikitsDaeSolver)
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    unittest.main()
diff --git a/tests/unit/test_models/test_submodels/test_electrolyte_current.py b/tests/unit/test_models/test_submodels/test_electrolyte_current.py
new file mode 100644
index 000000000..a74eb0c25
--- /dev/null
+++ b/tests/unit/test_models/test_submodels/test_electrolyte_current.py
@@ -0,0 +1,33 @@
+#
+# Tests for the electrolyte current submodel
+#
+import pybamm
+from pybamm.solvers.scikits_ode_solver import scikits_odes_spec
+import unittest
+
+
+@unittest.skipIf(scikits_odes_spec is None, "scikits.odes not installed")
+class TestMacInnesStefanMaxwell(unittest.TestCase):
+    def test_default_solver(self):
+        param = pybamm.standard_parameters_lithium_ion
+        model = pybamm.electrolyte_current.MacInnesStefanMaxwell(param)
+        self.assertIsInstance(model.default_solver, pybamm.ScikitsDaeSolver)
+
+
+@unittest.skipIf(scikits_odes_spec is None, "scikits.odes not installed")
+class TestMacInnesCapacitance(unittest.TestCase):
+    def test_default_solver(self):
+        param = pybamm.standard_parameters_lithium_ion
+        model = pybamm.electrolyte_current.MacInnesCapacitance(param, True)
+        self.assertIsInstance(model.default_solver, pybamm.ScikitsOdeSolver)
+        model = pybamm.electrolyte_current.MacInnesCapacitance(param, False)
+        self.assertIsInstance(model.default_solver, pybamm.ScikitsDaeSolver)
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    unittest.main()
diff --git a/tests/unit/test_spatial_methods/test_finite_volume.py b/tests/unit/test_spatial_methods/test_finite_volume.py
index fa503a1e7..ee2a4251e 100644
--- a/tests/unit/test_spatial_methods/test_finite_volume.py
+++ b/tests/unit/test_spatial_methods/test_finite_volume.py
@@ -195,10 +195,6 @@ class TestFiniteVolume(unittest.TestCase):
         )
 
         # test errors
-        with self.assertRaisesRegex(
-            TypeError, "discretised_symbol must be a StateVector or Concatenation"
-        ):
-            pybamm.FiniteVolume(mesh).add_ghost_nodes(var, pybamm.Scalar(1), None)
         bcs = {"left": (pybamm.Scalar(0), "x"), "right": (pybamm.Scalar(3), "Neumann")}
         with self.assertRaisesRegex(ValueError, "boundary condition must be"):
             pybamm.FiniteVolume(mesh).add_ghost_nodes(var, discretised_symbol, bcs)
@@ -493,24 +489,25 @@ class TestFiniteVolume(unittest.TestCase):
 
         np.testing.assert_array_equal(grad_eval, np.zeros([sec_pts, prim_pts + 1]))
 
-        # div: test on linear r^2
-        # div (grad r^2) = 6
-        const = 6 * np.ones(sec_pts * prim_pts)
-
+        # div
+        # div (grad r^2) = 6, N_left = N_right = 0
         N = pybamm.grad(var)
         div_eqn = pybamm.div(N)
+        bc_var = disc.process_symbol(
+            pybamm.SpatialVariable("x_n", domain="negative electrode")
+        )
         boundary_conditions = {
-            var.id: {
-                "left": (pybamm.Scalar(6), "Dirichlet"),
-                "right": (pybamm.Scalar(6), "Dirichlet"),
-            }
+            var.id: {"left": (bc_var, "Neumann"), "right": (bc_var, "Neumann")}
         }
         disc.bcs = boundary_conditions
-
         div_eqn_disc = disc.process_symbol(div_eqn)
+
+        const = 6 * np.ones(sec_pts * prim_pts)
         div_eval = div_eqn_disc.evaluate(None, const)
         div_eval = np.reshape(div_eval, [sec_pts, prim_pts])
-        np.testing.assert_array_almost_equal(div_eval, np.zeros([sec_pts, prim_pts]))
+        np.testing.assert_array_almost_equal(
+            div_eval[:, :-1], np.zeros([sec_pts, prim_pts - 1])
+        )
 
     def test_grad_div_shapes_Neumann_bcs(self):
         """Test grad and div with Neumann boundary conditions (applied by div on N)"""
@@ -959,6 +956,7 @@ class TestFiniteVolume(unittest.TestCase):
         grad_and_div = pybamm.div(pybamm.grad(phi))
         int_grad_phi = pybamm.IndefiniteIntegral(grad_and_div, x)
         disc.set_variable_slices([phi])
+        # Add boundary conditions so that the integrand has the right shape
         disc._bcs = {
             phi.id: {
                 "left": (pybamm.Scalar(0), "Dirichlet"),