Changed fundamental SEI variable from thickness to concentration

CHANGELOG.md

@@ -87,6 +87,8 @@
 - Fix to simulate c_rate steps with drive cycles ([#3186](https://github.com/pybamm-team/PyBaMM/pull/3186))
 - Always save last cycle in experiment, to fix issues with `starting_solution` and `last_state` ([#3177](https://github.com/pybamm-team/PyBaMM/pull/3177))
 - Fix simulations with `starting_solution` to work with `start_time` experiments ([#3177](https://github.com/pybamm-team/PyBaMM/pull/3177))
+- Fixed how `Loss of lithium due to loss of active material` variable is calculated ([#3171](https://github.com/pybamm-team/PyBaMM/pull/3171))
+- Rewritten SEI submodel so that concentration, not thickness, is the fundamental variable. This change ensures lithium is always conserved. ([#3171](https://github.com/pybamm-team/PyBaMM/pull/3171))
 - Fix SEI Example Notebook ([#3166](https://github.com/pybamm-team/PyBaMM/pull/3166))
 - Thevenin() model is now constructed with standard variables: `Time [s]`, `Time [min]`, `Time [h]` ([#3143](https://github.com/pybamm-team/PyBaMM/pull/3143))
 - Error generated when invalid parameter values are passed ([#3132](https://github.com/pybamm-team/PyBaMM/pull/3132))

docs/source/examples/notebooks/models/coupled-degradation.ipynb

@@ -57,7 +57,6 @@
     "        \"particle mechanics\": (\"swelling and cracking\", \"swelling only\"),\n",
     "        \"SEI on cracks\": \"true\",\n",
     "        \"loss of active material\": \"stress-driven\",\n",
-    "        \"calculate discharge energy\": \"true\",  # for compatibility with older PyBaMM versions\n",
     "    }\n",
     ")"
    ]

docs/source/examples/notebooks/models/using-submodels.ipynb

@@ -515,15 +515,27 @@
     "model.submodels[\"Negative sei\"] = pybamm.sei.NoSEI(\n",
     "    model.param, \"negative\", model.options\n",
     ")\n",
+    "model.submodels[\"Negative sei thickness\"] = pybamm.sei.SEIThickness(\n",
+    "    model.param, \"negative\", model.options\n",
+    ")\n",
     "model.submodels[\"Positive sei\"] = pybamm.sei.NoSEI(\n",
     "    model.param, \"positive\", model.options\n",
     ")\n",
+    "model.submodels[\"Positive sei thickness\"] = pybamm.sei.SEIThickness(\n",
+    "    model.param, \"positive\", model.options\n",
+    ")\n",
     "model.submodels[\"Negative sei on cracks\"] = pybamm.sei.NoSEI(\n",
     "    model.param, \"negative\", model.options, cracks=True\n",
     ")\n",
+    "model.submodels[\"Negative sei on cracks thickness\"] = pybamm.sei.SEIThickness(\n",
+    "    model.param, \"negative\", model.options, cracks=True\n",
+    ")\n",
     "model.submodels[\"Positive sei on cracks\"] = pybamm.sei.NoSEI(\n",
     "    model.param, \"positive\", model.options, cracks=True\n",
     ")\n",
+    "model.submodels[\"Positive sei on cracks thickness\"] = pybamm.sei.SEIThickness(\n",
+    "    model.param, \"positive\", model.options, cracks=True\n",
+    ")\n",
     "model.submodels[\"Negative lithium plating\"] = pybamm.lithium_plating.NoPlating(\n",
     "    model.param, \"Negative\"\n",
     ")\n",

examples/scripts/custom_model.py

@@ -75,9 +75,15 @@
     model.submodels[f"{domain} sei"] = pybamm.sei.NoSEI(
         model.param, domain, model.options
     )
+    model.submodels[f"{domain} sei thickness"] = pybamm.sei.SEIThickness(
+        model.param, domain, model.options
+    )
     model.submodels[f"{domain} sei on cracks"] = pybamm.sei.NoSEI(
         model.param, domain, model.options, cracks=True
     )
+    model.submodels[f"{domain} sei on cracks thickness"] = pybamm.sei.SEIThickness(
+        model.param, domain, model.options, cracks=True
+    )
     model.submodels[f"{domain} lithium plating"] = pybamm.lithium_plating.NoPlating(
         model.param, domain
     )

pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py

@@ -256,12 +256,6 @@ def set_open_circuit_potential_submodel(self):
 
     def set_sei_submodel(self):
         for domain in ["negative", "positive"]:
-            if self.options.electrode_types[domain] == "planar":
-                reaction_loc = "interface"
-            elif self.options["x-average side reactions"] == "true":
-                reaction_loc = "x-average"
-            else:
-                reaction_loc = "full electrode"
             sei_option = getattr(self.options, domain)["SEI"]
             phases = self.options.phases[domain]
             for phase in phases:
@@ -275,12 +269,16 @@ def set_sei_submodel(self):
                     submodel = pybamm.sei.SEIGrowth(
                         self.param,
                         domain,
-                        reaction_loc,
                         self.options,
                         phase,
                         cracks=False,
                     )
                 self.submodels[f"{domain} {phase} sei"] = submodel
+                self.submodels[
+                    f"{domain} {phase} sei thickness"
+                ] = pybamm.sei.SEIThickness(
+                    self.param, domain, self.options, phase, cracks=False
+                )
             if len(phases) > 1:
                 self.submodels[f"{domain} total sei"] = pybamm.sei.TotalSEI(
                     self.param, domain, self.options
@@ -304,19 +302,23 @@ def set_sei_on_cracks_submodel(self):
                             self.param, domain, self.options, phase, cracks=True
                         )
                     else:
-                        if self.options["x-average side reactions"] == "true":
-                            reaction_loc = "x-average"
-                        else:
-                            reaction_loc = "full electrode"
                         submodel = pybamm.sei.SEIGrowth(
                             self.param,
                             domain,
-                            reaction_loc,
                             self.options,
                             phase,
                             cracks=True,
                         )
                     self.submodels[f"{domain} {phase} sei on cracks"] = submodel
+                    self.submodels[
+                        f"{domain} {phase} sei on cracks thickness"
+                    ] = pybamm.sei.SEIThickness(
+                        self.param,
+                        domain,
+                        self.options,
+                        phase,
+                        cracks=True,
+                    )
                 if len(phases) > 1:
                     self.submodels[
                         f"{domain} total sei on cracks"

pybamm/models/submodels/active_material/loss_active_material.py

@@ -153,12 +153,12 @@ def set_rhs(self, variables):
             f"in {domain} electrode [mol]"
         ]
         # Multiply by mol.m-3 * m3 to get mol
-        c_s_av = variables[f"Average {domain} particle concentration [mol.m-3]"]
+        c_s_rav = variables[f"R-averaged {domain} particle concentration [mol.m-3]"]
         V = self.domain_param.L * self.param.A_cc
 
         self.rhs = {
             # minus sign because eps_solid is decreasing and LLI measures positive
-            lli_due_to_lam: -c_s_av * V * pybamm.x_average(deps_solid_dt),
+            lli_due_to_lam: -V * pybamm.x_average(c_s_rav * deps_solid_dt),
             eps_solid: deps_solid_dt,
         }
 

pybamm/models/submodels/interface/sei/__init__.py

@@ -3,3 +3,4 @@
 from .no_sei import NoSEI
 from .constant_sei import ConstantSEI
 from .sei_growth import SEIGrowth
+from .sei_thickness import SEIThickness