From afbea60eb31c395e3d6cafa32a647bcc10b6fb83 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Mon, 1 Aug 2022 12:10:13 -0400 Subject: [PATCH 1/8] #1986 remove zero electrolyte concentration cut-off, add some min/max to concentrations in functions --- .../lead_acid/basic_full.py | 4 +-- .../lithium_ion/basic_dfn.py | 28 ++----------------- .../lithium_ion/basic_dfn_half_cell.py | 2 +- .../base_electrolyte_conductivity.py | 15 ++-------- .../composite_conductivity.py | 1 + .../full_conductivity.py | 3 +- .../integrated_conductivity.py | 1 + .../full_surface_form_conductivity.py | 16 ++++------- .../base_electrolyte_diffusion.py | 10 ------- .../submodels/interface/base_interface.py | 4 --- pybamm/parameters/lead_acid_parameters.py | 7 +++++ pybamm/parameters/lithium_ion_parameters.py | 26 +++++++++++++---- pybamm/solvers/casadi_solver.py | 10 ++----- .../test_simulation_with_experiment.py | 3 ++ 14 files changed, 51 insertions(+), 79 deletions(-) diff --git a/pybamm/models/full_battery_models/lead_acid/basic_full.py b/pybamm/models/full_battery_models/lead_acid/basic_full.py index 60e1fe6ade..d5df7e17ef 100644 --- a/pybamm/models/full_battery_models/lead_acid/basic_full.py +++ b/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] = { diff --git a/pybamm/models/full_battery_models/lithium_ion/basic_dfn.py b/pybamm/models/full_battery_models/lithium_ion/basic_dfn.py index 5c6210bbb2..6133df35fb 100644 --- a/pybamm/models/full_battery_models/lithium_ion/basic_dfn.py +++ b/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), diff --git a/pybamm/models/full_battery_models/lithium_ion/basic_dfn_half_cell.py b/pybamm/models/full_battery_models/lithium_ion/basic_dfn_half_cell.py index 95ad8abf4a..fe1dc113de 100644 --- a/pybamm/models/full_battery_models/lithium_ion/basic_dfn_half_cell.py +++ b/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 diff --git a/pybamm/models/submodels/electrolyte_conductivity/base_electrolyte_conductivity.py b/pybamm/models/submodels/electrolyte_conductivity/base_electrolyte_conductivity.py index 483f274472..da27286549 100644 --- a/pybamm/models/submodels/electrolyte_conductivity/base_electrolyte_conductivity.py +++ b/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, ) diff --git a/pybamm/models/submodels/electrolyte_conductivity/composite_conductivity.py b/pybamm/models/submodels/electrolyte_conductivity/composite_conductivity.py index e433572d95..85b2e8c7d1 100644 --- a/pybamm/models/submodels/electrolyte_conductivity/composite_conductivity.py +++ b/pybamm/models/submodels/electrolyte_conductivity/composite_conductivity.py @@ -32,6 +32,7 @@ def __init__( def _higher_order_macinnes_function(self, x): "Function to differentiate between composite and first-order models" if self.higher_order_terms == "composite": + x = pybamm.maximum(x, 1e-15) return pybamm.log(x) elif self.higher_order_terms == "first-order": return x diff --git a/pybamm/models/submodels/electrolyte_conductivity/full_conductivity.py b/pybamm/models/submodels/electrolyte_conductivity/full_conductivity.py index 77058cc239..11d3eea25f 100644 --- a/pybamm/models/submodels/electrolyte_conductivity/full_conductivity.py +++ b/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 diff --git a/pybamm/models/submodels/electrolyte_conductivity/integrated_conductivity.py b/pybamm/models/submodels/electrolyte_conductivity/integrated_conductivity.py index 3819b4d7c2..e97843e3f7 100644 --- a/pybamm/models/submodels/electrolyte_conductivity/integrated_conductivity.py +++ b/pybamm/models/submodels/electrolyte_conductivity/integrated_conductivity.py @@ -34,6 +34,7 @@ def __init__(self, param, domain=None, options=None): pybamm.citations.register("BrosaPlanella2021") def _higher_order_macinnes_function(self, x): + x = pybamm.maximum(x, 1e-15) return pybamm.log(x) def get_coupled_variables(self, variables): diff --git a/pybamm/models/submodels/electrolyte_conductivity/surface_potential_form/full_surface_form_conductivity.py b/pybamm/models/submodels/electrolyte_conductivity/surface_potential_form/full_surface_form_conductivity.py index 34c96d4d79..d2780c8dc0 100644 --- a/pybamm/models/submodels/electrolyte_conductivity/surface_potential_form/full_surface_form_conductivity.py +++ b/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") diff --git a/pybamm/models/submodels/electrolyte_diffusion/base_electrolyte_diffusion.py b/pybamm/models/submodels/electrolyte_diffusion/base_electrolyte_diffusion.py index 75a3952ebd..ac8a9a7e19 100644 --- a/pybamm/models/submodels/electrolyte_diffusion/base_electrolyte_diffusion.py +++ b/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, - ) - ) diff --git a/pybamm/models/submodels/interface/base_interface.py b/pybamm/models/submodels/interface/base_interface.py index 7f3ac5fd62..e27b58b2e4 100644 --- a/pybamm/models/submodels/interface/base_interface.py +++ b/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) diff --git a/pybamm/parameters/lead_acid_parameters.py b/pybamm/parameters/lead_acid_parameters.py index a244dbaa97..1a8d2c8291 100644 --- a/pybamm/parameters/lead_acid_parameters.py +++ b/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 ( diff --git a/pybamm/parameters/lithium_ion_parameters.py b/pybamm/parameters/lithium_ion_parameters.py index 40d6b6d0cc..269713f854 100644 --- a/pybamm/parameters/lithium_ion_parameters.py +++ b/pybamm/parameters/lithium_ion_parameters.py @@ -199,11 +199,13 @@ def _set_dimensional_parameters(self): def D_e_dimensional(self, c_e, T): """Dimensional diffusivity in electrolyte""" + c_e = pybamm.maximum(c_e, 10) 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""" + c_e = pybamm.maximum(c_e, 10) inputs = {"Electrolyte concentration [mol.m-3]": c_e, "Temperature [K]": T} return pybamm.FunctionParameter("Electrolyte conductivity [S.m-1]", inputs) @@ -261,7 +263,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 +437,14 @@ 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 + """ + c_e = pybamm.maximum(c_e, 1e-2) + 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 +470,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 @@ -683,6 +693,9 @@ def D_dimensional(self, sto, T): def j0_dimensional(self, c_e, c_s_surf, T): """Dimensional exchange-current density [A.m-2]""" + # tol = 1e-3 + # c_e = c_e # pybamm.maximum(c_e, tol) + # c_s_surf = pybamm.maximum(pybamm.minimum(c_s_surf, 1 - tol), tol) inputs = { "Electrolyte concentration [mol.m-3]": c_e, f"{self.domain} particle surface concentration [mol.m-3]": c_s_surf, @@ -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,9 @@ def D(self, c_s, T): def j0(self, c_e, c_s_surf, T): """Dimensionless exchange-current density""" + tol = 1e-8 + c_e = pybamm.maximum(c_e, tol) + 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 +921,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} diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py index a78d89070c..f441385950 100644 --- a/pybamm/solvers/casadi_solver.py +++ b/pybamm/solvers/casadi_solver.py @@ -548,13 +548,7 @@ def create_integrator(self, model, inputs, t_eval=None, use_event_switch=False): algebraic = model.casadi_algebraic # When not in DEBUG mode (level=10), suppress warnings from CasADi - if ( - pybamm.logger.getEffectiveLevel() == 10 - or pybamm.settings.debug_mode is True - ): - show_eval_warnings = True - else: - show_eval_warnings = False + show_eval_warnings = False options = { **self.extra_options_setup, @@ -695,9 +689,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( diff --git a/tests/unit/test_experiments/test_simulation_with_experiment.py b/tests/unit/test_experiments/test_simulation_with_experiment.py index a8a01d4c69..e201e49a60 100644 --- a/tests/unit/test_experiments/test_simulation_with_experiment.py +++ b/tests/unit/test_experiments/test_simulation_with_experiment.py @@ -203,10 +203,12 @@ def test_run_experiment_breaks_early_error(self): model = pybamm.lithium_ion.DFN() parameter_values = pybamm.ParameterValues("Chen2020") + solver = pybamm.CasadiSolver(max_step_decrease_count=2) sim = pybamm.Simulation( model, experiment=experiment, parameter_values=parameter_values, + solver=solver, ) sol = sim.solve() self.assertEqual(len(sol.cycles), 1) @@ -220,6 +222,7 @@ def test_run_experiment_breaks_early_error(self): model, experiment=experiment, parameter_values=parameter_values, + solver=solver, ) sol = sim.solve() self.assertEqual(len(sol.cycles), 1) From 75d42b1e319d75d53561081ef1edce9bb130f6e8 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Mon, 1 Aug 2022 15:12:34 -0400 Subject: [PATCH 2/8] #1986 smaller current for LFP benchmarks --- .../parameters/lithium_ion/cells/A123_Lain2019/parameters.csv | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pybamm/input/parameters/lithium_ion/cells/A123_Lain2019/parameters.csv b/pybamm/input/parameters/lithium_ion/cells/A123_Lain2019/parameters.csv index 3544fe954b..5b8d6189f5 100644 --- a/pybamm/input/parameters/lithium_ion/cells/A123_Lain2019/parameters.csv +++ b/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, From 4704d183db69a49c0fda8d43794a1218ac5237c4 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Mon, 1 Aug 2022 15:14:17 -0400 Subject: [PATCH 3/8] #1986 make sure solver fails when supposed to fail --- pybamm/solvers/casadi_solver.py | 5 +---- .../test_experiments/test_simulation_with_experiment.py | 7 ++++++- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py index f441385950..050c8166fc 100644 --- a/pybamm/solvers/casadi_solver.py +++ b/pybamm/solvers/casadi_solver.py @@ -547,14 +547,11 @@ def create_integrator(self, model, inputs, t_eval=None, use_event_switch=False): rhs = model.casadi_rhs algebraic = model.casadi_algebraic - # When not in DEBUG mode (level=10), suppress warnings from CasADi - show_eval_warnings = False - options = { + "show_eval_warnings": False, **self.extra_options_setup, "reltol": self.rtol, "abstol": self.atol, - "show_eval_warnings": show_eval_warnings, } # set up and solve diff --git a/tests/unit/test_experiments/test_simulation_with_experiment.py b/tests/unit/test_experiments/test_simulation_with_experiment.py index e201e49a60..c2d20fd956 100644 --- a/tests/unit/test_experiments/test_simulation_with_experiment.py +++ b/tests/unit/test_experiments/test_simulation_with_experiment.py @@ -198,7 +198,12 @@ def test_run_experiment_breaks_early_infeasible(self): def test_run_experiment_breaks_early_error(self): experiment = pybamm.Experiment( - [("Rest for 10 minutes", "Discharge at 10 C for 1 minute")] + [ + ( + "Rest for 10 minutes", + "Discharge at 20 C for 10 minutes (10 minute period)", + ) + ] ) model = pybamm.lithium_ion.DFN() From c2260143228199ef2962a7702d13183c830924e5 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Mon, 1 Aug 2022 18:14:25 -0400 Subject: [PATCH 4/8] #1986 update callback --- pybamm/callbacks.py | 3 ++- pybamm/simulation.py | 7 ++----- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/pybamm/callbacks.py b/pybamm/callbacks.py index d536a0aae0..12f9830f04 100644 --- a/pybamm/callbacks.py +++ b/pybamm/callbacks.py @@ -222,7 +222,8 @@ def on_experiment_end(self, logs): self.logger.notice("Finish experiment simulation, took {}".format(elapsed_time)) def on_experiment_error(self, logs): - pass + error = logs["error"] + self.logger.error(f"Simulation errored: {error}") def on_experiment_infeasible(self, logs): termination = logs["termination"] diff --git a/pybamm/simulation.py b/pybamm/simulation.py index 02d5d45650..2cde917937 100644 --- a/pybamm/simulation.py +++ b/pybamm/simulation.py @@ -837,6 +837,7 @@ def solve( except pybamm.SolverError as e: logs["error"] = e callbacks.on_experiment_error(logs) + step_solution.termination = f"error: {e}" feasible = False # If none of the cycles worked, raise an error if cycle_num == 1 and step_num == 1: @@ -858,6 +859,7 @@ def solve( or step_solution.termination == "final time" or "[experiment]" in step_solution.termination ): + callbacks.on_experiment_infeasible(logs) feasible = False break @@ -911,11 +913,6 @@ def solve( if min_voltage <= voltage_stop[0]: break - # Break if the experiment is infeasible (or errored) - if feasible is False: - callbacks.on_experiment_infeasible(logs) - break - if self.solution is not None and len(all_cycle_solutions) > 0: self.solution.cycles = all_cycle_solutions self.solution.set_summary_variables(all_summary_variables) From 50c7c20ccda3ee11f207ce267bb0f98f0896d426 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Thu, 4 Aug 2022 12:37:30 -0400 Subject: [PATCH 5/8] #1986 make tolerances a setting --- .../composite_conductivity.py | 3 ++- .../integrated_conductivity.py | 3 ++- pybamm/parameters/lithium_ion_parameters.py | 17 +++++++++-------- pybamm/settings.py | 9 +++++++++ 4 files changed, 22 insertions(+), 10 deletions(-) diff --git a/pybamm/models/submodels/electrolyte_conductivity/composite_conductivity.py b/pybamm/models/submodels/electrolyte_conductivity/composite_conductivity.py index 85b2e8c7d1..abd5afb658 100644 --- a/pybamm/models/submodels/electrolyte_conductivity/composite_conductivity.py +++ b/pybamm/models/submodels/electrolyte_conductivity/composite_conductivity.py @@ -32,7 +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": - x = pybamm.maximum(x, 1e-15) + 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 diff --git a/pybamm/models/submodels/electrolyte_conductivity/integrated_conductivity.py b/pybamm/models/submodels/electrolyte_conductivity/integrated_conductivity.py index e97843e3f7..d520dc9409 100644 --- a/pybamm/models/submodels/electrolyte_conductivity/integrated_conductivity.py +++ b/pybamm/models/submodels/electrolyte_conductivity/integrated_conductivity.py @@ -34,7 +34,8 @@ def __init__(self, param, domain=None, options=None): pybamm.citations.register("BrosaPlanella2021") def _higher_order_macinnes_function(self, x): - x = pybamm.maximum(x, 1e-15) + tol = pybamm.settings.tolerances["macinnes__c_e"] + x = pybamm.maximum(x, tol) return pybamm.log(x) def get_coupled_variables(self, variables): diff --git a/pybamm/parameters/lithium_ion_parameters.py b/pybamm/parameters/lithium_ion_parameters.py index 269713f854..25ad42813e 100644 --- a/pybamm/parameters/lithium_ion_parameters.py +++ b/pybamm/parameters/lithium_ion_parameters.py @@ -199,13 +199,15 @@ def _set_dimensional_parameters(self): def D_e_dimensional(self, c_e, T): """Dimensional diffusivity in electrolyte""" - c_e = pybamm.maximum(c_e, 10) + 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""" - c_e = pybamm.maximum(c_e, 10) + 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) @@ -442,7 +444,8 @@ def chiT_over_c(self, c_e, T): chi * (1 + Theta * T) / c, as it appears in the electrolyte potential equation """ - c_e = pybamm.maximum(c_e, 1e-2) + 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): @@ -693,9 +696,6 @@ def D_dimensional(self, sto, T): def j0_dimensional(self, c_e, c_s_surf, T): """Dimensional exchange-current density [A.m-2]""" - # tol = 1e-3 - # c_e = c_e # pybamm.maximum(c_e, tol) - # c_s_surf = pybamm.maximum(pybamm.minimum(c_s_surf, 1 - tol), tol) inputs = { "Electrolyte concentration [mol.m-3]": c_e, f"{self.domain} particle surface concentration [mol.m-3]": c_s_surf, @@ -711,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) @@ -879,8 +879,9 @@ def D(self, c_s, T): def j0(self, c_e, c_s_surf, T): """Dimensionless exchange-current density""" - tol = 1e-8 + 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 diff --git a/pybamm/settings.py b/pybamm/settings.py index c3fd20bf73..b213eaaf38 100644 --- a/pybamm/settings.py +++ b/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): From b02977701f13a9170ab5b6f9e7598665be968655 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Thu, 4 Aug 2022 19:17:57 -0400 Subject: [PATCH 6/8] #1986 remove duplicate-ish test --- pybamm/solvers/casadi_solver.py | 2 +- .../test_simulation_with_experiment.py | 14 -------------- 2 files changed, 1 insertion(+), 15 deletions(-) diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py index 34ddc7b1bf..7baaa8691c 100644 --- a/pybamm/solvers/casadi_solver.py +++ b/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 diff --git a/tests/unit/test_experiments/test_simulation_with_experiment.py b/tests/unit/test_experiments/test_simulation_with_experiment.py index c6021e5e6a..f5062a69af 100644 --- a/tests/unit/test_experiments/test_simulation_with_experiment.py +++ b/tests/unit/test_experiments/test_simulation_with_experiment.py @@ -239,20 +239,6 @@ def test_run_experiment_breaks_early_error(self): self.assertEqual(len(sol.cycles), 1) self.assertEqual(len(sol.cycles[0].steps), 1) - # Different experiment setup style - experiment = pybamm.Experiment( - ["Rest for 10 minutes", "Discharge at 10 C for 1 minute"] - ) - sim = pybamm.Simulation( - model, - experiment=experiment, - parameter_values=parameter_values, - solver=solver, - ) - sol = sim.solve() - self.assertEqual(len(sol.cycles), 1) - self.assertEqual(len(sol.cycles[0].steps), 1) - # Different callback - this is for coverage on the `Callback` class sol = sim.solve(callbacks=pybamm.callbacks.Callback()) From 66def4fa69858eb7c6c0acfeb78e47e2c7c43fb5 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Thu, 4 Aug 2022 20:34:27 -0400 Subject: [PATCH 7/8] #1986 coverage --- .../test_simulation_with_experiment.py | 17 +++++++++++++++++ tests/unit/test_solvers/test_solution.py | 8 ++++++-- 2 files changed, 23 insertions(+), 2 deletions(-) diff --git a/tests/unit/test_experiments/test_simulation_with_experiment.py b/tests/unit/test_experiments/test_simulation_with_experiment.py index f5062a69af..e33c027378 100644 --- a/tests/unit/test_experiments/test_simulation_with_experiment.py +++ b/tests/unit/test_experiments/test_simulation_with_experiment.py @@ -239,6 +239,23 @@ def test_run_experiment_breaks_early_error(self): self.assertEqual(len(sol.cycles), 1) self.assertEqual(len(sol.cycles[0].steps), 1) + # Different experiment setup style + experiment = pybamm.Experiment( + [ + "Rest for 10 minutes", + "Discharge at 20 C for 10 minutes (10 minute period)", + ] + ) + sim = pybamm.Simulation( + model, + experiment=experiment, + parameter_values=parameter_values, + solver=solver, + ) + sol = sim.solve() + self.assertEqual(len(sol.cycles), 1) + self.assertEqual(len(sol.cycles[0].steps), 1) + # Different callback - this is for coverage on the `Callback` class sol = sim.solve(callbacks=pybamm.callbacks.Callback()) diff --git a/tests/unit/test_solvers/test_solution.py b/tests/unit/test_solvers/test_solution.py index 9fea56d51b..23c78f44a3 100644 --- a/tests/unit/test_solvers/test_solution.py +++ b/tests/unit/test_solvers/test_solution.py @@ -90,10 +90,14 @@ def test_add_solutions(self): sol3 = pybamm.Solution(t3, y3, pybamm.BaseModel(), {"a": 3}) self.assertEqual((sol_sum + sol3).all_ts, sol_sum.copy().all_ts) - # radd - sol4 = None + sol3 + # add None + sol4 = sol3 + None self.assertEqual(sol3.all_ys, sol4.all_ys) + # radd + sol5 = None + sol3 + self.assertEqual(sol3.all_ys, sol5.all_ys) + # radd failure with self.assertRaisesRegex( pybamm.SolverError, "Only a Solution or None can be added to a Solution" From 36870f3d405ec6069f4589a8338e7be166f92434 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Fri, 5 Aug 2022 09:43:12 -0400 Subject: [PATCH 8/8] #1986 changelog --- CHANGELOG.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2c890713a1..c247426f61 100644 --- a/CHANGELOG.md +++ b/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