Repository: tinosulzer/PyBaMM

Issue:
Bug when using model.rhs.evaluate() 
There is an issue using the `model.rhs.evaluate()` function which depends on the form of the flux `N_e` in the electrolyte model which is indicated by an arrow in the code below: 

```
       epsilon = pybamm.standard_parameters.epsilon_s  
        b = pybamm.standard_parameters.b
        delta = pybamm.standard_parameters.delta
        nu = pybamm.standard_parameters.nu
        t_plus = pybamm.standard_parameters.t_plus
        ce0 = pybamm.standard_parameters.ce0

        electrolyte_domain = ["whole cell"]
        c_e = pybamm.Variable("c_e", electrolyte_domain)

        N_e = -pybamm.grad(c_e)  <--------------------------------------------------(1)
        self.rhs = {
            c_e: -pybamm.div(N_e) / delta / epsilon + nu * (1 - t_plus) * G
        }  # / delta / epsilon + nu * (1 - t_plus) * G
        self.initial_conditions = {c_e: ce0}
        self.boundary_conditions = {
            N_e: {"left": pybamm.Scalar(0), "right": pybamm.Scalar(0)}
        } 
```
Within `test_electrolyte.py` I have
```
        G = pybamm.Scalar(1)
        model = pybamm.electrolyte.StefanMaxwellDiffusion(G)

        param = pybamm.ParameterValues(
            "input/parameters/lithium-ion/parameters/LCO.csv"
        )

        param.process_model(model)

        mesh = pybamm.FiniteVolumeMacroMesh(param, 2)
        disc = pybamm.FiniteVolumeDiscretisation(mesh)
        disc.process_model(model)

        y0 = model.initial_conditions
        model.rhs.evaluate(None, y0) <---------------------------------------(2)
```
However, if i change (1) to `N_e = epsilon * pybamm.grad(c_e)` , there is an error at (2). Link to the failed test output: https://gist.github.com/Scottmar93/b7bcbff685383d9a5dceaca630f3aa40

Discussion:

Patch:
diff --git a/.flake8 b/.flake8
index 70990935..9acc35b7 100644
--- a/.flake8
+++ b/.flake8
@@ -5,6 +5,12 @@ exclude=
     problems,
     __init__.py,
     venv,
+    bin,
+    etc,
+    lib,
+    lib64,
+    share,
+    pyvenv.cfg
 ignore=
     # False positive for white space before ':' on list slice
     # black should format these correctly
diff --git a/.gitignore b/.gitignore
index 65bb1425..4db0408e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -46,3 +46,9 @@ pyproject.toml
 
 # virtual enviroment
 venv/
+bin/
+etc/
+lib/
+lib64
+share/
+pyvenv.cfg
diff --git a/pybamm/models/submodels/electrolyte.py b/pybamm/models/submodels/electrolyte.py
index acad18a1..39609690 100644
--- a/pybamm/models/submodels/electrolyte.py
+++ b/pybamm/models/submodels/electrolyte.py
@@ -42,7 +42,7 @@ class StefanMaxwellDiffusion(pybamm.BaseModel):
 
         electrolyte_domain = ["whole cell"]
 
-        c_e = pybamm.Variable("c_e", domain=electrolyte_domain)
+        c_e = pybamm.Variable("c_e", electrolyte_domain)
 
         N_e = -(epsilon ** b) * pybamm.grad(c_e)
 
diff --git a/pybamm/parameters/parameter_values.py b/pybamm/parameters/parameter_values.py
index fe86c94d..74834bdd 100644
--- a/pybamm/parameters/parameter_values.py
+++ b/pybamm/parameters/parameter_values.py
@@ -95,11 +95,15 @@ class ParameterValues(dict):
             model.initial_conditions[variable] = self.process_symbol(equation)
 
         # Boundary conditions are dictionaries {"left": left bc, "right": right bc}
+        new_boundary_conditions = {}
         for variable, bcs in model.boundary_conditions.items():
+            processed_variable = self.process_symbol(variable)
+            new_boundary_conditions[processed_variable] = {}
             for side in ["left", "right"]:
