Merge pull request #1476 from pybamm-team/issue-1456-concentration

Issue 1456 concentration

Author: valentinsulzer

CHANGELOG.md

@@ -17,6 +17,8 @@
 ## Bug fixes
 
 ## Breaking changes
+
+-   Changed the variable in the full diffusion model from "Electrolyte concentration" to "Porosity times concentration" ([#1476](https://github.com/pybamm-team/PyBaMM/pull/1476))
 -   Renamed `lithium-ion` folder to `lithium_ion` and `lead-acid` folder to `lead_acid` in parameters ([#1464](https://github.com/pybamm-team/PyBaMM/pull/1464))
 
 # [v0.4.0](https://github.com/pybamm-team/PyBaMM/tree/v0.4.0) - 2021-03-28

examples/notebooks/initialize-model-with-solution.ipynb

@@ -23,6 +23,8 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
+      "\u001b[33mWARNING: You are using pip version 21.0.1; however, version 21.1 is available.\n",
+      "You should consider upgrading via the '/Users/vsulzer/Documents/Energy_storage/PyBaMM/.tox/dev/bin/python -m pip install --upgrade pip' command.\u001b[0m\n",
       "Note: you may need to restart the kernel to use updated packages.\n"
      ]
     }
@@ -159,7 +161,7 @@
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "0b5f9d809e244114b29f527062decbe9",
+       "model_id": "7ad44c260ec446e0bd04a2f7dfbf8c65",
        "version_major": 2,
        "version_minor": 0
       },
@@ -173,7 +175,7 @@
     {
      "data": {
       "text/plain": [
-       "<pybamm.plotting.quick_plot.QuickPlot at 0x7f3292dc9710>"
+       "<pybamm.plotting.quick_plot.QuickPlot at 0x1447078e0>"
       ]
      },
      "execution_count": 6,
@@ -240,7 +242,7 @@
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "d1e3967157584ff3802279544b0cc442",
+       "model_id": "aaf63750ed354d1987703fd2fdde83af",
        "version_major": 2,
        "version_minor": 0
       },
@@ -254,7 +256,7 @@
     {
      "data": {
       "text/plain": [
-       "<pybamm.plotting.quick_plot.QuickPlot at 0x7f32818d1990>"
+       "<pybamm.plotting.quick_plot.QuickPlot at 0x10437c370>"
       ]
      },
      "execution_count": 8,
@@ -299,6 +301,13 @@
    "source": [
     "pybamm.print_citations()"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
@@ -317,7 +326,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.4"
+   "version": "3.8.8"
   },
   "toc": {
    "base_numbering": 1,

pybamm/expression_tree/binary_operators.py

@@ -1192,7 +1192,10 @@ def simplified_matrix_multiplication(left, right):
         if right.right.evaluates_to_constant_number():
             r_left, r_right = right.orphans
             new_left = left / r_right
-            return new_left @ r_left
+            new_mul = new_left @ r_left
+            # Keep the domain of the old left
+            new_mul.copy_domains(left)
+            return new_mul
 
     # Simplify A @ (B @ c) to (A @ B) @ c if (A @ B) is constant
     # This is a common construction that appears from discretisation of spatial

pybamm/models/base_model.py

@@ -653,25 +653,6 @@ def check_ics_bcs(self):
                     """no initial condition given for variable '{}'""".format(var)
                 )
 
-        # Boundary conditions
-        for var, eqn in {**self.rhs, **self.algebraic}.items():
-            if eqn.has_symbol_of_classes(
-                (pybamm.Gradient, pybamm.Divergence)
-            ) and not eqn.has_symbol_of_classes(pybamm.Integral):
-                # I have relaxed this check for now so that the lumped temperature
-                # equation doesn't raise errors (this has and average in it)
-
-                # Variable must be in the boundary conditions
-                if not any(
-                    var.id == x.id
-                    for symbol in self.boundary_conditions.keys()
-                    for x in symbol.pre_order()
-                ):
-                    raise pybamm.ModelError(
-                        "no boundary condition given for "
-                        "variable '{}' with equation '{}'.".format(var, eqn)
-                    )
-
     def check_default_variables_dictionaries(self):
         """Check that the right variables are provided."""
         missing_vars = []

pybamm/models/standard_variables.py

@@ -31,6 +31,33 @@
     bounds=(0, np.inf),
 )
 
+# Electrolyte porosity times concentration
+eps_c_e_n = pybamm.Variable(
+    "Negative electrode porosity times concentration",
+    domain="negative electrode",
+    auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, np.inf),
+)
+eps_c_e_s = pybamm.Variable(
+    "Separator porosity times concentration",
+    domain="separator",
+    auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, np.inf),
+)
+eps_c_e_p = pybamm.Variable(
+    "Positive electrode porosity times concentration",
+    domain="positive electrode",
+    auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, np.inf),
+)
+eps_c_e = pybamm.concatenation(eps_c_e_n, eps_c_e_s, eps_c_e_p)
+
+eps_c_e_av = pybamm.Variable(
+    "X-averaged porosity times concentration",
+    domain="current collector",
+    bounds=(0, np.inf),
+)
+
 # Electrolyte potential
 phi_e_n = pybamm.Variable(
     "Negative electrolyte potential",

pybamm/models/submodels/electrolyte_diffusion/base_electrolyte_diffusion.py

@@ -78,6 +78,19 @@ def _get_standard_concentration_variables(self, c_e_n, c_e_s, c_e_p):
 
         return variables
 
+    def _get_standard_porosity_times_concentration_variables(
+        self, eps_c_e_n, eps_c_e_s, eps_c_e_p
+    ):
+        eps_c_e = pybamm.concatenation(eps_c_e_n, eps_c_e_s, eps_c_e_p)
+
+        variables = {
+            "Porosity times concentration": eps_c_e,
+            "Negative electrode porosity times concentration": eps_c_e_n,
+            "Separator porosity times concentration": eps_c_e_s,
+            "Positive electrode porosity times concentration": eps_c_e_p,
+        }
+        return variables
+
     def _get_standard_flux_variables(self, N_e):
         """
         A private function to obtain the standard variables which

pybamm/models/submodels/electrolyte_diffusion/constant_concentration.py

@@ -39,6 +39,19 @@ def get_fundamental_variables(self):
         return variables
 
     def get_coupled_variables(self, variables):
+        eps_n = variables["Negative electrode porosity"]
+        eps_s = variables["Separator porosity"]
+        eps_p = variables["Positive electrode porosity"]
+        c_e_n = variables["Negative electrolyte concentration"]
+        c_e_s = variables["Separator electrolyte concentration"]
+        c_e_p = variables["Positive electrolyte concentration"]
+
+        variables.update(
+            self._get_standard_porosity_times_concentration_variables(
+                eps_n * c_e_n, eps_s * c_e_s, eps_p * c_e_p
+            )
+        )
+
         c_e = variables["Electrolyte concentration"]
         eps = variables["Porosity"]
 

pybamm/models/submodels/electrolyte_diffusion/full_diffusion.py

@@ -25,17 +25,36 @@ def __init__(self, param):
         super().__init__(param)
 
     def get_fundamental_variables(self):
-        c_e_n = pybamm.standard_variables.c_e_n
-        c_e_s = pybamm.standard_variables.c_e_s
-        c_e_p = pybamm.standard_variables.c_e_p
+        eps_c_e_n = pybamm.standard_variables.eps_c_e_n
+        eps_c_e_s = pybamm.standard_variables.eps_c_e_s
+        eps_c_e_p = pybamm.standard_variables.eps_c_e_p
 
-        return self._get_standard_concentration_variables(c_e_n, c_e_s, c_e_p)
+        variables = self._get_standard_porosity_times_concentration_variables(
+            eps_c_e_n, eps_c_e_s, eps_c_e_p
+        )
+
+        return variables
 
     def get_coupled_variables(self, variables):
 
-        tor = variables["Electrolyte tortuosity"]
+        eps_n = variables["Negative electrode porosity"]
+        eps_s = variables["Separator porosity"]
+        eps_p = variables["Positive electrode porosity"]
+        eps_c_e_n = variables["Negative electrode porosity times concentration"]
+        eps_c_e_s = variables["Separator porosity times concentration"]
+        eps_c_e_p = variables["Positive electrode porosity times concentration"]
+
+        c_e_n = eps_c_e_n / eps_n
+        c_e_s = eps_c_e_s / eps_s
+        c_e_p = eps_c_e_p / eps_p
+
+        variables.update(
+            self._get_standard_concentration_variables(c_e_n, c_e_s, c_e_p)
+        )
+
         eps = variables["Porosity"]
         c_e = variables["Electrolyte concentration"]
+        tor = variables["Electrolyte tortuosity"]
         i_e = variables["Electrolyte current density"]
         v_box = variables["Volume-averaged velocity"]
         T = variables["Cell temperature"]
@@ -57,8 +76,7 @@ def set_rhs(self, variables):
 
         param = self.param
 
-        eps = variables["Porosity"]
-        deps_dt = variables["Porosity change"]
+        eps_c_e = variables["Porosity times concentration"]
         c_e = variables["Electrolyte concentration"]
         N_e = variables["Electrolyte flux"]
         div_Vbox = variables["Transverse volume-averaged acceleration"]
@@ -67,20 +85,16 @@ def set_rhs(self, variables):
         source_terms = sum_s_j / self.param.gamma_e
 
         self.rhs = {
-            c_e: (1 / eps)
-            * (
-                -pybamm.div(N_e) / param.C_e
-                + source_terms
-                - c_e * deps_dt
-                - c_e * div_Vbox
-            )
+            eps_c_e: -pybamm.div(N_e) / param.C_e + source_terms - c_e * div_Vbox
         }
 
     def set_initial_conditions(self, variables):
 
-        c_e = variables["Electrolyte concentration"]
+        eps_c_e = variables["Porosity times concentration"]
 
-        self.initial_conditions = {c_e: self.param.c_e_init}
+        self.initial_conditions = {
+            eps_c_e: self.param.epsilon_init * self.param.c_e_init
+        }
 
     def set_boundary_conditions(self, variables):
 
@@ -90,5 +104,5 @@ def set_boundary_conditions(self, variables):
             c_e: {
                 "left": (pybamm.Scalar(0), "Neumann"),
                 "right": (pybamm.Scalar(0), "Neumann"),
-            }
+            },
         }

tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_basic_models.py

@@ -25,7 +25,6 @@ def test_compare_dfns(self):
         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:

tests/integration/test_solvers/test_idaklu.py

@@ -31,7 +31,7 @@ def test_set_tol_by_variable(self):
         t_eval = np.linspace(0, 3600, 100)
         solver = pybamm.IDAKLUSolver()
 
-        variable_tols = {"Electrolyte concentration": 1e-3}
+        variable_tols = {"Porosity times concentration": 1e-3}
         solver.set_atol_by_variable(variable_tols, model)
 
         solver.solve(model, t_eval)

tests/unit/test_expression_tree/test_symbol.py

@@ -373,6 +373,9 @@ def test_symbol_visualise(self):
 
         variables = {
             "Porosity": param.epsilon,
+            "Negative electrode porosity": param.epsilon_n,
+            "Separator porosity": param.epsilon_s,
+            "Positive electrode porosity": param.epsilon_p,
             "Electrolyte tortuosity": param.epsilon ** 1.5,
             "Porosity change": zero_nsp,
             "Electrolyte current density": zero_nsp,

tests/unit/test_models/test_base_model.py

@@ -420,24 +420,6 @@ def test_check_well_posedness_initial_boundary_conditions(self):
         with self.assertRaisesRegex(pybamm.ModelError, "initial condition"):
             model.check_well_posedness()
 
-        # Model with bad boundary conditions - Dirichlet (expect assertion error)
-        d = pybamm.Variable("d", domain=whole_cell)
-        model.initial_conditions = {c: 3}
-        model.boundary_conditions = {
-            d: {"left": (0, "Dirichlet"), "right": (0, "Dirichlet")}
-        }
-        with self.assertRaisesRegex(pybamm.ModelError, "boundary condition"):
-            model.check_well_posedness()
-
-        # Model with bad boundary conditions - Neumann (expect assertion error)
-        d = pybamm.Variable("d", domain=whole_cell)
-        model.initial_conditions = {c: 3}
-        model.boundary_conditions = {
-            d: {"left": (0, "Neumann"), "right": (0, "Neumann")}
-        }
-        with self.assertRaisesRegex(pybamm.ModelError, "boundary condition"):
-            model.check_well_posedness()
-
         # Algebraic well-posed model
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         model = pybamm.BaseModel()
@@ -451,13 +433,6 @@ def test_check_well_posedness_initial_boundary_conditions(self):
         }
         model.check_well_posedness()
 
-        # Algebraic model with bad boundary conditions
-        model.boundary_conditions = {
-            d: {"left": (0, "Dirichlet"), "right": (0, "Dirichlet")}
-        }
-        with self.assertRaisesRegex(pybamm.ModelError, "boundary condition"):
-            model.check_well_posedness()
-
     def test_check_well_posedness_output_variables(self):
         model = pybamm.BaseModel()
         whole_cell = ["negative electrode", "separator", "positive electrode"]

tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spme.py

@@ -17,7 +17,7 @@ def test_external_variables(self):
         model = pybamm.lithium_ion.SPMe(model_options)
         self.assertEqual(
             model.external_variables[0].id,
-            model.variables["Electrolyte concentration"].id,
+            model.variables["Porosity times concentration"].id,
         )
 
         # a variable

tests/unit/test_models/test_submodels/test_electrolyte_diffusion/test_constant_concentration.py

@@ -11,7 +11,13 @@ class TestConstantConcentration(unittest.TestCase):
     def test_public_functions(self):
         param = pybamm.LithiumIonParameters()
         a = pybamm.Scalar(0)
-        variables = {"Porosity": a, "Electrolyte concentration": a}
+        variables = {
+            "Porosity": a,
+            "Negative electrode porosity": a,
+            "Separator porosity": a,
+            "Positive electrode porosity": a,
+            "Electrolyte concentration": a,
+        }
         submodel = pybamm.electrolyte_diffusion.ConstantConcentration(param)
         std_tests = tests.StandardSubModelTests(submodel, variables)
         std_tests.test_all()

tests/unit/test_models/test_submodels/test_electrolyte_diffusion/test_full_diffusion.py

@@ -18,6 +18,9 @@ def test_public_functions(self):
         )
         variables = {
             "Porosity": a,
+            "Negative electrode porosity": a,
+            "Separator porosity": a,
+            "Positive electrode porosity": a,
             "Electrolyte tortuosity": a,
             "Porosity change": a,
             "Volume-averaged velocity": a,

tests/unit/test_solvers/test_idaklu_solver.py

@@ -57,7 +57,7 @@ def test_set_atol(self):
         disc.process_model(model)
         solver = pybamm.IDAKLUSolver(root_method="lm")
 
-        variable_tols = {"Electrolyte concentration": 1e-3}
+        variable_tols = {"Porosity times concentration": 1e-3}
         solver.set_atol_by_variable(variable_tols, model)
 
     def test_failures(self):