Repository: pybamm-team/PyBaMM

Issue:
Create "complete" SPM and DFN classes
Create classes where the whole model can be seen in a single place

Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 834c1502f..8dfec92bd 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## Features
 
+-   Added some basic models (BasicSPM and BasicDFN) in order to clearly demonstrate the PyBaMM model structure for battery models ([#795](https://github.com/pybamm-team/PyBaMM/pull/795))
 -   Added the harmonic mean to the Finite Volume method, which is now used when computing fluxes ([#783](https://github.com/pybamm-team/PyBaMM/pull/783))
 -   Refactored `Solution` to make it a dictionary that contains all of the solution variables. This automatically creates `ProcessedVariable` objects when required, so that the solution can be obtained much more easily. ([#781](https://github.com/pybamm-team/PyBaMM/pull/781))
 -   Added notebook to explain broadcasts ([#776](https://github.com/pybamm-team/PyBaMM/pull/776))
diff --git a/docs/source/models/lithium_ion/dfn.rst b/docs/source/models/lithium_ion/dfn.rst
index 721c86786..9c3032760 100644
--- a/docs/source/models/lithium_ion/dfn.rst
+++ b/docs/source/models/lithium_ion/dfn.rst
@@ -3,3 +3,6 @@ Doyle-Fuller-Newman (DFN)
 
 .. autoclass:: pybamm.lithium_ion.DFN
     :members:
+
+.. autoclass:: pybamm.lithium_ion.BasicDFN
+    :members:
diff --git a/docs/source/models/lithium_ion/spm.rst b/docs/source/models/lithium_ion/spm.rst
index 5bf53a00d..fe891e98b 100644
--- a/docs/source/models/lithium_ion/spm.rst
+++ b/docs/source/models/lithium_ion/spm.rst
@@ -3,3 +3,6 @@ Single Particle Model (SPM)
 
 .. autoclass:: pybamm.lithium_ion.SPM
     :members:
+
+.. autoclass:: pybamm.lithium_ion.BasicSPM
+    :members:
diff --git a/examples/notebooks/README.md b/examples/notebooks/README.md
index 6a90a236e..fca28de81 100644
--- a/examples/notebooks/README.md
+++ b/examples/notebooks/README.md
@@ -48,7 +48,9 @@ The following notebooks are specific to different stages of the PyBaMM pipeline,
 
 ### Models
 
-The following models are implemented and can easily be used or [compared](./models/lead-acid.ipynb). We always welcome [new models](https://pybamm.readthedocs.io/en/latest/tutorials/add-model.html)!
+Several battery models are implemented and can easily be used or [compared](./models/lead-acid.ipynb). The notebooks below show the solution of each individual model. We always welcome [new models](https://pybamm.readthedocs.io/en/latest/tutorials/add-model.html)!
+
+Once you are comfortable with the expression tree structure, a good starting point to understand the models in PyBaMM is to take a look at the [basic SPM](https://github.com/pybamm-team/PyBaMM/blob/master/pybamm/models/full_battery_models/lithium_ion/basic_spm.py) and [basic DFN](https://github.com/pybamm-team/PyBaMM/blob/master/pybamm/models/full_battery_models/lithium_ion/basic_dfn.py), since these define the entire model (variables, equations, initial and boundary conditions, events) in a single class and so are easier to understand. However, we recommend that you subsequently use the full models as they offer much greater flexibility for coupling different physical effects and visualising a greater range of variables.
 
 #### Lithium-ion models
 
diff --git a/pybamm/models/full_battery_models/lithium_ion/__init__.py b/pybamm/models/full_battery_models/lithium_ion/__init__.py
index a807d22fc..167a38f3b 100644
--- a/pybamm/models/full_battery_models/lithium_ion/__init__.py
+++ b/pybamm/models/full_battery_models/lithium_ion/__init__.py
@@ -5,3 +5,5 @@ from .base_lithium_ion_model import BaseModel
 from .spm import SPM
 from .spme import SPMe
 from .dfn import DFN
+from .basic_dfn import BasicDFN
+from .basic_spm import BasicSPM
diff --git a/pybamm/models/full_battery_models/lithium_ion/basic_dfn.py b/pybamm/models/full_battery_models/lithium_ion/basic_dfn.py
new file mode 100644
index 000000000..f395877b0
--- /dev/null
+++ b/pybamm/models/full_battery_models/lithium_ion/basic_dfn.py
@@ -0,0 +1,271 @@
+#
+# Basic Doyle-Fuller-Newman (DFN) Model
+#
+import pybamm
+from .base_lithium_ion_model import BaseModel
+
+
+class BasicDFN(BaseModel):
+    """Doyle-Fuller-Newman (DFN) model of a lithium-ion battery, from [2]_.
+
+    This class differs from the :class:`pybamm.lithium_ion.DFN` model class in that it
+    shows the whole model in a single class. This comes at the cost of flexibility in
+    comparing different physical effects, and in general the main DFN class should be
+    used instead.
+
+    Parameters
+    ----------
+    name : str, optional
+        The name of the model.
+
+    References
+    ----------
+    .. [2] SG Marquis, V Sulzer, R Timms, CP Please and SJ Chapman. “An asymptotic
+           derivation of a single particle model with electrolyte”. In: arXiv preprint
+           arXiv:1905.12553 (2019).
+
+
+    **Extends:** :class:`pybamm.lithium_ion.BaseModel`
+    """
+
+    def __init__(self, name="Doyle-Fuller-Newman model"):
+        super().__init__({}, name)
+        # `param` is a class containing all the relevant parameters and functions for
+        # this model. These are purely symbolic at this stage, and will be set by the
+        # `ParameterValues` class when the model is processed.
+        param = self.param
+
+        ######################
+        # Variables
+        ######################
+        # Variables that depend on time only are created without a domain
+        Q = pybamm.Variable("Discharge capacity [A.h]")
+        # Variables that vary spatially are created with a domain
+        c_e_n = pybamm.Variable(
+            "Negative electrolyte concentration", domain="negative electrode",
+        )
+        c_e_s = pybamm.Variable(
+            "Separator electrolyte concentration", domain="separator",
+        )
+        c_e_p = pybamm.Variable(
+            "Positive electrolyte concentration", domain="positive electrode",
+        )
+        # Concatenations combine several variables into a single variable, to simplify
+        # implementing equations that hold over several domains
+        c_e = pybamm.Concatenation(c_e_n, c_e_s, c_e_p)
+
+        # Electrolyte potential
+        phi_e_n = pybamm.Variable(
+            "Negative electrolyte potential", domain="negative electrode",
+        )
+        phi_e_s = pybamm.Variable(
+            "Separator electrolyte potential", domain="separator",
+        )
+        phi_e_p = pybamm.Variable(
+            "Positive electrolyte potential", domain="positive electrode",
+        )
+        phi_e = pybamm.Concatenation(phi_e_n, phi_e_s, phi_e_p)
+
+        # Electrode potential
+        phi_s_n = pybamm.Variable(
+            "Negative electrode potential", domain="negative electrode",
+        )
+        phi_s_p = pybamm.Variable(
+            "Positive electrode potential", domain="positive electrode",
+        )
+        # Particle concentrations are variables on the particle domain, but also vary in
+        # the x-direction (electrode domain) and so must be provided with auxiliary
+        # domains
+        c_s_n = pybamm.Variable(
+            "Negative particle concentration",
+            domain="negative particle",
+            auxiliary_domains={"secondary": "negative electrode"},
+        )
+        c_s_p = pybamm.Variable(
+            "Positive particle concentration",
+            domain="positive particle",
+            auxiliary_domains={"secondary": "positive electrode"},
+        )
+
+        # Constant temperature
+        T = param.T_init
+
+        ######################
+        # Other set-up
+        ######################
+
+        # Current density
+        i_cell = param.current_with_time
+
+        # Porosity
+        # Primary broadcasts are used to broadcast scalar quantities across a domain
+        # into a vector of the right shape, for multiplying with other vectors
+        eps_n = pybamm.PrimaryBroadcast(param.epsilon_n, "negative electrode")
+        eps_s = pybamm.PrimaryBroadcast(param.epsilon_s, "separator")
+        eps_p = pybamm.PrimaryBroadcast(param.epsilon_p, "positive electrode")
+        eps = pybamm.Concatenation(eps_n, eps_s, eps_p)
+
+        # Tortuosity
+        tor = pybamm.Concatenation(
+            eps_n ** param.b_e_n, eps_s ** param.b_e_s, eps_p ** param.b_e_p
+        )
+
+        # Interfacial reactions
+        # Surf takes the surface value of a variable, i.e. its boundary value on the
+        # right side. This is also accessible via `boundary_value(x, "right")`, with
+        # "left" providing the boundary value of the left side
+        c_s_surf_n = pybamm.surf(c_s_n)
+        j0_n = (
+            param.m_n(T)
+            / param.C_r_n
+            * c_e_n ** (1 / 2)
+            * c_s_surf_n ** (1 / 2)
+            * (1 - c_s_surf_n) ** (1 / 2)
+        )
+        j_n = (
+            2
+            * j0_n
+            * pybamm.sinh(
+                param.ne_n / 2 * (phi_s_n - phi_e_n - param.U_n(c_s_surf_n, T))
+            )
+        )
+        c_s_surf_p = pybamm.surf(c_s_p)
+        j0_p = (
+            param.gamma_p
+            * param.m_p(T)
+            / param.C_r_p
+            * c_e_p ** (1 / 2)
+            * c_s_surf_p ** (1 / 2)
+            * (1 - c_s_surf_p) ** (1 / 2)
+        )
+        j_s = pybamm.PrimaryBroadcast(0, "separator")
+        j_p = (
+            2
+            * j0_p
+            * pybamm.sinh(
+                param.ne_p / 2 * (phi_s_p - phi_e_p - param.U_p(c_s_surf_p, T))
+            )
+        )
+        j = pybamm.Concatenation(j_n, j_s, j_p)
+
+        ######################
+        # State of Charge
+        ######################
+        I = param.dimensional_current_with_time
+        # The `rhs` dictionary contains differential equations, with the key being the
+        # variable in the d/dt
+        self.rhs[Q] = I * param.timescale / 3600
+        # Initial conditions must be provided for the ODEs
+        self.initial_conditions[Q] = pybamm.Scalar(0)
+
+        ######################
+        # Particles
+        ######################
+
+        # The div and grad operators will be converted to the appropriate matrix
+        # multiplication at the discretisation stage
+        N_s_n = -param.D_n(c_s_n, T) * pybamm.grad(c_s_n)
+        N_s_p = -param.D_p(c_s_p, T) * pybamm.grad(c_s_p)
+        self.rhs[c_s_n] = -(1 / param.C_n) * pybamm.div(N_s_n)
+        self.rhs[c_s_p] = -(1 / param.C_p) * pybamm.div(N_s_p)
+        # Boundary conditions must be provided for equations with spatial derivatives
+        self.boundary_conditions[c_s_n] = {
+            "left": (pybamm.Scalar(0), "Neumann"),
+            "right": (-param.C_n * j_n / param.a_n, "Neumann"),
+        }
+        self.boundary_conditions[c_s_p] = {
+            "left": (pybamm.Scalar(0), "Neumann"),
+            "right": (-param.C_p * j_p / param.a_p / param.gamma_p, "Neumann"),
+        }
+        self.initial_conditions[c_s_n] = param.c_n_init
+        self.initial_conditions[c_s_p] = param.c_p_init
+        # Events specify points at which a solution should terminate
+        self.events.update(
+            {
+                "Minimum negative particle surface concentration": (
+                    pybamm.min(c_s_surf_n) - 0.01
+                ),
+                "Maximum negative particle surface concentration": (1 - 0.01)
+                - pybamm.max(c_s_surf_n),
+                "Minimum positive particle surface concentration": (
+                    pybamm.min(c_s_surf_p) - 0.01
+                ),
+                "Maximum positive particle surface concentration": (1 - 0.01)
+                - pybamm.max(c_s_surf_p),
+            }
+        )
+        ######################
+        # Current in the solid
+        ######################
+        i_s_n = -param.sigma_n * (1 - eps_n) ** param.b_s_n * pybamm.grad(phi_s_n)
+        sigma_eff_p = param.sigma_p * (1 - eps_p) ** param.b_s_p
+        i_s_p = -sigma_eff_p * pybamm.grad(phi_s_p)
+        # The `algebraic` dictionary contains differential equations, with the key being
+        # the main scalar variable of interest in the equation
+        self.algebraic[phi_s_n] = pybamm.div(i_s_n) + j_n
+        self.algebraic[phi_s_p] = pybamm.div(i_s_p) + j_p
+        self.boundary_conditions[phi_s_n] = {
+            "left": (pybamm.Scalar(0), "Dirichlet"),
+            "right": (pybamm.Scalar(0), "Neumann"),
+        }
+        self.boundary_conditions[phi_s_p] = {
+            "left": (pybamm.Scalar(0), "Neumann"),
+            "right": (i_cell / pybamm.boundary_value(-sigma_eff_p, "right"), "Neumann"),
+        }
+        # Initial conditions must also be provided for algebraic equations, as an
+        # initial guess for a root-finding algorithm which calculates consistent initial
+        # conditions
+        self.initial_conditions[phi_s_n] = pybamm.Scalar(0)
+        self.initial_conditions[phi_s_p] = param.U_p(
+            param.c_p_init, param.T_init
+        ) - param.U_n(param.c_n_init, param.T_init)
+
+        ######################
+        # Current in the electrolyte
+        ######################
+        i_e = (param.kappa_e(c_e, T) * tor * param.gamma_e / param.C_e) * (
+            param.chi(c_e) * pybamm.grad(c_e) / c_e - pybamm.grad(phi_e)
+        )
+        self.algebraic[phi_e] = pybamm.div(i_e) - j
+        self.boundary_conditions[phi_e] = {
+            "left": (pybamm.Scalar(0), "Neumann"),
+            "right": (pybamm.Scalar(0), "Neumann"),
+        }
+        self.initial_conditions[phi_e] = -param.U_n(param.c_n_init, param.T_init)
+
+        ######################
+        # Electrolyte concentration
+        ######################
+        N_e = -tor * param.D_e(c_e, T) * pybamm.grad(c_e)
+        self.rhs[c_e] = (1 / eps) * (
+            -pybamm.div(N_e) / param.C_e + (1 - param.t_plus) * j / param.gamma_e
+        )
+        self.boundary_conditions[c_e] = {
+            "left": (pybamm.Scalar(0), "Neumann"),
+            "right": (pybamm.Scalar(0), "Neumann"),
+        }
+        self.initial_conditions[c_e] = param.c_e_init
+        self.events["Zero electrolyte concentration cut-off"] = pybamm.min(c_e) - 0.002
+
+        ######################
+        # (Some) variables
+        ######################
+        voltage = pybamm.boundary_value(phi_s_p, "right")
+        # The `variables` dictionary contains all variables that might be useful for
+        # visualising the solution of the model
+        self.variables = {
+            "Negative particle surface concentration": c_s_surf_n,
+            "Electrolyte concentration": c_e,
+            "Positive particle surface concentration": c_s_surf_p,
+            "Current [A]": I,
+            "Negative electrode potential": phi_s_n,
+            "Electrolyte potential": phi_e,
+            "Positive electrode potential": phi_s_p,
+            "Terminal voltage": voltage,
+        }
+        self.events["Minimum voltage"] = voltage - param.voltage_low_cut
+        self.events["Maximum voltage"] = voltage - param.voltage_high_cut
+
+    @property
+    def default_geometry(self):
+        return pybamm.Geometry("1D macro", "1+1D micro")
diff --git a/pybamm/models/full_battery_models/lithium_ion/basic_spm.py b/pybamm/models/full_battery_models/lithium_ion/basic_spm.py
new file mode 100644
index 000000000..b907f5659
--- /dev/null
+++ b/pybamm/models/full_battery_models/lithium_ion/basic_spm.py
@@ -0,0 +1,172 @@
+#
+# Basic Single Particle Model (SPM)
+#
+import pybamm
+from .base_lithium_ion_model import BaseModel
+
+
+class BasicSPM(BaseModel):
+    """Single Particle Model (SPM) model of a lithium-ion battery, from [2]_.
+
+    This class differs from the :class:`pybamm.lithium_ion.SPM` model class in that it
+    shows the whole model in a single class. This comes at the cost of flexibility in
+    combining different physical effects, and in general the main SPM class should be
+    used instead.
+
+    Parameters
+    ----------
+    name : str, optional
+        The name of the model.
+
+    References
+    ----------
+    .. [2] SG Marquis, V Sulzer, R Timms, CP Please and SJ Chapman. “An asymptotic
+           derivation of a single particle model with electrolyte”. In: arXiv preprint
+           arXiv:1905.12553 (2019).
+
+
+    **Extends:** :class:`pybamm.lithium_ion.BaseModel`
+    """
+
+    def __init__(self, name="Single Particle Model"):
+        super().__init__({}, name)
+        # `param` is a class containing all the relevant parameters and functions for
+        # this model. These are purely symbolic at this stage, and will be set by the
+        # `ParameterValues` class when the model is processed.
+        param = self.param
+
+        ######################
+        # Variables
+        ######################
+        # Variables that depend on time only are created without a domain
+        Q = pybamm.Variable("Discharge capacity [A.h]")
+        # Variables that vary spatially are created with a domain
+        c_s_n = pybamm.Variable(
+            "X-averaged negative particle concentration", domain="negative particle",
+        )
+        c_s_p = pybamm.Variable(
+            "X-averaged positive particle concentration", domain="positive particle",
+        )
+
+        # Constant temperature
+        T = param.T_init
+
+        ######################
+        # Other set-up
+        ######################
+
+        # Current density
+        i_cell = param.current_with_time
+        j_n = i_cell / param.l_n
+        j_p = -i_cell / param.l_p
+
+        ######################
+        # State of Charge
+        ######################
+        I = param.dimensional_current_with_time
+        # The `rhs` dictionary contains differential equations, with the key being the
+        # variable in the d/dt
+        self.rhs[Q] = I * param.timescale / 3600
+        # Initial conditions must be provided for the ODEs
+        self.initial_conditions[Q] = pybamm.Scalar(0)
+
+        ######################
+        # Particles
+        ######################
+
+        # The div and grad operators will be converted to the appropriate matrix
+        # multiplication at the discretisation stage
+        N_s_n = -param.D_n(c_s_n, T) * pybamm.grad(c_s_n)
+        N_s_p = -param.D_p(c_s_p, T) * pybamm.grad(c_s_p)
+        self.rhs[c_s_n] = -(1 / param.C_n) * pybamm.div(N_s_n)
+        self.rhs[c_s_p] = -(1 / param.C_p) * pybamm.div(N_s_p)
+        # Boundary conditions must be provided for equations with spatial derivatives
+        self.boundary_conditions[c_s_n] = {
+            "left": (pybamm.Scalar(0), "Neumann"),
+            "right": (-param.C_n * j_n / param.a_n, "Neumann"),
+        }
+        self.boundary_conditions[c_s_p] = {
+            "left": (pybamm.Scalar(0), "Neumann"),
+            "right": (-param.C_p * j_p / param.a_p / param.gamma_p, "Neumann"),
+        }
+        self.initial_conditions[c_s_n] = param.c_n_init
+        self.initial_conditions[c_s_p] = param.c_p_init
+        # Surf takes the surface value of a variable, i.e. its boundary value on the
+        # right side. This is also accessible via `boundary_value(x, "right")`, with
+        # "left" providing the boundary value of the left side
+        c_s_surf_n = pybamm.surf(c_s_n)
+        c_s_surf_p = pybamm.surf(c_s_p)
+        # Events specify points at which a solution should terminate
+        self.events.update(
+            {
+                "Minimum negative particle surface concentration": (
+                    pybamm.min(c_s_surf_n) - 0.01
+                ),
+                "Maximum negative particle surface concentration": (1 - 0.01)
+                - pybamm.max(c_s_surf_n),
+                "Minimum positive particle surface concentration": (
+                    pybamm.min(c_s_surf_p) - 0.01
+                ),
+                "Maximum positive particle surface concentration": (1 - 0.01)
+                - pybamm.max(c_s_surf_p),
+            }
+        )
+
+        # Note that the SPM does not have any algebraic equations, so the `algebraic`
+        # dictionary remains empty
+
+        ######################
+        # (Some) variables
+        ######################
+        # Interfacial reactions
+        j0_n = (
+            param.m_n(T)
+            / param.C_r_n
+            * 1 ** (1 / 2)
+            * c_s_surf_n ** (1 / 2)
+            * (1 - c_s_surf_n) ** (1 / 2)
+        )
+        j0_p = (
+            param.gamma_p
+            * param.m_p(T)
+            / param.C_r_p
+            * 1 ** (1 / 2)
+            * c_s_surf_p ** (1 / 2)
+            * (1 - c_s_surf_p) ** (1 / 2)
+        )
+        eta_n = (2 / param.ne_n) * pybamm.arcsinh(j_n / (2 * j0_n))
+        eta_p = (2 / param.ne_p) * pybamm.arcsinh(j_p / (2 * j0_p))
+        phi_s_n = 0
+        phi_e = -eta_n - param.U_n(c_s_surf_n, T)
+        phi_s_p = eta_p + phi_e + param.U_p(c_s_surf_p, T)
+        V = phi_s_p
+
+        whole_cell = ["negative electrode", "separator", "positive electrode"]
+        # The `variables` dictionary contains all variables that might be useful for
+        # visualising the solution of the model
+        # Primary broadcasts are used to broadcast scalar quantities across a domain
+        # into a vector of the right shape, for multiplying with other vectors
+        self.variables = {
+            "Negative particle surface concentration": pybamm.PrimaryBroadcast(
+                c_s_surf_n, "negative electrode"
+            ),
+            "Electrolyte concentration": pybamm.PrimaryBroadcast(1, whole_cell),
+            "Positive particle surface concentration": pybamm.PrimaryBroadcast(
+                c_s_surf_p, "positive electrode"
+            ),
+            "Current [A]": I,
+            "Negative electrode potential": pybamm.PrimaryBroadcast(
+                phi_s_n, "negative electrode"
+            ),
+            "Electrolyte potential": pybamm.PrimaryBroadcast(phi_e, whole_cell),
+            "Positive electrode potential": pybamm.PrimaryBroadcast(
+                phi_s_p, "positive electrode"
+            ),
+            "Terminal voltage": V,
+        }
+        self.events["Minimum voltage"] = V - param.voltage_low_cut
+        self.events["Maximum voltage"] = V - param.voltage_high_cut
+
+    @property
+    def default_geometry(self):
+        return pybamm.Geometry("1D macro", "1D micro")
diff --git a/pybamm/models/submodels/particle/base_particle.py b/pybamm/models/submodels/particle/base_particle.py
index 5360472aa..bbc3cc487 100644
--- a/pybamm/models/submodels/particle/base_particle.py
+++ b/pybamm/models/submodels/particle/base_particle.py
@@ -92,7 +92,7 @@ class BaseParticle(pybamm.BaseSubModel):
         tol = 0.01
 
         self.events[
-            "Minumum " + self.domain.lower() + " particle surface concentration"
+            "Minimum " + self.domain.lower() + " particle surface concentration"
         ] = (pybamm.min(c_s_surf) - tol)
 
         self.events[
diff --git a/pybamm/models/submodels/thermal/isothermal/isothermal.py b/pybamm/models/submodels/thermal/isothermal/isothermal.py
index 3c12d58d5..f4d6a7f44 100644
--- a/pybamm/models/submodels/thermal/isothermal/isothermal.py
+++ b/pybamm/models/submodels/thermal/isothermal/isothermal.py
@@ -23,7 +23,7 @@ class Isothermal(BaseThermal):
 
     def get_fundamental_variables(self):
 
-        T_x_av = pybamm.PrimaryBroadcast(0, "current collector")
+        T_x_av = pybamm.PrimaryBroadcast(self.param.T_init, "current collector")
         T_n = pybamm.PrimaryBroadcast(T_x_av, "negative electrode")
         T_s = pybamm.PrimaryBroadcast(T_x_av, "separator")
         T_p = pybamm.PrimaryBroadcast(T_x_av, "positive electrode")
diff --git a/pybamm/processed_variable.py b/pybamm/processed_variable.py
index 847fdef80..77aef5bfe 100644
--- a/pybamm/processed_variable.py
+++ b/pybamm/processed_variable.py
@@ -41,14 +41,12 @@ class ProcessedVariable(object):
             self.base_eval, self.known_evals[solution.t[0]] = base_variable.evaluate(
                 solution.t[0],
                 solution.y[:, 0],
-                {name: inp[0] for name, inp in solution.inputs.items()},
+                solution.inputs,
                 known_evals=self.known_evals[solution.t[0]],
             )
         else:
             self.base_eval = base_variable.evaluate(
-                solution.t[0],
-                solution.y[:, 0],
-                {name: inp[0] for name, inp in solution.inputs.items()},
+                solution.t[0], solution.y[:, 0], solution.inputs
             )
 
         # handle 2D (in space) finite element variables differently
@@ -92,14 +90,14 @@ class ProcessedVariable(object):
         # Evaluate the base_variable index-by-index
         for idx in range(len(self.t_sol)):
             t = self.t_sol[idx]
-            u = self.u_sol[:, idx]
-            inputs = {name: inp[0] for name, inp in self.inputs.items()}
             if self.known_evals:
                 entries[idx], self.known_evals[t] = self.base_variable.evaluate(
-                    t, u, inputs, known_evals=self.known_evals[t]
+                    t, self.u_sol[:, idx], self.inputs, known_evals=self.known_evals[t]
                 )
             else:
-                entries[idx] = self.base_variable.evaluate(t, u, inputs)
+                entries[idx] = self.base_variable.evaluate(
+                    t, self.u_sol[:, idx], self.inputs
+                )
 
         # No discretisation provided, or variable has no domain (function of t only)
         self._interpolation_function = interp.interp1d(
@@ -117,15 +115,14 @@ class ProcessedVariable(object):
         for idx in range(len(self.t_sol)):
             t = self.t_sol[idx]
             u = self.u_sol[:, idx]
-            inputs = {name: inp[0] for name, inp in self.inputs.items()}
             if self.known_evals:
                 eval_and_known_evals = self.base_variable.evaluate(
-                    t, u, inputs, known_evals=self.known_evals[t]
+                    t, u, self.inputs, known_evals=self.known_evals[t]
                 )
                 entries[:, idx] = eval_and_known_evals[0][:, 0]
                 self.known_evals[t] = eval_and_known_evals[1]
             else:
-                entries[:, idx] = self.base_variable.evaluate(t, u, inputs)[:, 0]
+                entries[:, idx] = self.base_variable.evaluate(t, u, self.inputs)[:, 0]
 
         # Process the discretisation to get x values
         nodes = self.mesh[0].nodes
@@ -221,10 +218,9 @@ class ProcessedVariable(object):
         for idx in range(len(self.t_sol)):
             t = self.t_sol[idx]
             u = self.u_sol[:, idx]
-            inputs = {name: inp[0] for name, inp in self.inputs.items()}
             if self.known_evals:
                 eval_and_known_evals = self.base_variable.evaluate(
-                    t, u, inputs, known_evals=self.known_evals[t]
+                    t, u, self.inputs, known_evals=self.known_evals[t]
                 )
                 entries[:, :, idx] = np.reshape(
                     eval_and_known_evals[0],
@@ -234,7 +230,7 @@ class ProcessedVariable(object):
                 self.known_evals[t] = eval_and_known_evals[1]
             else:
                 entries[:, :, idx] = np.reshape(
-                    self.base_variable.evaluate(t, u, inputs),
+                    self.base_variable.evaluate(t, u, self.inputs),
                     [first_dim_size, second_dim_size],
                     order="F",
                 )


Test patch:
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_basic_models.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_basic_models.py
new file mode 100644
index 000000000..dfb2f3c89
--- /dev/null
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_basic_models.py
@@ -0,0 +1,68 @@
+#
+# Compare basic models with full models
+#
+import pybamm
+
+import numpy as np
+import unittest
+
+
+class TestCompareBasicModels(unittest.TestCase):
+    def test_compare_dfns(self):
+        basic_dfn = pybamm.lithium_ion.BasicDFN()
+        dfn = pybamm.lithium_ion.DFN()
+
+        # Solve basic DFN
+        basic_sim = pybamm.Simulation(basic_dfn)
+        t_eval = np.linspace(0, 1)
+        basic_sim.solve(t_eval)
+        basic_sol = basic_sim.solution
+
+        # Solve main DFN
+        sim = pybamm.Simulation(dfn)
+        t_eval = np.linspace(0, 1)
+        sim.solve(t_eval)
+        sol = sim.solution
+
+        # Compare solution data
+        np.testing.assert_array_almost_equal(basic_sol.y, sol.y, decimal=4)
+        np.testing.assert_array_almost_equal(basic_sol.t, sol.t, decimal=4)
+        # Compare variables
+        for name in basic_dfn.variables:
+            np.testing.assert_array_almost_equal(
+                basic_sol[name].entries, sol[name].entries, decimal=4
+            )
+
+    def test_compare_spms(self):
+        basic_spm = pybamm.lithium_ion.BasicSPM()
+        spm = pybamm.lithium_ion.SPM()
+
+        # Solve basic SPM
+        basic_sim = pybamm.Simulation(basic_spm)
+        t_eval = np.linspace(0, 1)
+        basic_sim.solve(t_eval)
+        basic_sol = basic_sim.solution
+
+        # Solve main SPM
+        sim = pybamm.Simulation(spm)
+        t_eval = np.linspace(0, 1)
+        sim.solve(t_eval)
+        sol = sim.solution
+
+        # Compare solution data
+        np.testing.assert_array_almost_equal(basic_sol.y, sol.y)
+        np.testing.assert_array_almost_equal(basic_sol.t, sol.t)
+        # Compare variables
+        for name in basic_spm.variables:
+            np.testing.assert_array_almost_equal(
+                basic_sol[name].entries, sol[name].entries
+            )
+
+
+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_full_battery_models/test_lithium_ion/test_basic_models.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_basic_models.py
new file mode 100644
index 000000000..2b56b0bb2
--- /dev/null
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_basic_models.py
@@ -0,0 +1,35 @@
+#
+# Tests for the basic lithium-ion models
+#
+import pybamm
+import unittest
+
+
+class TestBasicModels(unittest.TestCase):
+    def test_dfn_well_posed(self):
+        model = pybamm.lithium_ion.BasicDFN()
+        model.check_well_posedness()
+
+    def test_dfn_default_geometry(self):
+        model = pybamm.lithium_ion.BasicDFN()
+        self.assertIsInstance(model.default_geometry, pybamm.Geometry)
+        self.assertTrue("secondary" in model.default_geometry["negative particle"])
+
+    def test_spm_well_posed(self):
+        model = pybamm.lithium_ion.BasicSPM()
+        model.check_well_posedness()
+
+    def test_spm_default_geometry(self):
+        model = pybamm.lithium_ion.BasicSPM()
+        self.assertIsInstance(model.default_geometry, pybamm.Geometry)
+        self.assertTrue("secondary" not in model.default_geometry["negative particle"])
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    pybamm.settings.debug_mode = True
+    unittest.main()
diff --git a/tests/unit/test_solvers/test_solution.py b/tests/unit/test_solvers/test_solution.py
index 7aa4db34f..32e81cfbc 100644
--- a/tests/unit/test_solvers/test_solution.py
+++ b/tests/unit/test_solvers/test_solution.py
@@ -95,30 +95,6 @@ class TestSolution(unittest.TestCase):
         self.assertEqual(solution.model.name, solution_load.model.name)
         np.testing.assert_array_equal(solution["c"].entries, solution_load["c"].entries)
 
-    def test_solution_evals_with_inputs(self):
-        model = pybamm.lithium_ion.SPM()
-        geometry = model.default_geometry
-        param = model.default_parameter_values
-        param.update({"Electrode height [m]": "[input]"})
-        param.process_model(model)
-        param.process_geometry(geometry)
-        var = pybamm.standard_spatial_vars
-        var_pts = {var.x_n: 5, var.x_s: 5, var.x_p: 5, var.r_n: 10, var.r_p: 10}
-        spatial_methods = model.default_spatial_methods
-        solver = model.default_solver
-        sim = pybamm.Simulation(
-            model=model,
-            geometry=geometry,
-            parameter_values=param,
-            var_pts=var_pts,
-            spatial_methods=spatial_methods,
-            solver=solver,
-        )
-        inputs = {"Electrode height [m]": 0.1}
-        sim.solve(t_eval=np.linspace(0, 0.01, 10), inputs=inputs)
-        time = sim.solution["Time [h]"](sim.solution.t)
-        self.assertEqual(len(time), 10)
-
 
 if __name__ == "__main__":
     print("Add -v for more debug output")