-                model.boundary_conditions[variable][side] = self.process_symbol(
+                new_boundary_conditions[processed_variable][side] = self.process_symbol(
                     bcs[side]
                 )
+        model.boundary_conditions = new_boundary_conditions
 
         for variable, equation in model.variables.items():
             model.variables[variable] = self.process_symbol(equation)


Test patch:
diff --git a/tests/test_models/test_electrolyte.py b/tests/test_models/test_electrolyte.py
index 8b8a2cd0..686f1096 100644
--- a/tests/test_models/test_electrolyte.py
+++ b/tests/test_models/test_electrolyte.py
@@ -6,6 +6,7 @@ from __future__ import print_function, unicode_literals
 import pybamm
 
 import unittest
+import numpy as np
 
 
 class TestStefanMaxwellDiffusion(unittest.TestCase):
@@ -33,8 +34,28 @@ class TestStefanMaxwellDiffusion(unittest.TestCase):
 
         param.process_model(model)
 
-        print(model.initial_conditions)
+        mesh = pybamm.FiniteVolumeMacroMesh(param, 2)
+        disc = pybamm.FiniteVolumeDiscretisation(mesh)
+        disc.process_model(model)
+
+    def test_solving(self):
+        G = pybamm.Scalar(0.001)
+        model = pybamm.electrolyte.StefanMaxwellDiffusion(G)
+
+        param = pybamm.ParameterValues(
+            "input/parameters/lithium-ion/parameters/LCO.csv"
+        )
+
+        param.process_model(model)
 
         mesh = pybamm.FiniteVolumeMacroMesh(param, 2)
         disc = pybamm.FiniteVolumeDiscretisation(mesh)
         disc.process_model(model)
+
+        y0 = model.concatenated_initial_conditions
+        np.testing.assert_array_equal(y0, np.ones_like(mesh["whole cell"].nodes))
+
+        # Solve
+        solver = pybamm.ScipySolver(tol=1e-8, method="RK45")
+        t_eval = mesh["time"]
+        solver.solve(model, t_eval)
diff --git a/tests/test_parameters/test_parameter_values.py b/tests/test_parameters/test_parameter_values.py
index 90321430..cc008556 100644
--- a/tests/test_parameters/test_parameter_values.py
+++ b/tests/test_parameters/test_parameter_values.py
@@ -156,18 +156,17 @@ class TestParameterValues(unittest.TestCase):
         self.assertTrue(isinstance(model.rhs[var].children[0], pybamm.Scalar))
         self.assertTrue(isinstance(model.rhs[var].children[1], pybamm.Gradient))
         self.assertEqual(model.rhs[var].children[0].value, 1)
-        # initial_conditions
+        # initial conditions
         self.assertTrue(isinstance(model.initial_conditions[var], pybamm.Scalar))
         self.assertEqual(model.initial_conditions[var].value, 2)
-        # boundary_conditions
-        self.assertTrue(
-            isinstance(model.boundary_conditions[var]["left"], pybamm.Scalar)
-        )
-        self.assertEqual(model.boundary_conditions[var]["left"].value, 3)
-        self.assertTrue(
-            isinstance(model.boundary_conditions[var]["right"], pybamm.Scalar)
-        )
-        self.assertEqual(model.boundary_conditions[var]["right"].value, 42)
+        # boundary conditions
+        bc_key = list(model.boundary_conditions.keys())[0]
+        self.assertTrue(isinstance(bc_key, pybamm.Variable))
+        bc_value = list(model.boundary_conditions.values())[0]
+        self.assertTrue(isinstance(bc_value["left"], pybamm.Scalar))
+        self.assertEqual(bc_value["left"].value, 3)
+        self.assertTrue(isinstance(bc_value["right"], pybamm.Scalar))
+        self.assertEqual(bc_value["right"].value, 42)
         # variables
         self.assertEqual(model.variables["var"].id, var.id)
         self.assertTrue(isinstance(model.variables["grad_var"], pybamm.Gradient))