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

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.
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"),

pip install -e .[dev,docs]
parse_log_pytest
pytest
flake8 sphinx pytest
apt-get update apt-get install -y gcc
3.6
requirements.txt
pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
alabaster==0.7.13
anyio==3.6.2
anytree==2.8.0
argon2-cffi==21.3.0
argon2-cffi-bindings==21.2.0
async-generator==1.10
attrs @ file:///opt/conda/conda-bld/attrs_1642510447205/work
autograd==1.6.2
Babel==2.11.0
backcall==0.2.0
bleach==4.1.0
certifi==2021.5.30
cffi==1.15.1
charset-normalizer==2.0.12
comm==0.1.4
contextvars==2.4
cycler==0.11.0
dataclasses==0.8
decorator==5.1.1
defusedxml==0.7.1
docutils==0.17.1
entrypoints==0.4
flake8==5.0.4
future==1.0.0
guzzle-sphinx-theme==0.7.11
idna==3.10
imagesize==1.4.1
immutables==0.19
importlib-metadata==4.2.0
iniconfig @ file:///home/linux1/recipes/ci/iniconfig_1610983019677/work
ipykernel==5.5.6
ipython==7.16.3
ipython-genutils==0.2.0
ipywidgets==7.8.5
jedi==0.17.2
Jinja2==3.0.3
json5==0.9.16
jsonschema==3.2.0
jupyter==1.1.1
jupyter-client==7.1.2
jupyter-console==6.4.3
jupyter-core==4.9.2
jupyter-server==1.13.1
jupyterlab==3.2.9
jupyterlab-pygments==0.1.2
jupyterlab-server==2.10.3
jupyterlab_widgets==1.1.11
kiwisolver==1.3.1
MarkupSafe==2.0.1
matplotlib==3.3.4
mccabe==0.7.0
mistune==0.8.4
more-itertools @ file:///tmp/build/80754af9/more-itertools_1637733554872/work
nbclassic==0.3.5
nbclient==0.5.9
nbconvert==6.0.7
nbformat==5.1.3
nest-asyncio==1.6.0
notebook==6.4.10
numpy==1.19.5
packaging @ file:///tmp/build/80754af9/packaging_1637314298585/work
pandas==1.1.5
pandocfilters==1.5.1
parso==0.7.1
pexpect==4.9.0
pickleshare==0.7.5
Pillow==8.4.0
pluggy @ file:///tmp/build/80754af9/pluggy_1615976315926/work
prometheus-client==0.17.1
prompt-toolkit==3.0.36
ptyprocess==0.7.0
py @ file:///opt/conda/conda-bld/py_1644396412707/work
-e git+https://github.com/pybamm-team/PyBaMM.git@ecbfa8446d7a770f03ad83ad47d6b40263e5fa54#egg=pybamm
pycodestyle==2.9.1
pycparser==2.21
pyflakes==2.5.0
Pygments==2.14.0
pyparsing @ file:///tmp/build/80754af9/pyparsing_1635766073266/work
pyrsistent==0.18.0
pytest==6.2.4
python-dateutil==2.9.0.post0
pytz==2025.2
pyzmq==25.1.2
requests==2.27.1
scipy==1.5.4
Send2Trash==1.8.3
six==1.17.0
sniffio==1.2.0
snowballstemmer==2.2.0
Sphinx==4.3.2
sphinxcontrib-applehelp==1.0.2
sphinxcontrib-devhelp==1.0.2
sphinxcontrib-htmlhelp==2.0.0
sphinxcontrib-jsmath==1.0.1
sphinxcontrib-qthelp==1.0.3
sphinxcontrib-serializinghtml==1.1.5
terminado==0.12.1
testpath==0.6.0
toml @ file:///tmp/build/80754af9/toml_1616166611790/work
tornado==6.1
traitlets==4.3.3
typing_extensions @ file:///opt/conda/conda-bld/typing_extensions_1647553014482/work
urllib3==1.26.20
wcwidth==0.2.13
webencodings==0.5.1
websocket-client==1.3.1
widgetsnbextension==3.6.10
zipp @ file:///tmp/build/80754af9/zipp_1633618647012/work

