Merge pull request #2207 from pybamm-team/issue-1986-electrolyte-cut-off

CHANGELOG.md

@@ -1,5 +1,9 @@
 # [Unreleased](https://github.com/pybamm-team/PyBaMM/)
 
+## Optimizations
+
+-   Added limits for variables in some functions to avoid division by zero, sqrt(negative number), etc ([#2213](https://github.com/pybamm-team/PyBaMM/pull/2213))
+
 # [v22.7](https://github.com/pybamm-team/PyBaMM/tree/v22.7) - 2022-07-31
 
 ## Features 

pybamm/input/parameters/lithium_ion/cells/A123_Lain2019/parameters.csv

@@ -34,5 +34,5 @@ Positive current collector thermal conductivity [W.m-1.K-1],237,CRC Handbook,alu
 ,,,
 # Electrical,,,
 Nominal cell capacity [A.h],1.1,Lain 2019,
-Current function [A],4.4,,
-Typical current [A],4.4,Lain 2019,
+Current function [A],1.1,,
+Typical current [A],1.1,Lain 2019,

pybamm/models/full_battery_models/lead_acid/basic_full.py

@@ -103,7 +103,7 @@ def __init__(self, name="Basic full model"):
 
         # transport_efficiency
         tor = pybamm.concatenation(
-            eps_n**param.n.b_e, eps_s**param.s.b_e, eps_p**param.p.b_e
+            eps_n ** param.n.b_e, eps_s ** param.s.b_e, eps_p ** param.p.b_e
         )
 
         # Interfacial reactions
@@ -177,7 +177,7 @@ def __init__(self, name="Basic full model"):
         # Current in the electrolyte
         ######################
         i_e = (param.kappa_e(c_e, T) * tor * param.gamma_e / param.C_e) * (
-            param.chi(c_e, T) * pybamm.grad(c_e) / c_e - pybamm.grad(phi_e)
+            param.chiT_over_c(c_e, T) * pybamm.grad(c_e) - pybamm.grad(phi_e)
         )
         self.algebraic[phi_e] = pybamm.div(i_e) - j
         self.boundary_conditions[phi_e] = {

pybamm/models/full_battery_models/lithium_ion/basic_dfn.py

@@ -188,25 +188,6 @@ def __init__(self, name="Doyle-Fuller-Newman model"):
         }
         self.initial_conditions[c_s_n] = param.n.c_init
         self.initial_conditions[c_s_p] = param.p.c_init
-        # Events specify points at which a solution should terminate
-        self.events += [
-            pybamm.Event(
-                "Minimum negative particle surface concentration",
-                pybamm.min(c_s_surf_n) - 0.01,
-            ),
-            pybamm.Event(
-                "Maximum negative particle surface concentration",
-                (1 - 0.01) - pybamm.max(c_s_surf_n),
-            ),
-            pybamm.Event(
-                "Minimum positive particle surface concentration",
-                pybamm.min(c_s_surf_p) - 0.01,
-            ),
-            pybamm.Event(
-                "Maximum positive particle surface concentration",
-                (1 - 0.01) - pybamm.max(c_s_surf_p),
-            ),
-        ]
         ######################
         # Current in the solid
         ######################
@@ -238,7 +219,7 @@ def __init__(self, name="Doyle-Fuller-Newman model"):
         # Current in the electrolyte
         ######################
         i_e = (param.kappa_e(c_e, T) * tor * param.gamma_e / param.C_e) * (
-            param.chi(c_e, T) * pybamm.grad(c_e) / c_e - pybamm.grad(phi_e)
+            param.chiT_over_c(c_e, T) * pybamm.grad(c_e) - pybamm.grad(phi_e)
         )
         self.algebraic[phi_e] = pybamm.div(i_e) - j
         self.boundary_conditions[phi_e] = {
@@ -260,11 +241,6 @@ def __init__(self, name="Doyle-Fuller-Newman model"):
             "right": (pybamm.Scalar(0), "Neumann"),
         }
         self.initial_conditions[c_e] = param.c_e_init
-        self.events.append(
-            pybamm.Event(
-                "Zero electrolyte concentration cut-off", pybamm.min(c_e) - 0.002
-            )
-        )
 
         ######################
         # (Some) variables
@@ -281,7 +257,9 @@ def __init__(self, name="Doyle-Fuller-Newman model"):
             "Electrolyte potential": phi_e,
             "Positive electrode potential": phi_s_p,
             "Terminal voltage": voltage,
+            "Terminal voltage [V]": voltage * param.potential_scale + param.ocv_ref,
         }
+        # Events specify points at which a solution should terminate
         self.events += [
             pybamm.Event("Minimum voltage", voltage - param.voltage_low_cut),
             pybamm.Event("Maximum voltage", voltage - param.voltage_high_cut),

pybamm/models/full_battery_models/lithium_ion/basic_dfn_half_cell.py

@@ -248,7 +248,7 @@ def __init__(self, options=None, name="Doyle-Fuller-Newman half cell model"):
         # Current in the electrolyte
         ######################
         i_e = (param.kappa_e(c_e, T) * tor * gamma_e / param.C_e) * (
-            param.chi(c_e, T) * pybamm.grad(c_e) / c_e - pybamm.grad(phi_e)
+            param.chiT_over_c(c_e, T) * pybamm.grad(c_e) - pybamm.grad(phi_e)
         )
         self.algebraic[phi_e] = pybamm.div(i_e) - j
 

pybamm/models/submodels/electrolyte_conductivity/base_electrolyte_conductivity.py

@@ -267,10 +267,7 @@ def _get_electrolyte_overpotentials(self, variables):
             c_e_n = variables["Negative electrolyte concentration"]
             T_n = variables["Negative electrode temperature"]
             indef_integral_n = pybamm.IndefiniteIntegral(
-                param.chi(c_e_n, T_n)
-                * (1 + param.Theta * T_n)
-                * pybamm.grad(c_e_n)
-                / c_e_n,
+                param.chiT_over_c(c_e_n, T_n) * pybamm.grad(c_e_n),
                 pybamm.standard_spatial_vars.x_n,
             )
 
@@ -284,17 +281,11 @@ def _get_electrolyte_overpotentials(self, variables):
 
         # concentration overpotential
         indef_integral_s = pybamm.IndefiniteIntegral(
-            param.chi(c_e_s, T_s)
-            * (1 + param.Theta * T_s)
-            * pybamm.grad(c_e_s)
-            / c_e_s,
+            param.chiT_over_c(c_e_s, T_s) * pybamm.grad(c_e_s),
             pybamm.standard_spatial_vars.x_s,
         )
         indef_integral_p = pybamm.IndefiniteIntegral(
-            param.chi(c_e_p, T_p)
-            * (1 + param.Theta * T_p)
-            * pybamm.grad(c_e_p)
-            / c_e_p,
+            param.chiT_over_c(c_e_p, T_p) * pybamm.grad(c_e_p),
             pybamm.standard_spatial_vars.x_p,
         )
 

pybamm/models/submodels/electrolyte_conductivity/composite_conductivity.py

@@ -32,6 +32,8 @@ def __init__(
     def _higher_order_macinnes_function(self, x):
         "Function to differentiate between composite and first-order models"
         if self.higher_order_terms == "composite":
+            tol = pybamm.settings.tolerances["macinnes__c_e"]
+            x = pybamm.maximum(x, tol)
             return pybamm.log(x)
         elif self.higher_order_terms == "first-order":
             return x

pybamm/models/submodels/electrolyte_conductivity/full_conductivity.py

@@ -43,8 +43,7 @@ def get_coupled_variables(self, variables):
         phi_e = variables["Electrolyte potential"]
 
         i_e = (param.kappa_e(c_e, T) * tor * param.gamma_e / param.C_e) * (
-            param.chi(c_e, T) * (1 + param.Theta * T) * pybamm.grad(c_e) / c_e
-            - pybamm.grad(phi_e)
+            param.chiT_over_c(c_e, T) * pybamm.grad(c_e) - pybamm.grad(phi_e)
         )
 
         # Override print_name

pybamm/models/submodels/electrolyte_conductivity/integrated_conductivity.py

@@ -34,6 +34,8 @@ def __init__(self, param, domain=None, options=None):
         pybamm.citations.register("BrosaPlanella2021")
 
     def _higher_order_macinnes_function(self, x):
+        tol = pybamm.settings.tolerances["macinnes__c_e"]
+        x = pybamm.maximum(x, tol)
         return pybamm.log(x)
 
     def get_coupled_variables(self, variables):

pybamm/models/submodels/electrolyte_conductivity/surface_potential_form/full_surface_form_conductivity.py

@@ -57,7 +57,7 @@ def get_coupled_variables(self, variables):
             T = variables[self.domain + " electrode temperature"]
 
             i_e = conductivity * (
-                ((1 + param.Theta * T) * param.chi(c_e, T) / c_e) * pybamm.grad(c_e)
+                param.chiT_over_c(c_e, T) * pybamm.grad(c_e)
                 + pybamm.grad(delta_phi)
                 + i_boundary_cc / sigma_eff
             )
@@ -79,11 +79,11 @@ def get_coupled_variables(self, variables):
             tor_s = variables["Separator porosity"]
             T = variables["Separator temperature"]
 
-            chi_e_s = param.chi(c_e_s, T)
+            chiT_over_c_e_s = param.chiT_over_c(c_e_s, T)
             kappa_s_eff = param.kappa_e(c_e_s, T) * tor_s
 
             phi_e = phi_e_n_s + pybamm.IndefiniteIntegral(
-                (1 + param.Theta * T) * chi_e_s / c_e_s * pybamm.grad(c_e_s)
+                chiT_over_c_e_s * pybamm.grad(c_e_s)
                 - param.C_e * i_boundary_cc / kappa_s_eff,
                 x_s,
             )
@@ -161,10 +161,7 @@ def set_boundary_conditions(self, variables):
             flux_left = -i_boundary_cc * pybamm.BoundaryValue(1 / sigma_eff, "left")
             flux_right = (
                 (i_boundary_cc / pybamm.BoundaryValue(conductivity, "right"))
-                - pybamm.BoundaryValue(
-                    (1 + param.Theta * T) * param.chi(c_e, T) / c_e, "right"
-                )
-                * c_e_flux
+                - pybamm.BoundaryValue(param.chiT_over_c(c_e, T), "right") * c_e_flux
                 - i_boundary_cc * pybamm.BoundaryValue(1 / sigma_eff, "right")
             )
 
@@ -178,10 +175,7 @@ def set_boundary_conditions(self, variables):
             c_e_flux = pybamm.BoundaryGradient(c_e, "left")
             flux_left = (
                 (i_boundary_cc / pybamm.BoundaryValue(conductivity, "left"))
-                - pybamm.BoundaryValue(
-                    (1 + param.Theta * T) * param.chi(c_e, T) / c_e, "left"
-                )
-                * c_e_flux
+                - pybamm.BoundaryValue(param.chiT_over_c(c_e, T), "left") * c_e_flux
                 - i_boundary_cc * pybamm.BoundaryValue(1 / sigma_eff, "left")
             )
             flux_right = -i_boundary_cc * pybamm.BoundaryValue(1 / sigma_eff, "right")

pybamm/models/submodels/electrolyte_diffusion/base_electrolyte_diffusion.py

@@ -156,13 +156,3 @@ def _get_total_concentration_electrolyte(self, eps_c_e):
         }
 
         return variables
-
-    def set_events(self, variables):
-        c_e = variables["Electrolyte concentration"]
-        self.events.append(
-            pybamm.Event(
-                "Zero electrolyte concentration cut-off",
-                pybamm.min(c_e) - 0.002,
-                pybamm.EventType.TERMINATION,
-            )
-        )

pybamm/models/submodels/interface/base_interface.py

@@ -106,10 +106,6 @@ def _get_exchange_current_density(self, variables):
                     c_e = c_e.orphans[0]
                     T = T.orphans[0]
 
-            tol = 1e-8
-            c_e = pybamm.maximum(tol, c_e)
-            c_s_surf = pybamm.maximum(tol, pybamm.minimum(c_s_surf, 1 - tol))
-
             j0 = (
                 domain_param.gamma
                 * domain_param.j0(c_e, c_s_surf, T)

pybamm/parameters/lead_acid_parameters.py

@@ -422,6 +422,13 @@ def kappa_e(self, c_e, T):
         kappa_scale = self.F ** 2 * self.D_e_typ * self.c_e_typ / (self.R * self.T_ref)
         return self.kappa_e_dimensional(c_e_dimensional, self.T_ref) / kappa_scale
 
+    def chiT_over_c(self, c_e, T):
+        """
+        chi * (1 + Theta * T) / c,
+        as it appears in the electrolyte potential equation
+        """
+        return self.chi(c_e, T) * (1 + self.Theta * T) / c_e
+
     def chi(self, c_e, T, c_ox=0, c_hy=0):
         """Thermodynamic factor"""
         return (

pybamm/parameters/lithium_ion_parameters.py

@@ -199,11 +199,15 @@ def _set_dimensional_parameters(self):
 
     def D_e_dimensional(self, c_e, T):
         """Dimensional diffusivity in electrolyte"""
+        tol = pybamm.settings.tolerances["D_e__c_e"]
+        c_e = pybamm.maximum(c_e, tol)
         inputs = {"Electrolyte concentration [mol.m-3]": c_e, "Temperature [K]": T}
         return pybamm.FunctionParameter("Electrolyte diffusivity [m2.s-1]", inputs)
 
     def kappa_e_dimensional(self, c_e, T):
         """Dimensional electrolyte conductivity"""
+        tol = pybamm.settings.tolerances["D_e__c_e"]
+        c_e = pybamm.maximum(c_e, tol)
         inputs = {"Electrolyte concentration [mol.m-3]": c_e, "Temperature [K]": T}
         return pybamm.FunctionParameter("Electrolyte conductivity [S.m-1]", inputs)
 
@@ -261,7 +265,7 @@ def _set_scales(self):
 
         # Electrolyte diffusion timescale
         self.D_e_typ = self.D_e_dimensional(self.c_e_typ, self.T_ref)
-        self.tau_diffusion_e = self.L_x**2 / self.D_e_typ
+        self.tau_diffusion_e = self.L_x ** 2 / self.D_e_typ
 
         # Thermal diffusion timescale
         self.tau_th_yz = self.therm.tau_th_yz
@@ -435,6 +439,15 @@ def chi(self, c_e, T):
         """
         return (2 * (1 - self.t_plus(c_e, T))) * (self.one_plus_dlnf_dlnc(c_e, T))
 
+    def chiT_over_c(self, c_e, T):
+        """
+        chi * (1 + Theta * T) / c,
+        as it appears in the electrolyte potential equation
+        """
+        tol = pybamm.settings.tolerances["chi__c_e"]
+        c_e = pybamm.maximum(c_e, tol)
+        return self.chi(c_e, T) * (1 + self.Theta * T) / c_e
+
     def t_plus(self, c_e, T):
         """Cation transference number (dimensionless)"""
         inputs = {
@@ -460,7 +473,7 @@ def D_e(self, c_e, T):
     def kappa_e(self, c_e, T):
         """Dimensionless electrolyte conductivity"""
         c_e_dimensional = c_e * self.c_e_typ
-        kappa_scale = self.F**2 * self.D_e_typ * self.c_e_typ / (self.R * self.T_ref)
+        kappa_scale = self.F ** 2 * self.D_e_typ * self.c_e_typ / (self.R * self.T_ref)
         T_dim = self.Delta_T * T + self.T_ref
         return self.kappa_e_dimensional(c_e_dimensional, T_dim) / kappa_scale
 
@@ -698,7 +711,7 @@ def U_dimensional(self, sto, T):
         # bound stoichiometry between tol and 1-tol. Adding 1/sto + 1/(sto-1) later
         # will ensure that ocp goes to +- infinity if sto goes into that region
         # anyway
-        tol = 1e-10
+        tol = pybamm.settings.tolerances["U__c_s"]
         sto = pybamm.maximum(pybamm.minimum(sto, 1 - tol), tol)
         inputs = {f"{self.domain} particle stoichiometry": sto}
         u_ref = pybamm.FunctionParameter(f"{self.domain} electrode OCP [V]", inputs)
@@ -757,7 +770,7 @@ def _set_scales(self):
         self.tau_r = main.F * self.c_max / (self.j0_ref_dimensional * self.a_typ)
         # Particle diffusion timescales
         self.D_typ_dim = self.D_dimensional(pybamm.Scalar(1), main.T_ref)
-        self.tau_diffusion = self.R_typ**2 / self.D_typ_dim
+        self.tau_diffusion = self.R_typ ** 2 / self.D_typ_dim
 
     def _set_dimensionless_parameters(self):
         main = self.main_param
@@ -807,7 +820,7 @@ def _set_dimensionless_parameters(self):
         self.sigma_cc = (
             self.sigma_cc_dimensional * main.potential_scale / main.i_typ / main.L_x
         )
-        self.sigma_cc_prime = self.sigma_cc * main.delta**2
+        self.sigma_cc_prime = self.sigma_cc * main.delta ** 2
         self.sigma_cc_dbl_prime = self.sigma_cc_prime * main.delta
 
         # Electrolyte Properties
@@ -866,6 +879,10 @@ def D(self, c_s, T):
 
     def j0(self, c_e, c_s_surf, T):
         """Dimensionless exchange-current density"""
+        tol = pybamm.settings.tolerances["j0__c_e"]
+        c_e = pybamm.maximum(c_e, tol)
+        tol = pybamm.settings.tolerances["j0__c_s"]
+        c_s_surf = pybamm.maximum(pybamm.minimum(c_s_surf, 1 - tol), tol)
         c_e_dim = c_e * self.main_param.c_e_typ
         c_s_surf_dim = c_s_surf * self.c_max
         T_dim = self.main_param.Delta_T * T + self.main_param.T_ref
@@ -905,7 +922,7 @@ def k_cr(self, T):
         Dimensionless cracking rate for the electrode;
         """
         T_dim = self.main_param.Delta_T * T + self.main_param.T_ref
-        delta_k_cr = self.E**self.m_cr * self.l_cr_0 ** (self.m_cr / 2 - 1)
+        delta_k_cr = self.E ** self.m_cr * self.l_cr_0 ** (self.m_cr / 2 - 1)
         return (
             pybamm.FunctionParameter(
                 f"{self.domain} electrode cracking rate", {"Temperature [K]": T_dim}

pybamm/settings.py

@@ -12,6 +12,15 @@ class Settings(object):
     _abs_smoothing = "exact"
     max_words_in_line = 4
     max_y_value = 1e5
+    tolerances = {
+        "D_e__c_e": 10,  # dimensional
+        "kappa_e__c_e": 10,  # dimensional
+        "chi__c_e": 1e-2,  # dimensionless
+        "U__c_s": 1e-10,  # dimensional
+        "j0__c_e": 1e-8,  # dimensionless
+        "j0__c_s": 1e-8,  # dimensionless
+        "macinnes__c_e": 1e-15,  # dimensionless
+    }
 
     @property
     def debug_mode(self):

pybamm/solvers/casadi_solver.py

@@ -267,7 +267,7 @@ def _integrate(self, model, t_eval, inputs_dict=None):
                             f"Maximum number of decreased steps occurred at t={t_dim}. "
                             "Try solving the model up to this time only or reducing "
                             f"dt_max (currently, dt_max={dt_max_dim}) and/or reducing "
-                            "the size of the time steps or period of the expeeriment."
+                            "the size of the time steps or period of the experiment."
                         )
                 # Check if the sign of an event changes, if so find an accurate
                 # termination point and exit
@@ -686,9 +686,11 @@ def _run_integrator(
                 inputs_with_tmin = casadi.vertcat(inputs, t_min)
                 # Call the integrator once, with the grid
                 timer = pybamm.Timer()
+                pybamm.logger.debug("Calling casadi integrator")
                 casadi_sol = integrator(
                     x0=y0_diff, z0=y0_alg, p=inputs_with_tmin, **self.extra_options_call
                 )
+                pybamm.logger.debug("Finished casadi integrator")
                 integration_time = timer.time()
                 y_sol = casadi.vertcat(casadi_sol["xf"], casadi_sol["zf"])
                 sol = pybamm.Solution(