name: PyBaMM
channels:
  - defaults
  - https://repo.anaconda.com/pkgs/main
  - https://repo.anaconda.com/pkgs/r
  - conda-forge
dependencies:
  - _libgcc_mutex=0.1=main
  - _openmp_mutex=5.1=1_gnu
  - attrs=21.4.0=pyhd3eb1b0_0
  - ca-certificates=2025.2.25=h06a4308_0
  - certifi=2021.5.30=py36h06a4308_0
  - importlib_metadata=4.8.1=hd3eb1b0_0
  - iniconfig=1.1.1=pyhd3eb1b0_0
  - ld_impl_linux-64=2.40=h12ee557_0
  - libffi=3.3=he6710b0_2
  - libgcc-ng=11.2.0=h1234567_1
  - libgomp=11.2.0=h1234567_1
  - libstdcxx-ng=11.2.0=h1234567_1
  - more-itertools=8.12.0=pyhd3eb1b0_0
  - ncurses=6.4=h6a678d5_0
  - openssl=1.1.1w=h7f8727e_0
  - packaging=21.3=pyhd3eb1b0_0
  - pip=21.2.2=py36h06a4308_0
  - pluggy=0.13.1=py36h06a4308_0
  - py=1.11.0=pyhd3eb1b0_0
  - pyparsing=3.0.4=pyhd3eb1b0_0
  - pytest=6.2.4=py36h06a4308_2
  - python=3.6.13=h12debd9_1
  - readline=8.2=h5eee18b_0
  - setuptools=58.0.4=py36h06a4308_0
  - sqlite=3.45.3=h5eee18b_0
  - tk=8.6.14=h39e8969_0
  - toml=0.10.2=pyhd3eb1b0_0
  - typing_extensions=4.1.1=pyh06a4308_0
  - wheel=0.37.1=pyhd3eb1b0_0
  - xz=5.6.4=h5eee18b_1
  - zipp=3.6.0=pyhd3eb1b0_0
  - zlib=1.2.13=h5eee18b_1
  - pip:
      - alabaster==0.7.13
      - anyio==3.6.2
      - anytree==2.8.0
      - argon2-cffi==21.3.0
      - argon2-cffi-bindings==21.2.0
      - async-generator==1.10
      - autograd==1.6.2
      - babel==2.11.0
      - backcall==0.2.0
      - bleach==4.1.0
      - cffi==1.15.1
      - charset-normalizer==2.0.12
      - comm==0.1.4
      - contextvars==2.4
      - cycler==0.11.0
      - dataclasses==0.8
      - decorator==5.1.1
      - defusedxml==0.7.1
      - docutils==0.17.1
      - entrypoints==0.4
      - flake8==5.0.4
      - future==1.0.0
      - guzzle-sphinx-theme==0.7.11
      - idna==3.10
      - imagesize==1.4.1
      - immutables==0.19
      - importlib-metadata==4.2.0
      - ipykernel==5.5.6
      - ipython==7.16.3
      - ipython-genutils==0.2.0
      - ipywidgets==7.8.5
      - jedi==0.17.2
      - jinja2==3.0.3
      - json5==0.9.16
      - jsonschema==3.2.0
      - jupyter==1.1.1
      - jupyter-client==7.1.2
      - jupyter-console==6.4.3
      - jupyter-core==4.9.2
      - jupyter-server==1.13.1
      - jupyterlab==3.2.9
      - jupyterlab-pygments==0.1.2
      - jupyterlab-server==2.10.3
      - jupyterlab-widgets==1.1.11
      - kiwisolver==1.3.1
      - markupsafe==2.0.1
      - matplotlib==3.3.4
      - mccabe==0.7.0
      - mistune==0.8.4
      - nbclassic==0.3.5
      - nbclient==0.5.9
      - nbconvert==6.0.7
      - nbformat==5.1.3
      - nest-asyncio==1.6.0
      - notebook==6.4.10
      - numpy==1.19.5
      - pandas==1.1.5
      - pandocfilters==1.5.1
      - parso==0.7.1
      - pexpect==4.9.0
      - pickleshare==0.7.5
      - pillow==8.4.0
      - prometheus-client==0.17.1
      - prompt-toolkit==3.0.36
      - ptyprocess==0.7.0
      - pycodestyle==2.9.1
      - pycparser==2.21
      - pyflakes==2.5.0
      - pygments==2.14.0
      - pyrsistent==0.18.0
      - python-dateutil==2.9.0.post0
      - pytz==2025.2
      - pyzmq==25.1.2
      - requests==2.27.1
      - scipy==1.5.4
      - send2trash==1.8.3
      - six==1.17.0
      - sniffio==1.2.0
      - snowballstemmer==2.2.0
      - sphinx==4.3.2
      - sphinxcontrib-applehelp==1.0.2
      - sphinxcontrib-devhelp==1.0.2
      - sphinxcontrib-htmlhelp==2.0.0
      - sphinxcontrib-jsmath==1.0.1
      - sphinxcontrib-qthelp==1.0.3
      - sphinxcontrib-serializinghtml==1.1.5
      - terminado==0.12.1
      - testpath==0.6.0
      - tornado==6.1
      - traitlets==4.3.3
      - urllib3==1.26.20
      - wcwidth==0.2.13
      - webencodings==0.5.1
      - websocket-client==1.3.1
      - widgetsnbextension==3.6.10
prefix: /opt/conda/envs/PyBaMM

tests/unit/test_expression_tree/test_jac.py::TestJacobian::test_index tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_index tests/unit/test_models/test_lead_acid/test_composite_capacitance.py::TestLeadAcidCompositeCapacitance::test_well_posed tests/unit/test_models/test_lead_acid/test_composite_capacitance.py::TestLeadAcidCompositeCapacitance::test_well_posed_no_capacitance tests/unit/test_models/test_lead_acid/test_loqs_capacitance.py::TestLeadAcidLOQSCapacitance::test_well_posed tests/unit/test_models/test_lead_acid/test_loqs_capacitance.py::TestLeadAcidLOQSCapacitance::test_well_posed_no_capacitance tests/unit/test_models/test_lead_acid/test_newman_tiedemann.py::TestLeadAcidNewmanTiedemann::test_well_posed tests/unit/test_models/test_lead_acid/test_newman_tiedemann_capacitance.py::TestLeadAcidNewmanTiedemannCapacitance::test_well_posed tests/unit/test_models/test_lead_acid/test_newman_tiedemann_capacitance.py::TestLeadAcidNewmanTiedemannCapacitance::test_well_posed_no_capacitance tests/unit/test_models/test_lithium_ion/test_lithium_ion_dfn.py::TestDFN::test_default_geometry tests/unit/test_models/test_lithium_ion/test_lithium_ion_dfn.py::TestDFN::test_well_posed
tests/unit/test_models/test_lead_acid/test_composite_capacitance.py::TestLeadAcidCompositeCapacitance::test_default_solver tests/unit/test_models/test_lead_acid/test_loqs_capacitance.py::TestLeadAcidLOQSCapacitance::test_default_solver tests/unit/test_models/test_lead_acid/test_newman_tiedemann.py::TestLeadAcidNewmanTiedemann::test_default_solver tests/unit/test_models/test_lead_acid/test_newman_tiedemann_capacitance.py::TestLeadAcidNewmanTiedemannCapacitance::test_default_solver tests/unit/test_models/test_lithium_ion/test_lithium_ion_dfn.py::TestDFN::test_default_solver
tests/unit/test_expression_tree/test_jac.py::TestJacobian::test_functions tests/unit/test_expression_tree/test_jac.py::TestJacobian::test_linear tests/unit/test_expression_tree/test_jac.py::TestJacobian::test_nonlinear tests/unit/test_expression_tree/test_unary_operators.py::test_const_function tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_absolute tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_average tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_boundary_operators tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_boundary_value tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_diff tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_function tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_gradient tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_id tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_integral tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_negation tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_printing tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_unary_operator tests/unit/test_models/test_base_model.py::TestBaseModel::test_algebraic_set_get tests/unit/test_models/test_base_model.py::TestBaseModel::test_boundary_conditions_set_get tests/unit/test_models/test_base_model.py::TestBaseModel::test_check_well_posedness_initial_boundary_conditions tests/unit/test_models/test_base_model.py::TestBaseModel::test_check_well_posedness_output_variables tests/unit/test_models/test_base_model.py::TestBaseModel::test_check_well_posedness_variables tests/unit/test_models/test_base_model.py::TestBaseModel::test_initial_conditions_set_get tests/unit/test_models/test_base_model.py::TestBaseModel::test_model_dict_behaviour tests/unit/test_models/test_base_model.py::TestBaseModel::test_rhs_set_get tests/unit/test_models/test_base_model.py::TestBaseModel::test_update tests/unit/test_models/test_base_model.py::TestBaseModel::test_variables_set_get tests/unit/test_models/test_base_model.py::TestStandardBatteryBaseModel::test_default_solver tests/unit/test_models/test_lead_acid/test_loqs.py::TestLeadAcidLOQS::test_default_geometry tests/unit/test_models/test_lead_acid/test_loqs.py::TestLeadAcidLOQS::test_default_spatial_methods tests/unit/test_models/test_lead_acid/test_loqs.py::TestLeadAcidLOQS::test_well_posed tests/unit/test_models/test_lithium_ion/test_lithium_ion_spm.py::TestSPM::test_default_geometry tests/unit/test_models/test_lithium_ion/test_lithium_ion_spm.py::TestSPM::test_well_posed tests/unit/test_models/test_lithium_ion/test_lithium_ion_spme.py::TestSPMe::test_default_geometry tests/unit/test_models/test_lithium_ion/test_lithium_ion_spme.py::TestSPMe::test_well_posed tests/unit/test_models/test_reaction_diffusion.py::TestReactionDiffusionModel::test_default_parameter_values tests/unit/test_models/test_reaction_diffusion.py::TestReactionDiffusionModel::test_well_posed tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_add_ghost_nodes tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_add_ghost_nodes_concatenation tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_compute_diffusivity tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_definite_integral tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_discretise_diffusivity_times_spatial_operator tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_discretise_spatial_variable tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_extrapolate_left_right tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_grad_div_shapes_Dirichlet_and_Neumann_bcs tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_grad_div_shapes_Dirichlet_bcs tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_grad_div_shapes_Neumann_bcs tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_grad_div_shapes_mixed_domain tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_indefinite_integral tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_jacobian tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_mass_matrix_shape tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_p2d_add_ghost_nodes tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_p2d_mass_matrix_shape tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_p2d_spherical_grad_div_shapes_Dirichlet_bcs tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_p2d_spherical_grad_div_shapes_Neumann_bcs tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_spherical_grad_div_shapes_Dirichlet_bcs tests/unit/test_spatial_methods/test_finite_volume.py::TestFiniteVolume::test_spherical_grad_div_shapes_Neumann_bcs
pybamm/expression_tree/unary_operators.py pybamm/models/lead_acid/newman_tiedemann.py pybamm/models/lead_acid/composite_capacitance.py pybamm/models/lead_acid/newman_tiedemann_capacitance.py pybamm/models/submodels/electrolyte_current.py pybamm/models/reaction_diffusion.py pybamm/discretisations/discretisation.py pybamm/models/lithium_ion/spme.py pybamm/spatial_methods/finite_volume.py pybamm/models/lithium_ion/spm.py pybamm/models/lead_acid/loqs_capacitance.py pybamm/models/submodels/electrode.py pybamm/models/base_models.py pybamm/models/lead_acid/loqs.py pybamm/models/lithium_ion/dfn.py
pybamm/expression_tree/unary_operators.py pybamm/models/lead_acid/newman_tiedemann.py pybamm/models/lead_acid/composite_capacitance.py pybamm/models/lead_acid/newman_tiedemann_capacitance.py pybamm/models/submodels/electrolyte_current.py pybamm/models/reaction_diffusion.py pybamm/discretisations/discretisation.py pybamm/models/lithium_ion/spme.py pybamm/spatial_methods/finite_volume.py pybamm/models/lithium_ion/spm.py pybamm/models/lead_acid/loqs_capacitance.py pybamm/models/submodels/electrode.py pybamm/models/base_models.py pybamm/models/lead_acid/loqs.py pybamm/models/lithium_ion/dfn.py